Pagani2016filter.H

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
     ██╗████████╗██╗  ██╗ █████╗  ██████╗ █████╗       ███████╗██╗   ██╗
     ██║╚══██╔══╝██║  ██║██╔══██╗██╔════╝██╔══██╗      ██╔════╝██║   ██║
     ██║   ██║   ███████║███████║██║     ███████║█████╗█████╗  ██║   ██║
     ██║   ██║   ██╔══██║██╔══██║██║     ██╔══██║╚════╝██╔══╝  ╚██╗ ██╔╝
     ██║   ██║   ██║  ██║██║  ██║╚██████╗██║  ██║      ██║      ╚████╔╝
     ╚═╝   ╚═╝   ╚═╝  ╚═╝╚═╝  ╚═╝ ╚═════╝╚═╝  ╚═╝      ╚═╝       ╚═══╝
 
 * In real Time Highly Advanced Computational Applications for Finite Volumes
 * Copyright (C) 2017 by the ITHACA-FV authors
-------------------------------------------------------------------------------
License
    This file is part of ITHACA-FV
    ITHACA-FV is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
    ITHACA-FV 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.
    You should have received a copy of the GNU Lesser General Public License
    along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
Class
    Pagani2016filter
Description
    A class that implements the Ensamble Kalman Filter in the version by Pagani et al. 2016
SourceFiles
    Pagani2016filter.C
\*---------------------------------------------------------------------------*/

 
#ifndef pagani2016filter_H
#define pagani2016filter_H
 
#include "ITHACAutilities.H"
#include "muq2ithaca.H"
#include "ensembleClass.H"
 
 
 
 
namespace ITHACAmuq
{
//--------------------------------------------------------------------------
class Pagani2016filter
{
        int Nseeds = 0;
 
        int timeSampI = 0;
        int stateSize = 0;
        int parameterSize = 0;
        int observationSize = 0;
 
        int timeStepI = 0;
        double startTime;
        double deltaTime;
        double endTime;
        int Ntimes = 0;
 
        bool initialStateFlag = 0;
        bool parameterPriorFlag = 0;
        bool modelErrorFlag = 0;
        bool measurementNoiseFlag = 0;

        Eigen::VectorXd initialStateMean;

        Eigen::MatrixXd initialStateCov;

        Eigen::VectorXd parameterPriorMean;

        Eigen::MatrixXd parameterPriorCov;

        Eigen::MatrixXd parameterMean;

        Eigen::MatrixXd stateMean;

        Eigen::VectorXd timeVector;

        int observationStart;

        int observationDelta;

        Eigen::VectorXi observationBoolVec;

        Eigen::MatrixXd trueObservations;

        std::shared_ptr<muq::Modeling::Gaussian> initialStateDensity;

        std::shared_ptr<muq::Modeling::Gaussian> parameterPriorDensity;

        Eigen::MatrixXd parameterObservation_crossCov;

        Eigen::MatrixXd stateObservation_crossCov;

        Eigen::MatrixXd observation_Cov;
 
 
        int innerLoopI = 0;
 
 
    public:
        // Constructors
        Pagani2016filter();
 
        Pagani2016filter(int _Nseeds);
 
        virtual ~Pagani2016filter() = default;

        std::shared_ptr<muq::Modeling::Gaussian> measNoiseDensity;

        std::shared_ptr<muq::Modeling::Gaussian> parameterErrorDensity;

        ensemble stateEns;

        ensemble oldStateEns;

        ensemble parameterEns;

        ensemble observationEns;

        Eigen::MatrixXd state_maxConf;

        Eigen::MatrixXd state_minConf;

        Eigen::MatrixXd parameter_maxConf;

        Eigen::MatrixXd parameter_minConf;
 
        // Functions
 
        //--------------------------------------------------------------------------
        int getNumberOfSamples();
 
        //--------------------------------------------------------------------------
        int getObservationSize();
 
        //--------------------------------------------------------------------------
        double getTime();
 
        //--------------------------------------------------------------------------
        double getTime(int _timeStepI);
 
        //--------------------------------------------------------------------------
        int getTimeStep();
 
        //--------------------------------------------------------------------------
        Eigen::VectorXd getTimeVector();
 
        //--------------------------------------------------------------------------
        Eigen::MatrixXd getStateMean();
 
        //--------------------------------------------------------------------------
        Eigen::MatrixXd getParameterMean();
 
        //--------------------------------------------------------------------------
        void computeStateMean(int _index);
 
        //--------------------------------------------------------------------------
        void computeParameterMean(int _index);
 
        //--------------------------------------------------------------------------
        int getObservationCounter();
 
        //--------------------------------------------------------------------------
        int getObservationTimestep(int _timeSampleI);
 
        //--------------------------------------------------------------------------
        void setTrueObservations(Eigen::MatrixXd _observations);
 
        //--------------------------------------------------------------------------
        void setTime(double _startTime, double _deltaTime, double _endTime);
 
        //--------------------------------------------------------------------------
        void setObservationTime(int _observationStart, int _observationDelta);
 
        //--------------------------------------------------------------------------
        void setMeasNoise(double cov);
 
        //--------------------------------------------------------------------------
        void setParameterError(double cov);
 
        //--------------------------------------------------------------------------
        void setInitialStateDensity(Eigen::VectorXd _mean, Eigen::MatrixXd _cov);
 
        //--------------------------------------------------------------------------
        void sampleInitialState();
 
        //--------------------------------------------------------------------------
        void sampleInitialParameter();
 
        //--------------------------------------------------------------------------
        void setParameterPriorDensity(Eigen::VectorXd _mean, Eigen::MatrixXd _cov);
 
        //--------------------------------------------------------------------------
        Eigen::MatrixXd ensembleFromDensity(std::shared_ptr<muq::Modeling::Gaussian>
                                            _density);
 
        //--------------------------------------------------------------------------
        virtual void stateProjection() = 0;
 
        //--------------------------------------------------------------------------
        virtual void observeState() = 0;
 
        //--------------------------------------------------------------------------
        void setObservationSize(int _size);
 
        //--------------------------------------------------------------------------
        void setStateSize(int _size);
 
        //--------------------------------------------------------------------------
        int getStateSize();
 
        //--------------------------------------------------------------------------
        void setParameterSize(int _size);
 
        //--------------------------------------------------------------------------
        int getParameterSize();
 
        //--------------------------------------------------------------------------
        void update();
 
        //--------------------------------------------------------------------------
        void run(word outputFolder);
};
}
#endif