text:
classdef Layer < handle
    %LAYER Base class for a network layer in a DagNN

    properties (Access = {?dagnn.DagNN, ?dagnn.Layer}, Hidden, Transient)
        net
        layerIndex
    end

    methods
        function outputs = forward(obj, inputs, params)
            %FORWARD Forward step
            %  OUTPUTS = FORWARD(OBJ, INPUTS, PARAMS) takes the layer object OBJ
            %  and cell arrays of inputs and parameters and produces a cell
            %  array of outputs evaluating the layer forward.
            outputs = {} ;
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutpus)
            %BACKWARD  Bacwkard step
            %  [DERINPUTS, DERPARAMS] = BACKWARD(OBJ, INPUTS, INPUTS, PARAMS,
            %  DEROUTPUTS) takes the layer object OBJ and cell arrays of
            %  inputs, parameters, and output derivatives and produces cell
            %  arrays of input and parameter derivatives evaluating the layer
            %  backward.
            derInputs = {} ;
            derOutputs = {} ;
        end

        function reset(obj)
            %RESET Restore internal state
            %  RESET(OBJ)resets the layer OBJ, clearing any internal state.
        end

        function params = initParams(obj)
            %INIT Initialize layer parameters
            %  PARAMS = INIT(OBJ) takes the layer OBJ and returns a cell
            %  array of layer parameters PARAMS with some initial
            %  (e.g. random) values.
            params = {} ;
        end

        function move(obj, device)
            %MOVE Move data to CPU or GPU
            %  MOVE(DESTINATION) moves the data associated to the layer object OBJ
            %  to either the 'gpu' or the 'cpu'. Note that variables and
            %  parameters are moved automatically by the DagNN object
            %  containing the layer, so this operation affects only data
            %  internal to the layer (e.g. the mask in dropout).
        end

        function forwardAdvanced(obj, layer)
            %FORWARDADVANCED  Advanced driver for forward computation
            %  FORWARDADVANCED(OBJ, LAYER) is the advanced interface to compute
            %  the forward step of the layer.
            %
            %  The advanced interface can be changed in order to extend DagNN
            %  non-trivially, or to optimise certain blocks.

            in = layer.inputIndexes ;
            out = layer.outputIndexes ;
            par = layer.paramIndexes ;
            net = obj.net ;

            inputs = {net.vars(in).value} ;

            % give up if any of the inputs is empty (this allows to run
            % subnetworks by specifying only some of the variables as input --
            % however it is somewhat dangerous as inputs could be legitimaly
            % empty)
            if any(cellfun(@isempty, inputs)), return ; end

            % clear inputs if not needed anymore
            for v = in
                net.numPendingVarRefs(v) = net.numPendingVarRefs(v) - 1 ;
                if net.numPendingVarRefs(v) == 0
                    if ~net.vars(v).precious & ~net.computingDerivative & net.conserveMemory
                        net.vars(v).value = [] ;
                    end
                end
            end

            %[net.vars(out).value] = deal([]) ;

            % call the simplified interface
            outputs = obj.forward(inputs, {net.params(par).value}) ;
            for oi = 1:numel(out)
                net.vars(out(oi)).value = outputs{oi};
            end
        end

        function backwardAdvanced(obj, layer)
            %BACKWARDADVANCED Advanced driver for backward computation
            %  BACKWARDADVANCED(OBJ, LAYER) is the advanced interface to compute
            %  the backward step of the layer.
            %
            %  The advanced interface can be changed in order to extend DagNN
            %  non-trivially, or to optimise certain blocks.
            in = layer.inputIndexes ;
            out = layer.outputIndexes ;
            par = layer.paramIndexes ;
            net = obj.net ;

            inputs = {net.vars(in).value} ;
            derOutputs = {net.vars(out).der} ;
            for i = 1:numel(derOutputs)
                if isempty(derOutputs{i}), return ; end
            end

            if net.conserveMemory
                % clear output variables (value and derivative)
                % unless precious
                for i = out
                    if net.vars(i).precious, continue ; end
                    net.vars(i).der = [] ;
                    net.vars(i).value = [] ;
                end
            end

            % compute derivatives of inputs and paramerters
            [derInputs, derParams] = obj.backward ...
                (inputs, {net.params(par).value}, derOutputs) ;
            if ~iscell(derInputs) || numel(derInputs) ~= numel(in)
                error('Invalid derivatives returned by layer "%s".', layer.name);
            end

            % accumuate derivatives
            for i = 1:numel(in)
                v = in(i) ;
                if net.numPendingVarRefs(v) == 0 || isempty(net.vars(v).der)
                    net.vars(v).der = derInputs{i} ;
                elseif ~isempty(derInputs{i})
                    net.vars(v).der = net.vars(v).der + derInputs{i} ;
                end
                net.numPendingVarRefs(v) = net.numPendingVarRefs(v) + 1 ;
            end

            for i = 1:numel(par)
                p = par(i) ;
                if (net.numPendingParamRefs(p) == 0 && ~net.accumulateParamDers) ...
                        || isempty(net.params(p).der)
                    net.params(p).der = derParams{i} ;
                else
                    net.params(p).der = vl_taccum(...
                        1, net.params(p).der, ...
                        1, derParams{i}) ;
                    %net.params(p).der = net.params(p).der + derParams{i} ;
                end
                net.numPendingParamRefs(p) = net.numPendingParamRefs(p) + 1 ;
                if net.numPendingParamRefs(p) == net.params(p).fanout
                    if ~isempty(net.parameterServer) && ~net.holdOn
                        net.parameterServer.pushWithIndex(p, net.params(p).der) ;
                        net.params(p).der = [] ;
                    end
                end
            end
        end

        function rfs = getReceptiveFields(obj)
            %GETRECEPTIVEFIELDS  Get receptive fields.
            %   RFS = GETRECEPTIVEFIELDS(OBJ) gets the receptive fields
            %   of each output varaibles in each input variable.
            %
            %   A *receptive field* is a structure with fields
            %
            %   - size: size of the receptive field
            %   - stride: stride of the receptive field
            %   - offset: offset of the receptive field
            %
            %   It should be interpreted as follows. Given a pixel of
            %   vertical coordinate u in an output variable OUT(y,...) , the first and last
            %   pixels affecting that pixel in an input variable IN(v,...) are:
            %
            %        v_first = stride(1) * (y - 1) + offset(1) - size(1)/2 + 1
            %        v_last  = stride(1) * (y - 1) + offset(1) + size(1)/2 + 1
            %
            %   RFS is a struct array of such structure, with one row for each
            %   input variable and one column for each output variable, expressing
            %   all possible combinations of inputs and outputs.
            rfs = [] ;
        end

        function outputSizes = getOutputSizes(obj, inputSizes)
            outputSizes = {} ;
        end

        function load(obj, varargin)
            %LOAD Initialize the layer from a paramter structure
            %  LOAD(OBJ, S) initializes the layer object OBJ from the parameter
            %  structure S.  It is the opposite of S = SAVE(OBJ).
            %
            %  LOAD(OBJ, OPT1, VAL1, OPT2, VAL2, ...) uses instead the
            %  option-value pairs to initialize the object properties.
            %
            %  LOAD(OBJ, {OPT1, VAL1, OPT2, VAL2, ...}) is an equivalent form
            %  to the previous call.
            s = dagnn.Layer.argsToStruct(varargin{:}) ;
            for f = fieldnames(s)'
                fc = char(f) ;
                if ~isprop(obj, fc)
                    error('No property `%s` for a layer of type `%s`.', ...
                        fc, class(obj));
                end;
                obj.(fc) = s.(fc) ;
            end
        end

        function s = save(obj)
            %SAVE Save the layer configuration to a parameter structure
            %  S = SAVE(OBJ) extracts all the properties of the layer object OBJ
            %  as a structure S. It is the oppostie of LOAD(OBJ, S).
            %
            %  By default, properties that are marked as transient,
            %  dependent, abstract, or private in the layer object are not
            %  saved.
            s = struct ;
            m = metaclass(obj) ;
            for p = m.PropertyList'
                if p.Transient || p.Dependent || p.Abstract, continue ; end
                s.(p.Name) = obj.(p.Name) ;
            end
        end

        function attach(obj, net, index)
            %ATTACH  Attach the layer to a DAG.
            %   ATTACH(OBJ, NET, INDEX). Override this function to
            %   configure parameters or variables.
            obj.net = net ;
            obj.layerIndex = index ;
            for input = net.layers(index).inputs
                net.addVar(char(input)) ;
            end
            for output = net.layers(index).outputs
                net.addVar(char(output)) ;
            end
            for param = net.layers(index).params
                net.addParam(char(param)) ;
                p = net.getParamIndex(char(param)) ;
            end
        end

    end % methods

    methods (Static)
        function s = argsToStruct(varargin)
            %ARGSTOSTRUCT  Convert varadic arguments to structure
            %  S = ARGSTOSTRCUT('opt1', val1, ....) converts the list of options
            %  into a structure.
            %
            %  S = ARGSTOSTRUCT({'opt1', val1, ...}) is equivalent to the
            %  previous call.
            %
            %  S = ARGSTOSTRUCT(S) where S is a structure returns S as is.

            if numel(varargin) == 1 && isstruct(varargin{1})
                s = varargin{1} ;
            else
                if numel(varargin) == 1 && iscell(varargin{1})
                    args = varargin{1} ;
                else
                    args = varargin ;
                end
                s = cell2struct(args(2:2:end),args(1:2:end),2) ;
            end
        end
    end

end % classdef

author: aimerykong


text:
% fkn_dbl_exp_fit - makes a double-exponential curve fit of F(k,t)/F(k,0)
%   [dy,y_fit]=fkn_dbl_exp_fit(p,t,Fr,k) where:
%
%   dy = the error in the fit
%   y_fit = the fitted values of F(k,t)/F(k,0)
%
%   p(1) = gamma2 (the fraction of component 2)
%   p(2) = D1 (the diffusion coefficient of component 1)
%   p(3) = D2 (the diffusion coefficient of component 2)
%   p(4) = gamma0 (the fraction of immobile molecules)
%   p(5:end) = amplitudes for each value of k (correcting for errors in
%       F(k,0))
%   t = the times for each frame
%   Fr = F(k,t)/F(k,0) (where F(k,t) is the Hankel transform of the
%       experimental data)
%   k = spatial frequencies

function [dy,y_fit]=fkn_dbl_exp_fit(p,t,Fr,k)

gamma2=p(1);
D1=p(2);
D2=p(3);
gamma0=p(4);
y_fit=zeros(size(Fr));
dy=zeros(size(Fr));

% Computes dy(k) from the curve fits to Fr for each value of k with gamma2,
% D1, D2 and gamma0 common for all k
for i=1:length(k)
     y_fit(i,:)=p(4+i)*((1-gamma2-gamma0)*exp(-4*pi^2*k(i)^2*D1*t)+...
         gamma2*exp(-4*pi^2*k(i)^2*D2*t)+gamma0);
     dy(i,:)=y_fit(i,:)-Fr(i,:);
end
dy=dy(:);   % Places all errors in a vector

author: Sable

repo: mcbench-benchmarks


text:
function script_create_MRCNN_VOC2007_2012

models_dir             = fullfile(pwd, 'models-exps');
net_def_file_directory = fullfile(pwd, 'model-defs');

% create the model directory
model_dir_dst          = fullfile(models_dir,    'MRCNN_VOC2007_2012');
mkdir_if_missing(model_dir_dst);

% path to the model mat file
model_mat_file  = fullfile(model_dir_dst, 'detection_model_svm.mat');

% path to the model caffe definition file
model_net_def_file = fullfile(model_dir_dst, 'deploy_svm.prototxt');
% source path to the model caffe definition file
model_net_def_file_scr = fullfile(net_def_file_directory, 'MRCNN_model_svm.prototxt');

assert(exist(model_net_def_file_scr,'file')>0);
copyfile(model_net_def_file_scr,model_net_def_file);
assert(exist(model_net_def_file,'file')>0);

region_dir_src              = {};
net_weigths_file_region_dst = {};
pooler_regions              = {};

% set the source directory of the region adaptation modules
region_dir_src{end+1} = fullfile(models_dir, 'vgg_R0010_voc2012_2007_EB_ZP');  % region  1
region_dir_src{end+1} = fullfile(models_dir, 'vgg_RHalf1_voc2012_2007_EB_ZP'); % region  2
region_dir_src{end+1} = fullfile(models_dir, 'vgg_RHalf2_voc2012_2007_EB_ZP'); % region  3
region_dir_src{end+1} = fullfile(models_dir, 'vgg_RHalf3_voc2012_2007_EB_ZP'); % region  4
region_dir_src{end+1} = fullfile(models_dir, 'vgg_RHalf4_voc2012_2007_EB_ZP'); % region  5
region_dir_src{end+1} = fullfile(models_dir, 'vgg_R0005_voc2012_2007_EB_ZP');  % region  6
region_dir_src{end+1} = fullfile(models_dir, 'vgg_R0308_voc2012_2007_EB_ZP');  % region  7
region_dir_src{end+1} = fullfile(models_dir, 'vgg_R0510_voc2012_2007_EB_ZP');  % region  8
region_dir_src{end+1} = fullfile(models_dir, 'vgg_R0815_voc2012_2007_EB_ZP');  % region  9
region_dir_src{end+1} = fullfile(models_dir, 'vgg_R1018_voc2012_2007_EB_ZP');  % region 10
num_regions = length(region_dir_src);

for r = 1:num_regions
    [~, region_name] = fileparts(region_dir_src{r});
    region_model_mat_file = fullfile(region_dir_src{r}, 'detection_model_softmax.mat');
    assert(exist(region_model_mat_file,'file')>0);

    ld = load(region_model_mat_file, 'model'); model_this = ld.model; clear ld;
    pooler_regions{r} = model_this.pooler;

    net_weigths_file_region_dst{r} = fullfile(model_dir_dst,sprintf('%s.caffemodel', region_name));

    net_def_file_region_dst_this1 = fullfile(net_def_file_directory, 'auxiliary_def_files',...
        sprintf('VGG_ILSVRC_16_layers_pascal_train_test_stream%d.prototxt',r));
    assert(exist(net_def_file_region_dst_this1, 'file')>0)

    caffe.set_mode_cpu();
    curr_dir = pwd;
    cd(region_dir_src{r});
    net_region_src         = caffe_load_model(model_this.net_def_file, model_this.net_weights_file);
    cd(curr_dir);
    net_region_dst         = caffe.Net(net_def_file_region_dst_this1, 'test');
    layers_region_src_this = {'fc6','fc7'};
    layers_region_dst_this = strcat(layers_region_src_this, sprintf('_s%d',r));
    layers_region_dst_this
    net_region_dst = caffe_copy_weights_from_net2net( net_region_dst, net_region_src, layers_region_dst_this, layers_region_src_this);
    net_region_dst.save(net_weigths_file_region_dst{r});
    caffe.reset_all();
end

% merge region types
pooler = pooler_regions{1};
for r = 1:num_regions
    fprintf('Region #%d - pooler: \n', r);
    disp(pooler_regions{r})
    pooler(r) = pooler_regions{r};
end

% convert paths from absolute to relative to the model directory
[a,b,c] = fileparts(model_net_def_file);
model_net_def_file = ['./',b,c];
for i = 1:length(net_weigths_file_region_dst)
    [a,b,c] = fileparts(net_weigths_file_region_dst{i});
    net_weigths_file_region_dst{i} = ['./',b,c];
end

% prepare and save the model structure of the multi-region cnn model

feat_blob_name   = {'fc_feat'}; % name of the output blob of the last hidden layer of the model
model_feat_cache = 'VGG_ILSVRC_16_layers'; % code-name of the activation maps that the VGG16 convolutional layers produce

model                  = struct;
model.net_def_file     = model_net_def_file;
model.net_weights_file = net_weigths_file_region_dst;
model.pooler           = pooler;
model.feat_blob_name   = feat_blob_name;
model.feat_cache       = model_feat_cache;
model.score_out_blob   = 'pascal_svm';
model.svm_layer_name   = 'pascal_svm';

VOCopts = initVOCOpts( '', '2007');
model.classes = VOCopts.classes;
fprintf('model:\n')
disp(model);
save(model_mat_file, 'model', '-v7.3');

end

author: gidariss

repo: mrcnn-object-detection


text:
% Reference performance
accTestRef = 0.773333;
allowedError = 0.001;
method = 'SVC1V1';

% Create the algorithm object
algorithmObj = SVC1V1();

% Clear parameter struct
clear param;

% Parameter C (Cost)
param.C = 10;

% Parameter k (kernel width)
param.k = 1;

% Run the algorithm
info = algorithmObj.fitpredict(train,test,param);

trainCM = confusionmat(info.predictedTrain,train.targets);
testCM = confusionmat(info.predictedTest,test.targets);

accTrain = CCR.calculateMetric(trainCM);
accTest  = CCR.calculateMetric(testCM);

% Report accuracy
fprintf('Performing test for %s\n', method);
fprintf('Accuracy Train %f, Accuracy Test %f\n',accTrain,accTest);

if abs(accTestRef-accTest)<allowedError
    fprintf('Test accuracy matches reference accuracy\n');
else
    warning('Test accuracy does NOT match reference accuracy');
end

author: ayrna

repo: orca

save_path: github-repos/MATLAB/ayrna-orca


text:
function showTreeHTML(tree)
% showTreeHTML(tree)
%
% tree = specifies the tree. It contains two fields:
%    tree.parent   this uses the Matlab structure for specifiying trees
%    (is a list of parents)
% tree.nodewords = is a cell array. Each element is the name associated
%    with each node.

pagename = 'test.html'

threshold = 0;

% Hearder
page = {};
page = addline(page, '<title>LabelMe tree</title>');
page = addline(page, '');
page = addline(page, '+ WordNet
');

%page = addline(page, ' Photoalbum
');
%page = addline(page, '

');
%page = addline(page, '');
% page = addline(page, drawlinehtml(10,10,40,30));
% page = addline(page, drawlinehtml(10,10,20,30));
% page = addline(page, drawlinehtml(10,10,50,30));
% page = addline(page, drawlinehtml(10,10,0,30));

page = addline(page, 'LabelMe by MIT.
');
page = addline(page, 'WordNet 2.1 Copyright 2005 by Princeton University.
');

page = addline(page, '
');
page = addline(page, 'Use the scroll bars to view different parts of the tree.
');
page = addline(page, 'If you use Firefox, use the CTRL+ and CTRL- to zoom in and out.
');
page = addline(page, 'Click on one word to see the object images.
');

p = tree.parent;
[x,y,h,s] = trimtreelayout(tree.parent);
x = round(x*1000);
y = round((1-y)*400);

x = x+rand(size(x))*2;
y = y+rand(size(y))*max(y)/50;
y = y-min(y(:))+1;

f = find(p~=0);
pp = p(f);

% figure
% hold on; axis('ij')
% for i = 1:length(f)
%     page = addline(page, drawlinehtml(x(f(i)), y(f(i)), x(pp(i)), y(pp(i))));
%     plot([x(f(i)) x(pp(i))], [y(f(i)) y(pp(i))])
% end

% figure
% trimtreeplot(tree.parent);
% hold on
for n = 1:length(x)
    if tree.counts(n)>threshold
        p =  tree.parent(n);
        if p>0
            page = addline(page, drawlinehtml(x(n), y(n), x(p), y(p)));
        end
    end
end

for n = 1:length(x)
    chain = tree.nodewords{n};
    w=getwords(chain); w = w{1};
    w(1) = upper(w(1));

    if tree.counts(n)>threshold
        page = addline(page, texthtml(x(n), y(n), chain));
        text (x(n), y(n), w, 'fontsize', 3)
    end
end

% Close web page

page = addline(page, '');

page = strrep(page, '  ', ' ');
page = strrep(page, '  ', ' ');

% write web page
fid = fopen(pagename,'w');
for i = 1:length(page);
    fprintf(fid, [page{i} '\n']);
end
fclose(fid);

function page = addline(page, line)
page = [page {line}];

function line = texthtml(x,y,s)
%line{1} = sprintf('
%s
', y-.5, double(x)-length(s)/4, s);

% Anyadir a la query los terminos de la rama. No usar solo la leave. Por
% ejemplo, en mouse, hacer la query con "mouse + device", al menos usar el
% siguente padre.

w = getwords(s);
w = w{1};
w(1) = upper(w(1));

sq = strrep(s, ',', '');
sq = strrep(sq, '-', ' ');

line = sprintf('
', y-.5, double(x)-length(w)/4);
line = [line sprintf('', sq, s)];
line = [line sprintf('%s
', w)];

function line = drawlinehtml(x1,y1,x2,y2)
cols = [1 2 4 8 16 32 64 128 256];

r = 128*abs(x2-x1)/abs(y2-y1);
bestratio = find(r-cols>0);
cols = [1 2 4 8 16 32 64 128 256 256];
if length(bestratio)==0
    bestratio = cols(1);
else
    bestratio = cols(bestratio(length(bestratio))+1);
end

if abs(x2-x1)<1
    x1 = x2+.2*sign(rand-.5);
    bestratio = 2;
end

if y2<y1
    tmp = x1; x1=x2; x2 = tmp;
    tmp = y1; y1=y2; y2 = tmp;
end
if x2>x1
    line = sprintf('
', y1, x1);
    line = [line sprintf('',  x2-x1, y2-y1, bestratio)];
    line = [line '
'];
else
    line = sprintf('
', y1, x2);
    line = [line sprintf('', (x1-x2), (y2-y1), bestratio)];
    line = [line '
'];
end

function [x,y,h,s] = trimtreelayout(parent,post)
%   TRIMTREELAYOUT A modified standard function TREELAYOUT.
%     Produces different heights for the leaves. They appear in
%     their respective layer instead of the deepest layer. The
%     format is the same as for TREELAYOUT, see below.
%
%   [x,y,h,s] = trimtreelayout(parent,post)
%
%
%   TREELAYOUT Lay out tree or forest.
%   [x,y,h,s] = treelayout(parent,post)
%       parent is the vector of parent pointers, with 0 for a root.
%       post is a postorder permutation on the tree nodes.
%       (If post is omitted we compute it here.)
%       x and y are vectors of coordinates in the unit square at which
%       to lay out the nodes of the tree to make a nice picture.
%       Optionally, h is the height of the tree and s is the
%       number of vertices in the top-level separator.
%
%   See also ETREE, TREEPLOT, ETREEPLOT, SYMBFACT.

%   Copyright 1984-2000 The MathWorks, Inc.
%   $Revision: 5.11 $  $Date: 2000/06/08 20:18:46 $
%
%   Modified by RG 01-Nov-1 to plot the leaves in the right layer
%

% This is based on the C code in sptrees.c by John Gilbert.
% Leaves are spaced evenly on the x axis, and internal
% nodes are centered over their descendant leaves with
% y coordinate proportional to height in the tree.

n = length(parent);

 pv = [];
 if (size(parent,1)>1), parent = parent(:)'; end
 if (nargin<2) & ~all(parent==0 | parent>(1:n))
     % This does not appear to be in the form generated by ETREE.
     if (any(parent>n) | any(parent<0) | any(parent~=floor(parent)) ...
	 | any(parent==[1:n]))
        error('Bad vector of parent pointers.');
     end
     [parent,pv] = fixparent(parent);
 end

if nargin < 2,

    % Create the adjacency matrix A of the given tree,
    % and get the postorder with another call to etree.

    j = find(parent);
    A = sparse (parent(j), j, 1, n, n);
    A = A + A' + speye(n,n);
    [ignore, post] = etree(A);
%    post
end;

% Add a dummy root node #n+1, and identify the leaves.

parent = rem(parent+n, n+1) + 1;  % change all 0s to n+1s
isaleaf = ones(1,n+1);
isaleaf(parent) = zeros(n,1);

% In postorder, compute heights and descendant leaf intervals.
% Space leaves evenly in x (in postorder).

xmin = n(1,ones(1,n+1)); % n+1 copies of n
xmax = zeros(1,n+1);
height = zeros(1,n+1);
nkids = zeros(1,n+1);
nleaves = 0;

for i = 1:n,
    node = post(i);
    if isaleaf(node),
        nleaves = nleaves+1;
        xmin(node) = nleaves;
        xmax(node) = nleaves;
    end;
    dad = parent(node);
    %RG
%    height(dad) = max (height(dad), height(node)+1);
    xmin(dad)   = min (xmin(dad),   xmin(node));
    xmax(dad)   = max (xmax(dad),   xmax(node));
    nkids(dad)  = nkids(dad)+1;
end;

% RG compute heights
% traverse the tree from the root downwards in a layer-manner
lay_ind = n+1;
par_ind = n+1;
while(1)
  lay_ind = find(ismember(parent, lay_ind));
  if isempty(lay_ind)
    break;
  end
  par_ind = [par_ind lay_ind];
  height(par_ind) = height(par_ind)+1;
end
height(1:n) = height(1:n)-1;

% Compute coordinates, leaving a little space on all sides.

treeht = height(n+1) - 1;
deltax = 1/(nleaves+1);
deltay = 1/(treeht+2);
x = deltax * (xmin+xmax)/2;
y = deltay * (height+1);

% Omit the dummy node.

x = x(1:n);
y = y(1:n);

% Return the height and top separator size.

h = treeht;
s = n+1 - max(find(nkids~=1));

if ~isempty(pv)
   x(pv) = x;
   y(pv) = y;
end

% ----------------------------
function [a,pv] = fixparent(parent)
%FIXPARENT  Fix order of parent vector
%   [A,PV] = FIXPARENT(B) takes a vector of parent nodes for an
%   elimination tree, and re-orders it to produce an equivalent vector
%   A in which parent nodes are always higher-numbered than child
%   nodes.  If B is an elimination tree produced by the TREE
%   functions, this step will not be necessary.  PV is a
%   permutation vector, so that A = B(PV);

n = length(parent);
a = parent;
a(a==0) = n+1;
pv = 1:n;

niter = 0;
while(1)
   k = find(a<(1:n));
   if isempty(k), break; end
   k = k(1);
   j = a(k);

   % Put node k before its parent node j
   a  = [ a(1:j-1)  a(k)  a(j:k-1)  a(k+1:end)];
   pv = [pv(1:j-1) pv(k) pv(j:k-1) pv(k+1:end)];
   t = (a >= j & a < k);
   a(a==k) = j;
   a(t) = a(t) + 1;

   niter = niter+1;
   if (niter>n*(n-1)/2), error('Bad vector of parent pointers.'); end
end

a(a>n) = 0;

function trimtreeplot(p,c,d)
% TRIMTREEPLOT A modified standard function TREEPLOT. Plots the
%   leaves differently from TREEPLOT. They appear in their
%   respective layer instead of the deepest layer so that the tree
%   appears as "trimmed". The format is the same as for TREEPLOT,
%   see below.
%
% trimtreeplot(p [,c,d])
%
% TREEPLOT Plot picture of tree.
%   TREEPLOT(p) plots a picture of a tree given a vector of
%   parent pointers, with p(i) == 0 for a root.
%
%   TREEPLOT(P,nodeSpec,edgeSpec) allows optional parameters nodeSpec
%   and edgeSpec to set the node or edge color, marker, and linestyle.
%   Use '' to omit one or both.
%
%   See also ETREE, TREELAYOUT, ETREEPLOT.

%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.8 $  $Date: 1997/11/21 23:44:59 $
%
%   Modified by RG 01-Nov-1 to use trimtreelayout instead of
%   TREELAYOUT

% RG use another treelayout
%[x,y,h]=treelayout(p);
[x,y,h]=trimtreelayout(p);

f = find(p~=0);
pp = p(f);
X = [x(f); x(pp); repmat(NaN,size(f))];
Y = [y(f); y(pp); repmat(NaN,size(f))];
X = X(:);
Y = Y(:);

if nargin == 1,
    n = length(p);
    if n < 500,
        plot (x, y, 'ro', X, Y, 'r-');
    else,
        plot (X, Y, 'r-');
    end;
else,
    [ignore, clen] = size(c);
    if nargin < 3,
        if clen > 1,
            d = [c(1:clen-1) '-'];
        else,
            d = 'r-';
        end;
    end;
    [ignore, dlen] = size(d);
    if clen>0 & dlen>0
        plot (x, y, c, X, Y, d);
    elseif clen>0,
        plot (x, y, c);
    elseif dlen>0,
        plot (X, Y, d);
    else
    end;
end;
xlabel(['height = ' int2str(h)]);
axis([0 1 0 1]);

author: CSAILVision

repo: LabelMeToolbox


text:
filename ='test2d_micro';%'RVE_Square_Triangle_FineFine';
ptype = 'MICRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
cost={'chomog_alphabeta','perimeterConstraint'};
weights=[1 0.1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';
initial_case = 'Vigdergauz';
vigdergauzType = 'VolumeAndRatio';
superEllipseRatio = 0.9;
volumeMicro = 0.61;

nsteps = 1;
Vfrac_final = 0.5;
Perimeter_target=1;
optimality_final =1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

%Micro
epsilon_isotropy_initial=1e-1;
epsilon_isotropy_final = 1e-3;
micro.alpha =[1 1 0]';
micro.beta =[1 1 0]';

% For all tests
plotting = false;
printing = false;
monitoring = false;
maxiter = 0;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%%Ex 2 if - elseif - else (This example is self-explanatory. Try to change the given value
% of num1 and observe the outcome.)

num1 = 4;
if (num1 >= 5)
 fprintf('%4i is greater than or equal to 5 \r', num1)
elseif (num1 > 1)
 fprintf('%4i is less than 5 but greater than 1 \r', num1)
elseif (num1 == 1)
 fprintf('%4i equals 1 \r', num1)
elseif (num1 > -3)
 fprintf('%4i is less than 1 but greater than -3 \r', num1)
else
 fprintf('%4i is less than or equal to -3 \r', num1)
end

author: TheAlgorithms

repo: MATLAB-Octave


text:
function o = mrAnatComputeMutualInfo(im1,im2,sampDensity,xform,cf,fwhm)
% Returns one of various mutual-info type measures of image similarity.
%
% o = mrAnatComputeMutualInfo(im1,im2,sampDensity,xform,cf,fwhm)
%
% NOTE: this code is copied from SPM5's spm_coreg.m! It was a sub-function
% in there, but we needed to call it separately. There are very trivial
% modifications in the calling convention and defaults. See spm_coreg for
% implementation notes and references.
%
% Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
% Written by John Ashburner.
%
% HISTORY:
% 2007.05.01: RFD- extracted code from spm_coreg.

if nargin<6, fwhm = [7 7];   end;
if nargin<5, cf   = 'nmi';    end;
if nargin<4, xform = eye(4);    end;
if nargin<3, sampDensity = [2 2 2];    end;

% Create the joint histogram
H = spm_hist2(im1, im2, xform, sampDensity);

% Smooth the histogram
lim  = ceil(2*fwhm);
krn1 = smoothing_kernel(fwhm(1),-lim(1):lim(1)) ; krn1 = krn1/sum(krn1); H = conv2(H,krn1);
krn2 = smoothing_kernel(fwhm(2),-lim(2):lim(2))'; krn2 = krn2/sum(krn2); H = conv2(H,krn2);

% Compute cost function from histogram
H  = H+eps;
sh = sum(H(:));
H  = H/sh;
s1 = sum(H,1);
s2 = sum(H,2);

switch lower(cf)
	case 'mi',
		% Mutual Information:
		H   = H.*log2(H./(s2*s1));
		mi  = sum(H(:));
		o   = -mi;
	case 'ecc',
		% Entropy Correlation Coefficient of:
		% Maes, Collignon, Vandermeulen, Marchal & Suetens (1997).
		% "Multimodality image registration by maximisation of mutual
		% information". IEEE Transactions on Medical Imaging 16(2):187-198
		H   = H.*log2(H./(s2*s1));
		mi  = sum(H(:));
		ecc = -2*mi/(sum(s1.*log2(s1))+sum(s2.*log2(s2)));
		o   = -ecc;
	case 'nmi',
		% Normalised Mutual Information of:
		% Studholme,  Hill & Hawkes (1998).
		% "A normalized entropy measure of 3-D medical image alignment".
		% in Proc. Medical Imaging 1998, vol. 3338, San Diego, CA, pp. 132-143.
		nmi = (sum(s1.*log2(s1))+sum(s2.*log2(s2)))/sum(sum(H.*log2(H)));
		o   = -nmi;
	case 'ncc',
		% Normalised Cross Correlation
		i     = 1:size(H,1);
		j     = 1:size(H,2);
		m1    = sum(s2.*i');
		m2    = sum(s1.*j);
		sig1  = sqrt(sum(s2.*(i'-m1).^2));
		sig2  = sqrt(sum(s1.*(j -m2).^2));
		[i,j] = ndgrid(i-m1,j-m2);
		ncc   = sum(sum(H.*i.*j))/(sig1*sig2);
		o     = -ncc;
	otherwise,
		error('Invalid cost function specified');
end;

return;

%_______________________________________________________________________
function krn = smoothing_kernel(fwhm,x)

% Variance from FWHM
s = (fwhm/sqrt(8*log(2)))^2+eps;

% The simple way to do it. Not good for small FWHM
% krn = (1/sqrt(2*pi*s))*exp(-(x.^2)/(2*s));

% For smoothing images, one should really convolve a Gaussian
% with a sinc function.  For smoothing histograms, the
% kernel should be a Gaussian convolved with the histogram
% basis function used. This function returns a Gaussian
% convolved with a triangular (1st degree B-spline) basis
% function.

% Gaussian convolved with 0th degree B-spline
% int(exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t= -0.5..0.5)
% w1  = 1/sqrt(2*s);
% krn = 0.5*(erf(w1*(x+0.5))-erf(w1*(x-0.5)));

% Gaussian convolved with 1st degree B-spline
%  int((1-t)*exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t= 0..1)
% +int((t+1)*exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t=-1..0)
w1  =  0.5*sqrt(2/s);
w2  = -0.5/s;
w3  = sqrt(s/2/pi);
krn = 0.5*(erf(w1*(x+1)).*(x+1) + erf(w1*(x-1)).*(x-1) - 2*erf(w1*x   ).* x)...
      +w3*(exp(w2*(x+1).^2)     + exp(w2*(x-1).^2)     - 2*exp(w2*x.^2));

krn(krn<0) = 0;
return;
%_______________________________________________________________________

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% SP_HCURL_ERROR: Evaluate the error in H(curl) norm.
%
%   [errhcurl, errl2, errcurl] = sp_hcurl_error (space, msh, u, uex, curluex);
%
% INPUT:
%
%    space:   object defining the space of discrete functions (see sp_multipatch)
%    msh:     object defining the domain partition and the quadrature rule (see msh_multipatch)
%    u:       vector of dof weights
%    uex:     function handle to evaluate the exact solution
%    curluex: function handle to evaluate the curl of the exact solution
%
% OUTPUT:
%
%     errhcurl: error in H(curl) norm
%     errl2:    error in L^2 norm
%     errcurl:  error of the curl in L^2 norm
%
% Copyright (C) 2015 Rafael Vazquez
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.

function [errhcurl, errl2, errcurl] = sp_hcurl_error (space, msh, u, uex, curluex)

  if (space.npatch ~= msh.npatch)
    error ('The number of patches does not coincide')
  end

  if (numel(u) ~= space.ndof)
    error ('Wrong size of the vector of degrees of freedom')
  end

  for iptc = 1:msh.npatch
    if (isempty (space.dofs_ornt))
      u_ptc = u(space.gnum{iptc});
    else
      u_ptc = u(space.gnum{iptc}) .* space.dofs_ornt{iptc}.';
    end
    [error_hcurl(iptc), error_l2(iptc), error_curl(iptc)] = ...
        sp_hcurl_error (space.sp_patch{iptc}, msh.msh_patch{iptc}, u_ptc, uex, curluex);
  end
  errl2 = sqrt (sum (error_l2 .* error_l2));
  errhcurl = sqrt (sum (error_hcurl .* error_hcurl));
  errcurl = sqrt (sum (error_curl .* error_curl));

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [cm] = yd2cm(yd)
% Convert length from yards to centimeters.
% Chad Greene 2012
cm = yd*91.44;


text:
function [ x, seed ] = reciprocal_sample ( a, b, seed )

%*****************************************************************************80
%
%% RECIPROCAL_SAMPLE samples the Reciprocal PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < A <= B.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ cdf, seed ] = r8_uniform_01 ( seed );

  x = b^cdf / a^( cdf - 1.0 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = franke ( x, y )

%*****************************************************************************80
%
%% FRANKE returns the value of the Franke function #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Richard Franke,
%    Scattered Data Interpolation: Tests of Some Methods,
%    Mathematics of Computation,
%    Volume 38, Number 157, January 1982, pages 181-200.
%
%  Parameters:
%
%    Input, real X, Y, the evalution points.
%
%    Output, real VALUE, the function values.
%
  value = ...
      0.75 * exp ( ...
        - ( ( 9.0 * x - 2.0 ).^2 ...
          + ( 9.0 * y - 2.0 ).^2 ) / 4.0 ) ...
    + 0.75 * exp ( ...
        - ( ( 9.0 * x + 1.0 ).^2 ) / 49.0 ...
          - ( 9.0 * y + 1.0 )      / 10.0 ) ...
    + 0.5  * exp ( ...
        - ( ( 9.0 * x - 7.0 ).^2 ...
          + ( 9.0 * y - 3.0 ).^2 ) / 4.0 ) ...
    - 0.2  * exp ( ...
          - ( 9.0 * x - 4.0 ).^2 ...
          - ( 9.0 * y - 7.0 ).^2 );

  return
end


text:
function obj = reparse_contiguous(obj, varargin)
% Re-construct list of contiguous voxels in an image based on in-image
% voxel coordinates.  Coordinates are taken from obj.volInfo.xyzlist.
%
% Results are saved in obj.volInfo.cluster.
%
% xyzlist can be generated from iimg_read_img, and is done automatically by
% object-oriented fMRI image classes (fmri_image, image_vector,
% statistic_image)
%
% If 'nonempty' is entered as an optional argument, will use only voxels
% that are non-zero, non-nan in all columns of obj.dat.
%
% :Usage:
% ::
%
%    obj = reparse_contiguous(obj, ['nonempty'])
%
% ..
%    Copyright tor wager, 2011
%
%    Programmers' notes:
%
%    Edited 1/27/13 by tor to use all columns when calculating 'nonempty'
%                also fixed bug - was not using 'nonempty' input in some
%                cases
% ..
%
% .cluster and .xyzlist should both always be length v in-mask voxels
% if 'nonempty' is entered, then .dat should be length v in-mask voxels too
%
% NOTES: this will only work if xyzlist in volInfo has only voxels for
% contiguous clusters. If the mask in volInfo contains the whole brain,
% use 'nonempty' to exclude empty data values when
% redefining clusters.

wh = true(obj.volInfo.n_inmask, 1); %size(obj.volInfo.cluster));   %obj.volInfo.wh_inmask;

if isempty(wh)

    disp('Warning: image_vector.reparse_contiguous: obj.volInfo.cluster is empty. Bad object...rebuilding.');

    % No clusters - rebuild cluster index
    obj.volInfo.cluster = spm_clusters(obj.volInfo.xyzlist')';

end

if isempty(obj.dat)

    disp('Warning: image_vector.reparse_contiguous: obj.dat is empty. Bad object...');
    disp('Stopping in debugger so you can check object.');
    keyboard

end

% restrict to voxels with actual data if desired
if any(strcmp(varargin, 'nonempty'))

    obj = replace_empty(obj,'voxels');  % changed to 'voxels' only. SG 2/23/18
    wh = all(obj.dat ~= 0, 2) & all(~isnan(obj.dat), 2); % this will expand wh to full in-mask dataset, .nvox

end

obj.volInfo(1).cluster = zeros(size(wh));

% .cluster can be either the size of a reduced, in-mask dataset after removing empties
% or the size of the full in-mask dataset that defined the image
% (volInfo.nvox).  We have to switch behavior according to which it is.

if size(obj.volInfo(1).cluster, 1) == size(obj.volInfo(1).xyzlist, 1)
    % if we have in-mask voxel list only in .cluster, .xyzlist

    if length(wh) ~= length(obj.volInfo(1).cluster)
        error('Bug: length(wh) does not match number of elements in volInfo.cluster.');
    end

    newcl = spm_clusters(obj.volInfo(1).xyzlist(wh, :)')';
    obj.volInfo(1).cluster = zeros(size(wh));
    obj.volInfo(1).cluster(wh) = newcl;

elseif size(obj.volInfo(1).xyzlist, 1) ~= length(wh)

    disp('obj.volInfo.xyzlist voxels:')
    size(obj.volInfo(1).xyzlist, 1)

    disp('Voxels in image list to re-parse:')
    length(wh)

    error('You have created an invalid image object somewhere along the line.');

else   % full in-mask -- assign to correct voxels

    obj.volInfo(1).cluster(wh) = spm_clusters(obj.volInfo(1).xyzlist(wh, :)')';

end

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [H, Ctour]= plot_scalp(mnt, w, varargin)
%PLOT_SCALP - Display a weight vector as scalp topography
%
%Description:
% This is the low level function for displaying a scalp topography.
% In many cases it is more wise to use one of the other plot_scalp*
% functions.
%
%Synopsis:
% H= plot_scalp(MNT, W, <OPT>)
%
%Input:
% MNT: An electrode montage, see mnt_setElectrodePositions
% W:   Vector to be displayed as scalp topography. The length of W must
%      concide with the length of MNT.clab or the number of non-NaN
%      entries of MNT.x, or OPT must include a field 'WClab'.
%      Warning: when you do not specify OPT.WClab you have to make sure
%      that the entries of W are matching with MNT.clab, or
%      MNT.clab(find(~isnan(MNT.x))).
% OPT: struct or property/value list of optional properties:
%  .CLim          - 'range', 'sym' (default), '0tomax', 'minto0',
%                   or [minVal maxVal]
%  .ScalePos      - Placement of the Colorbar 'horiz', 'vert' (default),
%                   or 'none'
%  .Contour       - Specifies at what Heights contour lines are drawn.
%                   If 'contour' is a vector, its entries define the
%                   Heights. If is a scalar it specifies
%                   (according to 'ContourPolicy') the number
%                   of or the spacing between contour levels. To display
%                   no contour lines set 'contour' to 0 (not []!).
%  .ContourPolicy - 'levels' (default): 'contour' specifies exactly the
%                   number of contour levels to be drawn, or
%                   'spacing': 'contour' specifies the spacing between two
%                   adjacent Height levels, or
%                   'choose': '.contour' specifies approximately the
%                   number of Height levels to be drawn, but the function
%                   'goodContourValues' is called to find nice values.
%  .Resolution    - default 40. Number of steps around circle used for
%                   plotting the scalp.
%  .ShowLabels    - Display channel names (1) or not (0), default 0.
%  .Shading         shading method for the pColor plot, default 'flat'.
%                   Use 'interp' to get nice, smooth plots. But saving
%                   needs more space.
%  .Extrapolation - Default value (1) extends the scalp plot to the
%                   peripheral areas where no channels are located.
%                   Value (0) turns off extrapolation.
%  .ExtrapolateToMean - Default value (1) paints peripheral area
%                   in Color of average (zero?) value. Needs .Extrapolation
%                   activated.
%  .ExtrapolateToZero - Value (1) paints peripheral area in "zero"-Color.
%                   Needs .Extrapolation activated.
%  .Renderer      - The function used for rendering the scalp map, 'pColor'
%                   or 'contourf' (default).
%  .ContourfLevels - number of levels for contourf function (default 100).
%  .Offset        - a vector of length 2  -  [x_Offset y_Offset]
%                   normally, the scalpplot is drawn centered at the origin,
%                   i.e. [x_Offset y_Offset] = [0 0] by default
%
%Output:
% H:     handle to several graphical objects
% Ctour: struct of contour information
%
%See also plot_scalpPatterns, plot_scalpEvolution.

% Author(s): Benjamin Blankertz, Aug 2000; Feb 2005, Matthias
% Added contourf: Matthias Treder 2010
% Added "ExtrapolateToZero" option: Simon Scholler, 2011
% "Offset" option added: Simon Scholler, 2011

props= {
        'CLim',                 'sym',             'CHAR|DOUBLE[2]';
        'Colormap',              get(gcf,'Colormap'), 'DOUBLE[- 3]';
        'ContourfLevels',        50,               'DOUBLE[1]';
        'ContourMargin',         0,                'DOUBLE[1]';
        'Contour',               5,                'DOUBLE';
        'ContourPolicy',         'levels',         'CHAR';
        'ContourLineprop',       {'LineWidth',1},  'PROPLIST';
        'ContourLabels',         0,                'BOOL';
        'Extrapolation',         1,                'BOOL';
        'ExtrapolateToMean',     0,                'BOOL';
        'ExtrapolateToZero',     0,                'BOOL';
        'Interpolation',         'linear',         'CHAR';
        'LineProperties',        {'k'},            'CELL';
        'MarkContour',           [],               'DOUBLE[0-2]';
        'MarkContourLineprop',   {'LineWidth',2},  'PROPLIST';
        'NewColormap',           0,                'BOOL';
        'Offset',                [0 0],            'DOUBLE[2]';
        'Resolution',            51,               'DOUBLE[1]';
        'Renderer',              'contourf',       'CHAR';
        'Shading',               'flat',           'CHAR';
        'ShrinkColorbar',        0,                'DOUBLE';
        'ScalePos',              'vert',           'CHAR';
        'TicksAtContourLevels',  1,                'BOOL';
        'WClab',                 {},               'CELL{CHAR}';
        };

props_scalpOutline = plot_scalpOutline;

if nargin==0,
  H= opt_catProps(props, props_scalpOutline); return
end

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props, props_scalpOutline);

opt_scalpOutline= opt_substruct(opt, props_scalpOutline(:,1));

if opt.NewColormap,
  acm= visutil_addColormap(opt.Colormap);
elseif isfield(opt, 'Colormap'),
  colormap(opt.Colormap);
end

if opt.Extrapolation,
  if isdefault.LineProperties,
    opt.LineProperties= {'Color','k', 'LineWidth',3};
  end
  if isdefault.Resolution,
    opt.Resolution= 101;
  end
end

w= w(:);
if ~isempty(opt.WClab),
  mnt= mnt_adaptMontage(mnt, opt.WClab);
  if length(mnt.clab)<length(opt.WClab),
    error('some channels of opt.WClab not found in montage');
  end
end
if length(w)==length(mnt.clab),
  DisplayChannels= find(~isnan(mnt.x));
  w= w(DisplayChannels);
else
  DisplayChannels= find(~isnan(mnt.x));
  if length(w)~=length(DisplayChannels),
    error(['length of w must match # of displayable channels, ' ...
           'i.e. ~isnan(mnt.x), in mnt']);
  end
end
xe= mnt.x(DisplayChannels);
ye= mnt.y(DisplayChannels);

H.ax= gca;
if isempty(opt.Resolution)
  oldUnits= get(H.ax, 'units');
  set(H.ax, 'units', 'normalized');
  pos= get(H.ax, 'position');
  set(H.ax, 'units', oldUnits);
  opt.Resolution= max(20, 60*pos(3));
end

% Allow radius of scalp data to go beyond scalp outline (>1)
maxrad = max(1,max(max(abs(mnt.x)),max(abs(mnt.y)))) + opt.ContourMargin;

% Extrapolation
if opt.Extrapolation,
  xx= linspace(-maxrad, maxrad, opt.Resolution);
  yy= linspace(-maxrad, maxrad, opt.Resolution)';
  if opt.ExtrapolateToMean
    xe_add = cos(linspace(0,2*pi,opt.Resolution))'*maxrad;
    ye_add = sin(linspace(0,2*pi,opt.Resolution))'*maxrad;
    w_add = ones(length(xe_add),1)*mean(w);
    xe = [xe;xe_add];
    ye = [ye;ye_add];
    w = [w;w_add];
  end
  if opt.ExtrapolateToZero
    xe_add = cos(linspace(0,2*pi,opt.Resolution))';
    ye_add = sin(linspace(0,2*pi,opt.Resolution))';
    xe = [xe;xe_add];
    ye = [ye;ye_add];
    w = [w; zeros(length(xe_add),1)];
  end

else
  xx= linspace(min(xe), max(xe), opt.Resolution);
  yy= linspace(min(ye), max(ye), opt.Resolution)';
end

if opt.Extrapolation,
  wstate= warning('off');
  [xg,yg,zg]= griddata(xe, ye, w, xx, yy, 'v4');
  warning(wstate);
  margin = maxrad +opt.ContourMargin;
  headmask= (sqrt(xg.^2+yg.^2)<=margin);
  imaskout= ~headmask;
  zg(imaskout)= NaN;

else
  if strcmp(opt.Interpolation, 'v4'),
    % get the convex hull from linear Interpolation
    [dummy1,dummy2,zconv]= griddata(xe, ye, w, xx, yy, 'linear');
    imaskout= isnan(zconv(:));
    [xg,yg,zg]= griddata(xe, ye, w, xx, yy, opt.Interpolation);
    zg(imaskout)= NaN;
  else
    [xg,yg,zg]= griddata(xe, ye, w, xx, yy, opt.Interpolation);
  end
end

xs= xg(1,2)-xg(1,1);
ys= yg(2,1)-yg(1,1);

% contour line coordinated
xgc= xg+opt.Offset(1);
ygc= yg+opt.Offset(2);
zgc= zg;
if ~isempty(strmatch(lower(opt.Shading), {'flat','faceted'})) ...
    && strcmp(opt.Renderer,'pColor'),
  % in shading FLAT last row/column is skipped, so add one
  xg= [xg-xs/2, xg(:,end)+xs/2];
  xg= [xg; xg(end,:)];
  yg= [yg, yg(:,end)]-ys/2;
  yg= [yg; yg(end,:)+ys];
  zg= [zg, zg(:,end)];
  zg= [zg; zg(end,:)];
end

% Render using pColor or contourf
xg= xg+opt.Offset(1);
yg= yg+opt.Offset(2);
if strcmp(opt.Renderer,'pColor')
  H.patch= pColor(xg, yg, zg);
else
  [dummy,H.patch]= contourf(xg, yg, zg, opt.ContourfLevels,'LineStyle','none');
  % *** Hack to enforce cdatamappig = scaled in Colorbarv6.m by introducing
  % a useless patch object
  hold on
  patch([0 0],[0 0],[1 2]);
  ccc = get(gca,'children');
  set(ccc(1),'Visible','off');
end

%
tight_caxis= [min(zg(:)) max(zg(:))];
if isequal(opt.CLim, 'sym'),
  zgMax= max(abs(tight_caxis));
  H.CLim= [-zgMax zgMax];
elseif isequal(opt.CLim, 'range'),
  H.CLim= tight_caxis;
elseif isequal(opt.CLim, '0tomax'),
  H.CLim= [0.0001*diff(tight_caxis) max(tight_caxis)];
elseif isequal(opt.CLim, 'minto0'),
  H.CLim= [min(tight_caxis) -0.0001*diff(tight_caxis)];
elseif isequal(opt.CLim, 'zerotomax'),
  H.CLim= [0 max(tight_caxis)];
elseif isequal(opt.CLim, 'mintozero'),
  H.CLim= [min(tight_caxis) 0];
else
  H.CLim= opt.CLim;
end
if diff(H.CLim)==0, H.CLim(2)= H.CLim(2)+eps; end
set(gca, 'CLim',H.CLim);

if strcmp(opt.Renderer,'pColor')
  shading(opt.Shading);
end

hold on;
if ~isequal(opt.Contour,0),
  if length(opt.Contour)>1,
    ctick= opt.Contour;
    v= ctick;
  else
    mi= min(H.CLim);
    ma= max(H.CLim);
    switch(opt.ContourPolicy),
     case {'levels','strict'},
      ctick= linspace(mi, ma, abs(opt.Contour)+2);
      v= ctick([2:end-1]);
     case 'withinrange',
      ctick= linspace(min(tight_caxis), max(tight_caxis), abs(opt.Contour)+2);
      v= ctick([2:end-1]);
     case 'spacing',
      mm= max(abs([mi ma]));
      v= 0:opt.Contour:mm;
      v= [fliplr(-v(2:end)), v];
      ctick= v(v>=mi & v<=ma);
      v(v<=mi | v>=ma)= [];
     case 'spacing_compatability',
      v= floor(mi):opt.Contour:ceil(ma);
      ctick= v(v>=mi & v<=ma);
      v(v<=mi | v>=ma)= [];
     case 'choose',
      ctick= visutil_goodContourValues(mi, ma, -abs(opt.Contour));
      v= ctick;
     otherwise
      error('ContourPolicy not known');
    end
  end
  v_tmp= v;
  if length(v)==1,
    v_tmp= [v v];
  end
  if isempty(v),
    H.contour= [];
  else
    [c,H.contour]= contour(xgc, ygc, zgc, v_tmp, 'k-');
    set(H.contour, opt.ContourLineprop{:});
    if opt.ContourLabels, %% & length(H.contour)>1,
      clabel(c,H.contour);
    end
  end
  H.contour_Heights= v;
else
  H.contour_Heights = [];   %% recently added
end
if ~isempty(opt.MarkContour),
  if length(opt.MarkContour)==1,
    v_tmp= [1 1]*opt.MarkContour;
  else
    v_tmp= opt.MarkContour;
  end
  [c,H.MarkContour]= contour(xgc, ygc, zgc, v_tmp, 'k-');
  set(H.MarkContour, opt.markContourLineprop{:});
end

% Scalp outline
H= plot_scalpOutline(mnt, opt_scalpOutline, 'H',H, 'DisplayChannels',DisplayChannels);

if strcmp(opt.ScalePos, 'none'),
  H.cb= [];
else
  H.cb= colorbar(opt.ScalePos);
  if opt.TicksAtContourLevels && opt.Contour,
    if strcmpi(opt.ScalePos, 'vert'),
      set(H.cb, 'yLim',H.CLim, 'yTick', ctick);
      if opt.ShrinkColorbar>0,
        cbpos= get(H.cb, 'Position');
        cbpos(2)= cbpos(2) + cbpos(4)*opt.ShrinkColorbar/2;
        cbpos(4)= cbpos(4) - cbpos(4)*opt.ShrinkColorbar;
        set(H.cb, 'Position',cbpos);
      end
    else
      set(H.cb, 'xLim',H.CLim, 'xTick', ctick);
      if opt.ShrinkColorbar>0,
        cbpos= get(H.cb, 'Position');
        cbpos(1)= cbpos(1) + cbpos(1)*opt.ShrinkColorbar/2;
        cbpos(3)= cbpos(3) - cbpos(3)*opt.ShrinkColorbar;
        set(H.cb, 'Position',cbpos);
      end
    end
  end
  if ~isempty(H.contour_Heights) && opt.TicksAtContourLevels,
    set(H.cb, 'YTick',H.contour_Heights);
  end
end
if opt.NewColormap,
  visutil_acmAdaptCLim(acm);
  set(H.cb, 'yLim',H.CLim); % otherwise ticks at the border of the
                            % Colorbar might get lost
end
axis('off');

if nargout==0,
  clear H;
end
if nargout>=2,
  if ~exist('c', 'var'),
    c= [];
  end
  Ctour= struct('xgrid',xg, 'ygrid',yg, 'zgrid',zg, ...
                'values',v, 'matrix',c);
end

author: bbci


text:
function x = sarkis_systematic_polar_encoder(u, info_bits, frozen_bits, N, lambda_offset, llr_layer_vec)
x = zeros(N, 1);
%First step encoding
x(info_bits) = u;
x = polar_encoder(x, lambda_offset, llr_layer_vec);
%Second step encoding
x(frozen_bits) = 0;
x = polar_encoder(x, lambda_offset, llr_layer_vec);
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
function [y_img,dice]=svm_test(img_dil,Out,s_test,ind)

% INPUT
% img_dil:3-D image (matrix) comprises region of interest and its uncertainty zone
% Out: classifier computed from svm_train.m
% s_test: similarity matrix of new subject in testing samples
% ind: index of all subcortical voxels in the 3-D atlas image (MNI152 space)

% OUTPUT
% y_img: probabilistic map
% dice: dice coefficient

Nxyz=size(img_dil);

Mdl=Out.Mdl;
ScoreMdl=Out.ScoreMdl;
y_gt=Out.y_gt;

ind1=find(img_dil(ind)==1); %class 1: other regions
ind2=find(img_dil(ind)==2); %class 2: region
ind_both=[ind1;ind2];
ind_out=setdiff(1:length(ind),ind_both);

s_test=s_test(ind_both,ind_out); %voxels x features x subjects
y_img=zeros(Nxyz);

y_pred=predict(Mdl,s_test);
[~,tmp]=predict(ScoreMdl,s_test);

%y_pred_score=tmp(:,2);
tmp_img=zeros(Nxyz);
tmp_img(ind(ind_both))=tmp(:,2);
y_img(:,:,:)=tmp_img;

%Dice coefficient
common=(y_pred==2 & y_gt==2);
a = sum(common(:));
b = sum(y_pred==2);
c = sum(y_gt==2);
dice=2*a/(b+c);

author: yetianmed

repo: subcortex

save_path: github-repos/MATLAB/yetianmed-subcortex


function statusAll = runMatlab2TikzTests(varargin)
%% This file runs the complete MATLAB2TIKZ test suite.
% It is mainly used for testing on a continuous integration server, but it can
% also be used on a development machine.

CI_MODE = strcmpi(getenv('CONTINUOUS_INTEGRATION'),'true') || strcmp(getenv('CI'),'true');
isJenkins = ~isempty(getenv('JENKINS_URL'));

%% Set path
addpath(fullfile(pwd,'..','src'));
addpath(fullfile(pwd,'suites'));

%% Select functions to run
suite = @ACID;
allTests = 1:numel(suite(0));

%% Prepare environment
if strcmpi(getEnvironment(), 'Octave')
  % Ensure that paging is disabled
  % https://www.gnu.org/software/octave/doc/interpreter/Paging-Screen-Output.html
  more off
end

%% Run tests
status = testHeadless('testFunctionIndices', allTests,...
                     'testsuite',           suite, varargin{:});

if isJenkins
    makeTapReport(status, 'stream', 'results.test.tap');
    makeTravisReport(status, 'stream', 'results.test.md');
end

nErrors = makeTravisReport(status);

%% Calculate exit code
if CI_MODE
    exit(nErrors);
end


text:
function varargout = binaryLaplaceGP(hyper, covfunc, lik, x, y, xstar)

% binaryLaplaceGP - Laplace's approximation for binary Gaussian process
% classification. Two modes are possible: training or testing: if no test
% cases are supplied, then the approximate negative log marginal likelihood
% and its partial derivatives wrt the hyperparameters is computed; this mode is
% used to fit the hyperparameters. If test cases are given, then the test set
% predictive probabilities are returned. The program is flexible in allowing
% several different likelihood functions and a multitude of covariance
% functions.
%
% usage: [nlZ, dnlZ     ] = binaryLaplaceGP(hyper, covfunc, lik, x, y);
%    or: [p, mu, s2, nlZ] = binaryLaplaceGP(hyper, covfunc, lik, x, y, xstar);
%
% where:
%
%   hyper    is a (column) vector of hyperparameters
%   covfunc  is the name of the covariance function (see below)
%   lik      is the name of the likelihood function (see below)
%   x        is a n by D matrix of training inputs
%   y        is a (column) vector (of size n) of binary +1/-1 targets
%   xstar    is a nn by D matrix of test inputs
%   nlZ      is the returned value of the negative log marginal likelihood
%   dnlZ     is a (column) vector of partial derivatives of the negative
%               log marginal likelihood wrt each log hyperparameter
%   p        is a (column) vector (of length nn) of predictive probabilities
%   mu       is a (column) vector (of length nn) of predictive latent means
%   s2       is a (column) vector (of length nn) of predictive latent variances
%
% The length of the vector of log hyperparameters depends on the covariance
% function, as specified by the "covfunc" input to the function, specifying the
% name of a covariance function. A number of different covariance function are
% implemented, and it is not difficult to add new ones. See "help covFunctions"
% for the details.
%
% The shape of the likelihood function is given by the "lik" input to the
% function, specifying the name of the likelihood function. The two implemented
% likelihood functions are:
%
%   logistic      the logistic function: 1/(1+exp(-x))
%   cumGauss      the cumulative Gaussian (error function)
%
% The function can conveniently be used with the "minimize" function to train
% a Gaussian process, eg:
%
% [hyper, fX, i] = minimize(hyper, 'binaryLaplaceGP', length, 'covSEiso',
%                                                             'logistic', x, y);
%
% Copyright (c) 2004, 2005, 2006, 2007 Carl Edward Rasmussen, 2010-07-21.

if nargin<5 || nargin>6
  disp('Usage: [nlZ, dnlZ     ] = binaryLaplaceGP(hyper, covfunc, lik, x, y);')
  disp('   or: [p, mu, s2, nlZ] = binaryLaplaceGP(hyper, covfunc, lik, x, y, xstar);')
  return
end

% Note, this function is just a wrapper provided for backward compatibility,
% the functionality is now provided by the more general gp function.
if strcmp(lik,'logistic')
  lik = @likLogistic;
elseif strcmpi(lik,'cumgauss')
  lik = @likErf;
else
  error('Allowable likelihoods: logistic and cumGauss.');
end
inf = @infLaplace; mean = @meanZero;

varargout = cell(nargout, 1);    % allocate the right number of output arguments
hyp.cov = hyper;
if nargin==5
  [varargout{:}] = gp(hyp,inf,mean,covfunc,lik,x,y);
  if nargout>1, varargout{2} = varargout{2}.cov; end
else
  [ymu,ys2,fmu,fs2] = gp(hyp,inf,mean,covfunc,lik,x,y,xstar);
  if nargout>0, varargout{1} = (1+ymu)/2; end
  if nargout>1, varargout{2} = fmu; end
  if nargout>2, varargout{3} = fs2; end
  if nargout>3, varargout{4} = gp(hyp,inf,mean,covfunc,lik,x,y); end       % nlZ
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function i4_max_test ( )

%*****************************************************************************80
%
%% I4_MAX_TEST tests I4_MAX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_MAX_TEST\n' );
  fprintf ( 1, '  I4_MAX returns the maximum of two I4''s.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '       A       B      C=I4_MAX(A,B)\n' );
  fprintf ( 1, '\n' );

  i4_lo = -100.0;
  i4_hi = +100.0;
  seed = 123456789;
  for i = 1 : 10
    [ a, seed ] = i4_uniform_ab ( i4_lo, i4_hi, seed );
    [ b, seed ] = i4_uniform_ab ( i4_lo, i4_hi, seed );
    c = i4_max ( a, b );
    fprintf ( 1, '  %8d  %8d  %8d\n', a, b, c );
  end

  return
end


text:
function h = mod(x, y)
%MOD   Modulus after division of two CHEBFUN objects.
%   MOD(X, Y) is X - n.*Y where n = floor(X./Y).
%
% See also FLOOR.

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

n = floor(x./y);
h = x - n.*y;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function h = times(f, g)
%.*   CHEBFUN3T pointwise multiplication.
%   F.*G multiplies CHEBFUN3T objects F and G. Alternatively F or G could be
%   a double.

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

if ( isa(f, 'chebfun3t') )     % CHEBFUN3T .* ???

    if ( isa(g, 'double') )    % CHEBFUN3T .* DOUBLE
        h = mtimes(f, g);
    elseif ( isa(g, 'chebfun3t') )
        sameDom = domainCheck(f, g);
        if ( sameDom )
            h = chebfun3t(@(x, y, z) feval(f, x, y, z).*feval(g, x, y, z), ...
                f.domain);
        else
           error('CHEBFUN:CHEBFUN3T:times:domain', 'Inconsistent domains');
       end
    else
        error('CHEBFUN:CHEBFUN3T:times:unknown', ...
            ['Undefined function ''times'' for input arguments of type %s ' ...
            'and %s.'], class(f), class(g));
    end

else
   h = times(g, f);
end

end


text:
function geometry_test0782 ( )

%*****************************************************************************80
%
%% TEST0782 tests POLYGON_AREA_3D_2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;
  test_num = 2;
  n_test = [ 4, 8 ];
  area_exact_test = [ 2.4494898, 6.0 ];

  v = [ ...
    1.0, 0.0, 0.0; ...
    2.0, 1.0, 1.0; ...
    1.0, 2.0, 1.0; ...
    0.0, 1.0, 0.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0782\n' );
  fprintf ( 1, '  For a polygon in 3D:\n' );
  fprintf ( 1, '  POLYGON_AREA_3D_2 computes the area;\n' );

  for test = 1 : test_num

    n = n_test(test);
    area_exact = area_exact_test(test);

    if ( test == 1 )

      v = [ ...
        1.0, 0.0, 0.0; ...
        2.0, 1.0, 1.0; ...
        1.0, 2.0, 1.0; ...
        0.0, 1.0, 0.0 ]';

    elseif ( test == 2 )

      v = [ ...
        0.00000,  0.00000,  0.00000; ...
        2.62679,  1.26009, -0.715657; ...
        1.48153,  3.97300, -0.142512; ...
        0.605932, 3.55297,  0.0960401; ...
        1.36944,  1.74437, -0.286056; ...
        0.493842, 1.32433, -0.0475041; ...
        0.112090, 2.22864,  0.143544; ...
       -0.763505, 1.80861,  0.382097 ]';

    end

    r8mat_transpose_print ( dim_num, n, v, '  The polygon vertices:' )

    area = polygon_area_3d_2 ( n, v );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Exact area is        %f\n', area_exact );
    fprintf ( 1, '  The computed area is %f\n', area );

  end

  return
end


text:
function splitNSxNTrode

% splitNSxNTrode
%
% Opens and splits an NSx file based on its NTrode groups. It depends on
% openCCF file.
%
% Use splitNSxNTrode
%
% This function does not take any inputs.
%
%   Example 1:
%   splitNSxNTrode;
%
%   In the example above, the user will be prompted to select a CCF file
%   first. The CCF contains the ntrode grouping infromation. Then the user
%   will be prompted to select a NSx file. The script will then split the
%   NSx file into smaller NSx files containing channels in given ntrode
%   groups. For example, if ntrode group one consists of channels 1,3,5,
%   and 12, then using splitNSxNTrode will split the file into a smaller
%   NSx file that contains those channels only. If there are multiple
%   ntrodes then the files will split into multiple smaller files, equal in
%   number of the ntrode groups.
%
%   Kian Torab
%   ktorab@blackrockmicro.com
%   Blackrock Microsystems
%   Version 1.0.1.0
%

% Validating input parameter
ccf = openCCF;
splitCount = length(ccf.NTrodeInfo.NTrodeID);

% Getting the file name
if ~ismac
    [fname, path] = getFile('*.ns*', 'Choose an NSx file...');
else
    [fname, path] = getFile('*.*', 'Choose an NSx file...');
end
if fname == 0
    disp('No file was selected.');
    if nargout
        clear variables;
    end
    return;
end
fext = fname(end-3:end);

% Loading the file
%% Reading Basic Header from file into NSx structure.
FID                       = fopen([path fname], 'r', 'ieee-le');
NSx.MetaTags.Filename     = fname;
NSx.MetaTags.FilePath     = path(1:end-1);
NSx.MetaTags.FileExt      = fext;
NSx.MetaTags.FileTypeID   = fread(FID, [1,8]   , '*char');
disp(['Splitting the NSx file in ' num2str(splitCount) ' pieces...']);
if strcmpi(NSx.MetaTags.FileTypeID, 'NEURALSG')
    disp('File type 2.1 is not yet implemented.');
    %NOT IMPLEMENTED YET
%     fseek(FID, 0, 'bof');
%     header = fread(FID, 314,'*uint8');
%     positionEOH = ftell(FID);
%     fseek(FID, 0, 'eof');
%     positionEOD = ftell(FID);
%     dataLength = positionEOD - positionEOH;
%     fseek(FID, 28, 'bof');
%     channelCount = fread(FID, 1       , 'uint32=>double');
elseif strcmpi(NSx.MetaTags.FileTypeID, 'NEURALCD')
    % Calculating different points in the file
    fseek(FID, 0, 'bof');
    basicHeader = fread(FID, 314, '*uint8');
    positionEOB = ftell(FID);
    positionEOE = typecast(basicHeader(11:14), 'uint32');
    fseek(FID, 0, 'eof');
    positionEOD = ftell(FID);

    % Calculating channelCount, data Length
    channelCount = typecast(basicHeader(311:314), 'uint32');
    dataLength = positionEOD - positionEOE - 9;

    % Reading extended header and calculating channel IDs
    fseek(FID, positionEOB, 'bof');
    extHeader = fread(FID, [(positionEOE-positionEOB)/channelCount, channelCount], '*uint8');

    % Reading the channel IDs. This only wokrs for as long as channel IDs
    % are 8-bit integers. If they change, a typecast of 3:4 is necessasry
    channelID = extHeader(3,1:channelCount);

    % Reading the number of packets
    fseek(FID, 28, 'bof');

    % Calculating the number of bytes in each segment
    channelBytes = (dataLength)/channelCount;

    % Reading the headers and the data header
    fseek(FID, 0, 'bof');
    fileHeader = fread(FID, positionEOE, 'char');
    dataHeader = fread(FID, 9, 'char');
	fseek(FID, positionEOE+9, 'bof');

    % Reading the data
    fprintf('\nReading the entire data file...\n');
    dataSegment = fread(FID, [channelCount, channelBytes], 'int16');

    for idx = 1:splitCount
        % Determining whether tetrode channel is recorded and valid in NSx
        tetrodeChannels = ccf.NTrodeInfo.NTrodeMembers{idx};
        for tetIDX = 1:length(tetrodeChannels)
            validChannel = find(channelID == tetrodeChannels(tetIDX));
            if isempty(validChannel)
                fprintf(2,'The tetrode channel %1.0f from tetrode group %d does not exist in the continuous file.\n',tetrodeChannels(tetIDX), idx);
                break;
            end
            validChannels(tetIDX) = validChannel;
        end
        if ~isempty(validChannel)
            % Opening a file for saving
            fprintf('Writing segment %d...\n', idx);
            FIDw = fopen([path fname(1:end-4) '-tet' sprintf('%03d', idx) fname(end-3:end)], 'w+', 'ieee-le');

            % Writing the segmented data into file
            basicHeader(end-3) = length(validChannels);
            fwrite(FIDw, basicHeader, 'char');
            fwrite(FIDw, extHeader(:,validChannels), 'char');
            fwrite(FIDw, dataHeader, 'char');
            fwrite(FIDw, dataSegment(validChannels,:), 'int16');
            fclose(FIDw);
        end
    end
    fprintf('\n');
else
    % Display error if non-compatible file is trying to open.
    disp('This version of splitNSx can only split File Specs 2.2 and 2.3');
    disp(['The selected file spec is ' NSx.MetaTags.FileSpec '.']);
    fclose(FID);
    clear variables;
    return;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [xtraj, ttraj, terminate_cond] = test_trajectory(start, stop, map, path, vis)
% TEST_TRAJECTORY simulates the robot from START to STOP following a PATH
% that's been planned for MAP.
% start - a 3d vector or a cell contains multiple 3d vectors
% stop  - a 3d vector or a cell contains multiple 3d vectors
% map   - map generated by your load_map
% path  - n x 3 matrix path planned by your dijkstra algorithm
% vis   - true for displaying visualization

%Controller and trajectory generator handles
controlhandle = @controller;
trajhandle    = @trajectory_generator;

% Make cell
if ~iscell(start), start = {start}; end
if ~iscell(stop),  stop  = {stop}; end
if ~iscell(path),  path  = {path} ;end

% Get nquad
nquad = length(start);

% Make column vector
for qn = 1:nquad
    start{qn} = start{qn}(:);
    stop{qn} = stop{qn}(:);
end

% Quadrotor model
params = crazyflie();

%% **************************** FIGURES *****************************
% Environment figure
if nargin < 5
    vis = true;
end

fprintf('Initializing figures...\n')
if vis
    h_fig = figure('Name', 'Environment');
else
    h_fig = figure('Name', 'Environment', 'Visible', 'Off');
end
if nquad == 1
    plot_path(map, path{1});
else
    % you could modify your plot_path to handle cell input for multiple robots
end
h_3d = gca;
drawnow;
xlabel('x [m]'); ylabel('y [m]'); zlabel('z [m]')
quadcolors = lines(nquad);
set(gcf,'Renderer','OpenGL')

%% *********************** INITIAL CONDITIONS ***********************
fprintf('Setting initial conditions...\n')
% Maximum time that the quadrotor is allowed to fly
time_tol = 30;          % maximum simulation time
starttime = 0;          % start of simulation in seconds
tstep     = 0.01;       % this determines the time step at which the solution is given
cstep     = 0.05;       % image capture time interval
nstep     = cstep/tstep;
time      = starttime;  % current time
max_iter  = time_tol / cstep;      % max iteration
for qn = 1:nquad
    % Get start and stop position
    x0{qn}    = init_state(start{qn}, 0);
    xtraj{qn} = zeros(max_iter*nstep, length(x0{qn}));
    ttraj{qn} = zeros(max_iter*nstep, 1);
end

% Maximum position error of the quadrotor at goal
pos_tol  = 0.05; % m
% Maximum speed of the quadrotor at goal
vel_tol  = 0.05; % m/s

x = x0;        % state

%% ************************* RUN SIMULATION *************************
OUTPUT_TO_VIDEO = 1;
if OUTPUT_TO_VIDEO == 1
    v = VideoWriter('map1.avi');
    open(v)
end

fprintf('Simulation Running....\n')
for iter = 1:max_iter
    timeint = time:tstep:time+cstep;
    tic;
    % Iterate over each quad
    for qn = 1:nquad
        % Initialize quad plot
        if iter == 1
            QP{qn} = QuadPlot(qn, x0{qn}, 0.1, 0.04, quadcolors(qn,:), max_iter, h_3d);
            desired_state = trajhandle(time, qn);
            QP{qn}.UpdateQuadPlot(x{qn}, [desired_state.pos; desired_state.vel], time);
            h_title = title(sprintf('iteration: %d, time: %4.2f', iter, time));
        end

        % Run simulation
        [tsave, xsave] = ode45(@(t,s) quadEOM(t, s, qn, controlhandle, trajhandle, params), timeint, x{qn});
        x{qn} = xsave(end, :)';
        % Save to traj
        xtraj{qn}((iter-1)*nstep+1:iter*nstep,:) = xsave(1:end-1,:);
        ttraj{qn}((iter-1)*nstep+1:iter*nstep)   = tsave(1:end-1);

        % Update quad plot
        desired_state = trajhandle(time + cstep, qn);
        QP{qn}.UpdateQuadPlot(x{qn}, [desired_state.pos; desired_state.vel], time + cstep);
        if OUTPUT_TO_VIDEO == 1
            im = frame2im(getframe(gcf));
            writeVideo(v,im);
        end
    end

    set(h_title, 'String', sprintf('iteration: %d, time: %4.2f', iter, time + cstep))
    time = time + cstep; % Update simulation time
    t = toc;

    % Pause to make real-time
    if (t < cstep)
        pause(cstep - t);
    end

    % Check termination criteria
    terminate_cond = terminate_check(x, time, stop, pos_tol, vel_tol, time_tol);
    if terminate_cond
        break
    end

end

fprintf('Simulation Finished....\n')

if OUTPUT_TO_VIDEO == 1
    close(v);
end

%% ************************* POST PROCESSING *************************
% Truncate xtraj and ttraj
for qn = 1:nquad
    xtraj{qn} = xtraj{qn}(1:iter*nstep,:);
    ttraj{qn} = ttraj{qn}(1:iter*nstep);
end

% Plot the saved position and velocity of each robot
if vis
    for qn = 1:nquad
        % Truncate saved variables
        QP{qn}.TruncateHist();
        % Plot position for each quad
        h_pos{qn} = figure('Name', ['Quad ' num2str(qn) ' : position']);
        plot_state(h_pos{qn}, QP{qn}.state_hist(1:3,:), QP{qn}.time_hist, 'pos', 'vic');
        plot_state(h_pos{qn}, QP{qn}.state_des_hist(1:3,:), QP{qn}.time_hist, 'pos', 'des');
        % Plot velocity for each quad
        h_vel{qn} = figure('Name', ['Quad ' num2str(qn) ' : velocity']);
        plot_state(h_vel{qn}, QP{qn}.state_hist(4:6,:), QP{qn}.time_hist, 'vel', 'vic');
        plot_state(h_vel{qn}, QP{qn}.state_des_hist(4:6,:), QP{qn}.time_hist, 'vel', 'des');
    end
end

end

author: yrlu

repo: quadrotor

save_path: github-repos/MATLAB/yrlu-quadrotor


text:
function exact = p28_exact ( )

%*****************************************************************************80
%
%% P28_EXACT returns the exact integral for problem 28.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  exact = ( 50.0 / 2501.0 ) * ( 1.0 - exp ( - 2.0 * pi ) );

  return
end


text:
%% ALL NONLINEAR VERSION
clc
clear
% Objective & Gradient
 fun = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
 grad = @(x)[2*x(1)-400*x(1)*(x(2)-x(1)^2)-2,200*x(2)-200*x(1)^2];

 % Constraints
A = [-1 1; 1 1];
rl = [-Inf;5];
ru = [-1;5];
lb = [0;0]; ub = [4;4];

% Nonlinear Constraint, Jacobian & Structure
 nlcon = @(x) A*x;
 nljac = @(x) sparse(A);
 jacstr = @() sparse(double(A~=0));

% Starting Guess
 x0 = [2;2];

% Build Function Structure
 funcs.objective = fun;
 funcs.gradient = grad;
 funcs.constraints = nlcon;
 funcs.jacobian = nljac;
 funcs.jacobianstructure = jacstr;
 funcs.hessian = @(x,sigma,lambda)sparse([[sigma*(1200*x(1)^2-400*x(2)+2),0];[-400*sigma*x(1),200*sigma]]);
 funcs.hessianstructure = @() sparse(tril(ones(2)));

% Build Options Structure
 opts.lb = lb;
 opts.ub = ub;
 opts.cl = rl;
 opts.cu = ru;
 opts.ipopt.linear_solver = 'pardiso';
%  opts.ipopt.hessian_approximation = 'limited-memory';

% Call IPOPT
[x,output] = ipopt(x0,funcs,opts)

%% ALL LINEAR VERSION
clc
clear all
% Objective & Gradient
 fun = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
 grad = @(x)[2*x(1)-400*x(1)*(x(2)-x(1)^2)-2,200*x(2)-200*x(1)^2];

 % Constraints
A = [-1 1; 1 1];
rl = [-Inf;5];
ru = [-1;5];
lb = [0;0]; ub = [4;4];

% Starting Guess
 x0 = [2;2];

% Build Function Structure
 funcs.objective = fun;
 funcs.gradient = grad;

% Build Options Structure
 opts.lb = lb;
 opts.ub = ub;
 opts.A = sparse(A);
 opts.rl = rl;
 opts.ru = ru;
 opts.ipopt.hessian_approximation = 'limited-memory';

% Call IPOPT
[x,output] = ipopt(x0,funcs,opts)

%% NL INEQ, LIN EQ
clc
clear all
% Objective & Gradient
 fun = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
 grad = @(x)[2*x(1)-400*x(1)*(x(2)-x(1)^2)-2,200*x(2)-200*x(1)^2];

 % Constraints
Aeq = [1 1];
rl = 5;
ru = 5;
lb = [0;0]; ub = [4;4];

% Nonlinear Constraint, Jacobian & Structure
 A = [-1 1];
 nlcon = @(x) [A]*x;
 nljac = @(x) sparse([A]);
 jacstr = @() sparse(double([A]~=0));
 cl = -Inf;
 cu = -1;

% Starting Guess
 x0 = [2;2];

% Build Function Structure
 funcs.objective = fun;
 funcs.gradient = grad;
 funcs.constraints = nlcon;
 funcs.jacobian = nljac;
 funcs.jacobianstructure = jacstr;

% Build Options Structure
 opts.lb = lb;
 opts.ub = ub;
 opts.A = sparse(Aeq);
 opts.rl = rl;
 opts.ru = ru;
 opts.cl = cl;
 opts.cu = cu;
 opts.ipopt.hessian_approximation = 'limited-memory';

% Call IPOPT
% for i = 1:1000
    [x,output] = ipopt(x0,funcs,opts)
% end

%% NL EQ, LIN INEQ
clc
clear all
% Objective & Gradient
 fun = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
 grad = @(x)[2*x(1)-400*x(1)*(x(2)-x(1)^2)-2,200*x(2)-200*x(1)^2];

 % Constraints
A = [-1 1];
rl = -Inf;
ru = -1;
lb = [0;0]; ub = [4;4];

% Nonlinear Constraint, Jacobian & Structure
 Aeq = [1 1];
 nlcon = @(x) [Aeq]*x;
 nljac = @(x) sparse([Aeq]);
 jacstr = @() double(sparse([Aeq]~=0));
 cl = 5;
 cu = 5;

% Starting Guess
 x0 = [2;2];

% Build Function Structure
 funcs.objective = fun;
 funcs.gradient = grad;
 funcs.constraints = nlcon;
 funcs.jacobian = nljac;
 funcs.jacobianstructure = jacstr;

% Build Options Structure
 opts.lb = lb;
 opts.ub = ub;
 opts.A = sparse(A);
 opts.rl = rl;
 opts.ru = ru;
 opts.cl = cl;
 opts.cu = cu;
 opts.ipopt.hessian_approximation = 'limited-memory';

% Call IPOPT
% for i = 1:1000
    [x,output] = ipopt(x0,funcs,opts)
% end

%% HESSIAN MODIFICATION FOR LAMBDA
clc
clear
f = @(x) x(1)*x(4)*sum(x(1:3)) + x(3);
g = @(x) [ x(1)*x(4) + x(4)*sum(x(1:3))
                x(1)*x(4)
                x(1)*x(4) + 1
                x(1)*sum(x(1:3)) ];
c = @(x) [ prod(x); sum(x.^2) ];
j = @(x) [ prod(x)./x'; 2*x' ];

cl = [25;40];
cu = [Inf;40];

H = @(x,sigma,lambda) sigma*[ 2*x(4)             0      0   0;
                              x(4)               0      0   0;
                              x(4)               0      0   0;
                              2*x(1)+x(2)+x(3)  x(1)  x(1)  0 ] + ...
                   lambda(1)*[   0          0         0         0;
                              x(3)*x(4)     0         0         0;
                              x(2)*x(4) x(1)*x(4)     0         0;
                              x(2)*x(3) x(1)*x(3) x(1)*x(2)     0  ] + ...
         		   lambda(2)*diag([2 2 2 2]);

lb = ones(4,1);
ub = 5*ones(4,1);

A = [1 0 0 0];
rl = 1.5;
ru = 1.5;

x0 = [1 5 5 1]';
opts = optiset('display','iter');

Opt = opti('fun',f,'grad',g,'hess',H,'nl',c,cl,cu,'nljac',j,'lin',A,rl,ru,'bounds',lb,ub,'options',opts)

[x,fval,exitflag,info]= solve(Opt,x0)

%%
clc
%Objective
obj = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
% Constraints
A = [-1 1];
b = -1;
Aeq = [1 1];
beq = 5;
lb = [0;0]; ub = [4;4];
% Solve
x0 = [2;2];
opts = optiset('solver','ipopt','display','iter');
Opt = opti('obj',obj,'ndec',2,'bounds',lb,ub,'ineq',A,b,'eq',Aeq,beq,'options',opts)
[x,fval,ef,info] = solve(Opt,x0)
%Plot
% plot(Opt,[],1)

%%
n = 2000;

% Objective & Gradient
fun = @(x) -sum(x);
grad = @(x) -ones(n,1);

% Linear constraints
A = randn(n,n);
rl = zeros(n,1);
ru = zeros(n,1);
%%
% Nonlinear Constraint, Jacobian & Structure
nlcon = @(x) x(:).^4;
nljac = @(x) sparse(diag(4*x(:).^3));
jacstr = @() speye(n);
hessian = @(x,sigma,lambda) sparse(diag(12*x(:).^2));
hessianstructure = @() speye(n);
cl = -Inf(n,1);
cu = ones(n,1);

% Starting Guess
x0 = 1*ones(n,1);
% Build Function Structure
funcs.objective = fun;
funcs.gradient = grad;
funcs.constraints = nlcon;
funcs.jacobian = nljac;
funcs.jacobianstructure = jacstr;
% funcs.hessian = hessian;
% funcs.hessianstructure = hessianstructure;

% Build Options Structure
lb = -10*ones(n,1);
ub = 40*ones(n,1);
opts = [];
opts.lb = lb;
opts.ub = ub;
opts.rl = rl;
opts.ru = ru;
opts.cl = cl;
opts.cu = cu;
opts.A = sparse(A);
opts.ipopt.hessian_approximation = 'limited-memory';
opts.ipopt.linear_solver = 'pardiso';
% Call IPOPT
[~,output] = ipopt(x0,funcs,opts);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function cellnorm = normCellVector(cellvec)
% this function calculates the second norm of a cell vector field
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Written by Ali A. Eftekhari
% See the license file

% check the size of the variable and the mesh dimension
d = cellvec.domain.dimension;

if (d ==1) || (d==1.5) || (d==1.8)
	cellnormval = abs(cellvec.xvalue);
elseif (d == 2) || (d == 2.5) || (d==2.8)
	cellnormval = realsqrt(cellvec.xvalue.*cellvec.xvalue+ ...
        cellvec.yvalue.*cellvec.yvalue);
elseif (d == 3) || (d==3.2)
    cellnormval = realsqrt(cellvec.xvalue.*cellvec.xvalue+ ...
        cellvec.yvalue.*cellvec.yvalue+ ...
        cellvec.zvalue.*cellvec.zvalue);
end
BC = createBC(cellvec.domain);
c=cellBoundary(cellnormval, BC);
cellnorm=CellVariable(cellvec.domain, c);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
% Juan Manuel Perez Rua

% I =input image
% algorithm =name of the algorithm
% varargin, threshold list, if needed. AN error message error appears
% if incorrect inputs are used.

%Example of usage: J = colorConstancy(I, 'modified white patch', 200);
function [OUT] = colorConstancy(I, algorithm, varargin)

    [m,n,~]=size(I);
    switch (algorithm)
        case 'grey world'
            Rmean      = sum(sum(I(:,:,1)))/(m*n);
            Gmean      = sum(sum(I(:,:,2)))/(m*n);
            Bmean      = sum(sum(I(:,:,3)))/(m*n);
            Avg        = mean([Rmean Gmean Bmean]);
            Kr         = Avg/Rmean;
            Kg         = Avg/Gmean;
            Kb         = Avg/Bmean;
            OUT(:,:,1) = Kr*double(I(:,:,1));
            OUT(:,:,2) = Kg*double(I(:,:,2));
            OUT(:,:,3) = Kb*double(I(:,:,3));
            OUT = uint8(OUT);
        case 'white patch'
            Kr = 255/max(max(double(I(:,:,1))));
            Kg = 255/max(max(double(I(:,:,2))));
            Kb = 255/max(max(double(I(:,:,3))));
            OUT(:,:,1) = Kr*double(I(:,:,1));
            OUT(:,:,2) = Kg*double(I(:,:,2));
            OUT(:,:,3) = Kb*double(I(:,:,3));
            OUT = uint8(OUT);
        case 'modified white patch'
            if (~isempty(varargin))
                th = varargin{1};
                R=I(:,:,1); Kr = 255/mean(R(R>th));
                G=I(:,:,2); Kg = 255/mean(G(G>th));
                B=I(:,:,3); Kb = 255/mean(B(B>th));
                OUT(:,:,1) = Kr*double(I(:,:,1));
                OUT(:,:,2) = Kg*double(I(:,:,2));
                OUT(:,:,3) = Kb*double(I(:,:,3));
                OUT = uint8(OUT);
            else
                OUT=I; disp('Modified white patch algorithm must have another parameter.');
            end
        case 'progressive'
            if (length(varargin)>1)
                h1 = varargin{1};
                h2 = varargin{2};
                imap = (double(I(:,:,1))+double(I(:,:,2))+double(I(:,:,3)))/3;
                R=I(:,:,1); G=I(:,:,2); B=I(:,:,3);

                Kr = zeros(size(imap));
                Kg = zeros(size(imap));
                Kb = zeros(size(imap));

                Kr( imap>=h1 )=255/mean(R(R>=h1));
                Kg( imap>=h1 )=255/mean(G(G>=h1));
                Kb( imap>=h1 )=255/mean(B(B>=h1));

                [m,n,~]=size(I);
                Rmean = sum(sum(R))/(m*n);
                Gmean = sum(sum(G))/(m*n);
                Bmean = sum(sum(B))/(m*n);
                Avg = mean([Rmean Gmean Bmean]);

                Kr( imap<=h2 ) = Avg/Rmean;
                Kg( imap<=h2 ) = Avg/Gmean;
                Kb( imap<=h2 ) = Avg/Bmean;

                deltha = imap/(h1-h2)-h2/(h1-h2);
                Kr( imap>=h2 & imap<=h1 ) = (1-deltha(imap>=h2 & imap<=h1))*...
                    255/mean(R(R>=h1)) + deltha(imap>=h2 & imap<=h1)*255/mean(R(R>=h1));
                Kg( imap>=h2 & imap<=h1 ) = (1-deltha(imap>=h2 & imap<=h1))*...
                    255/mean(G(G>=h1)) + deltha(imap>=h2 & imap<=h1)*255/mean(G(G>=h1));
                Kb( imap>=h2 & imap<=h1 ) = (1-deltha(imap>=h2 & imap<=h1))*...
                    255/mean(B(B>=h1)) + deltha(imap>=h2 & imap<=h1)*255/mean(B(B>=h1));
                OUT(:,:,1) = Kr.*double(I(:,:,1));
                OUT(:,:,2) = Kg.*double(I(:,:,2));
                OUT(:,:,3) = Kb.*double(I(:,:,3));
                OUT = uint8(OUT);
            else
                OUT=I; disp('Modified white patch algorithm must have another parameter.');
            end
        case 'single scale retinex'
            if (~isempty(varargin))
                c=varargin{1};
                [Y,X]=meshgrid(1:n,1:m);

                Fnok = exp(-((X.^2)+(Y.^2))./(c.^2));
                K = 1/(sum(sum(Fnok)));
                F = K.*Fnok;

                IR = double(I(:,:,1));
                IG = double(I(:,:,2));
                IB = double(I(:,:,3));
                FF  = fftshift(fft2(F));
                IFR = fftshift(fft2(IR)); IFR=FF.*IFR; IFR=real(ifft2(ifftshift(IFR)));
                IFG = fftshift(fft2(IG)); IFG=FF.*IFG; IFG=real(ifft2(ifftshift(IFG)));
                IFB = fftshift(fft2(IB)); IFB=FF.*IFB; IFB=real(ifft2(ifftshift(IFB)));

                RR = log10(double(IR))-log10(IFR);
                RG = log10(double(IG))-log10(IFG);
                RB = log10(double(IB))-log10(IFB);

                OUT(:,:,1)=uint8(255*RR/(max(max([RR RG RB]))));
                OUT(:,:,2)=uint8(255*RG/(max(max([RR RG RB]))));
                OUT(:,:,3)=uint8(255*RB/(max(max([RR RG RB]))));
            else
                OUT=I; disp('Single scale algorithm must have another parameter (c).');
            end
        case 'multi scale retinex'
            R1=colorConstancy(I,'single scale retinex',25);
            R2=colorConstancy(I,'single scale retinex',100);
            R3=colorConstancy(I,'single scale retinex',240);
            OUT=(1/3)*R1+(1/3)*R2+(1/3)*R3;
        case 'MSRCR'
            if (length(varargin)>1)
              alpha=varargin{1};
              betha=varargin{2};
              gain=varargin{3};
              Rmsr=colorConstancy(I, 'multi scale retinex');
              T(:,:,1)=I(:,:,1)+I(:,:,2)+I(:,:,3);
              T(:,:,2)=I(:,:,1)+I(:,:,2)+I(:,:,3);
              T(:,:,3)=I(:,:,1)+I(:,:,2)+I(:,:,3);
              C = log10(alpha*double(I))-log10(double(T));
              OUT = uint8(gain*(C.*double(Rmsr)+betha));
            else
                OUT=I; disp('MSRCR must have another parameter alpha, b and G.');
            end
        case 'ace'
            MR = zeros(m,n);
            MG = zeros(m,n);
            MB = zeros(m,n);
            RC = zeros(m,n,3);
            [Y,X]=meshgrid(1:n,1:m);
            cont=0;
            W=20;
            %STEP 1
            for i=1:m
               for j=1:n

                   cont=cont+1;
                   if (mod(cont,1000)==0) disp(['Ciclos: ' num2str(cont)]); end

                   MD = sqrt(((X-X(i,j)).^2)+((Y-Y(i,j)).^2));

                   MR(MD<=W & MD~=0)=sign(double(I(i,j,1))-double(I(MD<=W & MD~=0)));
                   MG(MD<=W & MD~=0)=sign(double(I(i,j,2))-double(I(MD<=W & MD~=0)));
                   MB(MD<=W & MD~=0)=sign(double(I(i,j,3))-double(I(MD<=W & MD~=0)));

                   T = zeros(size(m,n));
                   T(MD<=W & MD~=0)= MR(MD<=W & MD~=0)./MD(MD<=W & MD~=0);
                   RC(i,j,1)=sum(sum(T));

                   T = zeros(size(m,n));
                   T(MD<=W & MD~=0)= MG(MD<=W & MD~=0)./MD(MD<=W & MD~=0);
                   RC(i,j,2)=sum(sum(T));

                   T = zeros(size(m,n));
                   T(MD<=W & MD~=0)= MB(MD<=W & MD~=0)./MD(MD<=W & MD~=0);
                   RC(i,j,3)=sum(sum(T));
               end
            end
            %STEP 2
            mc = min(min(min(RC)));
            Mc = max(max(max(RC)));
            s = (255-Mc)/(-mc);
            OUT = uint8(round(127.5+s*RC));
        otherwise
            OUT=I; disp('Unknown alogrithm, please check name.');
    end
end

author: Sable

repo: mcbench-benchmarks


function [  ] = getColsCuts( visMatrixDir, rowsCutDir, saveFileDir )
%GETCOLSCUTS Summary of this function goes here
%   Detailed explanation goes here

    fileExt = '*.txt';
    files_pose = dir(fullfile(rowsCutDir, fileExt));

    rosCutLength= length(files_pose);

    for i=1:rosCutLength
        i
        fileName_pose = [rowsCutDir, num2str(i), '.txt'];
        file_pose = fopen(fileName_pose);
        pose_ID = fscanf(file_pose, '%d');

        points_for_section = [];  % all the ovserved Points
        for j = 1:length(pose_ID)
            fileName_point = [visMatrixDir, num2str(pose_ID(j)), '.txt'];
            file_point = fopen(fileName_point);
            point_ID = fscanf(file_pose, '%d');
            fclose(file_pose);

            %
            points_for_section = union(points_for_section, point_ID);
        end

        fileName_savePoint = [saveFileDir, num2str(i), '.txt'];
        dlmwrite(fileName_savePoint, points_for_section, 'precision', 10);

    end

end


text:
function [ phi, dphidx, dphidy ] = basis_11_t4 ( t, i, p )

%*****************************************************************************80
%
%% BASIS_11_T4: one basis at one point for a T4 element.
%
%  Discussion:
%
%    The T4 element is the cubic bubble triangle.
%
%    The routine is given the coordinates of the vertices of a triangle.
%    It works directly with these coordinates, and does not refer to a
%    reference element.
%
%    The sides of the triangle DO NOT have to lie along a coordinate
%    axis.
%
%    The routine evaluates the basis functions associated with each vertex,
%    and their derivatives with respect to X and Y.
%
%  Physical Element T4:
%
%            3
%           / \
%          /   \
%         /  4  \
%        /       \
%       1---------2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 March 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,4), the coordinates of the vertices
%    of the triangle, and the coordinates of the centroid.
%    It is common to list the first three points in counter clockwise
%    order.
%
%    Input, integer I, the index of the basis function.
%
%    Input, real P(2), the points where the basis function
%    is to be evaluated.
%
%    Output, real PHI, the value of the basis function
%    at the evaluation point.
%
%    Output, real DPHIDX, DPHIDY, the value of the
%    derivatives at the evaluation point.
%
%  Local parameters:
%
%    Local, real AREA, is (twice) the area of the triangle.
%
  area = t(1,1) * ( t(2,2) - t(2,3) ) ...
       + t(1,2) * ( t(2,3) - t(2,1) ) ...
       + t(1,3) * ( t(2,1) - t(2,2) );

  psi(1) =     (   ( t(1,3) - t(1,2) ) * ( p(2) - t(2,2) )     ...
                 - ( t(2,3) - t(2,2) ) * ( p(1) - t(1,2) ) );
  dpsidx(1) =    - ( t(2,3) - t(2,2) );
  dpsidy(1) =      ( t(1,3) - t(1,2) );

  psi(2) =     (   ( t(1,1) - t(1,3) ) * ( p(2) - t(2,3) )     ...
                 - ( t(2,1) - t(2,3) ) * ( p(1) - t(1,3) ) );
  dpsidx(2) =    - ( t(2,1) - t(2,3) );
  dpsidy(2) =      ( t(1,1) - t(1,3) );

  psi(3) =     (   ( t(1,2) - t(1,1) ) * ( p(2) - t(2,1) )     ...
                 - ( t(2,2) - t(2,1) ) * ( p(1) - t(1,1) ) );
  dpsidx(3) =    - ( t(2,2) - t(2,1) );
  dpsidy(3) =      ( t(1,2) - t(1,1) );
%
%  Normalize the first three functions.
%
  psi(1:3)    =    psi(1:3) / area;
  dpsidx(1:3) = dpsidx(1:3) / area;
  dpsidy(1:3) = dpsidy(1:3) / area;
%
%  Compute the cubic bubble function.
%
  psi(4) = 27.0D+00 * psi(1) * psi(2) * psi(3);

  dpsidx(4) = 27.0D+00 * ( ...
                dpsidx(1) *    psi(2) *    psi(3) ...
                +  psi(1) * dpsidx(2) *    psi(3) ...
                +  psi(1) *    psi(2) * dpsidx(3) );

  dpsidy(4) = 27.0D+00 * ( ...
                dpsidy(1) *    psi(2) *    psi(3) ...
                +  psi(1) * dpsidy(2) *    psi(3) ...
                +  psi(1) *    psi(2) * dpsidy(3) );
%
%  Subtract 1/3 of the cubic bubble function from each of the three linears.
%
  for j = 1 : 3
    psi(j)    =    psi(j) -    psi(4) / 3.0;
    dpsidx(j) = dpsidx(j) - dpsidx(4) / 3.0;
    dpsidy(j) = dpsidy(j) - dpsidy(4) / 3.0;
  end

  phi    = psi(i);
  dphidx = dpsidx(i);
  dphidy = dpsidy(i);

  return
end


text:
function runner_info=prepare_runner_test_simple(run_config)

class_info=run_config.class_info;

ds_config=[];

% set the input image scales, useful for multi-scale evaluation
run_config.input_img_scale=0.8;
run_config.gen_net_opts_fn=@gen_net_opts_model_type1;

run_config.run_evaonly=true;
ds_config.use_custom_data=true;
ds_config.use_dummy_gt=true;
run_config.use_dummy_gt=ds_config.use_dummy_gt;

ds_config.ds_name='tmp_data';
run_config.root_cache_dir=fullfile('../cache_data', 'tmp_runner_dir');
% mkdir_notexist(run_config.root_cache_dir);

run_config.model_name='tmp_model_name';

run_dir_name=fileparts(mfilename('fullpath'));
[~, run_dir_name]=fileparts(run_dir_name);
run_config.run_dir_name=run_dir_name;
run_config.run_file_name=mfilename();

train_opts=run_config.gen_net_opts_fn(run_config, class_info);

net_run_opts=train_opts.net_run_opts;
net_run_opts.batch_info_disp_fn=[];
net_run_opts.epoch_evaluate_fn=@(work_info, opts, imdb, net_config, work_info_epoch){};
train_opts.net_run_opts=net_run_opts;

ds_info=[];
ds_info.train_idxes=1;
ds_info.test_idxes=1;
imdb=my_gen_imdb(train_opts, ds_info);

data_norm_info=[];
data_norm_info.image_mean=128;

imdb.ref.data_norm_info=data_norm_info;

if run_config.use_gpu
	gpu_num=gpuDeviceCount;
	if gpu_num>=1
		gpuDevice(run_config.gpu_idx);
    else
        error('no gpu found!');
	end
end

[net_config, net_exp_info]=prepare_running_model(train_opts);

loss_group_info=net_config.ref.group_infos{end};
loss_group_info.forward_evaluate_fn=@tmp_loss_group_evaluate_runner;
net_config.ref.group_infos{end}=loss_group_info;

runner_info=[];
runner_info.train_opts=train_opts;
runner_info.imdb=imdb;
runner_info.net_config=net_config;
runner_info.net_exp_info=net_exp_info;
runner_info.class_info=class_info;
runner_info.ds_config=ds_config;

runner_info.run_task_fn=@run_task;

end

function ds_info=gen_runner_ds_info(ds_config, class_info, task_info)

    ds_info=[];

    img_dir=task_info.img_dir;
    img_files={task_info.img_filename};
    img_num=1;
    ds_info.img_idxes=uint32(1:img_num);
    ds_info.img_num=img_num;
    ds_info.train_idxes=uint32([]);
    ds_info.test_idxes=ds_info.img_idxes;

    img_names=cell(img_num, 1);
    for img_idx=1:img_num
        one_name=img_files{img_idx};
        [~, one_name]=fileparts(one_name);
        img_names{img_idx}=one_name;
    end

    ds_info.img_files=img_files;
    ds_info.img_names=img_names;
    ds_info.mask_files=[];
    ds_info.train_idxes=[];
    ds_info.test_idxes=uint32(1:img_num)';
    ds_info.name=ds_config.ds_name;

    ds_info.data_dir_idxes_img=repmat(uint8(1), [img_num, 1]);
    ds_info.data_dir_idxes_mask=[];
    ds_info.data_dirs={img_dir};

    ds_info.class_info=class_info;

end

function task_result=run_task(runner_info, task_info)

runner_info.imdb.ref.ds_info=gen_runner_ds_info(runner_info.ds_config, runner_info.class_info, task_info);
exp_info=my_net_tool(runner_info.train_opts, runner_info.imdb, runner_info.net_config, runner_info.net_exp_info);
work_info=exp_info.eva_val;

task_result=work_info.ref.task_result_imgs{1};

end

function tmp_loss_group_evaluate_runner(work_info_batch, group_idx)

    cnn_forward_evaluate_obj(work_info_batch, group_idx);
    cnn_softmax_evaluate(work_info_batch, group_idx);

    if ~work_info_batch.ref.run_eva
        return;
    end

    predict_info=gen_predict_info(work_info_batch, group_idx);
    do_full_eva(work_info_batch, predict_info);

end

function predict_info=gen_predict_info(work_info_batch, group_idx)

    prediction_info=work_info_batch.ref.prediction_info_groups{group_idx};
    one_output_info=prediction_info.output_infos{end};
    assert(check_valid_net_output(one_output_info));
    predict_info=one_output_info.mc_predict_info;

end

function do_full_eva(work_info_batch, predict_info)

    work_info=work_info_batch.ref.work_info;
    opts=work_info_batch.ref.train_opts;

    %assume no crop or data_aug
    batch_ds_info=work_info_batch.ref.ds_info;
    assert(isempty(batch_ds_info.aug_info));
    assert(isempty(batch_ds_info.crop_info));

    task_subidx=work_info_batch.ref.task_subidxes;
    img_idx=work_info_batch.ref.task_idxes;
    img_num=length(task_subidx);
    assert(img_num==1);

    seg_param=[];
    seg_param.eva_param=opts.eva_param;
    seg_param.img_idx=img_idx;

    batch_data=batch_ds_info.batch_data;
    seg_param.img_data_input=batch_data.img_data;
    seg_param.gt_mask_input=batch_data.label_data;
    seg_param.img_size_input=batch_data.img_size_origin;

    ds_info=work_info_batch.ref.imdb.ref.ds_info;
    result_info=do_gen_seg_result(seg_param, ds_info, predict_info);

    work_info.ref.task_result_imgs{task_subidx}=result_info;

end

function result_info=do_gen_seg_result(seg_param, ds_info, predict_info)

    img_size=seg_param.img_size_input;

    assert(~isempty(predict_info));

    score_map=single(predict_info.score_map);
    score_map_size=size(score_map);
    score_map_size=score_map_size(1:2);

    score_map_org=score_map;

    if any(img_size~=score_map_size)
        score_map=log(score_map);
        score_map=max(score_map, -20);
        score_map=my_resize(score_map, img_size);
        score_map=exp(score_map);
    end

    [~, predict_mask]=max(score_map,[],3);
    predict_mask=uint8(gather(predict_mask));

    result_info=gen_task_result(ds_info, seg_param, predict_mask, score_map_org);

end

function tmp_result_info=gen_task_result(ds_info, seg_param, predict_mask_net, score_map_org)

    img_idx=seg_param.img_idx;
    eva_param=seg_param.eva_param;
    class_info=eva_param.class_info;

    predict_mask_data=class_info.class_label_values(predict_mask_net);
    assert(isa(class_info.class_label_values, 'uint8'));
    assert(isa(predict_mask_data, 'uint8'));

    img_name=ds_info.img_names{img_idx};

    tmp_result_info=[];
    tmp_result_info.mask_data=predict_mask_data;
    tmp_result_info.img_size=size(predict_mask_data);
    tmp_result_info.class_info=class_info;
    tmp_result_info.img_name=img_name;
    tmp_result_info.score_map=im2uint8(score_map_org);

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function dy = tcell_rhs ( t, y )

%*****************************************************************************80
%
%% TCELL_RHS evaluates the right hand side for the ODE solver TESTGM.
%
%  Discussion:
%
%    This function evaluates the right hand side of an ODE for the
%    time evolution of the solution coefficients for a reduced order
%    model of the Navier Stokes equations.
%
%    The file MATRIX.OUT contains information about various coefficient
%    matrices needed to evaluate the right hand side, and is loaded by
%    this routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2004
%
%  Author:
%
%    Hyung-Chun Lee
%    Department of Mathematics
%    Ajou University, Korea
%
%  Parameters:
%
%    Input, real T, the current time.
%
%    Input, real Y(8), the current solution.
%
%    Output, real DY(8), the corresponding value of dY/dT.
%

%
%  Load the coefficient data from the file.
%
  load matrix.out
%
%  Distribute the data from the file.
%
  M = zeros(8,8);
  M = matrix(1:8,1:8);

  H = zeros(8,8);
  H = matrix(9:16,1:8);

  Q = zeros(8,8,8);
  for i = 1 : 8
    ii = 16+8*(i-1)+1;
    ij = 16+8*(i-1)+8;
    Q(:,:,i) = matrix(ii:ij,1:8);
  end

  N = zeros(8,1);
  N = matrix(81,1:8);

  R = zeros(8,8);
  R = matrix(82:89,1:8);

  C = zeros(8,1);
  C = matrix(90,1:8);

  MI = inv ( M );

  mu = 1;

  if ( t < 0.05 )
    alpha = 80.0 * t + 1.0;
  else
    alpha = -80.0 * ( t - 0.1 ) + 1.0;
  end

  beta = alpha / 3.0;

  if ( t < 0.05 )
    dbeta = 80.0 / 3.0;
  else
    dbeta = -80.0 / 3.0;
  end

  dy = zeros(8,1);
  for i = 1 : 8
    dy(i) = -mu * H(i,:) * y - y' * Q(:,:,i) * y - beta * R(i,:) * y ...
      - beta * ( beta - 1.0 ) * C(i) - dbeta * N(i);
  end

  dy = MI * dy;

  return
end