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

Post by keithpickering »

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
	read x(i)
next
data 12, 11, 14, 13   ' y data
for i=1 to 4
	read y(i)
next
data .1, .2, .3, .4   ' weight data
for 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

Post by keithpickering »

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

Post by JohnK »

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

Post by keithpickering »

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
	Read dataa(i)
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

Post by keithpickering »

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
   Read dataa(i)
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

Post by TJF »

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

Post by keithpickering »

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

GSL example: Roots of polynomials

Post by keithpickering »

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: 7396
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: GSL example: Roots of polynomials

Post by dodicat »

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

Post by Makoto WATANABE »

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.
Please consent to this.
Makoto WATANABE
Posts: 215
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

Re: GSL example usage

Post by Makoto WATANABE »

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

Post by srvaldez »

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: 7396
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: GSL example usage

Post by dodicat »

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
    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 ";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 



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 polynomial
dim 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: 2973
Joined: Sep 25, 2005 21:54

Re: GSL example usage

Post by srvaldez »

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:

Re: GSL example usage

Post by Makoto WATANABE »

Dear srvaldez ;
Thanks for your quick post.
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.
Post Reply