root/trunk/Library/Generator22D.fss @ 3917

(*******************************************************************************
  Generator-of-Generators Library is copyrighted free software by Kento Emoto 
  <emoto[at]ipl.t.u-tokyo.ac.jp> developed under the collaborative 
  research on Fortress Programming Language between Sun Microsystems, 
  Inc. and the University of Tokyo.  
  
  You can redistribute it and/or modify it under the following 
  BSD-style license or the Sun Contributor Agreement that Kento Emoto signed.
  
  
  Copyright (c) 2009 by Kento Emoto
  All rights reserved.
  
  Redistribution and use in source and binary forms, with or without 
  modification, are permitted provided that the following conditions 
  are met:
  
      * Redistributions of source code must retain the above copyright 
        notice, this list of conditions and the following disclaimer.
      * Redistributions in binary form must reproduce the above copyright 
        notice, this list of conditions and the following disclaimer 
        in the documentation and/or other materials provided with the 
        distribution.
      * Neither the name of Kento Emoto nor the names of its 
        contributors may be used to endorse or promote products derived 
        from this software without specific prior written permission.
  
  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
  THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 
  PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 
  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
  PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; 
  OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
  WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR 
  OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF 
  ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  
  
 ******************************************************************************)

component Generator22D

import List.{...}
import Generator2.{...}
export Generator22D

(*---------------------------------- Rows ----------------------------------*)
(* GoG rows generates all rows of the given 2D array. *)


rowsImpl[\E\](g) = 
do 
  bs = g.bounds()
  (il, jl) = bs.lower()
  (iu, ju) = bs.upper()
  array[\Generator[\E\]\](iu-il+1).fill(fn i =>g[(i,jl)#(1,ju-jl+1)])
end

(* Generator2 Rows *)
object Rows[\ E, nat b0, nat s0, nat b1, nat s1 \](g : Array2[\E, b0, s0, b1, s1\]) extends Generator2[\ E \]
  getter seed() : Generator[\ E \] = g
  toString() : String = rowsImpl[\ E \](g).toString()

  (* actual generation of the nested data structure *)
  generate[\ R \](r : Reduction[\ R \], body : Generator[\ E \] -> R) : R = do
        rowsImpl[\ E \](g).generate[\ R \](r, body)
  end
end

rows[\ E, nat b0, nat s0, nat b1, nat s1 \](g : Array2[\E, b0, s0, b1, s1\]) : Rows[\E, b0, s0, b1, s1\] = Rows[\E, b0, s0, b1, s1\](g)


(*---------------------------------- Cols ----------------------------------*)
(* GoG cols generates all columns of the given 2D array. *)

colsImpl[\E\](g) = 
do 
  bs = g.bounds()
  (il, jl) = bs.lower()
  (iu, ju) = bs.upper()
  array[\Generator[\E\]\](ju-jl+1).fill(fn j =>g[(il,j)#(iu-il+1,1)])
end

(* Generator2 Cols *)
object Cols[\ E, nat b0, nat s0, nat b1, nat s1 \](g : Array2[\E, b0, s0, b1, s1\]) extends Generator2[\ E \]
  getter seed() : Generator[\ E \] = g
  toString() : String = colsImpl[\ E \](g).toString()

  (* actual generation of the nested data structure *)
  generate[\ R \](r : Reduction[\ R \], body : Generator[\ E \] -> R) : R = do
        colsImpl[\ E \](g).generate[\ R \](r, body)
  end
end

cols[\ E, nat b0, nat s0, nat b1, nat s1 \](g : Array2[\E, b0, s0, b1, s1\]) : Cols[\E, b0, s0, b1, s1\] = Cols[\E, b0, s0, b1, s1\](g)

(*---------------------------------- Rects ----------------------------------*)

rectsImpl[\E\](g) = 
do 
  bs = g.bounds()
  (il, jl) = bs.lower()
  (iu, ju) = bs.upper()
  s0 = iu - il + 1
  s1 = ju - jl + 1
  ss0 = (s0 (s0+1))/2
  ss1 = (s1 (s1+1))/2
  ret = array[\Generator[\E\]\](ss0, ss1).fill(
fn (i,j) => do

est_n(u) = narrow(|\ SQRT(2 u) /|)
ss(n) = ((n-1) n) / 2
chk(u, n) = do d = u - ss n; d > 0 AND d <= n end
fst(a,b) = a
snd(a,b) = b

str_n(u) = do en = est_n u; cs = <|(chk(u,n),n) | n <- <|en-1, en, en+1|> |>; n = if fst cs[0] then snd cs[0] elif fst cs[1] then snd cs[1] elif fst cs[2] then snd cs[2] end; (n - 1, n - (u - ss n)) end

(i0, i1) : (ZZ32, ZZ32) = str_n(ss0 - i)
(j0, j1) : (ZZ32, ZZ32) = str_n(ss1 - j)

g[(iu-i0,ju-j0)#(i1+1,j1+1)]
end
)  

(* println(BIG //[r <- rows ret] (BIG |||[a <- r] a)) *)
 ret
end

(* Generator2 Rects *) 
(* wanted to be Array2? (GGs on 2D is preferred to be Array2?) *)
object Rects[\ E, nat b0, nat s0, nat b1, nat s1 \](g : Array2[\E, b0, s0, b1, s1\]) extends {Generator2[\ E \]}
  getter seed() : Generator[\ E \] = g
  toString() : String = rectsImpl[\ E \](g).toString()

  (* actual generation of the nested data structure *)
  generate[\ R \](r : Reduction[\ R \], body : Generator[\ E \] -> R) : R = do
        rectsImpl[\ E \](g).generate[\ R \](r, body)
  end

theorems[\R, L1, L2\]() =
<|[\((ActualReduction[\ R, L1 \], ActualReduction[\ R, L2 \], E -> R)->Boolean, (ActualReduction[\ R, L1 \], ActualReduction[\ R, L2 \], E -> R)->R)\]
(fn (r, q, f) => conditionDistributive[\ R, L1, L2 \](r, q, f), fn (r, q, f) => efficientImplDistributive[\ R, L1, L2 \](r, q, f))
|>
  
  (*--- conditions ---*)
  conditionDistributive[\ R, L1, L2 \](q : ActualReduction[\ R, L1 \], r : ActualReduction[\ R, L2 \], f : E -> R) : Boolean = distributes(q, r)

  (*--- wrapper methods for lifting/unlifting. *)
  efficientImplDistributive[\ R, L1, L2 \](q : ActualReduction[\ R, L1 \], r : ActualReduction[\ R, L2 \], f : E -> R) = do
    (* println("efficientImplDistributive") *)

    rp = r.distribute(q)
    i = rp.inner()
    o = rp.outer()
    (* we suppose that the operators i and o are lifted to L1.
       This lifting type should be replaced with the type L of ReductionPair[\R, L\] returned by r.distribute(q).
       But, we do not have any way to capture the type L....
     *)
    o.unlift(efficientImplDistributiveImpl[\L1\](o, i, fn a => i.lift(f(a))))
  end
  (*--- efficient implementations for some conditions ---*)
  efficientImplDistributiveImpl[\ R \](q : Reduction[\ R \], r : Reduction[\ R \], f : E -> R) = do
    (* defintions of operators and functions are found below. *)
    op = MSS2DReductionAbove[\R\](q, r)
    ot = MSS2DReductionBeside[\R\](q, r)


    tup = op.unlift(GeneratorATWrapper(seed()).generateAT[\AnyMaybe\](op, ot, fn a => op.lift(mssBody[\R\](f(a)))))
    tup.b
  end
  
end

rects[\ E, nat b0, nat s0, nat b1, nat s1 \](g : Array2[\E, b0, s0, b1, s1\]) : Rects[\E, b0, s0, b1, s1\] = Rects[\E, b0, s0, b1, s1\](g)



(*-------------------- for abide-tree reductions. -------------------------*)
trait GeneratorAT[\E\]
  generateAT[\R\](op:Reduction[\R\], ot:Reduction[\R\], body: E->R) : R
end

object GeneratorATWrapper[\E, nat b1, nat s1, nat b2, nat s2\](x : Array2[\E, b1, s1, b2, s2\]) extends GeneratorAT[\E\]
  generateAT[\R\](op:Reduction[\R\], ot:Reduction[\R\], body: E->R) : R = do
    loop'(lo1:E,hi1:E,lo2:E,hi2:E) =
      if lo1=hi1 AND lo2=hi2 then body(x[(lo1,lo2)])
      elif lo1=hi1 then
           split2 = partitionL((lo2 BITXOR hi2)+1)                      
           mid2   = hi2 BITAND (BITNOT (split2-1))                      
           ot.join(loop'(lo1,hi1,lo2,mid2-1),loop'(lo1,hi1,mid2,hi2))
      elif lo2=hi2 then
           split1 = partitionL((lo1 BITXOR hi1)+1)                      
           mid1   = hi1 BITAND (BITNOT (split1-1))                      
           op.join(loop'(lo1,mid1-1,lo2,hi2),loop'(mid1,hi1,lo2,hi2))
      else 
           split2 = partitionL((lo2 BITXOR hi2)+1)                      
           mid2   = hi2 BITAND (BITNOT (split2-1))                      
           split1 = partitionL((lo1 BITXOR hi1)+1)                      
           mid1   = hi1 BITAND (BITNOT (split1-1))                      
           op.join(ot.join(loop'(lo1,mid1-1,lo2,mid2-1),loop'(lo1,mid1-1,mid2,hi2)),
                   ot.join(loop'(mid1,hi1,lo2,mid2-1),loop'(mid1,hi1,mid2,hi2)))
      end
    loop'(b1,b1+s1-1,b2,b2+s2-1)
  end
  end

(* for nested abide-tree reductions *)
trait Generator2AT[\E\] extends { Generator2[\ E \], GeneratorAT[\GeneratorAT[\E\]\] }
  generate2AT[\R\](op:Reduction[\R\], ot:Reduction[\R\], om:Reduction[\R\], od:Reduction[\R\], body: E->R) : R
  theoremsAT[\R\]() : List[\((Reduction[\ R \], Reduction[\ R \], Reduction[\ R \], Reduction[\ R \], E -> R)->Boolean, 
                        (Reduction[\ R \], Reduction[\ R \], Reduction[\ R \], Reduction[\ R \], E -> R)->R)\] = <| |>
  naiveImplAT[\ R \](op:Reduction[\R\], ot:Reduction[\R\], om:Reduction[\R\], od:Reduction[\R\], f : E -> R) : R = 
                            generateAT[\R\](op, ot, (fn (x) => x.generateAT[\R\](om,od,f)))
end





(*------------------------ for optimization of rects ---------------------- *)

(* an exception for size-mismatching in array operations. *)
object SizeMismatchException end

(* zipwith operations *)
zipWith(f, a, b, c) = a.ivmap(fn (ij, v) => do (i, j) = ij; f (v, b[(i,j)], c[(i,j)]) end )
zipWith(f, a, b) = a.ivmap(fn (ij, v) => do (i, j) = ij; f (v, b[(i,j)]) end)

(* generalized matrix multiplication: op and ot are used instead of + and * *)
gemm[\T, nat s0, nat s1, nat s2 \](op:Reduction[\T\], ot:Reduction[\T\], sel:Array2[\T,0,s0,0,s1\], other: Array2[\T,0,s1,0,s2\]): Array2[\T, 0, s0, 0, s2\] = do
        res = array2[\T, s0, s2\]()
        mma(a:ZZ32,i:ZZ32,b:ZZ32,j:ZZ32,c:ZZ32,k:ZZ32):() =
            if k>=i AND k>=j then
              if k=1 then
                pr : T = ot.join(sel.get(a,b), other.get(b,c))
                (* If this were atomic, we could parallelize j-partition. *)
                res.put((a,c), op.join(res.get(a,c), pr))
              else
                (k0,k1) = partition(k)
                (mma(a,i,b,j,c,k0),mma(a,i,b,j,c+k0,k1))
              end
            elif j>=i then
                (j0,j1) = partition(j)
                mma(a,i,b,j0,c,k)
                mma(a,i,b+j0,j1,c,k)
            else
                (i0,i1) = partition(i)
                (mma(a,i0,b,j,c,k),mma(a+i0,i1,b,j,c,k))
            end
        mm(a:ZZ32,i:ZZ32,b:ZZ32,j:ZZ32,c:ZZ32,k:ZZ32):() =
            if k>=i AND k>=j then
              if k=1 then
                res.put((a,c), ot.join(sel.get(a,b), other.get(b,c)))
              else
                (k0,k1) = partition(k)
                (mm(a,i,b,j,c,k0),mm(a,i,b,j,c+k0,k1))
              end
            elif j>=i then
                (j0,j1) = partition(j)
                mm(a,i,b,j0,c,k)
                mma(a,i,b+j0,j1,c,k)
            else
                (i0,i1) = partition(i)
                (mm(a,i0,b,j,c,k),mm(a+i0,i1,b,j,c,k))
            end
        if s0=0 OR s1=0 OR s2=0 then
          res
        else
          mm(0,s0,0,s1,0,s2)
          res
        end
      end

(* concatenation of 2D arrays: [A] BSD [B] = [A B] *)
opr BSD[\T, nat s0, nat s1, nat s2\](x : Array2[\T, 0, s0, 0, s1\], y : Array2[\T, 0, s0, 0, s2\]) : Array2[\T, 0, s0, 0, s1 + s2\] = array2[\T, s0, s1 + s2\]().fill(fn (i, j) => if j < s1 then x[(i,j)] else y[(i, j - s1)] end)

(* concatenation of 2D arrays: [A] ABV [B] = [A
                                              B] *)
opr ABV[\T, nat s0, nat s1, nat s2\](x : Array2[\T, 0, s0, 0, s1\], y : Array2[\T, 0, s2, 0, s1\]) : Array2[\T, 0, s0 + s2, 0, s1\] = array2[\T, s0 + s2, s1\]().fill(fn (i, j) => if i < s0 then x[(i,j)] else y[(i-s0, j)] end)

(* construction of block-upper-triangular array: 
     tri(A, B, C, z) = [A B
                        Z C] , where Z is an array filled with z (zero).
*)
tri[\T, nat s0, nat s1, nat s2, nat s3\](x : Array2[\T, 0, s0, 0, s1\], u: Array2[\T, 0, s0, 0, s2\], y : Array2[\T, 0, s2, 0, s3\], z : T) : Array2[\T, 0, s0 + s2, 0, s1 + s3\] = array2[\T, s0 + s2, s1 + s3\]().fill(fn (i, j) => if i < s0 AND j < s1 then x[(i,j)] elif i < s0 AND j >= s1 then u[(i, j-s1)] elif i >= s0 AND j >= s1 then y[(i-s0, j-s1)] else z end)

(* generates with function an array of the same size as the given array. *)
createArray2[\T, nat s0, nat s1\](_:Array2[\T, 0, s0, 0, s1\], f:(ZZ32,ZZ32)->T) : Array2[\T, 0, s0, 0, s1\] = array2[\T, s0, s1\]().fill(f)

(* the lifted data structure for 
   oprimized computation of nested reduction with rects. *)
trait SomeMSS2DTuple[\T\] end
object MSS2DTuple[\T, nat s0, nat s1\](
b : T,
n : Array2[\T, 0, s1, 0, s1\],
s : Array2[\T, 0, s1, 0, s1\],
e : Array2[\T, 0, s0, 0, s0\],
w : Array2[\T, 0, s0, 0, s0\],
ne: Array2[\T, 0, s0, 0, s1\],
nw: Array2[\T, 0, s0, 0, s1\],
se: Array2[\T, 0, s0, 0, s1\],
sw: Array2[\T, 0, s0, 0, s1\],
c : Array2[\T, 0, s1, 0, s1\],
r : Array2[\T, 0, s0, 0, s0\]
) extends SomeMSS2DTuple[\T\]
  getter size() : (ZZ32, ZZ32) = (s0, s1)
  getter asString() : String = "(" b ", " (BIG ||| <|[\String\] "" n, "" s, "" e, "" w, "" ne, "" nw, "" se, "" sw, "" c, "" r|>) ")"
end

mssBody[\T\](v:T) : SomeMSS2DTuple[\T\] = do
  m = array2[\T,1,1\](v)
  MSS2DTuple[\T,1,1\](v, m, m, m, m, m, m, m, m, m, m)
end

__MSS2DTuple[\T, nat s0, nat s1\](
b : T,
n : Array2[\T, 0, s1, 0, s1\],
s : Array2[\T, 0, s1, 0, s1\],
e : Array2[\T, 0, s0, 0, s0\],
w : Array2[\T, 0, s0, 0, s0\],
ne: Array2[\T, 0, s0, 0, s1\],
nw: Array2[\T, 0, s0, 0, s1\],
se: Array2[\T, 0, s0, 0, s1\],
sw: Array2[\T, 0, s0, 0, s1\],
c : Array2[\T, 0, s1, 0, s1\],
r : Array2[\T, 0, s0, 0, s0\]
) : SomeMSS2DTuple[\T\] = MSS2DTuple[\T, s0, s1\](b, n, s, e, w, ne, nw, se, nw, c, r)

object MSS2DReductionAbove[\T\](op, ot) extends { AssociativeReduction[\SomeMSS2DTuple[\T\]\] }
  getter asString(): String = "MSS2DReductionAbove"
  simpleJoin(x: SomeMSS2DTuple[\T\], y: SomeMSS2DTuple[\T\]): SomeMSS2DTuple[\T\] = 
  do
    (s0x, s1x) = x.size()
    (s0y, s1y) = y.size()
    (* make sure s1x = s1y *)
    if (s1x =/= s1y) then
      throw SizeMismatchException
    end
    b = op.join(x.s.indices().generate[\T\](op, fn ij => ot.join(x.s[ij], y.n[ij])), op.join(x.b, y.b))
    n = createArray2(x.n, fn ij => op.join(x.n[ij], ot.join(x.c[ij], y.n[ij])))
    s = createArray2(x.s, fn ij => op.join(y.s[ij], ot.join(y.c[ij], x.s[ij])))
    e = tri(x.e, gemm(op, ot, x.se, y.ne.t()), y.e, op.empty())
    w = tri(x.w, gemm(op, ot, x.sw, y.nw.t()), y.w, op.empty())
    
    ne = (x.ne ABV (createArray2(y.ne, fn ij => do (i, j) = ij; ot.join(y.ne[ij], x.ne[(s0x-1, j)]) end)))
    nw = (x.nw ABV (createArray2(y.nw, fn ij => do (i, j) = ij; ot.join(y.nw[ij], x.nw[(s0x-1, j)]) end)))
    se = ((createArray2(x.se, fn ij => do (i, j) = ij; ot.join(x.se[ij], y.se[(0, j)]) end) ABV y.se))
    sw = ((createArray2(x.sw, fn ij => do (i, j) = ij; ot.join(x.sw[ij], y.sw[(0, j)]) end) ABV y.sw))
    
    c = createArray2(x.c, fn ij => ot.join(x.c[ij], y.c[ij]))
    r = tri(x.r, gemm(op, ot, x.r[(0,s0x-1)#(s0x, 1)], y.r[(0,0)#(1,s0y)]), y.r, op.empty())
    __MSS2DTuple(b, n, s, e, w, ne, nw, se, nw, c, r)
  end
end

object MSS2DReductionBeside[\T\](op, ot) extends { AssociativeReduction[\SomeMSS2DTuple[\T\]\] }
    getter asString(): String = "MSS2DReductionBeside"
    simpleJoin(x: SomeMSS2DTuple[\T\], y: SomeMSS2DTuple[\T\]): SomeMSS2DTuple[\T\] = 
    do
      (s0x, s1x) = x.size()
      (s0y, s1y) = y.size()
      (* make sure s0x = s0y *)
      if (s0x =/= s0y) then
        throw SizeMismatchException
      end
      b = op.join(x.e.indices().generate[\T\](op, fn ij => ot.join(x.e[ij], y.w[ij])), op.join(x.b, y.b))

      n = tri(x.n, gemm(op, ot, x.ne.t(), y.nw), y.n, op.empty())
      s = tri(x.s, gemm(op, ot, x.se.t(), y.sw), y.s, op.empty())
      e = createArray2(x.e, fn ij => op.join(y.e[ij], ot.join(y.r[ij], x.e[ij])))
      w = createArray2(x.w, fn ij => op.join(x.w[ij], ot.join(x.r[ij], y.w[ij])))

      ne = ((createArray2 (x.ne, fn ij => do (i, j) = ij; ot.join(x.ne[ij], y.ne[(i, 0)]) end) BSD y.ne))

      nw = (x.nw BSD (createArray2 (y.nw, fn ij => do (i, j) = ij; ot.join(y.nw[ij], x.nw[(i, s1x-1)]) end)))
      se = ((createArray2(x.se, fn ij => do (i, j) = ij; ot.join(x.se[ij], y.se[(i, 0)]) end) BSD y.se))
      sw = (x.sw BSD (createArray2(y.sw, fn ij => do (i, j) = ij; ot.join(y.sw[ij], x.sw[(i, s1x-1)]) end)))


      c = tri(x.c, gemm(op, ot, (x.c[(0,s1x-1)#(s1x,1)]), (y.c[(0,0)#(1,s1y)])), y.c, op.empty())
      r = createArray2(x.r,fn  ij => ot.join(x.r[ij], y.r[ij]))

      __MSS2DTuple(b, n, s, e, w, ne, nw, se, nw, c, r)      
    end
end



end