program prln;   { PRLN04.PAS }
{$R+}
{$N+}
{$I-}
{$D+,L+}

uses crt, dos;

type
  float = double;

const
  wmn = FALSE;    {wmn TRUE: mean income is used in calculations}
  NQuantiles = 100;         {# of quantiles, e.g. decile = 10}
  QuantD = 1/NQuantiles;   {quantile interval, e.g. decile = 0.1}
  epfx = 2;   {set epfx to:  }
    { 1 to assume income in open category}
    {   is constant and = lower bound of open category;}
    { > 1 (e.g. 1.1, 1.25, 2.0,..) to fix upper bound  }
    {   of open category at value epfx*(lower bound);  }
    { < 1 (e.g., 0) to not fix an upper bound}

type
  DisType = array [0..100] of float;   {index corresponds to lb of interval}
  QuantileType = array [0..NQuantiles] of float;
  String80 = string[80];
  BufType = array [1..32768] of Char;   {32K buffer for text files}
  PBuf = ^BufType;

   {Xlb contains LOWER BOUNDS of categories. Xlb[CatMax] is initialized to 0}
   {1970 data}  {1970 epfx = 2}
{  Xlb: distype = (0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000,
                  9000, 10000, 12000, 15000, 25000, 50000, 0);  }

   {1980 data}  {1980 epfx = 2}
{  Xlb: distype = (0, 2500, 5000, 7500, 10000, 12500, 15000, 17500, 20000, 22500,
                 25000, 27500, 30000, 35000, 40000, 50000, 75000, 0); }

   {1990 data}  {1990 epfx = 1.5}
{  Xlb: distype = (0, 5000, 10000, 12500, 15000, 17500, 20000, 22500, 25000,
                 27500, 30000, 32500, 35000, 37500, 40000, 42500, 45000,
                 47500, 50000, 55000, 60000, 75000, 100000, 125000, 150000, 0);}

var
  GiniR, GiniQ: float;
      { GiniR from original categories; GiniQ from quantiles }
  CatMax, i, j : integer;
  Dist, Xlb, p, L: DisType;
  QuantCml: QuantileType;
  ID1, ID2, NumRec: longint;
  InFile, LBFile, OutFile: text;
  InBuf, OutBuf: PBuf;

{**************  Utility procedures  ****************}
const
  IOerr: boolean = False;

const
  Null = #0;          { Ascii character codes }
  Bell = #7;
  Esc  = #27;
  Cr   = #13;

procedure Beep;
begin
  write(Bell);
end;

procedure Pause;
begin
  HighVideo;
  write('                           Press <Enter>');
  NormVideo;
  readln;
end;

procedure IOCheck;
var
  IOcode : integer;

procedure Error(Msg : String80);
begin
  Writeln;
  Beep;
  Writeln(Msg);
  Writeln;
end; { procedure Error }

begin { procedure IOCheck }
  IOcode := IOresult;
  IOerr := IOcode <> 0;
  if IOerr then
    case IOcode of
      2   : Error('File not found.');
      3   : Error('Path not found.');
      4   : Error('Too many open files.');
      5   : Error('File access denied.');
      6   : Error('Invalid file handle.');
      12  : Error('Invalid file access code.');
      15  : Error('Invalid drive number.');
      16  : Error('Cannot remove current directory.');
      17  : Error('Cannot rename across drives.');
      100 : Error('Disk read error.');
      101 : Error('Disk write error.');
      102 : Error('File not assigned.');
      103 : Error('File not open.');
      104 : Error('File not open for input.');
      105 : Error('File not open for output.');
      106 : Error('Invalid numeric format.');
      150 : Error('Disk is write-protected.');
      151 : Error('Unknown unit.');
      152 : Error('Drive not ready.');
      153 : Error('Unknown command.');
      154 : Error('CRC error in data.');
      155 : Error('Bad drive request structure length.');
      156 : Error('Disk seek error.');
      157 : Error('Unknown media type.');
      158 : Error('Sector not found.');
      159 : Error('Printer out of paper.');
      160 : Error('Device write fault.');
      161 : Error('Device read fault.');
      162 : Error('Hardware failure.');
    else
      begin
        Writeln;
        Writeln(Bell);
        Writeln('Unidentified error message = ', IOcode, '. See manual.');
        Writeln;
      end;
    end; { case }
end; { procedure IOCheck }

procedure GetInputFile(var InFile: text);
var
 filename: string;

begin
  InBuf := nil;
  writeln;
  repeat
    write('Input file name? ');
    readln(filename);
    assign(infile, filename);
    reset(infile);
    iocheck;
    if not ioerr then begin
      close(InFile);
      New (InBuf);
      reset (InFile);
      SetTextBuf (InFile, InBuf^);
    end;
  until not ioerr;

end;

procedure GetLBFile (var LBFile: text);
var
 filename: string;

begin
  writeln;
  repeat
    write('Lower Bounds file name? ');
    readln(filename);
    assign(LBFile, filename);
    reset(LBFile);
    iocheck;
    if not ioerr then begin
      close(LBFile);
      reset(LBFile);
    end;
  until not ioerr;
end;

procedure GetOutputFile(var OutFile: text);

{ This procedure determines whether output should  }
{ be sent to the screen, printer, or a disk file.  }
{ The variable OutFile is returned as the standard }
{ output channel.                                  }

var
  FileName : String;
  Ch       : char;

begin
  OutBuf := nil;
  writeln;
  write('Output to (S)creen, (F)ile, (P)rinter? ');
  Ch := ReadKey;
  if Ch = #0 then Ch := ReadKey;
  Ch := UpCase (Ch);
  writeln(Ch);

  case Ch of
             'S'  : begin
                      FileName  := 'CON';
                      Assign(OutFile, FileName);
                      Rewrite(OutFile);
                    end;

             'P'  : begin
                      FileName := 'PRN';
                      Assign(OutFile, FileName);
                      Rewrite(OutFile);
                    end;

             'F'  : begin
                      repeat
                        Ch := 'Y';
                        Writeln;
                        Write('Output file name? ');
                        Readln(FileName);
                        Assign(OutFile, FileName);
                        Reset(OutFile);
                        if IOresult = 0 then  { The file already exists. }
                        begin
                          Close(OutFile);
                          Writeln;
                          Write('File ', FileName, ' already exists. ');
                          Write('Write over it (Y/N)? ');
                          Ch := UpCase(ReadKey);
                          Writeln(Ch);
                        end;
                        if Ch = 'Y' then
                        begin
                          New(OutBuf);
                          Rewrite(OutFile);
                          SetTextBuf (OutFile, OutBuf^);
                          IOCheck;
                        end;
                      until((Ch = 'Y') and not(IOerr));
                    end;

             Esc  : Halt;     { Halt the program! }

    else begin
      FileName  := 'CON';
      Assign(OutFile, FileName);
      Rewrite(OutFile);
    end;

  end; { case }
end; { procedure GetOutputFile }

procedure GetLBs (var LBFile: text;
                  var Xlb: DisType;
                  var CatMax:integer);
var
  i: integer;

begin
  fillchar (Xlb, sizeof (Xlb), 0);
  i := 0;
  read (LBFile, Xlb [i]);
  while (not eoln(LBFile)) do begin
    inc (i);
    read (LBFile, Xlb [i]);
  end;
  close (LBFile);
  CatMax := succ (i);
  Xlb [CatMax] := 0;
end;

procedure GetDist (var ID1, ID2: longint;
                   CatMax: integer;
                   var Dist: DisType);
var
  i: integer;

begin
  ID1 := 0;
  ID2 := 0;
  fillchar (Dist, sizeof (Dist), 0);
  read (InFile, ID1, ID2);
  i := 0;
  read (InFile, Dist [i]);
  while (not eoln(InFile)) and (i<100) do begin
    inc (i);
    read (InFile, Dist [i]);
  end;
  readln (Infile);
  if (succ(i) <> CatMax) then begin
    writeln ('Incorrect number of items for ID: ', ID1:6, ID2:6);
    writeln ('Items read: ', succ (i):6);
    Pause;
    Halt;
  end;
end;


(* procedure WaitToGo;
{ Wait for the user to abort the program or continue }
const
  Esc = #27;
var
  Ch : char;
begin
  StatusLine('Esc aborts or press a key...');
  repeat until KeyPressed;
  Ch := ReadKey;
  if ch = #0 then ch := readkey;      { trap function keys }
  if Ch = Esc then
    Halt(0)                           { terminate program }
  else
    ClearDevice;                      { clear screen, go on with demo }
end; { WaitToGo } *)

function alphst (xl, nl, xu, nu: float): float;
  {estimates alpha of Pareto distribution for interval xl, xu}
begin
  alphst := (ln(nl/nu))/(ln(xu/xl));
end;

procedure GinParLin         (CatMax : integer;
                var Xlb, dist, p, L : distype;
                           var QTCml: QuantileType;
                            QuantMax: integer;
                         epfx, mean : float;
                                wmn : boolean;
                   var GiniR, GiniQ : float);
                  { GiniR from oRiginal categories, GiniQ from Quantiles }

var
  N, NCml, Xcml: distype;
  NTOT, agg, alpha,
  mdpt, EndPt, EPp, tnc, QD, QTFrac, QTX, QTN: float;
  i, jl, ju, jm, medloc, lindx, QI: integer;
  ascnd: boolean;

begin

  QD := 1/QuantMax;    {QD is quantile interval}

  { calculate N[i], the "reverse" cumulative distribution, i.e.  }
  { number of incomes GREATER than Xlb[i], so that N[CatMax-1]   }
  { is number of incomes in top category, N[0] is total number   }
  { of income receiving units  }

  N [CatMax-1] := dist[CatMax-1];
  for i := CatMax-2 downto 0 do
    N[i] := N[i+1] + dist[i];
  NTOT := N[0];      { NTOT = total number of income receiving units }

  { find index medloc of category in which median is located }

  mdpt := NTOT/2;
  ascnd := N[CatMax-1] > N[0];      { ascnd should normally be false }
  jl := -1;
  ju := CatMax;
  while ju - jl > 1 do begin
    jm := (ju+jl) div 2;
    if (mdpt > N[jm]) = ascnd then
      jl := jm
    else ju := jm
  end;
  medloc := jl;

(* writeln('medloc = ', medloc:6); readln; *)

  { calculate Ncml[i], the cumulative distribution of income,    }
  { i.e. number of incomes LESS THAN OR EQUAL TO Xlb[i], so that }
  { Ncml[0] = 0 and Ncml[CatMax] = total number of iru's NTOT    }
  { calculate Xcml[i], the cumulative sum of incomes less   }
  { than or equal to Xlb[i], so that Xcml[0] = 0 and        }
  { Xcml[CatMax] = total income of society.  The loop       }
  { estimates the aggregate income agg for each category    }

  Xcml[0] := 0;
  Ncml[0] := 0;
  p[0] := 0;
  QI := 0;   {index of quantile}

  for i := 0 to CatMax-1 do begin
    Ncml[i+1] := Ncml[i] + dist[i];
    p[i+1] := Ncml[i+1]/NTOT;
    if dist[i] = 0 then begin
      agg := 0;
      Xcml[i+1] := Xcml[i] + agg;
    end
    else begin

      { below median, estimate mean income as category midpoint }
      if i <= medloc then begin
        agg := dist[i]*(Xlb[i] + Xlb[i+1])/2;
        Xcml[i+1] := Xcml[i] + agg;
        while ((QI*QD)<p[i+1]) and not ((succ(QI)*QD)>p[i+1]) do begin
          inc(QI);
          QTFrac := (QI*QD - p[i])/(p[i+1] - p[i]);
          QTX := Xlb[i] + QTFrac*(Xlb[i+1]-Xlb[i]);
          QTCml[QI] := Xcml[i] + QTFrac*dist[i]*(Xlb[i]+QTX)/2;
        end;
      end
(* writeln('Xcml = ', i:6, Xcml[i+1]:16:0); readln; *)

      { above median, use Pareto interpolation }
      else if (i > medloc) and (i <= CatMax-2) then begin

        { if categories above are empty, cannot compute alpha, use midpoint }
        if N[i+1] = 0 then begin
          agg := dist[i]*(Xlb[i] + Xlb[i+1])/2;
          Xcml[i+1] := Xcml[i] + agg;
          while ((QI*QD)<p[i+1]) and not ((succ(QI)*QD)>p[i+1]) do begin
            inc(QI);
            QTFrac := (QI*QD - p[i])/(p[i+1] - p[i]);
            QTX := Xlb[i] + QTFrac*(Xlb[i+1]-Xlb[i]);
            QTCml[QI] := Xcml[i] + QTFrac*dist[i]*(Xlb[i]+QTX)/2;
          end;
        end
        else begin
          alpha := alphst (Xlb[i], N[i], Xlb[i+1], N[i+1]);
          if alpha > 1 then begin
            agg := (alpha/(alpha-1))*(Xlb[i]*N[i]-Xlb[i+1]*N[i+1]);
            Xcml[i+1] := Xcml[i] + agg;
            while ((QI*QD)<p[i+1]) and not ((succ(QI)*QD)>p[i+1]) do begin
              inc(QI);
              QTN := NTOT*(1-QI*QD);
              QTX := Xlb[i]*exp(ln(N[i]/QTN)/alpha);
              QTCml[QI] := Xcml[i] + (alpha/(alpha-1))
                               *(Xlb[i]*N[i]-QTX*QTN);
            end;
          end

          { if alpha <= 1 use midpoint }
          else begin
            agg := dist[i]*(Xlb[i] + Xlb[i+1])/2;
            Xcml[i+1] := Xcml[i] + agg;
            while ((QI*QD)<p[i+1]) and not ((succ(QI)*QD)>p[i+1]) do begin
              inc(QI);
              QTFrac := (QI*QD - p[i])/(p[i+1] - p[i]);
              QTX := Xlb[i] + QTFrac*(Xlb[i+1]-Xlb[i]);
              QTCml[QI] := Xcml[i] + QTFrac*dist[i]*(Xlb[i]+QTX)/2;
            end;
          end
        end
      end
      else begin

        { dealing with open category, i = CatMax-1 }
(* writeln('Xcml = ', i:6, Xcml[i+1]:16:0); readln; *)
(*        if wmn then begin      { if using mean income }
          tnc := mean*NTOT;
          if Xcml[CatMax-1] < tnc then begin
            agg := tnc - Xcml[CatMax-1];
            CatAlpha[i] := 0;
            InterPol[i] := Other;
          end
          else begin
            agg := 0;
            CatAlpha[i] := 0;
            InterPol[i] := Other;
            writeln('Inconsistency: Xcml = ', Xcml[CatMax-1]:16:0,
                    'Total Income from Mean = ', tnc:16:0);
            readln;
          end
        end
        else begin             { not using mean income }  *)

          { first calculate alpha for category just below open one }
          lindx := CatMax-2;
          alpha := alphst (Xlb[lindx], N[lindx], Xlb[CatMax-1], N[CatMax-1]);

          { if alpha less than 1, try interpolating from lower categories  }
          { until alpha > 1, but stay above category containing median     }
          while (alpha <= 1) and ((lindx-1) > medloc) do begin
            lindx := lindx - 1;
            alpha := alphst (Xlb[lindx], N[lindx], Xlb[CatMax-1], N[CatMax-1]);
          end;
          if (alpha <= 1) or (epfx = 1) then begin

            { assume income constant = lower bound of open category }
            agg := Xlb[CatMax-1]*dist[CatMax-1];
            Xcml[i+1] := Xcml[i] + agg;
            while ((QI*QD)<p[i+1]) and not ((succ(QI)*QD)>p[i+1]) do begin
              inc(QI);
              QTFrac := (QI*QD - p[i])/(p[i+1] - p[i]);
              QTX := Xlb[CatMax-1];
              QTCml[QI] := Xcml[i] + QTFrac*dist[i]*QTX;
            end;

          end
          else begin          { alpha > 1 and epfx <> 1 }
            if epfx < 1 then begin

              { do not fix upper bound in open category }
              agg := (alpha/(alpha-1))*Xlb[CatMax-1]*dist[CatMax-1];
              Xcml[i+1] := Xcml[i] + agg;
              while ((QI*QD)<p[i+1]) and not ((succ(QI)*QD)>p[i+1]) do begin
                inc(QI);
                QTN := NTOT*(1-QI*QD);
                QTX := Xlb[CatMax-1]*exp(ln(N[CatMax-1]/QTN)/alpha);
                QTCml[QI] := Xcml[i] + (alpha/(alpha-1))
                                 *(Xlb[CatMax-1]*N[CatMax-1]-QTX*QTN);
              end;
            end
            else begin

              { epfx >1, with value such as 1.1, 1.5, 2.0. etc.  }
              { want to fix upper bound of open category at      }
              { value epfx*(lower bound of open category)        }
              EndPt := Xlb[CatMax-1] * epfx;

              { estimate # above EndPt and store in N[CatMax] }
              { calculate EPp cumulative proportion at EndPt  }
              N[CatMax] := N[CatMax-1]*(exp(alpha*ln(Xlb[CatMax-1]/EndPt)));
              EPp := 1 - N[CatMax]/NTOT;

              { estimate aggregate income above lb of open category }
              agg := (alpha/(alpha-1))
                    *(Xlb[CatMax-1]*N[CatMax-1] - EndPt*N[CatMax])
                    + EndPt*N[CatMax];
              Xcml[i+1] := Xcml[i] + agg;

              { calculate quantiles below EPp }
              while ((QI*QD)<EPp) and not ((succ(QI)*QD)>EPp) do begin
                inc(QI);
                QTN := NTOT*(1-QI*QD);
                QTX := Xlb[CatMax-1]*exp(ln(N[CatMax-1]/QTN)/alpha);
                QTCml[QI] := Xcml[i] + (alpha/(alpha-1))
                                 *(Xlb[CatMax-1]*N[CatMax-1]-QTX*QTN);
              end;

              {calculate quantiles above EPp }
              while ((QI*QD)<1.0) and not ((succ(QI)*QD)>1.0) do begin
                inc(QI);
                QTFrac := (QI*QD - EPp)/(1.0 - EPp);
                QTX := EndPt;
                QTCml[QI] := Xcml[i] + (alpha/(alpha-1))
                             *(Xlb[CatMax-1]*N[CatMax-1] - EndPt*N[CatMax])
                                     + QTFrac*N[CatMax]*QTX;
              end;
            end;
          end;
(*        end;   *)
      end;
    end;
  end;

(* writeln('Xcml = ', Xcml[CatMax]:16:0); readln; *)

  { Xcml[CatMax] now contains total aggregate income. }
  { calculate cumulative population shares p[i], }
  { construct the Lorenz curve (p[i], L[i]), }
  { and calculate the Gini coefficient  GiniR }

  p[0] := 0;
  p[CatMax] := 1;
  L[0] := 0;
  L[CatMax] := 1;
  for i := 1 to pred(CatMax) do begin
    p[i] := (NTOT - N[i])/NTOT;
    L[i] := Xcml[i]/Xcml[CatMax];
  end;

  GiniR := p[CatMax];    { GiniR calculated from original intervals }
  for i := 1 to CatMax-1 do
    GiniR := GiniR - L[i]*(p[i+1]-p[i-1]);
(* writeln('GiniR = ', GiniR:12:4);  readln; *)

  { Xcml[CatMax] contains total aggregate income. }
  { Calculate cumulative quantiles QTCml }
  { and calculate the Gini coefficient  GiniQ }

  QTCml[0] := 0;
  QTCml[QuantMax] := 1;
  for i := 1 to pred(QuantMax) do
      QTCml[i] := QTCml[i]/Xcml[CatMax];

  GiniQ := 0;
  if (QuantMax>3) then begin
    GiniQ := QTCml[0]*5/12 + QTCml[1]*13/12;  { Press et al., E 4.1.12 p121 }
    for i := 2 to QuantMax-2 do               { Simpson's formula }
      GiniQ := GiniQ + QTCml[i];
    GiniQ := GiniQ + QTCml[QuantMax-1]*13/12
                   + QTCml[QuantMax]*5/12;
    GiniQ := GiniQ*QD;
    GiniQ := (0.5 - GiniQ)/0.5
  end;

end;


procedure DistOut (ID: integer; Gini:float; QuantCml: QuantileType;
                   QuantMax: integer);
const
  C: char = ',';
  Q: char = '"';

var
  i: integer;

begin
  write (OutFile, ID:8, C, Gini:8:4, C, 0.0:8:2);
  for i := 1 to QuantMax do
    write (OutFile, C, (QuantCml[i]*100):8:2);
  writeln (OutFile);
end;

(* procedure ginout (DisMax, CatMax: integer;
                  id: idtype;
                  gini: vartype);
var
  i: integer;
begin
  for i := 1 to DisMax do
    writeln (OutFile, id[i]:8, gini[i]:10:4);
  close(outfile);
end; *)

begin   {main}
  GetInputFile (InFile);
  GetLBFile (LBFile);
  GetOutputFile (OutFile);
  GetLBs (LBFile, Xlb, CatMax);
  writeln ('Number of income categories: ', CatMax:6);
  writeln ('Lower bounds of income categories:');
  for i:=0 to CatMax do write (Xlb[i]:10:0);
  Pause;
  NumRec := 0;
  repeat
    inc (NumRec);
    GetDist (ID1, ID2, CatMax, Dist);

(*    writeln ('Record: ', NumRec:6, ID1:6, ID2:6);
    for i:= 0 to CatMax-1 do write
      (Dist[i]:10:0);
    writeln; *)

    GinParLin (CatMax, Xlb, Dist, p, L, QuantCml, NQuantiles,
                 epfx, 0, wmn, GiniR, GiniQ);
    writeln (OutFile, NumRec:6, ID1:6, ID2:6, GiniR*100:8:2, GiniQ*100:8:2);

(*    writeln ('GiniR: ', GiniR:8:4, '    GiniQ: ', GiniQ:8:4);
    if NumRec mod 20 = 0 then Pause; *)

  until eof (InFile) (*   or (NumRec>=10)  *) ;

(*    gini[i] :=
    DistOut (id[i], gini[i], QuantCml, NQuantiles);
    if i mod 20 = 0 then Pause; *)

  close (InFile);
  close (OutFile);
  if InBuf<>nil then dispose(InBuf);
  if OutBuf<>nil then dispose(OutBuf);
  Pause;
end.