Curve fitting with polynomials
Posted: Fri Mar 30, 2018 1:12 pm
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.
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 (936.15 KiB) Viewed 3336 times
)