% function result = one_stage_regression(ip,ih,iitr)
% ip - input parameters
%   ip.oid - origin id
%   ip.gm - vaue of ground motion
%   ip.m - magnitude
%   ip.dist - distance from station to origin [km]
% ih - initail h parameter
% iitr - number of minimization iterations
% WARNING the data should be sorted
function result = one_stage_regression(X,Y,ic,nev,iitr,x0)
    N = length(X) ;
    v = eye(N) ;
    options = optimset('MaxFunEvals', iitr) ;
    % Initial B parameters
%     fprintf('X: %d,%d; Y: %d,%d; v: %d,%d\n',size(X),size(Y),size(v))
%     fprintf('v\\X: %d,%d\n',size(v\X))
%     fprintf('(XT * (v d X)): %d,%d\n',size( X' * (v \ X)))
%     fprintf('(XT * (v d X)) d XT: %d,%d\n',size( (X' * (v \ X)) \ X' ))
%     fprintf('v d Y: %d,%d\n',size( (v \ Y) ))
%     b1 = (X' * (v \ X)) \ X'; 
%     b2 = v \ Y;
%     fprintf('b1: %d,%d; b2: %d,%d\n',size(b1),size(b2))
%     sum(isnan(b1))
%     sum(isnan(b2))
%     B = b1 * b2
    % Minimisation B and epsilon 
    if nargin < 6
        x0 = rand() * 0.29 + 0.01 ;
    end
    x = fminsearch(@(x)ssr_maxlikelihood(x,X,Y,ic,nev,N), x0, options) ;
    % Result computation
    ipsylon = x ;
    v = v_matrixssr_maxlikelihood(ipsylon,ic,nev,N) ;
    B = (X' * (v \ X)) \ X' * (v \ Y) ;
    tm = (Y - X * B) ;
    sigma_2 = tm' * (v \ tm) / N ;
    unbiasedsigma_2 = sigma_2 * N / (N - 1 - size(B,1));
    % result setting
    result.x = x ;
    result.b = B ;
    result.ipsylon = ipsylon ;
    result.sigma = sqrt(sigma_2) ;
    result.sigma_r = sqrt((1-ipsylon)*unbiasedsigma_2) ;
    result.sigma_e = sqrt(ipsylon*unbiasedsigma_2) ;
    result.unbiased_sigma = sqrt(sigma_2) ;
    result.S2 = sqrt(tm' * tm / N) ;
    result.SSE = sum(tm .^ 2);
    result.N = N ;
    result.e = tm ;
    fprintf('std %e\n',std(tm))
end

function min_lnL = ssr_maxlikelihood(x,X,Y,ic,nev,N)
    v = v_matrixssr_maxlikelihood(x,ic,nev,N) ;
    size(v);
    B = (X' * (v \ X)) \ X' * (v \ Y);
    tm = (Y - X*B) ;
    sigma_2 = tm' * (v \ tm) / N;
    det_v = det(v);
%    min_lnL = N/2*log(sigma_2) + 0.5 * log(abs(det_v));
    min_lnL = sigma_2 *N + 0.5 * log(abs(det_v));
    fprintf('ipsylon=%f(det=%e), minlnL=%f, sigma^2=%f, \n',x,det_v,min_lnL,sigma_2)

end

function v = v_matrixssr_maxlikelihood(ipsylon,ic,nev,N)
    v = zeros(N,N) ;
    for i = 1:nev
        indexes = (ic == i) ;
        v(indexes,indexes) = ipsylon ;
    end
    v = v + eye(N,N)*(1-ipsylon) ;    
end
