#!/usr/bin/env ruby
require("rbgsl")

FFF = GSL::Function.alloc { |t, params|
  a = params[0]; lambda = params[1]; b = params[2]
  a*Math::exp(-lambda*t) + b
}
  
procf = Proc.new { |x, t, y, sigma, f|
  a = x[0]; lambda = x[1]; b = x[2]
  FFF.set_params(x)
  n = t.size
  for i in 0...n do
    yi = FFF.eval(t[i])
    f[i] = (yi - y[i])/sigma[i]
  end
}

procdf = Proc.new { |x, t, y, sigma, jac|
  a = x[0]; lambda = x[1]
  n = t.size
  for i in 0...n do
    ti = t[i]
    si = sigma[i]
    ei = Math::exp(-lambda*ti)
    jac[i][0] = ei/si
    jac[i][1] = -ti*a*ei/si
    jac[i][2] = 1.0/si
  end
}

n = 20
np = 3

f = GSL::MultiFit::Function_fdf.alloc(procf, procdf, np)

r = GSL::Rng.alloc()
t = GSL::Vector.alloc(n)
y = GSL::Vector.alloc(n)
sigma = GSL::Vector.alloc(n)
File.open("expdata.dat", "w") do |fp|
  for i in 0...n do
    t[i] = i
    y[i] = 1.0 + 5*Math::exp(-0.1*t[i]) + r.gaussian(0.2)
    sigma[i] = 0.2
    fp.printf("%d %g %g\n", t[i], y[i], sigma[i])
#    fp.printf("%d %g\n", t[i], y[i])
  end
end

f.set_data(t, y, sigma)
x = GSL::Vector.alloc([1.0, 0.0, 0.0])

solver = GSL::MultiFit::FdfSolver.alloc(GSL::MultiFit::FdfSolver::LMSDER, n, np)
#solver = GSL::MultiFit::FdfSolver.alloc(GSL::MultiFit::FdfSolver::LMDER, n, np)
solver.set(f, x)

iter = 0
solver.print_state(iter)
begin
  iter += 1
  status = solver.iterate
  solver.print_state(iter)
  status = solver.test_delta(1e-4, 1e-4)
end while status == GSL::CONTINUE and iter < 500

covar = solver.covar(0.0)
position = solver.position

chi2 = GSL::pow_2(solver.f.dnrm2)
dof = n - np
printf("A      = %.5f +/- %.5f\n", position[0], Math::sqrt(chi2/dof*covar[0][0]))
printf("lambda = %.5f +/- %.5f\n", position[1], Math::sqrt(chi2/dof*covar[1][1]))
printf("b      = %.5f +/- %.5f\n", position[2], Math::sqrt(chi2/dof*covar[2][2]))

str = sprintf("%4.3f*exp(-%4.3f*x)+%4.3f", position[0], position[1], position[2])

FFF.set_params(position)
x = 0
File.open("fit.dat", "w") do |f|
  while x < 40
    f.printf("%e %e\n", x, FFF.eval(x))
    x += 0.1
  end
end

system("graph -T X -C -g 3 -L '#{str}' -X x -I e -m -1 -S 4 expdata.dat -I a -m 2 -S 0 fit.dat")
#File.delete("expdata.dat")
#File.delete("fit.dat")