libmir / mir-algorithm Goto Github PK

The windows function in mir.ndslice.topology can be used to create a positive band matrix. However, there is no specific iterator for band matrices. The compressed band matrix format would allow for space savings, similar to the way that the StairsIterator works.

The key use-case I can think of is interacting with lapack functions that work on symmetric positive definite band matrices.

Does Mir compatible with extern(C):? If I am trying to build hello-world project like:

import core.stdc.stdio;
extern(C):
void main()
{
}

that have mir dependence in dub.sdl I am getting linking error.

Just found this package, and I must say I love the concept of the ndslice; really raises the bar in terms of functional numerical computing, id say. Best concept since the ndarray itself!

One thing I was wondering about; I found this project while searching for numpy-D bindings; and I don't see any projects to that end yet; nor do I see, for example, matplotlib getting ported to D anytime soon. As much as I like D and the ndslice, I personally struggle to think of any projects that I work on that would get a net productivity gain from doing them in D instead of python.

So my question is; would you consider such functionality to fit within the scope of this project? It certainly would generate a lot of interest for this repo, if all the pieces were in place to make D-ndslice extensions to python code a frictionless affair.

for internal iterator

The dimensions of ndslice seem to be fixed regardless of what I am assigning.

	import mir.ndslice;

	auto matrix = slice!double(5, 2);
	matrix[] = 0;

	for(int j = 0; j < matrix.length!1; j++)
		for(int i = 0; i < matrix.length!0; i++)
			matrix[i,j] = matrix.length!0*j + i;

but, similarly,

	import mir.ndslice;

	auto matrix = slice!double(2, 5);
	matrix[] = 0;


	for(int j = 0; j < matrix.length!1; j++)
		for(int i = 0; i < matrix.length!0; i++)
			matrix[j,i] = matrix.length!0*j + i;

but both seem to identically write the matrix by skipping around rather than sequentially, which is better for larger datasets.

e.g., both write this way:

00000000 00000005 00000001 00000006 00000002 00000007 00000003 00000008 00000004 00000009

rather than sequentially, even though it should be sequential. Of course, one can't swap the dimensions and keep the same indexing as this results in an error.

Does mir do this intentionally and not worried about cache misses or is it a bug?

auto matrix = slice!double(2, 5);

is a 2 row by 5 column matrix, correct? Regardless, Should it not store the rows sequentially.

BTW, the point is not to write the correct squence to memory but to demonstrate how ndslice writes a sequential list of bytes to memory, in which case it "shuffles" them up. That is, it transposes the memory layout in to an sub-optimal form.

Consider the following simplified collision kernel taken from a lattice Boltzmann simulation:

import std.datetime;
import std.stdio;
import std.random;
import mir.ndslice;

static void collide_kernel(Tensor, Matrix)
(Tensor n, Matrix ux, Matrix uy, Matrix rho,
  immutable double[] ex, immutable double[] ey, immutable double[] w, 
  immutable int Q, immutable int NX, immutable int NY, immutable double omega) 
{
  foreach (i; 0..NX) {
    foreach (j; 0..NY) {
      double ux0 = 0.0;
      double uy0 = 0.0;
      double rho0 = 0.0;
      foreach (q; 0..Q) {
        double nq = n[i][j][q];
        rho0 += nq;
        ux0 += ex[q] * nq;
        uy0 += ey[q] * nq;
      }
      double rho_inv = 1. / rho0;
      ux0 *= rho_inv;
      uy0 *= rho_inv;
      double usqr = ux0 * ux0 + uy0 * uy0;
      foreach (q; 0..Q) {
        double eu = 3 * (ex[q] * ux0 + ey[q] * uy0);
        double neq = rho0 * w[q] * (1.0 + eu + 0.5*eu*eu - 1.5*usqr);
        n[i][j][q] = (1 - omega)*n[i][j][q] + omega*neq;
      }
      ux[i][j] = ux0;
      uy[i][j] = uy0;
      rho[i][j] = rho0;
    }
  }
}

int main(string[] args) {
  immutable uint NX = 1000;
  immutable uint NY = 1000;
  immutable uint Q  = 9;
  auto n = new double[NX*NY*Q].sliced(NX,NY,Q);
  auto ux = new double[NX*NY].sliced(NX,NY);
  auto uy = new double[NX*NY].sliced(NX,NY);
  auto rho = new double[NX*NY].sliced(NX,NY);
  immutable double[Q] ex = [0., 1., 0., -1., 0., 1., -1., -1., 1.];
  immutable double[Q] ey = [0., 0., 1., 0., -1., 1., 1., -1., -1.];
  immutable double[Q] w = [4. / 9., 
    1. / 9., 1. / 9., 1. / 9., 1. / 9., 
    1. / 36., 1. / 36., 1. / 36., 1. / 36.];
  
  Random gen;
  foreach (i; 0..NX)
    foreach (j; 0..NY) {
      ux[i][j] = 0;
      uy[i][j] = 0;
      rho[i][j] = 1;
      foreach (q; 0..Q)
        n[i][j][q] = uniform(0.0, 1.0, gen);
    }
  
  uint benchmarkCount = 100u;
  StopWatch sw;
  alias Tensor = typeof(n);
  alias Matrix = typeof(ux);
  sw.start();
  foreach (t; 0..benchmarkCount)
    collide_kernel!(Tensor, Matrix)(n, ux, uy, rho, ex, ey, w, Q, NX, NY, 0.6);
  sw.stop();
  writeln(n[0][0][0]);
  double t1 = sw.peek().usecs;
  writeln(t1 / 1000 / benchmarkCount);
  return 0;
}

when compiled with LDC 1.2.0 (via dub, using --build=release --compiler=ldc and dflags-ldc":["-mcpu=native"] in dub.json and using mir-algorithm v0.6.3) on a Haswell CPU the code yields a consistent timing of 55 ms for the collide_kernel call. The equivalent code in C++ (given in https://gist.github.com/dextorious/9a65a20e353542d6fb3a8d45c515bc18 ), compiled with clang v4.0.0 using -std=c++14 -O3 -march=native yields a consistent timing of 31 ms on the same CPU.

Inspecting the generated assembly reveals the following issues:

• Introducing the temporaries rho0, ux0, uy0 is necessary because, unlike clang, LDC doesn't seem convinced that the arrays do not alias. Notice that the C++ version does not use explicit __restrict__ annotations, so there might be room for improvement here. Without the temporaries, the D code runs in 90 ms.
• Despite the introduction of the temporaries, the inner loops are not fully unrolled and some vectorization opportunities (particularly on loads/stores) are missed. This is the major outstanding issue, because I currently don't know of a way to get to performance parity with C++.

Please let me know if I can supply any more useful information or test possible workarounds.

allow _ for 0..$

in a tensor library i've worked on I defined _ to mean 0..$
it worked as intended. Would be nice to support that in mir to make user code more readable.

implementation wise, it was something like:

enum _=R();

struct R {
  uint a0 = 0; ///from
  uint a1; ///to
  int stride = 1; ///stride
  bool is_a1_end = true; ///when is_a1_end=true, a1 should be $
  this(T1)(T1 stride) {
    this.stride = cast(int) stride;
  }

  this(T1, T2)(T1 a0, T2 a1) {
    this.a0 = cast(uint) a0;
    this.a1 = cast(uint) a1;
    is_a1_end = false;
  }

  this(T1, T2, T3)(T1 a0, T2 a1, T3 stride) {
    this(a0, a1);
    this.stride = cast(int) stride;
  }
}

and opIndex overloaded to take R arguments.

strided indexing

This allowed things like:
slice[0, _, R(2, $, 2), R(0, $, -1), ]

There is currently a kindOf to return the kind of a slice. The isSlice has a similar functionality to return packs. I suggest it be renamed packsOf and replaced with a new isSlice function that is an enum bool template that just tests if it is Slice.

const Slice!(Contiguous, [1], int*) s1 = [1, 2, 3].sliced(3);
const Slice!(Contiguous, [1], int*) s2 = [1, 2, 3].sliced(3);

writeln(s1 == s2);

... gives:

source/app.d(9,9): Error: template mir.ndslice.slice.Slice!(cast(SliceKind)2, [1LU], int*).Slice.opEquals cannot deduce function from argument types !()(const(Slice!(cast(SliceKind)2, [1LU], int*))) const, candidates are:
../../.dub/packages/mir-algorithm-0.0.29/mir-algorithm/source/mir/ndslice/slice.d(1305,10):        mir.ndslice.slice.Slice!(cast(SliceKind)2, [1LU], int*).Slice.opEquals(SliceKind rkind, ulong[] rpacks, IteratorR)(Slice!(rkind, rpacks, IteratorR) rslice) if (rpacks.sum == N)
../../.dub/packages/mir-algorithm-0.0.29/mir-algorithm/source/mir/ndslice/slice.d(1333,10):        mir.ndslice.slice.Slice!(cast(SliceKind)2, [1LU], int*).Slice.opEquals(T)(T[] arr)
dmd failed with exit code 1.

shell returned 2

In the main page of the documentation:
http://docs.algorithm.dlang.io/latest/mir_ndslice.html
The topology section does not mention indexed:
http://docs.algorithm.dlang.io/latest/mir_ndslice_topology.html#.indexed

Here is a minimalistic example to replicate the bug:

auto s = slice!float(9, 9);
auto w = s.windows(3, 3);
w[0].writeln;

Here is the compiler output:

dub run --compiler=ldc2
Performing "debug" build using ldc2 for x86_64.
mirplayground ~master: building configuration "application"...
../../../.dub/packages/mir-algorithm-0.6.5/mir-algorithm/source/mir/ndslice/slice.d(1584,58): Error: cannot implicitly convert expression (this._strides[1..3]) of type long[] to long[0]
../../../.dub/packages/mir-algorithm-0.6.5/mir-algorithm/source/mir/ndslice/slice.d(1558,20): Error: template instance mir.ndslice.slice.Slice!(cast(SliceKind)1, [2LU, 2LU], float*).Slice.opIndex!1LU error instantiating
source/app.d(103,6):        instantiated from here: opIndex!int
ldc2 failed with exit code 1.

mir.ndslice.algorithm.eachUploPair currently operates on i<=j. It would be helpful for there to be a switch at compile-time so that it operates i<j, i.e. so that it excludes the diagonal.

The following code does not compile.

struct Foo
{
	ContiguousSlice!(1, int) bar;

	int get(size_t i) const
	{
		return bar[i];
	}
}

The problem is that .opIndex cannot be called on a const object. There is an easy workaround, namely bar.toConst[i]. In other contexts, the toConst is called automatically (using some alias this construction). But here the compiler sees the non-const opIndex method and gives an error before even trying to use the alias this.

Calling toConst explicitly is easy to do but kinda unintuitive, because it is used here to remove a const attribute, not to add it. More precisely, it is used for the conversion from const Slice!int to Slice!(const int).

ndslice.iteration works with universal slices. Some functions accept dimensions in CT, so it can be determined if canonical slices can be used too.

ping @jmh530 @wilzbach

iota operates with integral and pointer types.
See IotaIterator for more details.
fiota may require slightly different API
Phobos iota`s floating point API is wrong

i had implemented slice indexing by struct field projection in my own tensor library. Would be great to support that in mir.

Here's the idea:

struct A{  int x;   string name; }
auto s=3.iota.map!(a=>A(a,"foo")).array.slice;

assert(s.fields.name==["foo", "foo"]);

// s.fields.x has a stride and backed by int*; it doesnt' allocate.
s.fields.x[]+=10;

assert(s[0]==10);

Note: this is more powerful than s.map!(a=>a.x) since the latter wouldn't be a slice backed by a int* (and in particular wouldn't contain lvalues).

implementation

N-dimensional field projection is just an indexing operation when treating A as void[0..A.sizeof]:
s.fields.x is:
s.field.Cast!(void[]).sliced(s.length, A.sizeof)[0..$, 0..A.init.x.sizeof].reinterpretCast!int;

with:

auto Cast(T, S)(S a){ return cast(T)a;}

Side question: is there an ndslice reinterpretCast ? it would be the slice analog of cast(T[]) a which just casts an array of S as an array of T without changing the bits.

NOTE

using directly opDispatch on Slice is too dangerous as it could lead to name clashes between Slice's existing fields/methods and the elements, hence we use an intermediate fields which then forwards to opDispatch. Also, not sure if fields is the best name, as it's completely different meaning than field.

another possible syntax:s.fields!"x" but maybe less nice-looking.

We need to rework std.exprimental.allocators to be betterC.

Title says it all, but here is the minimal example:

import mir.ndslice;
import mir.ndslice.sorting;

void main()
{
    float[] c;

    auto a = c.sliced;
    auto b = a.sort;
}

Terminal output of just running dub:

Performing "debug" build using dmd for x86_64.
mir-algorithm 0.6.7: target for configuration "library" is up to date.
mirplayground ~master: building configuration "application"...
Linking...
To force a rebuild of up-to-date targets, run again with --force.
Running ./mirplayground 
Program exited with code -11

(so it's a 0.6.7 mir-algorithm version)

And here's what lldb gives me on the breaking point:

Process 8356 stopped
* thread #1, queue = 'com.apple.main-thread', stop reason = EXC_BAD_ACCESS (code=1, address=0xfffffffffffffff8)
    frame #0: 0x0000000100001c20 mirplayground`_D3mir7ndslice7sorting46__T13quickSortImplS213app4mainFZ9__lambda1TPfZ13quickSortImplMFNaNbNiS3mir7ndslice5slice52__T5SliceVE3mir7ndslice5slice9SliceKindi2VAmA1i1TPfZ5SliceZv(slice=Slice @ 0x00007fff5fbff8d0) at sorting.d:230
   227 	                    //{
   228 	                        auto d = p;
   229 	                        import mir.functional: unref;
-> 230 	                        auto temp = unref(*d);
   231 	                        auto c = d;
   232 	                        ++c;
   233 	                        if (less(*c, temp))
(lldb)

@9il
what do i use for toArray?

auto s=slice!double(3);
double[]s2= s.toArray;

actually not super clear from docs

Also, from old experimental.nslice in https://dlang.org/phobos/std_experimental_ndslice.html:

T median(Range, T)(Slice!(2, Range) sl, T[] buf)
{
    import std.algorithm.sorting : topN;
    // copy sl to the buffer
    size_t n;
    foreach (r; sl)
        foreach (e; r)
            buf[n++] = e;
    buf[0 .. n].topN(n / 2);
    return buf[n / 2];
}

this snippet (copy sl to the buffer) seems to indicate no such utility to flatten elements into an array?

EDIT: OK i guess it's here:
http://docs.algorithm.dlang.io/latest/mir_ndslice_allocation.html

• somehow i find it hard to find, not sure if i'm the only one.

• IMO this should be more prominent since we want to make it easy to go from standard D to slices (plenty of examples) as well as opposite direction (very sparse in docs)

• not sure why median above doesn't use ndarray

• the docs (Creates a common n-dimensional array from a slice.) are not clear whether it allocates or merely returns a view over the data

The key part of the code is

    foreach(i, S; Slices[1 .. $])
    {
        static assert(isSlice!S == packs, "zip: all Slices must have the same shape packs");
        assert(slices[i]._lengths == slices[0]._lengths, "zip: all slices must have the same lengths");
        static if (sameStrides)
            assert(slices[i].unpack.strides == slices[0].unpack.strides, "zip: all slices must have the same strides");
    }

note that the last part of the foreach is Slices[1..$], but then slices is indexed at slices[i]. The first iteration of the loop is i=0, so it's just comparing slices[0] to itself. If Slices.length is 1, then this is all that occurs.

There are two options to fix, either change the slices[i] to slices[i+1] or to S.

In index.d, the numeric subheading has interpolation, math.sum, and math.numeric. However, the doc/Makefile has interpolation as its own separate package and those two math packages under the math heading. There is no numeric package here.

The consequence is that the dropdown on the left side does not match the ndslice overview page ordering.

I can make the fix, but I wanted to check whether the index.d should have the same structure as the doc/Makefile, or vice-versa.

and relaxation for canonical slices

    auto v = ndiota(3, 3).map!(pos => pos).slice;

... gives: Error: cannot access frame pointer of app.main.MapField!(ndIotaField!2LU, __lambda2).MapField.

This works though:

    auto v = ndiota(3, 3).slice.map!(pos => pos).slice;

mir-algorithm version is 0.0.24.

Here is the full dub output:

Performing "debug" build using dmd for x86_64.
mir-internal 0.0.2: target for configuration "library" is up to date.
mir-algorithm 0.0.24: target for configuration "library" is up to date.
mir-math 0.0.1: target for configuration "library" is up to date.
mirplayground ~master: building configuration "application"...
../../.dub/packages/mir-algorithm-0.0.24/mir-algorithm/source/mir/ndslice/field.d(66,37): Error: cannot access frame pointer of app.main.MapField!(ndIotaField!2LU, __lambda2).MapField
../../.dub/packages/mir-algorithm-0.0.24/mir-algorithm/source/mir/ndslice/iterator.d(667,31): Error: template instance app.main.mapField!(__lambda2, ndIotaField!2LU) error instantiating
../../.dub/packages/mir-algorithm-0.0.24/mir-algorithm/source/mir/ndslice/iterator.d(528,24):        instantiated from here: __map!(__lambda2)
../../.dub/packages/mir-algorithm-0.0.24/mir-algorithm/source/mir/ndslice/topology.d(2051,52):        instantiated from here: mapIterator!(__lambda2, FieldIterator!(ndIotaField!2LU))
source/app.d(9,26):        instantiated from here: map!(cast(SliceKind)2, [2LU], FieldIterator!(ndIotaField!2LU))
dmd failed with exit code 1.

shell returned 2

@9il, I tried sorting this out by myself, but I can't make heads or tails with it. That said, I think it would be really productive if you could make some sort of tutorial (wiki page maybe) about internals of ndslice, i.e. how stuff works. I'd really like not to bug you from time to time to fix a bug in ndslice, but rather to fix it myself, but I feel without this introductory tutorial from you, lot of time has to pass before I'm able to do that. I'm sure other people would also appreciate a tutorial like this.

stride should be typeof ptrdiff_t

Is it possible to "sweep" a vector across a 2-D array for instance subtracting a vector from each "column" of a 2-D array using Mir? The example in R would be:

x = matrix(1:50, nc = 5, byrow = TRUE); y = matrix(1:10, nc = 1); x; y
sweep(x, 1, y, FUN =  "-")

Hi,

I would like to be able to convert a mir slice object back into an array, I have looked in the documentation but cannot find a function that allows this - perhaps I have missed it.

I am making a little project with mir-algorithm (and other).

When building via dub with --build=release or --build=release-debug I get the following error:
(--build=debug and --build=unittest work well)

.dub/obj/MyApp.o: In function `_D3mir7ndslice5slice52__T5SliceVE3mir7ndslice5slice9SliceKindi0VAmA1i2TPdZ5Slice11__xopEqualsFKxS3mir7ndslice5slice52__T5SliceVE3mir7ndslice5slice9SliceKindi0VAmA1i2
TPdZ5SliceKxS3mir7ndslice5slice52__T5SliceVE3mir7ndslice5slice9SliceKindi0VAmA1i2TPdZ5SliceZb':                                                                                                       
.../.dub/packages/mir-algorithm-0.6.6/mir-algorithm/source/mir/ndslice/slice.d:579: undefined reference to `_D3mir7ndslice5slice52__T5SliceVE3mir7ndslice5slic
e9SliceKindi0VAmA1i2TPdZ5Slice55__T8opEqualsVE3mir7ndslice5slice9SliceKindi0VAmA1i2TPdZ8opEqualsMxFNaNbNiNexS3mir7ndslice5slice52__T5SliceVE3mir7ndslice5slice9SliceKindi0VAmA1i2TPdZ5SliceZb'        
collect2: error: ld returned 1 exit status
Error: /usr/bin/gcc failed with status: 1
ldc2 failed with exit code 1.

The strings (via core.demangle) translate to

bool mir.ndslice.slice.Slice!(0, [2], double*).Slice.__xopEquals(ref const(mir.ndslice.slice.Slice!(0, [2], double*).Slice), ref const(mir.ndslice.slice.Slice!(0, [2], double*).Slice))

and

const(pure nothrow @nogc @trusted bool function(const(mir.ndslice.slice.Slice!(0, [2], double*).Slice))) mir.ndslice.slice.Slice!(0, [2], double*).Slice.opEquals!(0, [2], double*).opEquals

package versions:

mir-algorithm 0.6.6

dub version:
DUB version 1.2.1, built on Feb 17 2017

ldc version:

LDC - the LLVM D compiler (1.3.0git-b82d87d):
based on DMD v2.073.2 and LLVM 3.9.0
built with LDC - the LLVM D compiler (0.17.5git-64a274a)
Default target: x86_64-unknown-linux-gnu
Host CPU: haswell

As said the build fails with dub build --build=release and dub build --build=release-debug
but succeeds with dub build --build=debug. Earlier versions have also succeeded in release mode.

As the last changes where rather involved and the linker error gave little hint, what in my code triggered the failure, I have no minimal failing code sample. If you could give me some hints, where to look, I would try to find it.

Univarsal should be Universal.

Also the long dash - short dash lines for shape topology and subspace topology go under the ndSlice box unlike the one between Univarsal(sic) and Dynamic.

The thought bubble should have a full stop after "Any pointer is an Iterator"

This is a bit of a bike-shed, but I think ContiguousTensor, CanonicalTensor, and UniversalTensor should be renamed ContiguousSlice, CanonicalTensor, and UniversalSlice, respectively.

My issue is that Slice is used throughout mir-algorithm and Tensor is really only used in this specific case. It would be more consistent to reuse the Slice terminology here.

Moreover, while a vector could be considered a rank-1 tensor and a matrix could be considered a rank-2 tensor, some people might be confused by the tensor terminology and think these aliases only apply to Slices with dimensions 3 or higher when the alias actually applies to all Slices.

@9il
git log -SbyElement shows eb7b255 but this doesn't seem to show the corresponding replacement

EDIT:
is that : a.byElement => a.iterator ?

http://docs.algorithm.dlang.io/latest/mir_ndslice_sorting.html

However the plain text version looks correct:

operation $(D less). Performs $(BIGOH ndslice.elementsCount) evaluations of less.
Unlike isSorted, $(LREF isStrictlyMonotonic) does not allow for equal values,
i.e. values for which both less(a, b) and less(b, a) are false.
With either function, the predicate must be a strict ordering just like with

https://github.com/libmir/mir-algorithm/edit/master/source/mir/ndslice/sorting.d

Any ideas?

Auto-deduction for tensors with smaller dimension count.

	auto matrix = slice!byte(3, 4);
	matrix[] = 0;
	matrix.diagonal[] = 1;

	auto row = matrix[2];
	row[3] = 6;

fails but changing byte to double works fine?!?! Same issue with short, etc.

for internal iterator, 1D only