Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

[diderot] View of /branches/ein16/src/compiler/ein/mkoperators.sml
 [diderot] / branches / ein16 / src / compiler / ein / mkoperators.sml

# View of /branches/ein16/src/compiler/ein/mkoperators.sml

Thu Aug 11 14:00:48 2016 UTC (2 years, 11 months ago) by cchiw
File size: 32542 byte(s)
results update
(* creates EIN operators
*
* This code is part of the Diderot Project (http://diderot-language.cs.uchicago.edu)
*
* COPYRIGHT (c) 2015 The University of Chicago
*)

structure MkOperators = struct

local

structure E = Ein
structure P=Printer
in

fun specialize(alpha,inc)=  List.tabulate(length(alpha), (fn(x)=>E.V (x+inc)))
fun sumIds(n,inc,alpha)=let
val vs=List.tabulate(n, (fn v=>E.V (v+inc)))
in ListPair.map  (fn(v,i)=>(v, 0, i-1))   (vs, alpha) end
fun sumIds2(n,i)=List.tabulate(n, (fn v=>(E.V v, 0, i)))
(*   val subst_flag() = 1(*here*)*)

val substFlag= ref true
val controls = [("substFlag",substFlag,"substFlag")]
fun subst_flag()=if (!substFlag) then 1 else 0

fun mkCx es =List.map (fn c => E.C (c,true)) es
fun mkCxSingle c = E.C (c,true)
fun mkSxSingle c = E.C (c,false)
params = [E.TEN(subst_flag(),[]), E.TEN(subst_flag(),[])] ,
index = [],
body = E.Opn(E.Add,[ E.Tensor(0, []), E.Tensor(1, [])])
}

(* Adding tensors : < X{\alpha} + Y_{\alpha}>_{\alpha} *)
val expindex= specialize(alpha,0)
in
E.EIN{
params = [E.TEN(subst_flag(),alpha), E.TEN(subst_flag(),alpha)],
index = alpha,
body = E.Opn(E.Add,[E.Tensor(0, expindex), E.Tensor(1, expindex)])
}
end

(* mkField functions*)
(*Adding Fields : < F{\alpha} + G_{\alpha}>_{\alpha} *)
val expindex= specialize(shape,0)
in E.EIN{
params = [E.FLD(dim),E.FLD(dim)],
index = shape,
}
end

(*Tensor and Fields*)
val expindex= specialize(shape,0)
in E.EIN{
params = [E.TEN(subst_flag(),shape),E.FLD(dim)],
index = shape,
}
end
(********************************* Sub **************************************)
(* Subtract Scalars*)
val subRR = E.EIN{
params = [E.TEN(subst_flag(),[]), E.TEN(subst_flag(),[])],
index = [],
body = E.Op2(E.Sub, E.Tensor(0, []), E.Tensor(1, []))
}

fun subTT alpha=let
val expindex= specialize(alpha,0)
in
E.EIN{
params = [E.TEN(subst_flag(),alpha), E.TEN(subst_flag(),alpha)],
index = alpha,
body = E.Op2(E.Sub,E.Tensor(0, expindex), E.Tensor(1, expindex))
}
end

fun subTF(dim,shape) =let
val expindex= specialize(shape,0)
in
E.EIN{
params = [E.TEN(subst_flag(),shape),E.FLD(dim)],
index = shape,
}
end

fun subFT(dim,shape) = let
val expindex= specialize(shape,0)
in
E.EIN{
params = [E.TEN(subst_flag(),shape),E.FLD(dim)],
index = shape,
body = E.Op2(E.Sub,E.Field(1, expindex),E.Lift(E.Tensor(0, expindex)))
}
end

fun subFF(dim,shape) =let
val expindex= specialize(shape,0)
in E.EIN{
params = [E.FLD(dim),E.FLD(dim)],
index = shape,
body = E.Op2(E.Sub,E.Field(0, expindex),E.Field(1, expindex))
}
end
(********************************** mul *************************************)
(* Product Scalars*)
val mulRR = E.EIN{
params =[E.TEN(subst_flag(),[]), E.TEN(subst_flag(),[])],
index = [],
body = E.Opn(E.Prod,[ E.Tensor(0, []), E.Tensor(1, [])])
}
(* scalar times tensor product: <s * T_{\alpha}>_{\alpha} *)
fun mulRT alpha = let
val expindex= specialize(alpha,0)
in
E.EIN{
params = [E.TEN(subst_flag(),[]), E.TEN(subst_flag(),alpha)],
index = alpha,
body = E.Opn(E.Prod,[E.Tensor(0, []),  E.Tensor(1, expindex)])
}
end
fun mulRF(dim,shape) =let
val  expindex= specialize(shape,0)
in E.EIN{
params = [E.TEN(subst_flag(),[]),E.FLD(dim)],
index = shape,
body = E.Opn(E.Prod,[E.Lift( E.Tensor(0,[])), E.Field(1,expindex)])
}
end
fun mulST(dim,shape) =let
val  expindex= specialize(shape,0)
in E.EIN{
params = [E.TEN(subst_flag(),shape),E.FLD(dim)],
index = shape,
body = E.Opn(E.Prod,[E.Lift( E.Tensor(0,expindex)), E.Field(1,[])])
}
end
fun mulSS dim = E.EIN{
params = [E.FLD(dim),E.FLD(dim)],
index = [],
body = E.Opn(E.Prod,[E.Field(0, []),E.Field(1, [])])
}

fun mulSF(dim,shape) =let
val  expindex= specialize(shape,0)
in E.EIN{
params = [E.FLD(dim),E.FLD(dim)],
index = shape,
body = E.Opn(E.Prod,[E.Field(0, []),E.Field(1, expindex)])
}
end

(************************************div ************************************)
(* Divide Scalars*)
val divRR = E.EIN{
params = [E.TEN(subst_flag(),[]), E.TEN(subst_flag(),[])],
index = [],
body = E.Op2(E.Div, E.Tensor(0, []), E.Tensor(1, []))
}

fun divTR alpha =let
val expindex= specialize(alpha,0)
in
E.EIN{
params = [E.TEN(subst_flag(),alpha), E.TEN(subst_flag(),[])],
index = alpha,
body = E.Op2(E.Div,E.Tensor(0, expindex), E.Tensor(1,[]))
}
end

fun divFR(dim,shape) = let
val expindex= specialize(shape,0)
in E.EIN{
params = [E.FLD(dim),E.TEN(subst_flag(),[])],
index = shape,
body = E.Op2(E.Div,E.Field(0, expindex),E.Lift(  E.Tensor(1, [])))
}
end

fun divTS(dim, shape) = let
val expindex = specialize(shape,0)
in E.EIN{
params = [E.TEN(subst_flag(), shape), E.FLD(dim)],
index = shape,
body = E.Op2(E.Div, E.Lift(E.Tensor(0, expindex)), E.Field(1, []))
}
end

fun divSS dim = E.EIN{
params = [E.FLD(dim),E.FLD(dim)],
index = [],
body = E.Op2(E.Div,E.Field(0, []),E.Field(1, []))
}

fun divFS(dim,shape) = let
val  expindex= specialize(shape,0)
in  E.EIN{
params = [E.FLD(dim),E.FLD(dim)],
index = shape,
body = E.Opn(E.Prod,[E.Field(0, expindex),E.Op2(E.Div,E.B(E.Const 1),E.Field(1, []))])
}
end
(************************************* Neg **********************************)
fun negTT alpha=let
val  expindex= specialize(alpha,0)
(*changed tensor lift variable here *)
in
E.EIN{
params = [E.TEN(subst_flag(),alpha)],index = alpha,
body = E.Op1(E.Neg,(E.Tensor(0, expindex)))
}
end

fun negFF(dim,shape) = let
val expindex = specialize(shape,0)
in E.EIN{
params = [E.FLD(dim)],index = shape,
body = E.Op1(E.Neg,E.Field(0, expindex))
}
end

(****************************** cross ***********************************)
(*crossProduct is on 3D vectors ..vec3 t8=t0 × t1; *)
val cross3TT = E.EIN{
params = [E.TEN(subst_flag(),[3]), E.TEN(subst_flag(),[3])],
index= [3],
body=E.Sum([(E. V 1,0,2),(E.V 2,0,2)],
E.Opn(E.Prod,[ E.G(E.Epsilon(E.V 0, E.V 1, E.V 2)),
E.Tensor(0, [E.V 1]),  E.Tensor(1, [E.V 2]) ]))
}

(*2-d cross product Eps_{ij}U_i V_j*)
val cross2TT = E.EIN{
params = [E.TEN(subst_flag(),[2]), E.TEN(subst_flag(),[2])],
index= [],
body=E.Sum([(E.V 0,0,1),(E.V 1,0,1)],
E.Opn(E.Prod,[ E.G(E.Eps2(E.V 0, E.V 1)), E.Tensor(0,[E.V 0]),  E.Tensor(1, [E.V 1])]))
}

(*Field Cross Product*)
val cross2FF = E.EIN{
params = [E.FLD(2), E.FLD(2)],index= [],
body=E.Sum([(E.V 0,0,1),(E.V 1,0,1)],
E.Opn(E.Prod,[E.G(E.Eps2(E.V 0, E.V 1)), E.Field(0, [E.V 0]),  E.Field(1, [E.V 1]) ]))
}

(*Field Cross Product*)
val cross3FF = E.EIN{
params = [E.FLD(3), E.FLD(3)],index= [3],
body=E.Sum([(E. V 1,0,2),(E.V 2,0,2)],
E.Opn(E.Prod,[ E.G(E.Epsilon(E.V 0, E.V 1, E.V 2)),
E.Field(0, [E.V 1]),  E.Field(1, [E.V 2]) ]))
}

val cross2FT = E.EIN{
params = [E.FLD(2), E.TEN(subst_flag(),[2])],index= [],
body=E.Sum([(E.V 0,0,1),(E.V 1,0,1)],
E.Opn(E.Prod,[E.G(E.Eps2(E.V 0, E.V 1)), E.Field(0, [E.V 0]),  E.Lift(E.Tensor(1, [E.V 1])) ]))
}

val cross3FT = E.EIN{
params = [E.FLD(3), E.TEN(subst_flag(),[3])],index= [3],
body=E.Sum([(E.V 1,0,2),(E.V 2,0,2)],
E.Opn(E.Prod,[ E.G(E.Epsilon(E.V 0, E.V 1, E.V 2)),
E.Field(0, [E.V 1]), E.Lift( E.Tensor(1, [E.V 2])) ]))
}

val cross2TF = E.EIN{
params = [E.TEN(subst_flag(),[2]), E.FLD(2)],index= [],
body=E.Sum([(E.V 0,0,1),(E.V 1,0,1)],
E.Opn(E.Prod,[E.G(E.Eps2(E.V 0, E.V 1)), E.Lift( E.Tensor(0, [E.V 0])),  E.Field(1, [E.V 1]) ]))
}

val cross3TF = E.EIN{
params = [E.TEN(subst_flag(),[3]), E.FLD(3)],index= [3],
body=E.Sum([(E.V 1,0,2),(E.V 2,0,2)],
E.Opn(E.Prod,[E.G(E.Epsilon(E.V 0, E.V 1, E.V 2)),
E.Lift(E.Tensor(0, [E.V 1])),E.Field(1, [E.V 2])]))
}

(********************  outerProduct ********************************)
(*  newbie*)
fun outerTT(alpha,beta) =let
val expindexA= specialize(alpha,0)
val expindexB= specialize(beta,length(alpha))
in
E.EIN{
params = [E.TEN(subst_flag(),alpha), E.TEN(subst_flag(),beta)],
index= alpha@beta,
body= E.Opn(E.Prod,[E.Tensor(0,  expindexA), E.Tensor(1,  expindexB)])
}
end

(*Assumes same dimension vector field *)
fun outerFF(dim,alpha,beta) =let
val expindexA= specialize(alpha,0)
val expindexB= specialize(beta,length(alpha))
in E.EIN{
params = [E.FLD(dim),E.FLD(dim)],
index = alpha@beta,
body = E.Opn(E.Prod,[E.Field(0, expindexA),E.Field(1, expindexB)])
}
end

fun outerTF(dim,alpha,beta) =let
val expindexA= specialize(alpha,0)
val expindexB= specialize(beta,length(alpha))
in E.EIN{
params = [E.TEN(subst_flag(),alpha),E.FLD(dim)],
index = alpha@beta,
body = E.Opn(E.Prod,[E.Lift(E.Tensor(0, expindexA)),E.Field(1, expindexB)])
}
end

fun outerFT(dim,alpha,beta) =let
val _ = print (String.concat["length(alpha):", Int.toString(length(alpha))])
val _ = print (String.concat["length(beta):", Int.toString(length(beta))])
val expindexA= specialize(alpha,0)
val expindexB= specialize(beta,length(alpha))
in E.EIN{
params = [E.FLD(dim),E.TEN(subst_flag(),beta)],
index = alpha@beta,
body = E.Opn(E.Prod,[E.Field(0, expindexA),E.Lift(E.Tensor(1, expindexB))])
}
end

(*************************** inner product **********************************)
(* generic inner product: <T_{\alpha i} * T_{i \beta}>_{\alpha \beta} *)
fun innerTT(shape1,i::beta) = let
val alpha= List.take(shape1,length(shape1)-1)
val expindexA= specialize(alpha,0)
val expindexB= specialize(beta,(length(alpha)))
val s'=E.V(length(alpha)+ length(beta))
val s''=[(s',0,i-1)]
in
E.EIN{
params = [E.TEN(subst_flag(),shape1) ,E.TEN(subst_flag(),i::beta)],
index = alpha@beta,
body = E.Sum(s'', E.Opn(E.Prod,[
E.Tensor(0, expindexA@[s']),   (* T_{\alpha i} *)
E.Tensor(1, [s']@expindexB )  (* T'_{i \beta} *)]))
}
end
| innerTT _ = raise Fail "Wrong shape for inner product"

(* generic inner product: <T_{\alpha i} * T_{i \beta}>_{\alpha \beta} *)
fun innerFF(shape1,dim,i::beta) = let
val alpha= List.take(shape1,length(shape1)-1)
val expindexA= specialize(alpha,0)
val expindexB= specialize(beta,(length(alpha)))
val sid=E.V(length(alpha)+ length(beta))
in E.EIN{
params = [E.FLD(dim) ,E.FLD(dim)],index = alpha@beta,
body = E.Sum([(sid,0,i-1)],
E.Opn(E.Prod,[
E.Field(0, expindexA@[sid]),   (* F_{\alpha i} *)
E.Field(1, [sid]@expindexB )  (* F'_{i \beta} *)
]))}
end
| innerFF _ = raise Fail "Wrong shape for innerProductField"

fun innerFT(shape1,dim,i::beta) = let
val alpha= List.take(shape1,length(shape1)-1)
val expindexA= specialize(alpha,0)
val expindexB= specialize(beta,(length(alpha)))
val sid=E.V(length(alpha)+ length(beta))
in E.EIN{
params = [E.FLD(dim),E.TEN(subst_flag(),i::beta)],
index = alpha@beta,
body = E.Sum([(sid,0,i-1)],
E.Opn(E.Prod,[
E.Field(0, expindexA@[sid]),   (* F_{\alpha i} *)
E.Lift(E.Tensor(1, [sid]@expindexB ))  (* F'_{i \beta} *)
]))}
end
| innerFT _ = raise Fail "Wrong shape for innerProductFieldTensor"

fun innerTF(shape1,dim,i::beta) = let
val alpha= List.take(shape1,length(shape1)-1)
val expindexA= specialize(alpha,0)
val expindexB= specialize(beta,(length(alpha)))
val sid=E.V(length(alpha)+ length(beta))
in E.EIN{
params = [E.TEN(subst_flag(),shape1),E.FLD(dim)],index = alpha@beta,
body = E.Sum([(sid,0,i-1)],
E.Opn(E.Prod,[E.Lift(E.Tensor(0, expindexA@[sid])),   (* F_{\alpha i} *)
E.Field(1, [sid]@expindexB )  (* F'_{i \beta} *)
]))}
end
| innerTF _ = raise Fail "Wrong shape for innerProductTensorField"

(*************************** colon product **********************************)
(*<T_{\alpha i j} * B{i j \beta }>_\alpha \beta*)
fun colonTT(shape1,i::j::beta) = let
val alpha= List.take(shape1,length(shape1)-2)
val expindexA= specialize(alpha,0)
val expindexB= specialize(beta,(length(alpha)))
val sumi=length(alpha)+ length(beta)
val s'=[E.V sumi,E.V(sumi+1)]
val sx=[(E.V sumi,0,i-1),(E.V(sumi+1),0,j-1)]
in
E.EIN{
params = [E.TEN(subst_flag(),shape1),E.TEN(subst_flag(),i::j::beta)],
index = alpha@beta,
body = E.Sum(sx,E.Opn(E.Prod,[ E.Tensor(0, expindexA@s'), E.Tensor(1,s'@expindexB)]))
}
end

(*<F_{\alpha i j} * G_{i j \beta }>_\alpha \beta*)
fun colonFF(dim,shape1,i::j::beta) = let
val alpha= List.take(shape1,length(shape1)-2)
val expindexA= specialize(alpha,0)
val expindexB= specialize(beta,(length(alpha)))
val sumi=length(alpha)+ length(beta)
val s'=[E.V sumi,E.V(sumi+1)]
val sx=[(E.V sumi,0,i-1),(E.V(sumi+1),0,j-1)]
in
E.EIN{
params = [E.FLD dim,E.FLD dim],
index = alpha@beta,
body = E.Sum(sx,E.Opn(E.Prod,[E.Field(0, expindexA@s'), E.Field(1,s'@expindexB)]))
}
end

(*<F_{\alpha i j} * T_{i j \beta }>_\alpha \beta*)
fun colonFT(dim,shape1,i::j::beta) = let
val alpha= List.take(shape1,length(shape1)-2)
val expindexA= specialize(alpha,0)
val expindexB= specialize(beta,(length(alpha)))
val sumi=length(alpha)+ length(beta)
val s'=[E.V sumi,E.V(sumi+1)]
val sx=[(E.V sumi,0,i-1),(E.V(sumi+1),0,j-1)]
in
E.EIN{
params = [E.FLD dim,E.TEN(subst_flag(),shape1)],
index = alpha@beta,
body = E.Sum(sx,E.Opn(E.Prod,[E.Field(0, expindexA@s'), E.Lift(E.Tensor(1,s'@expindexB))]))
}
end

(*<T_{\alpha i j} * G{i j \beta }>_\alpha \beta*)
fun colonTF(dim,shape1,i::j::beta) = let
val alpha= List.take(shape1,length(shape1)-2)
val expindexA= specialize(alpha,0)
val expindexB= specialize(beta,(length(alpha)))
val sumi=length(alpha)+ length(beta)
val s'=[E.V sumi,E.V(sumi+1)]
val sx=[(E.V sumi,0,i-1),(E.V(sumi+1),0,j-1)]
in
E.EIN{
params = [E.TEN(subst_flag(),i::j::beta),E.FLD dim],
index = alpha@beta,
body = E.Sum(sx,E.Opn(E.Prod,[E.Lift(E.Tensor(0, expindexA@s')), E.Field(1,s'@expindexB)]))
}
end

(******************** Norm  ********************************)

(*n.t.s tensor norms/normalize use high-il ops
*we can norm any size and use einapp
* normalize only implemented for vectors and use mid il-op
*)
(*get norm, but without the sqrt
* implemented as a summation over a modulate
*)
fun magnitudeTT alpha =let
val expindex= specialize(alpha,0)
val sx= sumIds(length(alpha),0,alpha)
val e= E.EIN{
params = [E.TEN(subst_flag(),alpha), E.TEN(subst_flag(),alpha)],
index = [],
body = E.Sum(sx,
E.Opn(E.Prod,[E.Tensor(0, expindex), E.Tensor(1, expindex)]))
}
(*val _=print(String.concat["\ntesting mag tt",P.printerE e ])*)
in e
end

fun magnitudeFF(dim,[]) =E.EIN{params = [E.FLD(dim)],index = [],body =
E.Op1(E.Sqrt, E.Opn(E.Prod,[E.Field(0, []),E.Field(0,[])]))}
| magnitudeFF(dim,alpha) =let
val expindex= specialize(alpha,0)
val sx= sumIds(length(alpha),0,alpha)
val e=
E.EIN{
params = [E.FLD(dim)],
index = [],
body = E.Op1(E.Sqrt,E.Sum(sx, E.Opn(E.Prod,[E.Field(0, expindex),E.Field(0,expindex)])))
}
(*val _=print(String.concat["\ntesting mag FF",P.printerE e ])*)
in e
end

(*
fun magnitudeFFPow(dim,[]) =E.EIN{params = [E.FLD(dim),E.FLD(dim)],index = [],body =  E.Field(0, [])}
| magnitudeFF(dim,alpha) =let
val expindex= specialize(alpha,0)
val sx= sumIds(length(alpha),0,alpha)
val e=
E.EIN{
params = [E.FLD(dim),E.FLD(dim)],
index = [],
body =  E.PowEmb(E.Field(0, expindex),sx,2)
}
val _=print(String.concat["\ntesting mag FF",P.printerE e ])
in e
end
*)

fun normalizeTT [] =E.EIN{params = [E.TEN(subst_flag(),[])],index = [],body =  E.Tensor(0, [])}
| normalizeTT alpha =let
val expindex= specialize(alpha,0)
(*shift indices in the inner product*)
val expindexDot= specialize(alpha,length(alpha))
val sx= sumIds(length(alpha),length(alpha),alpha)
val f=E.Tensor(0, expindex)
val g=E.Tensor(1, expindexDot)
val e=
E.EIN{
params = [E.TEN(subst_flag() ,alpha) ,E.TEN(subst_flag(),alpha)],
index = alpha,
body = E.Opn(E.Prod,[f,E.Op2(E.Div,E.B(E.Const 1),E.Op1(E.Sqrt,E.Sum(sx, E.Opn(E.Prod,[g,g]))))])
}
in e
end

fun normalizeFF(dim,[]) =E.EIN{params = [E.FLD(dim)],index = [],body =  E.Field(0, [])}
| normalizeFF(dim,alpha) =let
val expindex= specialize(alpha,0)
(*shift indices in the inner product*)
val expindexDot= specialize(alpha,length(alpha))
val i=List.hd alpha
val sx= sumIds(length(alpha),length(alpha),alpha)
val f=E.Field(0, expindex)
val g=E.Field(1, expindexDot)
val e=
E.EIN{
params = [E.FLD(dim) ,E.FLD(dim)],
(* index = [dim],*)  index = alpha,
body = E.Opn(E.Prod,[f,E.Op2(E.Div,E.B(E.Const 1),E.Op1(E.Sqrt,E.Sum(sx, E.Opn(E.Prod,[g,g]))))])
}
(*val _=print(String.concat["\ntesting normalize FF",P.printerE e ])*)
in e
end

(************************* trace *************************)
(* Trace: <M_{i, i}>  This one Sx represents both i's*)
fun traceT dim = E.EIN{
params = [E.TEN(1,[dim,dim])],index = [],
body = E.Sum([(E.V 0,0,dim-1)],E.Tensor(0, [E.V 0, E.V 0]))
}

(* Trace: <Sigma_i F_{\alpha i, i}>  This one Sx represents both i's*)
fun traceF(dim, d2, alpha) =let
val expindex = specialize(alpha,0)
val s = E.V(length(alpha))
in
E.EIN{
params = [E.FLD(dim)],
index = alpha,
body = E.Sum([(s, 0, d2-1)], E.Field(0, expindex@[s,s]))
}
end

(************************* tranpose *************************)
fun transposeT alpha =E.EIN{
params = [E.TEN(subst_flag(),alpha)],
index= List.rev alpha,
body= E.Tensor(0, [E.V 1,E.V 0])
}

(*Transpose Field F_{ji}*)
fun transposeF(dim,i,j) =E.EIN{
params = [E.FLD(dim)], index= [i,j],
body= E.Field(0, [E.V 1,E.V 0])
}
(************************* determinant *************************)
val detT2 =E.EIN{
params = [E.TEN(0,[2,2])],
index= [],
body= E.Op2(E.Sub,
E.Opn(E.Prod,[E.Tensor(0, mkCx [0,0]),E.Tensor(0, mkCx [1,1])]),
E.Opn(E.Prod,[E.Tensor(0, mkCx [0,1]),E.Tensor(0, mkCx [1,0])])
)
}
(*
val detT3  = let
val a=E.Tensor(0,  mkCx [0,0])
val b=E.Tensor(0,  mkCx [0,1])
val c=E.Tensor(0,  mkCx [0,2])
val d=E.Tensor(0,  mkCx [1,0])
val e=E.Tensor(0,  mkCx [1,1])
val f=E.Tensor(0,  mkCx [1,2])
val g=E.Tensor(0,  mkCx [2,0])
val h=E.Tensor(0,  mkCx [2,1])
val i=E.Tensor(0,  mkCx [2,2])
in E.EIN{
params = [E.TEN(0,[3,3])],
index= [],
body=E.Op2(E.Sub,
}
end
*)
val detT3=
E.EIN{ params = [E.TEN(0,[3,3])], index= [],
body=
E.Sum([(E.V 0,0,2)],
E.Opn(E.Prod,[E.Tensor(0,[mkCxSingle 0, E. V 0]),
E.Sum([(E.V 1,0,2)],
E.Opn(E.Prod,[E.Tensor(0,[mkCxSingle 1, E. V 1]),
E.Sum([(E.V 2,0,2)], E.Opn(E.Prod,[E.G(E.Epsilon(E.V 0, E.V 1, E.V 2)),E.Tensor(0,[mkCxSingle 2, E. V 2])]))]))])
)
}

fun detF2(dim) = E.EIN{
params = [E.FLD dim], index= [],
body= E.Op2(E.Sub,
E.Opn(E.Prod,[E.Field(0, mkCx [0,0]),E.Field(0, mkCx [1,1])]),
E.Opn(E.Prod,[E.Field(0, mkCx [0,1]),E.Field(0, mkCx [1,0])]))
}

val detF3_full = let
val a=E.Field(0, mkCx [0,0])
val b=E.Field(0, mkCx [0,1])
val c=E.Field(0, mkCx [0,2])
val d=E.Field(0, mkCx [1,0])
val e=E.Field(0, mkCx [1,1])
val f=E.Field(0, mkCx [1,2])
val g=E.Field(0, mkCx [2,0])
val h=E.Field(0, mkCx [2,1])
val i=E.Field(0, mkCx [2,2])
in E.EIN{
params = [E.FLD(3)], index= [],
body=E.Op2(E.Sub,
}
end

val detF3A=
E.EIN{ params = [E.FLD(3)], index= [],
body=   E.Sum([(E.V 0,0,2),(E.V 1,0,2),(E.V 2,0,2)],
E.Opn(E.Prod,[ E.G(E.Epsilon(E.V 0, E.V 1, E.V 2)),
E.Field(0,[mkCxSingle 0, E. V 0]),
E.Field(0,[mkCxSingle 1, E. V 1]),
E.Field(0,[mkCxSingle 2, E. V 2])]))
}

fun detF3 (dim) =
E.EIN{ params = [E.FLD(dim)], index= [],
body=
E.Sum([(E.V 0,0,2)],
E.Opn(E.Prod,[E.Field(0,[mkCxSingle 0, E. V 0]),
E.Sum([(E.V 1,0,2)],
E.Opn(E.Prod,[E.Field(0,[mkCxSingle 1, E. V 1]),
E.Sum([(E.V 2,0,2)], E.Opn(E.Prod,[E.G(E.Epsilon(E.V 0, E.V 1, E.V 2)),E.Field(0,[mkCxSingle 2, E. V 2])]))]))])
)
}

fun mkField(ix, jx)= E.Field(0, [E.C ix, E.C jx])
fun invF2(dim) = let
val f00 = mkField(0,0)
val f11 = mkField(1,1) - mkField(0,1)*mkField(1,0);
return [[m[1,1],-m[0,1]],[-m[1,0],m[0,0]]]/d;

(************************* Pow **************************)
fun powF (dim,n) = E.EIN{params = [E.FLD(dim)],index = [], body = E.Op1(E.PowInt n, E.Field(0,[]))}
fun expF dim = E.EIN{params = [E.FLD(dim)],index = [], body = E.Op1(E.Exp,E.Field(0,[]))}
val expT = E.EIN{params = [E.TEN(0,[])],index = [], body = E.Op1(E.Exp,E.Tensor(0,[]))}
(************************* Lifted single-argument math functions *************************)
local
fun liftFn f dim = E.EIN{
params = [E.FLD dim],
index = [],
body = E.Op1(f, E.Field(0, []))
}
fun op1Fn op0= E.EIN{
params = [E.TEN(0,[])],
index = [],
body = E.Op1(op0, E.Field(0, []))
}
in
fun powF (dim, n) = E.EIN{
params = [E.FLD dim],
index = [], body = E.Op1(E.PowInt n, E.Field(0, []))
}
val sqrtF = liftFn E.Sqrt
val cosF = liftFn E.Cosine
val acosF = liftFn E.ArcCosine
val sinF = liftFn E.Sine
val asinF = liftFn E.ArcSine
val tanF = liftFn E.Tangent
val atanF = liftFn E.ArcTangent

val sqrtR =op1Fn E.Sqrt
val cosR =op1Fn E.Cosine
val acosR = op1Fn E.ArcCosine
val sinR = op1Fn E.Sine
val asinR = op1Fn E.ArcSine
val tanR = op1Fn E.Tangent
val atanR = op1Fn E.ArcTangent
end (* local *)

(*************************  modulate *************************)
fun modulateTT shape = let
val expindex = specialize(shape, 0)
in
E.EIN{
params = [E.TEN(subst_flag(), shape), E.TEN(subst_flag(), shape)],
index = shape,
body = E.Opn(E.Prod,[E.Tensor(0, expindex), E.Tensor(1, expindex)])
}
end

fun modulateFF(shape, dim)  = let
val expindex = specialize(shape, 0)
in
E.EIN{
params = [E.FLD(dim), E.FLD(dim)],
index = shape,
body = E.Opn(E.Prod,[E.Field(0, expindex), E.Field(1, expindex)])
}
end
fun modulateTF(shape, dim)  = let
val expindex = specialize(shape, 0)
in
E.EIN{
params = [E.TEN(subst_flag(), shape), E.FLD(dim)],
index = shape,
body = E.Opn(E.Prod,[E.Lift(E.Tensor(0, expindex)), E.Field(1, expindex)])
}
end
fun modulateFT(shape, dim)  = let
val expindex = specialize(shape, 0)
in
E.EIN{
params = [E.FLD(dim), E.TEN(subst_flag(), shape)],
index = shape,
body = E.Opn(E.Prod,[E.Field(0, expindex), E.Lift(E.Tensor(1, expindex))])
}
end

(************************* other tensor ops *************************)
fun identity dim =E.EIN{
params = [],index = [dim,dim], body = E.G(E.Delta(E.V(0), E.V(1)))
}

fun sliceT (mask, const, rstTy, argTy) = let
fun fg e = String.concat(List.map (fn e => Int.toString(e)^",") e)
val _ = print(String.concat["mask:",fg const, "argTy:", fg argTy,"rstTy:", fg rstTy])

fun iter ([], _, cnt) = []
| iter(true::es, c::cs, cnt)=(mkSxSingle c)::iter(es, cs, cnt)
|iter(false::es, cs, cnt)= (E.V cnt)::iter(es, cs, cnt+1)
val ix = iter(mask, const, 0)
val ein0 = E.EIN{params = [E.TEN(subst_flag(), argTy)], index = rstTy, body = E.Tensor(0, ix)}

in ein0 end

fun sliceF (mask, const, rstTy, dim) = let
fun iter ([], _, cnt) = []
| iter(true::es, c::cs, cnt)=(mkSxSingle c)::iter(es, cs, cnt)
| iter(false::es, cs, cnt)=(E.V cnt)::iter(es, cs, cnt+1)
val ix = iter(mask, const,0 )
val ein0 =  E.EIN{params = [E.FLD(dim)], index = rstTy, body = E.Field(0, ix)}

in ein0 end

(******************** other field ops  ********************************)
(*FLD here is bounded to image field, and dimension of h*)
fun conv(dim,shape) =let
val expindex= specialize(shape,0)
in E.EIN{
params = [E.IMG(dim,shape),E.KRN],
index = shape,
body= E.Conv(0,expindex,1,[])
}
end

(* Probe: <F(x)>_{\alpha}   *)
fun probe(alpha,dim) = let
val  expindex= specialize(alpha,0)
in E.EIN{
params = [E.FLD(dim),E.TEN(0,[])], index= alpha,
body= E.Probe(E.Field(0, expindex), E.Tensor(1,[]))
}
end
(***************************** derivative ****************************)
(*\EinExp{\sum_{ij}\mathcal{E}_{ij} \frac{ F_j}{\partial x_i}*)
val curl2d=E.EIN{
params = [E.FLD 2],
index = [],
body = E.Sum([(E.V 0,0,1), (E.V 1,0,1)],
E.Opn(E.Prod,[E.G(E.Eps2(E.V 0, E.V 1)),
E.Apply(E.Partial([E.V 0]), E.Field(0,[E.V 1]))]))
}

val curl3d=E.EIN{
params = [E.TEN(1,[3])],
index = [3],
body = E.Sum([(E.V 1,0,2), (E.V 2,0,2)],E.Opn(E.Prod,[E.G(E.Epsilon(E.V 0, E.V 1, E.V 2)),
E.Apply( E.Partial([E.V 1]), E.Field(0,[E.V 2])) ]))
}

(*< d F /  d_i>_i  *)
val a=List.hd(alpha)
val  expindex= specialize(alpha,0)
in E.EIN{
params = [E.FLD(a)],
index =alpha,
body = E.Apply(E.Partial(expindex),E.Field(0,[]))
}
end

fun gradConstant (alpha as a::_) =
E.EIN{
params = [E.FLD a],
index = [],
body = E.Apply(E.Partial [(E.C (0, true))], E.Field(0, []))
}

(*< Sigma d F_alpha /  d x_i>ALpha  i CHANGE HERE *)
fun dotimes(dim,alpha)= let
val n=length(alpha)
fun expIndex(n,inc)=List.tabulate(n, (fn(x)=>E.V (x+inc)))
val i'=expIndex(n,0)
in E.EIN{
params = [E.FLD(dim)], index =alpha@[dim],
body = E.Apply(E.Partial [E.V n] ,E.Field(0,i'))
}
end

(*  <d F_i /d_i> *)
fun divergence(dim,alpha)=let
val  expindex= specialize(alpha,0)
val sumI=length(alpha)
val sumIndex=E.V(sumI)
val sumIndexL=[sumIndex]
val S=expindex@sumIndexL
in E.EIN{
params = [E.FLD(dim)],
index = alpha,
body = E.Sum([(sumIndex,0,dim-1)],
E.Apply(E.Partial(sumIndexL),E.Field(0,S)))
}
end

(*newbies*)
val divCR = E.EIN{
params = [E.TEN(subst_flag(),[])],
index = [],
body = E.Op2(E.Div, E.B(E.Const 1), E.Tensor(0, []))
}

end; (* local *)

end (* local *)