/* --------------------------------------------------------------------- */
/*                                                                       */
/*           BAYESIAN MULTIMAPPER / INBRED LINE CROSSES  Version 1.0     */
/*                                                                       */
/*           (c) Copyright 1996-2006 by Mikko J. Sillanpaa               */
/*                                                                       */
/*            Rolf Nevanlinna Institute                                  */
/*            FIN-00014 University of Helsinki, FINLAND                  */
/*                                                                       */
/*            E-mail : mjs@rolf.helsinki.fi                              */
/*                                                                       */
/*                  last updated : February 2002                         */
/*                                                                       */
/*                                                                       */
/*                                                                       */
/*            Environmental covariates                                   */
/*                                                                       */
/*                                                                       */
/* --------------------------------------------------------------------- */
/* This file is part of the Multimapper software.                        */ 
/* The Multimapper software has been written only for scientific use.    */ 
/* There is absolutely no warranty in the Multimapper software.          */  
/* You are using the Multimapper software completely at your own risk.   */
/* Beware! this is the version 1.0, the Multimapper software might       */
/* contain some errors or bugs!                                          */
/*                                                                       */
/* --------------------------------------------------------------------- */
/*                                                                       */
/*   file name:       metropolis.c                                       */
/*                                                                       */
/*   This is a program that runs Metropolis-Hastings sampler             */
/*   with a reversible jumps for an inbred line cross                    */
/*   offspring data.                                                     */
/*                                                                       */
/*                                                                       */
/* --------------------------------------------------------------------- */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define _XOPEN_SOURCE 
#include <math.h>
#include <string.h>
#include "data_structure.h"
#include "conversions.h"

#define pi M_PI

/* --------------------------------------------------------------------- */

   double log_likelihood (int N_bc, int Num_ind, regpa *reg,
                      int take[MAX_N_BG_CONT], individuals *first, 
		      int trait,int genotype[3][MAX_N_INDIVID],
                      int bgc_table[MAX_N_INDIVID][MAX_N_BG_CONT],
                      int N_qtl, int num_cov, int cov_index[MAX_N_COV],
                      int cov_type [MAX_N_COV])

/* --------------------------------------------------------------------- */
{
  int i,j,k,index2,q;
  double l1,residual;  
 
  l1=Num_ind*log(1.0/sqrt(2.0*pi*(reg->sigma))); 

  for (i=0;i<Num_ind;i++) 
    {
      residual=first->phenotype[trait]-reg->a;

      for (k=0;k<num_cov;k++)
	if (cov_type[k]==0)
	  residual-=(reg->cov[k][0])*(first->phenotype[cov_index[k]]);
        else
	  residual-=(reg->cov[k][(int) first->phenotype[cov_index[k]]]);
      
      for (q=0;q<N_qtl;q++)  
         residual-=reg->b[q][genotype[q][i]];

      for (j=0;j<N_bc;j++)        
	{ 
   	  index2=bgc_table[i][j];
	  if (take[j]) residual-=reg->bc[j][index2];
	}        
       l1=l1+((-1.0/(2.0*(reg->sigma)))*residual*residual); 
      first=first->next;
    }   
  return(l1); 
}


/* --------------------------------------------------------------------- */

  double ind_log_likelihood (int N_bc, regpa *reg,
			     int take[MAX_N_BG_CONT],double phenotype[10],
                             int genotype[3][MAX_N_INDIVID],
                             int bgc_table[MAX_N_INDIVID][MAX_N_BG_CONT],
                             int individual, int N_qtl,int trait,
			     int num_cov, int cov_index [MAX_N_COV],
			     int cov_type [MAX_N_COV])

/* --------------------------------------------------------------------- */

{
  int j,k,index2,q;
  double l1,residual;  
 
  l1=log(1.0/sqrt(2.0*pi*(reg->sigma))); 

  residual=phenotype[trait]-reg->a;
  
  for (k=0;k<num_cov;k++)
    if (cov_type[k]==0)
      residual-=(reg->cov[k][0])*(phenotype[cov_index[k]]);
    else
      residual-=(reg->cov[k][(int) phenotype[cov_index[k]]]);
  
  for (q=0;q<N_qtl;q++)  
    residual-=reg->b[q][genotype[q][individual]];
      
  for (j=0;j<N_bc;j++)        
    { 
      index2=bgc_table[individual][j];
      if (take[j]) residual-=reg->bc[j][index2];
    }        
  l1=l1+((-1.0/(2.0*(reg->sigma)))*residual*residual); 

  return(l1); 

}

/* --------------------------------------------------------------------- */

void make_take (int take[MAX_N_BG_CONT],int bgc[MAX_N_BG_CONT][2],
                int N_bc, int cur_chr)

/* --------------------------------------------------------------------- */
{
  int i,current;

  current=cur_chr+1;
  for (i=0;i<N_bc;i++)
    if (bgc[i][0]==current) 
      take[i]=0;
    else take[i]=1; 
}


/* --------------------------------------------------------------------- */

void BCross_mrk_prior(int N_ind,int N_markers, 
                      double w[MAX_N_MARKERS-1],
                      double mc_prior[MAX_N_INDIVID],
                      int gene[MAX_N_INDIVID][MAX_N_MARKERS])

/* --------------------------------------------------------------------- */
{
  int i,j;
  double total;

  for (i=0;i<N_ind;i++)  
    { 
      total=1.0;   
      for (j=0;j<(N_markers-1);j++)     
        {
	  if (gene[i][j]!=gene[i][j+1])
	    total *=w[j]; 
          else
	    total *=(1.0-w[j]);            
	}
      mc_prior[i]=total;
    }
}

/* --------------------------------------------------------------------- */

void F2_mrk_prior(int N_ind,int N_markers, 
		  double w[MAX_N_MARKERS-1],
		  double mc_prior[MAX_N_INDIVID],
		  int gene[MAX_N_INDIVID][MAX_N_MARKERS])

/* --------------------------------------------------------------------- */
{
 int i,j;
 short same, homoz_j,homoz_j1;
 double total;

  for (i=0;i<N_ind;i++)  
    { 
      total=1.0;   
      for (j=0;j<(N_markers-1);j++)     
        {
	  same=(gene[i][j]==gene[i][j+1]);
	  homoz_j=(gene[i][j]<1000);
	  homoz_j1=(gene[i][j+1]<1000);
	  
	  if ((same)&&(homoz_j))
	    total *=((1.0-w[j])*(1.0-w[j]));
	  else if ((!same)&&(homoz_j)&&(homoz_j1))
	    total *=(w[j]*w[j]);
	  else if ((same)&&(!homoz_j))
	    total *=(1.0-((2.0*w[j])*(1.0-w[j])));
	  else if ((!same)&&(homoz_j1))
             total *=(w[j]*(1.0-w[j]));	    
	   else 
	    total *=((2.0*w[j])*(1.0-w[j]));	    

	}
      mc_prior[i]=total;
    }
}

/* --------------------------------------------------------------------- */

int flanking_objects (int f, QTL_type q[3], double *R1, double *R2,
		      int N_qtl, int N_ind,int *left_i,int *right_i)
		    
/* --------------------------------------------------------------------- */
/*                                                                       */
/*  returns  0 : no flanking QTLs                                        */
/*           1 : only left flanking QTL                                  */
/*           2 : only right flanking QTL                                 */
/*           3 : both left and right flanking QTL                        */
/*                                                                       */
/* --------------------------------------------------------------------- */
{
  int i, left, right;
  int minl_i, minr_i;
  double min_l, min_r;

   left=0;
   right=0;
   min_l=99999.0;
   min_r=99999.0;

   for (i=0;i<N_qtl;i++)
     if (i<f)
       {
	 if (((q[f].ML==q[i].ML) && (q[f].MR==q[i].MR)) &&
	     (q[i].l_x < q[f].l_x))
	   {
	     left++;
	     if ((q[f].l_x-q[i].l_x)<min_l)
	       {
		 min_l=(q[f].l_x-q[i].l_x);	     
		 minl_i=i;
	       }
	   }

	 if (((q[f].ML==q[i].ML) && (q[f].MR==q[i].MR)) &&
	     (q[f].l_x < q[i].l_x))
	   {
	     right++;
	     if ((q[i].l_x-q[f].l_x)<min_r)
	       {
		 min_r=(q[i].l_x-q[f].l_x);	     
		 minr_i=i;
	       }


	   }
       }

  if ((left+right)==0) 
    return(0);
  else
    {
      if (left>0)
	{
	  *R1=(q[f].r1-q[minl_i].r1)/(1.0-(2.0*q[minl_i].r1));  
	  *left_i=minl_i;
	}

      if (right>0)
	{
	  *R2=(q[f].r2-q[minr_i].r2)/(1.0-(2.0*q[minr_i].r2));  
	  *right_i=minr_i;
	}

      if ((left>0) && (right>0))
	return(3);
      else 
	{
	  if (left>0)
	    return(1);
	  else
	    return(2);
	}
    }

}




/* --------------------------------------------------------------------- */

   void get_regression_proposals(int f, int N_gen, int N_bc, double range[5], 
                                 regpa *regress, int num_cov,
                                 int cov_type[MAX_N_COV])

/* --------------------------------------------------------------------- */
/*  vector regress need to contain last round's accepted parameter       */
/*  values and range contain window sizes for parameters.                */
/*                                                                       */
/* --------------------------------------------------------------------- */
{
  double half,half2,ran1,ran2;
  int i,j,k;

  half=range[2]*0.5;         /* X x 0.5= X / 2.0 */    
  half2=range[3]*0.5;
  ran1=range[2];
  ran2=range[3];

if (f==0)
  {

(*regress).a=(*regress).a-(range[0]*0.5)+((rand()*range[0])/RAND_MAX);       
(*regress).sigma=(*regress).sigma-(range[1]*0.5)+((rand()*range[1])/RAND_MAX);

for (i=0;i<num_cov;i++)
  if (cov_type[i]==0)
    (*regress).cov[i][0]=(*regress).cov[i][0]-
                         (range[0]*0.5)+((rand()*range[0])/RAND_MAX);
  else
    {
      (*regress).cov[i][0]=0.0;
    for (k=1;k<cov_type[i];k++)
      (*regress).cov[i][k]=(*regress).cov[i][k]-
                         (range[0]*0.5)+((rand()*range[0])/RAND_MAX);
    }  

   }
 else if (f==1)
   {

  for(i=0;i<N_gen;i++)
    {
      (*regress).b[0][i]=(*regress).b[0][i]-half+((rand()*ran1)/RAND_MAX);    
      (*regress).b[1][i]=(*regress).b[1][i]-half+((rand()*ran1)/RAND_MAX);    
      (*regress).b[2][i]=(*regress).b[2][i]-half+((rand()*ran1)/RAND_MAX);    
    for(j=0;j<N_bc;j++)
     (*regress).bc[j][i]=(*regress).bc[j][i]-half2+((rand()*ran2)/RAND_MAX);
    }
    

   }

 }


/* --------------------------------------------------------------------- */

   void get_marker_proposal(int N_ind,int N_markers, int N_gen,
                            int gene [MAX_N_INDIVID][MAX_N_MARKERS],
                            int g_possible[MAX_N_GENOTYP],
                            individuals *start, int cur_chr)

/* --------------------------------------------------------------------- */
{
  int i,j,index,proposal,pro,X1,X2;
  geno_type *ptr;


  cur_chr++;
  for (i=0;i<N_ind;i++)
    {    
     ptr=start->start; 
     while(ptr->chr->chromosome_number!=cur_chr) ptr=ptr->next;
     for (j=0;j<N_markers;j++)     
	{
        proposal=ptr->genotype; 
        
        if (proposal<0) 
	  {
	    if (proposal==-1000)
	      {
		index=rand()%N_gen;
		proposal=g_possible[index];                
	      }
	    else
	      {
                while(proposal<0) 
		  {
		    index=rand()%N_gen;
		    pro=g_possible[index];
		    if (pro==-proposal) proposal=pro; 
		    else 
		      {
			X1=(int) pro/1000;
	                X2=pro%1000;
			if (X1==-proposal || X2==-proposal) proposal=pro;
		      }
		  }
	      }      
	  } 
 	gene[i][j]=proposal;	  
        ptr=ptr->next;
      }
    start=start->next;
    }
    
}

/* --------------------------------------------------------------------- */

void change_coding(int N_ind, 
                   int convert [MAX_N_GENOTYP*2][2],
                   individuals *start)

/* --------------------------------------------------------------------- */
/* This function is called only by function readin_gpossible             */
/* --------------------------------------------------------------------- */


{
  int i,k;
  geno_type *ptr;

  for (i=0;i<N_ind;i++)
    {    
    ptr=start->start; 
    while (ptr!=NULL)     
      {
        k=0;
        while(ptr->genotype!=convert[k][0]) k++;                  
	ptr->genotype=convert[k][1];       
        ptr=ptr->next;
       }
    start=start->next;
    }
    
}

/* --------------------------------------------------------------------- */

   void make_genotype(int N_ind, int N_gen, 
                      double p[MAX_N_INDIVID][MAX_N_GENOTYP],
                      int genotype[MAX_N_INDIVID])

/* --------------------------------------------------------------------- */
{
  int i,j,s; 
  double k,sum,add_in[MAX_N_GENOTYP];  

  for (i=0;i<N_ind;i++)
    {
      sum=0.0;     
      for (j=0;j<N_gen;j++)
	{
	  sum+=p[i][j];
	  add_in[j]=sum;  
	}
      k=(((double)rand())/RAND_MAX);    
      s=0;
      while (k>add_in[s]) s++;
      genotype[i]=s;
    }
}


/* --------------------------------------------------------------------- */

   void change_cur (int N_bc, int N_ind, int cur_chr, 
		   int bgc [MAX_N_BG_CONT][2],
		   int bgc_new [MAX_N_INDIVID] [MAX_N_BG_CONT],
		   int bgc_acc [MAX_N_INDIVID] [MAX_N_BG_CONT])

/* --------------------------------------------------------------------- */
{
  int i,j;

  cur_chr++;
  for (i=0;i<N_bc;i++)
      if (bgc[i][0]==cur_chr)   
	{ 
	  for (j=0;j<N_ind;j++)
	    bgc_new[j][i]=bgc_acc[j][i];
	}
}

/* --------------------------------------------------------------------- */

   void update_bgc_new( int N_ind,int N_bc, int cur_chr, 
			int bgc [MAX_N_BG_CONT][2], 
                        int bgc_table1 [MAX_N_INDIVID] [MAX_N_BG_CONT],
                        int gene [MAX_N_INDIVID][MAX_N_MARKERS],
			int a_gene [MAX_N_INDIVID][MAX_N_MARKERS],
                        int g_possible[MAX_N_GENOTYP],
			int operation)

/* --------------------------------------------------------------------- 

    if operation == 1 then convert bgc genotype update to index form   
    
    if operation == 0 then update accepted bgc genotypes to a_gene matrix  



   --------------------------------------------------------------------- */

{
  int i,j,index2,indexi,genocode;

    cur_chr++;
    for (i=0;i<N_bc;i++)
        if (bgc[i][0]==cur_chr) 
	  {
	    index2=(bgc[i][1])-1;

          if (operation==1)
	    {   
	      for (j=0;j<N_ind;j++)
		{
		  genocode=gene[j][index2];
		  indexi=0;
		  while (genocode!=g_possible[indexi]) indexi++; 
		  bgc_table1[j][i]=indexi;
		}
	    }
	  else	   
	    {
	      for (j=0;j<N_ind;j++)
		a_gene[j][index2]=gene[j][index2];
	    }
	  
	  } 
  
}

/* --------------------------------------------------------------------- */

   void make_cro_status (int N_bc, int cur_chr, int Num_marks,
			 int bgc [MAX_N_BG_CONT][2],
			 int cro_status[MAX_N_MARKERS])

/* --------------------------------------------------------------------- */

/*  cro_status vector have zeroes at place of bg and one otherwise       */

/* --------------------------------------------------------------------- */

{
  int i,j;

  for (j=0;j<Num_marks;j++)
    cro_status[j]=1;

  cur_chr++;
   for (i=0;i<N_bc;i++)
     if (bgc[i][0]==cur_chr) 
       cro_status[((bgc[i][1])-1)]=0;      
}

/* --------------------------------------------------------------------- */

 void bg_proposal (int N_ind,int N_bc, int N_gen, 
                   int bgc [MAX_N_BG_CONT][2], 
		   int bgc_table1 [MAX_N_INDIVID] [MAX_N_BG_CONT],
		   int bgc_table2 [MAX_N_INDIVID] [MAX_N_BG_CONT],
                   int g_possible[MAX_N_GENOTYP])

/* --------------------------------------------------------------------- */
{
  int i,j,proposal,index,pro,X1,X2,indexi;
 

  for (i=0;i<N_ind;i++)
    for (j=0;j<N_bc;j++)
      {
        proposal=bgc_table1[i][j];

        if (proposal>=0)
          {
            indexi=0; 
	    while (proposal!=g_possible[indexi]) indexi++;             
	   }
        else 
	  {
	    if (proposal==-1000) 
              {
		indexi=rand()%N_gen;		
              }
            else
	      {
		while(proposal<0)
		  {
		    index=rand()%N_gen;
		    pro=g_possible[index];
		    if (pro==-proposal)
                      {
                        proposal=pro;
                        indexi=index;
                      } 
		    else 
		      {
			X1=(int) pro/1000;
	                X2=pro%1000;
			if (X1==-proposal || X2==-proposal)
                          {
                            proposal=pro;
                            indexi=index;
                          }    
		      }

                   } /* while */

              }  /* else */             

          } /* if */

	bgc_table2[i][j]=indexi;

     } /* for */
}  

/* --------------------------------------------------------------------- */

   void make_bgc_table(int N_ind,int N_bc, int bgc [MAX_N_BG_CONT][2], 
		       int bgc_table [MAX_N_INDIVID] [MAX_N_BG_CONT],
		       individuals *start)

/* --------------------------------------------------------------------- */
{
  int i,j;
  geno_type *ptr;

  for (i=0;i<N_ind;i++)
    {

      for (j=0;j<N_bc;j++)
	{ 
          ptr=start->start;
	  while (bgc[j][0]!=ptr->chr->chromosome_number) ptr=ptr->next;
	  while (bgc[j][1]!=ptr->marker->marker) ptr=ptr->next;        
	  bgc_table[i][j]=ptr->genotype;
	}
      start=start->next;
    } 
}

/* --------------------------------------------------------------------- */

   void readin_randomwalk(int *cur_chr, int *cur_trait, int *maxround, 
                          double *y, double *window, double range[5],
                          double start_point[3], int *permission,
                          int *design, int *num_cov, int cov_index[MAX_N_COV],
                          int cov_type[MAX_N_COV])

/* --------------------------------------------------------------------- */
{
 FILE *filep;
 char line[250],temp[250];
 int i;

  filep=fopen("randomwalk.file","r");
  if (filep==NULL)
    {
      fprintf(stderr,"Cannot find file named randomwalk.file !!!!\n");
      exit(1);
    }

  for (i=0;i<10;i++)
    {
      fgets(line, sizeof(line),filep);
      strncpy(temp,line,9); 
      temp[9]='\0';
      if (i==0) *cur_chr=atoi(temp)-1;
      if (i==1) *cur_trait=atoi(temp)-1;
      if (i==2) *maxround=atoi(temp);
      if (i==3) {
	start_point[0]=atof(temp);                
	strncpy(temp,line+10,9); 
	temp[9]='\0';
        start_point[1]=atof(temp);                
	strncpy(temp,line+20,9); 
	temp[9]='\0';
        start_point[2]=atof(temp);                
     }

      if (i==4) *y=atof(temp); 
      if (i==5) *window=atof(temp); 
      if (i>5) range[i-6]=atof(temp);
    }
      fgets(line, sizeof(line),filep);
      strncpy(temp,line,9); 
      temp[9]='\0';
      *permission=atoi(temp);

      fgets(line, sizeof(line),filep);
      strncpy(temp,line,9); 
      temp[9]='\0';
      *design=atoi(temp);

      fgets(line, sizeof(line),filep);
      strncpy(temp,line,4); 
      temp[4]='\0';
      *num_cov=atoi(temp);     

      for (i=0;i<*num_cov;i++)
	{
	  strncpy(temp,line+5+(i*4),3); 
	  temp[3]='\0';
	  cov_index[i]=atoi(temp);
	}

      for (i=0;i<*num_cov;i++)
	{
	  strncpy(temp,line+5+((*num_cov)*4)+(i*4),3); 
	  temp[3]='\0';
	  cov_type[i]=atoi(temp);	  
	}

  fclose(filep);
  printf("now randomwalk.file succesfully read in ... \n");
}

/* --------------------------------------------------------------------- */
   
   void readin_priorlimits(double limit[4][2],double lx_limit[2])

/* --------------------------------------------------------------------- */
{
  FILE *filep;
  char line[250],temp[250];
  int i;

  filep=fopen("priorlimits.file","r");
  if (filep==NULL)
    {
      fprintf(stderr,"Cannot find file named priorlimits.file !!!!\n");
      exit(1);
    }

  for (i=0;i<5;i++)
    {
      fgets(line, sizeof(line),filep);
      strncpy(temp,line,9); 
      temp[9]='\0';
      if (i<4) limit[i][0]=atof(temp);
      else lx_limit[0]=atof(temp);
      strncpy(temp,line+10,9); 
      temp[9]='\0';
      if (i<4) limit[i][1]=atof(temp);
      else lx_limit[1]=atof(temp);      
    }
  fclose(filep);
  printf("now priorlimits.file succesfully read in ... \n");
}


/* --------------------------------------------------------------------- */

   void readin_bg_controls(int N_bc,int bgc[MAX_N_BG_CONT][2])

/* --------------------------------------------------------------------- */
{
  FILE *filep;
  char line[250],temp[250];
  int i;

  filep=fopen("bg_controls.file","r");
  if (filep==NULL)
    {
      fprintf(stderr,"Cannot find file named bg_controls.file !!!!\n");
      exit(1);
    }

  for (i=0;i<N_bc;i++)
     { 
       fgets(line, sizeof(line),filep);
       strncpy(temp,line,4); 
       temp[4]='\0';
       bgc[i][0]=atoi(temp); 
       strncpy(temp,line+5,4); 
       temp[4]='\0';
       bgc[i][1]=atoi(temp); 
     }
  fclose(filep); 
  printf("now bg_controls.file succesfully read in ... \n");
}

/* --------------------------------------------------------------------- */
    
   void readin_g_possible(int g_possible[MAX_N_GENOTYP],int *N_gen,
                          int N_ind, individuals *start)

/* --------------------------------------------------------------------- */
/*                                                                       */
/*  Read in codesystem file : codesystem.file and                        */
/*  make g_possible vector containing all possible genotypes.            */
/*                                                                       */
/*    codes for homozygotes      :   1,...,N_homoz                       */
/*                                                                       */
/*    code for heterozygote (ij) :   i*1000+j                            */
/*                                                                       */
/* --------------------------------------------------------------------- */
{
  FILE *filep;
  char line[250],temp[250],temp_1,temp_2;
  int i,temp2,out,out2,convert[MAX_N_GENOTYP*2][2];

  filep=fopen("codesystem.file","r");
  if (filep==NULL)
    {
      fprintf(stderr,"Cannot find file named codesystem.file !!!!\n");
      exit(1);
    }
  out=0;
  *N_gen=0;
  while (!out)
    {
      fgets(line, sizeof(line),filep);
      strncpy(temp,line,2); 
      temp[2]='\0';
      if (strchr(temp,'-')!=NULL) out=1;
      else
	{
	  strncpy(temp,line+3,4);
	  temp[4]='\0';
          convert[*N_gen][0]=atoi(temp);

	  strncpy(temp,line+8,6);
	  temp[6]='\0';
          temp2=atoi(temp);
	  g_possible[*N_gen]=temp2;
          convert[*N_gen][1]=temp2;
          (*N_gen)++;
	}
    }
  out2=0;
  i=0;
  while (!out2)
    {
      temp_1=line[0];
      temp_2=line[1];
      if (temp_1=='-' && temp_2=='-') out2=1;
      strncpy(temp,line+3,4);
      temp[4]='\0';
      convert[*N_gen+i][0]=atoi(temp);
      strncpy(temp,line+8,6);
      temp[6]='\0';
      convert[(*N_gen)+i][1]=atoi(temp);
      i++;
      fgets(line, sizeof(line),filep);
    }

  fclose(filep);
  printf("now codesystem.file succesfully read in ... \n");


/* for debugging ------------------------ */
/*  
  for (j=0;j<*N_gen;j++) 
    printf("convert[%d,0.1]= %d, %d \n",j,convert[j][0],convert[j][1]);
  for (j=0;j<i;j++) 
    printf("convert[%d,0.1]= %d, %d \n",
    j+(*N_gen),convert[j+(*N_gen)][0],convert[j+(*N_gen)][1]); */

/* --------------------------------------- */

  change_coding(N_ind,convert,start);
}


/* --------------------------------------------------------------------- */

   void make_p_ratio(int N_ind, int genotype[MAX_N_INDIVID],
                     int genotype1[MAX_N_INDIVID],
                     double p[MAX_N_INDIVID][MAX_N_GENOTYP], 
                     double p1[MAX_N_INDIVID][MAX_N_GENOTYP],
                     double *ratio)

/* --------------------------------------------------------------------- */
{
  int i; 
  double prod;

  prod=1.0;
  for (i=0;i<N_ind;i++)
        prod=prod*(p[i][genotype[i]]/p1[i][genotype1[i]]);
  *ratio=prod;
}


/* --------------------------------------------------------------------- */

void cp_field(int f, int N_qtl, QTL_type q_new[3],
	      QTL_type q[3], QTL_type q_propo)

/* --------------------------------------------------------------------- */
{
  int i;

  for (i=0;i<N_qtl;i++)	  
    if (i!=f)
      {
	q_new[i].l_x=q[i].l_x;
	q_new[i].ML=q[i].ML;
	q_new[i].MR=q[i].MR;
	q_new[i].g_l=q[i].g_l;               
	q_new[i].r1=q[i].r1;   
	q_new[i].r2=q[i].r2;                  	
      }
    else
      {
	q_new[f].l_x=q_propo.l_x;
	q_new[f].ML=q_propo.ML;
	q_new[f].MR=q_propo.MR;
	q_new[f].g_l=q_propo.g_l;               
	q_new[f].r1=q_propo.r1;   
	q_new[f].r2=q_propo.r2;                  	
      }

}

/* --------------------------------------------------------------------- */
   void x_prior_distribution (int N_ind, QTL_type q, int N_gen, 
			      int g_possible[MAX_N_GENOTYP],
                              double p[MAX_N_INDIVID][MAX_N_GENOTYP], 
                              int gene[MAX_N_INDIVID][MAX_N_MARKERS],
                              int design, int flank_type,
			      int left_i, int right_i, double R1, double R2,
			      int QTL_geno [3][MAX_N_INDIVID])

/* --------------------------------------------------------------------- */
/*                                                                       */
/*  Returns matrix p which contains prior distribution for different     */
/*  genotypes of QTL locus conditional to each individuals complete      */
/*  flanking marker genotypes.                                           */
/*                                                                       */
/* --------------------------------------------------------------------- */

{
  int i,j,gen,X1,X2,ML1,ML2,MR1,MR2,al1,al2,alleleshare;
  int COMB,k,Q1,Q2,M1,M2,QTL_heterozygote,RFM_heterozygote;  
  double p1,p2,p3,p4,total[MAX_N_INDIVID];

   if (design==1)
     {
/* --------------------------------------------------------*/              
/*          INBRED LINE CROSS , BACKCROSS (P X F1)         */
/* --------------------------------------------------------*/
    
   for (j=0;j<N_ind;j++)
     total[j]=0.0;

    for (i=0;i<N_gen;i++)
      {
       
	gen=g_possible[i];
	if (gen<1000) 
	  {
	    X1=gen;
	    X2=gen;
	  }
	else 
	  {
	    X1=(int) gen/1000;
	    X2=gen%1000;
	  }

	for (j=0;j<N_ind;j++)
          {
	    
       if (flank_type==0)       /* no flanking QTLs                 */
	 {
	    Q1=gene[j][q.ML];
	    Q2=gene[j][q.MR];
	    R1=q.r1;
	    R2=q.r2;
	 }
	else 
	  if (flank_type==1 )   /* only left flanking QTL           */
	    {
	      Q1=QTL_geno[left_i][j];
	      Q2=gene[j][q.MR];
	      R2=q.r2;
	    }
	  else 
	    if (flank_type==2)  /* only right flanking QTL          */
	      {
		Q1=gene[j][q.ML];
		Q2=QTL_geno[right_i][j];
		R1=q.r1;
	      }
	    else                /* both left and right flanking QTL */
	      {
		Q1=QTL_geno[left_i][j];
		Q2=QTL_geno[right_i][j];
	      }

       
       if (Q1<1000)
	 {
	   ML1=Q1;
	   ML2=Q1;
	 }
       else 
	 {
	   ML1=(int) Q1/1000;
	   ML2=Q1%1000;
	 }
	    
       if (Q2<1000)
	 {
	   MR1=Q2;
	   MR2=Q2;
	 }
       else 
	 {
	   MR1=(int) Q2/1000;
	   MR2=Q2%1000;
	 }
       
       al1=(X1==ML1)+(X2==ML2);
       al2=(X1==ML2)+(X2==ML1);

       if (al1>al2) alleleshare=al1;
       else alleleshare=al2;

       if (alleleshare==2) p1=1.0-R1;
       else p1=R1;

       al1=(X1==MR1)+(X2==MR2);
       al2=(X1==MR2)+(X2==MR1);

       if (al1>al2) alleleshare=al1;
       else alleleshare=al2;

       if (alleleshare==2) p2=1.0-R2;
       else p2=R2;

       p[j][i]=p1*p2;
     
       total[j]=total[j]+p[j][i];
	  
	  } /* for j */ 
      } /* for i */
     
    for (j=0;j<N_ind;j++)
     for (i=0;i<N_gen;i++)
       p[j][i]=p[j][i]/total[j];
 
 
     
 /*   printf("r1= %f \n",q.r1);
    printf("r2= %f \n",q.r2);
    printf("total[1]= %f \n",total[1]);
    printf("p[1][0]= %f \n",p[1][0]);
    printf("p[1][1]= %f\n",p[1][1]);
    printf("gene[1][q.ML]= %d\n",gene[1][q.ML]);
    printf("gene[1][q.MR]= %d\n",gene[1][q.MR]);

    printf("total[15]= %f \n",total[15]);
    printf("p[15][0]= %f \n",p[15][0]);
    printf("p[15][1]= %f \n",p[15][1]);
    printf("gene[15][q.ML]= %d\n",gene[15][q.ML]);
    printf("gene[15][q.MR]= %d\n",gene[15][q.MR]);

    printf("total[19]= %f \n",total[19]);
    printf("p[19][0]= %f \n",p[19][0]);
    printf("p[19][1]= %f \n",p[19][1]);
    printf("gene[19][q.ML]= %d\n",gene[19][q.ML]);
    printf("gene[19][q.MR]= %d\n",gene[19][q.MR]);   */

     }
else if (design==2) 
  {
/* --------------------------------------------------------*/
/*     INBRED LINE CROSS , F2 INTERCROSS (F1 x F1)         */
/* --------------------------------------------------------*/
    for (j=0;j<N_ind;j++)
      total[j]=0.0;

    for (i=0;i<N_gen;i++)
      {
       
	gen=g_possible[i];
	if (gen<1000) 
	  {
	    X1=gen;
	    X2=gen;
	  }
	else 
	  {
	    X1=(int) gen/1000;
	    X2=gen%1000;
	  }

	for (j=0;j<N_ind;j++)
          {

	    
       if (flank_type==0)       /* no flanking QTLs                 */
	 {
	    al1=gene[j][q.ML];
	    al2=gene[j][q.MR];
	    R1=q.r1;
	    R2=q.r2;
	 }
	else 
	  if (flank_type==1 )   /* only left flanking QTL           */
	    {
	      al1=QTL_geno[left_i][j];
	      al2=gene[j][q.MR];
	      R2=q.r2;
	    }
	  else 
	    if (flank_type==2)  /* only right flanking QTL          */
	      {
		al1=gene[j][q.ML];
		al2=QTL_geno[right_i][j];
		R1=q.r1;
	      }
	    else                /* both left and right flanking QTL */
	      {
		al1=QTL_geno[left_i][j];
		al2=QTL_geno[right_i][j];
	      }


            if (al1<1000)
	      {
                ML1=al1;
                ML2=al1;
	      }
            else 
              {
		ML1=(int) al1/1000;
		ML2=al1%1000;
              }
	    
	    if (al2<1000)
	      {
                MR1=al2;
                MR2=al2;
	      }
            else 
              {
		MR1=(int) al2/1000;
		MR2=al2%1000;
              }
                   
            QTL_heterozygote=(X1!=X2);   

            RFM_heterozygote=(MR1!=MR2);

            if (QTL_heterozygote) COMB=2;
            else COMB=1;

            if (RFM_heterozygote) COMB=COMB*2;

            p[j][i]=0.0; 

            Q1=X1;
            Q2=X2;
            M1=MR1;
            M2=MR2;

            for (k=0;k<COMB;k++)
	      {
		
		if ((QTL_heterozygote)&&((k==1)||(k==3))) 
                  {
		    al1=Q1;
                    Q1=Q2;
                    Q2=al1;		    
		  }
              
                if (((!QTL_heterozygote)&&(k==1))|| 
                    ((RFM_heterozygote)&&(k==2))) 
		    {
		     al1=M1;
		     M1=M2;
                     M2=al1; 
		    }

		if (Q1==ML1) p1=(1.0-R1);
		else p1=R1;
		
		if (Q1==M1) p2=(1.0-R2);
		else p2=R2;
	    
		if (Q2==ML2) p3=(1.0-R1);
		else p3=R1;
	    
		if (Q2==M2) p4=(1.0-R2);
		else p4=R2;

		p[j][i]+=(p1*p2*p3*p4);
	      }

	    total[j]=total[j]+p[j][i];
	  }
      }

    for (j=0;j<N_ind;j++)
      for (i=0;i<N_gen;i++)
	p[j][i]=p[j][i]/total[j];

  }

}

/* --------------------------------------------------------------------- */
 
    void vec_cp(int N_ind,int acc_genotype[MAX_N_INDIVID],
                int genotype[MAX_N_INDIVID])

/* --------------------------------------------------------------------- */
{
  int i;
  
  for(i=0;i<N_ind;i++)
    acc_genotype[i]=genotype[i];
}

/* --------------------------------------------------------------------- */

   void mat_cp(int N_ind, int N_markers,
               int acc_gene[MAX_N_INDIVID][MAX_N_MARKERS], 
               int gene[MAX_N_INDIVID][MAX_N_MARKERS])

/* --------------------------------------------------------------------- */
{
  int i,j;

  for(i=0;i<N_markers;i++)
	for (j=0;j<N_ind;j++)
	  acc_gene[j][i]=gene[j][i];

}

/* --------------------------------------------------------------------- */
 
 void cp_regress(int N_gen, int N_bc, regpa *target, regpa *source, int N_qtl,
                 int num_cov, int cov_type[MAX_N_COV])

/* --------------------------------------------------------------------- */
{
  int i,j,k;

  (*target).a=source->a;
  (*target).sigma=source->sigma;
  for(i=0;i<N_gen;i++)
     {
      for(j=0;j<N_qtl;j++)
       (*target).b[j][i]=source->b[j][i];

       for(j=0;j<N_bc;j++) 
          (*target).bc[j][i]=source->bc[j][i];
     }

  for (i=0;i<num_cov;i++)
    {
      (*target).cov[i][0]=source->cov[i][0];
      for (k=1;k<cov_type[i];k++)
	(*target).cov[i][k]=source->cov[i][k];
    }

}

/* --------------------------------------------------------------------- */

   void update_bgc(int N_ind, int N_bc, 
                   int target[MAX_N_INDIVID][MAX_N_BG_CONT],
                   int source[MAX_N_INDIVID][MAX_N_BG_CONT])

/* --------------------------------------------------------------------- */
{
  int i,j;
  for(i=0;i<N_ind;i++)
    for(j=0;j<N_bc;j++)
      target[i][j]=source[i][j];

}

/* --------------------------------------------------------------------- */

   void print_start_of_file( FILE *fp, int cur_chr, int max_round,
                             double limit [5][2], double range[5],
			     double window,double y)

/* --------------------------------------------------------------------- */
{

  fprintf(fp,"-------------------------------------------------------\n");
  fprintf(fp,"For chromosome : %9d\n",(cur_chr));
  fprintf(fp,"Number of Metropolis-Hastings rounds        : %9d\n",max_round);
  fprintf(fp,"proposal probability for add/delete a QTL : 1/%f\n",y);
  fprintf(fp,"-------------------------------------------------------\n");
  fprintf(fp,"priorlimits : \n");
  fprintf(fp,"for a       = [%9.2f %9.2f]\n",limit[0][0],limit[0][1]);
  fprintf(fp,"    sigma^2 = [%9.2f %9.2f]\n",limit[1][0],limit[1][1]);
  fprintf(fp,"    b       = [%9.2f %9.2f]\n",limit[2][0],limit[2][1]);
  fprintf(fp,"    bc      = [%9.2f %9.2f]\n",limit[3][0],limit[3][1]);
  fprintf(fp,"-------------------------------------------------------\n");
  fprintf(fp,"proposal ranges : \n");
  fprintf(fp,"for l_x            = %9.2f\n",window);
  fprintf(fp,"    a, covariates  = %9.2f\n",range[0]);
  fprintf(fp,"    sigma^2        = %9.2f\n",range[1]);
  fprintf(fp,"    b              = %9.2f\n",range[2]);
  fprintf(fp,"    bc             = %9.2f\n",range[3]);
  fprintf(fp,"-------------------------------------------------------\n");

}

/* --------------------------------------------------------------------- */

void metropolis_hastings (chromosome_type *first, int Num_chr,
                          int Num_markers[MAX_N_CHROMOS], individuals *start,
                          int N_ind, int N_bc)

/* --------------------------------------------------------------------- */
{
  FILE *fp,*fp0,*fp1,*fp2,*fp3,*fp5, *fr0, *fr1, *fr2, *fr3, *fr20;   
  int i,j,m,g_possible[MAX_N_GENOTYP],N_gen,fail,round,rejected[3],
      max_round; /*,s */ 
  int acc_geno[3][MAX_N_INDIVID],genotype[3][MAX_N_INDIVID],possi,in;
  int cur_chr,trait,Num_marks,take[MAX_N_BG_CONT];
  int num_reject[3][1],num_rej[2], /*last_stage,*/ design;
  int bgc[MAX_N_BG_CONT][2], bgc_table[MAX_N_INDIVID][MAX_N_BG_CONT],
      bgc_new[MAX_N_INDIVID][MAX_N_BG_CONT],cro_status[MAX_N_MARKERS],
      bgc_acc[MAX_N_INDIVID][MAX_N_BG_CONT],permission,mup,mup2,mup3[3],
      prop_mup,f; 
  int a_gene[MAX_N_INDIVID][MAX_N_MARKERS],gene[MAX_N_INDIVID][MAX_N_MARKERS];
  double cumul[MAX_N_MARKERS],lenght,location,acc_loca[3],window,y,l1,l0,
         range[5];
  double p[MAX_N_INDIVID][MAX_N_GENOTYP],limit[4][2];
  double pp[3][MAX_N_INDIVID][MAX_N_GENOTYP],
         pp1[3][MAX_N_INDIVID][MAX_N_GENOTYP];
  double p1[MAX_N_INDIVID][MAX_N_GENOTYP],ratio,dist,temp_d; 
  double lx_limit[2],chance,l01,l11,given_limit;
  markers *ptr, *ptr_table[MAX_N_MARKERS];  
  QTL_type q[3],q_propo, q_new[3];
  regpa regress,reg_acc;
  individuals *start_1;
  double p_del,p_add,lambda; /* mem[3][MAX_N_GENOTYP]; */
  int N_qtl,add_one,del_one,move_all,which_one,how_many;
  double mc_prior[MAX_N_INDIVID],acc_prior[MAX_N_INDIVID];
  double wg[MAX_N_MARKERS-1],mean_effect;

  int num_cov, cov_index[MAX_N_COV], cov_type[MAX_N_COV];

  int flank_type, left_i, right_i, calc_al1,calc_al2,calc_add,calc_del;
  double R1,R2;
  double true_inter;

         
  /*  acc_loca[1]=50.0; acc_loca[2]=80.0; */

  lambda=2;
  N_qtl=3;

  permission=0;
  round=0;
  num_reject[0][0]=0;  
  num_reject[1][0]=0;  
  num_reject[2][0]=0;
  how_many=0;  
  /*  num_reject[0][1]=0; */  
  num_rej[0]=0;
  num_rej[1]=0;
  mup=0;
  mup2=0;  
  mup3[0]=0;
  mup3[1]=0;
  mup3[2]=0;
  prop_mup=0;
  calc_al1=0;
  calc_al2=0; 
  calc_add=0; 
  calc_del=0;

/*------ read in codesystem.file  -------------- */   
  readin_g_possible(g_possible,&N_gen,N_ind,start);      

  readin_priorlimits(limit,lx_limit);    /* read in priorlimits.file */
  readin_bg_controls(N_bc,bgc);              /* read in bg_controls.file */

  make_bgc_table(N_ind,N_bc,bgc,bgc_table,start);

/*------ read in randomwalk.file -------------- */
readin_randomwalk(&cur_chr,&trait,&max_round,&y,&window,range,acc_loca,
		  &permission,&design, &num_cov, cov_index, cov_type); 

Num_marks=Num_markers[cur_chr]; 

printf("acc loca 0= %f\n",acc_loca[0]);
printf("acc loca 1= %f\n",acc_loca[1]);
printf("acc loca 2= %f\n",acc_loca[2]);


  given_limit=window; 

 fp2=fopen("MH_output.logi","w");   

/* initial values :  */

  get_marker_proposal (N_ind,Num_marks,N_gen,a_gene,g_possible,start,cur_chr);
  

  reg_acc.a=(limit[0][0]+limit[0][1])*0.5;
  reg_acc.sigma=(limit[1][0]+limit[1][1])*0.5;
  for (i=0;i<N_gen;i++)
    {
      reg_acc.b[0][i]=(limit[2][0]+limit[2][1])*0.5;
      reg_acc.b[1][i]=(limit[2][0]+limit[2][1])*0.5;
      reg_acc.b[2][i]=(limit[2][0]+limit[2][1])*0.5;
    }
  for (i=0;i<N_gen;i++)
    for (j=0;j<N_bc;j++)
      reg_acc.bc[j][i]=(limit[3][0]+limit[3][1])*0.5;      

  for (i=0;i<num_cov;i++)
    for (j=0;j<MAX_N_CLASSES;j++)
      reg_acc.cov[i][j]=0.0;
     
  get_regression_proposals(0,N_gen,N_bc,range,&reg_acc, num_cov, cov_type);
  get_regression_proposals(1,N_gen,N_bc,range,&reg_acc, num_cov, cov_type);
    
  bg_proposal (N_ind,N_bc,N_gen,bgc,bgc_table,bgc_acc,g_possible);

  make_cro_status(N_bc,cur_chr,Num_marks,bgc,cro_status);

  make_take(take,bgc,N_bc,cur_chr);  
 		 
/*
       for (i=0;i<20;i++)
        for (j=0;j<N_bc;j++)
	  {
	  if (i==0) printf(" bgc[%d,0. .1]= %d  %d\n",j,bgc[j][0],bgc[j][1]);
	  printf(" bgc_table[%d,%d]= %d\n",i,j,bgc_table[i][j]);      
          printf(" bgc_acc[%d,%d]= %d\n",i,j,bgc_acc[i][j]);      
	}
*/       
       
    fprintf(fp2," Starting values for regression parameters:\n");
    fprintf(fp2,"--------------------------------------------\n");  
    fprintf(fp2," reg_acc.a     %f\n",reg_acc.a);
    fprintf(fp2," reg_acc.sigma %f\n",reg_acc.sigma);
     for (i=0;i<N_gen;i++)
       {
      fprintf(fp2," reg_acc.b[0][%d] %f\n",i,reg_acc.b[0][i]);
      fprintf(fp2," reg_acc.b[1][%d] %f\n",i,reg_acc.b[1][i]);
      fprintf(fp2," reg_acc.b[2][%d] %f\n",i,reg_acc.b[2][i]);

      for (j=0;j<N_bc;j++)
        fprintf(fp2," reg_acc.bc[%d][%d] %f\n",j,i,reg_acc.bc[j][i]);
       }
    fprintf(fp2,"--------------------------------------------\n");  

  fprintf(fp2," Number of individuals %d\n",N_ind);
  fprintf(fp2," Number of BG controls %d\n",N_bc);
  if (design==1)
    fprintf(fp2," Design 1 => Backcross (%d genotypes)\n",N_gen);
  else if (design==2)
    fprintf(fp2," Design 2 => F2 (%d genotypes)\n",N_gen);
  else
    fprintf(fp2," Design %d => UNKNOWN DESIGN (%d genotypes)\n",design,N_gen);

  fprintf(fp2,"----------------------------\n");
  fprintf(fp2," Background controls :\n");

  for (i=0;i<N_bc;i++)
     fprintf(fp2,"CR  %d  MARKER  %d\n",bgc[i][0],bgc[i][1]);

 fprintf(fp2,"Number of covariates = %d\n",num_cov);
 for (i=0;i<num_cov;i++)
   if (cov_type[i]==0)
     	fprintf(fp2,"  Type of cov # %d = Regression                  id=%d\n",i,cov_index[i]);
      else 
        fprintf(fp2,"  Type of cov # %d = Classification (%d classes)  id=%d\n",i,cov_type[i],cov_index[i]);

  fprintf(fp2,"----------------------------\n");
  fprintf(fp2," possible genotypes :\n");
  possi=0;
  for (i=0;i<N_gen;i++)
    {
      fprintf(fp2,"%d\n",g_possible[i]);
      if (g_possible[i]>1000) possi=1;
    }
  if (possi==0) 
    { 
     fprintf(fp2," !!!!!! you have only made homozygote genotypes           !!!!\n");
     fprintf(fp2," !!!!!! please check that your codesystem.file is correct !!!!\n");
     }
  fprintf(fp2,"----------------------------\n");


  srand(clock());                              /* shuffle rand generator */
  fp=fopen("MH_output__lx0","w");
  fp0=fopen("MH_output__lx1","w");
  fp1=fopen("MH_output__lx2","w");

  fp3=fopen("MH_output__logL","w");

  /*  fp4=fopen("MH_output__regr","w"); */

  fp5=fopen("MH_output__Nqtl","w"); 

  fr0=fopen("MH_output__reg_sigma_a","w"); 
  fr1=fopen("MH_output__reg0","w"); 
  fr2=fopen("MH_output__reg1","w"); 
  fr3=fopen("MH_output__reg2","w"); 

  fr20=fopen("MH_output__env","w");
/* 
  print_start_of_file(fp,(cur_chr+1),max_round,limit,range,window);
  print_start_of_file(fp3,(cur_chr+1),max_round,limit,range,window);
  */

   print_start_of_file(fp2,(cur_chr+1),max_round,limit,range,window,y);

/*     print_start_of_file(fp4,(cur_chr+1),max_round,limit,range,window); */


  
  fflush(NULL);                          /* for debugging */

/* ------------------------------------------------------------------- */
/*                                                                     */
/*  Here we calculate current chromosomal lenght                       */
/*                                                                     */
/* ------------------------------------------------------------------- */     

  lenght=0.0;
  while(cur_chr!=(first->chromosome_number-1)) first=first->next;
  ptr=first->start;
  i=0;
  cumul[0]=0.0;
  while (ptr!=NULL) {
    lenght+=ptr->r; 
    ptr_table[i]=ptr;
    cumul[++i]=lenght;
    ptr=ptr->next; 
  }
  cumul[i+1]=0.0;
  fprintf(fp2,"Lenght of chromosome #%d is %f cM\n",cur_chr+1,lenght);  
 
  if (window>=lenght) 
    { 
     printf("------------------------------------------------\n"); 
     printf("ERROR : CHECK YOUR PROPOSAL RANGES!!!!!!!\n"); 
     printf("PROPOSAL RANGE FOR QTL LOCATION PARAMETERS (%f)\n",window); 
     printf("IS HIGHER OR EQUIVALENT TO THE LENGTH OF\n"); 
     printf("THE CHROMOSOME #%d (%f)\n",cur_chr+1,lenght); 
     printf("------------------------------------------------\n"); 
     exit(1); 
    } 

  for (i=0;i<3;i++)
    if(acc_loca[i]>lenght)
      {
     printf("------------------------------------------------\n");
     printf("ERROR: INITIAL LOCATION OUT OF MARKER MAP!!!!\n");
     printf("PLEASE MODIFY INITIAL LOCATIONS IN FILE : randomwalk.file\n");
     printf("Lenght of chromosome #%d is %f cM\n",cur_chr+1,lenght);  
     printf("------------------------------------------------\n");
     exit(1);
      }

  for (i=0;i<12;i++)
    printf("cumul[%d]=%f\n",i,cumul[i]);

   
  if (lx_limit[0]<0) 
    lx_limit[0]=0; 

  if (lx_limit[1]>lenght)
    lx_limit[1]=lenght;
     

/* ------------------------------------------------------------------- */

  /* acc_loca=lenght*0.5; */                          /* starting value */
          


   p_add=1.0/y;
   p_del=p_add;   

  
   /*  for initialization */
      
   for (i=0;i<(Num_marks-1);i++)
     {
       wg[i]=Haldane_to_recomb(ptr_table[i]->r);

         if (wg[i]==0.0) 
	 {
        printf("------------------------------------------------\n");
        printf("ERROR : DISTANCE BETWEEN TWO MARKERS IS 0.0 cM !!!!!\n");
	printf("PLEASE DELETE ONE MARKER OR REPLACE ZERO DISTANCE\n");
        printf("WITH SOME SMALL NON-ZERO VALUE IN <filename>.map\n");
	printf("THIS ERROR OCCURRED IN CHROMOSOME %d\n",cur_chr+1);
	printf("------------------------------------------------\n");
	exit(1);
	 }
     }

   
   if (design==1) 
     BCross_mrk_prior(N_ind,Num_marks,wg,acc_prior,a_gene);
   else
     F2_mrk_prior(N_ind,Num_marks,wg,acc_prior,a_gene);

  /*----------------( M-H iteration starts from here )----------------*/ 

   while (round<max_round)
     {

  if (round%1000==0) printf("round # %d\n",round);



/* ------------------------------------------------------------------- */
/*                                                                     */
/*  Here we choice one of 3 operations                                 */
/*                                                                     */
/*   1) delete one of existing QTLs, if N_qtl>0                        */
/*   2) add one new QTL, if N_qtl<3                                    */
/*   3) propose new location(s) of current QTL(s)                      */
/*                                                                     */
/* ------------------------------------------------------------------- */

  
    chance=(((double)rand())/RAND_MAX);
    add_one=(chance<p_add);
    del_one=((!add_one)&&(chance<(p_add+p_del)));
    move_all=(((!add_one)&&(!del_one))||((add_one) && (N_qtl==3))||
               ((del_one) && (N_qtl==0)));


   if ((!move_all) && (round!=0))
     {
      if (add_one)
        {
          /* add one QTL to the model*/

	  location=((rand()*lenght)/RAND_MAX);     /* proposal */

    	  in=0;
	  while (cumul[in]<location) in++;
	  if (cumul[in]>location) in--;
          if (location==lenght) in--;       
     
	  q_propo.l_x=location;
	  q_propo.ML=in;
	  q_propo.MR=in+1;
  	  q_propo.g_l=ptr_table[in];               
          dist=location-cumul[in];
	  q_propo.r1=Haldane_to_recomb(dist);

	  dist=Haldane_to_recomb(q_propo.g_l->r);
          q_propo.r2=(dist-q_propo.r1)/(1.0-(2.0*q_propo.r1));
  	  
          if (location==lenght) 
      	    {
               temp_d=q_propo.r1;
               q_propo.r1=q_propo.r2; 
               q_propo.r2=temp_d;      
	    }

    	cp_field(N_qtl,(N_qtl+1),q_new,q,q_propo);

	flank_type=flanking_objects(N_qtl,q_new,&R1,&R2,(N_qtl+1),
				    N_ind, &left_i, &right_i);

	x_prior_distribution(N_ind,q_propo,N_gen,g_possible,p,a_gene,design,
			     flank_type,left_i,right_i,R1,R2,acc_geno); 

        make_genotype(N_ind,N_gen,p,acc_geno[N_qtl]);      
       
        for (i=0;i<N_gen;i++)
	   { 	    		 
	     reg_acc.b[N_qtl][i]=limit[2][0]+
	       ((double)(rand()*(limit[2][1]-limit[2][0]))/RAND_MAX); 
	   }

       l0=log_likelihood(N_bc,N_ind,&reg_acc,take,start,
			 trait,acc_geno,bgc_acc,(N_qtl+1),
			 num_cov,cov_index,cov_type);     

        l1=log_likelihood(N_bc,N_ind,&reg_acc,take,start,
			  trait,acc_geno,bgc_acc,N_qtl,
			  num_cov,cov_index,cov_type);             

        ratio=(lambda*p_del)/((N_qtl+1)*(N_qtl+1)*p_add);

        
        chance=(((double)rand())/RAND_MAX);           	


	calc_al1++;
	/*        if (chance<=exp(l0-l1)*ratio) */             
       if (log(chance)<=(l0-l1+log(ratio)))             
          {            
	   q[N_qtl].l_x=q_propo.l_x;
           q[N_qtl].ML=q_propo.ML;
	   q[N_qtl].MR=q_propo.MR;
	   q[N_qtl].g_l=q_propo.g_l;               
	   q[N_qtl].r1=q_propo.r1;   
	   q[N_qtl].r2=q_propo.r2;          
	   acc_loca[N_qtl]=location;            
                   
           N_qtl++;
	  
	  }
	else calc_add++;

        }
       else 
	 {
	   /* delete one */
                                    
        cp_regress(N_gen,N_bc,&regress,&reg_acc,N_qtl,
		   num_cov,cov_type);    /* regress=reg_acc */

        which_one=(rand()%N_qtl);
           
	for (f=0;f<which_one;f++)
          {
           for (i=0;i<N_ind;i++)
	    genotype[f][i]=acc_geno[f][i];
    
           for (i=0;i<N_gen;i++)
	     regress.b[f][i]=reg_acc.b[f][i];
	  }

	for (f=which_one;f<(N_qtl-1);f++)
	  {
	    for (i=0;i<N_ind;i++)
	      genotype[f][i]=acc_geno[f+1][i];    

	    for (i=0;i<N_gen;i++)
	      regress.b[f][i]=reg_acc.b[f+1][i];
	  }

        
        l0=log_likelihood(N_bc,N_ind,&regress,take,start,
			  trait,genotype,bgc_acc,(N_qtl-1),
			  num_cov,cov_index,cov_type);        

        l1=log_likelihood(N_bc,N_ind,&reg_acc,take,start,
			  trait,acc_geno,bgc_acc,N_qtl,
			  num_cov,cov_index,cov_type);  

        ratio=(N_qtl*N_qtl*p_add)/(lambda*p_del);
	

        chance=(((double)rand())/RAND_MAX);           	

	calc_al2++;
	/*   if (chance<=exp(l0-l1)*ratio) */             

	if (log(chance)<=(l0-l1+log(ratio)))             
          {       
	  
	    for (f=which_one;f<(N_qtl-1);f++)            
	      {
                q[f]=q[f+1];
		acc_loca[f]=acc_loca[f+1];

                for(i=0;i<N_gen;i++)                   
		  reg_acc.b[f][i]=reg_acc.b[f+1][i];                  

		for (i=0;i<N_ind;i++)
		  acc_geno[f][i]=acc_geno[f+1][i];    
              }
           
           N_qtl--;
	
	  }
	 else calc_del++;

	 } /* delete one */
      }      
    else 
     {
        
/* ------------------------------------------------------------------- */
/*                                                                     */
/*  Here we prose new value for QTL location                           */
/*  using random walk and                                              */
/*  then check Hastings ratio (prior ratio).                           */
/*                                                                     */
/* ------------------------------------------------------------------- */ 
    
  rejected[0]=0;
  rejected[1]=0;
  rejected[2]=0;
 
  /*  num_reject[1]++; */

   how_many++;

   for (f=0;f<N_qtl;f++)
     { 
   location=acc_loca[f]-(window*0.5)+((rand()*window)/RAND_MAX);/* proposal */

      if ((location<lx_limit[0]) || (location>lx_limit[1]))
	/*||
	  ((f==0) && (acc_loca[1]-location<given_limit)) ||
          ((f==1) && (acc_loca[2]-location<given_limit)) || 
          ((f==1) && (location-acc_loca[0]<given_limit)) || 
          ((f==2) && (location-acc_loca[1]<given_limit)) ) */
         {

	   if (round!=0) 
	     { 
	       rejected[f]=1; 
	       num_reject[f][0]++;
	     } 
	   else 
	     location=acc_loca[f]; 
	 }
        
      if (!rejected[f])
	{

	  in=0;
	  while (cumul[in]<location) in++;
	  if (cumul[in]>location) in--;
          if (location==lenght) in--;       
     
	  q_propo.l_x=location;  
          q_propo.ML=in;
          q_propo.MR=in+1;
  	  q_propo.g_l=ptr_table[in];               
/* 	  q_propo.g_r=ptr_table[in+1];  */              
          dist=location-cumul[in];
	  q_propo.r1=Haldane_to_recomb(dist);

  /*  from Haldane map function X2=(X3-X1)/(1-2X1), because X3=X1+X2-2X1X2  */

	  dist=Haldane_to_recomb(q_propo.g_l->r);
          q_propo.r2=(dist-q_propo.r1)/(1.0-(2.0*q_propo.r1));
	  
          if (location==lenght) 
      	    {
               temp_d=q_propo.r1;
               q_propo.r1=q_propo.r2; 
               q_propo.r2=temp_d;      
	      }
         
/* ------------------------------------------------------------------- */
/*  Here is core for updating QTL location                             */
/* ------------------------------------------------------------------- */

    cp_field(f,N_qtl,q_new,q,q_propo);
    flank_type=flanking_objects(f,q_new,&R1,&R2,N_qtl,
				N_ind, &left_i, &right_i);

    x_prior_distribution(N_ind,q_propo,N_gen,g_possible,p,a_gene,design,
			 flank_type,left_i,right_i,R1,R2,acc_geno); 
  
 /*   x_prior_distribution(N_ind,q_propo,N_gen,g_possible,p,a_gene,design);*/
   
    if (round!=0)
      {

	flank_type=flanking_objects(f,q,&R1,&R2,N_qtl,
				    N_ind, &left_i, &right_i);

	x_prior_distribution(N_ind,q[f],N_gen,g_possible,p1,a_gene,design,
			     flank_type,left_i,right_i,R1,R2,acc_geno); 


 /*	x_prior_distribution(N_ind,q[f],N_gen,g_possible,p1,a_gene,design); */
     

	/* tarkista acc_geno:n parametrivalityksen oikeellisuus */

	make_p_ratio(N_ind,acc_geno[f],acc_geno[f],p,p1,&ratio);
      }

    chance=(((double)rand())/RAND_MAX);           
    
    if ((round==0)||(chance<=ratio))
      {
        q[f].l_x=q_propo.l_x;
        q[f].ML=q_propo.ML;
        q[f].MR=q_propo.MR;
        q[f].g_l=q_propo.g_l;               
	q[f].r1=q_propo.r1;   
        q[f].r2=q_propo.r2;          
     	acc_loca[f]=location;        
	/*  make_take(take,bgc,q.ML,q.MR,N_bc,cur_chr);  */ 
      }
    else num_reject[f][0]++;


	}
 
     
    
/* ------------------------------------------------------------------- */
/*                                                                     */
/*  Here we update QTL genotypes individually                          */
/*                                                                     */
/* ------------------------------------------------------------------- */

 /* x_prior_distribution(N_ind,q,N_gen,g_possible,p,a_gene,design); */


      flank_type=flanking_objects(f,q,&R1,&R2,N_qtl,
				  N_ind, &left_i, &right_i);

      x_prior_distribution(N_ind,q[f],N_gen,g_possible,p,a_gene,design,
			   flank_type,left_i,right_i,R1,R2,acc_geno); 
  

    make_genotype(N_ind,N_gen,p,genotype[f]);      

     
    if ((round==0) && (f==0)) 
       {
        make_genotype(N_ind,N_gen,p,acc_geno[0]);      
        make_genotype(N_ind,N_gen,p,acc_geno[1]);      
        make_genotype(N_ind,N_gen,p,acc_geno[2]);      
       }

    start_1=start;

    for (i=0;i<N_ind;i++)
       {

         if (f!=0) genotype[0][i]=acc_geno[0][i]; 
         if (f!=1) genotype[1][i]=acc_geno[1][i]; 
         if (f!=2) genotype[2][i]=acc_geno[2][i]; 

	 l01=ind_log_likelihood (N_bc,&reg_acc,take,start_1->phenotype,
				 genotype,bgc_acc,i,N_qtl,trait,
				 num_cov,cov_index,cov_type);
    
	 l11=ind_log_likelihood (N_bc,&reg_acc,take,start_1->phenotype,
				 acc_geno,bgc_acc,i,N_qtl,trait,
				 num_cov,cov_index,cov_type);
     
        chance=(((double)rand())/RAND_MAX); 

        if ((round==0) || (chance<=exp(l01-l11)))             
           acc_geno[f][i]=genotype[f][i];                  
        else mup3[f]++;

        start_1=start_1->next;
       }
	            
     }
     
}
   
/* ------------------------------------------------------------------- */
/*                                                                     */
/*  Here we update marker genotype matrix                              */
/*  in current chromosome under the search.                            */
/*                                                                     */
/* ------------------------------------------------------------------- */     

   get_marker_proposal (N_ind,Num_marks,N_gen,gene,g_possible,start,cur_chr);

     if (design==1)
        BCross_mrk_prior(N_ind, Num_marks,wg,mc_prior,gene); 
     else
        F2_mrk_prior(N_ind, Num_marks,wg,mc_prior,gene); 

     
   for (f=0;f<N_qtl;f++)   
     {

       flank_type=flanking_objects(f,q,&R1,&R2,N_qtl,
				   N_ind, &left_i, &right_i);

       x_prior_distribution(N_ind,q[f],N_gen,g_possible,pp[f],gene,design,
			    flank_type,left_i,right_i,R1,R2,acc_geno); 
  
       x_prior_distribution(N_ind,q[f],N_gen,g_possible,pp1[f],a_gene,design,
			    flank_type,left_i,right_i,R1,R2,acc_geno); 
  

 /* x_prior_distribution(N_ind,q[f],N_gen,g_possible,pp[f],gene,design);   */  
 /* x_prior_distribution(N_ind,q[f],N_gen,g_possible,pp1[f],a_gene,design);*/
 
     }     
   
       
   for (i=0;i<N_ind;i++) 
     {         
       ratio=(mc_prior[i])/(acc_prior[i]);
       
       for (f=0;f<N_qtl;f++)   
	 ratio=ratio*(pp[f][i][acc_geno[f][i]]/pp1[f][i][acc_geno[f][i]]);  

       chance=(((double)rand())/RAND_MAX);
          
       if (chance<=ratio)  
	 {   
	   for(j=0;j<Num_marks;j++)
	     a_gene[i][j]=gene[i][j];       
	   
	   acc_prior[i]=mc_prior[i];
	 }
       else prop_mup++;      
     }     


/* ------------------------------------------------------------------- */
/*                                                                     */
/*  Here we update regression parameters and check acceptance          */
/*                                                                     */
/*                                                                     */
/*  prior is omitted from acceptance ratio, because it is uniform      */
/*  distribution and therefore it cancels out from the equation.       */
/*                                                                     */
/* ------------------------------------------------------------------- */

   for (f=0;f<2;f++)
     {

     cp_regress(N_gen,N_bc,&regress,&reg_acc,N_qtl, num_cov,
		cov_type);  /* regress=reg_acc */

     get_regression_proposals(f,N_gen,N_bc,range,&regress,num_cov,
			      cov_type);
     fail=0;

     if (f==0)
       {
     if ((regress.a < limit[0][0]) || (regress.a > limit[0][1])) 
       fail=1;
     if ((regress.sigma <limit[1][0]) || (regress.sigma > limit[1][1])) 
       fail=1; 
    
     for (i=0;i<num_cov;i++)
       {
	 if ((regress.cov[i][0] < limit[0][0]) || 
	    (regress.cov[i][0] > limit[0][1])) 
	      fail=1;
       
	 for (m=1;m<cov_type[i];m++)
	   if ((regress.cov[i][m] < limit[0][0]) ||
               (regress.cov[i][m] > limit[0][1]))
                  fail=1;

       }
       }

    if (f==1)
      {
     for (i=0;i<N_gen;i++)
       {     
         for(j=0;j<N_qtl;j++)
      if ((regress.b[j][i] < limit[2][0]) || (regress.b[j][i] > limit[2][1])) 
	fail=1;

	 for (j=0;j<N_bc;j++)
	   if ((regress.bc[j][i] < limit[3][0]) || 
	       (regress.bc[j][i] > limit[3][1])) 
	       fail=1; 
       }	

      }
     
     if (!fail)
       {
      l0=log_likelihood(N_bc,N_ind,&regress,take,start,trait,acc_geno,
			bgc_acc,N_qtl,num_cov,cov_index,cov_type);
      l1=log_likelihood(N_bc,N_ind,&reg_acc,take,start,trait,acc_geno,
			bgc_acc,N_qtl,num_cov,cov_index,cov_type);

      if ((((double)rand())/RAND_MAX)<=exp(l0-l1))/* accept with prob (l0/l1)*/
	{
         cp_regress(N_gen,N_bc,&reg_acc,&regress,N_qtl,num_cov,cov_type); 
	                                  /* reg_acc=regress */
	 /* l1=l0; */
	}
       else num_rej[f]++;
     }
     else num_rej[f]++;

   }

/* ------------------------------------------------------------------- */
/*                                                                     */
/*  Here we update BG control parameters locating in other chromosomes */
/*  and check an acceptance                                            */
/*                                                                     */
/* ------------------------------------------------------------------- */
       
      bg_proposal (N_ind,N_bc,N_gen,bgc,bgc_table,bgc_new,g_possible);
      change_cur(N_bc,N_ind,cur_chr,bgc,bgc_new,bgc_acc);

/* l0=log_likelihood(N_bc,N_ind,&reg_acc,take,start,trait,acc_geno,bgc_new); */
      
      start_1=start;

      for (i=0;i<N_ind;i++)
        {
      l01=ind_log_likelihood (N_bc,&reg_acc,take,start_1->phenotype,
			      acc_geno,bgc_new,i,N_qtl,trait,num_cov,cov_index,
			      cov_type);

      l11=ind_log_likelihood (N_bc,&reg_acc,take,start_1->phenotype,
			      acc_geno,bgc_acc,i,N_qtl,trait,num_cov,cov_index,
			      cov_type);

         chance=(((double)rand())/RAND_MAX); 
	  
         if (chance<=exp(l01-l11))          
	   {
	     for(j=0;j<N_bc;j++)
	       bgc_acc[i][j]=bgc_new[i][j];
	   }
         else mup2++;
                      
         start_1=start_1->next;
        }
       
       
/* accept with prob (l0/l1) */

	/*      if ((((double)rand())/RAND_MAX)<=exp(l0-l1)) 
	{
          update_bgc(N_ind,N_bc,bgc_acc,bgc_new);
	  l1=l0;
        }
        else mup2++; */

 
  if (N_qtl>0) 
    fprintf(fp,"%12.7f\n",acc_loca[0]);
  if (N_qtl>1)
    fprintf(fp0,"%12.7f\n",acc_loca[1]);
  if (N_qtl>2) 
    fprintf(fp1,"%12.7f\n",acc_loca[2]);
  
   fprintf(fp5,"%8d\n",N_qtl);

if (permission)
  {
      true_inter=reg_acc.a;
      for (i=0;i<N_qtl;i++)
	true_inter+=reg_acc.b[i][0];
      for (i=0;i<N_bc;i++)
	if (take[i]) true_inter+=reg_acc.bc[i][0];

      fprintf(fr0,"%12.7f ",reg_acc.a);
      fprintf(fr0,"%12.7f ",true_inter);
      fprintf(fr0,"%12.7f \n",reg_acc.sigma);


      for (i=0;i<num_cov;i++)
	{
	  fprintf(fr20,"%12.7f ",reg_acc.cov[i][0]);
	  for (j=1;j<cov_type[i];j++)
	    fprintf(fr20,"%12.7f ",reg_acc.cov[i][j]);		    
	  fprintf(fr20,"\n");
	}

    
      if (N_qtl>0)
	{
	  for (i=0;i<N_gen;i++)
	    fprintf(fr1,"%12.7f ",reg_acc.b[0][i]);
	  fprintf(fr1,"\n");
	}

      if (N_qtl>1)
	{
	  for (i=0;i<N_gen;i++)
	    fprintf(fr2,"%12.7f ",reg_acc.b[1][i]);
	  fprintf(fr2,"\n");
	}

      if (N_qtl>2)
	{
	  for (i=0;i<N_gen;i++)
	    fprintf(fr3,"%12.7f ",reg_acc.b[2][i]);
	  fprintf(fr3,"\n");
	}

    }

  /*     fprintf(fp,"%8d %12.7f\n",round,acc_loca); */       
  /*        fprintf(fp2,"%8d ",round);
	 for (i=0;i<N_ind;i++)
	    fprintf(fp2,"%8d ",acc_geno[i]);
         fprintf(fp2,"\n");   */


/*         fprintf(fp3,"%8d %12.7f %12.7f ",round,reg_acc.a,reg_acc.sigma);
         for (i=0;i<N_gen;i++)
	   {
	   fprintf(fp3,"%12.7f ",reg_acc.b[i]);
           for (j=0;j<N_bc;j++) 
	     fprintf(fp3,"%12.7f ",reg_acc.bc[j][i]);
  	 }
         fprintf(fp3,"\n");  
	 */

/*       fprintf(fp4,"%8d ",round);
       for (i=0;i<N_ind;i++)
           for (j=0;j<Num_marks;j++)
	     fprintf(fp4,"%8d ",a_gene[i][j]);
       fprintf(fp4,"\n");  */

     round++;
      }

 for (f=0;f<3;f++)
  { 
 fprintf(fp2,"Rejection rate %d for QTL location: %f \n",
	 f,(double) num_reject[f][0]/how_many);

 fprintf(fp2,"Rejection rate %d for QTL genotypes: %f\n",
	 f,(double) mup3[f]/(how_many*N_ind));

   }


 fprintf(fp2,"Rejection rate for marker update in current chr.: %f\n",
	 (double) prop_mup/(max_round*N_ind));

 fprintf(fp2,"Rejection rate for regression parameters : \n");

 fprintf(fp2,"  block 0 = (a, env. covariates, sigma)     : %f\n",    
           	  (double) num_rej[0]/max_round);

 fprintf(fp2,"  block 1 = (all QTL and BGC coefficients   : %f\n",
	          (double) num_rej[1]/max_round);

/* fprintf(fp2,"Rejection rate for BG marker update in current chr.: %f\n",
   (double) mup/(max_round*N_ind)); */

 fprintf(fp2,"Rejection rate for BG update in other chrom.: %f\n",
	 (double) mup2/(max_round*N_ind));

 fprintf(fp2,"Rejection rate for add step: %f\n",
	 (double) calc_add/calc_al1);

 fprintf(fp2,"Rejection rate for del step: %f\n",
	 (double) calc_del/calc_al2);

fclose(fp);
fclose(fp0);
fclose(fp1);
fclose(fp2); 
fclose(fp3); 
/* fclose(fp4); */
fclose(fp5); 

fclose(fr0); 
fclose(fr1); 
fclose(fr2); 
fclose(fr3); 
fclose(fr20);

}