/* This file contains the main routines in */
/* a modified version of the ILINK program*/
/* The modifications are described in the papers: */
/* R. W. Cottingham, Jr., R. M. Idury, and A. A. Schaffer, */
/* Faster Sequential Genetic Linkage Computations */
/* American Journal of Human Genetics, 53(1993), pp. 252-263*/
/* and A. A. Schaffer, S. K. Gupta, K. Shriram, and R. W. Cottingham, Jr.,*/
/* Avoiding Recomputation in Linkage Analysis */
/* Human Heredity 44(1994), pp. 225-237.*/

#include "commondefs.h"
#if defined(GENERR)
#include "err.h"
#endif
#if !defined(DOS)
#include "checkpointdefs.h"
#endif
#include "gemdefs.h"
#include "ildefs.h"
#ifndef LESSMEMORY
#include "moddefs.h"
#endif

#if defined(GENERR)
/* Meg: Begin */
void indexx(n,arrin,indx)
     int n,indx[];
     double arrin[];
{
  int l,j,ir,indxt,i;
  double q;

  for (j=1;j<=n;j++) indx[j]=j;
  l=(n >> 1) + 1;
  ir=n;
  for (;;) {
    if (l > 1)
      q=arrin[(indxt=indx[--l])];
    else {
      q=arrin[(indxt=indx[ir])];
      indx[ir]=indx[1];
      if (--ir == 1) {
	indx[1]=indxt;
	return;
      }
    }
    i=l;
    j=l << 1;
    while (j <= ir) {
      if (j < ir && arrin[indx[j]] < arrin[indx[j+1]]) j++;
      if (q < arrin[indx[j]]) {
	indx[i]=indx[j];
	j += (i=j);
      }
      else j=ir+1;
    }
    indx[i]=indxt;
  }
}
/* Meg: End */
#endif

static Void outcontrol(z)
double *z;
{
  int i, j, FORLIM;

  j = 0;
  FORLIM = n;
  for (i = 0; i < FORLIM; i++) {
    fprintf(outfile, "% .5le", z[i]);
    j++;
    if (j == 4) {
      putc('\n', outfile);
      j = 0;
    }
  }
  if (j != 0)
    putc('\n', outfile);
}


double mapfunction(theta1, theta2)
double theta1, theta2;
{
  /*User defined function giving recombination between
    flanking markers as a function of recombination
    between adjacent markers*/
  return ((theta1 + theta2) / (1 + 4 * theta1 * theta2));
}


double getdist(theta)
double *theta;
{
  if (*theta < 0.5)
    return (log(1.0 - 2.0 * *theta) / -2.0);
  else
    return 10.0;
}


double invdist(dist)
double *dist;
{
  if (*dist != 10.0)
    return ((1.0 - exp(-2 * *dist)) / 2.0);
  else
    return 0.5;
}

/*saveparams saves the parameter values at the beginning of
gforward or gcentral, so we capture the optimal thetas*/
Local Void saveparams()
{
  int i, FORLIM;

  FORLIM = nall;
  for (i=0; i < FORLIM; i++){
       outsavex[i] = x[i];
  }
}

static Void fun(f, x)
double *f;
double *x;
{
  int i, k, FORLIM;

  k = 0;
  FORLIM = nall;
  /*Test code inserted by Alex
  for (i = 0; i < FORLIM; i++) {
    printf("%lf ", x[i]);
  } 
  printf("\n");*/
  for (i = 0; i < FORLIM; i++) {
    if (itp[i] == 1) {
      k++;
      xall[i] = x[k - 1];
    }
  }
  setparam();
  recombination();
  if (penalty) {
    *f = 1.5 * penlike;
    return;
  }
  FORLIM = totperson;
  for (i = 1; i <= FORLIM; i++)
    person[i]->done = false;
  FORLIM = totperson;
  for (i = 1; i <= FORLIM; i++)
    person[i]->firstpass = true;
  alike = 0.0;
  thisc = minint;
  FORLIM = nuped;
  for (i = 1; i <= FORLIM; i++) {
    likelihood(i, proband[i - 1]);
    likebyped[i-1] = like; /*Store likelihood for later*/
    alike += like;
  }
  if (firsttime) {
    if (thisc < maxcensor)
      printf("Maxcensor can be reduced to %12ld\n", thisc);
    else {
      if (thisc > maxcensor)
	printf("You may gain efficiency by increasing maxcensor\n");
    }
  }
  firsttime = false;
  *f = -2 * alike;
  penlike = *f;
  firsttime = false;
/*  exit(EXIT_SUCCESS); /*exitramana*/
}


/* Local variables for outf: */
struct LOC_outf {
  double lods;
  double thisval[maxped];
} ;

Local Void getlods(f, x, LINK)
double *f;
double *x;
struct LOC_outf *LINK;
{
  int i, k, FORLIM;

  firstapprox = true;
  k = 0;
  FORLIM = nall;
  for (i = 0; i < FORLIM; i++) {
    if (itp[i] == 1) {
      k++;
      xall[i] = x[k - 1];
    }
  }
  FORLIM = nuped;
  for (i = 0; i < FORLIM; i++)
    LINK->thisval[i] = 0.0;
  setparam();
  if (!penalty) {
    if (byfamily) {
      for (i = 1; i <= 35; i++)
	putc('-', final);
      fprintf(final, "\nPEDIGREE | LN LIKE | LOG 10 LIKE\n");
      for (i = 1; i <= 35; i++)
	putc('-', final);
      putc('\n', final);
    }
    if (byfamily)   {  /*Then need to redo likelihoods on a by-family basis*/
      FORLIM = totperson;
      for (i = 1; i <= FORLIM; i++)
        person[i]->done = false;
      FORLIM = totperson;
      for (i = 1; i <= FORLIM; i++)
        person[i]->firstpass = true;
      recombination();
      thisc = minint;
      FORLIM = nuped;
      for (i = 1; i <= FORLIM; i++) {
	likelihood(i, proband[i - 1]);
	LINK->thisval[i - 1] = like;
      }
    }
  }				/* Shriram: reindented */
  firsttime = true;
  FORLIM = mlocus;
  for (i = 0; i < FORLIM; i++)
    maletheta->theta[i] = 0.5;
  FORLIM = mlocus;
  for (i = 0; i < FORLIM; i++)
    femaletheta->theta[i] = 0.5;
  if (penalty) {
    *f = 1.5 * penlike;
    return;
  }
  dolod = true;
  FORLIM = totperson;
  for (i = 1; i <= FORLIM; i++)
    person[i]->done = false;
  FORLIM = totperson;
  for (i = 1; i <= FORLIM; i++)
    person[i]->firstpass = true;
  recombination();
  alike = 0.0;
  thisc = minint;
  FORLIM = nuped;
  for (i = 0; i < FORLIM; i++) {
    likelihood(i + 1, proband[i]);
    alike += like;
    if (byfamily) {
      fprintf(final, "%9ld %9ld %12.6f",
	      i + 1, proband[i]->ped, LINK->thisval[i] - like);
      fprintf(final, "%12.6f", (LINK->thisval[i] - like) / log10_);
      if (like == zerolike)
	fprintf(final, " inconsistent data\n");
      else
	putc('\n', final);
      if (byfamily)  /*Need to compute family lodscore*/
        if (like != zerolike)
	  LINK->thisval[i] -= like;
    }
  }
  *f = -2 * alike;
  penlike = *f;
  dolod = false;
}

Local Void streamout(LINK)
struct LOC_outf *LINK;
{
  int i, j, k, l;
  double deriv[3];
  thetarray oldtheta;
  double dist;
  int FORLIM;
  locusvalues *WITH;
  int FORLIM1, FORLIM2;
  thetavalues *WITH1;
  double TEMP;

  inconsistent = false;
  setparam();
  recombination();
  firstapprox = true;
  /*Start stream*/
  fprintf(stream, "@\n");
  fprintf(stream, "ILINK\n");

  /*Gemini information*/
  fprintf(stream, "% .5le %12ld % .5le %12ld %12ld %12ld\n",
	  f, nit, ptg, idg, nall, n);
  /*Valid or not*/
  if (nit > 1 && fabs(ptg) < 0.001 || n == 0)
    fprintf(stream, " 2\n");
  else {
    if (nit > 1 && fabs(ptg) < 0.1)
      fprintf(stream, " 1\n");
    else
      fprintf(stream, " 0\n");
  }
  FORLIM = nall;
  for (i = 0; i < FORLIM; i++)
    fprintf(stream, " % .5le\n", xall[i]);
  putc('\n', stream);
  FORLIM = nall;
  for (i = 0; i < FORLIM; i++)
    fprintf(stream, " %2ld\n", itp[i]);
  putc('\n', stream);
  FORLIM = n;
  for (i = 0; i < FORLIM; i++)
    fprintf(stream, " % .5le\n", g[i]);
  putc('\n', stream);
  /*Variance-covariance calculated if ivar=1 and icall=0*/
  if (ivar == 1 && icall == 0)
    fprintf(stream, "1\n");
  else
    fprintf(stream, "0\n");

  if (inconsistent)
    fprintf(stream, "1\n");
  else
    fprintf(stream, "0\n");

  /*Control parameters for sexlink,interference and sex difference*/
  if (sexlink)
    fprintf(stream, "1\n");
  else
    fprintf(stream, "0\n");
  if (interfer)
    fprintf(stream, "1 ");
  else
    fprintf(stream, "0 ");
  if (mapping)
    fprintf(stream, "1\n");
  else
    fprintf(stream, "0\n");
  if (sexdif && readfemale)
    fprintf(stream, "2 ");
  else {
    if (sexdif)
      fprintf(stream, "1\n");
    else
      fprintf(stream, "0\n");
  }

  /*Genetic information*/
  fprintf(stream, "%12ld\n", mlocus);
  FORLIM = mlocus;
  for (i = 1; i <= FORLIM; i++) {
    j = 0;
    do {
      j++;
    } while (order[j - 1] != i);
    fprintf(stream, " %12ld", j);
  }
  putc('\n', stream);

  /*Iterated locus*/
  if (itsys != 0)
    fprintf(stream, "%12ld\n", order[itsys - 1]);
  else
    fprintf(stream, " 0\n");
  if (itsys != 0) {
    WITH = thislocus[itsys - 1];
    if (disequi) {
      fprintf(stream, "1 %12ld %12d\n", WITH->nallele, nuhap);
      FORLIM = nuhap;
      for (i = 0; i < FORLIM; i++)
	fprintf(stream, "%9.6f", hapfreq->genarray[i]);
    } else {
      fprintf(stream, "0 %12ld\n", WITH->nallele);
      FORLIM = WITH->nallele;
      for (i = 0; i < FORLIM; i++)
	fprintf(stream, "%9.6f", WITH->freq[i]);
    }
    putc('\n', stream);

    WITH = thislocus[itsys - 1];
    if (WITH->which == binary_ && WITH->format == 3)
      fprintf(stream, " 3\n");
    else {
      if (WITH->which == binary_ && WITH->format == 2)
	fprintf(stream, " 2\n");
      else {
	if (WITH->which == quantitative) {
	  fprintf(stream, "0 %12ld\n", WITH->UU.U1.ntrait);
	  invert(WITH->UU.U1.vmat, WITH->UU.U1.ntrait, &WITH->UU.U1.det);
	  FORLIM = WITH->UU.U1.ntrait;
	  for (i = 0; i < FORLIM; i++) {
	    FORLIM1 = WITH->nallele;
	    for (j = 1; j <= FORLIM1; j++) {
	      FORLIM2 = WITH->nallele;
	      for (k = j - 1; k < FORLIM2; k++)
		fprintf(stream, "%6.3f ", WITH->UU.U1.pm[j][k][i]);
	    }
	    putc('\n', stream);
	  }
	  FORLIM = WITH->UU.U1.ntrait;
	  for (i = 0; i < FORLIM; i++) {
	    FORLIM1 = WITH->UU.U1.ntrait;
	    for (j = i; j < FORLIM1; j++)
	      fprintf(stream, "%6.3f ", WITH->UU.U1.vmat[i][j]);
	    putc('\n', stream);
	  }
	} else {
	  if (WITH->which == affection) {
	    fprintf(stream, "1 %12ld\n", WITH->UU.U0.nclass);
	    FORLIM = WITH->UU.U0.nclass;
	    for (l = 0; l < FORLIM; l++) {
	      FORLIM1 = WITH->nallele;
	      for (i = 1; i <= FORLIM1; i++) {
		FORLIM2 = WITH->nallele;
		for (j = i - 1; j < FORLIM2; j++)
		  fprintf(stream, " %5.3f", WITH->UU.U0.pen[i][j][2][l]);
	      }
	      putc('\n', stream);
	      if (sexlink) {
		FORLIM1 = WITH->nallele;
		for (i = 0; i < FORLIM1; i++)
		  fprintf(stream, "%5.3f ", WITH->UU.U0.pen[0][i][2][l]);
	      }
	      putc('\n', stream);
	    }
	  }
	}
      }
    }
  }
  if (interfer && !mapping) {
    WITH1 = maletheta;
    memcpy(oldtheta, WITH1->theta, sizeof(thetarray));
    WITH1->theta[0] =
      (oldtheta[0] + oldtheta[mlocus - 1] - oldtheta[mlocus - 2]) / 2.0;
    WITH1->theta[mlocus - 2] =
      (oldtheta[mlocus - 2] + oldtheta[mlocus - 1] - oldtheta[0]) / 2.0;
    WITH1->theta[mlocus - 1] =
      (oldtheta[0] + oldtheta[mlocus - 2] - oldtheta[mlocus - 1]) / 2.0;
    FORLIM = mlocus;
    for (i = 0; i < FORLIM; i++)
      fprintf(stream, " %5.3f", WITH1->theta[i]);
    putc('\n', stream);
    memcpy(WITH1->theta, oldtheta, sizeof(thetarray));
    if (sexdif) {
      WITH1 = femaletheta;
      memcpy(oldtheta, WITH1->theta, sizeof(thetarray));
      WITH1->theta[0] =
	(oldtheta[0] + oldtheta[mlocus - 1] - oldtheta[mlocus - 2]) / 2.0;
      WITH1->theta[mlocus - 2] =
	(oldtheta[mlocus - 2] + oldtheta[mlocus - 1] - oldtheta[0]) / 2.0;
      WITH1->theta[mlocus - 1] =
	(oldtheta[0] + oldtheta[mlocus - 2] - oldtheta[mlocus - 1]) / 2.0;
      FORLIM = mlocus;
      for (i = 0; i < FORLIM; i++)
	fprintf(stream, " %5.3f", WITH1->theta[i]);
      putc('\n', stream);
      memcpy(WITH1->theta, oldtheta, sizeof(thetarray));
    }
    if (ivar == 1 && icall == 0) {
      for (i = 0; i <= 2; i++) {
	TEMP = 1 + exp(xall[nall - i - 1]);
	deriv[i] = exp(xall[nall - i - 1]) / (TEMP * TEMP);
      }
      k = -1;
      for (i = 0; i <= 2; i++) {
	if (itp[nall - i - 1] == 1) {
	  k++;
	  l = -1;
	  for (j = 0; j <= 2; j++) {
	    if (itp[nall - j - 1] == 1) {
	      l++;
	      bmat[n - k - 1][n - l - 1] *= deriv[j] * deriv[i];
	    }
	  }
	  if (l == -1)
	    l = 0;
	  for (j = n - l - 1; j >= 0; j--) {
	    bmat[j][n - k - 1] *= deriv[i];
	    bmat[n - k - 1][j] = bmat[j][n - k - 1];
	  }
	}
      }
    }
  }

  WITH1 = maletheta;
  FORLIM = mlocus - 2;
  /*Recombination*/
  for (i = 0; i <= FORLIM; i++)
    fprintf(stream, " %5.3f", WITH1->theta[i]);
  if (interfer)
    fprintf(stream, " %5.3f\n", maletheta->theta[mlocus - 1]);
  else
    putc('\n', stream);
  if (sexdif) {
    WITH1 = femaletheta;
    FORLIM = mlocus - 2;
    for (i = 0; i <= FORLIM; i++)
      fprintf(stream, " %5.3f", WITH1->theta[i]);
    if (interfer)
      fprintf(stream, " %5.3f\n", femaletheta->theta[mlocus - 1]);
    else
      putc('\n', stream);
    if (readfemale) {
      FORLIM = mlocus - 2;
      for (i = 0; i <= FORLIM; i++) {
	dist = getdist(&maletheta->theta[i]);
	if (dist != 0.0)
	  dist = getdist(&femaletheta->theta[i]) / dist;
	else
	  dist = 0.0;
	fprintf(stream, " %5.3f\n", dist);
      }
      if (interfer) {
	dist = getdist(&maletheta->theta[mlocus - 1]);
	if (dist != 0.0)
	  dist = getdist(&femaletheta->theta[mlocus - 1]) / dist;
	else
	  dist = 0.0;
	fprintf(stream, " %5.3f\n", dist);
      }
    } else
      fprintf(stream, "%5.3f\n", distratio);
  }

  /*Valid parameters or not*/
  /*Lod scores are valid only if penalty is false*/
  if (penalty)
    fprintf(stream, "1\n");
  else
    fprintf(stream, "0\n");
  /*Number of pedigrees*/
  fprintf(stream, "%12ld\n", nuped);
  FORLIM = nuped;
  for (i = 0; i < FORLIM; i++)
    fprintf(stream, "%12ld % .5le % .5le\n",
	    proband[i]->ped, outsavelike[i], outsavelike[i] / log10_);
  /*Total lod-score*/
  fprintf(stream, "% .5le\n", (LINK->lods - f) / (2 * log10_));
  if (ivar == 1 && icall == 0) {
    FORLIM = n;
    for (i = 0; i < FORLIM; i++) {
      FORLIM1 = n;
      for (j = i; j < FORLIM1; j++)
	fprintf(stream, " % .5le\n", 2.0 * bmat[i][j]);
      putc('\n', stream);
    }
  }
  /*Stop stream*/
  fprintf(stream, "~\n");

  /*sexdif=true and readfemale=false*/
}

static Void outf()
{
  struct LOC_outf V;
  int i, j, k, l;
  double deriv[3];
  thetarray oldtheta;
  double dist;
  int FORLIM;
  locusvalues *WITH;
  int FORLIM1, FORLIM2;
  thetavalues *WITH1;
  double TEMP; 

#if !defined(DOS)
  /* Shriram: begin */
  if ( normalRun == checkpointStatus )
  {
    funCallPath = fCP_outf ;
    performCheckpoint ( functionLocation , funCallPath , 0 ) ;

    fclose ( final ) ;
    copyFile ( "final.dat" , OutfILINKFinalDat ) ;
#ifdef vms
    final = fopen ( "final.dat" , "a", "ctx=rec","shr=get,put,upd" ) ;
#else
    final = fopen ( "final.dat" , "a" ) ;
#endif

    if ( dostream )
    {
      fclose ( stream ) ;
      copyFile ( "stream.dat" , OutfILINKStreamDat ) ;
#ifdef vms
      stream = fopen ( "stream.dat" , "a", "ctx=rec","shr=get,put,upd" ) ;
#else
      stream = fopen ( "stream.dat" , "a" ) ;
#endif
    }
  }
  else
  {
    if ( functionLocation == ckptTuple . ckptLocation )
    {
      recoverCheckpoint ( NULL ) ;

      fclose ( final ) ;
      copyFile ( OutfILINKFinalDat , "final.dat" ) ;
#ifdef vms
      final = fopen ( "final.dat" , "a", "ctx=rec","shr=get,put,upd" ) ;
#else
      final = fopen ( "final.dat" , "a" ) ;
#endif
      
      if ( dostream )
      {
        fclose ( stream ) ;
	copyFile ( OutfILINKStreamDat , "stream.dat" ) ;
#ifdef vms
        stream = fopen ( "stream.dat" , "a", "ctx=rec","shr=get,put,upd" ) ;
#else
        stream = fopen ( "stream.dat" , "a" ) ;
#endif
      }
    }

  }
  /* Shriram: end */
#endif

  firstapprox = true;
  lasttime = true;
  /*recover optimal x and f values, put x value into xall*/
  f = outsavefvalue;
  FORLIM = nall;
  k = 0;
  for (i = 0; i < FORLIM; i++) {
    if (itp[i] == 1) {
      k++;
      xall[i] = outsavex[k - 1];
      x[i] = outsavex[k - 1];  /* Meg */
    }
  }

  /*recover xall values into maletheta and femaletheta*/
  k = 0;
  FORLIM = mlocus - 1;
  for (i=0; i < FORLIM; i++){
    k++;
    maletheta->theta[i] = xall[k-1];
  }
  if (interfer && !(mapping)) {
    k++;
    maletheta->theta[mlocus - 1] = xall[k -1];
  }
  if (sexdif) {
    if (readfemale) {
    FORLIM = mlocus -1;
    for (i=0; i < FORLIM; i++) {
      k++;
      femaletheta->theta[i] = xall[k-1];
    }
      if (interfer && !(mapping)) {
        k++;
        femaletheta->theta[mlocus - 1] = xall[k -1];
      }
    }
    else {
      k++;
      distratio = xall[k-1];
    }
  }
  fprintf(final, "CHROMOSOME ORDER OF LOCI : \n");
  FORLIM = mlocus;
  for (i = 1; i <= FORLIM; i++) {
    j = 0;
    do {
      j++;
    } while (order[j - 1] != i);
    fprintf(final, "%3ld", j);
  }
  putc('\n', final);
  if (itsys != 0) {
    fprintf(final, "****************** FINAL VALUES **********************\n");
    fprintf(final, "PROVIDED FOR LOCUS %3ld (CHROMOSOME ORDER)\n",
	    order[itsys - 1]);
    fprintf(final, "******************************************************\n");
    WITH = thislocus[itsys - 1];
    if (disequi) {
      fprintf(final, "HAPLOTYPE FREQUENCIES:\n");
      FORLIM = nuhap;
      for (i = 0; i < FORLIM; i++)
	fprintf(final, "%9.6f", hapfreq->genarray[i]);
    } else {
      fprintf(final, "GENE FREQUENCIES :\n");
      FORLIM = WITH->nallele;
      for (i = 0; i < FORLIM; i++)
	fprintf(final, "%9.6f", WITH->freq[i]);
    }
    putc('\n', final);
    if (WITH->which == quantitative) {
      fprintf(final, "VALUES FOR GENOTYPE MEANS:\n");
      invert(WITH->UU.U1.vmat, WITH->UU.U1.ntrait, &WITH->UU.U1.det);
      FORLIM = WITH->UU.U1.ntrait;
      for (i = 0; i < FORLIM; i++) {
	FORLIM1 = WITH->nallele;
	for (j = 1; j <= FORLIM1; j++) {
	  FORLIM2 = WITH->nallele;
	  for (k = j - 1; k < FORLIM2; k++)
	    fprintf(final, "%6.3f ", WITH->UU.U1.pm[j][k][i]);
	}
	putc('\n', final);
      }
      fprintf(final, "COVARIANCE MATRIX:\n");
      FORLIM = WITH->UU.U1.ntrait;
      for (i = 1; i <= FORLIM; i++) {
	for (j = 0; j < i; j++)
	  fprintf(final, "%6.3f ", WITH->UU.U1.vmat[i - 1][j]);
	putc('\n', final);
      }
    } else {
      if (WITH->which == affection) {
	fprintf(final, "PENETRANCES:\n");
	FORLIM = WITH->UU.U0.nclass;
	for (l = 0; l < FORLIM; l++) {
	  FORLIM1 = WITH->nallele;
	  for (i = 1; i <= FORLIM1; i++) {
	    FORLIM2 = WITH->nallele;
	    for (j = i - 1; j < FORLIM2; j++)
	      fprintf(final, "%5.3f ", WITH->UU.U0.pen[i][j][2][l]);
	  }
	  putc('\n', final);
	  if (sexlink) {
	    FORLIM1 = WITH->nallele;
	    for (i = 0; i < FORLIM1; i++)
	      fprintf(final, "%5.3f ", WITH->UU.U0.pen[0][i][2][l]);
	  }
	  putc('\n', final);
	}
      }
    }
  }
  fprintf(final, "******************************************************\n");
  if (interfer && !mapping) {
    fprintf(final, "P VALUES:\n");
    WITH1 = maletheta;
    memcpy(oldtheta, WITH1->theta, sizeof(thetarray));
    WITH1->theta[0] =
      (oldtheta[0] + oldtheta[mlocus - 1] - oldtheta[mlocus - 2]) / 2.0;
    WITH1->theta[mlocus - 2] =
      (oldtheta[mlocus - 2] + oldtheta[mlocus - 1] - oldtheta[0]) / 2.0;
    WITH1->theta[mlocus - 1] =
      (oldtheta[0] + oldtheta[mlocus - 2] - oldtheta[mlocus - 1]) / 2.0;
    FORLIM = mlocus;
    for (i = 0; i < FORLIM; i++)
      fprintf(final, " %5.3f", WITH1->theta[i]);
    putc('\n', final);
    memcpy(WITH1->theta, oldtheta, sizeof(thetarray));
    if (sexdif) {
      WITH1 = femaletheta;
      memcpy(oldtheta, WITH1->theta, sizeof(thetarray));
      WITH1->theta[0] =
	(oldtheta[0] + oldtheta[mlocus - 1] - oldtheta[mlocus - 2]) / 2.0;
      WITH1->theta[mlocus - 2] =
	(oldtheta[mlocus - 2] + oldtheta[mlocus - 1] - oldtheta[0]) / 2.0;
      WITH1->theta[mlocus - 1] =
	(oldtheta[0] + oldtheta[mlocus - 2] - oldtheta[mlocus - 1]) / 2.0;
      fprintf(final, "FEMALE: \n");
      FORLIM = mlocus;
      for (i = 0; i < FORLIM; i++)
	fprintf(final, " %5.3f", WITH1->theta[i]);
      putc('\n', final);
      memcpy(WITH1->theta, oldtheta, sizeof(thetarray));
    }
    if (ivar == 1 && icall == 0) {
      for (i = 0; i <= 2; i++) {
	TEMP = 1 + exp(xall[nall - i - 1]);
	deriv[i] = exp(xall[nall - i - 1]) / (TEMP * TEMP);
      }
      k = -1;
      for (i = 0; i <= 2; i++) {
	if (itp[nall - i - 1] == 1) {
	  k++;
	  l = -1;
	  for (j = 0; j <= 2; j++) {
	    if (itp[nall - j - 1] == 1) {
	      l++;
	      bmat[n - k - 1][n - l - 1] *= deriv[j] * deriv[i];
	    }
	  }
	  if (l == -1)
	    l = 0;
	  for (j = n - l - 1; j >= 0; j--) {
	    bmat[j][n - k - 1] *= deriv[i];
	    bmat[n - k - 1][j] = bmat[j][n - k - 1];
	  }
	}
      }
    }
  }
  fprintf(final, "THETAS:\n");
  WITH1 = maletheta;
  FORLIM = mlocus - 2;
  
  for (i = 0; i <= FORLIM; i++)
    fprintf(final, " %5.3f", WITH1->theta[i]);
    
  if (interfer)
    fprintf(final, " %5.3f\n", maletheta->theta[mlocus - 1]);
  else
    putc('\n', final);
  if (sexdif) {
    fprintf(final, "FEMALE:\n");
    WITH1 = femaletheta;
    FORLIM = mlocus - 2;
    for (i = 0; i <= FORLIM; i++)
      fprintf(final, " %5.3f", WITH1->theta[i]);
    if (interfer)
      fprintf(final, " %5.3f\n", femaletheta->theta[mlocus - 1]);
    else
      putc('\n', final);
    if (readfemale) {
      fprintf(final, "FEMALE/MALE DIST RATIO :\n");
      FORLIM = mlocus - 2;
      for (i = 0; i <= FORLIM; i++) {
	dist = getdist(&maletheta->theta[i]);
	if (dist != 0.0)
	  dist = getdist(&femaletheta->theta[i]) / dist;
	else
	  dist = 0.0;
	fprintf(final, " %5.3f", dist);
      }
      if (interfer) {
	dist = getdist(&maletheta->theta[mlocus - 1]);
	if (dist != 0.0)
	  dist = getdist(&femaletheta->theta[mlocus - 1]) / dist;
	else
	  dist = 0.0;
	fprintf(final, " %5.3f", dist);
      }
      putc('\n', final);
    } else {
      fprintf(final, "CONSTANT FEMALE/MALE DIST RATIO :\n");
      fprintf(final, "%5.3f\n", distratio);
    }
  }
  fprintf(final, "******************************************************\n");
  fprintf(final, "-2 LN(LIKE) = % .5le\n", f);
  getlods(&V.lods, x, &V);
  if (mlocus != 2)
    fprintf(final, "OTTS GENERALIZED LOD SCORE =% .5le\n",
	    (V.lods - f) / (2 * log10_));
  else
    fprintf(final, "LOD SCORE =% .5le\n", (V.lods - f) / (2 * log10_));
  fprintf(final, "NUMBER OF ITERATIONS = %5ld\n", nit);
  fprintf(final, "NUMBER OF FUNCTION EVALUATIONS = %5ld\n", nfe);
  fprintf(final, "PTG = % .5le\n", ptg);
  idg++;
  fprintf(outfile, "EXIT CONDITION%12ld\n", idg);
  switch (idg) {

  case 1:
    fprintf(outfile, "Maximum possible accuracy reached\n");
    break;

  case 2:
    fprintf(outfile, "Search direction no longer downhill\n");
    break;

  case 3:
    fprintf(outfile,
	    "Accumulation of rounding error prevents further progress\n");
    break;

  case 4:
    fprintf(outfile,
      "All significant differences lost through cancellation in conditioning\n");
    break;

  case 5:
    fprintf(outfile, "Specified tolerance on normalized gradient met\n");
    break;

  case 6:
    fprintf(outfile, "Specified tolerance on gradient met\n");
    break;

  case 7:
    fprintf(outfile, "Maximum number of iterations reached\n");
    break;

  case 8:
    fprintf(outfile, "Excessive cancellation in gradient\n");
    break;
  }
  if (ivar == 1 && icall == 0) {
    fprintf(final, "VARIANCE-COVARIANCE OF THE ESTIMATES\n");
    if (interfer && !mapping)
      fprintf(final, "VALUES GIVEN FOR P VALUES\n");
    FORLIM = n;
    for (i = 0; i < FORLIM; i++) {
      FORLIM1 = n;
      for (j = 0; j < FORLIM1; j++)
	fprintf(final, "% .5le ", 2.0 * bmat[i][j]);
      putc('\n', final);
    }
  }
  fprintf(final, "******************************************************\n");
  fprintf(final, "******************************************************\n");
  if (dostream)
    streamout(&V);
}


/* Local variables for step: */
struct LOC_step {			/* If this structure is changed,    */
  enum {				/* some code may have to be altered */
    go_, exit1, exit2, exit3		/* in performCheckpoint() and       */
  } cont;				/* recoverCheckpoint ().            */
} ;					/*                        --Shriram */

Local Void firststep(LINK)
struct LOC_step *LINK;
{
  int i, FORLIM;

#if !defined(DOS)
  /* Shriram: begin */
  if ( normalRun == checkpointStatus )
  {
    funCallPath = fCP_gem_iter_st_first ;
    performCheckpoint ( functionLocation , funCallPath , LINK -> cont ) ;
  }
  else
  {
    if ( ( functionLocation == ckptTuple . ckptLocation )
        && ( fCP_gem_iter_st_first == ckptTuple . ckptAttribute ) )
    {
      recoverCheckpoint ( & LINK -> cont ) ;
      funCallPath = ckptTuple . ckptAttribute ;
    }
  }
  /* Shriram: end */
#endif

  sumt = 0.0;
  idg = 0;
  FORLIM = n;
  for (i = 0; i < FORLIM; i++) {
   /*printf("first %lf %lf %lf %lf\n", xt[i], x[i], t, p[i]);*/ /*Alex*/
    xt[i] = x[i] + t * p[i];
  }
  fun(&ft, xt);
  /*save likelihood values for later*/
  for (i=0; i < nuped; i++)
    outsavelike[i] = likebyped[i];
  nfe++;
}

Local Void decreaset(LINK)
struct LOC_step *LINK;
{
  int i, FORLIM;

#if !defined(DOS)
  /* Shriram: begin */
  if ( normalRun == checkpointStatus )
  {
    funCallPath = fCP_gem_iter_st_decr ;
    performCheckpoint ( functionLocation , funCallPath , LINK -> cont ) ;
  }
  else
  {
    if ( ( functionLocation == ckptTuple . ckptLocation )
        && ( fCP_gem_iter_st_decr == ckptTuple . ckptAttribute ) )
    {
      recoverCheckpoint ( & LINK -> cont ) ;
      funCallPath = ckptTuple . ckptAttribute ;
    }
  }
  /* Shriram: end */
#endif

  /*Modified by M. Lathrop 29/04/86 to trap tmin=0 problem*/
  if (tmin < small)
    tmin = small;
  LINK->cont = go_;
  while (LINK->cont == go_) {
    if (f - fsave == 0.0 && idif == 2)
      idg = 2;
    if (t < tmin) {
      LINK->cont = exit3;
      break;
    }
    t = 0.5 * t;
    FORLIM = n;
    for (i = 0; i < FORLIM; i++)
      xt[i] = x[i] + t * p[i];
    fun(&ft2, xt);
    nfe++;
    if (ft2 < f) {
      sumt += t;
      f = ft2;
      idg = 0;
      memcpy(x, xt, sizeof(vector));
      /*save likelihood values for later*/
      for (i=0; i < nuped; i++)
        outsavelike[i] = likebyped[i];
      if (sumt < tmin)
	LINK->cont = exit3;
      else
	LINK->cont = exit2;
      continue;
    }
    if (ft + f - ft2 - ft2 <= 0.0)
      scal = 0.1;
    else {
      scal = 1.0 + 0.5 * (f - ft) / (f + ft - ft2 - ft2);
      if (scal < 0.1)
	scal = 0.1;
    }
    t = scal * t;
    FORLIM = n;
    for (i = 0; i < FORLIM; i++)
      xt[i] = x[i] + t * p[i];
    fun(&ft, xt);
    nfe++;
    if (f <= ft)
      continue;
    sumt += t;
    idg = 0;
    f = ft;
    memcpy(x, xt, sizeof(vector));
    /*save likelihood values for later*/
    for (i=0; i < nuped; i++)
      outsavelike[i] = likebyped[i];
    if (t < tmin)
      LINK->cont = exit3;
    else
      LINK->cont = exit2;
  }
}

Local Void increaset(LINK)
struct LOC_step *LINK;
{
  int i, j, FORLIM;

#if !defined(DOS)
  /* Shriram: begin */
  if ( normalRun == checkpointStatus )
  {
    funCallPath = fCP_gem_iter_st_incr ;
    performCheckpoint ( functionLocation , funCallPath , LINK -> cont ) ;
  }
  else
  {
    if ( ( functionLocation == ckptTuple . ckptLocation )
        && ( fCP_gem_iter_st_incr == ckptTuple . ckptAttribute ) )
    {
      recoverCheckpoint ( & LINK -> cont ) ;
      funCallPath = ckptTuple . ckptAttribute ;
    }
  }
  /* Shriram: end */
#endif

  sumt += t;
  LINK->cont = go_;
  while (LINK->cont == go_) {
    twot = t + t;
    if (ibnd > 0 && tbnd >= 0.0 && twot > tbnd) {
      f = ft;
      memcpy(x, xt, sizeof(vector));
      fprintf(outfile, "****** ACTIVE BOUNDARY CONSTRAINT *****\n");
      LINK->cont = exit1;
      continue;
    }
    memcpy(x, xt, sizeof(vector));
    FORLIM = n;
    for (i = 0; i < FORLIM; i++)
      xt[i] = x[i] + t * p[i];
    fun(&f2t, xt);
    nfe++;
    if (f2t > ft) {
      f = ft;
      LINK->cont = exit2;
      break;
    }
    if (f2t + f - ft - ft < 0.0 || ft - f2t + curv * (ft - f) >= 0.0) {
      sumt += t;
      t += t;
      ft = f2t;
    } else {
      sumt += t;
      f = f2t;
      memcpy(x, xt, sizeof(vector));
      /*save likelihood values for later*/
      for (j=0; j < nuped; j++)
        outsavelike[i] = likebyped[i];
      LINK->cont = exit2;
    }
  }
}

static Void step()
{
  struct LOC_step V;

#if !defined(DOS)
  /* Shriram: begin */
  int    ckptFunCallPath = ckptTuple . ckptAttribute ;

  if ( ( checkpointedRun == checkpointStatus ) &&
      ( functionLocation == ckptTuple . ckptLocation ) )
  {
    if ( fCP_gem_iter_st_first == ckptFunCallPath )
      goto FunRecover3 ;
    else if ( fCP_gem_iter_st_decr == ckptFunCallPath )
      goto FunRecover4 ;
    else if ( fCP_gem_iter_st_incr == ckptFunCallPath )
      goto FunRecover5 ;
  }
  /* Shriram: end */
#endif
  firstapprox = true;
 FunRecover3:
  firststep(&V);
  if (f > ft)
 FunRecover5:
    increaset(&V);
  else
 FunRecover4:
    decreaset(&V);
  if (V.cont == exit2)
    t = sumt;
  if (V.cont == exit3) {
    iret = 1;
    return;
  }
  memcpy(gsave, g, sizeof(vector));
  isw = 0;
  iret = 2;
}


/* Local variables for update: */
struct LOC_update {
  int j, jm1;
  boolean cont, dfp;
} ;

Local Void prep(LINK)
struct LOC_update *LINK;
{
  int i, k, FORLIM;

  iswup = 1;
  wtil[0] = -gsave[0];
  ztil[0] = y[0];
  FORLIM = n;
  for (k = 2; k <= FORLIM; k++) {
    wtil[k - 1] = -gsave[k - 1];
    ztil[k - 1] = y[k - 1];
    for (i = 0; i <= k - 2; i++) {
      tl[i][k - 1] = tl[k - 1][i];
      wtil[k - 1] -= tl[k - 1][i] * wtil[i];
      ztil[k - 1] -= tl[k - 1][i] * ztil[i];
    }
  }
  sb = 0.0;
  sc = 0.0;
  sd = 0.0;
  FORLIM = n;
  for (i = 0; i < FORLIM; i++) {
    sb += t * ztil[i] * wtil[i] / d[i];
    sc += t * t * wtil[i] * wtil[i] / d[i];
    sd += ztil[i] * ztil[i] / d[i];
  }
}

Local Void sr1update(LINK)
struct LOC_update *LINK;
{
  alpha = -1.0;
  if (sc - sb == 0.0 || sb - sd == 0.0 || sd - 2.0 * sb + sc == 0.0) {
    idg = 3;
    return;
  }
  sa = 1.0 / (sqrt(fabs(sc - sb)) * sqrt(fabs(sb - sd)));
  thet1 = ((sd - sb) * sa + 1.0) / (2.0 * sb - sd - sc);
  thet2 = sa + (sa * (sb - sc) + 1.0) / (sd - 2.0 * sb + sc);
}

Local Void sr2update(LINK)
struct LOC_update *LINK;
{
  aa = sb / sc - 2.0 * (sd / sb) + sd / sc;
  bb = sb / sc - 1.0;
  cc = 1.0 - sb / sd;
  del2 = bb * bb - aa * cc;
  LINK->dfp = true;
  if (del2 > 0.00000001) {
    LINK->dfp = false;
    del = sqrt(del2);
    alph1 = (del - bb) / aa;
    alph2 = (-bb - del) / aa;
    if (fabs(alph1) < fabs(alph2))
      alpha = alph1;
    else
      alpha = alph2;
    if (fabs(alpha) < 0.00001)
      LINK->dfp = true;
    else {
      sa = (alpha + 1.0) * (alpha + 1.0) + sc / sb -
	   alpha * alpha * (sc / sb) * (sd / sb) - 1.0 + alpha * alpha * sd / sb;
      if (sa <= 0.0)
	sa = 0.0;
      else {
	sa = sqrt(sa);
	sa = 1.0 / (sa * sb);
      }
      rdiv = 1.0 / (alpha * alpha * sd + 2.0 * alpha * sb + sc);
      thet1 = -(sa * alpha * (alpha * sd + sb) + 1.0) * rdiv;
      thet2 = sa + (alpha * sa * (sc + alpha * sb) - alpha) * rdiv;
    }
  }
  if (!LINK->dfp)
    return;
  alpha = 0.0;
  sa = 1.0 / (sqrt(sb) * sqrt(sc));
  thet1 = -1.0 / sc;
  thet2 = sa;
}

Local Void getwzs(LINK)
struct LOC_update *LINK;
{
  int i, k, FORLIM;

  FORLIM = n;
  for (i = 0; i < FORLIM; i++) {
    w[i] = t * wtil[i] + alpha * ztil[i];
    z[i] = t * thet1 * wtil[i] + thet2 * ztil[i];
  }
  s[n - 1] = 0.0;
  FORLIM = n;
  for (k = 1; k < FORLIM; k++) {
    LINK->j = n - k + 1;
    LINK->jm1 = LINK->j - 1;
    s[LINK->jm1 - 1] = s[LINK->j - 1] +
		       w[LINK->j - 1] * w[LINK->j - 1] / d[LINK->j - 1];
  }
  nu = 1.0;
  eta = 0.0;
}

Local Void recur(LINK)
struct LOC_update *LINK;
{
  int i, k, FORLIM;

  LINK->cont = true;
  while (LINK->cont) {
    LINK->cont = false;
    if (iswup < 2) {
      FORLIM = n;
      for (i = 0; i < FORLIM; i++)
	wtjp1[i] = -gsave[i];
      memcpy(ztjp1, y, sizeof(vector));
    } else {
      FORLIM = n;
      for (i = 0; i < FORLIM; i++) {
	wtil[i] = alpha * y[i] - t * g[i];
	ztil[i] = thet2 * y[i] - t * thet1 * g[i];
      }
    }
    LINK->j = 0;
    lambj2 = 0.0;
    while (LINK->j < n - 1) {
      LINK->j++;
      if (iswup < 2) {
	FORLIM = n;
	for (k = LINK->j; k < FORLIM; k++) {
	  wtjp1[k] -= wtil[LINK->j - 1] * tl[k][LINK->j - 1];
	  ztjp1[k] -= ztil[LINK->j - 1] * tl[k][LINK->j - 1];
	}
      } else {
	FORLIM = n;
	for (k = LINK->j; k < FORLIM; k++) {
	  wtjp1[k] = wtil[k] - w[LINK->j - 1] * tl[k][LINK->j - 1];
	  ztjp1[k] = ztil[k] - z[LINK->j - 1] * tl[k][LINK->j - 1];
	}
      }
      aj = nu * z[LINK->j - 1] - eta * w[LINK->j - 1];
      thj = 1.0 + aj * w[LINK->j - 1] / d[LINK->j - 1];
      lambj2 = thj * thj + aj * aj * s[LINK->j - 1] / d[LINK->j - 1];
      if (iswup < 2) {
	if (lambj2 > 10.0) {
	  LINK->cont = true;
	  iswup = 2;
	  FORLIM = n;
	  for (k = 2; k <= FORLIM; k++) {
	    for (i = 0; i <= k - 2; i++)
	      tl[k - 1][i] = tl[i][k - 1];
	  }
	  LINK->j = n;
	}
      }
      if (LINK->cont)
	continue;
      dp[LINK->j - 1] = d[LINK->j - 1] * lambj2;
      lambj = sqrt(lambj2);
      if (thj > 0.0)
	lambj = -lambj;
      muj = thj - lambj;
      bj = thj * w[LINK->j - 1] + aj * s[LINK->j - 1];
      gamlj = bj * nu / (lambj2 * d[LINK->j - 1]);
      betlj = (aj - bj * eta) / (lambj2 * d[LINK->j - 1]);
      nu = -(nu / lambj);
      eta = -((eta + aj * aj / (muj * d[LINK->j - 1])) / lambj);
      if (iswup < 2) {
	FORLIM = n;
	for (k = LINK->j; k < FORLIM; k++)
	  tl[k][LINK->j - 1] += t * (betlj + thet1 * gamlj) * wtjp1[k] +
				(alpha * betlj + thet2 * gamlj) * ztjp1[k];
      } else {
	FORLIM = n;
	for (k = LINK->j; k < FORLIM; k++) {
	  tl[k][LINK->j - 1] = tl[k][LINK->j - 1] / lambj2 + betlj * wtil[k] +
			       gamlj * ztil[k];
	  wtil[k] = wtjp1[k];
	  ztil[k] = ztjp1[k];
	}
      }
    }
  }
  aj = nu * z[n - 1] - eta * w[n - 1];
  lambj = 1.0 + aj * w[n - 1] / d[n - 1];
  dp[n - 1] = d[n - 1] * lambj * lambj;
  memcpy(d, dp, sizeof(vector));
}


/*UPDATE*/
static Void update()
{
  struct LOC_update V;

  prep(&V);
  if (sb < sc)
    fbcd = sb;
  else
    fbcd = sc;
  if (sd < fbcd)
    fbcd = sd;
  if (fbcd > small) {
    fbcd = 2.0 * sc * (sd / sb) / (sc + sd);
    if (fbcd < 1.0)
      sr1update(&V);
    else
      sr2update(&V);
  } else
    sr1update(&V);
  if (idg != 3) {
    getwzs(&V);
    recur(&V);
  }
}


Local Void bldlt(b)
vector *b;
{
  int i, j, ic, k;
  double su, tlic, ff, hh, temp, temp1;
  matrix s;
  vector tu, xvec;
  int FORLIM, FORLIM1, FORLIM2;

  firstapprox = true;
  ff = f;
  if (icall > 0 && ihess <= 0) {
    memcpy(xvec, xall, sizeof(vector));
    hh = 10 * h;
    FORLIM = n;
    for (i = 0; i < FORLIM; i++) {
      temp = x[i];
      x[i] += hh;
      fun(&f, x);
      tu[i] = f;
      for (j = 0; j < i + 1; j++) {
	temp1 = x[j];
	x[j] += hh;
	fun(&f, x);
	b[i][j] = f;
	b[j][i] = b[i][j];
	x[j] = temp1;
      }
      x[i] = temp;
    }
    memcpy(xall, xvec, sizeof(vector));
    FORLIM = n;
    for (i = 0; i < FORLIM; i++) {
      for (j = 0; j < i + 1; j++) {
	b[i][j] = (ff + b[i][j] - tu[i] - tu[j]) / (hh * hh);
	b[j][i] = b[i][j];
      }
    }
  }
  ic = 1;
  while (b[ic - 1][ic - 1] > 0.0 && ic <= n) {
    temp = b[ic - 1][ic - 1];
    d[ic - 1] = b[ic - 1][ic - 1];
    tl[ic - 1][ic - 1] = 1.0;
    if (ic != n) {
      FORLIM = n;
      for (k = ic; k < FORLIM; k++)
	tl[k][ic - 1] = b[k][ic - 1] / temp;
      FORLIM = n;
      for (i = ic; i < FORLIM; i++) {
	tlic = tl[i][ic - 1];
	FORLIM1 = n;
	for (k = i; k < FORLIM1; k++)
	  b[k][i] -= tlic * tl[k][ic - 1] * temp;
      }
    }
    ic++;
  }
  if (ic > n) {
    icall--;
    fprintf(outfile, "FACTORIZATION SUCCEEDED\n");
  } else {
    icall++;
    fprintf(outfile, "FACTORIZATION FAILED\n");
  }
  if (icall != 0)
    return;
  s[0][0] = 1.0;
  FORLIM = n;
  for (i = 2; i <= FORLIM; i++) {
    for (k = 0; k <= i - 2; k++) {
      su = 0.0;
      for (j = k; j <= i - 2; j++)
	su += tl[i - 1][j] * s[j][k];
      s[i - 1][k] = -su;
    }
    s[i - 1][i - 1] = 1.0;
  }
  FORLIM = n;
  for (i = 0; i < FORLIM; i++) {
    FORLIM1 = n;
    for (j = i; j < FORLIM1; j++) {
      su = 0.0;
      FORLIM2 = n;
      for (k = j; k < FORLIM2; k++)
	su += s[k][i] * s[k][j] / d[k];
      b[i][j] = su;
      b[j][i] = su;
    }
  }
  fprintf(outfile, "B-MATRIX\n");
  FORLIM = n;
  for (i = 1; i <= FORLIM; i++) {
    for (j = 0; j < i; j++)
      fprintf(outfile, "% .5le", b[i - 1][j]);
    putc('\n', outfile);
  }
  if (ivar != 1)
    return;
  k = 0;
  FORLIM = nall;
  for (i = 0; i < FORLIM; i++) {
    se[i] = 0.0;
    if (itp[i] == 1) {
      k++;
      se[i] = sqrt(b[k - 1][k - 1]);
    }
  }
}

Local Void inib()
{
  int i, j;
  FILE *in1;
  int FORLIM;

  in1 = NULL;
  if (icall != 1 && ihess > 1) {
    if (in1 != NULL)
      in1 = freopen("in1.dat", "r", in1);
    else
      in1 = fopen("in1.dat", "r");
    if (in1 == NULL)
      exit(FileNotFound);
    FORLIM = n;
    for (i = 1; i <= FORLIM; i++) {
      for (j = 0; j < i; j++)
	fscanf(in1, "%lf", &bmat[i - 1][j]);
    }
    FORLIM = n;
    for (i = 0; i < FORLIM; i++) {
      for (j = 0; j < i + 1; j++)
	bmat[j][i] = bmat[i][j];
    }
  }
  bldlt(bmat);
  if (icall == 1) {
    FORLIM = n;
    for (i = 0; i < FORLIM; i++) {
      for (j = 0; j < i + 1; j++)
	tl[i][j] = 0.0;
      tl[i][i] = 1.0;
      d[i] = 1.0;
    }
  }
  if (in1 != NULL)
    fclose(in1);
}

Local Void initialize()
{
  int i, j, FORLIM;

#if !defined(DOS)
  /* Shriram: begin */
  if ( normalRun == checkpointStatus )
  {
    funCallPath = fCP_gem_init ;
    performCheckpoint ( functionLocation , funCallPath , 0 ) ;
  }
  else
  {
    if ( ( functionLocation == ckptTuple . ckptLocation )
        && ( fCP_gem_init == ckptTuple . ckptAttribute ) )
    {
      recoverCheckpoint ( NULL ) ;
      funCallPath = ckptTuple . ckptAttribute ;
    }
  }
  /* Shriram: end */
#endif
  firstapprox = true;
  fprintf(outfile, "DIFFER INTER = % .5le TRUNC UPPER = % .5le\n", h, trupb);
  nit = 0;
  idg = 0;
  idif = 1;
  t = 0.1;
  tmin = 0.0;
  fsmf = 0.0;
  fun(&f, x);
  /*save likelihood values for later*/
  for (i=0; i < nuped; i++)
    outsavelike[i] = likebyped[i];
  nfe = 1;
  FORLIM = n;
  for (i = 0; i < FORLIM; i++) {
    for (j = 0; j < i + 1; j++)
      tl[i][j] = 0.0;
    tl[i][i] = 1.0;
    d[i] = 1.0;
  }
  if (ihess <= 0)
    return;
  icall = 0;
  inib();
  icall = 1;
  t = 1.0;
}

Local Void gforward()
{
  int i, k, FORLIM;
  /*The following declaration was introduced by A. A. Schaffer 7/27/94
  to fix a bug where the wrong xall was being used to move to the
  next theta*/
  double xallsave[maxn]; /*Save value of xall to recover it later*/

  firstapprox = true;
  FORLIM = n;
  outsavefvalue = f;
  saveparams();
  k = 0;
  for (i = 0; i < maxn; i++) {
    if (itp[i] == 1) {
      k++;
      xallsave[i] = outsavex[k - 1];
    }
    else
      xallsave[i] = xall[i];
  }
  for (i = 0; i < FORLIM; i++) {
#if !defined(DOS)
    /* Shriram: begin */
    if ( normalRun == checkpointStatus )
    {
      if ( fCP_gem_iter_ == funCallPath )
        funCallPath = fCP_gem_iter_gfor ;
      else if ( fCP_gem_ == funCallPath )
        funCallPath = fCP_gem_gfor ;
      performCheckpoint ( functionLocation , funCallPath , i ) ;
    }
    else if ( ( functionLocation == ckptTuple . ckptLocation ) &&
             ( ( fCP_gem_iter_gfor == ckptTuple . ckptAttribute ) ||
              ( fCP_gem_gfor == ckptTuple . ckptAttribute ) ) )
    {
      recoverCheckpoint ( & i ) ;
      FORLIM = n ;
      if ( i >= FORLIM )
        break ;
      funCallPath = ckptTuple . ckptAttribute ;
    }
    /* Shriram: end */
#endif
    hx = h;
    if (ihx == 1) {
      if (fabs(x[i]) < 0.1)
	hx = h * 0.1;
      else
	hx = h * fabs(x[i]);
    }
    xsave = x[i];
    x[i] += hx;
    if (continue_ != quit) {
      fun(&fxph, x);
      savedf[i] = fxph;
    }
    else
      fxph = savedf[i];
    firstapprox = false;
    g[i] = (fxph - f) / hx;
    x[i] = xsave;
  }
  for (i = 0; i < maxn; i++)
    xall[i] = xallsave[i];
  if (continue_ != quit)
    nfe += n;
}

Local Void gcentral()
{
  int i, k, FORLIM;
  double xallsave[maxn]; /*see similar declaration in gforward, above*/

  firstapprox = true;
  FORLIM = n;
  outsavefvalue = f;
  saveparams();
  k = 0;
  for (i = 0; i < maxn; i++) {
    if (itp[i] == 1) {
      k++;
      xallsave[i] = outsavex[k - 1];
    }
    else
      xallsave[i] = xall[i];
  }
  for (i = 0; i < FORLIM; i++) {
#if !defined (DOS)
    /* Shriram: begin */
    if ( normalRun == checkpointStatus ) 
    {
      if ( fCP_gem_iter_gcen1_ == funCallPath )
        funCallPath = fCP_gem_iter_gcen1 ;
      else if ( fCP_gem_iter_gcen2_ == funCallPath )
        funCallPath = fCP_gem_iter_gcen2 ;
      performCheckpoint ( functionLocation , funCallPath, i ) ;
    }
    else if ( ( functionLocation == ckptTuple . ckptLocation ) &&
             ( ( fCP_gem_iter_gcen1 == ckptTuple . ckptAttribute ) ||
              ( fCP_gem_iter_gcen2 == ckptTuple . ckptAttribute ) ) )
    {
      recoverCheckpoint ( & i ) ;
      FORLIM = n ;
      if ( i >= FORLIM )
        break ;
      funCallPath = ckptTuple . ckptAttribute ;
    }
    /* Shriram: end */
#endif
    hx = h;
    if (ihx == 1) {
      if (fabs(x[i]) < 0.1)
	hx = h * 0.1;
      else
	hx = h * fabs(x[i]);
    }
    xsave = x[i];
    x[i] += hx;
    if (continue_ != quit) {
       fun(&fxph, x);
       savedf[2*i] = fxph;
    }
    else
      fxph = savedf[2*i];
    firstapprox = false;
    x[i] = xsave - hx;
    if (continue_ != quit) {
       fun(&fxmh, x);
       savedf[2 * i + 1] = fxmh;
    }
    else
      fxmh = savedf[2 * i + 1];
    g[i] = (fxph - fxmh) / (hx + hx);
    x[i] = xsave;
  }
  for (i = 0; i < maxn; i++)
     xall[i] = xallsave[i];
  if (continue_ != quit)
    nfe += 2 * n ;			/* Shriram */
}

Local Void getp()
{
  int i, nmj, j, FORLIM, FORLIM1;

  nit++;
  fsav2 = fsave;
  fsave = f;
  fprintf(outfile, "\nITERATION %5ld T = %10.3f NFE = %5ld F = % .5le\n",
	  nit, t, nfe, f);
  printf("ITERATION %5ld T = %10.3f NFE = %5ld F = % .5le\n", nit, t, nfe, f);
  if (nit > maxit) {
    idg = 6;
    return;
  }
  if (n < 20) {
    fprintf(outfile, "X= ");
    outcontrol(x);
    fprintf(outfile, "G= ");
    outcontrol(g);
  }
  if (n == 1)
    p[0] = -(g[0] / d[0]);
  else {
    v[0] = -g[0];
    FORLIM = n;
    for (i = 2; i <= FORLIM; i++) {
      v[i - 1] = -g[i - 1];
      for (j = 0; j <= i - 2; j++)
	v[i - 1] -= tl[i - 1][j] * v[j];
    }
    p[n - 1] = v[n - 1] / d[n - 1];
    FORLIM = n;
    for (j = 1; j < FORLIM; j++) {
      nmj = n - j;
      p[nmj - 1] = v[nmj - 1] / d[nmj - 1];
      FORLIM1 = n;
      for (i = nmj; i < FORLIM1; i++)
	p[nmj - 1] -= tl[i][nmj - 1] * p[i];
    }
  }
  fprintf(outfile, "P= ");
  outcontrol(p);
  if (ibnd == 1) {
    FORLIM = n;
    for (i = 0; i < FORLIM; i++) {
      if (p[i] == 0.0)
	idg = 7;
    }
  }
  if (idg == 7)
    return;
  ptg = 0.0;
  FORLIM = n;
  for (i = 0; i < FORLIM; i++)
    ptg += p[i] * g[i];
  if (ptg >= 0.0) {
    idg = 1;
    return;
  }
  if (fabs(ptg) < tol) {
    idg = 4;
    fsmf = fsav2 - f;
    fprintf(outfile, "FSMF = % .5le PTG = % .5le TMIN = % .5le\n",
	    fsmf, ptg, tmin);
    return;
  }
  xpm = xpmcon;
  if (nit == 1)
    return;
  FORLIM = n;
  for (i = 0; i < FORLIM; i++) {
    if (ihx == 1) {
      xp = fabs(x[i]);
      if (xp < 0.1)
	xp = 0.1;
      xp /= fabs(p[i]);
      if (xp < xpm)
	xpm = xp;
    } else {
      if (xpm > fabs(1.0 / p[i]))
	xpm = fabs(1.0 / p[i]);
    }
  }
  tmin = 0.5 * xpm * h / trupb;
  if (idif == 2)
    t = 1.0;
  else {
    t = -2.0 * fsmf / ptg;
    if (t <= 0.0)
      t = 1.0;
    else {
      if (1.0 < t)
	t = 1.0;
    }
  }
  fprintf(outfile, "FSMF = % .5le PTG =% .5le TMIN=% .5le\n", fsmf, ptg, tmin);
  fprintf(outfile, "INITIAL T = % .5le\n", t);
}

Local Void getytp()
{
  int i, FORLIM;

  ytp = 0.0;
  fsmf = fsave - f;
  FORLIM = n;
  for (i = 0; i < FORLIM; i++) {
    y[i] = g[i] - gsave[i];
    ytp += y[i] * p[i];
  }
}

Local Void chkbnd()
{
  /*This procedure modified by M. Lathrop 29/04/86
    with the introduction of continue and repeat*/
  int i, j, k, ik, ii, jk;
  boolean continue_;
  int FORLIM, FORLIM1;

  do {
    clb = clbcon;
    for (j = 1; j <= 2; j++) {
      k = 0;
      FORLIM1 = nall;
      for (i = 1; i <= FORLIM1; i++) {
	if (itp[i - 1] == 1) {
	  k++;
	  check = fabs(x[k - 1] - bnd[i - 1][j - 1]);
	  if (check <= clb) {
	    clb = check;
	    ik = k;
	    ii = i;
	    jk = j;
	  }
	}
      }
    }
    if (clb < 0.1 * h) {
      ihess = 0;
      if (jk == 2)
	eh = -h - h;
      else
	eh = h + h;
      x[ik - 1] = bnd[ii - 1][jk - 1] + eh;
      xall[ii - 1] = x[ik - 1];
      itp[ii - 1] = 0;
      /*slide down saved x*/
      for(k = ik; k < nall; k++)
        outsavex[k-1] = outsavex[k];
      k = 0;
      FORLIM1 = nall;
      for (i = 0; i < FORLIM1; i++) {
	if (itp[i] == 1) {
	  k++;
	  x[k - 1] = xall[i];
	}
      }
      n = k;
      tol = tolconst;
      /*      tol:=tolconst*sqrt(n);*/
      continue_ = true;
      fprintf(outfile, "******* A VARIABLE WAS SET TO A BOUND ***********\n");
    } else {
      tbnd = tbndcon;
      for (j = 0; j <= 1; j++) {
	k = 0;
	FORLIM = nall;
	for (i = 0; i < FORLIM; i++) {
	  if (itp[i] == 1) {
	    k++;
	    teq = (bnd[i][j] - x[k - 1]) / p[k - 1];
	    if (teq < tbnd && teq > 0.0)
	      tbnd = teq;
	  }
	}
      }
      continue_ = false;
      if (t * (2.0 + h) >= tbnd)
	t = tbnd * (0.5 - h);
      fprintf(outfile, "TBND = % .5le RESET T = % .5le\n", tbnd, t);
    }
  } while (continue_);
}

Local Void iterate()
{
#if !defined(DOS)
  /* Shriram: begin */
  int  ckptFunCallPath = ckptTuple . ckptAttribute ;

  if ( normalRun == checkpointStatus )
    funCallPath = fCP_gem_iter_ ;
  /* Shriram: end */
#endif
  while (continue_ == go) {
    active = false;
#if !defined(DOS)
    /* Shriram: begin */
    if ( checkpointedRun == checkpointStatus )
    {
      if ( iterationLocation == ckptTuple . ckptLocation )
        checkpointStatus = normalRun ;
      else      /* checkpointedRun AND functionLocation */
      {
        if ( ( fCP_gem_iter_st_first == ckptFunCallPath ) ||
            ( fCP_gem_iter_st_decr == ckptFunCallPath ) ||
            ( fCP_gem_iter_st_incr == ckptFunCallPath ) )
          goto FunRecover345 ;
        else if ( fCP_gem_iter_gcen1 == ckptFunCallPath )
          goto FunRecover6 ;
        else if ( fCP_gem_iter_gfor == ckptFunCallPath )
          goto FunRecover7 ;
        else if ( fCP_gem_iter_gcen2 == ckptFunCallPath )
          goto FunRecover9 ;
      }
    }
    else
      performCheckpoint ( iterationLocation , nit , 0 ) ;
#endif
    getp();
    /* Shriram: end */
    if (idg != 0)
      continue_ = quit;
    else  {   /*idg == 0*/          /*Shriram*/
      iret = 2;
      if (ibnd == 1)
	chkbnd();
      if (n == 0)
	continue_ = quit;
      if (!active)
      {
#if !defined(DOS)
      FunRecover345:
#endif
	step();
#if !defined(DOS)
	funCallPath = fCP_gem_iter_ ;
#endif
      }
    }
    /* THE NEXT IS NOT TRUE IF ACTIVE */
    if (iret != 2) {
      if (idif != 1) {
	continue_ = quit;
	break;
      }
      idif = 2;
      isw = 1;
#if !defined(DOS)
    FunRecover9:
      funCallPath = fCP_gem_iter_gcen2_ ;
#endif
      gcentral();
#if !defined(DOS)
      funCallPath = fCP_gem_iter_ ;
#endif
      fprintf(outfile, "***** SWITCHING TO CENTRAL DIFFERERENCE*****\n");
      continue;
    }
    if (active) {
      continue_ = restart;
      break;
    }
    if (idif == 1)
    {
#if !defined(DOS)
    FunRecover7:
      funCallPath = fCP_gem_iter_gfor ;
#endif
      gforward();
#if !defined(DOS)
      funCallPath = fCP_gem_iter_ ;
#endif
    }
    else
    {
#if !defined(DOS)
    FunRecover6:
      funCallPath = fCP_gem_iter_gcen1_ ;
#endif
      gcentral();
#if !defined(DOS)
      funCallPath = fCP_gem_iter_ ;
#endif
    }
    if (isw != 0)
      continue;
    getytp();
    if (ytp <= 0.0)
      continue;
    if (n == 1)
      d[0] = -(y[0] * d[0] / (t * gsave[0]));
    else
      update();
    if (idg == 3)
      continue_ = quit;
  }
}


/*START GEMINI*/
static Void gemini()
{
#if !defined(DOS)
  /* Shriram: begin */
  int  ckptFunCallPath = ckptTuple . ckptAttribute ;
  
  if ( normalRun == checkpointStatus )
    funCallPath = fCP_gem_ ;
  /* Shriram: end */
#endif

  continue_ = go;
  while (continue_ != quit) {
#if !defined(DOS)
    /* Shriram: begin */
    if ( ( checkpointedRun == checkpointStatus ) &&
         ( iterationLocation == ckptTuple . ckptLocation ) )
      recoverCheckpoint ( NULL ) ;
    else
    {
      if ( checkpointedRun == checkpointStatus )
      {
        /* must be a functionLocation checkpoint, but which? */
        if ( fCP_gem_init == ckptFunCallPath )
          goto FunRecover2 ;
        else if ( ( fCP_gem_iter_st_first == ckptFunCallPath ) ||
                  ( fCP_gem_iter_st_decr == ckptFunCallPath ) ||
                  ( fCP_gem_iter_st_incr == ckptFunCallPath ) ||
                  ( fCP_gem_iter_gcen1 == ckptFunCallPath ) ||
                  ( fCP_gem_iter_gcen2 == ckptFunCallPath ) ||
                  ( fCP_gem_iter_gfor == ckptFunCallPath ) )
          goto FunRecover34567 ;
        else if ( fCP_gem_gfor == ckptFunCallPath )
          goto FunRecover8 ;
      }
    FunRecover2:
#endif
      initialize();
#if !defined(DOS)
      funCallPath = fCP_gem_ ;
    FunRecover8:
#endif
      isw = 0;
      gforward();
#if !defined(DOS)
      funCallPath = fCP_gem_ ;
    }
#endif
    /* Shriram: end */
    continue_ = go;
#if !defined(DOS)
  FunRecover34567:
#endif
    iterate();
#if !defined(DOS)
    funCallPath = fCP_gem_ ;
#endif
  }
  if (ivar == 1)
    inib();
}


static Void initilink()
{
  printf("Program ILINK version%6.2f (1-Feb-1991)\n\n", version);
  printf("FASTLINK version%6.2f (21-Mar-1994)\n\n", fastversion);
  printf("The program constants are set to the following maxima:\n");
  printf("%6ld loci in mapping problem\n", (int)maxlocus);
  printf("%6ld alleles at a single locus\n", (int)maxall);
/*  printf("%6ld recombination probabilities (maxneed)\n", (int)maxneed);*/
  printf("%6ld maximum of censoring array (maxcensor)\n", maxcensor);
  printf("%6ld haplotypes = n1 x n2 x ... where ni = current # alleles at locus i\n",
	 (int)maxhap);
  printf("%6ld joint genotypes for a female\n", (int)maxfem);
  printf("%6ld joint genotypes for a male\n", (int)maxmal);
  printf("%6ld individuals in all pedigrees combined\n", (int)maxind);
  printf("%6ld pedigrees (maxped\n", (int)maxped);
  printf("%6ld quantitative factor(s) at a single locus\n", (int)maxtrait);
  printf("%6ld liability classes\n", (int)maxliab);
  printf("%6ld binary codes at a single locus\n", (int)maxfact);
  printf("%8.2f base scaling factor for likelihood (scale)\n", scale);
  printf("%8.2f scale multiplier for each locus (scalemult)\n", scalemult);
  printf("%8.5f frequency for elimination of heterozygotes (minfreq)\n",
	 minfreq);
  if (minfreq != 0.0) {
    printf("IMPORTANT : RECOMPILE THIS PROGRAM WITH MINFREQ=0.0\n");
    printf("FOR THE ANALYSIS OF RECESSIVE TRAITS\n");
  }
  putchar('\n');
}


main(argc, argv)
int argc;
char *argv[];
{
  /* Shriram: begin */
  int           statReturnCode ;
  struct stat   statBuffer ;
  char          dateTimeStamp [ DateTimeStampStringLength ] ;
  FILE          * scriptCheckpoint ;
  int           scriptRun , scriptToSleep ;
  FILE          * fileTester ;
  /* Shriram: end */

#if defined(GENERR)
  /* Meg: Begin*/
  int ii,jj,kk,ll,mm,indx[maxlocus*maxind],tempind,templocus;
  double savef, orderlhdiff[maxlocus*maxind],numknown; 
  FILE *errout;

  NewPene = false;

  numknown=0.0;
  for (jj=0; jj<maxlocus; jj++)
    for (ii=0; ii<maxind; ii++)
      personunk[jj][ii]=1.0;

  /* Meg: End*/
#endif


#ifdef vms  
  final = fopen("final.dat", "w", "ctx=rec","shr=get,put,upd" ) ;     
#else
  final = fopen("final.dat", "w");      /* If you change things in        */
#endif
  if (final == NULL)                    /* connexion with either final    */
    exit(FileNotFound);               /* or stream, be sure to make     */
  if (dostream) {                       /* the appropriate modifications  */
#ifdef vms
    stream = fopen("stream.dat", "w", "ctx=rec","shr=get,put,upd" ) ;  
#else
    stream = fopen("stream.dat", "w");  /* to the checkpointing code near */
#endif
    if (stream == NULL)                 /* the beginning of outf().       */
      exit(FileNotFound);             /*                      --Shriram */
  }
  else
    stream = NULL;
  
  /* Shriram: begin */

  /* We perform the script-level checkpointing work here because we want
     the files final and stream, which will be used by external scripts,
     to have been created at this point. */

#if !defined(DOS)
#ifdef vms
  scriptCheckpoint = fopen ( ScriptILINKCheckpointFilename , "r", "ctx=rec","shr=get,put,upd" ) ;
#else
  scriptCheckpoint = fopen ( ScriptILINKCheckpointFilename , "r" ) ;
#endif
  if ( NULL != scriptCheckpoint )
  {
    fscanf ( scriptCheckpoint , "%d %d" , & scriptRun , & scriptToSleep ) ;
    fclose ( scriptCheckpoint ) ;
    if ( 0 != scriptToSleep )
    {
#ifdef vms
      scriptCheckpoint = fopen ( ScriptILINKCheckpointFilename , "w", "ctx=rec","shr=get,put,upd" ) ;
#else
      scriptCheckpoint = fopen ( ScriptILINKCheckpointFilename , "w" ) ;
#endif
      scriptRun ++ ;  scriptToSleep -- ;
      fprintf ( scriptCheckpoint , "%d %d\n" , scriptRun , scriptToSleep ) ;
      fclose ( scriptCheckpoint ) ;
      printf ( "Recovering to checkpoint: %d more run(s)\n" ,
               scriptToSleep + 1 ) ;
      if ( NULL != final )
        fclose ( final ) ;
      if ( NULL != stream )
        fclose ( stream ) ;

      if ( 0 == scriptToSleep )
      {
#ifdef vms
        fileTester = fopen ( ScriptILINKFinalOut , "r", "ctx=rec","shr=get,put,upd" ) ;
#else
        fileTester = fopen ( ScriptILINKFinalOut , "r" ) ;
#endif
        if ( NULL != fileTester )
        {
          fclose ( fileTester ) ;
	  copyFile ( ScriptILINKFinalOut , "final.out" ) ;
        }
        if ( dostream )
        {
#ifdef vms
          fileTester = fopen ( ScriptILINKStreamOut , "r", "ctx=rec","shr=get,put,upd" ) ;
#else
          fileTester = fopen ( ScriptILINKStreamOut , "r" ) ;
#endif
          if ( NULL != fileTester )
          {
            fclose ( fileTester ) ;
	    copyFile ( ScriptILINKStreamOut , "stream.out" ) ;
          }
        }
      }
      exit ( EXIT_FAILURE ) ;
    }
  }
#endif

  /* Shriram: end */

  datafile = fopen("datafile.dat", "r");
  if (datafile == NULL)
    exit(FileNotFound);
  ipedfile = fopen("ipedfile.dat", "r");
  if (ipedfile == NULL)
    exit(FileNotFound);
  speedfile = fopen("speedfile.dat", "r");
  if (speedfile == NULL)
    exit(FileNotFound);
#ifdef vms
  outfile = fopen("outfile.dat", "w", "ctx=rec","shr=get,put,upd" ) ;
#else
  outfile = fopen("outfile.dat", "w");
#endif
  if (outfile == NULL)
    exit(FileNotFound);

  initilink();
  ihess = 0;
  ibnd = 1;
  icall = 1;
  ivar = 0;
  ihx = 1;
  inconsistent = false;
  inputdata();
  
#if defined(GENERR)
  /* Meg */
  if (mutsys==1)
    {
      printf("\nError:  Mutation data not compatible with error checking option. \nThe program will exit politely to allow you to correct the problem.\n\n");
      exit(EXIT_FAILURE);
    }
  if (sexlink==1) 
    {
      printf("\nError:  Sexlinked data not compatible with error checking option. \nThe program will exit politely to allow you to correct the problem.\n\n");
      exit(EXIT_FAILURE);
    }
 if (sexdif==1) 
    {
      printf("\nError:  Different male and female recombination fractions are not compatible with error checking option.\nThe program will exit politely to allow you to correct the problem.\n\n");
      exit(EXIT_FAILURE);
    }
#endif

  if (datafile != NULL)
    fclose(datafile);
  datafile = NULL;
  if (ipedfile != NULL)
    fclose(ipedfile);
  ipedfile = NULL;
  if (speedfile != NULL)
    fclose(speedfile);
  speedfile = NULL;
#if !defined(DOS)  
  /* Shriram: begin */
  
  /* We will use stat() to determine the size of the checkpointing
     file.  If stat() reports any of the following: the file cannot
     be accessed; there is a problem with the filename or with the
     buffer; the file does not exist; or the file is shorter than the
     length of the date/time stamp we put in it, then we assume that
     no checkpointing has taken place.
     If we don't ensure the file has contents now, we would be doing
     this later, and might need to return here, etc. */

  statBuffer.st_size = 0 ;
  statReturnCode = stat ( CheckpointILINKFilename , & statBuffer ) ;
/*  if ( ( statReturnCode == EACCES ) || ( statReturnCode == EFAULT ) 
      || ( statReturnCode == ENOENT )
      || ( statBuffer.st_size < DateTimeStampStringLength ) )*/

  if ( statReturnCode != STAT_SUCCESS )
    checkpointStatus = normalRun ;
  else
  {

#ifdef vms
    checkpointDatafile = fopen ( CheckpointILINKFilename , "r", "ctx=rec","shr=get,put,upd" ) ;
#else
    checkpointDatafile = fopen ( CheckpointILINKFilename , "r" ) ;
#endif
    checkpointStatus = checkpointedRun ;
    puts ( "NOTE: attempting to continue previous (unfinished) run" ) ;
    fgets ( dateTimeStamp , DateTimeStampStringLength , checkpointDatafile ) ;
    printf ( "      from %s" , dateTimeStamp ) ;
    getCkptTuple ( ) ;
  }
#endif

  /* Shriram: end */

#if !defined(DOS)
  if ( checkpointStatus != checkpointedRun )	/* Shriram */
  {
#endif
    h = sqrt(exp(-nbit * log(2.0)));
    tol = tolconst;
    /* tol:=tolconst*sqrt(n);*/
    trupb = sqrt(h);
    maxit = n * 50;
    firsttime = true;
#if !defined(DOS)
  }
  firsttime = ( checkpointStatus == normalRun ) ;	/* K. Shriram */
#endif
  lasttime = false;
  dolod = false;
  censorstruct = (censorrec *)Malloc(sizeof(censorrec));
  if (censorstruct == NULL)
    malloc_err("censorstruct");
  gennustruct = (gennurec *)Malloc(sizeof(gennurec));
  if (gennustruct == NULL)
    malloc_err("gennustruct");
  if (approximate) {
    approxstruct = (approxrec *)Malloc(sizeof(approxrec));
    if (approxstruct == NULL)
      malloc_err("approxstruct");
  }
  firstapprox = true;
  if (DIAGNOSTIC)   /* A. A. Schaffer*/
    check_constants();
  getlocations();
  getgeneindices(); /* R. M. Idury */
#if !defined(DOS)
  /* Shriram: begin */
  if ( checkpointedRun == checkpointStatus )
  {
   
#if defined(GENERR)
    /* Meg */
    if ( ( functionLocation == ckptTuple . ckptLocation ) && ( fCP_error == ckptTuple . ckptAttribute ))
      goto FunRecoverErr ;
    /* Meg */
#endif
    if ( ( functionLocation == ckptTuple . ckptLocation ) &&
        ( fCP_outf == ckptTuple . ckptAttribute ))
      goto FunRecover1 ;
  }
  /* Shriram: end */
#endif
  gemini();
 FunRecover1:
  outf();

  if (final != NULL)
    fclose(final);
  if (dostream) {
    if (stream != NULL)
      fclose(stream);
  }

#if defined(GENERR)
  /* Meg: Begin */
  savef = f;
   
#if !defined(DOS)
  if ( checkpointStatus == normalRun )
    {
#endif
    
      for (ii=0; ii<maxlocus; ii++) 
	for (jj=0; jj<maxind; jj++)
	  lhdiff[ii][jj] = savef;
        
     
#if !defined(DOS)
    }

  FunRecoverErr :

  funCallPath = fCP_error;
#endif
 
  ErrorRate = rate;

  for (ii=0; ii<maxlocus; ii++) 
    for (jj=0; jj<=maxall; jj++) 
      for (kk=0; kk<maxall; kk++) 	
	for (ll=0; ll<3; ll++) 	 
	  for (mm=0;mm<maxliab; mm++) 	 
	    altlocuspen[ii][jj][kk][ll][mm] = 0.0;
 
  putlocuspen(ErrorRate);
  NewPene = true;


 #if !defined(DOS)
  
  if ( checkpointedRun == checkpointStatus )
    jplace = ckptTuple . jplace ;
  else
#endif
    jplace = 0 ;
    
    
  for (/*do nothing */; jplace < mlocus; jplace++) 
    {
#if !defined(DOS)
      if ( checkpointedRun == checkpointStatus )
	iplace = ckptTuple . iplace ;
      else
#endif
      iplace = 0;
	for (/*do nothing */; iplace < totperson; iplace++) 
	{
#if !defined(DOS)
	  ckptTuple . jplace  = jplace ;
	  ckptTuple . iplace  = iplace ;

	  if ( normalRun == checkpointStatus )
	    {
	     /* funCallPath = fCP_error ;*/
	      performCheckpoint ( functionLocation , funCallPath , 0 ) ;
	    }
	  else
	    {
	      if ( ( functionLocation == ckptTuple . ckptLocation )
		  && ( fCP_error == ckptTuple . ckptAttribute ) )
		{
		  recoverCheckpoint ( NULL ) ;
		  funCallPath = ckptTuple . ckptAttribute ;
		}
	    }
	 
#endif
	  if (numknown==0.0)
	    {
	      for (jj=0; jj<mlocus; jj++)
		for (ii=0; ii<totperson; ii++)
		  numknown += personunk[jj][ii];

	      fflush(stdout);
	      printf("\n\nGenoCheck will test %4.0lf genotypes for errors.\n\n",numknown);
	    }
 
	  if (personunk[jplace][iplace]==1.0)
	    {
	      fun(&savef, x);
	      lhdiff[jplace][iplace] -= savef;
	      nrep++;
	      if (nrep % 10 == 0)
		putchar('*');
	      else
		putchar('.');
	      fflush(stdout);
	    }
	  else
	    lhdiff[jplace][iplace] =0.0;
	}
    }
  printf("\n\n");

  errout= fopen("PosError","w"); 
  
  if (prettyout == 0)
    {
      for (ii=0; ii<mlocus; ii++) 
	{
	  for (jj=0; jj<totperson; jj++)
	    fprintf(errout,"%12.9lf ",lhdiff[ii][jj]);
	  
	  fprintf(errout,"\n");
	}
      fprintf(errout,"\n");
    }
  else  /* ordered output */
    {

      fprintf(errout,"\nGenoCheck, version 1.0 \n\n");
      fprintf(errout,"The error index for each genotype is separated by pedigree\n");
      fprintf(errout,"and then by locus within each pedigree.  Large error index\n");
      fprintf(errout,"indicates likely to contain error.\n");

      for (ll=1; ll<=nuped; ll++)
	{
	  fprintf(errout,"\n\n#####################################################\nPedigree %d:\n#####################################################\n",person[pedind[ll]-1]->ped);
	  for (ii=0; ii<mlocus; ii++)                             
	    { 
	      totalgenotypes=pedind[ll]-pedind[ll-1];
	      kk=0;
	      for (jj=pedind[ll-1]; jj<pedind[ll]; jj++)
		{
		  if (personunk[ii][jj]==0.0)
		    {
		      lhdiff[ii][jj]=-100.00;
		      totalgenotypes -=1;
		    }
		
		  orderlhdiff[kk++]=lhdiff[ii][jj];
		}
	      indexx(pedind[ll]-pedind[ll-1],orderlhdiff-1,indx-1);
	      fprintf(errout,"\n\nList of individuals for locus %d ordered by decreasing\n",ii+1);
	      fprintf(errout,"error index.\n\n");
	      fprintf(errout,"Individual     Error Index\n\n");  

	      for (kk=(pedind[ll]-pedind[ll-1]-1); kk>=0; kk--)
		{
		  tempind = (indx[kk]-1)-(((indx[kk]-1)/(pedind[ll]-pedind[ll-1]))*(pedind[ll]-pedind[ll-1]))+1;
		
		  if (personunk[ii][pedind[ll-1]+tempind-1]>0)
		    fprintf(errout," %3d           %9.6lf\n",tempind, lhdiff[ii][pedind[ll-1]+tempind-1]);
		}
	    }
	}
    }
/*

   totalgenotypes = mlocus*totperson;
   kk=0;
   for (ii=0; ii<mlocus; ii++) 
   for (jj=0; jj<totperson; jj++)
   {
   if (personunk[ii][jj]==0.0)
   {
   lhdiff[ii][jj]=-100.00;
   totalgenotypes -=1;
   }
   orderlhdiff[kk++]=lhdiff[ii][jj];
   }

   indexx(totperson*mlocus,orderlhdiff-1,indx-1);
      
   fprintf(errout,"\nList of individuals and loci ordered by decreasing\n");
   fprintf(errout,"error index.  Large index indicates likely to contain error.\n\n");
   fprintf(errout,"Locus     Individual     Error Index\n\n");  
    
   for (kk=(totperson*mlocus-1); kk>=0; kk--)
   {
   tempind = (indx[kk]-1)-(((indx[kk]-1)/totperson)*totperson)+1;
   templocus = (indx[kk]-1)/totperson+1;

   if (personunk[templocus-1][tempind-1]>0)
   fprintf(errout,"%3d         %3d           %9.6lf\n",templocus, tempind, lhdiff[templocus-1][tempind-1]);
   }
   }
   fprintf(errout,"\nTotal number of genotypes analyzed:  %d\n",totalgenotypes);
*/
  fclose(errout);
  /* Meg: End */
#endif

  /* Shriram: begin */
  
  /* This is the time at which to update the script-level checkpointing
     routine, since we have closed final and stream; any disasters that
     occur after this are "safe". */

#if !defined(DOS)
#ifdef vms
  scriptCheckpoint = fopen ( ScriptILINKCheckpointFilename , "r+", "ctx=rec","shr=get,put,upd" ) ;
#else
  scriptCheckpoint = fopen ( ScriptILINKCheckpointFilename , "r+" ) ;
#endif
  if ( NULL != scriptCheckpoint )
  {
    rewind ( scriptCheckpoint ) ;
    scriptRun ++ ;
    fprintf ( scriptCheckpoint , "%d 0\n" , scriptRun ) ;
    fclose ( scriptCheckpoint ) ;

    copyFile ( "final.out" , ScriptILINKFinalOut ) ;
    appendFile ( "final.dat" , ScriptILINKFinalOut ) ;
    
    if ( dostream )
    {
      copyFile ( "stream.out" , ScriptILINKStreamOut ) ; 
      appendFile ( "stream.dat" , ScriptILINKStreamOut ) ; 
    }
  }
#endif

#if !defined(DOS)  
  if ( NULL != checkpointDatafile )
    fclose ( checkpointDatafile ) ;
  
  unlink ( CheckpointILINKFileBackup ) ;
  unlink ( CheckpointILINKFilename ) ;

  unlink ( OutfILINKFinalDat ) ;
  unlink ( OutfILINKStreamDat ) ;
#endif

  /* Shriram: end */
  
  if (outfile != NULL)
    fclose(outfile);
  if (ipedfile != NULL)
    fclose(ipedfile);
  if (datafile != NULL)
    fclose(datafile);
  if (speedfile != NULL)
    fclose(speedfile);
  exit(EXIT_SUCCESS);
}

/* End. */