Python numpy multidimensional arrays with slicing in C99 (single header)

Python numpy multidimensional arrays with slicing in C99 (single header)

9 comments · 87% upvoted · source · other discussions

comments

operamint [OP]

Some notes:

Spans are views only. E.g transposing a 2D span does not change the underlying data array.

This header should cover many needs when dealing with dynamic multidimensional arrays in C as it is very fast and compact, but the underlying data must be allocated and freed independently, and care must be taken that spans does not outlive the lifetime of the underlying array.

Add: To understand how single- and multi-dimensional array slicing works, here is a short tutorial. The behaviour in cspan is pretty much identical apart from the syntax. Negative indices are not supported.

​

Slicing Expression Python numpy C cspan
Select a single row / column / element. Reduces output span rank. a2[m][:] cspan_slice(S, &a2, {m}, {c_ALL})
Select elements with an index starting at m and ending at n-1. a[m:n] cspan_slice(S, &a, {m,n})
Select all elements in a given axis a[:] or a[0:-1] cspan_slice(S, &a, {c_ALL})
Select elements starting with index 0 and up to element with index n-1 a[:n] cspan_slice(S, &a, {0,n})
Select elements starting with index m and up to the last element a[m:] cspan_slice(S, &a, {m,c_END})
Select elements with index m through n (exclusive), with an increment k a[m:n:k] cspan_slice(S, &a, {m,n,k})
dainasol

This is "just" the array, right? (Not to undervalue the project by any means but to clarify)

I've been wondering if the C code of numpy could be taken out of the project and use that directly, but numpy is like 2/3 python and 1/3 C so idk if that's viable

operamint [OP]

Ah, I get it. You mean all the support functions associated, like max() mean() stddev() etc... Sure, but there are many libraries which have these.

What cspan brings to the table is very good ergonomics combined with type safety and a pretty much optimal implementation of handling dynamic multidimensional array views, slicing and iteration.

dainasol

Understood. But I think numpy does bring a lot of value besides the array because every function relies on SIMD and they handle all the intricacies to make it cross platform. I don't know how difficult this is though (I'm a noob in C) but I think SIMD varies from one chip to another (?).

EDIT: All functions is an exaggeration because things like cumsum don't parallelize

nerd4code

Pretty much every aspect of a chip can vary, and newer things vary more. SIMD support in a processor varies from none at all, to ugly (e.g., MMX—register overlays were de rigeur for a bit because it saved on die space & connectivity), to kinda okay, to all the thing is capable of (e.g., GPUs).

Styles of expansion and fractioning vary (e.g., x86 supports fixed-width vectors but AVX supports masking and “swizzling” more-or-less on-the-fly; other arches support counted ops or special encoding structures for fast loops), operations vary (e.g., AVX-pre-512 can do full-width packed byte ops, but most AVX512 extensions only support byte access to memory by narrowing from/widening to 32-bit integers in-register), data types operated on vary (integer, floating-point, fixed-point, _Accum, bytestrings, or bitstrings), vector width and pipeline/bus throughput vary.

SIMD impls exist with 2D and 3D register sets, enabling tiled matrix, image, or voxel ops to be performed internally. Most use 1D registers.

Register count and sharing vary. Some vector ops will downmodulate your CPU’s clock rate to prevent overheating or undeserved computational efficiency, so you have to be burning through RAM to make use of them worthwhile. Some vector ops will interfere with floating-point or ALU usage. Some vector ops are implemented by ganged coprocessors, others by separate VPUs.

You can even do SIMD shit in integer registers to some extent, which is honestly not that hard for ASCII since holding bit 7 at 0 is a form of carry inhibition (which is offered directly by some, not by most even though it’s a super-cheap way to get packed integer add/sub). Pipelining and superscalar trickery can give you enough low-latency ILP that vector ops might not always be faster, or maybe using a mix of SIMD and non-SIMD insns is best.

All of these processors and techniques are usable from C by a number of means. If you’re feeding more than choppy little operator functions to the compiler, a modern optimizer (compiler-specific) will generally at least take a swipe at SIMDizing. You can force it to try with pragmas (of many fascinating, highly-touchy varieties like OpenMP, optimizer hints, or offloading suggestions) or attributes (highly target-, compiler-, and version-specific), built-in functions (ditto), library functions (ditto, +runtime-specific +vendor-specific), or inline assembler (oddly less compiler-/version-specific in the GNUish subset of compilers). There are also offloading frameworks for purer SIMD (most things that can run kernels on a GPU can potentially be used), and SIMD techniques can be applied across networks as well.

So a big part of something like PyTorch (which is mostly a wretchedly huge amalgamation of secondary libraries) is an implementation of every sufficiently-fine-grained operation that might need to happen in bulk, for most modern compilers, processors, etc., and a bunch of glue so it can translate the equivalent of a Python-embedded DSL form of a program specification into orchestration of all the necessary routines on all necessary devices—in sum it’s a compiler-interpreter-runtime via Python frontend.

It’s not dissimilar to the gulch that static dataflow projects have ended up in since time immemorial, but the algorithms are regular and bulky enough, and data scientists blithe enough, that it actually works in its setting.

operamint [OP]

I don't get what you mean by "just" the array? This literally does the same as the numpy C code, except managing the underlaying data array. You create a multi-dimensional array view using a C array as "source", and can access it.

#include <stdio.h>
#include "cspan.h"
using_cspan3(Span, float);
...
float data[3*5*6];
for (int i=0; i<3*5*6; ++i) data[i] = i + 1;

Span3 md = cspan_md(data, 3, 5, 6); // md[3][5][6] view of data
Span2 mat = cspan_submd3(&md, 0); // md[0] => mat 5x6 matrix

cspan_transpose(&mat); // switch row/col view => 6x5

// create a view of a sub-3x3 matrix of mat
mat = cspan_slice(Span2, &mat, {2,5}, {2,5});

for (int i=0; i<mat.shape[0]; ++i) {
  for (int j=0; j<mat.shape[1]; ++j)
    printf(" %g", *cspan_at(&mat, i,j));
  puts("");
}
// or iterate it flat
puts(""); c_foreach (i, Span2, mat) printf(" %g", *i.ref);
puts("");

Outputs:

 15 21 27
 16 22 28
 17 23 29

 15 21 27 16 22 28 17 23 29
dainasol

Numpy has a lot of functions, what I wanted to confirm is that this is an implementation of the arrays themselves exclusively and not, say, np.random.normal, np.any, np.tile, etc. I put "just" in double quotes because even if np.array is only one part it's still an important one and it's cool that you made this

MarcoxD

Just a heads-up, when you slice up to -1, e.g. arr[:-1], it selects the remaining of the array, but not the last element (the slicing end is exclusive). Here is an example:

In [1]: import numpy as np
In [2]: arr = np.arange(5)
In [3]: arr
Out[3]: array([0, 1, 2, 3, 4])
In [4]: arr[:]
Out[4]: array([0, 1, 2, 3, 4])
In [5]: arr[:-1]
Out[5]: array([0, 1, 2, 3])
operamint [OP]

Thanks. I don't want to support negative values (starting from end like in python) so this is intentional, but unfortunate. I will change to require use a maxint constant c_END instead.