## GSL example usage

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
keithpickering
Posts: 9
Joined: Feb 03, 2012 5:37

### GSL example usage

Thought I would post this example of using GSL within FB, including the proper use of pointers and arrays. This example is a pretty much direct translation of the C example for using the gsl_fit header. KP.

Code: Select all

' Test program for linear regression, using gsl_fit.bi functions.
'
#include once "gsl/gsl_cblas.bi"
#include "gsl/gsl_fit.bi"
dim x(4) as double
dim y(4) as double
dim w(4) as double
dim as double c0, c1, cov00, cov01, cov11, chisq
dim as integer i, n=4, s=1
data 1970, 1980, 1990, 2000  ' x data
for i=1 to 4
next
data 12, 11, 14, 13   ' y data
for i=1 to 4
next
data .1, .2, .3, .4   ' weight data
for i=1 to 4
next
'
' Standard linear fit:
' Note 1: We're passing pointers in many cases instead of actual
' variables. See gsl_fit.bi for which parms are pointers.
' Note 2: When passing a pointer to an array, we specify the first
' desired element of the array only. The function knows (from parm n)
' how many array elements there need to be in the regression.
'
gsl_fit_linear (@x(1),s,@y(1),s,n,@c0,@c1,@cov00,@cov01,@cov11,@chisq)
print using "best fit: y = ####.###### + ####.###### x"; c0, c1
print " covariance matrix:"
print using " ######.####  ######.#### "; cov00, cov01
print using " ######.####  ######.#### "; cov01, cov11
print using " chisq = ####.#### "; chisq
'
' Weighted linear fit, with data weights in w() array.
'
gsl_fit_wlinear (@x(1),s,@w(1),s,@y(1),s,n,@c0,@c1,@cov00,@cov01,@cov11,@chisq)
print
print "Weighted:"
print
print using "best fit: y = ####.###### + ####.###### x"; c0, c1
print " covariance matrix:"
print using " ######.####  ######.#### "; cov00, cov01
print using " ######.####  ######.#### "; cov01, cov11
print using " chisq = ####.#### "; chisq
' You should get a slope of .06 and intercept of -106.6 for both weighted
' and non-weighted versions.


keithpickering
Posts: 9
Joined: Feb 03, 2012 5:37

### GSL example: Fast Fourier Transform

Here's an example of how to use the FFT functions in the GSL package. Note that the actual transform functions in FB require you to pass a pointer to the data array rather than the array itself, which you do by pointing to the first array element. It's always best to eyeball the gsl_whatever.bi to check for datatypes you need to pass.

Code: Select all

' Test program for Fast Fourier Transform, using GSL functions.
' This is pretty much a direct translation of the C example found at
' See that page for more details.
'
#include once "gsl/gsl_cblas.bi"
#include "gsl/gsl_errno.bi"
#include "gsl/gsl_fft_real.bi"
#include "gsl/gsl_fft_halfcomplex.bi"

dim as integer i, n=100
dim as double dataa(n)

dim  real as gsl_fft_real_wavetable ptr
dim  hc   as gsl_fft_halfcomplex_wavetable ptr
dim  work as gsl_fft_real_workspace ptr

' Fill the data array and print it out.

for i=1 to n
dataa(i)=2
next

for i=n/3 to (2*n/3)
dataa(i)=1
next

for i=1 to n
print using "###:  ###.##############"; i, dataa(i)
next

' Allocate space.

work = gsl_fft_real_workspace_alloc(n)
real = gsl_fft_real_wavetable_alloc(n)

' Do the forward transform on the dataa array

gsl_fft_real_transform(@dataa(1), 1, n, real, work)

' Check the results

for i=1 to n
print using "###:  ###.##############"; i, dataa(i)
next

' Free wavetable, re-alloc for halfcomplex,
' and do the reverse transform on the resulting dataa array

gsl_fft_real_wavetable_free(real)
hc=gsl_fft_halfcomplex_wavetable_alloc(n)
gsl_fft_halfcomplex_inverse (@dataa(1), 1, n, hc, work)
gsl_fft_halfcomplex_wavetable_free(hc)

' Check the results again; should reconstruct original array within
' rounding error.

for i=1 to n
print using "###:  ###.##############"; i, dataa(i)
next

gsl_fft_real_workspace_free(work)

end

JohnK
Posts: 279
Joined: Sep 01, 2005 5:20
Location: Earth, usually
Contact:

### Re: GSL example: Fast Fourier Transform

Thanks for the GSL examples. Nice to know how to use this great library.
keithpickering
Posts: 9
Joined: Feb 03, 2012 5:37

### GSL example: basic statistics and sort

Continuing with GSL examples. In this one, we create an array, run some basic statistics on it, then sort the array and run a few other statistics that require a sorted input. Once again, note that we use a pointer to the dataa() array. In this example, I'm running in Windows instead of Linux, so I include the Sleep statement to keep the input visible.

Code: Select all

' FreeBasic example program for GSL basic statistical functions, with sort.
' This is pretty much a copy of the C example from the GSL page at
' http://www.gnu.org/software/gsl/manual/html_node/Example-statistical-programs.html
' See that page for more details.
'
#Include "gsl/gsl_blas.bi"
#Include "gsl/gsl_sort.bi"
#Include "gsl/gsl_statistics.bi"

Dim As Integer i
Dim As Double dataa(5)
Dim As Double median, upperq, lowerq, mean, variance, largest, smallest

Data 17.2, 18.1, 16.5, 18.3, 12.6
For i=1 To 5
Next

mean = gsl_stats_mean(@dataa(1), 1, 5)
variance = gsl_stats_variance(@dataa(1), 1, 5)
largest  = gsl_stats_max(@dataa(1), 1, 5)
smallest = gsl_stats_min(@dataa(1), 1, 5)

Print Using "The sample mean is ##.###";mean
Print Using "The estimated variance is ##.###";variance
Print Using "The largest value is ##.###";largest
Print Using "The smallest value is ##.###";smallest

Print Using "Original dataset: ##.#, ##.#, ##.#, ##.#, ##.#";dataa(1),dataa(2),dataa(3),dataa(4),dataa(5)

' Now sort the data, so we can run more statistics

gsl_sort(@dataa(1),1,5)

Print Using "Sorted dataset: ##.#, ##.#, ##.#, ##.#, ##.#";dataa(1),dataa(2),dataa(3),dataa(4),dataa(5)

' These functions require sorted input data

median=gsl_stats_median_from_sorted_data(@dataa(1),1,5)
upperq=gsl_stats_quantile_from_sorted_data(@dataa(1),1,5,0.75)
lowerq=gsl_stats_quantile_from_sorted_data(@dataa(1),1,5,0.25)

Print Using "The median is ##.###";median
Print Using "The upper quartile is ##.###";upperq
Print Using "The lower quartile is ##.###";lowerq

Sleep

end

keithpickering
Posts: 9
Joined: Feb 03, 2012 5:37

### Re: GSL example: basic statistics and sort

Perhaps even better is to use a variable to store the size of the array. You can then use the variable (instead of hard coding) to pass to the various functions. I should have thought of this before, but I was following the C example too slavishly. :)

Code: Select all

' FreeBasic example program for GSL basic statistical functions, with sort.
' This is pretty much a copy of the C example from the GSL page at
' http://www.gnu.org/software/gsl/manual/html_node/Example-statistical-programs.html
' See that page for more details.
'
#Include "gsl/gsl_blas.bi"
#Include "gsl/gsl_sort.bi"
#Include "gsl/gsl_statistics.bi"

Dim As Integer i, n=5       ' n is array size
Dim As Double dataa(n)
Dim As Double median, upperq, lowerq, mean, variance, largest, smallest

Data 17.2, 18.1, 16.5, 18.3, 12.6
For i=1 To n
Next

mean = gsl_stats_mean(@dataa(1), 1, n)
variance = gsl_stats_variance(@dataa(1), 1, n)
largest  = gsl_stats_max(@dataa(1), 1, n)
smallest = gsl_stats_min(@dataa(1), 1, n)

Print Using "The sample mean is ##.###";mean
Print Using "The estimated variance is ##.###";variance
Print Using "The largest value is ##.###";largest
Print Using "The smallest value is ##.###";smallest

Print Using "Original dataset: ##.#, ##.#, ##.#, ##.#, ##.#";dataa(1),dataa(2),dataa(3),dataa(4),dataa(5)

' Now sort the data, so we can run more statistics

gsl_sort(@dataa(1),1,n)

Print Using "Sorted dataset: ##.#, ##.#, ##.#, ##.#, ##.#";dataa(1),dataa(2),dataa(3),dataa(4),dataa(5)

' These functions require sorted input data

median=gsl_stats_median_from_sorted_data(@dataa(1),1,n)
upperq=gsl_stats_quantile_from_sorted_data(@dataa(1),1,n,0.75)
lowerq=gsl_stats_quantile_from_sorted_data(@dataa(1),1,n,0.25)

Print Using "The median is ##.###";median
Print Using "The upper quartile is ##.###";upperq
Print Using "The lower quartile is ##.###";lowerq

Sleep

end

TJF
Posts: 3758
Joined: Dec 06, 2009 22:27
Location: N47°, E15°
Contact:

### Re: GSL example: basic statistics and sort

The C source
• double dataa[5];
creates an array of 5 elements indexed 0 to 4. Its different in FB:
• DIM AS DOUBLE dataa(5)
creates an array of 6 elements indexed 0 to 5. Your example doesn't use the first element in this array. It's not a bug, just a waste of memory.

Instead of the CONST n you may use the UBOUND function:

Code: Select all

DIM AS DOUBLE dataa(...) = {17.2, 18.1, 16.5, 18.3, 12.6}
FOR i AS INTEGER = 0 TO UBOUND(dataa)
? dataa(i)
NEXT

mean = gsl_stats_mean(@dataa(0), 0, UBOUND(dataa))
...

BTW: nice examples!
keithpickering
Posts: 9
Joined: Feb 03, 2012 5:37

### Re: GSL example: basic statistics and sort

Thanks for the tip, we can all learn something new. :)
keithpickering
Posts: 9
Joined: Feb 03, 2012 5:37

### GSL example: Roots of polynomials

Here's an example program that uses the GSL library to find roots of polynomials of arbitrary power. The resultant roots are complex numbers, so if you're looking for the real roots, you need to look for those roots that have an imaginary part of zero.

Code: Select all

' Example program to find roots of polynomials.
' This is basicly a translation from the C program at
' http://www.gnu.org/software/gsl/manual/html_node/Roots-of-Polynomials-Examples.html
' See that page for more details.

' To demonstrate the use of the general polynomial solver we will take
' the polynomial P(x) = x^5 - 1 which has the following roots:
'     1, e^{2\pi i /5}, e^{4\pi i /5}, e^{6\pi i /5}, e^{8\pi i /5}

#include "gsl/gsl_cblas.bi"
#include "gsl/gsl_poly.bi"

dim as integer i, n=6  ' n is the highest power of x, plus 1.

'  coefficients of P(x) =  -1 + x^5
'  so all lesser powers of x have coefficients of zero. Put coefficients
'  into an array starting with the zeroth power and going up.
dim as double a(6)={-1,0,0,0,0,1}
'  array z will hold the roots as complex numbers.
'  Size accordingly: n-1 roots, with 2 doubles each.
dim as double z(2*(n-1))
'  allocate worspace
dim w as gsl_poly_complex_workspace ptr=gsl_poly_complex_workspace_alloc(6)

gsl_poly_complex_solve(@a(0),n,w,@z(0))

gsl_poly_complex_workspace_free(w)

print "       Real part           Imaginary part"
for i=0 to 4
print using "z# = ##.###############  ##.###############";i,z(2*i),z(2*i+1)
next
sleep
'
end


dodicat
Posts: 7395
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: GSL example: Roots of polynomials

Hi.
I got the two dll's by google, and it runs fine.
Thank you.
Here's my own polynomial root finder/plotter, gives the same answer for x^5-1.
(WINDOWS)
http://www.freebasic.net/forum/viewtopi ... al#p147077
Makoto WATANABE
Posts: 215
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

### Re: GSL example usage

Dear keithpickering ;

Thanks for posting these fine examples.
I would like to translate your examples into Japanese and introduce them to Japanese people on my website.
Makoto WATANABE
Posts: 215
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

### Re: GSL example usage

In the last example "Example program to find roots of polynomials", the following errors are displayed.

FindRootsOfPolynomials.o:fake:(.text+0x144): undefined reference to gsl_poly_complex_workspace_alloc'
FindRootsOfPolynomials.o:fake:(.text+0x168): undefined reference to gsl_poly_complex_solve'
FindRootsOfPolynomials.o:fake:(.text+0x176): undefined reference to gsl_poly_complex_workspace_free'

Please let me know if there is a correction method.
srvaldez
Posts: 2973
Joined: Sep 25, 2005 21:54

### Re: GSL example usage

hello Makoto WATANABE
add #inclib "gsl" after the #include's, for example

Code: Select all

#include "gsl/gsl_cblas.bi"
#include "gsl/gsl_poly.bi"
#inclib "gsl"

Last edited by srvaldez on Feb 08, 2017 14:03, edited 1 time in total.
dodicat
Posts: 7395
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: GSL example usage

Here is a rough and ready method to get the roots of a polynomial (real or complex coefficients).
You can tweak two parameters if necessary.
It is not 100% guaranteed for every polynomial, but then again neither is the eigenvalue method, which I assume gsl uses.

Code: Select all



Type complex
As Double re,im
End Type

Operator +(Byval n1 As complex,Byval n2 As complex) As complex
Return Type<complex>(n1.re+n2.re,n1.im+n2.im)
End Operator

Operator -(Byval n1 As complex,Byval n2 As complex) As complex
Return Type<complex>(n1.re-n2.re,n1.im-n2.im)
End Operator

Operator *(Byval n1 As complex,Byval n2 As complex) As complex
Return Type<complex>(n1.re*n2.re - n1.im*n2.im,n1.im*n2.re + n1.re*n2.im)
End Operator

Operator /(Byval n1 As complex,Byval n2 As complex) As complex
Dim As Double d = n2.re*n2.re+n2.im*n2.im
Return Type<complex>((n1.re*n2.re+n1.im*n2.im)/d,(n1.im*n2.re - n1.re*n2.im)/d)
End Operator

Operator *(Byval n1 As Double,Byval n2 As complex) As complex
Return Type<complex>(n1*n2.re,n1*n2.im)
End Operator

Operator =(Byval n1 As complex,Byval n2 As complex) As Integer
Return n1.re=n2.re And n1.im=n2.im
End Operator

#define Cabs(z)  sqr(z.re*z.re+z.im*z.im)
#define conjugate(z) type<complex>(z.re,-z.im)

Function poly(coff() As complex,Byval z As complex) As complex
Dim As Integer deg=Ubound(coff)
Dim As complex p=coff(deg)
For i As Integer=deg-1 To 0 Step -1
p=p*z+coff(i)
Next i
Return p
End Function

Sub differentiate(pol() As complex,Dpol() As complex)
Redim Dpol(Ubound(pol)-1)
For count As Integer=0 To Ubound(Dpol)
Dpol(count)=(count+1)*pol(count+1)
Next count
End Sub

Dim Shared As Integer E
Function Newton1(f0() As complex,f1() As complex,Byval start As complex,Byval iterations As Integer) As complex
Dim As complex x1=start,x2
E=iterations
For z As Integer=0 To iterations
x2=x1-poly(f0(),x1)/poly(f1(),x1)
If x1=x2 Then E=z: Exit For
x1=x2
Next z
Return x2
End Function

Function Newton(f0() As complex,f1() As complex,f2() As complex,Byval start As complex,Byval iterations As Integer) As complex
Dim As complex x1=start,x2
E=iterations
Dim As complex numerator,denominator,d0,d1,d2

For z As Integer=0 To iterations
d0=poly(f0(),x1)
d1=poly(f1(),x1)
d2=poly(f2(),x1)
numerator=d0*d1
denominator=d1*d1-.5*d0*d2
x2=x1-numerator/denominator
If x1=x2 Then E=z: Exit For
x1=x2
Next z
Return x2
End Function

Sub get_roots(pol() As complex,Byval iterate As Integer=1,Byval range As Single=1,roots() As complex)
#macro sort(a)
Scope
Dim As Integer n=Ubound((a))
For p1 As Integer  = 1 To n - 1
For p2 As Integer  = p1 + 1 To n
If a(p1).re>=a(p2).re Then Swap a(p1),a(p2)
Next p2
Next p1
End Scope
#endmacro

Dim As Integer real=1
Dim As complex cj
For z As Integer=0 To Ubound(pol)
If pol(z).im<>0 Then real=0:Exit For
Next z
Dim As Single mag
For z As Integer=0 To Ubound(pol)
mag=mag+Abs(pol(z).re)
Next z
mag=mag/Ubound(pol)
Dim As complex start,root,er,centre
Redim As complex lastroot()
Dim As Integer count,flag,count2,deg=Ubound(pol)/10+1,ecount,rootflag=1
Dim As Double stepsize=1,small=1/(10^deg)
small=.00000001
start:
Randomize
'do
For x As Double=centre.re-range To centre.re+range Step stepsize
For y As Double=centre.im-range To centre.im+range Step stepsize
start=Type<complex>(x,y)
er=poly(pol(),root)
If cabs(er)> small  Then Goto jump '1e-6
If Abs(root.im)<1e-12 Then root.im=0
flag=1
If e=iterate Then
ecount=ecount+1
centre=Type<complex>(10*Rnd*1-10*Rnd*1,10*Rnd*1-10*Rnd*1)
End If
count=count+1
For z2 As Integer=1 To Ubound(lastroot)
If cabs((root-lastroot(z2)))<1e-10 Then flag=0:Exit For
Next z2
If flag =1 Then
count2=count2+1
Redim Preserve roots(count2)
roots(count2)=root
Redim Preserve lastroot(count2)
lastroot(count2)=root

If Ubound(roots)= Ubound(pol) Then Goto oot
centre=Type<complex>(root.re,root.im)

End If
nxt:
jump:
Next y
Next x
If range>29 Then  'arbitary get out
Beep
Print "ERROR -- Range "; range;"  Iteratations ";iterate
sort(roots)
Exit Sub
End If
If range>23 Then rootflag=2
If range>20 Then

centre=Type<complex>(10*Rnd*1-10*Rnd*1,10*Rnd*1-10*Rnd*10)
End If
If range>10 Then

If Ubound(roots) Then centre=Type<complex>(roots(Ubound(roots)-1).re,roots(Ubound(roots)-1).im)
End If

If range>5 Then
stepsize=.9
End If

If Ubound(roots)< Ubound(pol) Then
range=range+1:iterate=iterate+10:Goto start
End If

oot:
sort(roots)
' Print "Range ";range;"  Iteratations ";iterate;" r ";rootflag
End Sub

Sub printroots(pol() As complex,roots() As complex)
Dim As complex er
Dim As Double er2
Print "ROOTS";
Print Tab(5);"Real";Tab(30);"Imag";Tab(57);"Error"
For z As Double=1 To Ubound(roots)
er=poly(pol(),roots(z)):er2=cabs(er)
Print Tab(0);z & ") ";Tab(5);roots(z).re;Tab(30);roots(z).im & "i";Tab(55); Using "###.###";er2;
Next z
print
Print string(60,"_")
print
End Sub

Dim As complex pol1(5)={ (-1),(0),(0),(0),(0),(1)} '-1 + x^5

dim as complex pol2(...)={ (-1,1),(8.8,-7),(3,-3)}

dim as complex pol3(...)={(0),(5),(0),(-20),(0),(16)} ' 5*x-20*x^3 +16*x^5

Sleep


srvaldez
Posts: 2973
Joined: Sep 25, 2005 21:54

### Re: GSL example usage

here's the integration example from the gsl manual

Code: Select all

#include once "gsl/gsl_integration.bi"
#inclib "gsl"

function f cdecl (byval x as double, byval params as any ptr) as double
dim alph as double = *cptr(double ptr, params)
f = log(alph * x) / sqrt(x)
end function

dim w as gsl_integration_workspace ptr = gsl_integration_workspace_alloc(1000)
dim result as double
dim errr as double
dim expected as double = -4.0
dim alph as double = 1.0
dim func as gsl_function

func.function = @f
func.params = @alph
'gsl_integration_qags (const gsl_function * f, double a, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)
gsl_integration_qags(@func, 0, 1, 0, 1e-7, 1000, w, @result, @errr)
print using "result          = ##.###############"; result
print using "exact result    = ##.###############"; expected
print using "estimated error = ##.###############"; errr
print using "actual error    = ##.###############"; result - expected
print "intervals       = "; w->size
gsl_integration_workspace_free(w)

sleep
`
Makoto WATANABE
Posts: 215
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

Dear srvaldez ;