/* This file contains the low-level probability routines in */
/* a modified version of the ILINK program*/

#include "commondefs.h"
#if defined(GENERR)
#include "err.h"
#endif

/* Local variables for getvect: */
struct LOC_getvect {
  struct LOC_likelihood *LINK;
  thisperson *p;
  hapvector hap1, hap2;
} ;

/*Local Void getgene PP((int syste, double val, struct LOC_getvect *LINK));
Local Void ugetgene PP((int syste, double val, struct LOC_getvect *LINK));*/

/*static void getgene (int syste, double val, struct LOC_getvect *LINK);
static void ugetgene (int syste, double val, struct LOC_getvect *LINK);*/

void getgene ();
void ugetgene ();

Local double quanfun(phen, thislocus, i, j, mean, LINK)
phenotype *phen;
locusvalues *thislocus;
int i, j;
double *mean;
struct LOC_getvect *LINK;
{
  double val;
  int it, jt, FORLIM, FORLIM1;

  val = 1.0;
  if (phen->missing)
    return val;
  val = 0.0;
  FORLIM = thislocus->UU.U1.ntrait;
  for (it = 0; it < FORLIM; it++) {
    FORLIM1 = thislocus->UU.U1.ntrait;
    for (jt = 0; jt < FORLIM1; jt++) {
      if (i == j)
	val += thislocus->UU.U1.vmat[it]
	       [jt] * (phen->x[jt] - mean[jt]) *
	       (phen->x[it] - mean[it]);
      else
	val += thislocus->UU.U1.conmat * thislocus->UU.U1.vmat[it]
	       [jt] * (phen->x[jt] - mean[jt]) *
	       (phen->x[it] - mean[it]);
    }
  }
  val = thislocus->UU.U1.det * exp(-val * 0.5);
  if (i != j)
    val *= thislocus->UU.U1.contrait;
  return val;
}  /*quanfun*/


#if defined(GENERR)
/* Meg: Begin */

void incompene(numa,err,longpene)
     int numa;
     double err, longpene[maxliab*maxliab];
{
  int numcom, i, j, count;
  double pene[maxliab][maxliab];

  numcom = numa*(numa+1)/2;

  for (i=0; i<=(numcom-1); i++)
    for (j=0; j<= (numcom-1); j++)
      pene[i][j] = 0.0;
  
  for (i=0; i<(maxliab*maxliab); i++)
    longpene[i] = 0.0;

  for (i=0; i<=(numcom-1); i++)
    {
      for (j=0; j<=(numcom-1); j++)
	{
	  pene[i][j]= err/((double) numcom);
	  if (i==j) 
	    pene[i][j] += 1-err;
	}
    }

  count = 0;

  for (i=0; i<numcom; i++)
    for (j=0; j<numcom; j++)
      {
	longpene[count] = pene[i][j];
	count += 1;
      }
}

void putlocuspen(err)
     double err;
{
  double longpene[maxliab*maxliab];
  int i,j,k,l,n,count,numcom, numa;
  void incompene();
 
  for (n=1; n<=mlocus; n++)
    {
      numa = thislocus[n-1]->nallele;
      numcom = numa*(numa+1)/2;

      incompene(numa, err, longpene);
           
      count = 0;
      for (l=0; l<numcom; l++) {
	for (i=1; i<=numa; i++) {
	  for (j=(i-1); j<numa; j++) {
	    altlocuspen[n-1][i][j][2][l] = longpene[count];
	    count += 1;
	    altlocuspen[n-1][i][j][1][l] = 1-altlocuspen[n-1][i][j][2][l];
	    altlocuspen[n-1][i][j][0][l] = 1.0;
	    for (k=0; k<=2; k++)
	      altlocuspen[n-1][j+1][i-1][k][l] = altlocuspen[n-1][i][j][k][l];
	  }
	}

	for (i=0; i<numa; i++)
	  altlocuspen[n-1][0][i][0][l] = 1.0;
      }
    }
}
/* Meg: End */
#endif

Local Void getval(syste, i, j, val, LINK)
int syste, i, j;
double *val;
struct LOC_getvect *LINK;
{
  locusvalues *WITH;
  phenotype *WITH1;

  WITH = thislocus[syste - 1];
  switch (WITH->which) {

  case quantitative:
    *val *= quanfun(LINK->p->phen[syste - 1], thislocus[syste - 1], i, j,
		    WITH->UU.U1.pm[i][j - 1], LINK);
    break;

  case affection:
    WITH1 = LINK->p->phen[syste - 1];
    
#if defined(GENERR)
    /* Meg: Begin */

    if ((LINK->p->phen[syste-1]->liability==1)&&(LINK->p->phen[syste-1]->aff==0)&&(firsttime==true))
      personunk[syste-1][LINK->p->id-1]=0;

    /*Meg: End */
#endif

#if defined(GENERR)
    /* Meg: Begin */
    if ((NewPene==true)&&((LINK->p->id-1)==iplace)&&(jplace==(syste-1))) 
      *val *= altlocuspen[syste - 1][i][j - 1][WITH1->aff][WITH1->liability - 1];
    else
      *val *= WITH->UU.U0.pen[i][j - 1][WITH1->aff]
	  [WITH1->liability - 1];
    /* Meg: End */
#else
    *val *= WITH->UU.U0.pen[i][j - 1][WITH1->aff][WITH1->liability - 1];
#endif
    

    break;
  }
}  /*getval*/

#include "comgetvect.c"