unarmath.pas

神经网络用于X荧光分析的源码。输入元素峰的强度

    Y : Double;
begin
  Result := True;
  Step := (StopY - StartY) / ChanCount;
  Y := StartY + Step / 2;
  I := 0;
  repeat
    DistrX[I] := 0;
    DistrY[I] := Y + Step / 2;
    for J := StartPosition to StopPosition do
      if (ResidsValues[J] >= Y) and (ResidsValues[J] < Y + Step) then
        DistrX[I] := DistrX[I] + 1;
    Y := Y + Step;
    Inc(I);
  until (Y >= StopY) or (I >= ChanCount);
end;

//------------------------------------------------------------------------------
// Calculate Max of float array (returned value) in StartPosition - StopPosition
// additional is position of Max in array calculated
// by equal values the lowest position will be prefeared.
function MaxFloatArray(const StartPosition, StopPosition : LongInt; // > 0
                       var MaxPosition                   : LongInt;
                       const FloatArray                  : array of Double //start by 0
                       )                                 : Double;
var I : Word; // counter
begin
  if (Length(FloatArray) < StopPosition)
  then begin
    Result := 0;
    MessageBeep(0);
    Exit;
  end;

  Result := FloatArray[StartPosition];
  MaxPosition := StartPosition;
  for I := StartPosition + 1 to StopPosition do
  begin
    if FloatArray[I] > Result then
    begin
      Result := FloatArray[I];
      MaxPosition := I;
    end
  end;
end;

//------------------------------------------------------------------------------
// Calculate Min of float array (returned value) in StartPosition - StopPosition
// additional is position of Min in array calculated
// by equal values the lowest position will be prefeared.
function MinFloatArray(const StartPosition, StopPosition : LongInt; // > 0
                       var MaxPosition                   : LongInt;
                       const FloatArray                  : array of Double //start by 0
                       )                                 : Double;
var I : Word; // counter
begin
  if (Length(FloatArray) < StopPosition)
  then begin
    Result := 0;
    MessageBeep(0);
    Exit;
  end;

  Result := FloatArray[StartPosition];
  MaxPosition := StartPosition;
  for I := StartPosition + 1 to StopPosition do
  begin
    if FloatArray[I] < Result then
    begin
      Result := FloatArray[I];
      MaxPosition := I;
    end
  end;
end;

//------------------------------------------------------------------------------
// Weight factors
function CalcWeigtArray(const PoisonWeighted : Boolean;
                        const StartPosition, StopPosition : LongInt; // > 0
                        const MeasValues_Y  : array of Double;
                        var   WeightArray : array of Double
                        ) : Boolean;
var I : LongInt; //counter
begin
  Result := True;
  // Verify dimensions
  if (Length(MeasValues_Y) < StopPosition + 1)
  or (Length(WeightArray) < StopPosition + 1)
  then begin
    Result := False;
    MessageBeep(0);
    Exit;
  end;

  for I := StartPosition to StopPosition do
    WeightArray[I] := 1;

  if PoisonWeighted then
    for I := StartPosition to StopPosition do
      if MeasValues_Y[I] <> 0 then
        WeightArray[I] := 1.0 / MeasValues_Y[I]
//      else WeightArray[I] := 0.0; //?????????????
end;

//------------------------------------------------------------------------------

function Standard_Deviation;
const MaxValue : Extended = 1.0E300; //double limit
var   I        : Word;
      Resids   : array of Double;
      ChannelCount : LongInt;
begin
  // Verify dimensions
  ChannelCount := StopPosition - StartPosition + 1;
  if (Length(MeasValues_Y)  < StopPosition + 1)
  or (Length(FittedValues) < StopPosition + 1 )
  or (Length(WeightArray)  < StopPosition + 1 )
  or (ChannelCount <= 0)
  then begin
    Result := -1; //Array dimension error
    MessageBeep(0);
    Exit;
  end;

  SetLength(Resids,StopPosition + 1);

  try
    Result := 0;
    for I := StartPosition to StopPosition do
      Resids[I]:= FittedValues[I] - MeasValues_Y[I];
    for I := StartPosition to StopPosition do
    begin
      Result := Result + Resids[I] * Resids[I] * WeightArray[I];
      if ABS(Result) > MaxValue then
      begin
        Result := MaxValue;
        Break;
      end;
    end;
    //Calculate normalized standard deviation
    if ChannelCount > FreedomDegree + 1 then
      Result := Result / (ChannelCount - FreedomDegree - 1);

//    Result := SQRT(Result);
  except
    Result := -2; //Math error
  end;
end;

//------------------------------------------------------------------------------
// Calculats wighted residuals
function Residuals
                  (const PoisonWeighted              : Boolean;
                   const StartPosition, StopPosition : LongInt; // > 0
                   const MeasValues_Y, FittedValues    : array of Double;
                   var ResidsValues                  : array of Double
                    )                                : Boolean;
var I : Word; // counter
    WeightArray : array of Double;
//    MaxPos : LongInt;
begin
  Result := True;
  // Verify dimensions
  if (Length(MeasValues_Y) < StopPosition  + 1)
  or (Length(FittedValues) < StopPosition + 1 )
  or (Length(ResidsValues) < StopPosition + 1 )
  then begin
    Result := False;
    MessageBeep(0);
    Exit;
  end;

  SetLength(WeightArray,StopPosition + 1);

  try
    // Set weight
    if not CalcWeigtArray(PoisonWeighted,StartPosition,StopPosition,
           MeasValues_Y,WeightArray)
    then begin
      Result := False;
      MessageBeep(0);      
      Exit;
    end;
    // Residuals
    for I := StartPosition to StopPosition do
      ResidsValues[I] := (MeasValues_Y[I] - FittedValues[I]) * Power(WeightArray[I],0.5);

  except
    Result := False;
  end;
end;

//------------------------------------------------------------------------------

function Autocorrelation
                  (const PoisonWeighted              : Boolean;
                   const StartPosition, StopPosition : LongInt; // > 0
                   const MeasValues_Y, FittedValues    : array of Double;
                   var   AutoCorrValues              : array of Double;
                   var   Count                       : LongInt
                   )                                 : Boolean;
var I,J,M,N        : LongInt; // counter
    WeightArray    : array of Double;
    MyBuff         : array of Double;
//    MaxPos         : LongInt;
    Sum, Buff      : Double;
begin
//Muss geprueft werden!!!!!!!!!!!!!!!
  Result := True;
  // Count
  N :=  StopPosition - StartPosition + 1;
  Count := N div 2 - 1;
  M := N - Count;
  // Verify dimensions
  if Count <= 1 then //??
  begin
    Result := False;
    MessageBeep(0);
    Exit;
  end;

  if (Length(MeasValues_Y) < StopPosition  + 1)
  or (Length(FittedValues) < StopPosition + 1)
  or (Length(AutoCorrValues) < StopPosition + 1)
  then begin
    Result := False;
    MessageBeep(0);
    Exit;
  end;

  SetLength(WeightArray,StopPosition + 1);
  SetLength(MyBuff,StopPosition + 1);

  try
    // Set weight
    if not CalcWeigtArray(PoisonWeighted,StartPosition,StopPosition,
           MeasValues_Y,WeightArray)
    then begin
      Result := False;
      MessageBeep(0);
      Exit;
    end;

    // Autocorrelation
    Sum := 0;
    J := 0;

    for I := StartPosition to StopPosition do
    begin
       MyBuff[j] := MeasValues_Y[I] - FittedValues[I]; //???
       Sum := Sum + WeightArray[I] * MyBuff[J] * MyBuff[J];
       Inc(J);
    end;

    Buff := Sum / N;

    for J := 0 to Count - 1 do
    begin
      Sum := 0;
      for I := 1 to M do
        Sum  := Sum + MyBuff[I] * MyBuff[I + J] *
                Sqrt(WeightArray[I + StartPosition])*
                Sqrt(WeightArray[I + J + StartPosition]);
      if M <> 0 then Sum := Sum / M;

      if Buff <> 0 then AutoCorrValues[J] := Sum / Buff
                   else AutoCorrValues[J] := 0;
     end;

  except
    Result := False;
  end;

end;

//------------------------------------------------------------------------------
//                     FITTING
//------------------------------------------------------------------------------

// Back-Matrix, by Gaus-Jordan Method from 'Numerical Resepies in Pascal...'
//------------------------------------------------------------------------------

function Gauss_Jordan(
                      N               : Byte;         //Dimension
                      var Matrix      : TFloatMatrix; //
                      var Vector      : TFloatVector
                      )               : LongInt;
var
  big,dum,pivinv        : Double;
  I,Icol,Irow,J,K,L,LL  : Integer;
  Indxc,Indxr,Ipiv      : array of LongInt;
begin
   Result := 0;
   SetLength(indxc,N);
   SetLength(indxr,N);
   SetLength(ipiv,N);
   SetLength(Matrix,N,N); //????????????????????????????
   SetLength(Vector,N);
   try
     Icol := 0;
     Irow := 0;
     for J := 0 to N - 1 do ipiv[J] := 0;
     for I := 0 to N - 1 do
     begin
       Big := 0.0;
       for J := 0 to N - 1 do
       begin
         if ipiv[J] <> 1 then
         begin
           for K := 0 to N - 1 do
           begin
             if ipiv[K] = 0 then
             begin
               if abs(Matrix[J,K]) >= big then
               begin
                  big := abs(Matrix[J,K]);
                  irow := J;
                  icol := K;
               end
               else
                  if ipiv[K] > 1 then Inc(Result);
             end;
           end;
         end;
       end;

       ipiv[icol] := ipiv[icol] + 1;

       if irow <> icol then
       begin
         for L := 0 to N - 1 do
         begin
            dum := Matrix[irow,L];
            Matrix[icol,L] := dum;
          end;
          dum := Vector[irow]; //
          Vector[irow] :=  Vector[icol];
          Vector[icol] := dum;
       end;
       indxr[I] := irow;
       indxc[I] := icol;
       if Matrix[icol,icol] = 0.0 then
       begin
          Inc(Result);
          Matrix[icol,icol] := 1.;
       end;
       //start AR
       if Matrix[icol,icol] < 1E-200 then Matrix[icol,icol] := 1E-200;
       if Matrix[icol,icol] > 1E200 then Matrix[icol,icol] := 1E200;
       //end AR
       pivinv := 1.0 / Matrix[icol,icol];
       Matrix[icol,icol] := 1.0;
       for L := 0 to N - 1 do
         Matrix[icol,L] := Matrix[icol,L] * pivinv;
       Vector[icol] := Vector[icol] * pivinv;
       for LL := 0 to N - 1 do
       begin
         if LL <> icol then
         begin
           dum := Matrix[LL,icol];
           Matrix[LL,icol] := 0.0;
           for L := 0 to N - 1 do Matrix[LL,L] := Matrix[LL,L] - Matrix[icol,L] * dum;
           Vector[LL] := Vector[LL] - Vector[icol] * dum;
         end;
       end;
     end;

     for L := N - 1 downto 0 do
     begin
       if indxr[L] <> indxc[L] then
       begin
         for K := 0 to N - 1 do
         begin
            dum := Matrix[K,indxr[L]];
            Matrix[K,indxr[L]] := Matrix[K,indxc[L]];
            Matrix[K,indxc[L]] := dum;
         end;
       end;
     end;

   except
     Inc(Result);
   end;
end;

PROCEDURE Fitting;
{////////////////////////////////////////////////////////////////////}
{//////              MARQUARDT-LEVENBERG ALGORITHM         //////////}
{////////////////////////////////////////////////////////////////////}
const
  StepDevider = 10;
  Tau         = 0.001; //zu vermeiden divisiov by zero
var
  WorkArr         : array of array of Double; //deriv. matrix
  WorkVect        : array of Double; //derivatives save
  WeightArray     : array of Double;
  A,C             : TFloatMatrix; //matrix calculation Gauss-Jordan
  Save_XX,                  // var parametres, save
  Diag, G,
  DeltaX          : TFloatVector; //for fitted params
  J               : Word; //counter
  Save_Step,
  F,F1,F2,F3,
  Save_StdDev     : Double; //StdDev save variable
  W               : Double; //power, step by iter type iii
  Step            : Double;
//------------------------------------------------------------------------------
// Fitting function as parameter
function Fit_Function : Double;
var I : LongInt;
begin
  Result := 0;

  case Fitting_Formula_No of
  1 :  if not Gauss(StartPosition, StopPosition, FittedParams, FittedValues) then
       begin
          Inc(error_Count);