gsl-shell.git - gsl-shell

diff --git a/gslext.lua b/gslext.lua
index ef77265b..ff0d696e 100644
--- a/gslext.lua
+++ b/gslext.lua

@@ -2,5 +2,6 @@

require('base')

require('matrix')

+require('misc')

require('integ')

require('draw')

diff --git a/lm_test_enso.lua b/lm_test_enso.lua
index 15dd25b6..761eb717 100644
--- a/lm_test_enso.lua
+++ b/lm_test_enso.lua

@@ -1,5 +1,3 @@

-local template = require 'template'

-

local sin, cos, exp = math.sin, math.cos, math.exp

local pi = math.pi

@@ -258,21 +256,21 @@ local function enso_model_f(x, t)

return y

end

-lm = template.load('num/lmfit.lua.in', {N= enso_N, P= enso_P})

+local s = gsl.nlinfit {n= enso_N, p= enso_P}

-lm.set(enso_fdf, enso_x0)

+s:set(enso_fdf, enso_x0)

-print(gsl.tr(lm.x))

+print(gsl.tr(s.x))

for i=1, 80 do

- lm.iterate()

- print('ITER=', i, ': ', gsl.tr(lm.x))

- if lm.test(0, 1e-7) then print('solution found'); break end

+ s:iterate()

+ print('ITER=', i, ': ', gsl.tr(s.x), s.chisq)

+ if s:test(0, 1e-7) then print('solution found'); break end

end

p = graph.plot()

pts = graph.ipath(gsl.sequence(function(i) return i, enso_F[i] end, enso_N))

-fitln = graph.fxline(function(t) return enso_model_f(lm.x, t) end, 0, 168, 512)

+fitln = graph.fxline(function(t) return enso_model_f(s.x, t) end, 0, 168, 512)

p:addline(pts, 'blue', {{'marker', size=4}})

p:addline(fitln)

p.title = 'ENSO non-linear fit NIST test'

diff --git a/lm_test_thurber.lua b/lm_test_thurber.lua
index 60b5a711..461b6516 100644
--- a/lm_test_thurber.lua
+++ b/lm_test_thurber.lua

@@ -1,5 +1,3 @@

-local template = require 'template'

-

local sin, cos, exp = math.sin, math.cos, math.exp

local pi = math.pi

@@ -127,21 +125,21 @@ local function thurber_model_f(x, t)

return num / den

end

-lm = template.load('num/lmfit.lua.in', {N= thurber_N, P= thurber_P})

+local s = gsl.nlinfit {n= thurber_N, p= thurber_P}

-lm.set(thurber_fdf, thurber_x0)

+s:set(thurber_fdf, thurber_x0)

-print(gsl.tr(lm.x))

+print(gsl.tr(s.x), s.chisq)

for i=1, 200 do

- lm.iterate()

- print('ITER=', i, ': ', gsl.tr(lm.x))

- if lm.test(0, 1e-7) then print('solution found'); break end

+ s:iterate()

+ print('ITER=', i, ': ', gsl.tr(s.x), s.chisq)

+ if s:test(0, 1e-7) then print('solution found'); break end

end

p = graph.plot()

pts = graph.ipath(gsl.sequence(function(i) return thurber_t[i], thurber_F[i] end, thurber_N))

-fitln = graph.fxline(function(t) return thurber_model_f(lm.x, t) end, -3.1, 2.2, 512)

+fitln = graph.fxline(function(t) return thurber_model_f(s.x, t) end, -3.1, 2.2, 512)

p:addline(pts, 'blue', {{'marker', size=4}})

p:addline(fitln)

p.title = 'Thurber non-linear fit NIST test'

diff --git a/lua-lmtest.lua b/lua-lmtest.lua
index fe4fdd4a..46c8944d 100644
--- a/lua-lmtest.lua
+++ b/lua-lmtest.lua

@@ -26,13 +26,12 @@ r:set(0)

yrf = gsl.new(n, 1, function(i) return A * exp(-lambda*(i-1)) + b + gsl.rnd.gaussian(r, 0.1) end)

sigrf = gsl.new(n, 1, function() return 0.1 end)

-local template = require 'template'

-lm = template.load('num/lmfit.lua.in', {N= n, P= 3})

+local s = gsl.nlinfit {n= n, p= 3}

-lm.set(fdf, gsl.vector {1, 0, 0})

-print(gsl.tr(lm.x), gsl.prod(lm.f, lm.f))

+s:set(fdf, gsl.vector {1, 0, 0})

+print(gsl.tr(s.x), s.chisq)

for i=1, 10 do

- lm.iterate()

- print('ITER=', i, ': ', gsl.tr(lm.x), sqrt(gsl.prod(lm.f, lm.f)[1]))

+ s:iterate()

+ print('ITER=', i, ': ', gsl.tr(s.x), s.chisq)

end

diff --git a/matrix.lua b/matrix.lua
index 8d7e6319..d4f4e08b 100644
--- a/matrix.lua
+++ b/matrix.lua

@@ -235,74 +235,3 @@ function gsl.linfit(gener, x, y, w)

end

return f, c

end

-

-local function nlinfitwrapf(fmodel, x, y)

- local n = gsl.dim(y)

- return function(p, f, J)

- for k=1, n do

- local ym = fmodel(p, x[k], J and J:row(k))

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

- end

- end

-

-end

-

-function gsl.nlinfit(f, x, y, p0)

- local N = gsl.dim(y)

- local P = gsl.dim(p0)

- local s = gsl.nlfsolver {fdf= nlinfitwrapf(f, x, y), n= N, p0= p0}

- s:run()

- return function(x) return f(s.p, x) end, s.p

-end

-

-function gsl.ode(spec)

- local required = {N= 'number', eps_abs= 'number'}

- local defaults = {eps_rel = 0, a_y = 1, a_dydt = 0}

- local is_known = {rkf45= true, rk8pd= true}

-

- for k, tp in pairs(required) do

- if type(spec[k]) ~= tp then

- error(string.format('parameter %s should be a %s', k, tp))

- end

- end

- for k, v in pairs(defaults) do

- if not spec[k] then spec[k] = v end

- end

-

- local method = spec.method and spec.method or 'rkf45'

- if not is_known[method] then error('unknown ode method: ' .. method) end

- spec.method = nil

-

- local ode = template.load(string.format('num/%s.lua.in', method), spec)

-

- local ode_methods = {

- evolve = function(s, f, t)

- s._sync = false

- return ode.evolve(s._state, f, t)

- end,

- init = function(s, t0, h0, f, ...)

- s._sync = false

- return ode.init(s._state, t0, h0, f, ...)

- end

- }

-

- local ODE = {__index= function(s, k)

- if k == 't' then return s._state.t end

- if k == 'y' then

- if not s._sync then

- for k=1, s.dim do

- s._y[k] = s._state.y[k-1]

- end

- s._sync = true

- end

- return s._y

- end

- return ode_methods[k]

- end}

-

- local dim = spec.N

- local solver = {_state = ode.new(), dim= dim, _y = gsl.new(dim,1)}

- setmetatable(solver, ODE)

-

- return solver

-end