## 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 doubledim y(4) as doubledim w(4) as doubledim as double c0, c1, cov00, cov01, cov11, chisqdim as integer i, n=4, s=1data 1970, 1980, 1990, 2000  ' x datafor i=1 to 4   read x(i)nextdata 12, 11, 14, 13   ' y datafor i=1 to 4   read y(i)nextdata .1, .2, .3, .4   ' weight datafor i=1 to 4   read w(i)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' http://www.gnu.org/software/gsl/manual/html_node/Mixed_002dradix-FFT-routines-for-real-data.html' 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=100dim as double dataa(n)dim  real as gsl_fft_real_wavetable ptrdim  hc   as gsl_fft_halfcomplex_wavetable ptrdim  work as gsl_fft_real_workspace ptr' Fill the data array and print it out.for i=1 to n   dataa(i)=2nextfor i=n/3 to (2*n/3)   dataa(i)=1nextfor 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 arraygsl_fft_real_transform(@dataa(1), 1, n, real, work)' Check the resultsfor 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 arraygsl_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)nextgsl_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 iDim As Double dataa(5)Dim As Double median, upperq, lowerq, mean, variance, largest, smallestData 17.2, 18.1, 16.5, 18.3, 12.6For i=1 To 5   Read dataa(i)Nextmean = 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 ##.###";meanPrint Using "The estimated variance is ##.###";variancePrint Using "The largest value is ##.###";largestPrint Using "The smallest value is ##.###";smallestPrint Using "Original dataset: ##.#, ##.#, ##.#, ##.#, ##.#";dataa(1),dataa(2),dataa(3),dataa(4),dataa(5)' Now sort the data, so we can run more statisticsgsl_sort(@dataa(1),1,5)Print Using "Sorted dataset: ##.#, ##.#, ##.#, ##.#, ##.#";dataa(1),dataa(2),dataa(3),dataa(4),dataa(5)' These functions require sorted input datamedian=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 ##.###";medianPrint Using "The upper quartile is ##.###";upperqPrint Using "The lower quartile is ##.###";lowerqSleepend
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 sizeDim As Double dataa(n)Dim As Double median, upperq, lowerq, mean, variance, largest, smallestData 17.2, 18.1, 16.5, 18.3, 12.6For i=1 To n   Read dataa(i)Nextmean = 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 ##.###";meanPrint Using "The estimated variance is ##.###";variancePrint Using "The largest value is ##.###";largestPrint Using "The smallest value is ##.###";smallestPrint Using "Original dataset: ##.#, ##.#, ##.#, ##.#, ##.#";dataa(1),dataa(2),dataa(3),dataa(4),dataa(5)' Now sort the data, so we can run more statisticsgsl_sort(@dataa(1),1,n)Print Using "Sorted dataset: ##.#, ##.#, ##.#, ##.#, ##.#";dataa(1),dataa(2),dataa(3),dataa(4),dataa(5)' These functions require sorted input datamedian=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 ##.###";medianPrint Using "The upper quartile is ##.###";upperqPrint Using "The lower quartile is ##.###";lowerqSleepend
TJF
Posts: 3698
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)NEXTmean = 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 worspacedim 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)nextsleep' end
dodicat
Posts: 7194
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)
viewtopic.php?f=7&t=16763&p=147077&hilit=polynomial#p147077
Makoto WATANABE
Posts: 202
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: 202
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: 2888
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: 7194
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,imEnd TypeOperator +(Byval n1 As complex,Byval n2 As complex) As complexReturn Type<complex>(n1.re+n2.re,n1.im+n2.im)End OperatorOperator -(Byval n1 As complex,Byval n2 As complex) As complexReturn Type<complex>(n1.re-n2.re,n1.im-n2.im)End OperatorOperator *(Byval n1 As complex,Byval n2 As complex) As complexReturn Type<complex>(n1.re*n2.re - n1.im*n2.im,n1.im*n2.re + n1.re*n2.im)End OperatorOperator /(Byval n1 As complex,Byval n2 As complex) As complexDim As Double d = n2.re*n2.re+n2.im*n2.imReturn Type<complex>((n1.re*n2.re+n1.im*n2.im)/d,(n1.im*n2.re - n1.re*n2.im)/d)End OperatorOperator *(Byval n1 As Double,Byval n2 As complex) As complexReturn Type<complex>(n1*n2.re,n1*n2.im)End OperatorOperator =(Byval n1 As complex,Byval n2 As complex) As IntegerReturn n1.re=n2.re And n1.im=n2.imEnd 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 pEnd FunctionSub 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 countEnd SubDim Shared As Integer EFunction 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 x2End FunctionFunction 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 x2End FunctionSub 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    Dim As complex grad(Ubound(pol))    Redim As complex lastroot()    differentiate(pol(),grad())    Dim As complex grad2(Ubound(grad))    differentiate(grad(),grad2())    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)            If rootflag=1 Then   root=newton(pol(),grad(),grad2(),start,iterate)            If rootflag=2 Then root=newton1(pol(),grad(),start,iterate)            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 ";rootflagEnd SubSub 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,"_")    printEnd Sub Redim As complex answer()Dim As complex pol1(5)={ (-1),(0),(0),(0),(0),(1)} '-1 + x^5  get_roots(pol1(),,,answer())printroots(pol1(),answer())'complex quadratic polynomialdim as complex pol2(...)={ (-1,1),(8.8,-7),(3,-3)} get_roots(pol2(),,,answer())printroots(pol2(),answer())dim as complex pol3(...)={(0),(5),(0),(-20),(0),(16)} ' 5*x-20*x^3 +16*x^5 get_roots(pol3(),,,answer())printroots(pol3(),answer())Sleep 
srvaldez
Posts: 2888
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 functiondim w as gsl_integration_workspace ptr = gsl_integration_workspace_alloc(1000)dim result as doubledim errr as doubledim expected as double = -4.0dim alph as double = 1.0dim func as gsl_functionfunc.function = @ffunc.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          = ##.###############"; resultprint using "exact result    = ##.###############"; expectedprint using "estimated error = ##.###############"; errrprint using "actual error    = ##.###############"; result - expectedprint "intervals       = "; w->sizegsl_integration_workspace_free(w)sleep`
Makoto WATANABE
Posts: 202
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

### Re: GSL example usage

Dear srvaldez ;
I was able to find the roots of polynomial successfully.
Thank you for always teaching precisely.

Dear dodicat ;
Thank you for introducing a very interesting program.
I will study along with your "polynomial root finder / plotter".
Then I will translate your programs into Japanese and introduce them to Japanese people on my website.