example_num_opt_SHS: Numerical optimization example for mean

function [estlambda, lambda] = example_num_opt_SHS
%
% Numerical optimization example for mean-field model of SHS
%

%% For testing, generate some random trajectories

% Sample once every dt seconds
dt = 5;
t = 0:dt:120;

% Noise scale
noisesd = 0.01;

% Arbitrary rates that we hopefully can fit later
% (scale to accomodate sampling interval)
lambda = [ 0 2 3 ; 5 0 9 ; 2 7 0 ];
lambda = lambda/max(max(lambda)) * 1/(2*dt);

% Graph laplacian (to generate fake data)
laplacian = lambda.*(1-eye(size(lambda))) - (1-eye(size(lambda)))*lambda.*eye(size(lambda));

% Fake data, with random initial condition and noise added at every step

% Initial condition for mean-field trajectory is randomized
meandata = zeros(length(lambda), length(t));
meandata(:, 1) = randperm(length(lambda))' + noisesd*randn(length(lambda),1);

% Generate mean-field trajectory based on Laplacian and some noise
for i = 2:length(t)
    meandata(:, i) = expm(laplacian*t(i))*meandata(:, 1) + noisesd*randn(length(lambda),1);
end

%% Use a least-squares fit to get lambda back out from meandata

% fmincon (non-linear optimization with constraints) arguments:
%    1. Function handle defining the metric to minimize
%    2. Initial condition (all of the off-diagonal elements of lambda matrix)
%    3. Matrix "A" in the inequality constraint A*x <= b
%    4. Vector "b" in the inequality constraint A*x <= b
% In the linear inequality constraint, A*x <= b, "x" is a feasible solution
% to the optimization problem.

% The output vector is the optimal off-diagonal elements of the rate matrix
estlv = fmincon( @(lv)( lsmetric(lv, meandata, t) ), ...
            randperm(length(lambda)*(length(lambda) - 1)), ...
            [ -eye(length(lambda)*(length(lambda) - 1)) ; ...
              eye(length(lambda)*(length(lambda) - 1)) ], ...
            [ zeros(length(lambda)*(length(lambda) - 1), 1) ;
              (1/(2*5))*ones(length(lambda)*(length(lambda) - 1), 1) ] );

% These lines turn the estlv vector into an estimated rate matrix estlambda
estlambda = zeros( length(lambda) );
offdiagonals = setdiff(1:(length(lambda))^2, ...
                       (1:length(lambda))+length(lambda)*((1:length(lambda))-1));
estlambda(offdiagonals) = estlv;

end

%% Metric for least-squares fit

% This metric sums up the squared differences between the data (i.e.,
% meandata at times t) and a trajectory defined by "lambdavec"

function lssum = lsmetric( lambdavec, meandata, t )
    % This maps the number of elements in lambdavec to the number of states
    num_states = max( roots( [1 -1 -length(lambdavec)] ) );

    % This generates a rate matrix whose diagonal entries are zero and all
    % off-diagonal entires come from lambdavec
    lambda = zeros( num_states );
    offdiagonals = setdiff(1:(num_states^2), ...
                           (1:num_states)+num_states*((1:num_states)-1));
    lambda(offdiagonals) = lambdavec;

    % This generates the Laplacian from the rate matrix
    laplacian = lambda.*(1-eye(size(lambda))) - (1-eye(size(lambda)))*lambda.*eye(size(lambda));

    % This sums up the squared differences over time
    lssum = 0;
    for i = 2:length(t)
        lssum = lssum + norm( expm(laplacian*t(i))*meandata(:,1) - meandata(:,i) )^2;
    end

end