root/trunk/blade/BladeSimplify.d

//  Written in the D programming language 1.0
/**
* Simplify a vector expression.
* Part of BLADE : Basic Linear Algebra D Expressions
*
*  The following transformations are performed:
* (A) Remove all duplicate symbols.
* (B) Create 'compound' variables, consisting of multiple symbols:
*    - Combine all scalars into a single scalar.
*    - Rank reduction: Combine slicing operations into vectors
* (C) Arithmetic transformations.
*    - Use slicing distributive law: Given A[B..C] for expressions A,B,C
*      where B and C are both rank 0, and A is rank 1 or more, the slice can
*      be moved to every vector inside A.
*    - Use associativity of *: A*(B*C[]) == (A*B)*C[] (Not strictly true for
*      floating point; results may differ by 1ulp,
*       eg (1.3L*3.1L)*4.7L < 1.3L*(3.1L*4.7L)
*      Note that floating point addition is not associative at all).
*    - Remove unary minus where possible, eg A-(-B) => A+B, abs(-A) => abs(A).
*    - Use associativity of * in intrinsics:
*         sum(A*V) => A*sum(V), abs(A*B) => abs(A)*abs(B)
* (D) Expression standardisation
*    - Move multiplies to left: Convert A[]*B into B*A[] (assumes * is commutative,
*      not valid for quaternions).
*    - Convert C-A*B into C+(-A)*B whenever possible.
*
* Author:
*   Don Clugston.
* License:
*   Public domain.
*/

module blade.BladeSimplify;

public import blade.SyntaxTree : AbstractSyntaxTree, Symbol;
private import blade.BladeVisitor;
private import blade.BladeRank : exprRank, subexprRank, getRankErrorText;


/// A simplified vector expression
struct RevisedExpression {
    char [] braceExpression; // the expression with compounds in braces
    char [] expression; // the revised expression using new variable names
                        // (so, for example, B+=(D-F) becomes A+=(B-C) ).
    Symbol[] symbolTable; // the original symbol table, with all the old names
    char [][] compounds; // the compound variables, defined using the old names
    char [] rank;   // rank of all symbols (including original & compounds)
    char [] mapping; // mapping from new names onto old names.
    char [] errorMessage; // or null if no errors.
}

// TODO: ".re", ".im" should join this list.
bool isBladeIntrinsic(char [] str)
{
    return str=="dot" || str=="sum" || str=="max" || str=="min"
           || str=="abs" || str=="sqrt" || str=="prod";
}


// Given an abstract syntax tree, returns a RevisedExpression.
RevisedExpression simplifySyntaxTree(AbstractSyntaxTree tree)
{
    char [] ranks;
    char [] err="";
    for (int i=0; i<tree.symbolTable.length; ++i) {
        ranks~=tree.symbolTable[i].rank;
        if (tree.symbolTable[i].type=="" && !isBladeIntrinsic(tree.symbolTable[i].value)) err ~= " " ~ tree.symbolTable[i].value;
    }
    // Check for undefined symbols
    if (err.length > 0)
        return RevisedExpression(tree.expression, "", tree.symbolTable, [""], "","", "Undefined symbols:" ~ err);
    else {
        // Remove duplicate symbols, convert intrinsics
        char [] expr2 = removeDuplicates(tree);
        // Check for rank errors
        int wholerank = exprRank(expr2, ranks);
        if (wholerank<0)
            return RevisedExpression(expr2, "", tree.symbolTable, [""], "","", getRankErrorText(wholerank));
        // Perform scalar foldings and dimension reduction
        char [] expr3 = foldScalars(foldIndices(expr2, ranks), ranks);
        return remapCompounds(expr3, ranks, tree.symbolTable);
    }
}

private:

/* Adjust the expression to remove all references to duplicated symbol table
 * entries. (Duplicates can occur as a result of resolving aliases or constants).
 * Also move intrinsics in-line.
*/
char [] removeDuplicates(AbstractSyntaxTree tree)
{
    int numdups=0;
    char [] mapping = ""; // The new letter which this symbol should become,
                          // or '!' if it is an intrinsic
    for (int i=0; i<tree.symbolTable.length; ++i) {
        char c = 'A'+i;
        if (isBladeIntrinsic(tree.symbolTable[i].value)) {
            ++numdups;
            c = '!';
        } else {
            for (int j=0; j<i; ++j) {
                if (tree.symbolTable[i].value==tree.symbolTable[j].value) {
                    ++numdups;
                    c = ('A'+j);
                    break;
                }
            }
        }
        mapping ~= c;
    }
    if (numdups==0) return tree.expression;
    char [] e = "";
    for (int i=0; i<tree.expression.length;++i) {
        char c = tree.expression[i];
        if (c>='A' && c<='Z') {
            if (mapping[c-'A']=='!') e~=tree.symbolTable[c-'A'].value;
            else e~= mapping[c-'A'];
        } else e~=c;
    }
    return e;
}

unittest {
    AbstractSyntaxTree t = AbstractSyntaxTree("A+(B*C)", [Symbol("int", "125", 0),
    Symbol("int", "7", 0), Symbol("int", "125", 0)]);
    assert(removeDuplicates(t)=="A+(B*A)");
}

// Determine rank of a multidimensional index
int indexRank(char [] s)
{
   int r=0;
   int numbrack=0;
   int paren = 0;
   for(int i=1; i<s.length; ++i) {
        if (s[i]=='(') ++paren;
        else if (s[i]==')') --paren;
        if (paren==0 && s[i]==']') { numbrack--; }
        if (paren==0 && s[i]=='[') {
            if (numbrack==0) ++r;
            numbrack++;
        }
        if (paren==0 && numbrack==1 && s[i]==',') ++r; // commas increase the rank
        if (paren==0 && numbrack==1 && s[i]=='.' && s[i-1]=='.') {
            // if it's a slice, it does not increase the rank
             r--;
        }
   }
   return r;
}

public:
// Remove everything inside {} from expr, and create new variables for it.
RevisedExpression remapCompounds(char [] expr, char [] rank, Symbol[] symTable)
{
    char [][] comp;
    char [] used = ""; // which of the old variables are used; gives the new mapping
                  // of the name, or - if not used.
    for (int i=0; i<rank.length; ++i) used ~= "-";
    char [] r;
    char next = cast(char)('A' + rank.length);
    char [] e = "";
    for (int i=0; i<expr.length; ++i) {
        if (expr[i]=='{') {
            int k;
            int bracecount=1;
            for (k=i+1; bracecount>0; ++k) {
                if (expr[k]=='{') ++bracecount;
                if (expr[k]=='}') --bracecount;
            }
            --k;
            char [] newexpr = expr[i+1..k]; // strip off the {}
            int newi = k;
            if (i>0 && k<expr.length-1 && expr[i-1]=='(' && expr[k+1]==')') {
                e = e[0..$-1]; // remove the last '('
                newi=k+1; // don't add ')'
            }
            // Check for a duplicate
            int z;
            for (z=0; z<comp.length && comp[z]!=newexpr; ++z) {}
            if (z==comp.length) {
                e ~= next;
                ++next;
                comp ~= expr[i+1..k]; // strip off the {}
                if (expr[k-1]==']') {
                    // it's a vector/matrix of some kind, with rank reduced
                    // by indices. Can't just use exprRank, because the []
                    // aren't wrapped by ().
                    r ~= (rank[expr[i+1]-'A'] - indexRank(expr[i+1..k]));
                } else {
                    // it's a scalar expression. Note that it could involve
                    // a vector expression.
                    r~='0';
                }
            } else e ~= cast(char)('A'+z+rank.length);
            i = newi;
        } else {
            e ~= expr[i];
            if (expr[i]>='A' && expr[i]<='Z') used[expr[i]-'A']=expr[i];
        }
    }
    // Create a mapping from old to new variable names

    char [] old_ranks = "";
    char [] mapping="";
    char knt = 'A';
    for (int i=0; i<used.length; ++i) {
        if (used[i]!='-') {
            mapping ~= ('A'+i);
            old_ranks ~= rank[i];
            used[i] = knt;
            ++knt;
        }
    }
    for (int i=0; i<r.length; ++i) {
        mapping ~= ('A'+used.length+i);
    }
   // and set the expression to use the new names
    char [] f = "";
    for (int i=0; i<e.length; ++i) {
        char c = e[i];
        if (c>='A' && c<='Z') {
            if ((c-'A')<rank.length) f ~= used[c-'A'];
            else f~= (c-'A') - rank.length + knt;
        } else f ~= c;
    }
    return RevisedExpression(expr, f, symTable, comp, old_ranks~r, mapping);
}

private:
RevisedExpression simplifyVectorExpression(char [] expr, char [] rank, Symbol[] symTable)
{
    return remapCompounds(foldScalars(foldIndices(expr, rank), rank), rank, symTable);
}

unittest {
    RevisedExpression e = simplifyVectorExpression("A+=(((D[B])*C)[B])", "2004",[]);
    assert(e.rank=="202");
    assert(e.compounds[0]=="D[B,B]");
    assert(e.mapping=="ACE");
    assert(e.expression== "A+=(C*B)");
}

// -----------------------------------------------------------

// Given an array of slices/indices where
//  slicing[0]= start of slice, or index
//  slicing[1]= end of slice, or "" if it's an index,
// combine everything into a single slicing expression.
char [] createMultiSlice(char [][2][] slicing)
{
    char [] s="[";
    for (int i=0; i<slicing.length;++i) {
        if (i>0) s~= ",";
        s ~= slicing[i][0];
        if (slicing[i][1].length>0) s~=".." ~ slicing[i][1];
    }
    return s~"]";
}

// Combines all the indexing and slicing operations together (dimension reduction).
// Multiplication of sliced matrices and/or vectors is dimensionally
// reduced where possible (may even be converted into dot product).
// Returns the new expression. This eliminates all unnecessary slice operations.
// Furthermore, *any* value followed by '[' should be used as a new compound.
struct IndexFoldingVisitor {
    alias typeof(*this) This;
    alias char [] ReturnType;
    char [] rank;
    char [] dollar;
    char [][2][] slicing; // the indexing which applies to this complete expr.
static:
    ReturnType onVisitSymbol(This this_, char[] sym) {
       if (this_.slicing.length==0) {
           if (sym=="$") {
               return this_.dollar;
           }
           else return sym;
       } else {
           assert(sym!="$" && this_.rank[sym[0]-'A']>'0', "Rank error " ~ sym);
           // Note: Later, we'll want this to be a new terminal.
           return sym ~ createMultiSlice(this_.slicing);
       }
    }
    ReturnType onVisitFunction(This this_, char [] func, char [][] args) {
        // Intrinsics have no effect on slicing
        char [] result = "";
        for (int i=0; i<args.length; ++i) {
            if (i>0) result ~= ",";
            result ~= wrapInParens(doVisit(this_,args[i]));
        }
        return func ~ "(" ~ result ~ ")";
    }
    ReturnType onVisitPrefix(This this_, char [] op, char [] expr) {
        assert(this_.slicing.length==0, "BLADE ICE");
        return op ~ wrapInParens(doVisit(this_, expr));
    }
    ReturnType onVisitPostfix(This this_, char [] op, char [] expr) {
        assert(this_.slicing.length==0, "BLADE ICE");
        return wrapInParens(doVisit(this_, expr)) ~ op;
    }
    // Includes multi-dimensional slicing and indexing.
    ReturnType onVisitIndex(This this_, char [] base, char [][2][] slices) {
        if (slices.length==0) { // []  -- has no effect.
             return doVisit(this_, base);
         }
//        printf("  %.*s --> %.*s %.*s\n", base, slices[0][0], slices[0][1]);
        if (this_.slicing.length >0 && slices[$-1][1].length>0) {
            // the new dimension block ends with a slice. This needs to be combined
            // with the earliest existing dimension.
            // * If the existing dimension is an index,
            //   it might contain a dollar, which we need to replace.
            // * If the existing dimension is a slice, the two slices will combine.
            //
            // The items inside the slice are top-level, ie have no slice or dollar.
            char [] a = wrapInParens(doVisit(IndexFoldingVisitor(this_.rank,"$",[]), slices[$-1][0]));
            char [] b = wrapInParens(doVisit(IndexFoldingVisitor(this_.rank,"$",[]), slices[$-1][1]));
            char [] dollr = b ~ "-" ~ a;
            char [] c = wrapInParens(doVisit(IndexFoldingVisitor(this_.rank,dollr,[]), this_.slicing[0][0]));
            char [][2][] newslice=[];
             if (this_.slicing[0][1].length>0) { // slicing a slice
               if (b=="$" && this_.slicing[0][1]=="$") {
                   // very common special case, where both are sliced from end
                    newslice ~= ["(" ~ a ~ "+" ~ c ~")","$"];
                } else {
                   char [] d = wrapInParens(doVisit(IndexFoldingVisitor(this_.rank,dollr,[]), this_.slicing[0][1]));
                    newslice ~= ["(" ~ a ~ "+" ~ c~ ")", "(" ~ a ~ "+" ~ d~ ")"];
                }
            } else {
                newslice ~= [a ~ "+" ~ c, ""];
            }
            if (slices.length>1) {
                // append other slices, if any.
                return doVisit(IndexFoldingVisitor(this_.rank, "$", slices[0..$-1] ~ newslice ~ this_.slicing[1..$]), base);
            } else {
                return doVisit(IndexFoldingVisitor(this_.rank, "$",newslice ~ this_.slicing[1..$]), base);
            }
        } else { // just append them.
            return doVisit(IndexFoldingVisitor(this_.rank, "$", slices ~ this_.slicing), base);
        }
    }
    ReturnType onVisitBinaryOp(This this_, char [] op, char [] left, char [] right) {
        int lrank = subexprRank(left, this_.rank);
        int rrank = subexprRank(right, this_.rank);
        char [] first="";
        char [] second="";
        if ((op=="*" || op=="*=") && this_.slicing.length>0) {
            // If one of these is a matrix, the slicing gets interesting...
            // .. extremely so for slicing of matrix chain multiplication.
            if (lrank==0) {
                // All dimensions apply to right operand
                first = wrapInParens(doVisit(IndexFoldingVisitor(this_.rank, this_.dollar,[]), left));
                second = wrapInParens(doVisit(this_, right));
            } else if (rrank==0) {
                // All dimensions apply to the left operand
                first = wrapInParens(doVisit(this_, left));
                second = wrapInParens(doVisit(IndexFoldingVisitor(this_.rank, this_.dollar,[]), right));
            } else {
                assert(lrank>0 && rrank>0 && lrank<=2 && rrank<=2, "BLADE ICE: Tensor*tensor is unsupported");
                bool isDotProduct = false; // was it reduced to a dot product?

                // In the case of chained matrix multiplies, we can end up with an empty slice.
                if (this_.slicing.length>0 && this_.slicing[$-1][0]=="") {
                    this_.slicing=this_.slicing[0..$-1];
                }
                if (lrank==2) {
                    // First dimension applies to rows of the left operand
                    // If it's a slice, it will be a strided slice -- unless
                    // it comes from another matrix multiply, in which case the
                    // stride will drop out. (A[x]*B is strided).
                    char [][2][] newslice=[];
                    newslice ~= this_.slicing[0];
                    newslice ~= ["",""];
                    first = wrapInParens(doVisit(IndexFoldingVisitor(this_.rank, this_.dollar,newslice), left));
                } else {
                    assert(this_.slicing.length==1, "BLADE ICE: Rank error");
                    first = wrapInParens(doVisit(IndexFoldingVisitor(this_.rank, this_.dollar,[]), left));
                }
                if (lrank==2 && rrank==2) {
                    // Matrix * matrix
                    if (this_.slicing.length>1) {
                        // Second dimension applies to the right operand.
                        second = wrapInParens(doVisit(IndexFoldingVisitor(this_.rank, this_.dollar,this_.slicing[1..$]), right));
                        if (this_.slicing[0][1].length==0 && this_.slicing[1][1].length==0) {
                            // It's indices in both cases -- so it's a dot product.
                            isDotProduct = true;
                        }
                    } else {
                        second = wrapInParens(doVisit(IndexFoldingVisitor(this_.rank, this_.dollar,[]), right));
                    }
                } else if (lrank==1 && rrank==2) {
                    // vector * matrix, Matrix uses all the slicing
                    second = wrapInParens(doVisit(this_, right));
                    if (this_.slicing[0][1].length==0)  isDotProduct = true;
                } else if (lrank==2 && rrank==1) {
                    // matrix * vector, vector is unsliced.
                    second = wrapInParens(doVisit(IndexFoldingVisitor(this_.rank, this_.dollar,[]), right));
                    if (this_.slicing[0][1].length==0)  isDotProduct = true;
                } else assert(0, "BLADE ICE");
                if (isDotProduct) {
                  return "dot(" ~ first ~ "," ~ second ~ ")";
                }
            }
        } else { // not a multiplication
            return wrapInParens(doVisit(this_, left)) ~ op ~ wrapInParens(doVisit(this_, right));
        }
        return first ~ op ~ second;
    }
}

char [] foldIndices(char [] expr, char [] ranks)
{
    return beginVisit(IndexFoldingVisitor(ranks,"$",[]), expr);
}

unittest {
    assert(foldIndices("((A[C..D])+B)[($-E)]", "21000")=="(A[C+((D-C)-E)])+(B[($-E)])");
    assert(foldIndices("(A[C])[D]", "3100")=="A[C,D]");
    assert(foldIndices("(A[B..C])[D]", "3000")=="A[B+D]");
    assert(foldIndices("(A[B])[C..D]", "3000")=="A[B,C..D]");
    assert(foldIndices("((A[$])[(C-$)])[D]", "3000")=="A[$,(C-$),D]");
    assert(foldIndices("(A[B..$])[C..$]", "3000")=="A[(B+C)..$]");
    assert(foldIndices("((A[])[C..$])[]", "3000")=="A[C..$]");
    assert(foldIndices("((A[])[(B[C])..$])[]", "3100")=="A[(B[C])..$]");
    assert(foldIndices("A[,B..$,C]", "300")=="A[,B..$,C]");
    // Multidimensional slicing
    assert(foldIndices("(C*((A*B)[C]))[D]", "2200")=="C*dot((A[C,]),(B[D]))");
    assert(foldIndices("(A*B)[C..D,D]", "2200")=="(A[C..D,])*(B[D])");
    assert(foldIndices("(A*B)[C..D]", "2200")=="(A[C..D,])*B");
    assert(foldIndices("(A*B)[C..D]", "2100")=="(A[C..D,])*B");
    assert(foldIndices("(A*B)[C..D]", "1200")=="A*(B[C..D])");
    assert(foldIndices("(A*B)[C]", "120")=="dot(A,(B[C]))");

    assert(foldIndices("((A*B)*C)[D]", "2220")=="((A[D,])*B)*C");
    assert(foldIndices("((A+B)*C)[D]", "2220")=="((A[D,])+(B[D,]))*C");
    assert(foldIndices("((D*A)*B)[C]", "2100")=="dot((D*(A[C,])),B)");
    assert(foldIndices("(((A*B)*C)[D..E])[D]", "12200")=="dot((A*B),(C[D+D]))");
    assert(foldIndices("A+=(((D[B])*C)[B])", "2004")=="A+=((D[B,B])*C)");
    assert(foldIndices("dot(A,(A*dot(A,A)))","1")=="dot(A,(A*dot(A,A)))");
}

struct ScalarFold
{
    char [] expr;       // vector or matrix expression; empty for a pure scalar expression
    char [] multiplier; // scalar multiply of the entire expression. or "-" for unary minus
}