/*This file contains modified versions of the routines to do single
/* likelihood computations for the LINKAGE package. The modifications herein
/* are mostly needed to handle loops better and are described in the
/* paper: 
/* A. A. Schaffer, S. K. Gupta, K. Shriram, R. W. Cottingham, Jr. 
/* Avoiding Recomputation in Linkage Analysis*/
/* Human Heredity 44(1994), pp. 225-237. */

#include "commondefs.h"
#ifndef LESSMEMORY
#include "moddefs.h"
#endif

Void invert(m, n, det)
double (*m)[maxtrait];
int n;
double *det;
{
  covmatrix v;
  double val;
  int i, j, k;

  *det = 1.0;
  for (i = 1; i <= n; i++) {
    val = m[i - 1][i - 1];
    for (k = 0; k <= i - 2; k++)
      val -= v[k][i - 1] * v[k][i - 1];
    *det *= val;
    v[i - 1][i - 1] = sqrt(val);
    for (j = i; j < n; j++) {
      val = m[i - 1][j];
      for (k = 0; k <= i - 2; k++)
	val -= v[k][i - 1] * v[k][j];
      v[i - 1][j] = val / v[i - 1][i - 1];
      v[j][i - 1] = 0.0;
    }
  }
  for (i = 1; i <= n; i++) {
    m[i - 1][i - 1] = 1 / v[i - 1][i - 1];
    for (j = i + 1; j <= n; j++) {
      val = 0.0;
      for (k = i - 1; k <= j - 2; k++)
	val -= m[k][i - 1] * v[k][j - 1];
      m[j - 1][i - 1] = val / v[j - 1][j - 1];
      m[i - 1][j - 1] = 0.0;
    }
  }
  for (i = 0; i < n; i++) {
    for (j = 0; j < i + 1; j++) {
      val = 0.0;
      for (k = j; k < n; k++)
	val += m[k][i] * m[k][j];
      v[i][j] = val;
      v[j][i] = val;
    }
  }
  memcpy(m, v, sizeof(covmatrix));
}

/* static void collapsedown (thisperson *p, struct LOC_likelihood *LINK);*/
static void collapsedown ();
static void precollapsedown();

Void cleanup(p, LINK)
thisperson **p;
struct LOC_likelihood *LINK;
{
  thisperson *WITH;

  WITH = *p;
  if ((!(WITH->loopneeded)) || looplastgen) {   /*changed by A. A. Schaffer*/
    Free(WITH->gen);
    WITH->gen = NULL;
  }
  WITH->newgenexists = false;
}  /*cleanup*/

Local Void prob(p, LINK)
thisperson **p;
struct LOC_likelihood *LINK;
{
  int i,j;
  thisperson *WITH;
  thisarray *WITH1;
  unsigned char *tempflag1; /*Next two declarations by R. M. Idury*/
  double *tempwith1; 
  int FORLIM;

  WITH = *p;
  if (WITH->newgenexists) return;  /*Fixed by A. A. Schaffer*/
  if (WITH->gen == NULL) {
    WITH->gen = (thisarray *) Malloc(sizeof(thisarray));
    if (WITH->gen == NULL)
      malloc_err("gen field");
    WITH->newgenexists = true; /*Inserted by A. A. Schaffer*/
    /*Code to initialze genarray and sparseflag put in by R. M. Idury*/
    tempwith1 = WITH->gen->genarray;
    tempflag1 = WITH->gen->sparseflag;
    for(i=0;i<fgeno;i++) {
      tempflag1[i] = 0;
      tempwith1[i] = 0.0;
    }
    if (WITH->inloop > 0) {
      if (looppers[LINK->thisped - 1][WITH->inloop - 1][0] == *p) {
	tempwith1[LINK->loopgen[WITH->inloop - 1] - 1] = 
          LINK->holdpoint[WITH->inloop - 1]->
	    genarray[LINK->loopgen[WITH->inloop - 1] - 1];
	tempflag1[LINK->loopgen[WITH->inloop - 1] - 1] = 1;
      }
      else {
	tempwith1[LINK->loopgen[WITH->inloop - 1] - 1] = 1.0;
	tempflag1[LINK->loopgen[WITH->inloop - 1] - 1] = 1;
      }
      return;
    }
    getvect(*p, LINK);
  }
  if ((*p)->pa == NULL && (WITH->inloop==0))
    LINK->nuscale++;
}  /*prob*/


/*pseudoprob does just the scaling function of prob */
/*This routine written by A. A. Schaffer*/
Local Void pseudoprob(p, LINK)
thisperson **p;
struct LOC_likelihood *LINK;
{
  thisperson *WITH;

  WITH = *p;
  if (!(WITH->newgenexists))
    if ((*p)->pa == NULL) {
      LINK->nuscale++;
      if (WITH->inloop > 0) {
	LINK->nuscale--;
      }
    }
}  /*pseudoprob*/

Void initseg(LINK)
struct LOC_seg *LINK;
{
  if ((*LINK->p)->male) {
    LINK->nfirst = mgeno;
    LINK->nsecond = fgeno;
    LINK->firstsex = maletheta;
    LINK->secondsex = femaletheta;
    LINK->pf = mutmale;
    LINK->ps = mutfemale;
  } else {
    LINK->nfirst = fgeno;
    LINK->nsecond = mgeno;
    LINK->firstsex = femaletheta;
    LINK->secondsex = maletheta;
    LINK->pf = mutfemale;
    LINK->ps = mutmale;
  }
  prob(&LINK->father, LINK->LINK);
  prob(&LINK->mother, LINK->LINK);
  LINK->child = LINK->father->foff;
  nchild = 0;
  do {
    prob(&LINK->child, LINK->LINK);
    if (LINK->child->ma == LINK->mother && (!LINK->child->up)) {
      nchild++;
      childarray[nchild - 1] = LINK->child;  /*Line added by A. A. Schaffer*/
      thischild[nchild - 1] = LINK->child->gen;
      malechild[nchild - 1] = LINK->child->male;
    }
    LINK->child = LINK->child->nextpa;
  } while (LINK->child != NULL);   /*initseg*/
  if (nchild > maxchild) {
    fprintf(stderr, "\nA nuclear family has more children than maxchild");
    fprintf(stderr, "\nThe program will exit politely to allow you to");
    fprintf(stderr, "\nincrease maxchild in the file commondefs.h");
    exit(EXIT_FAILURE);
  }
}


/*pseudoseg simulates the calls that initseg would make to prob by
using pseudoprob just to get scaling right*/
/*This routine written by A. A. Schaffer*/
Local Void pseudoseg(LINK)
struct LOC_seg *LINK;
{
  pseudoprob(&LINK->father, LINK->LINK);
  pseudoprob(&LINK->mother, LINK->LINK);
  LINK->child = LINK->father->foff;
  nchild = 0;
  do {
    pseudoprob(&LINK->child, LINK->LINK);
    if (LINK->child->ma == LINK->mother && (!LINK->child->up)) {
      nchild++;
    }
    LINK->child = LINK->child->nextpa;
  } while (LINK->child != NULL);   
  LINK->LINK->nuscale++;  /*does extra scaling skipped in exitseg*/
} /*pseudoseg*/


Void exitseg(LINK)
struct LOC_seg *LINK;
{
  LINK->LINK->nuscale++;
  LINK->child = LINK->father->foff;
  do {
    if (LINK->child->ma == LINK->mother && (!LINK->child->up))
      cleanup(&LINK->child, LINK->LINK);
    LINK->child = LINK->child->nextpa;
  } while (LINK->child != NULL);   /*exitseg*/
}


/*computenumhaps computes the weighted sum over all genotypes of
how many haplotypes can be passed on to a child.
The following explanation assumes some basic understanding of
combinatorics at the undergraduate level. The total number of
haplotypes that can be passed on is
a weighted sum of the number of genotypes, where
different genotypes get different weights,  We first
classify the genotypes by heterozygosity pattern, where 0
means homozygous and 1 means heterozygous. There are mlocus loci
and therefore 2**mlocus different heterozygosity patterns.
The number of genotypes of a given heterozygosity pattern can be computed
by considering one locus at a time and multiplying the contributions
of all the loci, since they are independent.
A homozygous locus with A alleles contributes a factor of A.
The first heterozygous locus contributes a factor of A *(A-1)/2 if
it has A alleles.
Any other heterozygous locus with B alleles, contributes a factor of B*(B-1).

The weight of a genotype is 1 if it is homozygous; otherwise it
is 2**(h-1), where h is the number of heterozygous loci in the genotype.*/

int computenumhaps()
{
  int numpatterns;
  hapvector *patternmatrix;
  int p; /*index over heterozygosity patterns*/
  int l; /*index over loci*/
  int value; /*accumulator for return value*/
  int term; /*term for this heterozygosity pattern*/
  boolean hetyet; /*have we seen a heterozygous locus yet*/ 

  numpatterns = 2;
  for (i=2; i<=mlocus;i++)
  numpatterns*=2;
  patternmatrix = (hapvector *) malloc(numpatterns*sizeof(hapvector));
  if (patternmatrix == NULL)
    malloc_err("patternmatrix");
  for(l=0; l < mlocus; l++)
    patternmatrix[0][l]=0;
  for(p=1; p< numpatterns; p++) {
    l = 0;
    do{
      patternmatrix[p][l] = !patternmatrix[p-1][l];
      l++;
    } while(!patternmatrix[p][l-1]);
    do{
      patternmatrix[p][l] = patternmatrix[p-1][l];
      l++;
    } while(l < mlocus);
  }

  value = 2; /*offset by 1 because 0 index has special meaning*/
  for(p=0; p < numpatterns; p++){
    term = 1;
    hetyet = 0;
    for(l=0; l < mlocus; l++){
      if (!patternmatrix[p][l]) {
        term *= thislocus[l]->nallele;
      }
      else 
        if (!hetyet) {
          hetyet = 1;
          term *= (thislocus[l]->nallele * (thislocus[l]->nallele -1))/2;
        }
        else 
          term *= thislocus[l]->nallele * (thislocus[l]->nallele - 1) * 2;
    }      
    value+=term;
  }
  return(value);
}

/*getgeneindices is a preprocessing routine that compute the encodings
of all the haplotypes and genotypes and precomputes for each
genotype a list of the haplotypes it can pass on in one array
insted of using indirect addressing by genotype*/

Void getgeneindices()
{
  int locfstart, locfend, locfseg;
  int f1, f2, i, j, first;
  int finalfence;

/*find number of recombination probabilities*/  
  nuneed = 7;
  for (i= 3; i<= mlocus; i++)
    nuneed = 5 * nuneed - 3;
#ifndef LESSMEMORY
  segprob2 = (double *) malloc(nuneed * nuneed * sizeof(double));
  if (segprob2 == NULL) {
    fprintf(stderr, "\nYou probably need to use the slower version of this program");
    malloc_err("segprob2");
  }
#endif
/*find size of isozygote class and it's square for joint classes*/
  maxclasssize = 2;
  for (i = 3; i<= mlocus; i++)
    maxclasssize *= 2;
  maxisozygclass = maxclasssize * maxclasssize;
  segval = (double *) malloc(maxisozygclass * sizeof(double));
  if (segval == NULL)
    malloc_err("segval");
  tempseg = (double *) malloc(maxchild *maxisozygclass * sizeof(double));
  if (tempseg == NULL)
    malloc_err("tempseg");
  segindex = (unsigned *) malloc(4 * maxisozygclass * sizeof(unsigned));
  if (segindex == NULL)
    malloc_err("segindex");
#ifndef LESSMEMORY
  indpool = (unsigned *) malloc(2 * maxclasssize * fgeno * sizeof (unsigned));
  if (indpool == NULL) {
    fprintf(stderr, "\nYou probably need to use the slower version of this program");
    malloc_err("indpool");
  }
  invpool = (unsigned *) malloc(2 * maxclasssize * fgeno * sizeof (unsigned));
  if (invpool == NULL) {
    fprintf(stderr, "\nYou probably need to use the slower version of this program");
    malloc_err("invpool");
  }
  nextpool = (unsigned *) malloc(2 * maxclasssize * fgeno * sizeof (unsigned));
  if (nextpool == NULL) {
    fprintf(stderr, "\nYou probably need to use the slower version of this program");
    malloc_err("nextpool");
  }
  partialprob = (pairsums *) malloc(maxclasssize * sizeof(pairsums));
  if (partialprob == NULL) {
    fprintf(stderr, "\nYou probably need to use the slower version of this program");
    malloc_err("partialprob");
  }
#endif
#ifdef LESSMEMORY
  tempseg2 = (double *) malloc(maxchild *maxisozygclass * sizeof(double));
  if (tempseg2 == NULL)
    malloc_err("tempseg2");
#endif

  nuprobclass = 2;
  for(i = 2; i <= mlocus; i++)
    nuprobclass = 3 * nuprobclass - 1;
  finalfence = computenumhaps();
  haps1 = (unsigned short *) malloc(finalfence * sizeof(unsigned short));
  if (haps1 == NULL)
    malloc_err("haps1");
  haps2 = (unsigned short *) malloc(finalfence * sizeof(unsigned short));
  if (haps2 == NULL)
    malloc_err("haps2");
  hind = (unsigned int *) malloc(finalfence * sizeof(unsigned int));
  if (hind == NULL)
    malloc_err("hind");
  currentfence = 1;
  for (first = 0; first < fgeno; first++) {
    base[first] = currentfence;
    locfstart = probstart[first];
    locfend = probend[first];
    locfseg = segstart[first];
    for (i = locfstart-1; i < locfend; i++) {
      f1 = invgenenum1[locfseg - 1];
      haps1[currentfence] = f1;
      f2 = invgenenum2[locfseg - 1];
      haps2[currentfence] = f2;
      hind[currentfence++] = i;
      locfseg++;
    }
    fence[first] = currentfence;
  } 


  for(i = 1; i <= nuhap; i++) {
    for(j = i; j <= nuhap; j++) {
      invgenenum1[gennustruct->genenumber[i-1][j-1]-1] = i;
      invgenenum2[gennustruct->genenumber[i-1][j-1]-1] = j;
    }
  }
}  /* getgeneindices */

/*Routine determines in which way to do the next nuclear family genarray
update. Modified by R. M. Idury to distinguish between old versions
of segup (called oldsegup) and new versions of segdown and segup.
Modified by A. A. Schaffer to avoid doing the update in
certain loop traversal pedigree situations */
Void seg(p_, q_, r_, peel, LINK)
thisperson **p_, **q_, **r_;
direction peel;
struct LOC_likelihood *LINK;
{
  struct LOC_seg V;
  boolean phaseunkn;


  V.LINK = LINK;
  V.p = p_;
  V.q = q_;
  V.r = r_;
  phaseunkn = ((*V.q)->pa == NULL && (*V.q)->firstpass &&
	       (*V.q)->inloop == 0 && !disequi);
  if ((*V.p)->male) {
    V.father = *V.p;
    V.mother = *V.q;
  } else {
    V.father = *V.q;
    V.mother = *V.p;
  }
  if (peel == peelup) {
    if (((*V.p)->pa == NULL) && phaseunkn && approximate) {
      if (firstapprox && firsttime) {
	switch (thispath) {

	case auto_:
	  segtop(&V);
	  break;

	case mauto:
	  segtop(&V);
	  break;

	case sex:
          segsextop(&V);
	  break;

	case msex:
	  segsextop(&V);
	  break;
	}
      } else {  /*first approximate not first time*/
	if (firstapprox) {
	  switch (thispath) {

	  case auto_:
	    segctop(&V);
	    break;

	  case mauto:
	    segctop(&V);
	    break;

	  case sex:
	    segsexctop(&V);
	    break;

	  case msex:
	    segsexctop(&V);
	    break;
	  }
	} else {  /*approximate*/
	  switch (thispath) {

	  case auto_:
	    segcapprox(&V);
	    break;

	  case mauto:
	    segcapprox(&V);
	    break;

	  case sex:
	    segsexctop(&V);
	    break;

	  case msex:
	    segsexctop(&V);
	    break;
	  }
	}
      }
    } else {  /*do not approximate*/
      if (phaseunkn) {
	if (firsttime) {
	  switch (thispath) {

	  case auto_:
            if (((*V.p)->loopdepend) || loopfirstgen)  
	      segup(&V);
	    else
	      pseudoseg(&V);
	    break;

	  case mauto:
	    segtop(&V);
	    break;

	  case sex:
            if (((*V.p)->loopdepend) || loopfirstgen)
	      segsexup(&V);
            else
              pseudoseg(&V);
	    break;

	  case msex:
	    segsextop(&V);
	    break;
	  }
	} else {  /*not firsttime*/
	  switch (thispath) {

	  case auto_:
            if (((*V.p)->loopdepend) || loopfirstgen)  
	      segup(&V);
	    else
	      pseudoseg(&V);
	    break;

	  case mauto:
	    segctop(&V);
	    break;

	  case sex:
            if (((*V.p)->loopdepend) || loopfirstgen)
	      segsexup(&V);
            else
              pseudoseg(&V);
	    break;

	  case msex:
	    segsexctop(&V);
	    break;
	  }
	}
      } else {  /*phaseinfo*/
	switch (thispath) {

	case auto_:
          if (((*V.p)->loopdepend) || loopfirstgen)  
	    segup(&V);
          else
            pseudoseg(&V);
	  break;

	case mauto:
	  oldsegup(&V);
	  break;

	case sex:
          if (((*V.p)->loopdepend) || loopfirstgen)
	    segsexup(&V);
          else
            pseudoseg(&V);
	  break;

	case msex:
	  oldsegsexup(&V);
	  break;
	}
      }
    }
  } else {  /*not peelup*/
    switch (thispath) {

    case auto_:
      if (((*V.r)->loopdepend) || loopfirstgen)
        segdown(&V);
      else
        pseudoseg(&V);
      break;

    case mauto:
      msegdown(&V);
      break;

    case sex:
      if (((*V.r)->loopdepend) || loopfirstgen)
        segsexdown(&V);
      else
        pseudoseg(&V);
      break;

    case msex:
      msegsexdown(&V);
      break;
    }
  }
  (*V.q)->firstpass = false;
  (*V.p)->firstpass = false;
}  /*seg*/

Local Void collapseup(p, LINK)
thisperson *p;
struct LOC_likelihood *LINK;
{
  thisperson *q, *child, *nextchild;
  boolean down;

  p->done = true;
  if (p->foff == NULL)
    return;
  down = false;
  child = p->foff;
  while (child != NULL) {
    down = false;
    if (p->male)
      q = child->ma;
    else
      q = child->pa;
    if (!q->done) {
      collapsedown(q, LINK);
      nextchild = child;
      while (nextchild != NULL) {
	if (nextchild->pa == q || nextchild->ma == q) {
	  if (!nextchild->up)
	    collapseup(nextchild, LINK);
	  else
	    down = true;
	}
	if (p->male)
	  nextchild = nextchild->nextpa;
	else
	  nextchild = nextchild->nextma;
      }
      if (q->multi)
	collapseup(q, LINK);
      if (!down)
	seg(&p, &q, &child, peelup, LINK);
      else
	collapsedown(p, LINK);
    }
    if (p->male)
      child = child->nextpa;
    else
      child = child->nextma;
  }
}  /*collapseup*/


Local Void precollapseup(p)
thisperson *p;
{
  thisperson *q, *child, *nextchild, *childonlist;
  boolean down;
  p->done = true;
  if (p->foff == NULL)
    return;
  down = false;
  child = p->foff;
  while (child != NULL) {
    down = false;
    if (p->male)
      q = child->ma;
    else
      q = child->pa;
    if (!q->done) {
      precollapsedown(q);
      nextchild = child;
      while (nextchild != NULL) {
	if (nextchild->pa == q || nextchild->ma == q) {
	  if (!nextchild->up)
	    precollapseup(nextchild);
	  else
	    down = true;
	}
	if (p->male)
	  nextchild = nextchild->nextpa;
	else
	  nextchild = nextchild->nextma;
      }
      if (q->multi)
	precollapseup(q);
      if (!down) {
        childonlist = p->foff;
        while (childonlist != NULL) {
          if (childonlist->loopdepend) {
            p->loopdepend = true;
            p->loopneeded = false;
          }
	  if (p->male)
	    childonlist = childonlist->nextpa;
	  else
	    childonlist = childonlist->nextma; 
	}
      }
      else
	precollapsedown(p);
    }
    if (p->male)
      child = child->nextpa;
    else
      child = child->nextma;
  }
}  /*precollapseup*/



Local Void collapsedown(p, LINK)
thisperson *p;
struct LOC_likelihood *LINK;
{
  if (p->pa == NULL)
    return;
  p->up = true;
  collapseup(p->pa, LINK);
  seg(&p->pa, &p->ma, &p, peeldown, LINK);
}  /*collapsedown*/


Local Void precollapsedown(p, LINK)
thisperson *p;
struct LOC_likelihood *LINK;
{
  if (p->pa == NULL)
    return;
  p->up = true;
  precollapseup(p->pa);
  if (p->pa->loopdepend || p->ma->loopdepend) {
    p->loopdepend = true;
    p->loopneeded = false;
  }
}  /*precollapsedown*/


Local Void riskcumul(LINK)
struct LOC_likelihood *LINK;
{
  int i;
  thisarray *WITH;
  int FORLIM;

  WITH = LINK->proband->gen;
  if (sexlink && LINK->proband->male) {
    FORLIM = mgeno;
    for (i = 0; i < FORLIM; i++) {
      if (riskmale[i])
	LINK->hetero += WITH->genarray[i];
    }
    return;
  }
  FORLIM = fgeno;
  for (i = 0; i < FORLIM; i++) {
    if (risk2[i])
      LINK->homo += WITH->genarray[i];
    else if (risk1[i])
      LINK->hetero += WITH->genarray[i];
  }
}  /*riskcumul*/


Local Void riskcalc(LINK)
struct LOC_likelihood *LINK;
{
  double normal;

  LINK->homo /= like;
  LINK->hetero /= like;
  normal = 1 - LINK->homo - LINK->hetero;
  fprintf(outfile, "RISK FOR PERSON %6ld IN PEDIGREE %7ld\n",
	  LINK->proband->id, LINK->proband->ped);
  if (!LINK->proband->male || !sexlink)
    fprintf(outfile, "HOMOZYGOTE CARRIER   : %8.5f\n", LINK->homo);
  if (!LINK->proband->male || !sexlink)
    fprintf(outfile, "HETEROZYGOTE CARRIER : %8.5f\n", LINK->hetero);
  else
    fprintf(outfile, "MALE CARRIER         : %8.5f\n", LINK->hetero);
  fprintf(outfile, "NORMAL               : %8.5f\n", normal);
  printf("RISK FOR PERSON %6ld IN PEDIGREE %7ld\n",
	 LINK->proband->id, LINK->proband->ped);
  if (!LINK->proband->male || !sexlink)
    printf("HOMOZYGOTE CARRIER   : %8.5f\n", LINK->homo);
  if (!LINK->proband->male || !sexlink)
    printf("HETEROZYGOTE CARRIER : %8.5f\n", LINK->hetero);
  else
    printf("MALE CARRIER         : %8.5f\n", LINK->hetero);
  printf("NORMAL               : %8.5f\n", normal);
}  /*riskcalc*/

/*pollutedescendants finds which people descended from startper by
depth first search and ensures that they will be reevaluated for
each genotype of the loopbreaker */
/*This routine written by A. A. Schaffer*/

Local Void pollutedescendants(startper)
thisperson *startper;
{
  thisperson *nextchild;

  startper->loopdepend=true;
  nextchild = startper->foff;
  while(nextchild != NULL) {
    if (nextchild->loopdepend == false) 
      pollutedescendants(nextchild);
    nextchild = nextchild->nextpa;
  }
  nextchild = startper->foff;
  while(nextchild != NULL) {
    if (nextchild->loopdepend == false) 
      pollutedescendants(nextchild);
    nextchild = nextchild->nextma;
  }
}

Void likelihood(thisped_, proband_)
int thisped_;
thisperson *proband_;
{
  int ti, tj; /* put in by R. M. Idury */
  struct LOC_likelihood V;
  int loopmax[maxloop];
  double tmplike;
  int i, j;
  boolean gocalc, alldone;
  thisperson *WITH;
  thisarray *WITH1;
  int FORLIM1;
  int lastgen; /*last genotype that first loop breaker can have*/

  V.thisped = thisped_;
  V.proband = proband_;
  if (!informative[V.thisped - 1]) {
    like = 0.0;
    return;
  }
  V.homo = 0.0;
  V.hetero = 0.0;
  tmplike = 0.0;
  alldone = false;
  V.nuscale = 0;
  /* Next two pieces of code added by A. A. Schaffer*/
  /* initialize loop structure; all descendants of place where
     first loop is broken will need to be update for each genotype*/
  for (i = 1; i <= totperson; i++) {
    person[i]->loopdepend = false;
    person[i]->loopneeded = false;
  }
  if (looppers[V.thisped-1][0][0] != NULL){
    pollutedescendants(looppers[V.thisped-1][0][0]);
    pollutedescendants(looppers[V.thisped-1][0][1]);
    precollapseup(V.proband);
    precollapsedown(V.proband);
  }
  for (i = 0; i < maxloop; i++) {
    V.loopgen[i] = 1;
    loopmax[i] = 1;
    V.holdpoint[i] = NULL;
    if (looppers[V.thisped - 1][i][0] != NULL) {
      WITH = looppers[V.thisped - 1][i][0];
      WITH->gen = (thisarray *)Malloc(sizeof(thisarray));
      WITH->newgenexists = true; /*Added by A. A. Schaffer*/
      WITH1 = WITH->gen;
      FORLIM1 = fgeno;
      for (j = 0; j < FORLIM1; j++)
	WITH1->genarray[j] = 0.0;
      getvect(looppers[V.thisped - 1][i][0], &V);
      if (looppers[V.thisped - 1][i][0]->pa == NULL)
	V.nuscale++;
      V.holdpoint[i] = WITH->gen;
      WITH->gen = NULL;
      WITH->newgenexists = false;  /*Added by A. A. Schaffer*/
      if (WITH->male)
	loopmax[i] = mgeno;
      else
	loopmax[i] = fgeno;
    }
  }
  V.loopgen[0] = 0;
  /*find last nonzero genotype in first loop genarray*/
  /*Added by A. A. Schaffer*/
  if (looppers[V.thisped-1][0][0] != NULL){
    for (j = loopmax[0]-1; j>=0; j--){
      if (V.holdpoint[0]->genarray[j] > 0.0){
        lastgen =j;
        break;
      }
    }
  }
  /*Next line added by A. A. Schaffer*/
  loopfirstgen = 1;
  do {
    i = 1;
    /*Next line added by A. A. Schaffer*/
    looplastgen = 0;
    do {
      V.loopgen[i - 1]++;
      if (V.loopgen[i - 1] > loopmax[i - 1]) {
	V.loopgen[i - 1] = 1;
        /*Added by A. A. Schaffer*/
        if (i == 1)
          loopfirstgen = 1;
      }
      else
	i = maxloop;
      i++;
    } while (i <= maxloop);
    gocalc = true;
    for (i = 0; i < maxloop; i++) {
      /*ML change*/
      if (V.holdpoint[i] != NULL) {
	if (V.holdpoint[i]->genarray[V.loopgen[i] - 1] == 0.0)
	  gocalc = false;
        else               /*else clause added by A. A. Schaffer*/
          if ((i == 0) && ((V.loopgen[i] - 1) == lastgen))
            looplastgen = 1;
      }
    }
    if (gocalc) {
      FORLIM1 = totperson;
      for (i = 1; i <= FORLIM1; i++) {
	WITH = person[i];
        /*next four lines written by A. A. Schaffer*/
        if (loopfirstgen || !(WITH->loopneeded)) {
          WITH->gen = NULL;
        }
        WITH->newgenexists = false;
	WITH->done = false;
	WITH->up = false;
      }

/* accumulate segprob entries */
/* next test written by R. M. Idury*/
#ifndef LESSMEMORY
      for(ti = 0; ti < nuneed; ti++)
	for(tj = 0; tj < nuneed; tj++)
	  segprob2[ti*nuneed+tj] = maletheta->segprob[ti] * femaletheta->segprob[tj];
      getprobtable();
#endif
      collapseup(V.proband, &V);
      collapsedown(V.proband, &V);
      like = 0.0;
      WITH1 = V.proband->gen;
      if (V.proband->male) {
	FORLIM1 = mgeno;
	for (i = 0; i < FORLIM1; i++)
	  like += WITH1->genarray[i];
      } else {
	FORLIM1 = fgeno;
	for (i = 0; i < FORLIM1; i++)
	  like += WITH1->genarray[i];
      }

      /*Added by Alex*/
      loopfirstgen = 0;
      tmplike += like;
      if (risk && like != 0.0)
	riskcumul(&V);
      FORLIM1 = totperson;
      for (i = 1; i <= FORLIM1; i++) {
	if (person[i]->gen != NULL)
	  cleanup(&person[i], &V);
      }
    }
    alldone = true;
    for (i = 0; i < maxloop; i++)
      alldone = (alldone && V.loopgen[i] == loopmax[i]);
  } while (!alldone);
  like = tmplike;
  if (risk && like != 0.0)
    riskcalc(&V);
  if (like == 0.0)
    like = zerolike;
  else
    like = log(like) - V.nuscale * log(segscale);
  for (i = 0; i < maxloop; i++) {
    if (V.holdpoint[i] != NULL) {
      Free(V.holdpoint[i]);
      V.holdpoint[i] = NULL;
    }
  }
}  /*likelihood*/