InraPorc

unit UFCalibrage;

interface

uses

  Windows, Forms, Classes, Controls, StdCtrls, Buttons, dmath;

type

  TFCalibrage = class(TForm)

    BBOk: TBitBtn;

    LTitle: TLabel;

    LMethod: TLabel;

    LTime: TLabel;

    ETime: TEdit;

    LPV: TLabel;

    LP2: TLabel;

    EPVRSD: TEdit;

    EP2RSD: TEdit;

    EPVR2: TEdit;

    LInitial: TLabel;

    LFinal: TLabel;

    ECoefEntGestInitial: TEdit;

    LStDev: TLabel;

    LCoefEntGest: TLabel;

    LCoefEntLact: TLabel;

    LRatioLipProt: TLabel;

    LInfo: TLabel;

    ECoefEntLactInitial: TEdit;

    ERatioLipProtInitial: TEdit;

    ECoefEntGestFinal: TEdit;

    ECoefEntLactFinal: TEdit;

    ERatioLipProtFinal: TEdit;

    ECoefEntGestStDev: TEdit;

    ECoefEntLactStDev: TEdit;

    ERatioLipProtStDev: TEdit;

    EMethod: TEdit;

    GBParam: TGroupBox;

    GBVariables: TGroupBox;

    EP2R2: TEdit;

    EPVVar: TEdit;

    EP2Var: TEdit;

    LR2: TLabel;

    LRSD: TLabel;

    LVar: TLabel;

    Lw: TLabel;

    EwP2: TEdit;

    EwPds: TEdit;

    procedure FormCreate(Sender: TObject);

    procedure FormDestroy(Sender: TObject);

    procedure FormShow(Sender: TObject);

  private

    { D?clarations priv?es }

    Pds, P2: TVector;

    X, G: TVector;

    H_inv: TMatrix;

    F_min, Det: Float;

    Freq, Start, Stop: Int64;

    MPds, MP2: Extended;

  public

    { D?clarations publiques }

  end;

var

  FCalibrage: TFCalibrage;

implementation

uses

  Graphics, Math, SysUtils, gnugettext, UVariables, UStrings, UFProfilT,

  UCalcSimulT;

{$R *.dfm}

{ Calibration }

const

  NVar = 3;

  NData = NB_CYCLES + NB_CYCLES - 1;

  MaxIter = 100;

  Tol = 1.0E-6;

var

  wPds, wP2: Extended;

function Func(X: TVector): Float;

var

  cycle, jour: Integer;

  GainProt: double;

  PdsSail, P2Sail, PdsMB, P2MB: array[1..NB_CYCLES] of double;

begin

  PProfilT.CoefEntretienGest := X[1];

  PProfilT.CoefEntretienLact := X[2];

  PProfilT.RatioLipProt := X[3];

  // calcul des CoefNR

  for cycle := 1 to NB_CYCLES do

    PProfilT.CoefNR[cycle] := 1;

  CalcSimulT(-1, -1, 1, 1, 8, 0, PResSimulT);

  GainProt := 0;

  for jour := 1 to PResSimulT.NbJSim do

    if PResSimulT.TabResult[3, jour] = 0

    then // Gestation

      GainProt := GainProt + PResSimulT.TabResult[63, jour]

    else

      if (PResSimulT.TabResult[3, jour] = 1) and (PResSimulT.TabResult[3, jour - 1] = 0)

      then // Mise bas

      begin

        cycle := Round(PResSimulT.TabResult[2, jour]);

        with PProfilT.Truies[cycle] do

          PProfilT.CoefNR[cycle] := Math.Max((0.178 * (PdsApMB - PdsSail) - 0.333 * (P2MB - P2Sail)) / (GainProt / 1000), 0);

        GainProt := 0;

      end;

  // ?valuation du cas

  CalcSimulT(-1, -1, 1, 1, 8, 0, PResSimulT);

  FProfilT.PBCoefEntGest.AsFloat := PProfilT.CoefEntretienGest;

  FProfilT.PBCoefEntLact.AsFloat := PProfilT.CoefEntretienLact;

  FProfilT.PBRatioLipProt.AsFloat := 1 - PProfilT.RatioLipProt;

  FProfilT.AffGraphCal;

  if FProfilT.GaugeCalibr.Progress < FProfilT.GaugeCalibr.MaxValue

  then

    FProfilT.GaugeCalibr.AddProgress(1)

  else

    FProfilT.GaugeCalibr.Progress := 1;

  Application.ProcessMessages;

  // poids et ?paisseur de lard apr?s gestation (mise bas)

  for jour := 2 to PResSimulT.NbJSim do

    if (PResSimulT.TabResult[3, jour] = 1) and (PResSimulT.TabResult[3, jour - 1] = 0)

    then

    begin

      cycle := Round(PResSimulT.TabResult[2, jour]);

      PdsMB[cycle] := PResSimulT.TabResult[60, jour - 1] + PResSimulT.TabResult[61, jour];

      P2MB[cycle] := PResSimulT.TabResult[67, jour - 1] + PResSimulT.TabResult[68, jour];

    end ;

  // poids et ?paisseur de lard apr?s lactation (saillie)

  for jour := 2 to PResSimulT.NbJSim do

    if (PResSimulT.TabResult[3, jour] = 0) and (PResSimulT.TabResult[3, jour - 1] = 2)

    then

      begin

        cycle := Round(PResSimulT.TabResult[2, jour]);

        PdsSail[cycle] := PResSimulT.TabResult[60, jour - 1] + PResSimulT.TabResult[61, jour];

        P2Sail[cycle] := PResSimulT.TabResult[67, jour - 1] + PResSimulT.TabResult[68, jour];

      end;

  // ?cart par rapport aux valeurs du profil

  Result := 0;

  for cycle := 1 to NB_CYCLES do

    Result := Result + Power(PdsMB[cycle] - PProfilT.Truies[cycle].PdsApMB, 2) * wPds;

  for cycle := 2 to NB_CYCLES do

    Result := Result + Power(PdsSail[cycle] - PProfilT.Truies[cycle].PdsSail, 2) * wPds;

  for cycle := 1 to NB_CYCLES do

    Result := Result + Power(P2MB[cycle] - PProfilT.Truies[cycle].P2MB, 2) * wP2;

  for cycle := 2 to NB_CYCLES do

    Result := Result + Power(P2Sail[cycle] - PProfilT.Truies[cycle].P2Sail, 2) * wP2;

end;

// Newton-Raphson method

// procedure HessGrad to compute the gradient G and the hessian H of the function at point X

{$I numhess.inc}

// Broyden-Fletcher-Goldfarb-Shanno method

// procedure Gradient to compute the gradient G of the function at point X

{$I numgrad.inc}

{ TFCalibrage }

procedure TFCalibrage.FormCreate(Sender: TObject);

var

  cycle: Integer;

  ETPds, ETP2: Extended;

begin

  if Screen.Fonts.IndexOf('Arial Unicode MS') <> -1

  then

    Font.Name := 'Arial Unicode MS';

  TranslateComponent(Self);

  // Initialisation des structures

  New(PResSimulT);

  DimVector(X, NVar);

  X[1] := PProfilT.CoefEntretienGest;

  X[2] := PProfilT.CoefEntretienLact;

  X[3] := PProfilT.RatioLipProt;

  DimVector(G, NVar);

  DimMatrix(H_inv, NVar, NVar);

  // Poids : moyenne et ?cart type

  DimVector(Pds, 2 * NB_CYCLES);

  for cycle := 2 to NB_CYCLES do

    Pds[cycle] := PProfilT.Truies[cycle].PdsSail;

  for cycle := 1 to NB_CYCLES do

    Pds[NB_CYCLES + cycle] := PProfilT.Truies[cycle].PdsApMB;

  MPds := dmath.Mean(Pds, 1, NData);

  ETPds := dmath.StDev(Pds, 1, NData, MPds);

  // Epaisseurs de lard : moyenne et ?cart type

  DimVector(P2, 2 * NB_CYCLES);

  for cycle := 2 to NB_CYCLES do

    P2[cycle] := PProfilT.Truies[cycle].P2Sail;

  for cycle := 1 to NB_CYCLES do

    P2[NB_CYCLES + cycle] := PProfilT.Truies[cycle].P2MB;

  MP2 := dmath.Mean(P2, 1, NData);

  ETP2 := dmath.StDev(P2, 1, NData, MP2);

  // Pond?rations

  wPds := 1 ;

  wP2 := Power(ETPds, 2) / Power(ETP2, 2);

//  wP2 := StrToFloat(InputBox('Question', 'Pond?ration relative', FloatToStrF(wP2, ffFixed, 15, 2)));

end;

procedure TFCalibrage.FormDestroy(Sender: TObject);

begin

  Dispose(PResSimulT);

end;

procedure TFCalibrage.FormShow(Sender: TObject);

var

  cycle, jour: Integer;

  PdsSail, P2Sail, PdsMB, P2MB: array[1..NB_CYCLES] of double;

  {MPds,} SStPds, SSrPds, {MP2,} SStP2, SSrP2: Float;

  {wM, wSStPds, wSSePds,} wSSrPds, {wSStP2, wSSeP2,} wSSrP2: Float;

begin

  ECoefEntGestInitial.Text := FProfilT.PBCoefEntGest.Text;

  ECoefEntLactInitial.Text := FProfilT.PBCoefEntLact.Text;

  ERatioLipProtInitial.Text := FProfilT.PBRatioLipProt.Text;

//  Screen.Cursor := crHourGlass;

  FProfilT.GaugeCalibr.Progress := 0;

  FProfilT.GaugeCalibr.Left := (FProfilT.GraphCalibr.Width - FProfilT.GaugeCalibr.Width) div 2;

  FProfilT.GaugeCalibr.Visible := True;

  Application.ProcessMessages;

  try

    if not QueryPerformanceFrequency(Freq) then Freq := 1;

    QueryPerformanceCounter(Start);

    case FProfilT.CBMethod.ItemIndex of

      0: Newton(Func, HessGrad, X, 1, NVar, MaxIter, Tol, F_min, G, H_inv, Det);

      1: Marquardt(Func, HessGrad, X, 1, NVar, MaxIter, Tol, F_min, G, H_inv, Det);

      2: BFGS(Func, Gradient, X, 1, NVar, MaxIter, Tol, F_min, G, H_inv);

    end;

  finally

    QueryPerformanceCounter(Stop);

//    Screen.Cursor:= crDefault;

    FProfilT.GaugeCalibr.Visible := False;

    Application.ProcessMessages;

  end;

  // Affichage du r?sultat obtenu

  PProfilT.CoefEntretienGest := X[1];

  PProfilT.CoefEntretienLact := X[2];

  PProfilT.RatioLipProt := X[3];

  CalcSimulT(-1, -1, 1, 1, 8, 0, PResSimulT);

  FProfilT.AffGraphCal ;

  FProfilT.PBCoefEntGest.AsFloat := PProfilT.CoefEntretienGest;

  FProfilT.PBCoefEntLact.AsFloat := PProfilT.CoefEntretienLact;

  FProfilT.PBRatioLipProt.AsFloat := 1 -PProfilT.RatioLipProt ;

  Application.ProcessMessages;

  EMethod.Text := FProfilT.CBMethod.Text;

  ETime.Text := FloatToStrF((Stop - Start) / Freq, ffFixed, 15, 3);

  ECoefEntGestFinal.Text := FProfilT.PBCoefEntGest.Text;

  if (FProfilT.PBCoefEntGest.AsFloat < FProfilT.PBCoefEntGest.MinValue)

  or (FProfilT.PBCoefEntGest.AsFloat > FProfilT.PBCoefEntGest.MaxValue)

  then

    ECoefEntGestFinal.Font.Color := clRed;

  ECoefEntLactFinal.Text := FProfilT.PBCoefEntLact.Text;

  if (FProfilT.PBCoefEntLact.AsFloat < FProfilT.PBCoefEntLact.MinValue)

  or (FProfilT.PBCoefEntLact.AsFloat > FProfilT.PBCoefEntLact.MaxValue)

  then

    ECoefEntLactFinal.Font.Color := clRed;

  ERatioLipProtFinal.Text := FProfilT.PBRatioLipProt.Text;

  if (FProfilT.PBRatioLipProt.AsFloat < FProfilT.PBRatioLipProt.MinValue)

  or (FProfilT.PBRatioLipProt.AsFloat > FProfilT.PBRatioLipProt.MaxValue)

  then

    ERatioLipProtFinal.Font.Color := clRed;

  // poids et ?paisseur de lard apr?s gestation (mise bas)

  for jour := 2 to PResSimulT.NbJSim do

    if (PResSimulT.TabResult[3, jour] = 1) and (PResSimulT.TabResult[3, jour - 1] = 0)

    then

    begin

      cycle := Round(PResSimulT.TabResult[2, jour]);

      PdsMB[cycle] := PResSimulT.TabResult[60, jour - 1] + PResSimulT.TabResult[61, jour];

      P2MB[cycle] := PResSimulT.TabResult[67, jour - 1] + PResSimulT.TabResult[68, jour];

    end ;

  // poids et ?paisseur de lard apr?s lactation (saillie)

  for jour := 2 to PResSimulT.NbJSim do

    if (PResSimulT.TabResult[3, jour] = 0) and (PResSimulT.TabResult[3, jour - 1] = 2)

    then

    begin

      cycle := Round(PResSimulT.TabResult[2, jour]);

      PdsSail[cycle] := PResSimulT.TabResult[60, jour - 1] + PResSimulT.TabResult[61, jour];

      P2Sail[cycle] := PResSimulT.TabResult[67, jour - 1] + PResSimulT.TabResult[68, jour];

    end;

{

  // Moyennes non pond?r?es

  MPds := 0;

  for cycle := 1 to NB_CYCLES do

    MPds := MPds + PProfilT.Truies[cycle].PdsApMB;

  for cycle := 2 to NB_CYCLES do

    MPds := MPds + PProfilT.Truies[cycle].PdsSail;

  MPds := MPds / NData;

  MP2 := 0;

  for cycle := 1 to NB_CYCLES do

    MP2 := MP2 + PProfilT.Truies[cycle].P2MB;

  for cycle := 2 to NB_CYCLES do

    MP2 := MP2 + PProfilT.Truies[cycle].P2Sail;

  MP2 := MP2 / NData;

}

  // Somme carr?e des ?carts totaux (non pond?r?s)

  SStPds := 0;

  for cycle := 1 to NB_CYCLES do

    SStPds := SStPds + Power(PProfilT.Truies[cycle].PdsApMB - MPds, 2);

  for cycle := 2 to NB_CYCLES do

    SStPds := SStPds + Power(PProfilT.Truies[cycle].PdsSail - MPds, 2);

  SStP2 := 0;

  for cycle := 1 to NB_CYCLES do

    SStP2 := SStP2 + Power(PProfilT.Truies[cycle].P2MB - MP2, 2);

  for cycle := 2 to NB_CYCLES do

    SStP2 := SStP2 + Power(PProfilT.Truies[cycle].P2Sail - MP2, 2);

  // Somme carr?e des ?carts des r?sidus (non pond?r?s)

  SSrPds := 0;

  for cycle := 1 to NB_CYCLES do

    SSrPds := SSrPds + Power(PProfilT.Truies[cycle].PdsApMB - PdsMB[cycle], 2);

  for cycle := 2 to NB_CYCLES do

    SSrPds := SSrPds + Power(PProfilT.Truies[cycle].PdsSail - PdsSail[cycle], 2);

  SSrP2 := 0;

  for cycle := 1 to NB_CYCLES do

    SSrP2 := SSrP2 + Power(PProfilT.Truies[cycle].P2MB - P2MB[cycle], 2);

  for cycle := 2 to NB_CYCLES do

    SSrP2 := SSrP2 + Power(PProfilT.Truies[cycle].P2Sail - P2Sail[cycle], 2);

 {

  // Moyenne pond?r?e

  wM := 0;

  for cycle := 1 to NB_CYCLES do

    wM := wM + PProfilT.Truies[cycle].PdsApMB * wPds;

  for cycle := 2 to NB_CYCLES do

    wM := wM + PProfilT.Truies[cycle].PdsSail * wPds;

  for cycle := 1 to NB_CYCLES do

    wM := wM + PProfilT.Truies[cycle].P2MB * wP2;

  for cycle := 2 to NB_CYCLES do

    wM := wM + PProfilT.Truies[cycle].P2Sail * wP2;

  wM := wM / (NData * wPds + NData * wP2);

  // Somme carr?e des ?carts totaux (pond?r?s)

  wSStPds := 0;

  for cycle := 1 to NB_CYCLES do

    wSStPds := wSStPds + Power(PProfilT.Truies[cycle].PdsApMB - wM, 2) * wPds;

  for cycle := 2 to NB_CYCLES do

    wSStPds := wSStPds + Power(PProfilT.Truies[cycle].PdsSail - wM, 2) * wPds;

  wSStP2 := 0;

  for cycle := 1 to NB_CYCLES do

    wSStP2 := wSStP2 + Power(PProfilT.Truies[cycle].P2MB - wM, 2) * wP2;

  for cycle := 2 to NB_CYCLES do

    wSStP2 := wSStP2 + Power(PProfilT.Truies[cycle].P2Sail - wM, 2) * wP2;

  // Somme carr?e des ?carts de pr?diction (pond?r?s)

  wSSePds := 0;

  for cycle := 1 to NB_CYCLES do

    wSSePds := wSSePds + Power(PdsMB[cycle] - wM, 2) * wPds;

  for cycle := 2 to NB_CYCLES do

    wSSePds := wSSePds + Power(PdsSail[cycle] - wM, 2) * wPds;

  wSSeP2 := 0;

  for cycle := 1 to NB_CYCLES do

    wSSeP2 := wSSeP2 + Power(P2MB[cycle] - wM, 2) * wP2;

  for cycle := 2 to NB_CYCLES do

    wSSeP2 := wSSeP2 + Power(P2Sail[cycle] - wM, 2) * wP2;

}

  // Somme carr?e des ?carts des r?sidus (pond?r?s)

  wSSrPds := 0;

  for cycle := 1 to NB_CYCLES do

    wSSrPds := wSSrPds + Power(PProfilT.Truies[cycle].PdsApMB - PdsMB[cycle], 2) * wPds;

  for cycle := 2 to NB_CYCLES do

    wSSrPds := wSSrPds + Power(PProfilT.Truies[cycle].PdsSail - PdsSail[cycle], 2) * wPds;

  wSSrP2 := 0;

  for cycle := 1 to NB_CYCLES do

    wSSrP2 := wSSrP2 + Power(PProfilT.Truies[cycle].P2MB - P2MB[cycle], 2) * wP2;

  for cycle := 2 to NB_CYCLES do

    wSSrP2 := wSSrP2 + Power(PProfilT.Truies[cycle].P2Sail - P2Sail[cycle], 2) * wP2;

  // Affichage des statistiques

  EwPds.Text := FloatToStrF(wPds, ffFixed, 15, 0);

  EwP2.Text := FloatToStrF(wP2, ffFixed, 15, 0);

  EPVRSD.Text := Format('%1.1f %s', [Sqrt(SSrPds / (NData - Nvar / 2)), StrKg]);

  EP2RSD.Text := Format('%1.2f %s', [Sqrt(SSrP2 / (NData - Nvar / 2)), StrMM]);

  EPVR2.Text := Format('%1.2f %%', [100 - SSrPds / SStPds * 100]);

  EP2R2.Text := Format('%1.2f %%', [100 - SSrP2 / SStP2 * 100]);

  EPVVar.Text := Format('%1.2f %%', [wSSrPds / (wSSrPds + wSSrP2) * 100]);

  EP2Var.Text := Format('%1.2f %%', [wSSrP2 / (wSSrPds + wSSrP2) * 100]);

  ECoefEntGestStDev.Text := FloatToStrF(Sqrt((wSSrPds + wSSrP2) / (2 * NData - Nvar) * H_inv[1, 1]), ffFixed, 15, 5);

  ECoefEntLactStDev.Text := FloatToStrF(Sqrt((wSSrPds + wSSrP2) / (2 * NData - Nvar) * H_inv[2, 2]), ffFixed, 15, 5);

  ERatioLipProtStDev.Text := FloatToStrF(Sqrt((wSSrPds + wSSrP2) / (2 * NData - Nvar) * H_inv[3, 3]), ffFixed, 15, 5);

  case MathErr of

    OptNonConv: LInfo.Caption := _('Non-convergence');

    OptSing: LInfo.Caption := _('Singular Hessian matrix');

    OptBigLambda: LInfo.Caption := _('Marquardt parameter too high');

    else

      if (FProfilT.PBCoefEntGest.AsFloat < FProfilT.PBCoefEntGest.MinValue)

      or (FProfilT.PBCoefEntGest.AsFloat > FProfilT.PBCoefEntGest.MaxValue)

      or (FProfilT.PBCoefEntLact.AsFloat < FProfilT.PBCoefEntLact.MinValue)

      or (FProfilT.PBCoefEntLact.AsFloat > FProfilT.PBCoefEntLact.MaxValue)

      or (FProfilT.PBRatioLipProt.AsFloat < FProfilT.PBRatioLipProt.MinValue)

      or (FProfilT.PBRatioLipProt.AsFloat > FProfilT.PBRatioLipProt.MaxValue)

      then

        LInfo.Caption := _('Suspicious result')

      else

        LInfo.Visible := False;

  end;

end;

end.