/* This file contains a modified version of the LINKMAP program */
/* The modifications are described in the papers: */
/* R. W. Cottignham, 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"
#include "gemdefs.h"
#include "lidefs.h"
#if !defined(DOS)
#include "checkpointdefs.h"
#endif
#ifndef LESSMEMORY
#include "moddefs.h"
#endif


#if PARALLEL  /* cgh */
#include "compar.h"          /* parallel support code */
#endif  /* if PARALLEL -- cgh */

extern void checkzero();
#if !defined(KNR_PROTO)
extern boolean zerotest(double);
#else
extern boolean zerotest();
#endif  /* !defined(KNR_PROTO) */


#if PARALLEL  /* cgh */

/* Calculate the number of calls to iterpeds() by simulating the loop
   in main where it is called.  We utilize temporary copies of
   volatile global variables to ensure nothing external is affected. */
void simIpedLoop(status)
loopStatus status;  
{
  int i, j;
  int currentIter;  /* current number of simulated iteration */

  strBuff* strmBuff;             /* buffer for stream output */
  char tempBuff[TEMPBUFF_SIZE];  /* temporary buffer */
  
  /* initialize temporaries */
  thetarray tMaletheta;
  thetarray tFemaletheta;
  boolean tContinue_ = continue_;
  double tThisdist, tThisdistm, tThisdistf;
  tThisdist = thisdist;
  tThisdistm = thisdistm;
  tThisdistf = thisdistf;

  for (i = 0; i < mlocus - 1; i++)
    tMaletheta[i] = maletheta->theta[i];
  if (readfemale)
    for (i = 0; i < mlocus - 1; i++)
      tFemaletheta[i] = femaletheta->theta[i];

  currentIter = 0;
  if (status == countIpeds)
    /* initialize zero theta counter */
    numZeroThetas = 0;
  
  /* simulate loop in main that calls iterpeds */
  do {
    if (status == computeThetas) {
      /* save the computed theta values */
      gZeroTheta[currentIter] = false;
      for (i = 0; i < mlocus - 1; i++) {
	gmaletheta[currentIter][i] = tMaletheta[i];
	if (zerotest(tMaletheta[i]))
	  /* if this thetarray has a 0, record it */
	  gZeroTheta[currentIter] = true;
      }
      if (readfemale)
	for (i = 0; i < mlocus - 1; i++)
	  gfemaletheta[currentIter][i] = tFemaletheta[i];
      /* buffer output as well */
      /* stream buffer entry for this loop iteration */
      strmBuff = streamBuff->ipeds[currentIter];

      if (dostream) {
	if (whichvary == 1)
	  tThisdist = tThisdistm - getdist(tMaletheta);
	else 
	  tThisdist = tThisdistm + getdist(&tMaletheta[whichvary - 2]);
	if (normalRun == checkpointStatus) {
	  sprintf(tempBuff, "% .5e\n", tThisdist);
	  append(strmBuff, tempBuff);
	}
	if (sexdif) {
	  if (whichvary == 1)
	    tThisdist = tThisdistf - getdist(tFemaletheta);
	  else
	    tThisdist = tThisdistf + getdist(&tFemaletheta[whichvary - 2]);
	  if (normalRun == checkpointStatus) {
	    sprintf(tempBuff, "% .5e\n", tThisdist);
	    append(strmBuff, tempBuff);
	  }
	}
	for (i = 1; i < mlocus; i++)
	  if (normalRun == checkpointStatus) {
	    sprintf(tempBuff, "% .5e\n", tMaletheta[i - 1]);
	    append(strmBuff, tempBuff);
	  }
	if (normalRun == checkpointStatus)
	  append(strmBuff, "\n");
	if (sexdif) {
	  for (i = 1; i < mlocus; i++)
	    if (normalRun == checkpointStatus) {
	      sprintf(tempBuff, "% .5e\n", tFemaletheta[i - 1]);
	      append(strmBuff, tempBuff);
	    }
	}
	if (normalRun == checkpointStatus)
	  append(strmBuff, "\n");
      }  /* if dostream */
    } else {
      /* we are counting iterations -- but we also want to count
	 the number of zero thetas */
      for (i = 0; i < mlocus - 1; i++) {
	if (zerotest(tMaletheta[i])) {
	  numZeroThetas++;
	  break;
	}
      }
    }
    
    /* this is where iterpeds() is called */
    currentIter++;
    
    if (whichvary == 1) {
      tMaletheta[0] -= thetainc;
      if (zerotest(tMaletheta[0]))
	tMaletheta[0] = 0.0;
      /* don't care about this when counting iterations */
      if (status == computeThetas)
	if (readfemale) {
	  tFemaletheta[0] = distratio * getdist(tMaletheta);
	  tFemaletheta[0] = invdist(tFemaletheta);
	  if (zerotest(tFemaletheta[0]))
	    tFemaletheta[0] = 0.0;
	}
    } else if (whichvary == mlocus) {
      tMaletheta[mlocus - 2] += thetainc;
      if (zerotest(tMaletheta[mlocus - 2]))
	tMaletheta[mlocus - 2] = 0.0;
      /* don't care about this when counting iterations */
      if (status == computeThetas)
	if (readfemale) {
	  tFemaletheta[mlocus - 2] = distratio * getdist(&tMaletheta[mlocus - 2]);
	  tFemaletheta[mlocus - 2] = invdist(&tFemaletheta[mlocus - 2]);
	  if (zerotest(tFemaletheta[mlocus - 2]))
	    tFemaletheta[mlocus - 2] = 0.0;
	}
    } else {
      tMaletheta[whichvary - 2] += thetainc;
      if (zerotest(tMaletheta[whichvary - 2]))
	tMaletheta[whichvary - 2] = 0.0;
      tMaletheta[whichvary - 1] = (thetatot - tMaletheta[whichvary - 2]) /
	(1.0 - 2.0 * tMaletheta[whichvary - 2]);
      if (zerotest(tMaletheta[whichvary - 1]))
	tMaletheta[whichvary - 1] = 0.0;
      /* don't care about this when counting iterations */
      if (status == computeThetas)
	if (readfemale) {
	  tFemaletheta[whichvary - 2] =
	    distratio * getdist(&tMaletheta[whichvary - 2]);
	  tFemaletheta[whichvary - 2] =
	    invdist(&tFemaletheta[whichvary - 2]);
	  if (zerotest(tFemaletheta[whichvary - 2]))
	    tFemaletheta[whichvary - 2] = 0.0;
	  tFemaletheta[whichvary - 1] =
	    distratio * getdist(&tMaletheta[whichvary - 1]);
	  tFemaletheta[whichvary - 1] =
	    invdist(&tFemaletheta[whichvary - 1]);
	  if (zerotest(tFemaletheta[whichvary - 1]))
	    tFemaletheta[whichvary - 1] = 0.0;
	}
    }
    if (whichvary == 1)
      tContinue_ = (tMaletheta[0] >= finaltheta && tMaletheta[0] >= 0.0);
    else {
      if (whichvary == mlocus)
	tContinue_ = (tMaletheta[whichvary - 2] <= finaltheta &&
		      tMaletheta[whichvary - 2] <= 0.5);
      else
	tContinue_ = (tMaletheta[whichvary - 2] <= finaltheta &&
		      tMaletheta[whichvary - 1] >= 0.0);
    }
  } while (tContinue_);
  /* if it was our job to count the calls to iterpeds(), set the
     global numIpeds variable */
  if (status == countIpeds) {
    numIpeds = currentIter;
  } else {
    /* calculate the number of nonzero theta vectors */
    for (i = 0, numNonzeroThetas = 0; i < numIpeds; i++)
      if (gZeroTheta[i] == false)
	/* increment count of nonzero thetas, and save ordering
	   of nozero thetas */
	absoluteThetaIndex[numNonzeroThetas++] = i;
    /* save ordering of theta vectors that contain zeros */
    for (i = 0, j = numNonzeroThetas; i < numIpeds; i++)
      if (gZeroTheta[i] == true)
	absoluteThetaIndex[j++] = i;
  }
  return;
}

#else  /* if PARALLEL -- cgh */

/* The following routine iterates over the different pedigrees and
   handles output. */
void iterpeds()
{
  int i;
  short thisped;
  int locus_loop;  /* dwix */

  /* dwix: begin */
  /* I have improved the performace of this function by letting it */
  /* save it's results on runs which will be recomputed next */
  /* iteration of the linkmap program in the script.  The */
  /* recomputations will not be performed, but rather the saved data */
  /* is used.  I use the files s_screen, s_outfile, and s_stream to */
  /* save the output to the screen, outfile, and stream respectively. */
  /* I have mirrored every output to the screen with an output to */
  /* s_screen, every output to outfile to s_outfile, and every output */
  /* to stream to s_stream.   On the restore, I put the contents of */
  /* each file back to its respective place. */
  int srun; /* saved run */
  int rrun; /* restored run */

#if !defined(vms) && !defined(DOS) && (!PARALLEL) && (MULTI_LINKMAP)
  int    ppid;       /* the parent process ID of this process */
  int    omlocus;    /* the mlocus of the run LINKSAVE was from */
  int    count;      /* a counter variable for loops */
  double otheta;     /* the thetas of the other run */

  char   instr[DEFAULT_STR_LENGTH]; /* str to read files with */
#endif /* !defined(vms) && !defined(DOS) && !PARALLEL && MULTI_LINKMAP */

  FILE   *s_screen;  /* file to save screen output */
  FILE   *s_outfile; /* file to save outfile output */
  FILE   *s_stream;  /* file to save stream output */
  FILE   *s_save;    /* file with information regarding to save */

  locus_loop = mlocus - 2; /* common number for loop iterations */
  srun = false; /* is this a saved run? */
  rrun = false; /* is this a restored run */

  s_screen = NULL;  /* init all FILE* to NULL */
  s_outfile = NULL;
  s_stream = NULL;
  s_save = NULL;
  /* dwix: end */

  tlike = 0.0;
  alike = 0.0;
  for (i = 1; i <= totperson; i++)
    person[i]->done = false;
  for (i = 1; i <= totperson; i++)
    person[i]->firstpass = true;
  thisc = minint;
  recombination();
  checkzero();

#if !defined(vms) && !defined(DOS) && (!PARALLEL) && (MULTI_LINKMAP)
  /* dwix: begin */
  /* check if it is a run we need to save */
  if (maletheta->theta[whichvary-1]==0.0) srun = true;
  /* if it's a run we also need to restore, then we should not save */
  if ((whichvary > 1) && (maletheta->theta[whichvary-2]==0.0))
    srun = false;
  
  if (srun) {
    /* if run to save */
    unlink(LINKSAVE);    /* delete all the old save files */
    unlink(SCREENSAVE);
    unlink(OUTFILESAVE);
    unlink(STREAMSAVE);
    s_screen = fopen(SCREENSAVE,"w");  /* open the files for output */
    s_outfile = fopen(OUTFILESAVE,"w");
    s_stream = fopen(STREAMSAVE,"w");
    if ((s_screen==NULL) || (s_outfile==NULL) || (s_stream==NULL)) {
      /* dwix err msg */
      fprintf(stderr,"ERROR opening the save files.\n");
      srun = false;
      fclose(s_screen);
      fclose(s_outfile);
      fclose(s_stream);
    }
  }

  /* check if it is a run to restore */
  if ((whichvary > 1) && (maletheta->theta[whichvary-2] == 0.0))
    rrun = true;
  if (rrun) {
    /* if run to restore */
    s_save = fopen(LINKSAVE,"r");   /* try to open the head save file */
    if (s_save==NULL)
      rrun = false; /* if error no restore */
  }
  if (rrun) {
    /* still to restore */
    fscanf(s_save,"%d",&ppid);            /* get ppid from the file */
    /* check it against ours */
    if (ppid != getppid())
      rrun = false;
    fscanf(s_save,"%d",&omlocus);         /* get mlocus from other run */
    if (mlocus != omlocus)
      /* check it against ours */
      rrun = false;
    if (rrun) {
      /* still valid for restore -- do loop for loci */
      for (count=1;count<mlocus;count++) {
	/* read theta from other run */
	if (fscanf(s_save,"%lf",&otheta) != 1) {
	  /* scan conversion failed */
	  fprintf(stderr,"Warning reading LINKSAVE... Do not panic..\n");
	  rrun = false;  /* cgh */
	}
	/* cgh -- equality test didn't work due to rounding errors */
	if (fabs(otheta - maletheta->theta[count]) > epsilon)
	  rrun = false;
      }
    }
    fclose(s_save);                     /* and close the file */
  }
  if (rrun) {
     /* if still to restore */
    s_screen = fopen(SCREENSAVE,"r");   /* open the files for restore */
    s_outfile = fopen(OUTFILESAVE,"r");
    s_stream = fopen(STREAMSAVE,"r");
    if ((s_screen==NULL) || (s_outfile==NULL) || (s_stream==NULL)) {
      /* dwix err msg */
      fprintf(stderr,"ERROR opening the save files.\n");
      rrun = false;
      fclose(s_screen);
      fclose(s_outfile);
      fclose(s_stream);
    }
  }
  if (rrun) {
    /* if still to restore */
    /* output the screen file's data to screen */
    fgets(instr,DEFAULT_STR_LENGTH,s_screen);
    while (!feof(s_screen)) {
      printf("%s",instr);
      fgets(instr,DEFAULT_STR_LENGTH,s_screen);
    }
    
    /* output the outfile file's data to outfile */
    fgets(instr,DEFAULT_STR_LENGTH,s_outfile);
    while (!feof(s_outfile)) {
      fprintf(outfile,"%s",instr);
      fgets(instr,DEFAULT_STR_LENGTH,s_outfile);
    }

    /* output the stream file's data to stream */
    fgets(instr,DEFAULT_STR_LENGTH,s_stream);
    while (!feof(s_stream)) {
      fprintf(stream,"%s",instr);
      fgets(instr,DEFAULT_STR_LENGTH,s_stream);
    }

    /* close and delete all the files */
    fclose(s_screen);
    fclose(s_outfile);
    fclose(s_stream);
    unlink(LINKSAVE);
    unlink(SCREENSAVE);
    unlink(OUTFILESAVE);
    unlink(STREAMSAVE);
    return;
  }
  /* dwix: end */
#endif /* !defined(vms) && !defined(DOS) && !PARALLEL && MULTI_LINKMAP */
  
  fprintf(outfile, "%s\n%s\n", LINE, LINE);
  if (srun) 
    fprintf(s_outfile, "%s\n%s\n", LINE, LINE);

  if (sexdif) {
    fprintf(outfile, "MALE THETAS   ");
    if (srun)
      fprintf(s_outfile, "MALE THETAS   ");
  } else {
    fprintf(outfile, "THETAS ");
    if (srun)
      fprintf(s_outfile, "THETAS ");
  }
  
  for (i = 0; i <= locus_loop; i++) {
    fprintf(outfile, "%6.3f", maletheta->theta[i]);
    if (srun)
      fprintf(s_outfile, "%6.3f", maletheta->theta[i]);
  }
  putc('\n', outfile);
  if (srun)
    putc('\n', s_outfile);

  if (sexdif) {
    fprintf(outfile, "FEMALE THETAS ");
    if (srun)
      fprintf(s_outfile, "FEMALE THETAS ");
    for (i = 0; i <= locus_loop; i++) {
      fprintf(outfile, "%6.3f", femaletheta->theta[i]);
      if (srun)
	fprintf(s_outfile, "%6.3f", femaletheta->theta[i]);
    }
    putc('\n', outfile);
    if (srun)
      putc('\n', s_outfile);
  }
  
  fprintf(outfile, "%s\n", LINE);
  fprintf(outfile, "PEDIGREE |  LN LIKE  | LOG 10 LIKE\n");
  fprintf(outfile, "%s\n", LINE);
  if (srun) {
    fprintf(s_outfile, "%s\n", LINE);
    fprintf(s_outfile, "PEDIGREE |  LN LIKE  | LOG 10 LIKE\n");
    fprintf(s_outfile, "%s\n", LINE);
  }

  printf("%s\n%s\n", LINE, LINE);
  if (srun)
    fprintf(s_screen, "%s\n%s\n", LINE, LINE);

  if (sexdif) {
    printf("MALE THETAS   ");
    if (srun)
      fprintf(s_screen,"MALE THETAS   ");
  } else {
    printf("THETAS ");
    if (srun)
      fprintf(s_screen,"THETAS ");
  }
  
  for (i = 0; i <= locus_loop; i++) {
    printf("%6.3f", maletheta->theta[i]);
    if (srun)
      fprintf(s_screen,"%6.3f", maletheta->theta[i]);
  }
  if (interfer) {
    printf("%6.3f", maletheta->theta[mlocus - 1]);
    if (srun)
      fprintf(s_screen,"%6.3f", maletheta->theta[mlocus - 1]);
  }
  putchar('\n');
  if (srun)
    putc('\n',s_screen);

  if (sexdif) {
    printf("FEMALE THETAS ");
    if (srun) fprintf(s_screen,"FEMALE THETAS ");
    for (i = 0; i <= locus_loop; i++) {
      printf("%6.3f", femaletheta->theta[i]);
      if (srun)
	fprintf(s_screen,"%6.3f", femaletheta->theta[i]);
    }
    if (interfer) {
      printf("%6.3f", femaletheta->theta[mlocus - 1]);
      if (srun)
	fprintf(s_screen,"%6.3f", femaletheta->theta[mlocus - 1]);
    }
    putchar('\n');
    if (srun)
      putc('\n',s_screen);
  }
  
  printf("%s\n", LINE);
  printf("PEDIGREE |  LN LIKE  | LOG 10 LIKE\n");
  printf("%s\n", LINE);
  if (srun) {
    fprintf(s_screen, "%s\n", LINE);
    fprintf(s_screen,"PEDIGREE |  LN LIKE  | LOG 10 LIKE\n");
    fprintf(s_screen, "%s\n", LINE);
  }

#if LOOPSPEED
  open_loop_file();
#endif  
  for (thisped = 0; thisped < nuped; thisped++) {
#if LOOPSPEED
    read_loop_file(thisped + 1);
#endif
    likelihood((int)(thisped + 1), proband[thisped]);
    if (byfamily) {
      fprintf(outfile, "%9d %12.6f ", proband[thisped]->ped, like);
      if (srun)
	fprintf(s_outfile, "%9d %12.6f ", proband[thisped]->ped,like);
    }
    if (dostream) {
      fprintf(stream, "%12d % .5e ", proband[thisped]->ped, like);
      if (srun)
	fprintf(s_stream, "%12d % .5e ", proband[thisped]->ped, like);
    }
    printf("%9d %12.6f ", proband[thisped]->ped, like);
    if (srun)
      fprintf(s_screen,"%9d %12.6f ", proband[thisped]->ped, like);
    alike += like;
    like /= log10_;
    if (byfamily) {
      fprintf(outfile, "%12.6f\n", like);
      if (srun)
	fprintf(s_outfile, "%12.6f\n", like);
    }
    if (dostream) {
      fprintf(stream, "%12.6f\n", like);
      if (srun)
	fprintf(s_stream, "%12.6f\n", like);
    }
    printf("%12.6f\n", like);
    if (srun)
      fprintf(s_screen,"%12.6f\n", like);
    tlike += like;
  }
#if LOOPSPEED
  close_loop_file();
#endif
  fprintf(outfile, "%s\n", LINE);
  fprintf(outfile, "TOTALS    %12.6f %12.6f\n", alike, tlike);
  if (srun) {
    fprintf(s_outfile, "%s\n", LINE);
    fprintf(s_outfile, "TOTALS    %12.6f %12.6f\n", alike, tlike);
  }

  printf("%s\n", LINE);
  printf("TOTALS    %12.6f %12.6f\n", alike, tlike);
  if (srun) {
    fprintf(s_screen, "%s\n", LINE);
    fprintf(s_screen,"TOTALS    %12.6f %12.6f\n", alike, tlike);
  }
  if (dostream) {
    fprintf(stream, "% .5e % .5e\n", -2 * alike, tlike);
    if (srun)
      fprintf(s_stream, "% .5e % .5e\n", -2 * alike, tlike);
  }
  alike = -2 * alike;
  fprintf(outfile, "-2 LN(LIKE) = % .5e\n", alike);
  if (srun)
    fprintf(s_outfile, "-2 LN(LIKE) = % .5e\n", alike);
  printf("-2 LN(LIKE) = % .5e\n", alike);
  if (srun)
    fprintf(s_screen,"-2 LN(LIKE) = % .5e\n", alike);
  if (firsttime) {
    if (thisc < maxcensor) {
      printf("Maxcensor can be reduced to %12d\n", thisc);
      if (srun)
	fprintf(s_screen,"Maxcensor can be reduced to %12d\n", thisc);
    } else {
      if (thisc > maxcensor) {
	printf("You may gain efficiency by increasing maxcensor\n");
	if (srun)
	  fprintf(s_screen, "You may gain efficiency by increasing maxcensor\n");
      }
    }
  }
  
  /*  exit(EXIT_SUCCESS); */ /* exitramana */
  firsttime = false;

#if !defined(vms) && !defined(DOS) && (!PARALLEL) && (MULTI_LINKMAP)
  if (srun) {
    fclose(s_screen);
    fclose(s_outfile);
    fclose(s_stream);
    s_save = fopen(LINKSAVE,"w");
    if (s_save==NULL) {
      /* dwix err msg */
      fprintf(stderr,"ERROR: could not open LINKSAVE file.\n");
      exit(EXIT_FAILURE);
    }
    fprintf(s_save,"%d",getppid());
    fprintf(s_save," %d",mlocus);
    for (count = 1; count < mlocus; count++)
      fprintf(s_save," %lf",maletheta->theta[count]);
    fclose(s_save);
  }

  /* need to get rid of files because this might be the last theta */
  if ((maletheta->theta[whichvary-1] != 0.0) && (whichvary > 1) &&
      (maletheta->theta[whichvary-2] != 0.0)) {
    unlink(LINKSAVE);
    unlink(SCREENSAVE);
    unlink(OUTFILESAVE);
    unlink(STREAMSAVE);
  }
#endif /* !defined(vms) && !defined(DOS) && !PARALLEL && MULTI_LINKMAP */

}  /* iterpeds */


#endif  /* if PARALLEL -- cgh */


void preIpedLoop()
{
  if (dostream) {
    if (checkpointStatus !=checkpointedRun ) {
      fprintf(stream, "@\n");
      fprintf(stream, "LINKMAP\n");
      fprintf(stream, "%12d %12d %12d ", mlocus, whichvary, gridsize);
      if (sexlink)
	fprintf(stream, " 1 ");
      else
	fprintf(stream, " 0 ");
      if (sexdif && readfemale)
	fprintf(stream, "2 ");
      else {
	if (sexdif)
	  fprintf(stream, "1\n");
	else 
	  fprintf(stream, "0\n");
      }
      fprintf(stream, "%12d\n", nuped);
      for (i = 1; i <= mlocus; i++) {
	j = 0;
	do {
	  j++;
	} while (j != order[i - 1]);
	fprintf(stream, "%12d\n", j);
      }
      putc('\n', stream);
    }
  
    /* This call to recombination() is needed to make sure that the
       proper theta values are printed to stream.  recombination()
       is called again before liklihood() is ever called from within
       iterpeds() (in sequential LINKMAP), or multiPedLike() (in
       parallel LINKMAP).  Therefore, if dostream is false, this
       call can be safely omitted.
       In parallel LINKMAP, we don't want the barrier inside of
       recombination() to get called, and we don't want
       g{male,female}segprob assigned to.  This code in
       recombination() is only executed when the boolean flag infun
       is true, so we set it to false before the call. -- cgh */
#if PARALLEL  /* cgh */
    infun = false;
#endif  /* if PARALLEL -- cgh */
    recombination();

    if (normalRun == checkpointStatus) {
      for (i = 1; i < mlocus; i++)
	fprintf(stream, "% .5e\n", maletheta->theta[i - 1]);
      putc('\n', stream);
      if (sexdif) {
	for (i = 1; i < mlocus; i++)
	  fprintf(stream, "% .5e\n", femaletheta->theta[i - 1]);
      }
      putc('\n', stream);
    }
    
    thisdistm = 0.0;
    if (normalRun == checkpointStatus)
      fprintf(stream, "% .5e\n", thisdistm);
    if (whichvary == 1)
      j = 2;
    else
      j = 1;
    for (i = j; i < mlocus; i++) {
      thisdistm += getdist(&maletheta->theta[i - 1]);
      if ((i + 1 != whichvary) && (normalRun == checkpointStatus))
	fprintf(stream, "% .5e\n", thisdistm);
    }
    if (normalRun == checkpointStatus)
      putc('\n', stream);
    if (sexdif) {
      thisdistf = 0.0;
      fprintf(stream, "% .5e\n", thisdistf);
      if (whichvary == 1)
	j = 2;
      else
	j = 1;
      for (i = j; i < mlocus; i++) {
	thisdistf += getdist(&femaletheta->theta[i - 1]);
	if ((i + 1 != whichvary) && (normalRun == checkpointStatus))
	  fprintf(stream, "% .5e\n", thisdistf);
      }
      if (normalRun == checkpointStatus)
	putc('\n', stream);
    }
    thisdistm = 0.0;
    if (whichvary != 1) {
      for (i = 2; i <= whichvary; i++)
	thisdistm += getdist(&maletheta->theta[i - 2]);
    }
    thisdistf = 0.0;
    if (sexdif) {
      if (whichvary != 1) {
	for (i = 2; i <= whichvary; i++)
	  thisdistf += getdist(&femaletheta->theta[i - 2]);
      }
    }
  }
  
  if (whichvary != mlocus && whichvary != 1) {
    thetatot = maletheta->theta[whichvary - 2] + maletheta->theta[whichvary - 1] -
	       2 * maletheta->theta[whichvary - 2] * maletheta->theta[whichvary - 1];
    thetainc = thetatot / gridsize;
    if (readfemale) {
      distratio = femaletheta->theta[whichvary - 2] + femaletheta->theta[whichvary - 1] -
		  2 * femaletheta->theta[whichvary - 2] * femaletheta->theta[whichvary - 1];
      distratio = getdist(&distratio) / getdist(&maletheta->theta[whichvary - 1]);
    }
  } else {
    thetainc = 0.5 / gridsize;
    if (readfemale) {
      if (whichvary == 1) {
	distratio = femaletheta->theta[1];
	distratio = getdist(&distratio) / getdist(&maletheta->theta[1]);
      } else {
	distratio = femaletheta->theta[mlocus - 3];
	distratio = getdist(&distratio) / getdist(&maletheta->theta[mlocus - 3]);
      }
    }
  }
}


#if !PARALLEL
void ipedLoop()
{
#if !defined(DOS)  
  iterpeds_counter = 0;  /*set up for checkpointing*/
#endif  /* !defined(DOS) */
  do {    
    if (dostream) {
      if (whichvary == 1)
	thisdist = thisdistm - getdist(maletheta->theta);
      else
	thisdist = thisdistm + getdist(&maletheta->theta[whichvary - 2]);
      if (normalRun == checkpointStatus)
	fprintf(stream, "% .5e\n", thisdist);
      if (sexdif) {
	if (whichvary == 1)
	  thisdist = thisdistf - getdist(femaletheta->theta);
	else
	  thisdist = thisdistf + getdist(&femaletheta->theta[whichvary - 2]);
	if (normalRun == checkpointStatus)
	  fprintf(stream, "% .5e\n", thisdist);
      }
      for (i = 1; i < mlocus; i++)
	if (normalRun == checkpointStatus)
	  fprintf(stream, "% .5e\n", maletheta->theta[i - 1]);
      if (normalRun == checkpointStatus)
	putc('\n', stream);
      if (sexdif) {
	for (i = 1; i < mlocus; i++)
	  if (normalRun == checkpointStatus)
	    fprintf(stream, "% .5e\n", femaletheta->theta[i - 1]);
      }
      if (normalRun == checkpointStatus)
	putc('\n', stream);
    }
#if !defined(DOS)
    if ((checkpoint_counter == 0) && 
	(iterped_call_before == checkpoint_place))
	checkpointStatus = normalRun;
#endif  /* if !defined(DOS) */
    if (normalRun == checkpointStatus) {
      ensureWrite(outfile);
      ensureWrite(stream);
#if !defined(DOS)
      performCheckpoint(functionLocation, iterpeds_counter, iterped_call_before);
      copyFile("outfile.dat", MainLINKMAPOutfileDat);
      if (dostream)
	copyFile("stream.dat", MainLINKMAPStreamDat);
#endif  /* if !defined(DOS) */
    }
#if !defined(DOS)
    if ((checkpoint_counter == 0) && 
	(iterped_call_before == checkpoint_place)){
#endif  /* if !defined(DOS) */
      
       iterpeds();
       ensureWrite(outfile);
       ensureWrite(stream);
#if !defined(DOS)
       performCheckpoint(functionLocation, iterpeds_counter, iterped_call_after);
     }
    if ((checkpoint_counter == 0) && (iterped_call_after == checkpoint_place)){
      checkpoint_place = iterped_call_before;
      checkpointStatus = normalRun;
    }
    if (normalRun == checkpointStatus) {
      copyFile("outfile.dat", MainLINKMAPOutfileDat);
      if (dostream)
	copyFile("stream.dat", MainLINKMAPStreamDat);
    }
#endif  /* if !defined(DOS) */

    if (whichvary == 1) {
      maletheta->theta[0] -= thetainc;
      if (zerotest(maletheta->theta[0]))
	maletheta->theta[0] = 0.0;
      if (readfemale) {
	femaletheta->theta[0] = distratio * getdist(maletheta->theta);
	femaletheta->theta[0] = invdist(femaletheta->theta);
        if (zerotest(femaletheta->theta[0]))
	  femaletheta->theta[0] = 0.0;
      }
    } else if (whichvary == mlocus) {
      maletheta->theta[mlocus - 2] += thetainc;
      if (zerotest(maletheta->theta[mlocus - 2]))
	maletheta->theta[mlocus - 2] = 0.0;
      if (readfemale) {
	femaletheta->theta[mlocus - 2] = distratio * getdist(&maletheta->theta[mlocus - 2]);
	femaletheta->theta[mlocus - 2] = invdist(&femaletheta->theta[mlocus - 2]);
	if (zerotest(femaletheta->theta[mlocus - 2]))
	  femaletheta->theta[mlocus - 2] = 0.0;
      }
    } else {
      maletheta->theta[whichvary - 2] += thetainc;
      if (zerotest(maletheta->theta[whichvary - 2]))
	maletheta->theta[whichvary - 2] = 0.0;
      maletheta->theta[whichvary - 1] = (thetatot - maletheta->theta[whichvary - 2]) /
				   (1.0 - 2.0 * maletheta->theta[whichvary - 2]);
      if (zerotest(maletheta->theta[whichvary - 1]))
	maletheta->theta[whichvary - 1] = 0.0;
      if (readfemale) {
	femaletheta->theta[whichvary - 2] =
	  distratio * getdist(&maletheta->theta[whichvary - 2]);
	femaletheta->theta[whichvary - 2] = invdist(&femaletheta->theta[whichvary - 2]);
	if (zerotest(femaletheta->theta[whichvary - 2]))
	  femaletheta->theta[whichvary - 2] = 0.0;
	femaletheta->theta[whichvary - 1] =
	  distratio * getdist(&maletheta->theta[whichvary - 1]);
	femaletheta->theta[whichvary - 1] = invdist(&femaletheta->theta[whichvary - 1]);
	if (zerotest(femaletheta->theta[whichvary - 1]))
	  femaletheta->theta[whichvary - 1] = 0.0;
      }
    }
    if (whichvary == 1)
      continue_ = (maletheta->theta[0] >= finaltheta && maletheta->theta[0] >= 0.0);
    else {
      if (whichvary == mlocus)
	continue_ = (maletheta->theta[whichvary - 2] <= finaltheta &&
		     maletheta->theta[whichvary - 2] <= 0.5);
      else
	continue_ = (maletheta->theta[whichvary - 2] <= finaltheta &&
		     maletheta->theta[whichvary - 1] >= 0.0);
    }
#if !defined(DOS)
    if (checkpoint_counter > 0) 
      checkpoint_counter--;
    iterpeds_counter++;
#endif  /* !defined(DOS) */
  } while (continue_);
}
#endif  /* if !PARALLEL */

/* End. */