OVERVIEW

This directory contains facilities for dynamic parsing, optimization,
differentiation, and evaluation of expressions.  This capability has many uses,
such as specifying analytic test cases from the commandline.  It is not intended
to be a full-featured language.  It supports branches, but it does not support
loops.  The language itself is described in LANGUAGE_DESCRIPTION.txt.

The main class is PROGRAM_COMPILER, which does most of the work.  This class
gathers the string containing the program and information about what the input
and output variables are.  The class then parses the string and generates an
internal representation, which is rather similar to machine code.  After
parsing, this code is converted to SSA form.  From here, the code can be
manipulated in a few ways, including optimization, differentiation, appending
instructions, appending, or duplication.  After any manipulations are completed,
the result should be optimized and finalized.

Finalizing the PROGRAM_COMPILER initializes a PROGRAM object.  After this point,
the PROGRAM_COMPILER object is no longer useful and should be destroyed.  The
PROGRAM object is a list of instructions similar to assembly language.  The
instructions are executed in order to evaluate the expression, but the PROGRAM
object itself is not changed.

The PROGRAM_CONTEXT is used as a staging area for executing a PROGRAM.  It
stores the inputs to the program, which must be populated before executing the
program.  It stores the registers that are used internally when executing the
program, and it stores the outputs that will be populated during program
execution.  This object is reusable.

EXAMPLE USAGE

  PROGRAM_COMPILER<T> pc;
  PROGRAM<T> prog;
  PROGRAM_CONTEXT<T> context;
  MEM_IN in_x=pc.Add_Input("x");
  MEM_IN in_t=pc.Add_Input("t");
  MEM_OUT out_u=pc.Add_Output("u");
  pc.Parse("u=cos(x)*t;");
  pc.Optimize();
  pc.Finalize(prog);
  context.Initialize(prog);
  context.Set(in_x,1.3);
  context.Set(in_t,2.7);
  prog.Execute(context);
  T u=context.Get(out_u);

BASIC USAGE

The Simple_Compile routine provides a concise wrapper around the most basic
usage of these facilities.  It follows the bare minimum steps that all users
must follow to set up a program.

1. Create a PROGRAM_COMPILER<T> pc.  This is the object that will do all of the
interesting work.  The constructor takes no arguments and just sets up a blank
slate.  This is a one-time-use class.  After compiling one expression, it must
be thrown away.  If you need to compile additional expressions, create a new
PROGRAM_COMPILER object for each.

2. Register your inputs and outputs by calling pc.Add_Input and pc.Add_Output.
The first argument is the name of the variable.  The version that does not take
an integer registers a scalar.  The version taking one integer registers a
vector of that size.  The version taking two integers registers a matrix.  The
object being registered is allocated a register (or contiguous block of
registers in the case of vectors or matrices).  These routines return a MEM_IN
or MEM_OUT object, which stores the index of the first register in the object.
These objects are ELEMENT_ID objects and as such support many of the usual
operations (i++, i+3, i<j, etc.) needed to manipulate them conveniently while
avoiding accidental misuse.

Each of the Add_Input and Add_Output routines takes an optional MEM_IN or
MEM_OUT argument. If this is set, then the object is registered starting at the
specified index, aliasing the objects already stored there.  Note the extra
argument may only be used to register new names to existing places in memory;
the entire memory block must already be allocated.  This allows a name like "X"
to refer to an input vector, while also allowing its components to be
conveniently referred to as "x", "y", and "z".

Indices are registered consecutively; input and output are indices into separate
arrays.  The first input will be assigned at MEM_IN(0), and the first output
will be assigned at MEM_OUT(0).

Zero-sized matrices or vectors are permitted, but they are not assigned any
space, and the index that is returned will be used to store something else.

3. Call pc.Parse, which compiles the string provided into its internal
representation.  The Parse routine is not threadsafe.  Only one Parse routine
may execute at a time, even if those routines are called on different
PROGRAM_COMPILER objects.  The parsing is done using LEX/YACC, which maintain
global state during the parsing process.

4. At this point, the PROGRAM_COMPILER object is in its working state.  It is
intended to be manipulate while in this state in various ways.  These will be
discussed later.  Although not required, you should generally call pc.Optimize
before moving on to the next step; this makes execution faster.

5. Call pc.Finalize to fill in the PROGRAM object.  At this point, the
PROGRAM_COMPILER object should be discarded; it cannot be reused.

6. Set up the PROGRAM_CONTEXT object by calling context.Initialize(prog); This
allocates space in the PROGRAM_CONTEXT object for inputs, outputs, and space for
intermediate computations.  If you want to execute the program in parallel, you
can initialize many PROGRAM_CONTEXT objects from the same PROGRAM.

7. Fill in the inputs in the PROGRAM_CONTEXT using the Set routines.  Use the
indices you got from the corresponding calls to Add_Input.

8. Call Execute to run the program.

9. Extract the outputs from the PROGRAM_CONTEXT using the Get routines.  Use the
indices you got from the corresponding calls to Add_Output.

The steps 7-9 can be repeated to run the program many times.  PROGRAM_CONTEXT is
reusable, and PROGRAM should be treated as constant after it has been
initialized.  Program execution is thread-safe, provided each execution has its
own PROGRAM_CONTEXT object.  The PROGRAM object can be shared.

DEBUGGING

These classes have a few debugging facilities to help track down issues.  Both
PROGRAM and PROGRAM_COMPILER can be printed using their Print() functions; these
dump out something that looks rather like machine language.  These routines are
also helpful to understand how various manipulations work. PROGRAM_COMPILER has
a debug flag and a verbose flag.  The debug flag causes (expensive) assertions
to be performed frequently.  This is very helpful in tracking down bugs, but it
is very expensive.  The verbose flag prints out the program after many stages in
the process.




    PROGRAM_COMPILER(const PROGRAM_COMPILER& pc);

    // Differentiation; Note that these routines do not create new output
    // variables, but you can do so with Make_Output.
    REG Diff(REG func_var,REG diff_var);
    ARRAY<REG,MEM_OUT> Diff(REG diff_var); // func_var is all outputs
    ARRAY<REG> Diff(ARRAY_VIEW<const REG> func_vars,REG diff_var);
    template<class R,class T_ARRAY,class I>
    ARRAY<REG,I> Diff(const ARRAY_BASE<R,T_ARRAY,I>& func_vars,REG diff_var)
    {
        if constexpr(is_same<R,REG>::value && HAS_POINTER<T_ARRAY>::value)
            return ARRAY<REG,I>(Diff(ARRAY_VIEW<const REG>(func_vars,force_idx),diff_var),force_idx);
        else
        {
            ARRAY<REG> fv(func_vars,cast_data,force_idx);
            return ARRAY<REG,I>(Diff(ARRAY_VIEW<const REG>(fv),diff_var),force_idx);
        }
    }
    ARRAY<REG> Diff(ARRAY_VIEW<const MEM_OUT> func_vars,REG diff_var);
    ARRAY<REG,MEM_IN> Gradient(REG func_var); // diff_var is all inputs
    // ret(out_var)(in_var)
    ARRAY<ARRAY<REG,MEM_IN>> Gradient(ARRAY_VIEW<const REG> func_vars); // diff_var is all inputs
    REG Directional_Derivative(REG func_var,ARRAY_VIEW<const REG,MEM_IN> dir);
    ARRAY<REG,MEM_OUT> Directional_Derivative(ARRAY_VIEW<const REG,MEM_IN> dir); // func_var is all outputs
    
    void Parse(const std::string& str);
    void Print() const;
    void Optimize();
    MEM_CONST Add_Constant(T x);
    void Finalize(PROGRAM<T>& prog);

    void Shift_Registers();
    void Rename_Variable(REG var_old,REG var_new);

    MEM_OUT Append(const PROGRAM_COMPILER<T>& pc);
    MEM_REG Num_Reg() const {return num_reg;}
    MEM_OUT Num_Out() const {return num_out;}
    MEM_IN Num_In() const {return num_in;}
    MEM_REG New_Reg() {return {num_reg++};}
    MEM_OUT New_Out() {return {num_out++};}
    
    // Makes new output, copies var to it.
    MEM_OUT Make_Output(REG var);
    template<class A,class I>
    auto Make_Output(const ARRAY<A,I>& var)
    {
        ARRAY<decltype(Make_Output(var(I{}))),I> a(var.m);
        for(I i{0};i<var.m;i++) a(i)=Make_Output(var(i));
        return a;
    }
    void Test_Invariants() const;

    // Turns outputs into regular registers.  Returns first such register and
    // the number of former outputs
    std::pair<MEM_REG,int> Move_Outputs_To_Reg();

    void Register_Default_Func();
    void Register_Func(const std::string& name,int num_in,DATA(*func)(PROGRAM_COMPILER<T>*pc,const ARRAY<DATA>& in));
    REG Append_Instruction(const OPERATION& type,REG src0,REG src1={},REG src2={})
    {return Before_Instruction(End_It(),type,src0,src1,src2);}

    // make_true and make_false should add instructions (with
    // Append_Instruction) to compute quantities.  The returned arrays must be
    // the same length.  By emitting instructions on demand, this routine allows
    // lazy execution of that code in the resulting program.
    ARRAY<REG> Append_If_Else(REG cond,std::function<ARRAY<REG>()> make_true,std::function<ARRAY<REG>()> make_false);





1. 
}


PROGRAM_COMPILER

namespace PhysBAM{
template<class T> ANALYTIC_PROGRAM_HELPER<T>::
ANALYTIC_PROGRAM_HELPER(const std::string& str,const std::string& out_name,
    int dim,int m,int n,bool want_ddx)
    :dim(dim)
{
    PHYSBAM_ASSERT(dim<=3);
    PROGRAM_COMPILER<T> pc;
    const char* axes[]={"x","y","z"};
    idx_f.Resize(dim);
    pc.Add_Input("X",dim);
    ARRAY<MEM_IN> X(dim);
    for(int i=0;i<dim;i++){
        X(i)=MEM_IN(i);
        pc.Add_Input(axes[i],MEM_IN(i));}
    MEM_IN t=pc.Add_Input("t");

    int dofs=0;
    if(m<0)
    {
        pc.Add_Output(out_name);
        dofs=1;
    }
    else if(n<0)
    {
        dofs=m;
        pc.Add_Output(out_name,m);
    }
    else
    {
        dofs=m*n;
        pc.Add_Output(out_name,m,n);
    }
    idx_f.Resize(dofs);
    for(int i=0;i<dofs;i++) idx_f(i)=MEM_OUT(i);
    pc.Parse(str);
    pc.Optimize();

    PROGRAM_COMPILER<T> pc_dt(pc),pc_dx(pc);
    idx_dt=pc_dt.Make_Output(pc_dt.Diff(idx_f,t));
    pc_dt.Optimize();

    if(dofs>1)
    {
        for(int i=0;i<dim;i++)
            idx_dx.Append(pc_dx.Make_Output(pc_dx.Diff(idx_f,X(i))));
    }
    else
    {
        auto grad=pc.Gradient(MEM_OUT{0});
        idx_dx.Resize(dim);
        for(int i=0;i<dim;i++)
            idx_dx(i)=pc.Make_Output(grad(MEM_IN{i}));
    }
    pc_dx.Optimize();

    if(want_ddx)
    {
        PROGRAM_COMPILER<T> pc_ddx(pc_dx);
        for(int i=0;i<dim;i++)
            idx_ddx.Append(pc_ddx.Make_Output(pc_ddx.Diff(idx_dx,X(i))));
        pc_ddx.Optimize();
        pc_ddx.Finalize(prog_ddx);
        context_ddx.Initialize(prog_ddx);
    }

    pc.Finalize(prog);
    pc_dt.Finalize(prog_dt);
    pc_dx.Finalize(prog_dx);
    context.Initialize(prog);
    context_dt.Initialize(prog_dt);
    context_dx.Initialize(prog_dx);
}
//#####################################################################
// Function f
//#####################################################################
template<class T> void ANALYTIC_PROGRAM_HELPER<T>::
f(ARRAY_VIEW<const T> X,T t) const
{
    assert(X.m==dim);
    for(int i=0;i<dim;i++) context.data_in(MEM_IN(i))=X(i);
    context.data_in(MEM_IN(dim))=t;
    prog.Execute(context);
}
//#####################################################################
// Function dt
//#####################################################################
template<class T> void ANALYTIC_PROGRAM_HELPER<T>::
dt(ARRAY_VIEW<const T> X,T t) const
{
    assert(X.m==dim);
    for(int i=0;i<dim;i++) context_dt.data_in(MEM_IN(i))=X(i);
    context_dt.data_in(MEM_IN(dim))=t;
    prog_dt.Execute(context_dt);
}
//#####################################################################
// Function dX
//#####################################################################
template<class T> void ANALYTIC_PROGRAM_HELPER<T>::
dX(ARRAY_VIEW<const T> X,T t) const
{
    assert(X.m==dim);
    for(int i=0;i<dim;i++) context_dx.data_in(MEM_IN(i))=X(i);
    context_dx.data_in(MEM_IN(dim))=t;
    prog_dx.Execute(context_dx);
}
//#####################################################################
// Function ddX
//#####################################################################
template<class T> void ANALYTIC_PROGRAM_HELPER<T>::
ddX(ARRAY_VIEW<const T> X,T t) const
{
    assert(X.m==dim);
    for(int i=0;i<dim;i++) context_ddx.data_in(MEM_IN(i))=X(i);
    context_ddx.data_in(MEM_IN(dim))=t;
    prog_ddx.Execute(context_ddx);
}
template class ANALYTIC_PROGRAM_HELPER<float>;
template class ANALYTIC_PROGRAM_HELPER<double>;
}



