23. search Example

% cat exp.ss

// creates noisy exponential data values y = 2*exp(0.5*x),
// then finds minimum max-absolute-error approximation z = f1*exp(f2*x)

srand 12345;	// just to make the docs reproducible,
		// remove to use default srand(time())

#define M 11	// number of data points

#define xCell(c,r) c ## r
#define Range(c1,r1,c2,r2) xCell(c1,r1):xCell(c2,r2)

#define X Range(a,1,a,M)	// X,Y input data points
#define Y Range(b,1,b,M)

#define Z Range(c,1,c,M)	// approximation Y values
#define Err  Range(d,1,d,M)	// absolute error
#define Coef f1:f2		// exponential coefficients

fill X 0, 4/(M-1); a0:f0 = { "X", "Y", "Est", "Err", "maxErr", "Coef" };

Y = { 2*exp(0.5*a1) + 4*drand()-2 }; // y = 2*exp(0.5*x) + noise

Z = { $f$1*exp($f$2*a1) }; Err  = { fabs(R[]C[-2]-R[]C[-1]) };

e1 = max(Err); // e1 is the max-absolute-error to be minimized

{Coef} = search(e1,3,0.3); // start search with f1=3, f2=0.3

// evaluate Y just once here,
// otherwise search() will evaluate it multiple times
//
eval Y;

// eval symbols will evaluate search() which will evaluate everything else
//
eval symbols; plot "exp.out" Range(a,1,c,M);

format "%10.6f"; format A "%5.2f"; print all;

% cat exp.sh
#! /bin/sh

SS -t exp -H exp.ss > exp.html

printf "set term gif\nset output\nplot 'exp.out' using 1:2 with points notitle,\
 'exp.out' using 1:3 with lines notitle\n" | gnuplot > exp.gif

rm exp.out

# the points are the original noisy data,
# the line is the exponential approximation to the data

% ./exp.sh

exp.html