/***************************************************************************
 *   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 "mcmcthread.h"

QMutex MCMCThread::mutex;

MCMCThread::MCMCThread ( QObject *parent )
    :QThread ( parent )
{
  /* Create mcmc run */
  mcmc = new MCMC();
}
//

MCMCThread::~MCMCThread()
{
  if (isRunning())
    {
      terminate();
      wait();
    }

  delete mcmc;
}

void MCMCThread::run()
{
  /* Initialization */
  mcmc->init();

  /* Pilot run */
  for ( int pilot=0;pilot<MCMC::nb_pilot;pilot++ )
    {
      /* Display pilot run number */
      emit statusChanged ( mcmc->idf,tr ( "Pilot#%1" ).arg ( pilot+1 ) );

      /* MCMC updates */
      for ( int iter=1;iter<=MCMC::pilot_length;iter++ )
        {
          setTerminationEnabled ( false );
          mcmc->update();
          setTerminationEnabled ( true);
        }

      /* Display acceptance rates */
      emit acceptanceChanged ( mcmc->idf,MCMC::pilot_length,mcmc->nm_acc,mcmc->m_acc,mcmc->Mt_acc,mcmc->F_acc,mcmc->p_acc,mcmc->theta_acc,mcmc->pglob_acc,mcmc->psi_acc,mcmc->RJ_acc );

      /* Adjust incremental values */

      /* Non-migrant proportions */
      if ( mcmc->nm_acc <= 0.25*MCMC::pilot_length*MCMC::genodata->get_nbpop() )
        mcmc->e_nm /= 1.2;
      if ( mcmc->nm_acc >= 0.45*MCMC::pilot_length*MCMC::genodata->get_nbpop() )
        mcmc->e_nm *= 1.2;

      /* Immigration rates */
      if ( mcmc->m_acc <= 0.25 * MCMC::pilot_length * MCMC::genodata->get_nbpop () )
        mcmc->e_m /= 1.2;
      if ( mcmc->m_acc >= 0.45 * MCMC::pilot_length * MCMC::genodata->get_nbpop () )
        mcmc->e_m *= 1.2;

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

      /* Allele frequencies */
      if ( mcmc->p_acc <= 0.25*MCMC::pilot_length*MCMC::genodata->get_nbloci() *MCMC::genodata->get_nbpop() )
        mcmc->e_p /= 1.2;
      if ( mcmc->p_acc >= 0.45*MCMC::pilot_length*MCMC::genodata->get_nbloci() *MCMC::genodata->get_nbpop() )
        mcmc->e_p *= 1.2;

      /* F model parameters */
      if ( MCMC::useFmodel )
        {
          if ( mcmc->theta_acc <= 0.25 * MCMC::pilot_length * MCMC::genodata->get_nbpop () )
            mcmc->s2_theta /= 1.5;
          if ( mcmc->theta_acc >= 0.45 * MCMC::pilot_length * MCMC::genodata->get_nbpop () )
            mcmc->s2_theta *= 1.5;

          if ( mcmc->pglob_acc <= 0.25*MCMC::pilot_length*MCMC::genodata->get_nbloci() )
            mcmc->e_p /= 1.2;
          if ( mcmc->pglob_acc >= 0.45*MCMC::pilot_length*MCMC::genodata->get_nbloci() )
            mcmc->e_p *= 1.2;
        }

      /* Inbreeding coefficients */
      if ( mcmc->F_acc <= 0.25 * MCMC::pilot_length * MCMC::genodata->get_nbpop () )
        mcmc->e_F *= 0.5;
      if ( mcmc->F_acc >= 0.45 * MCMC::pilot_length * MCMC::genodata->get_nbpop () )
        mcmc->e_F *= 1.5;

      /* Reset acceptance rates */
      mcmc->nm_acc = 0;
      mcmc->m_acc = 0;
      mcmc->p_acc = 0;
      mcmc->Mt_acc = 0;
      mcmc->F_acc = 0;
      mcmc->psi_acc = 0;
      mcmc->theta_acc = 0;
      mcmc->pglob_acc = 0;
    }

  /* MH algorithm */
  int nbiter = MCMC::burnin + MCMC::samplesize * MCMC::thining;
  int stp = 1;
  /* Display MCMC state */
  emit statusChanged ( mcmc->idf,tr ( "Burnin" ) );
  for ( int iter = 1; iter <= nbiter; iter++ )
    {
      setTerminationEnabled ( false );

      /* Within model moves */
      mcmc->update();

      /* Between model move */
      if ( MCMC::useRJMCMC )
        mcmc->RJupdate ();

      /* Display MCMC state */
      if ( iter == MCMC::burnin+1 )
        emit statusChanged ( mcmc->idf,tr ( "Sampling" ) );

      /* Display progress */
      if ( 100 * iter == stp * 5 * nbiter )
        {
          stp++;
          emit acceptanceChanged ( mcmc->idf,iter,mcmc->nm_acc,mcmc->m_acc,mcmc->Mt_acc,mcmc->F_acc,mcmc->p_acc,mcmc->theta_acc,mcmc->pglob_acc,mcmc->psi_acc,mcmc->RJ_acc );
          mutex.lock();
          emit progressChanged ( mcmc->idf,iter*100.0/nbiter );
          mutex.unlock();
        }

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

          if ( MCMC::useRJMCMC )
            mcmc->visit[mcmc->model]++;

          mcmc->write_iter ( iter );
          emit samplesizeChanged ( mcmc->idf, mcmc->d_assign);

          for ( int h = 0; h < MCMC::genodata->get_nbindiv (); h++ )
            mcmc->post_assign[h][mcmc->Mt[h][1]][mcmc->Mt[h][0]] = ++mcmc->post_assign[h][mcmc->Mt[h][0]][mcmc->Mt[h][1]];

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

              for ( int l = 0; l < MCMC::genodata->get_nbpop (); l++ )
                {
                  mcmc->m_mean[q][l] += mcmc->m[q][l];
                  mcmc->m_sd[q][l] += mcmc->m[q][l] * mcmc->m[q][l];
                }
              for ( int j = 0; j < MCMC::genodata->get_nbloci (); j++ )
                for ( int i = 1; i < MCMC::genodata->get_nballel (); i++ )
                  {
                    mcmc->p_mean[q][j][i] += mcmc->p[q][j][i];
                    mcmc->p_sd[q][j][i] += mcmc->p[q][j][i] * mcmc->p[q][j][i];
                  }
            }
        }
      setTerminationEnabled ( true );
    }
  /* Display MCMC state */
  emit statusChanged ( mcmc->idf,tr ( "Finished" ) );
  mcmc->close_outputfiles();
  mcmc->print();
}

void MCMCThread::kill()
{}
