unit Hist;

interface

uses GLType, MathObj, Math,Classes,SysUtils, MVTU_TLB;

type
 //Свойства гистограммы
 THistProp=record
  Size:^integer;     //Размер серии
  CalcMode:^integer; //1-расчёт по сериям , 0-расчёт по всей выборке
  Tau:^double;       //Период квантования
  DelTrend:^integer; //Удаление линейного тренда 0-не удалять тренд, 1-удалять тренд
  Min:^double;       //Минимальное значение
  Max:^double;       //Максимальное значение
  Col:^integer;      //Число интервалов
  AutoRange:^integer;//Автоматическая установка границ интервалов в каждой серии 0-вручную, 1-автоматически
  OutMode:^integer;  //0-относительная частота попаданий
  end;               //1-число попаданий
                     //2-плотность вероятности

 PHistProp=^THistProp;

 THistNumericStates=record
  Position:integer;
  DMin,DMax,D0Min,D0Max,
  N,CN2,CNY,M,SmSqr:double;
  Auto:boolean;
  end;

 THistState=record
  Data,Sum,REZX,REZY:TExtArray;
  NumStates:THistNumericStates;
 end;
 PHistState=^THistState;

 function TStatHist(at,adt:RealType;var time,dtime:RealType;var AU,AY:TPtrExt;var AX,ADX:array of RealType;
   var CU,CY : TIntArray;var Prop:TPtrArray;var Vars : PHistState;RAction:Integer):Integer;export;
 function TStatHistRst(at,adt:RealType;var AU,AY:TPtrExt;var AX,ADX:PExtArr;
 var Prop:TPtrArray;var Vars : PHistState;Action:Integer;var F:File):Integer;export;

implementation

function TStatHistConvert(Action:  integer;
                          Vars:    Pointer;
                          Prop:    TPtrArray;
                          PropStr: TStringList;
                          BlockId: integer;
                          var MVTU: IMVTU_Server):integer;
  const RecName:PChar = 'Гистограмма распределения';
  var
      Res: integer;
begin
  case Action of
    //Возвращаем ссылку на имя записи в базе МВТУ-4
    cnv_GetRecName: Result:=integer(RecName);
    //Конвертируем и устанавливаем параметры при помощи интерфейса МВТУ
    cnv_Convert: with PHistProp(Prop.Arr)^ do try
                   //Присваиваем параметры для блока
                   MVTU.SetBlockProp(BlockId,'size',PChar(PropStr[0]),Res);
                   MVTU.SetBlockProp(BlockId,'calcmode',PChar(IntToStr(CalcMode^)),Res);
                   MVTU.SetBlockProp(BlockId,'tau',PChar(PropStr[2]),Res);
                   MVTU.SetBlockProp(BlockId,'deltrend',PChar(IntToStr(DelTrend^)),Res);
                   MVTU.SetBlockProp(BlockId,'min',PChar(PropStr[4]),Res);
                   MVTU.SetBlockProp(BlockId,'max',PChar(PropStr[5]),Res);
                   MVTU.SetBlockProp(BlockId,'col',PChar(PropStr[6]),Res);
                   MVTU.SetBlockProp(BlockId,'autorange',PChar(IntToStr(AutoRange^)),Res);
                   MVTU.SetBlockProp(BlockId,'outmode',PChar(IntToStr(OutMode^)),Res);
                 finally
                 end;
  end;
end;

function TStatHist;
 var
  i,j:integer;
  h,S,K,a,b:double;
 label L1,L2,L3,L4;
begin
 Result:=0;
 with PHistProp(Prop.Arr)^ ,Vars^,Vars^.NumStates do case RAction of
  f_GetConvertFuncAdr: Result:=integer(@TStatHistConvert);
  f_GetCount   :  begin
    CY.arr^[0]:=Col^;
    CY.arr^[1]:=Col^;
   end;
  f_InitTime   :  time:=at;
  f_InitState  :  begin
    Position:=0;
    DMin:=Min^;
    DMAx:=Max^;
    Auto:=(AutoRange^=1)and(CalcMode^=0);
    AY.Ptr(0).FillArray(0);
    AY.Ptr(1).FillArray(0);
    Sum.FillArray(0);
    N:=0;CN2:=0;CNY:=0;M:=0;SmSqr:=0;
    goto L1
   end;
  f_InitMem    :  begin
    Data.ChangeCount(Size^);
    Sum.ChangeCount(Col^);
    REZX.ChangeCount(Col^);
    REZY.ChangeCount(Col^);
   end;
  f_Create     :  begin
   Vars:=New(PHistState);
   with PHistState(Vars)^ do begin
    Data:=TExtArray.Create(1);
    Sum:=TExtArray.Create(1);
    REZX:=TExtArray.Create(1);
    REZY:=TExtArray.Create(1);
   end;
   end;
  f_Free       :  begin
    Data.Free;
    Sum.Free;
    REZX.Free;
    REZY.Free;
    Dispose(Vars)
   end;
  f_RestoreOuts: begin
    Move(REZX.arr^,AY.Ptr(0).arr^,REZX.Count*SOfR);
    Move(REZY.arr^,AY.Ptr(1).arr^,REZY.Count*SOfR);
   end;
  f_GoodStep   : if time-at<=adt/2 then begin
 L1:if CalcMode^=0 then begin
   if Auto then goto L2;
   N:=N+1;
   CN2:=CN2+N*N;
   CNY:=CNY+N*AU.Ptr(0).arr^[0];
   M:=M+AU.Ptr(0).arr^[0];
   if (DelTrend^=1) or (OutMode^=3) then begin
    if DelTrend^<>0 then b:=(CNY-0.5*M*N)/(CN2-N*sqr(0.5*N)) else b:=0;
    a:=M/N-b*0.5*N;
    S:=AU.Ptr(0).arr^[0]-(b*N+a);
   end else S:=AU.Ptr(0).arr^[0];
   h:=(DMax-DMin)/Col^;
   case OutMode^ of
    0: K:=1/N;
    2: K:=1/(N*h);
    3: begin
        SmSqr:=SmSqr+S*S;
        if N=1 then exit;
        K:=sqrt(SmSqr/(N-1));
        a:=1/K;
        K:=K/(N*h);
        for i:=0 to Col^-1 do AY.Ptr(0).arr^[i]:=(DMin+h*(i+0.5))*a;
        goto L3;
       end;
   else
    K:=1
   end;
   for i:=0 to Col^-1 do AY.Ptr(0).arr^[i]:=DMin+h*(i+0.5);
L3:if (S>=DMin)and(S<=DMax) then begin
    if S=DMax then i:=Col^-1 else i:=trunc((S-DMin)/h);
    Sum.arr^[i]:=Sum.arr^[i]+1;
   end;
   for j:=0 to Col^-1 do AY.Ptr(1).arr^[j]:=Sum.arr^[j]*K;
   Move(AY.Ptr(0).arr^,REZX.arr^,REZX.Count*SOfR);
   Move(AY.Ptr(1).arr^,REZY.arr^,REZY.Count*SOfR);
   time:=at+Tau^
  end else begin  //Расчёт по отдельным сериям
 L2:Data.arr^[Position]:=AU.Ptr(0).arr^[0];
   if Position=Size^-1 then begin
    //Вычисление коэффициентов линии регрессии и удаление линейного тренда
    N:=0;CN2:=0;CNY:=0;M:=0;SmSqr:=0;
    for i:=0 to Size^-1 do begin
      N:=i+1;
      CN2:=CN2+N*N;
      CNY:=CNY+N*Data.arr^[i];
      M:=M+Data.arr^[i];
    end;
    if (DelTrend^=1) or (OutMode^=3) then begin
     if DelTrend^=1 then b:=(CNY-0.5*M*N)/(CN2-N*sqr(0.5*N)) else b:=0;
     a:=M/N-b*0.5*N;
     end else begin b:=0;a:=0 end;
     //Расчёт границ интервалов
    if AutoRange^=1 then begin
    DMin:=MaxDouble;
    DMax:=-MaxDouble;
    for i:=0 to Size^-1 do begin
      Data.arr^[i]:=Data.arr^[i]-(b*(i+1)+a);
      h:=Data.arr^[i];
      if h>DMax then DMax:=h;
      if h<DMin then DMin:=h;
    end;
    end else begin
     DMin:=Min^;
     DMax:=Max^
    end;
     //Ширина интервала гистограммы
    h:=(DMax-DMin)/Col^;
    if h=0 then h:=0.000001;
     //Расчёт гистограммы
    case OutMode^ of
     0: K:=1/Size^;
     2: K:=1/(Size^*h);
     3: begin
         for i:=0 to Size^-1 do SmSqr:=SmSqr+sqr(Data.arr^[i]);
         K:=sqrt(SmSqr/(Size^-1));
         a:=1/K;
         K:=K/(Size^*h);
         for i:=0 to Col^-1 do AY.Ptr(0).arr^[i]:=(DMin+h*(i+0.5))*a;
         goto L4;
        end;
    else
     K:=1
    end;
    for i:=0 to Col^-1 do AY.Ptr(0).arr^[i]:=DMin+h*(i+0.5);
L4: AY.Ptr(1).FillArray(0);
    for i:=0 to Size^-1 do
     if (Data.arr^[i]>=Dmin)and(Data.arr^[i]<=DMax) then begin
      if Data.arr^[i]>=DMax then j:=Col^-1 else j:=trunc((Data.arr^[i]-DMin)/h);
      if Auto then Sum.arr^[j]:=Sum.arr^[j]+1;
      AY.Ptr(1).arr^[j]:=AY.Ptr(1).arr^[j]+K
    end;
    Move(AY.Ptr(0).arr^,REZX.arr^,REZX.Count*SOfR);
    Move(AY.Ptr(1).arr^,REZY.arr^,REZY.Count*SOfR);
   end;
     //Переустановка стека
    if Position>=Size^-1 then
     if Auto then Auto:=false else Position:=0
     else inc(Position);
    time:=time+Tau^
  end
 end   //END Case
 end   //END If
end;   //END Procedure

function TStatHistRst;
 var CNT:integer;
begin
 Result:=0;
 with PHistProp(Prop.arr)^,Vars^ do
 case Action of
 f_ReadRez   : try
                BlockRead(F,REZX.arr^,REZX.Count*SOfR);
                BlockRead(F,REZY.arr^,REZY.Count*SOfR);
               except
                Result:=er_ReadFile
               end;
 f_WriteRez  : try
                BlockWrite(F,REZX.arr^,REZX.Count*SOfR);
                BlockWrite(F,REZY.arr^,REZY.Count*SOfR);
               except
                Result:=er_WriteFile
               end;
 f_ReadRst   : try
                BlockRead(F,NumStates,SizeOf(NumStates));
                BlockRead(F,Sum.arr^,Sum.Count*SOfR);
                BlockRead(F,REZX.arr^,REZX.Count*SOfR);
                BlockRead(F,REZY.arr^,REZY.Count*SOfR);
                BlockRead(F,CNT,SOfI);
                Data.FillArray(0);
                if CNT>Size^ then Data.ChangeCount(CNT);
                BlockRead(F,Data.arr^,CNT*SOfR)
               except
                Result:=er_ReadFile
               end;
 f_WriteRst  : try
                BlockWrite(F,NumStates,SizeOf(NumStates));
                BlockWrite(F,Sum.arr^,Sum.Count*SOfR);
                BlockWrite(F,REZX.arr^,REZX.Count*SOfR);
                BlockWrite(F,REZY.arr^,REZY.Count*SOfR);
                BlockWrite(F,Size^,SOfI);
                BlockWrite(F,Data.arr^,Size^*SOfR);
               except
                Result:=er_WriteFile
               end;
 end
end;


end.