/***************************************************************************
 *   Copyright (C) 2007 by Faubet Pierre   *
 *   pierre.faubet@e.ujf-grenoble.fr   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program 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 General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

// #include <QWriteLocker>

#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_randist.h>

#include <iostream>
#include <sstream>
#include <time.h>

#include "mcmc.h"
#include "BIMrErrors.h"

// QReadWriteLock MCMC::lock ;

/* Static member initialization */
int MCMC::nb_mcmc = 0;

/* MCMC settings */
int MCMC::nb_pilot = 20;
int MCMC::pilot_length = 1000;
int MCMC::burnin = 10000;
int MCMC::samplesize = 10000;
int MCMC::thining = 100;

/* Enable MCMC options */
bool MCMC::useReg = true;
bool MCMC::useRJMCMC = true;
bool MCMC::useFmodel = true;
bool MCMC::useinter=true;

/* Ouput options */
bool MCMC::printLglik = true;
bool MCMC::printPsi = true;
bool MCMC::printMig = true;
bool MCMC::printReg = true;
bool MCMC::printP = true;
bool MCMC::printF = true;
bool MCMC::printFst = true;
bool MCMC::printPglob = true;
bool MCMC::printSd = true;
bool MCMC::printpopq = true;
bool MCMC::printindivq = true;
bool MCMC::printmtild = true;

/* Priors */
double MCMC::atau = 1;
double MCMC::btau = 1;
double MCMC::s2_alpha = 4;
double MCMC::Psi = 1;
double MCMC::lambda = 1;
double MCMC::omega = 1;
double MCMC::xi = 1;

/* Proposal */
double MCMC::e_F_init = 0.1;
double MCMC::e_p_init = 0.1;
double MCMC::e_m_init = 0.1;
double MCMC::e_nm_init = 0.1;
double MCMC::s2_psi_init = 5;
double MCMC::s2_theta_init = 5;

/* Data */
multilocus_genotypes *MCMC::genodata = NULL;
environmental_factors *MCMC::factors = NULL;

/* Input files */
string MCMC::MGinputfile=string();
string MCMC::EFinputfile=string();

MCMC::MCMC ()
{

  /* Create run */
  nb_mcmc++;
  idf = nb_mcmc;
  clog << "Create MCMC #" << idf << endl;

  stringstream ss;
  string suffix;
  ss << idf;
  ss >> suffix;

  if (!genodata)
    {
      cerr << "Multilocus genotypes missing" << endl;
      EXCEPTION_INFOS();
      throw MCMCError();
    }

  fileprefix = genodata->get_inputfilename() + "_run" + suffix;

  /* Result file */
  string resultfilename = fileprefix +".out";
  resultfile.open( resultfilename.data ());
  if (resultfile.fail())
    {
      EXCEPTION_INFOS();
      throw OpenFileError();
    }

  /* Random generator */
  time_t ti = time ( NULL );
  struct tm *tb = localtime ( &ti );
  unsigned long init[5] = { idf, tb->tm_yday, tb->tm_sec, tb->tm_min, tb->tm_hour};
  int length = 5;

  resultfile << "Chain #" << idf << endl;
  resultfile << "Seed for random number generator:";
  unsigned long myseed=0;
  for ( int s = 0; s < length; s++ )
    {
      myseed += init[s];
      resultfile << " " << init[s];
    }
  resultfile << endl;

  const gsl_rng_type * rngT;
  gsl_rng_env_setup();
  rngT = gsl_rng_default;
  rng = gsl_rng_alloc (rngT);
  gsl_rng_set (rng,myseed);

  /* Initial incremental values */
  e_F = e_F_init;
  e_p = e_p_init;
  e_m = e_m_init;
  e_nm = e_nm_init;
  s2_psi = s2_psi_init;
  s2_theta = s2_theta_init;

  /* Write data into the output file */
  if ( useReg )
    {
      nb_model = factors->get_nb_model();
      visit = new double[nb_model];
      if ( !visit )
        {
          EXCEPTION_INFOS();
          throw MemError();
        }
      for ( int mo = 0; mo < nb_model; mo++ )
        visit[mo] = 0;
    }

  d_assign=0;
  /* Create raw mcmc files */
  create_output_files ();

  clog << "Create MCMC #" << idf << " ... (done)" << endl;
}

MCMC::~MCMC ()
{
  close_outputfiles();

  gsl_rng_free(rng);

  for (int q=0;q<genodata->get_nbpop();q++)
    {
      delete [] m[q];
      delete [] m_mean[q];
      delete [] m_sd[q];
    }
  delete [] m;
  delete [] m_sd;
  delete [] m_mean;

  delete [] F_mean;
  delete [] F_sd;

  for (int q=0;q<genodata->get_nbpop();q++)
    {
      for (int j=0;j<genodata->get_nbloci();j++)
        {
          delete [] p_mean[q][j];
          delete [] p_sd[q][j];
        }
      delete [] p_mean[q];
      delete [] p_sd[q];
    }
  delete [] p_mean;
  delete [] p_sd;

  if (useFmodel)
    {
      for (int j=0;j<genodata->get_nbloci();j++)
        delete [] glob_all_frq[j];
      delete [] glob_all_frq;
      delete [] theta;
    }

  if (useReg)
    {
      for (int q=0;q<genodata->get_nbpop();q++)
        {
          delete [] psi[q];
          delete [] mu[q];
        }
      delete [] psi;
      delete [] mu;
      delete [] alpha;

      delete [] visit;
    }

  if (!genodata->get_nodata())
    {
      for (int q=0;q<genodata->get_nbpop();q++)
        {
          for (int l=0;l<genodata->get_nbpop();l++)
            delete [] ancescount[q][l];
          delete [] ancescount[q];
        }
      delete [] ancescount;

      for (int h=0;h<genodata->get_nbindiv();h++)
        {
          for (int q=0;q<genodata->get_nbpop();q++)
            delete [] post_assign[h][q];
          delete [] post_assign[h];
          delete [] Mt[h];
        }
      delete [] post_assign;
      delete [] Mt;
    }

  delete [] F;
  delete [] freqmin;

  for (int q=0;q<genodata->get_nbpop();q++)
    {
      for (int j=0;j<genodata->get_nbloci();j++)
        delete [] p[q][j];
      delete [] p[q];
    }
  delete [] p;

  delete [] tmp_ind_lglik;
  delete [] ind_lglik;
  delete [] pop_lglik;

  nb_mcmc--;
  if (!nb_mcmc)
    {
      if (!EFinputfile.empty())
        delete factors;
      delete genodata;
    }
}

void MCMC::create_output_files ()
{
  if ( printMig )
    create_mcmc_outputfile ();

  if ( printmtild )
    create_mtild_outputfile ();

  if ( useFmodel )
    {
      if ( printPglob )
        create_frq_outputfile ();
      if ( printFst )
        create_fst_outputfile ();
    }

  if ( printF )
    create_F_outputfile();

  if ( printLglik )
    create_lglk_file ();

  if ( printP )
    create_p_file ();

  if ( useReg && printPsi )
    create_psi_outputfile();

  if ( useReg && printReg )
    create_reg_outputfile ();
}

void MCMC::init ()
{
  clog << "Initialization of MCMC #" << idf << endl;

  if ( useFmodel )
    {
      init_glob_all_frq ();
      init_theta ();
    }

  if ( useReg )
    {
      model = nb_model - 1;
      nb_reg_param = factors->get_size();
      init_alpha ();
      init_s2 ();
      init_psi ();
    }

  init_m ();

  if ( !genodata->get_nodata () )
    init_Mt ();

  init_p ();
  init_F ();
  init_lglik ();

  init_stat ();

  RJ_acc = 0;

  clog << "Initialization of MCMC #" << idf << " ... (done)" << endl;
}

void MCMC::update ()
{
  if ( useFmodel )
    {
      update_glob_all_frq ();
      update_theta ();
    }

  update_p ();
  update_F ();

  if ( !genodata->get_nodata () )
    update_Mt ();

  if ( genodata->get_nbpop() > 2 )
    {
      update_non_migrant_proportions ();
      update_m ();
    }
  else
    update_m_2pop();

  if ( useReg )
    {
      update_psi ();
      update_s2 ();
      lgMV = MVregupdate ();
    }
}

void MCMC::pilot_run()
{
  clog << "Pilot runs for MCMC #" << idf << endl;
  for ( int rep = 0; rep < nb_pilot; rep++ )
    {
      for ( int iter=1;iter<=pilot_length;iter++ )
        update();

      /* Adjust incremental values */
      if ( nm_acc <= 0.25 * pilot_length * genodata->get_nbpop () )
        e_nm /= 1.2;
      if ( nm_acc >= 0.45 * pilot_length * genodata->get_nbpop () )
        e_nm *= 1.2;

      if ( m_acc <= 0.25 * pilot_length * genodata->get_nbpop () )
        e_m /= 1.2;
      if ( m_acc >= 0.45 * pilot_length * genodata->get_nbpop () )
        e_m *= 1.2;

      if ( useReg )
        {
          if ( psi_acc <= 0.25 * pilot_length * ( genodata->get_nbpop () - 1 ) * genodata->get_nbpop () )
            s2_psi /= 1.5;
          if ( psi_acc >= 0.45 * pilot_length * genodata->get_nbpop () * genodata->get_nbpop () )
            s2_psi *= 1.5;
        }

      if ( useFmodel )
        {
          if ( theta_acc <= 0.25 * pilot_length * genodata->get_nbpop () )
            s2_theta /= 1.5;
          if ( theta_acc >= 0.45 * pilot_length * genodata->get_nbpop () )
            s2_theta *= 1.5;
        }

      if ( F_acc <= 0.25 * pilot_length * genodata->get_nbpop () )
        e_F *= 0.5;
      if ( F_acc >= 0.45 * pilot_length * genodata->get_nbpop () )
        e_F *= 1.5;

      /* Print acceptance rates */
//       clog << "Acceptance rates (%)" << endl;
//       if ( useFmodel )
//         {
//           clog << "ptild: " << pglob_acc * 100 / ( pilot_length * genodata->get_nbloci () ) << endl;
//           clog << "theta: " << theta_acc * 100 / ( pilot_length * genodata->get_nbpop () ) << endl;
//           theta_acc = 0;
//           pglob_acc = 0;
//         }
//       clog << "p: " << p_acc * 100 / ( pilot_length * genodata->get_nbpop () * genodata->get_nbloci () ) << endl;
//       clog << "F: " << F_acc * 100 / ( pilot_length * genodata->get_nbpop () ) << endl;
//       clog << "Mt: " << Mt_acc * 100 / ( pilot_length * genodata->get_nbindiv () ) << endl;
//       clog << "m: " << m_acc * 100 / ( pilot_length * genodata->get_nbpop () ) << endl;
//       clog << "nm: " << nm_acc * 100 / ( pilot_length * genodata->get_nbpop () ) << endl;
//       if ( useReg )
//         {
//           clog << "psi: " << psi_acc * 100 / ( pilot_length * genodata->get_nbpop () * ( genodata->get_nbpop () - 1 ) ) << endl;
//           psi_acc = 0;
//         }

      /* Reset acceptance rates */
      nm_acc = 0;
      m_acc = 0;
      p_acc = 0;
      Mt_acc = 0;
      F_acc = 0;
    }
  clog << "Pilot runs for MCMC #" << idf << " ... (done)" << endl;
}

void MCMC::run ()
{
  /* Print proposal */
  write_parameters();

  /* Regression parameter initialization */
  if ( useReg )
    {
      model = nb_model - 1;
      nb_reg_param = factors->get_size();

      for ( int r = 1; r < factors->get_size(); r++ )
        if ( ! ( ( model >> ( r - 1 ) ) % 2 ) )
          alpha[r] = 0;
    }

  int nbiter = burnin + samplesize * thining;
  clog << "Metropolis-Hasting algorithm for MCMC #" << idf << endl;
  int stp = 1;
  for ( int iter = 1; iter <= nbiter; iter++ )
    {
      /* Within model moves */
      update ();

      /* Between model move */
      if ( useRJMCMC )
        RJupdate ();

      /* Print iteration */
      if ( 100 * iter == stp * 5 * nbiter )
        {
          clog << "Acceptance rates for chain #" << idf << " (%)" << endl;
          clog << "Iteration #" << iter << endl;
          if ( useFmodel )
            {
              clog << "ptild: " << pglob_acc * 100 / ( iter * genodata->get_nbloci () ) << endl;
              clog << "theta: " << theta_acc * 100 / ( iter * genodata->get_nbpop () ) << endl;
            }
          clog << "p: " << p_acc * 100 / ( iter * genodata->get_nbpop () * genodata->get_nbloci () ) << endl;
          clog << "F: " << F_acc * 100 / ( iter * genodata->get_nbpop () ) << endl;
          clog << "Mt: " << Mt_acc * 100 / ( iter * genodata->get_nbindiv () ) << endl;
          clog << "m: " << m_acc * 100 / ( iter * genodata->get_nbpop () ) << endl;
          clog << "nm: " << nm_acc * 100 / ( iter * genodata->get_nbpop () ) << endl;
          if ( useReg )
            clog << "psi: " << psi_acc * 100 / ( iter * ( genodata->get_nbpop () - 1 ) * genodata->get_nbpop () ) << endl;
          if ( useRJMCMC )
            clog << "Model: " << RJ_acc * 100 / iter << endl;
          stp++;
        }

      /* Sampling */
      if ( ( iter > burnin ) && ( ! ( ( iter - burnin ) % thining ) ) )
        {
          d_assign += -2*lglikMt;

          /* The number visits in current model */
          if ( useRJMCMC )
            visit[model]++;

          /* Write raw mcmc files */
          write_iter ( iter );

          /* Posterior statistics */
          for ( int h = 0; h < genodata->get_nbindiv (); h++ )
            post_assign[h][Mt[h][1]][Mt[h][0]] = ++post_assign[h][Mt[h][0]][Mt[h][1]];

          for ( int q = 0; q < genodata->get_nbpop (); q++ )
            {
              F_mean[q] += F[q];
              F_sd[q] += F[q] * F[q];

              for ( int l = 0; l < genodata->get_nbpop (); l++ )
                {
                  m_mean[q][l] += m[q][l];
                  m_sd[q][l] += m[q][l] * m[q][l];
                }
              for ( int j = 0; j < genodata->get_nbloci (); j++ )
                for ( int i = 1; i < genodata->get_nballel (); i++ )
                  {
                    p_mean[q][j][i] += p[q][j][i];
                    p_sd[q][j][i] += p[q][j][i] * p[q][j][i];
                  }
            }
        }
    }
  close_outputfiles();
  /* End of MH algorithm */
  clog << "Metropolis-Hasting algorithm for MCMC #" << idf << " ... (done)" << endl;
}

void MCMC::print ()
{
//   QWriteLocker locker ( &lock );

  print_acceptance_rates ( burnin + samplesize * thining );

  poststats ();
}

inline void MCMC::write_iter ( int iteration )
{
//   QWriteLocker locker ( &lock );

  /* Write m and psi */
  if ( printMig )
    {
      mcmcfile << iteration;
      for ( int q = 0; q < genodata->get_nbpop (); q++ )
        for ( int l = 0; l < genodata->get_nbpop (); l++ )
          mcmcfile << '\t' << m[q][l];
      mcmcfile << endl;
    }

  if ( printmtild )
    {
      mtildfile << iteration;
      for ( int q = 0; q < genodata->get_nbpop (); q++ )
        for ( int l = 0; l < genodata->get_nbpop (); l++ )
          for (int t = l;t<genodata->get_nbpop (); t++)
            mtildfile << '\t' << (t == l ? 1 : 2)*m[q][l]*m[q][t];
      mtildfile << endl;
    }

  if ( useReg && printPsi )
    {
      psifile << iteration;
      for ( int q = 0; q < genodata->get_nbpop (); q++ )
        for ( int l = 0; l < genodata->get_nbpop (); l++ )
          psifile << '\t' << psi[q][l];
      psifile << endl;
    }

  /* Write regression parameters */
  if ( useReg && printReg )
    {
      regfile << iteration;
      regfile << '\t' << model;
      regfile << '\t' << s2;
      for ( int r = 0; r < factors->get_size(); r++ )
        regfile << '\t' << alpha[r];
      regfile << endl;
    }

  /* Write loglik */
  if ( printLglik )
    {
      double loglik = lglikgeno + lglikMt;

      // f(m|psi)
      double loglik_m = 0;
      if ( useReg )
        {
          for ( int q = 0; q < genodata->get_nbpop (); q++ )
            {
              loglik_m += - ( genodata->get_nbpop () - 2 ) * log ( 1 - m[q][q] ) + gsl_sf_lngamma ( psi[q][q] );
              for ( int l = 0; l < genodata->get_nbpop (); l++ )
                if ( l != q )
                  loglik_m += ( psi[q][l] - 1 ) * ( log ( m[q][l] ) - log ( 1 - m[q][q] ) ) - gsl_sf_lngamma ( psi[q][l] );
            }
        }
      else
        {
          for ( int q = 0; q < genodata->get_nbpop (); q++ )
            {
              for ( int l = 0; l < genodata->get_nbpop (); l++ )
                {
                  loglik_m += ( Psi - 1 ) * log ( m[q][l] );
                }
              loglik_m += gsl_sf_lngamma ( genodata->get_nbpop() *Psi ) - genodata->get_nbpop() *gsl_sf_lngamma ( Psi );
            }
        }
      loglik += loglik_m;

      // f(psi|alpha,s2,G)
      double loglik_psi = 0;
      if ( useReg )
        for ( int q = 0; q < genodata->get_nbpop (); q++ )
          for ( int l = 0; l < genodata->get_nbpop (); l++ )
            if ( l != q )
              loglik_psi -= log ( psi[q][l] ) + 0.5 * log ( 2 * M_PI * s2 ) + 0.5 * ( log ( psi[q][l] ) - mu[q][l] ) * ( log ( psi[q][l] ) - mu[q][l] ) / s2;
      loglik += loglik_psi;

      // f(alpha)
      double loglik_alpha = 0;
      if ( useReg )
        {
          loglik_alpha -= 0.5 * log ( 2 * M_PI * s2_alpha ) + alpha[0] * alpha[0] / ( 2*s2_alpha );
          for ( int r = 1; r < factors->get_size(); r++ )
            if ( ( model >> ( r - 1 ) ) % 2 )
              loglik_alpha -= 0.5 * log ( 2 * M_PI * s2_alpha ) + alpha[r] * alpha[r] / ( 2*s2_alpha );
        }
      loglik += loglik_alpha;

      // f(s2)
      double loglik_s2 = 0;
      if ( useReg )
        loglik_s2 += atau*log ( btau ) - gsl_sf_lngamma ( atau ) - ( atau+1 ) *log ( s2 ) - btau/s2 ;
      loglik += loglik_s2;

      // f(p|theta,ptild)
      double loglik_p = 0;
      if ( useFmodel )
        for ( int q = 0; q < genodata->get_nbpop (); q++ )
          for ( int j = 0; j < genodata->get_nbloci (); j++ )
            {
              loglik_p += gsl_sf_lngamma ( theta[q] );
              for ( int k = 1; k < genodata->get_nballel (); k++ )
                if ( genodata->get_is_at_locus ( j, k ) && ( p[q][j][k] != 0 ) )
                  loglik_p += ( theta[q] * glob_all_frq[j][k] - 1 ) * log ( p[q][j][k] ) - gsl_sf_lngamma ( theta[q] * glob_all_frq[j][k] );
            }
      loglik += loglik_p;

      // f(theta)
      double loglik_theta = 0;
      if ( useFmodel )
        for ( int q = 0; q < genodata->get_nbpop (); q++ )
          loglik_theta -= log ( theta[q] ) + 0.5 * log ( 2 * M_PI * xi ) + 0.5 * ( log ( theta[q] ) - omega ) * ( log ( theta[q] ) - omega ) / ( xi );
      loglik += loglik_theta;

      loglikfile << iteration;
      loglikfile << '\t' << loglik;
      loglikfile << '\t' << lglikgeno;
      loglikfile << '\t' << lglikMt;
      loglikfile << '\t' << loglik_m;
      if ( useReg )
        {
          loglikfile << '\t' << loglik_psi;
          loglikfile << '\t' << loglik_alpha;
          loglikfile << '\t' << loglik_s2;
        }
      loglikfile << '\t' << loglik_p;
      if ( useFmodel )
        loglikfile << '\t' << loglik_theta;
      loglikfile << endl;
    }

  /* Write F-stat */
  if ( printF )
    {
      Fisfile << iteration;
      for ( int q = 0; q < genodata->get_nbpop (); q++ )
        Fisfile << '\t' << F[q];
      Fisfile << endl;
    }

  if ( useFmodel && printFst )
    {
      fstfile << iteration;
      for ( int q = 0; q < genodata->get_nbpop (); q++ )
        fstfile << '\t' << 1 / ( 1 + theta[q] );
      fstfile << endl;
    }

  /* Write ancestral allele frequencies */
  if ( useFmodel && printPglob )
    {
      frqfile << iteration;
      for ( int j = 0; j < genodata->get_nbloci (); j++ )
        for ( int k = 1; k < genodata->get_nballel (); k++ )
          if ( genodata->get_is_at_locus ( j, k ) )
            frqfile << '\t' << glob_all_frq[j][k];
      frqfile << endl ;
    }

  /* Write population allele frequencies */
  if ( printP )
    {
      pfile << iteration;
      for ( int q = 0; q < genodata->get_nbpop (); q++ )
        for ( int j = 0; j < genodata->get_nbloci (); j++ )
          for ( int k = 1; k < genodata->get_nballel (); k++ )
            if ( genodata->get_is_at_locus ( j, k ) )
              pfile << '\t' << p[q][j][k];
      pfile << endl;
    }
}

inline void MCMC::create_lglk_file ()
{
  string lglkfilename = fileprefix + "_lglik.txt";
  loglikfile.open( lglkfilename.data () );

  loglikfile << "iter";
  loglikfile << '\t' << "Likelihood";
  loglikfile << '\t' << "Pr(X|S;M,F,p)";
  loglikfile << '\t' << "Pr(M|m)";
  loglikfile << '\t' << "f(m|psi)";
  if ( useReg )
    {
      loglikfile << '\t' << "f(psi|alpha,s2,G)";
      loglikfile << '\t' << "f(alpha)";
      loglikfile << '\t' << "f(s2)";
    }
  loglikfile << '\t' << "f(p" << (useFmodel ? "|theta,ptild" : "") << ")";
  if ( useFmodel )
    loglikfile << '\t' << "f(theta)";
  loglikfile << endl;
}

inline void MCMC::create_p_file ()
{
  string pfilename = fileprefix + "_popfrq.txt";
  pfile.open( pfilename.data () );

  pfile << "iter";
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    for ( int j = 0; j < genodata->get_nbloci (); j++ )
      for ( int k = 1; k < genodata->get_nballel (); k++ )
        if ( genodata->get_is_at_locus ( j, k ) )
          pfile << '\t' << "p(" << q + 1 << ',' << j + 1 << ',' << k << ")";
  pfile << endl;
}

inline void MCMC::create_frq_outputfile ()
{
  string frqfilename = fileprefix + "_ancfrq.txt";
  frqfile.open( frqfilename.data () );

  frqfile << "iter";
  for ( int j = 0; j < genodata->get_nbloci (); j++ )
    for ( int k = 1; k < genodata->get_nballel (); k++ )
      if ( genodata->get_is_at_locus ( j, k ) )
        frqfile << '\t' << "p(" << j + 1 << ',' << k << ")";
  frqfile << endl;
}

inline void MCMC::create_F_outputfile ()
{
  string Ffilename = fileprefix + "_F.txt";
  Fisfile.open( Ffilename.data () );

  Fisfile << "iter";
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    Fisfile << '\t' << "F" << q + 1;
  Fisfile << endl;
}

inline void MCMC::create_fst_outputfile ()
{
  string fstfilename = fileprefix + "_Fst.txt";
  fstfile.open( fstfilename.data () );

  fstfile << "iter";
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    fstfile << '\t' << "Fst" << q + 1;
  fstfile << endl;
}

inline void MCMC::create_mcmc_outputfile ()
{
  string mcmcfilename = fileprefix + "_mig.txt";
  mcmcfile.open( mcmcfilename.data () );

  mcmcfile << "iter";
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    for ( int l = 0; l < genodata->get_nbpop (); l++ )
      mcmcfile << '\t' << "m" << q + 1 << l + 1;
  mcmcfile << endl ;
}

inline void MCMC::create_mtild_outputfile ()
{
  string mtildfilename = fileprefix + "_mtild.txt";
  mtildfile.open( mtildfilename.data () );

  mtildfile << "iter";
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    for ( int l = 0; l < genodata->get_nbpop (); l++ )
      for (int t = l;t<genodata->get_nbpop (); t++)
        mtildfile << '\t' << "m" << q + 1 << l + 1 << t+1;
  mtildfile << endl ;
}

inline void MCMC::create_psi_outputfile ()
{
  string psifilename = fileprefix + "_psi.txt";
  psifile.open( psifilename.data () );

  psifile << "iter";
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    for ( int l = 0; l < genodata->get_nbpop (); l++ )
      psifile << '\t' << "psi" << q + 1 << l + 1;
  psifile << endl ;
}

inline void MCMC::create_prop_file ()
{
  string ancfilename = fileprefix + "_assign.txt";
  propfile.open( ancfilename.data () );

  propfile << "iter";
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    for ( int l = 0; l < genodata->get_nbpop (); l++ )
      for ( int t = l; t < genodata->get_nbpop (); t++ )
        propfile << '\t' << "n" << q + 1 << l + 1 << t + 1;
  propfile << endl;
}

inline void MCMC::create_reg_outputfile ()
{
  string regfilename = fileprefix + "_reg.txt";
  regfile.open( regfilename.data () );

  regfile << "iter";
  regfile << '\t' << "Model";
  regfile << '\t' << "Sigma2";
  for ( int r = 0; r < factors->get_size(); r++ )
    regfile << '\t' << "Alpha" << r;
  regfile << endl;
}

inline void MCMC::write_parameters ()
{
  resultfile << "MCMC run paremeters" ;
  resultfile << "Burnin: " << burnin ;
  resultfile << "Sample size: " << samplesize ;
  resultfile << "Thining: " << thining ;
  resultfile << "e_p: " << e_p ;
  resultfile << "e_F: " << e_F ;
  resultfile << "e_m: " << e_m ;
  resultfile << "e_nm: " << e_nm ;
  if ( useReg )
    resultfile << "s2_psi: " << s2_psi ;
  if ( useFmodel )
    resultfile << "s2_theta: " << s2_theta ;
}

inline void MCMC::print_acceptance_rates ( int nbit )
{
  resultfile << endl;
  resultfile << "Acceptance rates for chain #" << idf << " (%)" << endl;
  if ( useFmodel )
    {
      resultfile << "ptild: " << pglob_acc * 100 / ( nbit * genodata->get_nbloci () ) << endl;
      resultfile << "theta: " << theta_acc * 100 / ( nbit * genodata->get_nbpop () ) << endl;
    }
  resultfile << "p: " << p_acc * 100 / ( nbit * genodata->get_nbpop () * genodata->get_nbloci () ) << endl;
  resultfile << "F: " << F_acc * 100 / ( nbit * genodata->get_nbpop () ) << endl;
  resultfile << "Mt: " << Mt_acc * 100 / ( nbit * genodata->get_nbindiv () ) << endl ;
  resultfile << "nm: " << nm_acc * 100 / ( nbit * genodata->get_nbpop () ) << endl;
  resultfile << "m: " << m_acc * 100 / ( nbit * genodata->get_nbpop () ) << endl;
  if ( useReg )
    {
      resultfile << "psi: " << psi_acc * 100 / ( nbit * ( genodata->get_nbpop () - 1 ) * genodata->get_nbpop () ) << endl;
    }
  if ( useRJMCMC )
    resultfile << "Model: " << RJ_acc * 100 / nbit << endl;
}


inline void MCMC::poststats ()
{
  if ( useRJMCMC )
    {
      resultfile << endl;
      resultfile << "Posterior model probabilities" << endl;
      for ( model = 0; model < nb_model; model++ )
        if ( factors->is_a_valid_model (model) )
          resultfile << "#" << model << '\t' << visit[model] * 100 / samplesize << endl;
    }
  resultfile << endl;
  resultfile << "Migration rates" << endl;
  resultfile << "Into/From";
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    resultfile << '\t' << q + 1;
  resultfile << endl;

  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    {
      resultfile << q + 1;
      for ( int l = 0; l < genodata->get_nbpop (); l++ )
        {
          m_mean[q][l] /= ( ( double ) samplesize );
          m_sd[q][l] = sqrt ( m_sd[q][l] / ( ( double ) samplesize ) - m_mean[q][l] * m_mean[q][l] );
          resultfile << '\t' << m_mean[q][l];
        }
      resultfile << endl;
    }

  if ( printSd )
    {
      resultfile << "Standard deviation"<< endl ;
      for ( int q = 0; q < genodata->get_nbpop (); q++ )
        {
          resultfile << q + 1;
          for ( int l = 0; l < genodata->get_nbpop (); l++ )
            resultfile << '\t' << m_sd[q][l];
          resultfile << endl;
        }
    }
  resultfile << endl;
  resultfile << "Inbreeding coefficients" << endl ;
  resultfile << "Population" << '\t' << "Mean";
  if ( printSd )
    resultfile  << '\t' << "Sd" ;
  resultfile  << endl;
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    {
      F_mean[q] /= ( ( double ) samplesize );
      F_sd[q] = sqrt ( F_sd[q] / ( ( double ) samplesize ) - F_mean[q] * F_mean[q] );
      resultfile << q + 1 << '\t' << F_mean[q];
      if ( printSd )
        resultfile << '\t' << F_sd[q] ;
      resultfile << endl;
    }

  resultfile << endl ;
  resultfile << "Allele frequencies" << '\t' << "Mean";
  if ( printSd )
    resultfile << '\t' << "Sd";
  resultfile << endl;
  for ( int l = 0; l < genodata->get_nbpop (); l++ )
    {
      resultfile << "In population " << l + 1 << endl;
      for ( int j = 0; j < genodata->get_nbloci (); j++ )
        {
          resultfile << "At locus " << j + 1 << endl;
          for ( int i = 1; i < genodata->get_nballel (); i++ )
            {
              p_mean[l][j][i] /= ( ( double ) samplesize );
              p_sd[l][j][i] = sqrt ( p_sd[l][j][i] / ( ( double ) samplesize ) - p_mean[l][j][i] * p_mean[l][j][i] );
              resultfile << i << '\t' << p_mean[l][j][i];
              if ( printSd )
                resultfile << '\t' << p_sd[l][j][i];
              resultfile << endl;
            }
        }
      resultfile << endl ;
    }

  if ( !genodata->get_nodata ())
    {
      int ***cum_post_assign;
      cum_post_assign = new int **[genodata->get_nbpop () ];
      if ( cum_post_assign )
        for ( int q = 0; q < genodata->get_nbpop (); q++ )
          {
            cum_post_assign[q] = new int *[genodata->get_nbpop () ];
            if ( cum_post_assign[q] )
              for ( int l = 0; l < genodata->get_nbpop (); l++ )
                {
                  cum_post_assign[q][l] = new int[genodata->get_nbpop () ];
                  if ( cum_post_assign[q][l] )
                    for ( int k = 0; k < genodata->get_nbpop (); k++ )
                      cum_post_assign[q][l][k] = 0;
                  else
                    exit ( 0 );
                }
            else
              exit ( 0 );
          }
      else
        exit ( 0 );

      resultfile << endl ;

      string indivqfilename = fileprefix + ".indivq";
      ofstream indivqfile;
      if (printindivq)
        {
          indivqfile.open(indivqfilename.data());
          if (!indivqfile)
            {
              cerr << "Can't open output file" << endl;
              exit(0);
            }
        }
      for ( int h = 0; h < genodata->get_nbindiv (); h++ )
        {
          if (printindivq)
            {
              indivqfile << h + 1 << '\t' << "Ind#" << h << '\t' << "(0)";
              indivqfile << '\t' << genodata->get_Sh ( h ) + 1 << '\t' << ":";
            }
          for ( int q = 0; q < genodata->get_nbpop (); q++ )
            {
              if (printindivq)
                {
                  indivqfile << '\t';
                  indivqfile.setf(ios::fixed );
                  indivqfile.precision(3);
                  indivqfile << ( ( double ) post_assign[h][q][q] ) / ( ( double ) samplesize );
                }
              cum_post_assign[genodata->get_Sh ( h ) ][q][q] += post_assign[h][q][q];
            }
          for ( int q = 0; q < genodata->get_nbpop () - 1; q++ )
            for ( int l = q + 1; l < genodata->get_nbpop (); l++ )
              {
                if (printindivq)
                  {
                    indivqfile << '\t';
                    indivqfile.setf(ios::fixed );
                    indivqfile.precision(3);
                    indivqfile << ( ( double ) post_assign[h][q][l] ) / ( ( double ) samplesize );
                  }
                cum_post_assign[genodata->get_Sh ( h ) ][q][l] += post_assign[h][q][l];
              }
          if (printindivq)
            indivqfile << endl ;
        }
      if (printindivq)
        indivqfile.close();

      string popqfilename = fileprefix + ".popq";
      ofstream popqfile;
      if (printpopq)
        {
          popqfile.open(popqfilename.data());
          if (!popqfile)
            {
              cerr << "Can't open output file" << endl;
              exit(0);
            }
          for ( int pop = 0; pop < genodata->get_nbpop (); pop++ )
            {
              popqfile << pop + 1 << ":";
              for ( int q = 0; q < genodata->get_nbpop (); q++ )
                {
                  popqfile << '\t';
                  popqfile.setf(ios::fixed );
                  popqfile.precision(3);
                  popqfile << ( ( double ) cum_post_assign[pop][q][q] ) / ( ( double ) samplesize * genodata->get_Nq ( pop ) );
                }
              for ( int q = 0; q < genodata->get_nbpop () - 1; q++ )
                for ( int l = q + 1; l < genodata->get_nbpop (); l++ )
                  {
                    popqfile << '\t';
                    popqfile.setf(ios::fixed );
                    popqfile.precision(3);
                    popqfile << ( ( double ) cum_post_assign[pop][q][l] ) / ( ( double ) samplesize * genodata->get_Nq ( pop ) );
                  }
              popqfile << '\t' << genodata->get_Nq ( pop ) << endl;

              for (int q=0;q<genodata->get_nbpop();q++)
                delete [] cum_post_assign[pop][q];
              delete [] cum_post_assign[pop];
            }
          delete [] cum_post_assign;
          popqfile.close();
        }
    }
}

inline void MCMC::init_stat ()
{
  m_mean = new double *[genodata->get_nbpop () ];
  if ( m_mean )
    for ( int q = 0; q < genodata->get_nbpop (); q++ )
      {
        m_mean[q] = new double[genodata->get_nbpop () ];
        if ( m_mean[q] )
          for ( int l = 0; l < genodata->get_nbpop (); l++ )
            m_mean[q][l] = 0;
        else
          exit ( 0 );
      }
  else
    exit ( 0 );

  m_sd = new double *[genodata->get_nbpop () ];
  if ( m_sd )
    for ( int q = 0; q < genodata->get_nbpop (); q++ )
      {
        m_sd[q] = new double[genodata->get_nbpop () ];
        if ( m_mean[q] )
          for ( int l = 0; l < genodata->get_nbpop (); l++ )
            m_sd[q][l] = 0;
        else
          exit ( 0 );
      }
  else
    exit ( 0 );

  post_assign = new int **[genodata->get_nbindiv () ];
  if ( post_assign )
    for ( int h = 0; h < genodata->get_nbindiv (); h++ )
      {
        post_assign[h] = new int *[genodata->get_nbpop () ];
        if ( post_assign[h] )
          for ( int q = 0; q < genodata->get_nbpop (); q++ )
            {
              post_assign[h][q] = new int[genodata->get_nbpop () ];
              if ( post_assign[h][q] )
                for ( int g = 0; g < genodata->get_nbpop (); g++ )
                  post_assign[h][q][g] = 0;
              else
                exit ( 0 );
            }
        else
          exit ( 0 );
      }
  else
    exit ( 0 );

  p_mean = new double **[genodata->get_nbpop () ];
  if ( p_mean )
    for ( int q = 0; q < genodata->get_nbpop (); q++ )
      {
        p_mean[q] = new double *[genodata->get_nbloci () ];
        if ( p_mean[q] )
          for ( int j = 0; j < genodata->get_nbloci (); j++ )
            {
              p_mean[q][j] = new double[genodata->get_nballel () ];
              if ( !p_mean[q][j] )
                exit ( 0 );
              else
                for ( int i = 0; i < genodata->get_nballel (); i++ )
                  p_mean[q][j][i] = 0;
            }
        else
          exit ( 0 );
      }
  else
    exit ( 0 );

  p_sd = new double **[genodata->get_nbpop () ];
  if ( p_sd )
    for ( int q = 0; q < genodata->get_nbpop (); q++ )
      {
        p_sd[q] = new double *[genodata->get_nbloci () ];
        if ( p_sd[q] )
          for ( int j = 0; j < genodata->get_nbloci (); j++ )
            {
              p_sd[q][j] = new double[genodata->get_nballel () ];
              if ( !p_sd[q][j] )
                exit ( 0 );
              else
                for ( int i = 0; i < genodata->get_nballel (); i++ )
                  p_sd[q][j][i] = 0;
            }
        else
          exit ( 0 );
      }
  else
    exit ( 0 );

  F_mean = new double[genodata->get_nbpop () ];
  if ( F_mean )
    for ( int q = 0; q < genodata->get_nbpop (); q++ )
      F_mean[q] = 0;
  else
    exit ( 0 );

  F_sd = new double[genodata->get_nbpop () ];
  if ( F_sd )
    for ( int q = 0; q < genodata->get_nbpop (); q++ )
      F_sd[q] = 0;
  else
    exit ( 0 );

}

inline void MCMC::init_alpha ()
{
  /* Regression parameters */
  alpha = new double[factors->get_size()];
  if ( !alpha )
    exit ( 0 );
  for ( int r = 0; r < factors->get_size(); r++ )
//     alpha[r] = generator.rnorm ( 0,1 );
    alpha[r]=gsl_ran_gaussian (rng,1);

  mu = new double *[genodata->get_nbpop () ];
  if ( !mu )
    exit ( 0 );
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    {
      mu[q] = new double[genodata->get_nbpop () ];
      if ( !mu[q] )
        exit ( 0 );
      for ( int l = 0; l < genodata->get_nbpop (); l++ )
        mu[q][l] = 0;
    }

  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    for ( int l = 0; l < genodata->get_nbpop (); l++ )
      if ( q != l )
        for ( int r = 0; r < factors->get_size(); r++ )
          mu[q][l] += alpha[r] * factors->get_G(r,q,l);
}

double MCMC::update_alpha ()
{
  /* For migration rates */
  double moy = 0;
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    for ( int l = 0; l < genodata->get_nbpop (); l++ )
      if ( l != q )
        moy += log ( psi[q][l] );

  double v2 = 1 / ( genodata->get_nbpop () * ( genodata->get_nbpop () - 1 ) / s2 + 1 / s2_alpha );
  moy *= v2 / s2;
//   alpha[0] = generator.rnorm ( moy, v2 );
  alpha[0] = moy + sqrt(v2)*gsl_ran_gaussian(rng,1);

  double loglikmod = -0.5 * log ( 2 * M_PI * v2 ) - ( alpha[0] - moy ) * ( alpha[0] - moy ) / ( 2 * v2 );

  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    for ( int l = 0; l < genodata->get_nbpop (); l++ )
      if ( l != q )
        mu[q][l] = alpha[0] * factors->get_G(0,q,l);

  for ( int r = 1; r < factors->get_size(); r++ )
    alpha[r] = 0;

  return loglikmod;
}

inline void MCMC::init_s2 ()
{
//   s2 = 1 / generator.rstgam ( atau + genodata->get_nbpop () * ( genodata->get_nbpop () - 1 ) / 2, 1/btau );
  s2= 1 / gsl_ran_gamma(rng,atau + genodata->get_nbpop () * ( genodata->get_nbpop () - 1 ) / 2,1/btau );
}

inline void MCMC::update_s2 ()
{
  double sum_ql = 0;
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    for ( int l = 0; l < genodata->get_nbpop (); l++ )
      if ( q != l )
        sum_ql += ( log ( psi[q][l] ) - mu[q][l] ) * ( log ( psi[q][l] ) - mu[q][l] );
//   s2 = 1 / generator.rstgam ( atau + genodata->get_nbpop () * ( genodata->get_nbpop () - 1 ) / 2.0, 1.0 / ( btau + sum_ql / 2 ) );
  s2 = 1 / gsl_ran_gamma(rng,atau + genodata->get_nbpop () * ( genodata->get_nbpop () - 1 ) / 2.0, 1.0 / ( btau + sum_ql / 2 ) );
}

inline void MCMC::init_psi ()
{

  psi_acc = 0;

  psi = new double *[genodata->get_nbpop () ];
  if ( psi )
    for ( int q = 0; q < genodata->get_nbpop (); q++ )
      {
        psi[q] = new double[genodata->get_nbpop () ];
        if ( !psi[q] )
          exit ( 0 );
        psi[q][q] = 0;
        for ( int l = 0; l < genodata->get_nbpop (); l++ )
          if ( l != q )
            {
//               psi[q][l] = exp ( generator.rnorm ( mu[q][l], s2 ) );
              psi[q][l] = exp ( mu[q][l] + sqrt(s2)*gsl_ran_gaussian(rng,1)  );
              psi[q][q] += psi[q][l];
            }
      }
  else
    exit ( 0 );
}

inline void MCMC::update_psi ()
{
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    for ( int l = 0; l < genodata->get_nbpop (); l++ )
      if ( l != q )
        {
          /* Save current value */
          double old_psi_ql = psi[q][l];

          /* Propose a new value */
          psi[q][l] = exp ( log ( psi[q][l] ) + sqrt(s2_psi)*gsl_ran_gaussian(rng,1) );

          /* Acceptance ratio */
          // Proposal
          double A = log ( psi[q][l] ) - log ( old_psi_ql );

          // f(m|psi)
          A += gsl_sf_lngamma ( psi[q][q] + psi[q][l] - old_psi_ql ) + gsl_sf_lngamma ( old_psi_ql );
          A -= ( gsl_sf_lngamma ( psi[q][q] ) + gsl_sf_lngamma ( psi[q][l] ) );
          A += ( psi[q][l] - old_psi_ql ) * ( log ( m[q][l] ) - log ( 1 - m[q][q] ) );

          // f(psi|alpha,s2,G)
          A += log ( old_psi_ql ) - log ( psi[q][l] );
          A -= ( log ( psi[q][l] ) - log ( old_psi_ql ) ) * ( log ( psi[q][l] ) + log ( old_psi_ql ) - 2 * mu[q][l] ) / ( 2 * s2 );

          /* Accept or reject proposed value */
          if ( log ( gsl_ran_flat(rng,0,1) ) < A )
            {
              psi_acc++;
              psi[q][q] += psi[q][l] - old_psi_ql;
            }
          else
            {
              psi[q][l] = old_psi_ql;
            }
        }
}

inline void MCMC::init_m ()
{

  m_acc = 0;
  nm_acc = 0;

  m = new double *[genodata->get_nbpop () ];
  if ( m )
    for ( int q = 0; q < genodata->get_nbpop (); q++ )
      {
        m[q] = new double[genodata->get_nbpop () ];
        if ( !m[q] )
          exit ( 0 );
      }
  else
    exit ( 0 );

  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    {
      m[q][q] = gsl_ran_flat(rng,2.0 / 3.0, 1 );
      int i = 0;
      for ( int l = 0; l < genodata->get_nbpop (); l++ )
        if ( l != q )
          {
            m[q][l] = ( 1 - m[q][q] ) / ( genodata->get_nbpop () - 1 );
            i++;
          }
    }
}

inline void MCMC::update_non_migrant_proportions ()
{
  double *old_m_q = new double[genodata->get_nbpop () ];
  if ( !old_m_q )
    exit ( 0 );

  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    {
      /* Save current values */
      for ( int l = 0; l < genodata->get_nbpop (); l++ )
        old_m_q[l] = m[q][l];

      double upper = min ( m[q][q] + e_nm, 1.0 );
      double lower = max ( m[q][q] - e_nm, 0.0 );

      /* Proposed value */
      m[q][q] = gsl_ran_flat(rng,lower,upper);
      for ( int l = 0; l < genodata->get_nbpop (); l++ )
        if ( l != q )
          m[q][l] *= ( 1 - m[q][q] ) / ( 1 - old_m_q[q] );

      /* New bounds for proposal density */
      double new_upper = min ( m[q][q] + e_nm, 1.0 );
      double new_lower = max ( m[q][q] - e_nm, 0.0 );

      // Proposal ratio
      double A = ( log ( upper - lower ) - log ( new_upper - new_lower ) );

      // f(m|psi)
      A += ( genodata->get_nbpop () - 2 ) * ( log ( 1 - old_m_q[q] ) - log ( 1 - m[q][q] ) );

      // Pr(M,t|m)
      double old_pop_lglik = pop_lglik[q];
      pop_loglik ( q );
      A += pop_lglik[q] - old_pop_lglik;

      if ( log ( gsl_ran_flat(rng,0,1) ) < A )
        {
          nm_acc++;
          lglikMt += pop_lglik[q] - old_pop_lglik;
          if ( lglikMt > 0 )
            {
              EXCEPTION_INFOS();
              exit ( 0 );
            }
        }
      else
        {
          pop_lglik[q] = old_pop_lglik;
          for ( int l = 0; l < genodata->get_nbpop (); l++ )
            m[q][l] = old_m_q[l];
        }
    }
  delete [] old_m_q;
}

inline void MCMC::update_migration_rate_matrix()
{
  for ( int q=0;q<genodata->get_nbpop();q++ )
    {
      /* Choose two migrant sources */
      int l1 = gsl_rng_uniform_int(rng, genodata->get_nbpop () );
      int l2;
      do
        l2 = gsl_rng_uniform_int(rng, genodata->get_nbpop () );
      while ( l1 == l2 );

      /* Save old values */
      double old_m_ql1 = m[q][l1];
      double old_m_ql2 = m[q][l2];

      /* Bounds of proposal density */
      double upper = min ( m[q][l1] + e_m, m[q][l1] + m[q][l2] );
      double lower = max ( m[q][l1] - e_m, 0.0 );

      /* Proposed values */
      m[q][l1] = gsl_ran_flat(rng,lower,upper);
      m[q][l2] = old_m_ql1 + old_m_ql2 - m[q][l1];

      /* New bounds for proposal density */
      double new_upper = min ( m[q][l1] + e_m, m[q][l1] + m[q][l2] );
      double new_lower = max ( m[q][l1] - e_m, 0.0 );

      /* Calculate the log of acceptance value */
      double A;

      /* 1) Ratio of proposal densities */
      A = ( log ( upper - lower ) - log ( new_upper - new_lower ) );

      /* 2) Ratio of M,t given m prior density */
      double lglikMtratio = 0;
      if ( !genodata->get_nodata () )
        {
          int sum1 = ancescount[q][l1][l1];
          int sum2 = ancescount[q][l2][l2];
          for ( int g = 0; g < genodata->get_nbpop (); g++ )
            {
              sum1 += ancescount[q][l1][g];
              sum2 += ancescount[q][l2][g];
            }

          lglikMtratio += sum1 * ( log ( m[q][l1] ) - log ( old_m_ql1 ) );
          lglikMtratio += sum2 * ( log ( m[q][l2] ) - log ( old_m_ql2 ) );

          A += lglikMtratio;
        }

      /* 3) Prior */
      A += ( Psi - 1 ) * ( ( log ( m[q][l1] ) - log ( old_m_ql1 ) ) + ( log ( m[q][l2] ) - log ( old_m_ql2 ) ) );

      if ( log ( gsl_ran_flat(rng,0,1) ) < A )
        {
          m_acc++;

          if ( !genodata->get_nodata () )
            {
              pop_lglik[q] += lglikMtratio;
              lglikMt += lglikMtratio;
            }
        }
      else
        {
          m[q][l1] = old_m_ql1;
          m[q][l2] = old_m_ql2;
        }
    }
}

inline void MCMC::update_m_2pop()
{
  for ( int q=0;q<genodata->get_nbpop();q++ )
    {
      /* The other population */
      int l = ( q+1 ) % 2;

      /* Save current value */
      double old_m_qq = m[q][q];

      /* Bounds of proposal density */
      double upper = min ( m[q][q] + e_m, 1.0 );
      double lower = max ( m[q][q] - e_m, 0.0 );

      /* Prospose new values */
      m[q][q] = gsl_ran_flat(rng,lower,upper);

      /* New bounds for proposal density */
      double new_upper = min ( m[q][q] + e_m, 1.0 );
      double new_lower = max ( m[q][q] - e_m, 0.0 );

      /* Calculate the log of acceptance value */
      double A;

      /* 1) Ratio of proposal densities */
      A = ( log ( upper - lower ) - log ( new_upper - new_lower ) );

      /* 2) Ratio of M,t given m prior density */
      double lglikMtratio = 0;
      if ( !genodata->get_nodata () )
        {
          int sum1 = 2*ancescount[q][q][q] + ancescount[q][q][l];
          int sum2 = 2*ancescount[q][l][l] + ancescount[q][q][l];

          lglikMtratio += sum1 * ( log ( m[q][q] ) - log ( old_m_qq ) );
          lglikMtratio += sum2 * ( log ( 1-m[q][q] ) - log ( 1-old_m_qq ) );

          A += lglikMtratio;
        }

      /* 3) Prior */
      A += ( Psi-1 ) * ( log ( m[q][q] ) - log ( old_m_qq ) + log ( 1-m[q][q] ) - log ( 1-old_m_qq ) );

      /* Accept or reject proposed values */
      if ( log ( gsl_ran_flat(rng,0,1) ) < A )
        {
          m_acc++;

          if ( !genodata->get_nodata () )
            {
              pop_lglik[q] += lglikMtratio;
              lglikMt += lglikMtratio;
            }
        }
      else
        m[q][q] = old_m_qq;

      m[q][l] =  1-m[q][q];
    }
}

inline void MCMC::update_m ()
{
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    {
      /* Choose two migrant sources */
      int l1;
      do
        l1 = gsl_rng_uniform_int(rng, genodata->get_nbpop () );
      while ( l1 == q );
      int l2;
      do
        l2 = gsl_rng_uniform_int(rng, genodata->get_nbpop () );
      while ( ( l1 == l2 ) || ( l2 == q ) );

      /* Save old values */
      double old_m_ql1 = m[q][l1];
      double old_m_ql2 = m[q][l2];

      /* Bounds of proposal density */
      double upper = min ( m[q][l1] + e_m, m[q][l1] + m[q][l2] );
      double lower = max ( m[q][l1] - e_m, 0.0 );

      /* Proposed values */
      m[q][l1] = gsl_ran_flat(rng,lower,upper);
      m[q][l2] = old_m_ql1 + old_m_ql2 - m[q][l1];

      /* New bounds for proposal density */
      double new_upper = min ( m[q][l1] + e_m, m[q][l1] + m[q][l2] );
      double new_lower = max ( m[q][l1] - e_m, 0.0 );

      /* Calculate the log of acceptance value */
      double A;

      /* 1) Ratio of proposal densities */
      A = ( log ( upper - lower ) - log ( new_upper - new_lower ) );

      /* 2) Ratio of M,t given m prior density */
      double lglikMtratio = 0;
      if ( !genodata->get_nodata () )
        {
          int sum1 = ancescount[q][l1][l1];
          int sum2 = ancescount[q][l2][l2];
          for ( int g = 0; g < genodata->get_nbpop (); g++ )
            {
              sum1 += ancescount[q][l1][g];
              sum2 += ancescount[q][l2][g];
            }

          lglikMtratio += sum1 * ( log ( m[q][l1] ) - log ( old_m_ql1 ) );
          lglikMtratio += sum2 * ( log ( m[q][l2] ) - log ( old_m_ql2 ) );

          A += lglikMtratio;
        }

      /* 3) Prior */
      if ( useReg )
        A += ( psi[q][l1] - 1 ) * ( log ( m[q][l1] ) - log ( old_m_ql1 ) ) + ( psi[q][l2] - 1 ) * ( log ( m[q][l2] ) - log ( old_m_ql2 ) );
      else
        A+= ( Psi-1 ) * ( log ( m[q][l1] ) - log ( old_m_ql1 ) +log ( m[q][l2] ) - log ( old_m_ql2 ) );

      if ( log ( gsl_ran_flat(rng,0,1) ) < A )
        {
          m_acc++;

          if ( !genodata->get_nodata () )
            {
              pop_lglik[q] += lglikMtratio;
              lglikMt += lglikMtratio;
            }
        }
      else
        {
          m[q][l1] = old_m_ql1;
          m[q][l2] = old_m_ql2;
        }
    }
}

inline void MCMC::init_Mt ()
{

  Mt_acc = 0;

  Mt = new int *[genodata->get_nbindiv () ];
  if ( Mt )
    for ( int h = 0; h < genodata->get_nbindiv (); h++ )
      {
        Mt[h] = new int[2];
        if ( !Mt[h] )
          exit ( 0 );
      }
  else
    exit ( 0 );

  ancescount = new int **[genodata->get_nbpop () ];
  if ( ancescount )
    for ( int q = 0; q < genodata->get_nbpop (); q++ )
      {
        ancescount[q] = new int *[genodata->get_nbpop () ];
        if ( ancescount[q] )
          for ( int l = 0; l < genodata->get_nbpop (); l++ )
            {
              ancescount[q][l] = new int[genodata->get_nbpop () ];
              if ( !ancescount[q][l] )
                exit ( 0 );
              else
                for ( int g = 0; g < genodata->get_nbpop (); g++ )
                  ancescount[q][l][g] = 0;
            }
        else
          exit ( 0 );
      }
  else
    exit ( 0 );

  for ( int h = 0; h < genodata->get_nbindiv (); h++ )
    {
      if ( gsl_ran_flat(rng,0,1) <= 0.7 )
        {
          Mt[h][0] = genodata->get_Sh ( h );
          Mt[h][1] = Mt[h][0];
        }
      else
        {
          do
            Mt[h][0] = gsl_rng_uniform_int(rng, genodata->get_nbpop () );
          while ( Mt[h][0] == genodata->get_Sh ( h ) );

          Mt[h][1] = gsl_rng_uniform_int(rng, genodata->get_nbpop () );
        }
      ancescount[genodata->get_Sh ( h ) ][Mt[h][1]][Mt[h][0]] = ++ancescount[genodata->get_Sh ( h ) ][Mt[h][0]][Mt[h][1]];
    }
}

inline void MCMC::update_Mt ()
{
  for ( int h = 0; h < genodata->get_nbindiv (); h++ )
    {
      /* Save current value of individual ancestry */
      int old_M = Mt[h][0];
      int old_g = Mt[h][1];

      /* Proposed ancestry */
      do
        {
          Mt[h][0] = gsl_rng_uniform_int(rng, genodata->get_nbpop () );
          Mt[h][1] = gsl_rng_uniform_int(rng, genodata->get_nbpop () );
        }
      while ( ( Mt[h][0] == old_M ) && ( Mt[h][1] == old_g ) );

      /* Update counter */
      ancescount[genodata->get_Sh ( h ) ][Mt[h][1]][Mt[h][0]] = ++ancescount[genodata->get_Sh ( h ) ][Mt[h][0]][Mt[h][1]];

      /* Calculate the log of acceptance value */
      double A;

      /* 1) Proposal */
      A = log ( ancescount[genodata->get_Sh ( h ) ][Mt[h][0]][Mt[h][1]] ) - log ( ancescount[genodata->get_Sh ( h ) ][old_M][old_g] );

      /* 2) X,S given M,t,F,p */
      double lglikgenoratio = 0;
      if ( !genodata->get_missingdata () )
        {
          tmp_ind_lglik[h] = ind_lglik[h];
          indiv_loglik ( h );
          lglikgenoratio = ( ind_lglik[h] - tmp_ind_lglik[h] );
          A += lglikgenoratio;
        }

      /* 3) M,t given m */
      double lglikMtratio = 0;
      lglikMtratio += log ( ancescount[genodata->get_Sh ( h ) ][old_M][old_g] ) - log ( ancescount[genodata->get_Sh ( h ) ][Mt[h][0]][Mt[h][1]] );
      lglikMtratio += ( Mt[h][0] == Mt[h][1] ? 0 : log ( 2 ) ) + log ( m[genodata->get_Sh ( h ) ][Mt[h][0]] ) + log ( m[genodata->get_Sh ( h ) ][Mt[h][1]] );
      lglikMtratio -= ( old_M == old_g ? 0 : log ( 2 ) ) + log ( m[genodata->get_Sh ( h ) ][old_M] ) + log ( m[genodata->get_Sh ( h ) ][old_g] );
      A += lglikMtratio;

      /* Accept or reject proposed value */
      if ( log ( gsl_ran_flat(rng,0,1) ) < A )
        {
          Mt_acc++;

          ancescount[genodata->get_Sh ( h ) ][old_g][old_M] = --ancescount[genodata->get_Sh ( h ) ][old_M][old_g];

          /* Update log(Pr(M,t|m) */
          pop_lglik[genodata->get_Sh ( h ) ] += lglikMtratio;
          lglikMt += lglikMtratio;

          /* Update Pr(X|S;M,t,F,p) */
          if ( !genodata->get_missingdata () )
            lglikgeno += lglikgenoratio;
        }
      else
        {
          ancescount[genodata->get_Sh ( h ) ][Mt[h][0]][Mt[h][1]] = --ancescount[genodata->get_Sh ( h ) ][Mt[h][1]][Mt[h][0]];
          Mt[h][0] = old_M;
          Mt[h][1] = old_g;
          if ( !genodata->get_missingdata () )
            ind_lglik[h] = tmp_ind_lglik[h];
        }
    }
}

inline void MCMC::init_p ()
{

  p_acc = 0;

  freqmin = new double[genodata->get_nbpop () ];
  if ( !freqmin )
    exit ( 0 );

  p = new double **[genodata->get_nbpop () ];
  if ( p )
    for ( int q = 0; q < genodata->get_nbpop (); q++ )
      {
        p[q] = new double *[genodata->get_nbloci () ];
        if ( p[q] )
          for ( int j = 0; j < genodata->get_nbloci (); j++ )
            {
              p[q][j] = new double[genodata->get_nballel () ];
              if ( !p[q][j] )
                exit ( 0 );
            }
        else
          exit ( 0 );
      }
  else
    exit ( 0 );

  double *vect = new double[genodata->get_nballel () ];
  if ( !vect )
    exit ( 0 );
  vect[0] = 0;
  for ( int l = 0; l < genodata->get_nbpop (); l++ )
    {
      freqmin[l] = 1;
      for ( int j = 0; j < genodata->get_nbloci (); j++ )
        {
          for ( int i = 1; i < genodata->get_nballel (); i++ )
            vect[i] = ( genodata->get_is_at_locus ( j, i ) ? genodata->get_allelcount ( l, j, i ) + 1 : 0 );

//           generator.rdirichlet ( p[l][j], vect, genodata->get_nballel () );
          gsl_ran_dirichlet(rng,genodata->get_nballel (),vect,p[l][j]);

          for ( int i = 1; i < genodata->get_nballel (); i++ )
            if ( genodata->get_is_at_locus ( j, i ) && ( p[l][j][i] < freqmin[l] ) )
              freqmin[l] = p[l][j][i];
        }
    }
  delete [] vect;
}

inline void MCMC::update_p ()
{
  for ( int l = 0; l < genodata->get_nbpop (); l++ )
    for ( int j = 0; j < genodata->get_nbloci (); j++ )
      {
        /* Choose an allele i */
        int i;
        do
          i = gsl_rng_uniform_int(rng, genodata->get_nballel () );
        while ( !genodata->get_is_at_locus ( j, i ) || ( i == 0 ) );

        /* Choose another allele k */
        int k;
        do
          k = gsl_rng_uniform_int(rng, genodata->get_nballel () );
        while ( ( k == i ) || ( !genodata->get_is_at_locus ( j, k ) ) || ( k == 0 ) );

        /* Store current values */
        double old_p_lji = p[l][j][i];
        double old_p_ljk = p[l][j][k];

        double p_min = max ( 0.0, F[l] / ( F[l] - 1 ) );

        /* Boundaries of proposed values */
        double upper = min ( p[l][j][i] + p[l][j][k] - p_min, p[l][j][i] + e_p );
        double lower = max ( p_min, p[l][j][i] - e_p );

        /* Generate proposed values */
        p[l][j][i] = gsl_ran_flat(rng,lower,upper);
        p[l][j][k] = old_p_lji + old_p_ljk - p[l][j][i];

        /* New boundaries */
        double new_upper = min ( p[l][j][i] + p[l][j][k] - p_min, p[l][j][i] + e_p );
        double new_lower = max ( p_min, p[l][j][i] - e_p );

        /* Compute the log of the acceptance probability */
        double A;

        /* 1) Ratio of proposal densities */
        A = log ( upper - lower ) - log ( new_upper - new_lower );

        /* 2) Likelihood ratio Pr(X|S;M,t,F,p) */
        double old_lglikgeno = lglikgeno;
        if ( !genodata->get_nodata () && !genodata->get_missingdata () )
          {
            for ( int h = 0; h < genodata->get_nbindiv (); h++ )
              {
                /* Assignments and genotypes */
                if ( ( genodata->get_Xhji ( h,j,0 ) ==i ) || ( genodata->get_Xhji ( h,j,0 ) ==k ) || ( genodata->get_Xhji ( h,j,1 ) ==i ) || ( genodata->get_Xhji ( h,j,1 ) ==k ) )
                  if ( ( Mt[h][0] == l ) || ( Mt[h][1] == l ) )
                    {
                      tmp_ind_lglik[h] = ind_lglik[h];
                      double test = ind_lglik[h];

                      if ( Mt[h][0] == Mt[h][1] )
                        test += log ( Phi ( h, j ) );

                      /* Genotype homozygote (i,i) or (k,k) */
                      if ( ( ( genodata->get_Xhji ( h,j,0 ) ==i ) && ( genodata->get_Xhji ( h,j,1 ) ==i ) ) || ( ( genodata->get_Xhji ( h,j,0 ) ==k ) && ( genodata->get_Xhji ( h,j,1 ) ==k ) ) )
                        {
                          int all = ( genodata->get_Xhji ( h, j, 0 ) == i ? i : k );
                          double freq = ( genodata->get_Xhji ( h, j, 0 ) == i ? old_p_lji : old_p_ljk );
                          switch ( Mt[h][0] == Mt[h][1] )
                            {
                            case true:
                              test -= log ( ( 1 - F[l] ) * pow ( freq, 2 ) + F[l] * freq );
                              break;
                            case false:
                              test += log ( p[l][j][all] );
                              test -= log ( freq );
                              break;
                            default:
                              EXCEPTION_INFOS();
                              exit ( 0 );
                            }
                        }

                      /* Genotype heterozygote (i,k) */
                      else
                        if ( ( ( genodata->get_Xhji ( h,j,0 ) ==i ) && ( genodata->get_Xhji ( h,j,1 ) ==k ) ) || ( ( genodata->get_Xhji ( h,j,0 ) ==k ) && ( genodata->get_Xhji ( h,j,1 ) ==i ) ) )
                          {
                            int pop;
                            switch ( Mt[h][0] == Mt[h][1] )
                              {
                              case true:
                                test -= log ( 2 * ( 1 - F[l] ) * old_p_lji * old_p_ljk );
                                break;
                              case false:
                                pop = ( Mt[h][0] == l ? Mt[h][1] : Mt[h][0] );
                                test += log ( phi ( h, j ) );
                                test -= log ( old_p_lji * p[pop][j][k] + old_p_ljk * p[pop][j][i] );
                                break;
                              default:
                                EXCEPTION_INFOS();
                                exit ( 0 );
                              }
                          }

                      /* Genotype heterozygote (i,.) */
                        else
                          if ( ( ( genodata->get_Xhji ( h, j, 0 ) == i ) && ( genodata->get_Xhji ( h, j, 1 ) != k ) && ( genodata->get_Xhji ( h, j, 1 ) != i ) ) || ( ( genodata->get_Xhji ( h, j, 0 ) != k ) && ( genodata->get_Xhji ( h, j, 0 ) != i ) && ( genodata->get_Xhji ( h, j, 1 ) == i ) ) )
                            {
                              int pop;
                              /* allele */
                              int all =
                                ( genodata->get_Xhji ( h, j, 0 ) == i ? genodata->get_Xhji ( h, j, 1 ) : genodata->get_Xhji ( h, j, 0 ) );
                              switch ( Mt[h][0] == Mt[h][1] )
                                {
                                case true:
                                  test -= log ( 2 * ( 1 - F[l] ) * old_p_lji * p[l][j][all] );
                                  break;
                                case false:
                                  pop = ( Mt[h][0] == l ? Mt[h][1] : Mt[h][0] );
                                  test += log ( phi ( h, j ) );
                                  test -= log ( old_p_lji * p[pop][j][all] + p[l][j][all] * p[pop][j][i] );
                                  break;
                                default:
                                  EXCEPTION_INFOS();
                                  exit ( 0 );
                                }
                            }

                      /* Genotype heterozygote (k,.) */
                          else
                            if ( ( ( genodata->get_Xhji ( h, j, 0 ) != i ) && ( genodata->get_Xhji ( h, j, 0 ) != k ) && ( genodata->get_Xhji ( h, j, 1 ) == k ) ) || ( ( genodata->get_Xhji ( h, j, 0 ) == k ) && ( genodata->get_Xhji ( h, j, 1 ) != i ) && ( genodata->get_Xhji ( h, j, 1 ) != k ) ) )
                              {
                                int pop;
                                int all =
                                  ( genodata->get_Xhji ( h, j, 0 ) == k ? genodata->get_Xhji ( h, j, 1 ) : genodata->get_Xhji ( h, j, 0 ) );
                                switch ( Mt[h][0] == Mt[h][1] )
                                  {
                                  case true:
                                    test -= log ( 2 * ( 1 - F[l] ) * old_p_ljk * p[l][j][all] );
                                    break;
                                  case false:
                                    pop = ( Mt[h][0] == l ? Mt[h][1] : Mt[h][0] );
                                    test += log ( phi ( h, j ) );
                                    test -= log ( old_p_ljk * p[pop][j][all] + p[l][j][all] * p[pop][j][k] );
                                    break;
                                  default:
                                    EXCEPTION_INFOS();
                                    exit ( 0 );
                                  }
                              }
                            else
                              {
                                EXCEPTION_INFOS();
                                exit ( 0 );
                              }

                      ind_lglik[h] = test;
                      lglikgeno += ind_lglik[h] - tmp_ind_lglik[h];
                      A += ( ind_lglik[h] - tmp_ind_lglik[h] );
                    }
              }
          }

        /* Prior */
        if ( useFmodel )
          {
            A += ( theta[l] * glob_all_frq[j][i]-1 ) * ( log ( p[l][j][i] ) - log ( old_p_lji ) );
            A += ( theta[l] * glob_all_frq[j][k]-1 ) * ( log ( p[l][j][k] ) - log ( old_p_ljk ) );
          }
        else
          A += ( lambda-1 ) * ( ( log ( p[l][j][i] ) - log ( old_p_lji ) ) + ( log ( p[l][j][k] ) - log ( old_p_ljk ) ) );

        /* Accept or reject proposed values */
        if ( log ( gsl_ran_flat(rng,0,1) ) < A )
          {
            p_acc++;
            freqmin[l] = min ( freqmin[l], min ( p[l][j][i], p[l][j][k] ) );
          }
        else
          {
            /* Restaure old values */
            p[l][j][i] = old_p_lji;
            p[l][j][k] = old_p_ljk;
            if ( !genodata->get_nodata () && !genodata->get_missingdata () )
              {
                lglikgeno = old_lglikgeno;
                for ( int h = 0; h < genodata->get_nbindiv (); h++ )
                  if ( ( genodata->get_Xhji ( h,j,0 ) ==i ) || ( genodata->get_Xhji ( h,j,0 ) ==k ) || ( genodata->get_Xhji ( h,j,1 ) ==i ) || ( genodata->get_Xhji ( h,j,1 ) ==k ) )
                    if ( ( Mt[h][0] == l ) || ( Mt[h][1] == l ) )
                      ind_lglik[h] = tmp_ind_lglik[h];
              }
          }
      }
}

inline void MCMC::init_F ()
{
  F_acc = 0;

  F = new double[genodata->get_nbpop () ];
  if ( !F )
    exit ( 0 );

  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    {
      double F_min = max ( -1.0, - ( freqmin[q] ) / ( 1 - freqmin[q] ) );
//       F[q] = generator.runiform ( F_min, 1 );
      F[q]=gsl_ran_flat(rng,F_min,1);
    }
}

inline void MCMC::update_F ()
{
  for ( int l=0;l<genodata->get_nbpop ();l++ )
    {
      /* Save current value */
      double old_F = F[l];

      double F_min = max ( -1.0, - ( freqmin[l] ) / ( 1 - freqmin[l] ) );

      /* Boundaries of proposal interval */
      double upper = min ( F[l] + e_F, 1.0 );
      double lower = max ( F[l] - e_F, F_min );

      /* Choose new F[l] uniformly in (lower,upper) */
      F[l] = gsl_ran_flat(rng,lower,upper);

      double new_upper = min ( F[l] + e_F, 1.0 );
      double new_lower = max ( F[l] - e_F, F_min );

      /* Calculate log of acceptance value */
      double A;

      /* 1) Ratio of proposal densities */
      A = log ( upper - lower ) - log ( new_upper - new_lower );

      /* 2) Ratio of individuals'likelihoods */
      double old_lglikgeno = lglikgeno;
      if ( !genodata->get_nodata () && !genodata->get_missingdata () )
        for ( int h = 0; h < genodata->get_nbindiv (); h++ )
          if ( ( Mt[h][0] == l ) && ( Mt[h][1] == l ) )
            {
              tmp_ind_lglik[h] = ind_lglik[h];
              indiv_loglik ( h );
              A += ( ind_lglik[h] - tmp_ind_lglik[h] );

              /* New Pr(X|S;M,t,F,p) */
              lglikgeno += ( ind_lglik[h] - tmp_ind_lglik[h] );
            }

      /* 3) Prior */
      A += 0;

      /* Accept or reject proposed value */
      if ( log ( gsl_ran_flat(rng,0,1) ) < A )
        F_acc++;
      else
        {
          /* Restaure old values */
          F[l] = old_F;

          if ( !genodata->get_nodata () && !genodata->get_missingdata () )
            {
              lglikgeno = old_lglikgeno;
              for ( int h = 0; h < genodata->get_nbindiv (); h++ )
                if ( ( Mt[h][0] == l ) && ( Mt[h][1] == l ) )
                  ind_lglik[h] = tmp_ind_lglik[h];
            }
        }
    }
}

inline void MCMC::init_lglik ()
{

  ind_lglik = new double[genodata->get_nbindiv () ];
  if ( !ind_lglik )
    exit ( 0 );

  tmp_ind_lglik = new double[genodata->get_nbindiv () ];
  if ( !tmp_ind_lglik )
    exit ( 0 );

  /* Assignments likelihood */
  pop_lglik = new double[genodata->get_nbpop () ];
  if ( !pop_lglik )
    exit ( 0 );

  if ( genodata->get_nodata () || genodata->get_missingdata () )
    {
      for ( int h = 0; h < genodata->get_nbindiv (); h++ )
        ind_lglik[h] = 0;
      lglikgeno = 0;
    }
  else
    {
      for ( int h = 0; h < genodata->get_nbindiv (); h++ )
        indiv_loglik ( h );
      loglikgeno ();
    }

  if ( genodata->get_nodata () )
    {
      for ( int q = 0; q < genodata->get_nbpop (); q++ )
        pop_lglik[q] = 0;
      lglikMt = 0;
    }
  else
    {
      for ( int q = 0; q < genodata->get_nbpop (); q++ )
        pop_loglik ( q );
      loglikMt ();
    }
}

inline void MCMC::loglikgeno ()
{
  lglikgeno = 0;
  for ( int h = 0; h < genodata->get_nbindiv (); h++ )
    lglikgeno += ind_lglik[h];
}

void MCMC::indiv_loglik ( int h )
{
  double loc_j = 0;

  /* Individual without admixture */
  if ( Mt[h][0] == Mt[h][1] )
    for ( int j = 0; j < genodata->get_nbloci (); j++ )
      {
        if ( ( genodata->get_Xhji ( h,j,0 ) != 0 ) && ( genodata->get_Xhji ( h,j,1 ) != 0 ) )
          loc_j += log ( Phi ( h, j ) );
      }

  /* Hybrids */
  else if ( Mt[h][0] != Mt[h][1] )
    for ( int j = 0; j < genodata->get_nbloci (); j++ )
      if ( ( genodata->get_Xhji ( h,j,0 ) != 0 ) && ( genodata->get_Xhji ( h,j,1 ) != 0 ) )
        loc_j += log ( phi ( h, j ) );

  if ( loc_j > 0 )
    {
      EXCEPTION_INFOS();
      exit ( 0 );
    }

  ind_lglik[h] = loc_j;
}

double MCMC::Phi ( int h, int j )
{
  if ( genodata->get_Xhji ( h,j,0 ) == genodata->get_Xhji ( h,j,1 ) )
    return ( ( 1-F[Mt[h][0]] ) *pow ( p[Mt[h][0]][j][genodata->get_Xhji ( h,j,0 ) ],2 ) +F[Mt[h][0]]*p[Mt[h][0]][j][genodata->get_Xhji ( h,j,0 ) ] );
  else
    return ( 2* ( 1-F[Mt[h][0]] ) *p[Mt[h][0]][j][genodata->get_Xhji ( h,j,0 ) ]*p[Mt[h][0]][j][genodata->get_Xhji ( h,j,1 ) ] );
}

double MCMC::phi ( int h, int j )
{
  if ( genodata->get_Xhji ( h,j,0 ) == genodata->get_Xhji ( h,j,1 ) )
    return ( p[Mt[h][0]][j][genodata->get_Xhji ( h,j,0 ) ]*p[Mt[h][1]][j][genodata->get_Xhji ( h,j,0 ) ] );
  else
    return ( p[Mt[h][0]][j][genodata->get_Xhji ( h,j,0 ) ]*p[Mt[h][1]][j][genodata->get_Xhji ( h,j,1 ) ]+p[Mt[h][0]][j][genodata->get_Xhji ( h,j,1 ) ]*p[Mt[h][1]][j][genodata->get_Xhji ( h,j,0 ) ] );
}

inline void MCMC::loglikMt ()
{
  lglikMt = 0;
  /* Cycle over populations */
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    lglikMt += pop_lglik[q];
}

void MCMC::pop_loglik ( int q )
{
  pop_lglik[q] = gsl_sf_lnfact( genodata->get_Nq ( q ) );

  /* Cycle over ancestry states */
  for ( int l = 0; l < genodata->get_nbpop (); l++ )
    for ( int g = 0; g <= l; g++ )
      if ( m[q][l] * m[q][g] != 0 )
        pop_lglik[q] += ancescount[q][l][g]* ( ( l == g ? 0 : log ( 2 ) ) +log ( m[q][l] ) +log ( m[q][g] ) )- gsl_sf_lnfact( ancescount[q][l][g] );

  if ( ( pop_lglik[q] > 0 ) || ( isnan ( pop_lglik[q] ) != 0 ) )
    {
      EXCEPTION_INFOS();
      for ( int l = 0; l < genodata->get_nbpop (); l++ )
        cerr << "m" << q + 1 << l + 1 << "=" << m[q][l] << endl;

      for ( int l = 0; l < genodata->get_nbpop (); l++ )
        for ( int g = 0; g <= l; g++ )
          cerr << "n_" << q + 1 << ',' << l + 1 << ',' << g + 1 << '=' << ancescount[q][l][g] << endl;

      exit ( 0 );
    }
}

void MCMC::RJupdate ()
{
  // Save the number of regressors in the current model
  int old_nb_reg = nb_reg_param;

  // Save current model
  int old_model = model;

  // Save current regression parameters
  double *old_alpha = new double[factors->get_size()];
  if ( !old_alpha )
    exit ( 0 );
  for ( int r = 0; r < factors->get_size(); r++ )
    old_alpha[r] = alpha[r];

  // Propose a new model
  do
    model = gsl_rng_uniform_int(rng, nb_model );
  while ( ( model == old_model ) || !factors->is_a_valid_model (model) );

  // Save current lognormal mean
  double **old_mu = new double *[genodata->get_nbpop () ];
  if ( !old_mu )
    exit ( 0 );
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    {
      old_mu[q] = new double[genodata->get_nbpop () ];
      if ( !old_mu[q] )
        exit ( 0 );
      else
        for ( int l = 0; l < genodata->get_nbpop (); l++ )
          old_mu[q][l] = mu[q][l];
    }

  //Generate new parameters
  double newlglikmod = MVregupdate ();

  double A = lgMV - newlglikmod;

  // f(psi|alpha,s2,G)
  double sum = 0;
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    for ( int l = 0; l < genodata->get_nbpop (); l++ )
      if ( l != q )
        sum += ( mu[q][l] - old_mu[q][l] ) * ( 2 * log ( psi[q][l] ) - mu[q][l] - old_mu[q][l] );
  A += sum / ( 2 * s2 );

  // f(alpha)
  sum = 0;
  for ( int r = 0; r < factors->get_size(); r++ )
    sum += old_alpha[r] * old_alpha[r] - alpha[r] * alpha[r];
  A += sum / ( 2 * s2_alpha );
  A -= 0.5* ( nb_reg_param-old_nb_reg ) *log ( 2*M_PI*s2_alpha );

  /* Accept or reject new model */
  if ( log ( gsl_ran_flat(rng,0,1) ) < A )
    {
      RJ_acc++;
      lgMV = newlglikmod;
    }
  else
    {
      nb_reg_param = old_nb_reg;
      model = old_model;
      for ( int q = 0; q < genodata->get_nbpop (); q++ )
        for ( int l = 0; l < genodata->get_nbpop (); l++ )
          if ( l != q )
            mu[q][l] = old_mu[q][l];

      for ( int r = 0; r < factors->get_size(); r++ )
        alpha[r] = old_alpha[r];
    }

  for (int q=0;q<genodata->get_nbpop();q++)
    delete [] old_mu[q];
  delete [] old_mu;
  delete [] old_alpha;
}

inline void MCMC::init_glob_all_frq ()
{
  pglob_acc = 0;

  glob_all_frq = new double *[genodata->get_nbloci () ];
  if ( !glob_all_frq )
    exit ( 0 );
  for ( int j = 0; j < genodata->get_nbloci (); j++ )
    {
      glob_all_frq[j] = new double[genodata->get_nballel () ];
      if ( !glob_all_frq[j] )
        exit ( 0 );
    }

  double *vect = new double[genodata->get_nballel () ];
  if ( !vect )
    exit ( 0 );
  vect[0] = 0;
  for ( int j = 0; j < genodata->get_nbloci (); j++ )
    {
      for ( int k = 1; k < genodata->get_nballel (); k++ )
        {
          vect[k] = ( genodata->get_is_at_locus ( j, k ) ? 1 : 0 );
          for ( int q = 0; q < genodata->get_nbpop (); q++ )
            vect[k] += ( genodata->get_is_at_locus ( j, k ) ? genodata->get_allelcount ( q, j, k ) : 0 );
        }
//       generator.rdirichlet ( glob_all_frq[j], vect, genodata->get_nballel () );
      gsl_ran_dirichlet(rng,genodata->get_nballel (),vect,glob_all_frq[j]);
    }
    delete [] vect;
}

inline void MCMC::update_glob_all_frq ()
{
  for ( int j = 0; j < genodata->get_nbloci (); j++ )
    {
      /* Choose an allele i */
      int i;
      do
        i = gsl_rng_uniform_int(rng, genodata->get_nballel () );
      while ( !genodata->get_is_at_locus ( j, i ) || ( i == 0 ) );

      /* Choose another allele k */
      int k;
      do
        k = gsl_rng_uniform_int(rng, genodata->get_nballel () );
      while ( ( k == i ) || ( !genodata->get_is_at_locus ( j, k ) ) || ( k == 0 ) );

      /* Store current values */
      double old_p_ji = glob_all_frq[j][i];
      double old_p_jk = glob_all_frq[j][k];

      /* Boundaries of proposed values */
      double upper = min ( glob_all_frq[j][i] + glob_all_frq[j][k], glob_all_frq[j][i] + e_p );
      double lower = max ( 0.0, glob_all_frq[j][i] - e_p );

      /* Generate proposed values */
      glob_all_frq[j][i] = gsl_ran_flat(rng,lower,upper);
      glob_all_frq[j][k] = old_p_ji + old_p_jk - glob_all_frq[j][i];

      /* New boundaries */
      double new_upper = min ( glob_all_frq[j][i] + glob_all_frq[j][k],glob_all_frq[j][i] + e_p );
      double new_lower = max ( 0.0, glob_all_frq[j][i] - e_p );

      /* Compute the log of the acceptance probability */
      double A;

      /* 1) Ratio of proposal densities */
      A = log ( upper - lower ) - log ( new_upper - new_lower );

      /* Pr(p|theta,glob_all_frq) */
      for ( int q = 0; q < genodata->get_nbpop (); q++ )
        {
          A += theta[q] * ( glob_all_frq[j][i] - old_p_ji ) * log ( p[q][j][i] );
          A += gsl_sf_lngamma ( theta[q] * old_p_ji ) - gsl_sf_lngamma ( theta[q] * glob_all_frq[j][i] );

          A += theta[q] * ( glob_all_frq[j][k] - old_p_jk ) * log ( p[q][j][k] );
          A += gsl_sf_lngamma ( theta[q] * old_p_jk ) - gsl_sf_lngamma ( theta[q] * glob_all_frq[j][k] );
        }

      /* Prior */
      A += ( lambda-1 ) * ( log ( glob_all_frq[j][k] ) - log ( old_p_jk ) +log ( glob_all_frq[j][i] ) - log ( old_p_ji ) );

      if ( log ( gsl_ran_flat(rng,0,1) ) < A )
        pglob_acc++;
      else
        {
          /* Restaure old values */
          glob_all_frq[j][i] = old_p_ji;
          glob_all_frq[j][k] = old_p_jk;
        }
    }
}

inline void MCMC::init_theta ()
{
  theta_acc = 0;

  theta = new double[genodata->get_nbpop () ];
  if ( !theta )
    exit ( 0 );

  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    theta[q] = 99;
}

inline void MCMC::update_theta ()
{
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    {
      // Save current value
      double old_theta_q = theta[q];

      // Propose a new value
      theta[q] *= exp ( sqrt(s2_theta)*gsl_ran_gaussian(rng,1) );

      // Acceptance ratio
      double A = log ( theta[q] ) - log ( old_theta_q );

      A += log ( old_theta_q ) - log ( theta[q] );
      A -= ( log ( theta[q] ) - log ( old_theta_q ) ) * ( log ( theta[q] ) + log ( old_theta_q ) - 2*omega ) / ( 2*xi );

      for ( int j = 0; j < genodata->get_nbloci (); j++ )
        for ( int k = 1; k < genodata->get_nballel (); k++ )
          if ( genodata->get_is_at_locus ( j, k ) )
            {
              // P(p|glob_all_frq,theta)
              A += gsl_sf_lngamma ( old_theta_q * glob_all_frq[j][k] ) - gsl_sf_lngamma ( theta[q] * glob_all_frq[j][k] );
              A += ( theta[q] - old_theta_q ) * glob_all_frq[j][k] * log ( p[q][j][k] );

            }
      A += genodata->get_nbloci () * ( gsl_sf_lngamma ( theta[q] ) - gsl_sf_lngamma ( old_theta_q ) );

      if ( log ( gsl_ran_flat(rng,0,1) ) < A )
        theta_acc++;
      else
        theta[q] = old_theta_q;
    }
}

double MCMC::MVregupdate ()
{
  int i = 0; /* row index */
  int j = 1; /* col index */

  // Univariate update
  if ( model == 0 )
    {
      nb_reg_param = 1;
      return update_alpha ();
    }

  // Multivariate update
  int nbfact = 1;
  for ( int r = 1; r < factors->get_size(); r++ )
    nbfact += ( ( model >> ( r - 1 ) ) % 2 );
  nb_reg_param = nbfact;

  /* s2_p*I Identity matrix */
  gsl_matrix *invVarCovar = gsl_matrix_alloc ( nbfact, nbfact );
  gsl_matrix_set_identity ( invVarCovar ); // I
  gsl_matrix_scale ( invVarCovar, 1 / s2_alpha ); // (1/s2_p)*I

  /* Matrix of covariates Z_m */
  gsl_matrix *Z_m = gsl_matrix_alloc ( genodata->get_nbpop () * ( genodata->get_nbpop () - 1 ), nbfact );
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    for ( int l = 0; l < genodata->get_nbpop (); l++ )
      if ( l != q )
        {
          gsl_matrix_set ( Z_m, i, 0, factors->get_G(0,q,l) );
          for ( int r = 1; r < factors->get_size(); r++ )
            if ( ( model >> ( r - 1 ) ) % 2 )
              {
                gsl_matrix_set ( Z_m, i, j, factors->get_G(r,q,l) );
                j++;
              }
          j = 1;
          i++;
        }

  /* Inverse of the variance/covariance matrix */
  gsl_matrix *tmp1 = gsl_matrix_alloc ( nbfact, nbfact );
  gsl_linalg_matmult_mod ( Z_m, GSL_LINALG_MOD_TRANSPOSE, Z_m, GSL_LINALG_MOD_NONE, tmp1 ); // tr(Z_m)*Z_m
  gsl_matrix_scale ( tmp1, 1 / s2 ); // (1/s2)*tr(Z_m)*Z_m
  gsl_matrix_add ( invVarCovar, tmp1 ); // inv(Sigma) = (1/s2)*tr(Z_m)*Z_m + (1/s2_p)*I
  gsl_matrix_memcpy ( tmp1, invVarCovar );

  /* Variance/covariance matrix */
  gsl_matrix *VarCovar = gsl_matrix_alloc ( nbfact, nbfact );
  gsl_permutation *perm = gsl_permutation_alloc ( nbfact ); /* Permutation */
  int signum; /* Sign of the permuation */
  gsl_linalg_LU_decomp ( tmp1, perm, &signum ); // LU decomposition
  double det = gsl_linalg_LU_det ( tmp1, signum );
  gsl_linalg_LU_invert ( tmp1, perm, VarCovar ); // Compute inverse from LU decomposition

  gsl_matrix_free ( tmp1 );
  gsl_permutation_free ( perm );

  /* Matrix with log(psi_ql) (Dirichlet prior for migration rates) */
  gsl_matrix *logPsi = gsl_matrix_alloc ( genodata->get_nbpop () * ( genodata->get_nbpop () - 1 ), 1 );
  i = 0;
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    for ( int l = 0; l < genodata->get_nbpop (); l++ )
      if ( l != q )
        {
          gsl_matrix_set ( logPsi, i, 0, log ( psi[q][l] ) );
          i++;
        }

  /* Mean vector */
  gsl_matrix *Mean = gsl_matrix_alloc ( nbfact, 1 );
  gsl_matrix *tmp2 = gsl_matrix_alloc ( nbfact, 1 );
  gsl_linalg_matmult_mod ( Z_m, GSL_LINALG_MOD_TRANSPOSE, logPsi, GSL_LINALG_MOD_NONE, tmp2 ); // tr(Z_m)*logPsi
  gsl_linalg_matmult ( VarCovar, tmp2, Mean ); // VarCovar*tr(Z_m)*logPsi
  gsl_matrix_scale ( Mean, 1 / s2 ); // (1/s2)*VarCovar*tr(Z_m)*logPsi

  // Free memory
  gsl_matrix_free ( logPsi );
  gsl_matrix_free ( Z_m );
  gsl_matrix_free ( tmp2 );

  gsl_matrix *deviate = gsl_matrix_alloc ( nbfact, 1 );
  MVNormaldev ( Mean, VarCovar, deviate );

  gsl_matrix *X = gsl_matrix_alloc ( nbfact, 1 );
  gsl_matrix_memcpy ( X, deviate );
  gsl_matrix_sub ( X, Mean );
  gsl_matrix *tmp3 = gsl_matrix_alloc ( nbfact, 1 );
  gsl_linalg_matmult ( invVarCovar, X, tmp3 );
  gsl_matrix *mat = gsl_matrix_alloc ( 1, 1 );
  gsl_linalg_matmult_mod ( X, GSL_LINALG_MOD_TRANSPOSE, tmp3, GSL_LINALG_MOD_NONE, mat );

  double loglikmod = -0.5*nbfact*log (2*M_PI) + 0.5*log(det) - gsl_matrix_get(mat,0,0)/2.0;

  gsl_matrix_free ( mat );
  gsl_matrix_free ( X );
  gsl_matrix_free ( tmp3 );
  gsl_matrix_free ( Mean );
  gsl_matrix_free ( VarCovar );
  gsl_matrix_free ( invVarCovar );

  alpha[0] = gsl_matrix_get ( deviate, 0, 0 );
  int k = 1;
  for ( int r = 1; r < factors->get_size(); r++ )
    if ( ( model >> ( r - 1 ) ) % 2 )
      {
        alpha[r] = gsl_matrix_get ( deviate, k, 0 );
        k++;
      }
    else
      alpha[r] = 0;

  gsl_matrix_free ( deviate );

  /* Compute mu_qls */
  for ( int q = 0; q < genodata->get_nbpop (); q++ )
    for ( int l = 0; l < genodata->get_nbpop (); l++ )
      if ( l != q )
        {
          mu[q][l] = 0;
          for ( int r = 0; r < factors->get_size(); r++ )
            mu[q][l] += alpha[r] * factors->get_G(r,q,l);
        }

  return loglikmod;
}

inline void MCMC::MVNormaldev ( const gsl_matrix * omega, const gsl_matrix * varcovar, gsl_matrix * vect )
{
  // Cholesky decomposition of the Var/Covar matrix
  gsl_matrix *L = gsl_matrix_alloc ( varcovar->size1, varcovar->size2 );
  gsl_matrix_memcpy ( L, varcovar );
  gsl_linalg_cholesky_decomp ( L );
  for ( size_t i = 0; i < L->size1; i++ )
    for ( size_t j = i + 1; j < L->size2; j++ )
      gsl_matrix_set ( L, i, j, 0 );

  // Independent standard normal variates
  gsl_matrix *Z = gsl_matrix_alloc ( omega->size1, 1 );
  for ( size_t i = 0; i < Z->size1; i++ )
    gsl_matrix_set ( Z, i, 0, gsl_ran_gaussian(rng, 1 ) );

  gsl_linalg_matmult ( L, Z, vect ); // L*Z
  gsl_matrix_add ( vect, omega ); // omega + L*Z

  gsl_matrix_free ( L );
  gsl_matrix_free ( Z );
}

void MCMC::setMCMCData()
{
  setMCMCData(MGinputfile,EFinputfile,useinter);
}

void MCMC::setMCMCData(const string& MGfilename,const string& EFfilename,bool interaction)
{
  try
    {
      if (!genodata)
        genodata = new multilocus_genotypes (MGfilename);

      if (!EFfilename.empty())
        {
          if (!factors)
            factors = new environmental_factors(EFfilename,interaction);
        }
      else
        {
          useReg=false;
          useRJMCMC=false;
        }
      MGinputfile=MGfilename;
      EFinputfile=EFfilename;
    }
  catch (exception &e)
    {
      cerr << "Error: " << e.what() << endl;
      cerr << "Exit" << endl;
      exit(-1);
    }
}

void MCMC::readMCMCSettings(const string& inputfilename)
{
  clog << "Read MCMC settings" << endl;

  ifstream infile;
  infile.open(inputfilename.data());
  if (infile.fail())
    {
      EXCEPTION_INFOS();
      throw OpenFileError();
    }

  string line;
  while (!infile.eof())
    {
      getline(infile,line,'=');

      if ((line.length() >= 10) && (line.substr(line.length()-10,line.length())=="[genodata]"))
        {
          getline(infile,line);
          if (MGinputfile.empty())
            MGinputfile=line;
        }
      else if ((line.length() >= 9) && (line.substr(line.length()-9,line.length())=="[factors]"))
        {
          getline(infile,line);
          if (EFinputfile.empty())
            EFinputfile=line;
        }
      else if ((line.length() >= 13) && (line.substr(line.length()-13,line.length())=="[interaction]"))
        {
          getline(infile,line);
          istringstream info(line);
          if (!factors)
            info >> useinter;
        }
      else if ((line.length() >= 8) && (line.substr(line.length()-8,line.length())=="[burnin]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> burnin;
        }
      else if ((line.length() >= 12) && (line.substr(line.length()-12,line.length())=="[samplesize]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> samplesize;
        }
      else if ((line.length() >= 9) && (line.substr(line.length()-9,line.length())=="[thining]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> thining;
        }
      else if ((line.length() >= 12) && (line.substr(line.length()-12,line.length())=="[regression]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> useReg;
        }
      else if ((line.length() >= 8) && (line.substr(line.length()-8,line.length())=="[RJMCMC]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> useRJMCMC;
        }
      else if ((line.length() >= 8) && (line.substr(line.length()-8,line.length())=="[Fmodel]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> useFmodel;
        }
      else if ((line.length() >= 7) && (line.substr(line.length()-7,line.length())=="[lglik]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> printLglik;
        }
      else if ((line.length() >= 5) && (line.substr(line.length()-5,line.length())=="[psi]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> printPsi;
        }
      else if ((line.length() >= 5) && (line.substr(line.length()-5,line.length())=="[mig]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> printMig;
        }
      else if ((line.length() >= 7) && (line.substr(line.length()-7,line.length())=="[mtild]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> printmtild;
        }
      else if ((line.length() >= 5) && (line.substr(line.length()-5,line.length())=="[reg]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> printReg;
        }
      else if ((line.length() >= 8) && (line.substr(line.length()-8,line.length())=="[popfrq]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> printP;
        }
      else if ((line.length() >= 3) && (line.substr(line.length()-3,line.length())=="[F]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> printF;
        }
      else if ((line.length() >= 5) && (line.substr(line.length()-5,line.length())=="[Fst]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> printFst;
        }
      else if ((line.length() >= 8) && (line.substr(line.length()-8,line.length())=="[ancfrq]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> printPglob;
        }
      else if ((line.length() >= 4) && (line.substr(line.length()-4,line.length())=="[sd]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> printSd;
        }
      else if ((line.length() >= 8) && (line.substr(line.length()-8,line.length())=="[indivq]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> printindivq;
        }
      else if ((line.length() >= 6) && (line.substr(line.length()-6,line.length())=="[popq]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> printpopq;
        }
      else if ((line.length() >= 13) && (line.substr(line.length()-13,line.length())=="[pilotlength]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> pilot_length;
        }
      else if ((line.length() >= 11) && (line.substr(line.length()-11,line.length())=="[pilotruns]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> nb_pilot;
        }
      else if ((line.length() >= 7) && (line.substr(line.length()-7,line.length())=="[a_tau]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> atau;
        }
      else if ((line.length() >= 7) && (line.substr(line.length()-7,line.length())=="[b_tau]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> btau;
        }
      else if ((line.length() >= 10) && (line.substr(line.length()-10,line.length())=="[s2_alpha]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> s2_alpha;
        }
      else if ((line.length() >= 5) && (line.substr(line.length()-5,line.length())=="[Psi]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> Psi;
        }
      else if ((line.length() >= 8) && (line.substr(line.length()-8,line.length())=="[lambda]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> lambda;
        }
      else if ((line.length() >= 7) && (line.substr(line.length()-7,line.length())=="[omega]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> omega;
        }
      else if ((line.length() >= 4) && (line.substr(line.length()-4,line.length())=="[xi]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> xi;
        }
      else if ((line.length() >= 10) && (line.substr(line.length()-10,line.length())=="[s2_theta]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> s2_theta_init;
        }
      else if ((line.length() >= 8) && (line.substr(line.length()-8,line.length())=="[s2_psi]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> s2_psi_init;
        }
      else if ((line.length() >= 5) && (line.substr(line.length()-5,line.length())=="[e_p]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> e_p_init;
        }
      else if ((line.length() >= 5) && (line.substr(line.length()-5,line.length())=="[e_F]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> e_F_init;
        }
      else if ((line.length() >= 5) && (line.substr(line.length()-5,line.length())=="[e_m]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> e_m_init;
        }
      else if ((line.length() >= 6) && (line.substr(line.length()-6,line.length())=="[e_nm]"))
        {
          getline(infile,line);
          istringstream info(line);
          info >> e_nm_init;
        }
    }

  if (MGinputfile.empty())
    {
      cerr << "You must provide at least multilocus genotype input file name" << endl;
      EXCEPTION_INFOS();
      throw FileFormatError();
    }

  clog << "Read MCMC settings ... (done)" << endl;
}

void MCMC::writeMCMCParameters(const string& outputfilename)
{
  ofstream out;
  out.open(outputfilename.data());
  if (out.fail())
    {
      EXCEPTION_INFOS();
      throw OpenFileError();
    }

  out << "[genodata]=" << MGinputfile << endl;

  if (factors)
    {
      out << "[factors]=" << EFinputfile << endl;
      out << "[interaction]=" << factors->get_with_interaction() << endl;
      out << "[nbfactors]=" << factors->get_nb_factors() << endl;
    }

  /* Run settings */
  out << "[burnin]=" << burnin << endl;
  out << "[samplesize]=" << samplesize << endl;
  out << "[thining]=" << thining << endl;
  out << "[runs]=" << nb_mcmc << endl;

  /* Advanced options */
  out << "[regression]=" << useReg << endl;
  out << "[RJMCMC]=" << useRJMCMC << endl;
  out << "[Fmodel]=" << useFmodel << endl;

  /* Output options */
  out << "[lglik]=" << printLglik << endl;
  out << "[psi]=" << printPsi << endl;
  out << "[mig]=" << printMig << endl;
  out << "[mtild]=" << printmtild << endl;
  out << "[reg]=" << printReg << endl;
  out << "[popfrq]=" << printP << endl;
  out << "[F]=" << printF << endl;
  out << "[Fst]=" << printFst << endl;
  out << "[ancfrq]=" << printPglob << endl;
  out << "[sd]=" << printSd << endl;
  out << "[indivq]=" << printindivq << endl;
  out << "[popq]=" << printpopq << endl;

  /* Pilot runs */
  out << "[pilotlength]=" << pilot_length << endl;
  out << "[pilotruns]=" << nb_pilot << endl;

  /* Parameters for prior distributions */
  out << "[a_tau]=" << atau << endl;
  out << "[b_tau]=" << btau << endl;
  out << "[s2_alpha]=" << s2_alpha << endl;
  out << "[Psi]=" << Psi << endl;
  out << "[lambda]=" << lambda << endl;
  out << "[omega]=" << omega << endl;
  out << "[xi]=" << xi << endl;

  /* Proposal parameters */
  out << "[s2_theta]=" << s2_theta_init << endl;
  out << "[s2_psi]=" << s2_psi_init << endl;
  out << "[e_p]=" << e_p_init << endl;
  out << "[e_F]=" << e_F_init << endl;
  out << "[e_m]=" << e_m_init << endl;
  out << "[e_nm]=" << e_nm_init << endl;

  out.close();
}

void MCMC::close_outputfiles()
{
  if ( mcmcfile.is_open() )
    mcmcfile.close();

  if ( mtildfile.is_open() )
    mtildfile.close();

  if ( useFmodel )
    {
      if ( frqfile.is_open())
        frqfile.close();
      if ( fstfile.is_open())
        fstfile.close();
    }

  if ( Fisfile.is_open() )
    Fisfile.close();

  if ( loglikfile.is_open())
    loglikfile.close();

  if ( pfile.is_open())
    pfile.close();

  if ( useReg && psifile.is_open() )
    psifile.close();

  if ( useReg && regfile.is_open() )
    regfile.close();

}
