#!/usr/bin/python
import numpy

data = numpy.loadtxt("baccam.txt")

def eclipse(x,t)
	b=par(1); k=par(2); d=par(3); p=par(4); c=par(5);
	V = 1; T = 2; E = 3; I = 4; D = 5;

	xdot = zeros(5,1);
	xdot(V) = p*x(I) - c*x(V);
	xdot(T) = -b*x(T)*x(V);
	xdot(E) = b*x(T)*x(V) - k*x(E);
	xdot(I) = k*x(E) - d*x(I);
	xdot(D) = d*x(I);

def solver(st,pp)
	global par;
	par = 10.^pp;
	# x0 = [ V0; T0; E0; I0; D0 ]
	x0 = [ par(end); 4.0e8; 0; 0; 0 ];
	y = lsode("eclipse", x0, [0:7] );
	y = log10(y(2:end,1));
	# This line is to deal with data below the detection limit.
	# It should be commented out if there is no such data.
	y = max(y,0.5*ones(7,1));

function ssr = wrapper(pp)
	global ldata;
	y = solver([],pp);
	10.^pp
	ssr = sum( (y-log10(ldata(:,2))).^2 )
endfunction;

pars = []; ssrs = [];
for i = 1:6 # For each patient data set
	# par = [ b, k, d, p, c, V0 ]
	p = log10([ 3.2e-5 4.0 5.2 4e-2 5.2 0.1 ]); # log ensures pars always >= 0
	# Selecting the 6 lines corresponding to patient "i"
	ldata = data(7*i-6:7*i,:);
	# Fitting the data
	[bfpar,ssr] = nelder_mead_min("wrapper",p);
	# Storing best pars and ssr for each patient in array
	pars = [pars; 10.^bfpar ]
	ssrs = [ssrs; ssr ]
end;
printf("%g %g %g %g %g %g\n",pars');
printf("ssrs = [ ");printf("%g ",ssrs);printf("]");

# If you want to save the resulting best
# fits to file put (1), otherwise put (0)
if( 0 )
	times = [0:0.1:7]';
	for pt = 1:6
		par = pars(i,:);
		x0 = [ par(end); 4.0e8; 0; 0; 0 ];
		fd = fopen(["eclipse-pt" num2str(pt) ".dat"], "w");
		fprintf( fd, "# Eclipse model\n\n# Patient %d\n", pt );
		fprintf( fd, "# ssr = %g\n# b k d p c V0\n", ssrs(pt) );
		fprintf( fd, "# par = [ %g %g %g %g %g %g ];\n", par );
		yd = lsode("eclipse", x0, times );
		fprintf( fd, "\n# t(d) V T E I D\n" );
		fprintf( fd, "%g %g %g %g %g %g\n", [times, yd]' );
		fclose(fd);
	end;
end;
