Power-Programmierung/CD1.mdf/pascal/library/dos/simplex/simplex2.pas

Pascal/Delphi Source File | 1986-04-15 | 23.2 KB | 924 lines

Overlay Procedure Simplex;
                  { Curve fitting with the Simplex Algorithm}
                  { m is the number of parameters to fit, n       }
                  { is m + 1, nvpp = number of vars per data point}
Const
  n_max              = 11;   { m + 1                             }
  nvpp_max           = 10;   { maximum no. of variables per data point}
  maxiter            = 100;  { maximum no. of iterations         }
  alfa               = 1.0;  { Reflection coefficient, > 0       }
  beta               = 0.5;  { Contraction coefficient, in(0..1) }
  gamma              = 2.0;  { expansion coefficient, > 1        }
  LW                 = 5;
  MaxNumLength       = 10;
  MaxExpression      = 100;
  MaxToken           = 100;
  MaxPriority        = 6;
  MaxParam           = 50;
  NameLength         = 7;
  FirstUnary         = 4;
  LastUnary          = 15;
  FirstBinary        = 16;
  LastBinary         = 22;
  LastOperand        = 24;
  MaxString          = 200;

Type
  vector             = array[1..n_max] of real;
  index              = 0..255;
  PriorRange         = 1..MaxPriority;
  Name               = String[NameLength];
  Token              = 0..MaxToken;
  TokenKind          = (Operand,UnaryOp,BinaryOp,EndExpression,LeftParen,
                        Variable,RightParen);
  DeftToken          = Record
                         NM      : Name;
                       Case K: TokenKind of
                         Variable,
                         Operand    :  (Val : Real);
                         UnaryOp,
                         BinaryOp   :  (Pri : PriorRange);
                         EndExpression,
                         LeftParen,
                         RightParen : ()
                       End;
  Expression         = array[1..MaxExpression] of Token;
  ExprIndex          = 0..MaxExpression;
  Param              = 0..MaxParam;
  Message            = String[255];


Var
  vars          : array[1..nvpp_max] of index;
  done          : boolean;
  Nvars         : index;
  H,L           : array[1..n_max] of index;
  niter,m,n,
  I,J           : integer;
  next, center,
  mean, error,maxerr,
  p,q,step      : vector;
  simp          : array[1..n_max] of vector;
  Root2,Z       : real;
  PostFix       : Expression;
  Lexicon       : Array[Token] of DeftToken;
  Nparameter    : Param;
  Parameter     : Array[Param] of Token;

  InString      : String[MaxString];

{$V-}
Procedure Errors(Line: Message);
Begin
  TextBackground(4) ; textcolor(0);
  Writeln(' ERROR!! ');
  TextBackground(0); TextColor(15);
  Writeln(Line);
  Halt;
End;
{$V+}

Procedure DefineTokens;
Begin
  Lexicon[1].NM  :=  '';
  Lexicon[1].K   := EndExpression;
  Lexicon[2].NM  := '(';
  Lexicon[2].K   := LeftParen;
  Lexicon[3].NM  := ')';
  Lexicon[3].K   := RightParen;
  Lexicon[4].NM  := '\';
  Lexicon[4].K   := UnaryOp;
  Lexicon[4].Pri := 6;
  Lexicon[5].NM  := 'ABS';
  Lexicon[5].K   := UnaryOp;
  Lexicon[5].Pri := 6;
  Lexicon[6].NM  := 'SQR';
  Lexicon[6].K   := UnaryOp;
  Lexicon[6].Pri := 6;
  Lexicon[7].NM  := 'SQRT';
  Lexicon[7].K   := UnaryOp;
  Lexicon[7].Pri := 6;
  Lexicon[8].NM  := 'EXP';
  Lexicon[8].K   := UnaryOp;
  Lexicon[8].Pri := 6;
  Lexicon[9].NM  := 'LN';
  Lexicon[9].K   := UnaryOp;
  Lexicon[9].Pri := 6;
  Lexicon[10].NM  := 'LOG';
  Lexicon[10].K   := UnaryOp;
  Lexicon[10].Pri := 6;
  Lexicon[11].NM  := 'SIN';
  Lexicon[11].K   := UnaryOp;
  Lexicon[11].Pri := 6;
  Lexicon[12].NM  := 'COS';
  Lexicon[12].K   := UnaryOp;
  Lexicon[12].Pri := 6;
  Lexicon[13].NM  := 'ARCTAN';
  Lexicon[13].K   := UnaryOp;
  Lexicon[13].Pri := 6;
  Lexicon[14].NM  := 'ROUND';
  Lexicon[14].K   := UnaryOp;
  Lexicon[14].Pri := 6;
  Lexicon[15].NM  := 'TRUNC';
  Lexicon[15].K   := UnaryOp;
  Lexicon[15].Pri := 6;
  Lexicon[16].NM  := '+';
  Lexicon[16].K   := BinaryOp;
  Lexicon[16].Pri := 4;
  Lexicon[17].NM  := '-';
  Lexicon[17].K   := BinaryOp;
  Lexicon[17].Pri := 4;
  Lexicon[18].NM  := '*';
  Lexicon[18].K   := BinaryOp;
  Lexicon[18].Pri := 5;
  Lexicon[19].NM  := '/';
  Lexicon[19].K   := BinaryOp;
  Lexicon[19].Pri := 5;
  Lexicon[20].NM  := 'DIV';
  Lexicon[20].K   := BinaryOp;
  Lexicon[20].Pri := 5;
  Lexicon[21].NM  := 'MOD';
  Lexicon[21].K   := BinaryOp;
  Lexicon[21].Pri := 5;
  Lexicon[22].NM  := '^';
  Lexicon[22].K   := BinaryOp;
  Lexicon[22].Pri := 6;
  Lexicon[23].NM  := 'PI';
  Lexicon[23].K   := Operand;
  Lexicon[23].Val := 3.14159;
  Lexicon[24].NM  := 'E';
  Lexicon[24].K   := Operand;
  Lexicon[24].Val := 2.71828;
End;

Function Kind(T: Token) : TokenKind;
Begin
  Kind := Lexicon[T].K
End;


Procedure StartExpression;
Const
  HashSize        = 101;

Type
  Address         = 0..HashSize;
  IndexName       = 0..NameLength;
  IndexString     = 0..MaxString;

Var
  InFix           : Expression;
  H               : Array[Address] of Token;
  TokenCount      : Token;

Function Hash(X: Name) : Address;
Var
  A     : Integer;
  Ch    : Char;
  Found : Boolean;
Begin
  If Length(X) <= 0 then
    Errors('Hash attempted with empty name')
  Else Begin
    Ch := X[1];
    A  := Abs( Ord(Ch)) mod HashSize;
    Repeat
      If H[A] = 0 then
        Found := True
      Else if Lexicon[H[A]].NM = X then
        Found := True
      Else Begin
        If Length(X) > 1 then
          Begin
            Ch := X[2];
            A  := A + Abs(Ord(Ch));
          End
        Else
          A := A + 29;
        If A > HashSize then A := A mod HashSize;
        Found := False
      End
    Until Found;
    Hash := A
  End
End;

Procedure MakeHashTable;
Var
  A     : Address;
  T     : Token;

Begin
  For A := 0 to HashSize do
    H[A] := 0;
  For T := 2 to LastOperand do
     H[Hash(Lexicon[T].NM) ] := T;
End;



Procedure ReadExpression;
Var
  ExprLength     : ExprIndex;
  Position       : IndexString;
  ParenCount     : Integer;
  Digit,Alphabet,
  Lower          : Set of Char;

Procedure PutToken(X: Token);
Begin
  ExprLength := Succ(ExprLength);
  InFix[ExprLength] := X;
End;


Function Leading: Boolean;
Var K : TokenKind;
Begin
  If ExprLength = 0 then
    Leading := True
  Else begin
    K := Kind(InFix[ExprLength]);
    Leading := (K = LeftParen) or (K = UnaryOp) or (K = BinaryOp)
  End;
End;

Procedure FindWord;
Var
  Word           : Name;
  T              : Token;
  I              : IndexName;
  NewPosition    : IndexString;

Begin
  NewPosition := Position + 1;
  While InString[NewPosition] in (Alphabet + Digit) do
    NewPosition := Newposition + 1;
  If NewPosition - Position <= NameLength then
    Word := Copy(InString,Position,NewPosition - Position)
  Else Begin
    Word := Copy(Instring,Position,NameLength);
    Writeln('Warning: the name ',Word,' has been truncated');
  End;
  For I := 1 to Length(Word) do
    If Word[I] in Lower then Word[I] := UpCase(Word[I]);
  T := H[Hash(word)];
  If T <> 0 then
    If Leading then
      If Kind(T) = BinaryOp then
        Errors('Binary Operator in illegal Postion')
      Else
        PutToken(T)
    Else
      If Kind(T) <> BinaryOp then
        Errors('Binary Operator Expected')
      Else
        PutToken(T)
  Else
    If TokenCount >= MaxToken then
      Errors('Too many distinct variables and constants')
    Else if not Leading then
      Errors('Operand immediately follows ) or another Operand')
    Else begin
      TokenCount := TokenCount + 1;
      H[Hash(Word)] := TokenCount;
      Lexicon[TokenCount].NM := Word;
      Lexicon[TokenCount].K  := Operand;
      If Nparameter >= MaxParam then
        Errors('Too many parameters')
      Else begin
        Nparameter := Nparameter + 1;
        Parameter[Nparameter] := TokenCount;
        PutToken(TokenCount);
      End;
    End;
  Position := NewPosition;
End;

Procedure FindNumber;
Var
  NumberName,X          : String[MaxNumLength];
  DecPoint,NewPosition  : IndexString;
  R                     : Real;
  Code                  : Integer;

Begin
  If Not Leading then
    Errors('Constant in Illegal Position')
  Else if TokenCount >= MaxToken then
    Errors('Too many Constants and Variables')
  Else Begin
    NewPosition := Position;
    While Instring[NewPosition] in (digit + ['.','e','E']) do
      NewPosition := NewPosition + 1;
    X := Copy(Instring,Position,NewPosition - Position);
    If Length(X) <= Namelength then
      NumberName := X
    Else
      NumberName := Copy(X,1,NameLength);
    Val(X,R,Code);
    TokenCount   := TokenCount + 1;
    Lexicon[TokenCount].NM  := NumberName;
    Lexicon[TokenCount].K   := Operand;
    Lexicon[TokenCount].Val := R;
    PutToken(TokenCount);
    Position := NewPosition;
  End;
End;

Procedure FindSymbol;
Var
  X  : Name;
  T  : Token;

Begin
  X := Copy(InString,Position,1);
  T := H[Hash(X)];
  If T = 0 then
    Begin
      Writeln(X,' ',T);
      Errors('Unrecognized symbol in Expression')
    End
  Else if Leading then
    If Kind(T) = RightParen then
      Errors('Illegal place for closing parenthesis')
    Else if Kind(T) = BinaryOp then
      If X = '+' then begin end
      Else if X = '-' then
      Begin
        X := '\';
        T := H[Hash(X)];
        PutToken(T);
      End
      Else
        Errors('Binary Operator in Illegal Position')
    Else
      PutToken(T)
  Else
    If (Kind(T) = RightParen) or (Kind(T) = BinaryOp) then
      PutToken(T)
    Else
      Errors('Binary operator or ) expected');
  If Kind(T) = LeftParen then
    ParenCount := ParenCount + 1
  Else if Kind(T) = RightParen then
    Begin
      ParenCount := ParenCount - 1;
      If Parencount < 0 then Errors('More right than left parentheses')
    End;
    Position := Position + 1
End;

Begin
  TokenCount := LastOperand;
  Writeln(' Type in the Expression to evaluate on the following line: ');
  Write(' Y = ');Readln(Instring);
  Instring   := Instring + ' ';
  ExprLength := 0;  NParameter := 0; ParenCount := 0;
  Position   := 1;
  Lower      := ['a'..'z'];
  Alphabet   := Lower + ['A'..'Z'];
  Digit      := ['0'..'9'];
  While Position <= Length(Instring) do
    If Instring[Position] = ' ' then
      Position := Position + 1
    Else if instring[Position] in alphabet then
      FindWord
    Else if instring[Position] in Digit + ['.'] then
      FindNumber
    Else
      FindSymbol;
  If ParenCount <> 0 then
    Errors('Numbers of Left and Right Parentheses are not Equal');
  If Leading then
    Errors('Input Expression is Incomplete.');
  PutToken(1)
End;

Procedure Translate;
Const
  MaxStack  = 100;

Var
  Stack     : Array[1..MaxStack] of Token;
  Nstack    : 0..MaxStack;
  T,X       : Token;
  EndRight  : Boolean;
  InFixCount,PostFixCount    : 0..100;

Procedure Push(X : Token);
Begin
  Nstack := Nstack + 1;
  Stack[Nstack] := X;
End;

Procedure Pop(Var X : Token);
Begin
  X := Stack[Nstack];
  Nstack := Nstack - 1;
End;

Procedure GetToken(Var X: Token);
Begin
  X := InFix[InFixCount];
  InFixCount := Succ(InFixCount);
End;

Procedure PutToken(X: Token);
Begin
  PostFix[PostFixCount] := X;
  PostFixCount := Succ(PostFixCount);
End;

Function Priority(T : Token): PriorRange;
Begin
  Priority := Lexicon[T].Pri;
End;

Begin
  Nstack := 0;
  InFixCount := 1;  PostFixCount := 1;
  Repeat
    GetToken(T);
    Case Kind(T) of
      Operand    : PutToken(T);
      LeftParen  : Push(T);
      RightParen : Begin
                     Pop(T);
                     While Kind(T) <> LeftParen do
                     Begin
                       PutToken(T);
                       Pop(T)
                     End
                   End;
      UnaryOp,
      BinaryOp   : Begin
                     Repeat
                       If Nstack = 0 then
                         EndRight := True
                       Else If Kind(Stack[Nstack]) = LeftParen then
                         EndRight := True
                       Else If Priority(Stack[Nstack]) < Priority(T) then
                         EndRight := True
                       Else Begin
                         EndRight := False;
                         Pop(X);
                         PutToken(X);
                       End
                     Until EndRight;
                     Push(T);
                   End;
      EndExpression :
                  While Nstack > 0 do
                    Begin
                      Pop(X);
                      PutToken(X);
                    End;
    End;
  Until Kind(T) = EndExpression;
  PutToken(T);
End;

Begin
  MakeHashTable;
  ReadExpression;
  Translate
End;

Procedure ReadParameters;
Var
  I,J   : Param;

Begin
  If Nparameter > 0 then
    Begin
      ClrScr;
      Writeln('The following parameters are present in your equation.');
      For I := 1 to Nparameter do
        With Lexicon[Parameter[I]] do
        Begin
          Writeln(I,':  ',NM:Namelength);
        End;
      Write(^J' How many of these parameters are variables? ');
      Readln(Nvars);
      m := nparameter - nvars;
      n := m + 1;
      Writeln(^J' Enter the parameter number for each variable : '^J);
      For I := 1 to Nvars do
        Begin
          Write(' Variable ',I,': ');
          Readln(J);
          Lexicon[Parameter[J]].K := Variable;
        End;
      ClrScr;
      J := 0;
      For I := 1 to Nparameter do
        With Lexicon[Parameter[I]] do
          Begin
            If  K = Variable then
              Begin
                J := J + 1;
                Writeln(NM,'  corresponds to which column of');
                Write(' data in your data set? ');
                Readln(Vars[J]);
              End;
          End;
      ClrScr;
      Writeln('Enter initial guesses for the following coefficients ');
      J := 0;
      For I := 1 to Nparameter do
        With Lexicon[Parameter[I]] do
        Begin
          If K = Operand then
            Begin
              J := J + 1;
              Write(NM:Namelength,'?  ');
              Readln(Simp[1,J]);
              Step[J] := Simp[1,J]/5.0;
            End;
        End;
      Write('Which column of your data is the dependant (Y) variable? ');
      Readln(IY);
    End;
End;

Procedure EvaluatePostFix(Var Result: Real);
Const
  MaxStack   = 100;
Var
  Stack      : Array[1..MaxStack] of Real;
  Nstack     : 0..MaxStack;
  T          : Token;
  A,B        : Real;
  PostFixCount : 1..MaxExpression;

Function GetValue(T: Token): Real;
Begin
  If (Kind(T) <> Operand) and (Kind(T) <> Variable) then
    Errors('Attempt to get Value for non-Operand')
  Else
    GetValue := Lexicon[T].Val
End;

Procedure Push(V: Real);
Begin
  If Nstack >= MaxStack then Errors('Stack OverFlow')
  Else Begin
    Nstack := Nstack + 1;
    Stack[Nstack] := V
  End
End;

Procedure Pop(Var V: Real);
Begin
  If Nstack <= 0 then Errors('Stack UnderFlow')
  Else Begin
    V := Stack[Nstack];
    Nstack := Nstack - 1
  End
End;


Function Exponent(X,Y: Real): Real;
Const
  Epsilon    = 0.000001;

Var
  I: Integer;
  P: Real;

Begin
  If Abs(Y - Round(Y)) < Epsilon then
  Begin
    P := 1.0;
    If Y >= 0.0 then
      For I := 1 to Round(Y) do P := P*X
    Else if X = 0.0 then Errors('Negative Power of 0.0')
    Else
      For I := -1 downto Round(Y) do
        P := P/X;
    Exponent := P;
  End
  Else if X > 0.0 then
    Exponent := Exp(Y*Ln(X))
  Else if Abs(X) < Epsilon then
    Exponent := 0.0
  Else
    Errors('Attempt to take negative number to non-integer power')
End;

Function DoBinary(T: Token; X,Y: Real): Real;
Begin
  If (T < FirstBinary) or (T > LastBinary) then
    Errors('Binary operator code out of range')
  Else Case T of
    16: DoBinary := X + Y;
    17: DoBinary := X - Y;
    18: DoBinary := X * Y;
    19: If Y = 0 then Errors('Division by 0') else DoBinary := X/Y;
    20: If Round(Y) = 0 then Errors('Integer Division by 0') else
        DoBinary := Round(X) div Round(Y);
    21: If Round(Y) = 0 then Errors('Attempt to use 0 Modulus') else
          DoBinary := Round(X) mod Round(Y);
    22: DoBinary := Exponent(X,Y);
  End;
End;

Function DoUnary(T: Token; X: Real): Real;
Begin
  If (T < FirstUnary) or (T > LastUnary) then
    Errors('Unary Operator code out of Range')
  Else Case T of
    4:  DoUnary := -1.0*X;
    5:  DoUnary := Abs(X);
    6:  DoUnary := Sqr(X);
    7:  DoUnary := Sqrt(X);
    8:  DoUnary := Exp(X);
    9:  DoUnary := Ln(X);
    10: DoUnary := Ln(X)/2.302585093;
    11: DoUnary := Sin(X);
    12: DoUnary := Cos(X);
    13: DoUnary := ArcTan(X);
    14: DoUnary := Round(X);
    15: DoUnary := Trunc(X);
  End;
End;

Procedure GetToken(Var T: Token);
Begin
  T := PostFix[PostFixCount];
  PostFixCount := Succ(PostFixCount);
End;


Begin
  PostFixCount := 1;
  Nstack := 0;
  Repeat
    GetToken(T);
    Case Kind(T) of
      Operand,Variable:
                Begin
                  Push(GetValue(T));
                End;
      UnaryOp :
                Begin
                  Pop(A);
                  Push(DoUnary(T,A))
                End;
     BinaryOp :
                Begin
                  Pop(B);
                  Pop(A);
                  Push(DoBinary(T,A,B))
                End;
     EndExpression :
                If Nstack = 1 then
                   Pop(Result)
                 Else
                   Errors('Stack Problem');
    End;
  Until Kind(T) = EndExpression;
End;

Function f (II: Integer ; X: vector) : real;
Var
  I,J,L      : Index;
begin
  J := 0;
  L := 0;
    For I := 1 to Nparameter do
      With Lexicon[Parameter[I]] do
        Begin
          If K = Variable then
            Begin
              J := J + 1;
              Val := GetX(II,Vars[J]);
            End
          Else
            Begin
              L := L + 1;
              Val := X[L];
            End;
        End;
    EvaluatePostfix(Z);
    F := Z;
  end;


Procedure Sum_of_residuals(var X : vector);
Var
  I   : Integer;
Begin
  x[n] := 0.0;
  For I := 1 to Nobs do
    Begin
      x[n] := x[n] + sqr(f(I,X) - GetX(I,IY));
    End;
End;

Procedure enter;
  Begin
    Write(IO,' SIMPLEX optimization curve fitting ');
    Writeln(IO,^J^J'Evaluating the Expression: ');
    Writeln(IO,^J' ',InString);
    Writeln(IO,^J' max number of iterations = ',maxiter:5);
    Write(IO,' start coord.: ');
    For I := 1 to M do
      Begin
        If (I mod LW) = 0 then writeln(IO);
        Write(IO,Simp[1,I]);
      End;
    Writeln(IO);
    Write(IO,' Start Steps: ');
    For I := 1 to M do
      Begin
        If (I mod LW) = 0 then writeln(IO);
        Write(IO,Step[I]);
      End;
    Writeln(IO);
    Write(IO,' Maximum Errors: ');
    For I := 1 to N do
      Begin
        maxerr[I] := 1.0E-04;
        If (I mod LW) = 0 then writeln(IO);
        Write(IO,maxerr[I])
      End;
    Writeln(IO);
    Nobs := Nobs;
  End;
Procedure Report;
Var
  Y, DY, Sigma     : real;
Begin
  GotoXY(1,24);Writeln;
  Writeln(IO,' Program exited after ',niter:5,' iteration ');
  Writeln(IO,' The final SIMPLEX is');
  For J := 1 to N do
    Begin
      For I := 1 to N do
        Begin
          If (I mod LW) = 0 then Writeln(IO);
          Write(IO,Simp[j,i]:10,'  ')
        End;
      Writeln(IO);
    End;
  Writeln(IO,' The mean is');
  For I := 1 to N do
    Begin
      If (I mod LW) = 0 then Writeln;
      Write(IO,Mean[I]);
    End;
  Writeln(IO);
  Writeln(IO,' The estimated fractional error is');
  For I := 1 to N do
    Begin
      If (I mod LW) = 0 then Writeln;
      Write(IO,error[I]);
    End;
  Writeln(IO);
  Sigma := 0.0;
  For I := 1 to Nobs do
    Begin
      Y := f(I,mean);
      dY := GetX(I,IY) - Y;
      PutX(Y,I,FC);
      PutX(dY,I,RC);
      Sigma := Sigma + Sqr(dY);
    End;
  Sigma := Sqrt(Sigma / Nobs);
  Writeln(IO,' The standard deviation is ',Sigma);
  Sigma := Sigma / Sqrt(Nobs - M);
  Writeln(IO,' The estimated error of the function is', Sigma);
End;

Procedure First;
  Begin
    Writeln(' Starting Simplex');
    For J := 1 to N do
      Begin
        Write(' Simp[',J:1,']');
        For I := 1 to N do
          Begin
            If (I mod LW) = 0 then Writeln;
            Write(Simp[J,I]);
          End;
        Writeln;
      End;
      Writeln;
  GotoXY(50,1);
  TextBackground(15); TextColor(0);
  Write('   Iteration No.              ');
  TextBackground(0); TextColor(0);
End;
Procedure New_vertex;
  Begin
    GotoXY(65,1);Write(niter:4);
    For I := 1 to N do
      Begin
        Simp[H[N],I] := Next[I];
      End;
    Writeln;
  End;
Procedure Order;
Var
  I, J            : Index;

  Begin
    For J := 1 to N do
      Begin
        For I := 1 to N do
          Begin
            If Simp[I,J] < Simp[L[J],J] then L[J] := I;
            If Simp[I,J] > Simp[H[J],J] then H[J] := I;
          End;
      End;
  End;
{****************************************************************************}
{Start of Main                                                               }
Begin
  Root2 := Sqrt(2.00);
  ClrScr;
  DefineTokens;
  StartExpression;
  ReadParameters;
  Enter;
  Sum_of_residuals(Simp[1]);
  For I := 1 to M do
    Begin
      P[I] := Step[I]*(Sqrt(N) + M - 1)/(M * Root2);
      Q[I] := Step[I]*(Sqrt(N) -1)/(M * Root2);
    End;
  For I := 2 to N do
    Begin
      For J := 1 to M do Simp[I,J] := Simp[1,J] + Q[J];
      Simp[I,I-1] := Simp[1,I-1] + P[I-1];
      Sum_of_residuals(Simp[I]);
    End;
  For I := 1 to N do
    Begin;
      L[I] := 1; H[I] := 1;
    End;
  Order;
  First;
  Niter := 0;
  Repeat;
  Done := True;
  Niter := Succ(Niter);
  For I := 1 to N do Center[I] := 0.0;
  For I := 1 to N do
     If I <> H[N] then For J := 1 to M do Center[J] := Center[J] + Simp[I,J];
  For I := 1 to N do
    Begin
      Center[I] := Center[I]/ M;
      Next[I] := (1.0 + Alfa)* Center[I] - Alfa*Simp[H[N],I]
    End;
  Sum_of_residuals(next);
  If Next[N] <= Simp[L[N],N] then
    Begin
      New_vertex;
      For I := 1 to M do
         Next[I] := Gamma*Simp[H[N],I] + (1.0 - Gamma)*Center[I];
      Sum_of_Residuals(next);
      If Next[N] <= Simp[L[N],N] then New_Vertex;
    End
  Else
    Begin
      If Next[N] <= Simp[H[N],N] then
        New_vertex
     Else
       Begin
         For I := 1 to M do
           Next[I] := Beta*Simp[H[N],I] + (1.0 - Beta)*Center[I];
         Sum_of_residuals(next);
         If Next[N] <= Simp[H[N],N] then
           New_vertex
         Else
           Begin
             For I := 1 to N do
               Begin
                 For J := 1 to M do
                   Simp[I,J] := (Simp[I,J] + Simp[L[N],J])*Beta;
                 Sum_of_residuals(Simp[I]);
               End
           End
       End
    End;
    Order;
    For J := 1 to N do
      Begin
        Error[J] := (Simp[H[J],J] - Simp[L[J],J])/Simp[H[J],J];
        If done then
          If error[J] > Maxerr[J] then  Done := false;
      End;
    Until (done or (Niter = Maxiter));
    For I := 1 to N do
      Begin
        Mean[I] := 0.0;
        For J := 1 to N do Mean[I] := Mean[I] + Simp[J,I];
        Mean[I] := Mean[I]/N;
      End;
    Report;
  Pause;
End;