unit MassScale;
interface

uses
     Peak,IniFiles;

const
     cNoiseScale=10000;

     cMinPeakCount=0;
     cMaxPeakCount=1000;

     cDefaultPeakCount=200;
     cDefaultMinMass=0;
     cDefaultMaxMass=400*cMassScale;
     cDefaultMinAmplitude=0;
     cDefaultMaxAmplitude=cSignalScale;
     cDefaultSigma=cMassScale div 100;
     cDefaultSigmaSigma=0;
     cDefaultNoise=round(cNoiseScale);

     cDefaultNameSpectrFile='Peaks.txt';
     cDefaultNamePeakFile='Peaks.ini';

     
     
type
     tPeakCount=cMinPeakCount..cMaxPeakCount;
     tPeaks=array[tPeakCount] of tPeak;
     tpPeaks=^tPeaks;

     tSignalFile=record
          m:tMass;
          s:tSignal;
     end;

     tSignalsFile=record
          Data: array of tSignalFile;
          Count:cardinal;
     end;

     tPeakFile=TIniFile;

     tMassScale=class
     private
          mutex          :THandle;

          SignalsFile    :tSignalsFile;
          PeakFile       :tPeakFile;

          lPeaks         :tpPeaks;
          PeakCount      :tPeakCount;
          AllocatePeak   :Cardinal;
          MinMass        :tMass;
          MaxMass        :tMass;
          Sigma          :tMass;
          SigmaSigma     :tMass;
          MinAmplitude   :tSignal;
          MaxAmplitude   :tSignal;

          Noise          :tSignal;

          function Signal0(m:tMass):tSignal;
          function Signal2(m:tMass):tSignal;
          procedure Init0;
          procedure Init2;
          procedure Init1;

          procedure Allocate;
          procedure UnAllocate;
          function Size:cardinal;
          procedure Init;
          function FindPeakPlace(m:tMass; MaxNum:tPeakCount):tPeakCount;
          procedure AddPeakToPlace(NumPeak:tPeakCount; MaxNum:tPeakCount);
          function FindNearestPeak(m:tMass):tPeakCount;

          function SignalFromFile(m:tMass):tSignal;
          function FindMass(m:tMass; MaxNum:cardinal):Cardinal;

          function TryLock:boolean;
          function Lock:boolean;
          procedure UnLock;
     public
          constructor Create;
          destructor Destroy; override;
          procedure ReInit;

          function Signal(m:tMass):tSignal;
          function GetNoise(Time:Cardinal;sgnl:tSignal):tSignal;
     end;

implementation

uses
     Windows, uMainForm, StdCtrls,
     sysutils,Classes,
     Hardware, Math;

//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
function XRnd(l,r:cardinal):cardinal;
var
     tmp:cardinal;
begin
     if l>r then begin
          tmp:=l; l:=r; r:=tmp;
     end;
     if l=r then begin
          Result:=l;
       end else begin
          Result:=Random(r-l)+l;
     end;
end;


//*****************************************************************************
//
//*****************************************************************************


//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
constructor tMassScale.Create;
begin
     mutex:=CreateMutex(nil, false, nil);
     AllocatePeak:=0;
     MinMass:=Params.MinMass;
     MaxMass:=Params.MaxMass;
     Sigma:=Params.Sigma;
     SigmaSigma:=Params.SigmaSigma;
     MinAmplitude:=Params.MinAmplitude;
     MaxAmplitude:=Params.MaxAmplitude;


     SetLength(SignalsFile.Data ,0);
     SignalsFile.Count:=0;
     UnLock;

     Init;
end;



//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
destructor tMassScale.Destroy;
begin
     UnAllocate;
end;



//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
function tMassScale.Size:cardinal;
begin
     Result:=PeakCount*SizeOf(tPeak);
end;



//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
procedure tMassScale.Allocate;
begin
     if Assigned(lPeaks) then begin
          UnAllocate;
     end;
     if PeakCount>0 then begin
          AllocatePeak:=Size;
          GetMem(lPeaks,AllocatePeak);
     end;
end;



//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
procedure tMassScale.UnAllocate;
begin
     if Assigned(lPeaks) then begin
          FreeMem(lPeaks,AllocatePeak);
          AllocatePeak:=0;
     end;
     lPeaks:=nil;
     if SignalsFile.Count>0 then begin
          SetLength(SignalsFile.Data,0);
     end;
end;



//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
procedure tMassScale.ReInit;
begin
     UnAllocate;
     Init;
end;



//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
procedure tMassScale.Init;
begin
if Lock then begin
                         // в зависимости от выбранного режима генерации пиков
     case Params.ModeGenPeak of
                                   //случайная генерация спектра
          0:begin
               Init0;
          end;
                                   // чтение параметров пиков из файла
          1:begin
               Init1;
          end;
                                   // чтение реального пика из файла
          2:begin
               Init2;
          end;
       else begin
          Init0;
          Params.ModeGenPeak:=0;
       end;
     end;
end;
UnLock;

end;


//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
function tMassScale.FindPeakPlace(m:tMass; MaxNum:tPeakCount):tPeakCount;
var
     i,il,ir:tPeakCount;
begin
     il:=Low(i);
     ir:=MaxNum;
     while (il<ir) do begin
          i:=(il+ir) div 2;
          if lPeaks^[i].CenterPeak >m then begin
               ir:=i;
            end else begin
               il:=Succ(i);
          end;
     end;
     Result:=ir;
end;


//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
procedure tMassScale.AddPeakToPlace(NumPeak:tPeakCount; MaxNum:tPeakCount);
begin
     if NumPeak<MaxNum then begin
          System.Move(lPeaks^[NumPeak],
		     lPeaks^[Succ(NumPeak)],
		     SizeOf(lPeaks^[NumPeak])*(MaxNum-NumPeak)
          );
     end;

end;



//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
function tMassScale.Signal(m:tMass):tSignal;
begin
     if Lock then begin
          case Params.ModeGenPeak of
               0,1:begin
                    Result:=Signal0(m);
               end;
               2:begin
                    Result:=Signal2(m);
               end;
          else
               Result:=0;
          end;
     end;
     UnLock;
end;



//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
function tMassScale.FindNearestPeak(m:tMass):tPeakCount;
begin
  if (m<lPeaks^[0].CenterPeak) then begin
    Result:=Low(tPeakCount);
  end else if (m>lPeaks^[Pred(PeakCount)].CenterPeak) then begin
    Result:=Pred(PeakCount);
  end else begin
    Result:=FindPeakPlace(m, Pred(PeakCount));
  end;
end;





//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
function tMassScale.FindMass(m:tMass; MaxNum:cardinal):Cardinal;
var
     i,il,ir:cardinal;
begin
     il:=0;
     ir:=MaxNum;
     while (il<ir) do begin
          i:=(il+ir) div 2;
          if SignalsFile.Data[i].m >m then begin
               ir:=i;
            end else begin
               il:=Succ(i);
          end;
     end;
     Result:=ir;
end;



//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
function tMassScale.SignalFromFile(m:tMass) :tSignal;
var
     lm:tMass;
     x1,x2:tMass;
     y1,y2:Extended;
begin
     lm:=FindMass(m, SignalsFile.Count);
     if lm>0 then begin
          x1:=SignalsFile.Data[Round(lm)-1].m;
          y1:=SignalsFile.Data[Round(lm)-1].s;
       end else begin
          x1:=m;
          y1:=0;
     end;
     if lm=SignalsFile.Count-1 then begin
          Result:=0;
       end else begin
          Result:=(m-x1)/(x2-x1)*(y2-y1)+y1;
     end;
end;




//-----------------------------------------------------------------------------
// подготовка данных спектра при выборе режима генерации спектра случайным образом
//-----------------------------------------------------------------------------
procedure tMassScale.Init0;
var
     i:tPeakCount;
     m:tMass;
     NumPeak:tPeakCount;
     f:text;
begin
     PeakCount:=Params.PeakCount;
     Allocate;
     Noise:=Params.Noise0;
     if Assigned(lPeaks) and (PeakCount>0) then begin
          for i:=Low(i) to Pred(Low(i)+PeakCount) do begin
               m:=xRnd(Params.MinMass,Params.MaxMass);
               NumPeak:=FindPeakPlace(m,i);
               AddPeakToPlace(NumPeak,i);
               lPeaks^[NumPeak]:=(tPeak.Create);
               with lPeaks^[NumPeak] do begin
                    CenterPeak:=m;
                    AmplitudePeak:=XRnd(round(Params.MinAmplitude),round(Params.MaxAmplitude));
                    if Params.SigmaSigma>0 then
                         SigmaPeak:=Params.Sigma+(Random(2*Params.SigmaSigma)-Params.SigmaSigma)
                    else
                         SigmaPeak:=Params.Sigma;
               end;
          end;
     end;
end;


//-----------------------------------------------------------------------------
// подготовка данных спектра при выборе режима чтени параметров пиков из файла
//-----------------------------------------------------------------------------
procedure tMassScale.Init1;
var
     pPeaksRead:tStrings;
     pPeaks:tStrings;
     i:integer;
     m:tMass;
     NumPeak:tPeakCount;
     tmp:string;
     t:integer;
begin
 try
     pPeaksRead:=tStringList.Create;
     pPeaks:=tStringList.Create;
//     pPeaksRead.Items.Clear;
//     pPeaks.Clear;
     If FileExists(Params.NamePeakFile) then begin
          Noise:=Params.Noise1;
          PeakCount:=0;
          PeakFile:=tPeakFile.Create(Params.NamePeakFile);
          PeakFile.ReadSections(pPeaksRead);
          for i:=0 to pPeaksRead.Count-1 do begin
               if Pos('PeakNum',pPeaksRead[i])=1 then begin
                    Inc(PeakCount);
                    pPeaks.Add(pPeaksRead.Strings[i]);
               end;
          end;
          Allocate;
          for i:=0 to PeakCount-1 do begin
               m:=PeakFile.ReadInteger(pPeaks.Strings[i],'Mass',xRnd(cDefaultMinMass,cDefaultMaxMass));
               NumPeak:=FindPeakPlace(m,i);
               AddPeakToPlace(NumPeak,i);
               lPeaks^[NumPeak]:=(tPeak.Create);
               with lPeaks^[NumPeak] do begin
                    CenterPeak:=m;
                    AmplitudePeak:=PeakFile.ReadInteger(pPeaks.Strings[i], 'Amplitude', Round(XRnd(Round(cDefaultMinAmplitude),round(cDefaultMaxAmplitude))));
                    SigmaPeak:=PeakFile.ReadInteger(pPeaks.Strings[i], 'Sigma', cDefaultSigma);
               end;
          end;
          PeakFile.Free;
       end else begin
          Params.ModeGenPeak:=0;
          Init0;
     end;
 except
     Params.ModeGenPeak:=0;
     Init0;
 end;
     pPeaksRead.Free;
     pPeaks.Free;
end;



//-----------------------------------------------------------------------------
// подготовка данных спектра при выборе режима генерации спекра чтением реального из файла
//-----------------------------------------------------------------------------
procedure tMassScale.Init2;
var
     F:text;
     tmp:Extended;
begin
 try
     if FileExists(Params.NameSpectrFile) then begin
          AssignFile(F,Params.NameSpectrFile);
          Reset(f);
          while not(Eof(f)) do begin
               Inc(SignalsFile.Count);
               SetLength(SignalsFile.Data, SignalsFile.Count);
               Read(F,tmp);
               SignalsFile.Data[SignalsFile.Count-1].m:=Round(tmp*cMassScale);
               Read(F,tmp);
               SignalsFile.Data[SignalsFile.Count-1].s:=tmp;
               if (SignalsFile.Data[SignalsFile.Count-1].m=0) and (SignalsFile.Data[SignalsFile.Count-1].s=0) then begin
                    Dec(SignalsFile.Count);
                    SetLength(SignalsFile.Data, SignalsFile.Count);
               end;
          end;
          tmp:=1;
      end else begin
          Params.ModeGenPeak:=0;
          Init0;
     end;
 except
     Params.ModeGenPeak:=0;
     Init0;
 end;
end;




//-----------------------------------------------------------------------------
// сигнал в случае режима генерации спектра случайным образом
//-----------------------------------------------------------------------------
function tMassScale.Signal0(m:tMass):tSignal;
var
     i:tPeakCount;
begin
     i:=FindNearestPeak(m);
     Result:=lPeaks^[i].Signal(m);
     if i<(PeakCount-1) then begin
          Result:=Result+lPeaks^[i+1].Signal(m);
     end;
     if i>0 then begin
          Result:=Result+lPeaks^[i-1].Signal(m);
     end;
end;


//-----------------------------------------------------------------------------
// сигнал в случае режима генерации спектра чтением реального спектра из файла
//-----------------------------------------------------------------------------
function tMassScale.Signal2(m:tMass):tSignal;
var
     lm:tMass;
     t:tMass;
     x1,x2:tMass;
     y1,y2:Extended;
     f:textfile;
begin

     if (m<SignalsFile.Data[0].m) or (m>SignalsFile.Data[SignalsFile.Count-1].m) then begin
          Result:=0;
          exit;
     end;

     lm:=FindMass(m, SignalsFile.Count);
     if lm>0 then begin
          x1:=SignalsFile.Data[lm-1].m;
          y1:=SignalsFile.Data[lm-1].s;
       end else begin
          x1:=m;
          y1:=0;
     end;
     x2:=SignalsFile.Data[lm].m;
     y2:=SignalsFile.Data[lm].s;
     if lm=(SignalsFile.Count-1) then begin
          Result:=0;
       end else begin
          Result:=(m-x1)/(x2-x1)*(y2-y1)+y1;
     end;
end;


//-----------------------------------------------------------------------------
function tMassScale.TryLock:boolean;
begin
     Result:=WaitForSingleObject(mutex,3000)=Wait_OBJECT_0;

end;

//-----------------------------------------------------------------------------
function tMassScale.Lock:boolean;
begin
     Result:=WaitForSingleObject(mutex,infinite)=Wait_OBJECT_0;

end;

//-----------------------------------------------------------------------------
procedure tMassScale.UnLock;
begin
     ReleaseMutex(mutex)
end;

//-----------------------------------------------------------------------------
function tMassScale.GetNoise(Time:Cardinal;sgnl:tSignal):tSignal;
begin
     if Params.ModeGenPeak=2 then
          Result:=0
     else
//          Result:=((Noise/cNoiseScale)*random)/sqrt(Time);
          Result:=0.02+((Noise/cNoiseScale)*random)/sqrt(Time);
end;


end.