ggeD.ggeD source code

1 /*
2 Copyright (c) 2022 Zenw
3 Released under the MIT license
4 https://opensource.org/licenses/mit-license.php
5 */
6 
7 module ggeD.ggeD;
8 import std;
9 import ggeD.indexVec;
10 import mir.ndslice;
11 import ggeD.einsum;
12 import mir.ndslice.topology : iota;
13 
14 
15 auto gged(T,R,Args...)(R value,Args xyz) if(isInputRange!(Unqual!R))
16 {
17     auto gg = gged!T(xyz);
18     ulong idx = 0;
19     foreach(ijk;gg.index)
20     {
21         gg[ijk] = value.front;
22         value.popFront();
23     }
24     return gg;
25 }
26 
27 auto gged(T,size_t N)(ulong[N] xyz)  
28 {
29     return GgedStruct!(T*,N,false,mir_slice_kind.contiguous)(slice!(T)(xyz));
30 }
31 auto gged(T,Args...)(Args xyz)  if(allSameType!(Args) && isIntegral!(Args[0])) 
32 {
33     return GgedStruct!(T*,Args.length,false,mir_slice_kind.contiguous)(slice!(T)(xyz));
34 }
35 auto gged(X,ulong Y,SliceKind Z)(Slice!(X,Y,Z) slice_) 
36 {
37     return GgedStruct!(X,Y,false, Z)(slice_);
38 }
39 T gged(T)(T value) if(!__traits(isSame, TemplateOf!(T), Slice))
40 {
41     return value;
42 }
43 
44 private struct IndexLoop(GGED) 
45 {
46     private alias MakeSerialIndex(X) =  SerialIndex!X;
47     GGED gg;
48     alias IndexTypes =GGED.IndexTypes;
49     alias TypeSerialIndex = staticMap!(MakeSerialIndex,IndexTypes);
50 
51     int opApply(int delegate(TypeSerialIndex) fun) {
52         mixin(genLoop!(false,0));
53         return 1;
54     }
55     alias Rank = GGED.Rank;
56     static if(Rank > 1)
57     {
58         int opApply(int delegate(IndexVec!(IndexTypes)) fun) {
59             mixin(genLoop!(true,0));
60             return 1;
61         }
62         static foreach(N ; 1 .. Rank -1)
63         {
64             int opApply(int delegate(IndexVec!(IndexTypes[0..Rank-N]),IndexTypes[Rank-N .. $]) fun) {
65                 mixin(genLoop!(true,N));
66                 return 1;
67             }
68         }
69     }
70     static string genLoop(bool vec,ulong N)(){
71         string result;
72         static foreach(idx;0..Rank){
73             result ~= "foreach(_" ~idx.to!string~ ";0 .. gg._slice.shape["~idx.to!string~"])\n";
74         }
75         result ~= "fun(";
76 
77        	static if(vec) result ~= "IndexVec!(IndexTypes[0..Rank-"~N.to!string~"])(";
78         static foreach(idx;0..Rank-N)
79 	        result ~= "SerialIndex!(IndexTypes["~idx.to!string~"])(gg.invIndex!"~idx.to!string~"(_" ~idx.to!string ~"),gg._slice.shape["~idx.to!string~"],gg.offsets["~idx.to!string~"]),";
80        	static if(vec) result ~= "),";
81 
82         static foreach(idx;Rank - N .. Rank)
83             result ~=  "SerialIndex!(IndexTypes["~idx.to!string~"])(gg.invIndex!"~idx.to!string~"(_" ~idx.to!string ~"),gg._slice.shape["~idx.to!string~"],gg.offsets["~idx.to!string~"]),";
84 
85         result ~= ");";
86         return result;
87     }
88 }
89 
90 @nogc nothrow auto setOffset(GGED,ulong Rank)(GGED gged,double[Rank] offsets) if(is(GGED == GgedStruct!(T,Rank,isOffsetIndex,kind),T,bool isOffsetIndex,SliceKind kind))
91 {
92     return GgedStruct!(GGED.TypePointer,GGED.Rank,true,GGED.Kind)(gged._slice,offsets.tupleof);
93 }
94 
95 @nogc nothrow auto setOffset(R,ulong Rank)(R value,double[Rank] offsets) if(!is(R == GgedStruct!(T,Rank,isOffsetIndex,kind),T,bool isOffsetIndex,SliceKind kind) )
96 {
97     return value;
98 }
99 
100 
101 @nogc nothrow auto setOffset(GGED,Args...)(GGED gged,Args offsets) 
102 {
103     static if(is(GGED == GgedStruct!(T,Rank,isOffsetIndex,kind),T,ulong Rank,bool isOffsetIndex,SliceKind kind) && Args.length == Rank && allSatisfy!(isFloatingPoint,Args))
104     return GgedStruct!(GGED.TypePointer,GGED.Rank,true,GGED.Kind)(gged._slice,offsets);
105     else
106     return gged;
107 }
108 
109  auto sum(T,ulong N,bool offsetflag,SliceKind kind)(GgedStruct!(T,N,offsetflag,kind) g)
110 {
111     auto result = 0.;
112     foreach(ijk;g.index)
113     {
114         result += g[ijk];
115     }
116     return result;
117 }
118 
119 alias Gged(T,ulong Rank, bool isOffsetIndex = false,SliceKind kind = SliceKind.contiguous) = GgedStruct!(T*,Rank,isOffsetIndex,kind);
120 
121 struct GgedStruct(T,ulong RANK, bool isOffsetIndex = false,SliceKind kind = SliceKind.contiguous)
122 {
123 	import mir.ndslice;
124     alias SliceType = Slice!(T, Rank, kind);
125     SliceType _slice;
126     alias Kind = kind;
127 
128     static if(isOffsetIndex)
129         alias IndexTypes = Repeat!(RANK,double);
130     else
131         alias IndexTypes = Repeat!(RANK,size_t);
132 
133     IndexTypes offsets = 0;
134 
135     private size_t getIndex(ulong n,Arg)(Arg idx)
136     {
137         return cast(size_t)(idx-offsets[n]);
138     }
139 
140     auto dup()
141     {
142         auto d = gged!(Type)(shape).setOffset(offsets);
143         foreach(ijk;this.index)
144         {
145             d[ijk] = this[ijk];
146         }
147         return d;
148     }
149 
150     private IndexTypes[n] invIndex(ulong n)(size_t idx)
151     {
152         return cast(IndexTypes[n])(idx+offsets[n]);
153     }
154     private auto getIndexes(Args...)(Args idxs) if(Args.length == RANK)
155     {
156         template Filting(X)
157         {
158             static if(isContiguousSlice!X) alias Filting = X;
159             else alias Filting = size_t;
160         }
161         struct Dummy{
162             staticMap!(Filting,Args) value;
163             alias value this;
164         }
165         Dummy result;
166         static foreach(i;0..RANK)
167         {
168             static if(isContiguousSlice!(Args[i])) 
169                 result[i] =  idxs[i];
170             else
171                 result[i] = getIndex!i(idxs[i]);
172         }
173         return result;
174     }
175     
176     
177     static if(isPointer!T)
178     {
179     	alias Type = PointerTarget!T;
180     }
181     else static if(isArray!T)
182     {
183         alias Type = ElementType!T;
184     }
185     else
186     {
187         alias Type = ReturnType!(T.opUnary!"*");
188     } 
189 	alias TypePointer = T;
190 	alias Rank = Alias!(RANK);
191 
192     alias _slice this;
193     
194     @nogc nothrow auto shape() => _slice.shape; 
195     
196     @nogc nothrow auto index(){
197         return IndexLoop!(typeof(this))(this);
198     }
199     template PickupArray(Ns...)
200         {
201             alias PickupArray = AliasSeq!();
202             static foreach(i; Ns)
203             {
204                 PickupArray = AliasSeq!(PickupArray,offsets[i]);
205             }
206         }
207 
208     /// Ns次元だけを走査するindexでかえす
209     @nogc nothrow auto index(Ns...)(){
210         template allInDim(ulong N){alias allInDim= Alias!(N<Rank);}
211         static assert(allSatisfy!(allInDim,Ns),"all Ns should be lower than Rank");
212         auto newone = this.pick!(Ns)(offsets).setOffset(PickupArray!(Ns));
213         return IndexLoop!(typeof(newone))(newone);
214     }
215     
216     auto toString()=>_slice.to!string;
217 
218     /// pick!(Args...)(ijk)
219     /// Argsにない次元の要素がijkの位置を含むArgs.length次元のgged配列を返す
220     /// pick!0(ijk)では、ijkを含む0次元目への1次配列になる
221     template pick(ns...)
222     {
223         @nogc nothrow auto pick(Args...)(Args value)
224         {
225             mixin(genPick([ns]));
226         }
227     }
228         
229     private static string genPick(ulong[] ns)
230     {
231         string result =  "return this[";
232         foreach(i ; 0..Rank)
233         {
234             if(!ns.canFind(i))
235                 result ~= "value["~i.to!string~"],";
236             else
237                 result ~= "offsets["~i.to!string~"] .. $,";
238         }
239         result ~= "];";
240         return result;
241     }
242     @nogc nothrow auto opSlice(X,Y)(X start, Y end) if(RANK == 1 && is(X == SerialIndex!(IndexTypes[0])) && is(Y == SerialIndex!(IndexTypes[0])))
243     {
244         return gged(_slice.opSlice(start,end));
245     }
246     @nogc nothrow auto opSlice(size_t dim,X,Y)(X start, Y end) 
247     {
248         return _slice.opSlice!dim(getIndex!dim(start),getIndex!dim(end));
249     }
250     @nogc nothrow auto opIndex(IndexVec!(IndexTypes) args){
251         return gged(_slice.opIndex(getIndexes(args.idx).value ));
252     }
253     static foreach(N; 1 .. Rank-1)
254     {
255         @nogc nothrow auto opIndex(Args...)(IndexVec!(IndexTypes[0..Rank-N]) args1,Args args2) if(Args.length == N) {
256             return gged(_slice.opIndex(getIndexes(args1.idx,args2).value));
257         }
258     }
259     @nogc nothrow auto opIndex(Args...)(Args args) if(Args.length == Rank) {
260         return gged(_slice.opIndex(getIndexes(args).value));
261     }
262     import std.traits;
263     static if( __traits(compiles, () {T itr; itr[0] = Type.init; } ))
264     {
265         @nogc nothrow auto opIndexAssign(AssignType)(AssignType value,IndexVec!(IndexTypes) args){
266             _slice.opIndexAssign(value,getIndexes(args.idx).value);
267         }
268         
269         static foreach(N; 1 .. Rank-1)
270         {
271             @nogc nothrow auto opIndexAssign(AssignType,Args...)(AssignType value,IndexVec!(IndexTypes[0..Rank-N])  args1,Args args2) if(Args.length == N) {
272                 _slice.opIndexAssign(value,getIndexes(args1.idx,args2).value);
273             }
274         }
275         @nogc nothrow auto opIndexAssign(AssignType,Args...)(AssignType value,Args args) if(Args.length == Rank)  {
276             _slice.opIndexAssign(value,getIndexes(args).value);
277         }
278     }
279     @nogc nothrow auto opDollar(ulong dim)(){
280         return invIndex!dim(_slice.opDollar!dim);
281     }
282     @nogc nothrow auto opDispatch(string idx)() if(idx.replace("_","").length == Rank)
283     {
284         return Leaf!("",idx.replace("_",""),typeof(this))(this);
285     }
286     auto opEquals(RHS)(RHS rhs)
287     {
288         static if(is(RHS == GgedStruct!(T2,RANK,isOffsetIndex2,kind2),T2,bool isOffsetIndex2,SliceKind kind2))
289         {
290             bool result = true;
291             foreach(ijk ; index)
292             {
293                 result &= isClose(rhs[ijk] , this[ijk] );
294             }
295             return result;
296         }
297         else
298         {
299             return _slice == rhs;
300         }
301     }
302 
303     auto opBinary(string op,GGED)(GGED rhs) 
304     {
305         static if(is(GGED==GgedStruct!(T2,Rank,isOffsetIndex2,kind2),T2,bool isOffsetIndex2,SliceKind kind2))
306             return _slice.opBinary!(op,T2,Rank,kind2)(rhs._slice).gged.setOffset(offsets);
307         else
308             return _slice.opBinary!(op,GGED)(rhs).gged.setOffset(offsets);
309     }
310     
311     auto opUnary(string op)()
312         if (op == "*" || op == "~" || op == "-" || op == "+")
313     {
314         import mir.ndslice.topology: map;
315         static if (op == "+")
316             return this;
317         else
318             return _slice.opUnary!op.gged.setOffset(offsets);
319     }
320 }