Programmer 7500/MAX_PROGRAMMERS.iso/INFO/TURBOPAS/NONLIN.ZIP/SIMPLEX.PAS

Pascal/Delphi Source File | 1985-12-28 | 21.3 KB | 709 lines

program simp;   {curve fitting with the simplex algorithm}
const

{********************************************************************}
{THESE PARAMETERS MUST BE CHANGED FOR EACH DIFFERENT MODEL TO BE FIT}
     m=2;                    {number of parameters to fit}
     n=3;                    {n=m+1}
     nvpp=3;                 {total # of vars per data point}
{*********************************************************************}

        variance = false;       { bootstrap ?}
        debug = false;
        mnp=900;                {maximum # of points -- depends on memory available}
        no_simplex=100;         {number of times bootstrap performed}
        alfa=1.0;               {reflection coefficient, > 0}
        beta= 0.5;              {contraction coeff, 0 to 1 }
        gamma = 2.0;            {expansion coeff, > 1}
        lw = 5;                 {width of line in data fields +1}
        page=12;
        root2=1.414214;

type
        vector=array [1..n] of real;
        datarow=array [1..nvpp] of real;
        index=0..mnp;
        input = array[1..mnp] of datarow;

var
        residual,plots: boolean;        {flags for output files and plots}
        ok,done:boolean;                {flag for simplex}
        i,j:index;
        h,l:array[1..n] of index;       {number high/low paramts}
        np,                             {number of records read}
        maxiter,                        {max # iterations}
        niter:integer;                  {# of iterations}
        next,                           {new vertex to be tested}
        center,                         {center of hyperplane described
                        by all vertexes of the simplex excluding the worst}
        mean, error,
        maxerr,                         {maximum error accepted}
        p,q,                            {to compute the first simplex}
        step :vector;                   {input starting steps}
        simp:array[1..n] of vector;     {the simplex}
        data1,data: input;              {the data}
        fname: string[14];              {input file name}
        nerr,nfit:string[14];           {file names for fitted and residual files}
        ferr,ffit,din,
        dout,fin,fout: text;            {file variables}
        printer: boolean;               {printer flag}
        response:char;
        title:string[80];               {for title}
        position:integer;               {new position from 'sample' used in bootstrap}
        stat:boolean;                   {for bootstrap}
        bootstrap:integer;              {"    "}
        sample_no:string[19];           {"    "}

{**********************************************************************}


{*********************************************************************}
{$i nonfunc.pas}  {has function definition}
{*********************************************************************}



procedure sum_of_residuals (var x:vector;matrix:input);
                {computes the sum of the squares of the residuals}
                {x(1..m) passes parameters. Result returned in x(n)}
var
        i:index;


begin
        x[n]:=0.0;
        if debug then writeln ('np:= ',np);
        for i:=1 to np do
                begin
                if debug then writeln (matrix[i,1],matrix[i,2]);
                x[n]:=x[n] + sqr (f(x,matrix[i])-matrix[i,nvpp]);
                if debug then writeln ('Sum of residuals--x[',i,']:= ',x[n]);
                end
        end;    {sum of residuals}

{*******************************************************************}

procedure enter;
                {enters data from disk file fname}

var
   f1:text;            {parameters for errors, starting points, incr}

begin                                                   {enter}
        assign (f1, 'a:parms.dat');
        reset (f1);
        writeln(dout);
        writeln(dout, 'Simplex optimization program              ');
        writeln(dout,' Taken from Byte May 1984--written at FSU');
        writeln(dout,'Adapted by Tom Dean, Dept. For. Res.,Utah State University');
        writeln(dout);
        writeln(dout,' accessing file ',fname);
        read (f1,maxiter);
        writeln(dout, 'maximum # of iterations is: ',maxiter:5);
        write(dout,' start coord.: ');
        for i:= 1 to m do
                begin
                read(f1,simp[1,i]);
                if (i mod lw)=0 then writeln (dout);
                write(dout,simp[1,i])
                end;
        writeln(dout);
        write(dout, ' start steps: ');
        for i:=1 to m do
                begin
                read(f1,step[i]);
                if (i mod lw)=0 then writeln(dout);
                write(dout,step[i])
                end;
        writeln(dout);
        write(dout, ' max. errors: ');
        for i:= 1 to n do
                begin
                read (f1,maxerr[i]);
                if (i mod lw)=0 then writeln(dout);
                write(dout,maxerr[i])
                end;
        close (f1);
        writeln(dout);
        writeln(dout,' data: ');
        writeln(dout,' x':14,'y':14);
        np:=0;
        while not EOF(din) do
                begin
                np:=succ(np);
                for j:= 1 to (nvpp-1) do read (din, data[np,j]);
                readln (din, data[np,nvpp]);
                end;
        for j:=1 to 4 do
            begin
            for i:=1 to nvpp do write (dout, data[j,i]);
            writeln (dout);
            end;
        writeln(dout);
        writeln(dout, 'Total # of records read= ',np);
        writeln(dout);
end;{enter}


{**********************************************************************}

procedure report;
type
    plotinp    = array [1..mnp] of real;
    variable   = string[3];
var
        sumy, ssy, sstot, r2,
        y,dy,
        sigma:real;
        k:integer;
        ind,res,pred,obs         : plotinp;
        xchar                    : variable;


procedure plot (varx,vary:variable; m:integer;  x,y,yhat:plotinp);
{if m neg then yhat plotted and its symbol takes precedence
* = > 10 pts at a point}

const
     y_sym  = 'o';
     yhat_sym = 'p';
     over_sym = '*';
     xlen = 51;
     ylen = 25;

type
    points     = array [1..xlen, 1..ylen] of char;
    counters =integer;

var
   np,i,j,k,l                : counters;
   done, plot_pred           : boolean;
   graph                     : points;
   xscale,yscale,yhscale,
   xnext,ynext               : real;
   xmin,xmax,ymin,ymax,
   yhmin,yhmax, temp         : real;
   symbol,tempsym            : char;
   mx,my,myh                 : real;
   x1,y1                     : integer;


begin
{find out if yhat to be plotted}
if m < 0 then
begin
     plot_pred:=false;
     m:=-m;
end
else plot_pred:=true;

fillchar(graph, sizeof(graph), $20);      {initialize graph to blanks}

{find max and mins}
xmin:=x[1];
xmax:=x[2];
ymin:=y[1];
ymax:=y[2];
if Plot_pred then
begin
     yhmin:=yhat[1];
     yhmax:=yhat[2];
end;
for i:=2 to m do
begin
     if xmin > x[i] then xmin:=x[i];
     if xmax < x[i] then xmax:=x[i];
     if ymin > y[i] then ymin:=y[i];
     if ymax < y[i] then ymax:=y[i];
     if plot_pred then
     begin
          if yhmin > yhat[i] then yhmin:=yhat[i];
          if yhmax < yhat[i] then yhmax:=yhat[i];
     end;
end;

{calc scales increments}
xscale:=(xmax-xmin)/(xlen-1);
yscale:=(ymax-ymin)/(ylen-1);
if plot_pred then
begin
     yhscale:=(yhmax-yhmin)/(ylen-1);
     if yhscale > yscale then
     begin
          yscale:=yhscale;     {if yscale larger then use these max mins for stan}
          ymax:=yhmax;
          ymin:=yhmin;
     end;
end;


{get ready to fill graph -- decide on symbols}
if plot_pred then symbol:=y_sym  else symbol:='1';  { char 31h}
{Standardize data to fit into graph space}
mx:=(49)/(xmax-xmin);
my:=(24)/(ymax-ymin);
for i:=1 to m do
begin
     x1:=round(1+mx*(x[i]-xmin));
     y1:=round(1+my*(y[i]-ymin));
     if plot_pred then graph[x1,y1]:=symbol
     else
     begin
          tempsym:=graph[x1,y1];
          case tempsym of
               ' ' : graph[x1,y1]:= symbol;
               '9' : graph[x1,y1]:='*'
               else graph[x1,y1]:=succ(tempsym);
          end;
     end;
end;

{overlay predicted symbol on graph space}
if plot_pred then
begin
     symbol:=yhat_sym;
     for i:= 1 to m do
     begin
          x1:=round(1+mx*(x[i]-xmin));
          y1:=round(1+my*(yhat[i]-ymin));
          graph[x1,y1]:=symbol;
     end;
end;


{output the plot}
write (dout,'  Variable ',vary);
if plot_pred then
   writeln(dout,'   p = predicted   o = observed')
   else writeln (dout);
write(dout,'           ');
writeln(dout,'-+---------+---------+---------+---------+---------+');
ynext:=ymax;
for j:= ylen downto 1 do
begin
     if (j mod 5) = 0 then
     begin
          write(dout,ynext:10,' +');
          ynext:=ynext-yscale
     end
     else
     begin
          write (dout,'           |');
          ynext:=ynext-yscale;
     end;
     for i:= 1 to xlen-1 do write (dout,graph[i,j]);
     writeln (dout,graph[xlen,j]);
end;
write(dout,ynext:10,' ');
writeln(dout,'-+---------+---------+---------+---------+---------+');
xnext:=xmin;
write (dout,'       ',xnext:8);
for i:= 1 to 5 do
begin
     xnext:=xnext+xscale;
     write (dout,'   ',xnext:8);
end;
writeln(dout);writeln(dout);
writeln (dout,'                             Variable ',varx);
writeln(dout);
end; { PLOT }
{************************************************************}



begin
        writeln(dout,' program exited after',niter:5,' iterations');
        writeln(dout);
        writeln(dout,' 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(dout);
                        write(dout,simp[j,i]:20)
                        end;
                writeln(dout);
                end;
        writeln(dout);
        writeln(dout,' the mean is');
        for i:=1 to n do
                begin
                if (i mod lw)=0 then writeln(dout);
                write(dout,mean[i]);
                end;
        writeln(dout);
        writeln(dout);
        writeln(dout,' the estimated fractional error is');
        for i:= 1 to n do
                begin
                if(i mod lw)=0 then writeln(dout);
                write(dout, error[i])
                end;
        writeln(dout);
        writeln(dout);
        writeln (dout);
        sumy:=0;
        ssy:=0;
        for i:= 1 to np do
            begin
            ssy:=ssy + sqr(data[i,nvpp]);
            sumy:=sumy + data[i,nvpp];
            end;
        sstot:= ssy-sqr(sumy)/np;
        r2:=(sstot-mean[n])/sstot;
        writeln (dout, 'The approximate R-square is: ',r2);
        writeln(dout);
        sigma:=0.0;
        for i:=1 to np do
            begin
             y:= f(mean, data[i]);
             dy:=data[i,nvpp]-y;
             sigma:=sigma+ sqr(dy);
             end;
        sigma:=sqrt(sigma/np);
        writeln(dout,' the standard deviation is',sigma);
        writeln(dout);
        sigma:=sigma/sqrt(np-m);
        write(dout,' the estimated error of the');
        writeln(dout,' function is: ',sigma,^L,^G);

    if plots then
    begin

        { Plot residuals by each independent variable}

        for j:= 1 to nvpp-1 do
        begin
             for i:= 1 to np do
             begin
                  ind[i]:= data[i,j];
                  res[i]:=data[i,nvpp]-f(mean,data[i]);
             end;
             str(j,xchar);
             plot (xchar,'res',-np,ind,res,res);
             writeln(dout);writeln(dout);
        end;
        {this is a special plot for a specific function}
        {*********************************************}
        for i:= 1 to np do
        begin
             ind[i]:=data[i,1]*data[i,2];
             res[i]:=data[i,nvpp]-f(mean,data[i]);
        end;
        plot('1*2','res',-np,ind,res,res);
        {*********************************************}

        for j:=1 to nvpp-1 do
        begin
             for i:=1 to np do
             begin
                  ind[i]:=data[i,j];
                  obs[i]:=data[i,nvpp];
                  pred[i]:=f(mean,data[i]);
             end;
             str(j,xchar);
             plot(xchar,'3',np,ind,obs,pred);
         end;
         {************************************************
         this is a special graph for a particular function I use}
         for i:=1 to np do
         begin
              ind[i]:=data[i,1]*data[i,2];
              obs[i]:=data[i,nvpp];
              pred[i]:=f(mean,data[i]);
         end;
         plot ('1*2','3',np,ind,obs,pred);
         {************************************************}
    end;



end;    {report}

{*********************************************************************}

procedure first;
        begin
        writeln(dout,' starting simplex');
        for j:=1 to n do
                begin
                write(dout,' simp[',j:1,']');
                for i:=1 to n do
                        begin
                        if(i mod lw)=0 then writeln(dout);
                        write(dout,simp[j,i])
                        end;
                writeln(dout)
                end;
        writeln(dout)
end;    {first}

{*******************************************************************}

procedure new_vertex;           {next in place of the worst vertex}
begin
        if not stat then write(dout,' ---',niter:4);
        for     i:=1 to n do
                begin
                simp[h[n],i]:=next[i];
                if not stat then write(next[i])
                end;
        if not stat then writeln(dout);
end;    {new vertex}

{********************************************************************}

procedure order;                {gives high.low in each parameter}
                                {in simp. CAUTION:NOT INITIALIZED}
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;    {order}

{*******************************************************************}
procedure simplex ( matrix:input);
begin
sum_of_residuals(simp[1],matrix);      {first vertex}

        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],matrix)
                end;

        for i:=1 to n do
                begin
                l[i]:=1;h[i]:=1
                end;

        order;

        if not stat then
        begin
        first;
        if printer then
           begin
           assign(dout,'lst:');
           rewrite(dout);
           first;
           assign(dout,'con:');
           rewrite(dout);
           end;
        end;

        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,matrix);

        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,matrix);
                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,matrix);
                        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],matrix)
                                        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;
end;
{****************************************************************************}

begin           {simplex main program}
repeat
        write ('Name of data file: ');
        readln(fname);          {input file }
{$i-}   assign(din,fname);      {fname is on disk}
        reset (din);
{$i+}   ok:=(ioresult=0);
        if not ok then
        begin
             writeln ('Cannot find file ',fname);
             writeln;
        end;
until ok;

        write(con, 'Plot residuals and fitted data? (Y/N): ');
        read (kbd,response);
        writeln(upcase(response));
        if upcase(response) = 'Y'  then plots:=true else plots:=false;

        write (con, 'Save residuals and fitted data in files? (Y/N): ');
        read(kbd,response);
        writeln (upcase(response));
        if upcase(response)= 'Y' then residual:= true else residual:=false;

        write ('Should output go to printer? (Y/N): ');
        read (kbd,response);
        writeln (upcase(response));
        if upcase(response) = 'Y' then printer:=true else printer:=false;

        if printer then
           begin
           lowvideo;
           writeln('Make sure printer is on.');
           normvideo;
           write('Enter descriptive title for run');
           writeln('    Blank line terminates input');
           repeat
                 readln (title);
                 writeln(lst,title);
           until title = '';
           end;
        assign (dout,'con:');
        rewrite(dout);
        enter;                  {get the data}
        if printer then
           begin
           assign(dout,'lst:');
           rewrite(dout);
           reset (din);
           enter;
           assign(dout,'con:');
           rewrite(dout);
           end;
stat:=false;
simplex(data);

        report;
        if printer then
           begin
           assign(dout,'lst:');
           rewrite(dout);
           report;
           end;

if residual then
begin
     i:=0;
     repeat i:= i +1 until  (fname[i] = '.');
     nerr:= copy(fname,1,(i-1)) + '.err';
     nfit:= copy(fname,1,(i-1)) + '.fit';
     writeln ;
     writeln ('Residual and fit files written to ',nerr,' and ',nfit);
     assign (ferr,nerr);
     assign (ffit,nfit);
     rewrite(ferr);
     rewrite (ffit);
     for i:= 1 to np do
     begin
          for j:= 1 to nvpp-1 do
          begin
               write (ferr, data[i,j]);
               write (ffit, data[i,j]);
          end;
          writeln (ferr, data[i,nvpp]-f(mean, data[i]));
          writeln (ffit, f(mean, data[i]));
     end;
     close (ferr);
     close (ffit);
end;
{ Beginning of bootstrap routine ++ Beware this can take forever without a
math coprocessor}
if variance then begin
printer:=false;
stat:=true;
bootstrap:=1;

assign(fin,'sample');
reset(fin);
assign(fout,'b:bootstra');
rewrite(fout);

repeat
      readln(fin);
      readln(fin,sample_no);
      for i:=1 to np do
      begin
           read(fin,position);
           data1[i]:=data[position];
      end;
      simplex(data1);
      write(fout, bootstrap,' ');
      write(con, bootstrap,' ');
      for i:=1 to n do write (fout,mean[i]);
      writeln(fout);
      for i:=1 to n do write (con, mean[i]);
      writeln(con);
      readln(fin);
      bootstrap:=bootstrap+1;
until bootstrap = no_simplex+1;
close (fin);
close(fout);
writeln (con, ^g,^g ,'          +++ DONE+++');
end;{if variance}

end.