Curve Fitting 3
 Volume Number: 1 Issue Number: 13 Column Tag: Forth Forum

# Curve Fitting, Part III

By Jörg Langowski, Chemical Engineer, Grenoble, France, MacTutor Editorial Board

Last time I promised you to add some input/output to the curve fitter program so that it becomes a useful application. We'll do this here today.

A lot of basic graphic routines which you would have to write yourself on other machines are included in the Quickdraw ROM. This makes life very easy; we can use the screen like we would use a plotter (of course, we can do much better than that, but this is what we are dealing with today). What is not included in the ROM, of course, are utilities that plot coordinate axes, draw lines through pairs of x/y points, add symbols etc.

Most computer systems that are used for numerical calculation purposes contain a set of standard plotting routines, mostly callable from FORTRAN, that will help you generating a 'nice' graphical output with not too much effort. This column gives you some similar routines in FORTH; I have tried to stick closely to the Calcomp plotting routines, since they are something of a standard.

Our objective is to produce something like the curve in Fig. 1. The data points that the model is fitted to will be displayed using some symbols, while the fitted curve is plotted as a smooth line. X- and Y-axes are added to the plot.

The numbers that are to be plotted will be given in floating point format (in this case, single precision arrays as defined previously). So first of all we have to know how to scale the data to the integer values of the bit map that we are going to plot to.

The Calcomp convention is, for an array of N floating point numbers, to store the minimum of these numbers in cell N+1 of the array and the difference between the minimum and the maximum in cell N+2. These two numbers are then rounded to one significant figure. Given this information, one can write routines that automatically draw an axis that spans the range of the data values.

The rounding to one significant figure is done by the word next.int, which given the address of an extended number on the stack returns with the address of the rounded number (which is a local variable of next.int).

fscale finds the minimum and maximum values of an array and stores the scaling numbers above the last data point after rounding them to one significant digit.

xaxis and yaxis draw coordinate axes in x- and y-direction. Input parameters are (#ticks\length\array\npts), #ticks gives the number of ticks to be drawn (to the bottom of the x- and to the left of the y-axis), length is the total axis length in points, array the address of the scaled data array and npts the number of points in array.

Given two scaled arrays, line draws a line through the (x y) pairs defined by the two arrays. The symbol defined by the global symbol is plotted at each of the data points. Symbols could be e.g. '+' or '*'. Setting symbol to zero will draw no symbols (how could you have guessed). Depending on the setting of the global flag connecting, the symbols are either connected by a line or not.

For all the plotting routines cartesian has to be switched on and the origin to be defined within the output window by xyoffset. init.plot is used to do this job.

The main curve fitter program now consists of the following:

- the data arrays xdat and ydat are initialized with the simulated 'experimental data' (init). In this month's example, we use a sum of two exponentials as our model; this gives us five parameters to fit.

- a routine is called in which one can deliberately change the values of the parameters (for this we need floating point input, see below), so that one can see how the correct curve is fitted through the data starting from wrong parameters. One may also change the total number of parameters to be fitted; for instance, fitting three parameters only will force a single exponential fit to the data points. In this case, parameters 4 and 5 will have to be set to zero.

- the main loop calculates theoretical function values from xdat and the parameters and stores them in zdat. The data arrays are scaled and the plot displayed (like in Fig. 1); then one iteration of the fitting routine is taken, the parameters changed accordingly, displayed and the process repeated until parameter changes are below 1 part in 104.

## Floating point input

As I mentioned, the program would not be very useful without floating point input routines. In order to be independent of any particular Forth implementation, I am including a simple floating point input routine here which accepts a string and converts it to the decimal string format used by the SANE conversion routines.

Any such routine would mainly consist of taking successive characters from the input string, generating a decimal mantissa and keeping track of the number of digits behind the decimal point, then looking for an 'e' or 'E' to indicate an exponent and converting the exponent finally. In that aspect, the routine written here is very similar to the MacForth floating input routine. There is one difference, in that numbers without an exponent will be accepted and converted to floating point numbers, too. Therefore, you won't be able to use this routine as an extension of the Forth interpreter, as MacForth level 2 does. But in entering a lot of data, it can be very tedious always having to type an exponent.

fnumber takes a string as its input parameter and leaves on the stack:

```( addr of float\true ) if a valid floating point number was read from
the string,
( false ) if the conversion was not successful.
```

In the latter case, an error message is printed.

input.float reads a string from the keyboard into pad, then calls fnumber to convert it.

With the additions from this column, the curve fitter should finally be a utility that is useful to you (in case you have any curves to fit). Changing the function to be fitted is easy, and you might even install a make switch for vectored execution (V1#7) so that you can easily switch between different functions within one program.

Data input still has to be done manually, number by number. However, input.float may easily be extended to read input from a file. To transfer data to/from other applications through the scrap requires some more work; I'll deal with that problem soon.

## Feedback dept. Re: finding object code of unnamed tokens

The procedure to decompile the object code of an axed token is the following (V1#2):

- convert the token to an absolute address, using token>addr. If the word contained there is \$4e4f (TRAP \$F), the next words will be Forth code. Start decompiling at the following word. Of course, different versions of MacForth will give different addresses for the individual words due to different dictionary arrangements; but this procedure should work for any version of Level 1 or Level 2.

## Re: MacModula floating point

Shortly after my comment on errors in MacModula 2's 32-bit floating point arithmetic, I received a letter from the author of those routines, Daan Strebe:

" A few weeks ago I followed several bug reports to find two fundamental errors in the floating-point routines, one a conceptual error in the rounding and the other a misunderstanding about the 68000 processor instruction set. These two errors affected the multiply routine most, although the others were somewhat affected also. I revamped those routines and then ran a complex set of comprehensive tests, and the results allow me to state with a margin of confidence that the floating point in the latest rev is now not only slightly faster than it was, but also that the only errors in the basic routines (not necessarily including those in the math library) are from rounding. The rounding as it is now implemented yields 3 downward, 4 upward, and 1 non-round per 8 random floating-point operations. This should be satisfactory for most users; full IEEE rounding would considerably decrease the speed. I believe the compromise was successful. "

So everything should be fixed now. Let's hope that the new update will be distributed soon (it probably is when this goes to the printer's).

```Listing 1: Plotting routines and floating point input for the curve fitter
program

( Additions to the curve fitting program)
( for graphical output and floating point input.)
( © 1985 J. Langowski for MacTutor)
( Again, only the parts that have been changed)
( or added with respect to the last two columns )
( are printed here)

: s> 1008 fp68k @sr 10 and ;

: numstring create 24 allot ;
numstring zzs1 ( internal conversion string )

: dec. ( float\format# -- )
zzformat ! zzformat swap zzs1 b2d
zzs1 dup w@ 255 > if ." -" else ."  " then
dup 4+ count over 1 type ." ."
swap 1+ swap 1- type ( mantissa )
2+ w@ ( get exponent )
1 w* zzformat @ + 1-
." E" 0 .r ;

create xdat 400 allot   create ydat 400 allot
create zdat 400 allot   create residues 400 allot
100 10 matrix derivmat    10 11 matrix resmat
5 constant npars        80 constant npts
create par 5 10 * allot
create delta 5 4*  allot

float logten  ten logten x2x  logten lnx
: log10x dup lnx logten swap f/ ;

( define function )                               ( 090485 jl )
float temp ( local to func)

func ( x -- f[x] = par[1] + par[2] * exp[par[3]*x]
+ par[4] * exp[par[5]*x] )
dup temp x2x
par 40 + temp f* temp expx  par 30 + temp f*
par 20 + over f* dup expx  par 10 + over f*
par over f+  temp over f+ ;
axe temp
: >fa1  fa1 s2x ;
: init_pars   one par x2x  two par 10+ x2x
-one par 20 + x2x  two par 20 + f/
one par 30 + x2x
-one par 40 + x2x ten par 40 + f/ ;
init_pars

: one_iter
make_derivmat  residuals
make_resmat     delta 0 0 resmat npars gauss ;

: new_pars 16 ( true if no significant changes)
npars 0 do par i 10 * +
delta i 4* +  over  s+
delta i 4* + fa1 s2x  fa1 f/
fa1 fabs  eps fa1 f> and loop ;

80 ' npts !    5 ' npars !

: init ( initialize data arrays)
npts 0 do i sp@ fa1 in2x 4*
xdat over + fa1 swap x2s
ydat over + fa1 func  ranf fa2 x2x
ten fa2 f/  fa2 over f+  swap x2s
drop  loop ;

( plotting routines)                        ( 092385 jl )
float 1/2  1 sp@ 1/2 in2x drop  2 sp@ 1/2 in/ drop
float small  ten small x2x  -200 sp@ fa1 in2x drop
fa1 small x^y
( anything smaller than small will be zero)
float sc.aux  float sc.exp
: next.int ( float -- rounded to 1 dec. place )
dup small f>
if dup sc.aux x2x  sc.aux  dup fabs  dup log10x
dup 1/2 swap f- frti
ten sc.exp x2x  sc.aux sc.exp x^y   sc.aux x2x
sc.exp sc.aux f/  sc.aux frti  sc.exp sc.aux f*
sc.aux
else drop zero then  ;

float xlow  float xhi  float sc.factor
: fscale
( array \ n -- | start, scale -> array[n+1..n+2] )
over xlow s2x  over xhi s2x
over over 1 do
dup i 4* + dup xlow s> not if dup xlow s2x then
dup xhi  s> if xhi s2x  else drop then  loop
xlow xhi f-
over 4*  +  xlow next.int swap x2s
1+ 4* +  xhi  next.int swap x2s  ;

: .fscale ( array \ n -- )
4* + dup fa1 s2x ." min = " fa1 7 dec.
4+ fa1 s2x ." , scale = " fa1 7 dec. cr ;

( xtick, ytick )                                  ( 092385 jl )

: xtick ( rx x -- ) dup 0 move.to
dup -5 draw.to  dup 15 - -16 move.to
get.textsize >r 9 textsize swap 2 dec. r> textsize
0 move.to ;

: ytick ( ry y -- ) -45 over 4- move.to
get.textsize >r 9 textsize swap 2 dec. r> textsize
0 over move.to    -5 over draw.to
0 over 6- move.to  0 swap draw.to ;

( xaxis)                                           ( 092385 jl )
float start  float del
: xaxis
( #ticks\length\array\npts -- | starts at origin )
0 0 move.to
4* + dup  start s2x  4+ delta s2x
over /  ( #ticks\length/tick -- )
over 0 do dup i 1+ * dup 0 draw.to
delta fa1 x2x i 1+ sp@ fa1 in* drop
3 pick sp@ fa1 in/ drop
start fa1 f+    fa1 swap xtick loop
drop drop ;

( yaxis)                                          ( 092385 jl )
: yaxis
( #ticks\length\array\npts -- | starts at origin )
0 0 move.to
4* + dup  start s2x  4+ delta s2x
over /  ( #ticks\length/tick -- )
over 0 do dup i 1+ * 0 over draw.to
delta fa1 x2x i 1+ sp@ fa1 in* drop
3 pick sp@ fa1 in/ drop
start fa1 f+    fa1 swap ytick loop
drop drop ;

( line, variables )                               ( 092585 jl )
variable xline.length variable yline.length
variable xpos  variable ypos  variable symbol
variable connecting  connecting off
float xstart float ystart float xdel float ydel

( line)                                           ( 092585 jl )
: line ( xarray\yarray\npts\xlength\ylength -- )
yline.length !  xline.length !   0 0 move.to
over over 4* + dup ystart s2x 4+ ydel s2x
3 pick over 4* + dup xstart s2x 4+ xdel s2x
0 do over i 4* + fa1 s2x  xstart fa1 f-  xdel fa1 f/
dup  i 4* + fa2 s2x  ystart fa2 f-  ydel fa2 f/
xline.length fa1 in*  fa1 xpos x2in
yline.length fa2 in*  fa2 ypos x2in
xpos @  ypos @  over over
connecting @ if draw.to else move.to then
-4 -4 rmove symbol @ emit move.to
loop   drop drop ;

( calc.theor, disp.theor, disp.data, scale.all)
( 092585 jl )
: calc.theor
npts 0 do xdat i 4* + fa1 s2x fa1 func
zdat i 4* + x2s  loop ;
: disp.theor
0 symbol !  connecting on
xdat zdat npts 400 250 line ;
: disp.data
43 symbol !  connecting off
xdat ydat npts 400 250 line ;
: scale.all
xdat npts fscale  ydat npts fscale
ydat npts 4* + @  zdat npts 4* + !
ydat npts 1+ 4* + @  zdat npts 1+ 4* + !  ;

( disp.axes)                                      ( 092585 jl )
: disp.axes
5 400 xdat npts xaxis  5 250 ydat npts yaxis ;
: init.plot
page 50 255 xyoffset cartesian on ;

( floating point input)                           ( 092385 jl )
: fsign zzs1 ;  : fexpo zzs1 2+ ;
: fmant zzs1 4+ ;  variable frac   float float.out
variable decimals
: input.sign 0 fsign !
dup c@ case 45 of -1 fsign w! 1+ endof
43 of             1+ endof endcase ;
: input.mantissa   0 decimals !   0 frac !
fmant 21 +  fmant 1+ do  i fmant 1+ -  fmant c!
dup c@ dup
case 48 57 range.of ic! frac @ decimals +!
1+ 1 endof
46 of   drop 1 frac !    1+ 0 endof
20 swap endcase  +loop drop ;
: input.exponent
dup c@ dup bl = not
if dup 69 =  swap 101 =  or not
if  drop 0
else dup 1+ c@ 43 = 1 and +
number decimals @ - fexpo w! -1
then
else  decimals @ -1 * fexpo w! -1  2drop
then  ;

: fnumber zzs1 24 blanks
if fsign float.out d2b float.out -1
else ." floating input error!" cr 0 then ;

: get.pars
npars 0 do
begin cr ." par[" i 1+ . ." ] = " input.float until
par i 10 * + x2x loop
cr begin cr ." total # of parameters to be fitted "
5 input.number until
' npars !  cr ;

: .pars npars 0 do
." par[" i 1+ . ." ] = " i 10 * par + 7 dec. cr loop ;

( main curve fitter program)                 ( 092685 jl )
: fit.curve  page
init_pars ." initializing data arrays..." cr init
5 ' npars ! init.plot get.pars .pars
." fitting " npars . ." parameters to "
npts . ." data points "
calc.theor scale.all page
disp.data disp.axes disp.theor
begin  one_iter new_pars not while
page .pars
calc.theor page disp.data disp.axes disp.theor
repeat
300 300 move.to ;
```

Community Search:
MacTech Search:

Capto 1.2.9 - \$29.99
Capto (was Voila) is an easy-to-use app that takes capturing, recording, video and image editing to the next level. With an intelligent file manager and quick sharing options, Capto is perfect for... Read more
Opera 51.0.2830.40 - High-performance We...
Opera is a fast and secure browser trusted by millions of users. With the intuitive interface, Speed Dial and visual bookmarks for organizing favorite sites, news feature with fresh, relevant content... Read more
GarageSale 7.0.13 - Create outstanding e...
GarageSale is a slick, full-featured client application for the eBay online auction system. Create and manage your auctions with ease. With GarageSale, you can create, edit, track, and manage... Read more
1Password is a password manager that uniquely brings you both security and convenience. It is the only program that provides anti-phishing protection and goes beyond password management by adding Web... Read more
Evernote 7.0.1 - Create searchable notes...
Evernote allows you to easily capture information in any environment using whatever device or platform you find most convenient, and makes this information accessible and searchable at anytime, from... Read more
MacUpdate Desktop 6.2.0 - \$20.00
MacUpdate Desktop brings seamless 1-click app installs and version updates to your Mac. With a free MacUpdate account and MacUpdate Desktop 6, Mac users can now install almost any Mac app on... Read more
HoudahSpot is a versatile desktop search tool. Use HoudahSpot to locate hard-to-find files and keep frequently used files within reach. HoudahSpot will immediately feel familiar. It works just the... Read more
EtreCheck 4.0.4 - For troubleshooting yo...
EtreCheck is an app that displays the important details of your system configuration and allow you to copy that information to the Clipboard. It is meant to be used with Apple Support Communities to... Read more
WhatsApp 0.2.8361 - Desktop client for W...
WhatsApp is the desktop client for WhatsApp Messenger, a cross-platform mobile messaging app which allows you to exchange messages without having to pay for SMS. WhatsApp Messenger is available for... Read more
iClock is a menu-bar replacement for Apple's default clock but with 100x features. Have your Apple or Google calendar in the menubar. Have the day, date, and time in different fonts and colors in the... Read more

## Latest Forum Discussions

The best games like Florence
Florence is a great little game about relationships that we absolutely adored. The only problem with it is it's over a little too soon. If you want some other games with some emotional range like Florence, check out these suggestions: [Read more] | Read more »
Angry Birds Champions adds cash prizes t...
Collaborating with developer Rovio Entertainment, GSN Games has released a twist on the Angry Birds formula. Angry Birds Champions features the same bird-flinging gameplay, but now you can catapult Red and co for cash. | Read more »
Around the Empire: What have you missed...
148Apps is part of a family. A big family of sites that make sure you're always up to date with all the portable gaming news. Just like a real family, I guess. I don't know, my mum never told me anything about Candy Crush to be fair. [Read more] | Read more »
The Battle of Polytopia Guide - Tips for...
The addition of multiplayer to The Battle of Polytopia has catapulted the game from a fun enough time waster to a fully-fledged 4X experience on your phone. We've been playing quite a few matches over the past week or so, and we've put together a... | Read more »
All the best games on sale for iPhone an...
Hi there, and welcome to our round up of all the best games that are on sale for iOS at the moment. It's not a vintage week in terms of numbers, but I'm pretty sure that every single one of these is worth picking up if you haven't already played... | Read more »
Disc Drivin' 2 Guide - Tips for the...
We're all still playing quite a bit of Disc Drivin' 2 over here at 148Apps, and we've gotten pretty good at it. Now that we've spent some more time with the game and unlocked more powerups, check out some of these more advanced tips: | Read more »
Alto's Odyssey Guide - How to Tackl...
Alto’s Odyssey is a completely stunning and serene runner, but it can also be a bit tricky. Check out these to try and keep your cool while playing this endless runner: Don’t focus too much on tasks [Read more] | Read more »
Here's everything you need to know...
Alto's Odyssey is a really, really good game. If you don't believe me, you should definitely check out our review by clicking this link right here. It takes the ideas from the original Alto's Adventure, then subtly builds on them, creating... | Read more »
Alto's Odyssey (Games)
Alto's Odyssey 1.0.1 Device: iOS Universal Category: Games Price: \$4.99, Version: 1.0.1 (iTunes) Description: Just beyond the horizon sits a majestic desert, vast and unexplored. Join Alto and his friends and set off on an endless... | Read more »
Vainglory 5v5: Everything you need to kn...

## Price Scanner via MacPrices.net

Sale! Amazon offers 13″ 2.3GHz MacBook Pros f...
Amazon has 2017 13″ 2.3GHz Apple MacBook Pros on sale today for \$151-\$150 off MSRP, each including free shipping: – 13″ 2.3GHz/128GB Space Gray MacBook Pro (MPXQ2LL/A): \$1148 \$151 off MSRP – 13″ 2.... Read more
Apple AirPods in stock today for \$159, free s...
Adorama reports stock of Apple AirPods today for \$159 including free shipping, plus pay no sales tax outside of NY & NJ. See our Apple AirPod Price Tracker for the latest prices and stock status... Read more
Saturday Sale: Amazon offers 12″ 1.3GHz MacBo...
Amazon has Silver and Gold 2017 12″ 1.3GHz Retina MacBooks on sale for \$250 off MSRP. Shipping is free: – 12″ 1.3GHz Silver MacBook: \$1349.99 \$250 off MSRP – 12″ 1.3GHz Gold MacBook: \$1349.99 \$250... Read more
Use your Apple Education discount and save up...
Purchase a new Mac using Apple’s Education discount, and take up to \$400 off MSRP. All teachers, students, and staff of any educational institution with a .edu email address qualify for the discount... Read more
Apple Canada offers 2017 21″ and 27″ iMacs fo...
Canadian shoppers can save up to \$470 on the purchase of a 2017 current-generation 21″ or 27″ iMac with Certified Refurbished models at Apple Canada. Apple’s refurbished prices are the lowest... Read more
9″ iPads available online at Walmart for \$50...
Walmart has 9.7″ Apple iPads on sale for \$50 off MSRP for a limited time. Sale prices are for online orders only, in-store prices may vary: – 9″ 32GB iPad: \$279.99 \$50 off – 9″ 128GB iPad: \$379.99 \$... Read more
15″ Apple MacBook Pros, Certified Refurbished...
Save \$360-\$420 on the purchase of a 2017 15″ MacBook Pro with Certified Refurbished models at Apple. Apple’s refurbished prices are the lowest available for each model from any reseller. An standard... Read more
Amazon restocks MacBook Pros with models avai...
Amazon has restocked 15″ and 13″ Apple MacBook Pros with models on sale for up to \$251 off MSRP. Shipping is free. Note that stock of some Macs may come and go (and some sell out quickly), so check... Read more
Lowest price of the year: 15″ 2.8GHz Apple Ma...
Amazon has the 2017 Space Gray 15″ 2.8GHz MacBook Pro on sale today for \$251 off MSRP. Shipping is free: – 15″ 2.8GHz Touch Bar MacBook Pro Space Gray (MPTR2LL/A): \$2148, \$251 off MSRP Their price is... Read more
Apple restocks full line of Certified Refurbi...
Apple has restocked a full line of Apple Certified Refurbished 2017 13″ MacBook Pros for \$200-\$300 off MSRP. A standard Apple one-year warranty is included with each MacBook, and shipping is free.... Read more

## Jobs Board

*Apple* Media Products Commerce Engineering...
# Apple Media Products Commerce Engineering Manager Job Number: 56207285 Santa Clara Valley, California, United States Posted: 26-Jan-2018 Weekly Hours: 40.00 **Job Read more
Digital Platforms Lead, Today at *Apple* -...
# Digital Platforms Lead, Today at Apple Job Number: 56178747 Santa Clara Valley, California, United States Posted: 23-Feb-2018 Weekly Hours: 40.00 **Job Summary** Read more
*Apple* Retail - Multiple Positions - Apple,...
Job Description:SalesSpecialist - Retail Customer Service and SalesTransform Apple Store visitors into loyal Apple customers. When customers enter the store, Read more
*Apple* Retail - Multiple Positions - Apple,...
Job Description: Sales Specialist - Retail Customer Service and Sales Transform Apple Store visitors into loyal Apple customers. When customers enter the store, Read more
*Apple* Retail - Multiple Positions - Apple,...
Job Description:SalesSpecialist - Retail Customer Service and SalesTransform Apple Store visitors into loyal Apple customers. When customers enter the store, Read more