sources-request@panda.UUCP (02/23/86)
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."