
function [deltamin,negsumlogL,logL,G,H, checkOutput] = OPminFundTheorem(delta,y,x,beta,cri0,cri1,stepsize,maxiter,h,dh,sw)

% The function OPminFundTheorem minimizes the log-likelhood function of the
% ordinal probit model with respect to the cut-points 'delta' using the
% Fundamental Theorem of Integral Calculus.
% OUTPUT
% -------------------------------------------------------------------------
% dmin      : cutpoints that minimize the log-likelihood fucntion
% sumlogL   : sum of log-likelihood values
% logL     : logLikelihood values for each observation
% G         : gradient vector, nxk matrix. i-th row gives the score for the
%             i-th unit
% H         : hessian
% INPUT
%--------------------------------------------------------------------------
% deltam    : cutpoint starting value
% y         : dependent variable
% x         : covariate (nxk)
% beta      : a column vector of coeffcients (kx1)
% cri0      : initial criterion, cri0=1.
% cri1      : = sum(abs(sum(G))) lies between (0.001 to 0.0001), criterion to stop iteration.
% stepsize  : lies between (0.1, 1)
% maxiter   : maximum number of iteration
% h         : change in value of each delta, holding other delta's
%             constant, for first derivatives
% dh        : change in each value of delta, holding other deltas constant
%             for second derivaties. Normally dh < h, d=0.0002 and
%             dh=0.00002
% sw        : iteration to switch from BHHH to inv(-H)
%--------------------------------------------------------------------------
% First derivative from first principle
% dy/dx=[f(x+h)-f(x-h)]/2h.
%
% Second derivative from First principle
% f'(x-h)=(f(x)-f(x-h))/h
% f''(x)= [{f(x+h)-f(x)/h} - (f(x)-f(x-h)/h)]/h
%       = [f(x+h)+f(x-h)-2 f(x)]/h^2
%
% cross partial derivatives
% f(x) = [f(x+dh,y)-f(x-dh,y)]/2dh
% f(x,y)=[{f(x+dh,y+dh) - f(x+dh,y-dh)/2dh} - {f(x-dh,y+dh) -
% f(x-dh,y-dh)/2dh}]/2dh
% = 0.25* [{f(x+dh,y+dh)-f(x+dh,y-dh)} -{  f(x-dh,y+dh)-f(x-dh,y-dh)}]/dh2
%
% Code wtitten by Mohammad Arshad Rahman, Assistant Professor, IIT Kanpur. 
% Email at marshad@iitk.ac.in.


n           = length(y);        % number of observation
% k         = size(x,2);        % number of covariates
d           = length(delta);    % number of transformed cutpoints
storevn     = zeros(n,d);
storevp     = zeros(n,d);
checkOutput = zeros(n,3*d+1);

% Stuff for optimizer
cri     = cri0;
der     = zeros(n,d);       % storage matrix for derivatives
dh2     = dh^2;

% disp(' iter  logLikelihood  crit')

jj   = 0;
while (cri > cri1) && (jj < maxiter)
       jj = jj + 1;
       [nlogL, ~] = NegloglikenSum(delta,y,x,beta);
       v0        = -nlogL;
       delta0    = delta;
            
            % Calculation for first derivative
            for i = 1:d
                delta(i)  = delta(i) - h;
                [nlogL, ~] = NegloglikenSum(delta,y,x,beta);
                vn        = -nlogL;                
                delta     = delta0;
                
                storevn(:,i) = vn;
                
                delta(i)  = delta(i) + h;
                [nlogL, ~] = NegloglikenSum(delta,y,x,beta);
                vp        = -nlogL;
                delta     = delta0;
                
                storevp(:,i)   = vp; 
                
                der(:,i)  = (0.5*(vp - vn))/h;          
            end
            %--------------------------------------
            % Debugging checks
            % checkOutput = [y storevn storevp der];
            % sumder = sum(der); 
            % sumder
            % min([vp vn]);
            %[vp(1:100) vn(1:100) der(1:100,:)];
            %--------------------------------------
            % Calculation for second derivative
            hess          = zeros(d,d);
            i             = 1;
            j             = 1;
            
            while j <= d
                while i <= j
                    if i == j
                        delta(i)  = delta(i) + dh;
                        [nlogL, ~]  = NegloglikenSum(delta,y,x,beta);
                        vp2        = -nlogL;                        
                        delta     = delta0;
                        
                        delta(i)  = delta(i) - dh;
                        [nlogL, ~]  = NegloglikenSum(delta,y,x,beta);
                        vn2        = -nlogL;                        
                        delta     = delta0;
                        
                        hess(i,j) = sum( (vp2+vn2-2*v0)/dh2);
                    else
                        c         = [i,j];
                        delta(c)  = delta(c) + dh;
                        [nlogL, ~] = NegloglikenSum(delta,y,x,beta);
                        vpp       = -nlogL;
                        delta     = delta0;
                        
                        delta(c)  = delta(c) - dh;
                        [nlogL, ~] = NegloglikenSum(delta,y,x,beta);
                        vnn       = -nlogL;
                        delta     = delta0;
                        
                        delta(c)  = delta(c) + [dh; -dh];
                        [nlogL, ~]  = NegloglikenSum(delta,y,x,beta);
                        vpn       = -nlogL;
                        delta     = delta0;
                        
                        delta(c)  = delta(c) + [-dh; dh];
                        [nlogL, ~] = NegloglikenSum(delta,y,x,beta);
                        vnp       = -nlogL;
                        delta     = delta0;

                        hess(i,j) = 0.25*sum( (vpp+vnn-vpn-vnp)/dh2 );
                        % disp('calculating Hessian')
                    end
                    i = i + 1;
                end
                i = 1;
                j = j+1;
            end
            
            hess = eye(d).*hess + (1 - eye(d)).*(hess+hess');
            % This operation adds up the diagonal and off-diagonal
            % elements.
            
            cri   = sum(abs(sum(der)));
            
            % For the direction vector, we want to use the BHHH early in
            % the game and th Hessian a bit later, since the hessian can be
            % erratic in early iterations. 
            
            ddeltaBHHH  = (der'*der)\(sum(der)');  % pinv(der'*der)*(sum(der)')
            ddeltahess  = (-hess)\(sum(der)');     % pinv(hess)*(sum(der)');    
            
            % Note: Directional derivative of a function f(v) in the
            % direction u is df(v)/dv.u
            
            % nifty one-liner to switch from BHHH to hess at "sw"
            % iterations
            
            ddelta = (1-min(1,max(0,jj-sw)))*ddeltaBHHH + min(1,max(0,jj-sw))*ddeltahess;
            
            % Next compute update using the step-size as coefficient on
            % ddelta
            delta = delta + stepsize*ddelta;
            
            % iter = [jj sum(nlogL) cri];
            % disp('iter')
            
            if jj==maxiter
                disp('Maximum number of iterations reached')
            end
end

deltamin   = delta;
[negloglike, negsumlogL] = NegloglikenSum(deltamin,y,x,beta);
logL       = - negloglike;
G          = der;
H          = hess;   
                    
end
       

