---

StochasticSensitivityAnalysis.h

Go to the documentation of this file.

// StochasticSensitivityAnalysis.h: interface for the CStochasticSensitivityAnalysis class.
//

#if !defined(AFX_STOCHASTICSENSITIVITYANALYSIS_H__6F1FADC6_28C2_4963_AFA8_607527338B17__INCLUDED_)
#define AFX_STOCHASTICSENSITIVITYANALYSIS_H__6F1FADC6_28C2_4963_AFA8_607527338B17__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

#ifdef _WIN32
#include <mkl.h>
#include <strstream>
#else
#include "lapack.h"
#include <algorithm>
using namespace std;
#define __min min
#define __max max
#endif

#include <iostream>
#include <fstream>
#include <vector>
#include <iomanip>
#include <sstream>
#include "SubjectObserver.h"
#include "NLLSMinimizable.h"
#include "Rand.h"
#include "ParameterFilter.h"
#include "EnsembleMember.h"
using namespace std;

class CStochasticSensitivityAnalysis : public Subject  
{
public:
        CStochasticSensitivityAnalysis();
        CStochasticSensitivityAnalysis(bool readFlag, bool hessianFlag, int nEnsembleSize, int seed, CParameterFilter *pFilter, double funcAccuracy, double temperature=1.0);
        virtual ~CStochasticSensitivityAnalysis();
        // computes the Levenberg-Marquardt approximation to the Hessian
        void ComputeApproximateHessian(double *pParameters, NLLSMinimizable *pNLLS);
        // computes the true Hessian
        void ComputeTrueHessian(double *pParameters, NLLSMinimizable *pNLLS);
        // reads the Hessian in from a file
        void ReadHessian();
        void DoAnalysis(double *pFittedParameters, Minimizable *pMinimizable);
        double SampleGaussian(double *pParDelta);
        double SampleWithoutExplicitHessian(double *pParDelta);
        void ConstructTrialMove(double *pParTrial, double *pParCurrent, double *pParDelta);
        void AddParametersToEnsemble(double *pParSave, double *resid, double quadratic, double cost);
        bool AcceptMove(double quadratic, double deltaCost);
        double GetRWValue() const {return m_dRWValue;}
        std::vector<CEnsembleMember *> *GetEnsemble() {return &m_vpEnsemble;}
        const double **GetEigenvectorMatrix() const {return (const double **)m_pdMatrixV;}
        double GetRescale() {return m_dRescale;}
        double *GetSingularValues() {return m_pdSingularValues;}
private:
        // Set to true to decompose the Hessian 
        bool m_bSVDFlag;
        // Set to true to use the true Hessian, false for LM Hessian
        bool m_bHessianFlag;
        // Set to true to read the (previously computed) Hessian in from a file
        bool m_bReadFlag;
        // Number of residuals in cost function
        int m_iNResiduals;
        // Number of parameters in cost function
        int m_iNParameters;
        // Number of samples
        int m_iNEnsembleSize;
        int m_iRNGSeed;
        // Thermodynamic temperature
        double m_dTemperature;
        // control parameter for stepsize to sample efficiently
        double m_dRescale;
        // Acceptance Ratio for the duration of the sampling
        double m_dAcceptanceRatio;
        // This is needed in order to compute the optimal stepsize for finite differencing
        double m_dFuncAccuracy;
        // Size of parameter fluctuations as estimated by a simple 1D RW
        double m_dRWValue;
        // Used in the computation of RW estimate
        double m_dMNormSquared;
        // Vector of singular values of the Hessian
        double *m_pdSingularValues;
        // Matrix of Eigenvectors (columns are eigenvectors of H)
        double **m_pdMatrixV;
        // The Hessian
        double **m_pdHessian;
        // Pseudo-Square root of the Hessian
        double **m_pdSamplingMatrix;
        // The Jacobian (Square root of the LM Hessian)
        double **m_pdJacobian;
        CParameterFilter *m_pFilter;
        Rand *m_pRNG;
        std::vector<CEnsembleMember *> m_vpEnsemble;
};

#endif // !defined(AFX_STOCHASTICSENSITIVITYANALYSIS_H__6F1FADC6_28C2_4963_AFA8_607527338B17__INCLUDED_)

---