/* This file contains input code shared by modified versions of the */
/* ILINK,  LINKMAP, and MLINK programs */
/* 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.  */
/* Most of the  code  in this file is not changed from LINKAGE 5.1 */

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

/* Local variables for getlocations: */
struct LOC_getlocations {
  long ngene, nseg, here, there, start, nhet, thisseg;
  boolean rarepresent, riskhom, riskhet;
  hapvector hap1, hap2;
  boolean thishet[maxlocus];
} ;

Local boolean checkrare(LINK)
struct LOC_getlocations *LINK;
{
  long i;
  boolean check;
  long FORLIM;
  locusvalues *WITH;

  check = false;
  FORLIM = mlocus;
  for (i = 1; i <= FORLIM; i++) {
    if (nohom[i - 1]) {
      WITH = thislocus[i - 1];
      if ((WITH->freq[LINK->hap1[i] - 1]< minfreq) ||
	  (WITH->freq[LINK->hap2[i] - 1] < minfreq))
	check = true;
    }
  }
  return check;
}

Local Void checkrisk(riskhet, riskhom, LINK)
boolean *riskhet, *riskhom;
struct LOC_getlocations *LINK;
{
  *riskhet = false;
  *riskhom = false;
  if (LINK->hap1[risksys - 1] == riskall &&
      LINK->hap2[risksys - 1] == riskall)
    *riskhom = true;
  else {
    if ((LINK->hap1[risksys - 1] != riskall &&
	 LINK->hap2[risksys - 1] == riskall) ||
	(LINK->hap2[risksys - 1] != riskall &&
	 LINK->hap1[risksys - 1] == riskall))
      *riskhet = true;
  }
}

Local long gethapn(hap, LINK)
hapvector hap;
struct LOC_getlocations *LINK;
{
  long i, n, FORLIM;

  n = 1;
  FORLIM = mlocus;
  for (i = 1; i <= FORLIM; i++)
    n += increment[i - 1] * (hap[i - 1] - 1);
  return n;
}

/* Local variables for domalerisk: */
struct LOC_domalerisk {
  struct LOC_getlocations *LINK;
} ;

Local Void setrisk(LINK)
struct LOC_domalerisk *LINK;
{
  long n;

  n = gethapn(LINK->LINK->hap1, LINK->LINK);
  if (LINK->LINK->hap1[risksys - 1] == riskall)
    riskmale[n - 1] = true;
  else
    riskmale[n - 1] = false;
}

Local Void getriskhap(system, LINK)
long system;
struct LOC_domalerisk *LINK;
{
  long i;
  locusvalues *WITH;
  long FORLIM;

  WITH = thislocus[system - 1];
  FORLIM = WITH->nallele;
  for (i = 1; i <= FORLIM; i++) {
    LINK->LINK->hap1[system - 1] = i;
    if (system != mlocus)
      getriskhap(system + 1, LINK);
    else
      setrisk(LINK);
  }
}

Local Void domalerisk(LINK)
struct LOC_getlocations *LINK;
{
  struct LOC_domalerisk V;

  V.LINK = LINK;
  getriskhap(1L, &V);
}

/* Local variables for domutation: */
struct LOC_domutation {
  struct LOC_getlocations *LINK;
} ;

Local Void setmutation(LINK)
struct LOC_domutation *LINK;
{
  long i, n;

  n = gethapn(LINK->LINK->hap1, LINK->LINK);
  if (LINK->LINK->hap1[mutsys - 1] == thislocus[mutsys - 1]->nallele) {
    muthap[n - 1] = n;
    return;
  }
  i = LINK->LINK->hap1[mutsys - 1];
  LINK->LINK->hap1[mutsys - 1] = thislocus[mutsys - 1]->nallele;
  muthap[n - 1] = gethapn(LINK->LINK->hap1, LINK->LINK);
  LINK->LINK->hap1[mutsys - 1] = i;
}

Local Void getmuthap(system, LINK)
long system;
struct LOC_domutation *LINK;
{
  long i;
  locusvalues *WITH;
  long FORLIM;

  WITH = thislocus[system - 1];
  FORLIM = WITH->nallele;
  for (i = 1; i <= FORLIM; i++) {
    LINK->LINK->hap1[system - 1] = i;
    if (system != mlocus)
      getmuthap(system + 1, LINK);
    else
      setmutation(LINK);
  }
}

Local Void domutation(LINK)
struct LOC_getlocations *LINK;
{
  struct LOC_domutation V;

  V.LINK = LINK;
  getmuthap(1L, &V);
}

Local Void setnumbers(LINK)
struct LOC_getlocations *LINK;
{
  long nhap1, nhap2;

  LINK->ngene++;

  segstart[LINK->ngene - 1] = LINK->here + 1;
  probstart[LINK->ngene - 1] = LINK->there + 1;
  probend[LINK->ngene - 1] = LINK->there + LINK->nseg;

  LINK->there += LINK->nseg;

  nhap1 = gethapn(LINK->hap1, LINK);
  nhap2 = gethapn(LINK->hap2, LINK);
  gennustruct->genenumber[nhap1 - 1][nhap2 - 1] = LINK->ngene;
  gennustruct->genenumber[nhap2 - 1][nhap1 - 1] = LINK->ngene;

  if (minfreq != 0.0) {
    if (LINK->rarepresent)
      rare[LINK->ngene - 1] = true;
    else
      rare[LINK->ngene - 1] = false;
  } else
    rare[LINK->ngene - 1] = false;
  if (risk) {
    risk1[LINK->ngene - 1] = LINK->riskhet;
    risk2[LINK->ngene - 1] = LINK->riskhom;
  }

  LINK->thisseg++;
  invgenenum1[LINK->thisseg - 1] = nhap1;
  invgenenum2[LINK->thisseg - 1] = nhap2;
}

Local Void hapscr(system, nscramble, LINK)
long system, nscramble;
struct LOC_getlocations *LINK;
{
  long i, j;

  if (LINK->thishet[system - 1])
    nscramble++;
  if (system != mlocus)
    hapscr(system + 1, nscramble, LINK);
  else
    setnumbers(LINK);
  if (nscramble <= 1)
    return;
  if (LINK->hap1[system - 1] == LINK->hap2[system -1])
    return;
  i = LINK->hap1[system - 1];
  j = LINK->hap2[system - 1];
  LINK->hap1[system - 1] = j;
  LINK->hap2[system - 1] = i;
  if (system != mlocus)
    hapscr(system + 1, nscramble, LINK);
  else
    setnumbers(LINK);
  LINK->hap1[system - 1] = i;
  LINK->hap2[system - 1] = j;
}

Local Void sethap(system, LINK)
long system;
struct LOC_getlocations *LINK;
{
  long i, j;
  locusvalues *WITH;
  long FORLIM, FORLIM1;

  WITH = thislocus[system - 1];
  if (LINK->thishet[system - 1]) {
    FORLIM = WITH->nallele;
    for (i = 1; i < FORLIM; i++) {
      LINK->hap1[system - 1] = i;
      FORLIM1 = WITH->nallele;
      for (j = i + 1; j <= FORLIM1; j++) {
	LINK->hap2[system - 1] = j;
	if (system != mlocus)
	  sethap(system + 1, LINK);
	else {
	  LINK->rarepresent = checkrare(LINK);
	  if (risk)
	    checkrisk(&LINK->riskhet, &LINK->riskhom, LINK);
	  LINK->there = LINK->start;
	  LINK->thisseg = LINK->here;
	  hapscr(1L, 0L, LINK);
	  LINK->here += LINK->nseg;
	}
      }
    }
    return;
  }
  FORLIM = WITH->nallele;
  for (i = 1; i <= FORLIM; i++) {
    LINK->hap1[system - 1] = i;
    LINK->hap2[system - 1] = i;
    if (system != mlocus)
      sethap(system + 1, LINK);
    else {
      LINK->rarepresent = checkrare(LINK);
      if (risk)
	checkrisk(&LINK->riskhet, &LINK->riskhom, LINK);
      LINK->thisseg = LINK->here;
      LINK->there = LINK->start;
      hapscr(1L, 0L, LINK);
      LINK->here += LINK->nseg;
    }
  }
}

Local Void starthap(LINK)
struct LOC_getlocations *LINK;
{
  long i, FORLIM;

  LINK->nseg = 1;
  FORLIM = LINK->nhet;
  for (i = 2; i <= FORLIM; i++)
    LINK->nseg *= 2;
  sethap(1L, LINK);
  LINK->start = LINK->there;
}

Local Void gethet1(system, LINK)
long system;
struct LOC_getlocations *LINK;
{
  LINK->thishet[system - 1] = false;
  if (system != mlocus)
    gethet1(system + 1, LINK);
  else
    starthap(LINK);
  LINK->thishet[system - 1] = true;
  LINK->nhet++;
  if (system != mlocus)
    gethet1(system + 1, LINK);
  else
    starthap(LINK);
  LINK->nhet--;
}


Void getlocations()
{
  struct LOC_getlocations V;

  V.nhet = 0;
  V.here = 0;
  V.there = 0;
  V.ngene = 0;
  V.start = 0;
  gethet1(1L, &V);
  if (mutsys != 0)
    domutation(&V);
  if (sexlink && risk)
    domalerisk(&V);
}


Void inputerror(nerror, par1, par2, LINK)
long nerror, par1, par2;
struct LOC_inputdata *LINK;
{
  printf("Fatal error detected in procedure inputdata\n");
  switch (nerror) {

  case 0:
    printf("Number of loci %2ld exceeds the constant maxlocus\n", par1);
    break;

  case 1:
    printf("Number of loci read %2ld. Less than minimum of 1\n", par1);
    break;

  case 2:
    printf(
      "Error detected reading loci order. Locus number %2ld in position %2ld exceeds number of loci\n",
      par2, par1);
    break;

  case 3:
    printf(
      "Error detected reading loci order. Illegal locus number %2ld in position %2ld\n",
      par2, par1);
    break;

  case 4:
    printf(
      "Error detected reading loci order. Locus number repeated in positions %2ld and %2ld\n",
      par1, par2);
    break;

  case 5:
    printf(
      "Error detected reading locus description. Illegal locus type %2ld for locus %2ld\n",
      par2, par1);
    break;

  case 6:
    printf(
      "Error detected reading locus description for system %2ld. Number of alleles  %2ld exceeds maxall\n",
      par1, par1);
    break;

  case 7:
    printf(
      "Error detected reading locus description for system %2ld. Illegal number of alleles  %2ld\n",
      par1, par2);
    break;

  case 8:
    printf(
      "Error detected reading locus description for system %2ld. Number of factors  %2ld exceeds maxfact\n",
      par1, par2);
    break;

  case 9:
    printf(
      "Error detected reading locus description for system %2ld. Illegal number of factors  %2ld\n",
      par1, par2);
    break;

  case 10:
    printf(
      "Error detected reading locus description for system %2ld. Alleles not codominant\n",
      par1);
    break;

  case 11:
    printf("Error detected reading pedigree record %2ld. Illegal code for sex %2ld\n",
	   par1, par2);
    break;

  case 12:
    printf(
      "Error detected reading pedigree record at pedigree%2ld. Maximum number of pedigree records exceeded\n",
      par1);
    break;

  case 13:
    printf(
      "Error detected reading pedigree record %2ld. Maximum number of individuals exceeded\n",
      par1);
    break;

  case 14:
    printf(
      "Error detected reading pedigree record %2ld. Illegal binary factor code %2ld\n",
      par1, par2);
    break;

  case 15:
    printf(
      "Error detected reading pedigree record %2ld. No allelic pair for genotype\n",
      par1);
    break;

  case 16:
    printf(
      "Error detected reading pedigree record %2ld. Allele number %2ld exceeds maxall\n",
      par1, par2);
    break;

  case 17:
    printf(
      "Error detected reading pedigree record %2ld. Illegal allele number %2ld\n",
      par1, par2);
    break;

  case 18:
    printf("Number of systems after factorization (%3ld) exceeds maxsystem\n",
	   par1);
    break;

  case 19:
    printf("Number of systems after factorization (%3ld) less than minimum of 1\n",
	   par1);
    break;

  case 20:
    printf("Number of recombination types (%3ld) exceeds maxrectype\n", par1);
    break;

  case 21:
    printf("Number of recombination types (%3ld) less than minimum of 1\n",
	   par1);
    break;

  case 22:
    printf(
      "End of file detected in tempdat by procedure readthg before all data found\n");
    break;

  case 23:
    printf(
      "Error detected reading iterated locus in datafile. Value (%3ld) greater than nlocus\n",
      par1);
    break;

  case 24:
    printf(
      "Error detected reading iterated locus in datafile. Illegal value (%3ld)\n",
      par1);
    break;

  case 25:
    printf("Number of iterated parameters greater then maxn\n");
    break;

  case 26:
    printf(
      "Error detected reading pedigree record %2ld. Liability class (%2ld) exceeds nclass\n",
      par1, par2);
    break;

  case 27:
    printf(
      "Error detected reading pedigree record %2ld. Illegal liability class (%2ld)\n",
      par1, par2);
    break;

  case 28:
    printf(
      "Error detected reading locus description for system%2ld. Liability classes (%3ld) exceed maxliab\n",
      par1, par2);
    break;

  case 29:
    printf(
      "Error detected reading locus description for system%2ld. Illegal number of liability classes (%3ld)\n",
      par1, par2);
    break;

  case 30:
    printf(
      "Error detected reading locus description for system%2ld. Penetrance out of range\n",
      par1);
    break;

  case 31:
    printf(
      "Error detected reading locus description for system%2ld. Number of traits (%3ld) exceeds maxtrait\n",
      par1, par2);
    break;

  case 32:
    printf(
      "Error detected reading locus description for system%2ld. Number of traits out of range (%3ld)\n",
      par1, par2);
    break;

  case 33:
    printf(
      "Error detected reading locus description for system%2ld. Variance must be positive\n",
      par1);
    break;

  case 34:
    printf(
      "Error detected reading locus description for system%2ld. Variance multiplier must be positive\n",
      par1);
    break;

  case 35:
    printf(
      "Error detected reading locus description for system%2ld. Risk allele %3ld) exceeds nallele\n",
      par1, par2);
    break;

  case 36:
    printf(
      "Error detected reading locus description for system%2ld. Illegal risk allele (%3ld)\n",
      par1, par2);
    break;

  case 37:
    printf("Error detected reading datafile. Risk locus %3ld) exceeds nlocus\n",
	   par2);
    break;

  case 38:
    printf("Error detected reading datafile. Illegal value for risk locus %3ld)\n",
	   par2);
    break;

  case 39:
    printf("Error detected reading datafile. Mutation locus %3ld) exceeds nlocus\n",
	   par2);
    break;

  case 40:
    printf(
      "Error detected reading datafile. Illegal value for mutation locus %3ld)\n",
      par2);
    break;

  case 41:
    printf(
      "Error detected reading datafile. Linkage disequilbirium is not allowed with this program\n");
    break;

  case 42:
    printf("Locus %5ld in lod score list exceeds nlocus %5ld\n", par1, par2);
    break;

  case 43:
    printf("Illegal locus number %5ld in lod score list\n", par1);
    break;
  }
}

Void inputwarning(nwarning, par1, par2, LINK)
long nwarning, par1, par2;
struct LOC_inputdata *LINK;
{
  printf("Warning number from procedure inputdata\n");
  switch (nwarning) {

  case 0:
    printf("Illegal sex difference parameter %2ld Parameter should be 0, 1, or 2\n",
	   par1);
    break;

  case 1:
    printf("Illegal interference parameter %2ld Lack of interference assumed\n",
	   par1);
    break;

  case 2:
    printf(
      "Illegal sex difference parameter %2ld Parameter must be 0 with sex-linked data\n",
      par1);
    break;

  case 3:
    printf(
      "Non-standard affection status%4ld interpreted as normal in pedigree record%5ld\n",
      par2, par1);
    break;
  }
}

Void readspeed(LINK)
struct LOC_inputdata *LINK;
{
  long i, a, b, sys;
  char ch;
  information *WITH;
  long FORLIM;

  while (!P_eof(speedfile)) {
    ch = getc(speedfile);
    if (ch == '\n')
      ch = ' ';
    if (ch == 'i' || ch == 'I') {
      fscanf(speedfile, "%c%ld", &ch, &i);
      if (ch == '\n')
	ch = ' ';
      person[i]->unknown = true;
      person[i]->store = (information *)Malloc(sizeof(information));
      if ((person[i]->store) == NULL)
        malloc_err("store field");
      WITH = person[i]->store;
      FORLIM = mlocus;
      for (sys = 0; sys < FORLIM; sys++) {
	for (a = 0; a < maxall; a++) {
	  for (b = 0; b < maxall; b++)
	    WITH->possible[sys][a][b] = false;
	}
      }
    } else {
      WITH = person[i]->store;
      fscanf(speedfile, "%ld%ld%ld", &sys, &a, &b);
      if (sys <= mlocus)
	WITH->possible[order[sys - 1] - 1][a - 1][b - 1] = true;
    }
    fscanf(speedfile, "%*[^\n]");
    getc(speedfile);
  }
}

/* Local variables for readped: */
struct LOC_readped {
  struct LOC_inputdata *LINK;
  long sequence;
  long startped[maxped], endped[maxped];
} ;

/* Local variables for getphenotype: */
struct LOC_getphenotype {
  struct LOC_readped *LINK;
  thisperson **p;
  long system;
} ;

Local Void readbin(phen, thislocus, LINK)
phenotype **phen;
locusvalues *thislocus;
struct LOC_getphenotype *LINK;
{
  long i, j;
  phenotype *WITH;
  long FORLIM;

  WITH = *phen;
  WITH->which = binary_;
  WITH->phenf = 0;
  FORLIM = LINK->LINK->LINK->nfactor[LINK->system - 1];
  for (i = 1; i <= FORLIM; i++) {
    fscanf(ipedfile, "%ld", &j);
    if (j != 0 && j != 1)
      inputerror(14L, (*LINK->p)->id, j, LINK->LINK->LINK);
    if (j == 1)
      WITH->phenf = ((long)WITH->phenf) | (1 << ((long)i));
  }
}

Local Void readnumber(phen, thislocus, LINK)
phenotype **phen;
locusvalues *thislocus;
struct LOC_getphenotype *LINK;
{
  long i, j;
  phenotype *WITH;

  WITH = *phen;
  WITH->which = binary_;
  WITH->phenf = 0;
  for (i = 1; i <= 2; i++) {
    fscanf(ipedfile, "%ld", &j);
    if (j > maxall)
      inputerror(16L, (*LINK->p)->id, j, LINK->LINK->LINK);
    if (j < 0)
      inputerror(17L, (*LINK->p)->id, j, LINK->LINK->LINK);
    if (j != 0)
      WITH->phenf = ((long)WITH->phenf) | (1 << ((long)j));
  }
}

Local Void readaff(phen, thislocus, LINK)
phenotype **phen;
locusvalues *thislocus;
struct LOC_getphenotype *LINK;
{
  long thisval;
  phenotype *WITH;

  WITH = *phen;
  WITH->which = affection;
  fscanf(ipedfile, "%ld", &thisval);
  if (thisval == missaff)
    WITH->aff = 0;
  else {
    if (thisval == affval)
      WITH->aff = 2;
    else {
      if (thisval != 1)
	inputwarning(3L, (*LINK->p)->id, thisval, LINK->LINK->LINK);
      WITH->aff = 1;
    }
  }
  if (thislocus->UU.U0.nclass == 1)
    WITH->liability = 1;
  else
    fscanf(ipedfile, "%ld", &(*phen)->liability);
  if (WITH->liability > thislocus->UU.U0.nclass)
    inputerror(26L, (*LINK->p)->id, WITH->liability, LINK->LINK->LINK);
  if (WITH->liability <= 0)
    inputerror(27L, (*LINK->p)->id, WITH->liability, LINK->LINK->LINK);
}

Local Void readquan(phen, thislocus, LINK)
phenotype **phen;
locusvalues *thislocus;
struct LOC_getphenotype *LINK;
{
  long i;
  double xval;
  phenotype *WITH;
  long FORLIM;

  WITH = *phen;
  if (!sexlink || !(*LINK->p)->male) {
    WITH->which = quantitative;
    FORLIM = thislocus->UU.U1.ntrait;
    for (i = 0; i < FORLIM; i++)
      fscanf(ipedfile, "%lf", &(*phen)->x[i]);
    WITH->missing = true;
    FORLIM = thislocus->UU.U1.ntrait;
    for (i = 0; i < FORLIM; i++) {
      if (WITH->x[i] != missval)
	WITH->missing = false;
    }
    return;
  }
  WITH->which = affection;
  fscanf(ipedfile, "%lf", &xval);
  if (xval == missval)
    WITH->aff = missaff;
  else {
    if (xval == affall)
      WITH->aff = affall;
    else
      WITH->aff = -11;
  }
  WITH->liability = 1;
  FORLIM = thislocus->UU.U1.ntrait;
  for (i = 2; i <= FORLIM; i++)
    fscanf(ipedfile, "%lf", &xval);
}

Local Void getphenotype(p_, LINK)
thisperson **p_;
struct LOC_readped *LINK;
{
  struct LOC_getphenotype V;
  long thisread;
  thisperson *WITH;
  long FORLIM;

  V.LINK = LINK;
  V.p = p_;
  WITH = *V.p;
  FORLIM = mlocus;
  for (thisread = 0; thisread < FORLIM; thisread++) {
    V.system = order[thisread];
    WITH->phen[V.system - 1] = NULL;
    if (thislocus[V.system - 1]->which != null_) {
      WITH->phen[V.system - 1] = (phenotype *)Malloc(sizeof(phenotype));
      if ((WITH->phen[V.system - 1]) == NULL)
        malloc_err("phen field");
    }    
    switch (thislocus[V.system - 1]->which) {

    case quantitative:
      readquan(&WITH->phen[V.system - 1], thislocus[V.system - 1], &V);
      break;

    case affection:
      readaff(&WITH->phen[V.system - 1], thislocus[V.system - 1], &V);
      break;

    case binary_:
      if (thislocus[V.system - 1]->format == 3)
	readnumber(&WITH->phen[V.system - 1], thislocus[V.system - 1], &V);
      else
	readbin(&WITH->phen[V.system - 1], thislocus[V.system - 1], &V);
      break;
    }
    if (lastpriv == V.system) {
      WITH->privphen = (phenotype *)Malloc(sizeof(phenotype));
      if (WITH->privphen == NULL)
        malloc_err("privphen field");
      switch (thislocus[V.system - 1]->privlocus->which) {

      case quantitative:
	readquan(&WITH->privphen, thislocus[V.system - 1]->privlocus, &V);
	break;

      case affection:
	readaff(&WITH->privphen, thislocus[V.system - 1]->privlocus, &V);
	break;
      }
    }
  }
}


Local Void getinformative(LINK)
struct LOC_readped *LINK;
{
  long i, j, k, l, m, count, nchild;
  thisperson *child;
  long FORLIM, FORLIM1, FORLIM2;
  phenotype *WITH;
  locusvalues *WITH1;
  long FORLIM3;

  if (fitmodel || risk) {
    FORLIM = nuped;
    for (i = 0; i < FORLIM; i++)
      informative[i] = true;
    return;
  }
  FORLIM = nuped;
  for (i = 0; i < FORLIM; i++) {
    informative[i] = false;
    FORLIM1 = LINK->endped[i];
    for (j = LINK->startped[i]; j <= FORLIM1; j++) {
      if (person[j]->foff != NULL) {
	nchild = 0;
	child = person[j]->foff;
	do {
	  nchild++;
	  if (person[j]->male)
	    child = child->nextpa;
	  else
	    child = child->nextma;
	} while (child != NULL);
	count = 0;
	if (nchild > 1 || nchild == 1 && person[j]->pa != NULL) {
	  FORLIM2 = mlocus;
	  for (k = 0; k < FORLIM2; k++) {
	    WITH = person[j]->phen[k];
	    WITH1 = thislocus[k];
	    if (WITH1->which != binary_)
	      count++;
	    else {
	      if (WITH->phenf == 0)
		count++;
	      else {
		l = 0;
		FORLIM3 = WITH1->nallele;
		for (m = 1; m <= FORLIM3; m++) {
		  if ((unsigned long)m < 32 &&
		      ((1L << m) & WITH->phenf) != 0)
		    l++;
		}
		if (l > 1)
		  count++;
	      }
	    }
	  }
	}
	if (count > 1)
	  informative[i] = true;
      }
    }
  }
}



Local Void getind(id, LINK)
long *id;
struct LOC_readped *LINK;
{
  fscanf(ipedfile, "%ld", id);
  if (*id == 0)
    return;
  *id += LINK->sequence;
  if (*id > maxind)
    inputerror(13L, *id, *id, LINK->LINK);
  if (person[*id] == NULL) {
    person[*id] = (thisperson *)Malloc(sizeof(thisperson));
    if (person[*id] == NULL)
      malloc_err("person");
  }
}  /*getind*/


Local Void multimarriage(p, LINK)
thisperson **p;
struct LOC_readped *LINK;
{
  thisperson *q, *child, *WITH;

  if ((*p)->foff == NULL) {
    (*p)->multi = false;
    return;
  }
  WITH = *p;
  if (WITH->male)
    q = WITH->foff->ma;
  else
    q = WITH->foff->pa;
  child = WITH->foff;
  (*p)->multi = false;
  do {
    if (WITH->male) {
      WITH->multi = (q == child->ma);
      child = child->nextpa;
    } else {
      WITH->multi = (q == child->pa);
      child = child->nextma;
    }
  } while (!(child == NULL || WITH->multi));
}  /*multimarriage*/


Void readped(LINK)
struct LOC_inputdata *LINK;
{
  struct LOC_readped V;
  /*profield reads column 9 of pedin.dat which indicates proband or
    loopbreaker status*/
  long i, newid, sex_, profield, newped, nuperson, thisone, thisped;
  thisperson *holdloop;
  long FORLIM;
  thisperson *WITH;

  V.LINK = LINK;
  for (i = 0; i <= maxind; i++)
    person[i] = NULL;
  V.sequence = 0;
  nuperson = 0;
  nuped = 1;
  
#if defined(GENERR)
  pedind[0]=0;
#endif

  for (i = 0; i < maxloop; i++) {
    looppers[nuped - 1][i][0] = NULL;
    looppers[nuped - 1][i][1] = NULL;
  }
  proband[nuped - 1] = NULL;
  fscanf(ipedfile, "%ld", &newped);
  thisped = newped;
  V.startped[0] = 1;
  while (!P_eof(ipedfile)) { 
    nuperson++;
    getind(&thisone, &V);
    if (proband[nuped - 1] == NULL)
      proband[nuped - 1] = person[thisone];
    WITH = person[thisone];
    WITH->ped = thisped;
    WITH->id = thisone;
    getind(&newid, &V);
    WITH->pa = person[newid];
    getind(&newid, &V);
    WITH->ma = person[newid];
    getind(&newid, &V);
    WITH->foff = person[newid];
    getind(&newid, &V);
    WITH->nextpa = person[newid];
    getind(&newid, &V);
    WITH->nextma = person[newid];
    fscanf(ipedfile, "%ld", &sex_);
    if (sex_ != 1 && sex_ != 2)
      inputerror(11L, WITH->id, sex_, LINK);
    if (sex_ == 1)
      WITH->male = true;
    else
      WITH->male = false;
    WITH->unknown = false;
    WITH->inloop = 0;
    WITH->loopdepend = false; /*A. A. Schaffer added this line and next*/
    WITH->loopneeded = false;
    fscanf(ipedfile, "%ld", &profield);
    /*Diagnostic added by A. A. Schaffer 7/25/94*/
    if ((profield - 1) > maxloop) {
      fprintf(stderr, "\nYour pedigree has more loops than allowed by the constant maxloop");
      fprintf(stderr, "\nYou must increase the constant maxloop defined in commondefs.h and recompile");
      fprintf(stderr, "\nYou are encouraged to read the loops.ps document distributed with FASTLINK");
      fprintf(stderr, "\nThe program will exit politely to allow you to correct the problem\n");
      exit(EXIT_FAILURE);
    }
    if (profield == 1)
      proband[nuped - 1] = person[thisone];
    else if (profield > 1 && profield - 1 <= maxloop) {
      if (looppers[nuped - 1][profield - 2][1] == NULL)
	looppers[nuped - 1][profield - 2][1] = person[thisone];
      else
	looppers[nuped - 1][profield - 2][0] = person[thisone];
    }
    getphenotype(&person[thisone], &V);
    fscanf(ipedfile, "%*[^\n]");
    getc(ipedfile);
    if (!P_eof(ipedfile))
      fscanf(ipedfile, "%ld", &newped);
    if (newped == thisped)
      continue;
    V.sequence += nuperson;
    V.endped[nuped - 1] = V.sequence;

#if defined(GENERR)
    pedind[nuped]=pedind[nuped-1]+nuperson;
#endif

    nuperson = 0;
    nuped++;
    if (nuped > maxped)
      inputerror(12L, newped, nuped, LINK);
    V.startped[nuped - 1] = V.sequence + 1;
    for (i = 0; i < maxloop; i++) {
      looppers[nuped - 1][i][0] = NULL;
      looppers[nuped - 1][i][1] = NULL;
    }
    proband[nuped - 1] = NULL;
    thisped = newped;
  }

#if defined(GENERR)
    pedind[nuped]=pedind[nuped-1]+nuperson;
#endif

  totperson = V.sequence + nuperson;
  V.endped[nuped - 1] = totperson;
  FORLIM = nuped;
  for (newped = 0; newped < FORLIM; newped++) {
    if (looppers[newped][0][1] != NULL && looppers[newped][0][0] == NULL)
      looppers[newped][0][0] = proband[newped];
    for (i = 0; i < maxloop; i++) {
      if (looppers[newped][i][0] == NULL)
	looppers[newped][i][1] = NULL;
      else {
	looppers[newped][i][0]->inloop = i + 1;
	looppers[newped][i][1]->inloop = i + 1;
	if (looppers[newped][i][0]->pa == NULL &&
	    looppers[newped][i][1]->pa != NULL) {
	  holdloop = looppers[newped][i][0];
	  looppers[newped][i][0] = looppers[newped][i][1];
	  looppers[newped][i][1] = holdloop;
	}
      }
    }
  }
  FORLIM = totperson;
  for (thisone = 1; thisone <= FORLIM; thisone++)
    multimarriage(&person[thisone], &V);
  getinformative(&V);
}  /*readped*/

/* Local variables for getlocus: */
struct LOC_getlocus {
  struct LOC_readloci *LINK;
  long system;
} ;

Local Void getbin(locus, system, LINK)
locusvalues **locus;
long *system;
struct LOC_getlocus *LINK;
{
  long i, j, k;
  locusvalues *WITH;
  long FORLIM, FORLIM1;

  fscanf(datafile, "%ld%*[^\n]", &LINK->LINK->LINK->nfactor[*system - 1]);
  getc(datafile);
  if (LINK->LINK->LINK->nfactor[*system - 1] > maxfact)
    inputerror(8L, *system, LINK->LINK->LINK->nfactor[*system - 1],
	       LINK->LINK->LINK);
  if (LINK->LINK->LINK->nfactor[*system - 1] <= 0)
    inputerror(9L, *system, LINK->LINK->LINK->nfactor[*system - 1],
	       LINK->LINK->LINK);
  WITH = *locus;
  FORLIM = WITH->nallele;
  for (i = 0; i < FORLIM; i++)
    WITH->UU.allele[i] = 0;
  FORLIM = WITH->nallele;
  for (i = 0; i < FORLIM; i++) {
    FORLIM1 = LINK->LINK->LINK->nfactor[*system - 1];
    for (j = 1; j <= FORLIM1; j++) {
      fscanf(datafile, "%ld", &k);
      if (k == 1)
	WITH->UU.allele[i] = ((long)WITH->UU.allele[i]) | (1 << ((long)j));
    }
  }
  fscanf(datafile, "%*[^\n]");
  getc(datafile);
}

Local Void getnumber(locus, system, LINK)
locusvalues **locus;
long *system;
struct LOC_getlocus *LINK;
{
  long i;
  locusvalues *WITH;
  long FORLIM;

  WITH = *locus;
  FORLIM = WITH->nallele;
  for (i = 1; i <= FORLIM; i++)
    WITH->UU.allele[i - 1] = 1L << ((int)i);
}


Local Void getpen(locus, LINK)
locusvalues **locus;
struct LOC_getlocus *LINK;
{
  long i, j, k, l;
  locusvalues *WITH;
  long FORLIM, FORLIM1, FORLIM2;

  WITH = *locus;
  fscanf(datafile, "%ld%*[^\n]", &(*locus)->UU.U0.nclass);
  getc(datafile);
  if (WITH->UU.U0.nclass > maxliab)
    inputerror(28L, LINK->system, WITH->UU.U0.nclass, LINK->LINK->LINK);
  if (WITH->UU.U0.nclass <= 0)
    inputerror(29L, LINK->system, WITH->UU.U0.nclass, LINK->LINK->LINK);
  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++) {
	fscanf(datafile, "%lf", &(*locus)->UU.U0.pen[i][j][2][l]);
	if ((unsigned)WITH->UU.U0.pen[i][j][2][l] > 1.0)
	  inputerror(30L, LINK->system, LINK->system, LINK->LINK->LINK);
	WITH->UU.U0.pen[i][j][1][l] = 1 - WITH->UU.U0.pen[i][j][2][l];
	WITH->UU.U0.pen[i][j][0][l] = 1.0;
	for (k = 0; k <= 2; k++)
	  WITH->UU.U0.pen[j + 1][i - 1][k][l] = WITH->UU.U0.pen[i][j][k][l];
      }
    }
    fscanf(datafile, "%*[^\n]");
    getc(datafile);
    FORLIM1 = WITH->nallele;
    for (i = 0; i < FORLIM1; i++)
      WITH->UU.U0.pen[0][i][0][l] = 1.0;
    if (sexlink) {
      FORLIM1 = WITH->nallele;
      for (i = 0; i < FORLIM1; i++) {
	fscanf(datafile, "%lf", &(*locus)->UU.U0.pen[0][i][2][l]);
	if ((unsigned)WITH->UU.U0.pen[0][i][2][l] > 1.0)
	  inputerror(30L, LINK->system, LINK->system, LINK->LINK->LINK);
      }
      FORLIM1 = WITH->nallele;
      for (i = 0; i < FORLIM1; i++)
	WITH->UU.U0.pen[0][i][1][l] = 1.0 - WITH->UU.U0.pen[0][i][2][l];
      fscanf(datafile, "%*[^\n]");
      getc(datafile);
    }
  }
}


Local Void getquan(locus, privelege, LINK)
locusvalues **locus;
boolean privelege;
struct LOC_getlocus *LINK;
{
  /*Get information of a quantitative locus. Privelege says whether it is
    a priveleged locus or not*/
  long i, j, k;
  locusvalues *WITH;
  long FORLIM, FORLIM1, FORLIM2;

  WITH = *locus;
  fscanf(datafile, "%ld%*[^\n]", &(*locus)->UU.U1.ntrait);
  getc(datafile);
  if (WITH->UU.U1.ntrait > maxtrait)
    inputerror(31L, LINK->system, WITH->UU.U1.ntrait, LINK->LINK->LINK);
  if (WITH->UU.U1.ntrait <= 0)
    inputerror(32L, LINK->system, WITH->UU.U0.nclass, LINK->LINK->LINK);
  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; k <= FORLIM2; k++) {
	fscanf(datafile, "%lf", &(*locus)->UU.U1.pm[j][k - 1][i]);
	WITH->UU.U1.pm[k][j - 1][i] = WITH->UU.U1.pm[j][k - 1][i];
      }
    }
  }
  fscanf(datafile, "%*[^\n]");
  getc(datafile);
  if (privelege && LINK->LINK->nupriv != lastpriv)
    return;
  FORLIM = WITH->UU.U1.ntrait;
  for (i = 0; i < FORLIM; i++) {
    FORLIM1 = WITH->UU.U1.ntrait;
    for (j = i; j < FORLIM1; j++) {
      fscanf(datafile, "%lf", &(*locus)->UU.U1.vmat[i][j]);
      if (i + 1 == j + 1 && WITH->UU.U1.vmat[i][j] <= 0.0)
	inputerror(33L, LINK->system, LINK->system, LINK->LINK->LINK);
      WITH->UU.U1.vmat[j][i] = WITH->UU.U1.vmat[i][j];
    }
  }
  fscanf(datafile, "%*[^\n]");
  getc(datafile);
  invert(WITH->UU.U1.vmat, WITH->UU.U1.ntrait, &WITH->UU.U1.det);
  WITH->UU.U1.det = 1 / sqrt(WITH->UU.U1.det);
  fscanf(datafile, "%lf%*[^\n]", &(*locus)->UU.U1.conmat);
  getc(datafile);
  if (WITH->UU.U1.conmat <= 0)
    inputerror(34L, LINK->system, LINK->system, LINK->LINK->LINK);
  WITH->UU.U1.conmat = 1 / WITH->UU.U1.conmat;
  WITH->UU.U1.contrait = 1.0;
  FORLIM = WITH->UU.U1.ntrait;
  for (i = 1; i <= FORLIM; i++)
    WITH->UU.U1.contrait *= WITH->UU.U1.conmat;
  WITH->UU.U1.contrait = sqrt(WITH->UU.U1.contrait);
}

Void getlocus(system_, LINK)
long system_;
struct LOC_readloci *LINK;
{
  struct LOC_getlocus V;
  long j;
  locusvalues *WITH, *WITH1;
  long FORLIM;
  int TEMP;

  V.LINK = LINK;
  V.system = system_;
  thislocus[V.system - 1] = (locusvalues *)Malloc(sizeof(locusvalues));
  if (thislocus[V.system - 1] == NULL)
    malloc_err("entry in thislocus");
  WITH = thislocus[V.system - 1];
  WITH->privlocus = NULL;
  fscanf(datafile, "%ld%ld", &LINK->whichtype,
	 &thislocus[V.system - 1]->nallele);
  
#if defined(GENERR)
  if (LINK->whichtype != 1)
    {
      printf("\nError:  Data must be in affection status format in order to perform \nerror checking.  The program <toaff> has been provided to convert any \ndata format to the affection status format.  The program will exit \npolitely to allow you to correct the problem.\n\n");
      exit(EXIT_FAILURE);
    }
#endif

  if (LINK->whichtype < 0 && LINK->whichtype > 4)
    inputerror(5L, V.system, LINK->whichtype, LINK->LINK);
  if (WITH->nallele > maxall)
    inputerror(6L, V.system, WITH->nallele, LINK->LINK);
  if (WITH->nallele <= 0)
    inputerror(7L, V.system, WITH->nallele, LINK->LINK);
  switch (LINK->whichtype) {

  case 0:
    WITH->which = quantitative;
    break;

  case 1:
    WITH->which = affection;
    break;

  case 2:
    WITH->which = binary_;
    break;

  case 3:
    WITH->which = binary_;
    break;
  }
  WITH->format = LINK->whichtype;
  if (lastpriv == 0) {
    fscanf(datafile, "%*[^\n]");
    getc(datafile);
  } else {
    fscanf(datafile, "%ld%*[^\n]", &LINK->whichtype);
    getc(datafile);
    if (LINK->whichtype == 0 || LINK->whichtype == 1) {
      LINK->nupriv++;
      WITH->privlocus = (locusvalues *)Malloc(sizeof(locusvalues));
      if (WITH->privlocus == NULL)
        malloc_err("privlocus field");
      WITH->privlocus->nallele = WITH->nallele;
      WITH1 = WITH->privlocus;
      switch (LINK->whichtype) {

      case 0:
	WITH1->which = quantitative;
	break;

      case 1:
	WITH1->which = affection;
	break;
      }
    }
  }
  if (!disequi) {
    FORLIM = WITH->nallele;
    for (j = 0; j < FORLIM; j++)
      fscanf(datafile, "%lf", &thislocus[V.system - 1]->freq[j]);
    fscanf(datafile, "%*[^\n]");
    getc(datafile);
  }
  switch (WITH->which) {

  case binary_:
    if (WITH->format == 3)
      getnumber(&thislocus[V.system - 1], &V.system, &V);
    else
      getbin(&thislocus[V.system - 1], &V.system, &V);
    break;

  case affection:
    getpen(&thislocus[V.system - 1], &V);
    break;

  case quantitative:
    getquan(&thislocus[V.system - 1], false, &V);
    break;
  }
  if (WITH->privlocus != NULL) {
    switch (WITH->privlocus->which) {

    case affection:
      getpen(&WITH->privlocus, &V);
      break;

    case quantitative:
      getquan(&WITH->privlocus, true, &V);
      break;
    }
  }
  if (LINK->nupriv == lastpriv && lastpriv != 0)
    lastpriv = V.system;
  if (!risk)
    return;
  if (V.system == risksys) {
    fscanf(datafile, "%d%*[^\n]", &TEMP);
    getc(datafile);
    riskall = TEMP;
  }
  if (riskall > thislocus[V.system - 1]->nallele)
    inputerror(35L, V.system, (int)riskall, LINK->LINK);
  if (riskall < 0)
    inputerror(36L, V.system, (int)riskall, LINK->LINK);
}


Void gettheta(sex_, LINK)
thetavalues **sex_;
struct LOC_readloci *LINK;
{
  thetarray oldtheta;
  long i;
  thetavalues *WITH;
  long FORLIM;

  *sex_ = (thetavalues *)Malloc(sizeof(thetavalues));
  nuneed = 7;
  for(i = 2; i < mlocus; i++)
    nuneed = 5 * nuneed - 3;
  if (*sex_ == NULL)
    malloc_err("item of type thetavalues");
   /*Next line added by A. A. Schaffer*/
  (*sex_)->segprob = (double*) malloc(nuneed * sizeof(double));
  if ((*sex_)->segprob == NULL)
    malloc_err("a segprob array in procedure gettheta");
  WITH = *sex_;
  if (*sex_ == maletheta || readfemale) {
    FORLIM = mlocus - 2;
    for (i = 0; i <= FORLIM; i++)
      fscanf(datafile, "%lf", &(*sex_)->theta[i]);
    if (interfer && !mapping)
      fscanf(datafile, "%lf", &(*sex_)->theta[mlocus - 1]);
    fscanf(datafile, "%*[^\n]");
    getc(datafile);
  } else {
    fscanf(datafile, "%lf%*[^\n]", &distratio);
    getc(datafile);
  }
  /*     FOR j:=1 TO maxneed DO segprob[j]:=0.0;*/
  if (!interfer || mapping)
    return;
  memcpy(oldtheta, WITH->theta, sizeof(thetarray));
  WITH->theta[0] = (oldtheta[0] + oldtheta[mlocus - 1] - oldtheta[mlocus - 2]) / 2.0;
  WITH->theta[mlocus - 2] =
    (oldtheta[mlocus - 2] + oldtheta[mlocus - 1] - oldtheta[0]) / 2.0;
  WITH->theta[mlocus - 1] =
    (oldtheta[0] + oldtheta[mlocus - 2] - oldtheta[mlocus - 1]) / 2.0;
  FORLIM = mlocus;
  for (i = 0; i < FORLIM; i++) {   /*=ln(1/0.0001-1.0)*/
    if (WITH->theta[i] > 0.0) {
      if (WITH->theta[i] < 0.999)
	WITH->theta[i] = log(1.0 / WITH->theta[i] - 1.0);
      else
	WITH->theta[i] = -6.9;   /*=ln(1/0.999-1.0)*/
    } else
      WITH->theta[i] = 9.21;
  }
}