InitialConditions.m

Go to the documentation of this file.

namespace models{
namespace motorunit{
namespace experiments{


/* (Autoinserted by mtoc++)
 * This source code has been filtered by the mtoc++ executable,
 * which generates code that can be processed by the doxygen documentation tool.
 *
 * On the other hand, it can neither be interpreted by MATLAB, nor can it be compiled with a C++ compiler.
 * Except for the comments, the function bodies of your M-file functions are untouched.
 * Consequently, the FILTER_SOURCE_FILES doxygen switch (default in our Doxyfile.template) will produce
 * attached source files that are highly readable by humans.
 *
 * Additionally, links in the doxygen generated documentation to the source code of functions and class members refer to
 * the correct locations in the source code browser.
 * However, the line numbers most likely do not correspond to the line numbers in the original MATLAB source files.
 */

function  InitialConditions() {



/*  Create with no dynamic initial conditions as this script is intended to
 * compute them :-) */
file = fullfile(fileparts(mfilename(" fullpath "))," ic_singletwitch.mat ");
filex0 = fullfile(fileparts(mfilename(" fullpath "))," .. "," x0coeff ");
if exist(file," file ") == 2
    load(file);
else
    m = models.motorunit.Shorten(false,true);
    m.UseNoise= false;
    s = sampling.ManualSampler;
    p = linspace(0,1,200);
    p = [p; ones(size(p))*1];
    s.Samples= p;

    m.off1_createParamSamples;
    m.T= 5000;
    m.dt= .1;

    n = size(p,2);
    maxvals = zeros(1,n);
    maxidx = [];
    endvals = zeros(m.System.f.xDim,n);
    ctimes = zeros(1,n);
    pi = ProcessIndicator(" Running %d simulations with one twitch ",n/(max(1,matlabpool(" size "))),false,n);
    /* for k=1:n */
    parfor k=1:n
        [t, y, ctimes(k), x] = m.simulate(p(:,k),1);/* #ok */

        endvals(:,k) = x(:,end);
        pi.step;/* #ok */

    end
    /*  Collect all trajectories into the local/main model */
    m.Data.TrajectoryData.consolidate(m);
    pi.stop;
    m.save;
    save(file, " m ", " endvals ", " p ", " ctimes ");
end

/* % Fit every dimension by polynomial */
deg = 8;
dim = size(endvals,1);
coeff = zeros(dim,deg+1);
err = zeros(dim);
a = dscomponents.AffineInitialValue;
pi = ProcessIndicator(" Fitting polynomials ",dim);
for k=1:dim
    coeff(k,:) = polyfit(p(1,:),endvals(k,:),deg);
    err(k) = Norm.Linf((endvals(k,:) - polyval(coeff(k,:),p(1,:)))^t);
    a.addMatrix(sprintf(" polyval([%s],mu(1)) ",sprintf(" %g  ",coeff(k,:))),full(sparse(k,1,1,dim,1)));
    pi.step;
end
pi.stop;

/* % Error checks */
endvals_afflin = zeros(size(endvals));
for k=1:n
    endvals_afflin(:,k) = a.evaluate(p(1,k));
end
semilogy(p(1,:),Norm.L2(endvals_afflin-endvals)./Norm.L2(endvals));
title(sprintf(" Relative errors of initial values for %d increasing mu_1 ",n));

save(file, " a ", " endvals_afflin ", " err ", " coeff ", " -APPEND ");
save(filex0, " coeff ");
}
};
};
};