#include <stdio.h>
#include <math.h>
#define TRUE 1
#define FALSE 0

void weights();
void gamma_find();
double nblik();
int inverse(double matrix[3][3]);
double f_prime(double gamma);
double dlgama(double z);
double ztail(double z);




int numdos, mumdos;
double xbar[11];
float dose[11];
int kount[11][11];
char backgr[11][11];
int n[11], nn[11];
double weight[11];
double muhat[11];
double adtvcnx;
double adtvcn;
double alpha;
double beta;
double gamma;
double c;
int model;
float outparm[30];
double stats[11];
double dose_last;
double xmut;

FILE *fopen(), *ip, *op;


main ()

{

char line[80], ch;
char raw_backgr[11][11], reject_dose[11], reject_kount[11][11];
float raw_dose[11];
int raw_kount[11][11];
int raw_n[11];
long kntsum;
int i, j, k, l, imany, jstat, numtox;
double xsum, dsum, ddsum, dxsum, xxbar;
double c_all, w, cc[11], x, x2, cnew, z[11];
int nsum, jmax;
int model2_flag, finish_up, nodel;
double fnew;
int ioffset, ierror, dataleft;
#define GCUT 9.5
int peak_dose, ndf, isign;
double peaktest, tvalue, sigma1, test;
double covar[3][3];
static double peak[8] = {0.0, 2.06, 2.16, 2.23, 2.29, 2.34, 2.38, 2.42};
static double sigmul[23] = {0.000, 2.539, 2.576, 2.609, 2.638, 2.665, 2.690,
                            2.713, 2.734, 2.754, 2.773, 2.791, 2.807, 2.823,
                            2.838, 2.852, 2.865, 2.878, 2.891, 2.902, 2.914,
                            2.925, 2.935 };

void get_initial();
double nblik();
void iterate();
void weights();
void parms();
void gamma_find();
double f_prime();
double addconx();
double addcon();
double cmputc();
double dlgama();
void output();
double chicdf();
double ztail();
void fisher();
void muhatc();
int inverse();

ip = fopen("salm.in","r");
if (!ip) {
  printf("Error on salm.in");
  exit (1);
}

op = fopen("salm.out","w");
if (!op) {
  printf("Error on salm.out");
  exit (1);
}
/*
fprintf(stdout,"\7\7\7"
               "                          #############################\n"
               "                          ####                     ####\n"
               "                          #### Salmonella Analysis ####\n"
               "                          #### á release  ver. 1.1 ####\n"
               "                          ####  November 11, 1990  ####\n"
               "                          ####                     ####\n"
               "                          #### Report any problems ####\n"
               "                          #### to Kenneth J. Risko ####\n"
               "                          ####                     ####\n"
               "                          #############################\n"
               "\7\7\7");

fprintf(op,"                          #############################\n"
           "                          ####                     ####\n"
           "                          #### Salmonella Analysis ####\n"
           "                          #### á release  ver. 1.1 ####\n"
           "                          ####  November 11, 1990  ####\n"
           "                          ####                     ####\n"
           "                          #### Report any problems ####\n"
           "                          #### to Kenneth J. Risko ####\n"
           "                          ####                     ####\n"
           "                          #############################\n");


fprintf(stdout,"Press any key to continue");
fscanf(stdin,"%c",&ch);
*/

while(fgets(line,80,ip)) {

  fprintf(stdout,"\n\n");
  fprintf(op,"\n\n");

  i = 0;
  while((ch = line[i++]) != '\n') {
    fprintf(stdout,"%c",ch);
    fprintf(op,"%c",ch);
  }
  fprintf(stdout,"%c",ch);
  fprintf(op,"%c",ch);

  fscanf(ip,"%d",&numdos);
  fprintf(op,"number of doses = %2d\n",numdos);
  fprintf(stdout,"number of doses = %2d\n",numdos);
  for(i=1;i<=numdos;i++) {
    fscanf(ip,"%f,%d",&raw_dose[i],&raw_n[i]);
    fprintf(op,"%2d,%7.2f,%2d",i,raw_dose[i],raw_n[i]);
    fprintf(stdout,"%2d,%7.2f,%2d",i,raw_dose[i],raw_n[i]);
    reject_dose[i] = ' ';
    for(j=1;j<=raw_n[i];j++) {
      fscanf(ip,",%d,%c",&raw_kount[i][j],&raw_backgr[i][j]);
      fprintf(op,",%4d(%c)",raw_kount[i][j],raw_backgr[i][j]);
      fprintf(stdout,",%4d(%c)",raw_kount[i][j],raw_backgr[i][j]);
      reject_kount[i][j] = ' ';
    }
    fprintf(op,"\n");
    fprintf(stdout,"\n");
  }

/***** CLEANSE DATA *****/

  for(i=1;i<=numdos;i++) {
    for(j=1;j<=raw_n[i];j++) {
      if(raw_backgr[i][j] != ' ' &&
         raw_backgr[i][j] != 'P' &&
         raw_backgr[i][j] != 'T' &&
         raw_backgr[i][j] != 'S'   ) {
           reject_kount[i][j] = raw_backgr[i][j];
           fprintf(op,"Rejecting plate count %4d at dose %7.2f because "
             "background code = %c\n",raw_kount[i][j],raw_dose[i],
             raw_backgr[i][j]);
           fprintf(stdout,"Rejecting plate count %4d at dose %7.2f because "
             "background code = %c\n",raw_kount[i][j],raw_dose[i],
             raw_backgr[i][j]);
      }
    }
  }

  fprintf(op,"\n");
  fprintf(stdout,"\n");

  for(i=1;i<=numdos;i++) {
    numtox = 0;
    for(j=1;j<=raw_n[i];j++) {
      if(raw_backgr[i][j] == 'T') {
        ++numtox;
        reject_kount[i][j] = 'T';
        fprintf(op,"Rejecting plate count %4d at dose %7.2f because "
          "background code = %c\n",raw_kount[i][j],raw_dose[i],raw_backgr[i][j]);
        fprintf(stdout,"Rejecting plate count %4d at dose %7.2f because "
          "background code = %c\n",raw_kount[i][j],raw_dose[i],raw_backgr[i][j]);
      }
    }
    if(3*numtox > raw_n[i]) {
      for(j=i;j<=numdos;j++) {
        reject_dose[j] = 'T';
        fprintf(op,"Rejecting dose %7.2f because of toxicity\n",raw_dose[j]);
        fprintf(stdout,"Rejecting dose %7.2f because of toxicity\n",raw_dose[j]);
      }
      break;
    }
  }

  numdos = i - 1;

  for(i=1;i<=numdos;i++) {
    kntsum = 0;
    dataleft = FALSE;
    for(j=1;j<=raw_n[i];j++) {
      if(reject_kount[i][j] == ' ') {
        dataleft = TRUE;
        kntsum += raw_kount[i][j];
      }
    }
    if(kntsum == 0 && dataleft) {
      for(j=i;j<=numdos;j++) {
        reject_dose[j] = '0';
        fprintf(op,"Rejecting dose %7.2f because of zero mean "
          "at or before dose %6f\n",raw_dose[j],raw_dose[j]);
        fprintf(stdout,"Rejecting dose %7.2f because of zero mean "
          "at or before dose %6f\n",raw_dose[j],raw_dose[j]);
      }
      break;
    }
  }

  numdos = i - 1;

  if(numdos == 0) return;

  for(i=numdos;i>=1;i--) {
    dose[i] = raw_dose[i];
    n[i] = 0;
    for(j=1;j<=raw_n[i];j++) {
      if(reject_kount[i][j] == ' ') {
        n[i] += 1;
        kount[i][n[i]] = raw_kount[i][j];
        backgr[i][n[i]] = raw_backgr[i][j];
      }
    }
    if(n[i] == 0) {
      if(i != numdos) {
        for(k=i;k<=numdos-1;k++) {
          dose[k] = dose[k+1];
          n[k] = n[k+1];
          for(l=1;l<=n[k];l++) {
            kount[k][l] = kount[k+1][l];
            backgr[k][l] = backgr[k+1][l];
          }

        }
      }
    --numdos;
    }
  }

/***** REMOVE OUTLIERS *****/

  c_all = 0.0;
  w = 0.0;
  nsum = 0;
  for(i=1;i<=numdos;i++) {
    nsum += n[i];
    cc[i] = 0.0;
    xbar[i] = 0.0;
    for(j=1;j<=n[i];j++) xbar[i] += kount[i][j];
    xbar[i] /= n[i];
    if(xbar[i] != 0.0 && n[i] > 1) {
      x = xbar[i];
      x2 = x*x;
      w = w + x2;
      for(j=1;j<=n[i];j++) cc[i] += kount[i][j]*((double) kount[i][j] - 1.0)/x2;
      cc[i] = (cc[i] + 1.0/x - n[i])/(n[i] - 1);
      c_all += cc[i]*x2;
    }
  }
  c_all /= w;

  for(i=1;i<=numdos;i++) {
    if(n[i] >= 3) {
      cnew = (w*c_all - xbar[i]*xbar[i]*cc[i])/(w - xbar[i]*xbar[i]);
      if(cnew < 0.0) cnew = 0.0;
      jmax = 1;
      for(j=1;j<=n[i];j++) {
        z[j] = (xbar[i] == 0.0) ? 0.0 : fabs(kount[i][j] - xbar[i])/
                                  sqrt(2.0*xbar[i]*(1.0 + cnew*xbar[i])/3.0);
        if(z[j] > z[jmax]) jmax = j;
      }
      if(z[jmax] >= sigmul[nsum-8]) {
        fprintf(op,"Plate count %4d at dose %7.2f determined to "
          "an outlier\n",kount[i][jmax],dose[i]);
        fprintf(stdout,"Plate count %4d at dose %7.2f determined to "
          "an outlier\n",kount[i][jmax],dose[i]);
        --n[i];
        for(k=jmax;k<=n[i];k++) {
          kount[i][k] = kount[i][k+1];
          backgr[i][k] = backgr[i][k+1];
        }
      }
    }
  }

  for(i=1;i<=numdos;i++) {
	xbar[i] = 0.0;
	for(j=1;j<=n[i];j++) xbar[i] += kount[i][j];
	xbar[i] = xbar[i]/n[i];
    if(xbar[i] == 0.0) {
      for(j=i;j<=numdos;j++) {
        fprintf(op,"Rejecting dose %7.2f because of zero mean (after outlier"
          " deletion) at or before dose %6f\n",raw_dose[j],raw_dose[j]);
        fprintf(stdout,"Rejecting dose %7.2f because of zero mean (after outlier"
          " deletion) at or before dose %6f\n",raw_dose[j],raw_dose[j]);
      }
      break;
    }
  }

  numdos = i - 1;

  if(numdos < 4) {
    fprintf(stdout,"Too few remaining doses (%2d) to perform analysis\n",numdos);
    return(-1);
  }

  imany = 0;
  for(j=1;j<=numdos;j++) imany += (n[j] >= 2) ? 1 : 0;
  if(imany < 3) {
	fprintf(stdout,"Too few doses with replicates (%2d) to perform analysis\n",imany);
	return(-2);
  }

  xsum  = 0.0;
  nsum  = 0;
  dsum  = 0.0;
  ddsum = 0.0;
  dxsum = 0.0;
  for(j=1;j<=numdos;j++) {
	nsum  += n[j];
	dsum  += n[j]*dose[j];
	ddsum += n[j]*dose[j]*dose[j];
	dxsum += n[j]*dose[j]*xbar[j];
	xsum  += n[j]*xbar[j];
  }
  xxbar = xsum/nsum;

  c = cmputc();

  peak_dose = 1;
  for(i=2;i<=numdos;i++) peak_dose = (xbar[i] >= xbar[peak_dose]) ?
                                     i : peak_dose;
  stats[1] = -peak_dose;
  for(i=peak_dose+1;i<=numdos;i++) {
    peaktest = (xbar[peak_dose] - xbar[i])/
               sqrt(xbar[peak_dose]*(1.0 + c*xbar[peak_dose])/n[peak_dose] +
                    xbar[i]*(1.0 + c*xbar[i])/n[i]);
    if(peaktest >= 2.576) {
      for(j=i+1;j<=numdos;j++) {
        fprintf(op,"Rejecting dose %2d because depressed mean response of "
          "%6.2f is substantially below peak response of %6.2f\n",j,xbar[j],
          xbar[peak_dose]);
        fprintf(stdout,"Rejecting dose %2d because depressed mean response of "
          "%6.2f is substantially below peak response of %6.2f\n",j,xbar[j],
          xbar[peak_dose]);
      }
      break;
    }
  }
  numdos = (i < numdos) ? i : numdos;

  if(numdos <= 3) {
    fprintf(stdout,"Too few remaining doses (%2d) to perform analysis\n",numdos);
    return(-3);
  }

  fprintf(op,"\n *** ANALYSED DATA ***\n\n");
  fprintf(op," No.   Dose                 Plates\n"
                 " ---   ----             --------------\n");
  fprintf(stdout,"\n *** ANALYSED DATA ***\n\n");
  fprintf(stdout," No.   Dose                 Plates\n"
                 " ---   ----             --------------\n");
  for(i=1;i<=numdos;i++) {
    fprintf(op," %2d  %7.2f  ",i,dose[i]);
    fprintf(stdout," %2d  %7.2f  ",i,dose[i]);
    for(j=1;j<=n[i];j++) fprintf(op,"%5d (%c) ",kount[i][j],backgr[i][j]);
    for(j=1;j<=n[i];j++) fprintf(stdout,"%5d (%c) ",kount[i][j],backgr[i][j]);
    fprintf(op,"\n");
    fprintf(stdout,"\n");
  }

  fprintf(stdout,"Running");

  for(i=1;i<=26;i++) outparm[i] = -1.0;

  c = cmputc();

  dose_last = dose[numdos];
  adtvcnx = addconx();
  adtvcn = addcon(xsum,nsum);
  for(i=1;i<=numdos;i++) muhat[i] = xxbar;
  outparm[13] = xxbar;
  outparm[14] = 0.0;
  outparm[15] = 0.0;
  outparm[16] = nblik() - adtvcn;	

  for(i=1;i<=numdos;i++) muhat[i] = xbar[i];
  outparm[25] = c;
  outparm[26] = nblik() - adtvcn;

  model2_flag = FALSE;
  model = 0;
  weights();
  gamma = 0.0;
  numdos -= 1;
  parms();
  numdos += 1;
  xmut = alpha + beta*dose_last;
  gamma = -log(xbar[numdos]/xmut)/dose_last;
  if(gamma < 0.0) gamma = 0.0;
  parms();
  weights();
  finish_up = FALSE;

  while(!finish_up) {
    iterate(&fnew);
    ioffset = 4*model;
    outparm[ioffset+1] = alpha;
    outparm[ioffset+2] = beta;
    outparm[ioffset+3] = gamma;
    outparm[ioffset+4] = fnew;
    if(model == 5) {
      finish_up = TRUE;
    } else {
      if(model == 4) {
        model = 1;
        get_initial(model);
      } else {
        if(gamma < 0.0) {
          model2_flag = TRUE;
          for(model=0;model<=2;model++) {
	        ioffset = 4*model;
	        outparm[ioffset+1] = alpha;
	        outparm[ioffset+2] = beta;
	        outparm[ioffset+3] = 0.0;
	        outparm[ioffset+4] = fnew;
          }
          for(model=4;model<=5;model++) {
	        ioffset = 4*model;
	        for(jstat=1;jstat<=4;jstat++) outparm[ioffset+jstat] = outparm[12+jstat];
          }
          output();
          return;
        }
        model = model + 4;
        if(beta <= 0.0 || outparm[13] <= (dxsum/dsum)) {
          ioffset = 4*model;
          outparm[ioffset+1] = (beta <= 0.0) ? alpha : outparm[13];
          outparm[ioffset+2] = 0.0;
          outparm[ioffset+3] = (beta <= 0.0) ? gamma : 0.0;
          outparm[ioffset+4] = (beta <= 0.0) ? fnew : outparm[16];
          if(model == 5) {
            finish_up = TRUE;
          } else {
            model = 1;
            get_initial(model);
          }
        } else {
          if(model == 5) {
            get_initial(model);
          } else {
            double ysum = 0.0;
            double dysum = 0.0;
            double y;
            double determ;
            for(j=1;j<=numdos;j++) {
	          for(i=1;i<=n[j];i++) {
	            y = kount[j][i];
	            if(y < 1.0) y = 1.0;
	            y = log(y);
	            ysum += y;
	            dysum += dose[j]*y;
              }
            }
            determ = nsum*ddsum - dsum*dsum;
            gamma = (dsum*ysum - nsum*dysum)/determ;
            if(gamma < 0.0) gamma = 0.0;
            alpha = exp((ysum + gamma*dsum)/nsum);
            beta = 0.0;
            for(j=1;j<=numdos;j++) muhat[j] = alpha*exp(-gamma*dose[j]);
          }
          weights();
        }
      }
    }
  }

  if(!model2_flag) {
    if(outparm[13] >= (dxsum/dsum)) {
      outparm[9] = outparm[13];
      outparm[10] = 0.0;
      outparm[11] = 0.0;
      outparm[12] = outparm[16];
    } else {
      double determ;
      determ = nsum*ddsum - dsum*dsum;
      beta = (nsum*dxsum - dsum*xsum)/determ;
      alpha = (xsum - beta*dsum)/nsum;
      gamma = 0.0;
      for(j=1;j<=numdos;j++) muhat[j] = alpha + beta*dose[j];
      weights();
      model = 2;
      iterate(&fnew);
      outparm[9] = alpha;
      outparm[10] = beta;
      outparm[11] = 0.0;
      outparm[12] = fnew;
    }
  }

  fprintf(stdout,"\n\n");
  c = outparm[25];
  model = (outparm[4] < outparm[8]) ? 0 : 1;
  ioffset = 4*model;
  if(outparm[12] - outparm[ioffset+4] < 0.5*GCUT) {
    model = 2;
    tvalue = 2.0*(outparm[16] - outparm[12]);
    tvalue = (tvalue > 0.0) ? sqrt(tvalue) : 0.0;
  } else {
    tvalue = 2.0*(outparm[ioffset+20] - outparm[ioffset+4]);
    tvalue = (tvalue > 0.0) ? sqrt(tvalue) : 0.0;
    sigma1 = xbar[1]*(1.0 + c*xbar[1])/n[1];
    isign = -1;
    for(i=2;i<=numdos;i++) {
      test = sqrt((1.0 + (double) n[1]/n[i])*sigma1)*peak[numdos-3] + xbar[1];
      if(xbar[i] > test) isign = 1;
    }
    tvalue = isign*tvalue;
  }
  stats[1] = model;
  ioffset = 4*model;
  stats[2] = outparm[ioffset+1];
  stats[3] = outparm[ioffset+2];
  stats[4] = outparm[ioffset+3];
  stats[5] = c;
  stats[6] = outparm[ioffset+4];
  stats[7] = 2.0*(outparm[ioffset+4] - outparm[26]);
  ndf = numdos - ((model ==2) ? 2 : 3);
  stats[8] = chicdf(stats[7],ndf);
  stats[9] = tvalue;
  stats[10] = ztail(fabs(stats[9]));

  for(nodel=0;nodel<=2;nodel++) {
    model = nodel;
    if(outparm[4*model+3] == 0.0) model = 2;
    ioffset = 4*model;
    alpha = outparm[ioffset+1];
    beta  = outparm[ioffset+2];
    gamma = outparm[ioffset+3];
    muhatc();
    weights();
    fisher(covar);
    ierror = inverse(covar);
    tvalue = 0.0;
    if(!ierror) tvalue = beta/sqrt(covar[1][1]);
    outparm[27+nodel] = tvalue;
  }

  output();

  fgets(line,80,ip);

}

return;

}


void output()

{

int i, j, model, maxn;
double fit, xbar_max;

  model = stats[1];
  switch (model) {
    case 0:
      fprintf(stdout,"\n\nFitted model = (%6f + %6f x dose) x "
                     "exp( -%6f x dose)\n",stats[2],stats[3],stats[4]);
      fprintf(op,"\n\nFitted model = (%6f + %6f x dose) x "
                     "exp( -%6f x dose)\n",stats[2],stats[3],stats[4]);
      break;
    case 1:
      fprintf(stdout,"\n\nFitted model = (%6f + %6f x dose) x "
                     "(2 - exp( %6f x dose))\n",stats[2],stats[3],stats[4]);
      fprintf(op,"\n\nFitted model = (%6f + %6f x dose) x "
                     "(2 - exp( %6f x dose))\n",stats[2],stats[3],stats[4]);
      break;
    case 2:
      fprintf(stdout,"\n\nFitted model = (%6f + %6f x dose)\n",
                      stats[2],stats[3]);
      fprintf(op,"\n\nFitted model = (%6f + %6f x dose)\n",
                      stats[2],stats[3]);
      break;
    default:
      break;
  }

  fprintf(stdout,"\n *** ANALYSED DATA ***\n\n");
  fprintf(stdout," No.   Dose          Plates"
                 "         Mean        Model\n"
                 " ---   ----      --------------"
                 "   --------    ---------\n");
  fprintf(op,"\n *** ANALYSED DATA ***\n\n");
  fprintf(op," No.   Dose          Plates"
                 "         Mean        Model\n"
                 " ---   ----      --------------"
                 "   --------    ---------\n");

  xbar_max = 0.0;
  maxn = 0;
  for(i=1;i<=numdos;i++) {
    xbar_max = (xbar[i] > xbar_max) ? xbar[i] : xbar_max;
    maxn = (n[i] > maxn) ? n[i] : maxn;
    switch (model) {
      case 0:
        fit = (stats[2] + stats[3]*dose[i])*exp(-stats[4]*dose[i]);
        break;
      case 1:
        fit = (stats[2] + stats[3]*dose[i])*(2-exp(stats[4]*dose[i]));
        break;
      case 2:
        fit = (stats[2] + stats[3]*dose[i]);
        break;
      default:
        fit = 0.0;
    }
    fprintf(op," %2d  %7.2f  ",i,dose[i]);
    fprintf(stdout," %2d  %7.2f  ",i,dose[i]);
    for(j=1;j<=n[i];j++) fprintf(op,"%5d ",kount[i][j]);
    for(j=1;j<=n[i];j++) fprintf(stdout,"%5d ",kount[i][j]);
    for(j=n[i]+1;j<=maxn;j++) fprintf(op,"      ");
    for(j=n[i]+1;j<=maxn;j++) fprintf(stdout,"      ");
    fprintf(op,"  %6.1f       %6.1f\n",xbar[i],fit);
    fprintf(stdout,"  %6.1f       %6.1f\n",xbar[i],fit);
  }

  fprintf(stdout,"\nc = %6.2f\n",stats[5]);
  fprintf(op,"\nc = %6.2f\n",stats[5]);

  fprintf(stdout,"\nTest for GOODNESS OF FIT p-value: %5.3f\n",stats[8]);
  fprintf(op,"\nTest for GOODNESS OF FIT p-value: %5.3f\n",stats[8]);

  fprintf(stdout,"\nTest for MUTAGENICITY p-value: %5.3f\n",stats[10]);
  fprintf(op,"\nTest for MUTAGENICITY p-value: %5.3f\n",stats[10]);

  if(stats[9] < 0.0) {
    fprintf(stdout,"\n***  NOTE: Peak mean response of %6.2f not "
                   "significanltly elevated above control mean response "
                   " of %6.2f  ***\n",xbar_max,xbar[1]);
    fprintf(op,"\n***  NOTE: Peak mean response of %6.2f not "
               "significanltly elevated above control mean response "
               " of %6.2f\n  ***",xbar_max,xbar[1]);
  }

}


void get_initial(model)

int model;

{

int j, ioffset;

  ioffset = 4*model;
  alpha = outparm[ioffset - 3];
  beta = outparm[ioffset - 2];
  gamma = log(2.0 - exp(-outparm[ioffset - 1]*dose_last))/dose_last;
  for(j=1;j<=numdos;j++) muhat[j] = (alpha + beta*dose[j])*
                                        (2.0 - exp(gamma*dose[j]));
  weights();
  return;

}


void iterate(fnew)

double *fnew;

{

  int loop;
  double fold;

  gamma_find();
  *fnew = nblik() - adtvcn;

/*  iterate and reweight  */

  for(loop=1;loop<=10;loop++) {
    fprintf(stdout,".");
	weights();
	fold = *fnew;
 	gamma_find();
 	*fnew = nblik() - adtvcn;
 	if(fabs(fold-*fnew) <= 1.0e-4) break;
  }
  return;

}


void parms()

{

  double sumwdx, sumwd, sumwx, sumw, sumwdd;
  double tox[11], tj, wt2;
  int j;

  sumwdx = 0.0;
  sumwd = 0.0;
  sumwx = 0.0;
  sumw = 0.0;
  sumwdd = 0.0;
  for(j=1;j<=numdos;j++) {
    tox[j] = (model != 0 && model != 4) ? 2.0 - exp(gamma*dose[j])
                                        : exp(-gamma*dose[j]);
    tj = tox[j];
    wt2 = weight[j]*tj*tj;
    sumwx  += wt2         *(xbar[j]/tj);
    sumw   += wt2;
    if(model < 4) {
 	  sumwdx += wt2*dose[j]*(xbar[j]/tj);
	  sumwd  += wt2*dose[j];
	  sumwdd += wt2*dose[j]*dose[j];
	}
  }
  beta = 0.0;
  if(model < 4) {
	beta = (sumwdx - sumwd*sumwx/sumw)/(sumwdd - sumwd*sumwd/sumw);
	if(beta < 0.0 || model >= 3) beta = 0.0;
  }
  alpha = (sumwx - beta*sumwd)/sumw;
  if(alpha <= 0.0) {
	alpha = 1.0;
	beta = (sumwdx - sumwd)/sumwdd;
  }
  for(j=1;j<=numdos;j++) muhat[j] = (alpha + beta*dose[j])*tox[j];

  return;

}


void gamma_find()

{

static double gamma_int[2][10] =
     { {6.908, 2.303, 1.609, 1.204, 0.916, 0.693, 0.511, 0.357, 0.223, 0.105 },
       {0.693, 0.642, 0.588, 0.531, 0.470, 0.405, 0.336, 0.262, 0.182, 0.095 } };
int mod_type;
double tox_at;
int q1, q2, q3, did_break;
int istep, igam, loop;
double gamma_low, gamma_high, gamma_at, g_diff, g2, g1;
double g1lik, g2lik;
int ig0, ig1, ig2, ig3;

  dose_last = dose[numdos];
  mod_type = 0;
  if(model == 1 || model == 5) mod_type = 1;

/*  check for zero gamma and/or model 2  */

  if(f_prime(0.0) >= 0.0 || model == 2) {
	gamma = 0.0;
	parms;
	return;
  }

  tox_at = exp(-gamma*dose_last);
  if(mod_type == 1) tox_at = 2.0 - 1.0/tox_at;
  q1 = (f_prime(gamma) <= 0.0);
  q2 = (gamma*dose_last >= gamma_int[mod_type][0]);
  istep = (f_prime(gamma) <= 0.0) ? -1 : 1;
  if(q1 && q2) {
	gamma = gamma_int[mod_type][0]/dose_last;
	parms();
	return;
  }
  gamma_low = 0.0;
  gamma_high = gamma;
  q3 = (tox_at < 0.9);
  if((!q1 && q2) || (q1 && q3)) {
    igam = (gamma == 0.0) ? 10.0 : 10.0*tox_at + 1;
	if(tox_at < 0.001) igam--;
	if(istep == 1) igam++;
	if(istep > 0) gamma_low  = gamma;
	if(istep < 0) gamma_high = gamma;
	did_break = FALSE;
	while((istep == -1) || (igam <= 10)) {
	  if((istep != 1) && (igam < 1)) {
	    gamma = gamma_int[mod_type][0]/dose_last;
	    parms();
	    return;
	  }
      gamma_at = gamma_int[mod_type][igam-1]/dose_last;
	  if(istep*f_prime(gamma_at) >= 0) {
	    igam += istep;
	    if(istep > 0) gamma_low  = gamma_at;
	    if(istep < 0) gamma_high = gamma_at;
	  } else {
        if(istep > 0) gamma_high = gamma_at;
	    if(istep < 0) gamma_low  = gamma_at;
	    did_break = TRUE;
	    break;  /* This will force an exit from the while loop */
	  }
	}
	if(!did_break) {
	  gamma_low = 0.0;
	  gamma_high = gamma_int[mod_type][9]/dose_last;
	}
  }
  g_diff = (gamma_high - gamma_low)/1597.0;
  ig0 = 0;
  ig1 = 610;
  ig2 = 987;
  ig3 = 1597;
  g2 = gamma_low + ig2*g_diff;
  gamma = g2;
  parms();
  g2lik = nblik();
  g1 = gamma_low + ig1*g_diff;
  gamma = g1;
  parms();
  g1lik = nblik();
  for(loop=1;loop<=13;loop++) {
	if(g2lik >= g1lik) {
	  ig3 = ig2;
	  ig2 = ig1;
	  g2lik = g1lik;
	  ig1 = ig0 + ig3 - ig2;
	  g1 = gamma_low + ig1*g_diff;
	  gamma = g1;
	  parms();
	  g1lik = nblik();
	} else {
      ig0 = ig1;
	  ig1 = ig2;
	  g1lik = g2lik;
	  ig2 = ig0 + ig3 - ig1;
	  g2 = gamma_low + ig2*g_diff;
	  gamma = g2;
	  parms();
	  g2lik = nblik();
	}
  }
  gamma = g2;
  if(g1lik < g2lik) gamma = g1;
  parms();
  return;
}


double f_prime(gamma)

double gamma;

{

double sumwdx, sumwd, sumwx, sumw, sumwdd;
double sumwdx_prime, sumwd_prime, sumwx_prime, sumw_prime, sumwdd_prime;
double tox[11], tj;
double tox_prime[11], tj_prime;
double muhat_prime[11];
double mut[11];
int j;
int flag_alpha;
double alpha_prime, beta_prime, f_prim;
double beta_num, beta_denom;


  sumwdx = 0.0;
  sumwd = 0.0;
  sumwx = 0.0;
  sumw = 0.0;
  sumwdd = 0.0;
  sumwdx_prime = 0.0;
  sumwd_prime = 0.0;
  sumwx_prime = 0.0;
  sumw_prime = 0.0;
  sumwdd_prime = 0.0;

  for(j=1;j<=numdos;j++) {
    if(model != 0 && model != 4) {
  	tox[j] = 2.0 - exp(gamma*dose[j]);
  	tox_prime[j] = dose[j]*(tox[j] - 2.0);
    } else {
      tox[j] = exp(-gamma*dose[j]);
  	tox_prime[j] = -dose[j]*tox[j];
    }
    tj = tox[j];
    tj_prime = tox_prime[j];
    sumwx += (weight[j]*tj*tj)         *(xbar[j]/tj);
    sumw  += (weight[j]*tj*tj);
    sumwx_prime += weight[j]*tj_prime*xbar[j];
    sumw_prime  += weight[j]*tj*tj_prime*2.0;
    if(model < 4) {
   	  sumwdx += (weight[j]*tj*tj)*dose[j]*(xbar[j]/tj);
	  sumwd  += (weight[j]*tj*tj)*dose[j];
  	  sumwdd += (weight[j]*tj*tj)*dose[j]*dose[j];
	  sumwdx_prime += weight[j]*tj_prime*dose[j]*xbar[j];
	  sumwd_prime  += weight[j]*tj*tj_prime*dose[j]*2.0;
	  sumwdd_prime += weight[j]*tj*tj_prime*dose[j]*dose[j]*2.0;
    }
  }

  beta = 0.0;
  if(model < 4) {
	beta_num = sumw*sumwdx - sumwd*sumwx;
	beta_denom = sumw*sumwdd - sumwd*sumwd;
	beta = beta_num/beta_denom;
	if(beta < 0.0 || model >= 3) beta = 0.0;
  }
  alpha = (sumwx - beta*sumwd)/sumw;

  flag_alpha = FALSE;
  if(alpha <= 0.0) {
	flag_alpha = TRUE;
	alpha = 1.0;
	beta = (sumwdx - sumwd)/sumwdd;
  }
  beta_prime = 0.0;
  if(beta != 0.0 && !flag_alpha) {
    beta_prime = ((beta_denom*(sumw*sumwdx_prime +
              sumw_prime*sumwdx - sumwd*sumwx_prime - sumwd_prime*sumwx) -
              beta_num*(sumw*sumwdd_prime + sumw_prime*sumwdd -
              2.0*sumwd*sumwd_prime))/beta_denom)/beta_denom;
  }
  if(flag_alpha) {
	beta_prime = (sumwdd*(sumwdx_prime - sumwd_prime) -
	(sumwdx - sumwd)*sumwdd_prime)/(sumwdd*sumwdd);
  }
  alpha_prime = 0.0;
  if(flag_alpha) {
	alpha_prime = (sumwx_prime - (beta*sumwd_prime + sumwd*beta_prime))/sumw
	- (alpha*sumw_prime)/sumw;
  }
  f_prim = 0.0;
  for(j=1;j<=numdos;j++) {
	mut[j] = alpha + beta*dose[j];
	muhat[j] = mut[j]*tox[j];
	muhat_prime[j] = mut[j]*tox_prime[j] +
	                 tox[j]*(alpha_prime + beta_prime*dose[j]);
	f_prim -= weight[j]*muhat_prime[j]*(xbar[j] - muhat[j]);
  }

  return(f_prim);

}


void weights()

{

  int j;

  for(j=1;j<=numdos;j++) weight[j] = (muhat[j] > 0.0) ?
                                     n[j]/(muhat[j]*(1.0 + c*muhat[j])) :
                                     1.0;
}


double addconx()

{

  double sum;
  int i, j;

  sum = 0.0;
  for(j=1;j<=numdos;j++) {
	for(i=1;i<=n[j];i++) sum -= dlgama((double) kount[j][i] + 1.0);
  }
  return(sum);
}


double nblik()

{

  int j;
  double cj, sum, ratioj;

/* check for violated constraints.  if found, impose penalty. */

  sum = 0.0;
  for(j=1;j<=numdos;j++) {
	if(muhat[j] > 0.0) {
	  cj = 1.0 + c*muhat[j];
	  ratioj = (c != 0.0) ? log(cj)/c : muhat[j];
	  sum -= n[j]*(xbar[j]*log(muhat[j]/cj) - ratioj);
	} else {
	  if(xbar[j] != 0.0) return(-7777.0);
	}
  }
  return(sum);
}


double addcon(xsum,nsum)

double xsum;
int nsum;

{

  double sum, cinv;
  int i, j;

  sum = adtvcnx;
  if(c == 0.0) return(sum);
  cinv = 1.0/c;
  for(j=1;j<=numdos;j++) {
	for(i=1;i<=n[j];i++) sum += dlgama((double) kount[j][i] + cinv);
  }
  sum -= (nsum*dlgama(cinv) - xsum*log(c));
  return(sum);
}


double cmputc()

{

  double c;
  double ci;
  double z;
  double z2;
  double z2sum;
  int i, j;

  c = 0.0;
  z2sum = 0.0;
  for(i=1;i<=numdos;i++) {
	if(xbar[i] != 0.0 && n[i] > 1) {
	  ci = 0.0;
	  z = xbar[i];
	  z2 = z*z;
	  z2sum += z2;
	  for(j=1;j<=n[i];j++) ci += (z2 != 0.0) ?
	                             kount[i][j]*((double) kount[i][j] - 1)/z2 : 0;
	  ci = (ci + 1.0/z - n[i])/(n[i]-1);
	  c += ci*z2;
	}
  }
  c = c/z2sum;
  return((c < 1.0e-8) ? 0.0 : c);
}


double dlgama(z)

double z;

{

  double zfac;
  double s1;
  double s2;
  double x;
  double dlgam;
  int i;
  static double a[7] = {0.833333333333333e-1,  0.333333333333333e-1,
	                    0.252380952380952e0,  0.525606469002695e0,
	                    0.101152306812684e1,  0.151747364915329e1,
	                    0.226948897420496e1 };
  static double dcon = 0.918938533204673e0;
  static double plgam3[8] = {-2.407986980173375e8, -2.448321769032882e8,
	                          4.818077102773628e8,  1.119385354299855e8,
	                         -8.731675438238387e7, -1.971830115860921e7,
	                         -9.827102281420492e5, -1.037701651732974e4 };
  static double plgam2[8] = {-1.513831834115067e3, -1.508630228766725e4,
	                         -2.064829420532528e4,  1.204317380987164e4,
	                          1.955360554063045e4,  5.268983255914981e3,
	                          3.775106797972170e2,  5.155057617640817e0 };
  static double plgam1[8] = {-2.208843997216182e0, -7.741064071332953e1,
	                         -5.176383498023218e2, -9.222613728801522e2,
	                         -2.617218583856145e2,  2.431752435244210e2,
	                          8.568982062831317e1,  4.120843185847770e0 };
  static double qlgam3[8] = {-7.902611114187634e7, -4.353707148043742e8,
	                         -4.044359282914355e8, -1.119254116263318e8,
	                         -1.048577583049937e7, -3.114062847340679e5,
	                         -2.015275195500483e3,  1.000000000000000e0 };
  static double qlgam2[8] = { 6.983274140573510e2,  1.440209037170085e4,
	                          5.262286383841199e4,  5.712025539602503e4,
	                          2.202956214415664e4,  3.039903041439440e3,
	                          1.289093189012958e2,  1.000000000000000e0 };
  static double qlgam1[8] = { 4.097792921092615e-1, 2.443519662506312e1,
	                          2.623083470269460e2,  8.460755362020782e2,
	                          9.513235976797060e2,  3.778372484823942e2,
	                          4.564677187585908e1,  1.000000000000000e0 };

/* the interval [12,infinity] . . . */

  if(z >= 12.0) {
    dlgam = (z - 0.5)*log(z) + dcon - z;
	x = a[6]/z;
	for(i=5;i>=0;i--) x = a[i]/(z + x);
	return(dlgam + x);

  } else if(z >= 4.0) {

/* the interval [4.0,12.0] . . . */

    s1 = plgam3[7]*z + plgam3[6];
	s2 = qlgam3[7]*z + qlgam3[6];
	for(i=5;i>=0;i--) {
	  s1 = s1*z + plgam3[i];
 	  s2 = s2*z + qlgam3[i];
	}
	return(s1/s2);

  } else if(z >= 1.5) {

/* the interval [1.5,4.0] . . . */

	s1 = plgam2[7]*z + plgam2[6];
	s2 = qlgam2[7]*z + qlgam2[6];
	for(i=5;i>=0;i--) {
	  s1 = s1*z + plgam2[i];
	  s2 = s2*z + qlgam2[i];
	}
	return((z - 2.0)*(s1/s2));
	
  } else {
	
/* the interval [0.0,1.5]: for z in [0.0,0.5] we evaluate instead the
   function at z+1 and the use the recursion mentioned above . . .  */

	zfac = z + (z < 0.50) ? 1.0 : 0.0;
	s1 = plgam1[7]*zfac + plgam1[6];
	s2 = qlgam1[7]*zfac + qlgam1[6];
	for(i=5;i>=0;i--) {
	  s1 = s1*zfac + plgam1[i];
	  s2 = s2*zfac + qlgam1[i];
	}
	return((z >= 0.5) ? (z - 1.0)*(s1/s2) : -log(z) + z*(s1/s2));
  }
	
}

double chicdf(chi,idegre)

double chi;
int idegre;
{

  #define PIE 3.1415926536

  int i, index;
  double product, sum, rootch, atemp, zxsum, dlnzx, chisum;
  double return_value;

  if(idegre > 100) {
    return_value = ztail(sqrt(2.0*chi) - sqrt((double) (2.0*idegre-1)));
    return(return_value);
  }
  product = 1.0;
  sum = 0.0;
  rootch = sqrt(chi);
  atemp = sqrt(2.0*PIE);
  if((idegre%2) != 0) {
    zxsum = 0.0;
    if(idegre == 1) {
      return_value = 2.0*(ztail(rootch));
      return(return_value);
    }
    dlnzx = -0.5*chi - log(atemp);
    for(i=2;i<=idegre-1;i+=2) {
      product *= (double) (i-1);
      sum += pow(rootch,(double) (i-1))/product;
    }
    sum = dlnzx + log(sum);
    zxsum = (sum > -85.0) ? exp(sum) : 0.0;
    return_value = 2.0*(ztail(rootch) + zxsum);
    return(return_value);
  }
  index = (idegre-2)/2;
  dlnzx = -0.5*chi;
  chisum = 1.0;
  for(i=1;i<=index;i++) {
    product *= (double) i;
    chisum += pow((0.5*chi),(double) i)/product;
  }
  sum = dlnzx + log(chisum);
  zxsum = (sum > -85.0) ? exp(sum) : 0.0;
  return_value = zxsum;
  return(zxsum);
}

double ztail(z)

double z;

{

  double a[6],t,x,sum;
  int i,isign;

  a[0] = 0.0;
  a[1] = 0.254829592;
  a[2] = -0.284496736;
  a[3] = 1.421413741;
  a[4] = -1.453152027;
  a[5] = 1.061405429;

  isign = (z < 0.0) ? -1 : 1;
  x = fabs(z)/sqrt((double) 2.0);
  if(x > 6.0) {
    return((isign == 1) ? 0.0 : 1.0);
  }
  t = 1.0/(1.0 + 0.3275911*x);
  sum = 0.0;
  for(i=1;i<=5;i++) {
    sum += a[i]*pow(t,(double) i);
  }
  return(0.5 - isign*0.5*(1.0 - sum*exp(-x*x)));
}


void fisher(covar)

double covar[3][3];

{

int i, j, l;
double tox, prime[3];

  covar[0][0] = 0.0;
  covar[1][0] = 0.0;	
  covar[1][1] = 0.0;
  covar[2][0] = 0.0;
  covar[2][1] = 0.0;
  covar[2][2] = 0.0;

  for(j=1;j<=numdos;j++) {
    xmut = alpha + beta*dose[j];
    tox = 1.0;
    if(model == 0 || model == 4) {
      tox = exp(-gamma*dose[j]);
    } else if(model == 1 || model == 5) {
      tox = 2.0 - exp(gamma*dose[j]);
    }

    prime[0] = tox;
    prime[1] = dose[j]*tox;
    prime[2] = 1.0;
    if(model == 0 || model == 4) {
      prime[2] = -xmut*dose[j]*tox;
    } else if(model == 1 || model == 5) {
      prime[2] = -xmut*dose[j]*(2.0 - tox);
    }

    for(i=0;i<=2;i++) {
      for(l=0;l<=i;l++) {
        covar[i][l] += weight[j]*prime[i]*prime[l];
      }
    }
  }

  if(model <= 1) return;
  if(model <= 3) {
    covar[2][0] = 0.0;
    covar[2][1] = 0.0;
    covar[2][2] = 1.0;
    if(model == 2) return;
  }
  covar[1][0] = 0.0;
  covar[1][1] = 1.0;
  covar[2][1] = 0.0;
  return;
}



void muhatc()

{

int j;
double tox;

  for(j=1;j<=numdos;j++) {
    tox = 1.0;
    if(model == 0 || model == 4) {
      tox = exp(-gamma*dose[j]);
    } else if(model == 1 || model == 5) {
      tox = 2.0 - exp(gamma*dose[j]);
    }
    muhat[j] = (alpha + beta*dose[j])*tox;
  }

  return;

}


int inverse(matrix)

double matrix[3][3];

{

  double x[3][3], y[3][3], temp;

  if(matrix[0][0] <= 0.0) return(1);

  x[0][0] = sqrt(matrix[0][0]);
  x[1][0] = matrix[1][0]/x[0][0];
  temp = matrix[1][1] - x[1][0]*x[1][0];
  if(temp <= 0.0) return(2);

  x[1][1] = sqrt(temp);
  x[2][0] = matrix[2][0]/x[0][0];
  x[2][1] = (matrix[2][1] - x[2][0]*x[1][0])/x[1][1];
  temp = matrix[2][2] - (x[2][0]*x[2][0] + x[2][1]*x[2][1]);
  if(temp <= 0.0) return(3);

  x[2][2] = sqrt(temp);

  y[0][0] = 1.0/x[0][0];
  y[1][0] = -x[1][0]/(x[0][0]*x[1][1]);
  y[1][1] = 1.0/x[1][1];
  y[2][0] = (x[1][0]*x[2][1] - x[2][0]*x[1][1])/(x[0][0]*x[1][1]*x[2][2]);
  y[2][1] = -x[2][1]/(x[1][1]*x[2][2]);
  y[2][2] = 1.0/x[2][2];

  matrix[0][0] = y[0][0]*y[0][0] + y[1][0]*y[1][0] + y[2][0]*y[2][0];
  matrix[1][0] = y[1][0]*y[1][1] + y[2][0]*y[2][1];
  matrix[1][1] = y[1][1]*y[1][1] + y[2][1]*y[2][1];
  matrix[2][0] = y[2][0]*y[2][2];
  matrix[2][1] = y[2][1]*y[2][2];
  matrix[2][2] = y[2][2]*y[2][2];
  return(0);
}
#include <stdio.h>
#include <math.h>
#define TRUE 1
#define FALSE 0

void weights();
void gamma_find();
double nblik();
int inverse(double matrix[3][3]);
double f_prime(double gamma);
double dlgama(double z);
double ztail(double z);




int numdos, mumdos;
double xbar[11];
float dose[11];
int kount[11][11];
char backgr[11][11];
int n[11], nn[11];
double weight[11];
double muhat[11];
double adtvcnx;
double adtvcn;
double alpha;
double beta;
double gamma;
double c;
int model;
float outparm[30];
double stats[11];
double dose_last;
double xmut;

FILE *fopen(), *ip, *op;


main ()

{

char line[80], ch;
char raw_backgr[11][11], reject_dose[11], reject_kount[11][11];
float raw_dose[11];
int raw_kount[11][11];
int raw_n[11];
long kntsum;
int i, j, k, l, imany, jstat, numtox;
double xsum, dsum, ddsum, dxsum, xxbar;
double c_all, w, cc[11], x, x2, cnew, z[11];
int nsum, jmax;
int model2_flag, finish_up, nodel;
double fnew;
int ioffset, ierror, dataleft;
#define GCUT 9.5
int peak_dose, ndf, isign;
double peaktest, tvalue, sigma1, test;
double covar[3][3];
static double peak[8] = {0.0, 2.06, 2.16, 2.23, 2.29, 2.34, 2.38, 2.42};
static double sigmul[23] = {0.000, 2.539, 2.576, 2.609, 2.638, 2.665, 2.690,
                            2.713, 2.734, 2.754, 2.773, 2.791, 2.807, 2.823,
                            2.838, 2.852, 2.865, 2.878, 2.891, 2.902, 2.914,
                            2.925, 2.935 };

void get_initial();
double nblik();
void iterate();
void weights();
void parms();
void gamma_find();
double f_prime();
double addconx();
double addcon();
double cmputc();
double dlgama();
void output();
double chicdf();
double ztail();
void fisher();
void muhatc();
int inverse();

ip = fopen("salm.in","r");
if (!ip) {
  printf("Error on salm.in");
  exit (1);
}

op = fopen("salm.out","w");
if (!op) {
  printf("Error on salm.out");
  exit (1);
}
/*
fprintf(stdout,"\7\7\7"
               "                          #############################\n"
               "                          ####                     ####\n"
               "                          #### Salmonella Analysis ####\n"
               "                          #### á release  ver. 1.1 ####\n"
               "                          ####  November 11, 1990  ####\n"
               "                          ####                     ####\n"
               "                          #### Report any problems ####\n"
               "                          #### to Kenneth J. Risko ####\n"
               "                          ####                     ####\n"
               "                          #############################\n"
               "\7\7\7");

fprintf(op,"                          #############################\n"
           "                          ####                     ####\n"
           "                          #### Salmonella Analysis ####\n"
           "                          #### á release  ver. 1.1 ####\n"
           "                          ####  November 11, 1990  ####\n"
           "                          ####                     ####\n"
           "                          #### Report any problems ####\n"
           "                          #### to Kenneth J. Risko ####\n"
           "                          ####                     ####\n"
           "                          #############################\n");


fprintf(stdout,"Press any key to continue");
fscanf(stdin,"%c",&ch);
*/

while(fgets(line,80,ip)) {

  fprintf(stdout,"\n\n");
  fprintf(op,"\n\n");

  i = 0;
  while((ch = line[i++]) != '\n') {
    fprintf(stdout,"%c",ch);
    fprintf(op,"%c",ch);
  }
  fprintf(stdout,"%c",ch);
  fprintf(op,"%c",ch);

  fscanf(ip,"%d",&numdos);
  fprintf(op,"number of doses = %2d\n",numdos);
  fprintf(stdout,"number of doses = %2d\n",numdos);
  for(i=1;i<=numdos;i++) {
    fscanf(ip,"%f,%d",&raw_dose[i],&raw_n[i]);
    fprintf(op,"%2d,%7.2f,%2d",i,raw_dose[i],raw_n[i]);
    fprintf(stdout,"%2d,%7.2f,%2d",i,raw_dose[i],raw_n[i]);
    reject_dose[i] = ' ';
    for(j=1;j<=raw_n[i];j++) {
      fscanf(ip,",%d,%c",&raw_kount[i][j],&raw_backgr[i][j]);
      fprintf(op,",%4d(%c)",raw_kount[i][j],raw_backgr[i][j]);
      fprintf(stdout,",%4d(%c)",raw_kount[i][j],raw_backgr[i][j]);
      reject_kount[i][j] = ' ';
    }
    fprintf(op,"\n");
    fprintf(stdout,"\n");
  }

/***** CLEANSE DATA *****/

  for(i=1;i<=numdos;i++) {
    for(j=1;j<=raw_n[i];j++) {
      if(raw_backgr[i][j] != ' ' &&
         raw_backgr[i][j] != 'P' &&
         raw_backgr[i][j] != 'T' &&
         raw_backgr[i][j] != 'S'   ) {
           reject_kount[i][j] = raw_backgr[i][j];
           fprintf(op,"Rejecting plate count %4d at dose %7.2f because "
             "background code = %c\n",raw_kount[i][j],raw_dose[i],
             raw_backgr[i][j]);
           fprintf(stdout,"Rejecting plate count %4d at dose %7.2f because "
             "background code = %c\n",raw_kount[i][j],raw_dose[i],
             raw_backgr[i][j]);
      }
    }
  }

  fprintf(op,"\n");
  fprintf(stdout,"\n");

  for(i=1;i<=numdos;i++) {
    numtox = 0;
    for(j=1;j<=raw_n[i];j++) {
      if(raw_backgr[i][j] == 'T') {
        ++numtox;
        reject_kount[i][j] = 'T';
        fprintf(op,"Rejecting plate count %4d at dose %7.2f because "
          "background code = %c\n",raw_kount[i][j],raw_dose[i],raw_backgr[i][j]);
        fprintf(stdout,"Rejecting plate count %4d at dose %7.2f because "
          "background code = %c\n",raw_kount[i][j],raw_dose[i],raw_backgr[i][j]);
      }
    }
    if(3*numtox > raw_n[i]) {
      for(j=i;j<=numdos;j++) {
        reject_dose[j] = 'T';
        fprintf(op,"Rejecting dose %7.2f because of toxicity\n",raw_dose[j]);
        fprintf(stdout,"Rejecting dose %7.2f because of toxicity\n",raw_dose[j]);
      }
      break;
    }
  }

  numdos = i - 1;

  for(i=1;i<=numdos;i++) {
    kntsum = 0;
    dataleft = FALSE;
    for(j=1;j<=raw_n[i];j++) {
      if(reject_kount[i][j] == ' ') {
        dataleft = TRUE;
        kntsum += raw_kount[i][j];
      }
    }
    if(kntsum == 0 && dataleft) {
      for(j=i;j<=numdos;j++) {
        reject_dose[j] = '0';
        fprintf(op,"Rejecting dose %7.2f because of zero mean "
          "at or before dose %6f\n",raw_dose[j],raw_dose[j]);
        fprintf(stdout,"Rejecting dose %7.2f because of zero mean "
          "at or before dose %6f\n",raw_dose[j],raw_dose[j]);
      }
      break;
    }
  }

  numdos = i - 1;

  if(numdos == 0) return;

  for(i=numdos;i>=1;i--) {
    dose[i] = raw_dose[i];
    n[i] = 0;
    for(j=1;j<=raw_n[i];j++) {
      if(reject_kount[i][j] == ' ') {
        n[i] += 1;
        kount[i][n[i]] = raw_kount[i][j];
        backgr[i][n[i]] = raw_backgr[i][j];
      }
    }
    if(n[i] == 0) {
      if(i != numdos) {
        for(k=i;k<=numdos-1;k++) {
          dose[k] = dose[k+1];
          n[k] = n[k+1];
          for(l=1;l<=n[k];l++) {
            kount[k][l] = kount[k+1][l];
            backgr[k][l] = backgr[k+1][l];
          }

        }
      }
    --numdos;
    }
  }

/***** REMOVE OUTLIERS *****/

  c_all = 0.0;
  w = 0.0;
  nsum = 0;
  for(i=1;i<=numdos;i++) {
    nsum += n[i];
    cc[i] = 0.0;
    xbar[i] = 0.0;
    for(j=1;j<=n[i];j++) xbar[i] += kount[i][j];
    xbar[i] /= n[i];
    if(xbar[i] != 0.0 && n[i] > 1) {
      x = xbar[i];
      x2 = x*x;
      w = w + x2;
      for(j=1;j<=n[i];j++) cc[i] += kount[i][j]*((double) kount[i][j] - 1.0)/x2;
      cc[i] = (cc[i] + 1.0/x - n[i])/(n[i] - 1);
      c_all += cc[i]*x2;
    }
  }
  c_all /= w;

  for(i=1;i<=numdos;i++) {
    if(n[i] >= 3) {
      cnew = (w*c_all - xbar[i]*xbar[i]*cc[i])/(w - xbar[i]*xbar[i]);
      if(cnew < 0.0) cnew = 0.0;
      jmax = 1;
      for(j=1;j<=n[i];j++) {
        z[j] = (xbar[i] == 0.0) ? 0.0 : fabs(kount[i][j] - xbar[i])/
                                  sqrt(2.0*xbar[i]*(1.0 + cnew*xbar[i])/3.0);
        if(z[j] > z[jmax]) jmax = j;
      }
      if(z[jmax] >= sigmul[nsum-8]) {
        fprintf(op,"Plate count %4d at dose %7.2f determined to "
          "an outlier\n",kount[i][jmax],dose[i]);
        fprintf(stdout,"Plate count %4d at dose %7.2f determined to "
          "an outlier\n",kount[i][jmax],dose[i]);
        --n[i];
        for(k=jmax;k<=n[i];k++) {
          kount[i][k] = kount[i][k+1];
          backgr[i][k] = backgr[i][k+1];
        }
      }
    }
  }

  for(i=1;i<=numdos;i++) {
	xbar[i] = 0.0;
	for(j=1;j<=n[i];j++) xbar[i] += kount[i][j];
	xbar[i] = xbar[i]/n[i];
    if(xbar[i] == 0.0) {
      for(j=i;j<=numdos;j++) {
        fprintf(op,"Rejecting dose %7.2f because of zero mean (after outlier"
          " deletion) at or before dose %6f\n",raw_dose[j],raw_dose[j]);
        fprintf(stdout,"Rejecting dose %7.2f because of zero mean (after outlier"
          " deletion) at or before dose %6f\n",raw_dose[j],raw_dose[j]);
      }
      break;
    }
  }

  numdos = i - 1;

  if(numdos < 4) {
    fprintf(stdout,"Too few remaining doses (%2d) to perform analysis\n",numdos);
    return(-1);
  }

  imany = 0;
  for(j=1;j<=numdos;j++) imany += (n[j] >= 2) ? 1 : 0;
  if(imany < 3) {
	fprintf(stdout,"Too few doses with replicates (%2d) to perform analysis\n",imany);
	return(-2);
  }

  xsum  = 0.0;
  nsum  = 0;
  dsum  = 0.0;
  ddsum = 0.0;
  dxsum = 0.0;
  for(j=1;j<=numdos;j++) {
	nsum  += n[j];
	dsum  += n[j]*dose[j];
	ddsum += n[j]*dose[j]*dose[j];
	dxsum += n[j]*dose[j]*xbar[j];
	xsum  += n[j]*xbar[j];
  }
  xxbar = xsum/nsum;

  c = cmputc();

  peak_dose = 1;
  for(i=2;i<=numdos;i++) peak_dose = (xbar[i] >= xbar[peak_dose]) ?
                                     i : peak_dose;
  stats[1] = -peak_dose;
  for(i=peak_dose+1;i<=numdos;i++) {
    peaktest = (xbar[peak_dose] - xbar[i])/
               sqrt(xbar[peak_dose]*(1.0 + c*xbar[peak_dose])/n[peak_dose] +
                    xbar[i]*(1.0 + c*xbar[i])/n[i]);
    if(peaktest >= 2.576) {
      for(j=i+1;j<=numdos;j++) {
        fprintf(op,"Rejecting dose %2d because depressed mean response of "
          "%6.2f is substantially below peak response of %6.2f\n",j,xbar[j],
          xbar[peak_dose]);
        fprintf(stdout,"Rejecting dose %2d because depressed mean response of "
          "%6.2f is substantially below peak response of %6.2f\n",j,xbar[j],
          xbar[peak_dose]);
      }
      break;
    }
  }
  numdos = (i < numdos) ? i : numdos;

  if(numdos <= 3) {
    fprintf(stdout,"Too few remaining doses (%2d) to perform analysis\n",numdos);
    return(-3);
  }

  fprintf(op,"\n *** ANALYSED DATA ***\n\n");
  fprintf(op," No.   Dose                 Plates\n"
                 " ---   ----             --------------\n");
  fprintf(stdout,"\n *** ANALYSED DATA ***\n\n");
  fprintf(stdout," No.   Dose                 Plates\n"
                 " ---   ----             --------------\n");
  for(i=1;i<=numdos;i++) {
    fprintf(op," %2d  %7.2f  ",i,dose[i]);
    fprintf(stdout," %2d  %7.2f  ",i,dose[i]);
    for(j=1;j<=n[i];j++) fprintf(op,"%5d (%c) ",kount[i][j],backgr[i][j]);
    for(j=1;j<=n[i];j++) fprintf(stdout,"%5d (%c) ",kount[i][j],backgr[i][j]);
    fprintf(op,"\n");
    fprintf(stdout,"\n");
  }

  fprintf(stdout,"Running");

  for(i=1;i<=26;i++) outparm[i] = -1.0;

  c = cmputc();

  dose_last = dose[numdos];
  adtvcnx = addconx();
  adtvcn = addcon(xsum,nsum);
  for(i=1;i<=numdos;i++) muhat[i] = xxbar;
  outparm[13] = xxbar;
  outparm[14] = 0.0;
  outparm[15] = 0.0;
  outparm[16] = nblik() - adtvcn;	

  for(i=1;i<=numdos;i++) muhat[i] = xbar[i];
  outparm[25] = c;
  outparm[26] = nblik() - adtvcn;

  model2_flag = FALSE;
  model = 0;
  weights();
  gamma = 0.0;
  numdos -= 1;
  parms();
  numdos += 1;
  xmut = alpha + beta*dose_last;
  gamma = -log(xbar[numdos]/xmut)/dose_last;
  if(gamma < 0.0) gamma = 0.0;
  parms();
  weights();
  finish_up = FALSE;

  while(!finish_up) {
    iterate(&fnew);
    ioffset = 4*model;
    outparm[ioffset+1] = alpha;
    outparm[ioffset+2] = beta;
    outparm[ioffset+3] = gamma;
    outparm[ioffset+4] = fnew;
    if(model == 5) {
      finish_up = TRUE;
    } else {
      if(model == 4) {
        model = 1;
        get_initial(model);
      } else {
        if(gamma < 0.0) {
          model2_flag = TRUE;
          for(model=0;model<=2;model++) {
	        ioffset = 4*model;
	        outparm[ioffset+1] = alpha;
	        outparm[ioffset+2] = beta;
	        outparm[ioffset+3] = 0.0;
	        outparm[ioffset+4] = fnew;
          }
          for(model=4;model<=5;model++) {
	        ioffset = 4*model;
	        for(jstat=1;jstat<=4;jstat++) outparm[ioffset+jstat] = outparm[12+jstat];
          }
          output();
          return;
        }
        model = model + 4;
        if(beta <= 0.0 || outparm[13] <= (dxsum/dsum)) {
          ioffset = 4*model;
          outparm[ioffset+1] = (beta <= 0.0) ? alpha : outparm[13];
          outparm[ioffset+2] = 0.0;
          outparm[ioffset+3] = (beta <= 0.0) ? gamma : 0.0;
          outparm[ioffset+4] = (beta <= 0.0) ? fnew : outparm[16];
          if(model == 5) {
            finish_up = TRUE;
          } else {
            model = 1;
            get_initial(model);
          }
        } else {
          if(model == 5) {
            get_initial(model);
          } else {
            double ysum = 0.0;
            double dysum = 0.0;
            double y;
            double determ;
            for(j=1;j<=numdos;j++) {
	          for(i=1;i<=n[j];i++) {
	            y = kount[j][i];
	            if(y < 1.0) y = 1.0;
	            y = log(y);
	            ysum += y;
	            dysum += dose[j]*y;
              }
            }
            determ = nsum*ddsum - dsum*dsum;
            gamma = (dsum*ysum - nsum*dysum)/determ;
            if(gamma < 0.0) gamma = 0.0;
            alpha = exp((ysum + gamma*dsum)/nsum);
            beta = 0.0;
            for(j=1;j<=numdos;j++) muhat[j] = alpha*exp(-gamma*dose[j]);
          }
          weights();
        }
      }
    }
  }

  if(!model2_flag) {
    if(outparm[13] >= (dxsum/dsum)) {
      outparm[9] = outparm[13];
      outparm[10] = 0.0;
      outparm[11] = 0.0;
      outparm[12] = outparm[16];
    } else {
      double determ;
      determ = nsum*ddsum - dsum*dsum;
      beta = (nsum*dxsum - dsum*xsum)/determ;
      alpha = (xsum - beta*dsum)/nsum;
      gamma = 0.0;
      for(j=1;j<=numdos;j++) muhat[j] = alpha + beta*dose[j];
      weights();
      model = 2;
      iterate(&fnew);
      outparm[9] = alpha;
      outparm[10] = beta;
      outparm[11] = 0.0;
      outparm[12] = fnew;
    }
  }

  fprintf(stdout,"\n\n");
  c = outparm[25];
  model = (outparm[4] < outparm[8]) ? 0 : 1;
  ioffset = 4*model;
  if(outparm[12] - outparm[ioffset+4] < 0.5*GCUT) {
    model = 2;
    tvalue = 2.0*(outparm[16] - outparm[12]);
    tvalue = (tvalue > 0.0) ? sqrt(tvalue) : 0.0;
  } else {
    tvalue = 2.0*(outparm[ioffset+20] - outparm[ioffset+4]);
    tvalue = (tvalue > 0.0) ? sqrt(tvalue) : 0.0;
    sigma1 = xbar[1]*(1.0 + c*xbar[1])/n[1];
    isign = -1;
    for(i=2;i<=numdos;i++) {
      test = sqrt((1.0 + (double) n[1]/n[i])*sigma1)*peak[numdos-3] + xbar[1];
      if(xbar[i] > test) isign = 1;
    }
    tvalue = isign*tvalue;
  }
  stats[1] = model;
  ioffset = 4*model;
  stats[2] = outparm[ioffset+1];
  stats[3] = outparm[ioffset+2];
  stats[4] = outparm[ioffset+3];
  stats[5] = c;
  stats[6] = outparm[ioffset+4];
  stats[7] = 2.0*(outparm[ioffset+4] - outparm[26]);
  ndf = numdos - ((model ==2) ? 2 : 3);
  stats[8] = chicdf(stats[7],ndf);
  stats[9] = tvalue;
  stats[10] = ztail(fabs(stats[9]));

  for(nodel=0;nodel<=2;nodel++) {
    model = nodel;
    if(outparm[4*model+3] == 0.0) model = 2;
    ioffset = 4*model;
    alpha = outparm[ioffset+1];
    beta  = outparm[ioffset+2];
    gamma = outparm[ioffset+3];
    muhatc();
    weights();
    fisher(covar);
    ierror = inverse(covar);
    tvalue = 0.0;
    if(!ierror) tvalue = beta/sqrt(covar[1][1]);
    outparm[27+nodel] = tvalue;
  }

  output();

  fgets(line,80,ip);

}

return;

}


void output()

{

int i, j, model, maxn;
double fit, xbar_max;

  model = stats[1];
  switch (model) {
    case 0:
      fprintf(stdout,"\n\nFitted model = (%6f + %6f x dose) x "
                     "exp( -%6f x dose)\n",stats[2],stats[3],stats[4]);
      fprintf(op,"\n\nFitted model = (%6f + %6f x dose) x "
                     "exp( -%6f x dose)\n",stats[2],stats[3],stats[4]);
      break;
    case 1:
      fprintf(stdout,"\n\nFitted model = (%6f + %6f x dose) x "
                     "(2 - exp( %6f x dose))\n",stats[2],stats[3],stats[4]);
      fprintf(op,"\n\nFitted model = (%6f + %6f x dose) x "
                     "(2 - exp( %6f x dose))\n",stats[2],stats[3],stats[4]);
      break;
    case 2:
      fprintf(stdout,"\n\nFitted model = (%6f + %6f x dose)\n",
                      stats[2],stats[3]);
      fprintf(op,"\n\nFitted model = (%6f + %6f x dose)\n",
                      stats[2],stats[3]);
      break;
    default:
      break;
  }

  fprintf(stdout,"\n *** ANALYSED DATA ***\n\n");
  fprintf(stdout," No.   Dose          Plates"
                 "         Mean        Model\n"
                 " ---   ----      --------------"
                 "   --------    ---------\n");
  fprintf(op,"\n *** ANALYSED DATA ***\n\n");
  fprintf(op," No.   Dose          Plates"
                 "         Mean        Model\n"
                 " ---   ----      --------------"
                 "   --------    ---------\n");

  xbar_max = 0.0;
  maxn = 0;
  for(i=1;i<=numdos;i++) {
    xbar_max = (xbar[i] > xbar_max) ? xbar[i] : xbar_max;
    maxn = (n[i] > maxn) ? n[i] : maxn;
    switch (model) {
      case 0:
        fit = (stats[2] + stats[3]*dose[i])*exp(-stats[4]*dose[i]);
        break;
      case 1:
        fit = (stats[2] + stats[3]*dose[i])*(2-exp(stats[4]*dose[i]));
        break;
      case 2:
        fit = (stats[2] + stats[3]*dose[i]);
        break;
      default:
        fit = 0.0;
    }
    fprintf(op," %2d  %7.2f  ",i,dose[i]);
    fprintf(stdout," %2d  %7.2f  ",i,dose[i]);
    for(j=1;j<=n[i];j++) fprintf(op,"%5d ",kount[i][j]);
    for(j=1;j<=n[i];j++) fprintf(stdout,"%5d ",kount[i][j]);
    for(j=n[i]+1;j<=maxn;j++) fprintf(op,"      ");
    for(j=n[i]+1;j<=maxn;j++) fprintf(stdout,"      ");
    fprintf(op,"  %6.1f       %6.1f\n",xbar[i],fit);
    fprintf(stdout,"  %6.1f       %6.1f\n",xbar[i],fit);
  }

  fprintf(stdout,"\nc = %6.2f\n",stats[5]);
  fprintf(op,"\nc = %6.2f\n",stats[5]);

  fprintf(stdout,"\nTest for GOODNESS OF FIT p-value: %5.3f\n",stats[8]);
  fprintf(op,"\nTest for GOODNESS OF FIT p-value: %5.3f\n",stats[8]);

  fprintf(stdout,"\nTest for MUTAGENICITY p-value: %5.3f\n",stats[10]);
  fprintf(op,"\nTest for MUTAGENICITY p-value: %5.3f\n",stats[10]);

  if(stats[9] < 0.0) {
    fprintf(stdout,"\n***  NOTE: Peak mean response of %6.2f not "
                   "significanltly elevated above control mean response "
                   " of %6.2f  ***\n",xbar_max,xbar[1]);
    fprintf(op,"\n***  NOTE: Peak mean response of %6.2f not "
               "significanltly elevated above control mean response "
               " of %6.2f\n  ***",xbar_max,xbar[1]);
  }

}


void get_initial(model)

int model;

{

int j, ioffset;

  ioffset = 4*model;
  alpha = outparm[ioffset - 3];
  beta = outparm[ioffset - 2];
  gamma = log(2.0 - exp(-outparm[ioffset - 1]*dose_last))/dose_last;
  for(j=1;j<=numdos;j++) muhat[j] = (alpha + beta*dose[j])*
                                        (2.0 - exp(gamma*dose[j]));
  weights();
  return;

}


void iterate(fnew)

double *fnew;

{

  int loop;
  double fold;

  gamma_find();
  *fnew = nblik() - adtvcn;

/*  iterate and reweight  */

  for(loop=1;loop<=10;loop++) {
    fprintf(stdout,".");
	weights();
	fold = *fnew;
 	gamma_find();
 	*fnew = nblik() - adtvcn;
 	if(fabs(fold-*fnew) <= 1.0e-4) break;
  }
  return;

}


void parms()

{

  double sumwdx, sumwd, sumwx, sumw, sumwdd;
  double tox[11], tj, wt2;
  int j;

  sumwdx = 0.0;
  sumwd = 0.0;
  sumwx = 0.0;
  sumw = 0.0;
  sumwdd = 0.0;
  for(j=1;j<=numdos;j++) {
    tox[j] = (model != 0 && model != 4) ? 2.0 - exp(gamma*dose[j])
                                        : exp(-gamma*dose[j]);
    tj = tox[j];
    wt2 = weight[j]*tj*tj;
    sumwx  += wt2         *(xbar[j]/tj);
    sumw   += wt2;
    if(model < 4) {
 	  sumwdx += wt2*dose[j]*(xbar[j]/tj);
	  sumwd  += wt2*dose[j];
	  sumwdd += wt2*dose[j]*dose[j];
	}
  }
  beta = 0.0;
  if(model < 4) {
	beta = (sumwdx - sumwd*sumwx/sumw)/(sumwdd - sumwd*sumwd/sumw);
	if(beta < 0.0 || model >= 3) beta = 0.0;
  }
  alpha = (sumwx - beta*sumwd)/sumw;
  if(alpha <= 0.0) {
	alpha = 1.0;
	beta = (sumwdx - sumwd)/sumwdd;
  }
  for(j=1;j<=numdos;j++) muhat[j] = (alpha + beta*dose[j])*tox[j];

  return;

}


void gamma_find()

{

static double gamma_int[2][10] =
     { {6.908, 2.303, 1.609, 1.204, 0.916, 0.693, 0.511, 0.357, 0.223, 0.105 },
       {0.693, 0.642, 0.588, 0.531, 0.470, 0.405, 0.336, 0.262, 0.182, 0.095 } };
int mod_type;
double tox_at;
int q1, q2, q3, did_break;
int istep, igam, loop;
double gamma_low, gamma_high, gamma_at, g_diff, g2, g1;
double g1lik, g2lik;
int ig0, ig1, ig2, ig3;

  dose_last = dose[numdos];
  mod_type = 0;
  if(model == 1 || model == 5) mod_type = 1;

/*  check for zero gamma and/or model 2  */

  if(f_prime(0.0) >= 0.0 || model == 2) {
	gamma = 0.0;
	parms;
	return;
  }

  tox_at = exp(-gamma*dose_last);
  if(mod_type == 1) tox_at = 2.0 - 1.0/tox_at;
  q1 = (f_prime(gamma) <= 0.0);
  q2 = (gamma*dose_last >= gamma_int[mod_type][0]);
  istep = (f_prime(gamma) <= 0.0) ? -1 : 1;
  if(q1 && q2) {
	gamma = gamma_int[mod_type][0]/dose_last;
	parms();
	return;
  }
  gamma_low = 0.0;
  gamma_high = gamma;
  q3 = (tox_at < 0.9);
  if((!q1 && q2) || (q1 && q3)) {
    igam = (gamma == 0.0) ? 10.0 : 10.0*tox_at + 1;
	if(tox_at < 0.001) igam--;
	if(istep == 1) igam++;
	if(istep > 0) gamma_low  = gamma;
	if(istep < 0) gamma_high = gamma;
	did_break = FALSE;
	while((istep == -1) || (igam <= 10)) {
	  if((istep != 1) && (igam < 1)) {
	    gamma = gamma_int[mod_type][0]/dose_last;
	    parms();
	    return;
	  }
      gamma_at = gamma_int[mod_type][igam-1]/dose_last;
	  if(istep*f_prime(gamma_at) >= 0) {
	    igam += istep;
	    if(istep > 0) gamma_low  = gamma_at;
	    if(istep < 0) gamma_high = gamma_at;
	  } else {
        if(istep > 0) gamma_high = gamma_at;
	    if(istep < 0) gamma_low  = gamma_at;
	    did_break = TRUE;
	    break;  /* This will force an exit from the while loop */
	  }
	}
	if(!did_break) {
	  gamma_low = 0.0;
	  gamma_high = gamma_int[mod_type][9]/dose_last;
	}
  }
  g_diff = (gamma_high - gamma_low)/1597.0;
  ig0 = 0;
  ig1 = 610;
  ig2 = 987;
  ig3 = 1597;
  g2 = gamma_low + ig2*g_diff;
  gamma = g2;
  parms();
  g2lik = nblik();
  g1 = gamma_low + ig1*g_diff;
  gamma = g1;
  parms();
  g1lik = nblik();
  for(loop=1;loop<=13;loop++) {
	if(g2lik >= g1lik) {
	  ig3 = ig2;
	  ig2 = ig1;
	  g2lik = g1lik;
	  ig1 = ig0 + ig3 - ig2;
	  g1 = gamma_low + ig1*g_diff;
	  gamma = g1;
	  parms();
	  g1lik = nblik();
	} else {
      ig0 = ig1;
	  ig1 = ig2;
	  g1lik = g2lik;
	  ig2 = ig0 + ig3 - ig1;
	  g2 = gamma_low + ig2*g_diff;
	  gamma = g2;
	  parms();
	  g2lik = nblik();
	}
  }
  gamma = g2;
  if(g1lik < g2lik) gamma = g1;
  parms();
  return;
}


double f_prime(gamma)

double gamma;

{

double sumwdx, sumwd, sumwx, sumw, sumwdd;
double sumwdx_prime, sumwd_prime, sumwx_prime, sumw_prime, sumwdd_prime;
double tox[11], tj;
double tox_prime[11], tj_prime;
double muhat_prime[11];
double mut[11];
int j;
int flag_alpha;
double alpha_prime, beta_prime, f_prim;
double beta_num, beta_denom;


  sumwdx = 0.0;
  sumwd = 0.0;
  sumwx = 0.0;
  sumw = 0.0;
  sumwdd = 0.0;
  sumwdx_prime = 0.0;
  sumwd_prime = 0.0;
  sumwx_prime = 0.0;
  sumw_prime = 0.0;
  sumwdd_prime = 0.0;

  for(j=1;j<=numdos;j++) {
    if(model != 0 && model != 4) {
  	tox[j] = 2.0 - exp(gamma*dose[j]);
  	tox_prime[j] = dose[j]*(tox[j] - 2.0);
    } else {
      tox[j] = exp(-gamma*dose[j]);
  	tox_prime[j] = -dose[j]*tox[j];
    }
    tj = tox[j];
    tj_prime = tox_prime[j];
    sumwx += (weight[j]*tj*tj)         *(xbar[j]/tj);
    sumw  += (weight[j]*tj*tj);
    sumwx_prime += weight[j]*tj_prime*xbar[j];
    sumw_prime  += weight[j]*tj*tj_prime*2.0;
    if(model < 4) {
   	  sumwdx += (weight[j]*tj*tj)*dose[j]*(xbar[j]/tj);
	  sumwd  += (weight[j]*tj*tj)*dose[j];
  	  sumwdd += (weight[j]*tj*tj)*dose[j]*dose[j];
	  sumwdx_prime += weight[j]*tj_prime*dose[j]*xbar[j];
	  sumwd_prime  += weight[j]*tj*tj_prime*dose[j]*2.0;
	  sumwdd_prime += weight[j]*tj*tj_prime*dose[j]*dose[j]*2.0;
    }
  }

  beta = 0.0;
  if(model < 4) {
	beta_num = sumw*sumwdx - sumwd*sumwx;
	beta_denom = sumw*sumwdd - sumwd*sumwd;
	beta = beta_num/beta_denom;
	if(beta < 0.0 || model >= 3) beta = 0.0;
  }
  alpha = (sumwx - beta*sumwd)/sumw;

  flag_alpha = FALSE;
  if(alpha <= 0.0) {
	flag_alpha = TRUE;
	alpha = 1.0;
	beta = (sumwdx - sumwd)/sumwdd;
  }
  beta_prime = 0.0;
  if(beta != 0.0 && !flag_alpha) {
    beta_prime = ((beta_denom*(sumw*sumwdx_prime +
              sumw_prime*sumwdx - sumwd*sumwx_prime - sumwd_prime*sumwx) -
              beta_num*(sumw*sumwdd_prime + sumw_prime*sumwdd -
              2.0*sumwd*sumwd_prime))/beta_denom)/beta_denom;
  }
  if(flag_alpha) {
	beta_prime = (sumwdd*(sumwdx_prime - sumwd_prime) -
	(sumwdx - sumwd)*sumwdd_prime)/(sumwdd*sumwdd);
  }
  alpha_prime = 0.0;
  if(flag_alpha) {
	alpha_prime = (sumwx_prime - (beta*sumwd_prime + sumwd*beta_prime))/sumw
	- (alpha*sumw_prime)/sumw;
  }
  f_prim = 0.0;
  for(j=1;j<=numdos;j++) {
	mut[j] = alpha + beta*dose[j];
	muhat[j] = mut[j]*tox[j];
	muhat_prime[j] = mut[j]*tox_prime[j] +
	                 tox[j]*(alpha_prime + beta_prime*dose[j]);
	f_prim -= weight[j]*muhat_prime[j]*(xbar[j] - muhat[j]);
  }

  return(f_prim);

}


void weights()

{

  int j;

  for(j=1;j<=numdos;j++) weight[j] = (muhat[j] > 0.0) ?
                                     n[j]/(muhat[j]*(1.0 + c*muhat[j])) :
                                     1.0;
}


double addconx()

{

  double sum;
  int i, j;

  sum = 0.0;
  for(j=1;j<=numdos;j++) {
	for(i=1;i<=n[j];i++) sum -= dlgama((double) kount[j][i] + 1.0);
  }
  return(sum);
}


double nblik()

{

  int j;
  double cj, sum, ratioj;

/* check for violated constraints.  if found, impose penalty. */

  sum = 0.0;
  for(j=1;j<=numdos;j++) {
	if(muhat[j] > 0.0) {
	  cj = 1.0 + c*muhat[j];
	  ratioj = (c != 0.0) ? log(cj)/c : muhat[j];
	  sum -= n[j]*(xbar[j]*log(muhat[j]/cj) - ratioj);
	} else {
	  if(xbar[j] != 0.0) return(-7777.0);
	}
  }
  return(sum);
}


double addcon(xsum,nsum)

double xsum;
int nsum;

{

  double sum, cinv;
  int i, j;

  sum = adtvcnx;
  if(c == 0.0) return(sum);
  cinv = 1.0/c;
  for(j=1;j<=numdos;j++) {
	for(i=1;i<=n[j];i++) sum += dlgama((double) kount[j][i] + cinv);
  }
  sum -= (nsum*dlgama(cinv) - xsum*log(c));
  return(sum);
}


double cmputc()

{

  double c;
  double ci;
  double z;
  double z2;
  double z2sum;
  int i, j;

  c = 0.0;
  z2sum = 0.0;
  for(i=1;i<=numdos;i++) {
	if(xbar[i] != 0.0 && n[i] > 1) {
	  ci = 0.0;
	  z = xbar[i];
	  z2 = z*z;
	  z2sum += z2;
	  for(j=1;j<=n[i];j++) ci += (z2 != 0.0) ?
	                             kount[i][j]*((double) kount[i][j] - 1)/z2 : 0;
	  ci = (ci + 1.0/z - n[i])/(n[i]-1);
	  c += ci*z2;
	}
  }
  c = c/z2sum;
  return((c < 1.0e-8) ? 0.0 : c);
}


double dlgama(z)

double z;

{

  double zfac;
  double s1;
  double s2;
  double x;
  double dlgam;
  int i;
  static double a[7] = {0.833333333333333e-1,  0.333333333333333e-1,
	                    0.252380952380952e0,  0.525606469002695e0,
	                    0.101152306812684e1,  0.151747364915329e1,
	                    0.226948897420496e1 };
  static double dcon = 0.918938533204673e0;
  static double plgam3[8] = {-2.407986980173375e8, -2.448321769032882e8,
	                          4.818077102773628e8,  1.119385354299855e8,
	                         -8.731675438238387e7, -1.971830115860921e7,
	                         -9.827102281420492e5, -1.037701651732974e4 };
  static double plgam2[8] = {-1.513831834115067e3, -1.508630228766725e4,
	                         -2.064829420532528e4,  1.204317380987164e4,
	                          1.955360554063045e4,  5.268983255914981e3,
	                          3.775106797972170e2,  5.155057617640817e0 };
  static double plgam1[8] = {-2.208843997216182e0, -7.741064071332953e1,
	                         -5.176383498023218e2, -9.222613728801522e2,
	                         -2.617218583856145e2,  2.431752435244210e2,
	                          8.568982062831317e1,  4.120843185847770e0 };
  static double qlgam3[8] = {-7.902611114187634e7, -4.353707148043742e8,
	                         -4.044359282914355e8, -1.119254116263318e8,
	                         -1.048577583049937e7, -3.114062847340679e5,
	                         -2.015275195500483e3,  1.000000000000000e0 };
  static double qlgam2[8] = { 6.983274140573510e2,  1.440209037170085e4,
	                          5.262286383841199e4,  5.712025539602503e4,
	                          2.202956214415664e4,  3.039903041439440e3,
	                          1.289093189012958e2,  1.000000000000000e0 };
  static double qlgam1[8] = { 4.097792921092615e-1, 2.443519662506312e1,
	                          2.623083470269460e2,  8.460755362020782e2,
	                          9.513235976797060e2,  3.778372484823942e2,
	                          4.564677187585908e1,  1.000000000000000e0 };

/* the interval [12,infinity] . . . */

  if(z >= 12.0) {
    dlgam = (z - 0.5)*log(z) + dcon - z;
	x = a[6]/z;
	for(i=5;i>=0;i--) x = a[i]/(z + x);
	return(dlgam + x);

  } else if(z >= 4.0) {

/* the interval [4.0,12.0] . . . */

    s1 = plgam3[7]*z + plgam3[6];
	s2 = qlgam3[7]*z + qlgam3[6];
	for(i=5;i>=0;i--) {
	  s1 = s1*z + plgam3[i];
 	  s2 = s2*z + qlgam3[i];
	}
	return(s1/s2);

  } else if(z >= 1.5) {

/* the interval [1.5,4.0] . . . */

	s1 = plgam2[7]*z + plgam2[6];
	s2 = qlgam2[7]*z + qlgam2[6];
	for(i=5;i>=0;i--) {
	  s1 = s1*z + plgam2[i];
	  s2 = s2*z + qlgam2[i];
	}
	return((z - 2.0)*(s1/s2));
	
  } else {
	
/* the interval [0.0,1.5]: for z in [0.0,0.5] we evaluate instead the
   function at z+1 and the use the recursion mentioned above . . .  */

	zfac = z + (z < 0.50) ? 1.0 : 0.0;
	s1 = plgam1[7]*zfac + plgam1[6];
	s2 = qlgam1[7]*zfac + qlgam1[6];
	for(i=5;i>=0;i--) {
	  s1 = s1*zfac + plgam1[i];
	  s2 = s2*zfac + qlgam1[i];
	}
	return((z >= 0.5) ? (z - 1.0)*(s1/s2) : -log(z) + z*(s1/s2));
  }
	
}

double chicdf(chi,idegre)

double chi;
int idegre;
{

  #define PIE 3.1415926536

  int i, index;
  double product, sum, rootch, atemp, zxsum, dlnzx, chisum;
  double return_value;

  if(idegre > 100) {
    return_value = ztail(sqrt(2.0*chi) - sqrt((double) (2.0*idegre-1)));
    return(return_value);
  }
  product = 1.0;
  sum = 0.0;
  rootch = sqrt(chi);
  atemp = sqrt(2.0*PIE);
  if((idegre%2) != 0) {
    zxsum = 0.0;
    if(idegre == 1) {
      return_value = 2.0*(ztail(rootch));
      return(return_value);
    }
    dlnzx = -0.5*chi - log(atemp);
    for(i=2;i<=idegre-1;i+=2) {
      product *= (double) (i-1);
      sum += pow(rootch,(double) (i-1))/product;
    }
    sum = dlnzx + log(sum);
    zxsum = (sum > -85.0) ? exp(sum) : 0.0;
    return_value = 2.0*(ztail(rootch) + zxsum);
    return(return_value);
  }
  index = (idegre-2)/2;
  dlnzx = -0.5*chi;
  chisum = 1.0;
  for(i=1;i<=index;i++) {
    product *= (double) i;
    chisum += pow((0.5*chi),(double) i)/product;
  }
  sum = dlnzx + log(chisum);
  zxsum = (sum > -85.0) ? exp(sum) : 0.0;
  return_value = zxsum;
  return(zxsum);
}

double ztail(z)

double z;

{

  double a[6],t,x,sum;
  int i,isign;

  a[0] = 0.0;
  a[1] = 0.254829592;
  a[2] = -0.284496736;
  a[3] = 1.421413741;
  a[4] = -1.453152027;
  a[5] = 1.061405429;

  isign = (z < 0.0) ? -1 : 1;
  x = fabs(z)/sqrt((double) 2.0);
  if(x > 6.0) {
    return((isign == 1) ? 0.0 : 1.0);
  }
  t = 1.0/(1.0 + 0.3275911*x);
  sum = 0.0;
  for(i=1;i<=5;i++) {
    sum += a[i]*pow(t,(double) i);
  }
  return(0.5 - isign*0.5*(1.0 - sum*exp(-x*x)));
}


void fisher(covar)

double covar[3][3];

{

int i, j, l;
double tox, prime[3];

  covar[0][0] = 0.0;
  covar[1][0] = 0.0;	
  covar[1][1] = 0.0;
  covar[2][0] = 0.0;
  covar[2][1] = 0.0;
  covar[2][2] = 0.0;

  for(j=1;j<=numdos;j++) {
    xmut = alpha + beta*dose[j];
    tox = 1.0;
    if(model == 0 || model == 4) {
      tox = exp(-gamma*dose[j]);
    } else if(model == 1 || model == 5) {
      tox = 2.0 - exp(gamma*dose[j]);
    }

    prime[0] = tox;
    prime[1] = dose[j]*tox;
    prime[2] = 1.0;
    if(model == 0 || model == 4) {
      prime[2] = -xmut*dose[j]*tox;
    } else if(model == 1 || model == 5) {
      prime[2] = -xmut*dose[j]*(2.0 - tox);
    }

    for(i=0;i<=2;i++) {
      for(l=0;l<=i;l++) {
        covar[i][l] += weight[j]*prime[i]*prime[l];
      }
    }
  }

  if(model <= 1) return;
  if(model <= 3) {
    covar[2][0] = 0.0;
    covar[2][1] = 0.0;
    covar[2][2] = 1.0;
    if(model == 2) return;
  }
  covar[1][0] = 0.0;
  covar[1][1] = 1.0;
  covar[2][1] = 0.0;
  return;
}



void muhatc()

{

int j;
double tox;

  for(j=1;j<=numdos;j++) {
    tox = 1.0;
    if(model == 0 || model == 4) {
      tox = exp(-gamma*dose[j]);
    } else if(model == 1 || model == 5) {
      tox = 2.0 - exp(gamma*dose[j]);
    }
    muhat[j] = (alpha + beta*dose[j])*tox;
  }

  return;

}


int inverse(matrix)

double matrix[3][3];

{

  double x[3][3], y[3][3], temp;

  if(matrix[0][0] <= 0.0) return(1);

  x[0][0] = sqrt(matrix[0][0]);
  x[1][0] = matrix[1][0]/x[0][0];
  temp = matrix[1][1] - x[1][0]*x[1][0];
  if(temp <= 0.0) return(2);

  x[1][1] = sqrt(temp);
  x[2][0] = matrix[2][0]/x[0][0];
  x[2][1] = (matrix[2][1] - x[2][0]*x[1][0])/x[1][1];
  temp = matrix[2][2] - (x[2][0]*x[2][0] + x[2][1]*x[2][1]);
  if(temp <= 0.0) return(3);

  x[2][2] = sqrt(temp);

  y[0][0] = 1.0/x[0][0];
  y[1][0] = -x[1][0]/(x[0][0]*x[1][1]);
  y[1][1] = 1.0/x[1][1];
  y[2][0] = (x[1][0]*x[2][1] - x[2][0]*x[1][1])/(x[0][0]*x[1][1]*x[2][2]);
  y[2][1] = -x[2][1]/(x[1][1]*x[2][2]);
  y[2][2] = 1.0/x[2][2];

  matrix[0][0] = y[0][0]*y[0][0] + y[1][0]*y[1][0] + y[2][0]*y[2][0];
  matrix[1][0] = y[1][0]*y[1][1] + y[2][0]*y[2][1];
  matrix[1][1] = y[1][1]*y[1][1] + y[2][1]*y[2][1];
  matrix[2][0] = y[2][0]*y[2][2];
  matrix[2][1] = y[2][1]*y[2][2];
  matrix[2][2] = y[2][2]*y[2][2];
  return(0);
}