PostTheta.h

#ifndef POST_THETA_H
#define POST_THETA_H
 
#include <random>
#include <vector>
#include <iostream>
#include <fstream>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
 
#include <omp.h>
#include "mpi.h"
 
// std::setwd print out
#include <iomanip>
 
#include <Eigen/Dense>
#include <Eigen/Core>
#include <Eigen/SparseCholesky>
//#include <Eigen/CholmodSupport>
#include <unsupported/Eigen/KroneckerProduct>
 
//#include "../read_write_functions.h"
 
//#include "solver_cholmod.h" -> pardiso can do inversion now
#include "PardisoSolver.h"
#include "BTASolver.h"
//#include "BTASolver_dummy.h"
#include "EigenCholSolver.h"
 
//#include "Hyperparameters.h"
 
//#define SMART_GRAD
 
//#define PRINT_MSG
//#define PRINT_TIMES
//#define RECORD_TIMES
//#define DATA_SYNTHETIC
 
using namespace Eigen;
using namespace std;
 
// typedef Eigen::Matrix<double, Dynamic, Dynamic> Mat;
typedef Eigen::SparseMatrix<double> SpMat;
typedef Eigen::SparseMatrix<double, RowMajor> SpRmMat;
typedef Eigen::VectorXd Vect;
 
 
class PostTheta{
 
    private:
 
    int MPI_size;       
    int MPI_rank;       
    int threads_level1; 
    int threads_level2;
 
 
    int ns;             
    int nt;             
    int nss;            
    int nb;             
    int no;             
    int nst;            
    int nu;             
    int n;              
    size_t nnz_Qst;
    size_t nnz_Qs;
 
    int dim_th;         
    int dim_spatial_domain; 
    string manifold;    
    int dim_grad_loop;  
    int num_solvers;    
    Solver* solverQ;
    Solver* solverQst;
 
    int threadID_solverQst;
    int threadID_solverQ;
 
    string likelihood;  
    Vect extraCoeffVecLik;
 
    string solver_type;
 
    std::string prior;  
    int fct_count;      
    int iter_count;     
    int iter_acc;
    Vect y;         
    Vect theta_prior_param; 
    Vector3i dimList;
#if 0
    Hyperparameters* theta_prior_test;
    Hyperparameters* theta_test;
    //Hyperparameters theta_prior_test = new Hyperparameters(dim_spatial_domain, manifold, dimList, theta_prior_param, theta_prior_param);
#endif
    // either Ax or B used
    SpMat Ax;           
    MatrixXd B;         
    // used in spatial and spatial-temporal case
    SpMat c0;           
    SpMat g1;           
    SpMat g2;           
    // only used in spatial-temporal case
    SpMat g3;           
    SpMat M0;           
    SpMat M1;           
    SpMat M2;           
    SpMat Qb;           
    SpMat Qu;           
    SpMat Qst;
    SpMat Qs;
    SpMat Qx;           
    SpMat Qxy;          
    double yTy;         
    Vect BTy;           
    Vect AxTy;          
    SpMat AxTAx;        
    Vect mu_initial;
    Vect mu;            
    Vect mu_midpoint;    
    // new attempt, store previous mu for each stencil point 
    // figure out later what to do for Hessian ...
    MatrixXd mu_matrix; 
    Vect t_grad;        
    double min_f_theta; 
    double w_sum;       
    int no_f_eval;      
    ArrayXi task_to_rank_list_grad;
 
    MatrixXd G;         
#ifdef SMART_GRAD
    bool thetaDiff_initialized; 
    VectorXd theta_prev;
    MatrixXd ThetaDiff;
#endif
 
    const bool constr;      
    const MatrixXd Dx;         
    const MatrixXd Dxy;         
    //MatrixXd V;           /**< V = Q^-1 * t(D)                                */
    //MatrixXd W;           /**< W = A*V                                        */
    //Vect U;               /**< U = W^-1 * t(V)                                */
 
    const bool validate;
    const Vect w;
 
    bool printed_eps_flag = false;
 
#ifdef RECORD_TIMES
    std::string log_file_name;
    // to record times 
    double t_Ftheta_ext;
    double t_priorHyp;
    double t_priorLat;
    double t_priorLatAMat;
    double t_likel;
    double t_condLat;
    double t_condLatAMat;
    double t_thread_nom;
    double t_thread_denom;
#endif
 
    // these have to always be defined otherwise compile error.
    // always get measured but not written to file
    double t_priorLatChol;
    double t_condLatChol;
    double t_condLatSolve;
 
    double t_bfgs_iter;
 
    public:
    PostTheta(int ns, int nt, int nb, int no, 
        MatrixXd B, Vect y, 
        Vect theta_prior, 
        //Hyperparameters* theta_prior_test,
        Vect mu_initial, 
        string likelihood, Vect extraCoeffVecLik,
        string solver_type,
        const bool constr, const MatrixXd Dxy,
        const bool validate, const Vect w);
    PostTheta(int ns, int nt, int nb, int no, 
        SpMat Ax, Vect y, SpMat c0, SpMat g1, SpMat g2,
        Vect theta_prior, 
        Vect mu_initial, 
        string likelihood, Vect extraCoeffVecLik,
        string solver_type, 
        int dim_spatial_domain, string manifold,
        const bool constr, const MatrixXd Dx, const MatrixXd Dxy,
        const bool validate, const Vect w);
 
    PostTheta(int ns, int nt, int nb, int no, 
        SpMat Ax, Vect y, SpMat c0, SpMat g1, SpMat g2, SpMat g3, 
        SpMat M0, SpMat M1, SpMat M2, 
        Vect theta_prior, 
        Vect mu_initial, 
        string likelihood, Vect extraCoeffVecLik,
        string solver_type, 
        int dim_spatial_domain, string manifold,
        const bool constr, const MatrixXd Dx, const MatrixXd Dxy,
        const bool validate, const Vect w); 
 
    PostTheta(int ns, int nt, int nss, int nb, int no, 
        SpMat Ax, Vect y, SpMat c0, SpMat g1, SpMat g2, SpMat g3, 
        SpMat M0, SpMat M1, SpMat M2, 
        Vect theta_prior_param, 
        Vect mu_initial, 
        string likelihood, Vect extraCoeffVecLik,
        string solver_type, 
        int dim_spatial_domain, string manifold,
        const bool constr, const MatrixXd Dx, const MatrixXd Dxy, 
        const bool validate, const Vect w);
 
 
    double operator()(Vect& theta, Vect& grad);
 
#ifdef DATA_SYNTHETIC
    double compute_error_bfgs(Vect& theta);
#endif
 
    //Hyperparameters create_hp(Vect param, char scale);
 
    void computeG(Vect& theta);
 
    int get_fct_count();
 
    // ============================================================================================ //
    // CONVERT MODEL PARAMETRISATION TO INTERPRETABLE PARAMETRISATION & VICE VERSA
 
    void convert_theta2interpret(Vect& theta, Vect& theta_interpret);
 
    void convert_interpret2theta(Vect& theta_interpret, Vect& theta);
    
    void convert_theta2interpret_spatTemp(double lgamE, double lgamS, double lgamT, double& sigU, double& ranS, double& ranT);
    
    void convert_interpret2theta_spatTemp(double sigU, double ranS, double ranT, double& lgamE, double& lgamS, double& lgamT);
 
    void convert_interpret2theta_spat(double lranS, double lsigU, double& lgamS, double& lgamE);
 
    void convert_theta2interpret_spat(double lgamS, double lgamE, double& lranS, double& lsigU);
 
 
    // ============================================================================================ //
    // FUNCTIONS TO BE CALLED AFTER THE BFGS SOLVER CONVERGED
 
    void get_mu(Vect& theta, Vect& mu_);
 
    #if 0
    void get_mu(Vect& mu);
    #endif
 
    Vect get_grad();
 
    void get_Qprior(Vect theta, SpMat& Qprior);
 
    MatrixXd get_Covariance(Vect theta, double eps);
 
    MatrixXd get_Cov_interpret_param(Vect interpret_theta, double eps);
 
    double f_eval(Vect& theta);
 
    void get_marginals_f(Vect& theta, Vect& mu, Vect& vars);
 
    void get_fullFact_marginals_f(Vect& theta, SpMat& Qinv);
 
    void compute_fullInverseQ(Vect& theta, MatrixXd& Qinv);
 
    MatrixXd hess_eval(Vect theta, double eps);
 
    MatrixXd hess_eval_interpret_theta(Vect interpret_theta, double eps);
 
     void check_pos_def(MatrixXd &hess);
 
    // ============================================================================================ //
    // ALL FOLLOWING FUNCTIONS CONTRIBUT TO THE EVALUATION OF F(THETA) & GRADIENT
 
    double eval_post_theta(Vect& theta, Vect& mu);
 
    void eval_log_gaussian_prior_hp(Vect& theta_param, Vect& theta_prior_param, double& log_prior);
 
    void eval_log_pc_prior_hp(double& log_sum, Vect& lambda, Vect& theta_interpret);
 
    void update_mean_constr(const MatrixXd& D, Vect& e, Vect& sol, MatrixXd& V, MatrixXd& W, MatrixXd& U, Vect& updated_sol);
 
    void eval_log_dens_constr(Vect& x, Vect& mu, SpMat&Q, double& log_det_Q, const MatrixXd& D, MatrixXd& W, double& val_log_dens);
 
    void eval_log_prior_lat(Vect& theta, Vect& mu, double &val);
 
    void eval_likelihood(Vect& theta, Vect& mu, double &log_det, double &val);
    
    void construct_Q_spatial(Vect& theta, SpMat& Qs);
 
    void construct_Q_spat_temp(Vect& theta, SpMat& Qst);
 
    void construct_Qprior(Vect& theta, SpMat& Qx);
 
    void construct_Q(Vect& theta, Vect& mu, SpMat& Q);
 
    void construct_b(Vect& theta, Vect &rhs);
 
    //void update_mean_constr(MatrixXd& D, Vect& e, Vect& sol, MatrixXd& V, MatrixXd& W);
 
    void eval_denominator(Vect& theta, double& val, SpMat& Q, Vect& rhs, Vect& mu);
 
    // ============================================================================================ //
    // INNER ITERATION & everything that is needed for it
 
    double cond_LogPriorLat(SpMat& Qprior, Vect& x);
 
    double cond_LogPoisLik(Vect& eta);
 
    double cond_negLogPoisLik(Vect& eta);
 
    Vect grad_cond_negLogPoisLik(Vect& eta);
 
    Vect diagHess_cond_negLogPoisLik(Vect& eta);
 
 
    double cond_negLogPois(SpMat& Qprior, Vect& x);
 
    void link_f_sigmoid(Vect& x, Vect& sigmoidX);
 
    double cond_negLogBinomLik(Vect& eta);
 
    double cond_negLogBinom(SpMat& Qprior, Vect& x);
 
    double cond_negLogDist(SpMat &Qprior, Vect& x, function<double(Vect&, Vect&)> lik_func);
 
    void FD_gradient(Vect& eta, Vect& grad);
 
    void FD_diag_hessian(Vect& eta, Vect& diag_hess);
 
    void NewtonIter(Vect& theta, Vect& x, SpMat& Q, double& log_det);
 
 
    // measure times within each iterationstd::string file_name, int& iter_count, double& t_Ftheta_ext, double& t_priorHyp, 
    void record_times(std::string file_name, int iter_count, double t_Ftheta_ext, double t_thread_nom, double t_priorHyp, 
                                double t_priorLat, double t_priorLatAMat, double t_priorLatChol, double t_likel, 
                                double t_thread_denom, double t_condLat, double t_condLatAMat, double t_condLatChol, double t_condLatSolve);
 
 
    ~PostTheta();
 
};
 
#endif