home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / unix / volume4 / simplex < prev    next >
Text File  |  1986-11-30  |  11KB  |  462 lines
  1. Subject: Simplex Curve Fitting Algorithm in C
  2. Keywords: simplex C
  3. Newsgroups: mod.sources
  4. Approved: jpn@panda.UUCP
  5.  
  6. Mod.sources:  Volume 4, Issue 2
  7. Submitted by: panda!genrad!decvax!ihnp4!chinet!blm
  8.  
  9. This program originally appeared in the May 1984 issue of Byte Magazine.  It
  10. was originally written in Pascal by Marco Caceci and William Caceris at Florida
  11. State University.  I have translated it to 'C'.
  12.  
  13. This program is based upon the Simplex curve fitting algorithm.  For a detailed
  14. descripstion of this program and it's workings see the above mentioned article.
  15.  
  16. I acknowledge the work of Marco Caceci and William Caceris for writing the
  17. original Pascal program from which this is derived.  The original authors
  18. explicitly stated ''no copy-right''.
  19.  
  20. I've had some problems with accuracy.  I've checked and rechecked my source
  21. code and I can't find anything that would account for it.  I could very well be
  22. overlooking something, though.  I suppose it could have to do with differences
  23. in the precision of the floating-point libraries of the authors Pascal compiler
  24. (Pascal/Z Version 4.0 CP/M) and my C compiler (Xenix 3.0).  I welcome E-mail
  25. on this subject.  Nevertheless, the differences in values returned between
  26. the original and mine are comparatively small.
  27.  
  28. I hope this helps someone in Netland.
  29. -----
  30. The preceding announcement was conceived in my own mind and is not
  31. necessarily the opinion of whom ever I happen to work for at the time.
  32.  
  33. Name  : Brad L. McKinley --- (503) 889-4321
  34. USMail: M D R Professional Software, Inc., 915 SW 3rd Avenue, Ontario, OR 97914
  35. Usenet: ihnp4!chinet!blm OR ihnp4!chinet!mdr!blm
  36. CIS   : 70015,1205
  37. Source: BDY171
  38. "God created Arrakis to test the faithful." -- Maud'dib
  39.  
  40. --- cut here ------ cut here ------ cut here ------ cut here ---
  41. : This is a shar archive.  Extract with sh, not csh.
  42. : The rest of this file will extract:
  43. :      Makefile
  44. :      data              sample data file
  45. :      enter.c
  46. :      f.c               a typical function
  47. :      first.c
  48. :      main.c
  49. :      new_vertex.c
  50. :      order.c
  51. :      report.c
  52. :      sum_residual.c
  53. :
  54. echo extracting -- Makefile
  55. sed 's/^X//' > Makefile << E_O_F
  56. XBINARY = simplex
  57. X
  58. XSOURCES = main.c f.c sum_residual.c enter.c first.c new_vertex.c order.c report.c
  59. X
  60. XOBJECTS = main.o f.o sum_residual.o enter.o first.o new_vertex.o order.o report.o
  61. X
  62. XLIBRARIES = -lm -lc
  63. X
  64. XCFLAGS = -O
  65. XLDFLAGS = -n -s -x
  66. XLINTFLAGS =
  67. X
  68. X$(BINARY): $(OBJECTS)
  69. X    @echo "    loading $(BINARY)"
  70. X    @ld -o $@ $(LDFLAGS) /lib/crt0.o $(OBJECTS) $(LIBRARIES)
  71. X
  72. Xlint:
  73. X    lint $(LINTFLAGS) $(SOURCES)
  74. X
  75. X$(OBJECTS): simplex.h
  76. E_O_F
  77.  
  78. echo extracting -- data
  79. sed 's/^X//' > data << E_O_F
  80. X100
  81. X0.2 3
  82. X0.1 1
  83. X1e-6 1e-6 1e-6
  84. X1.68 0.172
  85. X3.33 0.250
  86. X5.00 0.286
  87. X6.67 0.303
  88. X10.0 0.334
  89. X20.0 0.384
  90. E_O_F
  91.  
  92. echo extracting -- simplex.h
  93. sed 's/^X//' > simplex.h << E_O_F
  94. X#include <math.h>
  95. X#include <stdio.h>
  96. X
  97. X#define M      2
  98. X#define NVPP   2
  99. X#define N      M+1
  100. X#define MNP    200
  101. X#define ALPHA  1.0
  102. X#define BETA   0.5
  103. X#define GAMMA  2.0
  104. X#define LW     5
  105. X#define ROOT2  1.414214
  106. X
  107. Xtypedef double vector[N];
  108. Xtypedef double datarow[NVPP];
  109. X
  110. Xextern int       h[N], l[N];
  111. Xextern int       np, maxiter, niter;
  112. Xextern vector    next, mean, error, maxerr, step, simp[N];
  113. Xextern datarow   data[MNP];
  114. Xextern FILE      *fpdata;
  115. X
  116. Xextern double    f();
  117. E_O_F
  118.  
  119. echo extracting -- enter.c
  120. sed 's/^X//' > enter.c << E_O_F
  121. X#include "simplex.h"
  122. X
  123. Xenter(fname)
  124. Xchar *fname;
  125. X{
  126. X     register int   i, j;
  127. X
  128. X     printf("SIMPLEX Optimization --- 'C' Version\n\n");
  129. X     printf("Accessing file: %s\n\n", fname);
  130. X
  131. X     fscanf(fpdata, "%d", &maxiter);
  132. X     printf("maximum number of iterations = %d\n\n", maxiter);
  133. X
  134. X     printf("Start Coordinates: ");
  135. X     for (i=0 ; i<M ; i++) {
  136. X          fscanf(fpdata, "%F", &simp[0][i]);
  137. X          if ((i+1) % LW == 0)
  138. X               printf("\n");
  139. X          printf(" %e", simp[0][i]);
  140. X     }
  141. X     printf("\n\n");
  142. X
  143. X     printf("Start Steps: ");
  144. X     for (i=0 ; i<M ; i++) {
  145. X          fscanf(fpdata, "%F", &step[i]);
  146. X          if ((i+1) % LW == 0)
  147. X               printf("\n");
  148. X          printf(" %e", step[i]);
  149. X     }
  150. X     printf("\n\n");
  151. X
  152. X     printf("Maximum Error: ");
  153. X     for (i=0 ; i<N ; i++) {
  154. X          fscanf(fpdata, "%F", &maxerr[i]);
  155. X          if ((i+1) % LW == 0)
  156. X               printf("\n");
  157. X          printf(" %e", maxerr[i]);
  158. X     }
  159. X     printf("\n\n");
  160. X
  161. X     printf("DATA\n");
  162. X     printf("      X             Y\n");
  163. X     np = 0;
  164. X     while (!feof(fpdata)) {
  165. X          for (j=0 ; j<NVPP ; j++) {
  166. X               if (fscanf(fpdata, "%F", &data[np][j]) == EOF)
  167. X                    break;
  168. X               printf("  %e", data[np][j]);
  169. X          }
  170. X          np++;
  171. X          printf("\n");
  172. X     }
  173. X     np--;
  174. X
  175. X}
  176. E_O_F
  177.  
  178. echo extracting -- f.c
  179. sed 's/^X//' > f.c << E_O_F
  180. X#include "simplex.h"
  181. X
  182. Xdouble f(x, d)
  183. Xvector    x;
  184. Xdatarow   d;
  185. X{
  186. X     return (x[0]*d[0]/(x[1]+d[0]));
  187. X}
  188. E_O_F
  189.  
  190. echo extracting -- first.c
  191. sed 's/^X//' > first.c << E_O_F
  192. X#include "simplex.h"
  193. X
  194. Xfirst()
  195. X{
  196. X     register int   i, j;
  197. X
  198. X     printf("Starting Simplex\n");
  199. X     for (j=0 ; j<N ; j++) {
  200. X          printf(" simp[%d]", j+1);
  201. X          for (i=0 ; i<N ; i++) {
  202. X               if ((i+1) % LW == 0)
  203. X                    printf("\n");
  204. X               printf("  %e", simp[j][i]);
  205. X          }
  206. X          printf("\n");
  207. X     }
  208. X     printf("\n");
  209. X}
  210. E_O_F
  211.  
  212. echo extracting -- main.c
  213. sed 's/^X//' > main.c << E_O_F
  214. X#include "simplex.h"
  215. X
  216. X#define until(x) while (!(x))
  217. X
  218. Xint       h[N], l[N];
  219. Xint       np, maxiter, niter;
  220. Xvector    next, mean, error, maxerr, step, simp[N];
  221. Xdatarow   data[MNP];
  222. XFILE      *fpdata;
  223. X
  224. Xmain(argc, argv)
  225. Xint   argc;
  226. Xchar  *argv[];
  227. X{
  228. X   register int   i, j, done;
  229. X   vector         center, p, q;
  230. X
  231. X   if (argc != 2) {
  232. X      fprintf(stderr, "usage: simplex file_name\n", argv[1]);
  233. X      exit(1);
  234. X   }
  235. X
  236. X   if ((fpdata = fopen(argv[1], "r")) == NULL) {
  237. X      fprintf(stderr, "simplex: can't open %s\n", argv[1]);
  238. X      exit(1);
  239. X   }
  240. X
  241. X   enter(argv[1]);
  242. X
  243. X   /* First Vertex */
  244. X   sum_residual(simp[0]);
  245. X
  246. X   /* Compute offset of Vertices */
  247. X   for (i=0 ; i<M ; i++) {
  248. X      p[i] = step[i] * (sqrt((double) N) + M - 1) / (M * ROOT2);
  249. X      q[i] = step[i] * (sqrt((double) N) - 1) / (M * ROOT2);
  250. X   }
  251. X
  252. X   /* All Vertices of the Starting Simplex */
  253. X   for (i=1 ; i<N ; i++) {
  254. X      for (j=0 ; j<M ; j++)
  255. X         simp[i][j] = simp[0][j] + q[j];
  256. X      simp[i][i-1] = simp[0][i-1] + p[i-1];
  257. X      sum_residual(simp[i]);
  258. X   }
  259. X
  260. X   /* Preset */
  261. X   for (i=0 ; i<N ; i++) {
  262. X      l[i] = 1;
  263. X      h[i] = 1;
  264. X   }
  265. X   order();
  266. X
  267. X   first();
  268. X
  269. X   niter = 0;
  270. X
  271. X   /* Iterate */
  272. X   do {
  273. X      /* Wish it were True */
  274. X      done = 1;
  275. X      niter++;
  276. X
  277. X      /* Compute Centroid...Excluding the Worst */
  278. X      for (i=0 ; i<N ; i++)
  279. X         center[i] = 0.0;
  280. X      for (i=0 ; i<N ; i++)
  281. X         if (i != h[N-1])
  282. X            for (j=0 ; j<M ; j++)
  283. X               center[j] += simp[i][j];
  284. X
  285. X      /* First Attempt to Reflect */
  286. X      for (i=0 ; i<N ; i++) {
  287. X         center[i] /= M;
  288. X         next[i] = (1.0+ALPHA) * center[i] - ALPHA * simp[h[N-1]][i];
  289. X      }
  290. X      sum_residuals(next);
  291. X
  292. X      if (next[N-1] <= simp[l[N-1]][N-1]) {
  293. X         new_vertex();
  294. X         for (i=0 ; i<M ; i++)
  295. X            next[i] = GAMMA * simp[h[N-1]][i] + (1.0-GAMMA) * center[i];
  296. X         sum_residual(next);
  297. X         if (next[N-1] <= simp[l[N-1]][N-1])
  298. X            new_vertex();
  299. X      }
  300. X      else {
  301. X         if (next[N-1] <= simp[h[N-1]][N-1])
  302. X            new_vertex();
  303. X         else {
  304. X            for (i=0 ; i<M ; i++)
  305. X               next[i] = BETA * simp[h[N-1]][i] + (1.0-BETA) * center[i];
  306. X            sum_residual(next);
  307. X            if (next[N-1] <= simp[h[N-1]][N-1])
  308. X               new_vertex();
  309. X            else {
  310. X               for (i=0 ; i<N ; i++) {
  311. X                  for (j=0 ; j<M ; j++)
  312. X                     simp[i][j] = BETA * (simp[i][j] + simp[l[N-1]][j]);
  313. X                  sum_residual(simp[i]);
  314. X               }
  315. X            }
  316. X         }
  317. X      }
  318. X
  319. X      order();
  320. X
  321. X      /* Check For Convergence */
  322. X      for (j=0 ; j<N ; j++) {
  323. X         error[j] = (simp[h[j]][j] - simp[l[j]][j]) / simp[h[j]][j];
  324. X         if (done)
  325. X            if (error[j] > maxerr[j])
  326. X               done = 0;
  327. X      }
  328. X
  329. X   } until(done || (niter == maxiter));
  330. X
  331. X   /* Average Each Parameter */
  332. X   for (i=0 ; i<N ; i++) {
  333. X      mean[i] = 0.0;
  334. X      for (j=0 ; j<N ; j++)
  335. X         mean[i] += simp[j][i];
  336. X      mean[i] /= N;
  337. X   }
  338. X
  339. X   report();
  340. X}
  341. E_O_F
  342.  
  343. echo extracting -- new_vertex.c
  344. sed 's/^X//' > new_vertex.c << E_O_F
  345. X#include "simplex.h"
  346. X
  347. Xnew_vertex()
  348. X{
  349. X     register int   i;
  350. X
  351. X     printf(" --- %4d", niter);
  352. X     for (i=0 ; i<N ; i++) {
  353. X          simp[h[N-1]][i] = next[i];
  354. X          printf("  %e", next[i]);
  355. X     }
  356. X     printf("\n");
  357. X}
  358. E_O_F
  359.  
  360. echo extracting -- order.c
  361. sed 's/^X//' > order.c << E_O_F
  362. X#include "simplex.h"
  363. X
  364. Xorder()
  365. X{
  366. X     register int   i, j;
  367. X
  368. X     for (j=0 ; j<N ; j++)
  369. X          for (i=0 ; i<N ; i++) {
  370. X               if (simp[i][j] < simp[l[j]][j])
  371. X                    l[j] = i;
  372. X               if (simp[i][j] > simp[h[j]][j])
  373. X                    h[j] = i;
  374. X          }
  375. X}
  376. E_O_F
  377.  
  378. echo extracting -- report.c
  379. sed 's/^X//' > report.c << E_O_F
  380. X#include "simplex.h"
  381. X
  382. Xreport()
  383. X{
  384. X     register int   i, j;
  385. X     double         y, dy, sigma;
  386. X
  387. X     printf("\nProgram exited after %d iterations.\n\n", niter);
  388. X
  389. X     printf("The final simplex is:\n");
  390. X     for (j=0 ; j<N ; j++) {
  391. X          for (i=0 ; i<N ; i++) {
  392. X               if ((i+1) % LW == 0)
  393. X                    printf("\n");
  394. X               printf("  %e", simp[j][i]);
  395. X          }
  396. X          printf("\n");
  397. X     }
  398. X     printf("\n\n");
  399. X
  400. X     printf("The mean is:");
  401. X     for (i=0 ; i<N ; i++) {
  402. X          if ((i+1) % LW == 0)
  403. X               printf("\n");
  404. X          printf("  %e", mean[i]);
  405. X     }
  406. X     printf("\n\n");
  407. X
  408. X     printf("The estimated fractional error is:");
  409. X     for (i=0 ; i<N ; i++) {
  410. X          if ((i+1) % LW == 0)
  411. X               printf("\n");
  412. X          printf("  %e", error[i]);
  413. X     }
  414. X     printf("\n\n");
  415. X
  416. X     printf("   #         X              Y             Y''             DY\n");
  417. X     sigma = 0.0;
  418. X     for (i=0 ; i<np ; i++) {
  419. X          y = f(mean, data[i]);
  420. X          dy = data[i][1] - y;
  421. X          sigma += (dy*dy);
  422. X          printf("%4d  ", i);
  423. X          printf("%13e  %13e  ", data[i][0], data[i][1]);
  424. X          printf("%13e  %13e\n", y, dy);
  425. X     }
  426. X     printf("\n");
  427. X     sigma = sqrt(sigma);
  428. X     printf("The standard deviation is %e\n\n", sigma);
  429. X     sigma /= sqrt((double) (np-M));
  430. X     printf("The estimated error of the function is %e\n\n", sigma);
  431. X}
  432. E_O_F
  433.  
  434. echo extracting -- sum_residual.c 
  435. sed 's/^X//' > sum_residual.c << E_O_F
  436. X#include "simplex.h"
  437. X
  438. X#define sqr(x) ((x) * (x))
  439. X
  440. Xsum_residual(x)
  441. Xvector    x;
  442. X{
  443. X     register int   i;
  444. X
  445. X     x[N-1] = 0.0;
  446. X
  447. X     for (i=0 ; i<np ; i++)
  448. X          x[N-1] += sqr(f(x, data[i]) - data[i][1]);
  449. X}
  450. E_O_F
  451.  
  452. exit
  453.  
  454. ---
  455. Name  : Brad L. McKinley --- (503) 889-4321
  456. Usenet: ihnp4!chinet!blm             US Mail: M D R Professional Software, Inc.
  457. CIS   : 70015,1205                            915 SW 3rd Avenue
  458. Source: BDY171                                Ontario, Oregon 97914
  459. "God created Arrakis to test the faithful."
  460.  
  461.  
  462.