function [logmax,dL]=nolin_ssts_mxd_fast_v2(x,yz,M1,expA2,ord,ff,Tz,g1)

% NOLIN_SSTS_MXD_FAST_V2    This optimisation script uses only Gradient 
% information to perform Gaussian Process prior model optimisation.
% Hessian information is supplied by user. The hyperparamters, g represents
% the length-scale, and V, the noise variance.
%
% [logmax,dL] = NOLIN_SSTS_MXD_FAST_V2(x,y,M1,expA2,ord,ff,Tz,g1)
%
% Inputs:
%           x       : Noise hyperparameter for GP, x = [a1 V]
%           y       : One-dimensional outcome column vector
%           M1      : Time-series input explantory variables
%           expA2   : Cov. function Qg = a2*exp{ -0.5*g2*(ZZi-ZZj)^2 }
%           ord     : Order of data
%           ff      : Fudge factor, e
%           Tz      : Reduced "y"
%           g1      : Lengthscale hyperparameter, gts
%
% Outputs:
%           logmax  = Log likelihood function value
%           dL      = Gradient vector value
%
% Dependencies: 
%   - get_eigen_block
%   - get_eigen_block_traceqq_a1
%   - get_eigen_block_traceqdq_a
%   - schur_invssts_c (C)
%   - schur_invssts_qq_c (C)
%   - toep2_genct (C)
%   - toep2_logdet (C)
%   - toep2_invcy (C)
%   - toep2_trace (C) {optional, for testing}
%   - toe_times_y (C)
%
%   **(C): defines C codes
%
%
% (C) Gaussian Process SSTS Toolbox 1.0
% (C) Copyright 2006-2007, Keith Neo <keith.neo@nuim.ie>
%                          http://www.hamilton.ie

a1=exp(x(1));
V=exp(x(2));
Vxp=V*(1+ff);

G=get_eigen_block(M1,[a1 g1 Vxp]);    % Deploy Schur algorithm here
% [invQ,logdetQftemp,trQf,invFz]=schur_invssts(G,ord,yz);     % Matlab code
[invQ,logdetQf,trQf,invFz]=schur_invssts_c(G,ord,yz);   % C code
% disp(sum(sum(abs(invQ)-abs(invQc))));
% disp(logdetQftempc/logdetQftemp-1);
% disp(trQf/trQfc-1);
% disp(norm(invFz-invFzc)/norm(invFz));
clear G

QG=diag_add(expA2,Vxp,length(Tz))-(V*V)*invQ;
MatGF=inv(QG);
Y=MatGF*(Tz-V*getord(ord,invFz));

ct=toep2_genct(M1,[a1 g1 Vxp]);       % Deploy Durbin-Levinson algorithm here
[logdetQf,ds,vn1]=toep2_logdet(ct);     % For use in derivative LL
% dt=[1 zeros(1,N1-1)];
% trQftemp=toep2_trace(ds,vn1,dt);
% disp(logdetQf/logdetQftemp-1);  %% Test for correctness of logdet(invQf)
% disp(trQf/trQftemp-1);          %% Test for correctness of trace(invQf)
X=invord(ord,Y,length(yz));
X=toep2_invcy(ct,X);
X=invFz-V*X;
invQy=[X;Y];

[cholC,pp]=chol(QG); % Calculation of determinant of state-space matrix
if pp==0, logdetQg=2*sum(log(diag(cholC)));
else logdetQg=sum(log(svd(QG)));
end
clear cholC

logdet=logdetQf+logdetQg;
logmax=logdet+[yz;Tz]'*invQy;


%% Computing the first order derivatives of log-likelihood function
ydQy3=(1+ff)*(X'*X+Y'*Y)+2*(Y'*getord(ord,X));
ct(1)=a1;       % variable for ct has changed
ydQy1=X'*toe_times_y(ct,X);
dt=-0.5*g1*(M1'-M1(1)).^2;
dt=dt.*exp(dt);     % Note: dt is less factor 'a1'
ydQy2=a1*(X'*toe_times_y(dt,X));

trQg=trace(MatGF);
trNQN=traceAxxB(invQ,MatGF);
trdQda=toep2_trace(ds,vn1,ct);
trdQdg=a1*toep2_trace(ds,vn1,dt);

coef=Vxp/a1;
G1=get_eigen_block_traceqq_a1(M1,[g1 coef],coef);   % Set fudge factor = coef
% invQQ=schur_invssts_qq(G1,ord);       % Matlab code
invQQ=schur_invssts_qq_c(G1,ord);       % C code
% disp(sum(sum(abs(invQQ)-abs(invQQc))));
invQQ=invQQ/(Vxp*a1)-invQ/Vxp;
trNQQN=traceAxxB(invQQ,MatGF);

G1=get_eigen_block_traceqdq_a(M1,[g1 coef],1);      % reformulated for lower rank
invQQ=schur_invssts_qq_c(G1,ord);   % re-using C code from invQ*invQ
invQQ=invQQ/a1-invQ;
trNQdPQN=traceAxxB(invQQ,MatGF);

dL=zeros(2,1);
dL(1)=trdQda+V*V*(trNQN-Vxp*trNQQN)-ydQy1;
% dL(2)=trdQdg+V*V*trNQdPQN-ydQy2;
dL(2)=V*((1+ff)*(trQf+trQg+V*V*trNQQN)-2*V*trNQN-ydQy3);



function z=getord(x,y)
if x==0
    z=y;
else
    x1=x+1;
    D0=ceil(size(y,1)/x1);
    z=zeros(D0,size(y,2));
    for k=1:D0
        pp=x1*k-x;
        z(k,:)=y(pp,:);
    end
end


function z=invord(x,y,N1)
Y0=size(y);
if x==0 && Y0(1)==N1
    z=y;
else
    yx1=Y0(1)*(x+1);
    z=zeros(N1,Y0(2));
    for k=1:Y0(2)
        tmpZ=[y(:,k)';zeros(x,Y0(1))];
        tmpZ=reshape(tmpZ,yx1,1);
        z(:,k)=tmpZ(1:N1);
    end
end


function Q=diag_add(Q,v,N0)
for pp=1:N0
    Q(pp,pp)=Q(pp,pp)+v;
end


function trace_AB=traceAxxB(A,B)
trace_AB=0;
for i=1:size(A,1)
    trace_AB=trace_AB+A(i,:)*B(:,i);
end