Curve fitting with polynomials

Post Reply
Henko
Posts: 814
Joined: Tue Apr 09, 2013 12:23 pm
My devices: iPhone,iPad
Windows
Location: Groningen, Netherlands
Flag: Netherlands

Curve fitting with polynomials

Post by Henko »

This is a remake of an curve fitting program which i made a couple of years ago. I used a better programming technique (event driven MAIN,like in Windows), and some presently available interface objects in SB. As a result, the user interface is far better than the earlier design.
After running the program, please hit the "help" button for operating info.

Code: Select all

' curve fitting with polynomials in 2D
'
np_max=200   ' maximum # of points, may be modified
m_max=9      ' maximum degree of fitting polynomial
f$="ls_data" ' file used to store and restore point data
option base 1
dim xy(np_max,2),xlist$(np_max),ylist$(np_max)
dim coeff(m_max+1),coef$(m_max+1)
'
prog_init()
do slowdown
  if b_p("new_x") then add_new_X()
  k=list_selected("xp") ! if k>0 and k<=.np then edit_X(k)
  k=list_selected("yp") ! if k>0 then edit_Y(k)
  if b_p("del_p") then del_point()
  if b_p("degree") then set_degree()
  if b_p("calc") then calc()
  if b_p("help") then help()
  if b_p("load") then load_data(.f$)
  if b_p("save") then save_data(.f$)
  if b_p("stop") then stop
until forever
end

def prog_init()
graphics ! graphics clear ! draw color 0,0,0
get screen size sw,sh ! cx=sw/2
.gx1=20 ! .gy1=500 ! .gx2=748 ! .gy2=908
.np=0 ! help_page()
c_list("xp","X-coordinates",.xlist$,.np,20,20,190,420)
c_list("yp","Y-coordinates",.ylist$,.np,210,20,380,420)
button "new_x" text "Add X-value" at 20,440 size 170,30
button "del_p" text "Delete point" at 210,440 size 170,30
field "points" text "    "&"0 Points" at 450,20 size 170,30 RO
field "points" font name "Courier-Bold"
button "calc" text "y = f(x)" at 450,440 size 290,30
button "degree" text "Degree" at 450,70 size 170,30
button "load" text "Load data" at 240,sh-50 size 120,30
button "save" text "Save data" at 400,sh-50 size 120,30
button "help" text "Help" at 20,sh-50 size 120,30
button "stop" text "Quit" at sw-140,sh-50 size 120,30
button "calc" hide
init_numpad(450,120,60,.7,.7,.7,1)
end def

def add_new_X()
if .np=.np_max then return
.np+=1 ! .xy(.np,2)=0 ! edit_X(.np)
field "points" text "    " & .np & " Points"
list "xp" select .np ! list "xp" select -1
list "yp" select .np ! list "yp" select -1
end def

def edit_X(k)
.xy(k,1)=numpad(0,0)
tablesort (.np,2,.xy,1)
disp_lists() ! plot()
end def

def edit_Y(k)
.xy(k,2)=numpad(0,0)
for i=1 to .np ! .ylist$(i)=i & ":  " &.xy(i,2) ! next i
c_list("yp","Y-coordinates",.ylist$,.np,210,20,380,420)
plot()
end def

def del_point()
k=numpad(1,.np)
if k<.np then
  for i=k to .np-1
    .xy(i,1)=.xy(i+1,1) ! .xy(i,2)=.xy(i+1,2)
    next i
  end if
.np-=1
field "points" text "    " & .np & " Points"
disp_lists() ! plot()
end def

def disp_lists()
for i=1 to .np
  .xlist$(i)=i & ":   " &.xy(i,1) ! .ylist$(i)=i & ":   " &.xy(i,2)
  next i
c_list("xp","X-coordinates",.xlist$,.np,20,20,190,420)
c_list("yp","Y-coordinates",.ylist$,.np,210,20,380,420)
end def

def stat()
.xmin=1E+8 ! .xmax=-1E+8 ! .ymin=1E+8 ! .ymax=-1E+8
for i=1 to .np
  x=.xy(i,1) ! .xmin=min(.xmin,x) ! .xmax=max(.xmax,x)
  y=.xy(i,2) ! .ymin=min(.ymin,y) ! .ymax=max(.ymax,y)
  next i
end def

def set_degree()
dim x(.np),y(.np)
if .np<2 then return
.m=numpad(0,min(.np-1,.m_max))
button "degree" text "Degree = "&.m
for i=1 to .np ! x(i)=.xy(i,1) ! y(i)=.xy(i,2) ! next i
ls (.np,.m,x,y,.coeff)
button "calc" text "y = f(x)" ! button "calc" show
c_table() ! plot()
end def

def calc()
x=numpad(.xmin,.xmax)
y=poly(.m,.coeff,x)
button "calc" text "x="&x&"   =>   "&"y="&y
end def

def load_data(f$)
i=0
do
 i+=1 ! file f$ input .xy(i,1),.xy(i,2)
 until not data_exist(f$)
.np=i ! file f$ reset
field "points" text "    " & .np & " Points"
disp_lists() ! plot()
end def

def save_data(f$)
file f$ delete
for i=1 to .np ! file f$ print .xy(i,1),.xy(i,2) ! next i
end def

def plot()
stat() ! if .xmin=.xmax or .ymin=.ymax or .np<2 then return
if .xmin>0 then
  sc_x=720/.xmax ! y_as=.gx1+4 ! goto lab1
  end if
if .xmax<0 then
  sc_x=720/-.xmin ! y_as=.gx2-4 ! goto lab1
  end if
sc_x=720/(.xmax-.xmin) ! y_as=.gx1+4-.xmin*sc_x
lab1:
if .ymin>0 then
  sc_y=400/.ymax ! x_as=.gy2-4 ! goto lab2
  end if
if .ymax<0 then
  sc_y=400/-.ymin ! x_as=.gy1+4 ! goto lab2
  end if
sc_y=400/(.ymax-.ymin) ! x_as=.gy2-4+.ymin*sc_y
lab2:
refresh off
fill color 0.88,0.85,0.65 ! fill rect .gx1,.gy1 to .gx2,.gy2
draw rect .gx1,.gy1 to .gx2,.gy2
draw color .7,.7,.7 ! draw size 1    ' draw grid
for x=.gx1+4 to .gx2-4 step 20
  draw line x,.gy1+4 to x,.gy2-4
  next x
for y=.gy1+4 to .gy2-4 step 20
  draw line .gx1+4,y to .gx2-4,y
  next y
draw color 0,0,1 ! draw size 3
draw line y_as,.gy1+4 to y_as,.gy2-4   ' draw axes
draw line .gx1+4,x_as to .gx2-4,x_as
draw color 0,0,0
for i=1 to .np
  x=xp(.xy(i,1)) ! y=yp(.xy(i,2))
  draw circle x,y size 2
  next i
if .m>0 then
  draw color 1,0,0 ! draw size 2
  draw to y_as+sc_x*.xmin,x_as-sc_y*poly(.m,.coeff,.xmin)
  dx=(.xmax-.xmin)/48
  for x=.xmin+dx to .xmax step dx
    draw line to y_as+sc_x*x,x_as-sc_y*poly(.m,.coeff,x)
    next x
  end if
refresh
draw color 0,0,0
end def

def xp(x)=plot.y_as+plot.sc_x*x
def yp(y)=plot.x_as-plot.sc_y*y

def c_table()
for i=1 to .m+1 ! .coef$(i)="a("&(i-1)&") = "&.coeff(i) ! next i
c_list("cf","Coefficients",.coef$,.m+1,450,120,620,420)
end def

def poly (n,coef(),x)
res=coef(n+1)
for i=n to 1 step-1
  res=res*x+coef(i)
next i
return res
end def

def n2a$(num,lang,dec)
b$="               "
fac=10^dec
num$=floor(fac*num)/fac
tot=lang-len(num$)
if tot<1 then tot=1
rem test ("", tot)
a$=substr$(b$,1,tot) & num$
n2a$=a$
end def


def help()
n$="help" ! page n$ set ! page n$ alpha 0 ! page n$ show
for i=0 to 1 step 0.01 ! page n$ alpha i ! pause 0.01 ! next i
do slowdown ! until b_p("h_stop")
for i=1 to 0 step -0.01 ! page n$ alpha i ! pause 0.01 ! next i
page n$ hide ! page "" set
end def

def help_page()
xs=40 ! ys=200 ! w=680 ! h=680
h$="Curve fitting with polynomials in 2 dimensions."&chr$(10)&"A set of x-y points can be fitted with a polynomial of a suitable degree."&chr$(10)&"In this program a maximum number of 200 point is set, but this can easily be changed in the dim function."&chr$(10)&"The maximum degree of the polynomial is set at 9, as a higher degrees seldom add to the quality of the fit."&chr$(10)&"x- and y-coordinates are kept in separate lists with their point number (for synchronizing purpose). Both lists are kept in sequence of the x-coordinates."&chr$(10)&"Points are added by touching the ""new x"" button. the corresponding y-value is defaulted to zero."&chr$(10)&"Existing x- and y-values may be edited bij touching them in the appropriate list."&chr$(10)&"Points may be deleted by touching the ""delete point"" button and entering the point number."&chr$(10)&"The collection of points can be saved to a file and loaded later on by touching the ""save"" and ""load"" buttons."&chr$(10)&"A graph will appear as soon as two points are entered with different x- and y-values. Each modification in the point data is right away reflected in the graph."&chr$(10)&"Also, when more than one point is present, a fit command may be issued by touching the ""degree"" button and entering the degree of the fitting polynomial."&chr$(10)&"The coefficients of the polynomial are then shown and the calculated curve is drawn in the graph."&chr$(10)&"An y-value may be calculated for a given x-value by using the ""y=f(x)"" button."
n$="help" ! page n$ set ! page n$ frame xs,ys,w,h
field "hh" text h$ at 0,0 size w,h ML RO
field "hh" font color 0,0,1 ! field "hh" font size 20
button "h_stop" text "Close" at w/2-30,h-40 size 60,25
page n$ hide
page "" set
end def

    ' {LS-utils}

' id$ = object name
' cont$ = array met elementen
' size = aantal elementen in de list
'
def c_list(id$,title$,cont$(),size,xt,yt,xb,yb)
size=max(1,size)
dim temp$(size)
for i=1 to size ! temp$(i)=cont$(i) ! next i
list id$ text temp$ at xt+2,yt+32 size xb-xt-4,yb-yt-34
draw size 3
draw rect xt,yt to xb,yb ! draw line xt,yt+30 to xb,yt+30
draw color 0,0,1 ! draw text title$ at xt+5,yt+5
draw color 0,0,0
end def

' numeric table sorting function (option base 1 assumed)
' r = number of rows
' c = number of columns
' in(,) = table to be sorted
' col = column to be used for sorting
'
def tablesort (r,c,in(,),col)
dim sortcol(r),index(r),out(r,c)
for i=1 to r ! sortcol(i)=in(i,col) ! next i
sort sortcol as index
for i=1 to r ! for j=1 to c
  out(i,j)=in(index(i),j)
  next j ! next i
for i=1 to r ! for j=1 to c ! in(i,j)=out(i,j) ! next j ! next i
return
end def

' numerical keypad object
' 
' produce a simple keypad to quickly enter a number in an app
' upon entry, the keypad disappears
' initialize once, multiple use after
' left upper corner is placed at "xtop,ytop"
' "bs" is the button size (keypad becomes 4.3 times larger)
' size of number is accepted between "minval" and "maxval"
' if both "minval" and "maxval" are zero, then no restrictions
' max number of tokens in the number is 10 (minus and dot included)
' works for option base 0 and 1
'
def init_numpad(xtop,ytop,bs,R,G,B,alpha)
name$="numpad" ! cn=10
page name$ set 
page name$ frame xtop,ytop,0,0
set buttons custom
if bs<20 then bs=20
sp=4 ! th=.5*bs+4 ! ww=4*bs+5*sp ! hh=th+4*bs+6*sp
fsize=.5*bs
draw font size fsize ! set buttons font size fsize
draw color 1,1,1 ! fill color .7,.7,.7
button "rec" title "" at 0,0 size ww,hh
button "res" title "" at 0,0 size ww,th+4
fill color R,G,B ! fill alpha alpha
button "0" title "0" at sp,th+3*bs+5*sp size bs,bs
for k=1 to 9
  x=(k-1)%3 ! y=2-floor((k-1)/3)
  button k title k at (x+1)*sp+x*bs,th+y*bs+(y+2)*sp size bs,bs
  next k
button "-" title "-" at 2*sp+bs,th+3*bs+5*sp size bs,bs
button "." title "." at 3*sp+2*bs,th+3*bs+5*sp size bs,bs
button "Cl" title "C" at 4*sp+3*bs,th+2*sp size bs,bs
button "del" title "<-" at 4*sp+3*bs,th+bs+3*sp size bs,bs
button "ok" title "ok" at 4*sp+3*bs,th+2*bs+4*sp size bs,2*bs+sp
page name$ hide
page name$ frame xtop,ytop,ww,hh
set buttons default ! set buttons font size 20
draw font size 20 ! draw color 0,0,0
end def

def numpad(minval,maxval)
page "numpad" set ! page "numpad" show
a$="" ! pflag=0 ! sflag=0 ! ob=1-option_base()
nump1:
if b_p("ok") then
  number=val(a$) ! a$="" ! button "res" text ""
  if minval<>0 or maxval<>0 then
    if number<minval or number>maxval then
      button "res" text "range error"
      pflag=0 ! a$="" ! pause 1
      button "res" text ""
      goto nump1
      end if
    end if
  page "numpad" hide ! page "" set
  return number
  end if
if b_p("Cl") then
  a$ = "" ! pflag=0 ! sflag=0 ! goto nump3
  end if
if b_p("del") and len(a$) then
  ll=len(a$) ! if substr$(a$,ll-ob,ll-ob)="." then pflag=0
  a$ = left$(a$,ll-1) ! sflag=0 ! goto nump3
  end if
if b_p("-") then
  a$ = "-" ! pflag=0 ! sflag=0 ! goto nump3
  end if
if b_p(".") and not pflag and not sflag then
  a$ &= "." ! pflag=1 ! goto nump3
  end if
for k=0 to 9
  t$=k
  if b_p(t$) and not sflag then
    a$ &= t$ ! goto nump3
    end if
  next k
goto nump1
nump3:
if len(a$)>10 then ! sflag=1 ! goto nump1 ! end if
button "res" text a$
goto nump1
end def

' this is a m-degree polynomial least squares fit of a number of
' point in 2 dimensional space
' there are n points, with x- and y-coordinates in vectors x() and y()
' m is the degree of the polynomial (take 1 for straight line fit,
' m=2 for parabola, and so on)
' the coefficients of the best fit polynomial are returned in vector c()
' f.i. for m=2 : y = c(1) + c(2)*x + c(3)*x^2
'
def ls (n,m,x(),y(),c())
dim a1(20,20),a2(20,20),a3(20,20),rl(20)
m=m+1
for i=1 to n ! a1(i,1)=1
  for j=2 to m ! a1(i,j)=a1(i,j-1)*x(i) ! next j
  next i
mat_trans (n,m,a1,a2)
mat_mul (m,n,a2,a1,a3)
mat_vec(m,n,a2,y,rl)
mat_inv (m,a3,a1)
mat_vec(m,n,a1,rl,c)
end def

def vec_uit (n$,n,v(),x,y,l,d)
ys=y-30
if n$<>"" then
  draw text n$ at x,y
  ys=ys+30
end if
for i=1 to n
  a$=n2a$(i,3,0) & " " & n2a$(v(i),l,d)
  draw text a$ at x-12,ys+25*i
next i
end def

def mat_trans (n,m,mat(,),matt(,))
for i=1 to n
  for j=1 to m
    matt(j,i)=mat(i,j)
  next j
next i
end def

def mat_mul (n,m,mata(,),matb(,),matc(,))
for i=1 to n
  for j=1 to n
    tot=0
    for k=1 to m
      tot=tot+mata(i,k)*matb(k,j)
    next k
  matc(i,j)=tot
  next j
next i
end def

def mat_vec (n,m,mat(,),vin(),vuit())
for i=1 to n
  tot=0
  for j=1 to m
    tot=tot+mat( i,j)*vin(j)
  next j
  vuit(i)=tot
next i
end def

' produce the inverse of square matrix a() giving matrix ainv()
' returns 1 on success, 0 if matrix is singular
'
def mat_inv (nvar,a(,),ainv(,))
dim w(nvar,2*nvar)                      
for i=1 to nvar                 
  for j=1 to nvar ! w(i,j)=a(i,j) ! w(i,j+nvar)=0  ! next j
  w(i,i+nvar)=1
  next i
for piv=1 to nvar
  fac=w(piv,piv) ! if fac=0 then return 0
  for j=piv to piv+nvar ! w(piv,j)=w(piv,j)/fac ! next j
  for i=1 to nvar
    if i<>piv then
      fac=w(i,piv)
      for j=piv to piv+nvar ! w(i,j)=w(i,j)-fac*w(piv,j) ! next j
      endif
    next i
  next piv
for i=1 to nvar
  for j=1 to nvar ! ainv(i,j)=w(i,j+nvar) ! next j
  next i
return 1
end def

def b_p(a$) = button_pressed(a$)
IMG_1487.PNG
IMG_1487.PNG (936.15 KiB) Viewed 3334 times

User avatar
rbytes
Posts: 1338
Joined: Sun May 31, 2015 12:11 am
My devices: iPhone 11 Pro Max
iPad Pro 11
MacBook
Dell Inspiron laptop
CHUWI Plus 10 convertible Windows/Android tablet
Location: Calgary, Canada
Flag: Canada
Contact:

Re: Curve fitting with polynomials

Post by rbytes »

This works very well. Thanks for posting. :D
I am exceeding the chart range by quite a lot. When I change values, my old curves that extend beyond the chart don't get erased. Maybe it would be possible to rescale the chart so that the curves always stay inside it? Or just a screen clear when the values are changed..
Attachments
1F930ECA-94B6-4D98-AE42-A81C7C622DD5.png
1F930ECA-94B6-4D98-AE42-A81C7C622DD5.png (1005.84 KiB) Viewed 3325 times
The only thing that gets me down is gravity...

Henko
Posts: 814
Joined: Tue Apr 09, 2013 12:23 pm
My devices: iPhone,iPad
Windows
Location: Groningen, Netherlands
Flag: Netherlands

Re: Curve fitting with polynomials

Post by Henko »

Hi Richard,
I had already noticed that problem, and i will add a clipping facility to the curve drawing procedure.
Your example is extreme with a clear reason: the number of calculated coefficients equals the number of points. This means that the polynomial is forced to exactly cover the points. This leads to extreme values between the points.
Normally, the number of points ("measurements") is far, far greater than the degree of the fitting polynomial.
I already have added some little mods to the program, tomorrow i'll add the clipping thing and upload the prog.

User avatar
rbytes
Posts: 1338
Joined: Sun May 31, 2015 12:11 am
My devices: iPhone 11 Pro Max
iPad Pro 11
MacBook
Dell Inspiron laptop
CHUWI Plus 10 convertible Windows/Android tablet
Location: Calgary, Canada
Flag: Canada
Contact:

Re: Curve fitting with polynomials

Post by rbytes »

Thanks. I will start adding more points, because I see your point! :lol:
The only thing that gets me down is gravity...

Henko
Posts: 814
Joined: Tue Apr 09, 2013 12:23 pm
My devices: iPhone,iPad
Windows
Location: Groningen, Netherlands
Flag: Netherlands

Re: Curve fitting with polynomials

Post by Henko »

Final (?) version with the graph staying within the boundaries. Not by clipping, but by resizing (as suggested by mr. rbytes ;) )

Code: Select all

' curve fitting with polynomials in 2D
' by Henko, version 31-03-2018
'
np_max=200   ' maximum # of points, may be modified
m_max=9      ' maximum degree of fitting polynomial
f$="ls_data" ' file used to store and restore point data
option base 1
dim xy(np_max,2),xlist$(np_max),ylist$(np_max)
dim coeff(m_max+1),coef$(m_max+1)
'
prog_init()
do slowdown
  if b_p("new_x") then add_new_X()
  k=list_selected("xp") ! if k>0 and k<=.np then edit_X(k)
  k=list_selected("yp") ! if k>0 then edit_Y(k)
  if b_p("del_p") then del_point()
  if b_p("degree") then set_degree()
  if b_p("calc") then calc()
  if b_p("help") then help()
  if b_p("load") then load_data(f$)
  if b_p("save") then save_data(f$)
  if b_p("stop") then stop
until forever
end

def prog_init()
graphics ! graphics clear ! draw color 0,0,0
get screen size sw,sh ! cx=sw/2
.gx1=20 ! .gy1=500 ! .gx2=748 ! .gy2=908
.np=0 ! help_page()
c_list("xp","X-coordinates",.xlist$,.np,20,20,190,420)
c_list("yp","Y-coordinates",.ylist$,.np,210,20,380,420)
button "new_x" text "Add X-value" at 20,440 size 170,30
button "del_p" text "Delete point" at 210,440 size 170,30
field "points" text "    "&"0 Points" at 450,20 size 190,30 RO
field "points" font name "Courier-Bold"
button "calc" text "y = f(x)" at 450,440 size 290,30
button "degree" text "Degree" at 450,70 size 190,30
button "load" text "Load data" at 240,sh-50 size 120,30
button "save" text "Save data" at 400,sh-50 size 120,30
button "help" text "Help" at 20,sh-50 size 120,30
button "stop" text "Quit" at sw-140,sh-50 size 120,30
button "calc" hide
init_numpad(450,120,60,.7,.7,.7,1)
end def

def add_new_X()
if .np=.np_max then return
.np+=1 ! .xy(.np,2)=0 ! edit_X(.np)
field "points" text "    " & .np & " Points"
list "xp" select .np ! list "xp" select -1
list "yp" select .np ! list "yp" select -1
end def

def edit_X(k)
.xy(k,1)=numpad(0,0)
tablesort (.np,2,.xy,1)
disp_lists() ! plot()
end def

def edit_Y(k)
.xy(k,2)=numpad(0,0)
for i=1 to .np ! .ylist$(i)=i & ":  " &.xy(i,2) ! next i
c_list("yp","Y-coordinates",.ylist$,.np,210,20,380,420)
plot()
end def

def del_point()
k=numpad(1,.np)
if k<.np then
  for i=k to .np-1
    .xy(i,1)=.xy(i+1,1) ! .xy(i,2)=.xy(i+1,2)
    next i
  end if
.np-=1
field "points" text "    " & .np & " Points"
disp_lists() ! plot()
end def

def disp_lists()
for i=1 to .np
  .xlist$(i)=i & ":   " &.xy(i,1) ! .ylist$(i)=i & ":   " &.xy(i,2)
  next i
c_list("xp","X-coordinates",.xlist$,.np,20,20,190,420)
c_list("yp","Y-coordinates",.ylist$,.np,210,20,380,420)
end def

def stat()
.xmin=1E+8 ! .xmax=-1E+8 ! .ymin=1E+8 ! .ymax=-1E+8
for i=1 to .np
  x=.xy(i,1) ! .xmin=min(.xmin,x) ! .xmax=max(.xmax,x)
  y=.xy(i,2) ! .ymin=min(.ymin,y) ! .ymax=max(.ymax,y)
  next i
end def

def set_degree()
if .np<2 then return
.m=numpad(0,min(.np-1,.m_max))
button "degree" text "Degree = "&.m
button "calc" text "y = f(x)" ! button "calc" show
plot()
end def

def calc()
x=numpad(.xmin,.xmax)
y=poly(.m,.coeff,x)
button "calc" text "x="&x&"   =>   "&"y="&y
end def

def load_data(f$)
i=0
do
 i+=1 ! file f$ input .xy(i,1),.xy(i,2)
 until not data_exist(f$)
.np=i ! file f$ reset
field "points" text "    " & .np & " Points"
disp_lists() ! plot()
end def

def save_data(f$)
file f$ delete
for i=1 to .np ! file f$ print .xy(i,1),.xy(i,2) ! next i
end def

def plot()
dim xls(.np),yls(.np)
stat() ! if .xmin=.xmax or .ymin=.ymax or .np<2 then return
if .m>0 then
  for i=1 to .np ! xls(i)=.xy(i,1) ! yls(i)=.xy(i,2) ! next i
  ls (.np,.m,xls,yls,.coeff) ! c_table()
  dx=(.xmax-.xmin)/48
  for x=.xmin to .xmax step dx
    y=poly(.m,.coeff,x)
    .ymin=min(.ymin,y) ! .ymax=max(.ymax,y)
    next x
  end if
if .xmin>0 then
  sc_x=720/.xmax ! y_as=.gx1+4 ! goto lab1
  end if
if .xmax<0 then
  sc_x=720/-.xmin ! y_as=.gx2-4 ! goto lab1
  end if
sc_x=720/(.xmax-.xmin) ! y_as=.gx1+4-.xmin*sc_x
lab1:
if .ymin>0 then
  sc_y=400/.ymax ! x_as=.gy2-4 ! goto lab2
  end if
if .ymax<0 then
  sc_y=400/-.ymin ! x_as=.gy1+4 ! goto lab2
  end if
sc_y=400/(.ymax-.ymin) ! x_as=.gy2-4+.ymin*sc_y
lab2:
refresh off
fill color 0.88,0.85,0.65 ! fill rect .gx1,.gy1 to .gx2,.gy2
draw rect .gx1,.gy1 to .gx2,.gy2
draw color .7,.7,.7 ! draw size 1    ' draw grid
for x=.gx1+4 to .gx2-4 step 20
  draw line x,.gy1+4 to x,.gy2-4
  next x
for y=.gy1+4 to .gy2-4 step 20
  draw line .gx1+4,y to .gx2-4,y
  next y
draw color 0,0,1 ! draw size 3
draw line y_as,.gy1+4 to y_as,.gy2-4   ' draw axes
draw line .gx1+4,x_as to .gx2-4,x_as
draw color 0,0,0
for i=1 to .np
  x=xp(.xy(i,1)) ! y=yp(.xy(i,2))
  draw circle x,y size 2
  next i
if .m>0 then
  draw color 1,0,0 ! draw size 2
  y=yp(poly(.m,.coeff,.xmin)) ! y=max(y,.gy1+4) ! y=min(y,.gy2-4)
  draw to xp(.xmin),y
  for x=.xmin+dx to .xmax step dx
    y=yp(poly(.m,.coeff,x)) ! y=max(y,.gy1+4) ! y=min(y,.gy2-4)
    draw line to xp(x),y
    next x
  end if
refresh
draw color 0,0,0
end def

def xp(x)=plot.y_as+plot.sc_x*x
def yp(y)=plot.x_as-plot.sc_y*y

def c_table()
for i=1 to .m+1 ! .coef$(i)="a("&(i-1)&") = "&.coeff(i) ! next i
c_list("cf","Coefficients",.coef$,.m+1,450,120,640,420)
end def

def poly (n,coef(),x)
res=coef(n+1)
for i=n to 1 step-1
  res=res*x+coef(i)
next i
return res
end def

def n2a$(num,lang,dec)
b$="               "
fac=10^dec
num$=floor(fac*num)/fac
tot=lang-len(num$)
if tot<1 then tot=1
rem test ("", tot)
a$=substr$(b$,1,tot) & num$
n2a$=a$
end def

def help()
n$="help" ! page n$ set ! page n$ alpha 0 ! page n$ show
for i=0 to 1 step 0.01 ! page n$ alpha i ! pause 0.01 ! next i
do slowdown ! until b_p("h_stop")
for i=1 to 0 step -0.01 ! page n$ alpha i ! pause 0.01 ! next i
page n$ hide ! page "" set
end def

def help_page()
xs=40 ! ys=200 ! w=680 ! h=680
h$="Curve fitting with polynomials in 2 dimensions."&chr$(10)&"A set of x-y points can be fitted with a polynomial of a suitable degree."&chr$(10)&"In this program a maximum number of 200 point is set, but this can easily be changed in the dim function."&chr$(10)&"The maximum degree of the polynomial is set at 9, as a higher degrees seldom add to the quality of the fit."&chr$(10)&"x- and y-coordinates are kept in separate lists with their point number (for synchronizing purpose). Both lists are kept in sequence of the x-coordinates."&chr$(10)&"Points are added by touching the ""add x-value"" button. the corresponding y-value is defaulted to zero."&chr$(10)&"Existing x- and y-values may be edited bij touching them in the appropriate list."&chr$(10)&"Points may be deleted by touching the ""delete point"" button and entering the point number."&chr$(10)&"The collection of points can be saved to a file and loaded later on by touching the ""save"" and ""load"" buttons."&chr$(10)&"A graph will appear as soon as two points are entered with different x- and y-values. Each modification in the point data is right away reflected in the graph."&chr$(10)&"Also, when more than one point is present, a fit command may be issued by touching the ""degree"" button and entering the degree of the fitting polynomial."&chr$(10)&"The coefficients of the polynomial are then shown and the calculated curve is drawn in the graph."&chr$(10)&"An y-value may be calculated for a given x-value by using the ""y=f(x)"" button."
n$="help" ! page n$ set ! page n$ frame xs,ys,w,h
field "hh" text h$ at 0,0 size w,h ML RO
field "hh" font color 0,0,1 ! field "hh" font size 20
button "h_stop" text "Close" at w/2-30,h-40 size 60,25
page n$ hide
page "" set
end def

    ' {LS-utils}

' id$ = object name
' cont$ = array met elementen
' size = aantal elementen in de list
'
def c_list(id$,title$,cont$(),size,xt,yt,xb,yb)
size=max(1,size)
dim temp$(size)
for i=1 to size ! temp$(i)=cont$(i) ! next i
list id$ text temp$ at xt+2,yt+32 size xb-xt-4,yb-yt-34
draw size 3
draw rect xt,yt to xb,yb ! draw line xt,yt+30 to xb,yt+30
draw color 0,0,1 ! draw text title$ at xt+5,yt+5
draw color 0,0,0
end def

' numeric table sorting function (option base 1 assumed)
' r = number of rows
' c = number of columns
' in(,) = table to be sorted
' col = column to be used for sorting
'
def tablesort (r,c,in(,),col)
dim sortcol(r),index(r),out(r,c)
for i=1 to r ! sortcol(i)=in(i,col) ! next i
sort sortcol as index
for i=1 to r ! for j=1 to c
  out(i,j)=in(index(i),j)
  next j ! next i
for i=1 to r ! for j=1 to c ! in(i,j)=out(i,j) ! next j ! next i
return
end def

' numerical keypad object
' 
' produce a simple keypad to quickly enter a number in an app
' upon entry, the keypad disappears
' initialize once, multiple use after
' left upper corner is placed at "xtop,ytop"
' "bs" is the button size (keypad becomes 4.3 times larger)
' size of number is accepted between "minval" and "maxval"
' if both "minval" and "maxval" are zero, then no restrictions
' max number of tokens in the number is 10 (minus and dot included)
' works for option base 0 and 1
'
def init_numpad(xtop,ytop,bs,R,G,B,alpha)
name$="numpad" ! cn=10
page name$ set 
page name$ frame xtop,ytop,0,0
set buttons custom
if bs<20 then bs=20
sp=4 ! th=.5*bs+4 ! ww=4*bs+5*sp ! hh=th+4*bs+6*sp
fsize=.5*bs
draw font size fsize ! set buttons font size fsize
draw color 1,1,1 ! fill color .7,.7,.7
button "rec" title "" at 0,0 size ww,hh
button "res" title "" at 0,0 size ww,th+4
fill color R,G,B ! fill alpha alpha
button "0" title "0" at sp,th+3*bs+5*sp size bs,bs
for k=1 to 9
  x=(k-1)%3 ! y=2-floor((k-1)/3)
  button k title k at (x+1)*sp+x*bs,th+y*bs+(y+2)*sp size bs,bs
  next k
button "-" title "-" at 2*sp+bs,th+3*bs+5*sp size bs,bs
button "." title "." at 3*sp+2*bs,th+3*bs+5*sp size bs,bs
button "Cl" title "C" at 4*sp+3*bs,th+2*sp size bs,bs
button "del" title "<-" at 4*sp+3*bs,th+bs+3*sp size bs,bs
button "ok" title "ok" at 4*sp+3*bs,th+2*bs+4*sp size bs,2*bs+sp
page name$ hide
page name$ frame xtop,ytop,ww,hh
set buttons default ! set buttons font size 20
draw font size 20 ! draw color 0,0,0
end def

def numpad(minval,maxval)
page "numpad" set ! page "numpad" show
a$="" ! pflag=0 ! sflag=0 ! ob=1-option_base()
nump1:
if b_p("ok") then
  number=val(a$) ! a$="" ! button "res" text ""
  if minval<>0 or maxval<>0 then
    if number<minval or number>maxval then
      button "res" text "range error"
      pflag=0 ! a$="" ! pause 1
      button "res" text ""
      goto nump1
      end if
    end if
  page "numpad" hide ! page "" set
  return number
  end if
if b_p("Cl") then
  a$ = "" ! pflag=0 ! sflag=0 ! goto nump3
  end if
if b_p("del") and len(a$) then
  ll=len(a$) ! if substr$(a$,ll-ob,ll-ob)="." then pflag=0
  a$ = left$(a$,ll-1) ! sflag=0 ! goto nump3
  end if
if b_p("-") then
  a$ = "-" ! pflag=0 ! sflag=0 ! goto nump3
  end if
if b_p(".") and not pflag and not sflag then
  a$ &= "." ! pflag=1 ! goto nump3
  end if
for k=0 to 9
  t$=k
  if b_p(t$) and not sflag then
    a$ &= t$ ! goto nump3
    end if
  next k
goto nump1
nump3:
if len(a$)>10 then ! sflag=1 ! goto nump1 ! end if
button "res" text a$
goto nump1
end def

' this is a m-degree polynomial least squares fit of a number of
' point in 2 dimensional space
' there are n points, with x- and y-coordinates in vectors x() and y()
' m is the degree of the polynomial (take 1 for straight line fit,
' m=2 for parabola, and so on)
' the coefficients of the best fit polynomial are returned in vector c()
' f.i. for m=2 : y = c(1) + c(2)*x + c(3)*x^2
'
def ls (n,m,x(),y(),c())
dim a1(20,20),a2(20,20),a3(20,20),rl(20)
m=m+1
for i=1 to n ! a1(i,1)=1
  for j=2 to m ! a1(i,j)=a1(i,j-1)*x(i) ! next j
  next i
mat_trans (n,m,a1,a2)
mat_mul (m,n,a2,a1,a3)
mat_vec(m,n,a2,y,rl)
mat_inv (m,a3,a1)
mat_vec(m,n,a1,rl,c)
end def

def vec_uit (n$,n,v(),x,y,l,d)
ys=y-30
if n$<>"" then
  draw text n$ at x,y
  ys=ys+30
end if
for i=1 to n
  a$=n2a$(i,3,0) & " " & n2a$(v(i),l,d)
  draw text a$ at x-12,ys+25*i
next i
end def

def mat_trans (n,m,mat(,),matt(,))
for i=1 to n
  for j=1 to m
    matt(j,i)=mat(i,j)
  next j
next i
end def

def mat_mul (n,m,mata(,),matb(,),matc(,))
for i=1 to n
  for j=1 to n
    tot=0
    for k=1 to m
      tot=tot+mata(i,k)*matb(k,j)
    next k
  matc(i,j)=tot
  next j
next i
end def

def mat_vec (n,m,mat(,),vin(),vuit())
for i=1 to n
  tot=0
  for j=1 to m
    tot=tot+mat( i,j)*vin(j)
  next j
  vuit(i)=tot
next i
end def

' produce the inverse of square matrix a() giving matrix ainv()
' returns 1 on success, 0 if matrix is singular
'
def mat_inv (nvar,a(,),ainv(,))
dim w(nvar,2*nvar)                      
for i=1 to nvar                 
  for j=1 to nvar ! w(i,j)=a(i,j) ! w(i,j+nvar)=0  ! next j
  w(i,i+nvar)=1
  next i
for piv=1 to nvar
  fac=w(piv,piv) ! if fac=0 then return 0
  for j=piv to piv+nvar ! w(piv,j)=w(piv,j)/fac ! next j
  for i=1 to nvar
    if i<>piv then
      fac=w(i,piv)
      for j=piv to piv+nvar ! w(i,j)=w(i,j)-fac*w(piv,j) ! next j
      endif
    next i
  next piv
for i=1 to nvar
  for j=1 to nvar ! ainv(i,j)=w(i,j+nvar) ! next j
  next i
return 1
end def

def b_p(a$) = button_pressed(a$)

User avatar
rbytes
Posts: 1338
Joined: Sun May 31, 2015 12:11 am
My devices: iPhone 11 Pro Max
iPad Pro 11
MacBook
Dell Inspiron laptop
CHUWI Plus 10 convertible Windows/Android tablet
Location: Calgary, Canada
Flag: Canada
Contact:

Re: Curve fitting with polynomials

Post by rbytes »

Fascinating to watch the graph adjust itself to new coordinates. Well done!
The only thing that gets me down is gravity...

Post Reply