modellingbase.h

/******************************************************************************
 *   Copyright (C) 2005-2024 by the GIMLi development team                    *
 *   Carsten Rücker carsten@resistivity.net                                   *
 *                                                                            *
 *   Licensed under the Apache License, Version 2.0 (the "License");          *
 *   you may not use this file except in compliance with the License.         *
 *   You may obtain a copy of the License at                                  *
 *                                                                            *
 *       http://www.apache.org/licenses/LICENSE-2.0                           *
 *                                                                            *
 *   Unless required by applicable law or agreed to in writing, software      *
 *   distributed under the License is distributed on an "AS IS" BASIS,        *
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
 *   See the License for the specific language governing permissions and      *
 *   limitations under the License.                                           *
 *                                                                            *
 ******************************************************************************/
 
#ifndef _GIMLI_MODELLINGBASE__H
#define _GIMLI_MODELLINGBASE__H
 
#include "gimli.h"
#include "matrix.h"
//#include "blockmatrix.h"
 
namespace GIMLI{
 
//class H2SparseMapMatrix;

class DLLEXPORT ModellingBase{
public:
 
    ModellingBase(bool verbose=false) ;
 
    ModellingBase(DataContainer & dataContainer, bool verbose=false);
 
    ModellingBase(const Mesh & mesh, bool verbose=false);
 
    ModellingBase(const Mesh & mesh, DataContainer & dataContainer, bool verbose=false);
 
    virtual ~ModellingBase();

    void setVerbose(bool verbose);

    inline bool verbose() const {return verbose_;}
 
    virtual RVector response(const RVector & model) {
        throwError(WHERE_AM_I + " you need to implement a response "
        "function.");
        return RVector(0);
    }

    virtual RVector response_mt(const RVector & model, Index i=0) const {
        throwError(WHERE_AM_I + " if you want to use read only response "
        "function to use in multi threading environment .. you need to "
        " implement me");
        return RVector(0);
    }
 
    inline RVector operator() (const RVector & model){ return response(model); }

    void setData(DataContainer & data);

    DataContainer & data() const;

    virtual RVector createDefaultStartModel() { return RVector(0); }

    virtual RVector startModel();
 
    virtual void setStartModel(const RVector & startModel);

    void setMesh(const Mesh & mesh, bool ignoreRegionManager=false);
 
    inline Mesh * mesh() { return mesh_; }
 
    void createRefinedForwardMesh(bool refine=true, bool pRefine=false);

    void deleteMesh();

    virtual void setJacobian(MatrixBase * J);

    virtual void createJacobian(const RVector & model);

    virtual void createJacobian(const RVector & model, const RVector & resp);

    virtual void createJacobian_mt(const RVector & model, const RVector & resp);

    virtual void initJacobian();

    MatrixBase * jacobian() { return jacobian_; }

    MatrixBase * jacobian() const { return jacobian_; }

    virtual RMatrix & jacobianRef() const {
        if (! jacobian_) {
            throwError(WHERE_AM_I + " Jacobian matrix is not initialized.");
        }
        return *dynamic_cast< RMatrix * >(jacobian_);
    }
 
    virtual RMatrix & jacobianRef() {
        if (! jacobian_) {
            throwError(WHERE_AM_I + " Jacobian matrix is not initialized.");
        }
        return *dynamic_cast< RMatrix * >(jacobian_);
    }

    virtual void clearJacobian(){ jacobian_->clear(); }

//     inline void setConstraintMatrix(const RSparseMapMatrix & C){ C_ = C; }
//     inline const RSparseMapMatrix & constraintMatrix() const { return C_; }
//         /*! Clear Constraints matrix */
    virtual void setConstraints(MatrixBase * C);
 
    virtual void clearConstraints();
 
    virtual void initConstraints();
 
    virtual void createConstraints();
 
    virtual MatrixBase * constraints();
 
    virtual MatrixBase * constraints() const;
 
    virtual RSparseMapMatrix & constraintsRef() const;
 
    virtual RSparseMapMatrix & constraintsRef();
 
    const RMatrix & solution() const { return solutions_; }
 
    void mapModel(const RVector & model, double background=0);

    RVector createMappedModel(const RVector & model, double background=-9e99) const;
 
    void setRegionManager(RegionManager * reg);
 
    const RegionManager & regionManager() const;
 
    RegionManager & regionManager();
 
    RegionManager & regionManagerRef(){ return regionManager(); }
 
    bool verbose() { return verbose_; }

    Region * region(int marker);
 
 
    RVector createStartModel();

    RVector createStartVector();
 
    void initRegionManager();

    void setThreadCount(Index nThreads);

    Index threadCount();

    void setMultiThreadJacobian(Index nThreads);

    inline Index multiThreadJacobian() const { return nThreadsJacobian_; }
 
 
protected:
 
    virtual void init_();
 
    virtual void deleteMeshDependency_(){}
 
    virtual void updateMeshDependency_(){}
 
    virtual void updateDataDependency_(){}
 
    void setMesh_(const Mesh & mesh, bool update=true);
 
    Mesh                    * mesh_;
 
    DataContainer           * dataContainer_;
 
    MatrixBase              * jacobian_;
    bool                    ownJacobian_;
 
    MatrixBase              * constraints_;
    bool                    ownConstraints_;
 
    RMatrix                 solutions_;
 
    RVector                 startModel_;
 
    bool                    verbose_;
 
    bool                    regionManagerInUse_;
    bool                    ownRegionManager_;
 
    Index                   nThreads_;
    Index                   nThreadsJacobian_;
 
private:
    RegionManager            * regionManager_;
 
};

class DLLEXPORT LinearModelling : public ModellingBase {
public:
    LinearModelling(Mesh & mesh, MatrixBase & A, bool verbose=false)
        : ModellingBase(mesh, verbose){
            setJacobian(&A);
        }
 
    LinearModelling(Mesh & mesh, MatrixBase * A, bool verbose=false)
        : ModellingBase(mesh, verbose){
            setJacobian(A);
        }
 
    LinearModelling(MatrixBase & A, bool verbose=false);
 
    virtual ~LinearModelling() { }

    RVector response(const RVector & model);
 
    virtual void createJacobian(const RVector & model);
 
    RVector createDefaultStartModel();
 
protected:
 
};
 
 
 
} // namespace GIMLI
 
#endif // MODELLINGBASE__H