ginac_tools.cpp

Go to the documentation of this file.

// Copyright (C) 2006-2009 Kent-Andre Mardal and Simula Research Laboratory
//
// This file is part of SyFi.
//
// SyFi is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 2 of the License, or
// (at your option) any later version.
//
// SyFi is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with SyFi. If not, see <http://www.gnu.org/licenses/>.

#include "ginac_tools.h"
#include "symbol_factory.h"

#include <ginac/ginac.h>

#include <iostream>
#include <sstream>
#include <fstream>
#include <stdexcept>

using namespace std;
using namespace GiNaC;

namespace SyFi
{

        GiNaC::lst cross(GiNaC::lst& v1, GiNaC::lst& v2)
        {
                GiNaC::lst ret;
                if ( v1.nops() != v2.nops() )
                {
                        cout <<"incompatible vectors "<<endl;
                        cout <<"v1.nops() "<<v1.nops();
                        cout <<"  v2.nops() "<<v2.nops()<<endl; ;
                        return GiNaC::lst();
                }
                ret.append(  v1.op(1)*v2.op(2) - v1.op(2)*v2.op(1));
                ret.append(- v1.op(0)*v2.op(2) + v1.op(2)*v2.op(0));
                ret.append(  v1.op(0)*v2.op(1) - v1.op(1)*v2.op(0));
                return ret;
        }

        GiNaC::ex inner(GiNaC::ex a, GiNaC::ex b, bool transposed)
        {
                if (GiNaC::is_a<GiNaC::matrix>(a) && GiNaC::is_a<GiNaC::matrix>(b))
                {
                        GiNaC::matrix ma = GiNaC::ex_to<GiNaC::matrix>(a);
                        GiNaC::matrix mb = GiNaC::ex_to<GiNaC::matrix>(b);
                        if ( !transposed )
                        {
                                if (ma.cols() != mb.cols() || ma.rows() != mb.rows() )
                                {
                                        cout <<"Incompatible matrices "<<endl;
                                        cout <<"a.cols() "<<ma.cols()<<endl;
                                        cout <<"a.rows() "<<ma.rows()<<endl;
                                        cout <<"b.cols() "<<mb.cols()<<endl;
                                        cout <<"b.rows() "<<mb.rows()<<endl;
                                        cout <<"a="<<a<<endl;
                                        cout <<"b="<<b<<endl;
                                        throw std::runtime_error("Incompatible matrices.");
                                }

                                GiNaC::ex ret;
                                for (unsigned int i=0; i<ma.rows(); i++)
                                {
                                        for (unsigned int j=0; j<ma.cols(); j++)
                                        {
                                                ret += ma(i,j)*mb(i,j);
                                        }
                                }
                                return ret;
                        }
                        else
                        {
                                if (ma.cols() != mb.rows() || ma.rows() != mb.cols() )
                                {
                                        cout <<"Incompatible matrices "<<endl;
                                        cout <<"a.cols() "<<ma.cols()<<endl;
                                        cout <<"a.rows() "<<ma.rows()<<endl;
                                        cout <<"b.cols() "<<mb.cols()<<endl;
                                        cout <<"b.rows() "<<mb.rows()<<endl;
                                        cout <<"a="<<a<<endl;
                                        cout <<"b="<<b<<endl;
                                        throw std::runtime_error("Incompatible matrices.");
                                }

                                GiNaC::ex ret;
                                for (unsigned int i=0; i<ma.rows(); i++)
                                {
                                        for (unsigned int j=0; j<ma.cols(); j++)
                                        {
                                                ret += ma(i,j)*mb(j,i);
                                        }
                                }
                                return ret;
                        }
                }
                else if (GiNaC::is_a<GiNaC::lst>(a)
                        && GiNaC::is_a<GiNaC::lst>(b))
                {
                        return inner(GiNaC::ex_to<GiNaC::lst>(a), GiNaC::ex_to<GiNaC::lst>(b));
                }
                else
                {
                        return a*b;
                }
        }

        GiNaC::ex inner(GiNaC::lst v1, GiNaC::lst v2)
        {
                GiNaC::ex ret;

                if ( v1.nops() != v2.nops() )
                {
                        cout <<"incompatible vectors "<<endl;
                        cout <<"v1.nops() "<<v1.nops();
                        cout <<"  v2.nops() "<<v2.nops()<<endl; ;
                        return 0;
                }
                for (unsigned i = 0; i <= v1.nops()-1 ; ++i)
                {
                        if ( GiNaC::is_a<GiNaC::lst>(v1.op(i)) &&
                                GiNaC::is_a<GiNaC::lst>(v2.op(i)) )
                        {
                                ret += inner(GiNaC::ex_to<GiNaC::lst>(v1.op(i)),
                                        GiNaC::ex_to<GiNaC::lst>(v2.op(i)));
                        }
                        else
                        {
                                ret += v1.op(i)*v2.op(i);
                        }
                }
                return ret;
        }

        GiNaC::ex inner(GiNaC::exvector& v1, GiNaC::exvector& v2)
        {
                GiNaC::ex ret;
                for (unsigned int i=0; i< v1.size(); i++)
                {
                        ret += v1[i]*v2[i];
                }
                return ret;
        }

        GiNaC::lst matvec(GiNaC::matrix& M, GiNaC::lst& x)
        {
                GiNaC::lst ret;
                int nr = M.rows();
                int nc = M.cols();
                for (int i = 0; i < nr; i++)
                {
                        GiNaC::ex tmp;
                        for (int j = 0; j < nc; j++)
                        {
                                tmp = tmp +  M(i,j)*(x.op(j));
                        }
                        ret.append(tmp);
                }
                return ret;
        }

        GiNaC::ex matvec(GiNaC::ex A, GiNaC::ex x)
        {
                ex sol;

                if (GiNaC::is_a<GiNaC::matrix>(A) && GiNaC::is_a<GiNaC::matrix>(x))
                {
                        GiNaC::matrix AA = GiNaC::ex_to<GiNaC::matrix>(A);
                        GiNaC::matrix xx = GiNaC::ex_to<GiNaC::matrix>(x);
                        sol = AA.mul(xx);
                }
                else
                {
                        throw std::runtime_error("Invalid argument types, need matrices");
                }
                return sol;
        }

        GiNaC::lst ex2equations(GiNaC::ex rel)
        {
                GiNaC::ex lhs = rel.lhs();
                GiNaC::ex rhs = rel.rhs();

                GiNaC::ex l;
                GiNaC::ex r;

                GiNaC::lst eqs;

                for (int i=lhs.ldegree(x); i<=lhs.degree(x); ++i)
                {
                        for (int j=lhs.ldegree(y); j<=lhs.degree(y); ++j)
                        {
                                for (int k=lhs.ldegree(z); k<=lhs.degree(z); ++k)
                                {
                                        l = lhs.coeff(x,i).coeff(y, j).coeff(z,k);
                                        r = rhs.coeff(x,i).coeff(y, j).coeff(z,k);
                                        //      if (! (l == 0 && r == 0 ) )  eqs.append(l == r); OLD VERSION
                                        if ( (l != 0 && (r == 0 || r == 1) ) )  eqs.append(l == r);
                                }
                        }
                }
                eqs.sort();
                return eqs;
        }

        GiNaC::lst collapse(GiNaC::lst l)
        {
                GiNaC::lst lc;
                GiNaC::lst::const_iterator iter1, iter2;

                for (iter1 = l.begin(); iter1 != l.end(); ++iter1)
                {
                        if (GiNaC::is_a<GiNaC::lst>(*iter1))
                        {
                                for (iter2 = GiNaC::ex_to<GiNaC::lst>(*iter1).begin(); iter2 != GiNaC::ex_to<GiNaC::lst>(*iter1).end(); ++iter2)
                                {
                                        lc.append(*iter2);
                                }
                        }
                        else
                        {
                                lc.append(*iter1);
                        }
                }
                lc.sort();
                lc.unique();
                return lc;
        }

        GiNaC::matrix equations2matrix(const GiNaC::ex &eqns, const GiNaC::ex &symbols)
        {

                GiNaC::matrix sys(eqns.nops(),symbols.nops());
                GiNaC::matrix rhs(eqns.nops(),1);
                GiNaC::matrix vars(symbols.nops(),1);

                for (size_t r=0; r<eqns.nops(); r++)
                {
                                                                 // lhs-rhs==0
                        const GiNaC::ex eq = eqns.op(r).op(0)-eqns.op(r).op(1);
                        GiNaC::ex linpart = eq;
                        for (size_t c=0; c<symbols.nops(); c++)
                        {
                                const GiNaC::ex co = eq.coeff(GiNaC::ex_to<GiNaC::symbol>(symbols.op(c)),1);
                                linpart -= co*symbols.op(c);
                                sys(r,c) = co;
                        }
                        linpart = linpart.expand();
                        rhs(r,0) = -linpart;
                }
                return sys;
        }

        void matrix_from_equations(const GiNaC::ex &eqns, const GiNaC::ex &symbols, GiNaC::matrix& A, GiNaC::matrix& b)
        {
                // build matrix from equation system
                GiNaC::matrix sys(eqns.nops(),symbols.nops());
                GiNaC::matrix rhs(eqns.nops(),1);
                GiNaC::matrix vars(symbols.nops(),1);

                for (size_t r=0; r<eqns.nops(); r++)
                {
                                                                 // lhs-rhs==0
                        const GiNaC::ex eq = eqns.op(r).op(0)-eqns.op(r).op(1);
                        GiNaC::ex linpart = eq;
                        for (size_t c=0; c<symbols.nops(); c++)
                        {
                                const GiNaC::ex co = eq.coeff(GiNaC::ex_to<GiNaC::symbol>(symbols.op(c)),1);
                                linpart -= co*symbols.op(c);
                                sys(r,c) = co;
                        }
                        linpart = linpart.expand();
                        rhs(r,0) = -linpart;
                }
                A = sys;
                b = rhs;
        }

        GiNaC::ex lst_to_matrix2(const GiNaC::lst& l)
        {
                GiNaC::lst::const_iterator itr, itc;

                // Find number of rows and columns
                size_t rows = l.nops(), cols = 0;
                for (itr = l.begin(); itr != l.end(); ++itr)
                {
                        if (!GiNaC::is_a<GiNaC::lst>(*itr))
                                //              throw (std::invalid_argument("lst_to_matrix: argument must be a list of lists"));
                                cols = 1;
                        if (itr->nops() > cols)
                                cols = itr->nops();
                }
                // Allocate and fill matrix
                GiNaC::matrix &M = *new GiNaC::matrix(rows, cols);
                M.setflag(GiNaC::status_flags::dynallocated);

                unsigned i;
                for (itr = l.begin(), i = 0; itr != l.end(); ++itr, ++i)
                {
                        unsigned j;
                        if (cols == 1)
                        {
                                M(i, 0) = *itr;
                        }
                        else
                        {
                                for (itc = GiNaC::ex_to<GiNaC::lst>(*itr).begin(), j = 0; itc != GiNaC::ex_to<GiNaC::lst>(*itr).end(); ++itc, ++j)
                                        M(i, j) = *itc;
                        }
                }
                return M;
        }

        GiNaC::lst matrix_to_lst2(const GiNaC::ex& m)
        {
                if (GiNaC::is_a<GiNaC::matrix>(m))
                {
                        GiNaC::matrix A = GiNaC::ex_to<GiNaC::matrix>(m);
                        int cols = A.cols();
                        int rows = A.rows();

                        GiNaC::lst ret;
                        if ( cols == 1 )
                        {
                                for (unsigned int i=0; i<=A.rows()-1; i++)
                                {
                                        ret.append(A(i,0));
                                }
                        }
                        else if ( rows == 1 )
                        {
                                for (unsigned int i=0; i<=A.cols()-1; i++)
                                {
                                        ret.append(A(0,i));
                                }
                        }
                        else
                        {
                                for (unsigned int i=0; i<=A.rows()-1; i++)
                                {
                                        GiNaC::lst rl;
                                        for (unsigned int j=0; j<=A.cols()-1; j++)
                                        {
                                                rl.append(A(i,j));
                                        }
                                        ret.append(rl);
                                }
                        }
                        return ret;
                }
                else
                {
                        return GiNaC::lst();
                }
        }

        GiNaC::lst lst_equals(GiNaC::ex a, GiNaC::ex b)
        {
                GiNaC::lst ret;
                if ( (GiNaC::is_a<GiNaC::lst>(a)) && (GiNaC::is_a<GiNaC::lst>(b)) /*&& (a.nops() == b.nops())*/ ) {
                for (unsigned int i=0; i<= a.nops()-1; i++)
                {
                        ret.append(b.op(i) == a.op(i));
                }
        }
        else if ( !(GiNaC::is_a<GiNaC::lst>(a)) && !(GiNaC::is_a<GiNaC::lst>(b)))
        {
                ret.append(b == a);
        }
        else if ( !(GiNaC::is_a<GiNaC::lst>(a)) && (GiNaC::is_a<GiNaC::lst>(b)))
        {
                ret.append(b.op(0) == a);
        }
        else
        {
                throw(std::invalid_argument("Make sure that the lists a and b are comparable."));
        }
        return ret;
}


int find(GiNaC::ex e, GiNaC::lst list)
{
        for (unsigned int i=0; i< list.nops(); i++)
        {
                if ( e == list.op(i) ) return i;
        }
        return -1;
}


void visitor_subst_pow(GiNaC::ex e, GiNaC::exmap& map, ex_int_map& intmap, string a)
{
        static int i=0;
        if (map.find(e) != map.end())
        {
                intmap[e] = intmap[e]+1;
                return;
        }
        if (GiNaC::is_exactly_a<GiNaC::power>(e))
        {
                std::ostringstream s;
                s <<a<<i++;
                map[e] = GiNaC::symbol(s.str());
                intmap[e] = 0;
                for (unsigned int i=0; i< e.nops(); i++)
                {
                        GiNaC::ex e2 = e.op(i);
                        //       cout <<"power e "<<e2<<endl;
                        visitor_subst_pow(e2,map,intmap, a);
                }
        }
        else if (GiNaC::is_a<GiNaC::function>(e))
        {
                std::ostringstream s;
                s <<a<<i++;
                map[e] = GiNaC::symbol(s.str());
                intmap[e] = 0;
                for (unsigned int i=0; i< e.nops(); i++)
                {
                        GiNaC::ex e2 = e.op(i);
                        //       cout <<"function e "<<e2<<endl;
                        visitor_subst_pow(e2,map,intmap, a);
                }
        }
        else if (GiNaC::is_a<GiNaC::mul>(e))
        {
                if (e.nops() > 4 && e.nops() < 10 )
                {
                        std::ostringstream s;
                        s <<a<<i++;
                        map[e] = GiNaC::symbol(s.str());
                        intmap[e] = 0;
                }

                for (unsigned int i=0; i< e.nops(); i++)
                {
                        GiNaC::ex e2 = e.op(i);
                        visitor_subst_pow(e2,map,intmap, a);
                }
        }
        else if (GiNaC::is_a<GiNaC::add>(e))
        {
                for (unsigned int i=0; i< e.nops(); i++)
                {
                        GiNaC::ex e2 = e.op(i);
                        visitor_subst_pow(e2,map,intmap,a);
                }
        }
}


void check_visitor(GiNaC::ex e, GiNaC::lst& exlist)
{
        if (find(e, exlist) >= 0) return;

        //  cout <<"ex e "<<e<<endl;
        if (GiNaC::is_a<GiNaC::numeric>(e))
        {
        }
        else if (GiNaC::is_a<GiNaC::add>(e) )
        {
                //    cout <<"e "<<e <<endl;
                //    cout <<"e.nops() "<<e.nops() <<endl;
                if (e.nops() > 4 && e.nops() < 10 ) exlist.append(e);
                for (unsigned int i=0; i< e.nops(); i++)
                {
                        GiNaC::ex e2 = e.op(i);
                        //       cout <<"add e "<<e2<<endl;
                        //       exlist.append(e2);
                        check_visitor(e2,exlist);
                }
        }
        else if (GiNaC::is_a<GiNaC::mul>(e))
        {
                for (unsigned int i=0; i< e.nops(); i++)
                {
                        GiNaC::ex e2 = e.op(i);
                        //       cout <<"mul e "<<e2<<endl;
                        exlist.append(e2);
                        check_visitor(e2,exlist);
                }
        }
        else if (GiNaC::is_a<GiNaC::lst>(e))
        {
                for (unsigned int i=0; i< e.nops(); i++)
                {
                        GiNaC::ex e2 = e.op(i);
                        //       cout <<"GiNaC::lst e "<<e2<<endl;
                        //       exlist.append(e2);
                        check_visitor(e2,exlist);
                }
        }
        else if (GiNaC::is_exactly_a<GiNaC::power>(e))
        {
                exlist.append(e);
                for (unsigned int i=0; i< e.nops(); i++)
                {
                        GiNaC::ex e2 = e.op(i);
                        //       cout <<"power e "<<e2<<endl;
                        check_visitor(e2,exlist);
                }
        }
        else if (GiNaC::is_a<GiNaC::function>(e))
        {
                exlist.append(e);
                for (unsigned int i=0; i< e.nops(); i++)
                {
                        GiNaC::ex e2 = e.op(i);
                        //       cout <<"function e "<<e2<<endl;
                        check_visitor(e2,exlist);
                }
        }

        else
        {
                //       exlist.append(e);
                //    cout <<"atom e "<<e<<endl;
        }

        exlist.sort();
        exlist.unique();
}


GiNaC::ex homogenous_pol(unsigned int order, unsigned int nsd, const string a)
{
        using SyFi::x;
        using SyFi::y;
        using SyFi::z;

        if ( nsd == 1)
        {
                GiNaC::symbol a0(istr(a,0));
                return GiNaC::lst(a0*pow(x,order), a0, pow(x,order));
        }
        else if ( nsd == 2 )
        {
                GiNaC::ex variables = get_symbolic_matrix(1,order+1, a);
                GiNaC::lst basis;
                GiNaC::ex ret;
                for (unsigned int i=0; i<= order; i++)
                {
                        basis.append(pow(x,i)*pow(y,order-i));
                        ret += variables.op(i)*basis.op(i);
                }
                return GiNaC::lst(ret, matrix_to_lst2(variables), basis);
        }
        else if ( nsd == 3 )
        {
                GiNaC::lst basis;
                for (unsigned int i=0; i<= order; i++)
                {
                        for (unsigned int j=0; j<= order; j++)
                        {
                                for (unsigned int k=0; k<= order; k++)
                                {
                                        if ( i + j + k == order )
                                        {
                                                basis.append(pow(x,i)*pow(y,j)*pow(z,k));
                                        }
                                }
                        }
                }
                GiNaC::ex variables = get_symbolic_matrix(1,basis.nops(), a);
                GiNaC::ex ret;
                for (unsigned int i=0; i<basis.nops(); i++)
                {
                        ret += variables.op(i)*basis.op(i);
                }
                return GiNaC::lst(ret, matrix_to_lst2(variables), basis);
        }
        throw std::runtime_error("Homogenous polynomials only implemented in 1D, 2D and 3D");
}


GiNaC::lst homogenous_polv(unsigned int no_fields, unsigned int order, unsigned int nsd, const string a)
{
        GiNaC::lst ret1;                         // contains the polynom
        GiNaC::lst ret2;                         // contains the coefficients
        GiNaC::lst ret3;                         // constains the basis functions
        GiNaC::lst basis_tmp;
        for (unsigned int i=0; i< no_fields; i++)
        {
                GiNaC::lst basis;
                std::ostringstream s;
                s <<a<<""<<i<<"_";
                GiNaC::ex polspace = homogenous_pol(order, nsd, s.str());
                ret1.append(polspace.op(0));
                ret2.append(polspace.op(1));
                basis_tmp = GiNaC::ex_to<GiNaC::lst>(polspace.op(2));
                for (GiNaC::lst::const_iterator basis_iterator = basis_tmp.begin();
                        basis_iterator != basis_tmp.end(); ++basis_iterator)
                {
                        GiNaC::lst tmp_lst;
                        for (unsigned int d=1; d<=no_fields; d++) tmp_lst.append(0);
                        tmp_lst.let_op(i) = (*basis_iterator);
                        ret3.append(tmp_lst);
                }
        }
        return GiNaC::lst(ret1,ret2,ret3);
}


GiNaC::ex pol(unsigned int order, unsigned int nsd, const string a)
{
        using SyFi::x;
        using SyFi::y;
        using SyFi::z;

        GiNaC::ex ret;                           // ex to return
        int dof;                                         // degrees of freedom
        GiNaC::ex A;                             // ex holding the coefficients a_0 .. a_dof
        GiNaC::lst basis;

        if (nsd == 1)
        {
                /* 1D:
                 * P^n = a_0 + a_1*x + .... + a_n*x^n
                 * dof : n+1
                 */
                dof = order+1;
                A = get_symbolic_matrix(1,dof, a);
                int o=0;
                for (GiNaC::const_iterator i = A.begin(); i != A.end(); ++i)
                {
                        ret += (*i)*pow(x,o);
                        basis.append(pow(x,o));
                        o++;
                }
        }
        else if ( nsd == 2)
        {

                /* 2D: structure of coefficients (a_i)
                 * [ a_0      a_1 x     a_3 x^2     a_6 x^3
                 * [ a_2 y    a_4 xy    a_7 x^2y
                 * [ a_5 y^2  a_8 xy^2
                 * [ a_9 y^3
                 */
                dof = (order+1)*(order+2)/2;
                A = get_symbolic_matrix(1, dof , a);

                size_t i=0;
                for (unsigned int o = 0; o <= order; o++)
                {
                        for (unsigned int d = 0; d <= o; d++)
                        {
                                ret += A.op(i)*pow(y,d)*pow(x,o-d);
                                basis.append(pow(y,d)*pow(x,o-d));
                                i++;
                        }
                }
        }
        else if (nsd == 3)
        {

                /* Similar structure as in 2D, but
                 * structured as a tetraheder, i.e.,
                 *   a_o + a_1 x + a_2 y + a_3 z
                 * + a_4 x^2 + a_5 xy +
                 */
                dof = 0;
                for (unsigned int j=0; j<= order; j++)
                {
                        dof += (j+1)*(j+2)/2;
                }
                A = get_symbolic_matrix(1, dof , a);

                size_t i=0;
                for (unsigned int o = 0; o <= order; o++)
                {
                        for (unsigned int d = 0; d <= o; d++)
                        {
                                for (unsigned int f = 0; f <= o; f++)
                                {
                                        if ( int(o)-int(d)-int(f) >= 0)
                                        {
                                                ret += A.op(i)*pow(y,f)*pow(z,d)*pow(x,o-d-f);
                                                basis.append(pow(y,f)*pow(z,d)*pow(x,o-d-f));
                                                i++;
                                        }
                                }
                        }
                }
        }
        return GiNaC::lst(ret,matrix_to_lst2(A), basis);
}


GiNaC::lst polv(unsigned int no_fields, unsigned int order, unsigned int nsd, const string a)
{
        GiNaC::lst ret1;                         // contains the polynom
        GiNaC::lst ret2;                         // contains the coefficients
        GiNaC::lst ret3;                         // constains the basis functions
        GiNaC::lst basis_tmp;
        for (unsigned int i=0; i< no_fields; i++)
        {
                GiNaC::lst basis;
                std::ostringstream s;
                s <<a<<""<<i<<"_";
                GiNaC::ex polspace = pol(order, nsd, s.str());
                ret1.append(polspace.op(0));
                ret2.append(polspace.op(1));
                basis_tmp = GiNaC::ex_to<GiNaC::lst>(polspace.op(2));
                for (GiNaC::lst::const_iterator basis_iterator = basis_tmp.begin();
                        basis_iterator != basis_tmp.end(); ++basis_iterator)
                {
                        GiNaC::lst tmp_lst;
                        for (unsigned int d=1; d<=no_fields; d++) tmp_lst.append(0);
                        tmp_lst.let_op(i) = (*basis_iterator);
                        ret3.append(tmp_lst);
                }
        }
        return GiNaC::lst(ret1,ret2,ret3);

        /* Old Code:
           GiNaC::lst ret;
           for (int i=1; i<= nsd; i++) {
           std::ostringstream s;
           s <<a<<"^"<<i<<"_";
           GiNaC::ex p = pol(order, nsd, s.str());
        ret.append(p);
        }
        return ret;
        */
}


GiNaC::ex polb(unsigned int order, unsigned int nsd, const string a)
{
        using SyFi::x;
        using SyFi::y;
        using SyFi::z;

        GiNaC::ex ret;                           // ex to return
        int dof;                                         // degrees of freedom
        GiNaC::ex A;                             // ex holding the coefficients a_0 .. a_dof
        GiNaC::lst basis;

        if (nsd == 1)
        {
                /* 1D:
                 * P^n = a_0 + a_1*x + .... + a_n*x^n
                 * dof : n+1
                 */
                dof = order+1;
                A = get_symbolic_matrix(1,dof, a);
                int o=0;
                for (GiNaC::const_iterator i = A.begin(); i != A.end(); ++i)
                {
                        ret += (*i)*pow(x,o);
                        basis.append(pow(x,o));
                        o++;
                }
        }
        else if ( nsd == 2)
        {

                /* 2D: structure of coefficients (a_i)
                 * [ a_0      a_1 x     a_3 x^2     a_6 x^3
                 * [ a_2 y    a_4 xy    a_7 x^2y
                 * [ a_5 y^2  a_8 xy^2
                 * [ a_9 y^3
                 */

                dof = (order+1)*(order+1);
                A = get_symbolic_matrix(1, dof , a);

                size_t i=0;
                for (unsigned int o = 0; o <= order; o++)
                {
                        for (unsigned int d = 0; d <= order; d++)
                        {
                                ret += A.op(i)*pow(y,d)*pow(x,o);
                                basis.append(pow(y,d)*pow(x,o));
                                i++;
                        }
                }
        }
        else if (nsd == 3)
        {

                /* Similar structure as in 2D, but
                 * structured as a tetraheder, i.e.,
                 *   a_o + a_1 x + a_2 y + a_3 z
                 * + a_4 x^2 + a_5 xy +
                 */
                dof = (order+1)*(order+1)*(order+1);
                A = get_symbolic_matrix(1, dof , a);

                size_t i=0;
                for (unsigned int o = 0; o <= order; o++)
                {
                        for (unsigned int d = 0; d <= order; d++)
                        {
                                for (unsigned int f = 0; f <= order; f++)
                                {
                                        ret += A.op(i)*pow(y,f)*pow(z,d)*pow(x,o);
                                        basis.append(pow(y,f)*pow(z,d)*pow(x,o));
                                        i++;
                                }
                        }
                }
        }

        return GiNaC::lst(ret,matrix_to_lst2(A), basis);
}


GiNaC::lst coeffs(GiNaC::lst pols)
{
        GiNaC::lst cc;
        GiNaC::lst tmp;
        for (unsigned int i=0; i<= pols.nops()-1; i++)
        {
                tmp = coeffs(pols.op(i));
                cc = collapse(GiNaC::lst(cc, tmp));
        }
        return cc;
}


GiNaC::lst coeffs(GiNaC::ex pol)
{
        using SyFi::x;
        using SyFi::y;
        using SyFi::z;

        GiNaC::lst cc;
        GiNaC::ex c, b;
        for (int i=pol.ldegree(x); i<=pol.degree(x); ++i)
        {
                for (int j=pol.ldegree(y); j<=pol.degree(y); ++j)
                {
                        for (int k=pol.ldegree(z); k<=pol.degree(z); ++k)
                        {
                                c = pol.coeff(x,i).coeff(y, j).coeff(z,k);
                                if ( c != 0 ) cc.append(c);
                        }
                }
        }
        return cc;
}


GiNaC::exvector coeff(GiNaC::ex pol)
{
        using SyFi::x;
        using SyFi::y;
        using SyFi::z;

        GiNaC::exvector cc;
        GiNaC::ex c, b;
        for (int i=pol.ldegree(x); i<=pol.degree(x); ++i)
        {
                for (int j=pol.ldegree(y); j<=pol.degree(y); ++j)
                {
                        for (int k=pol.ldegree(z); k<=pol.degree(z); ++k)
                        {
                                c = pol.coeff(x,i).coeff(y, j).coeff(z,k);
                                if ( c != 0 ) cc.insert(cc.begin(),c);
                        }
                }
        }
        return cc;
}


GiNaC::exmap pol2basisandcoeff(GiNaC::ex e, GiNaC::ex s)
{
        if (GiNaC::is_a<GiNaC::symbol>(s))
        {
                GiNaC::symbol ss = GiNaC::ex_to<GiNaC::symbol>(s);
                e = expand(e);
                GiNaC::ex c;
                GiNaC::ex b;
                GiNaC::exmap map;
                for (int i=e.ldegree(ss); i<=e.degree(ss); ++i)
                {
                        c = e.coeff(ss,i);
                        b = pow(ss,i);
                        map[b] = c;
                }
                return map;
        }
        else
        {
                throw(std::invalid_argument("The second argument must be a symbol."));
        }
}


GiNaC::exmap pol2basisandcoeff(GiNaC::ex e)
{
        using SyFi::x;
        using SyFi::y;
        using SyFi::z;

        e = expand(e);
        GiNaC::ex c;
        GiNaC::ex b;
        GiNaC::exmap map;
        for (int i=e.ldegree(x); i<=e.degree(x); ++i)
        {
                for (int j=e.ldegree(y); j<=e.degree(y); ++j)
                {
                        for (int k=e.ldegree(z); k<=e.degree(z); ++k)
                        {
                                c = e.coeff(x,i).coeff(y, j).coeff(z,k);
                                b = pow(x,i)*pow(y,j)*pow(z,k);
                                map[b] = c;
                        }
                }
        }
        return map;
}


GiNaC::ex legendre1D(const GiNaC::symbol x, unsigned int n)
{
        GiNaC::ex P;
        // Rodrigue's formula for Legendre polynomial of 1D
        // The interval [-1, 1]
        P=1/(pow(2,n)*GiNaC::factorial(n))*GiNaC::diff(GiNaC::pow((x*x-1),n),x,n);
        // -----------------
        // The interval [0,1]
        //  GiNaC::ex xx = 2*x - 1;
        // P=1/(pow(2,2*n)*GiNaC::factorial(n))*GiNaC::diff(GiNaC::pow((xx*xx-1),n),x,n);
        return P;
}


GiNaC::ex legendre(unsigned int order, unsigned int nsd, const string s)
{
        using SyFi::x;
        using SyFi::y;
        using SyFi::z;

        // The Legendre polynomials to be used in FiniteElement
        GiNaC::ex leg;
        GiNaC::ex A;
        GiNaC::lst basis;
        int dof;

        GiNaC::ex b;

        // 1D
        if(nsd == 1)
        {
                dof = order+1;
                A = get_symbolic_matrix(1,dof,s);
                int o=0;
                for(GiNaC::const_iterator i = A.begin(); i!=A.end(); ++i)
                {
                        b= legendre1D(x,o);
                        leg+= (*i)*b;
                        basis.append(b);
                        o++;
                }
        }
        // 2D
        /*
        else if(nsd == 2){  // NB: Only for tensor products on TRIANGLES (not boxes)
                / * 2D: structure of coefficients (a_i)
                 * [ a_0           a_1 P_1(x)           a_3 P_2(x)        a_6 P_3(x)
                 * [ a_2 P_1(y)    a_4 P_1(x)*P_1(y)    a_7 P_2(x)*P_1(y)
                 * [ a_5 P_2(y)    a_8 P_1(x)*P_2(y)
        * [ a_9 P_3(y)
        * /
        dof = (order+1)*(order+2)/2;
        A = get_symbolic_matrix(1,dof,s);
        size_t i=0;
        for (int o = 0; o <= order; o++) {
        for (int d = 0; d <= o; d++) {
        b = legendre1D(y,d)*legendre1D(x,o-d);
        leg += A.op(i)*b;
        basis.append(b);
        i++;

        }
        }
        }
        */
        else if(nsd == 2)                        // NB: Only for tensor products on rectangles
        {
                dof = (order+1)*(order+1);
                A = get_symbolic_matrix(1,dof,s);
                size_t i=0;
                for (unsigned int o = 0; o <= order; o++)
                {
                        for (unsigned int d = 0; d <= order; d++)
                        {
                                b = legendre1D(y,d)*legendre1D(x,o);
                                leg += A.op(i)*b;
                                basis.append(b);
                                i++;

                        }
                }
        }

        /* tetrahedron
        else if(nsd==3){
                dof = 0;
                for (int j=0; j<= order; j++) {
                        dof += (j+1)*(j+2)/2;
                }
        A = get_symbolic_matrix(1, dof , s);

        size_t i=0;
        for (int o = 0; o <= order; o++) {
        for (int d = 0; d <= o; d++) {
        for (int f = 0; f <= o; f++) {
        if ( o-d-f >= 0) {
        b = legendre1D(y,f)*legendre1D(z,d)*legendre1D(x,o-d-f);
        leg += A.op(i)*b;
        basis.append(b);
        i++;
        }
        }
        }
        }
        }
        */

        else if(nsd==3)
        {
                dof = (order+1)*(order+1)*(order+1);
                A = get_symbolic_matrix(1, dof , s);

                size_t i=0;
                for (unsigned int o = 0; o <= order; o++)
                {
                        for (unsigned int d = 0; d <= order; d++)
                        {
                                for (unsigned int f = 0; f <= order; f++)
                                {
                                        b = legendre1D(y,f)*legendre1D(z,d)*legendre1D(x,o);
                                        leg += A.op(i)*b;
                                        basis.append(b);
                                        i++;
                                }
                        }
                }
        }
        return GiNaC::lst(leg,matrix_to_lst2(A), basis);
}


GiNaC::lst legendrev(unsigned int no_fields, unsigned int order, unsigned int nsd, const string a)
{
        GiNaC::lst ret1;                         // contains the polynom
        GiNaC::lst ret2;                         // contains the coefficients
        GiNaC::lst ret3;                         // constains the basis functions
        GiNaC::lst basis_tmp;
        for (unsigned int i=1; i<= no_fields; i++)
        {
                GiNaC::lst basis;
                std::ostringstream s;
                s <<a<<""<<i<<"_";
                GiNaC::ex polspace = legendre(order, nsd, s.str());
                ret1.append(polspace.op(0));
                ret2.append(polspace.op(1));
                basis_tmp = GiNaC::ex_to<GiNaC::lst>(polspace.op(2));
                for (GiNaC::lst::const_iterator basis_iterator = basis_tmp.begin();
                        basis_iterator != basis_tmp.end(); ++basis_iterator)
                {
                        GiNaC::lst tmp_lst;
                        for (unsigned int d=1; d<=no_fields; d++) tmp_lst.append(0);
                        tmp_lst.let_op(i-1) = (*basis_iterator);
                        ret3.append(tmp_lst);
                }
        }
        return GiNaC::lst(ret1,ret2,ret3);
}


bool compare(const ex & e, const string & s)
{
        ostringstream ss;
        ss << e;
        return ss.str() == s;
}


void EQUAL_OR_DIE(const ex & e, const string & s)
{
        if (!compare(e, s))
        {
                ostringstream os;
                os << "ERROR: expression e: " <<e<<" is not equal to "<<s<<endl;
                throw runtime_error(os.str());
        }
}


// =========== expression inspection utilities

class SymbolMapBuilderVisitor:
public visitor,
        public symbol::visitor
{
        public:
                SymbolMapBuilderVisitor(GiNaC::exmap & em): symbolmap(em) {}

                GiNaC::exmap & symbolmap;

                void visit(const symbol & s)
                {
                        GiNaC::exmap::iterator it = symbolmap.find(s);
                        if(it != symbolmap.end())
                        {
                                it->second++;
                        }
                        else
                        {
                                //GiNaC::exmap::value_type p(ex(s), ex(s));
                                //symbolmap.insert(p);
                                symbolmap[ex(s)] = ex(s);
                        }
                }
};

class SymbolCounterVisitor:
public visitor,
        public symbol::visitor,
        public basic::visitor
{
        public:
                exhashmap<int> symbolcount;

                void visit(const basic & s)
                {
                        std::cout << "visiting basic " << std::endl;
                }

                void visit(const symbol & s)
                {
                        ex e = s;
                        std::cout << "visiting symbol " << e << std::endl;
                        exhashmap<int>::iterator it = symbolcount.find(s);
                        if(it != symbolcount.end())
                        {
                                std::cout << "found symbol " << e << std::endl;
                                it->second++;
                        }
                        else
                        {
                                std::cout << "adding symbol " << e << std::endl;
                                pair<ex,int> p;
                                p.first = ex(s);
                                p.second = 1;
                                symbolcount.insert(p);
                        }
                }
};

exhashmap<int> count_symbols(const ex & e)
{
        SymbolCounterVisitor v;
        e.traverse(v);
        return v.symbolcount;
}


/*
GiNaC::exvector get_symbols3(const GiNaC::ex & e)
{
        // Implemented directly to avoid copying map:
        SymbolCounterVisitor v;
        e.traverse(v);
exhashmap<int> & sc = v.symbolcount;

int i = 0;
GiNaC::exvector l(sc.size());
for(exhashmap<int>::iterator it=sc.begin(); it!=sc.end(); it++)
{
//l.push_back(it->first);
l[i] = it->first;
i++;
}
return l;
}

std::list<GiNaC::ex> get_symbols2(const ex & e)
{
// Implemented directly to avoid copying map:
SymbolCounterVisitor v;
e.traverse(v);
exhashmap<int> & sc = v.symbolcount;

list<ex> l;
for(exhashmap<int>::iterator it=sc.begin(); it!=sc.end(); it++)
{
l.push_back(it->first);
}
return l;
}

void get_symbols(const ex & e, GiNaC::exmap & em)
{
SymbolMapBuilderVisitor v(em);
e.traverse(v);
}
*/
ex extract_symbols(const ex & e)
{
        // Implemented directly to avoid copying map:
        SymbolCounterVisitor v;
        e.traverse(v);
        exhashmap<int> & sc = v.symbolcount;

        lst l;
        for(exhashmap<int>::iterator it=sc.begin(); it!=sc.end(); it++)
        {
                l.append(it->first);
                std::cout << (it->first) << std::endl;
        }
        ex ret = l;
        return ret;
}


// Collect all symbols of an expression
void collect_symbols(const GiNaC::ex & e, exset & v)
{
        if (GiNaC::is_a<GiNaC::symbol>(e))
        {
                v.insert(e);
        }
        else
        {
                for (size_t i=0; i<e.nops(); i++)
                {
                        collect_symbols(e.op(i), v);
                }
        }
}


GiNaC::exvector collect_symbols(const GiNaC::ex & e)
{
        exset s;
        collect_symbols(e, s);
        GiNaC::exvector v(s.size());
        for(exset::iterator i=s.begin(); i!= s.end(); i++)
        {
                v.push_back(*i);
        }
        return v;
}


bool compare_archives(const string & first, const string & second, std::ostream & os)
{
        bool ret = true;

        // read both archives
        archive a1, a2;
        ifstream if1(first.c_str()), if2(second.c_str());
        if1 >> a1;
        if2 >> a2;

        // compare size
        int n = a1.num_expressions();
        int n2 = a2.num_expressions();
        if(n != n2)
        {
                os << "Archives " << first << " and " << second
                        << " has a different number of expressions, " << n << " and " << n2 << "." << endl;
                os << "Comparing common expressions." << endl;
                ret = false;
        }

        // iterate over all expressions in first archive
        ex e1,e2;
        for(int i=0; i<n; i++)
        {
                lst syms;
                string exname;

                e1 = a1.unarchive_ex(syms, exname, i);

                syms = ex_to<lst>(extract_symbols(e1));
                //        os << "Comparing " << exname << " with symbols " << syms << endl;

                // is this in the second archive?
                try
                {
                        e2 = a2.unarchive_ex(syms, exname.c_str());

                        // got it, now compare
                        bool isequal = is_zero(e1-e2);
                        if(!isequal)
                        {
                                if(ret)
                                {
                                        os << "Archives " << first << " and " << second
                                                << " are not equal, details follow:" << endl;
                                }
                                os << "Expression with name " << exname << " is not equal:" << endl;
                                os << "First:  " << endl << e1 << endl;
                                os << "Second: " << endl << e2 << endl;
                                ret = false;
                        }
                }
                catch(...)
                {
                        os << "Expression " << exname << " is missing from " << second << "." << endl;
                        ret = false;
                }
        }

        return ret;
}


class ExStatsVisitor:
public visitor,
        public power::visitor,
        public mul::visitor,
        public add::visitor,
        public function::visitor
{
        public:
                ExStats es;

                void visit(const power & e)
                {
                        //cout << "pow" << endl;
                        ex b = e.op(1);
                        ex c = b.evalf();
                        if( is_a<numeric>(c) )
                        {
                                numeric nu = ex_to<numeric>(c);
                                int n = (int) to_double(nu);
                                if( is_zero(nu-n) )
                                {
                                        // is an integer power, interpret as a sequence of multiplications
                                        es.muls  += n-1;
                                        es.flops += n-1;
                                }
                        }

                        es.pows++;
                }

                void visit(const mul & e)
                {
                        //cout << "mul" << endl;
                        es.muls  += e.nops()-1;
                        es.flops += e.nops()-1;
                }

                void visit(const add & e)
                {
                        //cout << "add" << endl;
                        es.adds  += e.nops()-1;
                        es.flops += e.nops()-1;
                }

                void visit(const function & e)
                {
                        //cout << "func" << endl;
                        es.functions++;
                }
};

ExStats count_ops(const ex & e)
{
        //cout << "count_ops " << e << endl;
        //cout << "is an add: " << GiNaC::is<GiNaC::add>(e) << endl;
        //cout << "is a  mul: " << GiNaC::is<GiNaC::mul>(e) << endl;
        ExStatsVisitor v;
        e.traverse(v);
        return v.es;
}


// =========== expression manipulation utilities

// Returns in sel a list with symbol/expression pairs for all
// power replacement variables introduced. Works best on
// expressions in expanded form.
ex replace_powers(const ex & ein, const list<symbol> & symbols, list<symexpair> & sel, const string & tmpsymbolprefix)
{
        ex e = ein;
        // build power expressions
        list<symbol>::const_iterator it = symbols.begin();
        for(; it != symbols.end(); it++)
        {
                int deg      = e.degree(*it);
                if(deg > 0)
                {
                        symbol sym   = ex_to<symbol>(*it);
                        string sname = tmpsymbolprefix + sym.get_name();

                        // make list of new symbols
                        vector<symbol> symbols(deg);
                        symbols[0] = sym;
                        for(int i=1; i<deg; i++)
                        {
                                symbols[i] = get_symbol( sname + int2string(i+1) );
                                sel.push_back(make_pair(symbols[i], symbols[i-1]*sym));
                        }

                        // with highest order first, subs in e
                        ex prod = sym;
                        for(int i=deg-1; i>=1; i--)
                        {
                                e = e.subs(power(sym,i+1) == symbols[i], subs_options::algebraic);
                        }
                }
        }
        return e;
}


};                                                               // namespace SyFi