// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace Rcpp;

/*
## Introduction of the function ##
# The 'fn_GetBeta()' function computes the values of parameters that are identifiable in a case-control study without any additional assumption.

## Input Variables ##
# rctr/rcas: the control/case cell probabilities (Order {GE: 00, 01, 10, 11}). They can also be control/case cell counts if necessary.

## Output ##
# a vector of (\beta_{e}, \beta_{g}, \beta_{eg})
*/

// [[Rcpp::export]]
NumericVector fn_GetBeta(NumericVector rctr, NumericVector rcas)
{
  rctr = rctr/sum(rctr);
  rcas = rcas/sum(rcas);
  
  NumericVector lrctr = log(rctr), lrcas = log(rcas);
  
  double beta0 = lrcas[0]-lrctr[0];
  double beta1 = lrcas[1]-lrctr[1]-beta0;
  double beta2 = lrcas[2]-lrctr[2]-beta0;
  double beta3 = lrcas[3]-lrctr[3]-beta0-beta1-beta2;
  
  return NumericVector::create( Named("beta1") = beta1, 
                                Named("beta2") = beta2,
                                Named("beta3") = beta3);
}

/*
## Introduction of the function ##
# The 'fn_LambdaToTheta()' function computes the value(s) of \theta given the case-control cell probabilities and the G-E odds ratio in the source population. 

## Input Variables ##
# rctr/rcas: the control/case cell probabilities (Order {GE: 00, 01, 10, 11}). They can also be control/case cell counts if necessary.
# rho: the GE odds ratio in the source population, with the default being 1 (the GEI assumption)

## Output ##
# a list containing zero, one, or two values of \theta. 
*/
// [[Rcpp::export]]
List fn_LambdaToTheta(NumericVector rctr, NumericVector rcas, double rho=1)
{
  NumericVector rdif;

  rctr = rctr/sum(rctr);
  rcas = rcas/sum(rcas);
  rdif = rcas-rctr;
  
  double a = rdif[0]*rdif[3]-rho*rdif[1]*rdif[2];
  double b = rdif[0]*rctr[3]-rho*rdif[1]*rctr[2]+rdif[3]*rctr[0]-rho*rdif[2]*rctr[1];
  double c = rctr[0]*rctr[3]-rho*rctr[1]*rctr[2];
  double d = b*b-4*a*c;
  
  List ThetaList = List::create();
  
  if(a==0) 
  { 
    double theta = -c/b; 
    if((0<theta)&&(theta<1)) { ThetaList.push_back(theta); }
  }  
  else
  {
    if(d==0) 
    { 
      double theta = -b/(2*a); 
      if((0<theta)&&(theta<1)) { ThetaList.push_back(theta); }
    }  
    if(d>0)  
    { 
      double theta;
    
      theta = (-b+sqrt(d))/(2*a); 
      if((0<theta)&&(theta<1)) { ThetaList.push_back(theta); }
      
      theta = (-b-sqrt(d))/(2*a);
      if((0<theta)&&(theta<1)) { ThetaList.push_back(theta); } 
    }
  }
  
  return ThetaList;
}

/*
## Introduction of the function ##
# The 'fn_LambdaToPhi()' determines the value(s) of \phi given the case-control cell probabilities and the G-E odds ratio in the source population

## Input Variables ##
# rctr/rcas: the control/case cell probabilities (Order {GE: 00, 01, 10, 11}). They can also be control/case cell counts if necessary.
# rho: the GE odds ratio in the source population, with the default being 1 (the GEI assumption)

## Output ##
# a list containing zero, one, or two values of \phi
*/
// [[Rcpp::export]]
List fn_LambdaToPhi(NumericVector rctr, NumericVector rcas, double rho=1)
{
  rctr = rctr/sum(rctr);
  rcas = rcas/sum(rcas);
  
  NumericVector rdif = rcas-rctr;
  
  double a = rdif[0]*rdif[3]-rho*rdif[1]*rdif[2];
  double b = rdif[0]*rctr[3]-rho*rdif[1]*rctr[2]+rdif[3]*rctr[0]-rho*rdif[2]*rctr[1];
  double c = rctr[0]*rctr[3]-rho*rctr[1]*rctr[2];
  double d = b*b-4*a*c;
  
  
  List PhiList = List::create();
  
  if(a==0) 
  { 
    double theta = -c/b; 
    if((0<theta)&&(theta<1)) 
    {
      NumericVector pctr = rctr*(1-theta), pcas = rcas*theta;
      NumericVector lpctr = log(pctr), lpcas = log(pcas);
      
      double prvle = pctr[1]+pctr[3]+pcas[1]+pcas[3];
      double prvlg = pctr[2]+pctr[3]+pcas[2]+pcas[3];
      double beta0 = lpcas[0]-lpctr[0];
      double beta1 = lpcas[1]-lpctr[1]-beta0;
      double beta2 = lpcas[2]-lpctr[2]-beta0;
      double beta3 = lpcas[3]-lpctr[3]-beta0-beta1-beta2;
      
      NumericVector Phi = NumericVector::create(
                             Named("beta0") = beta0,
                             Named("beta1") = beta1,
                             Named("beta2") = beta2,
                             Named("beta3") = beta3,
                             Named("prvle") = prvle,
                             Named("prvlg") = prvlg,
                             Named("rho"  ) = rho);
                             
      PhiList.push_back(Phi);
    }  
  }  
  else
  {
    if(d==0) 
    { 
      double theta = -b/(2*a); 
      if((0<theta)&&(theta<1)) 
      {
        NumericVector pctr = rctr*(1-theta), pcas = rcas*theta;
        NumericVector lpctr = log(pctr), lpcas = log(pcas);
        
        double prvle = pctr[1]+pctr[3]+pcas[1]+pcas[3];
        double prvlg = pctr[2]+pctr[3]+pcas[2]+pcas[3];
        double beta0 = lpcas[0]-lpctr[0];
        double beta1 = lpcas[1]-lpctr[1]-beta0;
        double beta2 = lpcas[2]-lpctr[2]-beta0;
        double beta3 = lpcas[3]-lpctr[3]-beta0-beta1-beta2;
        
        NumericVector Phi = NumericVector::create(
                               Named("beta0") = beta0,
                               Named("beta1") = beta1,
                               Named("beta2") = beta2,
                               Named("beta3") = beta3,
                               Named("prvle") = prvle,
                               Named("prvlg") = prvlg,
                               Named("rho"  ) = rho);
                               
        PhiList.push_back(Phi);
      }  
    }  
    if(d>0)  
    { 
      double theta;
      
      theta = (-b+sqrt(d))/(2*a); 
      if((0<theta)&&(theta<1)) 
      {
        NumericVector pctr = rctr*(1-theta), pcas = rcas*theta;
        NumericVector lpctr = log(pctr), lpcas = log(pcas);
        
        double prvle = pctr[1]+pctr[3]+pcas[1]+pcas[3];
        double prvlg = pctr[2]+pctr[3]+pcas[2]+pcas[3];
        double beta0 = lpcas[0]-lpctr[0];
        double beta1 = lpcas[1]-lpctr[1]-beta0;
        double beta2 = lpcas[2]-lpctr[2]-beta0;
        double beta3 = lpcas[3]-lpctr[3]-beta0-beta1-beta2;
        
        NumericVector Phi = NumericVector::create(
                               Named("beta0") = beta0,
                               Named("beta1") = beta1,
                               Named("beta2") = beta2,
                               Named("beta3") = beta3,
                               Named("prvle") = prvle,
                               Named("prvlg") = prvlg,
                               Named("rho"  ) = rho);
                               
        PhiList.push_back(Phi);
      } 
      
      theta = (-b-sqrt(d))/(2*a);
      if((0<theta)&&(theta<1)) 
      {
        NumericVector pctr = rctr*(1-theta), pcas = rcas*theta;
        NumericVector lpctr = log(pctr), lpcas = log(pcas);
        
        double prvle = pctr[1]+pctr[3]+pcas[1]+pcas[3];
        double prvlg = pctr[2]+pctr[3]+pcas[2]+pcas[3];
        double beta0 = lpcas[0]-lpctr[0];
        double beta1 = lpcas[1]-lpctr[1]-beta0;
        double beta2 = lpcas[2]-lpctr[2]-beta0;
        double beta3 = lpcas[3]-lpctr[3]-beta0-beta1-beta2;
        
        NumericVector Phi = NumericVector::create(
                               Named("beta0") = beta0,
                               Named("beta1") = beta1,
                               Named("beta2") = beta2,
                               Named("beta3") = beta3,
                               Named("prvle") = prvle,
                               Named("prvlg") = prvlg,
                               Named("rho"  ) = rho);
                               
        PhiList.push_back(Phi);
      } 
    }
  }
  
  return PhiList;
}

/*    
## Introduction of the function ##
# The 'fn_PhiToCellprb()' function performs the transformation from \phi to cell probabilities (Pr[D,G,E]).

## Input Variables ##
# Betas: the vector of (\beta_{0}, \beta_{e}, \beta_{g}, \beta_{eg})
# prvle: the prevalence of environmental exposure Pr(E=1)
# prvlg: the prevalence of genetic variant Pr(G=1)
# rho:   the GE odds ratio in the source population, with the default of 1 (GEI assumption)

## Output ##
# a list with two elements:
#     pctr: the vector of the control cell probabilities, Pr[G,E,D=0].
#     pacs: the vector of the case cell probabilities, Pr[G,E,D=1].
*/
// [[Rcpp::export]]
List fn_PhiToCellprb(NumericVector Betas, double prvle, double prvlg, double rho=1)
{
  NumericVector pctr(4), pcas(4);

  double intTerm = sqrt(pow((1+(prvle+prvlg)*(rho-1)),2)-4*rho*(rho-1)*prvle*prvlg);
  double addTerm = (4*rho*prvle*prvlg-(rho-1)*pow((prvle+prvlg),2)-2*(prvle+prvlg))/(2+2*intTerm);
  double logitPr0 = Betas[0];
  double logitPr1 = Betas[0]+Betas[1];
  double logitPr2 = Betas[0]+Betas[2];
  double logitPr3 = Betas[0]+Betas[1]+Betas[2]+Betas[3];
    
  pctr[0] = (1+addTerm-(prvle+prvlg)/2)*(1/(1+exp(logitPr0)));
  pctr[1] = ( -addTerm+(prvle-prvlg)/2)*(1/(1+exp(logitPr1)));
  pctr[2] = ( -addTerm-(prvle-prvlg)/2)*(1/(1+exp(logitPr2)));
  pctr[3] = (  addTerm+(prvle+prvlg)/2)*(1/(1+exp(logitPr3)));
  
  pcas[0] = (1+addTerm-(prvle+prvlg)/2)*(exp(logitPr0)/(1+exp(logitPr0)));
  pcas[1] = ( -addTerm+(prvle-prvlg)/2)*(exp(logitPr1)/(1+exp(logitPr1)));
  pcas[2] = ( -addTerm-(prvle-prvlg)/2)*(exp(logitPr2)/(1+exp(logitPr2)));
  pcas[3] = (  addTerm+(prvle+prvlg)/2)*(exp(logitPr3)/(1+exp(logitPr3)));

  return List::create( Named("pctr") = pctr, Named("pcas") = pcas);
}

/*
## Introduction of the function ##
# The 'fn_LogCondtPri()' function computes the logarithm of the prior density for \phi|\rho, which is assumed to be independent of \rho in the manuscript.

## Input Variables ##
# Betas: the vector of (\beta_{0}, \beta_{e}, \beta_{g}, \beta_{eg})
# prvle: the prevalence of environmental exposure Pr(E=1)
# prvlg: the prevalence of genetic variant Pr(G=1)
# hypr1: the hyper-parameter, \sigma_{e}, for the prior of \beta_{e}
# hypr2: the hyper-parameter, \sigma_{g}, for the prior of \beta_{g}
# hypr3: the hyper-parameter, \sigma_{eg}, for the prior of \beta_{eg}

## Output ##
# log of the prior density for \phi|\rho
*/
// [[Rcpp::export]]
double fn_LogCondtPri(NumericVector Betas, double prvle, double prvlg,
                          double hypr1, double hypr2, double hypr3)
{
  double dns = ((exp(Betas[0]))/(pow((1+exp(Betas[0])),2)))*\
               (exp(-pow((Betas[1]/hypr1), 2)/2)/(hypr1*sqrt(2*M_PI)))*\
               (exp(-pow((Betas[2]/hypr2), 2)/2)/(hypr2*sqrt(2*M_PI)))*\
               (exp(-pow((Betas[3]/hypr3), 2)/2)/(hypr3*sqrt(2*M_PI)));
               
  return log(dns);
}

/*
## Introduction of the function ##
# The 'fn_LogDetJacobian()' function computes log of the absolute value of the Jacobian determinant as given in Appendix A of the manuscript.

## Input Variables ##
# Betas: the vector of (\beta_{0}, \beta_{e}, \beta_{g}, \beta_{eg})
# prvle: the prevalence of environmental exposure Pr(E=1)
# prvlg: the prevalence of genetic variant Pr(G=1)
# rho:   the GE odds ratio in the source population, with the default of 1 (GEI assumption)

## Output ##
# log of the absolute value of the Jacobian determinant
*/
// [[Rcpp::export]]
double fn_LogDetJacobian(NumericVector Betas, double prvle, double prvlg, double rho=1)
{
  List cellprb = fn_PhiToCellprb(Betas, prvle, prvlg, rho);
  
  NumericVector pctr = cellprb["pctr"];
  NumericVector pcas = cellprb["pcas"];
  
  double theta = sum(pcas);
  NumericVector peg = pctr+pcas;
  
  double d1 = pow(((prvle+prvlg)*(rho-1)+1),2)-4*rho*(rho-1)*prvle*prvlg;
  double d2 = (rho-1)*pow((prvle+prvlg),2)+2*(prvle+prvlg)-4*rho*prvle*prvlg;
  
  double ke = ((rho+1)*prvlg-(rho-1)*prvle-1)/(2*sqrt(d1));
  double kg = ((rho+1)*prvle-(rho-1)*prvlg-1)/(2*sqrt(d1));
  double kr = d2*(d2+2*prvle*prvlg-prvle-prvlg)/(2*sqrt(d1)*pow((1+sqrt(d1)),2))\
              -pow((prvle-prvlg),2)/(2+2*sqrt(d1));
  
  NumericVector u, v, w;
  u = NumericVector::create(-1,1,-1,1);
  v = NumericVector::create(-1,-1,1,1);
  w = NumericVector::create(1,-1,-1,1);
  
  NumericVector DevTheta(7);
  NumericMatrix DevPctr(4,7), DevPcas(4,7);
  
  DevPctr(_,0) = (pctr/peg)*(w*ke+u/2);
  DevPctr(_,1) = (pctr/peg)*(w*kg+v/2);
  DevPctr(_,2) = (pctr/peg)*(w*kr);
  DevPctr(_,3) = -pcas*pctr/peg;
  DevPctr(1,4) = DevPctr(1,3);
  DevPctr(3,4) = DevPctr(3,3);
  DevPctr(2,5) = DevPctr(2,3);
  DevPctr(3,5) = DevPctr(3,3);
  DevPctr(3,6) = DevPctr(3,3);
  
  DevPcas(_,0) = (pcas/peg)*(w*ke+u/2);
  DevPcas(_,1) = (pcas/peg)*(w*kg+v/2);
  DevPcas(_,2) = (pcas/peg)*(w*kr);
  DevPcas(_,3) = pcas*pctr/peg;
  DevPcas(1,4) = DevPcas(1,3);
  DevPcas(3,4) = DevPcas(3,3);
  DevPcas(2,5) = DevPcas(2,3);
  DevPcas(3,5) = DevPcas(3,3);
  DevPcas(3,6) = DevPcas(3,3);
  
  for(int i=0; i<7; i++) { DevTheta[i] = sum(DevPcas(_,i)); } 
  
  NumericMatrix Jcob(7,7);
  
  Jcob(0,_) = (DevPctr(0,_)*(1-theta)+DevTheta*pctr[0])/pow((1-theta),2);
  Jcob(1,_) = (DevPctr(1,_)*(1-theta)+DevTheta*pctr[1])/pow((1-theta),2);
  Jcob(2,_) = (DevPctr(2,_)*(1-theta)+DevTheta*pctr[2])/pow((1-theta),2);
  Jcob(3,_) = (DevPcas(0,_)*   theta -DevTheta*pcas[0])/pow(   theta ,2);
  Jcob(4,_) = (DevPcas(1,_)*   theta -DevTheta*pcas[1])/pow(   theta ,2);
  Jcob(5,_) = (DevPcas(2,_)*   theta -DevTheta*pcas[2])/pow(   theta ,2);
  Jcob(6,_) =                         DevTheta;
  
  return log(std::abs(arma::det(as<arma::mat>(Jcob))));
}

/*
## Introduction of the function ##
# The 'fn_LogInverseProp()' function computes the factor that accounts for the transformation from (\lambda, \rho) to \xi = (\lambda, \theta), on a logarithmic scale

## Input Variables ##
# rctr/rcas: the control/case cell probabilities (Order {GE: 00, 01, 10, 11}). 
# rho: the GE odds ratio in the source population, with the default being 1 (the GEI assumption)

## Output ##
# the factor that accounts for the transformation from (\lambda, \rho) to \xi = (\lambda, \theta), on a logarithmic scale
*/
// [[Rcpp::export]]
double fn_LogInverseProp(NumericVector rctr, NumericVector rcas, double rho=1)
{
  NumericVector rdif;

  rctr = rctr/sum(rctr);
  rcas = rcas/sum(rcas);
  rdif = rcas-rctr;
  
  double InvProp = 0;
  List ThetaList = fn_LambdaToTheta(rctr, rcas, rho);
  int  nsol = ThetaList.size();
  
  for(int i=0; i<nsol; i++)
  {
    NumericVector Theta = ThetaList[i];
    double theta = Theta[0];
    
    double tm1 = (2*rdif[0]*rdif[3]*theta+rdif[0]*rctr[3]+rctr[0]*rdif[3])*(rdif[1]*theta+rctr[1])*(rdif[2]*theta+rctr[2]);
    double tm2 = (2*rdif[1]*rdif[2]*theta+rdif[1]*rctr[2]+rctr[1]*rdif[2])*(rdif[0]*theta+rctr[0])*(rdif[3]*theta+rctr[3]);
    double div = pow(((rdif[1]*theta+rctr[1])*(rdif[2]*theta+rctr[2])),2);
    
    InvProp = InvProp+std::abs(div/(tm1-tm2));
  }
  
  return log(InvProp);
}
  

/*#################################################
  ###                MAIN FUNCTIONS             ###
  #################################################*/

/*
## Introduction of the function ##
# The 'fn_MainGeneral()' function samples from the general posterior distribution of \beta_{eg}, which assumes a non-informative conjugate prior.

## Input Variables ##
# dctr/dcas: the control/case cell counts. (Order {GE: 00, 01, 10, 11})
# n: the size of the sample from the general posterior distributon

## Output ##
# a list of two elements:
#     beta3GEN: a sample of size n from the general posterior distribution of \beta_{eg}
#     propNonGEI: the proportion of mass that the general posterior puts outside the GEI space
*/
// [[Rcpp::export]]
List fn_MainGeneral(NumericVector dctr, NumericVector dcas, int n)
{
  NumericVector beta3GEN(n), nonGEI(n);
  
  for(int i=0; i<n; i++)
  {
    NumericVector mdctr(4), mdcas(4); 
    
    mdctr[0] = R::rgamma(dctr[0]+1,1);
    mdctr[1] = R::rgamma(dctr[1]+1,1);
    mdctr[2] = R::rgamma(dctr[2]+1,1);
    mdctr[3] = R::rgamma(dctr[3]+1,1);
    mdctr    = mdctr/sum(mdctr);

    mdcas[0] = R::rgamma(dcas[0]+1,1);
    mdcas[1] = R::rgamma(dcas[1]+1,1);
    mdcas[2] = R::rgamma(dcas[2]+1,1);
    mdcas[3] = R::rgamma(dcas[3]+1,1);
    mdcas    = mdcas/sum(mdcas);
    
    beta3GEN[i] = fn_GetBeta(mdctr, mdcas)[2];
    
    List ThetaList = fn_LambdaToTheta(mdctr, mdcas, 1);
    if( ThetaList.size() == 0 ) { nonGEI[i] = 1; } else { nonGEI[i] = 0; } 
  }

  double propNonGEI = mean(nonGEI);
  return List::create( Named("beta3GEN") = beta3GEN,
                       Named("propNonGEI") = propNonGEI );
}

/*
## Introduction of the function ##
# The 'fn_VectorExpand()' function is an auxiliary function that appends some zeros at the end of the input vector.

## Input Variables ##
# x: the input vector
# n: number of zeros to be appended

## Output ##
# a vector of (x, z), where z is vector of n zeros
*/
// [[Rcpp::export]]
NumericVector fn_VectorExpand(const NumericVector& x, int m)
{
  int n = x.size();
  NumericVector y(n+m);
  for( int i=0; i<n; i++) y[i] = x[i];
  return y;
}

/*
## Introduction of the function ##
# The 'fn_VectorExpand()' function is an auxiliary function that appends some rows of zeros at the end of the input matrix.

## Input Variables ##
# X: the input matrix
# mrows: number of rows to be appended

## Output ##
# a matrix of rbind(X, Z), where Z is a matrix of mrows*ncol(X) with all zeros
*/
// [[Rcpp::export]]
NumericMatrix fn_MatrixExpand(const NumericMatrix& X, int mrows )
{
  int nrows = X.nrow(), ncols = X.ncol();
  NumericMatrix Y(nrows+mrows, ncols);
  for(int i=0; i<nrows; i++) { for(int j=0; j<ncols; j++) Y(i,j) = X(i,j); }
  return Y;
}

/*
## Introduction of the function ##
# The 'fn_MainRGEI()' function obtains an importance weighted sample to numerically represent the RGEI posterior distribution of \phi.
# Note that the GEI posterior distribution is just a special case where \sigma_{\rho} = 0.

## Input Variables ##
# dctr/dcas: control/case cell counts
# ESS: the desired effective sample size of the importance sample, with the default being 10000
# hypr1: the hyper-parameter, \sigma_{e}, for the prior of \beta_{e}
# hypr2: the hyper-parameter, \sigma_{g}, for the prior of \beta_{g}
# hypr3: the hyper-parameter, \sigma_{eg}, for the prior of \beta_{eg}
# hyprr: the hyper-parameter, \sigma_{\rho}, for the prior of \rho, with the default being 0 (GEI)

## Output ##
# a list of two elements:
#     Phi: A 7-column matrix, with each row being a sample point of \phi
#     Weight: the importance weight associated with the corresponding sample point
*/
// [[Rcpp::export]]
List fn_MainRGEI(NumericVector dctr, NumericVector dcas, int ESS=10000,
                 double hypr1=2, double hypr2=2, double hypr3=2, double hyprr=0)
{  
  int nlp = 0;
  double ess;

  NumericVector LogWeight(5000);
  NumericMatrix MatrixPhi(5000, 7);
  
  for(int i=0; i<5000; i++)
  {
    int nsol;
    double rho;
    NumericVector rctr(4), rcas(4);
    List PhiList;

    do {
      rho  = exp(R::rnorm(0, hyprr));
  
      rctr[0] = R::rgamma(dctr[0]+1,1);
      rctr[1] = R::rgamma(dctr[1]+1,1);
      rctr[2] = R::rgamma(dctr[2]+1,1);
      rctr[3] = R::rgamma(dctr[3]+1,1);
      rctr    = rctr/sum(rctr);
      
      rcas[0] = R::rgamma(dcas[0]+1,1);
      rcas[1] = R::rgamma(dcas[1]+1,1);
      rcas[2] = R::rgamma(dcas[2]+1,1);
      rcas[3] = R::rgamma(dcas[3]+1,1);
      rcas    = rcas/sum(rcas);
      
      PhiList = fn_LambdaToPhi(rctr, rcas, rho);
      nsol = PhiList.size(); 
      
    } while(nsol==0);
  
    int index = 0;
    if(nsol==2)
    {
      double tmp = R::runif(0,1);
      if(tmp > 0.5) index = 1;
    }
    
    NumericVector Phi = PhiList[index];
    
    for(int j=0; j<7; j++) MatrixPhi(nlp,j) = Phi[j];  
    
    NumericVector Betas(4);
    Betas[0] = Phi[0];
    Betas[1] = Phi[1];
    Betas[2] = Phi[2];
    Betas[3] = Phi[3];
    double prvle = Phi[4];
    double prvlg = Phi[5];
    
    LogWeight[nlp] = fn_LogCondtPri(Betas, prvle, prvlg, hypr1, hypr2, hypr3)-\
                     fn_LogDetJacobian(Betas, prvle, prvlg, rho)+\
                     fn_LogInverseProp(rctr, rcas, rho);  

    nlp++;
  }

  {
    NumericVector rWeight = exp(LogWeight-max(LogWeight)); 
    NumericVector  Weight = rWeight/sum(rWeight);
    ess = pow(sum(Weight),2)/sum(pow(Weight,2));
  }

  // dynamically increase the sample size until reach the desired effective sample size
  while (ess < ESS) 
  { 
    LogWeight = fn_VectorExpand(LogWeight, 5000);
    MatrixPhi = fn_MatrixExpand(MatrixPhi, 5000);
    
    for(int i=0; i<5000; i++)
    {
      int nsol;
      double rho;
      NumericVector rctr(4), rcas(4);
      List PhiList;
      
      do {
        rho  = exp(R::rnorm(0, hyprr));
  
        rctr[0] = R::rgamma(dctr[0]+1,1);
        rctr[1] = R::rgamma(dctr[1]+1,1);
        rctr[2] = R::rgamma(dctr[2]+1,1);
        rctr[3] = R::rgamma(dctr[3]+1,1);
        rctr    = rctr/sum(rctr);
        
        rcas[0] = R::rgamma(dcas[0]+1,1);
        rcas[1] = R::rgamma(dcas[1]+1,1);
        rcas[2] = R::rgamma(dcas[2]+1,1);
        rcas[3] = R::rgamma(dcas[3]+1,1);
        rcas    = rcas/sum(rcas);
        
        PhiList = fn_LambdaToPhi(rctr, rcas, rho);
        nsol = PhiList.size();
        
      } while(nsol==0);

      int index = 0;
      if(nsol==2)
      {
        double tmp = R::runif(0,1);
        if(tmp > 0.5) index = 1;
      }
      
      NumericVector Phi = PhiList[index];
      
      for(int j=0; j<7; j++) MatrixPhi(nlp,j) = Phi[j];
      
      NumericVector Betas(4);
      Betas[0] = Phi[0];
      Betas[1] = Phi[1];
      Betas[2] = Phi[2];
      Betas[3] = Phi[3];
      double prvle = Phi[4];
      double prvlg = Phi[5];

      LogWeight[nlp] = fn_LogCondtPri(Betas, prvle, prvlg, hypr1, hypr2, hypr3)-\
                       fn_LogDetJacobian(Betas, prvle, prvlg, rho)+\
                       fn_LogInverseProp(rctr, rcas, rho);
                       
      nlp++;

    } 
    
    NumericVector rWeight = exp(LogWeight-max(LogWeight));
    NumericVector  Weight = rWeight/sum(rWeight);
    ess = pow(sum(Weight),2)/sum(pow(Weight,2));

  }

  NumericVector rWeight = exp(LogWeight-max(LogWeight));
  NumericVector  Weight = rWeight/sum(rWeight);

  return List::create(
               Named("Phi") = MatrixPhi,
               Named("Weight") = Weight);
}