FreeBASIC GSL ライブラリの使用例

' Test program for linear regression, using gsl_fit.bi functions.
' gsl_fit.bi 関数を使った線形回帰のテスト・プログラム。

'GSL example usage
'https://www.freebasic.net/forum/viewtopic.php?f=7&t=19374
'by keithpickering ≫ Feb 03, 2012 5:42

#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: 多くの場合、実際の変数ではなく、ポインタを渡します。
' どの parms がポインタに対応するかは gsl_fit.bi を参照してください。
' Note 2: 配列にポインタを渡すときは、配列の最初の要素のみを指定します。 
' この関数は（parm n から）回帰に必要な配列要素の数を知っています。
'
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 " 共分散行列:"
   print using " ######.####  ######.#### "; cov00; cov01
   print using " ######.####  ######.#### "; cov01; cov11
   Print using " chisq = ####.#### "; chisq
'   
' ウェイト付き線形近似。データの重みは w() 配列です。
'
gsl_fit_wlinear (@x(1),s,@w(1),s,@y(1),s,n,@c0,@c1,@cov00,@cov01,@cov11,@chisq)
   Print
   Print "加重版:"
   Print
   Print using "best fit: y = ####.###### + ####.###### x"; c0; c1
   Print " 共分散行列:"
   print using " ######.####  ######.#### "; cov00; cov01
   print using " ######.####  ######.#### "; cov01; cov11
   print using " chisq = ####.#### "; chisq

' 加重版と非加重版の両方で、
' 傾きは .06 で、切片(intercept)は -106.6 になるはずです。

Sleep

' GSL 関数を使った高速フーリエ変換のテストプログラム。
' これは、下のサイトで見つけた C からの移植です
' http://www.gnu.org/software/gsl/manual/html_node/Mixed_002dradix-FFT-routines-for-real-data.html
' 詳細は上のページを参照下さい。

'Fast Fourier Transform
'https://www.freebasic.net/forum/viewtopic.php?f=7&t=19374#p170746
'by keithpickering ≫ Feb 12, 2012 7:18

#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

' データ配列を満たして、表示します。

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   

Print "データ配列を満たして、表示しました。"
Print "何かキーを押すと、データ配列の前向き変換を行って表示します。"
Sleep

' スペースを割り当てます。

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

' データ配列の前向き変換を行う

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

' 結果を確認する

For i=1 To n
   Print Using "###:  ###.##############"; i; dataa(i)
Next   

Print "データ配列を前向き変換した結果を表示しました。"
Print "何かキーを押すと、結果データ配列の逆変換を行って表示します。"
Sleep

' ウェーブテーブルを解放、ハーフコンプレックスの再割り当て、
' 結果データ配列の逆変換

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)

' 結果をもう一度確認してください。 
' 丸め誤差内で元の配列を再構築する必要があります。

For i=1 To n
   Print Using "###:  ###.##############"; i; dataa(i)
Next

gsl_fft_real_workspace_free(work)

Print "結果データ配列の逆変換を表示しました。"
Print "何かキーを押すと、終了します。"
Sleep
End

' GSL 基本統計関数のための FreeBasic サンプルプログラム、ソート付き.
' これは GSL ページの C の例のコピーです。
' http://www.gnu.org/software/gsl/manual/html_node/Example-statistical-programs.html
' 詳細は上のページを参照してください.
'
'basic statistics and sort
'https://www.freebasic.net/forum/viewtopic.php?f=7&t=19374#p170746
'by keithpickering ≫ Feb 13, 2012 15:57

#Include "gsl/gsl_blas.bi"
#Include "gsl/gsl_sort.bi"
#Include "gsl/gsl_statistics.bi"

Dim As Double dataa(...) = {17.2, 18.1, 16.5, 18.3, 12.6}
Dim As Double median, upperq, lowerq, mean, variance, largest, smallest
Dim As Integer count

count = UBound(dataa)+1  '-lang fb では、Option Base は使えない。配列の下限は 0 。 
Print "データ数：", count
Print

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

Print Using "サンプルの平均値は ##.###";mean
Print Using "推定分散は ##.###";variance
Print Using "最大値は ##.###";largest
Print Using "最小値は ##.###";smallest
Print

Print "元のデータセット　　　: ";
For i As Integer = 0 To UBound(dataa)
   Print Using " ##.# ";dataa(i);
Next
Print

' データを並べ替えます。もっと多くの統計計算ができます

gsl_sort(@dataa(0),1,count)

Print "並べ替えたデータセット: ";
For i As Integer = 0 To UBound(dataa)
   Print Using " ##.# ";dataa(i);
Next
Print
Print

' 下の関数は、ソートされた入力データが必要です

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

Print Using "中央値は ##.###";median
Print Using "上の四分位数は ##.###";upperq
Print Using "下の四分位数は ##.###";lowerq

Sleep

End

' 多項式の根を見つけるプログラム例
' これは基本的に C プログラムの移植です。
' http://www.gnu.org/software/gsl/manual/html_node/Roots-of-Polynomials-Examples.html
' 詳細は上のページを参照してください.

'Roots of polynomials
'https://www.freebasic.net/forum/viewtopic.php?f=7&t=19374#p170746
'by keithpickering ≫ Feb 15, 2012 18:48


' 一般多項式ソルバの使用法を示すために、
' 多項式 P(x) = x^5 - 1 を使います。この多項式は次の根を持ちます：

'     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"
#inclib "gsl"

Dim As Integer i, n=6  ' n は x の最大累乗に 1 を加えたものです。

'  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.
'  x のすべてのより小さい累乗はゼロの係数を持つ。 
'  係数を0次から始まり増加する配列に入れる。

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.
'  配列zは、根を複素数として保持します。
'  サイズはそれに応じて：n-1 根、それぞれ 2 倍。

Dim As Double z(2*(n-1))
'  ワークスペースを割り当てる
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 "         実数部                虚数部"
For i=0 To 4
   Print Using "z# = ##.##################  ##.##################";i;z(2*i);z(2*i+1)
Next
Sleep
' 
End

'これは gsl マニュアルの積分の例です
'http://www.gnu.org/software/gsl/manual/html_node/Numerical-integration-examples.html


'https://www.freebasic.net/forum/viewtopic.php?t=19374
'by srvaldez ≫ Feb 08, 2017 14:08

#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