time_extrapolation.h

/*
 * SPDX-FileCopyrightText: Copyright (c) 2014-2020:  G-CSC, Goethe University Frankfurt
 * SPDX-License-Identifier: LicenseRef-UG4-LGPL-3.0
 *
 * Author: Arne Naegel
 *
 * This file is part of UG4.
 *
 * UG4 is free software: you can redistribute it and/or modify it under the
 * terms of the GNU Lesser General Public License version 3 (as published by the
 * Free Software Foundation) with the following additional attribution
 * requirements (according to LGPL/GPL v3 §7):
 *
 * (1) The following notice must be displayed in the Appropriate Legal Notices
 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
 *
 * (2) The following notice must be displayed at a prominent place in the
 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
 *
 * (3) The following bibliography is recommended for citation and must be
 * preserved in all covered files:
 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
 *   parallel geometric multigrid solver on hierarchically distributed grids.
 *   Computing and visualization in science 16, 4 (2013), 151-164"
 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
 *   flexible software system for simulating pde based models on high performance
 *   computers. Computing and visualization in science 16, 4 (2013), 165-179"
 *
 * This program 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 Lesser General Public License for more details.
 */
 
#ifndef TIME_EXTRAPOLATION_H_
#define TIME_EXTRAPOLATION_H_
 
#include <vector>
#include <cmath>
 
// ug libraries
#include "common/common.h"
#include "common/util/smart_pointer.h"
 
#include "lib_algebra/lib_algebra.h"
#include "lib_algebra/operator/debug_writer.h"
 
#include "lib_disc/function_spaces/grid_function.h"
#include "lib_disc/function_spaces/grid_function_user_data.h"
#include "lib_disc/function_spaces/integrate.h"
#include "lib_disc/function_spaces/metric_spaces.h"
 
#include "lib_disc/time_disc/time_disc_interface.h"
#include "lib_disc/spatial_disc/user_data/linker/scale_add_linker.h"
 
 
 
namespace ug{
 
/*
 std::vector<size_t> steps {1, 2, 4};
 timex = new AitkenNevilleTimex(steps);
 
 //
 timex.set_solution(sol0, 0);  // single step
 timex.set_solution(sol1, 1);  // double step
 timex.set_solution(sol2, 2);  // four step
 
 timex.set_estimator(L2NormEstimator())
 timex.set_estimator(InfNormEstimator())
 timex.set_estimator(RelativeEstimator())
 //
 timex.apply()                // updating sol1 and sol2
 
 timex.get_error_estimate()
 
 */
 
namespace tools {
 

    inline void VecScaleAddWithNormRel(double &vUpdate, double alpha2, const double &vFine, double alpha3, const double &vCoarse, double &norm)
    {
        const double update = alpha2*vFine + alpha3*vCoarse;
        vUpdate = vUpdate + update;
        norm = std::max(norm, 0.5*fabs(update)/(1.0+fabs(vFine)+fabs(vCoarse)));
    }

    template<typename vector_t, template <class T> class TE_VEC>
    inline void VecScaleAddWithNormRel(TE_VEC<vector_t> &vUpdate, double alpha2, const TE_VEC<vector_t> &vFine, double alpha3, const TE_VEC<vector_t> &vCoarse, double &norm)
    {
        for(size_t i=0; i<vUpdate.size(); i++)
            VecScaleAddWithNormRel(vUpdate[i], alpha2, vFine[i], alpha3, vCoarse[i], norm);
    }

    inline void VecScaleAddWithNormInf(double &vUpdate, double alpha2, const double &vFine, double alpha3, const double &vCoarse, double &norm)
    {
        const double update = alpha2*vFine + alpha3*vCoarse;
        vUpdate = vUpdate + update;
        norm = std::max(norm, fabs(update));
    }

    template<typename vector_t, template <class T> class TE_VEC>
    inline void VecScaleAddWithNormInf(TE_VEC<vector_t> &vUpdate, double alpha2, const TE_VEC<vector_t> &vFine, double alpha3, const TE_VEC<vector_t> &vCoarse, double &norm, const int delta=1, const int offset=0)
    {
        // std::cerr << norm << " "<< delta << offset << std::endl;
        for(size_t i=offset; i<vUpdate.size(); i+=delta)
            VecScaleAddWithNormInf(vUpdate[i], alpha2, vFine[i], alpha3, vCoarse[i], norm);
 
        // std::cerr << norm << std::endl;
    }
 

    inline void VecScaleAddWithNorm2(double &vUpdate, double alpha2, const double &vFine, double alpha3, const double &vCoarse, double &norm)
    {
        const double update = alpha2*vFine + alpha3*vCoarse;
        vUpdate = vUpdate+ update;
        norm += update*update;
    }

    template<typename vector_t, template <class T> class TE_VEC>
    inline void VecScaleAddWithNorm2(TE_VEC<vector_t> &vUpdate, double alpha2, const TE_VEC<vector_t> &vFine, double alpha3, const TE_VEC<vector_t> &vCoarse, double &norm, const int delta=1, const int offset=0)
    {
        for(size_t i=offset; i<vUpdate.size(); i+=delta)
            VecScaleAddWithNorm2(vUpdate[i], alpha2, vFine[i], alpha3, vCoarse[i], norm);
    }
 
 
 
}


template <class TVector>
class ISubDiagErrorEst
{
public:
    // constructor
    ISubDiagErrorEst() : m_est(0.0) {};
 
    // destructor
    virtual ~ISubDiagErrorEst() {};

    virtual bool update(SmartPtr<TVector> vUpdate, number alpha2, SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse) = 0;

    number get_current_estimate() {return m_est; };
    void reset_estimate() {  m_est=0.0; };
 
    virtual std::string config_string() const {return "ISubDiagErrorEst";}
 
protected:
    number m_est;
};
 


template <class TVector>
class NormInfEstimator : public ISubDiagErrorEst<TVector>
{
protected:
    typedef ISubDiagErrorEst<TVector> base_type;
    int m_stride;
    int m_offset;
 
public:
    // constructor
    NormInfEstimator() :
    ISubDiagErrorEst<TVector>(), m_stride(1), m_offset(0) {};
 
 
    // apply
    bool update(SmartPtr<TVector> vUpdate, number alpha2, SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse)
    {
        const int delta = m_stride;
        const int offset = m_offset;
        base_type::m_est=0.0;
        tools::VecScaleAddWithNormInf(*vUpdate, alpha2, *vFine, -alpha2, *vCoarse, base_type::m_est, delta, offset);
        return true;
    }
 
    // offset (e.g., component to work on for systems w/ CPU1)
    void set_offset(int offset) {m_offset=offset;}
 
    // delta (e.g., total number components for systems w/ CPU1)
    void set_stride(int delta) {m_stride=delta;}
 
};
 

template <class TVector>
class Norm2Estimator : public ISubDiagErrorEst<TVector>
{
protected:
    typedef ISubDiagErrorEst<TVector> base_type;
    int m_stride;
    int m_offset;
 
public:
    // constructor
    Norm2Estimator() :
        ISubDiagErrorEst<TVector>(), m_stride(1), m_offset(0) {};
    Norm2Estimator(int stride) :
        ISubDiagErrorEst<TVector>(), m_stride(stride), m_offset(0) {};
    Norm2Estimator(int delta, int offset) :
        ISubDiagErrorEst<TVector>(), m_stride(delta), m_offset (offset)  {};
 
    // apply
    bool update(SmartPtr<TVector> vUpdate,  number alpha2,  SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse)
    {
        const int delta = m_stride;
        const int offset = m_offset;
        base_type::m_est=0.0;
        tools::VecScaleAddWithNorm2(*vUpdate, alpha2, *vFine, -alpha2, *vCoarse, base_type::m_est, delta, offset);
#ifdef UG_PARALLEL
        double locEst = base_type::m_est;
        double globEst =0.0;
        vUpdate->layouts()->proc_comm().allreduce(&locEst, &globEst, 1, PCL_DT_DOUBLE, PCL_RO_SUM);
        base_type::m_est = globEst;
#endif
        base_type::m_est = sqrt(base_type::m_est);
        return true;
    }
 
    // offset (e.g., component to work on for systems w/ CPU1)
    void set_offset(int offset) { m_offset=offset;}
 
    // delta (e.g., total number components for systems w/ CPU1)
    void set_stride(int delta) { m_stride=delta;}
 
};
 

template <class TVector>
class NormRelEstimator : public ISubDiagErrorEst<TVector>
{
protected:
    typedef ISubDiagErrorEst<TVector> base_type;
 
public:
    // constructor
    NormRelEstimator() : ISubDiagErrorEst<TVector>() {};
 
    // apply w/ rel norm
    bool update(SmartPtr<TVector> vUpdate, number alpha2,  SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse)
    {
        base_type::m_est=0;
        tools::VecScaleAddWithNormRel(*vUpdate, alpha2, *vFine, -alpha2, *vCoarse, base_type::m_est);
        return true;
    }
};
 
/*
template<class TVector>
class VectorDebugWritingEstimator
        : public VectorDebugWritingObject<TVector>
{
public:
    typedef TVector vector_type;
 
public:
    VectorDebugWritingEstimator()
    : VectorDebugWritingObject<vector_type>() {}
 
    VectorDebugWritingEstimator(SmartPtr<IVectorDebugWriter<vector_type> > spDebugWriter)
    : VectorDebugWritingObject<vector_type>(spDebugWriter) {}
 
    int get_call_id() { return m_dgbCall; }
    void inc_call_id() { m_dgbCall++; }
 
protected:
    int m_dgbCall;
};
*/
 
 
 

template <typename TGridFunction>
class SupErrorEvaluator
: public IComponentSpace<TGridFunction>
{
    public:
        typedef IComponentSpace<TGridFunction> base_type;
 
        SupErrorEvaluator(const char *fctNames) : base_type(fctNames) {};
        //SupErrorEvaluator(const char *fctNames, number scale) : base_type(fctNames, 1, scale) {};
        SupErrorEvaluator(const char *fctNames, const char* ssNames /*, number scale*/)
            : base_type(fctNames, ssNames, 1/*, scale*/) {};
        ~SupErrorEvaluator() {};
 
 
        using IComponentSpace<TGridFunction>::norm;
        using IComponentSpace<TGridFunction>::distance;
 
        double norm(SmartPtr<TGridFunction> uFine) 
        { return norm(*uFine); }
 
        double norm(TGridFunction& uFine)
        {
            // gather subsets in group
            SubsetGroup ssGrp(uFine.domain()->subset_handler());
            if (base_type::m_ssNames != NULL)
                ssGrp.add(TokenizeString(base_type::m_ssNames));
            else
                ssGrp.add_all();
 
            double maxVal = 0.0;
 
            // loop subsets
            for (size_t i = 0; i < ssGrp.size(); ++i)
            {
                // get subset index
                const int si = ssGrp[i];
 
                // loop elements of subset and dim
                switch (ssGrp.dim(i))
                {
                    case DIM_SUBSET_EMPTY_GRID: break;
                    case 0: maxVal = std::max(maxVal, findFctMaxOnSubset<Vertex>(uFine, si)); break;
                    case 1: maxVal = std::max(maxVal, findFctMaxOnSubset<Edge>(uFine, si)); break;
                    case 2: maxVal = std::max(maxVal, findFctMaxOnSubset<Face>(uFine, si)); break;
                    case 3: maxVal = std::max(maxVal, findFctMaxOnSubset<Volume>(uFine, si)); break;
                    default: UG_THROW("SupErrorEvaluator::norm: Dimension " << ssGrp.dim(i) << " not supported.");
                }
            }
 
            #ifdef UG_PARALLEL
                // max over processes
                if (pcl::NumProcs() > 1)
                {
                    pcl::ProcessCommunicator com;
                    number local = maxVal;
                    com.allreduce(&local, &maxVal, 1, PCL_DT_DOUBLE, PCL_RO_MAX);
                }
            #endif
 
            // return the result
            return maxVal;
        }
 
        double norm2(TGridFunction& uFine) 
        { double norm1 = norm(uFine); return norm1*norm1;}
 
        double distance(TGridFunction& uFine, TGridFunction& uCoarse) 
        {
            UG_COND_THROW(uFine.dof_distribution().get() != uCoarse.dof_distribution().get(),
                "Coarse and fine solutions do not have the same underlying dof distro.");
 
            SmartPtr<TGridFunction> uErr = uCoarse.clone();
            uErr->operator-=(uFine);
            return norm(*uErr);
        }
 
        double distance2(TGridFunction& uFine, TGridFunction& uCoarse)
        { double dist = distance(uFine, uCoarse); return dist*dist;}
 
    protected:
        template <typename TBaseElem>
        number findFctMaxOnSubset(const TGridFunction& u, int si) const
        {
            ConstSmartPtr<DoFDistribution> dd = u.dof_distribution();
 
            size_t fct;
            try {fct = dd->fct_id_by_name(base_type::m_fctNames.c_str());}
            UG_CATCH_THROW("Function index could not be determined.\n"
                "Bear in mind that only one function can be evaluated in this error evaluator.");
 
            number maxVal = 0.0;
            typename DoFDistribution::traits<TBaseElem>::const_iterator it = dd->begin<TBaseElem>(si);
            typename DoFDistribution::traits<TBaseElem>::const_iterator itEnd = dd->end<TBaseElem>(si);
            for (; it != itEnd; ++it)
            {
                std::vector<DoFIndex> vInd;
 
                // we compare against all indices on the element
                // which means most indices will be compared against quite often
                // but as this is a sup norm, this is not a problem (only in terms of performance)
                size_t nInd = dd->dof_indices(*it, fct, vInd, false, false);
                for (size_t i = 0; i < nInd; ++i)
                    maxVal = std::max(maxVal, fabs(DoFRef(u, vInd[i])));
            }
 
            return maxVal;
        }
};
 
 

template <typename TDataIn, typename TGridFunction>
class DeltaSquareIntegrand
    : public StdIntegrand<number, TGridFunction::dim, DeltaSquareIntegrand<TDataIn, TGridFunction> >
{
    public:
    //  world dimension of grid function
        static const int worldDim = TGridFunction::dim;
 
    //  data type
        typedef TDataIn data_type;
 
    private:
 
        static number product(const number &x, const number &y)
        {return x*y;}
 
        static number product(const MathVector<worldDim> &x, const MathVector<worldDim> &y)
        {return VecDot(x,y);}
 
 
    //  data to integrate
        SmartPtr<UserData<TDataIn, worldDim> > m_spData;
 
    //  grid function
        const TGridFunction* m_pGridFct1;
        const TGridFunction* m_pGridFct2;
 
    //  time
        number m_time;
 
    public:
        DeltaSquareIntegrand(SmartPtr<UserData<TDataIn, worldDim> > spData,
                        const TGridFunction* pGridFct1, const TGridFunction* pGridFct2,
                        number time)
        : m_spData(spData), m_pGridFct1(pGridFct1), m_pGridFct2(pGridFct2), m_time(time)
        {
            m_spData->set_function_pattern(pGridFct1->function_pattern());
        };

        DeltaSquareIntegrand(SmartPtr<UserData<TDataIn, worldDim> > spData, number time)
        : m_spData(spData), m_pGridFct1(NULL), m_pGridFct2(NULL), m_time(time)
        {
            if(m_spData->requires_grid_fct())
                UG_THROW("UserDataDeltaIntegrand: Missing GridFunction, but "
                        " data requires grid function.")
        };
 
 
 
 
        template <int elemDim>
        void get_values(TDataIn vValue[],
                      ConstSmartPtr<UserData<TDataIn, worldDim> > spData,
                      const TGridFunction& gridFct,
                      const MathVector<worldDim> vGlobIP[],
                      GridObject* pElem,
                      const MathVector<worldDim> vCornerCoords[],
                      const MathVector<elemDim> vLocIP[],
                      const MathMatrix<elemDim, worldDim> vJT[],
                      const size_t numIP)
        {
            //  collect local solution, if required
            if(spData->requires_grid_fct())
            {
            //  create storage
                LocalIndices ind;
                LocalVector u;
 
            //  get global indices
                gridFct.indices(pElem, ind);
 
            //  adapt local algebra
                u.resize(ind);
 
            //  read local values of u
                GetLocalVector(u, gridFct);
                std::cout << u << std::endl;
 
            //  compute data
                try{
                    (*spData)(vValue, vGlobIP, m_time, this->m_si, pElem,
                                vCornerCoords, vLocIP, numIP, &u, &vJT[0]);
                }
                UG_CATCH_THROW("UserDataDeltaIntegrand: Cannot evaluate data.");
            }
            else
            {
            //  compute data
                try{
                    (*spData)(vValue, vGlobIP, m_time, this->m_si, numIP);
                }
                UG_CATCH_THROW("UserDataDeltaIntegrand: Cannot evaluate data.");
            }
 
        };
 

            template <int elemDim>
            void evaluate(number vValue[],
                          const MathVector<worldDim> vGlobIP[],
                          GridObject* pElem,
                          const MathVector<worldDim> vCornerCoords[],
                          const MathVector<elemDim> vLocIP[],
                          const MathMatrix<elemDim, worldDim> vJT[],
                          const size_t numIP)
            {
                std::vector<TDataIn> v1(numIP);
 
                get_values<elemDim>(&v1[0], m_spData, *m_pGridFct1, vGlobIP, pElem, vCornerCoords, vLocIP, vJT, numIP);
                std::cout << "--- got v1!" << std::endl;
 
                if (m_pGridFct2 != NULL)
                {
                    std::vector<TDataIn> v2(numIP);
                /*  m_spGridFct->set(0.5);
                    m_spGridFct2->set(0.5);*/
                    get_values<elemDim>(&v2[0], m_spData, *m_pGridFct2, vGlobIP, pElem, vCornerCoords, vLocIP, vJT, numIP);
                    std::cout << "--- got v2!" << std::endl;
 
                    for (size_t ip=0; ip<numIP; ++ip)
                    {
 
                        std::cout << std::setprecision(12) << v1[ip] <<" "<<  std::setprecision(12) << v2[ip] << std::endl;
                        v1[ip] -= v2[ip];
                        vValue[ip] = this->product(v1[ip], v1[ip]);
 
                    }
                } else {
 
                    for (size_t ip=0; ip<numIP; ++ip)
                    { vValue[ip] = this->product(v1[ip], v1[ip]); }
 
                }
            };
};


template <typename TGridFunction, typename TDataInput>
class UserDataSpace : public IComponentSpace<TGridFunction>
{
public:
    typedef IComponentSpace<TGridFunction> base_type;
    typedef UserData<TDataInput, TGridFunction::dim> input_user_data_type;
 
    UserDataSpace(const char *fctNames) : base_type(fctNames) {};
    UserDataSpace(const char *fctNames, int order) : base_type(fctNames, order) {};
    ~UserDataSpace() {};
 
    void set_user_data(SmartPtr<input_user_data_type> spData)
    { m_userData = spData; }
 

    using IComponentSpace<TGridFunction>::norm;
    using IComponentSpace<TGridFunction>::distance;
 
    double norm2(TGridFunction& uFine) 
    {
        DeltaSquareIntegrand<TDataInput, TGridFunction> spIntegrand(m_userData, &uFine, NULL, 0.0);
        return IntegrateSubsets(spIntegrand, uFine, base_type::m_ssNames, base_type::m_quadorder, "best");
    }
 
    double distance2(TGridFunction& uFine, TGridFunction& uCoarse) 
    {
        DeltaSquareIntegrand<TDataInput, TGridFunction> spIntegrand(m_userData, &uFine, &uCoarse, 0.0);
        std::cerr << "uFine="<<(void*) (&uFine) << ", uCoarse="<< (void*) (&uCoarse) << std::endl;
        return IntegrateSubsets(spIntegrand, uFine, base_type::m_ssNames, base_type::m_quadorder, "best");
    }
 
protected:
    SmartPtr<input_user_data_type> m_userData;
};
 
 
/*
template <class TDomain, class TAlgebra>
class GridFunctionEstimator :
        public ISubDiagErrorEst<typename TAlgebra::vector_type>
{
protected:
    typedef typename TAlgebra::vector_type TVector;
    typedef GridFunction<TDomain, TAlgebra> grid_function_type;
    typedef IErrorEvaluator<grid_function_type> evaluator_type;
 
    std::vector<SmartPtr<evaluator_type> > m_evaluators;
    number m_refNormValue;
 
public:
    typedef ISubDiagErrorEst<TVector> base_type;
 
    // constructor
    ScaledGridFunctionEstimator() : m_refNormValue(0.0) {}
    ScaledGridFunctionEstimator(number ref) : m_refNormValue(ref) {}
 
    void add(SmartPtr<evaluator_type> eval)
    {
        m_evaluators.push_back(eval);
    }
 
    // apply w/ rel norm
    bool update(SmartPtr<TVector> vUpdate, number alpha,  SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse)
    {
        // typedef ScaleAddLinker<number, TDomain::dim, number> linker_type;
        typedef GridFunctionNumberData<TGridFunction> TNumberData;
 
        // try upcast
        SmartPtr<grid_function_type> uFine = vFine.template cast_dynamic<grid_function_type>();
        SmartPtr<grid_function_type> uCoarse = vCoarse.template cast_dynamic<grid_function_type>();
        if (uFine.invalid() || uCoarse.invalid()) return false;
 
        // error estimate
        if (m_refNormValue<=0.0)
        {
            // relative error estimator
            //number unorm = L2Norm(uFine, m_fctNames.c_str(), m_quadorder);
            //number enorm = alpha*L2Error(uFine, m_fctNames.c_str(), uCoarse, m_fctNames.c_str() ,m_quadorder);
            number unorm = 0.0;
            number enorm = 0.0;
            double est = 0.0;
            for (typename std::vector<evaluator_type>::iterator it = m_evaluators.begin(); it!= m_evaluators.end(); ++it)
            {
                enorm =  alpha * it->distance(uFine, uCoarse);
                unorm =  std::max(it->norm(uFine), 1e-10);
                est += (enorm*enorm)/(unorm*unorm);
                std::cerr << "unorm=" << unorm << "enorm=" << enorm << "est="<<est << std::endl;
            }
 
            base_type::m_est = sqrt(est)/m_evaluators.size();
            std::cerr << "eps="<< base_type::m_est << std::endl;
 
        }
        else
        {
            // weighted error estimator
            number enorm = 0.0;
            for (typename std::vector<evaluator_type>::iterator it = m_evaluators.begin(); it!= m_evaluators.end(); ++it)
            {
                enorm +=  alpha * it->distance(uFine, uCoarse);
            }
            base_type::m_est = enorm/m_refNormValue;
 
            std::cerr << "unorm (FIXED)=" << m_refNormValue << "enorm=" << enorm << "eps="<< base_type::m_est << std::endl;
 
        }
 
        // update
        VecScaleAdd(*vUpdate, 1.0+alpha, *vFine, -alpha, *vCoarse);
        return true;
    }
 
    void set_reference_norm(number norm)
    {m_refNormValue = norm; }
 

    std::string config_string() const
    {
        std::stringstream ss;
        ss << "GridFunctionEstimator:\n";
        for (typename std::vector<GridFunctionEvaluator>::const_iterator it = m_evaluators.begin(); it!= m_evaluators.end(); ++it)
        {
            ss << it->config_string();
        }
        return ss.str();
    }
};
*/


template <class TDomain, class TAlgebra>
class GridFunctionEstimator :
        public ISubDiagErrorEst<typename TAlgebra::vector_type>
{
protected:
    typedef typename TAlgebra::vector_type TVector;
public:
    typedef ISubDiagErrorEst<TVector> base_type;
    typedef GridFunction<TDomain, TAlgebra> grid_function_type;
    typedef IGridFunctionSpace<grid_function_type> subspace_type ;
    //typedef CompositeSpace<grid_function_type> composite_type;
 
protected:
    typedef std::pair<SmartPtr<subspace_type>, number> weighted_obj_type;
 
    number m_refNormValue;
    std::vector<weighted_obj_type> m_spWeightedSubspaces;
 
 
public:
    GridFunctionEstimator() : ISubDiagErrorEst<TVector>(), m_refNormValue(0.0)
    {}
 
 
    GridFunctionEstimator(double ref) : ISubDiagErrorEst<TVector>(), m_refNormValue(ref)
    {}

    void add(SmartPtr<subspace_type> spSubspace)
    {  m_spWeightedSubspaces.push_back(std::make_pair(spSubspace, 1.0)); }
 
    void add(SmartPtr<subspace_type> spSubspace, number sigma)
    {  m_spWeightedSubspaces.push_back(std::make_pair(spSubspace, sigma)); }
 

    bool update(SmartPtr<TVector> vUpdate, number alpha,  SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse)
    {
        // typedefs
        typedef ScaleAddLinker<number, TDomain::dim, number> linker_type;
        typedef GridFunction<TDomain, TAlgebra> grid_function_type;
 
        //  upcast to GridFunction
        SmartPtr<grid_function_type> uFine = vFine.template cast_dynamic<grid_function_type>();
        SmartPtr<grid_function_type> uCoarse = vCoarse.template cast_dynamic<grid_function_type>();
        if (uFine.invalid() || uCoarse.invalid()) return false;
 
        // error estimate
        if (m_refNormValue<=0.0)
        {
            // relative error estimator
            number unorm2 = 0.0;
            number enorm2 = 0.0;
            for (typename std::vector<weighted_obj_type>::iterator it = m_spWeightedSubspaces.begin(); it!= m_spWeightedSubspaces.end(); ++it)
            {
              const double sigma = it->second;
              const double norm2 = it->first->norm2(*uFine);
              const double dist2 = alpha * alpha * (it->first->distance2(*uFine, *uCoarse));
              unorm2 += sigma * norm2;
              enorm2 += sigma * dist2;
              
              UG_LOG("unorm=" << norm2 << "\tenorm=" << dist2 <<  "\tsigma=" << sigma << std::endl);            
            }
 
            base_type::m_est = sqrt(enorm2/unorm2);
            UG_LOG(">>> unorm2=" << unorm2 << "\tenorm2=" << enorm2 << "\teps="<< base_type::m_est << std::endl);
        }
        else
        {
            // weighted error estimator
            number enorm2 = 0.0;
            for (typename std::vector<weighted_obj_type>::iterator it = m_spWeightedSubspaces.begin(); it!= m_spWeightedSubspaces.end(); ++it)
            {
              const double sigma = it->second;
              const double dist2 = alpha * alpha * (it->first->distance2(*uFine, *uCoarse));
              enorm2 += sigma*dist2;
            }
 
            base_type::m_est = sqrt(enorm2)/m_refNormValue;
            UG_LOG("unorm (FIXED)=" << m_refNormValue << "enorm2=" << enorm2 << "eps="<< base_type::m_est << std::endl);
        }
 
        // update
        VecScaleAdd(*vUpdate, 1.0+alpha, *vFine, -alpha, *vCoarse);
        return true;
    }
 
    void set_reference_norm(number norm)
    {m_refNormValue = norm; }
 

    std::string config_string() const
    {
        std::stringstream ss;
        ss << "GridFunctionEstimator:\n";
        for (typename std::vector<weighted_obj_type>::const_iterator it = m_spWeightedSubspaces.begin(); it!= m_spWeightedSubspaces.end(); ++it)
        {
            ss << it->second << "*" << it->first->config_string();
        }
        return ss.str();
    }
};

template <class TDomain, class TAlgebra>
class ScaledGridFunctionEstimator :
        public ISubDiagErrorEst<typename TAlgebra::vector_type>
{
protected:
    typedef typename TAlgebra::vector_type TVector;
 
public:
    typedef ISubDiagErrorEst<TVector> base_type;
    typedef GridFunction<TDomain, TAlgebra> grid_function_type;
    typedef IComponentSpace<grid_function_type> subspace_type;
    typedef CompositeSpace<grid_function_type> composite_type;
 
    // constructor
    ScaledGridFunctionEstimator() : base_type() {}
 
    // add (single subspace)
    void add(SmartPtr<subspace_type> spSubspace)
    {  m_spSubspaces.push_back(spSubspace); }
 
    // add subspaces (from container)
    void add(SmartPtr<composite_type> spCompositeSpace)
    {
        typedef typename composite_type::weighted_obj_type weighted_obj_type;
        const std::vector<weighted_obj_type> &spaces = spCompositeSpace->get_subspaces();
        for (typename std::vector<weighted_obj_type>::const_iterator it = spaces.begin(); it != spaces.end(); ++it)
        {
            m_spSubspaces.push_back(it->first);
        }
    }
 
 
    // apply w/ rel norm
    bool update(SmartPtr<TVector> vUpdate, number alpha,  SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse)
    {
 
        // upcast
        SmartPtr<grid_function_type> uFine = vFine.template cast_dynamic<grid_function_type>();
        SmartPtr<grid_function_type> uCoarse = vCoarse.template cast_dynamic<grid_function_type>();
        if (uFine.invalid() || uCoarse.invalid()) return false;
 
        // error estimate
        number est = 0.0;
        for (typename std::vector<SmartPtr<subspace_type> >::iterator it = m_spSubspaces.begin();
                it!= m_spSubspaces.end(); ++it)
        {
            // use sub-diagonal error estimator (i.e. multiply with alpha)
            double enorm2 =  (alpha*alpha) * (*it)->distance2(*uFine, *uCoarse);
            double unorm2 =  std::max((*it)->norm2(*uFine), 1e-10*enorm2);
            est += (enorm2)/(unorm2);
            UG_LOGN("unorm2=" << unorm2 << "\tenorm2=" << enorm2 << "\tratio2="<< (enorm2)/(unorm2) << "est2=" << est);
        }
 
        base_type::m_est = sqrt(est)/m_spSubspaces.size();
        UG_LOGN("eps="<< base_type::m_est);
 
        // update
        VecScaleAdd(*vUpdate, 1.0+alpha, *vFine, -alpha, *vCoarse);
        return true;
    }
 

    std::string config_string() const
    {
        std::stringstream ss;
        ss << "ScaledGridFunctionEstimator:\n";
        for (typename std::vector<SmartPtr<subspace_type> >::const_iterator it = m_spSubspaces.begin();
                it!= m_spSubspaces.end(); ++it)
        {
            ss << (*it)->config_string();
        }
        return ss.str();
    }
 
protected:
 
    std::vector<SmartPtr<subspace_type> > m_spSubspaces;
 
};
 


template <class TDomain, class TAlgebra>
class CompositeGridFunctionEstimator :
        public ISubDiagErrorEst<typename TAlgebra::vector_type>
{
protected:
    typedef typename TAlgebra::vector_type TVector;
 
public:
    typedef ISubDiagErrorEst<TVector> base_type;
    typedef GridFunction<TDomain, TAlgebra> grid_function_type;
    typedef IComponentSpace<grid_function_type> subspace_type;
    typedef CompositeSpace<grid_function_type> composite_type;
 
    // constructor
    CompositeGridFunctionEstimator() : base_type(), m_strictRelativeError(0) {}
 
    void use_strict_relative_norms(int b)//bool b)
    { m_strictRelativeError  = b; }
 
    void use_strict_relative_norms(bool b){
        if(b){
            m_strictRelativeError=1;
        }else{
            m_strictRelativeError=0;
        }
    }
 
    // add (single subspace)
    void add(SmartPtr<subspace_type> spSubspace)
    {  m_spSubspaces.push_back(spSubspace); }
 
    // add subspaces (from container)
    void add(SmartPtr<composite_type> spCompositeSpace)
    {
        typedef typename composite_type::weighted_obj_type weighted_obj_type;
        const std::vector<weighted_obj_type> &spaces = spCompositeSpace->get_subspaces();
        for (typename std::vector<weighted_obj_type>::const_iterator it = spaces.begin(); it != spaces.end(); ++it)
        {
            m_spSubspaces.push_back(it->first);
        }
    }
 
 
    // apply w/ rel norm
    bool update(SmartPtr<TVector> vUpdate, number alpha,  SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse)
    {
 
        // upcast
        SmartPtr<grid_function_type> uFine = vFine.template cast_dynamic<grid_function_type>();
        SmartPtr<grid_function_type> uCoarse = vCoarse.template cast_dynamic<grid_function_type>();
        if (uFine.invalid() || uCoarse.invalid()) return false;
 
        size_t numFct=uFine->num_fct();
        std::vector<double> vcmp_e2, vcmp_u2;
        vcmp_e2.resize(numFct, 0);
        vcmp_u2.resize(numFct, 0);
 
        // error estimate
        double enorm2 = 0.0;
        double unorm2 = 0.0;
        const double SMALL = 1e-10;
 
        double cmp_rel=0.0;
        double cmp_rel4=0.0;
        double max_rel = 0.0;
        double max_rel_e4=0.0;
 
        for (typename std::vector<SmartPtr<subspace_type> >::iterator it = m_spSubspaces.begin();
                it!= m_spSubspaces.end(); ++it)
        {
            // use sub-diagonal error estimator (i.e. multiply with alpha)
            std::string cmp=(*it)->function_name();
            size_t cmp_id= uFine->fct_id_by_name(cmp.c_str());
            double cmp_e2 = (alpha*alpha) * (*it)->distance2(*uFine, *uCoarse);
            double cmp_u2 = (*it)->norm2(*uFine);
            vcmp_e2[cmp_id]+=cmp_e2;
            vcmp_u2[cmp_id]+=cmp_u2;
            cmp_rel=std::max(cmp_rel,(cmp_e2==0)?0:cmp_e2/(cmp_u2 + SMALL/ (1.0+cmp_e2+cmp_u2)));
            cmp_rel4=std::max(cmp_rel4,(cmp_e2==0)?0:cmp_e2/std::sqrt(cmp_u2*cmp_u2+SMALL*SMALL));
            UG_LOGN("cmp=" << cmp << "(id=" << cmp_id << "):\tui-2=" << cmp_u2 << "\tei-2=" <<  cmp_e2<<
                    "\tratio2="<< (vcmp_e2[cmp_id])/(vcmp_u2[cmp_id]) <<
                    "\tmax.rel(squared)="<< cmp_rel<<"\tmax.rel(ei4)="<<cmp_rel4
            );
        }
 
        for(size_t cmp_id=0; cmp_id<numFct; ++cmp_id){
            double cmp_e2=vcmp_e2[cmp_id];
            double cmp_u2=vcmp_u2[cmp_id];
            // |delta_i|/|u_i|
            max_rel = std::max(max_rel, (cmp_e2==0)?0:cmp_e2/
                    (cmp_u2 + SMALL/ (1.0+cmp_e2+cmp_u2)));
                    //std::max(cmp_u2, SMALL*cmp_e2));
            max_rel_e4=std::max(max_rel_e4, (cmp_e2==0)?0:cmp_e2/std::sqrt(cmp_u2*cmp_u2+SMALL*SMALL));
 
            enorm2 += cmp_e2;
            unorm2 += cmp_u2;
 
            UG_LOGN("summe by cmp_id=" << cmp_id << ":\tui-2=" << cmp_u2 << "\tei-2=" <<  cmp_e2<<
                    "\tratio2="<< (enorm2)/(unorm2) <<
                    "\tmax.rel(squared)="<< max_rel<<"\tmax.rel(ei4)="<<max_rel_e4);
        }
 
        //prevent division by zero
        /*base_type::m_est  = (m_strictRelativeError) ? sqrt(max_rel) :
                sqrt(enorm2/std::max(unorm2, SMALL*enorm2));*/
 
            switch(m_strictRelativeError)
                {
                       case 0:
                                base_type::m_est=sqrt(enorm2/std::max(unorm2, SMALL*enorm2)); /* absolute: */
                                break;
                       case 1:
                                base_type::m_est=sqrt(cmp_rel); /* relative max: max(e2/(u2+SMALL/(1+e2+u2))) for each cmp on each subset */
                                break;
                       case 2:
                                base_type::m_est=sqrt(cmp_rel4); /* relative max: max(e2/sqrt(u2*u2+SMALL*SMALL)) for each cmp on each subset*/
                                break;
                       case 3:
                                base_type::m_est=sqrt(max_rel); /* relative max based on cmp: max(sum_si e2/( sum_si u2+SMALL/(1+sum_si e2+ sum_si u2))) for each cmp on whole domain.*/
                                break;
                       case 4:
                                base_type::m_est=sqrt(max_rel_e4); /* relative max based on cmp: max(sum_si e2/sqrt( sum_si u2* sum_si u2+SMALL*SMALL)) for each cmp on whole domain.*/
                                 break;
                }
        UG_LOGN("eps="<< base_type::m_est<<" with strictRelativeError="<<m_strictRelativeError);
 
        // update
        VecScaleAdd(*vUpdate, 1.0+alpha, *vFine, -alpha, *vCoarse);
        return true;
    }
 

    std::string config_string() const
    {
        std::stringstream ss;
        ss << "CompositeGridFunctionEstimator:\n";
        for (typename std::vector<SmartPtr<subspace_type> >::const_iterator it = m_spSubspaces.begin();
                it!= m_spSubspaces.end(); ++it)
        {
            ss << (*it)->config_string();
        }
        return ss.str();
    }
 
protected:
 
    std::vector<SmartPtr<subspace_type> > m_spSubspaces;
    //bool m_strictRelativeError;
    int m_strictRelativeError;
};
 
 
 
template <typename TVector>
class AitkenNevilleTimex
{
    public:
        typedef TVector vector_type;
 
    public:

        AitkenNevilleTimex(std::vector<size_t> nsteps)
        : m_stepsize(0.0),
          m_subdiag (make_sp(new Norm2Estimator<TVector>())),
          m_num_steps(nsteps),
          m_solution(nsteps.size()),
          m_subdiag_error_est(nsteps.size(), INFINITY)
        {};
 
        AitkenNevilleTimex(std::vector<size_t> nsteps, SmartPtr<ISubDiagErrorEst<vector_type> > error)
        : m_stepsize(0.0),
          m_subdiag(error),
          m_num_steps(nsteps),
          m_solution(nsteps.size()),
          m_subdiag_error_est(nsteps.size(), INFINITY)
        {};
 
        virtual ~AitkenNevilleTimex() {}
 
        void set_global_stepsize(number H) {m_stepsize=H;}
        number get_global_stepsize() {return m_stepsize;}

         void set_solution(SmartPtr<vector_type> soli, int i)
         { m_solution[i] = soli; }

        SmartPtr<vector_type> get_solution(size_t i)
        { return m_solution[i]; }

        void set_error_estimate(SmartPtr<ISubDiagErrorEst<vector_type> > subdiag)
        {
            m_subdiag = subdiag;
        }
 
 
        const std::vector<number>& get_error_estimates() const
        {return m_subdiag_error_est; }

        number get_error_estimate(int k) const
        { return m_subdiag_error_est[k];}
 

        /*number get_error_estimate()
        {
            std::vector<number>::iterator best = std::min_element(m_subdiag_error_est.begin(), m_subdiag_error_est.end());
            return *best;
        }

        int get_best_index() const
        {
            std::vector<number>::iterator best = std::min_element(m_subdiag_error_est.begin(), m_subdiag_error_est.end());
            //std::cout << "min element at: " << std::distance(std::begin(m_subdiag_error_est), best);
            return std::distance(m_subdiag_error_est.begin(), best);
        }
 
         */

        void apply(size_t nstages, bool with_error=true)
        {
            UG_ASSERT(nstages <= m_solution.size(),
                     "Dimensions do not match:"  << nstages << ">" << m_solution.size());
 
            if (with_error)
            {
                // reset (for safety reasons...)
                for (size_t k=1; k<m_solution.size(); ++k)
                { m_subdiag_error_est[k] = INFINITY; }
            }
 
            //m_subdiag_error_est[0] = ;
            // process columns (left to right)
            for (size_t k=1; k<nstages; ++k)
            {
 
                // process rows (bottom-up -> recycling memory)
                for (size_t i=nstages-1; i>=k; --i)
                {
                    UG_ASSERT(m_solution[i].valid(), "Invalid SmarPtr!");
                    UG_ASSERT(m_solution[i-1].valid(), "Invalid SmarPtr!");
 
                    SmartPtr<vector_type> solcoarse = m_solution[i-1];
                    SmartPtr<vector_type> solfine = m_solution[i];
 
                    // (2^p -1)
                    // m_solution[i] += (1.0/scal)*(m_solution[i]- m_solution[i-1]);
                    const number scaling = ((1.0*m_num_steps[i])/(1.0*m_num_steps[i-k])-1.0);
                    UG_LOG("scaling="<<i << ","<< k <<
                            ": ns["<<i<<"]="<< m_num_steps[i] <<
                            "ns2["<<i-k<<"]=" <<  m_num_steps[i-k] <<"=" << scaling << std::endl);
 
                    if (with_error && (i==k))
                    {
                        // compute subdiagonal error estimate
                        m_subdiag->update(solfine, (1.0/scaling), solfine,  solcoarse);
                        number subdiag_error_est=m_subdiag->get_current_estimate();
 
                        //m_subdiag_error_est[k] = sqrt(subdiag_error_est*scaling);
                        m_subdiag_error_est[k]=subdiag_error_est; //*scaling;
 
                        UG_LOG(" ErrorEst["<< k<<"]=" << m_subdiag_error_est[k] << ";" << std::endl);
                    }
                    else
                    {
                        // standard case
                        VecScaleAdd(*solfine, (1.0+1.0/scaling), *solfine,
                                              -(1.0/scaling), *solcoarse);
 
                    }
                } // rows i
 
            } // columns k
 
        }

        void apply()
        {
            //UG_ASSERT(m_num_steps.size() == m_solution.size(), "Dimensions do not match");
            const size_t nstages = m_num_steps.size();
            apply(nstages);
        }
 
 
 
protected:
        number substep(size_t i) {return m_stepsize/m_num_steps[i];}
 
private:
        number m_stepsize;
        static const int m_order=1;

        SmartPtr<ISubDiagErrorEst<vector_type> > m_subdiag;

        std::vector<size_t> m_num_steps;

        std::vector<SmartPtr<vector_type> > m_solution;

        std::vector<number> m_subdiag_error_est;
 
 
};
 
 

/*
class TimeStepEstimator{
 
public:
    TimeStepEstimator(double tol) :
    m_rho(0.9), m_gamma(2.0), m_tol(tol)
    {
 
    }
 
    bool check(double eps, double &factor)
    {
        if (eps>=tol) return false;
 
        double val = (m_rho*m_tol)/eps;
        factor = pow(val, m_gamma);
 
    }
 
 
private:
    double m_rho;
    double m_gamma;
    double m_tol;
};
*/
 
}
#endif