function [bestP,min_cost] = recover 

clear all
warning off

options = optimset('TolX', 1e-8, 'TypicalX', [1e-10, 1e-10]);
m0 = .5;
n0 = m0 * 5e-2;
alpha = 0.99;  %pulled from  model2d    
Ptrue = 1e-3;
ktrue = 1e-6;
maxt=180;
t = linspace(0,maxt,4);
err = 1e-3;

m.true = (Ptrue/ktrue + (m0 - Ptrue/ktrue).*exp(-ktrue*t));
n.true = n0*exp(-ktrue*alpha*t);
data.m = (Ptrue/ktrue + (m0 - Ptrue/ktrue).*exp(-ktrue*t)).*(1 + randn(4,1)'*err);
data.n = n0*exp(-ktrue*alpha*t).*(1 + randn(4,1)'*err*0.01);

extra.m.m0 = m0;
extra.n.n0 = n0;
extra.m.data = data.m;
extra.n.data = data.n;
extra.t = t;
extra.alpha = alpha;

p0 = [0,0];

k0 = -1 * polyfit(t,log(data.n),1); %linear regression on 13C
extra.k0 = k0(1);
[bestP0,min_cost0]=fminsearch(@cost0,[0],options,extra);
bestP0(2) = extra.k0;
[bestP1,min_cost1]=fminsearch(@cost1,p0,options,extra);
[bestP2,min_cost2]=fminsearch(@cost2,p0,options,extra);

bestP = [Ptrue, ktrue, bestP0 , bestP1 , bestP2]';
min_cost = [min_cost0 , min_cost1 , min_cost2];

t2=linspace(0,maxt);
m0recov = bestP0(1)/bestP0(2) + (m0 - bestP0(1)/bestP0(2)).*exp(-bestP0(2)*t2);
n0recov = n0*exp(-bestP0(2)*alpha*t2);

m1recov = bestP1(1)/bestP1(2) + (m0 - bestP1(1)/bestP1(2)).*exp(-bestP1(2)*t2);
n1recov = n0*exp(-bestP1(2)*alpha*t2);

m2recov = bestP2(1)/bestP2(2) + (m0 - bestP2(1)/bestP2(2)).*exp(-bestP2(2)*t2);
n2recov = n0*exp(-bestP2(2)*alpha*t2);

figure;
plot(t,m.true,':',t,data.m, 'o', t2,m1recov,'-b',t2,m2recov,'-r',t2,m0recov,'-g');
legend('True','Data','L^1','L^2','LR');
title('12C recovery from data');
figure;
plot(t,n.true,':',t,data.n,'o',t2,n1recov,'-b',t2,n2recov,'-r',t2,n0recov,'-g');
legend('True','Data','L^1','L^2','LR');
title('13C recovery from data');
%-----------------------------
function errlr = cost0(param,extra);

P = param(1);
k = extra.k0;
t = extra.t;
data.m = extra.m.data;   %C12 data
m0 = extra.m.m0;        %C12 orig
data.n = extra.n.data;   %C13 data
data.n = extra.n.data;   %C13 data
n0 = extra.n.n0;        %C13 orig
alpha = extra.alpha;     %fractionization constant

M = P/k + (m0 - P/k).*exp(-k*t);
N = n0*exp(-k*alpha*t) ;
%AP =  N ./ (M + N) * 100;
M = M(:);
N = N(:);
m = data.m(:);
n = data.n(:);

errlr = .95*sum(abs(M-m)) + .99*sum(abs(N - n));

%-----------------------------
function l1error = cost1(param,extra);


P = param(1);
k = param(2);
t = extra.t;
data.m = extra.m.data;   %C12 data
m0 = extra.m.m0;        %C12 orig
data.n = extra.n.data;   %C13 data
data.n = extra.n.data;   %C13 data
n0 = extra.n.n0;        %C13 orig
alpha = extra.alpha;     %fractionization constant

M = P/k + (m0 - P/k).*exp(-k*t);
N = n0*exp(-k*alpha*t) ;
%AP =  N ./ (M + N) * 100;
M = M(:);
N = N(:);
m = data.m(:);
n = data.n(:);

l1error = .95*sum(abs(M-m)) + .99*sum(abs(N - n));

%-----------------------------
function l2error = cost2(param,extra);


P = param(1);
k = param(2);
t = extra.t;
data.m = extra.m.data;
data.n = extra.n.data;
m0 = extra.m.m0;
n0 = extra.n.n0;
alpha = extra.alpha;

M = P/k + (m0 - P/k).*exp(-k*t);
N = n0*exp(-k*alpha*t) ;

M = M(:);
N = N(:);

m = data.m(:);
n = data.n(:);

error.m = abs(M-m);
error.n = abs(N-n);

l2error = .95*(error.m'*error.m)^(1/2)+ .99*(error.n'*error.n)^(1/2);
