program MaximumLikelihood (input, output, sibships, Analysis);
	{2/1/90         Copyright 1990 Ruth G. Shaw
			Department of Botany and Plant Sciences
			University of California
			Riverside, CA 92521

        This program employs maximum likelihood or restricted
        maximum likelihood to estimate genetic parameters for pedigrees
        containing parents, offspring, halfsibs, and fullsibs.
        NOTE: THIS VERSION WILL RUN ON CRAY UNDER UNICOS} 

const r=2; {r=the max. number of traits}
      nvarcov = 3; {the number of (co)variances of a given type; r*(r+1)/2}
      maxnparm=9; {the max. number of params nvarcov * 3, generally}  
      f=20;  {the max number of families}
      sibs=13; {max # of individuals in a family, including unmeasured}
      rhsibs=18; {max # of recip half sibs in a family}
      rfsibs=6; {max # of recip full sibs in a family}
      missing=-98; {Value in the dataset indicating missing phenotypic data}
      pos=26; {maximum number of observations in a family}
      nmax=260; {maximum total # of individuals in the dataset, including
                 those missing data}
      delta = 0.01;  {Convergence criterion}
      maxround=10;
      df = 1;  {the max. number of fixed effects}
      p = 2;  {the maximum number of levels of fixed effects=r~nlevels}
      ceo = 1; {common environment option..4 for c.e, 1 for dominance}



 type list = array[1..df] of integer;
      Observation=record
                   ID:integer;
                   Father: integer;
                   Mother: integer;
                   NPS, NMS, FO, Famid: integer;
		   FixedVals: list;
                 end;  {Observation}
     Dataset=array[1..nmax] of Observation;
     FamZset = array[1..sibs,1..r] of real;
     FamZmat = array[1..f] of FamZset;
     slevelsbool = array[1..p] of boolean;
     Famset=array[1..sibs] of Observation;
     Famarray=array[1..f] of Famset;
     Paramvar=array[1..r,1..r] of real;
     Parammean=array[1..p] of real;
     Parmvec=array[1..maxnparm] of real;
     Parmbool=array[1..maxnparm] of boolean;
     info = array [1..maxnparm,1..maxnparm] of real;
     Vector=array[1..pos] of real;
     FamVector = array[1..f] of Vector;
     Covrow = array[1..pos] of real;
     Covmatrix =  array[1..pos] of Covrow;
     Covarray = array[1..f] of Covmatrix;
     Desmatrix = array [1..pos,1..p] of integer;
     Desarray = array [1..f] of Desmatrix;
     Measures = array[1..nmax,1..r] of real;
     Famphens = array[1..pos,1..p] of real;
     Size=array[1..f] of integer;
     Traitsquare = array[1..p,1..p] of real;
     items = array[1..df,1..p] of integer;
     arow = array[1..pos] of integer;
     amatrix = array[1..rhsibs] of arow;
     dmatrix = array[1..rfsibs] of arow;
     ematrix = array[1..1] of arow;
     Desset = array[1..maxnparm] of Famphens;
     tderivsa = array[1..nvarcov] of amatrix;
     tderivsd = array[1..nvarcov] of dmatrix;
     tderivse = array[1..nvarcov] of ematrix;
     derivsa = array[1..f] of tderivsa;
     derivsd = array[1..f] of tderivsd;
     derivse = array[1..f] of tderivse;
     Traitset = array[1..maxnparm] of Traitsquare;
     countrow = array[1..pos] of integer;
     tcountset = array[1..nvarcov] of countrow;
     countset = array[1..f] of tcountset;

var  PresentData: Dataset;
     Design : Desarray;
     CU : Desset;
     Ldet :real;
     LogLike, OldLogLike, Firstloglike :real;
     XVY,WtdVals, Means:Parammean;
     AddCov, DomCov, EnvCov :Paramvar;
     D, eD, SumD, Parms, OldParms, BParms, First:Parmvec;
     Zeroed3, Zeroed2, Zeroed : Parmbool;
     better, Done, REML, ok, Maxed, Trunked: boolean;
     bfl, TA, Trunked2, feasable : boolean;
     flagvec: slevelsbool;
     sibships: text;
     Analysis: text;
     nchar, ntrunc, slevels, nfxfact, nvc: integer;
     newfam, I, Iter: integer;
     Families: Famarray;
     FamZ: FamZmat;
     AB, VInvX: Famphens;
     Py: FamVector;
     Addcoeff, Domcoeff, Envcoeff, Matcoeff, Patcoeff, Kcoeff: real;
     M, famnum, vartypes : integer; {tot #, fam#, #types of comps estimable}
     Zmatrix:Measures;
     Q: Covarray;
     Vec, DiffVec, MeanVec:  Vector;
     famsize, vert: Size;
     Empty  : info;
     sBcov, XVXInv, XVX, Rfix : Traitsquare;
     UpCU, SumUpCU :  Traitset;   
     B, C, eB, Big, RBig : info;
     nfxlevels: list; {contains the true number of levels for each fixed factor}
     levels: items; {Contains the values of the levels as in the data}
     DataVec: FamVector;
     tdva, tdvd, tdve : countset;
     DVA : derivsa;
     DVD : derivsd;
     DVE : derivse;


procedure Getphtypes(var Data:Dataset; var M: integer; var nchar:integer;
                     var Zmatrix: Measures; var slevels: integer;
		     var nfxlevels: list);

          {read the data from the input file. Each line of data is a record
          for an individual in the pedigree, containing, in order, an  ID
          number, the sire's ID, the dam's ID, the ID of the next paternal sib,
          and that of the next maternal sib
          (the values of relatives not present in the pedigree are assumed
          to be set at 0), followed by phenotypic values for an arbitrary
          number of traits in a consistent order.  (The value of a
          trait that is missing for a given individual is assumed to be -99.
          Both parents of individuals for which phenotypic data is available
          must be included in the pedigree, regardless of whether
          phenotypic information is available for them.)}

   var  TempID, TempFather, TempMother: integer;
        Pheno: real;
        J, K, L, which, count :integer;
	already: boolean;
   begin
      slevels := 0;
      for K := 1 to nfxfact do begin
         nfxlevels[K] := 0;
         for L := 1 to p do
 	    levels[K,L] := -99;
	 end;
      M := 0;
      read (sibships,TempID);
      while not (TempID=0) do begin
         M :=M+1;
         if M > nmax 
         then begin
          writeln ('  The number of data records exceeds "nmax" ');
          halt;
         end;
         read (sibships, TempFather,TempMother);
         with Data[M] do begin
            ID := TempID ;
	    Father:=TempFather;
	    Mother:=TempMother;
            NPS := 0;
            NMS := 0;
            FO := 0;
            Famid := 0;
            for K := 1 to nfxfact do begin
               read (sibships,which);
               already := false;
               count := 0;
               for L := 1 to nfxlevels[K] do begin
                  count := count + 1;
                  if (which=levels[K,L])
                  then begin
                    already := true;
                    FixedVals[K] := count;
                  end;
                end; {for L}
               if not already
               then begin
                  nfxlevels[K]:= nfxlevels[K] + 1;
                  FixedVals[K]:= nfxlevels[K];
                  write (Analysis, which,'is fixed effect number',nfxlevels[K]);
                  writeln (Analysis, 'of factor ', K);
                  slevels := slevels + 1;
                  levels[K,nfxlevels[K]] := which;
                  end;
               end; {K}
            end; {with}
         for J := 1 to nchar do begin
            read (sibships, Pheno);
            Zmatrix[M,J]:=Pheno;
            end;
         readln (sibships);
         read (sibships,TempID);
      end; {while not}
   readln (sibships);
   if (nfxfact < 1)
   then slevels := 1;
   if (nfxfact > 1)
   then slevels := slevels - (nfxfact - 1); 
   slevels := nchar * slevels;
   writeln ('slevels ', slevels);
   end; {Getphtypes}

procedure FindGroupsN (var Data:Dataset; var M: integer);
     {The procedure fills in the pedigree information for
     each individual:  it finds first offspring, next maternal
     sib and next paternal sib.  If none exist, the value is 0.}

   var I, J, Dad, Mom, Me : integer;
       matched1, matched2 : boolean;
   begin
      for I := 1 to M do begin
         Dad := Data[I].Father;
         Mom := Data[I].Mother;
         Me := I;
         if Dad > 0 
	 then begin
           if Data[Dad].FO = 0
           then Data[Dad].FO := Me;
           if I < M then begin
  	   J := I;
	    matched1 := false;
	    repeat 
	       J := J + 1;
               if (Data[J].Father = Dad) or (Data[J].Mother = Dad)
	       then begin
	          Data[I].NPS := J;
	          matched1 := true;
	          end;
	    until matched1 or (J>M-1);
	    end;
	   end;
         if Mom > 0 
	 then begin
            if Data[Mom].FO = 0
            then Data[Mom].FO := Me;
            if I < M then begin
  	    J := I;
	    matched2 := false;
            repeat 
	       J := J + 1;
               if (Data[J].Father = Mom) or (Data[J].Mother = Mom)
	       then begin
	          Data[I].NMS := J;
	          matched2 := true;
	          end;
	    until matched2 or (J>M-1);
            end;
           end;
         end;
      end;

procedure Visit (var rel:Observation; family:integer);

      {Family designations are assigned to each individual, using all the
      information of relationships.}

   begin
      if rel.Famid = 0
      then begin
         rel.Famid := family;
         if (rel.FO > 0)
         then Visit (PresentData[rel.FO], family);
         if (rel.NPS > 0)
         then Visit (PresentData[rel.NPS], family);
         if (rel.NMS > 0)
         then Visit (PresentData[rel.NMS], family);
         if (rel.Father > 0)
         then Visit (PresentData[rel.Father], family);
         if (rel.Mother > 0)
         then Visit (PresentData[rel.Mother], family);
         end;
      end; {Visit}

procedure Check (I:integer; Data:Dataset; var y : integer);
   var Dad, Mom:integer;
   begin
      Dad := Data[I].Father;
      if Dad > 0
      then y := Data[Dad].Famid;
      if (y=0)
      then begin
         Mom := Data[I].Mother;
         if Mom > 0
         then y := Data[Mom].Famid;
	 end;
      end;

procedure PrintFams (var Data: Dataset; Zmatrix:Measures);
   {All the pedigree information is written to the Analysis file}
   var I, J: integer;
   begin
     writeln (Analysis,' ID   FA   MO   NPS  NMS FamID  DATA  ');
      for I := 1 to M do begin
         with Data[I] do begin
            write (Analysis, ID:5, Father:5, Mother:5, NPS:5);
            write (Analysis, NMS:5, Famid:5);
            end;
            for J:= 1 to nchar do
	      write (Analysis, Zmatrix[I,J]:8:2);
              writeln(Analysis);
      end;
   end; {PrintFams}

procedure Collecthem (var Data: Dataset; var Families: Famarray; 
                      var FamZ:FamZmat; var famsize: Size);

        {Two arrays are built for each family.  One contains the pedigree
        info. (Families) and the other contains the phenotypes (FamZ).}

   var I, J, K, L : integer;
   begin {Collecthem}
      for I := 1 to famnum do begin
         L:=1;
         for K := 1 to M do begin
            if Data[K].Famid=I
            then begin
               Families[I][L]:=Data[K];
               for J := 1 to nchar do
                   FamZ[I][L,J] := Zmatrix[K,J];
               L := L+1;
               end;
            end;
         famsize[I] := L-1;
         end;
       end; {Collecthem}

procedure CoeffMatrices (var Addcoeff, Domcoeff, Envcoeff : real;
                         I, L, K: integer);

          {The matrices of coefficients of the causal components of variance
          (Additive, Dominance, Environmental) are established for each
          family on the basis of the pedigree information.}

   var valA, valD, valE : real;
       SameFather, SameMother, Parent, Progen, Opposite1, Opposite2: boolean;
   begin {CoeffMatrices}
      if L=K
      then begin
         valA := 1.0;
         valD := 1.0;
         valE := 1.0;
         end  {Sameind}
      else begin
         Parent := false;
         Progen := false;
         SameFather := false;
         SameMother := false;
         valA := 0;
         valD := 0;
         valE := 0;
         if (Families[I][K].Father>0)
         then begin
 	    SameFather:=(Families[I][K].Father=Families[I][L].Father);
 	    Opposite1 :=(Families[I][K].Father=Families[I][L].Mother);
 	    end;
         if (Families[I][K].Mother>0)
         then begin
 	    SameMother:=(Families[I][K].Mother=Families[I][L].Mother);
 	    Opposite2 :=(Families[I][K].Mother=Families[I][L].Father);
 	    end;
         if SameMother or SameFather or Opposite1 or Opposite2
         then begin
            valA:=0.25;
            if (SameMother and SameFather) or (Opposite1 and Opposite2)
            then begin
               valA:= 0.5;
               valD:=0.25 * ceo; {common environment option}
               end;
            end {SameFather or SameMother or... or...}
         else begin
            Parent := (Families[I][L].Father=Families[I][K].ID) or
                            (Families[I][L].Mother=Families[I][K].ID);
            Progen := (Families[I][K].Father=Families[I][L].ID)
                            or (Families[I][K].Mother=Families[I][L].ID);
            if Parent or Progen
            then valA:=0.5;
            end; {else (not SameFather or SameMother)}
         end; {else (not Samind)}
      Addcoeff:=valA;
      Domcoeff:=valD;
      Envcoeff:=valE;
      end; {CoeffMatrices}

procedure Designit (FamZ: FamZmat; var vert:Size;
                    var Design:Desarray; var DataVec: FamVector);
       {Produces the arrays Design(vert x nchar) (design matrix for the fixed
       effects in each family) and DataVec(g*vert x 1) (The nonmissing
       observations in the family).}
    var I, J, K, L, g: integer;
        vertical, col, item: integer;
    begin
       for g := 1 to famnum do begin
          vertical:=0;
          for I := 1 to pos do begin
             for J := 1 to slevels do
               Design[g][I,J] := 0;
             end;
          for I:=1 to famsize[g] do begin
            for J:=1 to nchar do begin
               if not (FamZ[g][I,J]<missing)
                then begin
                  vertical:=vertical+1; 
                  if vertical > pos
                  then begin
                   writeln ('The number of observations for family ',g:5);
                   writeln ('exceeds "pos"');
                   halt;
                  end;
 		  DataVec[g][vertical]:=FamZ[g][I,J];
 		  Design[g][vertical,J] := 1;
 		  item := 0;
 		  for L := 1 to nfxfact do begin
 		     item := Families[g][I].FixedVals[L]-1;
 		     col := nchar;
		     for K := 1 to (L-1) do
 		        col:= col + nchar*(nfxlevels[K]-1);
                     if item > 0
 		     then begin
		        col := col + (item-1)*nchar + J;
 		        Design[g][vertical,col] := 1;
			end;
 		     end; {L}
 		  end; {if}
 	       end; {J}
             end; {I}
 	 vert[g] := vertical;
 	 end; {g}
       end;

procedure ZeroEmpty (var Empty: info);
  {The empty matrix is used for zeroing}
   var i,j: integer;
   begin
     for i := 1 to maxnparm do begin
        for j := 1 to maxnparm do
           Empty[i,j] := 0;
        end;
     end;

procedure EnterValues(var AddCov, DomCov, EnvCov :Paramvar);

         {values for the variance components are transferred
         into the matrices of variance and covariance components
         of the traits.}

   var I,J: integer;
       add, group, loc, top : integer;
   begin
      group := (nchar *(nchar+1)) div 2;
      top := 3*group;
      for I := 1 to top do
	 write (Analysis, Parms[I]:10:2);
      writeln (Analysis);
      add := 0;
      for I:=1 to nchar do begin
         for J:=I to nchar do begin
            add := add + 1;
            AddCov[I,J]:= Parms[add];
            AddCov[J,I] := AddCov[I,J];
            loc := group + add;
            EnvCov[I,J] := Parms[loc];;
            EnvCov[J,I] := EnvCov[I,J];
            loc := loc + group;
            DomCov[I,J]:= Parms[loc];
            DomCov[J,I]:=DomCov[I,J];
            end; {J loop}
         end;
      end;   {EnterValues}

procedure reMakeCovArray (var Q:Covarray; var Parms:Parmvec;
			      var DVA:derivsa; var DVD:derivsd;
			      var DVE:derivse; 
			      var tdva, tdvd, tdve :
			      countset);

var i,j,k,g,index,l : integer;
var fact: real;

      begin
       for g := 1 to famnum do begin
         for i := 1 to vert[g] do begin
            for k := 1 to vert[g] do
            Q[g][i][k] := 0.0;
       end; {i}
	   for j := 1 to nvc do begin
        if ((not Trunked) or (not Zeroed2[j])) then begin
         for i := 1 to vert[g] do begin
		 for l := 1 to tdva[g][j][i] do begin
                    index := DVA[g][j][l][i] mod 10000;
		    fact := trunc(DVA[g][j][l][i]/10000);
		    Q[g][i][index] := Q[g][i][index] + Parms[j]/fact;
		    end; {l}
               end; {i}
        end; {if not Trunked}
        if ((not Trunked) or (not Zeroed2[j+2*nvc])) then begin
         for i := 1 to vert[g] do begin
		 for l := 1 to tdvd[g][j][i] do begin
                    index := DVD[g][j][l][i] mod 10000;
		    fact := trunc(DVD[g][j][l][i]/10000);
		    Q[g][i][index] := Q[g][i][index] + Parms[j+2*nvc]/fact;
		    end; {l}
               end; {i}
        end; {if not Trunked}
        if ((not Trunked) or (not Zeroed2[j+nvc])) then begin
         for i := 1 to vert[g] do begin
		 for l := 1 to tdve[g][j][i] do begin
                    index := DVE[g][j][l][i];
		    Q[g][i][index] := Q[g][i][index] + Parms[j+nvc];
                    end; {l}
               end; {i}
        end; {if not Trunked}
        end; {j}
     end; {g}
   end; {remake}

procedure MakeCovArray (var EstCovArray:Covarray;
		        var DVA:derivsa; var DVD:derivsd; var DVE:derivse; 
			      var tdva, tdvd, tdve : countset);

      {The covariance matrix is made for each family.  The procedure
      produces a matrix that is (vert*nchar-missing) long, where vert
      is the # of measured individuals in the family, nchar is the number of
      traits, and missing is the # of data values that are missing.}
      {the following is experimental and good for only one family data}

   type sextet = array[1..6] of real;
        codeset= array[1..6] of integer;
   var G, I, J, K, L, skippedi, skippedj: integer;
       a, X, Y, s, t, u:integer;
       coeff: sextet;
       intcoeff: codeset;
   begin
      for G:= 1 to famnum do begin
      skippedi:=0;
         for I:= 1 to famsize[G] do begin
            for K:= 1 to nchar do begin
               if (FamZ[G][I,K]<missing)
               then skippedi := skippedi + 1
               else begin
                  X:= (I-1) * nchar + K - skippedi;
                  skippedj:=0;
                  for J:= 1 to famsize[G] do begin
                     for L:= 1 to nchar do begin
                        if (FamZ[G][J,L]<missing)
                        then skippedj := skippedj + 1
                        else begin
                           Y:= (J-1) * nchar + L - skippedj;
                                 CoeffMatrices(coeff[1],coeff[2],coeff[3],
				   G,I, J);
                           EstCovArray[G][X][Y]:=coeff[1] * AddCov[K,L] +
                                           coeff[2] * DomCov[K,L] +
                                           coeff[3] * EnvCov[K,L];
                        u := 0;
                        for s := 1 to nchar do begin
                           for t := s to nchar do begin
                              u := u + 1;
                              if (K=s) and (L=t) or ((K=t) and (L=s))
                              then begin
                                 for a := 1 to 2 do begin
				    if (coeff[a]>0)
				    then begin
                                      if coeff[a] <0.26
                                         then intcoeff[a] := 40000
                                         else if coeff[a] < 0.51
                                              then intcoeff[a] := 20000
                                              else if coeff[a] > 0.99
                                                  then intcoeff[a] := 10000;
                                      end
				    else intcoeff[a] := 0;
                                    end;           {a}
				 if (intcoeff[1]>0)
				 then begin
				    tdva[G][u][X] := tdva[G][u][X] + 1;
                                    DVA[G][u][tdva[G][u][X]][X]:=intcoeff[1]+Y;
				    end;
				 if (intcoeff[2]>0)
				 then begin
				    tdvd[G][u][X] := tdvd[G][u][X] + 1;
                                    DVD[G][u][tdvd[G][u][X]][X]:=intcoeff[2]+Y;
				    end;
				 if (coeff[3]>0)
				 then begin
				    tdve[G][u][X] := tdve[G][u][X] + 1;
                                    DVE[G][u][tdve[G][u][X]][X]:=Y;
				    end;
                                 end {if}
                              end;{t loop}
                             end; {s loop}
                           end; {else}
                        end; {L loop}
                     end; {J loop}
                  end; {else}
               end; {K loop}
            end; {I loop}
         end; {G loop}
end; {MakeCovMatrix}

procedure CheckValues (GenCov:Paramvar; var ok: boolean);
   {Checks to make sure the phenotypic variances are positive.}
   var Varsum{, CorrA,CorrE}: boolean;
       I{,J}: integer;
   begin
      for I:=1 to nchar do begin
         Varsum:=(GenCov[I,I]+DomCov [I,I]+EnvCov[I,I])>0; 
         if not Varsum
         then begin
            ok := false;
            writeln ({Analysis,}'The var. of trait',I,'is less than 0.');
            end;
         end;
      end; {CheckValues}

procedure ComputeLogLikelihood (var LogLike: real;
                                 var Q: Covarray; var Means: Parammean;
                                 var XVXInv: Traitsquare;  var Py : FamVector;
                                 var vert: Size;  var DataVec: FamVector;
				 var Design: Desarray; var ok:boolean);
   var g,I,J: integer;
       Matprod, SumLdet, XVXdet: real;
       SumS: real;

   procedure GetWtdVals (var WtdVals: Parammean; var VInvX: Famphens;
                         upper: integer; DataVec: Vector; Des: Desmatrix);
      {WtdVals: X'VInv y}

      var i, j, k :integer;
          tally : real;
      begin
         for i := 1 to upper do begin
            for j:= 1 to slevels do begin
               tally := 0;
               for k:= 1 to upper do
                  tally := Q[g][i][k] * Des[k,j] + tally;
               VInvX[i,j] := tally;
               end;
            end;
         for i := 1 to slevels do begin
            WtdVals[i] := 0;
            for j := 1 to upper do
               WtdVals[i] := VInvX[j,i] * DataVec[j] + WtdVals[i];
            end;
         end;

   procedure Vechs (var MeanVec: Vector; Means: Parammean;
                    Des: Desmatrix; upper: integer);
      {Vechs produces MeanVec, a vector of the means corresponding to the
      measures in DataVec.}
      var i, j : integer;
      begin
         for i := 1 to upper do begin
	    MeanVec[i] := 0.0;
            for j := 1 to slevels do begin
               if Des[i,j] = 1
               then MeanVec[i] := Means[j] + MeanVec[i];
               end;
            end;
         end;

procedure Standardize (DataVec: Vector; MeanVec: Vector; upper: integer;
                          var DiffVec: Vector);
      {DiffVec is the deviations of the measures from their means.}
      var I : integer;
      begin
         for I:=1 to upper do begin
            DiffVec[I]:=DataVec[I]-MeanVec[I];
            end;
         end;

procedure MatrixInversion (var A:Covmatrix; var Ldet: real; n: integer;
                           var ok: boolean);
   {Invert by triangularization. Return the inverse.}

   procedure Rootit (var A: Covmatrix; var Ldet: real; var ok: boolean);
     {The square root.}
      var i, j, k: integer;
          sum: real;
      begin
        if (A[1][1] < 0)
        then begin
	   ok := false;
           writeln (Analysis,'The variance matrix is not pos. def.');
	   end
        else begin
           A[1][1] := sqrt(A[1][1]);
           for i := 2 to n do
              A[i][1] := A[i][1]/A[1][1];
           for i := 2 to n do begin
              sum := 0.0;
              for j := 1 to (i-1) do
                 sum := sum +sqr(A[i][j]);
              if (A[i][i] < sum)
              then ok := false
              else begin
                 A[i][i] := sqrt(A[i][i]-sum);
                 if (i<n) then begin
                   for j := (i+1) to n do begin
                     sum := 0.0;
                     for k := 1 to (i-1) do
                       sum := sum + A[i][k] * A[j][k];
                       A[j][i] := (A[j][i]-sum)/ A[i][i];
                     end;
                   end;
                 end;
              end;
            end;
         Ldet := 0;
         for i := 1 to n do
            Ldet := Ldet + ln(A[i][i]);
         Ldet := 2*Ldet;
         end;

   procedure Invertit (var A: Covmatrix);
      {The inverse of the root, followed by L'*L.}
      var i, j, k : integer;
          one : real;
          sum : real;
      begin
         one := 1.0;
         for i := 1 to n do begin
            A[i][i] := one/A[i][i];
            for j := (i+1) to n do begin
               sum := 0.0;
               for k := i to (j-1) do
                  sum := sum + A[j][k]*A[k][i];
               A[j][i] := -sum/A[j][j];
               end;
            end;
         for i := 1 to n do begin
            for j := i to n do begin
               sum := 0.0;
               for k := j to n do
                  sum := sum + A[k][i] * A[k][j];
               A[i][j] := sum;
               end;
            end;
            for i := 1 to n do begin
             for j := i to n do begin
             A[j][i] := A[i][j];
             end;
            end;
         end;

      begin {MatrixInversion}
	 ok := true;
         Rootit(A, Ldet, ok);
         Invertit(A);
      end;

procedure sMatrixInversion (var A:Traitsquare; var Ldet: real; n: integer;
                           var ok: boolean);
   {Invert by triangularization. Return the inverse.}
   procedure sRootit (var A: Traitsquare; var Ldet: real; var ok: boolean);
     {The square root.}
      var i, j, k: integer;
          sum: real;
      begin
        if (A[1][1] < 0)
        then begin
	   ok := false;
           writeln (Analysis,'The XVX matrix is not pos. def.');
	   end
        else begin
           A[1][1] := sqrt(A[1][1]);
           for i := 2 to n do
              A[i][1] := A[i][1]/A[1][1];
           for i := 2 to n do begin
              sum := 0.0;
              for j := 1 to (i-1) do
                 sum := sum +sqr(A[i][j]);
              if (A[i][i] < sum)
              then ok := false
              else begin
                 A[i][i] := sqrt(A[i][i]-sum);
                 if (i<n) then begin
                   for j := (i+1) to n do begin
                     sum := 0.0;
                     for k := 1 to (i-1) do
                       sum := sum + A[i][k] * A[j][k];
                       A[j][i] := (A[j][i]-sum)/ A[i][i];
                     end;
                   end;
                 end;
              end;
            end;
         Ldet := 0;
         for i := 1 to n do
            Ldet := Ldet + ln(A[i][i]);
         Ldet := 2*Ldet;
         end;

   procedure sInvertit (var A: Traitsquare);
      {The inverse of the root, followed by L'*L.}
      var i, j, k : integer;
          one : real;
          sum : real;
      begin
         one := 1.0;
	 for i := 1 to n do
            sBcov[i] := A[i];
         for i := 1 to n do begin
            sBcov[i][i] := one/A[i][i];
            for j := (i+1) to n do begin
               sum := 0.0;
               for k := i to (j-1) do
                  sum := sum + A[j][k]*sBcov[k][i];
               sBcov[j][i] := -sum/A[j][j];
               end;
            end;
         for i := 1 to n do begin
            for j := i to n do begin
               sum := 0.0;
               for k := j to n do
                  sum := sum + sBcov[k][i] * sBcov[k][j];
               A[i][j] := sum;
               A[j][i] := sum;
               end;
            end;
         end;
      begin {sMatrixInversion}
	 ok := true;
         sRootit(A, Ldet, ok);
         sInvertit(A);
      end;

   procedure FormPy (var Q: Covmatrix; var Py:Vector; upper: integer);
          {Py is VInv * DiffVec.}
      var I, J: integer;
          tally: real;
      begin
         Vec := DiffVec;
         for I := 1 to upper do begin
            tally := 0;
            for J := 1 to upper do
               tally := Q[I][J]*Vec[J] + tally;
            Py[I] := tally;
            end;
         end;

   procedure Correction (var Rfix: Traitsquare; Des:Desmatrix);
      {Rfix is X' VInv X for a single family.}
      var j, k, l : integer;
      begin
         for j := 1 to slevels do begin
            for l:= 1 to vert[g] do begin
               VInvX[l,j] :=0;
               for k := 1 to vert[g]
               do VInvX[l,j] := Des[k,j]*Q[g][k][l] + VInvX[l,j];
               end; {l}
            end; {j}
         for j := 1 to slevels do begin
            for l := 1 to slevels do begin
               Rfix[j,l] :=0;
               for k := 1 to vert[g]
                 do Rfix[j,l] := VInvX[k,j] * Des[k,l] + Rfix[j,l];
               end; {l}
            end; {j}
      end; {Correction}

   procedure InvertXVX (var XVX: Traitsquare; var XVXInv: Traitsquare;
                        var bfl: boolean; var flagvec: slevelsbool;
                        var XVXdet: real; slevels: integer);

   {This procedure truncates the XVX matrix in case it's degenerate..
    a situation arising from the absence of observed individuals in
    some fixed factor at some level}

   var i, j, s, t: integer;
   var flag: boolean;
   
     begin
      if Iter = 1 
      then begin
      bfl := true;
      for i := 1 to slevels do begin
       flag := false;
       flagvec[i] := true;
       for j := 1 to slevels do
       if (abs(XVX[i,j]) > 1.0E-08 ) then flag := true;
       if flag = false
        then begin
        flagvec[i] := false;
        bfl := false;
        end;
       end; {i}
       end; {iter = 1} 
      if bfl = false 
      then begin
      s := 0;
      for i := 1 to slevels do begin
         t := 0;
         if flagvec[i]
         then begin
            s := s + 1;
            for j := 1 to slevels do begin
               if flagvec[j]
               then begin
                  t := t + 1;
                  XVX[s,t] := XVX[i,j];
                  end;
               end;
            end;
         end;
      end
   else s:= slevels;
   sMatrixInversion(XVX,XVXdet,s,ok);
      s := 0;
      for i := 1 to slevels do begin
         t := 0;
         if flagvec[i]
         then begin
            s := s + 1;
            for j := 1 to slevels do begin
               if flagvec[j]
               then begin
                  t := t + 1;
                  XVXInv[i,j] := XVX[s,t];
                  end
               else XVXInv[i,j] := 0.0;
               end;
            end
         else begin
            for j := 1 to slevels do
               XVXInv[i,j] := 0.0;
            end;
         end;
    end;  {procedure}

   procedure EstimateMeans (var Means: Parammean);
      {The means are estimated as (X'VInvX)Inv* X'VInvY.}
      var i, j : integer;
          tally : real;
      begin
         for i := 1 to slevels do begin
            tally := 0;
            for j:= 1 to slevels do
                  tally := XVXInv[i,j] * XVY[j] + tally;
               Means[i] := tally;
            end;
          end;
   begin {ComputeLogLikelihood}
      Matprod:=0;
      SumS:=0;
      SumLdet :=0;
      for I := 1 to slevels do begin
         for J := 1 to slevels do begin
            XVX[I,J] := 0;
            XVXInv[I,J] := 0;
            end;
         XVY[I] := 0;
         end;
      for g:=1 to famnum do begin
         if Iter > 1 then
         MatrixInversion (Q[g],Ldet, vert[g], ok)
         else
         Ldet := 0.0;
         SumLdet := Ldet + SumLdet;
         Correction (Rfix, Design[g]);
         for I:= 1 to slevels do begin
            for J := 1 to slevels do
               XVX[I,J] := Rfix[I,J] + XVX[I,J];
            end;
         GetWtdVals (WtdVals, VInvX, vert[g], DataVec[g],Design[g]);
         for I := 1 to slevels do
            XVY[I] := WtdVals[I] + XVY[I];
      end; {g}
      InvertXVX (XVX,XVXInv,bfl,flagvec,XVXdet,slevels);
      writeln ('done 2nd inversion');
      SumS := XVXdet;
      EstimateMeans (Means);
      for g := 1 to famnum do begin
         Vechs (MeanVec, Means, Design[g], vert[g]);
         Standardize (DataVec[g], MeanVec, vert[g], DiffVec);
         FormPy (Q[g], Py[g], vert[g]);
         for I:=1 to vert[g] do
	     Matprod:=DiffVec[I]*Py[g][I]+Matprod;
         end; {g}
      if REML
        then LogLike:=-(Matprod+SumLdet+SumS)/2
        else LogLike := -(Matprod + SumLdet)/2;
     writeln ('Matprod  SumLdet  ', Matprod:16:9, SumLdet:16:9);
     writeln (Analysis,'LogLike ', LogLike:15:4);
     writeln ('LogLike ', LogLike:15:4);
   end;{Compute LogLikelihood}

procedure CalcNewEsts (var Parms: Parmvec; var B: info;
                       var DVA: derivsa; var DVD: derivsd; var DVE: derivse;
                       var tdva, tdvd, tdve: countset);
   {New estimates are calculated using Fisher's scoring method. (See K.
   Meyer, 1983. J. Dairy Sci. 66: 1988-1997.  The inverse of the information
   matrix (B) and the first derivatives of the likelihood function with
   respect to each parameter (D) are obtained.  The new vector of estimates
   is Binv * y'P * Py.}
   var length, totvar : integer;
       relnum, fcount, fc, i, j, k, s, t,l, row, col: integer;
       tally : real;

procedure DVprods (var DVA:derivsa; var DVD:derivsd;
			      var DVE:derivse;
			      var tdva, tdvd, tdve: countset; 
			      length:integer; var Psub: Vector;
			      var D:Parmvec; var CU:Desset; var AB:
			      Famphens; var fc: integer);
      type shortarray=array[1..p] of real;
      var g, i, s, j : integer;
	  fact, index: integer;
	  tally :real;
          t:shortarray;
      begin
	 for g := 1 to totvar do
	    D[g] := 0.0;
         for g := 1 to nvc do begin
     if ((not Trunked) or (not Zeroed[g])) then begin
	    for i := 1 to length do begin
	       tally := 0.0;
	       for s := 1 to slevels do
                  t[s] := 0;
               for j := 1 to tdva[fc][g][i] do begin
		  index := DVA[fc][g][j][i] mod 10000;
		  fact := trunc(DVA[fc][g][j][i]/10000);
		  tally := tally + Psub[index]/fact;
		  if REML 
		  then begin
		    for s := 1 to slevels do
		       t[s] := t[s] + AB[index,s]/fact;
                    end;
                 end;
             D[g] := D[g] + Psub[i] * tally;
	     if REML
	     then begin
		for s := 1 to slevels do
   	           CU[g][i][s] := t[s];
                end;
	     end;
     end; {if not Trunked}
     if ((not Trunked) or (not Zeroed[g+nvc])) then begin
	    for i := 1 to length do begin
	       tally := 0.0;
	       for s := 1 to slevels do
                  t[s] := 0;
               for j := 1 to tdve[fc][g][i] do begin
		  index := DVE[fc][g][j][i] ;
		  tally := tally + Psub[index];
		  if REML 
		  then begin
		    for s := 1 to slevels do
		       t[s] := t[s] + AB[index,s];
                    end;
                 end;
             D[g+nvc] := D[g+nvc] + Psub[i] * tally;
	     if REML
	     then begin
		for s := 1 to slevels do
   	           CU[g+nvc][i][s] := t[s];
                end;
	     end;
     end; {if not Trunked}
     if ((not Trunked) or (not Zeroed[g+2*nvc])) then begin
	    for i := 1 to length do begin
	       tally := 0.0;
	       for s := 1 to slevels do
                  t[s] := 0;
               for j := 1 to tdvd[fc][g][i] do begin
		  index := DVD[fc][g][j][i] mod 10000;
		  fact := trunc(DVD[fc][g][j][i]/10000);
		  tally := tally + Psub[index]/fact;
		  if REML 
		  then begin
		    for s := 1 to slevels do
		       t[s] := t[s] + AB[index,s]/fact;
                    end;
                 end;
             D[g+2*nvc] := D[g+2*nvc] + Psub[i] * tally;
	     if REML
	     then begin
		for s := 1 to slevels do
   	           CU[g+2*nvc][i][s] := t[s];
                end;
	     end;
     end; {if not Trunked}
	 end; {g}
      end;{DVprods}

   procedure MakeA (var Des: Desmatrix; var VInv: Covmatrix;
		    var VInvX: Famphens);
      {A is VInvX.}
      var w,t,s : integer;
          tally : real;
      begin
         for w := 1 to length do begin
            for t := 1 to slevels do begin
               tally := 0;
               for s := 1 to length do
                  tally := VInv[w][s]*Des[s,t] + tally;
               VInvX[w,t] := tally;
               end;
            end;
         end;

   
procedure Root (var A: Traitsquare; n: integer);
     {The square root.}
      var i, j, k: integer;
          sum: real;
      begin
         A[1,1] := sqrt(A[1,1]);
         for i := 2 to n do
            A[i,1] := A[i,1]/A[1,1];
       for i := 2 to n do begin
            sum := 0.0;
            for j := 1 to (i-1) do
               sum := sum +sqr(A[i,j]);
            A[i,i] := sqrt(A[i,i]-sum);
            if (i<n) then begin
               for j := (i+1) to n do begin
                  sum := 0.0;
                  for k := 1 to (i-1) do
                     sum := sum + A[i,k] * A[j,k];
                 A[j,i] := (A[j,i]-sum)/ A[i,i];
                 end;
               end;
            end;
         end;

procedure MakeAB (var AB, VInvX: Famphens);
   var w, t, s: integer;
       tally: real;
   begin
      for w := 1 to vert[fcount] do begin
         for t:= 1 to slevels do begin
            tally := 0;
            for s := t to slevels do
               tally := VInvX[w,s]*XVXInv[s,t] + tally;
            AB[w,t] := tally;
            end;
         end;
      end;

   
   procedure ComputeBig (var Q: Covmatrix; var B: info; length, totvar, fc:integer);
      {The information matrix for a given family is generated.  Each element
      is the trace of (VInv dVdI VInv dVdJ).  These are summed over families
      in the main procedure to give the info matrix for the whole dataset.}
      type linemat=array[1..maxnparm,1..pos] of real;
      var g, h, i, j, k, l, m, fact, index : integer;
        tally : real;
	  QGI, QHI: linemat;
      begin
	for i := 1 to length do begin
	   for j := 1 to totvar do begin
	      for k := 1 to length do begin
                 QHI[j,k] := 0.0;
                 QGI[j,k] := 0.0;
		 end;
              end;
	   for j := 1 to nvc do begin
	      for k := 1 to length do begin
		 tally := 0.0;
        if ((not Trunked) or (not Zeroed[j])) then begin
		 for l := 1 to tdva[fc][j][k] do begin
                    index := DVA[fc][j][l][k] mod 10000;
		    fact := trunc(DVA[fc][j][l][k]/10000);
		    tally := tally + Q[i][index]/fact;
		    end; {l}
                 QGI[j,k] := tally;
        end; {if not Trunked}
        if ((not Trunked) or (not Zeroed[j+2*nvc])) then begin
		 tally := 0.0;
		 for l := 1 to tdvd[fc][j][k] do begin
                    index := DVD[fc][j][l][k] mod 10000;
		    fact := trunc(DVD[fc][j][l][k]/10000);
		    tally := tally + Q[i][index]/fact;
		    end; {l}
                 QGI[j+2*nvc,k] := tally;
        end; {if not Trunked}
        if ((not Trunked) or (not Zeroed[j+nvc])) then begin
		 for l := 1 to tdve[fc][j][k] do begin
                    index := DVE[fc][j][l][k];
                    QGI[j+nvc,k] := Q[i][index];
                    end; {l}
        end; {if not Trunked}
        end; {k}
        if ((not Trunked) or (not Zeroed[j])) then begin
		 for m := 1 to tdva[fc][j][i] do begin
		    index := DVA[fc][j][m][i] mod 10000;
		    fact := trunc(DVA[fc][j][m][i]/10000);
		    for k := 1 to length do
		       QHI[j,k] := QHI[j,k] + Q[k][index]/fact;
		    end; {m}
        end; {if not Trunked}
        if ((not Trunked) or (not Zeroed[j+2*nvc])) then begin
		 for m := 1 to tdvd[fc][j][i] do begin
		    index := DVD[fc][j][m][i] mod 10000;
		    fact := trunc(DVD[fc][j][m][i]/10000);
		    for k := 1 to length do
		    QHI[j+2*nvc,k] := QHI[j+2*nvc,k] + Q[k][index]/fact;
		    end; {m}
        end; {if not Trunked}
        if ((not Trunked) or (not Zeroed[j+nvc])) then begin
		 for m := 1 to tdve[fc][j][i] do begin
		    index := DVE[fc][j][m][i];
		    for k := 1 to length do
		       QHI[j+nvc,k] := Q[k][index];
		    end; {m}
        end; {if not Trunked}
		 end;  {j}
	      for g := 1 to totvar do begin
        if ((not Trunked) or (not Zeroed[g])) then begin
		 for h := 1 to totvar do begin
        if ((not Trunked) or (not Zeroed[h])) then begin
		    tally := 0.0;
		    for l := 1 to length do begin
		       tally := tally + QGI[g,l] * QHI[h,l];
		       end;
		    B[g,h] := B[g,h] + tally;
              end; {if no Trunked h}
                    end; {h}
              end; {if not Trunked g}
                 end; {g}
              end; {i}
           end; {computebig}

procedure Moreprods (var UpCU: Traitset; var CU: Desset; var RBig: info;
                     var V: Covmatrix);
   var i, w, j, s, q, t, n : integer;
       tally : real;
   begin
      n := vert[fcount];
      for i := 1 to totvar do begin
        if ((not Trunked) or (not Zeroed[i])) then begin
         for t := 1 to slevels do begin
            for w := 1 to slevels do begin
               tally := 0;
               for j := 1 to n do
                  tally := AB[j,t]*CU[i][j,w] + tally;
               UpCU[i][t,w] := tally;
               end;
            end;
         for s := 1 to n do begin
            for q := 1 to slevels do begin
               tally := 0;
               for t := 1 to n do
                  tally := V[s][t]*CU[i][t,q]+tally;
               VInvX[s,q] := tally;
               end;
            end;
        for j := i to totvar do begin
        if ((not Trunked) or (not Zeroed[j])) then begin
           tally := 0;
           for s := 1 to slevels do begin
              for q := 1 to n do
                 tally := CU[j][q,s]*VInvX[q,s] + tally;
              end;
           RBig[i,j] := -2*tally;
           RBig[j,i] := -2*tally;
           {This RBig is the second term in the expression for the info matrix
            for REML: see Meyer, K. 1983. eq. 9, substituting in the expression
            for P and multiplying out gives 4 terms.  The first is just the
            same as the info matrix for ML, the middle two can be combined and
            rearranged to give -2*tr(u'dV/dTj*VInv*dV/DTi*u), where
            uu'=VInvX(X'*VInv*X)X'VInv.  This term is the quantity caluclated
            this procedure.  See ComputeRBig for the last term.}
        end; {if not Trunked j}
           end; {j}
        end; {if not Trunked i}
        end; {i}
      end; {Moreprods}

procedure ComputeRBig (var RBig: info; var SumUpCU:Traitset);
   {Compute the last term in the info matrix of REML.  It is
   tr(u'dV/dTiu)(u'dV/dTju), where u is defined in Dvprods.
   see Meyer in J. Dairy Sci).}
   var tally : real;
       i, j, p, q : integer;
   begin
      for i := 1 to totvar do begin
        if ((not Trunked) or (not Zeroed[i])) then begin
         for j := i to totvar do begin
        if ((not Trunked) or (not Zeroed[j])) then begin
            tally := 0;
            for q := 1 to slevels do begin
               for p := 1 to slevels do
                  tally := SumUpCU[i][q,p]*SumUpCU[j][p,q] + tally;
               end;
            RBig[i,j] := tally;
            RBig[j,i] := tally;
        end; {if not Trunked j}
            end; {j}
        end; {if not Trunked i}
         end; {i}
     end; {ComputeRBig}

procedure MatrixInversion (var A:info ; dim: integer);
var i, j, k, lower, upper : integer;
      temp: real;
begin {Gauss-Jordan reduction -- invert
       matrix in place, overwriting previous
       contents of A.}
   lower := 1;
    upper := dim;
    for i := lower to upper do begin
      if A[i,i] = 0
      then writeln ('the info matrix is singular');
      temp := 1.0 / A[i, i];
      A[i, i] := 1.0;
      for j := lower to upper do
         A[i, j] := A[i, j] * temp;
      for j := lower to upper do
         if j <> i
         then begin
            temp := A[j, i];
            A[j,i] := 0.0;
            for k := lower to upper do
               A[j, k] := A[j, k] - temp * A[i, k];
            end;
      end;
   end; {invert}

   procedure Newtheta (D: Parmvec; var Theta: Parmvec; upper: integer;
         var Trunked2: boolean; var Zeroed2, Zeroed3: Parmbool);
      {takes the product of BInv and D to give new parameter estimates}
      var i, j, k, l, row, col, Mrow, Ml: integer;
      var Lam, Mlam: real;
      begin
         for i := 1 to upper do begin
            Theta[i] := 0;
            for j := 1 to upper do begin
               Theta[i] := B[i,j] * D[j] + Theta[i];
               end;
             end; {i}
          if feasable then begin
           l := 0;
           Mlam := 0.0;
           for i := 1 to vartypes do begin
             for j := 1 to r do begin
               for k := j to r do begin
                l := l+1;
                row := (i-1)*r + j;
                if k = j
                  then begin
               if Theta[l] < 0.0 
               then begin
               Lam := -Theta[l]/(OldParms[l]-Theta[l]);
               if Lam = 1 then begin
               Trunked2 := true;
               Zeroed2[l] := true;
               Zeroed3[row] := true;
               writeln ('param  ',l,Theta[l]:8:3,' constrained back to 0');
                Theta[l] := 0.0;
              end {if Lambda = 1..component went from 0 to negative value}
               else begin 
               if Mlam < Lam
                then begin
                Ml := l;
                Mrow := row;
                Mlam := Lam;
               end;{if Mlam}
               end; {else}
             end; {if theta}
            end; {if k=j}
           end; {k}
          end; {j}
         end; {i}
          if Mlam > 0.0
          then begin
               Trunked2 := true;
               writeln (' Mlam ',Mlam);
               Zeroed2[Ml] := true;
               Zeroed3[Mrow] := true;
               writeln ('param  ',Ml,Theta[Ml]:8:3,' constrained');
           if Mlam < 1.0
           then begin
         for i := 1 to upper do
          Theta[i] := Mlam * OldParms[i] + (1.0-Mlam) * Theta[i];
            Theta[Ml] := 0.0;
            end;
            end;
           l := 0;
            {now constrain covariances if necessary}
           for i := 1 to vartypes do begin
             for j := 1 to r do begin
               for k := j to r do begin
                l := l + 1;
                if k > j 
                then begin
                row := (i-1) * r + j;
                col := (i-1) * r + k;
                 if (not Zeroed[l] and (Zeroed3[row] or Zeroed3[col])) 
               then begin
               Zeroed2[l] := true;
               writeln ('param  ',l,Theta[l]:8:3,' constrained');
               Theta[l] := 0.0;
              end; {if Zeroed3}
              end; {if k>j}
             end; {k}
            end; {j}
           end; {i}
        end; {if feas}
          end;

 begin {CalcNewEsts}
    B := Empty;
    RBig := Empty;
    totvar := vartypes * nvc;
    for i := 1 to totvar do begin
       SumD[i] := 0;
       for j := 1 to slevels do begin
          for k := 1 to slevels do
             SumUpCU[i][j,k] := 0;
          end;
       end;
      XVX := XVXInv;
      if bfl = false 
      then begin
      s := 0;
      for i := 1 to slevels do begin
         t := 0;
         if flagvec[i]
         then begin
            s := s + 1;
            for j := 1 to slevels do begin
               if flagvec[j]
               then begin
                  t := t + 1;
                  XVX[s,t] := XVX[i,j];
                  end;
               end;
            end;
         end;
      end
   else s:= slevels;
    Root (XVX, s); 
      s := 0;
      for i := 1 to slevels do begin
         t := 0;
         if flagvec[i]
         then begin
            s := s + 1;
            for j := 1 to slevels do begin
               if flagvec[j]
               then begin
                  t := t + 1;
                  XVXInv[i,j] := XVX[s,t];
                  end
               else XVXInv[i,j] := 0.0;
               end;
            end
         else begin
            for j := 1 to slevels do
               XVXInv[i,j] := 0.0;
            end;
         end;
    for fcount := 1 to famnum do begin
       fc := fcount;
       length := vert[fcount];
       relnum := famsize[fcount];
       if REML
       then begin
  	  MakeA (Design[fcount],Q[fcount],VInvX);
          MakeAB (AB,VInvX);
          end;
       DVprods (DVA,DVD,DVE, tdva,tdvd,tdve,length,Py[fcount],
                     D,CU,AB,fc);
       Big := Empty;
       ComputeBig (Q[fcount],Big, length, totvar, fc);
       for i := 1 to totvar do begin
        if ((not Trunked) or (not Zeroed[i])) then begin
          SumD[i] := D[i] + SumD[i];
          for j := 1 to totvar do begin
        if ((not Trunked) or (not Zeroed[j])) then begin
             B[i,j] := Big[i,j] + B[i,j];
             end; {if not Trunked j}
            end;{j}
            end; {if not Trunked i}
          end; {i}
       if REML
       then begin
          Moreprods (UpCU, CU, RBig, Q[fcount]);
          for i := 1 to totvar do begin
        if ((not Trunked) or (not Zeroed[i])) then begin
             for j := 1 to slevels do begin
                for k := 1 to slevels do
                   SumUpCU[i][j,k] := SumUpCU[i][j,k] + UpCU[i][j,k];
                end;
             for j := 1 to totvar do begin
        if ((not Trunked) or (not Zeroed[j])) then begin
                B[i,j] := B[i,j] + RBig[i,j];
             end; {if not Trunked j}
             end; {j}
             end; {j}
          end; {if not Trunked i}
          end; {i}
       end; {fcount}
      if REML
      then begin
         ComputeRBig (RBig, SumUpCU);
	 writeln ('done computeRbig');
         for i := 1 to totvar do begin
        if ((not Trunked) or (not Zeroed[i])) then begin
             for j := 1 to totvar do begin
        if ((not Trunked) or (not Zeroed[j])) then begin
                B[i,j] := B[i,j] + RBig[i,j];
             end; {if j}
             end; {j}
             end; {j}
         end; {if not Trunked i}
         end;
 {build from B and D the vector of
  first derivitives eD.}
 writeln (' For the current iteration, the derivatives are:  ');
  for i := 1 to totvar do begin
     tally := 0.0;
    for j := 1 to totvar do begin
    eB [i,j] := -0.5 * B [i,j];
     tally := tally + eB[i,j] * Parms [j];
     end;
    eD [i] := 0.5 * SumD [i] + tally;
      writeln (i,eD[i]:8:4);
  end;
  if feasable
 then begin
  l := 0;  {release zeroed terms with positive derivatives }
  for i := 1 to vartypes do begin 
   for j := 1 to r do begin
     for k := j to r do begin
     l := l + 1;
     if j = k 
     then begin
     row := (i - 1)*r + j;
      if eD[l] > 0.0
       then begin
      if (Trunked2 and Zeroed2[l]) 
       then begin
       Zeroed2[l] := false;
       Zeroed3[row] := false;
       writeln ('Zeroed ',l,'released');
       end; {if Trunked2}
      end;   
      end;   
   end; {k}
   end; {j}
  end; {i}
  l := 0;  {release zeroed covariances if possible }
  for i := 1 to vartypes do begin 
   for j := 1 to r do begin
     for k := j to r do begin
     l := l + 1;
     if (Zeroed2[l] and (not Zeroed[l]))
      then begin
     if j < k 
     then begin
     row := (i - 1)*r + j;
     col := (i - 1)*r + k;
     if ((not Zeroed3[row]) and (not Zeroed3[col]))
      then begin
       writeln ('Zeroed ',l,'released');
       Zeroed2[l] := false;
      end;  {if}
    end; {j<k} 
      end;  {if Zeroed2} 
   end; {k}
   end; {j}
  end; {i}
 end; {if feas}
 if Trunked or Trunked2
   then begin
      s := 0;
      for i := 1 to totvar do begin
         t := 0;
         if not Zeroed2[i]
         then begin
            s := s + 1;
            for j := 1 to totvar do begin
               if not Zeroed2[j]
               then begin
                  t := t + 1;
                  B[s,t] := B[i,j];
                  end;
               end;
            end;
         end;
      end
   else s:= totvar;
   MatrixInversion (B, s);
   if Trunked or Trunked2
   then begin
      C := B;
      s := 0;
      for i := 1 to totvar do begin
         t := 0;
         if not Zeroed2[i]
         then begin
            s := s + 1;
            for j := 1 to totvar do begin
               if not Zeroed2[j]
               then begin
                  t := t + 1;
                  B[i,j] := C[s,t];
                  end
               else B[i,j] := 0.0;
               end;
            end
         else begin
            for j := 1 to totvar do
               B[i,j] := 0.0;
            end;
         end;
      end;
    writeln ('at newtheta');
    Newtheta (SumD, Parms, totvar,Trunked2,Zeroed2,Zeroed3);
   end;

procedure Test ( Old: real; Current: real; var Maxed, better: boolean);
   {Compares the absolute value of the difference between the current
   and last loglikelihood to delta, the stopping criterion, to determine
   whether a maximum has been reached.}
   var diff: real;
   begin
      diff := Old-Current;
      if (diff<0)
      then better := true;
      diff := abs(diff);
      if (diff - delta < 0) then Maxed := true
       else Maxed := false;
      end;

procedure Trunck (var Parms: Parmvec);
   var i: integer;
       totvar: integer;
   begin
      totvar := vartypes*((nchar*(nchar+1)) div 2);
      for i := 1 to totvar do
        if Zeroed[i]
        then Parms[i] := 0.0;
      end;

procedure Report;
   var i, j, k, l: integer;
       totvar: integer;
       count: integer;

   begin
    totvar := vartypes * ((nchar* (nchar + 1)) div 2);
    writeln (Analysis,'    At iteration ',Iter:4);
    writeln (Analysis);
    if Trunked
    then write (Analysis,'   *** The constrained analysis converged ')
    else write (Analysis,'   *** The unconstrained analysis converged ');
    writeln (Analysis,'with the following results ***');
    writeln (Analysis,'   The log likelihood is ', LogLike:8:4);
    writeln (Analysis,'   The mean of each trait is');
    for i := 1 to nchar do
       writeln (Analysis, '   ',Means[i]:14:6);
    writeln (Analysis,'   The effect of the fixed factors is');
    writeln (Analysis,'   (in the order given, levels within factors)');
     count := nchar;
      for i := 1 to nfxfact do begin
       for j := 1 to nfxlevels[i] do begin
        for k := 1 to nchar do begin
     if (j=1){modified by rgs 7/16/92}
     then begin
     writeln (Analysis,'Factor: ',i:4,' Level: ',j:4,' Label: ',levels[i,j]:4);
         writeln (Analysis,'has fixed factor effect  0');
         end
     else begin
        count := count + 1;
        writeln (Analysis,'Factor: ',i:4,' Level: ',j:4);
        writeln (Analysis,' Label: ',levels[i,j]:4);
         writeln (Analysis,'has fixed factor effect  ',Means[count]:14:6);
         end;
       end; {k}
     end; {j}
   end; {i}
    for i :=  1 to slevels do
        writeln (Analysis, '   ',Means[i]:10:5);
    writeln (Analysis);
    write (Analysis,'   The estimates of the');
    writeln (Analysis,' components are:');
    writeln (Analysis, '    Additive');
    writeln (Analysis,'              ');
    l := 1;
    for i := 1 to nchar do begin
       for j := i to nchar do begin
	  write (Analysis, BParms[l]:14:6);
          l := l+1;
          end;
       writeln(Analysis);
       write (Analysis,'              ');
       end;
    writeln(Analysis);
    writeln (Analysis, '   Environmental');
    writeln (Analysis,'              ');
    for i := 1 to nchar do begin
       for j := i to nchar do begin
          write (Analysis, BParms[l]:14:6);
          l := l+1;
          end;
       writeln(Analysis);
       write (Analysis,'              ');
       end;
    writeln(Analysis);
    if ceo = 1 then
    writeln (Analysis, '   Dominance');
    if ceo = 4 then
    writeln (Analysis, '   Common Environment');
    writeln (Analysis,'              ');
    for i := 1 to nchar do begin
       for j := i to nchar do begin
          write (Analysis, BParms[l]:14:6);
          l := l+1;
          end;
       writeln(Analysis);
       write (Analysis,'              ');
       end;
    writeln (Analysis,'    large-sample var-cov matrix of the estimates');
    for i := 1 to totvar do begin
        for j := i to totvar do
           write (Analysis,B[i,j]:14:6);
        writeln (Analysis);
        for j := 1 to i do
	   write (Analysis,'              ');
        end;
       writeln (Analysis);
       write (Analysis,'  The test statistic comparing two likelihoods is ');
       writeln (Analysis,'given by twice');
       write (Analysis, '  their difference and is compared to ');
       writeln (Analysis,'Chi-square with df given by the ');
       write (Analysis, '  number of parameters ');
       writeln (Analysis, 'specified by the hypothesis.');
    end;

procedure Cleanup (var hypo: Parmvec; var charnum: integer; 
                   var REML: boolean;
                   var ntrnc, nvc: integer; var tvec: Parmbool; 
		   var nfxfact: integer; var feasable: boolean);
   {reads and writes the parameters governing the simulation run}
   var Analtype, feas: integer;
      totvar, i, j, k, index : integer;
   begin
      for i := 1 to maxnparm do
	 tvec[i] := false;
      read (sibships, Analtype, charnum, nfxfact, feas);
      writeln (Analtype);
      writeln ('Read charnum. It is ', charnum:4);
      nvc := (charnum * (charnum+1)) div 2;
      totvar := 3 * nvc;
      writeln ('totvar: the total number of (co)variances to be estimated, is ', totvar:4);
      index := 0; {starting values for easy ride through first iteration}
      for k := 1 to 3 do begin
       for i := 1 to nchar do begin
        for j := i to nchar do begin
        index := index + 1;
        if (k = 2) and (i = j) then
         hypo[index] := 1.0
         else
         hypo[index] := 0.0;
       end;
      end;
     end;
      readln (sibships);
      writeln (Analysis,'  The suggested values for the components are ');
      for i := 1 to totvar do
         write (Analysis,'   ',hypo[i]:10:3);
      writeln (Analysis);
      if Analtype = 2
      then REML:=true
      else REML := false;
      if feas = 1
      then begin
      feasable := true;
      writeln ('Feasability constraints will be applied');
      end
      else feasable := false;
      if REML
      then writeln (Analysis,'  This is an REML analysis.')
      else writeln (Analysis,'  This is an ML analysis.');
      read (sibships, ntrnc);
      writeln ('The number of parameters to be constrained is ', ntrnc:4);
      for i := 1 to ntrnc do begin
         read (sibships, index);
         tvec[index] := true;
         writeln (' index of the component to be constrained to 0 ', index:6);
         end;
      readln (sibships);
      end;

   procedure SetNull (var Parms: Parmvec);
      begin
	 Parms := First;
         end;

procedure TryAgain (Old:Parmvec; var New: Parmvec);
   var i, top: integer;
       Temp: Parmvec;
   begin
      top := vartypes*((nchar*(nchar+1)) div 2);
      writeln ('In Tryagain');
      writeln (Analysis,'In Tryagain');
      for i := 1 to top do begin
	 write (Old[i]:14:4, New[i]:14:4);
         Temp[i] := (New[i] - Old[i])*(0.5);
	 New[i] := Old[i] + Temp[i];
	 writeln (New[i]:14:4);
         end;
      end;

begin {program MaximumLikelihood}
   reset(sibships,'sibships');
   rewrite (Analysis,'Analysis');
   Iter := 0;
   Done := false;
   Trunked := false;
   Trunked2 := false;
   Maxed := false;
   Cleanup (First, nchar, REML, ntrunc, nvc, Zeroed, nfxfact, feasable);
   if ntrunc>0 then Trunked := true;
   Zeroed2:=Zeroed;
   for I := 1 to maxnparm do
   Zeroed3[I] := false;
   Getphtypes (PresentData, M, nchar, Zmatrix, slevels, nfxlevels);
   FindGroupsN (PresentData, M);
   famnum:= 0;
   for I := 1 to M do  begin
      if PresentData[I].Famid=0
      then begin
	 Check (I, PresentData, newfam);
	 if (newfam > 0)
	 then Visit (PresentData[I], newfam)
	 else begin
	    famnum:= famnum+1;
	    Visit (PresentData[I], famnum);
	    end;
	 end;
      end;
   writeln (famnum:5,' families have been found.');
   if famnum > f then halt;
   PrintFams(PresentData, Zmatrix);
   Collecthem(PresentData,Families, FamZ, famsize);
   Designit (FamZ, vert, Design, DataVec);
   SetNull (Parms);
   OldParms := Parms;
   if Trunked
   then Trunck(Parms);
   ZeroEmpty (Empty);
   TA := false;
   repeat
      repeat
         EnterValues (AddCov, DomCov, EnvCov);
         if Iter = 0 then
         MakeCovArray (Q,DVA,DVD,DVE,
                 tdva,tdvd,tdve)
         else
         reMakeCovArray (Q,Parms,DVA,DVD,DVE,
                 tdva,tdvd,tdve);
         Iter := Iter + 1;
         vartypes := 3;
         better:=false;
	 ok := true;
         CheckValues (AddCov, ok);
         if ok
         then begin
	 writeln ('starting comploglike');
         ComputeLogLikelihood (LogLike, Q, Means, XVXInv, Py, vert, DataVec,
	                       Design, ok);
	 writeln ('done comploglike');
	    if Iter>1
	    then Test(OldLogLike, LogLike, Maxed, better)
            else better := true;
               if Iter=maxround
               then better:= true;
	       if not better and not Maxed
	       then begin
               TryAgain (OldParms, Parms);
               TA := true;
               end;
	    end
         else begin
               TryAgain (OldParms, Parms);
               TA := true; 
               end;
      until better or Maxed;
      if Maxed and not TA
      then begin
         BParms := Parms;
         Done := true;
         Report;
         end;
      OldLogLike := LogLike;
      OldParms := Parms;
      if not Done
      then begin
      CalcNewEsts (Parms,B,DVA,DVD,DVE, tdva,tdvd,tdve);
      TA := false;
      end;
      if Iter=maxround
      then begin
         writeln ('Exceeded max number of iterations without converging.');
         Done := true;
         end;
    until Done;
    BParms := Parms;
    {close(Analysis);}
    end {program MaximumLikelihood}.