function sensitivity
etohmodstruct;
bestfitparam;

k(1) = 1.769292045467269;
k(2) = 20.589996419308118;

%interpolate as chebfun for speedup
prior2 = chebfun(@(x)nctpdf(x,k(1),k(2)),domain(0,2e3));


param={'mu1','k1','alpha','g1','k2','g2','sigma1','sigma2'};

fitDataE = [ dE(:,1); dE(:,2) ].';
fitSTDE  = [ stdE(:,1); stdE(:,2) ].';

mini                        = 1e-5*ones(1,length(p));
maxi                        = 1e5*ones(1,length(p));

% Include priors?
%mStruct.p.mu1=1;
%mStruct.p.k1=2;
%mStruct.p.alpha=3;
%mStruct.p.g1=4;
%mStruct.p.Kg=5;
%mStruct.p.k2=6;
%mStruct.p.g2=7;
priors = 0;

fun = @(par)residETOHPROM( t, x0, par, fitDataE, fitSTDE, input, prior2 );


initPars                    = p;

opts = optimset( 'TolX', 1e-12, 'TolFun', 1e-6, 'MaxFunEvals', 1e4, 'MaxIter', 1e4, 'Jacobian', 'off', 'Display', 'off' );


[initPars, RN] = lsqnonlin( fun, initPars, mini, maxi, opts );
[initPars, RN] = lsqnonlin( fun, initPars, mini, maxi, opts );
[initPars, RN] = lsqnonlin( fun, initPars, mini, maxi, opts );
[initPars, RN] = lsqnonlin( fun, initPars, mini, maxi, opts );
[initPars, RN] = lsqnonlin( fun, initPars, mini, maxi, opts );
[initPars, RN] = lsqnonlin( fun, initPars, mini, maxi, opts );

par = initPars(1:6);


x0(1)       = par(3)/par(4);
x0(2) 		= par(5)/par(6)*x0(1);
x0(3)       = par(3)/par(4);
x0(4)       = par(5)/par(6)*x0(3);

[tp,dpdx] = ode23s(@(t,x)dxdp(t,x,par),[0 24],[x0';zeros(24,1)]);

for tx = 1:length(tp)
    DPDX = [dpdx(tx,5), dpdx(tx,9),  dpdx(tx,13), dpdx(tx,17), dpdx(tx,21), dpdx(tx,25);
            dpdx(tx,6), dpdx(tx,10), dpdx(tx,14), dpdx(tx,18), dpdx(tx,22), dpdx(tx,26);
            dpdx(tx,7), dpdx(tx,11), dpdx(tx,15), dpdx(tx,19), dpdx(tx,23), dpdx(tx,27);
            dpdx(tx,8), dpdx(tx,12), dpdx(tx,16), dpdx(tx,20), dpdx(tx,24), dpdx(tx,28)];
    NDPDX = DPDX.*(dpdx(tx,1:4)'*(1./[par]));
    abssens(tx,1:6) = DPDX(1,:);
    abssens(tx,7:12) = DPDX(2,:);
    abssens(tx,13:18) = DPDX(3,:);
    abssens(tx,19:24) = DPDX(4,:);
    relsens(tx,1:6) = NDPDX(1,:);
    relsens(tx,7:12) = NDPDX(2,:);
    relsens(tx,13:18) = NDPDX(3,:);
    relsens(tx,19:24) = NDPDX(4,:);
end

%mStruct.p.mu1=1;
%mStruct.p.k1=2;
%mStruct.p.alpha=3;
%mStruct.p.g1=4;
%mStruct.p.k2=5;
%mStruct.p.g2=6;

for m=[1,3]
    switch m
        case{1}
            disp = 'Induced';
        case{3}
            disp = 'Uninduced';
    end
    %subplot(2,4,m+1)
    figure(m)
    clf
    set(gca,'fontsize',16);
    set(gcf,'Position',[1 1 600 300])
    semilogy(tp,abs(abssens(:,1+m*6:6+m*6)),'LineWidth',3)
    title([disp ': Absolute Sensitivity'])
    legend('mu1','k1','alpha','g1','k2','g2','Location','NorthEastOutside')
    xlim([0,24])
    xlabel('Time [hours]')
    %subplot(2,4,m+5)
    figure(1+m)
    clf
    set(gca,'fontsize',16);
    set(gcf,'Position',[1 1 600 300])
    semilogy(tp,abs(relsens(:,1+m*6:6+m*6)),'LineWidth',3)
    title([disp ': Relative Sensitivity'])
    legend('mu1','k1','alpha','g1','k2','g2','Location','NorthEastOutside')
    xlim([0,24])
    xlabel('Time [hours]')
end

end

function [dxdp] = dxdp(t,x,p)

%M = reshape(x(5:16),4,3);
M = [x(5), x(9), x(13) x(17) x(21) x(25);
    x(6), x(10), x(14) x(18) x(22) x(26);
    x(7), x(11), x(15) x(19) x(23) x(27);
    x(8), x(12), x(16) x(20) x(24) x(28)];

H = [ 0                              1       1      -x(1)    0       0;
      p(5) *  x(1) - p(6) * x(2)     0       0       0       p(1)*x(1)   -p(1)*x(2);
      0                              0       1      -x(3)    0       0;
      p(5) *  x(3) - p(6) * x(4)     0       0       0      -p(1)*x(3)   -p(1)*x(2)];
    
S = [-p(4)  0       0       0;
      p(1)*p(5) -p(1)*p(6)    0       0;
      0     0      -p(4)    0;
      0     0       p(1)*p(5)   -p(1)*p(6)];

R = S*M + H;


dxdp = [
    p(3) + p(2) - p(4) * x(1);
    p(1)*(p(5) *  x(1) - p(6) * x(2));
    p(3) - p(4) * x(3);
    p(1)*(p(5) *  x(3) - p(6) * x(4));
    R(:,1);
    R(:,2);
    R(:,3);
    R(:,4);
    R(:,5);
    R(:,6);
    ];


end
