Fang2017filter.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
    Fang2017filter
Description
    A class that implements the Ensamble Kalman Filter in the version by Fang et al. 2017
SourceFiles
    Fang2017filter.C
\*---------------------------------------------------------------------------*/

 
#ifndef fang2017filter_H
#define fang2017filter_H
 
#include "ITHACAutilities.H"
#include "muq2ithaca.H"
#include "ensembleClass.H"
 
 
 
namespace ITHACAmuq
{
//--------------------------------------------------------------------------
class Fang2017filter
{
        int Nsamples = 0;
 
        int stateSize = 0;
        int parameterSize = 0;
        int observationSize = 0;
 
        int timeStepI = 0;
        double startTime;
        double deltaTime;
        double endTime;
        int Ntimes;
 
        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 observations;

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

        std::shared_ptr<muq::Modeling::Gaussian> parameterPriorDensity;
 
 
        int innerLoopI = 0;
 
 
    public:
        // Constructors
        Fang2017filter();
 
        Fang2017filter(int _Nsamples);
 
        virtual ~Fang2017filter() = default;

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

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

        ensemble stateEns;

        ensemble oldStateEns;

        ensemble parameterEns;

        ensemble observationEns;

        ensemble jointEns;

        Eigen::MatrixXd state_maxConf;

        Eigen::MatrixXd state_minConf;

        Eigen::MatrixXd parameter_maxConf;

        Eigen::MatrixXd parameter_minConf;
 
        // Functions
 
        //--------------------------------------------------------------------------
        int getNumberOfSamples();
 
        //--------------------------------------------------------------------------
        double getTime();
 
        //--------------------------------------------------------------------------
        double getTime(int _timeStepI);
 
        //--------------------------------------------------------------------------
        int getTimeStep();
 
        //--------------------------------------------------------------------------
        Eigen::VectorXd getTimeVector();
 
        //--------------------------------------------------------------------------
        Eigen::MatrixXd getStateMean();
 
        //--------------------------------------------------------------------------
        void setObservations(Eigen::MatrixXd _observations);
 
        //--------------------------------------------------------------------------
        void setTime(double _startTime, double _deltaTime, double _endTime);
 
        //--------------------------------------------------------------------------
        void setObservationTime(int _observationStart, int _observationDelta);
 
        //--------------------------------------------------------------------------
        void setModelError(double cov, bool univariate = 0);
 
        //--------------------------------------------------------------------------
        void setMeasNoise(double cov);
 
        //--------------------------------------------------------------------------
        void setInitialStateDensity(Eigen::VectorXd _mean, Eigen::MatrixXd _cov);
 
        //--------------------------------------------------------------------------
        void sampleInitialState();
 
        //--------------------------------------------------------------------------
        void setParameterPriorDensity(Eigen::VectorXd _mean, Eigen::MatrixXd _cov);
 
        //--------------------------------------------------------------------------
        void sampleParameterDist();
 
        //--------------------------------------------------------------------------
        Eigen::MatrixXd ensembleFromDensity(
            std::shared_ptr<muq::Modeling::Gaussian> _density);
 
        //--------------------------------------------------------------------------
        void buildJointEns();
 
        //--------------------------------------------------------------------------
        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 updateJointEns(Eigen::VectorXd _observation);
 
        //--------------------------------------------------------------------------
        void run(int innerLoopMax, word outputFolder);
};
}
#endif