TweetFollow Us on Twitter

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
   input.sign input.mantissa input.exponent
   if fsign float.out d2b float.out -1
   else ." floating input error!" cr 0 then ;

: input.float pad 22 blanks pad 22 expect pad fnumber ;

: 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:

Software Updates via MacUpdate

Trend Micro 11.0.2062 - An essential sec...
Trend Micro Antivirus provides essential security for macOS with real-time malware detection and mitigation in an affordable solution with a simple, intuitive interface. However, be aware that more... Read more
Backblaze 7.0.2.490 - Online backup serv...
Backblaze is an online backup service designed from the ground-up for the Mac. With unlimited storage available for $6 per month, as well as a free 15-day trial, peace of mind is within reach with... Read more
Pro Video Formats 2.2.1 - Updates for pr...
Pro Video Formats includes support for the following professional video codecs: Apple Intermediate Codec Apple ProRes AVC-Intra 50 / 100 / 200 / 4:4:4 / LT AVC-LongG XAVC XF-AVC DVCPRO HD HDV XDCAM... Read more
Boom 3D 1.3.11 - $19.99
Boom 3D is a revolutionary app with 3D Surround Sound and phenomenally rich and intense audio that is realistic and works on any headphones. Features 3D surround sound Built-in audio player... Read more
Final Cut Pro 10.5.2 - Professional vide...
Redesigned from the ground up, Final Cut Pro combines revolutionary video editing with a powerful media organization and incredible performance to let you create at the speed of thought.... Read more
Chromium 89.0.4389.72 - Fast and stable...
Chromium is an open-source browser project that aims to build a safer, faster, and more stable way for all Internet users to experience the web. Version 89.0.4389.72: Complete list of changes can... Read more
iMovie 10.2.3 - Edit personal videos and...
With a streamlined design and intuitive editing features, iMovie lets you create Hollywood-style trailers and beautiful movies like never before. Browse your video library, share favorite moments,... Read more
OmniPlan 4.2.5 - Professional-grade proj...
With OmniPlan, you can create logical, manageable project plans with Gantt charts, schedules, summaries, milestones, and critical paths. Break down the tasks needed to make your project a success,... Read more
Compressor 4.5.2 - Adds power and flexib...
Compressor adds power and flexibility to Final Cut Pro X export. Customize output settings, work faster with distributed encoding, and tap into a comprehensive set of delivery features. Features:... Read more
Motion 5.5.1 - Create and customize Fina...
Motion is designed for video editors, Motion 5 lets you customize Final Cut Pro titles, transitions, and effects. Or create your own dazzling animations in 2D or 3D space, with real-time feedback as... Read more

Latest Forum Discussions

See All

Frogger in Toy Town's latest update...
Konami Digital Entertainment has announced today that their Apple Arcade title Frogger in Toy Town has been updated to introduce a new ranked Endurance Mode. This new game variant's arrival is also accompanied by a few other changes. [Read more] | Read more »
Mitoza is surreal adventure toy you can...
The folks behind the Rusty Lake games have just put a new title onto the App Store. Second Maze, Rusty Lake's collaborative publishing brand, has just brought this 10 year old adventure game from Gal Mamalya to mobile. The best part about all of... | Read more »
Pocket Gamer Awards 2021: You have five...
Three weeks ago our sister site, Pocket Gamer, entered the voting stage for the upcoming Pocket Gamer Awards 2021 and is now in the final stretch. You only have a few hours left to vote for the games you’ve enjoyed on mobile in the past year, as... | Read more »
Patty Stack is a casual arcade game, ava...
Patty Stack is a casual arcade title that's available now for iOS and Android. It's the debut game from developer Feeka Games tasks players with making an increasingly giant burger tower. Think of it as Tower Bloxx but more edible. [Read more] | Read more »
Distract Yourself With These Great Mobil...
There’s a lot going on right now, and I don’t really feel like trying to write some kind of pithy intro for it. All I’ll say is lots of people have been coming together and helping each other in small ways, and I’m choosing to focus on that as I... | Read more »
Genshin Impact Guide - Gacha Strategy: W...
This is part 2 of our Genshin Impact gacha strategy guides. See part 1 here. You can check out more guides for Genshin Impact here. | Read more »
Slashy Camp is a new endless runner insp...
Blue Wizard Digital has released Slashy Camp onto iOS and Android after it spent a short amount of time in early access. [Read more] | Read more »
Kinder World is a relaxing game about lo...
Lumi Interactive is releasing a game called Kinder World later this year on iOS and Android, which is all about looking after houseplants. [Read more] | Read more »
Steam Link Spotlight - Fights in Tight S...
Steam Link Spotlight is a feature where we look at PC games that play exceptionally well using the Steam Link app. Our last entry was on Hades. Read about how it plays using Steam Link over here. | Read more »
Lyxo, the light-based puzzler for mobile...
Vienna-based independent game studio Emoak has just released its unique light-based puzzler for iOS and Android. Founded in 2014 by Tobias Sturn, the company is also the creative force behind the infinite climbing game Paper Climb, as well as the... | Read more »

Price Scanner via MacPrices.net

Weekend Sale: $100 off Apple iPad Magic Keybo...
Amazon has Apple iPad Magic Keyboards on sale for $100 off MSRP for a limited time. Amazon’s prices are the lowest available for iPad Magic Keyboard from any Apple reseller this weekend: – Magic... Read more
Gazelle now offering a full line of refurbish...
Gazelle is now offering a full range of discounted, refurbished, unlocked Apple iPhone 12 models starting at $649. iPhones are offered in Fair, Good, and Excellent conditions, and multiple colors are... Read more
These are the latest discounted iPhones Apple...
Apple has a range of Certified Refurbished iPhones available right now starting at only $339. Apple includes a standard one-year warranty, new outer shell, and shipping is free. According to Apple, “... Read more
Save up to $64 on new M1 MacBook Airs at Expe...
Apple reseller Expercom has 2020 13″ M1 MacBook Airs on sale for $51-$64 off Apple’s MSRP with prices starting at $947.96. In addition to their MacBook Air sale prices, take $50 off AppleCare+ when... Read more
Discounts available on 16″ MacBook Pros with...
Upgrade a 16″ 6-Core or 8-Core MacBook Pro from 16GB of standard RAM to 32GB at Adorama, and save $100-$210 over Apple’s price for this custom option: – 16″ 6-Core MacBook Pro/32GB RAM: $2699, save $... Read more
10.9″ iPad Airs on sale for $50-$70 off Apple...
Amazon has new 2020 10.9″ Apple WiFi iPad Airs in stock and on sale today for up to $70 off MSRP with prices starting at $549. Note that Amazon’s sale price might be restricted to certain colors (see... Read more
Apple restocks 2020 27″ 5K iMacs for up to $3...
After an initial offering in January, Apple has restocked a full line of Certified Refurbished 2020 27″ 5K iMacs starting at $1529 and up to $350 off original MSRP. Apple’s one-year warranty is... Read more
Sale! 16″ 8-Core MacBook Pro for $2449, $350...
Apple reseller Adorama has the 16″ 2.3GHz 8-Core Space Gray MacBook Pro in stock and on sale today for $2449 including free shipping. Their price is $350 off Apple’s MSRP for this model, and it’s the... Read more
Roundup of 13″ Multi-Core Intel MacBook Pro s...
Apple resellers are offering significant sales & deals this week on 2020 13″ MacBook Pros with 10th generation Intel CPUs. Take up to $250 off Apple’s MSRP, get free fast shipping, and/or pay no... Read more
64GB iPhone 8 Plus available for $379 at Appl...
Apple has the 64GB iPhone 8 Plus in Space Gray & Gold colors available for $379 today, Certified Refurbished. Each phone is unlocked and comes with Apple’s standard 1-year warranty and free... Read more

Jobs Board

Geek Squad Advanced Repair *Apple* Professi...
**795178BR** **Job Title:** Geek Squad Advanced Repair Apple Professional **Job Category:** Store Associates **Store Number or Department:** 001406-Allen Park-Store Read more
Geek Squad *Apple* Consultation Professiona...
**796549BR** **Job Title:** Geek Squad Apple Consultation Professional **Job Category:** Store Associates **Store Number or Department:** 001800-Hot Springs-Store Read more
*Apple* Mobility Specialist - Best Buy (Unit...
**796014BR** **Job Title:** Apple Mobility Specialist **Job Category:** Store Associates **Store Number or Department:** 001776-Woodmore Towne Centre-Store **Job Read more
Systems Architect, *Apple* Production Engin...
…package beginning on your first day? If so, we hope you'll keep reading! The Apple Sales Engineering and account team is looking for a stellar presales engineer with Read more
Systems Engineer, Webscale, *Apple* Retail,...
…beginning on your first day? If so, we hope you'll keep reading! The Apple Sales Engineering team is looking for a pre-sales engineer with Enterprise engineering Read more
All contents are Copyright 1984-2011 by Xplain Corporation. All rights reserved. Theme designed by Icreon.