text:
function mvnx = load_mvnx(filename)
% mvnx = load_mvnx(filename)
% loads a mvnx file
% filename is name of file (including path)
% mvnx is result struct containing all data of the mvnx file
% Xsens Technologies BV 28-05-2015

%% check filename
if isempty(strfind(filename,'mvnx'))
    filename = [filename '.mvnx'];
end
if ~exist(filename,'file')
    error([mfilename ':xsens:filename'],['No file with filename: ' filename ', file is not present or file has wrong format (function only reads .mvnx)'])
end
%% read data
fid = fopen(filename, 'r', 'n', 'UTF-8');
content = char(fread(fid,'char')');
fclose(fid);
%% get data into a cell array, one cell per line
lines = find(content == 10);
cellContent = cell(1,length(lines));

hWaitbar = waitbar(0, 'MVNX import: importing data...');
for n=1:length(lines)
    if n==1
        cellContent{n} = content(1:lines(n)-1);
    else
        cellContent{n} = content(lines(n-1)+1:lines(n)-1);
    end
end

% set file comments on one line
beginComment = find(cellfun(@(x) ~isempty(strfind(x, '<comment>')), cellContent));
endComment = find(cellfun(@(x) ~isempty(strfind(x, '</comment>')), cellContent));

if ~all(beginComment == endComment)
    for i=1:length(beginComment)
        cellContent{beginComment(i)} = [sprintf('%s ',cellContent{beginComment(i):endComment(i)-1}) cellContent{endComment(i)}];
        [cellContent{beginComment(i)+1:endComment(i)}] = deal('<!>');
    end
end
% look for comment lines and remove them
cellContent(cellfun(@(x) x(1)=='-' || (x(1)=='<' && (x(2)=='?' || x(2)=='!')),cellContent)) = [];
% get start of text and clean up some
index = cellfun(@(x) find(x==60,1),cellContent); % lines with only comments on them won't have this mark
openandclose = cellfun(@(x) sum(x==60)==2,cellContent); % lines with have an opening and closing statement in one line
cellContent = cellfun(@(x, y) x(y:end), cellContent, num2cell(index),'UniformOutput',false);

%% get start and end words
n = 0;clear value name
iWord = 0;
word = cell(1,length(cellContent));
wordindex = cell(1,length(cellContent));
wordvalue = cell(1,length(cellContent));
wordfields = cell(1,length(cellContent));
while n < length(cellContent)
    n=n+1;
    line = cellContent{n};oneword = false;
    hooks = find(line == 62);
    iWord = iWord+1;
    wordindex{iWord} = index(n);
    if any(line == 32) && ~openandclose(n)
        if ~isempty(hooks) && hooks(1) < find(line==32,1)
            word{iWord} = line(2:hooks(1)-1);
            iLine = hooks(1)+1;
        else
            word{iWord} = line(2:find(line==32,1)-1);
            iLine = find(line==32,1)+1;
        end
    elseif ~isempty(hooks) && length(line)>=8 && strcmp('comment',line(2:8))
        % add exception for comment
        word{iWord} = line(2:hooks(1)-1);
        iLine = hooks(1)+1;
    elseif openandclose(n)
        word{iWord} = line(2:find(line==62,1)-1);
        iLine = find(line==62,1)+1;
    else
        word{iWord} = line(2:end-1);
        oneword = true;
    end
    if word{iWord}(1) ~= '/'
        if ~oneword && ~openandclose(n)
            k = find(line == 34);
            k = reshape(k,2,length(k)/2)';
            l = [iLine find(line(iLine:end) == 61)+iLine-2];
            fieldname = cell(1,length(l)-1); value = cell(1,length(l)-1);
            if ~isempty(k)
                for il=1:size(k,1)
                    fieldname{il} = line(iLine:find(line(iLine:end) == 61,1)+iLine-2);
                    if size(k,1) > 1 && il < size(k,1)
                        a = strfind(line(iLine:end),'" ')+iLine+1;
                        iLine = a(1);
                    end
                    value{il} = line(k(il,1)+1:k(il,2)-1);
                end
            else
                value = []; fieldname =[];
                value = line(find(line == 62,1)+1:end);
            end
        elseif ~oneword && openandclose(n)
            value = []; fieldname =[];
            value = line(find(line == 62,1)+1:find(line==60,1,'last')-1);
        else
            value = NaN;fieldname = [];
        end
        wordvalue{iWord} = value;
        wordfields{iWord} = fieldname;
    end
end
isendword = cellfun(@(x) x(1) == '/',word);
endwords = cellfun(@(x) x(2:end),word(isendword),'UniformOutput',false);
kindofendwords = unique(endwords);
placeoffirststartword = zeros(1,length(kindofendwords));
placeofallendwords = cell(1,length(kindofendwords));
placeofallstartwords = cell(1,length(kindofendwords));
for n=1:length(kindofendwords)
    lengthWord = length(cell2mat(kindofendwords(n)));
    placeoffirststartword(n) = find(strncmp(word,kindofendwords(n),lengthWord),1);
    placeofallstartwords{n} = find(strncmp(word,kindofendwords(n),lengthWord));
    placeofallendwords{n} = find(strncmp(word,['/' kindofendwords{n}],lengthWord+1));
end
[a b] = sort(placeoffirststartword);
startwords = kindofendwords;
for n=1:length(startwords)
    obj.(startwords{n}).number = length(placeofallstartwords{n});
    obj.(startwords{n}).index = index(a(n));
    obj.(startwords{n}).count = 0;
end

%% get values
for n=1:length(wordvalue)
    if mod(n,5000) == 0
        waitbar((n/length(wordvalue))/2, hWaitbar);
    end

    if iscell(wordvalue{n})
        if length(wordvalue{n}) == 1
            B = [];
            try
                B = str2num(wordvalue{n}{1});
            end
            if ~isempty(B)
                wordvalue{n} = B;
            else
                wordvalue{n} = wordvalue{n}{1};
            end
        else
            for m=1:length(wordvalue{n})
                try
                    B = str2num(wordvalue{n}{m});
                    if ~isempty(B)
                        wordvalue{n}{m} = B;
                    end
                end
            end
        end
    else
        try
            B = str2num(wordvalue{n});
            if ~isempty(B)
                wordvalue{n} = B;
            end
        end
    end
end
%% put everything in struct
firstFramesFound = false;
nFramesToBeFound = 5;
nFramesFound = 0;
nFrames = length(cell2mat(placeofallendwords(strcmp(kindofendwords,'frame'))));
lastwordisendword = false;lastendword = '';
superstruct = struct;name = cell(1,max(index)+1);namecounter = ones(1,max(index));endplace = inf(1,max(index));
for iWord=1:length(word)
    if mod(iWord,5000) == 0
        waitbar(iWord/length(word) + 0.5, hWaitbar);
    end
    if word{iWord}(1)=='/'
        startword = word{iWord}(2:end);
        if any(strcmp(startword,name))
            name(find(strcmp(startword,name),1):end) = [];
        end
        if lastwordisendword
            obj.(lastendword).count = 0;
        end
        lastendword = startword;
        lastwordisendword = true;
    else
        lastwordisendword = false;
        try
            if ~isempty(wordfields{iWord})
                wordfields{iWord} = cellfun(@(x) x(x~=32), wordfields{iWord},'UniformOutput',false);
            end
            word(iWord) = checkText(word(iWord));
            wordfields{iWord} = checkText(wordfields{iWord});

            if any(strncmp(startwords,word{iWord},length(startwords)))
                if strcmp(word{iWord},'frame')
                    nFramesFound = nFramesFound + 1;
                    if ~firstFramesFound && nFramesFound >= nFramesToBeFound
                        firstFramesFound = true;
                        % preallocate frames in superstruct
                        % Assume remainder of the file contains frames,
                        % with exception of security code at the end
                        frames(1:nFrames,1) = superstruct.mvnx.subject.frames.frame(end);
                        frames(1:nFramesFound-1) = superstruct.mvnx.subject.frames.frame(1:nFramesFound-1);
                    end
                    if firstFramesFound
                        fields = wordfields{iWord};
                        for il = 1:length(fields)
                            frames(nFramesFound).(fields{il}) = wordvalue{iWord}{il};
                        end
                    end
                end
                if isfield(obj, word{iWord})
                    obj.(word{iWord}).count = obj.(word{iWord}).count + 1;
                else
                    obj.(word{iWord}).count = 1;
                end
                name{index(iWord)} = word{iWord};
                namecounter(index(iWord)) = obj.(word{iWord}).count;
                if ~firstFramesFound && (iscell(wordvalue{iWord}) || ~isempty(wordvalue{iWord}) && all(~isnan(wordvalue{iWord})))
                    superstruct = setvalue(superstruct,index(iWord),wordvalue{iWord},name(1:index(iWord)),wordfields{iWord},namecounter);
                end
            elseif iWord > 1 && strcmp(word{iWord},word{iWord-1})
                name{index(iWord)} = word{iWord};
                namecounter(index(iWord)) = namecounter(index(iWord))+1;
                superstruct = setvalue(superstruct,index(iWord),wordvalue{iWord},name(1:index(iWord)),wordfields{iWord},namecounter);
            else
                name{index(iWord)} = word{iWord};
                namecounter(index(iWord)) = 1;
                if firstFramesFound && ~strcmp(word{iWord},'securityCode')
                    frames(nFramesFound).(word{iWord}) = wordvalue{iWord};
                else
                    superstruct = setvalue(superstruct,index(iWord),wordvalue{iWord},name(1:index(iWord)),wordfields{iWord},namecounter);
                end
            end
        catch e
            disp(getReport(e))
        end
    end
end
if firstFramesFound
    superstruct.mvnx.subject.frames.frame = frames';
    if obj.frame.number > nFramesFound
        superstruct.mvnx.subject.frames.frame(nFramesFound+1:end) = [];
    end
end

%% get output
mvnx = superstruct.mvnx;

delete(hWaitbar);
end
function superstruct = setvalue(superstruct,index,value,name,fields,namecounter)
%superstruct = setvalue(superstruct,index,value,name,fields,namecounter)
% add values to a very big struct
if ~isempty(fields) && (~iscell(fields) || ~isempty(fields{1}))
    if length(fields) == 1
        name(end+1) = fields;
    else
        name{end+1} = fields;
    end
    index = index+1;
end
switch index
    case 1
        if ~iscell(name{1})
            superstruct.(name{1}) = value;
        else
            for il = 1:length(name{1})
                superstruct.(name{1}{il}) = value{il};
            end
        end
    case 2
        if ~iscell(name{2})
            superstruct.(name{1})(namecounter(1)).(name{2}) = value;
        else
            for il = 1:length(name{2})
                superstruct.(name{1})(namecounter(1)).(name{2}{il}) = value{il};
            end
        end
    case 3
        if ~iscell(name{3})
            superstruct.(name{1})(namecounter(1)).(name{2})(namecounter(2)).(name{3})...
                = value;
        else
            for il = 1:length(name{3})
                superstruct.(name{1})(namecounter(1)).(name{2})(namecounter(2)).(name{3}{il}) = value{il};
            end
        end
    case 4
        if ~iscell(name{4})
            superstruct.(name{1})(namecounter(1)).(name{2})(namecounter(2)).(name{3})(namecounter(3))....
                .(name{4}) = value;
        else
            for il = 1:length(name{4})
                superstruct.(name{1})(namecounter(1)).(name{2})(namecounter(2)).(name{3})(namecounter(3))...
                    .(name{4}{il}) = value{il};
            end
        end
    case 5
        if ~iscell(name{5})
            superstruct.(name{1})(namecounter(1)).(name{2})(namecounter(2)).(name{3})(namecounter(3))....
                .(name{4})(namecounter(4)).(name{5}) = value;
        else
            for il = 1:length(name{5})
                superstruct.(name{1})(namecounter(1)).(name{2})(namecounter(2)).(name{3})(namecounter(3))....
                    .(name{4})(namecounter(4)).(name{5}{il}) = value{il};
            end
        end
    case 6
        if ~iscell(name{6})
            superstruct.(name{1})(namecounter(1)).(name{2})(namecounter(2)).(name{3})(namecounter(3))....
                .(name{4})(namecounter(4)).(name{5})(namecounter(5)).(name{6}) = value;
        else
            for il = 1:length(name{6})
                superstruct.(name{1})(namecounter(1)).(name{2})(namecounter(2)).(name{3})(namecounter(3))....
                    .(name{4})(namecounter(4)).(name{5})(namecounter(5)).(name{6}{il}) = value{il};
            end
        end
    case 7
        if ~iscell(name{7})
            superstruct.(name{1})(namecounter(1)).(name{2})(namecounter(2)).(name{3})(namecounter(3))....
                .(name{4})(namecounter(4)).(name{5})(namecounter(5)).(name{6})(namecounter(6))...
                .(name{7}) = value;
        else
            for il = 1:length(name{7})
                superstruct.(name{1})(namecounter(1)).(name{2})(namecounter(2)).(name{3})(namecounter(3))....
                    .(name{4})(namecounter(4)).(name{5})(namecounter(5)).(name{6})(namecounter(6))...
                    .(name{7}{il}) = value{il};
            end
        end
    case 8
        if ~iscell(name{8})
            superstruct.(name{1})(namecounter(1)).(name{2})(namecounter(2)).(name{3})(namecounter(3))....
                .(name{4})(namecounter(4)).(name{5})(namecounter(5)).(name{6})(namecounter(6))...
                .(name{7})(namecounter(7)).(name{8}) = value;
        else
            for il = 1:length(name{7})
                superstruct.(name{1})(namecounter(1)).(name{2})(namecounter(2)).(name{3})(namecounter(3))....
                    .(name{4})(namecounter(4)).(name{5})(namecounter(5)).(name{6})(namecounter(6))...
                    .(name{7})(namecounter(7)).(name{8}{il}) = value{il};
            end
        end
end
end

function word = checkText(word)
%word = checkText(word)
% make sure the field name is just text and allowed symbols.
if ~isempty(word)
    if ~iscell(word)
        error('should get cells')
    end
    for i=1:length(word)
        if length(word{i})>1
            word{i}(word{i}=='!') = [];
            word{i}(word{i}==':') = '_';
            if word{i}(end) == '/'
                word{i}(end) = [];
            end
            while (word{i}(1) < 65 || word{i}(1) > 122 || (word{i}(1)> 90 && word{i}(1) < 97))
                word{i}(1)=[];
            end
        end
    end
end
end


text:
function [words, k] = compute_features(I, I_rgb, I_type, opts)
% Computes the features specified in spagglom_options.m. Variable. 'words'
% is a cell array, each cell of which is a matrix having the same size as
% the image. Matrix elements denote histogram bins, bin 0 meaning no-value.
% Output 'k' contains number of bins for each feature matrix.

% 'I' is the image in transformed format, or a duplicate of 'I_rgb',
% which must always be in RGB format. 'I_type' is a string specifying the
% format of 'I'.

features_num = 0;
k = [];
words = cell(0);

% denseSIFT with bag-of-words clustering
if opts.feature_dsift_bow
    features_num = features_num + 1;
    [words{features_num}, kt] = feature_dsift(I_rgb, opts); % dsift uses only the rgb image
    k(features_num) = kt;
end

% Cluster colors using k-means
if opts.feature_color_bow
    features_num = features_num + 1;
    [words{features_num}, kt] = feature_color(I, I_type, opts);
    k(features_num) = kt;
end

% Raw rgb histogram
if opts.feature_rgb_raw
    features_num = features_num + 1;
    [words{features_num}, kt] = feature_rgb_raw(I, I_type, opts);
    k(features_num) = kt;
end

% Texture feature used by van de Sande (may not be identical)
if opts.feature_grad_texture
    features_num = features_num + 1;
    [words{features_num}, kt] = feature_grad_texture2(I, opts); % This feature is calculated in three channels, so the 'words' variable unusually has size [3*h,w].
    k(features_num) = kt;
end

% Local binary patterns
if opts.feature_lbp
    features_num = features_num + 1;
    [words{features_num}, kt] = feature_lbp(I_rgb, opts); % rgb images only
    k(features_num) = kt;
end

author: Cloud-CV

repo: object-proposals


text:
function [controlFlux, objFlux] = runRobustnessAnalysis(model, controlRxn, nPoints, objRxn, plotRedCost)
% runRobustnessAnalysis
%   Performs robustness analysis for a reaction of interest and an objective
%   of interest. Modified from the COBRA robustnessAnalysis function.
%
% Input:
%   model           a model structure
%   controlRxn      reaction of interest whose value is to be controlled
%   nPoints         number of points to show on plot (opt, default 20)
%   objRxn          reaction identifier of objective to be maximized (opt,
%                   default it uses the objective defined in the model)
%   plotRedCost     logical whether reduced cost should also be plotted
%                   (opt, default false)
%
% Output:
%   controlFlux     flux values of the reaction of interest, ranging from
%                   its minimum to its maximum value
%   objFlux         optimal values of objective reaction at each control
%                   reaction flux value
%
% Modified from COBRA Toolbox robustnessAnalysis.m
%
% Usage: runRobustnessAnalysis(model, controlRxn, nPoints, objRxn)

if nargin < 3
    nPoints = 20;
end
if nargin < 4
    baseModel = model;
else
    baseModel = setParam(model,'obj',objRxn,1);
end
if nargin < 5
    plotRedCost = false;
end

if any(ismember(model.rxns,controlRxn))
    controlRxnIdx = getIndexes(model,controlRxn,'rxns');
    tmpModel = setParam(model,'obj',controlRxnIdx,1);
    solMax = solveLP(tmpModel);
    solMax = solMax.x(logical(tmpModel.c));
    tmpModel.c = -tmpModel.c;
    solMin = solveLP(tmpModel);
    solMin = solMin.x(logical(tmpModel.c));
else
    error('Control reaction does not exist!');
end

objFlux = zeros(nPoints,1);
redCost = zeros(nPoints,1);
controlFlux = linspace(solMin,solMax,nPoints)';

fprintf('Running robustness analysis...   0%% complete');
for i=1:length(controlFlux)
    progress=pad(num2str(floor(i/numel(controlFlux)*100)),3,'left');
    fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b%s%% complete',progress);
    modelControlled = setParam(baseModel,'eq',controlRxnIdx,controlFlux(i));
    solControlled = solveLP(modelControlled);
    objFlux(i) = solControlled.x(logical(modelControlled.c));
    redCost(i) = solControlled.rCost(controlRxnIdx);
end
fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\bCOMPLETE\n');

if plotRedCost
    yyaxis right
    plot(controlFlux,redCost)
    ylabel([strrep(controlRxn,'_','-') ' reduced cost']);
    yyaxis left
end
plot(controlFlux,objFlux)
xlabel(strrep(controlRxn,'_','-'));
ylabel('Objective');
end

author: SysBioChalmers

repo: RAVEN


text:
function u=transp(u0,lambda,dir,maxshift)
% Transports u0 a distance lambda and averages over grid (which is assumed
% to be equispaced) Boundary conditions are Neumann. It is easy to modify
% for other boundary conditions.
if nargin<3,
  dir=1;
end;
if (dir>2),
  error(strcat('You can easily extend this to',num2str(dir)','D yourself ...'));
end;
% Some things ensuring that you solve along the columns of u. Must
% transpose if u is a row vector or if dir = 2.
transpose=0;
S=size(u0);
N=S(1);
if ((N==1)||(dir==2)),
  u0=u0';
  transpose=1;
end
S=size(u0);
N=S(1); M=S(2);
Sl=size(lambda);
if ((Sl(1)==Sl(2))&&(Sl(1)==1)),
	lambda=lambda*ones(1,M);
elseif (Sl(2)==1)&&(Sl(1)>1),
	lambda=lambda';
end;
Sl=size(lambda);
if (Sl(2)~=M),
	error('The shift must have the right dimension');
end;
if nargin<4,
	maxshift=ceil(max(abs(lambda)))+1;
end;
shift=floor(lambda);
alpha=lambda-shift;
eshift=ones(maxshift,1);
uh=[eshift*u0(1,1:M); u0; eshift*u0(N,1:M)];
e=ones(N,1);
alpha=e*alpha;   % Making alpha a square matrix the size of u0
uh=colshift(uh,shift);
u=alpha.*uh(maxshift:maxshift+N-1,:)+(1-alpha).*uh(maxshift+1:maxshift+N,:);
if transpose,
	u=u';
end;

repo: Aero-matlab


text:
function hole_num = p00_hole_num ( test )

%*****************************************************************************80
%
%% P00_HOLE_NUM counts the holes in a problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TEST, the index of the test problem
%
%    Output, integer HOLE_NUM, the number of boundary segments.
%
  if ( test == 1 )
    hole_num = p01_hole_num ( );
  elseif ( test == 2 )
    hole_num = p02_hole_num ( );
  elseif ( test == 3 )
    hole_num = p03_hole_num ( );
  elseif ( test == 4 )
    hole_num = p04_hole_num ( );
  elseif ( test == 5 )
    hole_num = p05_hole_num ( );
  elseif ( test == 6 )
    hole_num = p06_hole_num ( );
  elseif ( test == 7 )
    hole_num = p07_hole_num ( );
  elseif ( test == 8 )
    hole_num = p08_hole_num ( );
  elseif ( test == 9 )
    hole_num = p09_hole_num ( );
  elseif ( test == 10 )
    hole_num = p10_hole_num ( );
  elseif ( test == 11 )
    hole_num = p11_hole_num ( );
  elseif ( test == 12 )
    hole_num = p12_hole_num ( );
  elseif ( test == 13 )
    hole_num = p13_hole_num ( );
  elseif ( test == 14 )
    hole_num = p14_hole_num ( );
  elseif ( test == 15 )
    hole_num = p15_hole_num ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_HOLE_NUM - Fatal error!\n' );
    fprintf ( 1, '  Input value of TEST is out of range.\n' );
    fprintf ( 1, '  TEST = %d\n', test );
    error ( 'P00_HOLE_NUM - Fatal error!' )
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function eg2OC_BVP
%EG2OC_BVP    Example 2 of optimal control tutorial.
%    This example is from D.E.Kirk's Optimal control theory: an introduction
%    example 6.2-2 on page 338-339
%    The problem is reformulated as a boundary value problem (BVP)
%    and solved with bvp4c solver

% Initial guess for the solution
solinit = bvpinit(linspace(0,0.78,50), [0 0 0.5 0.5]);
options = bvpset('Stats','on','RelTol',1e-1);
global R;
R = 0.1;
sol = bvp4c(@BVP_ode, @BVP_bc, solinit, options);
t = sol.x;
y = sol.y;

% Calculate u(t) from x1,x2,p1,p2
ut = (y(3,:).*(y(1,:) + 1/4))/(2*0.1);
n = length(t);
% Calculate the cost
J = 0.78*(y(1,:)*y(1,:)' + y(2,:)*y(2,:)' + ut*ut'*0.1)/n;

figure(1);
plot(t, y(1:2,:)','-');hold on;
plot(t,ut', 'r:')
text(.2,0.08,'x_1(t)');
text(.3,-.1,'x_2(t)');
text(.2,.4, 'u(t)');
s = strcat('Final cost is: J=',num2str(J));
text(.4,1,s);
xlabel('time');
ylabel('states');
hold off;
print -djpeg90 -r300 eg2_bvp.jpg

%------------------------------------------------
% ODE's for states and costates
function dydt = BVP_ode(t,y)
global R;
t1 = y(1)+.25;
t2 = y(2)+.5;
t3 = exp(25*y(1)/(y(2)+2));
t4 = 50/(y(1)+2)^2;
u =  y(3)*t1/(2*R);

dydt = [-2*t1+t2*t3-t2*u
        0.5-y(2)-t2*t3
        -2*y(1)+2*y(3)-y(3)*t2*t4*t3+y(3)*u+y(4)*t2*t4*t3
        -2*y(2)-y(3)*t3+y(4)*(1+t3)];

% -----------------------------------------------
% The boundary conditions:
% x1(0) = 0.05, x2(0) = 0, tf = 0.78, p1(tf) = 0, p2(tf) = 0;
function res = BVP_bc(ya,yb)
res = [ ya(1) - 0.05
        ya(2) - 0
        yb(3) - 0
        yb(4) - 0 ];

author: Sable

repo: mcbench-benchmarks


text:
function plot2dGaussian(Priors,Mus,Sigmas)
%PLOT2DGAUSSIAN Summary of this function goes here
%   Detailed explanation goes here

disp('in plot2dGaussians');

GMM=struct('Priors',Priors,'Mus',Mus,'Sigmas',Sigmas,'K',size(Mus,2));

spacing=1000;
xGMM=marginalizeGMM(GMM,1);
xs=linspace(min(xGMM.Mus'-2*sqrt(squeeze(xGMM.Sigmas))),max(xGMM.Mus'+2*sqrt(squeeze(xGMM.Sigmas))),spacing);
yGMM=marginalizeGMM(GMM,2);
ys=linspace(min(yGMM.Mus'-2*sqrt(squeeze(yGMM.Sigmas))),max(yGMM.Mus'+2*sqrt(squeeze(yGMM.Sigmas))),spacing);
pdraw=zeros(1,spacing^2);
[X Y]=meshgrid(xs,ys);
x=X(:);
y=Y(:);
vtests2=[x y]';
GMM.K
for k=1:GMM.K
    pdraw=pdraw+GMM.Priors(k).*gaussPDF(vtests2,GMM.Mus(:,k),GMM.Sigmas(:,:,k))';
end

ppdraw=reshape(pdraw,spacing,spacing);

contourf(xs,ys,ppdraw);
%pcolor(xs,ys,ppdraw);
%shading interp;

set(gca,'YDir','normal');

end

author: epfl-lasa


text:
function [Tcell, Fcell, Ycell, meanf, peakf] = spectrogram( w, varargin )
%SPECTROGRAM Plot an IceWeb-style spectrogram
%   spectrogram(w, s) Creates an IceWeb style spectrogram by wrapping the
%   function iceweb.spectrogram_iceweb(). If s is omitted it defaults to:
%       spectralobject(1024, 924, 10, [60 120])
%   where [60 120] are the dB limits that correspond to the color map used.
%   If this is a blank vector, [], the spectrogram is autoscaled in absolute
%   amplitude, not plotted in dB with absolute dB limits.
%
%   spectrogram(w, s, 'plot_metrics', 1) Will add frequency metrics computed
%   by the waveform/amplitude_spectrum function.
%
% For greater control, call
%   iceweb.spectrogram_iceweb() directly, or use spectralobject/specgram or
%   spectralobject/specgram2 (not clear how these differ). Note that
%   spectrogram_iceweb() is significantly faster.

% Glenn Thompson 2016/05/19 to provide a simple, fast way to geenrate nice
% spectrograms without having to delve into specgram and specgram2
figure
nfft = 1024;
overlap = 924;
fmax = 10;
dbLims = [60 120];
p = inputParser;
p.addRequired('w');
p.addParameter('spectralobject', spectralobject(nfft, overlap, fmax, dbLims));
p.addParameter('plot_metrics', 0, @isnumeric);
p.parse(w, varargin{:});
w=p.Results.w;
s=p.Results.spectralobject;

% if ~class(s, 'spectralobject')
%     disp('Oops, you did not give a valid spectralobject. Using default')
    if numel(w)>1
        w = reshape(w, numel(w), 1);
    end

    if p.Results.plot_metrics
        [result,Tcell,Fcell,Ycell, meanf, peakf] = iceweb.spectrogram_iceweb(s, w, 'plot_metrics',1);
    else
        [result,Tcell,Fcell,Ycell, meanf, peakf] = iceweb.spectrogram_iceweb(s, w);
    end

end

author: geoscience-community-codes

repo: GISMO


text:
function out = get(f, prop)
%GET   GET method for the CHEBTECH1 class
%   P = GET(F, PROP) returns the property P specified in the string PROP from
%   the CHEBTECH1 F.  The string PROP may be the name of a CHEBTECH1 property
%   (see the CHEBTECH and CHEBTECH1 classdef files for lists) or any of the
%   following strings:
%       'POINTS'          - 1st-kind Chebyshev grid corresponding to F.
%       'VALUES'          - Values of F at Chebyshev points.
%       'LVAL'            - Value of F at -1.
%       'RVAL'            - Value of F at +1.
%       'TECH'            - Handle to the CHEBTECH1 constructor. *
%       'VSCALE'          - Vertical scale of the CHEBTECH1.
%
% See also CHEBTECH, CHEBTECH1.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEVELOPER NOTE:
%  * Currently get(f, 'tech') returns a function handle to the tech constructor.
%    This may change in future to return instead an empty instance of the tech.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch prop
    case fieldnames(f)
        out = f.(prop);
    case 'points'
        out = f.points();
    case 'lval'
        out = lval(f);
    case 'rval'
        out = rval(f);
    case 'values'
        out = f.coeffs2vals(f.coeffs);
    case 'vscale'
        out = vscale(f);
    case 'tech'
        % TODO: Return function handle, or empty instance of the tech?
        out = @chebtech1;
    otherwise
        error('CHEBFUN:CHEBTECH1:get:propName', ...
            'Unknown property name ''%s'' for object of type CHEBTECH1.', prop);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% This code was originally written and distributed as part of the
% PASCAL VOC challenge
function ap = VOCap(rec,prec)

mrec=[0 ; rec ; 1];
mpre=[0 ; prec ; 0];
for i=numel(mpre)-1:-1:1
    mpre(i)=max(mpre(i),mpre(i+1));
end
i=find(mrec(2:end)~=mrec(1:end-1))+1;
ap=sum((mrec(i)-mrec(i-1)).*mpre(i));

author: zhoubolei

repo: CAM

save_path: github-repos/MATLAB/zhoubolei-CAM


text:
function DupEdit(topH)

% DupEdit(topH)
%
% Duplicate all the editable field contents in the present scan
% forward into all subsequent session fields. Input topH is the
% handle of the top-level edit figure.
%
% DBR 4/99
% DAR 1/22/07 - edited for compatibility with PC matlab 7.1

% Unpack the user data from the edit figure
uiData = get(topH, 'UserData');

iS = uiData.iScan; % Present scan index
nScans = length(uiData.session.functionals); % Number of scans in session

if iS == nScans, return; end % Nothing to copy
nFields = length(uiData.scanData); % Number of scan edit fields
for iScan=iS+1:nScans
    for iField=1:nFields
        editFlag = uiData.scanData(iField).edit;
        if editFlag
            hContent = uiData.scanData(iField).handle;
            fName = uiData.scanData(iField).field;
            origValue = uiData.original.functionals(iScan).(fName);
            curValue = get(hContent,'string');
            if isnumeric(origValue)
                curValue = str2num(curValue);
            end

            label = uiData.scanData(iField).label;
            uiData.session.functionals(iScan).(fName) = curValue;
        end
    end
end
set(topH, 'UserData', uiData);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function p4 = angle_half ( p1, p2, p3 )

%*****************************************************************************80
%
%% ANGLE_HALF finds half an angle.
%
%  Discussion:
%
%    The original angle is defined by the sequence of points P1, P2 and P3.
%
%    The point P4 is calculated so that:
%
%      (P1,P2,P4) = (P1,P2,P3) / 2
%
%        P1
%        /
%       /   P4
%      /  .
%     / .
%    P2--------->P3
%
%    Thanks to Cesar Fraga Bobis for pointing out a typographical error in
%    a previous version of this routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2,1), P2(2,1), P3(2,1), points defining the angle.
%
%    Input, real P4(2,1), a point defining the half angle.
%    The vector P4 - P2 will have unit norm.
%
  p4(1:2,1) = 0.5 * ( ...
      ( p1(1:2,1) - p2(1:2,1) ) / sqrt ( sum ( ( p1(1:2,1) - p2(1:2,1) ).^2 ) ) ...
    + ( p3(1:2,1) - p2(1:2,1) ) / sqrt ( sum ( ( p3(1:2,1) - p2(1:2,1) ).^2 ) ) );

  p4(1:2,1) = p2(1:2,1) + p4(1:2,1) / sqrt ( sum ( ( p4(1:2,1) ).^2 ) );

  return
end


text:
 % Data File DTPC3
 % Dynamics of a particle, moving
 % on a horizontal plane under action
 % of a body, bound up to the particle
 % by a non elastic thread. The body can
 % move in vertical direction.
   m    =  'm';  % mass of the particle
          % Projections of forces on axis [r] and [f]
   Fr   = '-G*(1+rtt/9.81)-9.81*mu*m*rt/sqrt(rt^2+r^2*ft^2)';
   Ff   = '-9.81*mu*m*r*ft/sqrt(rt^2+r^2*ft^2)';
   r0   = 1;     % Initial radius
   f0   = 0;     % Initial polar angle
   v0   = 10;    % Initial velocity
   alfa = pi/2;  % angle between v0 and polar axis
   Tend = 5;     % upper bound of integration
   eps  = 1e-10; % desirable accuracy
   np   = 3;     % number of parameters
   P{1} = 'm';   % mass of the particle
   P{2} = 'G';   % weight of the body
   P{3} = 'mu';  % coefficient of dry friction


text:
function outsig=operator(Op,insig);
%OPERATOR  Apply operator
%   Usage: c=operator(Op,f);
%
%   `c=operator(Op,f)` applies the operator *Op* to the input signal *f*.
%   The operator object *Op* must have been created using |operatornew|.
%
%   If *f* is a matrix, the transform will be applied along the columns
%   of *f*. If *f* is an N-D array, the transform will be applied along
%   the first non-singleton dimension.
%
%   See also: operatornew, ioperator, operatoradj

if nargin<2
  error('%s: Too few input parameters.',upper(mfilename));
end;

if ~isstruct(Op)
  error('%s: First agument must be a operator definition structure.',upper(mfilename));
end;

switch(Op.type)
  case 'framemul'
    outsig=framemul(insig,Op.Fa,Op.Fs,Op.s);
  case 'spread'
    outsig=spreadop(insig,Op.s);
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function varargout = decode(varargin)
% DECODE M-file for decode.fig
%      DECODE, by itself, creates a new DECODE or raises the existing
%      singleton*.
%
%      H = DECODE returns the handle to a new DECODE or the handle to
%      the existing singleton*.
%
%      DECODE('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in DECODE.M with the given input arguments.
%
%      DECODE('Property','Value',...) creates a new DECODE or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before decode_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to decode_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help decode

% Last Modified by GUIDE v2.5 27-Jul-2003 22:01:41

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @decode_OpeningFcn, ...
                   'gui_OutputFcn',  @decode_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before decode is made visible.
function decode_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to decode (see VARARGIN)

% Choose default command line output for decode
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes decode wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = decode_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in b1.
function b1_Callback(hObject, eventdata, handles)
% hObject    handle to b1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=697;f2=1209;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;
sound(y,fs)
subdecode;
% --- Executes on button press in b2.
function b2_Callback(hObject, eventdata, handles)
% hObject    handle to b2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=697;f2=1336;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;sound(y,fs)
subdecode;
% --- Executes on button press in A.
function A_Callback(hObject, eventdata, handles)
% hObject    handle to A (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=697;f2=1663;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;sound(y,fs)
subdecode;
% --- Executes on button press in b3.
function b3_Callback(hObject, eventdata, handles)
% hObject    handle to b3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=697;f2=1447;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;sound(y,fs)
subdecode;
% --- Executes on button press in b4.
function b4_Callback(hObject, eventdata, handles)
% hObject    handle to b4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=770;f2=1209;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;sound(y,fs)
subdecode;
% --- Executes on button press in b5.
function b5_Callback(hObject, eventdata, handles)
% hObject    handle to b5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=770;f2=1336;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;sound(y,fs)
subdecode;
% --- Executes on button press in B.
function B_Callback(hObject, eventdata, handles)
% hObject    handle to B (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=770;f2=1633;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;sound(y,fs)
subdecode;
% --- Executes on button press in b6.
function b6_Callback(hObject, eventdata, handles)
% hObject    handle to b6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=770;f2=1477;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;sound(y,fs)
subdecode;
%--- Executes on button press in b7.
function b7_Callback(hObject, eventdata, handles)
% hObject    handle to b7 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=852;f2=1209;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;sound(y,fs)
subdecode;
% --- Executes on button press in b8.
function b8_Callback(hObject, eventdata, handles)
% hObject    handle to b8 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=852;f2=1336;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;sound(y,fs)
subdecode;
% --- Executes on button press in C.
function C_Callback(hObject, eventdata, handles)
% hObject    handle to C (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=852;f2=1633;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;sound(y,fs)
subdecode;
% --- Executes on button press in b9.
function b9_Callback(hObject, eventdata, handles)
% hObject    handle to b9 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=852;f2=1477;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;sound(y,fs)
subdecode;
% --- Executes on button press in ba.
function ba_Callback(hObject, eventdata, handles)
% hObject    handle to ba (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=941;f2=1209;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;sound(y,fs)
subdecode;
% --- Executes on button press in b0.
function b0_Callback(hObject, eventdata, handles)
% hObject    handle to b0 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=941;f2=1336;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;sound(y,fs)
subdecode;
% --- Executes on button press in D.
function D_Callback(hObject, eventdata, handles)
% hObject    handle to D (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=941;f2=1633;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;sound(y,fs)
subdecode;
% --- Executes on button press in bn.
function bn_Callback(hObject, eventdata, handles)
% hObject    handle to bn (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
t=[0:0.000125:.05];
fs=8000;
f1=941;f2=1477;
y1=.25*sin(2*pi*f1*t);
y2=.25*sin(2*pi*f2*t);
y=y1+y2;sound(y,fs);
subdecode;

% --- Executes on button press in info.
function info_Callback(hObject, eventdata, handles)
% hObject    handle to info (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
msgbox('File was created by: Randolf C. Sequera                                       BSECE Adamson University Philippines','Info','warn')

% --- Executes on button press in close.
function close_Callback(hObject, eventdata, handles)
% hObject    handle to close (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
close;


text:
%  returns a 3-column matrix in which:
%  - The first column represents the rows
%  - The second column represents the columns
%  - The third column represents the values expander, such that vals(col3)
%  returns a vector of same length as the rows of matrix B
%
%  A is a cell array of cell arrays
%  - in each entry of A is a cell array
%    - where the first element is the row of the matrix
%    - the remaining entries are scalars or vectors that indicate the degree
%    of multiplicity of the corresponding derivative and the location of
%    those derivatives in the final matrix
%

author: jmaih


text:
function FEP_fluctuations
% This demonstration  uses an ensemble of particles with intrinsic (Lorentz
% attractor) dynamics and (Newtonian) short-range coupling.  The focus of
% this routine is to unpack the Bayesian perspective. We first simulate
% dynamics  to nonequilibrium steady-state, identify the Markov blanket and
% then examine the encoding of external states by internal states; in terms
% of their expected values.
%
% The crucial aspect of this implicit inference (and the basis of the free
% energy principle) is the existence of a conditional synchronisation
% manifold, when conditioning internal and external states on the Markov
% blanket. This provides the basis for a mapping between internal and
% external states that can be interpreted in terms of a probabilistic
% representation or inference.
%
% This Bayesian perspective is illustrated in terms of a mapping between
% the canonical modes of internal and external states (as approximated
% with a polynomial expansion). The canonical modes her are evaluated
% using an estimate of the conditional expectations  based upon the
% Euclidean proximity of Markov blanket states. The ensuing posterior over
% external states is than illustrated, in relation to the actual external
% states. We also  simulate event related potentials by identifying
% several points in time when the Markov blankets revisit the same
% neighbourhood. Finally, to illustrate the underlying dynamics, the
% Jacobians or coupling among internal and external states are
% presented; using different orders of coupling (i.e., degrees of
% separation)
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: FEP_fluctuations.m 7163 2017-09-04 09:12:50Z karl $

% default settings (GRAPHICS sets movies)
%--------------------------------------------------------------------------
rng('default')
GRAPHICS = 1;

% Demo of synchronization manifold using coupled Lorenz attractors
%==========================================================================
N    = 128;                         % number of (Lorenz) oscillators
T    = 2048;                        % number of time bins
dt   = 1/32;                        % time interval

% parameters
%--------------------------------------------------------------------------
P.k  = 1 - exp(-rand(1,N)*4);       % variations in temporal scale
P.d  = 1/8;                         % amplitude of random fluctuations

% states
%--------------------------------------------------------------------------
x.p  = randn(2,N)*4;                % microstates (position)
x.v  = zeros(2,N);                  % microstates (velocity)
x.q  = randn(3,N)/32;               % microstates (states)
u    = zeros(1,T);                  % exogenous fluctuations

% generate an dynamics from initial conditions
%==========================================================================
spm_figure('GetWin','Markov blanket');clf
if GRAPHICS
    subplot(2,2,1)
else
    subplot(2,1,1)
end
[Q,X,V,A,x] = spm_soup(x,u,P,T,dt,1);

% States
%--------------------------------------------------------------------------
% Q    - history of microstates (states)
% X    - history of microstates (position)
% V    - history of microstates (velocity)

for i = 1:size(X,3)
    S(:,:,i) = [Q(:,:,i);X(:,:,i);V(:,:,i)];
end

% Markov blanket - parents, children, and parents of children
%==========================================================================

% Adjacency matrix
%--------------------------------------------------------------------------
t     = (T - 256):T;                              % final time indices
L     = sparse(double(any(A(:,:,t),3)))';

% internal states (defined by principle eigenvector of Markov blanket)
%--------------------------------------------------------------------------
B     = double((L + L' + L'*L));
B     = B - diag(diag(B));
v     = spm_svd(B*B',1);
[v,j] = sort(abs(v(:,1)),'descend');

% get Markov blanket and divide into sensory and active states
%--------------------------------------------------------------------------
m     = j(1:8);                                   % internal cluster
mm    = sparse(m,1,1,N,1);                        % internal states
bb    = B*mm & (1 - mm);                          % Markov blanket
ee    = 1 - bb - mm;                              % external states
b     = find(bb);
e     = find(ee);
m     = find(mm);
s     = b(find( any(L(b,e),2)));
a     = b(find(~any(L(b,e),2)));

% adjacency matrix - with partition underneath (LL)
%--------------------------------------------------------------------------
k       = [e; s; a; m];
LL      = L;
LL(e,e) = LL(e,e) + 1/8;
LL(s,s) = LL(s,s) + 1/8;
LL(a,a) = LL(a,a) + 1/8;
LL(m,m) = LL(m,m) + 1/8;
LL      = LL(k,k);

% plot dynamics for the initial and subsequent time periods
%--------------------------------------------------------------------------
subplot(4,1,3)
r    = 1:512;
plot(r,squeeze(Q(1,e,r)),':c'), hold on
plot(r,squeeze(Q(1,m,r)),' b'), hold off
axis([r(1) r(end) -32 32])
xlabel('Time','FontSize',12)
title('Electrochemical dynamics','FontSize',16)

subplot(4,1,4)
r    = 1:T;
plot(r,squeeze(V(1,e,r)),':c'), hold on
plot(r,squeeze(V(1,m,r)),' b'), hold off
axis([r(1) r(end) -32 32])
xlabel('Time','FontSize',12)
title('Newtonian dymanics','FontSize',16)

% Markov blanket - self-assembly
%==========================================================================
subplot(2,2,1)
imagesc(1 - LL)
axis square
xlabel('Element','FontSize',12)
xlabel('Element','FontSize',12)
title('Adjacency matrix','FontSize',16)

% follow self-assembly
%--------------------------------------------------------------------------
clear M
for i = (T - 512):T

    % plot positions
    %----------------------------------------------------------------------
    subplot(2,2,2),set(gca,'color','w')

    px = ones(3,1)*X(1,:,i) + Q([1 2 3],:,i)/16;
    py = ones(3,1)*X(2,:,i) + Q([2 3 1],:,i)/16;
    plot(px,py,'.b','MarkerSize',8), hold on
    px = X(1,e,i); py = X(2,e,i);
    plot(px,py,'.c','MarkerSize',24)
    px = X(1,m,i); py = X(2,m,i);
    plot(px,py,'.b','MarkerSize',24)
    px = X(1,s,i); py = X(2,s,i);
    plot(px,py,'.m','MarkerSize',24)
    px = X(1,a,i); py = X(2,a,i);
    plot(px,py,'.r','MarkerSize',24)

    xlabel('Position','FontSize',12)
    ylabel('Position','FontSize',12)
    title('Markov Blanket','FontSize',16)
    axis([-1 1 -1 1]*8)
    axis square, hold off, drawnow

    % save
    %----------------------------------------------------------------------
    if i > (T - 128) && GRAPHICS
        M(i - T + 128) = getframe(gca);
    end

end

% set ButtonDownFcn
%--------------------------------------------------------------------------
if GRAPHICS
    h   = findobj(gca);
    set(h(1),'Userdata',{M,16})
    set(h(1),'ButtonDownFcn','spm_DEM_ButtonDownFcn')
    xlabel('Click for Movie','Color','r')
end

% illustrate the Bayesian perspective (predictability of external states)
%==========================================================================
spm_figure('GetWin','Bayesian perspective');clf

% establish a statistical dependency between internal (dynamic) states (XQ)
%--------------------------------------------------------------------------
T     = 512;                                   % length of timeseries
t     = size(X,3) - T - 2;
for i = 1:T
    Xe(i,:) = spm_vec(V(:,e,i + t));           % external states
    Xb(i,:) = spm_vec(S(:,[a;s],i + t));       % Markov blanket
    Xm(i,:) = spm_vec(Q(:,m,i + t));           % internal states
end
xe    = zeros(size(Xe));
xm    = zeros(size(Xm));
iC    = inv(cov(Xb));

% probabilistic proximity in the space of the Markov blanket
%--------------------------------------------------------------------------
for i = 1:T
    for j = 1:T
        r      = Xb(i,:) - Xb(j,:);
        w(i,j) = exp(-(r*iC*r')/128);
    end
end

% convert into proper probability distribution
%--------------------------------------------------------------------------
w = diag(sum(w,2))\w;

% mean
%--------------------------------------------------------------------------
for i = 1:T
    for j = 1:T
        xe(i,:) = xe(i,:) + w(i,j)*Xe(j,:);
        xm(i,:) = xm(i,:) + w(i,j)*Xm(j,:);
    end
end

% covariance (not used)
%--------------------------------------------------------------------------
% ce    = zeros(size(Xe,1),size(Xe,2),size(Xe,2));
% for i = 1:T
%     for j = 1:T
%         ce(i,:,:) = squeeze(ce(i,:,:)) + w(i,j)*(Xe(i,:) - xe(i,:))'*(Xe(j,:) - xe(j,:));
%     end
% end

% normalise and identify canonical eigenvariates
%--------------------------------------------------------------------------
xe    = spm_detrend(xe);
xm    = spm_detrend(xm);
CVA   = spm_cva(xe,xm);

% show results - canonical vectors over elements (mode M)
%--------------------------------------------------------------------------
subplot(3,2,1)

M     = 1;
Ve    = CVA.V(:,M);
Ve    = spm_unvec(Ve,V(:,e,1));
ve    = sum(Ve.^2);
ve    = ve/max(ve);
for k = 1:length(Ve)
    c = [0 1 1]*ve(k) + [1 1 1]*(1 - ve(k));
    plot(X(1,e(k),end),X(2,e(k),end),'.','MarkerSize',32,'Color',c), hold on
end

% overplot mode of motion
%--------------------------------------------------------------------------
quiver(X(1,e,end),X(2,e,end),Ve(1,:),Ve(2,:))

Vm    = CVA.W(:,M);
Vm    = spm_unvec(Vm,Q(:,m,1));
vm    = sum(Vm.^2);
vm    = vm/max(vm);
for k = 1:length(Vm)
    c = [0 0 1]*vm(k) + [1 1 1]*(1 - vm(k));
    plot(X(1,m(k),end),X(2,m(k),end),'.','MarkerSize',32,'Color',c), hold on
end
xlabel('Position', 'FontSize',12)
ylabel('Position','FontSize',12)
title('Canonical mode','FontSize',16)

% conditional synchronisation manifold (polynomial approximation)
%==========================================================================

% polynomial approximation
%--------------------------------------------------------------------------
xX    = xm*CVA.W(:,M);
XX    = [xX.^0 xX.^1 xX.^2 xX.^3 xX.^4 xX.^5];
bE    = pinv(XX)*Xe*CVA.V(:,M);
qE    = XX*bE;

% conditional expectation and variance
%--------------------------------------------------------------------------
% for i = 1:T
%     qC(i) = CVA.V(:,1)'*squeeze(ce(i,:,:))*CVA.V(:,1);
% end
qC    = ones(size(qE));
qC    = abs(var(Xe*CVA.V(:,1) - qE)*qC/mean(qC));

% show results - conditional synchronisation manifold
%--------------------------------------------------------------------------
subplot(3,2,2)
plot(CVA.w(:,1),Xe*CVA.V(:,1),'.c' ), hold on
plot(CVA.w(:,1),qE,'.b' ), hold off
xlabel('Internal mode', 'FontSize',12)
ylabel('External mode','FontSize',12)
title('Synchronisation manifold','FontSize',16), spm_axis tight

% show results - conditional distributions as a function of time
%--------------------------------------------------------------------------
subplot(3,1,2)
plot(Xe*CVA.V(:,1),'c' ), hold on
spm_plot_ci(qE',qC(:)'),  hold off
xlabel('Time', 'FontSize',12)
ylabel('External states','FontSize',12)
title('Inferred and real motion','FontSize',16), spm_axis tight

%  event related potentials
%==========================================================================

%  identify points of interest using the external canonical variate
%--------------------------------------------------------------------------
u     = CVA.v(:,1);
ue    = [];
for i = 1:8
    [d,j]   = max(u(33:end - 65));
    j       = j + 32;
    % j     = fix(rand*T);              % random times
    k       = (j - 32):(j + 64);        % perstimulus time (around j)
    u((j - 8):(j + 8)) = -Inf;          % eliminate from next max(u(:,1))
    ue(:,i) = Xe(k,:)*CVA.V(:,1);
    um(:,i) = Xm(k,:)*CVA.W(:,1);
    us(i)   = j;

end
j    = any(ue);
us   = us(j);
ue   = spm_detrend(ue(:,j));
um   = spm_detrend(um(:,j));

% plot points of interest on conditional density
%--------------------------------------------------------------------------
subplot(3,1,2)
uy    = get(gca,'Ylim'); hold on
for i = 1:length(us),plot([1 1]*us(i),uy,':'), end, hold off

%  show time locked (internal and external) fluctuations and their mean
%--------------------------------------------------------------------------
subplot(3,2,5)
pst   = (-32:64)*8;
plot(pst,ue,'c:',pst,um,'b:'),               hold on
plot(pst,mean(ue,2),'c',pst,mean(um,2),'b'), hold on
plot([0 0],get(gca,'YLim'),'--'),            hold off, axis square
xlabel('Time (milliseconds)', 'FontSize',12)
ylabel('Electrochemical response','FontSize',12)
title('Simulated ERP','FontSize',16)

% canonical correlations
%--------------------------------------------------------------------------
subplot(3,2,6)
bar(CVA.r,1/2,'c')
xlabel('Mode', 'FontSize',12)
ylabel('Correlation','FontSize',12)
title('Canonical correlations','FontSize',16), axis square

% Jacobian's and generalised synchronisation
%==========================================================================
spm_figure('GetWin','Jacobians');clf

% get Markov blanket indices the Jacobian
%--------------------------------------------------------------------------
xi    = spm_zeros(x); xi.v(:,e) = 1; iXe = find(spm_vec(xi));
xi    = spm_zeros(x); xi.q(:,m) = 1; iXm = find(spm_vec(xi));

% show results - Jacobians (of increasing order: 1 to n)
%--------------------------------------------------------------------------
[d,j] = max(CVA.v(:,1));
j     = j + (-8:8);
J     = spm_soup(Q(:,:,j),X(:,:,j),V(:,:,j),P);
J     = mean(J,3);
j     = [iXe;iXm];
% q   = 8;
% U   = blkdiag(CVA.V(:,1:q),CVA.W(:,1:q));  % eigenmodes (not used)

n     = 4;
for i = 1:n

    % all states
    %----------------------------------------------------------------------
    JJ    = J^i;
    subplot(n,2,(i - 1)*2 + 1)
    spy(abs(JJ) > 1e-2,'k')
    title(sprintf('%i-order coupling',i),'FontSize',16)
    xlabel('All states','FontSize',12)
    ylabel('All states','FontSize',12)

    % Internal and external states
    %----------------------------------------------------------------------
    JJ    = JJ(j,j);  % JJ     = pinv(U)*JJ(j,j)*U;
    subplot(n,2,(i - 1)*2 + 2)
    spy(abs(JJ) > 1e-2,'k')
    title(sprintf('%i-order coupling',i),'FontSize',16)
    xlabel('External and internal','FontSize',12)
    ylabel('External and internal','FontSize',12)

end

return

author: spm