Path: seismo!harvard!talcott!panda!sources-request
From: sources-request@panda.UUCP
Newsgroups: mod.sources
Subject: Simplex Curve Fitting Algorithm in C
Message-ID: <1430@panda.UUCP>
Date: 23 Feb 86 01:57:44 GMT
Sender: jpn@panda.UUCP
Lines: 454
Keywords: simplex C
Approved: jpn@panda.UUCP

Mod.sources:  Volume 4, Issue 2
Submitted by: panda!genrad!decvax!ihnp4!chinet!blm

This program originally appeared in the May 1984 issue of Byte Magazine.  It
was originally written in Pascal by Marco Caceci and William Caceris at Florida
State University.  I have translated it to 'C'.

This program is based upon the Simplex curve fitting algorithm.  For a detailed
descripstion of this program and it's workings see the above mentioned article.

I acknowledge the work of Marco Caceci and William Caceris for writing the
original Pascal program from which this is derived.  The original authors
explicitly stated ''no copy-right''.

I've had some problems with accuracy.  I've checked and rechecked my source
code and I can't find anything that would account for it.  I could very well be
overlooking something, though.  I suppose it could have to do with differences
in the precision of the floating-point libraries of the authors Pascal compiler
(Pascal/Z Version 4.0 CP/M) and my C compiler (Xenix 3.0).  I welcome E-mail
on this subject.  Nevertheless, the differences in values returned between
the original and mine are comparatively small.

I hope this helps someone in Netland.
-----
The preceding announcement was conceived in my own mind and is not
necessarily the opinion of whom ever I happen to work for at the time.

Name  : Brad L. McKinley --- (503) 889-4321
USMail: M D R Professional Software, Inc., 915 SW 3rd Avenue, Ontario, OR 97914
Usenet: ihnp4!chinet!blm OR ihnp4!chinet!mdr!blm
CIS   : 70015,1205
Source: BDY171
"God created Arrakis to test the faithful." -- Maud'dib

--- cut here ------ cut here ------ cut here ------ cut here ---
: This is a shar archive.  Extract with sh, not csh.
: The rest of this file will extract:
:      Makefile
:      data              sample data file
:      enter.c
:      f.c               a typical function
:      first.c
:      main.c
:      new_vertex.c
:      order.c
:      report.c
:      sum_residual.c
:
echo extracting -- Makefile
sed 's/^X//' > Makefile << E_O_F
XBINARY = simplex
X
XSOURCES = main.c f.c sum_residual.c enter.c first.c new_vertex.c order.c report.c
X
XOBJECTS = main.o f.o sum_residual.o enter.o first.o new_vertex.o order.o report.o
X
XLIBRARIES = -lm -lc
X
XCFLAGS = -O
XLDFLAGS = -n -s -x
XLINTFLAGS =
X
X$(BINARY): $(OBJECTS)
X	@echo "	loading $(BINARY)"
X	@ld -o $@ $(LDFLAGS) /lib/crt0.o $(OBJECTS) $(LIBRARIES)
X
Xlint:
X	lint $(LINTFLAGS) $(SOURCES)
X
X$(OBJECTS): simplex.h
E_O_F

echo extracting -- data
sed 's/^X//' > data << E_O_F
X100
X0.2 3
X0.1 1
X1e-6 1e-6 1e-6
X1.68 0.172
X3.33 0.250
X5.00 0.286
X6.67 0.303
X10.0 0.334
X20.0 0.384
E_O_F

echo extracting -- simplex.h
sed 's/^X//' > simplex.h << E_O_F
X#include <math.h>
X#include <stdio.h>
X
X#define M      2
X#define NVPP   2
X#define N      M+1
X#define MNP    200
X#define ALPHA  1.0
X#define BETA   0.5
X#define GAMMA  2.0
X#define LW     5
X#define ROOT2  1.414214
X
Xtypedef double vector[N];
Xtypedef double datarow[NVPP];
X
Xextern int       h[N], l[N];
Xextern int       np, maxiter, niter;
Xextern vector    next, mean, error, maxerr, step, simp[N];
Xextern datarow   data[MNP];
Xextern FILE      *fpdata;
X
Xextern double    f();
E_O_F

echo extracting -- enter.c
sed 's/^X//' > enter.c << E_O_F
X#include "simplex.h"
X
Xenter(fname)
Xchar *fname;
X{
X     register int   i, j;
X
X     printf("SIMPLEX Optimization --- 'C' Version\n\n");
X     printf("Accessing file: %s\n\n", fname);
X
X     fscanf(fpdata, "%d", &maxiter);
X     printf("maximum number of iterations = %d\n\n", maxiter);
X
X     printf("Start Coordinates: ");
X     for (i=0 ; i<M ; i++) {
X          fscanf(fpdata, "%F", &simp[0][i]);
X          if ((i+1) % LW == 0)
X               printf("\n");
X          printf(" %e", simp[0][i]);
X     }
X     printf("\n\n");
X
X     printf("Start Steps: ");
X     for (i=0 ; i<M ; i++) {
X          fscanf(fpdata, "%F", &step[i]);
X          if ((i+1) % LW == 0)
X               printf("\n");
X          printf(" %e", step[i]);
X     }
X     printf("\n\n");
X
X     printf("Maximum Error: ");
X     for (i=0 ; i<N ; i++) {
X          fscanf(fpdata, "%F", &maxerr[i]);
X          if ((i+1) % LW == 0)
X               printf("\n");
X          printf(" %e", maxerr[i]);
X     }
X     printf("\n\n");
X
X     printf("DATA\n");
X     printf("      X             Y\n");
X     np = 0;
X     while (!feof(fpdata)) {
X          for (j=0 ; j<NVPP ; j++) {
X               if (fscanf(fpdata, "%F", &data[np][j]) == EOF)
X                    break;
X               printf("  %e", data[np][j]);
X          }
X          np++;
X          printf("\n");
X     }
X     np--;
X
X}
E_O_F

echo extracting -- f.c
sed 's/^X//' > f.c << E_O_F
X#include "simplex.h"
X
Xdouble f(x, d)
Xvector    x;
Xdatarow   d;
X{
X     return (x[0]*d[0]/(x[1]+d[0]));
X}
E_O_F

echo extracting -- first.c
sed 's/^X//' > first.c << E_O_F
X#include "simplex.h"
X
Xfirst()
X{
X     register int   i, j;
X
X     printf("Starting Simplex\n");
X     for (j=0 ; j<N ; j++) {
X          printf(" simp[%d]", j+1);
X          for (i=0 ; i<N ; i++) {
X               if ((i+1) % LW == 0)
X                    printf("\n");
X               printf("  %e", simp[j][i]);
X          }
X          printf("\n");
X     }
X     printf("\n");
X}
E_O_F

echo extracting -- main.c
sed 's/^X//' > main.c << E_O_F
X#include "simplex.h"
X
X#define until(x) while (!(x))
X
Xint       h[N], l[N];
Xint       np, maxiter, niter;
Xvector    next, mean, error, maxerr, step, simp[N];
Xdatarow   data[MNP];
XFILE      *fpdata;
X
Xmain(argc, argv)
Xint   argc;
Xchar  *argv[];
X{
X   register int   i, j, done;
X   vector         center, p, q;
X
X   if (argc != 2) {
X      fprintf(stderr, "usage: simplex file_name\n", argv[1]);
X      exit(1);
X   }
X
X   if ((fpdata = fopen(argv[1], "r")) == NULL) {
X      fprintf(stderr, "simplex: can't open %s\n", argv[1]);
X      exit(1);
X   }
X
X   enter(argv[1]);
X
X   /* First Vertex */
X   sum_residual(simp[0]);
X
X   /* Compute offset of Vertices */
X   for (i=0 ; i<M ; i++) {
X      p[i] = step[i] * (sqrt((double) N) + M - 1) / (M * ROOT2);
X      q[i] = step[i] * (sqrt((double) N) - 1) / (M * ROOT2);
X   }
X
X   /* All Vertices of the Starting Simplex */
X   for (i=1 ; i<N ; i++) {
X      for (j=0 ; j<M ; j++)
X         simp[i][j] = simp[0][j] + q[j];
X      simp[i][i-1] = simp[0][i-1] + p[i-1];
X      sum_residual(simp[i]);
X   }
X
X   /* Preset */
X   for (i=0 ; i<N ; i++) {
X      l[i] = 1;
X      h[i] = 1;
X   }
X   order();
X
X   first();
X
X   niter = 0;
X
X   /* Iterate */
X   do {
X      /* Wish it were True */
X      done = 1;
X      niter++;
X
X      /* Compute Centroid...Excluding the Worst */
X      for (i=0 ; i<N ; i++)
X         center[i] = 0.0;
X      for (i=0 ; i<N ; i++)
X         if (i != h[N-1])
X            for (j=0 ; j<M ; j++)
X               center[j] += simp[i][j];
X
X      /* First Attempt to Reflect */
X      for (i=0 ; i<N ; i++) {
X         center[i] /= M;
X         next[i] = (1.0+ALPHA) * center[i] - ALPHA * simp[h[N-1]][i];
X      }
X      sum_residuals(next);
X
X      if (next[N-1] <= simp[l[N-1]][N-1]) {
X         new_vertex();
X         for (i=0 ; i<M ; i++)
X            next[i] = GAMMA * simp[h[N-1]][i] + (1.0-GAMMA) * center[i];
X         sum_residual(next);
X         if (next[N-1] <= simp[l[N-1]][N-1])
X            new_vertex();
X      }
X      else {
X         if (next[N-1] <= simp[h[N-1]][N-1])
X            new_vertex();
X         else {
X            for (i=0 ; i<M ; i++)
X               next[i] = BETA * simp[h[N-1]][i] + (1.0-BETA) * center[i];
X            sum_residual(next);
X            if (next[N-1] <= simp[h[N-1]][N-1])
X               new_vertex();
X            else {
X               for (i=0 ; i<N ; i++) {
X                  for (j=0 ; j<M ; j++)
X                     simp[i][j] = BETA * (simp[i][j] + simp[l[N-1]][j]);
X                  sum_residual(simp[i]);
X               }
X            }
X         }
X      }
X
X      order();
X
X      /* Check For Convergence */
X      for (j=0 ; j<N ; j++) {
X         error[j] = (simp[h[j]][j] - simp[l[j]][j]) / simp[h[j]][j];
X         if (done)
X            if (error[j] > maxerr[j])
X               done = 0;
X      }
X
X   } until(done || (niter == maxiter));
X
X   /* Average Each Parameter */
X   for (i=0 ; i<N ; i++) {
X      mean[i] = 0.0;
X      for (j=0 ; j<N ; j++)
X         mean[i] += simp[j][i];
X      mean[i] /= N;
X   }
X
X   report();
X}
E_O_F

echo extracting -- new_vertex.c
sed 's/^X//' > new_vertex.c << E_O_F
X#include "simplex.h"
X
Xnew_vertex()
X{
X     register int   i;
X
X     printf(" --- %4d", niter);
X     for (i=0 ; i<N ; i++) {
X          simp[h[N-1]][i] = next[i];
X          printf("  %e", next[i]);
X     }
X     printf("\n");
X}
E_O_F

echo extracting -- order.c
sed 's/^X//' > order.c << E_O_F
X#include "simplex.h"
X
Xorder()
X{
X     register int   i, j;
X
X     for (j=0 ; j<N ; j++)
X          for (i=0 ; i<N ; i++) {
X               if (simp[i][j] < simp[l[j]][j])
X                    l[j] = i;
X               if (simp[i][j] > simp[h[j]][j])
X                    h[j] = i;
X          }
X}
E_O_F

echo extracting -- report.c
sed 's/^X//' > report.c << E_O_F
X#include "simplex.h"
X
Xreport()
X{
X     register int   i, j;
X     double         y, dy, sigma;
X
X     printf("\nProgram exited after %d iterations.\n\n", niter);
X
X     printf("The final simplex is:\n");
X     for (j=0 ; j<N ; j++) {
X          for (i=0 ; i<N ; i++) {
X               if ((i+1) % LW == 0)
X                    printf("\n");
X               printf("  %e", simp[j][i]);
X          }
X          printf("\n");
X     }
X     printf("\n\n");
X
X     printf("The mean is:");
X     for (i=0 ; i<N ; i++) {
X          if ((i+1) % LW == 0)
X               printf("\n");
X          printf("  %e", mean[i]);
X     }
X     printf("\n\n");
X
X     printf("The estimated fractional error is:");
X     for (i=0 ; i<N ; i++) {
X          if ((i+1) % LW == 0)
X               printf("\n");
X          printf("  %e", error[i]);
X     }
X     printf("\n\n");
X
X     printf("   #         X              Y             Y''             DY\n");
X     sigma = 0.0;
X     for (i=0 ; i<np ; i++) {
X          y = f(mean, data[i]);
X          dy = data[i][1] - y;
X          sigma += (dy*dy);
X          printf("%4d  ", i);
X          printf("%13e  %13e  ", data[i][0], data[i][1]);
X          printf("%13e  %13e\n", y, dy);
X     }
X     printf("\n");
X     sigma = sqrt(sigma);
X     printf("The standard deviation is %e\n\n", sigma);
X     sigma /= sqrt((double) (np-M));
X     printf("The estimated error of the function is %e\n\n", sigma);
X}
E_O_F

echo extracting -- sum_residual.c 
sed 's/^X//' > sum_residual.c << E_O_F
X#include "simplex.h"
X
X#define sqr(x) ((x) * (x))
X
Xsum_residual(x)
Xvector    x;
X{
X     register int   i;
X
X     x[N-1] = 0.0;
X
X     for (i=0 ; i<np ; i++)
X          x[N-1] += sqr(f(x, data[i]) - data[i][1]);
X}
E_O_F

exit

---
Name  : Brad L. McKinley --- (503) 889-4321
Usenet: ihnp4!chinet!blm             US Mail: M D R Professional Software, Inc.
CIS   : 70015,1205                            915 SW 3rd Avenue
Source: BDY171                                Ontario, Oregon 97914
"God created Arrakis to test the faithful."