unit NelderMead;

interface

const
maxparameters = 5;
dSimplex = 0.5;
Cpar : array[1..maxparameters] of real = (10,-6,7,100,1);

type parameters = array[1..maxparameters] of real;


type TSimplex = class(Tobject)
       vertices: array[0..maxparameters] of parameters;
       values: array[0..maxparameters] of real;
       Margin:real;
       IHigh:integer;
       INearhigh:integer;
       ILow :integer;
       Centre: parameters;
       CentreValue:real;
       Reflection: parameters;
       ReflValue:real;
       Expansion:parameters;
       ExpaValue:real;
       Contraction:parameters;
       ContrValue:real;
       alpha:real;
       beta:real;
       gamma:real;
       delta:real;
       maxpar:integer;
       constructor New(a,b,c,d:real; npar:integer);
       procedure StartSimplex;
       procedure CalculateCentre;
       procedure CalculateReflection;
       procedure CalculateExpansion;
       procedure CalculateContraction;
       function CalculateValue(P:Parameters):real; virtual;
       procedure CalculateMargin;
       procedure DetermineHighAndLow;
       procedure Replace(n:integer;P:parameters; v:real);
       procedure Shrink;
       function Iteration:real;
       end;

type TBrassSimplex = class(TSimplex)
       function CalculateValue(P:Parameters):real; override;
       end;



implementation

uses unit2;


procedure TSimplex.Replace(n: Integer; P: parameters;v:real);
var i:integer;
begin
  for i:= 1 to maxpar do Vertices[n,i] := P[i];
  values[n] := v;
end;

constructor TSimplex.New(a,b,c,d: Real; npar:integer);
begin
  inherited create;
  alpha := a;
  beta := b;
  gamma := c;
  delta := d;
  maxpar := npar;
end;

function TSimplex.CalculateValue(P:Parameters):real;
var r:real;
    i:integer;
begin
  with form2
  do begin
     For i:= 1 to maxpar
     do if i = 1 then par[i] := P[i]/100 else par[i] := P[i];
     CalculateParameters;
     BuildLifeTable;
     r := fit;
     end;
  result := r;
end;

function TBrassSimplex.CalculateValue(P:Parameters):real;
var r, rFit, rtot, rdist :real;
    i:integer;
begin
  rFit := 0;
  rtot := 0;
  for i:= 1 to PrincetonNumber -1
  do if form2.logitArray1[i] <> unknown
     then begin
          r := SQR(P[1] + P[2] * form2.LogitArray2[i] - form2.LogitArray1[i]);
          rdist := form2.Lifetable[LTageDist,i];
          rFit := rFit + r * rdist ;
          rtot := rtot + rdist;
          end;
  rFit := SQRT(rFit/rtot);
  result := rFit;
end;

procedure TSimplex.StartSimplex;
var i,j:integer;
    p, q:real;
begin
   p := dSimplex * (sqrt(maxpar+1) + maxpar - 1)/(maxpar * sqrt(2));
   q := dSimplex * (sqrt(maxpar+1)      -1)/(maxpar * sqrt(2));
   for j:= 1 to maxpar do Vertices[0,j] := Cpar[j] - (p + (maxpar-1) * q)/(maxpar+1);
   for i:= 1 to maxpar
   do for j:= 1 to maxpar
      do if j=i
         then Vertices[i,j] := Vertices[0,j] + p
         else Vertices[i,j] := Vertices[0,j] + q;
   for i:= 0 to maxpar do values[i] := CalculateValue(Vertices[i]);
   DetermineHighAndLow;
   CalculateCentre;
end;

procedure TSimplex.CalculateMargin;
var i,j:integer;
    r:real;
begin
  Margin := 0;
  for i:= 0 to maxpar
  do begin
     r := 0;
     for j:= 1 to maxpar do r := r + sqr(Vertices[i,j] - Centre[j]);
     if sqrt(r) > Margin then Margin := sqrt(r);
     end;
end;

procedure TSimplex.CalculateCentre;
var i,j:integer;
begin
  for i:= 1 to maxpar
  do begin
     Centre[i] := 0;
     for j:= 0 to maxpar do if j <> Ihigh then Centre[i] := Centre[i] + Vertices[j,i]/maxpar;
     end;
  CentreValue := CalculateValue(Centre);
  CalculateMargin;
end;

procedure TSimplex.CalculateReflection;
var i:integer;
begin
  for i:= 1 to maxpar
  do Reflection[i] := (1 + alpha) * Centre[i] - alpha * Vertices[Ihigh,i];
  ReflValue := CalculateValue(Reflection);
end;

procedure TSimplex.CalculateExpansion;
var i:integer;
begin
  for i:= 1 to maxpar
  do Expansion[i] := gamma * Reflection[i] + (1-gamma) * Centre[i];
  ExpaValue := CalculateValue(Expansion);
end;

procedure TSimplex.CalculateContraction;
var i:integer;
begin
  for i:= 1 to maxpar
  do Contraction[i] := beta * Vertices[IHigh,i] + (1-beta) * Centre[i];
  ContrValue := CalculateValue(Contraction);
end;

procedure TSimplex.Shrink;
var i,j:integer;
begin
  for i:= 0 to maxpar
  do if i <> ILow
     then begin
          for j:= 1 to maxpar
          do Vertices[i,j] := delta * Vertices[i,j] + (1-delta) * Vertices[ILow,j];
          values[i] := CalculateValue(Vertices[i]);
          end;
end;

procedure TSimplex.DetermineHighAndLow;
var r:real;
    i:integer;
begin
  r := values[0];
  ILow := 0;
  for i:= 1 to maxpar
  do if values[i] < r
     then begin
          Ilow := i;
          r := values[i];
          end;
  r := values[0];
  IHigh := 0;
  for i:= 1 to maxpar
  do if values[i] > r
     then begin
          IHigh := i;
          r := values[i];
          end;
  r := values[ILow];
  INearHigh := ILow;
  for i:= 0 to maxpar
  do if ((values[i] > r) and  (i<> IHigh))
     then begin
          INearHigh := i;
          r := values[i];
          end;
end;


function TSimplex.Iteration:real;
begin
  CalculateReflection;
  If ReflValue < values[Ilow]
  then begin
       CalculateExpansion;
       If ExpaValue < values[Ilow]
       then Replace(IHigh,Expansion, ExpaValue)
       else Replace(IHigh,Reflection, ReflValue);
       end
  else if NOT(ReflValue > values[INearHigh])
       then Replace(IHigh,Reflection, ReflValue)
       else begin
            if NOT(ReflValue > values[IHigh])
            then Replace(IHigh,Reflection, ReflValue);
            CalculateContraction;
            if Not(ContrValue > values[IHigh])
            then Replace(IHigh,Contraction, ContrValue)
            else Shrink;
            end;
  DetermineHighAndLow;
  CalculateCentre;
  result := CentreValue;
end;




end.





unit Unit2;

interface

uses
  Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
  Vcl.Controls, Vcl.Forms, Vcl.Dialogs, Vcl.Menus, Vcl.StdCtrls, Vcl.ExtCtrls,
  Vcl.Grids, Vcl.Samples.Spin, Vcl.ComCtrls, NelderMead;

const
unknown = -9999;
tab = #9;
maxNumber = 1000000000000000;
IPPrinceton = 1;
IPFileSpaces = 2;
IPFileTabs = 3;
AppSiler = 1;
App4Par = 10;
App3Par = 11;
App2Par = 12;
AppBrassLogitL = 13;
AppBrassLogitQ = 15;
AimFit = 1;
AimCompare = 2;
FMwrmse = 1;
FMLikeliHood = 2;
LTAge = 0;
LTq = 1;
LTl = 2;
LTLL = 3;
LTm = 4;
LTT = 5;
LTE = 6;
LTAgeDist = 7;
LTDeathDist = 8;
maxAge = 120;
HMDAgegroupColumn = -4;
HMDAgeColumn = -3;
Princetoncolumn = -2;
PrincetonNumber = 22;

PrincetonAges: array[1..22] of integer = (0,1,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95,100);
HMDages: array[1..24] of integer = (0,1,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95,100,105,110);


 mortWestMales : array[1..26,1..22] of real = (
 (419.50,260.00,67.59,48.48,65.15,92.35,103.64,119.73,139.33,165.99,188.02,231.10,270.76,351.92,440.37,557.86,708.62,801.55,906.40,1000.00,1000.00,1000.00),
 (383.86,236.18,61.64,44.25,59.64,84.54, 94.77,109.44,127.47,152.16,173.13,213.72,252.18,329.36,414.95,529.54,677.50,781.41,893.80, 963.38,1000.00,1000.00),
 (351.74,214.72,56.28,40.43,54.67,77.50, 86.77,100.17,116.79,139.70,159.72,198.06,235.44,309.05,392.05,504.02,649.46,762.27,881.41, 957.75,1000.00,1000.00),
 (322.57,195.23,51.41,36.97,50.17,71.10, 79.51, 91.75,107.09,128.38,147.54,183.83,220.24,290.59,371.25,480.85,623.98,744.08,869.24, 952.01,1000.00,1000.00),
 (295.89,177.41,46.96,33.80,46.04,65.25, 72.87, 84.04, 98.21,118.03,136.40,170.82,206.33,273.71,352.22,459.65,600.69,726.75,857.32, 946.19,1000.00,1000.00),
 (271.35,161.01,42.86,30.88,42.25,59.87, 66.76, 76.96, 90.05,108.50,126.15,158.85,193.54,258.18,334.72,440.15,579.26,710.24,845.64, 940.30,1000.00,1000.00),
 (248.65,145.85,39.07,28.18,38.74,54.90, 61.11, 70.40, 82.50, 99.70,116.67,147.78,181.71,243.82,318.53,422.12,559.44,694.46,834.22, 934.37,1000.00,1000.00),
 (227.57,131.76,35.55,25.68,35.48,50.28, 55.86, 64.32, 75.49, 91.52,107.87,137.50,170.72,230.49,303.50,405.38,541.03,679.38,823.05, 928.41,1000.00,1000.00),
 (207.91,118.63,32.27,23.34,32.44,45.97, 50.97, 58.64, 68.95, 83.89, 99.66,127.92,160.48,218.05,289.48,389.76,523.87,664.95,812.13, 922.43,1000.00,1000.00),
 (189.51,106.34,29.20,21.16,29.60,41.93, 46.39, 53.33, 62.83, 76.75, 91.98,118.95,150.89,206.41,276.36,375.14,507.80,651.12,801.46, 916.45, 977.88,1000.00),
 (172.25, 94.80,26.32,19.11,26.93,38.15, 42.09, 48.34, 57.09, 70.05, 84.77,110.53,141.89,195.49,264.05,361.43,492.72,637.84,791.04, 910.48, 975.76,1000.00),
 (155.99, 83.94,23.61,17.18,24.42,34.59, 38.04, 43.65, 51.68, 63.74, 77.98,102.60,133.42,185.20,252.46,348.51,478.53,626.08,780.85, 904.51, 973.60,1000.00),
 (140.17, 71.40,20.71,15.00,22.03,31.25, 34.25, 39.26, 46.62, 57.79, 71.48, 94.72,125.03,175.14,241.17,336.18,464.93,612.63,770.73, 898.48, 971.35,1000.00),
 (125.13, 60.49,18.17,13.26,19.75,28.05, 30.61, 35.00, 41.83, 52.24, 65.78, 88.33,118.38,167.19,232.36,326.47,454.40,602.83,762.67, 893.58, 969.49,1000.00),
 (111.96, 51.45,15.99,11.75,17.74,25.16, 27.31, 31.18, 37.46, 47.30, 60.62, 82.49,112.25,159.82,224.14,317.37,444.50,593.49,754.89, 888.78, 967.63,1000.00),
 ( 99.16, 43.08,13.89,10.28,15.77,22.35, 24.11, 27.49, 33.21, 42.44, 55.48, 76.60,106.00,152.26,215.67,307.96,434.22,583.66,746.60, 883.60, 965.59,1000.00),
 ( 86.78, 35.31,11.88, 8.87,13.87,19.62, 21.01, 23.92, 29.08, 37.69, 50.37, 70.69, 99.67,144.57,207.00,298.27,423.60,573.37,737.82, 878.02, 963.35,1000.00),
 ( 74.83, 28.07, 9.96, 7.52,12.02,16.99, 18.03, 20.49, 25.10, 33.05, 45.32, 64.80, 93.29,136.78,198.17,288.36,412.70,562.65,728.56, 872.05, 960.91,1000.00),
 ( 63.34, 21.29, 8.14, 6.21,10.24,14.45, 15.18, 17.19, 21.27, 28.54, 40.37, 58.96, 86.91,128.95,189.25,278.30,401.59,551.58,718.87, 865.70, 958.27,1000.00),
 ( 52.33, 14.91, 6.39, 4.97, 8.53,12.01, 12.44, 14.04, 17.59, 24.18, 35.52, 53.20, 80.58,121.14,180.29,268.15,390.35,540.21,708.79, 858.99, 955.41,1000.00),
 ( 41.34,  9.70, 4.81, 3.85, 6.88, 9.63,  9.85, 11.15, 14.17, 20.24, 30.99, 47.81, 74.27,113.24,171.03,257.54,378.54,528.10,697.91, 851.62, 952.22,1000.00),
 ( 30.99,  6.18, 3.45, 2.84, 5.27, 7.34,  7.39,  8.35, 10.77, 15.91, 25.59, 40.91, 65.95,102.52,158.13,242.17,361.00,509.78,681.16, 840.03, 947.06,1000.00),
 ( 21.62,  3.47, 2.26, 1.92, 3.74, 5.18,  5.11,  5.75,  7.56, 11.69, 20.02, 33.50, 56.64, 90.23,142.99,223.92,339.84,487.14,659.99, 824.97, 940.12,1000.00),
 ( 13.48,  1.63, 1.29, 1.15, 2.39, 3.29,  3.15,  3.53,  4.76,  7.80, 14.51, 25.78, 46.39, 76.33,125.31,202.06,313.96,458.66,632.60, 804.79, 930.43,1000.00),
 (  7.11,  0.58, 0.61, 0.57, 1.30, 1.77,  1.63,  1.83,  2.55,  4.51,  9.39, 18.08, 35.40, 60.86,104.82,175.83,282.05,422.30,596.40, 776.95, 916.34,1000.00),



  mortWestFemales : array[1..26,1..22] of real = (
 (365.55,261.78,73.18,57.22,74.00,91.90,102.62,115.61,125.88,133.36,140.78,178.86,221.87,313.59,396.53,528.84,669.53,778.31,893.57,964.25,1000.00,1000.00),
 (333.99,237.60,66.57,52.05,67.44,83.85, 93.69,105.58,115.11,122.27,129.67,165.18,205.71,291.44,371.18,498.58,637.20,756.01,879.19,957.85,1000.00,1000.00),
 (305.56,215.82,60.61,47.38,61.53,76.60, 85.65, 96.54,105.41,112.27,119.67,152.85,191.16,271.49,348.35,471.31,608.08,734.85,865.02,951.26,1000.00,1000.00),
 (279.73,196.04,55.19,43.14,56.16,70.02, 78.34, 88.34, 96.60,103.19,110.58,141.66,177.94,253.36,327.61,446.55,581.63,714.73,851.11,944.52,1000.00,1000.00),
 (256.11,177.94,50.24,39.27,51.25,63.99, 71.66, 80.83, 88.54, 94.89,102.26,131.42,165.85,236.78,308.64,423.90,557.44,695.59,837.46,937.67,1000.00,1000.00),
 (234.38,161.30,45.69,35.70,46.73,58.45, 65.51, 73.93, 81.12, 87.25, 94.62,122.00,154.72,221.53,291.19,403.06,535.18,677.35,824.08,930.72,1000.00,1000.00),
 (214.29,145.90,41.48,32.40,42.56,53.33, 59.82, 67.54, 74.26, 80.18, 87.54,113.29,144.43,207.42,275.05,383.80,514.60,659.93,810.99,923.70,1000.00,1000.00),
 (195.62,131.61,37.56,29.34,38.68,48.57, 54.54, 61.61, 67.89, 73.62, 80.98,105.20,134.88,194.32,260.06,365.90,495.48,643.29,798.18,916.64, 978.87,1000.00),
 (178.22,118.27,33.92,26.48,35.06,44.13, 49.62, 56.08, 61.95, 67.50, 74.85, 97.65,125.97,182.11,246.08,349.21,477.65,627.36,785.64,909.54, 976.45,1000.00),
 (161.93,105.80,30.50,23.81,31.68,39.98, 45.01, 50.91, 56.40, 61.78, 69.12, 90.59,117.63,170.67,233.01,333.60,460.97,612.09,773.38,902.42, 973.95,1000.00),
 (146.64, 94.09,27.30,21.30,28.50,36.08, 40.68, 46.05, 51.18, 56.40, 63.74, 83.96,109.80,159.94,220.73,318.94,445.31,597.44,761.39,895.29, 971.38,1000.00),
 (132.25, 83.06,24.28,18.94,25.51,32.41, 36.61, 41.48, 46.27, 51.34, 58.67, 77.72,102.44,149.84,209.17,305.14,430.57,583.37,749.66,888.16, 968.73,1000.00),
 (118.79, 72.32,21.46,16.71,22.71,28.98, 32.80, 37.20, 41.67, 46.61, 53.93, 71.89, 95.54,140.39,198.36,292.23,416.78,569.95,738.29,881.11, 966.05,1000.00),
 (106.04, 60.22,18.52,14.29,20.25,26.28, 29.83, 33.66, 37.75, 42.50, 49.90, 66.89, 89.77,132.44,189.68,282.03,405.73,559.03,728.90,875.17, 963.74,1000.00),
 ( 93.94, 50.44,15.84,12.27,17.52,22.85, 26.07, 29.59, 33.56, 38.38, 45.93, 62.12, 84.21,124.95,181.22,272.06,395.30,548.57,719.79,869.32, 961.42,1000.00),
 ( 82.31, 41.57,13.37,10.38,14.96,19.63, 22.53, 25.73, 29.52, 34.35, 41.98, 57.34, 78.60,117.35,172.55,261.81,384.50,537.60,710.11,863.01, 958.88,1000.00),
 ( 71.16, 33.51,11.07, 8.62,12.58,16.63, 19.21, 22.07, 25.65, 30.42, 38.07, 52.56, 72.96,109.68,163.74,251.31,373.37,526.14,699.86,856.21, 956.08,1000.00),
 ( 60.51, 26.11, 8.93, 6.97,10.36,13.82, 16.08, 18.60, 21.94, 26.60, 34.22, 47.83, 67.33,101.98,154.84,240.64,361.98,514.23,689.08,848.95, 953.03,1000.00),
 ( 50.35, 19.28, 6.94, 5.42, 8.28,11.20, 13.15, 15.33, 18.41, 22.92, 30.45, 43.16, 61.75, 94.31,145.89,229.85,350.40,501.95,677.80,841.22, 949.71,1000.00),
 ( 40.69, 12.91, 5.07, 3.97, 6.33, 8.74, 10.39, 12.23, 15.04, 19.37, 26.77, 38.59, 56.24, 86.73,136.97,219.03,338.71,489.38,666.09,833.05, 946.14,1000.00),
 ( 31.16,  7.81, 3.42, 2.75, 4.42, 6.34,  7.69,  9.18, 11.70, 15.88, 23.07, 34.04, 50.61, 79.10,127.75,207.77,326.47,476.02,653.47,824.09, 942.13,1000.00),
 ( 22.82,  4.69, 2.23, 1.82, 2.96, 4.16,  5.16,  6.45,  8.63, 12.33, 18.92, 28.63, 43.68, 69.29,115.24,191.66,308.18,455.69,633.89,809.87, 935.59,1000.00),
 ( 15.30,  2.44, 1.29, 1.07, 1.79, 2.56,  3.25,  4.18,  5.86,  8.92, 14.67, 22.92, 36.16, 58.47,101.07,173.11,286.67,431.17,609.72,791.77, 926.96,1000.00),
 (  9.05,  1.04, 0.63, 0.53, 0.92, 1.35,  1.77,  2.36,  3.53,  5.83, 10.51, 17.13, 28.22, 46.80, 85.10,151.48,260.73,400.74,578.79,767.76, 915.01,1000.00),
 (  4.45,  0.32, 0.24, 0.21, 0.38, 0.57,  0.78,  1.09,  1.78,  3.28,  6.70, 11.55, 20.18, 34.63, 67.43,126.45,229.32,362.56,538.55,735.07, 897.85,1000.00),


type Pattern = array[1..maxage] of real;

type
  TForm2 = class(TForm)
    MainMenu1: TMainMenu;
    Panel1: TPanel;
    Label1: TLabel;
    Label2: TLabel;
    Label3: TLabel;
    Label4: TLabel;
    Label5: TLabel;
    Label6: TLabel;
    Label7: TLabel;
    Label8: TLabel;
    Edit1: TEdit;
    Edit2: TEdit;
    Edit3: TEdit;
    Edit4: TEdit;
    Edit5: TEdit;
    Hazardfunction1: TMenuItem;
    Siler1: TMenuItem;
    Edit6: TEdit;
    Panel2: TPanel;
    Label9: TLabel;
    ComboBox1: TComboBox;
    Label10: TLabel;
    StringGrid1: TStringGrid;
    OpenDialog1: TOpenDialog;
    file1: TMenuItem;
    OpenLifeTabledatafilewithspaties1: TMenuItem;
    ComboBox2: TComboBox;
    Label12: TLabel;
    SpinEdit2: TSpinEdit;
    ComboBox3: TComboBox;
    Panel3: TPanel;
    Memo1: TMemo;
    Memo2: TMemo;
    Princetonmodellifetable1: TMenuItem;
    Panel4: TPanel;
    ComboBox4: TComboBox;
    Label11: TLabel;
    Label13: TLabel;
    SpinEdit1: TSpinEdit;
    comparefit1: TMenuItem;
    Compare4parametermodel: TMenuItem;
    Label14: TLabel;
    Compare3paramterModel: TMenuItem;
    Compare2Parametermodel: TMenuItem;
    Panel5: TPanel;
    Model4Par: TMenuItem;
    Model3Par: TMenuItem;
    Model2Par: TMenuItem;
    BrasslogitSurvivor: TMenuItem;
    BraaslogitQprocedure1: TMenuItem;
    Panel6: TPanel;
    Label15: TLabel;
    SpinEdit3: TSpinEdit;
    CompareSiler1: TMenuItem;
    Button4: TButton;
    Button5: TButton;
    Label16: TLabel;
    procedure Siler1Click(Sender: TObject);
    procedure FormCreate(Sender: TObject);
    procedure OpenLifeTabledatafilewithspaties1Click(Sender: TObject);
    procedure SpinEdit2Change(Sender: TObject);
    procedure Princetonmodellifetable1Click(Sender: TObject);
    procedure Compare2parametermodel1Click(Sender: TObject);
    procedure Compare3paramterModelClick(Sender: TObject);
    procedure Compare4parametermodelClick(Sender: TObject);
    procedure Model4ParClick(Sender: TObject);
    procedure Model3ParClick(Sender: TObject);
    procedure Model2ParClick(Sender: TObject);
    procedure BrasslogitSurvivorClick(Sender: TObject);
    procedure BraaslogitQprocedure1Click(Sender: TObject);
    procedure CompareSiler1Click(Sender: TObject);
    procedure Memo1DblClick(Sender: TObject);
    procedure Button4Click(Sender: TObject);
    procedure Button5Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
    // data
    Inputtype:integer;
    ageColumn:integer;
    UsedColumn:integer;
    UsedVariable:integer;
    modelvisible:boolean;
    DataArray: Pattern;
    // model
    Aim:integer;
    ApplicationNumber :integer;
    FittingMeasure:integer;
    country:string;
    ready:boolean;
    ParNumber:integer;
    par:parameters;
    gamma:real;
    ksiPar:real;
    growthrate:real;
    growthStartAge:real;
    generationGap:real;
    bestFit:real;
    Margin:real;
    //lifetable
    LifeTable :array[0..8] of Pattern;
    lifeTableSize :integer;
    TableSize:integer;
    CDR:real;
    TotalPopulation:real;
    TotalDeaths:real;
    iteration:integer;
    //logit
    logitArray1:Pattern;
    logitArray2:Pattern;
    // log
    logfile:string;
    logdir:string;
    // data
    procedure ClearStringgrid1;
    procedure FillLifeTableWithAges;
    procedure FillStringgridwithSpaceString(s:string;var l:integer);
    // model
    procedure InitieerParameters;
    procedure CalculateParameters;
    function Survivorship(x: Real):real;
    function Schrijfreal(r:real;decimalen:integer):string;

    //lifetable
    procedure FillDataArray;
    procedure FillDataArrayWithPrincetonModel;
    procedure FillDataArrayWithStringgrid;
    procedure BuildLifeTable;
    procedure DrawLifeTable;
    {procedure SetNewStep(var rOld, rBest,Rstep:real; Rmin:real);}
    procedure WriteLifeTable;
    procedure WriteBestFit(s:string);
    procedure Compare;
    procedure SetIterationOptions(n:integer);
    procedure FindBestFitNelderMead;
    function fit:real;
    function fitWRMSE(n:integer;Pmodel,Pdata,Pweights:Pattern):real;
    function fitMaxLikeliHood:real;
    function LComb(n,k:integer):real;
    function LBinom(n,k:integer; q:real):real;

    //logit procedure
    function LogitSurvivalProcedure:string;
    function LogitQprocedure:string;
  end;


var
  Form2: TForm2;
  BAbort:boolean;

implementation

{$R *.dfm}

uses unit1;





function Tform2.LogitSurvivalProcedure:string;
var i, ia, ib, code:integer;
    y1, y2, aBest, bBest, rfit, rtot, rdist, r, rr, a, b:real;
    s, ss:string;
begin
  //BuildLogit
  for i:= 1 to 23
  do begin
     s := Stringgrid1.Cells[6,i];
     val(s,rr,code);
     if ((rr > 0) and (rr < 100000)) then begin
                    y1 := rr/100000;
                    LogitArray1[i] := 0.5 * Ln((1-y1)/y1);
                    case ApplicationNumber of
                     AppBrassLogitL: y2 := Eng1841[i]/100000;
                     end;
                    logitArray2[i] := 0.5 * Ln((1-y2)/y2);
                    end
     else logitArray1[i] := unknown;
     end;
  //Find Best fit logit
  Bestfit := maxNumber;
  for ia := -250 to 250
  do for ib := -250 to 250
     do begin
        rfit := 0;
        rtot := 0;
        a := ia * 0.01;
        b := ib * 0.01;
        for i:= 1 to 23
        do begin
           if logitArray1[i] <> unknown
           then begin
                r := SQR(a + b * LogitArray2[i] - LogitArray1[i]);
                rdist := Lifetable[LTageDist,i];
                rFit := rFit + r * rdist ;
                rtot := rtot + rdist;
                end;
           end;
        rFit := SQRT(rFit/rtot);
        if rFit < Bestfit
        then begin
             aBest := a;
             bBest := b;
             Bestfit := rfit;
             end;
        end;
  rtot := 0;
  rfit := 0;
  for i:= 1 to 23
  do  begin
      s := Stringgrid1.Cells[6,i];
      val(s,rr,code);
      if ((rr > 0) and (rr < 100000)) then begin
                     y1 := rr/100000;
                     y2 := 1/(1+exp(2*(aBest + bBest *logitArray2[i])));
                     r := SQR(y2 - y1);
                     rdist := Lifetable[LTageDist,i];
                     rFit := rFit + r * rdist ;
                     rtot := rtot + rdist;
                     end;
      end;
  rFit := SQRT(rFit/rtot);
  s := 'logit' + tab + schrijfreal(aBest,3) + tab + schrijfreal(bBest,3) + tab + schrijfreal(rfit,5);
  result := s;
end;



function Tform2.LogitQProcedure:string;
var i, ia, ib, k:integer;
    y1, y2, rfit, x, xx, xy, VarXX, VarXY, W, Wtot, Xmean, a, b:real;
    ss, sss:string;
    Pmodel, Pdata:Pattern;
    S:TBrassSimplex;
begin
  BestFit := maxNumber;
  S := TBrassSimplex.New(1,0.5,2,0.5,ParNumber);
  S.StartSimplex;

  //Build Logit for the Actual and the Reference LifeTable
  k := spinedit3.Value;
  for i:= 1 to PrincetonNumber -1
  do begin
     y1 :=  DataArray[i];
     if ((y1 > 0) and (y1 < 1)) then logitArray1[i] := 0.5 * Ln((1-y1)/y1) else logitArray1[i] := unknown;
     y2 := 0.5 * (mortWestfemales[k,i]/1000 + mortWestmales[k,i]/1000);
     If ((y2 > 0) and (y2 < 1)) then logitArray2[i] := 0.5 * Ln((1-y2)/y2)
                                else begin logitArray2[i] := unknown; logitArray1[i] := unknown; end;
     end;

  //parameter fit
  iteration:= 0;
  While (( iteration <5000) and (S.Margin > dSimplex/10000))
  do begin
     S.Iteration;
     inc(iteration,1);
     caption := 'iteration  ' + inttostr(iteration);
     end;
  BestFit := S.CalculateValue(S.vertices[S.ILow]);
  a := S.vertices[S.Ilow,1];
  b := S.vertices[S.Ilow,2];
  Margin := S.Margin;


    //present the best fitting life table
  sss := '';
  for i:= 1 to PrincetonNumber -1
  do begin
     if logitArray2[i] <> unknown
     then Pmodel[i] := 1/(1+exp(2*(a + b *logitArray2[i])))
     else Pmodel[i] := unknown;
     if logitArray1[i] <> unknown
     then Pdata[i] := DataArray[i]
     else Pdata[i] := unknown;
     sss := sss + tab + inttostr(i) + tab + schrijfreal(Pmodel[i],3);
     end;
  rfit := fitWRMSE(PrincetonNumber -1, Pmodel, Pdata, Lifetable[LTageDist]);
  ss := schrijfreal(a,3) + tab + schrijfreal(b,3) + tab + schrijfreal(rfit,5) + tab + schrijfreal(S.Margin,5) + tab + inttostr(iteration) {+ sss};
  S.Free; S:= nil;
  result := ss;
end;



function TForm2.Schrijfreal(r:real;decimalen:integer):string;
var s:string;
begin
 str(r:2:decimalen,s);
 result := s;
end;

procedure TForm2.ClearStringgrid1;
var i,j:integer;
begin
for i:= 0 to stringgrid1.colCount do stringgrid1.Cells[i,0] := inttostr(i);
for i:= 0 to stringgrid1.colcount
do for j:= 1 to stringgrid1.rowcount
   do stringgrid1.cells[i,j] := '';
end;


procedure TForm2.FormCreate(Sender: TObject);
var F:textfile;
    s:string;
    i:integer;
begin
  BAbort := false;
  stringgrid1.ColCount := 10;
  Stringgrid1.RowCount := 120;
  country := '';
  Aim := unknown;
  ApplicationNumber := unknown;
  lifeTableSize := 0;
  BestFit := maxNumber;
  stringgrid1.ShowHint := true;
  modelvisible := true;
  s := DateToStr(Date) + ' ' + Timetostr(Time);
  for i:= 1 to length(s) do if s[i] in ['/',':'] then s[i] := '_';
  logdir :=  getcurrentdir + '\logfiles\';
  logfile := logdir + 'log' + s +  '.txt';
  assignfile(F,logfile);
  rewrite(F);
  closefile(F);
end;


procedure TForm2.SpinEdit2Change(Sender: TObject);
begin
  If Spinedit2.Value < 0
  then spinedit2.Value := 0;
  If Spinedit2.Value > Stringgrid1.colcount
  then Spinedit2.value := stringgrid1.colcount;
end;

procedure TForm2.InitieerParameters;
var code:integer;
begin
  val(edit1.Text,Par[1],code);
  val(edit2.Text,Par[2],code);
  val(edit3.Text,Par[3],code);
  val(edit4.Text,Par[4],code);
  val(edit5.Text,Par[5],code);
  val(edit6.Text,growthrate,code);
  ksiPar := 1;
  GrowthStartAge :=0;
  GenerationGap := 0;
  If combobox1.itemindex = 1 then begin GrowthStartAge := 20; GenerationGap := 30; end;
  If Inputtype = IPPrinceton
  then begin
       UsedVariable := LTq;
       AgeColumn := PrincetonColumn;
       end
  else begin
       case combobox2.itemindex of
           0 : UsedVariable := LTq;
           1 : UsedVariable := LTl;
           2 : UsedVariable := LTm;
           3 : UsedVariable := LTE;
           4 : UsedVariable := LTAgeDist;
           5 : UsedVariable := LTDeathDist;
           else UsedVariable := unknown;
          end;
       UsedColumn := spinedit2.value;
       case combobox3.ItemIndex of
         0: AgeColumn := HMDAgegroupColumn;
         1: AgeColumn := HMDAgeColumn;
         2: AgeColumn := PrincetonColumn;
         3: AgeColumn := 0;
         4: AgeColumn := 1;
         5: AgeColumn := 2;
         6: AgeColumn := 3;
         else AgeColumn := unknown;
        end;
       end;
end;



procedure TForm2.FillStringgridWithSpaceString(s:string;var l:integer);
var i, n:integer;
    ss:string;
    b:boolean;
begin
  If s <> ''
  then begin
       inc(l,1);
       n:= 0;
       b := false;
       for i:= 1 to length(s)
       do begin
          if s[i] = ' '
          then b := true
          else begin
               if b
               then begin
                    Stringgrid1.cells[n,l] := ss;
                    inc(n,1);
                    ss := s[i];
                    b := false;
                    end
               else ss := ss + s[i];
               end;
          end;
       Stringgrid1.Cells[n,l] := ss;
       ss := '';
       end;
end;




procedure TForm2.OpenLifeTabledatafilewithspaties1Click(Sender: TObject);
var F:textfile;
    s:string;
    l:integer;
begin
  Panel4.Visible := false;
  Panel2.Visible := true;
  InputType := IPFileSpaces;
  opendialog1.filter := 'Textfiles (*txt)|*TXT';
  if opendialog1.execute
  then begin
       assignfile(F,opendialog1.filename);
       reset(F);
       ClearStringgrid1;
       l:= 0;
       while NOT(EOF(F))
       do begin
          readln(F,s);
          FillStringgridWithSpaceString(s,l);
          end;
       closefile(F);
       end;
end;

procedure TForm2.Princetonmodellifetable1Click(Sender: TObject);
begin
 Panel2.Visible := false;
 panel4.Visible := true;
 InputType := IPPrinceton;
end;

procedure Tform2.fillDataArrayWithPrincetonModel;
var k,i:integer;
begin
  k := spinedit1.Value;
  case combobox4.ItemIndex of
    0 : for i:= 1 to PrincetonNumber do DataArray[i] := mortWestfemales[k,i]/1000;
    1 : for i:= 1 to PrincetonNumber do DataArray[i] := mortWestmales[k,i]/1000;
    2 : for i:= 1 to PrincetonNumber do DataArray[i] := (mortWestfemales[k,i] + mortWestmales[k,i])/2000;
   end;
  UsedVariable := LTq;
  DataArray[LifeTableSize] := 1;
end;

procedure Tform2.FillDataArrayWithStringgrid;
var r, rr:real;
    i,n, code:integer;
    s:string;
begin
   //Fill Data Array with data
  If UsedColumn = unknown then BAbort := true;
  If UsedVariable = LTl then rr := 100000 else rr := 1;
  If Not(BAbort)
  then begin
       n:= 0;
       for i:= 1 to stringgrid1.RowCount
       do begin
          s := Stringgrid1.Cells[UsedColumn,i];
          If s <> ''
          then begin
               if s = '?'
               then begin
                    inc(n,1);
                    DataArray[n] := unknown;
                    end
               else begin
                    val(s, r, code);
                    if code = 0
                    then begin
                         inc(n,1);
                         DataArray[n] := r/rr;
                         end
                    else BAbort := true;
                    end;
               end;
          end;
       if LifeTableSize > n+1 then LifeTablesize := n+1;
       if UsedVariable = LTq then DataArray[LifeTableSize] := 1 else DataArray[LifeTableSize] := 0;
       end;
end;

procedure Tform2.FillLifeTableWithAges;
var i,n, code:integer;
    r:real;
    s:string;
begin
   //Fill LifeTable with Ages
  case AgeColumn of
   HMDAgegroupColumn: for n:= 1 to 24 do LifeTable[LTAge,n] := HMDAges[n];
   HMDAgeColumn : for n:= 1 to 111 do LifeTable[LTAge,n] := n-1;
   PrincetonColumn: for n:= 1 to PrincetonNumber do LifeTable[LTAge,n] := PrincetonAges[n];
   0..3: begin
         n:= 0;
         for i:= 1 to stringgrid1.RowCount
         do begin
            s := Stringgrid1.Cells[AgeColumn,i];
            If s <> ''
            then begin
                 val(s, r, code);
                 if code = 0
                 then begin
                      inc(n,1);
                      LifeTable[LTAge,n] := r;
                      end
                 else BAbort := true;
                 end;
            end;
         inc(n,1);
         end;
   else BAbort := true;
   end;
   LifeTableSize := n;
   LifeTable[LTAge,LifeTableSize] := maxAge;
end;

procedure Tform2.FillDataArray;
begin
  FillLifeTableWithAges;
  case InputType of
    IPPrinceton  : FillDataArrayWithPrincetonModel;
    IPfileSpaces : FillDataArrayWithStringgrid;
    IPfileTabs   : FilldataArrayWithStringgrid;
   end;
  If BAbort then messagedlg('Your Data contain non-numbers, or you have failed to choose your columns properly.',mtinformation,[mbOk],0);
end;



procedure TForm2.BuildLifeTable;
 var i:integer;
    Age, Age2, Lhulp, Lhulp2, LhulpHalf, Lhulpkwart, Lhulpdriekwart, rRepeatSimpson:real;
    s:string;
begin
  for i:= 1 to LifetableSize
  do begin
     Age := LifeTable[LTAge,i];
     LifeTable[LTl,i] := Survivorship(Age);
     end;
  LifeTable[LTAge,LifeTableSize+1] := maxAge + 10;
  LifeTable[LTl,LifeTableSize+1] := 0;

  TableSize := LifeTableSize;

  {bereken de q's en L's}
  for i:= 1 to lifeTableSize
  do begin
     Age := LifeTable[LTage,i];
     Age2 := LifeTable[LTage,i+1];
     Lhulp := LifeTable[LTl,i];
     Lhulp2 := LifeTable[LTl,i+1];
     If ((Lhulp = 0) and (TableSize = LifeTableSize)) then TableSize := i-1;
     If TableSize = LifeTableSize
     then begin
          LifeTable[LTl,i] := Lhulp;
          LifeTable[LTq,i] := 1 - Lhulp2/ Lhulp;
          {afschatting L met herhaalde Simpson-regel voor numeriek integreren}
          LhulpHalf := Survivorship((Age2 + Age)/2);
          LhulpKwart := Survivorship((3 * Age + Age2)/4);
          LhulpDrieKwart := Survivorship((Age + 3 * Age2)/4);
          rRepeatSimpson := (Age2-Age)/12 * (Lhulp + 4 * LhulpKwart + 2 * LhulpHalf + 4 * LhulpDrieKwart + Lhulp2);
          Lifetable[LTLL,i] := rRepeatSimpson;
          end;
     end;

  {bereken de m's}
  for i:= 1 to TableSize
  do LifeTable[LTm,i] := lifeTable[LTq,i] * LifeTable[LTl,i] / LifeTable[LTLL,i];

  {bereken de T's en E's}
  LifeTable[LTT,TableSize +1] := 0;
  for i:= TableSize downto 1
  do LifeTable[LTT,i] := LifeTable[LTT,i+1] + LifeTable[LTLL,i];
  for i:= 1 to TableSize
  do LifeTable[LTE,i] := LifeTable[LTT,i] / LifeTable[LTl,i];

  {bereken AgeDistribution + deathDistribution}
  for i:= 1 to TableSize
  do begin
     Age := LifeTable[LTage,i];
     Age2 := LifeTable[LTage,i+1];
     if age2 > GrowthStartAge
     then begin
          LifeTable[LTAgeDist,i] := LifeTable[LTLL,i] * exp(-growthrate* ((Age2 + Age)/2- GrowthStartAge));
          LifeTable[LTDeathDist,i] := lifeTable[LTq,i] * LifeTable[LTl,i] * exp(-growthrate * ((Age2 + Age)/2- GrowthStartAge));
          end
     else begin
          LifeTable[LTAgeDist,i] := LifeTable[LTLL,i] * exp(-growthrate* ((Age2 + Age)/2- GrowthStartAge + GenerationGap));
          LifeTable[LTDeathDist,i] := lifeTable[LTq,i] * LifeTable[LTl,i] * exp(-growthrate * ((Age2 + Age)/2- GrowthStartAge+ GenerationGap));
          end;
     end;
  TotalPopulation := 0;
  for i:= 1 to TableSize do TotalPopulation := TotalPopulation + LifeTable[LTAgeDist,i];
  for i:= 1 to TableSize do LifeTable[LTAgeDist,i] := LifeTable[LTAgeDist,i] / TotalPopulation * 100;
  TotalDeaths := 0;
  for i:= 1 to TableSize do TotalDeaths := TotalDeaths + LifeTable[LTDeathDist,i];
  for i:= 1 to TableSize do LifeTable[LTDeathDist,i] := LifeTable[LTDeathDist,i] / TotalDeaths * 100;
  CDR := TotalDeaths / TotalPopulation;
end;


procedure TForm2.DrawLifeTable;
var i, n1, m1, n2, m2, n3, m3, n4, m4, n5, m5, n6, m6, n7, m7:integer;
    t:real;
    color:Tcolor;
begin
  n1 := unknown;  m1:= unknown;
  n2 := unknown;  m2:= unknown;
  n3 := unknown;  m3:= unknown;
  n4 := unknown;  m4:= unknown;
  n5 := unknown;  m5:= unknown;
  n6 := unknown;  m6:= unknown;
  n7 := unknown;  m7:= unknown;
  for i:= 1 to TableSize -1
  do begin
     t := LifeTable[LTage,i];
     if modelvisible
     then begin
          form1.tekenpuntQx(LifeTable[LTq,i],t,n1,m1,modelcolor,vsize,vrond,true);
          form1.tekenpuntlx(LifeTable[LTl,i],t,n2,m2,modelcolor,vsize,vrond,true);
          form1.tekenpuntMx(lifeTable[LTm,i],t,n3,m3,modelcolor,vsize,vrond,true);
          form1.tekenpuntEx(LifeTable[LTE,i],t,n4,m4,modelcolor,vsize,vrond,true);
          form1.tekenpuntAgeDist(LifeTable[LTAgeDist,i],t,n5,m5,modelcolor,vsize,vrond,true);
          form1.tekenpuntDeathDist(LifeTable[LTDeathDist,i],t,n6,m6,modelcolor,vsize,vrond,true);
          end;
     color := clBlack;
     case UsedVariable of
      LTq :form1.tekenpuntQx(DataArray[i],t,n7,m7,color,vsize,vrecht,false);
      Ltl :form1.tekenpuntlx(DataArray[i],t,n7,m7,color,vsize,vrecht,false);
      Ltm :form1.tekenpuntMx(DataArray[i],t,n7,m7,color,vsize,vrecht,false);
      LTE :form1.tekenpuntEx(DataArray[i],t,n7,m7,color,vsize,vrecht,false);
      LTAgedist :form1.tekenpuntAgeDist(DataArray[i],t,n7,m7,color,vsize,vrecht,false);
      LTDeathdist :form1.tekenpuntDeathDist(DataArray[i],t,n7,m7,color,vsize,vrecht,false);
      end;
     end;
  if modelvisible
  then with form1 do  case UsedVariable of
                        LTq :legenda({wwidth + 8*puntgrootte} xmarge + 20 * puntgrootte,200,Qxgraph);
                        LTl :Legenda(wwidth + 8*puntgrootte,200,lxgraph);
                        LTm :Legenda(wwidth + 8*puntgrootte,200,Mxgraph);
                        LTE :Legenda(wwidth + 8*puntgrootte,200,Exgraph);
                        LTAgedist:legenda(wwidth+8*puntgrootte,200,AgeDistgraph);
                        LTDeathdist:Legenda(wwidth+8*puntgrootte,200,AgeDistgraph);
                       end;
end;

procedure TForm2.writeLifeTable;
var i:integer;
    s:string;
begin
  memo1.Clear;
  s := 'Age'+ tab + 'nQx' + tab + 'lx' + tab + 'nLx' + tab + 'nMx' +tab + 'Tx' + tab + 'Ex' + tab + 'Age dist' + tab + 'Death dist';
  memo1.Lines.Add(s);
  for i:= 1 to TableSize
  do begin
     s := schrijfreal(Lifetable[LTAge,i],0)     + tab + schrijfreal(LifeTable[Ltq,i],4)  + tab
        + schrijfreal(LifeTable[Ltl,i],4)       + tab + schrijfreal(LifeTable[LtLL,i],4) + tab
        + schrijfreal(LifeTable[LTm,i],4)       + tab +schrijfreal(LifeTable[LTT,i],0)   + tab
        + schrijfreal(LifeTable[LtE,i],1)       + tab + schrijfreal(LifeTable[LTAgeDist,i],2) + tab
        + schrijfreal(LifeTable[LTDeathDist,i],2);
     memo1.Lines.Add(s);
     end;
  memo1.Lines.Add('Crude Death Rate CDR = ' + schrijfreal(CDR,4)
    + '  Total Population = ' + schrijfreal(TotalPopulation,2) + '  Total Deaths = ' + schrijfreal(TotalDeaths,2));
end;


procedure TForm2.WriteBestFit(s:string);
var ss:string;
    F:textfile;
    i:integer;
begin
  If BAbort = false
  then begin
       BuildLifeTable;
       ss := country + tab +s +tab + schrijfreal(Par[1],3) + tab + schrijfreal(Par[2],3)
               + tab + schrijfreal(Par[3],3) + tab + schrijfreal(Par[4],3)
               + tab + schrijfreal(Par[5],3) + tab + schrijfreal(fit,4) + tab + schrijfreal(Margin,4)
               + tab + inttostr(iteration);
       case applicationNumber of
         AppBrassLogitL: ss := country + tab + s + tab + logitSurvivalprocedure;
         AppBrassLogitQ: ss := country + tab + s + tab + logitQprocedure;
        end;
       memo2.lines.add(ss);

       assignfile(F,form2.logfile);
       append(F);
       writeln(F,ss);
       Flush(F);
       closeFile(F);
       end;
end;


procedure TForm2.Compare;
var ss:string;
    i:integer;
begin
  BAbort := false;
  Initieerparameters;
  form1.InitiateGraphs;
  If NOT(BAbort) then FillDataArray;
  CalculateParameters;
  If NOT(BAbort)  then begin
                       BuildLifeTable;
                       ss := schrijfreal(Par[1],3) + tab + schrijfreal(Par[2],3)
                       + tab + schrijfreal(Par[3],3) + tab + schrijfreal(Par[4],3)
                       + tab + schrijfreal(Par[5],3) + tab + schrijfreal(fit,8);
                       memo2.lines.add(ss);
                       DrawLifeTable;
                       form1.TemporarySave;
                       form1.Visible := true;
                       end;
end;


procedure TForm2.Button4Click(Sender: TObject);
var i:integer;
begin
  case Aim of
    AimFit : FindBestFitNelderMead;
    AimCompare : begin
                 BAbort := false;
                 ready := false;
                 Initieerparameters;
                 form1.InitiateGraphs;
                 If NOT(BAbort) then FillDataArray;
                 CalculateParameters;
                 Iteration := unknown;
                 end;
   end;
  If BAbort
  then messagedlg('Procedure Aborted',mtinformation,[mbOk],0)
  else begin
       writeBestfit('');
       writeLifeTable;
       DrawLifeTable;
       form1.TemporarySave;
       form1.Visible := true;
       end;
end;


procedure TForm2.Button5Click(Sender: TObject);
var F:textfile;
    s, ss:string;
    l, ntot, i:integer;
    rtot :real;
begin
  country := '';
  InputType := IPFileSpaces;
  if ApplicationNumber = unknown
  then messagedlg('Choose your application first',mtinformation,[mbOk],0)
  else begin
       opendialog1.filter := 'Textfiles (*txt)|*TXT';
       if opendialog1.execute
       then begin
            assignfile(F,opendialog1.filename);
            reset(F);
            ClearStringgrid1;
            rtot := 0;
            ntot := 0;
            while NOT(EOF(F))
            do begin
               l:= 0;
               readln(F,s);
               while s <> ''
               do begin
                  FillStringgridWithSpaceString(s,l);
                  readln(F,s);
                  end;
               if stringgrid1.Cells[0,1] <> '' then country := stringgrid1.Cells[0,1];
               FindBestFitNelderMead;
               If BAbort
               then memo2.Lines.Add('aborted')
               else begin
                    inc(ntot,1);
                    rtot := rtot + BestFit;
                    WriteBestfit(Stringgrid1.Cells[0,2] + tab + Stringgrid1.Cells[1,2]);
                    end;
               end;
            closefile(F);
            ss := schrijfreal(rtot/ntot,8);
            messagedlg(Opendialog1.filename + ' is done !!!  mean error = ' + ss,mtinformation,[mbOk],0);
            end;
       end;
end;

procedure TForm2.FindBestFitNelderMead;
var i :integer;
    r, BestFithelp :real;
    S:TSimplex;
begin
  BAbort := false;
  ready := false;
  Initieerparameters;
  form1.InitiateGraphs;
  If NOT(BAbort) then FillDataArray;
  BestFithelp := maxNumber;
  BestFit := BestFithelp -1;
  S := TSimplex.New(1,0.5,2,0.5,ParNumber);
  S.StartSimplex;
  iteration:= 0;
  While (( iteration <5000) and (S.Margin > dSimplex/10000))
  do begin
     S.Iteration;
     inc(iteration,1);
     caption := 'iteration  ' + inttostr(iteration);
     end;
  BestFit := S.CalculateValue(S.vertices[S.ILow]);
  Margin := S.Margin;
  S.free; S:= nil;
end;


function TForm2.fit:real;
var r:real;
begin
 case FittingMeasure of
  FMwrmse : r := fitWRMSE(TableSize-1,Lifetable[UsedVariable],DataArray,Lifetable[LTageDist]);
  FMlikelihood: r:= fitMaxLikeliHood;
  end;
  result := r;
end;




function TForm2.fitWRMSE(n:integer;Pmodel,Pdata,Pweights:Pattern):real;
var VarXY, Wtot, VarXX, Xmean, W, xy, xx, x :real;
    i:integer;
begin
  VarXY := 0;
  VarXX := 0;
  Xmean := 0;
  Wtot := 0;
  for i:= 1 to n
  do If ((Pdata[i] <> unknown) and (Pmodel[i] <> unknown))
     then begin
          xy := SQR(Pmodel[i] - Pdata[i]);
          xx := SQR(Pdata[i]);
          x  := Pdata[i];
          W  := Pweights[i];
          VarXY := VarXY + xy * W ;
          VarXX := VarXX + xx * W ;
          Xmean := Xmean + x * W;
          Wtot  := Wtot + W;
          end;
  VarXX := VarXX/Wtot;
  VarXY := VarXY/Wtot;
  Xmean := Xmean/Wtot;
  result := SQRT(VarXY/(varXX - Xmean * Xmean));
end;

function TForm2.LComb(n,k:integer):real;
var r:real;
    i,l:integer;
begin
  r := 0;
  if k < n/2 then l:= n-k else l:= k;
  for i:= l+1 to n do r := r + ln(i);
  for i:= 1 to n-l do r := r - ln(i);
  result := r;
end;

function TForm2.LBinom(n,k:integer; q:real):real;
var r:real;
    n1,n2:real;
begin
   if ((n > 0) and (n > k -1) and (q < 1))
   then r := LComb(n,k) + k * Ln(q) + (n - k )* Ln(1-q)
   else r := 0;
   result := r;
end;

function TForm2.fitMaxLikeliHood;
var qi, r, rtot:real;
    Li, di, i:integer;
begin
  rtot := 0;
  Li := 100000;
  for i:= 1 to TableSize -1
  do if DataArray[i] <> unknown
     then begin
          di := Round(Li * DataArray[i]);
          qi := Lifetable[LTq,i];
          r  := - LBinom(li,di,qi);
          rtot := rtot + r;
          Li := Li - di;
          end;
  result := rtot;
end;


procedure Tform2.SetIterationOptions(n:integer);
begin
  Panel6.Visible := false;
  ApplicationNumber := n;
  FittingMeasure := FMwrmse;
  Label6.caption := inttostr(ParNumber) + '-Parameter Hazard function';
  Panel1.Visible := true;
end;

procedure TForm2.Siler1Click(Sender: TObject);
begin
  Aim := Aimfit;
  ParNumber := 5;
  SetIterationOptions(AppSiler);
end;

procedure TForm2.Compare4parametermodelClick(Sender: TObject);
begin
  Aim := AimCompare;
  ParNumber := 0;
  SetIterationOptions(App4Par);
end;

procedure TForm2.CompareSiler1Click(Sender: TObject);
begin
  Aim := AimCompare;
  ParNumber := 0;
  SetIterationOptions(AppSiler);
end;

procedure TForm2.Compare3paramterModelClick(Sender: TObject);
begin
  Aim := AimCompare;
  ParNumber := 0;
  SetIterationOptions(App3Par);
end;

procedure TForm2.Compare2parametermodel1Click(Sender: TObject);
begin
  Aim := AimCompare;
  ParNumber := 0;
  SetIterationOptions(App2Par);
end;

procedure TForm2.Model4ParClick(Sender: TObject);
begin
  Aim:= AimFit;
  ParNumber := 4;
  SetIterationOptions(App4Par);
end;



procedure TForm2.Model3ParClick(Sender: TObject);
begin
  Aim := AimFit;
  ParNumber := 3;
  SetIterationOptions(App3Par);
end;

procedure TForm2.Memo1DblClick(Sender: TObject);
begin
  if memo1.Lines[0] = 'Hugo'
  then begin
      panel2.Visible := true;
      panel5.Visible := true;
      end;
end;

procedure TForm2.Model2ParClick(Sender: TObject);
begin
  Aim := AimFit;
  ParNumber := 2;
  SetIterationOptions(App2Par);
end;


procedure TForm2.BraaslogitQprocedure1Click(Sender: TObject);
begin
   Aim := AimFit;
   ParNumber := 2;
   SetIterationOptions(AppBrassLogitQ);
   panel6.visible := true;
end;

procedure TForm2.BrasslogitSurvivorClick(Sender: TObject);
begin
  Aim := AimFit;
  ParNumber := 2;
  SetIterationOptions(AppBrassLogitL);
end;

function TForm2.Survivorship(x: Real):real;
var r0, r :real;
begin
  case ApplicationNumber of
    AppSiler:
                   begin   //Siler:
                   r0 := 1 / exp( 1/Par[3] * exp(Par[5])
                   -1/Par[1] * exp(-Par[1] * Par[4]));

                   r := r0 * exp( 1/Par[3] * exp(-Par[3]*x + Par[5])
                        - exp(Par[2])*x
                        -1/par[1] * exp(Par[1]*(x - Par[4])));
                   end;

    App4Par, App3Par, App2Par, AppBrassLogitL, AppBrassLogitQ:
                   begin   //power formula for immature term, logistic for elderly term:

                   r0 := 1 / exp( -ksiPar/(1-Par[3]) * exp((1-par[3])*ln(1.5))
                        - exp(Par[2])/par[5] * ln(1 + exp(-16))
                        -1/Par[1] * ln(1 + exp(-Par[1]*Par[4])));

                   r := r0 * exp(-ksiPar/(1-Par[3]) * exp((1-par[3])*ln(x+1.5))

                        - exp(Par[2])*x
                        - exp(Par[2])/par[5] * ln(1 + exp(par[5]*x - 16))

                        -1/Par[1] * ln(1 + exp(Par[1] * (x - Par[4]))));
                   end;


    end;
  result := r;
end;

procedure TForm2.CalculateParameters;
begin

  case ApplicationNumber of
    AppSiler               : begin
                             If par[1] < 0.01 then par[1] := 0.01;
                             If par[2] < -20 then par[2]:= -20;
                             If par[3] < 0.01 then par[3] := 0.01;
                             if par[3] > 30 then par[3] := 30;
                             if par[5] < -20 then par[5] := -20;
                             if par[5] > 5 then Par[5] := 5;
                             end;
    App4Par                : begin
                             Par[5] := 1;
                             If Par[3] < 1.01 then Par[3] := 1.01;
                             end;
    App3Par               : begin
                             Par[5] := 1;
                             Par[4] := 100;
                             If Par[3] < 1.01 then Par[3] := 1.01;
                             end;
    App2Par, AppBrassLogitL, AppBrassLogitQ
                           : begin
                             Par[5] := 1;
                             Par[4] := 100;
                             {Par[3] := 121.7 * Par[1]  -0.97 * Par[2]  -14.5;}
                             Par[3] := 119.3 * Par[1] -1.01 * Par[2] -14.5;
                             If Par[3] < 2.7 then Par[3] := 2.7;
                             end;
  end;
end;

end.



unit Unit3;

interface

uses
  Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
  Vcl.Controls, Vcl.Forms, Vcl.Dialogs, Vcl.StdCtrls, Vcl.Samples.Spin,
  Vcl.ExtCtrls;

type
  TForm3 = class(TForm)
    Panel1: TPanel;
    Label1: TLabel;
    SpinEdit1: TSpinEdit;
    CheckBox1: TCheckBox;
    CheckBox2: TCheckBox;
    procedure FormClose(Sender: TObject; var Action: TCloseAction);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form3: TForm3;

implementation

{$R *.dfm}

uses Unit1, Unit2;

procedure TForm3.FormClose(Sender: TObject; var Action: TCloseAction);
begin
  form1.fontsize := spinedit1.Value;
  if checkbox1.Checked then form1.raster := true else form1.raster := false;
  if checkbox2.Checked then form2.modelvisible := true else form2.modelvisible := false;
end;




end.



unit Unit1;

interface

uses
  Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
  Vcl.Controls, Vcl.Forms, Vcl.Dialogs, ExtCtrls, StdCtrls, ComCtrls, Vcl.Menus;

const
small = 0.0000001;
hheight = 400;
wwidth = 600;
vrond = 1;
vrecht=2;



type
  TForm1 = class(TForm)
    MainMenu1: TMainMenu;
    fil1: TMenuItem;
    save1: TMenuItem;
    PageControl1: TPageControl;
    TabSheet1: TTabSheet;
    TabSheet2: TTabSheet;
    TabSheet3: TTabSheet;
    Tabsheet4: TTabSheet;
    Tabsheet5: TTabSheet;
    TabSheet6: TTabSheet;
    ImageQx: TImage;
    Imagelx: TImage;
    ImageMx: TImage;
    ImageEx: TImage;
    ImageAgeDist: TImage;
    ImageDeathDist: TImage;
    procedure save1Click(Sender: TObject);
    procedure FormCreate(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
    deciaantal:integer;
    puntgrootte:integer;
    raster :boolean;
    fontsize:integer;
    xmarge:integer;
    qmax:real;
    lmax:real;
    Mmax:real;
    Emax:real;
    AgeDistmax:real;
    DeathDistmax:real;
    Qxgraph:Tbitmap;
    lxgraph:Tbitmap;
    Mxgraph:Tbitmap;
    Exgraph:Tbitmap;
    AgeDistgraph:Tbitmap;
    DeathDistgraph:Tbitmap;
    procedure Legenda(xx,yy:integer;var grafiek:Tbitmap);
    procedure TemporarySave;
    procedure InitiateGraphs;
    procedure TekenAssenstelsel(xxmin,xxmax,xxstap,yymin,yymax,yystap:real;ydeci:integer;var grafiek:Tbitmap;xxstring,yystring:string);
    procedure tekenpunt(var xn,yn,xnn,ynn:integer;kleur:Tcolor;grootte,vorm:integer;verbinding:boolean;var grafiek:Tbitmap);
    procedure tekenpuntQx(qq,xx:real;var qn,xn:integer;kleur:Tcolor;grootte,vorm:integer;verbinding:boolean);
    procedure tekenpuntlx(ll,xx:real;var ln,xn:integer;kleur:Tcolor;grootte,vorm:integer;verbinding:boolean);
    procedure tekenpuntMx(MM,xx:real;var Mn,xn:integer;kleur:Tcolor;grootte,vorm:integer;verbinding:boolean);
    procedure tekenpuntEx(EE,xx:real;var En,xn:integer;kleur:Tcolor;grootte,vorm:integer;verbinding:boolean);
    procedure tekenpuntAgeDist(AD,xx:real;var ADn,xn:integer;kleur:Tcolor;grootte,vorm:integer;verbinding:boolean);
    procedure tekenpuntDeathDist(DD,xx:real;var DDn,xn:integer;kleur:Tcolor;grootte,vorm:integer;verbinding:boolean);
  end;

function known(r:real):boolean;

var
  Form1: TForm1;
  vsize :integer;
  modelcolor :Tcolor;



implementation

{$R *.dfm}

uses unit2;

function known(r:real):boolean;
begin
  if ABS(r - unknown) < small
  then result := false
  else result := true;
end;

procedure TForm1.save1Click(Sender: TObject);
var s:string;
    i:integer;
begin
  s := DateToStr(Date) + ' ' + Timetostr(Time);
  for i:= 1 to length(s) do if s[i] in ['/',':'] then s[i] := '_';
  case pagecontrol1.ActivePageindex of
   0 : imageQx.Picture.SaveToFile(form2.logdir + 'graph_'+ tabsheet1.caption + '_'+ s + '.bmp');
   1 : imagelx.Picture.SaveToFile(form2.logdir + 'graph_'+ tabsheet2.caption + '_'+ s + '.bmp');
   2 : imageMx.Picture.SaveToFile(form2.logdir + 'graph_'+ tabsheet3.caption + '_'+ s + '.bmp');
   3 : imageEx.Picture.SaveToFile(form2.logdir + 'graph_'+ tabsheet4.Caption + '_'+ s + '.bmp');
   4 : imageAgeDist.Picture.SaveToFile(form2.logdir + 'graph_'+ tabsheet5.caption + '_'+ s + '.bmp');
   5 : imageDeathDist.Picture.SaveToFile(form2.logdir + 'graph_'+ tabsheet6.caption + '_'+ s + '.bmp');
   end;
end;

procedure TForm1.FormCreate(Sender: TObject);
begin
  deciaantal := 2;
  puntgrootte := 2;
  vsize := 3;
  modelcolor := clblack;
  fontsize := 12;
  raster := false;
  xmarge := 70;
end;


procedure Tform1.TemporarySave;
begin
  Qxgraph.savetofile(form2.logdir + 'Qx.bmp');                     Qxgraph.free;  Qxgraph := nil;
  lxgraph.savetofile(form2.logdir + 'lx.bmp');                     lxgraph.free;  lxgraph := nil;
  Mxgraph.SaveToFile(form2.logdir + 'Mx.bmp');                     Mxgraph.Free;  Mxgraph := nil;
  Exgraph.savetofile(form2.logdir + 'Ex.bmp');                     Exgraph.free;  Exgraph := nil;
  AgeDistgraph.savetofile(form2.logdir + 'AgeDist.bmp');           AgeDistgraph.free; AgeDistgraph := nil;
  DeathDistgraph.savetofile(form2.logdir + 'DeathDist.bmp');       DeathDistgraph.Free; DeathDistgraph := nil;
  form1.imageQx.picture.loadfromfile(form2.logdir + 'Qx.bmp');
  form1.imagelx.picture.loadfromfile(form2.logdir + 'lx.bmp');
  form1.imageMx.picture.loadfromfile(form2.logdir + 'Mx.bmp');
  form1.ImageEx.Picture.LoadFromFile(form2.logdir + 'Ex.bmp');
  form1.imageAgeDist.picture.loadfromfile(form2.logdir + 'AgeDist.bmp');
  form1.imageDeathDist.Picture.loadfromfile(form2.logdir + 'DeathDist.bmp');
end;




procedure Tform1.InitiateGraphs;
begin
  qmax := 1;
  lmax := 1;
  Mmax := 1;
  Emax := 100;
  AgeDistmax := 20;
  DeathDistmax := 20;
  Tekenassenstelsel(0,maxAge,10,0,qmax,qmax/10,2,Qxgraph,'age-group starting age x','mortality probability nQx');
  Tekenassenstelsel(0,maxAge,10,0,lmax,lmax/10,2,lxgraph,'age  x','survivorship lx');
  Tekenassenstelsel(0,maxAge,10,0,Mmax,Mmax/10,2,Mxgraph,'age-group starting age x','death rate nMx');
  Tekenassenstelsel(0,maxAge,10,0,Emax,Emax/10,0,Exgraph,'age  x','life expectancy E0x (yr)');
  Tekenassenstelsel(0,maxAge,10,0,AgeDistmax,AgeDistmax/10,0,AgeDistgraph,'age-group starting age x','Age Distribution (%)');
  Tekenassenstelsel(0,maxAge,10,0,DeathDistmax,DeathDistmax/10,0,DeathDistgraph,'age-group starting age x','Death Distribution (%)');
end;



procedure TForm1.TekenAssenstelsel(xxmin,xxmax,xxstap,yymin,yymax,yystap:real;Ydeci:integer;var grafiek:Tbitmap;xxstring,yystring:string);
var xx,yy,i,l,l2,l3,l4,imin,imax, hulpsize:integer;
    xi,yi:real;
    s:string;
begin
  grafiek := Tbitmap.create;
  grafiek.width := wwidth + 100;
  grafiek.height := hheight + 100;
  grafiek.canvas.font.size := fontsize;

  if raster
  then grafiek.canvas.pen.color := clsilver
  else grafiek.Canvas.Pen.Color := clblack;

  //plaats horizontale lijnen
  imin := 0; while imin*yystap > yymin do dec(imin,1);
  imax := 0; while imax*yystap < yymax do inc(imax,1);
  for i:= imin-5 to imax+5
  do begin
     yi := i*yystap;
     xx := round(-xxmin/(xxmax-xxmin) * wwidth) + xmarge;
     yy := round(-(yi-yymax)/(yymax-yymin)*hheight) + 50;
     if raster
     then begin
          grafiek.canvas.moveto(0,yy+1);
          grafiek.canvas.lineto(wwidth+100,yy+1);
          end
     else begin
          grafiek.Canvas.moveto(xx-1,yy+1);
          grafiek.Canvas.lineto(xx+10,yy+1);
          end;
     str(yi:3:ydeci,s);
     l := grafiek.Canvas.textwidth(s) + puntgrootte;
     if ((i > 0) and (xx > 0) and (xx < grafiek.width))
     then grafiek.canvas.textout(xx - l,yy-5,s)
     else grafiek.canvas.textout(0 - l,yy-5,s);
     end;

  // plaats verticale lijnen
  imin := 0; while imin*xxstap > xxmin do dec(imin,1);
  imax := 0; while imax*xxstap < xxmax do inc(imax,1);
  for i:= imin-5 to imax+5
  do begin
     xi := i*xxstap;
     xx := round((xi-xxmin)/(xxmax-xxmin) * wwidth) + xmarge;
     yy := round(yymax/(yymax-yymin)*hheight) + 50;
     if raster
     then begin
          grafiek.canvas.moveto(xx,0);
          grafiek.canvas.lineto(xx,hheight+100);
          end
     else begin
          grafiek.Canvas.MoveTo(xx,yy-10);
          grafiek.Canvas.LineTo(xx,yy);
          end;
     str(xi:3:0,s);
     l2 := grafiek.canvas.textwidth(s);
     if ((yy > 0) and (yy < grafiek.height))
     then grafiek.canvas.textout(xx - l2 div 2,yy+puntgrootte,s)
     else grafiek.canvas.textout(xx - l2 div 2,hheight,s);
     end;

  // plaats de twee assen
  grafiek.canvas.pen.color := clblack;
  xx := round(-xxmin/(xxmax-xxmin) * wwidth) + xmarge;
  if ((xx < 0) or (xx > grafiek.width)) then xx := 0;
  grafiek.canvas.moveto(xx,0);
  grafiek.canvas.lineto(xx,hheight+100);
  yy := round(yymax/(yymax-yymin)*hheight) + 50;
  if ((yy < 0) or (yy > grafiek.height)) then yy := hheight;
  grafiek.canvas.moveto(0,yy);
  grafiek.canvas.lineto(wwidth+100,yy);

  //plaats de grootheden bij de assen
  hulpsize := grafiek.canvas.font.size;
  grafiek.canvas.font.size := hulpsize + 2* puntgrootte;
  {grafiek.Canvas.Font.Style := [fsbold];}
  grafiek.Canvas.Font.Orientation := 900;
  l3 := grafiek.Canvas.TextHeight(yystring);
  grafiek.canvas.textout(xx -l - l3- puntgrootte - fontsize div 2,{hheight div 2} hheight - 50,yystring);
  grafiek.Canvas.Font.Orientation := 0;
  l4 := grafiek.Canvas.TextWidth(xxstring);
  if l4 > wwidth / 3
  then  grafiek.canvas.TextOut(wwidth div 3,yy + puntgrootte + l3 - fontsize div 2,xxstring)
  else  grafiek.canvas.TextOut(wwidth div 2,yy + puntgrootte + l3 - fontsize div 2,xxstring);
  grafiek.canvas.font.size := hulpsize;
  {grafiek.Canvas.Font.Style := [];}
end;

procedure Tform1.tekenpunt(var xn,yn,xnn,ynn:integer;kleur:Tcolor;grootte,vorm:integer;verbinding:boolean;var grafiek:Tbitmap);
begin
  if ((xnn > 0) and (xnn < grafiek.width) and (ynn > 0) and (ynn < grafiek.height))
  then begin
       {grafiek.canvas.pixels[xn,yn] := kleur; }
       grafiek.canvas.brush.color := kleur;
       grafiek.canvas.brush.style := bsclear;
       grafiek.canvas.pen.color := kleur;
       grafiek.Canvas.Pen.width := 2;
       If ((verbinding) and (xn > 0) and (xn < grafiek.width) and (yn > 0) and (yn < grafiek.height))
       then begin
            grafiek.canvas.moveto(xn,yn);
            grafiek.canvas.lineto(xnn,ynn);
            end;
       case vorm of
        vrond: grafiek.canvas.ellipse(xnn-grootte*puntgrootte,ynn-grootte*puntgrootte,xnn+grootte*puntgrootte,ynn+grootte*puntgrootte);
        vrecht: begin
                grafiek.Canvas.Brush.Style := bssolid;
                grafiek.Canvas.Brush.Color := kleur;
                grafiek.canvas.rectangle(xnn-grootte*puntgrootte,ynn-grootte*puntgrootte,xnn+grootte*puntgrootte,ynn+grootte*puntgrootte);
                end;
        end;
       xn := xnn;
       yn := ynn;
       end;
end;

procedure Tform1.tekenpuntQx(qq,xx:real;var qn,xn:integer;kleur:Tcolor;grootte,vorm:integer;verbinding:boolean);
var qnn,xnn:integer;
begin
  if known(xx) then xnn := round(xx/maxAge * wwidth) + xmarge else xnn := unknown;
  if known(qq) then qnn := round(-(qq-qmax)/qmax * hheight) + 50 else qnn := unknown;
  Tekenpunt(xn,qn,xnn,qnn,kleur,grootte,vorm,verbinding,Qxgraph);
end;

procedure Tform1.tekenpuntlx(ll,xx:real;var ln,xn:integer;kleur:Tcolor;grootte,vorm:integer;verbinding:boolean);
var lnn,xnn:integer;
begin
  if known(xx) then xnn := round(xx/maxAge * wwidth) + xmarge else xnn := unknown;
  if known(ll) then lnn := round(-(ll-lmax)/lmax * hheight) + 50 else lnn := unknown;
  Tekenpunt(xn,ln,xnn,lnn,kleur,grootte,vorm,verbinding,lxgraph);
end;

procedure Tform1.tekenpuntMx(MM,xx:real;var Mn,xn:integer;kleur:Tcolor;grootte,vorm:integer;verbinding:boolean);
var Mnn,xnn:integer;
begin
  if known(xx) then xnn := round(xx/maxAge * wwidth) + xmarge else xnn := unknown;
  if known(MM) then Mnn := round(-(MM-Mmax)/Mmax * hheight) + 50 else Mnn := unknown;
  Tekenpunt(xn,Mn,xnn,Mnn,kleur,grootte,vorm,verbinding,Mxgraph);
end;

procedure Tform1.tekenpuntEx(EE,xx:real;var En,xn:integer;kleur:Tcolor;grootte,vorm:integer;verbinding:boolean);
var Enn,xnn:integer;
begin
  if known(xx) then xnn := round(xx/maxAge * wwidth) + xmarge else xnn := unknown;
  if known(EE) then Enn := round(-(EE-Emax)/Emax * hheight) + 50 else Enn := unknown;
  Tekenpunt(xn,En,xnn,Enn,kleur,grootte,vorm,verbinding,Exgraph);
end;

procedure Tform1.tekenpuntAgeDist(AD,xx:real;var ADn,xn:integer;kleur:Tcolor;grootte,vorm:integer;verbinding:boolean);
var ADnn,xnn:integer;
begin
  if known(xx) then xnn := round(xx/maxAge * wwidth) + xmarge else xnn := unknown;
  if known(AD) then ADnn := round(-(AD-AgeDistmax)/AgeDistmax * hheight) + 50 else ADnn := unknown;
  Tekenpunt(xn,ADn,xnn,ADnn,kleur,grootte,vorm,verbinding,AgeDistgraph);
end;

procedure Tform1.tekenpuntDeathDist(DD,xx:real;var DDn,xn:integer;kleur:Tcolor;grootte,vorm:integer;verbinding:boolean);
var DDnn,xnn:integer;
begin
  if known(xx) then xnn := round(xx/maxAge * wwidth) + xmarge else xnn := unknown;
  if known(DD) then DDnn := round(-(DD-DeathDistmax)/DeathDistmax * hheight) + 50 else DDnn := unknown;
  Tekenpunt(xn,DDn,xnn,DDnn,kleur,grootte,vorm,verbinding,DeathDistgraph);
end;

   { procedure Tform1.tekenpuntQx(qq,xx:real;var qn,xn:integer;kleur:Tcolor;grootte:integer;verbinding:boolean);
var qnn,xnn:integer;
begin
  if known(xx) then xnn := round((xx-xmin)/(xmax-xmin) * wwidth) + 50 else xnn := unknown;
  if known(qq) then qnn := round(-(qq-qmax)/(qmax-qmin) * hheight) + 50 else qnn := unknown;
  Tekenpunt(xn,qn,xnn,qnn,kleur,grootte,verbinding,Qxgrafiek);
end;  }

  { procedure Tform1.tekenpuntQx(qq,xx:real;var qn,xn:integer;kleur:Tcolor;grootte:integer;verbinding:boolean);
var qnn,xnn:integer;
begin
  if known(xx) then xnn := round((xx-xmin)/(xmax-xmin) * wwidth) + 50 else xnn := unknown;
  if known(qq) then qnn := round(-(qq-qmax)/(qmax-qmin) * hheight) + 50 else qnn := unknown;
  Tekenpunt(xn,qn,xnn,qnn,kleur,grootte,verbinding,Qxgrafiek);
end;  }

procedure Tform1.Legenda(xx,yy:integer;var grafiek:Tbitmap);
begin
  tekenpunt(xx,yy,xx,yy,modelcolor,5,vrond,false,grafiek);
  grafiek.canvas.brush.color := clwhite;
  grafiek.canvas.pen.color := clblack;
  grafiek.Canvas.TextOut(xx + 8*puntgrootte,yy - 4 * puntgrootte,'model');
  inc(yy,12 * puntgrootte);
  tekenpunt(xx, yy,xx,yy,clblack,5,vrecht,false,grafiek);
  grafiek.canvas.brush.color := clwhite;
  grafiek.canvas.pen.color := clblack;
  grafiek.Canvas.TextOut(xx + 8*puntgrootte,yy - 4 * puntgrootte,'data');
end;


end.


object Form2: TForm2
  Left = 0
  Top = 0
  Caption = 'Form2'
  ClientHeight = 436
  ClientWidth = 1297
  Color = clBtnFace
  Font.Charset = DEFAULT_CHARSET
  Font.Color = clWindowText
  Font.Height = -11
  Font.Name = 'Tahoma'
  Font.Style = []
  Menu = MainMenu1
  OldCreateOrder = False
  OnCreate = FormCreate
  PixelsPerInch = 96
  TextHeight = 13
  object Panel1: TPanel
    Left = 495
    Top = 8
    Width = 226
    Height = 384
    TabOrder = 0
    Visible = False
    object Label1: TLabel
      Left = 21
      Top = 50
      Width = 14
      Height = 15
      Caption = 'P1'
      Font.Charset = DEFAULT_CHARSET
      Font.Color = clWindowText
      Font.Height = 15
      Font.Name = 'Tahoma'
      Font.Style = []
      ParentFont = False
    end
    object Label2: TLabel
      Left = 20
      Top = 81
      Width = 14
      Height = 15
      Caption = 'P2'
      Font.Charset = DEFAULT_CHARSET
      Font.Color = clWindowText
      Font.Height = 15
      Font.Name = 'Tahoma'
      Font.Style = []
      ParentFont = False
    end
    object Label3: TLabel
      Left = 21
      Top = 108
      Width = 14
      Height = 15
      Caption = 'P3'
      Font.Charset = DEFAULT_CHARSET
      Font.Color = clWindowText
      Font.Height = 15
      Font.Name = 'Tahoma'
      Font.Style = []
      ParentFont = False
    end
    object Label4: TLabel
      Left = 21
      Top = 135
      Width = 14
      Height = 15
      Caption = 'P4'
      Font.Charset = DEFAULT_CHARSET
      Font.Color = clWindowText
      Font.Height = 15
      Font.Name = 'Tahoma'
      Font.Style = []
      ParentFont = False
    end
    object Label5: TLabel
      Left = 20
      Top = 163
      Width = 14
      Height = 15
      Caption = 'P5'
      Font.Charset = DEFAULT_CHARSET
      Font.Color = clWindowText
      Font.Height = 15
      Font.Name = 'Tahoma'
      Font.Style = []
      ParentFont = False
    end
    object Label6: TLabel
      Left = 20
      Top = 9
      Width = 193
      Height = 20
      Caption = 'parameters hazard function'
      Font.Charset = DEFAULT_CHARSET
      Font.Color = clWindowText
      Font.Height = 20
      Font.Name = 'Tahoma'
      Font.Style = []
      ParentFont = False
    end
    object Label7: TLabel
      Left = 22
      Top = 326
      Width = 70
      Height = 15
      Caption = 'growth rate '
      Font.Charset = DEFAULT_CHARSET
      Font.Color = clWindowText
      Font.Height = 15
      Font.Name = 'Tahoma'
      Font.Style = []
      ParentFont = False
    end
    object Label8: TLabel
      Left = 20
      Top = 252
      Width = 149
      Height = 20
      Caption = 'Demographic growth'
      Font.Charset = DEFAULT_CHARSET
      Font.Color = clWindowText
      Font.Height = 20
      Font.Name = 'Tahoma'
      Font.Style = []
      ParentFont = False
    end
    object Label10: TLabel
      Left = 24
      Top = 286
      Width = 84
      Height = 15
      Caption = 'type of growth'
      Font.Charset = DEFAULT_CHARSET
      Font.Color = clWindowText
      Font.Height = 15
      Font.Name = 'Tahoma'
      Font.Style = []
      ParentFont = False
    end
    object Label14: TLabel
      Left = 24
      Top = 347
      Width = 181
      Height = 13
      Caption = '(for example 0.05 means 5% growth)'
    end
    object Edit1: TEdit
      Left = 40
      Top = 48
      Width = 57
      Height = 21
      TabOrder = 0
      Text = '0.1'
    end
    object Edit2: TEdit
      Left = 40
      Top = 75
      Width = 57
      Height = 21
      TabOrder = 1
      Text = '-6'
    end
    object Edit3: TEdit
      Left = 40
      Top = 102
      Width = 57
      Height = 21
      TabOrder = 2
      Text = '5'
    end
    object Edit4: TEdit
      Left = 40
      Top = 129
      Width = 57
      Height = 21
      TabOrder = 3
      Text = '100'
    end
    object Edit5: TEdit
      Left = 40
      Top = 157
      Width = 57
      Height = 21
      TabOrder = 4
      Text = '1'
    end
    object Edit6: TEdit
      Left = 98
      Top = 320
      Width = 43
      Height = 21
      TabOrder = 5
      Text = '0'
    end
    object ComboBox1: TComboBox
      Left = 114
      Top = 284
      Width = 97
      Height = 21
      ItemIndex = 0
      TabOrder = 6
      Text = 'birth growth'
      Items.Strings = (
        'birth growth'
        'migration growth')
    end
    object Panel6: TPanel
      Left = 0
      Top = 196
      Width = 193
      Height = 50
      TabOrder = 7
      Visible = False
      object Label15: TLabel
        Left = 18
        Top = 1
        Width = 66
        Height = 16
        Caption = 'Brass Logit '
        Font.Charset = DEFAULT_CHARSET
        Font.Color = clWindowText
        Font.Height = -13
        Font.Name = 'Tahoma'
        Font.Style = []
        ParentFont = False
      end
      object Label16: TLabel
        Left = 16
        Top = 23
        Width = 119
        Height = 13
        Caption = 'ref. Princeton West level'
      end
      object SpinEdit3: TSpinEdit
        Left = 141
        Top = 19
        Width = 39
        Height = 22
        MaxValue = 25
        MinValue = 1
        TabOrder = 0
        Value = 12
      end
    end
    object Button4: TButton
      Left = 136
      Top = 77
      Width = 75
      Height = 25
      Caption = 'OK'
      TabOrder = 8
      OnClick = Button4Click
    end
  end
  object Panel2: TPanel
    Left = 8
    Top = 7
    Width = 481
    Height = 394
    TabOrder = 1
    Visible = False
    object Label9: TLabel
      Left = 16
      Top = 14
      Width = 69
      Height = 20
      Caption = 'Life Table'
      Font.Charset = DEFAULT_CHARSET
      Font.Color = clWindowText
      Font.Height = 20
      Font.Name = 'Tahoma'
      Font.Style = []
      ParentFont = False
    end
    object Label12: TLabel
      Left = 169
      Top = 350
      Width = 111
      Height = 15
      Caption = 'taken from column: '
      Font.Charset = DEFAULT_CHARSET
      Font.Color = clWindowText
      Font.Height = 15
      Font.Name = 'Tahoma'
      Font.Style = []
      ParentFont = False
    end
    object StringGrid1: TStringGrid
      Left = 16
      Top = 40
      Width = 441
      Height = 233
      FixedCols = 0
      Options = [goFixedVertLine, goFixedHorzLine, goVertLine, goHorzLine, goRangeSelect, goEditing]
      TabOrder = 0
    end
    object ComboBox2: TComboBox
      Left = 18
      Top = 345
      Width = 145
      Height = 21
      ItemIndex = 0
      TabOrder = 1
      Text = 'mortality probability nQx'
      Items.Strings = (
        'mortality probability nQx'
        'survivorship l(x),     l(0)=100000'
        'hazard rate nMx'
        'life expectance e0    (in years)'
        'Age distribution      (in %)'
        'Deaths distribution      (in %)')
    end
    object SpinEdit2: TSpinEdit
      Left = 286
      Top = 345
      Width = 41
      Height = 22
      MaxValue = 0
      MinValue = 0
      TabOrder = 2
      Value = 4
      OnChange = SpinEdit2Change
    end
    object ComboBox3: TComboBox
      Left = 18
      Top = 297
      Width = 145
      Height = 21
      ItemIndex = 0
      TabOrder = 3
      Text = 'HMD age-groups'
      Items.Strings = (
        'HMD age-groups'
        'HMD age'
        'Princeton age-groups'
        'Ages taken from column 0'
        'Ages taken from column 1'
        'Ages taken from column 2'
        'Ages taken from column 3')
    end
    object Panel5: TPanel
      Left = 352
      Top = 289
      Width = 105
      Height = 96
      TabOrder = 4
      Visible = False
      object Button5: TButton
        Left = 16
        Top = 33
        Width = 75
        Height = 25
        Caption = 'Button5'
        TabOrder = 0
        OnClick = Button5Click
      end
    end
  end
  object Panel3: TPanel
    Left = 727
    Top = 8
    Width = 562
    Height = 393
    TabOrder = 2
    object Memo1: TMemo
      Left = 8
      Top = 16
      Width = 545
      Height = 245
      Lines.Strings = (
        'Memo1')
      TabOrder = 0
      OnDblClick = Memo1DblClick
    end
    object Memo2: TMemo
      Left = 8
      Top = 284
      Width = 545
      Height = 89
      Lines.Strings = (
        'Memo2')
      TabOrder = 1
    end
  end
  object Panel4: TPanel
    Left = 8
    Top = 8
    Width = 481
    Height = 78
    TabOrder = 3
    Visible = False
    object Label11: TLabel
      Left = 137
      Top = 42
      Width = 33
      Height = 18
      Caption = 'level:'
      Font.Charset = DEFAULT_CHARSET
      Font.Color = clWindowText
      Font.Height = -15
      Font.Name = 'Tahoma'
      Font.Style = []
      ParentFont = False
    end
    object Label13: TLabel
      Left = 24
      Top = 16
      Width = 115
      Height = 20
      Caption = 'Princeton model'
      Font.Charset = DEFAULT_CHARSET
      Font.Color = clWindowText
      Font.Height = 20
      Font.Name = 'Tahoma'
      Font.Style = []
      ParentFont = False
    end
    object ComboBox4: TComboBox
      Left = 25
      Top = 42
      Width = 97
      Height = 21
      ItemIndex = 2
      TabOrder = 0
      Text = 'West Combi males and females'
      Items.Strings = (
        'West females'
        'West males'
        'West Combi males and females')
    end
    object SpinEdit1: TSpinEdit
      Left = 176
      Top = 40
      Width = 41
      Height = 22
      MaxValue = 25
      MinValue = 1
      TabOrder = 1
      Value = 1
    end
  end
  object MainMenu1: TMainMenu
    Left = 472
    Top = 400
    object file1: TMenuItem
      Caption = 'data'
      object OpenLifeTabledatafilewithspaties1: TMenuItem
        Caption = 'Open Life Table, data separated by spaces'
        OnClick = OpenLifeTabledatafilewithspaties1Click
      end
      object Princetonmodellifetable1: TMenuItem
        Caption = 'Princeton model life table'
        OnClick = Princetonmodellifetable1Click
      end
    end
    object Hazardfunction1: TMenuItem
      Caption = 'best fit'
      object Siler1: TMenuItem
        Caption = 'Siler'
        OnClick = Siler1Click
      end
      object Model4Par: TMenuItem
        Caption = '4-parameter model'
        OnClick = Model4ParClick
      end
      object Model3Par: TMenuItem
        Caption = '3-parameter model'
        OnClick = Model3ParClick
      end
      object Model2Par: TMenuItem
        Caption = '2-parameter model'
        OnClick = Model2ParClick
      end
      object BrasslogitSurvivor: TMenuItem
        Caption = 'Brass logit Survivor 1841'
        OnClick = BrasslogitSurvivorClick
      end
      object BraaslogitQprocedure1: TMenuItem
        Caption = 'Braas logit Q procedure'
        OnClick = BraaslogitQprocedure1Click
      end
    end
    object comparefit1: TMenuItem
      Caption = 'compare'
      object CompareSiler1: TMenuItem
        Caption = 'Siler'
        OnClick = CompareSiler1Click
      end
      object Compare4parametermodel: TMenuItem
        Caption = '4-parameter model'
        OnClick = Compare4parametermodelClick
      end
      object Compare3paramterModel: TMenuItem
        Caption = '3-parameter model'
        OnClick = Compare3paramterModelClick
      end
      object Compare2Parametermodel: TMenuItem
        Caption = '2 parameter model'
        OnClick = Compare2parametermodel1Click
      end
    end
  end
  object OpenDialog1: TOpenDialog
    Left = 528
    Top = 400
  end
end

object Form1: TForm1
  Left = 0
  Top = 0
  Caption = 'Form1'
  ClientHeight = 592
  ClientWidth = 734
  Color = clBtnFace
  Font.Charset = DEFAULT_CHARSET
  Font.Color = clWindowText
  Font.Height = -11
  Font.Name = 'Tahoma'
  Font.Style = []
  Menu = MainMenu1
  OldCreateOrder = False
  OnCreate = FormCreate
  PixelsPerInch = 96
  TextHeight = 13
  object PageControl1: TPageControl
    Left = 8
    Top = 8
    Width = 700
    Height = 550
    ActivePage = Tabsheet4
    TabOrder = 0
    object TabSheet1: TTabSheet
      Caption = 'mortality probability nQ(x)'
      ExplicitLeft = 0
      ExplicitTop = 0
      ExplicitWidth = 0
      ExplicitHeight = 0
      object ImageQx: TImage
        Left = 0
        Top = 0
        Width = 692
        Height = 522
        Align = alClient
        ExplicitLeft = 104
        ExplicitTop = 32
        ExplicitWidth = 105
        ExplicitHeight = 105
      end
    end
    object TabSheet2: TTabSheet
      Caption = 'survivorship l(x)'
      ImageIndex = 1
      ExplicitLeft = 0
      ExplicitTop = 0
      ExplicitWidth = 0
      ExplicitHeight = 0
      object Imagelx: TImage
        Left = 0
        Top = 0
        Width = 692
        Height = 522
        Align = alClient
        ExplicitLeft = 184
        ExplicitTop = 88
        ExplicitWidth = 105
        ExplicitHeight = 105
      end
    end
    object TabSheet3: TTabSheet
      Caption = 'hazard rate nM(x)'
      ImageIndex = 2
      ExplicitLeft = 0
      ExplicitTop = 0
      ExplicitWidth = 0
      ExplicitHeight = 0
      object ImageMx: TImage
        Left = 0
        Top = 0
        Width = 692
        Height = 522
        Align = alClient
        ExplicitLeft = 512
        ExplicitTop = 88
        ExplicitWidth = 105
        ExplicitHeight = 105
      end
    end
    object Tabsheet4: TTabSheet
      Caption = 'Life Expectancy E0(x)'
      ImageIndex = 3
      object ImageEx: TImage
        Left = 0
        Top = 0
        Width = 692
        Height = 522
        Align = alClient
        ExplicitLeft = 264
        ExplicitTop = 200
        ExplicitWidth = 105
        ExplicitHeight = 105
      end
    end
    object Tabsheet5: TTabSheet
      Caption = 'Age distribution'
      ImageIndex = 4
      ExplicitLeft = 0
      ExplicitTop = 0
      ExplicitWidth = 0
      ExplicitHeight = 0
      object ImageAgeDist: TImage
        Left = 0
        Top = 0
        Width = 692
        Height = 522
        Align = alClient
        ExplicitLeft = 368
        ExplicitTop = 120
        ExplicitWidth = 105
        ExplicitHeight = 105
      end
    end
    object TabSheet6: TTabSheet
      Caption = 'Deaths Distribution'
      ImageIndex = 5
      ExplicitLeft = 0
      ExplicitTop = 0
      ExplicitWidth = 0
      ExplicitHeight = 0
      object ImageDeathDist: TImage
        Left = 0
        Top = 0
        Width = 692
        Height = 522
        Align = alClient
        ExplicitLeft = 128
        ExplicitTop = 104
        ExplicitWidth = 105
        ExplicitHeight = 105
      end
    end
  end
  object MainMenu1: TMainMenu
    Left = 744
    object fil1: TMenuItem
      Caption = 'file'
      object save1: TMenuItem
        Caption = 'save'
        OnClick = save1Click
      end
    end
  end
end


object Form3: TForm3
  Left = 0
  Top = 0
  Caption = 'Form3'
  ClientHeight = 202
  ClientWidth = 447
  Color = clBtnFace
  Font.Charset = DEFAULT_CHARSET
  Font.Color = clWindowText
  Font.Height = -11
  Font.Name = 'Tahoma'
  Font.Style = []
  OldCreateOrder = False
  OnClose = FormClose
  PixelsPerInch = 96
  TextHeight = 13
  object Panel1: TPanel
    Left = 16
    Top = 32
    Width = 321
    Height = 137
    TabOrder = 0
    object Label1: TLabel
      Left = 40
      Top = 16
      Width = 44
      Height = 13
      Caption = 'font size '
    end
    object SpinEdit1: TSpinEdit
      Left = 90
      Top = 13
      Width = 63
      Height = 22
      MaxValue = 0
      MinValue = 0
      TabOrder = 0
      Value = 12
    end
    object CheckBox1: TCheckBox
      Left = 56
      Top = 48
      Width = 97
      Height = 17
      Caption = 'raster visible'
      TabOrder = 1
    end
    object CheckBox2: TCheckBox
      Left = 56
      Top = 80
      Width = 97
      Height = 17
      Caption = 'model visible'
      Checked = True
      State = cbChecked
      TabOrder = 2
    end
  end
end


program Project2;

uses
  Vcl.Forms,
  Unit2 in 'Unit2.pas' {Form2},
  Unit1 in 'Unit1.pas' {Form1},
  Unit3 in 'Unit3.pas' {Form3},
  NelderMead in 'NelderMead.pas';

{$R *.res}

begin
  Application.Initialize;
  Application.MainFormOnTaskbar := True;
  Application.CreateForm(TForm2, Form2);
  Application.CreateForm(TForm1, Form1);
  Application.CreateForm(TForm3, Form3);
  Application.Run;
end.