program  EG21;
(* GEE2 for binary data. Correction to GEE1 variance estimator. *)
(*****************************)
%include QMATRIX; (* quick matrix routines *)
%include CMS;     (* filedef command *)
(****************************)
const
  VERSION1 = 'Extended Generalized Estimating Equations (v0.5)';
  VERSION2 = '       for Correlated Binary Data';
  VERSION3 = '        (c) Bahjat Qaqish, 1989, 1990, 1991';
  LRECL   = 133;   (* max record length on control file *)
  MAXINPX = 64;    (* dim(x) in input data file *)
  MAXDIMB = 20;    (* dim(beta) *)
  MAXROWSX=16000;   (* max size of expanded data set *)
  MAXCOLSX= MAXDIMB;  (* max # colums of the x matrix *)
  MAXNTYP = 12;    (* max # of primary types *)
  MAXXTYP = 90;    (* max # of primary+extended types *)
                   (* = 2*MAXNTYP + ch2(.MAXNTYP.)
                     A linear predictor spec exists for
                     each primary and extended type. *)
  HALTYE  = TRUE;
  HALTNO  = FALSE;
  MAXSIZ  = 12;  (* max cluster size *)
  MAXSIZ2 = 78;  (* MAXSIZ + ch2(.MAXSIZ.) *)

  MAXITERLO = 0;
  MAXITERHI = 100;  (* MAXITERLO <= maxiter <= MAXITERHI *)

  VMZP = 1; (* Zhao and Prentice formulae *)
  VMBQ = 2; (* exact BQ *)


type
  xmat_  = array(.1..MAXROWSX,1..MAXCOLSX.) of shortreal;
  xx_    = array(.1..MAXINPX.) of real; (* input x data *)
  ndxv_  = array(.1..MAXDIMB.) of integer;  (* index vector *)
  lrvec_  = array(.1..MAXROWSX.) of real;  (* long real vector *)
  livec_  = array(.1..MAXROWSX.) of integer; (* long integer vector *)
  sivec_  = array(.1..MAXSIZ.) of integer; (* short integer vector *)
  typtbl_ = array(.0..MAXNTYP,0..MAXNTYP.) of integer;
  reg1_   =  record    (* one regression equation *)
      t1, t2,          (* primary types t1 < t2 *)
      dim   : integer; (* dimension , 0 means lp := 0 *)
      xndx,            (* x-index for one type *)
      bndx  : ndxv_;    (* beta-index for onr type *)
      defined : BOOLEAN; (* has been defined? *)
    end;
  reg_      = array(.1..MAXXTYP.) of reg1_; (* all regressions *)
  solvpb_ = record   (* parameter block for solv3 and solv4 *)
      u1, u2, u3, u4,
      u12, u13, u14, u23, u24, u34,
      u123, u124, u134, u234 : real;
    end;
  string_ = string(133);
  used_ = array(.1..MAXDIMB.) of BOOLEAN;
  lbl_  = array(.1..MAXDIMB.) of string(25);
  uijka_  = array(.1..MAXSIZ,1..MAXSIZ,1..MAXSIZ.) of real;
  uijkla_ = array(.1..MAXSIZ,1..MAXSIZ,1..MAXSIZ,
                   1..MAXSIZ.) of real;
    (* The storage for the last 2 arrays is to be optimized later *)


var
  y    : lrvec_;   (* data 0/1 *)
  u    : lrvec_;   (* mean *)
  lp   : lrvec_;   (* linear predictor *)
  psi  : lrvec_;   (* odds ratios *)
  typ  : livec_;   (* type of this row *)
  id   : livec_;   (* cluster id *)
  siz  : livec_;   (* cluster size *)
  x    : xmat_;    (* regressors *)
  b    : vector;     (* parameters *)
  lbl  : lbl_;     (* labels for betas *)
  rvar : matrix;   (* robust variance matrix *)
  nvar : matrix;   (* naive  variance matrix *)
  rse  : vector;     (* robust s.e. *)
  nse  : vector;     (* naive  s.e. *)
  rz   : vector;     (* robust z-statistics *)
  s1,
  s1a,             (* the asymmetric version os s1 *)
  s3   : matrix;   (* pxp used in the gee algorithm *)
  s2   : vector;     (* px1 used in the gee algorithm *)
  used : used_;    (* beta used ? indicators *)
  reg  : reg_;     (* all regressions described here *)
  strt : livec_;   (* where cluster primary data starts *)
  fnsh : livec_;   (* where cluster primary data ends   *)
  ch2  : sivec_;   (* n choose 2 *)
  sizf : sivec_;   (* freqs on cluster size *)
  typtbl:typtbl_;  (* type translation *)
  dimx : integer;  (* dim of x in input data file *)
  dimb : integer;  (* dim of beta in the model *)
  dimb1: integer;  (* dimb1 = dim(b) for main effects (GEE1) *)
  geelvl: integer; (* 1=GEE1, 2=GEE2 *)
  nclus: integer;  (* # clusters *)
  ntyp : integer;  (* # primary types *)
  xtyp : integer;  (* # primary+extended types *)
  expn : integer;  (* size of the expanded data set *)
  maxiter:integer; (* convergence criterion *)
  iter :integer;   (* iteration count *)
  tol  : real;     (* convergence criterion *)
  converged : BOOLEAN;
  date, time: alfa;
  CPUtime   : real;
  vmethod   : integer; (* method for 3rd and 4th order moments *)
  monbflag  : BOOLEAN; (* monitor beta each iteration ? *)
  vshape    :integer; (* use a block diagonal v ? *)
  indebug   : BOOLEAN;
  debugmode : BOOLEAN; (* control of the debug procedure *)
  solvflag : BOOLEAN; (* control of the solv4 procedure *)
  probescore:BOOLEAN; (* probe the scores?  *)
  debugclus : integer; (* cluster for debugging *)
  datafile  : text; (* data *)
  probef    : text; (* data *)
  ctrlfile  : text; (* data *)
  uijka     : uijka_;  (* array for means uijk = E(Yijk) *)
  uijkla    : uijkla_; (* array for means uijkl = E(Yijkl) *)
  bias      : integer; (* stage in bias computation *)


(************************)
function sign(f : real) : integer;
Begin
  if f > 0 then sign := 1 else if f < 0 then sign := -1
  else sign := 0;
End;
(************************)
procedure RepCPUTime;
Begin (* RepCPUTime *)

  DateTime (date, time);
  writeln('Date: ', date, '    Time: ', time);

  CPUtime := (Clock-CPUtime)*1.0E-6;
  write  ('CPU time: ', CPUtime:10:3, ' seconds = ');
  writeln(CPUtime/60:10:3, ' minutes.');

End;  (* RepCPUTime *)
(************************)
procedure err (s : string_; haltflag : BOOLEAN);
  (* central error control: all errors routed here *)
Begin
  writeln(s);
  if haltflag then HALT;
End;
(************************)
procedure incr(var i,j : integer; n : integer);
  (* increment a product term index pair for cluster size n *)
  (* assumes 1 <= i < j <= n *)
Begin (* incr *)
  if (j = n) then
  begin
    i := succ(i);
    j := succ(i);
  end else j := succ(j);
End;  (* incr *)
(************************)
procedure swap(var i, j : integer);
var temp : integer;
Begin    (* swap *)
  temp := i; i := j; j := temp;
End; (* swap *)
(************************)
function ui (i, c : integer) : real;
 (* returns the mean for unit i in cluster c *)
Begin (* ui *)
  ui := u(.pred(strt(.c.)) + i.);
End;  (* ui *)
(************************)
function uij (i, j, c : integer) : real;
 (* returns the mean for the i,j product in cluster c *)
 (* assumes: i <> j *)
Begin (* uij *)
  if (i > j) then swap (i, j);
  uij := u(. fnsh(.c.) + (j-i) + (pred(i)*(2*siz(.c.)-i)) DIV 2 .);
End;  (* uij *)
(************************)
function uijk (i, j, k : integer) : real;
  (* returns E(Yijk) for current cluster *)
  (* assumes uijka(.i,j,k.) has been filled properly *)
Begin (* uijk *)
  (* first sort the indices so that i<j<k *)
  if (i > j) then swap(i, j);
  if (j > k) then swap(j, k);  (* k *)
  if (i > j) then swap(i, j);  (* i,j *)
  uijk := uijka(.i,j,k.);
End; (* uijk *)
(************************)
function uijkl (i, j, k, l : integer) : real;
  (* returns E(Yijkl) for current cluster *)
  (* assumes uijka(.i,j,k,l.) has been filled properly *)
Begin (* uijkl *)
  (* first sort the indices so that i<j<k<l *)
  if (i > j) then swap(i, j);
  if (j > k) then swap(j, k);
  if (k > l) then swap(k, l);  (* l *)
  if (i > j) then swap(i, j);
  if (j > k) then swap(j, k);  (* k *)
  if (i > j) then swap(i, j);  (* i, j *)
  uijkl := uijkla(.i,j,k,l.);
End; (* uijkl *)
(************************)
function covyiyjk(i, j, k, c : integer) : real;
  (* returns the COV (Yi,Yjk) for cluster c *)
  (* assumes: j < k *)
Begin (* covyiyjk *)
  if (i = j) or (i = k) then covyiyjk :=
     uij(j, k, c) * (1-ui(i, c)) else
  if (vmethod = VMZP) then
  (* Z & P formula *)
  covyiyjk := (uij(i,j,c)-ui(i,c)*ui(j,c)) * ui(k,c) +
              (uij(i,k,c)-ui(i,c)*ui(k,c)) * ui(j,c)   else
  covyiyjk := uijk (i,j,k) - ui(i,c)*uij(j,k,c);
End;  (* covyiyjk *)
(************************)
function covyijyik (i, j, k, c : integer) : real;
  (* returns the COV (Yij,Yik) for cluster c *)
  (* assumes i, j, and k are all distinct *)
Begin (* covyiyik *)
  if (vmethod = VMZP) then
  (* Z & P approx. *)
  covyijyik := (uij(i, k, c) - ui(i, c)*ui(k, c)) * ui(j, c) +
               (uij(i, j, c) - ui(i, c)*ui(j, c)) * ui(k, c) +
                uij(j, k, c) * ui(i, c) -
                uij(i, j, c) * uij(i, k, c)    else
  (* exact *)
  covyijyik :=  uijk(i, j, k) - uij(i, j, c)*uij(i, k, c);

End;  (* covyiyik *)
(************************)
function covyijykl(i, j, k, l, c : integer) : real;
(*
  : returns the COV (Yij,Ykl) for cluster c
  : should be called only for the off-diagonals of b22
  : assumes (i<j) and (i<=k<l)and not ((i=k) and (j=l))
  : 7 cases exist
    i < j < k < l         case 1
    i < j = k < l         case 2
    i < k < j < l         case 3
    i < k < j = l         case 4
    i < k < l < j         case 5
    i = k < j < l         case 6
    i = k < l < j         case 7
*)
Begin (* covyijykl *)

  if (j = k) then               (* case 2 *)
    covyijykl := covyijyik (j, i, l, c)  else

  if (j = l) then               (* case 4 *)
    covyijykl := covyijyik (l, i, k, c) else

  if (i = k) then               (* cases 6 and 7 *)
    covyijykl := covyijyik (i, j, l, c) else

  if (vmethod = VMZP) then
  (* Z & P formula:                cases 1, 3 and 5 *)
  covyijykl := (uij(i, k,c) - ui(i,c)*ui(k,c)) * ui(j,c) * ui(l,c) +
               (uij(i, l,c) - ui(i,c)*ui(l,c)) * ui(j,c) * ui(k,c) +
               (uij(j, k,c) - ui(j,c)*ui(k,c)) * ui(i,c) * ui(l,c) +
               (uij(j, l,c) - ui(j,c)*ui(l,c)) * ui(i,c) * ui(k,c)
  else   (* exact *)
  covyijykl := uijkl(i, j, k, l) - uij(i, j, c)*uij(k,l,c);

End;  (* covyijykl *)
(************************)
procedure calclp; (* linear predictor *)

var
  i, j  : integer;
  result : real;

Begin (* calclp *)
  with b do
  for i := 1 to  expn do
  begin
    result := 0;
    for j := 1 to dimb do
      if (x(.i,j.)  <> 0) then
        result := result + x(.i,j.) * ve(.j.);
    lp(.i.) := result;
  end;
End;  (* calclp *)
(************************)
function duijdui (ui, uj, psi : real) : real;
  (* returns derivative of uij with respect to ui *)
  (* d uij / d ui *)
  (* evaluated at (ui,uj,psi) *)
var  a, b, k, d : real;
Begin (* duijdui *)
  if (psi=1) then duijdui := uj else
  begin
    a := ui + uj;
    b := psi - 1;
    k := 1 + a*b;
    d := sqrt(sqr(k)-4*psi*b*ui*uj);
    duijdui := 0.5*(1-(k-2*psi*uj)/d);
  end;
End;  (* duijdui *)
(************************)
function duijdpsi (ui, uj, psi : real) : real;
  (* returns derivative of uij with respect to psi *)
  (* d uij / d psi *)
  (* evaluated at (ui,uj,psi) *)
var  a, b, k, d, e, f : real;
Begin (* duijdpsi *)
  if (psi=1) then duijdpsi := 0 else
  begin
    a := ui + uj;
    b := psi - 1;
    k := 1 + a*b;
    f :=  (sqrt(sqr(k)-4*psi*b*ui*uj));
    d := 1.0 / f;
    e := (2*k*a-4*ui*uj*(psi+b))* 0.5 * d;
    e := (a-e) * b;
    e := e + e;
    f := k - f;
    e := e - f - f;
    duijdpsi := e * 0.25 / sqr(b);
  end;
End;  (* duijdpsi *)
(************************)
function uijpsi (ui, uj, psi : real) : real;
  (* returns E(YiYj) from E(Yi), E(Yj), O.R.(i,j) *)
var  a, b, k, d : real;
Begin (* uijpsi *)
  if (psi=1) then uijpsi := ui*uj   else
  begin
    a := ui + uj;
    b := psi - 1;
    k := 1 + a*b;
    d := k-sqrt(sqr(k)-4*psi*b*ui*uj);
    uijpsi := 0.5 * d / b;
  end;
End;  (* uijpsi *)
(************************)
procedure calcpsi; (* calculate the odds ratios *)
const e = 1.0E-6;
var
  i : integer;
Begin (* calcpsi *)
  for i := 1 to expn  do
  if (lp(.i.) > -e) and (lp(.i.) < e) then psi(.i.) := 1
    else psi(.i.) := exp(lp(.i.));
End; (* calcpsi *)
(************************)
procedure calcu ; (* mean *)

var
  clus, i, j, k : integer;
  uf  : file of real;   (* mean *)

Begin (* calcu *)
  for clus := 1 to nclus do
  begin
    (* means of original data *)
    for i := strt(.clus.) to fnsh(.clus.) do
    begin
        u(.i.) := psi(.i.)/(1+psi(.i.));
    end;

    (* means of products *)
    i := fnsh(.clus.);
    for j := strt(.clus.) to pred(fnsh(.clus.)) do
    for k := succ(j) to fnsh(.clus.) do
    begin
      i := succ(i);
      u(.i.) := uijpsi (u(.j.), u(.k.), psi(.i.));
    end;
  end;


(* bias code *)
  if (bias = 1) then
  begin
    rewrite (uf);
    for i := 1 to expn  do
    write (uf,u(.i.));
    close (uf);
    HALT;
  end;
(* bias code *)

End; (* calcu *)
(************************)
procedure dump (solvpb : solvpb_);
Begin
  with solvpb do
  begin
    writeln ('dump: solvbp begin -------------');
    writeln (u1 : 8 : 5);
    writeln (u2 : 8 : 5);
    writeln (u3 : 8 : 5);
    writeln (u4 : 8 : 5);
    writeln (u12 : 8 : 5);
    writeln (u13 : 8 : 5);
    writeln (u14 : 8 : 5);
    writeln (u23 : 8 : 5);
    writeln (u24 : 8 : 5);
    writeln (u34 : 8 : 5);
    writeln (u123 : 8 : 5);
    writeln (u124 : 8 : 5);
    writeln (u134 : 8 : 5);
    writeln (u234 : 8 : 5);
    writeln ('dump: solvbp end ------------');
  end;
End;
(************************)
function solv3 (var solvpb : solvpb_; clus : integer) : real;

  (* solve the 2^3 table *)
const
  epsilon = 1.0E-6;
  maxiter = 100; (* max # iterations in Newton-Raphson *)

var
  iter : integer;
  lower, upper, dx, x, f, fp : real;
  solvflag : BOOLEAN;
(*---------------- nested proc -> ---------------------------------*)
procedure getffp ( x : real; var f, fp : real; var solvpb : solvpb_);
  (* x = u123 = c7 *)
  (* f=function, fp=derivative *)
var
  c0,c1,c2,c3,c4,c5,c6 : real;

Begin (* solv3.getffp *)
  with solvpb do
  begin
    c0 := 1 + u13 - u1 + u12 - x - u2 + u23 - u3;
    c1 := u1 - u12 + x - u13;
    c2 := u2 - u12 + x - u23;
    c3 := u12 - x;
    c4 := u3 - u13 - u23 + x;
    c5 := u13 - x;
    c6 := u23 - x;
  end;
  f := c0*c3*c5*c6-c1*c2*c4*x;
  fp := -(c3*c5*c6+c0*c5*c6+c0*c3*c6+c0*c3*c5
        +c2*c4*x +c1*c4*x +c1*c2*x +c1*c2*c4);
  (* <- not efficient but the ONLY safe way to do it *)

  if solvflag then  (* print all cells *)
  begin
    writeln('solv3:');
    writeln(c0 :8:5, '':1, c2 :8:5);
    writeln(c1 :8:5, '':1, c3 :8:5);
    writeln(c4 :8:5, '':1, c6 :8:5);
    writeln(c5 :8:5, '':1, x  :8:5);
  end;

End; (* solv3.getffp *)
(*------------- <- nested proc --------------------------------------*)


Begin (* solv3 *)

  with solvpb do
  begin      (* find bounds on u123 *)
    upper := min(u12, u13, u23);
    lower := max(0, u13+u23-u3, u12+u23-u2, u12+u13-u1);
  end;

  if (lower > upper) then
  begin
    writeln('solv3:');
    writeln('  id= ', id(.strt(.clus.).) );
    writeln('  lower = ', lower);
    writeln('  upper = ', upper);
    writeln('  lower > upper');
  end;


  solvflag := FALSE;
  x := 0.5 * (lower + upper);   (* starting value *)
  dx := 9999;
  iter := 0;

  while (iter < maxiter) and (abs(dx) > epsilon) do
  begin
    iter := succ(iter);
    getffp (x, f, fp, solvpb); (* f=function, fp=derivative *)

    if indebug   then     (* debugging code *)
    begin
      writeln('solv3: x = ', x:18:15);
      writeln('solv3:     f = ', f:18:15);
      writeln('solv3:     fp= ', fp:18:15);
    end;

    dx := - f / fp;
    x := x + dx;
  end;

  if (x <= lower) or (x >= upper) then    (* out of range *)
  begin
    writeln('solv3: id = ', id(.strt(.clus.).):6);
    writeln('       Out of range x = ', x:18:15);
    writeln('       Range = (', lower:18:15, ' , ', upper:18:15, ')');
    x := 0.5 * (lower + upper);     (* starting value *)
  end;

  if (abs(dx) > epsilon) then   (* ? no convergence *)
  begin
    writeln('solv3: id = ', id(.strt(.clus.).):5);
    writeln('       large dx = ', dx:18:15);
  end;

  if indebug   then      (* debugging code *)
  begin
    solvflag := TRUE; (* to print cells *)
    getffp (x, f, fp, solvpb);
  end;

  solv3 := x;  (* the solution *)

End; (* solv3 *)
(************************)
function solv4 (var solvpb : solvpb_; clus : integer) : real;

  (* solve the 2^4 table *)
const
  epsilon = 1.0E-6;
  maxiter = 100; (* max # iterations in Newton-Raphson *)

var
  iter : integer;
  lower, upper, dx, x, f, fp : real;
  solvflag : BOOLEAN;

(*---------------- nested proc -> ---------------------------------*)
procedure getffp ( x : real; var f, fp : real; var solvpb : solvpb_);
  (* x = u1234 = c15 *)
  (* f=function, fp=derivative *)
var
  c0,c1,c2,c3,c4,c5,c6,
  c7,c8,c9,c10,c11,c12,c13,c14,
  p1, p2, p3, p4,
  p1p, p2p, p3p, p4p : real;

Begin (* solv4.getffp *)
  with solvpb do
  begin
    c0 :=  1 + u34 - u1 + u13 + u12 - u124 + u14 - u123
             - u134 + x - u2 - u3 + u24 + u23 - u234 - u4;
    c1 := u1 - u13 - u12 + u124 - u14 + u123 + u134 - x;
    c2 := u2 + u124 - u12 + u123 - u24 + u234 - x - u23;
    c3 := u12 - u124 + x - u123;
    c4 := u3 - u34 - u13 + u123 - u23 + u134 + u234 - x;
    c5 := u13 - u123 - u134 + x;
    c6 := u23 - u123 - u234 + x;
    c7 := u123 - x;
    c8 := u4 - u34 - u14 - u24 + u124 + u134 + u234 - x;
    c9 := u14 - u124 + x - u134;
    c10 := u24 - u124 - u234 + x;
    c11 := u124 - x;
    c12 := u34 - u134 - u234 + x;
    c13 := u134 - x;
    c14 := u234 - x;
  end;

  p1 := c0*c3*c5*c6;
  p2 := c9*c10*c12*x;
  p3 := c1*c2*c4*c7;
  p4 := c8*c11*c13*c14;

  p1p := c3*c5*c6 + c0*c5*c6 + c0*c3*c6 + c0*c3*c5;
  p2p := c10*c12*x + c9*c12*x + c9*c10*x + c9*c10*c12;
  p3p := -(c2*c4*c7 + c1*c4*c7 + c1*c2*c7 + c1*c2*c4);
  p4p := -(c11*c13*c14 + c8*c13*c14 + c8*c11*c14 + c8*c11*c13);

  f := p1*p2 - p3*p4;
  fp := p1*p2p + p1p*p2 - (p3*p4p + p3p*p4);
  (* <- not efficient but the ONLY safe way to do it *)

  if solvflag then     (* debugging  code *)
  begin
    writeln('solv4:');
    writeln(c0 :8:5, '':1, c2 :8:5, '':1, c8 :8:5, '':1, c10:8:5);
    writeln(c1 :8:5, '':1, c3 :8:5, '':1, c9 :8:5, '':1, c11:8:5);
    writeln(c4 :8:5, '':1, c6 :8:5, '':1, c12:8:5, '':1, c14:8:5);
    writeln(c5 :8:5, '':1, c7 :8:5, '':1, c13:8:5, '':1, x  :8:5 );
  end;

End; (* solv4.getffp *)
(*------------- <- nested proc --------------------------------------*)

Begin (* solv4 *)

  with solvpb do   (* bounds on u1234 *)
  begin
    upper := min(u123, u124, u134, u234);
    lower := max(0, u134+u234-u34, u124+u234-u24,
                    u123+u234-u23, u124+u134-u14,
                    u123+u134-u13, u123+u124-u12);
  end;

  if (lower > upper) then
  begin
    writeln('solv4:');
    writeln('  id= ', id(.strt(.clus.).) );
    writeln('  lower = ', lower);
    writeln('  upper = ', upper);
    writeln('  lower > upper');
  end;

  solvflag := FALSE;
  x := 0.5 * (lower + upper);     (* starting value *)
  dx := 9999;
  iter := 0;

  while (iter < maxiter) and (abs(dx) > epsilon) do
  begin         (* Newton-Raphson *)
    iter := succ(iter);
    getffp (x, f, fp, solvpb);

    if indebug   then    (* debugging code *)
    begin
      writeln('solv4: x = ', x:18:15);
      writeln('solv4:      f = ', f:18:15);
      writeln('solv4:      fp= ', fp:18:15);
    end;

    dx := - f / fp;
    x := x + dx;

  end;   (* Newton-Raphson *)

  if (x <= lower) or (x >= upper) then    (* out of range *)
  begin
    writeln('solv4: id = ', id(.strt(.clus.).):6);
    writeln('       Out of range x = ', x:18:15);
    writeln('       Range = (', lower:18:15, ' , ', upper:18:15, ')');
    x := 0.5 * (lower + upper);     (* starting value *)
  end;

  if (abs(dx) > epsilon) then   (* ? no convergence *)
  begin
    writeln('solv4: id = ', id(.strt(.clus.).):5);
    writeln('       large dx = ', dx:18:15);
  end;

  if indebug   then      (* debugging *)
  begin
    solvflag := TRUE; (* to print cells *)
    getffp (x, f, fp, solvpb);
  end;

  solv4 := x;  (* the solution *)

End; (* solv4 *)
(************************)
procedure calcu4 (c : integer);
  (* compute 4th order moments for cluster c and puts results
   into the global array uijkla(.i,j,k,l.)  *)

var
  solvpb : solvpb_;  (* parameter block for proc solv3  *)
  n, i, j, k, l : integer;

Begin (* calcu4 *)
  n := siz(.c.);
  with solvpb do
  for i := 1 to (n-3) do
  begin
    u1 := ui (i, c);
    for j := succ(i) to (n-2) do
    begin
      u2 := ui (j, c);
      u12 := uij (i, j, c);
      for k := succ(j) to pred(n) do
      begin
        u3 := ui (k, c);
        u13 := uij (i, k, c);
        u23 := uij (j, k, c);
        u123 := uijka (.i, j, k.);
        for l := succ(k) to n do
        begin
          u4 := ui (l, c);
          u14 := uij (i, l, c);
          u24 := uij (j, l, c);
          u34 := uij (k, l, c);
          u124 := uijka (.i, j, l.);
          u134 := uijka (.i, k, l.);
          u234 := uijka (.j, k, l.);
          uijkla(.i, j, k, l.) := solv4 (solvpb, c);
        end
      end
    end
  end
End; (* calcu4 *)
(************************)
procedure calcu3 (c : integer);
  (* compute 3rd order moments for cluster c and puts results
   into the global array uijka(.i,j,k.)  *)

var
  solvpb : solvpb_;  (* parameter block for proc solv3  *)
  n, i, j, k : integer;

Begin (* calcu3 *)
  n := siz(.c.);
  with solvpb do
  for i := 1 to (n-2) do
  begin
    u1 := ui (i, c);
    for j := succ(i) to pred(n) do
    begin
      u2 := ui (j, c);
      u12 := uij (i, j, c);
      for k := succ(j) to n do
      begin
        u3 := ui (k, c);
        u13 := uij (i, k, c);
        u23 := uij (j, k, c);
        uijka(.i, j, k.) := solv3 (solvpb, c);
      end
    end
  end
End; (* calcu3 *)
(************************)
procedure calcb11 (var b : matrix;
                       c : integer);
  (* compute the upper half of the (1,1) part of
  the cov matrix for cluster c.
  b11 is symmteric *)
  (*       called only if b is full or block-diag *)
var
  i, j, ij, n, base : integer;

Begin (* calcb11 *)
  (* off-diagonals *)
  base := pred(strt(.c.));   (* 0-base *)
  ij := fnsh(.c.);             (* pointer to product terms *)
  n := siz(.c.);
  with b do
  for i := 1 to pred(n) do for j := succ(i) to n do
  begin
    ij := succ(ij);
    el(.i,j.) := u(.ij.) - u(.i+base.)*u(.j+base.);
  end;

  (* diagonals *)
  i := strt(.c.);   (* pointer to y *)
  with b do
  for j := 1 to n do
  begin
    el(.j,j.) := u(.i.) * (1-u(.i.));
    i := succ(i);
  end;
End;  (* calcb11 *)
(************************)
procedure calcb12 (var b : matrix;
                       c : integer);
  (* compute the (1,2) part of the cov matrix for cluster c *)
  (* note: b12 is asymmetric  n x (n+ch2(.n.)) *)
var
  row, col, k, l, n : integer;

Begin (* calcb12 *)
  n := siz(.c.);
  with b do
  for row :=  1 to n do
  begin
    k := 1;   (* (k,l) is an index pair *)
    l := 2;   (* first column in each row corres.  to (1,2) *)
    for col := succ(n) to (n+ch2(.n.)) do
    begin
      el(.row,col.) := covyiyjk(row,k,l,c);
      incr (k,l,n); (* move to the right *)
    end; (* for col *)
  end; (* for row *)
End;  (* calcb12 *)
(************************)
procedure calcb22 (var b : matrix;
                       c : integer);
  (* compute the upper half of (2,2) part of the cov matrix
   for cluster  c.
   b22 is symmetric. *)
var
  i, j, k, l, n, row, col, ij : integer;
  uij, ui : real;

Begin (* calcb22 *)
  n := siz(.c.);

  (* off diagonals *)
  i := 1; j := 2; (* row index-pair *)
  with b do
  for row :=  succ(n) to  pred((n+ch2(.n.))) do
  begin
    k := i;
    l := j;
    for col := succ(row) to (n+ch2(.n.)) do
    begin
      incr (k, l, n);  (* start to right of diagonal on each row *)
      el(.row,col.) := covyijykl(i,j,k,l,c);
    end;
    incr (i, j, n); (* next row index pair *)
  end;

  (* diagonals *)
  ij := fnsh(.c.);  (* pointer into product terms *)
  with b do
  for row := succ(n) to (n+ch2(.n.)) do
  begin
    ij := succ(ij);
    uij := u(.ij.);
    el(.row,row.) := uij*(1-uij);
  end;

End;  (* calcb22 *)
(***********************)
procedure calcg11 (var g : matrix;
                       c : integer);
  (* note: g is (n+nC2)x(n+nC2)  asymmetric *)
 (* compute the g11 part of g for cluster c *)
 (* g11 is diagonal  ui*(1-ui)   *)

var
  i, p : integer;
  ui : real;
Begin (* calcg11 *)
  p := strt(.c.);
  with g do
  for i := 1 to siz(.c.) do
  begin
    ui := u(.p.);
    el(.i,i.) := ui * (1-ui);
     p := succ(p);
  end;
End; (* calcg11 *)
(***********************)
procedure calcg21 (var g : matrix;
                       c : integer);
  (* note: g is (n+nC2)x(n+nC2)  asymmetric *)
  (* compute the g21 part of g for cluster c *)
  (* each row in g21 has exactly 2 nonzero entries *)
var
  row, i, j, ij, n, base : integer;
  ui, uj, psiij : real;

Begin (* calcg21 *)
  n := siz(.c.);
  base := pred(strt(.c.));
  ij := fnsh(.c.); (* pointer into product terms *)
  i := 1; j := 2; (* index pair *)
  with g do
  for row := succ(n) to (n+ch2(.n.)) do
  begin (* process one row *)
    ij := succ(ij);  (* next product term *)
    ui := u(.base+i.);
    uj := u(.base+j.);
    psiij := psi (.ij.);
    el(.row,i.) := ui*(1-ui)*duijdui(ui, uj, psiij);
    el(.row,j.) := uj*(1-uj)*duijdui(uj, ui, psiij);
    incr(i,j,n); (* next row *)
  end;
End; (* calcg21 *)
(***********************)
procedure calcg22 (var g : matrix;
                       c : integer);
  (* note: g is (n+nC2)x(n+nC2)  asymmetric *)
  (* compute the g22 part of g for cluster c *)
  (* g22 is diagonal *)
var
  row, i, j, ij : integer;

Begin (* calcg22 *)
  ij := fnsh(.c.);
  row := siz(.c.);
  with g do
  for i := strt(.c.) to pred(fnsh(.c.)) do
  for j := succ(i) to fnsh(.c.) do
  begin
    row := succ(row);  ij := succ(ij);
    el(.row,row.) := duijdpsi(u(.i.), u(.j.), psi(.ij.)) * psi(.ij.);
  end;
End; (* calcg22 *)
(***********************)
procedure calcg (var g : matrix;
                     c : integer);
  (* compute the G matrix for cluster c *)
  (* G = the derivative of the mean
         with respect to the linear predictor *)

var n : integer;
    ui : real;

Begin (* calcg *)
  n := siz(.c.);

  if (n=1) then    (* cluster size 1 *)
  begin
    matdim (g, 1, 1);
    ui := u(.strt(.c.).);
    g.el(.1,1.) := ui * (1-ui);
    RETURN;
  end;

  matdim (g, n+ch2(.n.), n+ch2(.n.));
  matfill (g, 0);
  calcg11 (g, c);
  calcg21 (g, c);
  calcg22 (g, c);
End;  (* calcg *)
(***********************)
procedure calcb (var b, v : matrix;
                     c : integer);
  (* compute the B matrix for cluster c *)
  (* v <- cov (y,w)  *)
  (* B <- inverse(v) *)
var
  i, row, n, m : integer;
  ui : real;
  mfile : file of matrix;

Begin (* calcb *)
  n := siz(.c.);
  m := n+ch2(.n.);
  matdim (v, m, m);
  matdim (b, m, m);

  if (n=1) then    (* cluster size 1, quick *)
  begin
    ui := u(.strt(.c.).);
    v.el(.1,1.) := ui * (1-ui);
    b.el(.1,1.) := 1.0 / v.el(.1,1.);
    RETURN;    (* done *)
  end;

  calcb11 (v, c);  (* same for GEE1 & GEE2 *)

  (* computing the other 2 blocks of v depends on the method being
     used to compute means of products of 3 and 4 variables.
    The dependence is built into the covy... functions
    so the general structure for computing v is the same *)

  if (vmethod = VMBQ) then  (* prepare high order means *)
  begin
    if (n >= 3) then calcu3 (c);  (* E (Yijk) *)
    if (n >= 4) then calcu4 (c);  (* E (Yijkl) *)
  end;

  if (geelvl = 2) then calcb12 (v, c)
    else blockout (v, 1, n, succ(n), m);

  calcb22 (v, c);

  MatU2L  (v); (* reflect upper half to lower half *)

  (* now: v == the complete covariance matrix V *)

  b := v;

  (* if we are called by proc debug we will not invert v *)
  if not indebug then
  begin
    if debugmode then
    begin
      writeln(StdErr, 'calcb: MatSInv(b) iter, c = ',
              iter:3, '':1, c:3);
      writeln(StdErr, 'n = ':5, n:3);
    end;
    MatSInv (b);
    if MatFail then
    begin
      rewrite (mfile);
      write(mfile, v);
      close(mfile);
      writeln('calcb: iter c id n = ',
              iter:3, '':1, c:3, '':1,
              id(.strt(.c.).):3, '':1, siz(.c.):3);
      err('calcb: MatSInv(b) failed.', HALTYE);
    end;
  end;

  (* now: b == the inverse of the complete covariance matrix v *)

End;  (* calcb *)
(************************)
procedure bldtyptbl (ntyp : integer);
(* build type table *)
var
  i, j, t : integer;

Begin (* bldtyptbl *)
  for i := 1 to ntyp do typtbl(.0,i.) := i;
  for i := 1 to ntyp do typtbl(.i,0.) := i;
  t := ntyp;
  for i := 1 to ntyp do
  for j := i to ntyp do
  begin
    t := succ (t);
    typtbl(.i,j.) := t;
    typtbl(.j,i.) := t;
  end;
End;  (* bldtyptbl *)
(************************)
procedure bldx (var row    : integer;
                  idd, t   : integer;
                    yy     : real;
                var xx     : xx_);
  (* build a row into y and x *)
var
  j : integer;
Begin (* bldx *)
  (* error checking first *)
  (* assume:   1 <= t <= MAXXTYP *)
  if not ((yy = 0) or (yy = 1)) then
    err('bldx: y must be 0/1.', HALTYE);
  if (row > MAXROWSX) then
    err('bldx: data set is too large.', HALTYE);

  (* build x from xx *)
  for j := 1 to dimb do x(.row,j.) := 0;
  with reg(.t.) do
  begin
    for j := 1 to dim do
      x(.row,bndx(.j.).) := xx(.xndx(.j.).);
    (* note: if dim=0 the x row <- 0's which causes lp <- 0 *)
  end;

  y(.row.) := yy;
  typ(.row.) := t;
  id(.row.) := idd;
  row := succ(row);   (* the only control of the row pointer *)
End; (* bldx *)
(************************)
procedure expand (c, id : integer;
                var xx     : xx_;
                var row    : integer);
  (* expand data of cluster c, start at row *)
var
  i, j, t : integer;
  yy   : real;
Begin (* expand *)
  if (siz(.c.) = 1) then RETURN;
  for i := strt(.c.) to pred(fnsh(.c.)) do
  for j := succ(i) to fnsh(.c.) do
  begin
    t := typtbl(.typ(.i.),typ(.j.).);
    yy := y(.i.) * y(.j.);
    bldx (row, id, t, yy, xx);
  end;
End;  (* expand *)
(************************)
procedure  chkt (t:integer); (* 1 <= t <= ntyp *)
Begin (* chkt *)
  if (t < 1) or (t > ntyp) then
  err('chkt: class out of range.', HALTYE);
End; (* chkt *)
(************************)
procedure readdata;

var
  i, j, clus, curid, newid, t : integer;
  xx : xx_;
  y : real;
Begin (* readdata *)

  (* initialize while loop *)
  i := 1; (* the row pointer *)
  clus := 1;
  strt(.1.) := 1;

  read (datafile, curid, t, y);
  for j := 1 to dimx do read (datafile, xx(.j.));
  readln (datafile);

  chkt (t); (* 1 <= t <= ntyp *)
  bldx (i, curid, t, y, xx);

  while not EOF (datafile) do
  begin
    read (datafile, newid, t, y);
    if (newid <> curid) then
    begin
      fnsh(.clus.) := pred (i);
      siz(.clus.) := succ(fnsh(.clus.)) - strt(.clus.);
      expand (clus, curid, xx, i); (* expand cluster *)
      clus := succ (clus);
      strt (.clus.) := i;
      curid := newid;
    end;

    for j := 1 to dimx do read (datafile, xx(.j.));
    readln (datafile);

    chkt (t); (* 1 <= t <= ntyp *)
    bldx (i, curid, t, y, xx);

  end;  (* while not EOF *)

  (* process last cluster *)
  fnsh(.clus.) := pred (i);
  siz(.clus.) := succ(fnsh(.clus.)) - strt(.clus.);
  nclus := clus;
  expand (clus, curid, xx, i); (* expand cluster *)
  expn := pred (i);

  writeln('# clusters: ', nclus:8);
  close (datafile);

End;  (* readdata *)
(************************)
procedure chktyp;
  (* check range of primary types *)
var
  clus, i : integer;
Begin (* chktyp *)
  for clus := 1 to nclus do
  for i := strt(.clus.) to fnsh(.clus.) do
  begin
    writeln('chktyp: i=',i,' typ(.i.)= ', typ(.i.));
    if (typ(.i.) > ntyp) or (typ(.i.) > MAXNTYP)
      or (typ(.i.) < 1) then
      err('chktyp: Invalid class.', HALTYE);
  end;
End;  (* chktyp *)
(************************)
procedure chksiz;   (* check and list cluster sizes *)
var
  i : integer;
Begin (* chksiz *)
  for i := 1 to nclus do
  if (siz(.i.) > MAXSIZ) then
    err('chksiz: Maximum cluster size exceeded.', HALTYE);

  (* do frequencies on cluster sizes *)
  for i := 1 to MAXSIZ do sizf(.i.) := 0; (* freqs *)
  for i := 1 to nclus do
    sizf(.siz(.i.).) := succ(sizf(.siz(.i.).));

  (* write frequencies *)
  writeln('cluster size':-12, '':10, '#clusters':-9);
  for i := 1 to MAXSIZ do if (sizf(.i.) > 0) then
    writeln(i:6, '':19, sizf(.i.):4);

End;  (* chksiz *)
(************************)
procedure bldch2;  (* duild the n choose 2 in array ch2 *)
var i : integer;
Begin (* bldch2 *)
  ch2(.1.) := 0;
  for i := 2 to MAXSIZ do ch2(.i.) :=  (sqr(i)-i) DIV 2;
End;  (* bldch2 *)
(************************)
procedure  chk1 (dimx, dimb, ntyp : integer;
                 var maxiter: integer;
                     tol : real);
Begin   (* chk1 *)
  if (maxiter < MAXITERLO) then maxiter := MAXITERLO else
  if (maxiter > MAXITERHI) then maxiter := MAXITERHI;

  writeln('dim(x):     ', dimx:8);
  writeln('dim(b):     ', dimb:8);
  writeln('# classes   ', ntyp:8);
  writeln('max #iter:  ', maxiter:8);
  writeln('tolerance:  ', tol:15:13);

  if  ((ntyp < 1) or  (ntyp > MAXNTYP)) then
  err('chk1: # classes out of range.', HALTYE);

  if (dimb < 1) or (dimb > MAXDIMB) then
  err('chk1: dim(b) out of range.', HALTYE);

  if (dimx < 1) or (dimx > MAXINPX) then
  err('chk1: dim(x) out of range.', HALTYE);

  xtyp := ntyp + ntyp + (sqr(ntyp)-ntyp) DIV 2;

End; (* chk1 *)
(************************)
function validptyp (i : integer) : BOOLEAN;
  (* is i a valid primary type ? *)
Begin
  validptyp := (i>=1) and (i<=ntyp)
End;
(**************************)
procedure parse (s : string_;
                 var reg1 : reg1_);
  (* s contains pairs of (beta index, x index) *)
var
  l, p, bn , xn : integer;
  s1 : alpha;
begin (* parse *)
  with reg1 do
  begin
    if defined then
    err('parse: Regression redefinition attempted.', HALTYE) else
    defined := TRUE;
    l := length (s);
    if (l=0) then err('parse: Incomplete record.', HALTYE);
    p := 1;
    dim := 0;
    repeat
      dim := succ (dim);

      Token (p, s, s1);
      ReadStr (Str(s1), bn);
      if (bn < 0) or (bn > dimb) then
      err('parse: beta index out of range.', HALTYE);
      if (bn = 0) then  (* lp = 0 *)
      begin
        dim := 0;
        RETURN; (* don't look further *)
      end;

      if (p>l) then err('parse: Incomplete record.', HALTYE);
      Token (p, s, s1);
      ReadStr (Str(s1), xn);
      if (xn < 1) or (xn > dimx) then
      err('parse: x index out of range.', HALTYE);

      if (dim > dimb) then err('parse: Too many arguments.', HALTYE);
      bndx(.dim.) := bn;
      xndx(.dim.) := xn;

    until (p > l);
  end; (* with *)
end; (* parse *)
(**************************)
procedure chktypdef;  (* all regs defined ? *)

var
  t : integer;
Begin (* chktyp *)
  for t  := 1 to xtyp do if not reg(.t.).defined then
  err('chktypdef: Incomplete regression definitions.', HALTYE);
End;  (* chktyp *)
(**************************)
procedure chkused;
  (* make sure all betas are used *)
var
  i, j, t : integer;
  s : string_;
Begin (* chkused *)
  for i := 1 to dimb do used(.i.) := FALSE;
  for t := 1 to xtyp do
    with reg(.t.) do
      for i := 1 to dim do
        used(.bndx(.i.).) := TRUE;
  for i := 1 to dimb do
    if not used(.i.) then
      err('chkused: Some betas are not used.', HALTYE);
End;  (* chkused *)
(**************************)
procedure parsendx; (* parse the indices from the input file *)

var
  t, t1, t2  : integer;
  s : string_;

Begin (* parsendx *)
  for t := 1 to xtyp do reg(.t.).defined := FALSE;

  for t := 1 to xtyp do
  begin
    readln(ctrlfile, t1, t2, s);
    if (t1 > t2) then swap(t1, t2);

    if (validptyp(t1) and validptyp(t2)) or
       ((t1=0) and validptyp(t2)) then
    begin
      reg(.typtbl(.t1,t2.).).t1 := t1;
      reg(.typtbl(.t1,t2.).).t2 := t2;
      parse (s, reg(.typtbl(.t1,t2.).))
    end else
    err('parsendx: Invalid class specified', HALTYE);
  end;

  chktypdef; (* all regs defined ? *)

  chkused; (* all betas used ?  *)

End;  (* parsendx *)
(************************)
procedure readb; (* read initial values of beta *)
var
  j : integer;
Begin (* readb *)
  b.rows := dimb;
  with b do
  for j := 1 to rows do read (ctrlfile, ve(.j.));
  readln (ctrlfile);
End;  (* readb *)
(************************)
procedure readlbl; (* read beta labels *)
var i : integer;
Begin (* readlbl *)
  for i := 1 to dimb do readln (ctrlfile, lbl(.i.));
End;  (* readlbl *)
(************************)
procedure initpsi; (* initialize psi to -9 *)
var i : integer;
Begin (* initpsi *)
  for i := 1 to expn do psi(.i.) := -9;
End;  (* initpsi *)
(************************)
procedure open4input; (* open control and data files *)
var
  cmd : string_;
  rc : integer;
Begin (*  open4input *)
  readln (ctrlfile, cmd);
  cmd := compress(cmd);
  writeln ('Data file: ', cmd);
  cmd := 'filedef datafile disk ' || cmd;
  cms (cmd, rc);
  if (rc <> 0) then
    err ('Invalid file specification: ' || cmd, HALTYE);
  reset (datafile);
End;  (*  open4input *)
(************************)
procedure repmethods; (* report methods used *)

Begin (* repmethods *)
  write('Method of computing 3rd and 4th order moments: ');
  if (vmethod=VMZP) then
    writeln('Zhao & Prentice approximation.') else
  if (vmethod=VMBQ) then
    writeln('The exact method.') else
  err('Invalid method for 3rd and 4th order moments.', HALTYE);

  if (geelvl = 1) then
    writeln('-- GEE1 --') else  writeln('-- GEE2 --');

  writeln;
End;  (* repmethods *)
(************************)
procedure  showmodel;  (* print model spec *)
var
  t,i : integer;

Begin (* showmodel *)
  writeln('The model:');
  for t := 1 to xtyp do
  with reg(.t.) do
  begin
    if (t1 <> 0) then write(t1:1,'.':1);
    write(t2:1,' : ':3);
    if (dim = 0) then write('0':1) else
    for i := 1 to dim do
    write('B':1,bndx(.i.):1,'':1, 'X':1,xndx(.i.):1, '  ':2);
    writeln;
  end;
  writeln;

  repmethods; (* report methods used *)

End;  (* showmodel *)
(************************)
procedure init;
var
  title : string_;
  i1, i2, i3 : integer;
  uf  : file of real;   (* mean *)

Begin (* init *)

  CPUtime := Clock;

  writeln (VERSION1);
  writeln (VERSION2);
  writeln (VERSION3);
  writeln;
  DateTime (date,  time);
  writeln('Date: ', date, '    Time: ', time);
  writeln;

  reset (ctrlfile);
  readln(ctrlfile, i1, i2, i3, bias);
  debugmode := (i1 = 1);
  debugclus := i2;
  probescore := (i3 = 1);

  readln (ctrlfile, title);
  writeln (title);
  readln (ctrlfile, title);
  writeln (title);

  open4input;  (* open data file for input *)

  readln (ctrlfile, ntyp);    (* # types *)
  readln (ctrlfile, dimx);    (* dim(x)  *)
  readln (ctrlfile, dimb, dimb1);    (* dim(b), dim(b) for GEE1 *)

  readln (ctrlfile, tol);
  readln (ctrlfile, maxiter);

  readln (ctrlfile, i1);
  monbflag := (i1=1);  (* report beta each iteration *);

  readln (ctrlfile, geelvl); (*  1=GEE1, 2=GEE2 *)
  if (geelvl  < 1) or (geelvl > 2) then
  err('init: GEE level must be either 1 or 2.', HALTYE);

  readln (ctrlfile, vmethod);
  if (vmethod < VMZP) or (vmethod > VMBQ) then
  err('init: Invalid method for 3rd and 4th order moments.', HALTYE);

  readln (ctrlfile); (* skip a line *)

  chk1 (dimx, dimb, ntyp, maxiter, tol); (* check & compute xtyp *)

  bldtyptbl (ntyp);      (* build translation table *)

  readlbl;          (* labels of betas *)

  readln (ctrlfile); (* skip a line *)

  readb;    (* initial values of beta *)

  readln (ctrlfile); (* skip a line *)

  bldch2;  (* n choose 2 in array ch2 *)

  parsendx; (* parse index tables that define all regressions *)

  close (ctrlfile);   (* control file done *)

  readdata; (* read actual data: id type y x *)

  (* bias code *)
  if (bias = 2) then
  begin
    reset (uf);
    for i1 := 1 to expn  do
    read  (uf,y(.i1.));
    close (uf);
  end;
  (* bias code *)

  chksiz;   (* check cluster sizes *)

  showmodel;

  initpsi; (* initialize psi to some constant, aids debugging *)

End; (* init *)
(************************)
procedure  dumpdata;  (* print all data *)
var
  i, j : integer;
Begin (* dumpdata *)
  for i := 1 to expn do
  begin
    write  ('dumpdata:'   , '':1,
            i       :3, '':1,
            id(.i.) :3, '':1,
            typ(.i.):3, '':1,
            y(.i.):3:0, '':1);
    for j := 1 to dimb do write (x(.i,j.):3:0, '':1);
    writeln;
  end;
End;  (* dumpdata *)
(************************)
procedure calcr (var v : vector ; c : integer);
(* v <- residuals  of cluster c *)
var
  i, j : integer;

Begin (* calcr*)

  v.rows :=  (siz(.c.)+ch2(.siz(.c.).));
  i := 1;
  with v do
  for j := strt(.c.) to (fnsh(.c.)+ch2(.siz(.c.).)) do
  begin
    ve (.i.) := y(.j.) - u(.j.);
    i := succ (i);
  end;
End;  (* calcr*)
(************************)
procedure calcxtgt (var xtgt, g : matrix;
                        c : integer);     (* xtgt <- X'G' *)
var
  i, j, k, kk, n : integer;
  w : real;

Begin (*  calcxtgt *)
  n := siz(.c.) + ch2(.siz(.c.).);
  MatDim (xtgt, dimb, n);

  with g do
  for i := 1 to n do for j := 1 to dimb do
  begin
    w := 0;
    kk := strt(.c.); (* start of data for cluster c *)
    for k := 1 to n do
    begin
      w := w + el(.i,k.)*x(.kk,j.);
      kk := succ(kk);
    end;
    xtgt.el(.j,i.) := w;
  end;

End;  (*  calcxtgt *)
(************************)
procedure debug (c : integer);
  (* inspect important data-structures *)

var
  vmat, bmat, gmat, xtgt: matrix;
  r : vector;
  i, j, k, l, n : integer;

Begin (* debug *)
  if not debugmode then
  begin
    indebug := FALSE;
    RETURN;
  end;

  indebug := TRUE; (* debug mode on *)

  writeln ('debug: cluster = ':-18, c:4);

  n := siz(.c.);
  calcr (r, c);

  writeln (
    '  id    t    y  lp         psi      u      (y-u)       r'
          );
  j := 1;
  for i := strt(.c.) to (fnsh(.c.) + ch2(.n.)) do
  begin
    writeln(id(.i.) :4, '':1,
           typ(.i.) :4, '':1,
            y (.i.) :4:0, '':1,
            lp(.i.) :8:5, '':1,
           psi(.i.) :8:5, '':1,
            u (.i.) :8:5, '':1,
    y(.i.)- u (.i.) :8:5, '':1,
          r.ve(.j.) :8:5, '':1);
    j := succ(j);
  end;

  calcb (bmat, vmat, c);
  calcg (gmat, c);
  calcxtgt (xtgt, gmat, c);     (* xtgt <- X'G' *)

  writeln('debug: ui = ');
  for i := 1 to n do
  writeln(i:3, '':1, ui(i,c):18:15);

  writeln('debug: uij = ');
  for i := 1 to pred(n) do
  for j := succ(i) to n do
  writeln(i:3, j:3, '':1, uij(i, j, c):18:15);

  writeln('debug: uijk = ');
  for i := 1 to n-2     do
  for j := succ(i) to n-1  do
  for k := succ(j) to n    do
  writeln(i:3, j:3, '':1, k:3, '':1, uijka(.i, j, k.):18:15);

  writeln('debug: uijkl = ');
  for i := 1 to n-3     do
  for j := succ(i) to n-2  do
  for k := succ(j) to n-1  do
  for l := succ(k) to n    do
  writeln(i:3, j:3, '':1, k:3, '':1, l:3, '':1,
          uijkla(.i, j, k, l.):18:15);

  writeln ('debug: Cov mtx = ':-18);
  MatWrite (output, bmat);

  writeln ('debug: G   mtx = ':-18);
  MatWrite (output, gmat);

  writeln ('debug: X''G''    = ':-18);
  MatWrite (output, xtgt);

  indebug := FALSE;      (* debug mode off *)
End; (* debug *)
(************************)
procedure probe (var score : vector);
  (* print the score vector if requested *)
Begin   (* probe *)
  if probescore then
    VecWrite (probef, score);    (* print score vector *)
End;    (* probe *)
(************************)
procedure AccGEE2 (var s1 : matrix;
                     var s2 : vector;
                     var s3 : matrix;
                     var s1a: matrix;
                     s3flag : BOOLEAN;
                         c  : integer);

var
  v, w1, w2, w3, gx   : matrix;
  v1, v2  : vector;

Begin (* AccGEE2 *)
  if (c=debugclus) then debug(c); (* dump a cluster's data *)

  calcr (v1, c);            (* v1 <- residuals  *)
  calcg (w3, c);            (* w3 <- G *)
  calcxtgt (w1, w3, c);     (* w1 <- X'G' *)
  if (geelvl = 1) then      (* GEE1: modify X'G' *)
  begin
    MatTran (gx,   w1);       (* the unabridged GX *)
    ZeroOut  (w1, 1, dimb1, succ(siz(.c.)), w1.cols);
  end;
  calcb (w2, v, c);         (* w2 <- B = inv(v), according to geelvl *)

  (* If GEE1:  w2 and w1 hold abridged versions of B and G'X'
           while gx holds the original GX  *)

  MatMul (w3, w1, w2);      (* w3 <- X'G'B *)
  MatxVec (v2, w3, v1);     (* v2 <- X'G'Br *)
  VecAdd (s2, s2, v2);      (* s2 <- s2 + X'G'Br *)
  MatTran (w2, w1);         (* w2 <- GX *)

  MatMul (w1, w3, w2);      (* w1 <- X'G'BGX *)
  MatAdd (s1, s1, w1);      (* s1 <- s1 + X'G'BGX *)

  (* s3 and s1a are needed only for variance estimation, but not during
     the iteration *)
  if s3flag then
  begin
    if (geelvl = 1) then      (* GEE1: compute s1a *)
    begin
      MatMul (w1, w3, gx);      (* w1 <- X'G'BGX *)
      MatAdd (s1a, s1a, w1);    (* s1a <- s1a + X'G'BGX *)
    end;
    AddOP (s3, v2);         (* s3 <- s3 + X'G'Brr'BGX *)
    probe (v2);             (* special code: dump score vector *)
  end;

End; (* AccGEE2 *)
(************************)
procedure calccorr (var v : matrix;
                    var s : vector);
  (* put correlations into lower half of v *)
  (* v = a variance matrix *)
  (* s = sqrt(diag(v)) *)
var
  i, j : integer;
Begin (* calccorr *)
  if (v.rows <> v.cols) or (v.rows <> s.rows) then
    err('calccorr: not conformable.', HALTYE);
  with v do
  for i := 2 to rows do
    for j := 1 to pred (i) do
      if (s.ve(.i.) = 0) or (s.ve(.j.) = 0) then
        el(.i,j.) := 999  else
        el(.i,j.) := el(.i,j.)/(s.ve(.i.)*s.ve(.j.));
End;  (* calccorr *)
(************************)
procedure calcvse;
    (* calc var & se & z *)

var
  w1, w2 : matrix;
  i : integer;

Begin (* calcvse *)

  (* notice: s1 has allready been inverted in situ *)
  nvar := s1;  (* naive var est *)

  if (geelvl = 2) then
  begin
    MatMul (w1, s1, s3);
    MatMul (rvar, w1, s1);       (* robust var est *)
  end else

  begin (* GEE1 *)
    Invert (w1, s1a);
    s1a := w1;
    MatTran (w1, s1a);
    MatMul (w2, s1a, s3);
    MatMul (rvar, w2, w1);
  end;

  nse.rows := dimb;
  rse.rows := dimb;
  for i := 1 to dimb do
  begin
    nse.ve(.i.) := sqrt(nvar.el(.i,i.));
    rse.ve(.i.) := sqrt(rvar.el(.i,i.));
    if (rse.ve(.i.) > 0) then
      rz.ve(.i.) := b.ve(.i.) / rse.ve(.i.) else
      rz.ve(.i.) := 0; (* don't divide by 0 *)
  end;

  calccorr (rvar, rse);   (* correlation in lower half *)
  calccorr (nvar, nse);   (* correlation in lower half *)

End; (* calcvse *)
(************************)
procedure monitorb;

Begin  (* monitorb *)
  if monbflag then
  begin
    writeln('monitorb: iter = ':-16,iter:4, '  b = ':-6);
    VecWrite (output, b);
  end;
End; (* monitorb *)
(************************)
procedure gee; (* the gee fitting algorithm *)

var
  db : vector;  (* change in beta in one iteration *)
  lastiter : BOOLEAN;  (* true in the last iteration only *)
  c, i   : integer;

Begin (* gee *)

  MatDim (s1, dimb, dimb);
  MatDim (s1a, dimb, dimb);
  MatDim (s3, dimb, dimb);
  MatFill (s1a, 0);  (* initialize only once *)
  MatFill (s3 , 0);  (* initialize only once *)
  s2.rows := dimb;
  db.rows := dimb;
  iter := 0;
  converged := FALSE;
  lastiter := FALSE; (* one extra iteration for s3 *)

  repeat
    monitorb; (* report b if requested *)
    calclp;
    calcpsi;
    calcu;
    MatFill (s1, 0);  (* initialize for each iteration *)
    VecFill (s2, 0);
    for c := 1 to nclus do
      AccGEE2 (s1, s2, s3, s1a, lastiter, c);
       (* accumulate s1, s2, s3, s1a *)

    MatSInv (s1);
    if MatFail then
      err('gee: MatSInv(s1) failed.', HALTYE);

    MatxVec  (db, s1, s2); (* db <- s1 * s2 *)
    VecAdd (b, b, db);     (* b  <- b + db *)

    iter := succ(iter);
    converged := (VecSumAbs (db) < tol);

    if converged and (not lastiter) then
    begin
      converged := FALSE; (* go once more *)
      lastiter  := TRUE;  (* stop next time through *)
    end else
    if lastiter then converged := TRUE;   (* let go *)
         (* must have converged last time through *)
  until converged or (iter > maxiter);

  if converged then  calcvse  (* calc var & se *)

End; (* gee *)
(************************)
procedure VecWrite3 (var f : text;
                   var lbl : lbl_;
                       v1, v2, v3, v4 : vector);
var
  i : integer;

Begin (* VecWrite3 *)
  write ('Beta:                    Estimate  ');
  writeln('  Robust s.e.     Naive s.e.     Robust z');
  for i := 1 to v1.rows do
  writeln (f, lbl(.i.)   :-16    , '':2,
              v1.ve(.i.) : 13 : 5, '':2,
              v2.ve(.i.) : 13 : 5, '':2,
              v3.ve(.i.) : 13 : 5, '':2,
              v4.ve(.i.) : 13 : 5, '':2  );
End;  (* VecWrite3 *)
(************************)
procedure report; (* report estimates *)
Begin (*  report *)

  if not converged then
  begin
    writeln ('No convergence after ', iter-1:5,' iterations.');
    writeln ('Last estimate of beta:');
    VecWrite (output, b);
    RETURN;
  end;

  writeln ('Converged after ', iter-1:5,' iterations.');

  repmethods; (* report methods used *)

  writeln ('Estimates:');
  VecWrite3 (output, lbl, b, rse, nse, rz);

  writeln('(Lower = Correlation,  Upper = Variance)');
  writeln ('Robust Correlation \ Variance Matrix Estimate:');
  MatWrite (output, rvar);
  writeln ('Naive  Correlation \ Variance Matrix Estimate:');
  MatWrite (output, nvar);


  RepCPUTime;   (* report CPU time *)

End;  (*  report  *)
(************************)


BEGIN

  init;     (* setup *)
  gee;      (* fit model, if converged get var & s.e. *)
  report;   (* the output *)

END.