gsl-shell.git - gsl-shell

diff --git a/agg-plot/plot.h b/agg-plot/plot.h
index 96a35570..e2b6467b 100644
--- a/agg-plot/plot.h
+++ b/agg-plot/plot.h

@@ -83,7 +83,7 @@ public:

const char *get_title() const { return m_title; };

bool use_units() const { return m_use_units; };

- void set_units(bool use_units) { m_use_units = use_units; };

+ void set_units(bool use_units);

void add(VertexSource* vs, agg::rgba8 color, bool outline = false)

{

@@ -366,4 +366,11 @@ void plot<VS,RM>::viewport_scale(agg::trans_affine& m)

trans_affine_compose (m, rsz);

}

+template<class VS, class RM>

+void plot<VS,RM>::set_units(bool use_units)

+{

+ m_use_units = use_units;

+ this->update_viewport_trans();

+}

+

#endif

diff --git a/agg-plot/units.cpp b/agg-plot/units.cpp
new file mode 100644
index 00000000..d2db402e
--- /dev/null
+++ b/agg-plot/units.cpp

@@ -0,0 +1,68 @@

+

+#include <stdio.h>

+#include <math.h>

+

+#include "utils.h"

+#include "units.h"

+

+void units::init(double yinf, double ysup, double spacefact)

+{

+ double del;

+

+ if (ysup == yinf)

+ ysup = yinf + 1.0;

+

+ del = (ysup - yinf) / spacefact;

+

+ order = (int) floor(log10(del));

+

+ double expf = pow(10, order);

+ double delr = del / expf;

+

+ if (5 <= delr)

+ m_major = 5;

+ else if (2 <= delr)

+ m_major = 2;

+ else

+ m_major = 1;

+

+ m_inf = (int) floor(yinf / (m_major * expf) + 1e-5);

+ m_sup = (int) ceil (ysup / (m_major * expf) - 1e-5);

+

+ nb_decimals = (order < 0 ? -order : 0);

+

+ dmajor = m_major * expf;

+}

+

+void units::mark_label (char *lab, unsigned size, int mark) const

+{

+ bool minus = (m_inf < 0);

+ int asup = (minus ? -m_inf : m_sup);

+ char fmt[16];

+

+ if (size < 16)

+ return;

+

+ if (nb_decimals == 0)

+ {

+ snprintf (lab, size, "%d", int(mark * dmajor));

+ lab[size-1] = '0円';

+ }

+ else

+ {

+ int dec = (nb_decimals < 10 ? nb_decimals : 9);

+ int base = (int) floor(asup * dmajor);

+ int space = dec + (base > 0 ? (int)log10(base): 0) + 1 \

+ + (minus ? 1 : 0) + 1;

+ snprintf (fmt, 16, "%%%i.%if", space, dec);

+ fmt[15] = '0円';

+ snprintf (lab, size, fmt, mark * dmajor);

+ lab[size-1] = '0円';

+ }

+}

+

+double units::mark_scale (double x)

+{

+ double xinf = m_inf * dmajor, xsup = m_sup * dmajor;

+ return (x - xinf) / (xsup - xinf);

+}

diff --git a/base.lua b/base.lua
new file mode 100644
index 00000000..45be36ee
--- /dev/null
+++ b/base.lua

@@ -0,0 +1,68 @@

+

+local cat = table.concat

+local insert = table.insert

+

+function divmod(n, p)

+ local r = n % p

+ return (n-r)/p, r

+end

+

+local function tos(t, maxdepth)

+ if type(t) == 'table' then

+ if maxdepth <= 0 then return '<table>' end

+

+ local ils = {}

+ for i, v in ipairs(t) do ils[i] = v end

+

+ local ls = {}

+ for i, v in ipairs(ils) do insert(ls, tos(v, maxdepth-1)) end

+ for k, v in pairs(t) do

+ if not ils[k] then insert(ls, k .. '= ' .. tos(v, maxdepth-1)) end

+ end

+

+ return '{' .. cat(ls, ', ') .. '}'

+ elseif type(t) == 'function' then

+ return '<function>'

+ elseif type(t) == 'userdata' then

+ local ftostr = getmetatable(t).__tostring

+ if ftostr then return ftostr(t) else

+ return string.format('<%s>', gsltype(t))

+ end

+ else

+ return tostring(t)

+ end

+end

+

+local function myprint(...)

+ for i, v in ipairs(arg) do

+ if i > 1 then io.write(', ') end

+ io.write(tos(v, 3))

+ end

+ io.write('\n')

+end

+

+print = myprint

+

+function ilist(f, a, b)

+ a, b = (b and a or 1), (b and b or a)

+ if not b or type(b) ~= 'number' then

+ error 'argument #2 should be an integer number'

+ end

+ local ls = {}

+ for k= a, b do ls[#ls+1] = f(k) end

+ return ls

+end

+

+function sequence(f, a, b)

+ a, b = (b and a or 1), (b and b or a)

+ if not b or type(b) ~= 'number' then

+ error 'argument #2 should be an integer number'

+ end

+ local k = a

+ return function()

+ if k <= b then

+ k = k+1

+ return f(k-1)

+ end

+ end

+end

diff --git a/doc/source/contour-plot-hyper.png b/doc/source/contour-plot-hyper.png
new file mode 100644
index 00000000..9bfabef9
--- /dev/null
+++ b/doc/source/contour-plot-hyper.png

diff --git a/doc/source/contour.rst b/doc/source/contour.rst
new file mode 100644
index 00000000..df2b8bd7
--- /dev/null
+++ b/doc/source/contour.rst

@@ -0,0 +1,23 @@

+.. highlight:: lua

+

+.. include:: <isogrk1.txt>

+

+Contour Plots

+=============

+

+Overview

+--------

+

+GSL shell offer a contour plot function to draw contour curve of bidimensional functions. The current algorthm works correctly only for continous functions and it may gives bad results if the function have discontinuities.

+

+Here an example of its utilisation to plot the function :math:`f(x,y) = x^2 - y^2`::

+

+ require 'contour'

+

+ contour(|x,y| x^2 - y^2, {-8, -8}, {8, 8})

+

+.. figure:: contour-plot-hyper.png

+

+.. function:: contour(f, {xmin, ymin}, {xmax, ymax}[, ngridx, ngridy, contours])

+

+ Plot a contour plot of the function ``f``.

diff --git a/doc/source/fft-example-power-spectrum.png b/doc/source/fft-example-power-spectrum.png
index 2a2255a8..adcc25f7 100644
--- a/doc/source/fft-example-power-spectrum.png
+++ b/doc/source/fft-example-power-spectrum.png

diff --git a/doc/source/fft.rst b/doc/source/fft.rst
index e33e2467..97082019 100644
--- a/doc/source/fft.rst
+++ b/doc/source/fft.rst

@@ -1,205 +1,206 @@

-

-.. include:: <isogrk1.txt>

-

-Fast Fourier Transform

-==============================

-

-Mathematical Definitions

-------------------------

-

-Fast Fourier Transforms are efficient algorithms for calculating the discrete fourier transform (DFT)

-

-.. math::

- x_j = \sum_{k=0}^{n-1} z_k \exp(-2\pi i j k / n)

-

-The DFT usually arises as an approximation to the continuous fourier

-transform when functions are sampled at discrete intervals in space or

-time. The naive evaluation of the discrete fourier transform is a matrix-vector multiplication :math:`W\vec{z}`. A general matrix-vector multiplication takes O(n\ :sup:`2`) operations for n data-points. Fast fourier transform algorithms use a divide-and-conquer strategy to factorize the matrix W into smaller sub-matrices, corresponding to the integer factors of the length n. If n can be factorized into a product of integers f\ :sub:`1` f\ :sub:`2` ... f\ :sub:`m` then the DFT can be computed in O(n |Sgr| f\ :sub:`i`) operations. For a radix-2 FFT this gives an operation count of O(n log\ :sub:`2` n).

-

-All the FFT functions offer two types of transform: forwards and inverse, based on the same mathematical definitions. The definition of the forward fourier transform, ``fft(z)``, is,

-

-.. math::

- x_j = \sum_{k=0}^{n-1} z_k \exp(-2\pi i j k / n)

-

-and the definition of the inverse fourier transform, ``fft_inv(x)``, is,

-

-.. math::

- z_j = {1 \over n} \sum_{k=0}^{n-1} x_k \exp(2\pi i j k / n).

-

-The factor of 1/n makes this a true inverse.

-

-In general there are two possible choices for the sign of the

-exponential in the transform/ inverse-transform pair. GSL follows the

-same convention as fftpack, using a negative exponential for the

-forward transform. The advantage of this convention is that the

-inverse transform recreates the original function with simple fourier

-synthesis. Numerical Recipes uses the opposite convention, a positive

-exponential in the forward transform.

-

-GSL Shell interface

--------------------

-

-When you perform a Fourier transform different routines can be used depending on the size N of the sample. GSL features two kinds of routines, radix 2 and mixed radix. The radix 2 routine can be used if the size of the sample is a power of 2 while mixed radix can be used for any size of the input data but will be "fast" only if the size can be factorized into the product of small prime numbers.

-

-GSL Shell will select automatically the appropriate routine for you in a transparent way. But you should nevertheless be aware that the calculation will be fast only if the size of data can be factorised in small primes. You should also know that the mixed radix algorithm requires additional memory spaces so if you need to minimise memory usage you should feed only data of power of two size.

-

-GSL shell manage automatically the additional memory of precomputed trigonometric tables that may be required by the routine. The strategy used will be very efficient if you will perform many times the Fourier Transform of data always of the *same* size. If you perform many Fourier transform always with different size you will incur in a performance penality. Please not anyway that *no additional workspaces* are needed with radix two transformations.

-

-Fourier Transform of Real Data

-------------------------------

-

-In principle the Fourier Transform of real data is exactly the same of the Fourier transform of complex data, we only make the difference between them because:

-

-- in GSL Shell real value matrix and complex value matrix are considered

- different types

-- specialised routines exists to work efficiently in the case of real value

- data

-

-Another important reason is *redundancy*, that is, because data are real the Fourier transform comes with a special symmetry:

-

-.. math::

- z_k = z_{N-k}^*

-

-so you don't need to store N complex values but only N/2.

-

-A sequence with this symmetry is called "conjugate-complex" or

-"half-complex". The results of a Fourier Transform of real data will be an half-complex array. In GSL Shell half-complex array are treated like special kind of arrays in order to be more efficient. To access element in half-complex array you can use the following functions:

-

-.. class:: half-complex array

-

- Store the half-complex sequence resulting from a Fourier transform on real data.

-

- .. method:: get(i)

-

- Returns the complex value of index ``i``. The index are used as follows:

-

- - 0, the value corresponding to zero frequency

- - n, with n > 0, the values corresponding to the frequency n. If the input

- is a temporal sequence of interval of index ``k`` then it corresponds

- to the term exp(-2 i k n / N)

- - -n, with n >0, the values of negative frequency -n. The resulting values

- would be always the complex conjugate of the corresponding positive

- frequency value

-

- .. method:: set(i, v)

-

- Set the complex element of index ``i`` to the value ``v``. The indexing

- rules are the same that for the method :func:`get`.

-

- .. attribute:: length

-

- The index of the maximum frequency. If N is the size of the

- real data it is always equal to N/2. It is useful since the domain of

- possible frequencies are:

-

- .. math::

- -N/2, -N/2+1, \dots, -1, 0, 1, \dots, N/2-1, N/2

-

-.. function:: fft(v)

-

- Perform the Fourier transform of a matrix column of *real*

- numbers. The transformation will be done *in place* so your data

- will be actually transformed in an half-complex array. The function

- does not return any value. If you need to preserve the original

- data you should make a copy of your vector by doing something like::

-

- local f = v:copy() -- we take a copy of the vector

- fft(f) -- fourier transform is made in place

-

- Please note that the value you obtain is not an ordinary matrix but

- an half-complex array. You can use half-complex array only with the

- methods :func:`get` and :func:`set`. If you want to have an

- ordinary matrix you can easily build it with the following

- instructions::

-

- -- we suppose that f is an half-complex array

- m = cnew(f.length+1, 1, |i,j| f:get(i-1))

-

-.. function:: fft_inv(hc)

-

- Perform the inverse Fourier transform of an half-complex array *in

- place*. As a result the input is transformed in a real valued

- column matrix. The factor resulting matrix does include the 1/N

- factor to ensure that the use of :func:`fft` and :func:`fft_inv`

- gives exactly the original data. Please note that some authors

- use a factor of :math:`1/\sqrt{N}` for both forward and inverse

- matrix.

-

- A tipical usage of :func:`fft_inv` is to revert the trasformation

- made with :func:`fft` but by doing some transformations of the

- way. So a typical usage path could be::

-

- -- we assume v is a column matrix with our data

- fft(v) -- fourier transform

-

- -- here we can manipulate the half-complex array with

- -- using the methods `get' and `set'

- some code here

-

- fft_inv(v) -- we perform the inverse fourier transform

- -- now v is again a column matrix of real numbers

-

-Fourier Transform of Complex Data

----------------------------------

-

-The Fourier transform of complex data is simpler then the one for real data because both the input and the output will be complex column matrix of the same size. The algorithm actually used will always be a mixed-radix algorithm and GSL Shell will take care of allocating the required resources. As for the real data the table allocation strategy is very efficient for the case when many fourier transforms are made always with the same size.

-

-The Fourier trasform is made using the function :func:`cfft` that provides both direct and inverse Fourier transform.

-

-.. function:: cfft(v[, sign])

-

- Perform the Fourier transform of a given complex column vector *in

- place*. The function does not return any value. If ``sign`` is

- positive a direct transformation will be done, this is the defaut

- case. If sign is negative the inverse transfrmations will be done

- and a factor 1/N will be used to scale the results.

-

-FFT example

------------

-

-In this example we will treat a square pulse in the temporal domain. To illustrate a typical example of FFT usage we perform the Fourier Transform of the signal and we cut the higher order frequencies. Than we perform the inverse transform and we compare the result with the original time signal.

-

-So, first we define our square pulse in the time domain. Actually it will be a matrix with just one column::

-

- n = 256

- ncut = 16

-

- -- we create a pulse signal in the time domain

- sq = new(n, 1, |i| i < n/3 and 0 or (i < 2*n/3 and 1 or 0))

-

-Than we create two new plots, one for the Fourier transform and one for the signal itself::

-

- pt, pf = plot(), plot()

- pt:add_line(ipath(irow(sq, function(m,i) return i, m:get(i,1) end)), 'black')

-

-Now we are ready to perform:

-

- - the Fourier transform

- - cut the higher frequencies

- - transform back the signal in the time domain

-

-and plot the results::

-

- -- a convenience function, the norm of a complex number

- cmod = |z| sqrt(z*conj(z))

-

- fft(sq)

- pf:add(ibars(iter(function(i) return i, cmod(sq:get(i)) end, 0, n/2)), 'black')

- for k=ncut, n/2 do sq:set(k,0) end

- fft_inv(sq)

-

- pt:add_line(ipath(irow(sq, function(m,i) return i, m:get(i,1) end)), 'red')

-

- pf:show()

- pt:show()

-

-.. figure:: fft-example-power-spectrum.png

-

- Fourier transform spectrum

-

-.. figure:: fft-example-time-signal.png

-

- Time signal before (black) and after (red) the transformation

-

-You can observe in the reconstructed signal (the red curve) that we obtain approximatively the square pulse but with a lot of oscillations. Of course this is an artifact of our transformations. The reason is that in order to reproduce prefectly a sharp signal we need also all the high frequencies of the Fourier transform.

+.. highlight:: lua

+

+.. include:: <isogrk1.txt>

+

+Fast Fourier Transform

+==============================

+

+Mathematical Definitions

+------------------------

+

+Fast Fourier Transforms are efficient algorithms for calculating the discrete fourier transform (DFT)

+

+.. math::

+ x_j = \sum_{k=0}^{n-1} z_k \exp(-2\pi i j k / n)

+

+The DFT usually arises as an approximation to the continuous fourier

+transform when functions are sampled at discrete intervals in space or

+time. The naive evaluation of the discrete fourier transform is a matrix-vector multiplication :math:`W\vec{z}`. A general matrix-vector multiplication takes O(n\ :sup:`2`) operations for n data-points. Fast fourier transform algorithms use a divide-and-conquer strategy to factorize the matrix W into smaller sub-matrices, corresponding to the integer factors of the length n. If n can be factorized into a product of integers f\ :sub:`1` f\ :sub:`2` ... f\ :sub:`m` then the DFT can be computed in O(n |Sgr| f\ :sub:`i`) operations. For a radix-2 FFT this gives an operation count of O(n log\ :sub:`2` n).

+

+All the FFT functions offer two types of transform: forwards and inverse, based on the same mathematical definitions. The definition of the forward fourier transform, ``fft(z)``, is,

+

+.. math::

+ x_j = \sum_{k=0}^{n-1} z_k \exp(-2\pi i j k / n)

+

+and the definition of the inverse fourier transform, ``fft_inv(x)``, is,

+

+.. math::

+ z_j = {1 \over n} \sum_{k=0}^{n-1} x_k \exp(2\pi i j k / n).

+

+The factor of 1/n makes this a true inverse.

+

+In general there are two possible choices for the sign of the

+exponential in the transform/ inverse-transform pair. GSL follows the

+same convention as fftpack, using a negative exponential for the

+forward transform. The advantage of this convention is that the

+inverse transform recreates the original function with simple fourier

+synthesis. Numerical Recipes uses the opposite convention, a positive

+exponential in the forward transform.

+

+GSL Shell interface

+-------------------

+

+When you perform a Fourier transform different routines can be used depending on the size N of the sample. GSL features two kinds of routines, radix 2 and mixed radix. The radix 2 routine can be used if the size of the sample is a power of 2 while mixed radix can be used for any size of the input data but will be "fast" only if the size can be factorized into the product of small prime numbers.

+

+GSL Shell will select automatically the appropriate routine for you in a transparent way. But you should nevertheless be aware that the calculation will be fast only if the size of data can be factorised in small primes. You should also know that the mixed radix algorithm requires additional memory spaces so if you need to minimise memory usage you should feed only data of power of two size.

+

+GSL shell manage automatically the additional memory of precomputed trigonometric tables that may be required by the routine. The strategy used will be very efficient if you will perform many times the Fourier Transform of data always of the *same* size. If you perform many Fourier transform always with different size you will incur in a performance penality. Please not anyway that *no additional workspaces* are needed with radix two transformations.

+

+Fourier Transform of Real Data

+------------------------------

+

+In principle the Fourier Transform of real data is exactly the same of the Fourier transform of complex data, we only make the difference between them because:

+

+- in GSL Shell real value matrix and complex value matrix are considered

+ different types

+- specialised routines exists to work efficiently in the case of real value

+ data

+

+Another important reason is *redundancy*, that is, because data are real the Fourier transform comes with a special symmetry:

+

+.. math::

+ z_k = z_{N-k}^*

+

+so you don't need to store N complex values but only N/2.

+

+A sequence with this symmetry is called "conjugate-complex" or

+"half-complex". The results of a Fourier Transform of real data will be an half-complex array. In GSL Shell half-complex array are treated like special kind of arrays in order to be more efficient. To access element in half-complex array you can use the following functions:

+

+.. class:: half-complex array

+

+ Store the half-complex sequence resulting from a Fourier transform on real data.

+

+ .. method:: get(i)

+

+ Returns the complex value of index ``i``. The index are used as follows:

+

+ - 0, the value corresponding to zero frequency

+ - n, with n > 0, the values corresponding to the frequency n. If the input

+ is a temporal sequence of interval of index ``k`` then it corresponds

+ to the term exp(-2 i k n / N)

+ - -n, with n >0, the values of negative frequency -n. The resulting values

+ would be always the complex conjugate of the corresponding positive

+ frequency value

+

+ .. method:: set(i, v)

+

+ Set the complex element of index ``i`` to the value ``v``. The indexing

+ rules are the same that for the method :func:`get`.

+

+ .. attribute:: length

+

+ The index of the maximum frequency. If N is the size of the

+ real data it is always equal to N/2. It is useful since the domain of

+ possible frequencies are:

+

+ .. math::

+ -N/2, -N/2+1, \dots, -1, 0, 1, \dots, N/2-1, N/2

+

+.. function:: fft(v)

+

+ Perform the Fourier transform of a matrix column of *real*

+ numbers. The transformation will be done *in place* so your data

+ will be actually transformed in an half-complex array. The function

+ does not return any value. If you need to preserve the original

+ data you should make a copy of your vector by doing something like::

+

+ local f = v:copy() -- we take a copy of the vector

+ fft(f) -- fourier transform is made in place

+

+ Please note that the value you obtain is not an ordinary matrix but

+ an half-complex array. You can use half-complex array only with the

+ methods :func:`get` and :func:`set`. If you want to have an

+ ordinary matrix you can easily build it with the following

+ instructions::

+

+ -- we suppose that f is an half-complex array

+ m = cnew(f.length+1, 1, |i,j| f:get(i-1))

+

+.. function:: fft_inv(hc)

+

+ Perform the inverse Fourier transform of an half-complex array *in

+ place*. As a result the input is transformed in a real valued

+ column matrix. The factor resulting matrix does include the 1/N

+ factor to ensure that the use of :func:`fft` and :func:`fft_inv`

+ gives exactly the original data. Please note that some authors

+ use a factor of :math:`1/\sqrt{N}` for both forward and inverse

+ matrix.

+

+ A tipical usage of :func:`fft_inv` is to revert the trasformation

+ made with :func:`fft` but by doing some transformations of the

+ way. So a typical usage path could be::

+

+ -- we assume v is a column matrix with our data

+ fft(v) -- fourier transform

+

+ -- here we can manipulate the half-complex array with

+ -- using the methods `get' and `set'

+ some code here

+

+ fft_inv(v) -- we perform the inverse fourier transform

+ -- now v is again a column matrix of real numbers

+

+Fourier Transform of Complex Data

+---------------------------------

+

+The Fourier transform of complex data is simpler then the one for real data because both the input and the output will be complex column matrix of the same size. The algorithm actually used will always be a mixed-radix algorithm and GSL Shell will take care of allocating the required resources. As for the real data the table allocation strategy is very efficient for the case when many fourier transforms are made always with the same size.

+

+The Fourier trasform is made using the function :func:`cfft` that provides both direct and inverse Fourier transform.

+

+.. function:: cfft(v[, sign])

+

+ Perform the Fourier transform of a given complex column vector *in

+ place*. The function does not return any value. If ``sign`` is

+ positive a direct transformation will be done, this is the defaut

+ case. If sign is negative the inverse transfrmations will be done

+ and a factor 1/N will be used to scale the results.

+

+FFT example

+-----------

+

+In this example we will treat a square pulse in the temporal domain. To illustrate a typical example of FFT usage we perform the Fourier Transform of the signal and we cut the higher order frequencies. Than we perform the inverse transform and we compare the result with the original time signal.

+

+So, first we define our square pulse in the time domain. Actually it will be a matrix with just one column::

+

+

+ local n, ncut = 256, 16

+

+ -- we create a pulse signal in the time domain

+ local sq = new(n, 1, |i| i < n/3 and 0 or (i < 2*n/3 and 1 or 0))

+

+Than we create two new plots, one for the Fourier transform and one for the signal itself::

+

+ local pt = plot('Original signal / reconstructed')

+ local pf = plot('FFT Power Spectrum')

+

+ pt:addline(ipath(sample(|i| sq[i], 1, n, n-1)), 'black')

+

+Now we are ready to perform:

+

+ - the Fourier transform

+ - cut the higher frequencies

+ - transform back the signal in the time domain

+

+and plot the results::

+

+ fft(sq)

+

+ pf:add(ibars(sample(|k| abs(sq:get(k)), 0, 60, 60)), 'black')

+

+ for k=ncut, n/2 do sq:set(k,0) end

+ fft_inv(sq)

+

+ pt:addline(ipath(sample(|i| sq[i], 1, n, n-1)), 'red')

+

+ pf:show()

+ pt:show()

+

+.. figure:: fft-example-power-spectrum.png

+

+ Fourier transform spectrum

+

+.. figure:: fft-example-time-signal.png

+

+ Time signal before (black) and after (red) the transformation

+

+You can observe in the reconstructed signal (the red curve) that we obtain approximatively the square pulse but with a lot of oscillations. Of course this is an artifact of our transformations. The reason is that in order to reproduce prefectly a sharp signal we need also all the high frequencies of the Fourier transform.

diff --git a/doc/source/graphics.rst b/doc/source/graphics.rst
index 49cd3362..8714c709 100644
--- a/doc/source/graphics.rst
+++ b/doc/source/graphics.rst

@@ -24,17 +24,18 @@ Let's start with a simple example, let us suppose that we want to plot the funct

where |agr| and |ohgr| are constants and t vary from 0 to t1. We can plot this function with GSL Shell with the following instructions::

- require 'draw'

-

function myplot(alpha, omega, t1)

- local p = plot() -- create a new plot, it is not shown for the moment

+ -- create a new plot, it is not shown for the moment

+ local p = plot('f(x) = exp(-a t) sin(w t)')

-- we create a line that corresponds to our function

local ln = fxline(|t| exp(-alpha*t)*sin(omega*t), 0, t1)

-- we add the line to the plot and show it

- p:add_line(ln, 'red')

+ p:addline(ln, 'red')

p:show()

+

+ return p

end

Then to plot something you have just to call the 'myplot' function. For example::

@@ -45,10 +46,15 @@ Then to plot something you have just to call the 'myplot' function. For example:

The function :func:`fxline` takes three arguments, the function to plot and the initial anf final values of the variable. By default the function will be sampled with 256 points but if you want you can provide a fourth arguments to give the number of sample points.

-TODO : Explain how 'plot' works.

+In this example we have used the :func:`plot` function to create a plot, the :func:`fxline` function to create the line to draw and the method :func:`addline` to add the line to the plot (in red). These three operations can be done with a single function, :func:`fxplot`. It works like that::

+

+ p = fxplot(|x| sin(x), 0, 8*pi)

-A simpler example

------------------

+where the first arguments is the function to plot and the following

+arguments are the extrema of variation of the indipendent variable x.

+

+Graphics primitives

+-------------------

In order to better understand the way GSL shell graphics works it is better to take a step back. Let use suppose that we want to plot an equilateral triangle. We can proceed as follows:

- define a 'path' that describe the countour that we want to plot

@@ -70,7 +76,7 @@ So to plot a triangle you can give the following instructions::

.. figure:: simpler-example-1.png

-Please not that we have used the :func:`add` method instead of :func:`add_line` to add the path.

+Please not that we have used the :func:`add` method instead of :func:`addline` to add the path.

Now let us suppose that we want to plot only the contour of the triangle with a line 10 pixel thick and with round edges. Then what you have to do is to supply to the :func:`add` method a third argument where you specify a ``stroke`` transformation::

@@ -105,20 +111,44 @@ Some transformations are naturally expressed as post-transforms because they doe

plot

----

-We have seen in the previous paragraph that you can add more graphical elements in a plot by using the methods :func:`add` and :func:`add_line`. The method :func:`add_line` is just a shortcut to add elements with a 'stroke' post transform with unitary width.

+We have seen in the previous paragraph that you can add more graphical elements in a plot by using the methods :func:`add` and :func:`addline`. The method :func:`addline` is just a shortcut to add elements with a 'stroke' post transform of unitary width.

If can add elements to a plot in any moments even when it is already shown. GSL Shell will automatically calculate the bounding box so that every elements is shown on the window.

.. class:: plot

- .. function:: plot([show_grid])

+ .. function:: plot([title])

- Create a new empty plot. By default ``show_grid`` is true and the plot is shown with axes, grids and marks.

+ Create a new empty plot with an optional title.

.. method:: add(obj, color[, post_trans, pre_trans])

- Add the :ref:`graphical object <graphics-objects>` ``obj`` to the plot with the given ``color``.

- The optional arguments ``post_trans`` and ``pre_trans`` should be a table of :ref:`graphical transformations <graphics-transforms>`.

+ Add the :ref:`graphical object <graphics-objects>` ``obj`` to

+ the plot with the given ``color``. The optional arguments

+ ``post_trans`` and ``pre_trans`` should be a table of

+ :ref:`graphical transformations <graphics-transforms>`.

+

+ .. method:: addline(obj, color[, post_trans, pre_trans])

+

+ Add the :ref:`graphical object <graphics-objects>` ``obj`` to

+ the plot by performing automatically a stroke of it. It is

+ useful because you often need to draw lines and not filled

+ polygons. It is equivalent to add a 'stroke' operations of

+ unitary size in the viewport coordinates system.

+

+ .. method:: update()

+

+ Updates the window that display the plot.

+

+ .. attribute:: units

+

+ A boolean value that define if the axis and grids should be

+ drawn or not. By default it is true.

+

+ .. attribute:: title

+

+ The title of the plot. You can change or set the title using

+ this attribute.

.. _graphics-objects:

@@ -178,7 +208,8 @@ Graphical Objects

.. method:: set_point(x, y)

- Set the position where the test is diplayed. It corresponds to the bottom left corner of the text.

+ Set the position where the test is diplayed. It corresponds to

+ the bottom left corner of the text.

.. method:: rotate(angle)

diff --git a/doc/source/index.rst b/doc/source/index.rst
index 47559aff..88692f3a 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst

@@ -1,47 +1,48 @@

-##################################

-Welcome to GSL shell documentation

-##################################

-

-GSL shell is an interactive command line interface that gives easy access

-to the GNU Scientific Library (GSL) collection of mathematical methods for

-numerical computations.

-

-GSL shell can be used interactively to perform calculations with matrices or

-vectors but it does allow also to write complex user defined functions with

-the Lua scripting interpreter.

-

-Lua is a very interesting and easy to learn scripting language that features

-advanced functionalities like closures and metamethods. Lua is very

-easy to learn and will give you the power of defining your own complex routines

-to use the GSL routines more effectively.

-

-GSL Shell is hosted at `Savannah <http://savannah.nongnu.org>`_, here is the `project page <https://savannah.nongnu.org/projects/gsl-shell/>`_.

-

-You can download the first developer preview release of GSL Shell in the `download page <https://savannah.nongnu.org/files/?group=gsl-shell>`_. You will find the source code, a binary package that require GSL 1.13 and readline5 and a DEBIAN package for ubuntu.

-

-GSL Shell is currently developed by Francesco Abbate. You can contact him at the followng address:

-

-.. image:: email-gslshell.png

-

-Contents:

-

-.. toctree::

- :maxdepth: 2

-

- news.rst

- intro.rst

- matrices.rst

- linalg.rst

- random.rst

- randist.rst

- pdf.rst

- cdf.rst

- nlinfit.rst

- fft.rst

- ode.rst

- integ.rst

- sf.rst

- minim.rst

- graphics.rst

- examples.rst

- acknowledgments.rst

+##################################

+Welcome to GSL shell documentation

+##################################

+

+GSL shell is an interactive command line interface that gives easy access

+to the GNU Scientific Library (GSL) collection of mathematical methods for

+numerical computations.

+

+GSL shell can be used interactively to perform calculations with matrices or

+vectors but it does allow also to write complex user defined functions with

+the Lua scripting interpreter.

+

+Lua is a very interesting and easy to learn scripting language that features

+advanced functionalities like closures and metamethods. Lua is very

+easy to learn and will give you the power of defining your own complex routines

+to use the GSL routines more effectively.

+

+GSL Shell is hosted at `Savannah <http://savannah.nongnu.org>`_, here is the `project page <https://savannah.nongnu.org/projects/gsl-shell/>`_.

+

+You can download the first developer preview release of GSL Shell in the `download page <https://savannah.nongnu.org/files/?group=gsl-shell>`_. You will find the source code, a binary package that require GSL 1.13 and readline5 and a DEBIAN package for ubuntu.

+

+GSL Shell is currently developed by Francesco Abbate. You can contact him at the followng address:

+

+.. image:: email-gslshell.png

+

+Contents:

+

+.. toctree::

+ :maxdepth: 2

+

+ news.rst

+ intro.rst

+ matrices.rst

+ linalg.rst

+ random.rst

+ randist.rst

+ pdf.rst

+ cdf.rst

+ nlinfit.rst

+ fft.rst

+ ode.rst

+ integ.rst

+ sf.rst

+ minim.rst

+ graphics.rst

+ contour.rst

+ examples.rst

+ acknowledgments.rst

diff --git a/doc/source/minim.rst b/doc/source/minim.rst
index 970f5c72..81619788 100644
--- a/doc/source/minim.rst
+++ b/doc/source/minim.rst

@@ -24,7 +24,7 @@ n-dimensional functions. The algorithms proceed from an initial guess

using a search algorithm which attempts to move in a downhill direction.

Algorithms making use of the gradient of the function perform a

-one-dimensional line minimisation along this direction until the lowest

+one-dimensional line minimization along this direction until the lowest

point is found to a suitable tolerance. The search direction is then

updated with local information from the function and its derivatives,

and the whole process repeated until the true n-dimensional minimum is

@@ -46,7 +46,7 @@ main phases of the iteration. The steps are,

* test S for convergence, and repeat iteration if necessary

Each iteration step consists either of an improvement to the

-line-minimisation in the current direction or an update to the search

+line-minimization in the current direction or an update to the search

direction itself. At any time you can query the minimizer object to ask about the current best estimates of the point or to have other informations about its current state.

Creating a minimizer

@@ -72,7 +72,7 @@ Let us suppose that we want to minimize the following function:

.. math::

f(x, y) = e^x \left( 4 x^2 + 2 y^2 + 4 x y + 2 y + 1 \right)

-If we want to use a minimisation algorithm with derivatives whe should provide the derivatives of the function. These are easily calculated:

+If we want to use a minimization algorithm with derivatives whe should provide the derivatives of the function. These are easily calculated:

.. math::

{\partial f \over \partial x} = e^x \left( 4 x^2 + 2 y^2 + 4 x y + 8 x + 6 y + 1 \right)

@@ -98,3 +98,88 @@ Once we have defined the function we can create the minimizer and set the starti

x0 = vector {-0.5, 1.0} -- the starting point

m = minimizer {fdf= fex, n= 2}

m:set(x0, vector {1, 1})

+

+Then you can use the method :func:`step` on the object ``m`` to make a single "step" with the minimization algorithm and we check the returned value (a string) to know if the algorithm succeded or if more iterations are needed::

+

+ while m:step() == 'continue' do

+ print(tr(m.x))

+ end

+ print('Found minimum at: ', tr(m.x))

+ print('with function value: ', m.value)

+

+Here a plot of the path followed by the minimizer to find the solution inside a contour plot of the functions.

+

+.. figure:: minimizer-example-screenshot.png

+

+

+Minimizer class definition

+--------------------------

+

+Actually there are two minimizer type, those that use the function with derivative and those that does not. For the GSL shell users this difference is not very important because the appropriate type of object is created automatically by thefunction :func:`minimizer` depending on its arguments. The two class also have the same methods, so once the minimizer is created, you can use it in the same way regardless of its nature.

+

+.. class:: fdfmultimin

+

+ Multiple parameters minimizer class that make use of the

+ derivatives to find the minimum.

+

+ .. method:: set(x0, step)

+

+ Set the minimizer current point to ``x0`` and the step size of the

+ search to ``step``. This latter parameter can be also a vector

+ but in this case the geometric mean will be used as a step size.

+

+ .. method:: step()

+

+ Advance the minimizer of a single step. It does return

+ ``continue`` if it did not reach the optimal point or

+ ``success`` if it does succeded to find the minimum.

+

+ .. method:: run()

+

+ Advance the minimizer until a minimum is found or raise an error

+ if it cannot locate any minimum.

+

+ .. attribute:: x

+

+ Return the current estimate of the search point.

+

+ .. attribute:: value

+

+ Return the value of the function at the current position.

+

+ .. attribute:: gradient

+

+ Return the gradient of the function at the current position.

+

+.. class:: fmultimin

+

+ Minimizer class that does not use the derivatives. Its interface is

+ the same of the :class:`fdfmultimin`.

+

+ .. method:: set(x0, step)

+

+ Set the minimizer current point to ``x0`` and the step size of

+ the search to vector ``step``. This latter parameter should be a

+ vector of the same size of ``x0`` and each of its components

+ rapresents the step size in the corresponding direction. If a

+ number is provided the step size are taken to its value for all

+ the directions.

+

+ .. method:: step()

+

+ Advance the minimizer of a single step. It does return

+ ``continue`` if it did not reach the optimal point or

+ ``success`` if it does succeded to find the minimum.

+

+ .. method:: run()

+

+ Advance the minimizer until a minimum is found or raise an error

+ if it cannot locate any minimum.

+

+ .. attribute:: x

+

+ Return the current estimate of the search point.

+

+ .. attribute:: value

+

+ Return the value of the function at the current position.

diff --git a/doc/source/minimizer-example-screenshot.png b/doc/source/minimizer-example-screenshot.png
new file mode 100644
index 00000000..36117319
--- /dev/null
+++ b/doc/source/minimizer-example-screenshot.png

diff --git a/doc/source/nlinfit-example-plot.png b/doc/source/nlinfit-example-plot.png
new file mode 100644
index 00000000..3020bb9f
--- /dev/null
+++ b/doc/source/nlinfit-example-plot.png

diff --git a/doc/source/nlinfit.rst b/doc/source/nlinfit.rst
index df92f087..1b88ba2b 100644
--- a/doc/source/nlinfit.rst
+++ b/doc/source/nlinfit.rst

@@ -60,23 +60,37 @@ In order to perform a non linear fitting with GSL Shell you should use a *solver

Here an example::

local n = 50

- local p = {a= -3.1, A= 1.55}

- local y = new(n, 1, function (i,j) return p.A * exp(p.a * (i-1)/n) end)

+ local px = vector {1.55, -1.1, 12.5}

+ local p0 = vector {2.5, -1.5, 5.3}

+ local xs = |i| (i-1)/n

+ local r = rng()

+

+ local fmodel = function(p, t, J)

+ local e, s = exp(p[2] * t), sin(p[3] * t)

+ if J then

+ J:set(1,1, e * s)

+ J:set(1,2, t * p[1] * e * s)

+ J:set(1,3, t * p[1] * e * cos(p[3] * t))

+ end

+ return p[1] * e * s

+ end

+

+ local y = new(n, 1, |i,j| fmodel(px, xs(i)) * (1 + rnd.gaussian(r, 0.1)))

local function expf(x, f, J)

- for k= 1, n do

- local t = (k-1) / n

- local A, a = x[1], x[2]

- local e = exp(a * t)

- if f then f:set(k, 1, A * e - y[k]) end

- if J then

- J:set(k, 1, e)

- J:set(k, 2, t * A * e)

- end

+ for k=1, n do

+ local ym = fmodel(x, xs(k), J and J:row(k))

+ if f then f:set(k, 1, ym - y[k]) end

end

end

- s = solver {fdf= expf, n= n, p= 2, x0= vector {3.5, -2.5}}

+ pl = plot('Non-linear fit / A * exp(a t) sin(w t)')

+ pl:addline(ipath(sequence(function(k) return xs(k), y[k] end, n)), 'blue',

+ {{'marker', size= 5}})

+

+ s = solver {fdf= expf, n= n, p= 3, x0= p0}

+

+ pl:addline(fxline(|x| fmodel(s.x, x), 0, xs(n)), 'red', {{'dash', a=7, b=3}})

repeat

print_state (s)

@@ -84,9 +98,12 @@ Here an example::

until status ~= 'continue'

print_state (s)

+ pl:addline(fxline(|x| fmodel(s.x, x), 0, xs(n)), 'red')

+ pl:show()

+

where the function ``print_state`` could be defined like::

- unction print_state(s)

+ function print_state(s)

print ("x: ", tr(s.x))

print ("chi square: ", prod(s.f, s.f)[1])

end

@@ -94,17 +111,28 @@ where the function ``print_state`` could be defined like::

The output you obtain is::

- x: 3.5, -2.5

- chi square: 42.641185290635

- x: 1.53499, -2.71813

- chi square: 0.14091408323034

- x: 1.54408, -3.05951

- chi square: 0.0008951390111855

- x: 1.54994, -3.09966

- chi square: 5.5968773129138e-08

- x: 1.55, -3.1

- chi square: 2.0533525307766e-16

-

+ x: , [ 2.5 -1.5 5.3 ]

+ chi square: , 61.909477682545

+ x: , [ 0.816847 -2.19811 5.30633 ]

+ chi square: , 24.637775808867

+ x: , [ 1.1919 -5.81962 5.71798 ]

+ chi square: , 20.698635047305

+ x: , [ 2.55001 -11.3184 11.1346 ]

+ chi square: , 15.387167949514

+ [ ... ]

+ x: , [ 1.58178 -1.22193 12.5912 ]

+ chi square: , 0.34786217353905

+ x: , [ 1.56791 -1.14061 12.5125 ]

+ chi square: , 0.30630801846857

+ x: , [ 1.56791 -1.14019 12.5156 ]

+ chi square: , 0.30626868109332

+ x: , [ 1.56791 -1.1402 12.5156 ]

+ chi square: , 0.3062686809533

+

+.. figure:: nlinfit-example-plot.png

+

+ Non-linear fit of function A exp(a t) sin(w t) with gaussian noise

+

Solver class definition

-----------------------

@@ -134,7 +162,7 @@ Solver class definition

.. method:: iterate()

Advance the solver of a single step. It does return ``continue`` if it

- did not reach the optimal point and ``terminated`` otherwise.

+ did not reach the optimal point and ``success`` otherwise.

.. method:: run([maxiter])

diff --git a/examples/fft.lua b/examples/fft.lua
index 1d3161ce..9d3dec01 100644
--- a/examples/fft.lua
+++ b/examples/fft.lua

@@ -6,15 +6,22 @@ function demo1()

local sq = new(n, 1, |i| i < n/3 and 0 or (i < 2*n/3 and 1 or 0))

- local p = plot()

- p:addline(ipath(sample(|i| sq[i], 1, n, n-1)), 'black')

+ local pt = plot('Original signal / reconstructed')

+ local pf = plot('FFT Power Spectrum')

+

+ pt:addline(ipath(sample(|i| sq[i], 1, n, n-1)), 'black')

fft(sq)

+

+ pf:add(ibars(sample(|k| abs(sq:get(k)), 0, 60, 60)), 'black')

+

for k=ncut, n/2 do sq:set(k,0) end

fft_inv(sq)

- p:addline(ipath(sample(|i| sq[i], 1, n, n-1)), 'red')

- p:show()

+ pt:addline(ipath(sample(|i| sq[i], 1, n, n-1)), 'red')

+

+ pf:show()

+ pt:show()

return p

end

diff --git a/examples/minimize.lua b/examples/minimize.lua
index 7d02c761..3fa4fd5a 100644
--- a/examples/minimize.lua
+++ b/examples/minimize.lua

@@ -45,7 +45,7 @@ function demo1()

local px = new(2,1)

p = contour(cook(frosenbrock), {-1.5, -0.2}, {1.5, 2}, 20, 20, 12)

- c=path(m.x[1], m.x[2])

+ c = path(m.x[1], m.x[2])

while m:step() == 'continue' do

c:line_to(m.x[1], m.x[2])

end

@@ -55,6 +55,7 @@ function demo1()

p:addline(c, 'black', {{'marker', size=5}})

p:addline(c, 'red')

+ p.title = 'Rosenbrock function minimisation'

p:show()

return p

end

@@ -65,7 +66,7 @@ function demo2()

m = minimizer {f= f, n= 2}

m:set(x0, vector {1, 1})

- p=contour(cook(f), {-2, -2}, {8, 2})

+ p=contour(cook(f), {-2, -3}, {8, 2})

c=path(m.x[1], m.x[2])

while m:step() == 'continue' do

c:line_to(m.x[1], m.x[2])

@@ -76,6 +77,7 @@ function demo2()

p:addline(c, 'black', {{'marker', size=5}})

p:addline(c, 'red')

+ p.title = 'Quadratic function minimisation'

p:show()

return p

end

@@ -86,7 +88,7 @@ function demo3()

m = minimizer {fdf= fex, n= 2}

m:set(x0, vector {1, 1})

- p=contour(cook(fex), {-2, -2.5}, {1, 1}, 30, 30, 16)

+ p=contour(cook(fex), {-2, -2.5}, {1, 1.5}, 30, 30, 22)

c=path(m.x[1], m.x[2])

while m:step() == 'continue' do

c:line_to(m.x[1], m.x[2])

@@ -97,6 +99,7 @@ function demo3()

p:addline(c, 'black', {{'marker', size=5}})

p:addline(c, 'green')

+ p.title = 'function minimisation: f(x,y) = 4 x^2 + 2 y^2 + 4 x y + 2 y + 1'

p:show()

return p

end