% main2.m Fit of the pigs data
figure(1)
% read data  matrix
ty; N = max(size(Y)); T = T'; Y = Y';
% plot given points
ax = ppoints([T,Y]);
% initial partition
xi = 120;   n1 = sum(T<xi);
% initialize  coefficients
%  a = [-7 9 0.01]';  b = [133 -550 -0.01]';
a = [-1  5   0.03]'; b = [1 -200 -0.02]';

% Levenberg Marquart regularization (1=yes)
levenberg = 1;  ee = 1e-2;
% plot initial approximation functions
pexps(xi,a,b,ax);
% initialize loop counter
k = 0;
% initialise corrections
dp = 1;

while (norm(dp) > 1e-5) & (k<40),
   k= k+1;
   Jf = [ones(n1,1) exp(a(3)*T(1:n1))...
           a(2)*T(1:n1).*exp(a(3)*T(1:n1))];
   Jg = [ones(N-n1,1) exp(b(3)*T(n1+1:N)) ...
           b(2)*T(n1+1:N).*exp(b(3)*T(n1+1:N))];

   f = a(1) + a(2)*exp(a(3)*T(1:n1));
   g = b(1) + b(2)*exp(b(3)*T(n1+1:N));

   JF = [dirsum(Jf,Jg), zeros(N,1)];

   % for Levenberg-Marquart add
   if levenberg,  JF = [JF;ee*eye(7)]; end

   gradfxi = [1 exp(a(3)*xi) a(2)*xi*exp(a(3)*xi)];
   gradgxi = [1 exp(b(3)*xi) b(2)*xi*exp(b(3)*xi)];

   gradfsxi = [0 a(3)*exp(a(3)*xi) ...
                 a(2)*exp(a(3)*xi)+a(2)*a(3)*xi*exp(a(3)*xi)];
   gradgsxi = [0 b(3)*exp(b(3)*xi)...
                 b(2)*exp(b(3)*xi)+b(2)*b(3)*xi*exp(b(3)*xi)];

   JG = [gradfxi -gradgxi ...
           a(2)*a(3)*exp(a(3)*xi)-b(2)*b(3)*exp(b(3)*xi)
         gradfsxi -gradgsxi ...
           a(2)*a(3)^2*exp(a(3)*xi)-b(2)*b(3)^2*exp(b(3)*xi)];

   mG = -[a(1)+a(2)*exp(a(3)*xi)-b(1)-b(2)*exp(b(3)*xi)
          a(2)*a(3)*exp(a(3)*xi)-b(2)*b(3)*exp(b(3)*xi)];

   % solve the constrained least squares problems for corrections
   if ~levenberg,
      dp = lsqlce(JF,JG,Y-[f;g],mG);
   else   % for Levenberg-Marquart use
      dp = lsqlce(JF,JG,[Y-[f;g];zeros(7,1)],mG);
   end

   %  update unknowns
   a = a + dp(1:3), b = b + dp(4:6), xi = xi + dp(7)
   % print norm of correction vector
   normcorrection = norm(dp)
   % plot currrent approximation functions
   pexps(xi,a,b,ax);
   % new partition of the data
   n1 = sum(T<xi)
   pause(0.5)
end

% plot and print final results
clf; ax = ppoints([T,Y]); pexps(xi,a,b,ax);
totalsteps = k, a, b, xi

