(* creates EIN operators 
 *
 * This code is part of the Diderot Project (http://diderot-language.cs.uchicago.edu)
 *
 * COPYRIGHT (c) 2015 The University of Chicago
 * All rights reserved.
 *)

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*)
    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)
    (******************************* Add *****************************************)
    val addRR = E.EIN{
        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} *)
    fun addTT 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.Opn(E.Add,[E.Tensor(0, expindex), E.Tensor(1, expindex)])
            }
        end
        
    (* mkField functions*)
    (*Adding Fields : < F{\alpha} + G_{\alpha}>_{\alpha} *)
    fun addFF(dim,shape) =let
        val expindex= specialize(shape,0)
        in E.EIN{
            params = [E.FLD(dim),E.FLD(dim)],
            index = shape,
            body = E.Opn(E.Add,[E.Field(0, expindex),E.Field(1, expindex)])
            }
        end
        
    (*Tensor and Fields*)
    fun addTF(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.Add,[E.Lift(E.Tensor(0, expindex)),E.Field(1, expindex)])
            }
        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,
                body = E.Opn(E.Add,[E.Lift(E.Tensor(0, expindex)),E.Op1(E.Neg,E.Field(1, expindex))])
            }
        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 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(0, 1, 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(0, 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(0, 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(0, 1, 2)),
                E.Field(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 expindexA= specialize(alpha,0)
        val expindexB= specialize(beta,length(alpha))
            in E.EIN{
            params = [E.FLD(dim),E.TEN(subst_flag,alpha)],
            index = alpha@beta,
            body = E.Opn(E.Prod,[E.Field(0, expindexA),E.Lift(E.Tensor(1, expindexB))])
            }
        end
*)
(********************  outerProduct ********************************)

(*Vector Examples : <T_i * T_j>_ij..t0⊗t1*)
    fun outerTT(dimA,dimB) =
    E.EIN{
    params = [E.TEN(subst_flag,[dimA]), E.TEN(subst_flag,[dimB])],
    index= [dimA,dimB],
    body= E.Opn(E.Prod,[E.Tensor(0, [E.V 0]), E.Tensor(1, [E.V 1])])
    }
    fun outerTM(i,[j,k]) =
    E.EIN{
    params = [E.TEN(subst_flag,[i]), E.TEN(subst_flag,[j,k])],
    index= [i,j,k],
    body= E.Opn(E.Prod,[E.Tensor(0, [E.V 0]), E.Tensor(1, [E.V 1,E.V 2])])
    }
    
    fun outerMT([i,j],k) =
    E.EIN{
    params = [E.TEN(subst_flag,[i,j]), E.TEN(subst_flag,[k])],
    index= [i,j,k],
    body= E.Opn(E.Prod,[E.Tensor(0, [E.V 0,E.V 1]), E.Tensor(1, [E.V 2])])
    }
    (*Assumes same dimension vector field *)
    fun outerFF(dim,i,j) =
    E.EIN{
    params = [E.FLD(dim),E.FLD(dim)],
    index = [i,j],
    body = E.Opn(E.Prod,[E.Field(0, [E.V 0]),E.Field(1, [E.V 1])])
    }

    (*************************** 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),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.Op1(E.Sqrt,E.Sum(sx, E.Opn(E.Prod,[E.Field(0, expindex),E.Field(1,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,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,dim-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,
                E.Opn(E.Add,[E.Opn(E.Prod,[a,e,i]),E.Opn(E.Prod,[b,f,g]),E.Opn(E.Prod,[c,d,h])]),
                E.Opn(E.Add,[E.Opn(E.Prod,[c,e,g]), E.Opn(E.Prod,[b,d,i]),E.Opn(E.Prod,[a,f,h])]))
            }
        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(0, 1, 2)),E.Tensor(0,[mkCxSingle 2, E. V 2])]))]))])
                            )
    }
    
    val detF2  =E.EIN{
        params = [E.FLD 2], 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,
            E.Opn(E.Add,[E.Opn(E.Prod,[a,e,i]),E.Opn(E.Prod,[b,f,g]),E.Opn(E.Prod,[c,d,h]) ]),
            E.Opn(E.Add,[E.Opn(E.Prod,[c,e,g]), E.Opn(E.Prod,[b,d,i]),E.Opn(E.Prod,[a,f,h])]))
        }
    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(0, 1, 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])]))
        }
        
    val detF3=
        E.EIN{ params = [E.FLD(3)], 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(0, 1, 2)),E.Field(0,[mkCxSingle 2, E. V 2])]))]))])
                    )
        }
        
        
    (************************* 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 *)
    
    (************************* other tensor ops *************************)
    fun modulate dim =E.EIN{
        params = [E.TEN(subst_flag,[dim]), E.TEN(subst_flag,[dim])],
        index = [dim],
        body = E.Opn(E.Prod,[E.Tensor(0, [E.V 0]), E.Tensor(1, [E.V 0])])
    }
     
    fun identity dim =E.EIN{
        params = [],index = [dim,dim], body = E.G(E.Delta(E.V(0), E.V(1)))
    }
        
    fun slice (argTy,mask,const,rstTy) =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)
        in
            E.EIN{  params = [E.TEN(subst_flag,argTy)],index = rstTy,body = E.Tensor(0, ix)}
        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(0, 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(0, 1, 2)),
             E.Apply( E.Partial([E.V 1]), E.Field(0,[E.V 2])) ]))
        }
             
    (*< d F /  d_i>_i  *)
    fun grad alpha=let
        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
            
    (*< 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 *)