% eda11_06
% stochastic gradient method, full error every
% 100 iterations
% non-linear least squares fit of a sinusoid
% of unknown frequency, using the stochastic
% gradient method.  The data are subsampled on
% each iteration.  A fixed step size alpha is
% used to adcance the solution anti-parallel to the
% gradient. The solution is considered good enough
% when its posterior variance is reduced to a
% specififed value.

clear all;

% load Black Rock Forest temperature data
D=load('brf_temp.txt');
t=D(:,1);
dobs=D(:,2);
Ns=size(D);

% clean the data by deleting outliers and zero (= no data) values
k = find( (dobs>-35) & (dobs<40) & (dobs~=0.0) );
t = D(k,1);
dobs = D(k,2);
N = length(dobs);

% subsampling information
decimation_factor = 1000;
Ns = floor(N/decimation_factor);

% some statistics
dmean = mean(dobs);
amp = 0.5*(max(dobs)-min(dobs));

% the solution considered good enough when
% its posterior variance is smaller than
low_enough_variance = 35;

% Annual period, frequency, angular frequency
Ta = 365;
fa = 1/Ta;
wa = 2*pi*fa;

% the model parameters m in this problem have been
% normalized so that we expect them to be of order
% 1.  The DC offset is parameterized as a fraction
% of the mean of the data, dmean*m(1). The amplitudes
% of the sine and cosine are parameterized as a fraction
% of the overall amplitude of the data, m(2)*amp and
%  m(3)*amp. The frequency is parameterized as a
% fraction of the annual frequecy, m(4)*wa.

% The resulting nonlinear theory and its derivative
% d = m1*dmean + m2*amp*sin(m4 wa t) + m3*amp*cos(m4 wa t)
% dd/dm1 =  dmean
% dd/dm2 =  amp*sin(m4 wa t)
% dd/dm3 =  amp*cos(m4 wa t)
% dd/dm4 =  amp*m2*wa*t*cos(m4 wa t) - amp*m3*wa*t*sin(m4 wa t)

% inital guesses
M=4;
m = zeros(M,1);
m(1) =  1;          % average level is close to the mean
m(2) =  0;          % phase with max in summer, so no sine term
m(3) = -0.8;        % cosine term is negative (hot in summer) and
                    % probably less than 1 sine data contains noise
m(4) =  1;          % frequency close to annual frequency

% evaluate initial guess
s = amp*sin(m(4)*wa*t);
c = amp*cos(m(4)*wa*t);
dpre = dmean*m(1) + m(2)*s + m(3)*c;

% plot observations
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis([0,5000,-30,50]);
plot(t,dobs,'k.','LineWidth',2);
hold on;
% plot(t,dpre,'-','LineWidth',2,'Color',[0.7,0.7,0.7]);
xlabel('time (days)');
ylabel('temperature (deg C)');

% trial solution
mg = m;

% the step size alpha is fixed
alpha = 1e-4;

% number of iterations is expected range of model
% parameters divided by step size
range = 1;
Niter = floor(range/alpha);

for k = [1:Niter]
    
    % subsample data
    ns = unidrnd(N,Ns,1);
    ts = t(ns);
    dobss = dobs(ns);
    
    % compute predicted data, the error and the
    % gradient of the error of the subsampled data
    s = amp*sin(mg(4)*wa*ts);
    c = amp*cos(mg(4)*wa*ts);
    dpreg = dmean*mg(1) + mg(2)*s + mg(3)*c;
    eg = dobss - dpreg;
    Eg = eg'*eg;
    if( k==1 )
        Estart = Eg;
    end
    G = [ dmean*ones(Ns,1), s, c, m(2)*wa*ts.*c-m(3)*wa*ts.*s ];
    dEdm = -2*G'*eg;
    % posterior variance of data
    sigma2d =  Eg/(Ns-M);
    
    % compute just the error for the full dataset
    % but only every 100 iterations
    if( mod(k-1, 100) == 0 )
        sfull = amp*sin(mg(4)*wa*t);
        cfull = amp*cos(mg(4)*wa*t);
        dpregfull = dmean*mg(1) + mg(2)*sfull + mg(3)*cfull;
        egfull = dobs - dpregfull;
        Egfull = egfull'*egfull;
        % posterior variance of data
        sigma2dfull =  Egfull/(N-M);
        if( k==1 )
            Estartfull = Egfull;
        end
        % stop when the veriance is low enough
        if( sigma2dfull <= low_enough_variance )
            break;
        end
    end
    
    % downslope unit vector
    v = -dEdm / sqrt(dEdm'*dEdm);
    % update solution
    mg = mg+alpha*v;

end

m = mg;
s = amp*sin(m(4)*wa*t);
c = amp*cos(m(4)*wa*t);
dpre = dmean*m(1) + m(2)*s + m(3)*c;

plot(t,dpre,'-','LineWidth',2,'Color',[0.9,0.9,0.9]);

ffinal = wa*m(4)/(2*pi);
Tfinal = (2*pi)/(wa*m(4));

fprintf('Iterations: %d\n', k );
fprintf('mean: %f amp: %f %f\n', m(1)*dmean, m(2)*amp, m(3)*amp );
fprintf('Estimated period %.3f error ratio %.3f\n', Tfinal, Egfull/Estartfull );

