/* This file contains the input routines in */
/* a modified version of the ILINK program*/
/* The modifications are described in the papers: */
/* R. W. Cotingham, 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 Genetic Linkage Analysis, */
/* Human Heredity 44(1994), pp. 225-237. */

#include "commondefs.h"
#include "gemdefs.h"
#include "ildefs.h"

Void setparam()
{
  long i, j, k, l, m;
  thetavalues *WITH;
  long FORLIM;
  locusvalues *WITH1;
  thisarray *WITH2;
  long FORLIM1, FORLIM2;

  penalty = false;
  k = 0;
  WITH = maletheta;
  FORLIM = mlocus - 2;
  for (i = 0; i <= FORLIM; i++) {
    k++;
    WITH->theta[i] = xall[k - 1];
  }
  if (interfer && !mapping) {
    k++;
    maletheta->theta[mlocus - 1] = xall[k - 1];
  }
  if (sexdif) {
    if (readfemale) {
      WITH = femaletheta;
      FORLIM = mlocus - 2;
      for (i = 0; i <= FORLIM; i++) {
	k++;
	WITH->theta[i] = xall[k - 1];
      }
      if (interfer && !mapping) {
	k++;
	femaletheta->theta[mlocus - 1] = xall[k - 1];
      }
    } else {
      k++;
      distratio = xall[k - 1];
    }
  }
  if (itsys == 0)
    return;
  WITH1 = thislocus[itsys - 1];
  if (disequi) {
    WITH2 = hapfreq;
    FORLIM = nuhap;
    for (i = 0; i < FORLIM; i++)
      WITH2->genarray[i] = xall[i + k];
    WITH2->genarray[nuhap - 1] = 1.0;
    FORLIM = nuhap - 2;
    for (i = 0; i <= FORLIM; i++)
      WITH2->genarray[nuhap - 1] -= WITH2->genarray[i];
    penalty = (WITH2->genarray[nuhap - 1] < 0.0);
    k += nuhap - 1;
  } else {
    FORLIM = WITH1->nallele - 2;
    for (i = 0; i <= FORLIM; i++)
      WITH1->freq[i] = xall[i + k];
    WITH1->freq[WITH1->nallele - 1] = 1.0;
    FORLIM = WITH1->nallele - 2;
    for (i = 0; i <= FORLIM; i++)
      WITH1->freq[WITH1->nallele - 1] -= WITH1->freq[i];
    penalty = (WITH1->freq[WITH1->nallele - 1] < 0.0);
    k += WITH1->nallele - 1;
  }
  if (WITH1->which == affection) {
    FORLIM = WITH1->UU.U0.nclass;
    for (l = 0; l < FORLIM; l++) {
      FORLIM1 = WITH1->nallele;
      for (i = 1; i <= FORLIM1; i++) {
	FORLIM2 = WITH1->nallele;
	for (j = i; j <= FORLIM2; j++) {
	  k++;
	  WITH1->UU.U0.pen[i][j - 1][2][l] = xall[k - 1];
	  WITH1->UU.U0.pen[i][j - 1][1][l] = 1.0 - xall[k - 1];
	  for (m = 0; m <= 2; m++)
	    WITH1->UU.U0.pen[j][i - 1][m][l] = WITH1->UU.U0.pen[i][j - 1][m]
	      [l];
	}
      }
      if (sexlink) {
	FORLIM1 = WITH1->nallele;
	for (i = 0; i < FORLIM1; i++) {
	  k++;
	  WITH1->UU.U0.pen[0][i][2][l] = xall[k - 1];
	  WITH1->UU.U0.pen[0][i][1][l] = 1.0 - xall[k - 1];
	}
      }
    }
    return;
  }
  if (WITH1->which != quantitative)
    return;
  FORLIM = WITH1->UU.U1.ntrait;
  /*USER MUST CHANGE FOR SYSTEMS WITH MORE THAN 2 ALLELES*/
  for (i = 0; i < FORLIM; i++) {
    WITH1->UU.U1.pm[1][0][i] = xall[k];
    WITH1->UU.U1.pm[2][1][i] = xall[k + 1] + xall[k];
    WITH1->UU.U1.pm[1][1][i] = xall[k + 2] * xall[k + 1] + xall[k];
    WITH1->UU.U1.pm[0][0][i] = xall[k];
    WITH1->UU.U1.pm[0][1][i] = WITH1->UU.U1.pm[1][1][i];
    WITH1->UU.U1.pm[2][0][i] = WITH1->UU.U1.pm[1][1][i];
    k += 3;
  }
  FORLIM = WITH1->UU.U1.ntrait;
  for (i = 0; i < FORLIM; i++) {
    FORLIM1 = WITH1->UU.U1.ntrait;
    for (j = i; j < FORLIM1; j++) {
      k++;
      WITH1->UU.U1.vmat[i][j] = xall[k - 1];
      WITH1->UU.U1.vmat[j][i] = xall[k - 1];
    }
  }
  invert(WITH1->UU.U1.vmat, WITH1->UU.U1.ntrait, &WITH1->UU.U1.det);
  WITH1->UU.U1.det = 1 / sqrt(WITH1->UU.U1.det);
}


typedef enum {
  nobnd, probnd, zerobnd, logbnd
} boundary;

Local Void setiterate(LINK)
struct LOC_inputdata *LINK;
{
  boundary bndtype[maxn];
  long i, j, k, l;
  thetavalues *WITH;
  long FORLIM;
  locusvalues *WITH1;
  long FORLIM1, FORLIM2;

  fscanf(datafile, "%ld%*[^\n]", &i);
  getc(datafile);
  if (i > mlocus)
    inputerror(23L, mlocus, i, LINK);
  if (i < 0)
    inputerror(24L, i, i, LINK);
  itsys = i;
  if (itsys != 0)
    itsys = order[itsys - 1];
  k = 0;
  WITH = maletheta;
  FORLIM = mlocus - 2;
  for (i = 0; i <= FORLIM; i++) {
    k++;
    if (k > maxn)
      inputerror(25L, (long)maxn, k, LINK);
    xall[k - 1] = WITH->theta[i];
    if (interfer && !mapping)
      bndtype[k - 1] = logbnd;
    else
      bndtype[k - 1] = probnd;
  }
  if (interfer && !mapping) {
    k++;
    if (k > maxn)
      inputerror(25L, (long)maxn, k, LINK);
    xall[k - 1] = maletheta->theta[mlocus - 1];
    bndtype[k - 1] = logbnd;
  }
  if (sexdif) {
    if (readfemale) {
      WITH = femaletheta;
      FORLIM = mlocus - 2;
      for (i = 0; i <= FORLIM; i++) {
	k++;
	if (k > maxn)
	  inputerror(25L, (long)maxn, k, LINK);
	xall[k - 1] = WITH->theta[i];
	if (interfer && !mapping)
	  bndtype[k - 1] = logbnd;
	else
	  bndtype[k - 1] = probnd;
      }
      if (interfer && !mapping) {
	k++;
	if (k > maxn)
	  inputerror(25L, (long)maxn, k, LINK);
	xall[k - 1] = femaletheta->theta[mlocus - 1];
	bndtype[k - 1] = logbnd;
      }
    } else {
      k++;
      if (k > maxn)
	inputerror(25L, (long)maxn, k, LINK);
      xall[k - 1] = distratio;
      bndtype[k - 1] = zerobnd;
    }
  }
  nall = k;
  if (itsys != 0) {
    WITH1 = thislocus[itsys - 1];
    if (disequi) {
      FORLIM = nuhap - 2;
      for (i = 0; i <= FORLIM; i++)
	xall[i + k] = hapfreq->genarray[i];
      nall = nuhap + k - 1;
    } else {
      FORLIM = WITH1->nallele - 2;
      for (i = 0; i <= FORLIM; i++)
	xall[i + k] = WITH1->freq[i];
      nall = WITH1->nallele + k - 1;
    }
    if (nall > maxn)
      inputerror(25L, (long)maxn, nall, LINK);
    FORLIM = nall;
    for (i = k; i < FORLIM; i++)
      bndtype[i] = probnd;
    k = nall;

    if (WITH1->which == quantitative) {
      nall += WITH1->UU.U1.ntrait * (WITH1->UU.U1.ntrait + 1) / 2 +
	      WITH1->UU.U1.ntrait * 3;

    } else {
      if (WITH1->which == affection) {
	nall += WITH1->UU.U0.nclass * WITH1->nallele * (WITH1->nallele + 1) / 2;
	if (sexlink)
	  nall += WITH1->nallele;
      }
    }
    if (nall > maxn)
      inputerror(25L, (long)maxn, nall, LINK);

    if (WITH1->which == affection) {
      FORLIM = WITH1->UU.U0.nclass;
      for (l = 0; l < FORLIM; l++) {
	FORLIM1 = WITH1->nallele;
	for (i = 1; i <= FORLIM1; i++) {
	  FORLIM2 = WITH1->nallele;
	  for (j = i - 1; j < FORLIM2; j++) {
	    k++;
	    if (k > maxn)
	      inputerror(25L, (long)maxn, k, LINK);
	    xall[k - 1] = WITH1->UU.U0.pen[i][j][2][l];
	    bndtype[k - 1] = probnd;
	  }
	}
	if (sexlink) {
	  FORLIM1 = WITH1->nallele;
	  for (i = 0; i < FORLIM1; i++) {
	    k++;
	    if (k > maxn)
	      inputerror(25L, (long)maxn, k, LINK);
	    xall[k - 1] = WITH1->UU.U0.pen[0][i][2][l];
	    bndtype[k - 1] = probnd;
	  }
	}
      }

    } else {
      if (WITH1->which == quantitative) {
	invert(WITH1->UU.U1.vmat, WITH1->UU.U1.ntrait, &WITH1->UU.U1.det);
	FORLIM = WITH1->UU.U1.ntrait;
	/*USER MUST MODIFY FOR A SYSTEM WILL MORE THAN 2 ALLELES */
	for (l = 0; l < FORLIM; l++) {
	  FORLIM1 = WITH1->nallele;
	  for (i = 1; i <= FORLIM1; i++) {
	    FORLIM2 = WITH1->nallele;
	    for (j = i; j <= FORLIM2; j++) {
	      k++;
	      if (i == 1 && j == 1)
		xall[k - 1] = WITH1->UU.U1.pm[i][j - 1][l];
	      else {
		if (i == 1)
		  xall[k - 1] = WITH1->UU.U1.pm[j][j - 1][l] - xall[k - 2];
		else
		  xall[k - 1] = (WITH1->UU.U1.pm[i - 1][j - 1]
				 [l] - xall[k - 3]) / xall[k - 2];
	      }
	      bndtype[k - 1] = nobnd;
	    }
	  }
	}
	FORLIM = WITH1->UU.U1.ntrait;
	for (i = 1; i <= FORLIM; i++) {
	  FORLIM1 = WITH1->UU.U1.ntrait;
	  for (j = i; j <= FORLIM1; j++) {
	    k++;
	    if (k > maxn)
	      inputerror(25L, (long)maxn, k, LINK);
	    xall[k - 1] = WITH1->UU.U1.vmat[i - 1][j - 1];
	    if (i == j)
	      bndtype[k - 1] = zerobnd;
	    else
	      bndtype[k - 1] = nobnd;
	  }
	}
      }
    }

  }
  FORLIM = nall;
  for (i = 0; i < FORLIM; i++)
    itp[i] = 0;
  i = 0;
  while ((i <= nall) & (!P_eoln(datafile))) {
    i++;
    fscanf(datafile, "%d", &itp[i - 1]);
  }
  n = 0;
  FORLIM = nall;
  for (i = 0; i < FORLIM; i++) {
    if (itp[i] == 1) {
      n++;
      x[n - 1] = xall[i];
    }
  }
  FORLIM = nall;
  for (i = 0; i < FORLIM; i++) {
    bnd[i][0] = -1000.0;
    bnd[i][1] = 1000.0;
    if (bndtype[i] == probnd) {
      bnd[i][0] = 0.0;
      bnd[i][1] = 1.0;
    } else {
      if (bndtype[i] == logbnd) {
	bnd[i][0] = -20.0;
	bnd[i][1] = 20.0;
      } else {
	if (bndtype[i] == zerobnd)
	  bnd[i][0] = 0.0;
      }
    }
  }
}

Local Void readloci(LINK)
struct LOC_inputdata *LINK;
{
  struct LOC_readloci V;
  long i, j, coupling, autosomal, independent, difference, FORLIM;
  locusvalues *WITH;


  V.LINK = LINK;
  lastpriv = 0;
  fscanf(datafile, "%ld%ld%ld%*[^\n]", &mlocus, &risksys, &autosomal);
  getc(datafile);

#if defined(GENERR)
  if (mlocus > maxlocus)
    {
      printf("\nError:  The value of maxlocus is set to %d.  If you wish to analyze \nmore loci, please change the value of maxlocus.  The program will \nexit politely to allow you to correct the problem.\n\n",maxlocus);
      exit(EXIT_FAILURE);
    }
#endif

  /*Replace the above line with the next when using epistasis*/
  /*readln(datafile,mlocus,risksys,autosomal,lastpriv);*/
  if (mlocus > maxlocus)
    inputerror(0L, mlocus, mlocus, LINK);
  if (mlocus <= 0)
    inputerror(1L, mlocus, mlocus, LINK);
  if (risksys > maxlocus)
    inputerror(37L, risksys, risksys, LINK);
  if (risksys < 0)
    inputerror(38L, risksys, risksys, LINK);
  risk = (risksys != 0);
  sexlink = (autosomal == 1);
  printf("YOU ARE USING LINKAGE (V%3.1f (1-Feb-1991)) WITH%3ld-POINT\n",
	 version, mlocus);
  printf("YOU ARE USING FASTLINK (V%3.1f (1-Sep-1994))", fastversion);
  if (sexlink)
    printf(" SEXLINKED DATA\n");
  else
    printf(" AUTOSOMAL DATA\n");

#if defined(GENERR)
  printf("YOU ARE USING GenoCheck (V1.0 (20-Sep-1995)).\n");
#endif

  fscanf(datafile, "%ld%lf%lf%ld%*[^\n]", &mutsys, &mutmale, &mutfemale,
	 &coupling);
  getc(datafile);
  if (mutsys > maxlocus)
    inputerror(39L, mutsys, mutsys, LINK);
  if (mutsys < 0)
    inputerror(40L, mutsys, mutsys, LINK);
  if (coupling == 1)
    disequi = true;
  else
    disequi = false;
  if (disequi) {
    hapfreq = (thisarray *)Malloc(sizeof(thisarray));
    if (hapfreq == NULL)
      malloc_err("hapfreq");
  }
  FORLIM = mlocus;
  for (i = 1; i <= FORLIM; i++) {
    fscanf(datafile, "%ld", &j);
    if (j > mlocus)
      inputerror(2L, i, j, LINK);
    if (j <= 0)
      inputerror(3L, i, j, LINK);
    order[j - 1] = i;
  }
  FORLIM = mlocus;
  for (i = 1; i <= FORLIM; i++) {
    for (j = 1; j < i; j++) {
      if (order[i - 1] == order[j - 1])
	inputerror(4L, i, j, LINK);
    }
  }
  fscanf(datafile, "%*[^\n]");
  getc(datafile);
  if (mutsys != 0)
    mutsys = order[mutsys - 1];
  if (risksys != 0)
    risksys = order[risksys - 1];
  V.nupriv = 0;
  FORLIM = mlocus;
  for (i = 0; i < FORLIM; i++)
    getlocus(order[i], &V);
  increment[mlocus - 1] = 1;
  for (i = mlocus - 1; i >= 1; i--)
    increment[i - 1] = increment[i] * thislocus[i]->nallele;
  fgeno = 1;
  FORLIM = mlocus;
  for (j = 0; j < FORLIM; j++)
    fgeno *= thislocus[j]->nallele;
  mgeno = fgeno;
  nuhap = fgeno;
  FORLIM = mlocus;
  for (i = 0; i < FORLIM; i++)
    nohom[i] = false;
  if (disequi) {
    FORLIM = mgeno;
    for (i = 0; i < FORLIM; i++)
      fscanf(datafile, "%lf", &hapfreq->genarray[i]);
    fscanf(datafile, "%*[^\n]");
    getc(datafile);
  } else {
    FORLIM = mlocus;
    for (i = 0; i < FORLIM; i++) {
      WITH = thislocus[i];
      if (WITH->which == affection || WITH->which == quantitative) {
	if (WITH->freq[affall - 1] < minfreq)
	  nohom[i] = true;
      }
    }
  }
  fgeno = fgeno * (fgeno + 1) / 2;
  if (!sexlink)
    mgeno = fgeno;
  fscanf(datafile, "%ld", &difference);
  if ((unsigned long)difference > 2) {
    inputwarning(0L, difference, difference, LINK);
    difference = 0;
  }
  sexdif = (difference != 0);
  readfemale = (difference == 2);
  fscanf(datafile, "%ld%*[^\n]", &independent);
  getc(datafile);
  if ((unsigned long)independent > 2) {
    inputwarning(1L, independent, independent, LINK);
    independent = 0;
  }
  interfer = (independent != 0);
  mapping = (independent == 2);
  gettheta(&maletheta, &V);
  if (sexdif)
    gettheta(&femaletheta, &V);
  else
    femaletheta = maletheta;
  if (sexlink && sexdif) {
    inputwarning(2L, difference, difference, LINK);
    sexdif = false;
    readfemale = false;
  }
  if (!sexlink) {
    if (mutsys == 0)
      thispath = auto_;
    else
      thispath = mauto;
  } else if (mutsys == 0)
    thispath = sex;
  else
    thispath = msex;
  segscale = scale;
  FORLIM = mlocus;
  for (i = 1; i <= FORLIM; i++)
    segscale *= scalemult;
}


Void inputdata()
{
  struct LOC_inputdata V;

  readloci(&V);
  setiterate(&V);
  readped(&V);
  readspeed(&V);
}