source: project/release/4/9ML-toolkit/trunk/examples/tensor.sml @ 28666

Last change on this file since 28666 was 28666, checked in by Ivan Raikov, 7 years ago

9ML-toolkit: bug fixes to indexing in tensor SML library

File size: 109.4 KB
Line 
1(* Obtained at http://www.arrakis.es/~worm/ *)
2
3signature MONO_VECTOR =
4  sig
5    type vector
6    type elem
7    val maxLen : int
8    val fromList : elem list -> vector
9    val tabulate : (int * (int -> elem)) -> vector
10    val length : vector -> int
11    val sub : (vector * int) -> elem
12    val extract : (vector * int * int option) -> vector
13    val concat : vector list -> vector
14    val mapi : ((int * elem) -> elem) -> (vector * int * int option) -> vector
15    val map : (elem -> elem) -> vector -> vector
16    val appi : ((int * elem) -> unit) -> (vector * int * int option) -> unit
17    val app : (elem -> unit) -> vector -> unit
18    val foldli : ((int * elem * 'a) -> 'a) -> 'a -> (vector * int * int option) -> 'a
19    val foldri : ((int * elem * 'a) -> 'a) -> 'a -> (vector * int * int option) -> 'a
20    val foldl : ((elem * 'a) -> 'a) -> 'a -> vector -> 'a
21    val foldr : ((elem * 'a) -> 'a) -> 'a -> vector -> 'a 
22  end
23
24(*
25 Copyright (c) Juan Jose Garcia Ripoll.
26 All rights reserved.
27
28 Refer to the COPYRIGHT file for license conditions
29*)
30
31(* COPYRIGHT
32 
33Redistribution and use in source and binary forms, with or
34without modification, are permitted provided that the following
35conditions are met:
36
371. Redistributions of source code must retain the above copyright
38   notice, this list of conditions and the following disclaimer.
39
402. Redistributions in binary form must reproduce the above
41   copyright notice, this list of conditions and the following
42   disclaimer in the documentation and/or other materials provided
43   with the distribution.
44
453. All advertising materials mentioning features or use of this
46   software must display the following acknowledgement:
47        This product includes software developed by Juan Jose
48        Garcia Ripoll.
49
504. The name of Juan Jose Garcia Ripoll may not be used to endorse
51   or promote products derived from this software without
52   specific prior written permission.
53
54THIS SOFTWARE IS PROVIDED BY JUAN JOSE GARCIA RIPOLL ``AS IS''
55AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
56TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
57PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL HE BE
58LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
59OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
60PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
61OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
62THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
63TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
64OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
65OF SUCH DAMAGE.
66*)
67
68structure Loop =
69    struct
70        fun all (a, b, f) =
71            if a > b then
72                true
73            else if f a then
74                all (a+1, b, f)
75            else
76                false
77
78        fun any (a, b, f) =
79            if a > b then
80                false
81            else if f a then
82                true
83            else
84                any (a+1, b, f)
85
86        fun app (a, b, f) =
87            if a < b then
88                (f a; app (a+1, b, f))
89            else
90                ()
91
92        fun app' (a, b, d, f) =
93            if a < b then
94                (f a; app' (a+d, b, d, f))
95            else
96                ()
97
98        fun appi' (a, b, d, f) =
99            if a < b then
100                (f a; appi' (a+d, b, d, f))
101            else
102                ()
103    end
104(*
105  INDEX         -Signature-
106
107  Indices are a enumerable finite set of data with an order and a map
108  to a continous nonnegative interval of integers.  In the sample
109  implementation, Index, each index is a list of integers,
110        [i1,...,in]
111  and each set of indices is defined by a shape, which has the same
112  shape of an index but with each integer incremented by one
113        shape = [k1,...,kn]
114        0 <= i1 < k1
115
116  type storage = RowMajor | ColumnMajor
117  order : storage
118        Identifies:
119                1) the underlying algorithms for this structure
120                2) the most significant index
121                3) the index that varies more slowly
122                4) the total order
123        RowMajor means that first index is most significant and varies
124        more slowly, while ColumnMajor means that last index is the most
125        significant and varies more slowly. For instance
126                RowMajor => [0,0]<[0,1]<[1,0]<[1,1] (C, C++, Pascal)
127                ColumnMajor => [0,0]>[1,0]>[0,1]>[1,1] (Fortran)
128  last shape
129  first shape
130        Returns the last/first index that belongs to the sed defined by
131        'shape'.
132  inBounds shape index
133        Checkes whether 'index' belongs to the set defined by 'shape'.
134  toInt shape index
135        As we said, indices can be sorted and mapped to a finite set of
136        integers. 'toInt' obtaines the integer number that corresponds to
137        a certain index.
138  indexer shape
139        It is equivalent to the partial evaluation 'toInt shape' but
140        optimized for 'shape'.
141
142  next shape index
143  prev shape index
144  next' shape index
145  prev' shape index
146        Obtain the following or previous index to the one we supply.
147        next and prev return an object of type 'index option' so that
148        if there is no such following/previous, the output is NONE.
149        On the other hand, next'/prev' raise an exception when the
150        output is not well defined and their output is always of type
151        index. next/prev/next'/prev' raise an exception if 'index'
152        does not belong to the set of 'shape'.
153
154  all shape f
155  any shape f
156  app shape f
157        Iterates 'f' over every index of the set defined by 'shape'.
158        'all' stops when 'f' first returns false, 'any' stops when
159        'f' first returns true and 'app' does not stop and discards the
160        output of 'f'.
161
162  compare(a,b)
163        Returns LESS/GREATER/EQUAL according to the total order which
164        is defined in the set of all indices.
165  <,>,eq,<=,>=,<>
166        Reduced comparisons which are defined in terms of 'compare'.
167
168  validShape t
169  validIndex t
170        Checks whether 't' conforms a valid shape or index.
171
172  iteri shape f
173*)
174
175signature INDEX =
176    sig
177        type t
178        type indexer = t -> int
179        datatype storage = RowMajor | ColumnMajor
180
181        exception Index
182        exception Shape
183
184        val order : storage
185        val toInt : t -> t -> int
186        val length : t -> int
187        val first : t -> t
188        val last : t -> t
189        val next : t -> t -> t option
190        val prev : t -> t -> t option
191        val next' : t -> t -> t
192        val prev' : t -> t -> t
193        val indexer : t -> (t -> int)
194
195        val inBounds : t -> t -> bool
196        val compare : t * t -> order
197        val < : t * t -> bool
198        val > : t * t -> bool
199        val eq : t * t -> bool
200        val <= : t * t -> bool
201        val >= : t * t -> bool
202        val <> : t * t -> bool
203        val - : t * t -> t
204
205        val validShape : t -> bool
206        val validIndex : t -> bool
207
208        val all : t -> (t -> bool) -> bool
209        val any : t -> (t -> bool) -> bool
210        val app : t -> (t -> unit) -> unit
211    end
212
213structure Index : INDEX =
214    struct
215        type t = int list
216        type indexer = t -> int
217        datatype storage = RowMajor | ColumnMajor
218
219        exception Index
220        exception Shape
221
222        val order = ColumnMajor
223
224        fun validShape shape = List.all (fn x => x > 0) shape
225
226        fun validIndex index = List.all (fn x => x >= 0) index
227
228        fun toInt shape index =
229            let fun loop ([], [], accum, p) = accum
230                  | loop ([], _, _, _) = raise Index
231                  | loop (_, [], _, _) = raise Index
232                  | loop (i::ri, l::rl, accum, p) =
233                if (i >= 0) andalso (i < l) then
234                    loop (ri, rl, i * p + accum, p * l)
235                else
236                    raise Index
237            in loop (index, shape, 0, 1)
238            end
239
240        (* ----- CACHED LINEAR INDEXER -----
241
242           An indexer is a function that takes a list of
243           indices, validates it and produces a nonnegative
244           integer number. In short, the indexer is the
245           mapper from indices to element positions in
246           arrays.
247
248           'indexer' builds such a mapper by optimizing
249           the most common cases, which are 1d and 2d
250           tensors.
251         *)
252    local
253        fun doindexer [] _ = raise Shape
254          | doindexer [a] [dx] =
255            let fun f [x] = if (x > 0) andalso (x < a)
256                            then x
257                            else raise Index
258                  | f _ = raise Index
259            in f end
260          | doindexer [a,b] [dx, dy] =
261                let fun f [x,y] = if ((x >= 0) andalso (x < a) andalso
262                                      (y >= 0) andalso (y < b))
263                                  then x + dy * y
264                                  else raise Index
265                      | f _ = raise Index
266                in f end
267          | doindexer [a,b,c] [dx,dy,dz] =
268                let fun f [x,y,z] = if ((x >= 0) andalso (x < a) andalso
269                                        (y >= 0) andalso (y < b) andalso
270                                        (z >= 0) andalso (z < c))
271                                    then x + dy * y + dz * z
272                                    else raise Index
273                      | f _ = raise Index
274                in f end
275          | doindexer shape memo =
276                let fun f [] [] accum [] = accum
277                      | f _  _  _ [] = raise Index
278                      | f (fact::rf) (ndx::ri) accum (dim::rd) =
279                        if (ndx >= 0) andalso (ndx < dim) then
280                            f rf ri (accum + ndx * fact) rd
281                        else
282                            raise Index
283                in f shape memo 0
284                end
285    in
286        fun indexer shape =
287            let fun memoize accum [] = []
288                  | memoize accum (dim::rd) =
289                accum :: (memoize (dim * accum) rd)
290            in if validShape shape then
291                   doindexer shape (memoize 1 shape)
292               else
293                   raise Shape
294            end
295    end
296
297        fun length shape =
298            let fun prod (a,b) =
299                if b < 0 then raise Shape else a * b
300            in foldl prod 1 shape
301            end
302
303        fun first shape = map (fn x => 0) shape
304
305        fun last [] = []
306          | last (size :: rest) = size - 1 :: last rest
307
308        fun next' [] [] = raise Subscript
309          | next' _ [] = raise Index
310          | next' [] _ = raise Index
311          | next' (dimension::restd) (index::resti) =
312            if (index + 1) < dimension then
313                (index + 1) :: resti
314            else
315                0 :: (next' restd resti)
316
317        fun prev' [] [] = raise Subscript
318          | prev' _ [] = raise Index
319          | prev' [] _ = raise Index
320          | prev' (dimension::restd) (index::resti) =
321            if (index > 0) then
322                (index - 1) :: resti
323            else
324                (dimension - 1) :: (prev' restd resti)
325
326        fun next shape index = (SOME (next' shape index)) handle
327            Subscript => NONE
328
329        fun prev shape index = (SOME (prev' shape index)) handle
330            Subscript => NONE
331
332        fun inBounds shape index =
333            ListPair.all (fn (x,y) => (x >= 0) andalso (x < y))
334            (index, shape)
335
336        fun compare ([],[]) = EQUAL
337          | compare (_, []) = raise Index
338          | compare ([],_) = raise Index
339          | compare (a::ra, b::rb) =
340            case Int.compare (a,b) of
341                EQUAL => compare (ra,rb)
342              | LESS => LESS
343              | GREATER => GREATER
344
345    local
346        fun iterator a inner =
347            let fun loop accum f =
348                let fun innerloop i =
349                    if i < a then
350                        if inner (i::accum) f then
351                            innerloop (i+1)
352                        else
353                            false
354                    else
355                        true
356                in innerloop 0
357                end
358            in loop
359            end
360        fun build_iterator [a] =
361            let fun loop accum f =
362                let fun innerloop i =
363                    if i < a then
364                        if f (i::accum) then
365                            innerloop (i+1)
366                        else
367                            false
368                    else
369                        true
370                in innerloop 0
371                end
372            in loop
373            end
374          | build_iterator (a::rest) = iterator a (build_iterator rest)
375    in
376        fun all shape = build_iterator shape []
377    end
378
379    local
380        fun iterator a inner =
381            let fun loop accum f =
382                let fun innerloop i =
383                    if i < a then
384                        if inner (i::accum) f then
385                            true
386                        else
387                            innerloop (i+1)
388                    else
389                        false
390                in innerloop 0
391                end
392            in loop
393            end
394        fun build_iterator [a] =
395            let fun loop accum f =
396                let fun innerloop i =
397                    if i < a then
398                        if f (i::accum) then
399                            true
400                        else
401                            innerloop (i+1)
402                    else
403                        false
404                in innerloop 0
405                end
406            in loop
407            end
408          | build_iterator (a::rest) = iterator a (build_iterator rest)
409    in
410        fun any shape = build_iterator shape []
411    end
412
413    local
414        fun iterator a inner =
415            let fun loop accum f =
416                let fun innerloop i =
417                    case i < a of
418                        true => (inner (i::accum) f; innerloop (i+1))
419                      | false => ()
420                in innerloop 0
421                end
422            in loop
423            end
424        fun build_iterator [a] =
425            let fun loop accum f =
426                let fun innerloop i =
427                    case i < a of
428                        true => (f (i::accum); innerloop (i+1))
429                      | false => ()
430                in innerloop 0
431                end
432            in loop
433            end
434          | build_iterator (a::rest) = iterator a (build_iterator rest)
435    in
436        fun app shape = build_iterator shape []
437    end
438
439        fun a < b = compare(a,b) = LESS
440        fun a > b = compare(a,b) = GREATER
441        fun eq (a, b) = compare(a,b) = EQUAL
442        fun a <> b = not (a = b)
443        fun a <= b = not (a > b)
444        fun a >= b = not (a < b)
445        fun a - b = ListPair.map Int.- (a,b)
446
447    end
448
449
450signature RANGE =
451    sig
452        structure Index : INDEX
453        type index = Index.t
454        type t
455
456        exception Range
457
458        val fromto : index -> index * index -> t
459        val fromto' : index -> index * index -> t
460        val upto : index -> index -> t
461        val below : index -> index -> t
462
463        val first : t -> index
464        val last : t -> index
465
466        val length : t -> int
467        val shape : t -> index
468        val inRange : t -> index -> bool
469        val next : t -> index -> index option
470        val prev : t -> index -> index option
471        val ranges : index -> ((index * index) list) -> t
472
473        val iteri : (index -> bool) -> t -> bool
474        val iteri2 : (index * index -> bool) -> (t * t) -> bool
475
476       
477    end
478
479structure Range : RANGE =
480struct
481
482    structure Index = Index
483    type index = Index.t
484
485    datatype t = RangeEmpty | RangeIn of index * index * index | RangeSet of index * ((index * index) list)
486
487    exception Range
488
489    local
490        fun next'' [] [] [] = raise Range
491          | next'' [] _  _  = raise Index.Index
492          | next'' _  [] _  = raise Index.Index
493          | next'' _  _  [] = raise Index.Index
494          | next'' (low::rl) (up::ru) (index::ri) =
495            if index < up then
496                index + 1 :: ri
497            else
498                low :: (next'' rl ru ri)
499
500        fun prev'' [] [] [] = raise Range
501          | prev'' [] _  _  = raise Index.Index
502          | prev'' _  [] _  = raise Index.Index
503          | prev'' _  _  [] = raise Index.Index
504          | prev'' (low::rl) (up::ru) (index::ri) =
505            if index > low then
506                index - 1 :: ri
507            else
508                up :: (prev'' rl ru ri)
509
510        (* Builds the simple loop
511           for i := first to last
512              if not (g (i::ndx)) then
513                  break
514           endfor;
515           *)
516
517        fun simple_loop (first : int) (last : int) =
518            let fun loop (ndx : index) (g: index -> bool) =
519                let fun innerloop i =
520                    if i > last then
521                        true
522                    else if g (i::ndx) then
523                        innerloop (i+1)
524                    else
525                        false
526                in innerloop first end
527            in loop end
528
529        (* Builds the nested loop
530           for i := first to last
531              if not (f (i:ndx) g) then
532                  break
533           endfor
534         *)
535        fun nested_loop f (first : int) (last : int) =
536            let fun loop (ndx: index) (g: index -> bool) =
537                let fun innerloop i =
538                     (if i > last then
539                        true
540                    else if (f (i::ndx) g) then
541                        innerloop (i+1)
542                    else
543                        false)
544                in 
545                     innerloop first
546                end
547            in loop end
548
549        fun build_iterator ([a] : index) ([b] : index) = 
550            simple_loop a b
551          | build_iterator (a::ra) (b::rb) =
552            nested_loop (build_iterator ra rb) a b
553
554        fun simple_loop2 (first : int) (last : int) (first' : int) (last' : int) =
555            let fun loop (ndx : index) (ndx' : index) (g: index * index -> bool) =
556                let fun innerloop (i,j) =
557                    if i > last andalso j > last' then
558                        true
559                    else if g (i::ndx,j::ndx') then
560                        innerloop (i+1,j+1)
561                    else
562                        false
563                in innerloop (first,first') end
564            in loop end
565
566        fun nested_loop2 f (first : int) (last : int) (first' : int) (last' : int) =
567            let fun loop (ndx: index) (ndx': index) (g: index * index -> bool) =
568                let fun innerloop (i,j) =
569                    (if i > last andalso j > last' then
570                        true
571                    else if (f (i::ndx) (j::ndx') g) then
572                        innerloop (i+1,j+1)
573                    else
574                        false)
575                in 
576                     innerloop (first,first')
577                end
578            in loop end
579
580        fun build_iterator2 ([a] : index) ([b] : index) ([a'] : index) ([b'] : index) = 
581            simple_loop2 a b a' b'
582          | build_iterator2 (a::ra) (b::rb) (a'::ra') (b'::rb') =
583            nested_loop2 (build_iterator2 ra rb ra' rb') a b a' b'
584
585    in
586
587        (* ----- CONSTRUCTORS ----- *)
588
589        fun fromto shape (lo, up) =
590            if (Index.validShape shape) andalso (Index.validIndex lo) andalso (Index.validIndex up) andalso
591               (Index.inBounds shape lo) andalso (Index.inBounds shape up) andalso Index.< (lo,up)
592            then
593                RangeIn(shape,lo,up)
594            else
595                RangeEmpty
596
597        fun fromto' shape (lo, up) = fromto shape (lo, (Index.prev' shape up))
598
599        fun upto shape index = fromto shape (Index.first index, index)         
600
601        fun below shape index = fromto' shape (Index.first index, index)
602
603        fun range_append ((lo,up),((lo',up')::ranges)) = 
604            if Index.< (up,lo') then ((lo,up)::(lo',up')::ranges) else
605               (if Index.> (up,up') then (lo',up')::(range_append ((lo,up),ranges)) else
606                 (if Index.> (lo,lo') then ((lo',up')::ranges) else ((lo,up')::ranges) ))
607          | range_append ((lo,up),[]) = [(lo,up)]
608
609
610        fun listWrite converter file x =
611        (List.app (fn x => (TextIO.output(file, "," ^ (converter x)))) x)
612
613        fun intListWrite file x = listWrite Int.toString file x
614
615        fun ranges' shape ((lo, up) :: rest) = 
616            (if (Index.validShape shape) andalso (Index.validIndex lo) andalso (Index.validIndex up) andalso
617               (Index.inBounds shape lo) andalso (Index.inBounds shape up) andalso Index.< (lo,up)
618               then range_append ((lo,up), (ranges' shape rest)) else (ranges' shape rest))
619            | ranges' shape ([]) = []
620
621        fun ranges shape xs = 
622            let val set = ranges' shape xs
623            in 
624                if List.null set then RangeEmpty else RangeSet (shape,set)
625            end
626
627        fun length RangeEmpty = 0
628          | length (RangeIn(shape,lo,up)) =
629            let fun diff (x,y) = (y-x+1) in
630                Index.length (ListPair.map diff (lo,up))
631            end
632          | length (RangeSet(shape,set)) =
633            let fun diff (x,y) = (y-x+1) in
634                foldl (fn ((lo,up),ax) => (Index.length (ListPair.map diff (lo,up))) + ax)
635                      0 set
636            end
637 
638        fun shape RangeEmpty = []
639          | shape (RangeIn(shape,lo,up)) =
640            let fun diff (x,y) = (y-x+1) in
641                ListPair.map diff (lo,up)
642            end
643          | shape (RangeSet(shape,set)) =
644            let fun diff (x,y) = (y-x+1) in
645                foldl (fn ((lo,up),ax) => ListPair.map (op +) ((ListPair.map diff (lo,up)), ax))
646                      (List.map (fn (x) => 0) shape) set
647            end
648 
649
650        fun first RangeEmpty = raise Range
651          | first (RangeIn(shape,lo,up)) = lo
652          | first (RangeSet(shape,(lo,up)::_)) = lo
653
654        fun last RangeEmpty = raise Range
655          | last (RangeIn(shape,lo,up)) = up
656          | last (RangeSet(shape,set)) = let val (lo,up) = List.last set in up end
657
658        (* ----- PREDICATES & OPERATIONS ----- *)
659
660        fun inRange RangeEmpty _ = false
661          | inRange (RangeIn(shape,lo,up)) ndx =
662            (Index.<=(lo,ndx)) andalso (Index.<=(ndx,up))
663          | inRange (RangeSet(shape,set)) ndx =
664            List.exists (fn ((lo,up)) => (Index.<=(lo,ndx)) andalso (Index.<=(ndx,up))) set
665
666        fun next RangeEmpty _ = NONE
667          | next (RangeIn(shape,lo,up)) index =
668            (SOME (next'' lo up index) handle Range => NONE)
669          | next (RangeSet(shape,set)) index =
670            let val m = List.find (fn ((lo,up)) => (Index.<=(lo,index)) andalso (Index.<=(index,up))) set
671            in
672               case m of NONE => NONE
673                       | SOME (lo,up) => (SOME (next'' lo up index) handle Range => NONE)
674            end
675
676        fun prev RangeEmpty _ = NONE
677          | prev (RangeIn(shape,lo,up)) index =
678            (SOME (prev'' lo up index) handle Range => NONE)
679          | prev (RangeSet(shape,set)) index =
680            let val m = List.find (fn ((lo,up)) => (Index.<=(lo,index)) andalso (Index.<=(index,up))) set
681            in
682               case m of NONE => NONE
683                       | SOME (lo,up) => (SOME (prev'' lo up index) handle Range => NONE)
684            end
685
686        fun next' RangeEmpty _ = raise Range
687          | next' (RangeIn(shape,lo,up)) index = next'' lo up index
688          | next' (RangeSet(shape,set)) index =
689            let val m = List.find (fn ((lo,up)) => (Index.<=(lo,index)) andalso (Index.<=(index,up))) set
690                val (lo,up) = valOf m
691            in
692                next'' lo up index
693            end
694
695        fun prev' RangeEmpty _ = raise Range
696          | prev' (RangeIn(shape,lo,up)) index = prev'' lo up index
697          | prev' (RangeSet(shape,set)) index =
698            let val m = List.find (fn ((lo,up)) => (Index.<=(lo,index)) andalso (Index.<=(index,up))) set
699                val (lo,up) = valOf m
700            in
701                prev'' lo up index
702            end
703
704        (* ----- ITERATION ----- *)
705
706        (* Builds an interator that applies 'f' sequentially to
707           all the indices in the range, *)
708        fun iteri f RangeEmpty = f []
709          | iteri (f: index -> bool) (RangeIn(shape,lo: index,up: index)) = 
710           (build_iterator lo up) [] f
711          | iteri (f: index -> bool) (RangeSet(shape,set)) = 
712           List.all (fn (lo,up) => (build_iterator lo up) [] f) set
713
714        (* Builds an interator that applies 'f' sequentially to
715           all the indices of the two ranges, *)
716        fun iteri2 f (RangeEmpty,RangeEmpty) = f ([],[])
717          | iteri2 (f: index * index -> bool) (RangeIn(shape,lo: index,up: index),RangeIn(shape',lo': index,up': index)) = 
718            if shape=shape' then (build_iterator2 lo up lo' up') [] [] f else raise Range
719          | iteri2 (f: index * index -> bool) (RangeSet(shape,set),RangeSet(shape',set')) = 
720           if shape=shape' then ListPair.all (fn ((lo,up),(lo',up')) => (build_iterator2 lo up lo' up') [] [] f) (set,set') else raise Range
721
722    end
723end
724
725
726(*
727 Copyright (c) Juan Jose Garcia Ripoll.
728 All rights reserved.
729 
730 Refer to the COPYRIGHT file for license conditions
731*)
732
733(*
734 TENSOR         - Signature -
735
736 Polymorphic tensors of any type. With 'tensor' we denote a (mutable)
737 array of any rank, with as many indices as one wishes, and that may
738 be traversed (map, fold, etc) according to any of those indices.
739
740 type 'a tensor
741        Polymorphic tensor whose elements are all of type 'a.
742 val storage = RowMajor | ColumnMajor
743        RowMajor = data is stored in consecutive cells, first index
744        varying fastest (FORTRAN convention)
745        ColumnMajor = data is stored in consecutive cells, last
746        index varying fastest (C,C++,Pascal,CommonLisp convention)
747 new ([i1,...,in],init)
748        Build a new tensor with n indices, each of sizes i1...in,
749        filled with 'init'.
750 fromArray (shape,data)
751 fromList (shape,data)
752        Use 'data' to fill a tensor of that shape. An exception is
753        raised if 'data' is too large or too small to properly
754        fill the vector. Later use of a 'data' array is disregarded
755        -- one must think that the tensor now owns the array.
756 length tensor
757 rank tensor
758 shape tensor
759        Return the number of elements, the number of indices and
760        the shape (size of each index) of the tensor.
761 toArray tensor
762        Return the data of the tensor in the form of an array.
763        Mutation of this array may lead to unexpected behavior.
764
765 sub (tensor,[i1,...,in])
766 update (tensor,[i1,...,in],new_value)
767        Access the element that is indexed by the numbers [i1,..,in]
768
769 app f a
770 appi f a
771        The same as 'map' and 'mapi' but the function 'f' outputs
772        nothing and no new array is produced, i.e. one only seeks
773        the side effect that 'f' may produce.
774 map2 operation a b
775        Apply function 'f' to pairs of elements of 'a' and 'b'
776        and build a new tensor with the output. Both operands
777        must have the same shape or an exception is raised.
778        The procedure is sequential, as specified by 'storage'.
779 foldl operation a n
780        Fold-left the elements of tensor 'a' along the n-th
781        index.
782 all test a
783 any test a
784        Folded boolean tests on the elements of the tensor.
785*)
786
787signature TENSOR =
788    sig
789        structure Array : ARRAY
790        structure Index : INDEX
791        type index = Index.t
792        type 'a tensor
793
794        val new : index * 'a -> 'a tensor
795        val tabulate : index * (index -> 'a) -> 'a tensor
796        val length : 'a tensor -> int
797        val rank : 'a tensor -> int
798        val shape : 'a tensor -> (index)
799        val reshape : index -> 'a tensor -> 'a tensor
800        val fromList : index * 'a list -> 'a tensor
801        val fromArray : index * 'a array -> 'a tensor
802        val toArray : 'a tensor -> 'a array
803
804        val sub : 'a tensor * index -> 'a
805        val update : 'a tensor * index * 'a -> unit
806        val map : ('a -> 'b) -> 'a tensor -> 'b tensor
807        val map2 : ('a * 'b -> 'c) -> 'a tensor -> 'b tensor -> 'c tensor
808        val app : ('a -> unit) -> 'a tensor -> unit
809        val appi : (int * 'a -> unit) -> 'a tensor -> unit
810        val foldl : ('c * 'a -> 'c) -> 'c -> 'a tensor -> int -> 'c tensor
811        val all : ('a -> bool) -> 'a tensor -> bool
812        val any : ('a -> bool) -> 'a tensor -> bool
813
814        val cat : int * 'a tensor * 'a tensor -> 'a tensor       
815    end
816
817
818signature TENSOR_SLICE =
819    sig
820        structure Tensor : TENSOR
821        structure Range : RANGE
822
823        type index = Tensor.Index.t
824        type range = Range.t
825        type 'a tensor = 'a Tensor.tensor
826        type 'a slice
827
828        val slice  : ((index * index) list) * 'a tensor -> 'a slice
829        val length : 'a slice -> int
830        val base   : 'a slice -> 'a tensor
831        val shape  : 'a slice -> (index)
832
833        val app : ('a -> unit) -> 'a slice -> unit
834        val map : ('a -> 'b) -> 'a slice -> 'b Vector.vector
835        val map2 : ('a * 'b -> 'c) -> 'a slice -> 'b slice -> 'c Vector.vector
836        val foldl  : ('a * 'b -> 'b) -> 'b -> 'a slice -> 'b
837
838    end
839
840(*
841 Copyright (c) Juan Jose Garcia Ripoll.
842 All rights reserved.
843 
844 Refer to the COPYRIGHT file for license conditions
845*)
846
847structure Tensor : TENSOR =
848    struct
849        structure Array = Array
850        structure Index = Index
851           
852        type index = Index.t
853        type 'a tensor = {shape : index, indexer : Index.indexer, data : 'a array}
854
855        exception Shape
856        exception Match
857        exception Index
858
859    local
860    (*----- LOCALS -----*)
861
862        fun make' (shape, data) =
863            {shape = shape, indexer = Index.indexer shape, data = data}
864
865        fun toInt {shape, indexer, data} index = indexer index
866
867        fun array_map f a =
868            let fun apply index = f(Array.sub(a,index)) in
869                Array.tabulate(Array.length a, apply)
870            end
871
872        fun splitList (l as (a::rest), place) =
873            let fun loop (left,here,right) 0 =  (List.rev left,here,right)
874                  | loop (_,_,[]) place = raise Index
875                  | loop (left,here,a::right) place = 
876                loop (here::left,a,right) (place-1)
877            in
878                if place <= 0 then
879                    loop ([],a,rest) (List.length rest - place)
880                else
881                    loop ([],a,rest) (place - 1)
882            end
883
884    in
885    (*----- STRUCTURAL OPERATIONS & QUERIES ------*)
886
887      fun cat (dim, x: 'a tensor, y: 'a tensor) =
888        (let val xshape = (#shape x)
889             val yshape = (#shape y)
890             val xdata  = (#data x)
891             val ydata  = (#data y)
892        in
893           if  not (length xshape  = length yshape) then
894           raise Shape
895           else
896                let 
897                   val (_,newshape)   = ListPair.foldl
898                                      (fn (x,y,(i,ax)) => if (dim = i) then (i+1,(x+y) :: ax) 
899                                                                       else if not (x=y) then raise Shape else (i+1,x :: ax))
900                                       (0,[]) (xshape, yshape)
901                   val newlength  = Index.length newshape
902                   val newdata    = Array.array(newlength,Array.sub(xdata,0))
903                in
904                    Array.copy {src=xdata,dst=newdata,di=0};
905                    Array.copy {src=ydata,dst=newdata,di=(Index.length xshape)};
906                    {shape = newshape,
907                     indexer = Index.indexer newshape,
908                     data = newdata}
909                end
910        end)
911
912        fun new (shape, init) =
913            if not (Index.validShape shape) then
914                raise Shape
915            else
916                let val length = Index.length shape in
917                    {shape = shape,
918                     indexer = Index.indexer shape,
919                     data = Array.array(length,init)}
920                end
921
922        fun toArray {shape, indexer, data} = data
923
924        fun length {shape, indexer, data} =  Array.length data
925
926        fun shape {shape, indexer, data} = shape
927
928        fun rank t = List.length (shape t)
929
930        fun reshape new_shape tensor =
931            if Index.validShape new_shape then
932                case (Index.length new_shape) = length tensor of
933                    true => make'(new_shape, toArray tensor)
934                  | false => raise Match
935            else
936                raise Shape
937
938        fun fromArray (s, a) =
939            case Index.validShape s andalso 
940                 ((Index.length s) = (Array.length a)) of
941                 true => make'(s, a)
942               | false => raise Shape
943
944        fun fromList (s, a) = fromArray (s, Array.fromList a)
945
946        fun tabulate (shape,f) =
947            if Index.validShape shape then
948                let val last = Index.last shape
949                    val length = Index.length shape
950                    val c = Array.array(length, f last)
951                    fun dotable (c, indices, i) =
952                        (Array.update(c, i, f indices);
953                         case i of
954                             0 => c
955                           | i => dotable(c, Index.prev' shape indices, i-1))
956                in
957                    make'(shape,dotable(c, Index.prev' shape last, length-1))
958                end
959            else
960                raise Shape
961
962        (*----- ELEMENTWISE OPERATIONS -----*)
963
964        fun sub (t, index) = Array.sub(#data t, toInt t index)
965
966        fun update (t, index, value) =
967            Array.update(toArray t, toInt t index, value)
968
969        fun map f {shape, indexer, data} =
970            {shape = shape, indexer = indexer, data = array_map f data}
971
972        fun map2 f t1 t2=
973            let val {shape, indexer, data} = t1
974                val {shape=shape2, indexer=indexer2, data=data2} = t2
975                fun apply i = f (Array.sub(data,i), Array.sub(data2,i))
976                val len = Array.length data
977            in
978                if Index.eq(shape, shape2) then
979                    {shape = shape,
980                     indexer = indexer,
981                     data = Array.tabulate(len, apply)}
982                else
983                    raise Match
984        end
985
986        fun appi f tensor = Array.appi f (toArray tensor)
987
988        fun app f tensor = Array.app f (toArray tensor)
989
990        fun all f tensor =
991            let val a = toArray tensor
992            in Loop.all(0, length tensor - 1, fn i =>
993                        f (Array.sub(a, i)))
994            end
995
996        fun any f tensor =
997            let val a = toArray tensor
998            in Loop.any(0, length tensor - 1, fn i =>
999                        f (Array.sub(a, i)))
1000            end
1001
1002        fun foldl f init {shape, indexer, data=a} index =
1003            let val (head,lk,tail) = splitList(shape, index)
1004                val li = Index.length head
1005                val lj = Index.length tail
1006                val c = Array.array(li * lj,init)
1007                fun loopi (0, _,  _)  = ()
1008                  | loopi (i, ia, ic) =
1009                    (Array.update(c, ic, f(Array.sub(c,ic), Array.sub(a,ia)));
1010                     loopi (i-1, ia+1, ic+1))
1011                fun loopk (0, ia, _)  = ia
1012                  | loopk (k, ia, ic) = (loopi (li, ia, ic);
1013                                         loopk (k-1, ia+li, ic))
1014                fun loopj (0, _,  _)  = ()
1015                  | loopj (j, ia, ic) = loopj (j-1, loopk(lk,ia,ic), ic+li)
1016            in
1017                loopj (lj, 0, 0);
1018                make'(head @ tail, c)
1019            end
1020
1021    end
1022    end (* Tensor *)
1023
1024
1025structure TensorSlice : TENSOR_SLICE =
1026    struct
1027        structure Tensor = Tensor
1028        structure Index  = Tensor.Index
1029        structure Range  = Range
1030           
1031        type index = Tensor.Index.t
1032        type range = Range.t
1033        type 'a tensor = 'a Tensor.tensor
1034
1035        type 'a slice = {range : range, shape: index, tensor : 'a tensor}
1036
1037        fun slice (rs,tensor) =
1038            let val r = (Range.ranges (Tensor.shape tensor) rs)
1039            in
1040                {range=r,
1041                 shape=(Range.shape r),
1042                 tensor=tensor}
1043            end
1044
1045        fun length ({range, shape, tensor}) = Range.length range
1046        fun base ({range, shape, tensor}) = tensor
1047        fun shape ({range, shape, tensor}) = Range.shape range
1048
1049        fun map f slice = 
1050        let
1051           val te  = #tensor slice
1052           val ra  = #range slice
1053           val fndx  = Range.first ra
1054           val arr = Array.array(length(slice),f (Tensor.sub(te,fndx)))
1055           val i   = ref 0
1056        in 
1057           Range.iteri (fn (ndx) => let val v = f (Tensor.sub (te,ndx)) in (Array.update (arr, !i, v); i := (!i + 1); true) end) ra;
1058           Array.vector arr
1059        end
1060
1061        fun app f (slice: 'a slice) = 
1062        let
1063           val te  = #tensor slice
1064           val ra  = #range slice
1065           val fndx  = Range.first ra
1066        in 
1067           Range.iteri (fn (ndx) => (f (Tensor.sub (te,ndx)); true)) ra; ()
1068        end
1069
1070        fun map2 f (sl1: 'a slice) (sl2: 'b slice) = 
1071        let
1072           val _    = if not ((#shape sl1) = (#shape sl2)) then raise Index.Shape else ()
1073           val te1  = #tensor sl1
1074           val te2  = #tensor sl2
1075           val ra1  = #range sl1
1076           val ra2  = #range sl2
1077           val fndx1  = Range.first ra1
1078           val fndx2  = Range.first ra2
1079           val arr   = Array.array(length(sl1),f (Tensor.sub(te1,fndx1),Tensor.sub(te2,fndx2)))
1080           val i     = ref 0
1081        in 
1082           Range.iteri2 (fn (ndx,ndx') => let val v = f (Tensor.sub (te1,ndx),Tensor.sub (te2,ndx')) in (Array.update (arr, !i, v); i := (!i + 1); true) end) (ra1,ra2);
1083           Array.vector arr
1084        end
1085
1086        fun foldl f init (slice: 'a slice) = 
1087        let
1088           val te     = #tensor slice
1089           val ra     = #range slice
1090           val ax     = ref init
1091        in 
1092           Range.iteri (fn (ndx) => let val ax' = f (Tensor.sub (te,ndx),!ax) in ((ax := ax'); true) end) ra;
1093           !ax
1094        end
1095
1096    end                               
1097
1098
1099(*
1100 Copyright (c) Juan Jose Garcia Ripoll.
1101 All rights reserved.
1102 
1103 Refer to the COPYRIGHT file for license conditions
1104*)
1105
1106(*
1107 MONO_TENSOR            - signature -
1108
1109 Monomorphic tensor of arbitrary data (not only numbers). Operations
1110 should be provided to run the data in several ways, according to one
1111 index.
1112
1113 type tensor
1114        The type of the tensor itself
1115 type elem
1116        The type of every element
1117 val storage = RowMajor | ColumnMajor
1118        RowMajor = data is stored in consecutive cells, first index
1119        varying fastest (FORTRAN convention)
1120        ColumnMajor = data is stored in consecutive cells, last
1121        index varying fastest (C,C++,Pascal,CommonLisp convention)
1122 new ([i1,...,in],init)
1123        Build a new tensor with n indices, each of sizes i1...in,
1124        filled with 'init'.
1125 fromArray (shape,data)
1126 fromList (shape,data)
1127        Use 'data' to fill a tensor of that shape. An exception is
1128        raised if 'data' is too large or too small to properly
1129        fill the vector. Later use of a 'data' array is disregarded
1130        -- one must think that the tensor now owns the array.
1131 length tensor
1132 rank tensor
1133 shape tensor
1134        Return the number of elements, the number of indices and
1135        the shape (size of each index) of the tensor.
1136 toArray tensor
1137        Return the data of the tensor in the form of an array.
1138        Mutation of this array may lead to unexpected behavior.
1139        The data in the array is stored according to `storage'.
1140
1141 sub (tensor,[i1,...,in])
1142 update (tensor,[i1,...,in],new_value)
1143        Access the element that is indexed by the numbers [i1,..,in]
1144
1145 map f a
1146 mapi f a
1147        Produce a new array by mapping the function sequentially
1148        as specified by 'storage', to each element of tensor 'a'.
1149        In 'mapi' the function receives a (indices,value) tuple,
1150        while in 'map' it only receives the value.
1151 app f a
1152 appi f a
1153        The same as 'map' and 'mapi' but the function 'f' outputs
1154        nothing and no new array is produced, i.e. one only seeks
1155        the side effect that 'f' may produce.
1156 map2 operation a b
1157        Apply function 'f' to pairs of elements of 'a' and 'b'
1158        and build a new tensor with the output. Both operands
1159        must have the same shape or an exception is raised.
1160        The procedure is sequential, as specified by 'storage'.
1161 foldl operation a n
1162        Fold-left the elements of tensor 'a' along the n-th
1163        index.
1164 all test a
1165 any test a
1166        Folded boolean tests on the elements of the tensor.
1167
1168 map', map2', foldl'
1169        Polymorphic versions of map, map2, foldl.
1170*)
1171
1172signature MONO_TENSOR =
1173    sig
1174        structure Array : MONO_ARRAY
1175        structure Index : INDEX
1176        type index = Index.t
1177        type elem
1178        type tensor
1179        type t = tensor
1180
1181        val new : index * elem -> tensor
1182        val tabulate : index * (index -> elem) -> tensor
1183        val length : tensor -> int
1184        val rank : tensor -> int
1185        val shape : tensor -> (index)
1186        val reshape : index -> tensor -> tensor
1187        val fromList : index * elem list -> tensor
1188        val fromArray : index * Array.array -> tensor
1189        val toArray : tensor -> Array.array
1190
1191        val sub : tensor * index -> elem
1192        val update : tensor * index * elem -> unit
1193        val map : (elem -> elem) -> tensor -> tensor
1194        val map2 : (elem * elem -> elem) -> tensor -> tensor -> tensor
1195        val app : (elem -> unit) -> tensor -> unit
1196        val appi : (int * elem -> unit) -> tensor -> unit
1197        val foldl : (elem * 'a -> 'a) -> 'a -> tensor -> tensor
1198        val foldln : (elem * elem -> elem) -> elem -> tensor -> int -> tensor
1199        val all : (elem -> bool) -> tensor -> bool
1200        val any : (elem -> bool) -> tensor -> bool
1201
1202        val map' : (elem -> 'a) -> tensor -> 'a Tensor.tensor
1203        val map2' : (elem * elem -> 'a) -> tensor -> tensor -> 'a Tensor.tensor
1204        val foldl' : ('a * elem -> 'a) -> 'a -> tensor -> int -> 'a Tensor.tensor
1205    end
1206
1207(*
1208 NUMBER         - Signature -
1209
1210 Guarantees a structure with a minimal number of mathematical operations
1211 so as to build an algebraic structure named Tensor.
1212 *)
1213
1214signature NUMBER =
1215    sig
1216        type t
1217        val zero : t
1218        val one  : t
1219        val ~ : t -> t
1220        val + : t * t -> t
1221        val - : t * t -> t
1222        val * : t * t -> t
1223        val / : t * t -> t
1224        val toString : t -> string
1225    end
1226
1227signature NUMBER =
1228    sig
1229        type t
1230        val zero : t
1231        val one : t
1232
1233        val + : t * t -> t
1234        val - : t * t -> t
1235        val * : t * t -> t
1236        val *+ : t * t * t -> t
1237        val *- : t * t * t -> t
1238        val ** : t * int -> t
1239
1240        val ~ : t -> t
1241        val abs : t -> t
1242        val signum : t -> t
1243
1244        val == : t * t -> bool
1245        val != : t * t -> bool
1246
1247        val toString : t -> string
1248        val fromInt : int -> t
1249        val scan : (char,'a) StringCvt.reader -> (t,'a) StringCvt.reader
1250    end
1251
1252signature INTEGRAL_NUMBER =
1253    sig
1254        include NUMBER
1255
1256        val quot : t * t -> t
1257        val rem  : t * t -> t
1258        val mod  : t * t -> t
1259        val div  : t * t -> t
1260
1261        val compare : t * t -> order
1262        val < : t * t -> bool
1263        val > : t * t -> bool
1264        val <= : t * t -> bool
1265        val >= : t * t -> bool
1266
1267        val max : t * t -> t
1268        val min : t * t -> t
1269    end
1270
1271signature FRACTIONAL_NUMBER =
1272    sig
1273        include NUMBER
1274
1275        val pi : t
1276        val e : t
1277
1278        val / : t * t -> t
1279        val recip : t -> t
1280
1281        val ln : t -> t
1282        val pow : t * t -> t
1283        val exp : t -> t
1284        val sqrt : t -> t
1285
1286        val cos : t -> t
1287        val sin : t -> t
1288        val tan : t -> t
1289        val sinh : t -> t
1290        val cosh : t -> t
1291        val tanh : t -> t
1292
1293        val acos : t -> t
1294        val asin : t -> t
1295        val atan : t -> t
1296        val asinh : t -> t
1297        val acosh : t -> t
1298        val atanh : t -> t
1299        val atan2 : t * t -> t
1300    end
1301
1302signature REAL_NUMBER =
1303    sig
1304        include FRACTIONAL_NUMBER
1305
1306        val compare : t * t -> order
1307        val < : t * t -> bool
1308        val > : t * t -> bool
1309        val <= : t * t -> bool
1310        val >= : t * t -> bool
1311
1312        val max : t * t -> t
1313        val min : t * t -> t
1314    end
1315
1316signature COMPLEX_NUMBER =
1317    sig
1318        include FRACTIONAL_NUMBER
1319
1320        structure Real : REAL_NUMBER
1321        type real = Real.t
1322
1323        val make : real * real -> t
1324        val split : t -> real * real
1325        val realPart : t -> real
1326        val imagPart : t -> real
1327        val abs2 : t -> real
1328    end
1329
1330structure INumber : INTEGRAL_NUMBER =
1331    struct
1332        open Int
1333        type t = Int.int
1334        val zero = 0
1335        val one = 1
1336
1337        infix **
1338        fun i ** n =
1339            let fun loop 0 = 1
1340                  | loop 1 = i
1341                  | loop n =
1342                let val x = loop (Int.div(n, 2))
1343                    val m = Int.mod(n, 2)
1344                in
1345                    if m = 0 then
1346                        x * x
1347                    else
1348                        x * x * i
1349                end
1350            in if n < 0
1351               then raise Domain
1352               else loop n
1353            end
1354
1355        fun signum i = case compare(i, 0) of
1356            GREATER => 1
1357          | EQUAL => 0
1358          | LESS => ~1
1359
1360        infix ==
1361        infix !=
1362        fun a == b = a = b
1363        fun a != b = (a <> b)
1364        fun *+(b,c,a) = b * c + a
1365        fun *-(b,c,a) = b * c - b
1366
1367        fun scan getc = Int.scan StringCvt.DEC getc
1368    end
1369
1370structure RNumber : REAL_NUMBER =
1371    struct
1372        open Real
1373        open Real.Math
1374        type t = Real.real
1375        val zero = 0.0
1376        val one = 1.0
1377
1378        fun signum x = case compare(x,0.0) of
1379            LESS => ~1.0
1380          | GREATER => 1.0
1381          | EQUAL => 0.0
1382
1383        fun recip x = 1.0 / x
1384
1385        infix **
1386        fun i ** n =
1387            let fun loop 0 = one
1388                  | loop 1 = i
1389                  | loop n =
1390                let val x = loop (Int.div(n, 2))
1391                    val m = Int.mod(n, 2)
1392                in
1393                    if m = 0 then
1394                        x * x
1395                    else
1396                        x * x * i
1397                end
1398            in if Int.<(n, 0)
1399               then raise Domain
1400               else loop n
1401            end
1402
1403        fun max (a, b) = if a < b then b else a
1404        fun min (a, b) = if a < b then a else b
1405
1406        fun asinh x = ln (x + sqrt(1.0 + x * x))
1407        fun acosh x = ln (x + (x + 1.0) * sqrt((x - 1.0)/(x + 1.0)))
1408        fun atanh x = ln ((1.0 + x) / sqrt(1.0 - x * x))
1409
1410    end
1411(*
1412 Complex(R)     - Functor -
1413
1414 Provides support for complex numbers based on tuples. Should be
1415 highly efficient as most operations can be inlined.
1416 *)
1417
1418structure CNumber : COMPLEX_NUMBER =
1419struct
1420        structure Real = RNumber
1421
1422        type t = Real.t * Real.t
1423        type real = Real.t
1424
1425        val zero = (0.0,0.0)
1426        val one = (1.0,0.0)
1427        val pi = (Real.pi, 0.0)
1428        val e = (Real.e, 0.0)
1429
1430        fun make (r,i) = (r,i) : t
1431        fun split z = z
1432        fun realPart (r,_) = r
1433        fun imagPart (_,i) = i
1434
1435        fun abs2 (r,i) = Real.+(Real.*(r,r),Real.*(i,i)) (* FIXME!!! *)
1436        fun arg (r,i) = Real.atan2(i,r)
1437        fun modulus z = Real.sqrt(abs2 z)
1438        fun abs z = (modulus z, 0.0)
1439        fun signum (z as (r,i)) =
1440            let val m = modulus z
1441            in (Real./(r,m), Real./(i,m))
1442            end
1443
1444        fun ~ (r1,i1) = (Real.~ r1, Real.~ i1)
1445        fun (r1,i1) + (r2,i2) = (Real.+(r1,r2), Real.+(i1,i2))
1446        fun (r1,i1) - (r2,i2) = (Real.-(r1,r2), Real.-(i1,i1))
1447        fun (r1,i1) * (r2,i2) = (Real.-(Real.*(r1,r2),Real.*(i1,i2)),
1448                                 Real.+(Real.*(r1,i2),Real.*(r2,i1)))
1449        fun (r1,i1) / (r2,i2) =
1450            let val modulus = abs2(r2,i2)
1451                val (nr,ni) = (r1,i1) * (r2,i2)
1452            in
1453                (Real./(nr,modulus), Real./(ni,modulus))
1454            end
1455        fun *+((r1,i1),(r2,i2),(r0,i0)) =
1456            (Real.*+(Real.~ i1, i2, Real.*+(r1,r2,r0)),
1457             Real.*+(r2, i2, Real.*+(r1,i2,i0)))
1458        fun *-((r1,i1),(r2,i2),(r0,i0)) =
1459            (Real.*+(Real.~ i1, i2, Real.*-(r1,r2,r0)),
1460             Real.*+(r2, i2, Real.*-(r1,i2,i0)))
1461
1462        infix **
1463        fun i ** n =
1464            let fun loop 0 = one
1465                  | loop 1 = i
1466                  | loop n =
1467                let val x = loop (Int.div(n, 2))
1468                    val m = Int.mod(n, 2)
1469                in
1470                    if m = 0 then
1471                        x * x
1472                    else
1473                        x * x * i
1474                end
1475            in if Int.<(n, 0)
1476                   then raise Domain
1477               else loop n
1478            end
1479
1480        fun recip (r1, i1) = 
1481            let val modulus = abs2(r1, i1)
1482            in (Real./(r1, modulus), Real./(Real.~ i1, modulus))
1483            end
1484        fun ==(z, w) = Real.==(realPart z, realPart w) andalso Real.==(imagPart z, imagPart w)
1485        fun !=(z, w) = Real.!=(realPart z, realPart w) andalso Real.!=(imagPart z, imagPart w)
1486        fun fromInt i = (Real.fromInt i, 0.0)
1487        fun toString (r,i) =
1488            String.concat ["(",Real.toString r,",",Real.toString i,")"]
1489
1490        fun exp (x, y) =
1491            let val expx = Real.exp x
1492            in (Real.*(x, (Real.cos y)), Real.*(x, (Real.sin y)))
1493            end
1494
1495    local
1496        val half = Real.recip (Real.fromInt 2)
1497    in
1498        fun sqrt (z as (x,y)) =
1499            if Real.==(x, 0.0) andalso Real.==(y, 0.0) then
1500                zero
1501            else
1502                let val m = Real.+(modulus z, Real.abs x)
1503                    val u' = Real.sqrt (Real.*(m, half))
1504                    val v' = Real./(Real.abs y , Real.+(u',u'))
1505                    val (u,v) = if Real.<(x, 0.0) then (v',u') else (u',v')
1506                in (u, if Real.<(y, 0.0) then Real.~ v else v)
1507                end
1508    end
1509        fun ln z = (Real.ln (modulus z), arg z)
1510
1511        fun pow (z, n) =
1512            let val l = ln z
1513            in exp (l * n)
1514            end
1515
1516        fun sin (x, y) = (Real.*(Real.sin x, Real.cosh y),
1517                          Real.*(Real.cos x, Real.sinh y))
1518        fun cos (x, y) = (Real.*(Real.cos x, Real.cosh y),
1519                          Real.~ (Real.*(Real.sin x, Real.sinh y)))
1520        fun tan (x, y) =
1521            let val (sx, cx) = (Real.sin x, Real.cos x)
1522                val (shy, chy) = (Real.sinh y, Real.cosh y)
1523                val a = (Real.*(sx, chy), Real.*(cx, shy))
1524                val b = (Real.*(cx, chy), Real.*(Real.~ sx, shy))
1525            in a / b
1526            end
1527
1528        fun sinh (x, y) = (Real.*(Real.cos y, Real.sinh x),
1529                           Real.*(Real.sin y, Real.cosh x))
1530        fun cosh (x, y) = (Real.*(Real.cos y, Real.cosh x),
1531                           Real.*(Real.sin y, Real.sinh x))
1532        fun tanh (x, y) =
1533            let val (sy, cy) = (Real.sin y, Real.cos y)
1534                val (shx, chx) = (Real.sinh x, Real.cosh x)
1535                val a = (Real.*(cy, shx), Real.*(sy, chx))
1536                val b = (Real.*(cy, chx), Real.*(sy, shx))
1537            in a / b
1538            end
1539
1540        fun asin (z as (x,y)) =
1541            let val w = sqrt (one - z * z)
1542                val (x',y') = ln ((Real.~ y, x) + w)
1543            in (y', Real.~ x')
1544            end
1545
1546        fun acos (z as (x,y)) = 
1547            let val (x', y') = sqrt (one + z * z)
1548                val (x'', y'') = ln (z + (Real.~ y', x'))
1549            in (y'', Real.~ x'')
1550            end
1551
1552        fun atan (z as (x,y)) =
1553            let val w = sqrt (one + z*z)
1554                val (x',y') = ln ((Real.-(1.0, y), x) / w)
1555            in (y', Real.~ x')
1556            end
1557
1558        fun atan2 (y, x) = atan(y / x)
1559
1560        fun asinh x = ln (x + sqrt(one + x * x))
1561        fun acosh x = ln (x + (x + one) * sqrt((x - one)/(x + one)))
1562        fun atanh x = ln ((one + x) / sqrt(one - x * x))
1563
1564        fun scan getc =
1565            let val scanner = Real.scan getc
1566            in fn stream => 
1567                  case scanner stream of
1568                      NONE => NONE
1569                    | SOME (a, rest) =>
1570                      case scanner rest of
1571                          NONE => NONE
1572                        | SOME (b, rest) => SOME (make(a,b), rest)
1573            end
1574
1575end (* ComplexNumber *)
1576
1577(*
1578 Copyright (c) Juan Jose Garcia Ripoll.
1579 All rights reserved.
1580 
1581 Refer to the COPYRIGHT file for license conditions
1582*)
1583
1584structure INumberArray =
1585    struct
1586        open Array
1587        type array = INumber.t array
1588        type vector = INumber.t vector
1589        type elem  = INumber.t
1590        structure Vector =
1591            struct
1592                open Vector
1593                type vector = INumber.t Vector.vector
1594                type elem = INumber.t
1595            end
1596        fun map f a = tabulate(length a, fn x => (f (sub(a,x))))
1597        fun mapi f a = tabulate(length a, fn x => (f (x,sub(a,x))))
1598        fun map2 f a b = tabulate(length a, fn x => (f(sub(a,x),sub(b,x))))
1599    end
1600
1601structure RNumberArray =
1602    struct
1603        open Real64Array
1604        val sub = Unsafe.Real64Array.sub
1605        val update = Unsafe.Real64Array.update
1606        fun map f a = tabulate(length a, fn x => (f (sub(a,x))))
1607        fun mapi f a = tabulate(length a, fn x => (f (x,sub(a,x))))
1608        fun map2 f a b = tabulate(length a, fn x => (f(sub(a,x),sub(b,x))))
1609    end
1610structure RNumber : REAL_NUMBER =
1611    struct
1612        open Real
1613        open Real.Math
1614        type t = Real.real
1615        val zero = 0.0
1616        val one = 1.0
1617        fun signum x = case compare(x,0.0) of
1618            LESS => ~1.0
1619          | GREATER => 1.0
1620          | EQUAL => 0.0
1621        fun recip x = 1.0 / x
1622        infix **
1623        fun i ** n =
1624            let fun loop 0 = one
1625                  | loop 1 = i
1626                  | loop n =
1627                let val x = loop (Int.div(n, 2))
1628                    val m = Int.mod(n, 2)
1629                in
1630                    if m = 0 then
1631                        x * x
1632                    else
1633                        x * x * i
1634                end
1635            in if Int.<(n, 0)
1636               then raise Domain
1637               else loop n
1638            end
1639        fun max (a, b) = if a < b then b else a
1640        fun min (a, b) = if a < b then a else b
1641        fun asinh x = ln (x + sqrt(1.0 + x * x))
1642        fun acosh x = ln (x + (x + 1.0) * sqrt((x - 1.0)/(x + 1.0)))
1643        fun atanh x = ln ((1.0 + x) / sqrt(1.0 - x * x))
1644    end
1645(*
1646 Complex(R)     - Functor -
1647 Provides support for complex numbers based on tuples. Should be
1648 highly efficient as most operations can be inlined.
1649 *)
1650structure CNumber : COMPLEX_NUMBER =
1651struct
1652        structure Real = RNumber
1653        type t = Real.t * Real.t
1654        type real = Real.t
1655        val zero = (0.0,0.0)
1656        val one = (1.0,0.0)
1657        val pi = (Real.pi, 0.0)
1658        val e = (Real.e, 0.0)
1659        fun make (r,i) = (r,i) : t
1660        fun split z = z
1661        fun realPart (r,_) = r
1662        fun imagPart (_,i) = i
1663        fun abs2 (r,i) = Real.+(Real.*(r,r),Real.*(i,i)) (* FIXME!!! *)
1664        fun arg (r,i) = Real.atan2(i,r)
1665        fun modulus z = Real.sqrt(abs2 z)
1666        fun abs z = (modulus z, 0.0)
1667        fun signum (z as (r,i)) =
1668            let val m = modulus z
1669            in (Real./(r,m), Real./(i,m))
1670            end
1671        fun ~ (r1,i1) = (Real.~ r1, Real.~ i1)
1672        fun (r1,i1) + (r2,i2) = (Real.+(r1,r2), Real.+(i1,i2))
1673        fun (r1,i1) - (r2,i2) = (Real.-(r1,r2), Real.-(i1,i1))
1674        fun (r1,i1) * (r2,i2) = (Real.-(Real.*(r1,r2),Real.*(i1,i2)),
1675                                 Real.+(Real.*(r1,i2),Real.*(r2,i1)))
1676        fun (r1,i1) / (r2,i2) =
1677            let val modulus = abs2(r2,i2)
1678                val (nr,ni) = (r1,i1) * (r2,i2)
1679            in
1680                (Real./(nr,modulus), Real./(ni,modulus))
1681            end
1682        fun *+((r1,i1),(r2,i2),(r0,i0)) =
1683            (Real.*+(Real.~ i1, i2, Real.*+(r1,r2,r0)),
1684             Real.*+(r2, i2, Real.*+(r1,i2,i0)))
1685        fun *-((r1,i1),(r2,i2),(r0,i0)) =
1686            (Real.*+(Real.~ i1, i2, Real.*-(r1,r2,r0)),
1687             Real.*+(r2, i2, Real.*-(r1,i2,i0)))
1688        infix **
1689        fun i ** n =
1690            let fun loop 0 = one
1691                  | loop 1 = i
1692                  | loop n =
1693                let val x = loop (Int.div(n, 2))
1694                    val m = Int.mod(n, 2)
1695                in
1696                    if m = 0 then
1697                        x * x
1698                    else
1699                        x * x * i
1700                end
1701            in if Int.<(n, 0)
1702                   then raise Domain
1703               else loop n
1704            end
1705        fun recip (r1, i1) = 
1706            let val modulus = abs2(r1, i1)
1707            in (Real./(r1, modulus), Real./(Real.~ i1, modulus))
1708            end
1709        fun ==(z, w) = Real.==(realPart z, realPart w) andalso Real.==(imagPart z, imagPart w)
1710        fun !=(z, w) = Real.!=(realPart z, realPart w) andalso Real.!=(imagPart z, imagPart w)
1711        fun fromInt i = (Real.fromInt i, 0.0)
1712        fun toString (r,i) =
1713            String.concat ["(",Real.toString r,",",Real.toString i,")"]
1714        fun exp (x, y) =
1715            let val expx = Real.exp x
1716            in (Real.*(x, (Real.cos y)), Real.*(x, (Real.sin y)))
1717            end
1718    local
1719        val half = Real.recip (Real.fromInt 2)
1720    in
1721        fun sqrt (z as (x,y)) =
1722            if Real.==(x, 0.0) andalso Real.==(y, 0.0) then
1723                zero
1724            else
1725                let val m = Real.+(modulus z, Real.abs x)
1726                    val u' = Real.sqrt (Real.*(m, half))
1727                    val v' = Real./(Real.abs y , Real.+(u',u'))
1728                    val (u,v) = if Real.<(x, 0.0) then (v',u') else (u',v')
1729                in (u, if Real.<(y, 0.0) then Real.~ v else v)
1730                end
1731    end
1732        fun ln z = (Real.ln (modulus z), arg z)
1733        fun pow (z, n) =
1734            let val l = ln z
1735            in exp (l * n)
1736            end
1737        fun sin (x, y) = (Real.*(Real.sin x, Real.cosh y),
1738                          Real.*(Real.cos x, Real.sinh y))
1739        fun cos (x, y) = (Real.*(Real.cos x, Real.cosh y),
1740                          Real.~ (Real.*(Real.sin x, Real.sinh y)))
1741        fun tan (x, y) =
1742            let val (sx, cx) = (Real.sin x, Real.cos x)
1743                val (shy, chy) = (Real.sinh y, Real.cosh y)
1744                val a = (Real.*(sx, chy), Real.*(cx, shy))
1745                val b = (Real.*(cx, chy), Real.*(Real.~ sx, shy))
1746            in a / b
1747            end
1748        fun sinh (x, y) = (Real.*(Real.cos y, Real.sinh x),
1749                           Real.*(Real.sin y, Real.cosh x))
1750        fun cosh (x, y) = (Real.*(Real.cos y, Real.cosh x),
1751                           Real.*(Real.sin y, Real.sinh x))
1752        fun tanh (x, y) =
1753            let val (sy, cy) = (Real.sin y, Real.cos y)
1754                val (shx, chx) = (Real.sinh x, Real.cosh x)
1755                val a = (Real.*(cy, shx), Real.*(sy, chx))
1756                val b = (Real.*(cy, chx), Real.*(sy, shx))
1757            in a / b
1758            end
1759        fun asin (z as (x,y)) =
1760            let val w = sqrt (one - z * z)
1761                val (x',y') = ln ((Real.~ y, x) + w)
1762            in (y', Real.~ x')
1763            end
1764        fun acos (z as (x,y)) = 
1765            let val (x', y') = sqrt (one + z * z)
1766                val (x'', y'') = ln (z + (Real.~ y', x'))
1767            in (y'', Real.~ x'')
1768            end
1769        fun atan (z as (x,y)) =
1770            let val w = sqrt (one + z*z)
1771                val (x',y') = ln ((Real.-(1.0, y), x) / w)
1772            in (y', Real.~ x')
1773            end
1774        fun atan2 (y, x) = atan(y / x)
1775        fun asinh x = ln (x + sqrt(one + x * x))
1776        fun acosh x = ln (x + (x + one) * sqrt((x - one)/(x + one)))
1777        fun atanh x = ln ((one + x) / sqrt(one - x * x))
1778        fun scan getc =
1779            let val scanner = Real.scan getc
1780            in fn stream => 
1781                  case scanner stream of
1782                      NONE => NONE
1783                    | SOME (a, rest) =>
1784                      case scanner rest of
1785                          NONE => NONE
1786                        | SOME (b, rest) => SOME (make(a,b), rest)
1787            end
1788end (* ComplexNumber *)
1789
1790(*
1791 Copyright (c) Juan Jose Garcia Ripoll.
1792 All rights reserved.
1793 Refer to the COPYRIGHT file for license conditions
1794*)
1795structure INumberArray =
1796    struct
1797        open Array
1798        type array = INumber.t array
1799        type vector = INumber.t vector
1800        type elem  = INumber.t
1801        structure Vector =
1802            struct
1803                open Vector
1804                type vector = INumber.t Vector.vector
1805                type elem = INumber.t
1806            end
1807        fun map f a = tabulate(length a, fn x => (f (sub(a,x))))
1808        fun mapi f a = tabulate(length a, fn x => (f (x,sub(a,x))))
1809        fun map2 f a b = tabulate(length a, fn x => (f(sub(a,x),sub(b,x))))
1810    end
1811structure RNumberArray =
1812    struct
1813        open Real64Array
1814        val sub = Unsafe.Real64Array.sub
1815        val update = Unsafe.Real64Array.update
1816        fun map f a = tabulate(length a, fn x => (f (sub(a,x))))
1817        fun mapi f a = tabulate(length a, fn x => (f (x,sub(a,x))))
1818        fun map2 f a b = tabulate(length a, fn x => (f(sub(a,x),sub(b,x))))
1819    end
1820(*--------------------- COMPLEX ARRAY -------------------------*)
1821structure BasicCNumberArray =
1822struct
1823        structure Complex : COMPLEX_NUMBER = CNumber
1824        structure Array : MONO_ARRAY = RNumberArray
1825        type elem = Complex.t
1826        type array = Array.array * Array.array
1827        val maxLen = Array.maxLen
1828        fun length (a,b) = Array.length a
1829        fun sub ((a,b),index) = Complex.make(Array.sub(a,index),Array.sub(b,index))
1830        fun update ((a,b),index,z) =
1831            let val (re,im) = Complex.split z in
1832                Array.update(a, index, re);
1833                Array.update(b, index, im)
1834            end
1835    local
1836        fun makeRange (a, start, NONE) = makeRange(a, start, SOME (length a - 1))
1837          | makeRange (a, start, SOME last) =
1838            let val len = length a
1839                val diff = last - start
1840            in
1841                if (start >= len) orelse (last >= len) then
1842                    raise Subscript
1843                else if diff < 0 then
1844                    (a, start, 0)
1845                else
1846                    (a, start, diff + 1)
1847            end
1848    in
1849        fun array (size,z:elem) =
1850            let val realsize = size * 2
1851                val r = Complex.realPart z
1852                val i = Complex.imagPart z in
1853                    (Array.array(size,r), Array.array(size,i))
1854            end
1855        fun zeroarray size =
1856            (Array.array(size,Complex.Real.zero),
1857             Array.array(size,Complex.Real.zero))
1858        fun tabulate (size,f) =
1859            let val a = array(size, Complex.zero)
1860                fun loop i =
1861                    case i = size of
1862                        true => a
1863                      | false => (update(a, i, f i); loop (i+1))
1864            in
1865                loop 0
1866            end
1867        fun fromList list =
1868            let val length = List.length list
1869                val a = zeroarray length
1870                fun loop (_, []) = a
1871                  | loop (i, z::rest) = (update(a, i, z);
1872                                         loop (i+1, rest))
1873            in
1874                loop(0,list)
1875            end
1876        fun extract range =
1877            let val (a, start, len) = makeRange range
1878                fun copy i = sub(a, i + start)
1879            in tabulate(len, copy)
1880            end
1881        fun concat array_list =
1882            let val total_length = foldl (op +) 0 (map length array_list)
1883                val a = array(total_length, Complex.zero)
1884                fun copy (_, []) = a
1885                  | copy (pos, v::rest) =
1886                    let fun loop i =
1887                        case i = 0 of
1888                            true => ()
1889                          | false => (update(a, i+pos, sub(v, i)); loop (i-1))
1890                    in (loop (length v - 1); copy(length v + pos, rest))
1891                    end
1892            in
1893                copy(0, array_list)
1894            end
1895        fun copy {src : array, si : int, len : int option, dst : array, di : int } =
1896            let val (a, ia, la) = makeRange (src, si, len)
1897                val (b, ib, lb) = makeRange (dst, di, len)
1898                fun copy i =
1899                    case i < 0 of
1900                        true => ()
1901                      | false => (update(b, i+ib, sub(a, i+ia)); copy (i-1))
1902            in copy (la - 1)
1903            end
1904        val copyVec = copy
1905        fun modifyi f range =
1906            let val (a, start, len) = makeRange range
1907                val last = start + len
1908                fun loop i =
1909                    case i >= last of
1910                        true => ()
1911                      | false => (update(a, i, f(i, sub(a,i))); loop (i+1))
1912            in loop start
1913            end
1914        fun modify f a =
1915            let val last = length a
1916                fun loop i =
1917                    case i >= last of
1918                        true => ()
1919                      | false => (update(a, i, f(sub(a,i))); loop (i+1))
1920            in loop 0
1921            end
1922        fun app f a =
1923            let val size = length a
1924                fun loop i =
1925                    case i = size of
1926                        true => ()
1927                      | false => (f(sub(a,i)); loop (i+1))
1928            in
1929                loop 0
1930            end
1931        fun appi f range =
1932            let val (a, start, len) = makeRange range
1933                val last = start + len
1934                fun loop i =
1935                    case i >= last of
1936                        true => ()
1937                      | false => (f(i, sub(a,i)); loop (i+1))
1938            in
1939                loop start
1940            end
1941        fun map f a =
1942            let val len = length a
1943                val c = zeroarray len
1944                fun loop ~1 = c
1945                  | loop i = (update(a, i, f(sub(a,i))); loop (i-1))
1946            in loop (len-1)
1947            end
1948        fun map2 f a b =
1949            let val len = length a
1950                val c = zeroarray len
1951                fun loop ~1 = c
1952                  | loop i = (update(c, i, f(sub(a,i),sub(b,i)));
1953                              loop (i-1))
1954            in loop (len-1)
1955            end
1956        fun mapi f range =
1957            let val (a, start, len) = makeRange range
1958                fun rule i = f (i+start, sub(a, i+start))
1959            in tabulate(len, rule)
1960            end
1961        fun foldli f init range =
1962            let val (a, start, len) = makeRange range
1963                val last = start + len - 1
1964                fun loop (i, accum) =
1965                    case i > last of
1966                        true => accum
1967                      | false => loop (i+1, f(i, sub(a,i), accum))
1968            in loop (start, init)
1969            end
1970        fun foldri f init range =
1971            let val (a, start, len) = makeRange range
1972                val last = start + len - 1
1973                fun loop (i, accum) =
1974                    case i < start of
1975                        true => accum
1976                      | false => loop (i-1, f(i, sub(a,i), accum))
1977            in loop (last, init)
1978            end
1979        fun foldl f init a = foldli (fn (_, a, x) => f(a,x)) init (a,0,NONE)
1980        fun foldr f init a = foldri (fn (_, x, a) => f(x,a)) init (a,0,NONE)
1981    end
1982end (* BasicCNumberArray *)
1983structure CNumberArray =
1984    struct
1985        structure Vector =
1986            struct
1987                open BasicCNumberArray
1988                type vector = array
1989            end : MONO_VECTOR
1990        type vector = Vector.vector
1991        open BasicCNumberArray
1992    end (* CNumberArray *)
1993structure ITensor =
1994    struct
1995        structure Number = INumber
1996        structure Array = INumberArray
1997(*
1998 Copyright (c) Juan Jose Garcia Ripoll.
1999 All rights reserved.
2000 Refer to the COPYRIGHT file for license conditions
2001*)
2002structure MonoTensor  =
2003    struct
2004(* PARAMETERS
2005        structure Array = Array
2006*)
2007        structure Index  = Index
2008        type elem = Array.elem
2009        type index = Index.t
2010        type tensor = {shape : index, indexer : Index.indexer, data : Array.array}
2011        type t = tensor
2012        exception Shape
2013        exception Match
2014        exception Index
2015    local
2016    (*----- LOCALS -----*)
2017        fun make' (shape, data) =
2018            {shape = shape, indexer = Index.indexer shape, data = data}
2019        fun toInt {shape, indexer, data} index = indexer index
2020        fun splitList (l as (a::rest), place) =
2021            let fun loop (left,here,right) 0 =  (List.rev left,here,right)
2022                  | loop (_,_,[]) place = raise Index
2023                  | loop (left,here,a::right) place = 
2024                loop (here::left,a,right) (place-1)
2025            in
2026                if place <= 0 then
2027                    loop ([],a,rest) (List.length rest - place)
2028                else
2029                    loop ([],a,rest) (place - 1)
2030            end
2031    in
2032    (*----- STRUCTURAL OPERATIONS & QUERIES ------*)
2033      fun cat (dim, x: tensor, y: tensor) =
2034        (let val xshape = (#shape x)
2035             val yshape = (#shape y)
2036             val xdata  = (#data x)
2037             val ydata  = (#data y)
2038        in
2039           if  not (length xshape  = length yshape) then
2040           raise Shape
2041           else
2042                let 
2043                   val (_,newshape)   = ListPair.foldl
2044                                      (fn (x,y,(i,ax)) => if (dim = i) then (i+1,(x+y) :: ax) 
2045                                                                       else if not (x=y) then raise Shape else (i+1,x :: ax))
2046                                       (0,[]) (xshape, yshape)
2047                   val newlength  = Index.length newshape
2048                   val newdata    = Array.array(newlength,Array.sub(xdata,0))
2049                in
2050                    Array.copy {src=xdata,dst=newdata,di=0};
2051                    Array.copy {src=ydata,dst=newdata,di=(Index.length xshape)};
2052                    {shape = newshape,
2053                     indexer = Index.indexer newshape,
2054                     data = newdata}
2055                end
2056        end)
2057
2058        fun new (shape, init) =
2059            if not (Index.validShape shape) then
2060                raise Shape
2061            else
2062                let val length = Index.length shape in
2063                    {shape = shape,
2064                     indexer = Index.indexer shape,
2065                     data = Array.array(length,init)}
2066                end
2067        fun toArray {shape, indexer, data} = data
2068        fun length {shape, indexer, data} =  Array.length data
2069        fun shape {shape, indexer, data} = shape
2070        fun rank t = List.length (shape t)
2071        fun reshape new_shape tensor =
2072            if Index.validShape new_shape then
2073                case (Index.length new_shape) = length tensor of
2074                    true => make'(new_shape, toArray tensor)
2075                  | false => raise Match
2076            else
2077                raise Shape
2078        fun fromArray (s, a) =
2079            case Index.validShape s andalso 
2080                 ((Index.length s) = (Array.length a)) of
2081                 true => make'(s, a)
2082               | false => raise Shape
2083        fun fromList (s, a) = fromArray (s, Array.fromList a)
2084        fun tabulate (shape,f) =
2085            if Index.validShape shape then
2086                let val last = Index.last shape
2087                    val length = Index.length shape
2088                    val c = Array.array(length, f last)
2089                    fun dotable (c, indices, i) =
2090                        (Array.update(c, i, f indices);
2091                         if i <= 1
2092                         then c
2093                         else dotable(c, Index.prev' shape indices, i-1))
2094                in make'(shape,dotable(c, Index.prev' shape last, length-2))
2095                end
2096            else
2097                raise Shape
2098        (*----- ELEMENTWISE OPERATIONS -----*)
2099        fun sub (t, index) = Array.sub(#data t, toInt t index)
2100        fun update (t, index, value) =
2101            Array.update(toArray t, toInt t index, value)
2102        fun map f {shape, indexer, data} =
2103            {shape = shape, indexer = indexer, data = Array.map f data}
2104        fun map2 f t1 t2=
2105            let val {shape=shape1, indexer=indexer1, data=data1} = t1
2106                val {shape=shape2, indexer=indexer2, data=data2} = t2
2107            in
2108                if Index.eq(shape1,shape2) then
2109                    {shape = shape1,
2110                     indexer = indexer1,
2111                     data = Array.map2 f data1 data2}
2112                else
2113                    raise Match
2114        end
2115        fun appi f tensor = Array.appi f (toArray tensor)
2116        fun app f tensor = Array.app f (toArray tensor)
2117        fun all f tensor =
2118            let val a = toArray tensor
2119            in Loop.all(0, length tensor - 1, fn i =>
2120                        f (Array.sub(a, i)))
2121            end
2122        fun any f tensor =
2123            let val a = toArray tensor
2124            in Loop.any(0, length tensor - 1, fn i =>
2125                        f (Array.sub(a, i)))
2126            end
2127        fun foldl f init tensor = Array.foldl f init (toArray tensor)
2128        fun foldln f init {shape, indexer, data=a} index =
2129            let val (head,lk,tail) = splitList(shape, index)
2130                val li = Index.length head
2131                val lj = Index.length tail
2132                val c = Array.array(li * lj,init)
2133                fun loopi (0, _,  _)  = ()
2134                  | loopi (i, ia, ic) =
2135                    (Array.update(c, ic, f(Array.sub(c,ic), Array.sub(a,ia)));
2136                     loopi (i-1, ia+1, ic+1))
2137                fun loopk (0, ia, _)  = ia
2138                  | loopk (k, ia, ic) = (loopi (li, ia, ic);
2139                                         loopk (k-1, ia+li, ic))
2140                fun loopj (0, _,  _)  = ()
2141                  | loopj (j, ia, ic) = loopj (j-1, loopk(lk,ia,ic), ic+li)
2142            in
2143                loopj (lj, 0, 0);
2144                make'(head @ tail, c)
2145            end
2146        (* --- POLYMORPHIC ELEMENTWISE OPERATIONS --- *)
2147        fun array_map' f a =
2148            let fun apply index = f(Array.sub(a,index)) in
2149                Tensor.Array.tabulate(Array.length a, apply)
2150            end
2151        fun map' f t = Tensor.fromArray(shape t, array_map' f (toArray t))
2152        fun map2' f t1 t2 =
2153            let val d1 = toArray t1
2154                val d2 = toArray t2
2155                fun apply i = f (Array.sub(d1,i), Array.sub(d2,i))
2156                val len = Array.length d1
2157            in
2158                if Index.eq(shape t1, shape t2) then
2159                    Tensor.fromArray(shape t1, Tensor.Array.tabulate(len,apply))
2160                else
2161                    raise Match
2162            end
2163        fun foldl' f init {shape, indexer, data=a} index =
2164            let val (head,lk,tail) = splitList(shape, index)
2165                val li = Index.length head
2166                val lj = Index.length tail
2167                val c = Tensor.Array.array(li * lj,init)
2168                fun loopi (0, _,  _)  = ()
2169                  | loopi (i, ia, ic) =
2170                    (Tensor.Array.update(c,ic,f(Tensor.Array.sub(c,ic),Array.sub(a,ia)));
2171                     loopi (i-1, ia+1, ic+1))
2172                fun loopk (0, ia, _)  = ia
2173                  | loopk (k, ia, ic) = (loopi (li, ia, ic);
2174                                         loopk (k-1, ia+li, ic))
2175                fun loopj (0, _,  _)  = ()
2176                  | loopj (j, ia, ic) = loopj (j-1, loopk(lk,ia,ic), ic+li)
2177            in
2178                loopj (lj, 0, 0);
2179                make'(head @ tail, c)
2180            end
2181    end
2182    end (* MonoTensor *)
2183        open MonoTensor
2184    local
2185        (*
2186         LEFT INDEX CONTRACTION:
2187         a = a(i1,i2,...,in)
2188         b = b(j1,j2,...,jn)
2189         c = c(i2,...,in,j2,...,jn)
2190         = sum(a(k,i2,...,jn)*b(k,j2,...jn)) forall k
2191         MEANINGFUL VARIABLES:
2192         lk = i1 = j1
2193         li = i2*...*in
2194         lj = j2*...*jn
2195         *)
2196        fun do_fold_first a b c lk lj li =
2197            let fun loopk (0, _,  _,  accum) = accum
2198                  | loopk (k, ia, ib, accum) =
2199                    let val delta = Number.*(Array.sub(a,ia),Array.sub(b,ib))
2200                    in loopk (k-1, ia+1, ib+1, Number.+(delta,accum))
2201                    end
2202                fun loopj (0, ib, ic) = c
2203                  | loopj (j, ib, ic) =
2204                    let fun loopi (0, ia, ic) = ic
2205                          | loopi (i, ia, ic) =
2206                        (Array.update(c, ic, loopk(lk, ia, ib, Number.zero));
2207                         loopi(i-1, ia+lk, ic+1))
2208                    in
2209                        loopj(j-1, ib+lk, loopi(li, 0, ic))
2210                    end
2211            in loopj(lj, 0, 0)
2212            end
2213    in
2214        fun +* ta tb =
2215            let val (rank_a,lk::rest_a,a) = (rank ta, shape ta, toArray ta)
2216                val (rank_b,lk2::rest_b,b) = (rank tb, shape tb, toArray tb)
2217            in if not(lk = lk2)
2218               then raise Match
2219               else let val li = Index.length rest_a
2220                        val lj = Index.length rest_b
2221                        val c = Array.array(li*lj,Number.zero)
2222                    in fromArray(rest_a @ rest_b,
2223                                 do_fold_first a b c lk li lj)
2224                    end
2225            end
2226    end
2227    local
2228        (*
2229         LAST INDEX CONTRACTION:
2230         a = a(i1,i2,...,in)
2231         b = b(j1,j2,...,jn)
2232         c = c(i2,...,in,j2,...,jn)
2233         = sum(mult(a(i1,i2,...,k),b(j1,j2,...,k))) forall k
2234         MEANINGFUL VARIABLES:
2235         lk = in = jn
2236         li = i1*...*i(n-1)
2237         lj = j1*...*j(n-1)
2238         *)
2239        fun do_fold_last a b c lk lj li =
2240            let fun loopi (0, ia, ic, fac) = ()
2241                  | loopi (i, ia, ic, fac) =
2242                    let val old = Array.sub(c,ic)
2243                        val inc = Number.*(Array.sub(a,ia),fac)
2244                    in
2245                        Array.update(c,ic,Number.+(old,inc));
2246                        loopi(i-1, ia+1, ic+1, fac)
2247                    end
2248                fun loopj (j, ib, ic) =
2249                    let fun loopk (0, ia, ib) = ()
2250                          | loopk (k, ia, ib) =
2251                            (loopi(li, ia, ic, Array.sub(b,ib));
2252                             loopk(k-1, ia+li, ib+lj))
2253                    in case j of
2254                           0 => c
2255                         | _ => (loopk(lk, 0, ib);
2256                                 loopj(j-1, ib+1, ic+li))
2257                    end (* loopj *)
2258            in
2259                loopj(lj, 0, 0)
2260            end
2261    in
2262        fun *+ ta tb  =
2263            let val (rank_a,shape_a,a) = (rank ta, shape ta, toArray ta)
2264                val (rank_b,shape_b,b) = (rank tb, shape tb, toArray tb)
2265                val (lk::rest_a) = List.rev shape_a
2266                val (lk2::rest_b) = List.rev shape_b
2267            in if not(lk = lk2)
2268               then raise Match
2269               else let val li = Index.length rest_a
2270                        val lj = Index.length rest_b
2271                        val c = Array.array(li*lj,Number.zero)
2272                    in fromArray(List.rev rest_a @ List.rev rest_b,
2273                                 do_fold_last a b c lk li lj)
2274                    end
2275            end
2276    end
2277        (* ALGEBRAIC OPERATIONS *)
2278        infix **
2279        infix ==
2280        infix !=
2281        fun a + b = map2 Number.+ a b
2282        fun a - b = map2 Number.- a b
2283        fun a * b = map2 Number.* a b
2284        fun a ** i = map (fn x => (Number.**(x,i))) a
2285        fun ~ a = map Number.~ a
2286        fun abs a = map Number.abs a
2287        fun signum a = map Number.signum a
2288        fun a == b = map2' Number.== a b
2289        fun a != b = map2' Number.!= a b
2290        fun toString a = raise Domain
2291        fun fromInt a = new([1], Number.fromInt a)
2292        (* TENSOR SPECIFIC OPERATIONS *)
2293        fun *> n = map (fn x => Number.*(n,x))
2294        fun normInf a =
2295            let fun accum (y,x) = Number.max(x,Number.abs y)
2296            in  foldl accum Number.zero a
2297            end
2298    end (* NumberTensor *)
2299
2300structure RTensor =
2301    struct
2302        structure Number = RNumber
2303        structure Array = RNumberArray
2304(*
2305 Copyright (c) Juan Jose Garcia Ripoll.
2306 All rights reserved.
2307 Refer to the COPYRIGHT file for license conditions
2308*)
2309structure MonoTensor  =
2310    struct
2311(* PARAMETERS
2312        structure Array = Array
2313*)
2314        structure Index  = Index
2315        type elem = Array.elem
2316        type index = Index.t
2317        type tensor = {shape : index, indexer : Index.indexer, data : Array.array}
2318        type t = tensor
2319        exception Shape
2320        exception Match
2321        exception Index
2322    local
2323    (*----- LOCALS -----*)
2324        fun make' (shape, data) =
2325            {shape = shape, indexer = Index.indexer shape, data = data}
2326        fun toInt {shape, indexer, data} index = indexer index
2327        fun splitList (l as (a::rest), place) =
2328            let fun loop (left,here,right) 0 =  (List.rev left,here,right)
2329                  | loop (_,_,[]) place = raise Index
2330                  | loop (left,here,a::right) place = 
2331                loop (here::left,a,right) (place-1)
2332            in
2333                if place <= 0 then
2334                    loop ([],a,rest) (List.length rest - place)
2335                else
2336                    loop ([],a,rest) (place - 1)
2337            end
2338    in
2339    (*----- STRUCTURAL OPERATIONS & QUERIES ------*)
2340      fun cat (dim, x: tensor, y: tensor) =
2341        (let val xshape = (#shape x)
2342             val yshape = (#shape y)
2343             val xdata  = (#data x)
2344             val ydata  = (#data y)
2345        in
2346           if  not (length xshape  = length yshape) then
2347           raise Shape
2348           else
2349                let 
2350                   val (_,newshape)   = ListPair.foldl
2351                                      (fn (x,y,(i,ax)) => if (dim = i) then (i+1,(x+y) :: ax) 
2352                                                                       else if not (x=y) then raise Shape else (i+1,x :: ax))
2353                                       (0,[]) (xshape, yshape)
2354                   val newlength  = Index.length newshape
2355                   val newdata    = Array.array(newlength,Array.sub(xdata,0))
2356                in
2357                    Array.copy {src=xdata,dst=newdata,di=0};
2358                    Array.copy {src=ydata,dst=newdata,di=(Index.length xshape)};
2359                    {shape = newshape,
2360                     indexer = Index.indexer newshape,
2361                     data = newdata}
2362                end
2363        end)
2364
2365        fun new (shape, init) =
2366            if not (Index.validShape shape) then
2367                raise Shape
2368            else
2369                let val length = Index.length shape in
2370                    {shape = shape,
2371                     indexer = Index.indexer shape,
2372                     data = Array.array(length,init)}
2373                end
2374        fun toArray {shape, indexer, data} = data
2375        fun length {shape, indexer, data} =  Array.length data
2376        fun shape {shape, indexer, data} = shape
2377        fun rank t = List.length (shape t)
2378        fun reshape new_shape tensor =
2379            if Index.validShape new_shape then
2380                case (Index.length new_shape) = length tensor of
2381                    true => make'(new_shape, toArray tensor)
2382                  | false => raise Match
2383            else
2384                raise Shape
2385        fun fromArray (s, a) =
2386            case Index.validShape s andalso 
2387                 ((Index.length s) = (Array.length a)) of
2388                 true => make'(s, a)
2389               | false => raise Shape
2390        fun fromList (s, a) = fromArray (s, Array.fromList a)
2391        fun tabulate (shape,f) =
2392            if Index.validShape shape then
2393                let val last = Index.last shape
2394                    val length = Index.length shape
2395                    val c = Array.array(length, f last)
2396                    fun dotable (c, indices, i) =
2397                        (Array.update(c, i, f indices);
2398                         if i <= 1
2399                         then c
2400                         else dotable(c, Index.prev' shape indices, i-1))
2401                in make'(shape,dotable(c, Index.prev' shape last, length-2))
2402                end
2403            else
2404                raise Shape
2405        (*----- ELEMENTWISE OPERATIONS -----*)
2406        fun sub (t, index) = Array.sub(#data t, toInt t index)
2407        fun update (t, index, value) =
2408            Array.update(toArray t, toInt t index, value)
2409        fun map f {shape, indexer, data} =
2410            {shape = shape, indexer = indexer, data = Array.map f data}
2411        fun map2 f t1 t2=
2412            let val {shape=shape1, indexer=indexer1, data=data1} = t1
2413                val {shape=shape2, indexer=indexer2, data=data2} = t2
2414            in
2415                if Index.eq(shape1,shape2) then
2416                    {shape = shape1,
2417                     indexer = indexer1,
2418                     data = Array.map2 f data1 data2}
2419                else
2420                    raise Match
2421        end
2422        fun appi f tensor = Array.appi f (toArray tensor)
2423        fun app f tensor = Array.app f (toArray tensor)
2424        fun all f tensor =
2425            let val a = toArray tensor
2426            in Loop.all(0, length tensor - 1, fn i =>
2427                        f (Array.sub(a, i)))
2428            end
2429        fun any f tensor =
2430            let val a = toArray tensor
2431            in Loop.any(0, length tensor - 1, fn i =>
2432                        f (Array.sub(a, i)))
2433            end
2434        fun foldl f init tensor = Array.foldl f init (toArray tensor)
2435        fun foldln f init {shape, indexer, data=a} index =
2436            let val (head,lk,tail) = splitList(shape, index)
2437                val li = Index.length head
2438                val lj = Index.length tail
2439                val c = Array.array(li * lj,init)
2440                fun loopi (0, _,  _)  = ()
2441                  | loopi (i, ia, ic) =
2442                    (Array.update(c, ic, f(Array.sub(c,ic), Array.sub(a,ia)));
2443                     loopi (i-1, ia+1, ic+1))
2444                fun loopk (0, ia, _)  = ia
2445                  | loopk (k, ia, ic) = (loopi (li, ia, ic);
2446                                         loopk (k-1, ia+li, ic))
2447                fun loopj (0, _,  _)  = ()
2448                  | loopj (j, ia, ic) = loopj (j-1, loopk(lk,ia,ic), ic+li)
2449            in
2450                loopj (lj, 0, 0);
2451                make'(head @ tail, c)
2452            end
2453        (* --- POLYMORPHIC ELEMENTWISE OPERATIONS --- *)
2454        fun array_map' f a =
2455            let fun apply index = f(Array.sub(a,index)) in
2456                Tensor.Array.tabulate(Array.length a, apply)
2457            end
2458        fun map' f t = Tensor.fromArray(shape t, array_map' f (toArray t))
2459        fun map2' f t1 t2 =
2460            let val d1 = toArray t1
2461                val d2 = toArray t2
2462                fun apply i = f (Array.sub(d1,i), Array.sub(d2,i))
2463                val len = Array.length d1
2464            in
2465                if Index.eq(shape t1, shape t2) then
2466                    Tensor.fromArray(shape t1, Tensor.Array.tabulate(len,apply))
2467                else
2468                    raise Match
2469            end
2470        fun foldl' f init {shape, indexer, data=a} index =
2471            let val (head,lk,tail) = splitList(shape, index)
2472                val li = Index.length head
2473                val lj = Index.length tail
2474                val c = Tensor.Array.array(li * lj,init)
2475                fun loopi (0, _,  _)  = ()
2476                  | loopi (i, ia, ic) =
2477                    (Tensor.Array.update(c,ic,f(Tensor.Array.sub(c,ic),Array.sub(a,ia)));
2478                     loopi (i-1, ia+1, ic+1))
2479                fun loopk (0, ia, _)  = ia
2480                  | loopk (k, ia, ic) = (loopi (li, ia, ic);
2481                                         loopk (k-1, ia+li, ic))
2482                fun loopj (0, _,  _)  = ()
2483                  | loopj (j, ia, ic) = loopj (j-1, loopk(lk,ia,ic), ic+li)
2484            in
2485                loopj (lj, 0, 0);
2486                make'(head @ tail, c)
2487            end
2488    end
2489    end (* MonoTensor *)
2490        open MonoTensor
2491    local
2492        (*
2493         LEFT INDEX CONTRACTION:
2494         a = a(i1,i2,...,in)
2495         b = b(j1,j2,...,jn)
2496         c = c(i2,...,in,j2,...,jn)
2497         = sum(a(k,i2,...,jn)*b(k,j2,...jn)) forall k
2498         MEANINGFUL VARIABLES:
2499         lk = i1 = j1
2500         li = i2*...*in
2501         lj = j2*...*jn
2502         *)
2503        fun do_fold_first a b c lk lj li =
2504            let fun loopk (0, _,  _,  accum) = accum
2505                  | loopk (k, ia, ib, accum) =
2506                    let val delta = Number.*(Array.sub(a,ia),Array.sub(b,ib))
2507                    in loopk (k-1, ia+1, ib+1, Number.+(delta,accum))
2508                    end
2509                fun loopj (0, ib, ic) = c
2510                  | loopj (j, ib, ic) =
2511                    let fun loopi (0, ia, ic) = ic
2512                          | loopi (i, ia, ic) =
2513                        (Array.update(c, ic, loopk(lk, ia, ib, Number.zero));
2514                         loopi(i-1, ia+lk, ic+1))
2515                    in
2516                        loopj(j-1, ib+lk, loopi(li, 0, ic))
2517                    end
2518            in loopj(lj, 0, 0)
2519            end
2520    in
2521        fun +* ta tb =
2522            let val (rank_a,lk::rest_a,a) = (rank ta, shape ta, toArray ta)
2523                val (rank_b,lk2::rest_b,b) = (rank tb, shape tb, toArray tb)
2524            in if not(lk = lk2)
2525               then raise Match
2526               else let val li = Index.length rest_a
2527                        val lj = Index.length rest_b
2528                        val c = Array.array(li*lj,Number.zero)
2529                    in fromArray(rest_a @ rest_b,
2530                                 do_fold_first a b c lk li lj)
2531                    end
2532            end
2533    end
2534    local
2535        (*
2536         LAST INDEX CONTRACTION:
2537         a = a(i1,i2,...,in)
2538         b = b(j1,j2,...,jn)
2539         c = c(i2,...,in,j2,...,jn)
2540         = sum(mult(a(i1,i2,...,k),b(j1,j2,...,k))) forall k
2541         MEANINGFUL VARIABLES:
2542         lk = in = jn
2543         li = i1*...*i(n-1)
2544         lj = j1*...*j(n-1)
2545         *)
2546        fun do_fold_last a b c lk lj li =
2547            let fun loopi (0, ia, ic, fac) = ()
2548                  | loopi (i, ia, ic, fac) =
2549                    let val old = Array.sub(c,ic)
2550                        val inc = Number.*(Array.sub(a,ia),fac)
2551                    in
2552                        Array.update(c,ic,Number.+(old,inc));
2553                        loopi(i-1, ia+1, ic+1, fac)
2554                    end
2555                fun loopj (j, ib, ic) =
2556                    let fun loopk (0, ia, ib) = ()
2557                          | loopk (k, ia, ib) =
2558                            (loopi(li, ia, ic, Array.sub(b,ib));
2559                             loopk(k-1, ia+li, ib+lj))
2560                    in case j of
2561                           0 => c
2562                         | _ => (loopk(lk, 0, ib);
2563                                 loopj(j-1, ib+1, ic+li))
2564                    end (* loopj *)
2565            in
2566                loopj(lj, 0, 0)
2567            end
2568    in
2569        fun *+ ta tb  =
2570            let val (rank_a,shape_a,a) = (rank ta, shape ta, toArray ta)
2571                val (rank_b,shape_b,b) = (rank tb, shape tb, toArray tb)
2572                val (lk::rest_a) = List.rev shape_a
2573                val (lk2::rest_b) = List.rev shape_b
2574            in if not(lk = lk2)
2575               then raise Match
2576               else let val li = Index.length rest_a
2577                        val lj = Index.length rest_b
2578                        val c = Array.array(li*lj,Number.zero)
2579                    in fromArray(List.rev rest_a @ List.rev rest_b,
2580                                 do_fold_last a b c lk li lj)
2581                    end
2582            end
2583    end
2584        (* ALGEBRAIC OPERATIONS *)
2585        infix **
2586        infix ==
2587        infix !=
2588        fun a + b = map2 Number.+ a b
2589        fun a - b = map2 Number.- a b
2590        fun a * b = map2 Number.* a b
2591        fun a ** i = map (fn x => (Number.**(x,i))) a
2592        fun ~ a = map Number.~ a
2593        fun abs a = map Number.abs a
2594        fun signum a = map Number.signum a
2595        fun a == b = map2' Number.== a b
2596        fun a != b = map2' Number.!= a b
2597        fun toString a = raise Domain
2598        fun fromInt a = new([1], Number.fromInt a)
2599        (* TENSOR SPECIFIC OPERATIONS *)
2600        fun *> n = map (fn x => Number.*(n,x))
2601        fun a / b = map2 Number./ a b
2602        fun recip a = map Number.recip a
2603        fun ln a = map Number.ln a
2604        fun pow (a, b) = map (fn x => (Number.pow(x,b))) a
2605        fun exp a = map Number.exp a
2606        fun sqrt a = map Number.sqrt a
2607        fun cos a = map Number.cos a
2608        fun sin a = map Number.sin a
2609        fun tan a = map Number.tan a
2610        fun sinh a = map Number.sinh a
2611        fun cosh a = map Number.cosh a
2612        fun tanh a = map Number.tanh a
2613        fun asin a = map Number.asin a
2614        fun acos a = map Number.acos a
2615        fun atan a = map Number.atan a
2616        fun asinh a = map Number.asinh a
2617        fun acosh a = map Number.acosh a
2618        fun atanh a = map Number.atanh a
2619        fun atan2 (a,b) = map2 Number.atan2 a b
2620        fun normInf a =
2621            let fun accum (y,x) = Number.max(x,Number.abs y)
2622            in  foldl accum Number.zero a
2623            end
2624        fun norm1 a =
2625            let fun accum (y,x) = Number.+(x,Number.abs y)
2626            in  foldl accum Number.zero a
2627            end
2628        fun norm2 a =
2629            let fun accum (y,x) = Number.+(x, Number.*(y,y))
2630            in Number.sqrt(foldl accum Number.zero a)
2631            end
2632    end (* RTensor *)
2633structure CTensor =
2634struct
2635    structure Number = CNumber
2636    structure Array = CNumberArray
2637(*
2638 Copyright (c) Juan Jose Garcia Ripoll.
2639 All rights reserved.
2640 Refer to the COPYRIGHT file for license conditions
2641*)
2642structure MonoTensor  =
2643    struct
2644(* PARAMETERS
2645        structure Array = Array
2646*)
2647        structure Index  = Index
2648        type elem = Array.elem
2649        type index = Index.t
2650        type tensor = {shape : index, indexer : Index.indexer, data : Array.array}
2651        type t = tensor
2652        exception Shape
2653        exception Match
2654        exception Index
2655    local
2656    (*----- LOCALS -----*)
2657        fun make' (shape, data) =
2658            {shape = shape, indexer = Index.indexer shape, data = data}
2659        fun toInt {shape, indexer, data} index = indexer index
2660        fun splitList (l as (a::rest), place) =
2661            let fun loop (left,here,right) 0 =  (List.rev left,here,right)
2662                  | loop (_,_,[]) place = raise Index
2663                  | loop (left,here,a::right) place = 
2664                loop (here::left,a,right) (place-1)
2665            in
2666                if place <= 0 then
2667                    loop ([],a,rest) (List.length rest - place)
2668                else
2669                    loop ([],a,rest) (place - 1)
2670            end
2671    in
2672    (*----- STRUCTURAL OPERATIONS & QUERIES ------*)
2673        fun new (shape, init) =
2674            if not (Index.validShape shape) then
2675                raise Shape
2676            else
2677                let val length = Index.length shape in
2678                    {shape = shape,
2679                     indexer = Index.indexer shape,
2680                     data = Array.array(length,init)}
2681                end
2682        fun toArray {shape, indexer, data} = data
2683        fun length {shape, indexer, data} =  Array.length data
2684        fun shape {shape, indexer, data} = shape
2685        fun rank t = List.length (shape t)
2686        fun reshape new_shape tensor =
2687            if Index.validShape new_shape then
2688                case (Index.length new_shape) = length tensor of
2689                    true => make'(new_shape, toArray tensor)
2690                  | false => raise Match
2691            else
2692                raise Shape
2693        fun fromArray (s, a) =
2694            case Index.validShape s andalso 
2695                 ((Index.length s) = (Array.length a)) of
2696                 true => make'(s, a)
2697               | false => raise Shape
2698        fun fromList (s, a) = fromArray (s, Array.fromList a)
2699        fun tabulate (shape,f) =
2700            if Index.validShape shape then
2701                let val last = Index.last shape
2702                    val length = Index.length shape
2703                    val c = Array.array(length, f last)
2704                    fun dotable (c, indices, i) =
2705                        (Array.update(c, i, f indices);
2706                         if i <= 1
2707                         then c
2708                         else dotable(c, Index.prev' shape indices, i-1))
2709                in make'(shape,dotable(c, Index.prev' shape last, length-2))
2710                end
2711            else
2712                raise Shape
2713        (*----- ELEMENTWISE OPERATIONS -----*)
2714        fun sub (t, index) = Array.sub(#data t, toInt t index)
2715        fun update (t, index, value) =
2716            Array.update(toArray t, toInt t index, value)
2717        fun map f {shape, indexer, data} =
2718            {shape = shape, indexer = indexer, data = Array.map f data}
2719        fun map2 f t1 t2=
2720            let val {shape=shape1, indexer=indexer1, data=data1} = t1
2721                val {shape=shape2, indexer=indexer2, data=data2} = t2
2722            in
2723                if Index.eq(shape1,shape2) then
2724                    {shape = shape1,
2725                     indexer = indexer1,
2726                     data = Array.map2 f data1 data2}
2727                else
2728                    raise Match
2729        end
2730        fun appi f tensor = Array.appi f (toArray tensor, 0, NONE)
2731        fun app f tensor = Array.app f (toArray tensor)
2732        fun all f tensor =
2733            let val a = toArray tensor
2734            in Loop.all(0, length tensor - 1, fn i =>
2735                        f (Array.sub(a, i)))
2736            end
2737        fun any f tensor =
2738            let val a = toArray tensor
2739            in Loop.any(0, length tensor - 1, fn i =>
2740                        f (Array.sub(a, i)))
2741            end
2742        fun foldl f init tensor = Array.foldl f init (toArray tensor)
2743        fun foldln f init {shape, indexer, data=a} index =
2744            let val (head,lk,tail) = splitList(shape, index)
2745                val li = Index.length head
2746                val lj = Index.length tail
2747                val c = Array.array(li * lj,init)
2748                fun loopi (0, _,  _)  = ()
2749                  | loopi (i, ia, ic) =
2750                    (Array.update(c, ic, f(Array.sub(c,ic), Array.sub(a,ia)));
2751                     loopi (i-1, ia+1, ic+1))
2752                fun loopk (0, ia, _)  = ia
2753                  | loopk (k, ia, ic) = (loopi (li, ia, ic);
2754                                         loopk (k-1, ia+li, ic))
2755                fun loopj (0, _,  _)  = ()
2756                  | loopj (j, ia, ic) = loopj (j-1, loopk(lk,ia,ic), ic+li)
2757            in
2758                loopj (lj, 0, 0);
2759                make'(head @ tail, c)
2760            end
2761        (* --- POLYMORPHIC ELEMENTWISE OPERATIONS --- *)
2762        fun array_map' f a =
2763            let fun apply index = f(Array.sub(a,index)) in
2764                Tensor.Array.tabulate(Array.length a, apply)
2765            end
2766        fun map' f t = Tensor.fromArray(shape t, array_map' f (toArray t))
2767        fun map2' f t1 t2 =
2768            let val d1 = toArray t1
2769                val d2 = toArray t2
2770                fun apply i = f (Array.sub(d1,i), Array.sub(d2,i))
2771                val len = Array.length d1
2772            in
2773                if Index.eq(shape t1, shape t2) then
2774                    Tensor.fromArray(shape t1, Tensor.Array.tabulate(len,apply))
2775                else
2776                    raise Match
2777            end
2778        fun foldl' f init {shape, indexer, data=a} index =
2779            let val (head,lk,tail) = splitList(shape, index)
2780                val li = Index.length head
2781                val lj = Index.length tail
2782                val c = Tensor.Array.array(li * lj,init)
2783                fun loopi (0, _,  _)  = ()
2784                  | loopi (i, ia, ic) =
2785                    (Tensor.Array.update(c,ic,f(Tensor.Array.sub(c,ic),Array.sub(a,ia)));
2786                     loopi (i-1, ia+1, ic+1))
2787                fun loopk (0, ia, _)  = ia
2788                  | loopk (k, ia, ic) = (loopi (li, ia, ic);
2789                                         loopk (k-1, ia+li, ic))
2790                fun loopj (0, _,  _)  = ()
2791                  | loopj (j, ia, ic) = loopj (j-1, loopk(lk,ia,ic), ic+li)
2792            in
2793                loopj (lj, 0, 0);
2794                make'(head @ tail, c)
2795            end
2796    end
2797    end (* MonoTensor *)
2798    open MonoTensor
2799    local
2800        (*
2801         LEFT INDEX CONTRACTION:
2802         a = a(i1,i2,...,in)
2803         b = b(j1,j2,...,jn)
2804         c = c(i2,...,in,j2,...,jn)
2805         = sum(a(k,i2,...,jn)*b(k,j2,...jn)) forall k
2806         MEANINGFUL VARIABLES:
2807         lk = i1 = j1
2808         li = i2*...*in
2809         lj = j2*...*jn
2810         *)
2811        fun do_fold_first a b c lk lj li =
2812            let fun loopk (0, _, _, r, i) = Number.make(r,i)
2813                  | loopk (k, ia, ib, r, i) =
2814                    let val (ar, ai) = Array.sub(a,ia)
2815                        val (br, bi) = Array.sub(b,ib)
2816                        val dr = ar * br - ai * bi
2817                        val di = ar * bi + ai * br
2818                    in loopk (k-1, ia+1, ib+1, r+dr, i+di)
2819                    end
2820                fun loopj (0, ib, ic) = c
2821                  | loopj (j, ib, ic) =
2822                    let fun loopi (0, ia, ic) = ic
2823                          | loopi (i, ia, ic) =
2824                            (Array.update(c, ic, loopk(lk, ia, ib, RNumber.zero, RNumber.zero));
2825                             loopi(i-1, ia+lk, ic+1))
2826                    in loopj(j-1, ib+lk, loopi(li, 0, ic))
2827                    end
2828            in loopj(lj, 0, 0)
2829            end
2830    in
2831        fun +* ta tb =
2832            let val (rank_a,lk::rest_a,a) = (rank ta, shape ta, toArray ta)
2833                val (rank_b,lk2::rest_b,b) = (rank tb, shape tb, toArray tb)
2834            in if not(lk = lk2)
2835               then raise Match
2836               else let val li = Index.length rest_a
2837                        val lj = Index.length rest_b
2838                        val c = Array.array(li*lj,Number.zero)
2839                    in fromArray(rest_a @ rest_b, do_fold_first a b c lk li lj)
2840                    end
2841            end
2842    end
2843    local
2844        (*
2845         LAST INDEX CONTRACTION:
2846         a = a(i1,i2,...,in)
2847         b = b(j1,j2,...,jn)
2848         c = c(i2,...,in,j2,...,jn)
2849         = sum(mult(a(i1,i2,...,k),b(j1,j2,...,k))) forall k
2850         MEANINGFUL VARIABLES:
2851         lk = in = jn
2852         li = i1*...*i(n-1)
2853         lj = j1*...*j(n-1)
2854         *)
2855        fun do_fold_last a b c lk lj li =
2856            let fun loopi(0, _, _, _, _) = ()
2857                  | loopi(i, ia, ic, br, bi) =
2858                    let val (cr,ci) = Array.sub(c,ic)
2859                        val (ar,ai) = Array.sub(a,ia)
2860                        val dr = (ar * br - ai * bi)
2861                        val di = (ar * bi + ai * br)
2862                    in
2863                        Array.update(c,ic,Number.make(cr+dr,ci+di));
2864                        loopi(i-1, ia+1, ic+1, br, bi)
2865                    end
2866                fun loopj(j, ib, ic) =
2867                    let fun loopk(0, _, _) = ()
2868                          | loopk(k, ia, ib) =
2869                            let val (br, bi) = Array.sub(b,ib)
2870                            in
2871                                loopi(li, ia, ic, br, bi);
2872                                loopk(k-1, ia+li, ib+lj)
2873                            end
2874                in case j of
2875                    0 => c
2876                  | _ => (loopk(lk, 0, ib);
2877                          loopj(j-1, ib+1, ic+li))
2878                end (* loopj *)
2879            in
2880                loopj(lj, 0, 0)
2881            end
2882    in
2883        fun *+ ta tb  =
2884            let val (rank_a,shape_a,a) = (rank ta, shape ta, toArray ta)
2885                val (rank_b,shape_b,b) = (rank tb, shape tb, toArray tb)
2886                val (lk::rest_a) = List.rev shape_a
2887                val (lk2::rest_b) = List.rev shape_b
2888            in
2889                if not(lk = lk2) then
2890                    raise Match
2891                else
2892                    let val li = Index.length rest_a
2893                        val lj = Index.length rest_b
2894                        val c = Array.array(li*lj,Number.zero)
2895                    in
2896                        fromArray(List.rev rest_a @ List.rev rest_b,
2897                                  do_fold_last a b c lk li lj)
2898                    end
2899            end
2900    end
2901    (* ALGEBRAIC OPERATIONS *)
2902    infix **
2903    infix ==
2904    infix !=
2905    fun a + b = map2 Number.+ a b
2906    fun a - b = map2 Number.- a b
2907    fun a * b = map2 Number.* a b
2908    fun a ** i = map (fn x => (Number.**(x,i))) a
2909    fun ~ a = map Number.~ a
2910    fun abs a = map Number.abs a
2911    fun signum a = map Number.signum a
2912    fun a == b = map2' Number.== a b
2913    fun a != b = map2' Number.!= a b
2914    fun toString a = raise Domain
2915    fun fromInt a = new([1], Number.fromInt a)
2916    (* TENSOR SPECIFIC OPERATIONS *)
2917    fun *> n = map (fn x => Number.*(n,x))
2918    fun a / b = map2 Number./ a b
2919    fun recip a = map Number.recip a
2920    fun ln a = map Number.ln a
2921    fun pow (a, b) = map (fn x => (Number.pow(x,b))) a
2922    fun exp a = map Number.exp a
2923    fun sqrt a = map Number.sqrt a
2924    fun cos a = map Number.cos a
2925    fun sin a = map Number.sin a
2926    fun tan a = map Number.tan a
2927    fun sinh a = map Number.sinh a
2928    fun cosh a = map Number.cosh a
2929    fun tanh a = map Number.tanh a
2930    fun asin a = map Number.asin a
2931    fun acos a = map Number.acos a
2932    fun atan a = map Number.atan a
2933    fun asinh a = map Number.asinh a
2934    fun acosh a = map Number.acosh a
2935    fun atanh a = map Number.atanh a
2936    fun atan2 (a,b) = map2 Number.atan2 a b
2937    fun normInf a =
2938        let fun accum (y,x) = RNumber.max(x, Number.realPart(Number.abs y))
2939        in  foldl accum RNumber.zero a
2940        end
2941    fun norm1 a =
2942        let fun accum (y,x) = RNumber.+(x, Number.realPart(Number.abs y))
2943        in  foldl accum RNumber.zero a
2944        end
2945    fun norm2 a =
2946        let fun accum (y,x) = RNumber.+(x, Number.abs2 y)
2947        in RNumber.sqrt(foldl accum RNumber.zero a)
2948        end
2949end (* CTensor *)
2950
2951
2952structure RTensorSlice =
2953    struct
2954        structure Tensor = RTensor
2955        structure Index  = Tensor.Index
2956        structure Range  = Range
2957           
2958        type index = Tensor.Index.t
2959        type range = Range.t
2960        type tensor = RTensor.tensor
2961
2962        type slice = {range : range, shape: index, tensor : tensor}
2963
2964        fun slice (rs,tensor) =
2965            let val r = (Range.ranges (Tensor.shape tensor) rs)
2966            in
2967                {range=r,
2968                 shape=(Range.shape r),
2969                 tensor=tensor}
2970            end
2971
2972        fun length ({range, shape, tensor}) = Range.length range
2973        fun base ({range, shape, tensor}) = tensor
2974        fun shape ({range, shape, tensor}) = Range.shape range
2975
2976        fun map f slice = 
2977        let
2978           val te  = #tensor slice
2979           val ra  = #range slice
2980           val fndx  = Range.first ra
2981           val arr = Array.array(length(slice),f (Tensor.sub(te,fndx)))
2982           val i   = ref 0
2983        in 
2984           Range.iteri (fn (ndx) => let val v = f (Tensor.sub (te,ndx)) in (Array.update (arr, !i, v); i := (!i + 1); true) end) ra;
2985           Array.vector arr
2986        end
2987
2988        fun listWrite converter file x =
2989        (List.app (fn x => (TextIO.output(file, "," ^ (converter x)))) x)
2990
2991        fun intListWrite file x = listWrite Int.toString file x
2992
2993        fun app f (slice: slice) = 
2994        let
2995           val te  = #tensor slice
2996           val ra  = #range slice
2997           val fndx  = Range.first ra
2998        in 
2999           Range.iteri (fn (ndx) => (f (Tensor.sub (te,ndx)); true)) ra; ()
3000        end
3001
3002        fun map2 f (sl1: slice) (sl2: slice) = 
3003        let
3004           val _    = if not ((#shape sl1) = (#shape sl2)) then raise Index.Shape else ()
3005           val te1  = #tensor sl1
3006           val te2  = #tensor sl2
3007           val ra1  = #range sl1
3008           val ra2  = #range sl2
3009           val fndx1  = Range.first ra1
3010           val fndx2  = Range.first ra2
3011           val arr   = Array.array(length(sl1),f (Tensor.sub(te1,fndx1),Tensor.sub(te2,fndx2)))
3012           val i     = ref 0
3013        in 
3014           Range.iteri2 (fn (ndx,ndx') => let val v = f (Tensor.sub (te1,ndx),Tensor.sub (te2,ndx')) in (Array.update (arr, !i, v); i := (!i + 1); true) end) (ra1,ra2);
3015           Array.vector arr
3016        end
3017
3018        fun foldl f init (slice: slice) = 
3019        let
3020           val te     = #tensor slice
3021           val ra     = #range slice
3022           val ax     = ref init
3023        in 
3024           Range.iteri (fn (ndx) => let val ax' = f (Tensor.sub (te,ndx),!ax) in ((ax := ax'); true) end) ra;
3025           !ax
3026        end
3027
3028
3029    end                               
3030
3031
3032structure TensorFile =
3033struct
3034
3035type file = TextIO.instream
3036
3037exception Data
3038
3039fun assert NONE = raise Data
3040  | assert (SOME a) = a
3041
3042(* ------------------ INPUT --------------------- *)
3043
3044fun intRead file = assert(TextIO.scanStream INumber.scan file)
3045fun realRead file = assert(TextIO.scanStream RNumber.scan file)
3046fun complexRead file = assert(TextIO.scanStream CNumber.scan file)
3047
3048fun listRead eltScan file =
3049    let val length = intRead file
3050        fun eltRead file = assert(TextIO.scanStream eltScan file)
3051        fun loop (0,accum) = accum
3052          | loop (i,accum) = loop(i-1, eltRead file :: accum)
3053    in
3054        if length < 0
3055        then raise Data
3056        else List.rev(loop(length,[]))
3057    end
3058
3059fun intListRead file = listRead INumber.scan file
3060fun realListRead file = listRead RNumber.scan file
3061fun complexListRead file = listRead CNumber.scan file
3062
3063fun intTensorRead file =
3064    let val shape = intListRead file
3065        val length = Index.length shape
3066        val first = intRead file
3067        val a = ITensor.Array.array(length, first)
3068        fun loop 0 = ITensor.fromArray(shape, a)
3069          | loop j = (ITensor.Array.update(a, length-j, intRead file);
3070                      loop (j-1))
3071    in loop (length - 1)
3072    end
3073
3074fun realTensorRead file =
3075    let val shape = intListRead file
3076        val length = Index.length shape
3077        val first = realRead file
3078        val a = RTensor.Array.array(length, first)
3079        fun loop 0 = RTensor.fromArray(shape, a)
3080          | loop j = (RTensor.Array.update(a, length-j, realRead file);
3081                      loop (j-1))
3082    in loop (length - 1)
3083    end
3084
3085fun complexTensorRead file =
3086    let val shape = intListRead file
3087        val length = Index.length shape
3088        val first = complexRead file
3089        val a = CTensor.Array.array(length, first)
3090        fun loop j = if j = length
3091                     then CTensor.fromArray(shape, a)
3092                     else (CTensor.Array.update(a, j, complexRead file);
3093                           loop (j+1))
3094    in loop 1
3095    end
3096
3097(* ------------------ OUTPUT -------------------- *)
3098fun linedOutput(file, x) = (TextIO.output(file, x); TextIO.output(file, "\n"))
3099
3100fun intWrite file x = linedOutput(file, INumber.toString x)
3101fun realWrite file x = linedOutput(file, RNumber.toString x)
3102fun complexWrite file x =
3103    let val (r,i) = CNumber.split x
3104    in linedOutput(file, concat [RNumber.toString r, " ", RNumber.toString i])
3105    end
3106
3107fun listWrite converter file x =
3108    (intWrite file (length x);
3109     List.app (fn x => (linedOutput(file, converter x))) x)
3110
3111fun intListWrite file x = listWrite INumber.toString file x
3112fun realListWrite file x = listWrite RNumber.toString file x
3113fun complexListWrite file x = listWrite CNumber.toString file x
3114
3115fun intTensorWrite file x = (intListWrite file (ITensor.shape x); ITensor.app (fn x => (intWrite file x)) x)
3116fun realTensorWrite file x = (intListWrite file (RTensor.shape x); RTensor.app (fn x => (realWrite file x)) x)
3117fun complexTensorWrite file x = (intListWrite file (CTensor.shape x); CTensor.app (fn x => (complexWrite file x)) x)
3118
3119fun realTensorSliceWrite file x = (intListWrite file (RTensorSlice.shape x); RTensorSlice.app (fn x => (realWrite file x)) x)
3120
3121end
3122
3123
3124structure RandomTensor =
3125struct
3126
3127fun realRandomTensor (xseed,yseed) shape =
3128    let 
3129        val length = Index.length shape
3130        val seed   = Random.rand (xseed,yseed)
3131        val a      = RTensor.Array.array(length, Random.randReal seed)
3132        fun loop 0 = RTensor.fromArray(shape, a)
3133          | loop j = (RTensor.Array.update(a, length-j, Random.randReal seed);
3134                      loop (j-1))
3135    in loop (length - 1)
3136    end
3137
3138fun intRandomTensor (xseed,yseed) shape =
3139    let 
3140        val length = Index.length shape
3141        val seed   = Random.rand (xseed,yseed)
3142        val a      = ITensor.Array.array(length, Random.randInt seed)
3143        fun loop 0 = ITensor.fromArray(shape, a)
3144          | loop j = (ITensor.Array.update(a, length-j, Random.randInt seed);
3145                      loop (j-1))
3146    in loop (length - 1)
3147    end
3148
3149
3150end
3151
3152
3153val Ne = 2
3154val Ni = 1
3155
3156val S0  = RTensor.fromList ([2,2],[1.0,2.0,3.0,4.0])
3157val _ = (print "S0 = "; TensorFile.realTensorWrite (TextIO.stdOut) S0)
3158val v = RTensor.sub (S0,[0,0])
3159val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3160val v = RTensor.sub (S0,[0,1])
3161val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3162val v = RTensor.sub (S0,[1,0])
3163val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3164val v = RTensor.sub (S0,[1,1])
3165val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3166
3167val SNe = (RTensor.*> 0.5 (RandomTensor.realRandomTensor (13,17) [(Ne+Ni),Ne]) )
3168val SNi = (RTensor.~ (RandomTensor.realRandomTensor (19,23) [(Ne+Ni),Ni]))
3169
3170val _ = TensorFile.realTensorWrite (TextIO.stdOut) SNe
3171
3172val v = RTensor.sub (SNe,[0,0])
3173val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3174val v = RTensor.sub (SNe,[0,1])
3175val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3176val v = RTensor.sub (SNe,[1,0])
3177val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3178val v = RTensor.sub (SNe,[1,1])
3179val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3180val v = RTensor.sub (SNe,[2,0])
3181val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3182val v = RTensor.sub (SNe,[2,1])
3183val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3184
3185val S  = RTensor.cat (1, SNe, SNi)
3186val _ = TensorFile.realTensorWrite (TextIO.stdOut) S
3187
3188val v = RTensor.sub (S,[0,0])
3189val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3190val v = RTensor.sub (S,[0,1])
3191val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3192val v = RTensor.sub (S,[0,2])
3193val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3194val v = RTensor.sub (S,[1,0])
3195val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3196val v = RTensor.sub (S,[1,1])
3197val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3198val v = RTensor.sub (S,[1,2])
3199val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3200val v = RTensor.sub (S,[2,0])
3201val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3202val v = RTensor.sub (S,[2,1])
3203val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3204val v = RTensor.sub (S,[2,2])
3205val _ = (print "v = "; TensorFile.realWrite (TextIO.stdOut) v)
3206
3207val S1 = RTensorSlice.slice ([([0,0],[2,0])],S)
3208val S2 = RTensorSlice.slice ([([0,1],[2,1])],S)
3209val S3 = RTensorSlice.slice ([([0,2],[2,2])],S)
3210
3211val _ = TensorFile.realTensorSliceWrite (TextIO.stdOut) S1
3212val _ = TensorFile.realTensorSliceWrite (TextIO.stdOut) S2
3213val _ = TensorFile.realTensorSliceWrite (TextIO.stdOut) S3
3214
3215
Note: See TracBrowser for help on using the repository browser.