Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4410 - (download) (annotate)
Fri Aug 12 18:28:32 2016 UTC (2 years, 11 months ago) by cchiw
File size: 34563 byte(s)
added inverse
(* 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*)*)
 
     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)
    (******************************* 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*)
fun halfT(n) = E.EIN{
params = [E.TEN(subst_flag(),[3,3])],
index = [n,n],
body = E.Op2(E.Div, E.Tensor(0, [E.V 0 , E.V 1]), E.B(E.Const 2))
}
fun halfF(dim,n) = E.EIN{
params = [E.FLD(dim)],
index = [n,n],
body = E.Op2(E.Div, E.Field(0, [E.V 0 , E.V 1]), E.B(E.Const 2))
}


(*scale by delta*)
fun scaleIdT(n) = E.EIN{
params = [E.TEN(subst_flag(),[])],
index = [n,n],
body = E.Opn(E.Prod, [E.Tensor(0, []), E.G(E.Delta(E.V 0, E.V 1))])
}

(*scale by delta*)
fun scaleIdF(dim,n) = E.EIN{
params = [E.FLD(dim)],
index = [n,n],
body = E.Opn(E.Prod, [E.Field(0, []), E.G(E.Delta(E.V 0, E.V 1))])
}



    (* 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,
                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(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,
            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(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])]))]))])
        )
        }
        
    (************************* inverse **************************)
    fun mkFieldCx(f, ix, jx)= f(0, [mkCxSingle ix, mkCxSingle jx])
    fun mkFieldVx(f, ix, jx)= f(0, [E.V ix, E.V jx])
    fun mkAdd(es)= E.Opn(E.Add, es)
    fun mkSub(e1, e2)= E.Op2(E.Sub, e1, e2)
    fun mkProd(es) = E.Opn(E.Prod, es)
    fun mkDeltaCV(ix, jx)= E.G(E.Delta(mkCxSingle ix, E.V jx))
    fun mkEpsVx(ix, jx)= E.G(E.Eps2(E.V ix, E.V jx))
    
    fun mkInvf2(f) = let
        val f00 = mkFieldCx(f, 0,0)
        val f11 = mkFieldCx(f, 1,1)
        val f01 = mkFieldCx(f, 0,1)
        val f10 = mkFieldCx(f, 1,0)
        val i = 0
        val j = 1
        val k = 2
        val fij =   mkFieldVx(f, i,j)
        val fkk =   mkFieldVx(f, k,k)
        (* numerator*)
        (*
        val e1 = mkProd[f11, mkDeltaCV(0,i), mkDeltaCV(0,j)]
        val e2 = mkProd[f00, mkDeltaCV(1,i), mkDeltaCV(1,j)]
        val e3 = mkProd[fij, mkEpsVx(i,j), mkEpsVx(j,i)]
        val en = mkAdd[e1,e2,e3]
        *)
        val en = mkSub(mkProd[E.Sum([(E.V k, 0, 1)], fkk), E.G(E.Delta(E.V i, E.V j))], fij)
        
        (* denominator *)
        val d1 =   mkProd[f00,f11]
        val d2 =   mkProd[f01,f10]
        val dn =    mkSub(d1, d2)
        in
            E.Op2(E.Div, en, dn)
        end
    
    
    fun invF2(dim) = let
        val e =  mkInvf2(E.Field)
        val shape = [2,2]
        in
            E.EIN{params = [E.FLD dim], index= shape, body = e}
        end

    val invT2  = let
        val e =  mkInvf2(E.Tensor)
        val shape = [2,2]
        in
            E.EIN{params = [E.TEN (subst_flag(), shape)], index= shape, body = e}
    end

    
    (************************* 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  *)
    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
        
    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 *)

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0