text:
function [hits,ninvalid] = hits(bmus, mmax, values)

%HITS Calculate number of occurances of each value.
%
% hits = hits(bmus,[mmax],[values])
%
%   h = hits(bmus);
%   h = hits(bmus,length(sM.codebook));
%
%  Input and output arguments ([]'s are optional):
%   bmus     (vector) BMU indeces (or other similar)
%   [mmax]   (scalar) maximum index, default value max(bmus)
%            (struct) map or topology struct from where the maximum
%                     index is acquired
%   [values] (vector) values associated with the data (default = 1)
%
%   hits     (vector) the number of occurances of each index
%                     (or if values are given, their sum for each index)
%   ninvalid (scalar) number of invalid indeces (NaN, Inf or
%                     <=0 or > mmax)
%
% See also SOM_HITS, SOM_BMUS.

% Copyright (c) 2002 by the SOM toolbox programming team.
% Contributed to SOM Toolbox by Juha Vesanto, April 24th, 2002
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta juuso 240402

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<2 || isempty(mmax),
  mmax = max(bmus);
elseif isstruct(mmax),
  switch mmax.type,
   case 'som_map',   mmax = prod(mmax.topol.msize);
   case 'som_topol', mmax = prod(mmax.msize);
   otherwise,
    error('Illegal struct for 2nd argument.')
  end
end

if nargin<3, values = 1; end

valid_bmus = find(isfinite(bmus) & bmus>0 & bmus<=mmax);
ninvalid = length(bmus)-length(valid_bmus);

bmus = bmus(valid_bmus);
if length(values)>length(bmus), values = values(valid_bmus); end
hits = full(sum(sparse(bmus,1:length(bmus),values,mmax,length(bmus)),2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: ilarinieminen

repo: SOM-Toolbox


text:
function prepareDummyTracksKatie(moviename)

expdir = fileparts(moviename);
[~,nframes,fid,~] = get_readframe_fcn(moviename);
if fid>0,
  fclose(fid);
end
T = struct;
T.timestamps = ones(1,nframes)*now;
T.trx.x = ones(1,nframes)*5;
T.trx.y = ones(1,nframes)*5;
T.trx.a = ones(1,nframes)*1;
T.trx.b = ones(1,nframes)*1;
T.trx.theta = zeros(1,nframes);
T.trx.id = 1;
T.trx.firstframe = 1;
T.trx.endframe = nframes;
T.trx.x_mm = T.trx.x;
T.trx.y_mm = T.trx.y;
T.trx.a_mm = T.trx.a;
T.trx.b_mm = T.trx.b;
T.trx.theta_mm = T.trx.theta;
T.trx.dt = 1/1300;
T.trx.fps = 1300;
T.trx.pxpermm = 1;
T.trx.sex = 'M';
T.trx.off = 0;
T.trx.nframes = nframes;

save(fullfile(expdir,'trx.mat'),'-struct','T');

author: kristinbranson

repo: JAABA


text:
%SHOW	show all parameters of LINK object
%
%	SHOW(link)

% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function show(l)

	llab = 6;
	for n =fieldnames(l)'
		v = getfield(l, char(n));
		name = char(n);
		spaces = char(' '*ones(1,llab-length(name)));
		val = num2str(v);
		label = [name spaces ' = '];
		if numrows(val) > 1,
			pad = {label; char(' '*ones(numrows(val)-1,1))};
		else
			pad = label;
		end
		disp([char(pad) val]);
	end

author: petercorke

repo: robotics-toolbox-matlab


text:
function list = find(varargin)
% XMLTREE/FIND Find elements in a tree with specified characteristics
% FORMAT list = find(varargin)
%
% tree  - XMLTree object
% xpath - string path with specific grammar (XPath)
% uid   - lists of root uid's
% parameter/value - pair of pattern
% list  - list of uid's of matched elements
%
%      list = find(tree,xpath)
%      list = find(tree,parameter,value[,parameter,value])
%      list = find(tree,uid,parameter,value[,parameter,value])
%
% Grammar for addressing parts of an XML document:
%        XML Path Language XPath (http://www.w3.org/TR/xpath)
% Example: /element1//element2[1]/element3[5]/element4
%__________________________________________________________________________
%
% Find elements in an XML tree with specified characteristics or given
% a path (using a subset of XPath language).
%__________________________________________________________________________
% Copyright (C) 2002-2011  http://www.artefact.tk/

% Guillaume Flandin
% $Id: find.m 4460 2011-09-05 14:52:16Z guillaume $

% TODO:
%   - clean up subroutines
%   - find should only be documented using XPath (other use is internal)
%   - handle '*', 'last()' in [] and '@'
%   - if a key is invalid, should rather return [] than error ?

if nargin==0
    error('[XMLTree] A tree must be provided');
elseif nargin==1
    list = 1:length(tree.tree);
    return
elseif mod(nargin,2)
    list = sub_find_subtree1(varargin{1}.tree,root(varargin{1}),varargin{2:end});
elseif isa(varargin{2},'double') && ...
       ndims(varargin{2}) == 2 && ...
       min(size(varargin{2})) == 1
    list = unique(sub_find_subtree1(varargin{1}.tree,varargin{2:end}));
elseif nargin==2 && ischar(varargin{2})
    list = sub_pathfinder(varargin{:});
else
   error('[XMLTree] Arguments must be parameter/value pairs.');
end

%==========================================================================
function list = sub_find_subtree1(varargin)
    list = [];
    for i=1:length(varargin{2})
        res = sub_find_subtree2(varargin{1},...
                varargin{2}(i),varargin{3:end});
        list = [list res];
    end

%==========================================================================
function list = sub_find_subtree2(varargin)
    uid = varargin{2};
    list = [];
    if sub_comp_element(varargin{1}{uid},varargin{3:end})
        list = [list varargin{1}{uid}.uid];
    end
    if isfield(varargin{1}{uid},'contents')
        list = [list sub_find_subtree1(varargin{1},...
                varargin{1}{uid}.contents,varargin{3:end})];
    end

%==========================================================================
function match = sub_comp_element(varargin)
match = 0;
try
    % v = getfield(varargin{1}, varargin{2}); % slow...
    for i=1:floor(nargin/2)
        v = subsref(varargin{1}, struct('type','.','subs',varargin{i+1}));
        if strcmp(v,varargin{i+2})
            match = 1;
        end
    end
catch
end

% More powerful but much slower
%match = 0;
%for i=1:length(floor(nargin/2)) % bug: remove length !!!
%   if isfield(varargin{1},varargin{i+1})
%     if ischar(getfield(varargin{1},varargin{i+1})) & ischar(varargin{i+2})
%      if strcmp(getfield(varargin{1},varargin{i+1}),char(varargin{i+2}))
%           match = 1;
%        end
%     elseif isa(getfield(varargin{1},varargin{i+1}),'double') & ...
%           isa(varargin{i+2},'double')
%        if getfield(varargin{1},varargin{i+1}) == varargin{i+2}
%           match = 1;
%       end
%    else
%       warning('Cannot compare different objects');
%     end
%   end
%end

%==========================================================================
function list = sub_pathfinder(tree,pth)
    %- Search for the delimiter '/' in the path
    i = strfind(pth,'/');
    %- Begin search by root
    list = root(tree);
    %- Walk through the tree
    j = 1;
    while j <= length(i)
        %- Look for recursion '//'
        if j<length(i) && i(j+1)==i(j)+1
            recursive = 1;
            j = j + 1;
        else
            recursive = 0;
        end
        %- Catch the current tag 'element[x]'
        if j ~= length(i)
            element = pth(i(j)+1:i(j+1)-1);
        else
            element = pth(i(j)+1:end);
        end
        %- Search for [] brackets
        k = xml_findstr(element,'[',1,1);
        %- If brackets are present in current element
        if ~isempty(k)
            l   = xml_findstr(element,']',1,1);
            val = str2num(element(k+1:l-1));
            element = element(1:k-1);
        end
        %- Use recursivity
        if recursive
            list = find(tree,list,'name',element);
        %- Just look at children
        else
            if i(j)==1 % if '/root/...' (list = root(tree) in that case)
                if sub_comp_element(tree.tree{list},'name',element)
                    % list = 1; % list still contains root(tree)
                else
                    list = [];
                    return;
                end
            else
                list = sub_findchild(tree,list,element);
            end
        end
        % If an element is specified using a key
        if ~isempty(k)
            if val < 1 || val > length(list)+1
                error('[XMLTree] Bad key in the path.');
            elseif val == length(list)+1
                list = [];
                return;
            end
            list = list(val);
        end
        if isempty(list), return; end
        j = j + 1;
    end

%==========================================================================
function list = sub_findchild(tree,listt,elmt)
    list = [];
    for a=1:length(listt)
        for b=1:length(tree.tree{listt(a)}.contents)
            if sub_comp_element(tree.tree{tree.tree{listt(a)}.contents(b)},'name',elmt)
                list = [list tree.tree{tree.tree{listt(a)}.contents(b)}.uid];
            end
        end
    end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
I = vl_test_pattern(1) ;
ur = 1:size(I,2) ;
vr = 1:size(I,1) ;

[u,v] = meshgrid(ur(1:5:end),vr(1:5:end)) ;

f = [u(:)';v(:)'] ;
K = size(f,2) ;
f = [f ; 2 * ones(1,K) ; 0 * ones(1,K)] ;

f = vl_sift(single(I), 'frames', f, 'orientations') ;

%f = diag([1 1 6 1]) * f ;

figure(1) ; clf ;
imagesc(I) ; colormap gray ; hold on ;
vl_plotframe(f,'color','k','linewidth',3) ;
vl_plotframe(f,'color','y','linewidth',2) ; axis equal ; axis off ;

vl_demo_print('sift_or') ;

author: jianxiongxiao

repo: ProfXkit


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the following exercises:
%   A) Exercise A: Use the Runge-Kutta 4 algorithm to integrate the
%   differential equation:
%            dy/dt = 2*t
%      Integrate from t0=0, to tfinal = 10 s.
%      Compare the solution with the algebraic integration.
%   B) Exercise B: Use the Runge-Kutta 4 algorithm to integrate the
%   second order differential equation:
%            d2y/dt^2+5*dy/dt + 4*y(t) = 0
%
%   C) Exercise C: Use the Runge-Kutta 4 algorithm to simulate the movement
%   of a 1 dof robot arm with friction under the effect of gravity and with
%   zero torque applied.
%
%   D) Exercise D: Use the Runge-Kutta 4 algorithm to simulate the movement
%   of a 2 dof robot arm with friction under the effect of gravity and with
%   zero torques applied.
%
% Help: Function prototype
% [y, t] = runge_kutta(f, y0, [t0 tfinal], timestep)
% where f is the function being integrated as dy/dt = f(t, y).
% y0 are the initial conditions
% t0: initial time
% tfinal: final time.
% h: time step for the calculations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function runge_kutta_exercises()
close all;

%uncomment to execute each of the exercises
%exerciseA()
%exerciseB()
%exerciseC()
exerciseD()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Integrate a simple time function. dy/dt = 2*t
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function exerciseA()
t0 = 0;
tfinal = 10;
% TODO: define the line function below.
[t, y] = runge_kutta(@line, 0, [t0 tfinal], 0.1);

%Compare with the integral of 2*t
%error = y-t.^2';
%mean(error)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use Runge-Kutta to integrate a second order equation of the form.
% d2y/dt^2+5*dy/dt + 4*y(t) = 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function exerciseB()
t0 = 0;
tfinal = 10;
[t, y] = runge_kutta(@second_order_system, [10 1], [t0 tfinal], 0.1);

%plot results
plot(t, y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now use Runge-kutta to integrate the movement of a 1 dof robot arm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function exerciseC()
%these variables are shared by the forward_dynamic_robot1 function defined
%below
global robot tau g
t0 = 0;
tfinal = 10;
robot = load_robot('example','1dofplanar')
robot.dynamics.friction=0
tau = [0];
g = [0 -9.81 0]';
[t, y] = runge_kutta(@forward_dynamic_robot1, [0 0]', [t0 tfinal], 0.01);

% Animate the movement. Change speed from 1-30-100-200
speed = 10
animate(robot,[y(1,1:speed:length(y))])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now use Runge-Kutta to simulate the movement of a 2 DOF robot arm.
function exerciseD()
global robot tau g
t0 = 0;
tfinal = 10;
robot = load_robot('example','2dofplanar')
robot.dynamics.friction=1
tau = [0 0];
g = [0 -9.81 0]';
[t, y] = runge_kutta(@forward_dynamic_robot2, [0 0 0 0]', [t0 tfinal], 0.001);
%speed, change to 10, 30, 50, 100
speed = 50
animate(robot,[y(1,1:speed:length(y)); y(2,1:speed:length(y))])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function.
% The function returns dy/dt = 2*t. Function called from exerciseA()
%
% Integrate a line in time. dy/dt = 2*t. Obviously, the integration should
% yield. y(t) = t^2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy = line(t, y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function to solve a second order differential equation.
% Called from function exerciseB()
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xd = second_order_system(t, y)

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function to simulate the movement of a 1 DOF robot
% Called from function exerciseC()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xd = forward_dynamic_robot1(t, y)
global tau g robot

qdd = accel(robot, y(1), y(2), tau, g);
%return qd, qdd
xd = [y(2); qdd];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function to simulate the movement of a 2 DOF robot
% Called from function exerciseD()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xd = forward_dynamic_robot2(t, y)
global tau g robot
%we must return the solution of
% [dx1/dt; dx2/dt]
t
qdd=forwarddynamics_2dofplanar(robot, y(1:2,1), y(3:4,1), tau', 9.81, [0 0 0 0 0 0]);
%return qd, qdd
xd = [y(3:4,1); qdd];

repo: ARTE


text:
% add the path of RBM code
addpath('..');
addpath('~/work/Algorithms/liblinear-1.7/matlab');

% load MNIST
load 'mnist_14x14.mat';

% shuffle the training data
X_labels = X_labels + 1;
X_test_labels = X_test_labels + 1;

perm_idx = randperm (size(X,1));

n_all = size(X, 1);
n_train = ceil(n_all * 3 / 4);
n_valid = floor(n_all /4);

X_valid = X(perm_idx(n_train+1:end), :);
X_valid_labels = X_labels(perm_idx(n_train+1:end));
X = X(perm_idx(1:n_train), :);
X_labels = X_labels(perm_idx(1:n_train));

layers = [size(X,2), 500, 500, 10];
n_layers = length(layers);
blayers = [1, 1, 1, 1];

use_tanh = 0;
do_pretrain = 1;

if do_pretrain
    Ds = cell(n_layers - 2, 1);
    H = X;
    H_valid = X_valid;

    for l = 1:n_layers-2
        % construct DAE and use default configurations
        D = default_dae (layers(l), layers(l+1));

        D.data.binary = blayers(l);
        D.hidden.binary = blayers(l+1);

        if use_tanh
            if l > 1
                D.visible.use_tanh = 1;
            end
            D.hidden.use_tanh = 1;
        else
            if D.data.binary
                mH = mean(H, 1)';
                D.vbias = min(max(log(mH./(1 - mH)), -4), 4);
            else
                D.vbias = mean(H, 1)';
            end
        end

        D.learning.lrate = 1e-1;
        D.learning.lrate0 = 5000;
        D.learning.minibatch_sz = 128;

        D.noise.drop = 0.2;
        D.noise.level = 0;

        %D.adagrad.use = 1;
        %D.adagrad.epsilon = 1e-8;
        D.adadelta.use = 1;
        D.adadelta.epsilon = 1e-8;
        D.adadelta.momentum = 0.99;

        D.valid_min_epochs = 10;

        if blayers(l+1)
            D.cae.cost = 0.01;
            %D.sparsity.target = 0.1;
            %D.sparsity.cost = 0.01;
        end

        D.iteration.n_epochs = 500;

        % save the intermediate data after every epoch
        D.hook.per_epoch = {@save_intermediate, {sprintf('dae_mnist_%d.mat', l)}};

        % print learining process
        D.verbose = 0;
        % display the progress
        D.debug.do_display = 0;

        % train RBM
        fprintf(1, 'Training DAE (%d)\n', l);
        tic;
        D = dae (D, H, H_valid, 0.1);
        fprintf(1, 'Training is done after %f seconds\n', toc);

        H = dae_get_hidden(H, D);
        H_valid = dae_get_hidden(H_valid, D);

        Ds{l} = D;
    end
end

M = default_mlp (layers);

M.output.binary = blayers(end);
M.hidden.use_tanh = use_tanh;

M.valid_min_epochs = 10;
M.dropout.use = 1;

M.hook.per_epoch = {@save_intermediate, {'mlp_mnist.mat'}};

M.learning.lrate = 1e-3;
M.learning.lrate0 = 5000;
M.learning.minibatch_sz = 128;

M.adadelta.use = 1;
M.adadelta.epsilon = 1e-8;
M.adadelta.momentum = 0.99;

M.noise.drop = 0;
M.noise.level = 0;

M.iteration.n_epochs = 100;

if do_pretrain
    for l = 1:n_layers-2
        M.biases{l+1} = Ds{l}.hbias;
        M.W{l} = Ds{l}.W;
    end
end

fprintf(1, 'Training MLP\n');
tic;
M = mlp (M, X, X_labels, X_valid, X_valid_labels, 0.1);
fprintf(1, 'Training is done after %f seconds\n', toc);

[pred] = mlp_classify (M, X_test);
n_correct = sum(X_test_labels == pred);

fprintf(2, 'Correctly classified test samples: %d/%d\n', n_correct, size(X_test, 1));

author: kyunghyuncho

repo: deepmat


text:
function [log_rad,angle]=base(m,n)
% Compute radius and angle matrix base for a given size of m,n
% Input:
%       m:  number of rows
%       n:  number of cols
% Output:
%   log_rad:    radius matrix, store the distance between a pixel and
%               the center of matrix
%     angle:    angle matrix, store the angle of a pixel

ctrm=ceil(m/2);
ctrn=ceil(n/2);

% x span from -m/2 to m/2 with interval 2/m
[x,y] = meshgrid(( (1:m)-ctrm)/(m/2),( (1:n)-ctrn)/(n/2));

% radius matrix
rad = sqrt(x.^2 + y.^2);
rad(m/2,n/2) =  rad(m/2,n/2-1);
log_rad  = log2(rad);

% angle matrix
angle = atan2(y,x);

author: Sable

repo: mcbench-benchmarks


text:
function modelNew = mergeModelFieldPositions(model,type,positions,mergeFunctions)
% USAGE:
%    [modelNew] = mergeModelFieldPositions(model,type,positions)
%
% INPUTS:
%    model:           The model with the fields to merge
%    type:            the field type to merge ( rxns, mets, comps or genes)
%    positions:       The positions in the given field type to merge either as indices or as logical array.
%
% OPTIONAL INPUTS:
%    mergeFunctions:  A cell array of fieldNames and functions that can be
%                     called on a larger array of the type used for the
%                     field. e.g.
%
%                      - {'metCharges',@(x) x(1); 'grRules', @(x) strjoin(x,' or ');
%                      - by default, all numeric fields are added up, all
%                        unique entries in cell arrays are concatenated with
%                        ';',
%                      - grRules and rules are assumed to be merged with
%                        and (i.e. a batch reaction is assumed to need all
%                        associated GPRs)
%
%
% OUTPUT:
%
%    modelNew:         merged model with the positions merged into one.
%                      Unique entries in Related String fields will be
%                      concatenated with ';' matrices will be summed up.
%
% .. Authors:
%                   - Thomas Pfau Sept 2017

modelNew = model;

if islogical(positions)
    positions = find(positions);
end

if numel(positions) <= 1
    %If there is less than two positions to merge, we don't do anything.
    return;
end

%basic functions for merger
basicFunctions = {'rules',@(x) {strjoin(cellfun(@(x) strcat('(',x,')'), setdiff(x,''),'UniformOutput',false),' & ')};...
                  'grRules', @(x) {strjoin(cellfun(@(x) strcat('(',x,')'), setdiff(x,''),'UniformOutput',false),' and ')}};
if ~exist('mergeFunctions','var')
    mergeFunctions = basicFunctions;
else
    toRemove = ismember(basicFunctions(:,1),mergeFunctions(:,1));
    if all(toRemove)
        mergeFunctions=mergeFunctions(:,:);
    else
        mergeFunctions=[basicFunctions(~toRemove,:);mergeFunctions(:,:)];
    end
end

posToKeep = positions(1);
posToMerge = positions(2:end);

%Store the original names.
origNames = modelNew.(type)(positions);

%get The fields and associated dimensions
[fields,dimensions] = getModelFieldsForType(modelNew, type, numel(model.(type)));

for i = 1:numel(fields)
    specialFieldPos = ismember(mergeFunctions(:,1),fields{i});
    if any(specialFieldPos)
        mergeFun = mergeFunctions{specialFieldPos,2};
        data = mergeFun(getSlice(modelNew.(fields{i}),positions,dimensions(i)));
        modelNew.(fields{i}) = setSlice(modelNew.(fields{i}),posToKeep,dimensions(i),data);
        continue;
    end

    %Lets assume, that we only have 2 dimensional fields.
    if isnumeric(modelNew.(fields{i})) || islogical(modelNew.(fields{i}))
        %There are exceptions.
        modelNew.(fields{i}) = setSlice(modelNew.(fields{i}),posToKeep,dimensions(i),sum(getSlice(modelNew.(fields{i}),positions,dimensions(i)),dimensions(i)));
    end
    % if its cell arrays, concatenate unique data, we will assume here, that
    % all cell arrays are string arrays and that there are no
    % multi-dimensional
    %There are two options for cell arrays currently. Either its all char
    %arrays, or its all cell arrays.
    if iscell(modelNew.(fields{i}))
        data = getSlice(modelNew.(fields{i}),positions,dimensions(i));
        if ischar([data{:}])
            newData = {strjoin(unique(data),';')};
            modelNew.(fields{i}) = setSlice(modelNew.(fields{i}),posToKeep,dimensions(i),newData);
        else
            dataconcatenation = unique([data{:}]);
            if iscell(dataconcatenation) % We have a field with Cell arrays of Cell arrays
                modelNew.(fields{i}) = setSlice(modelNew.(fields{i}),posToKeep,dimensions(i),{unique([data{:}])});
            else
                %Not sure which fields this would apply to, but I assume,
                %this will lead to an error....
                modelNew.(fields{i}) = setSlice(modelNew.(fields{i}),posToKeep,dimensions(i),unique([data{:}]));
            end
        end
    end
    if ischar(modelNew.(fields{i}))
        %This is a problem. as this is VERY dependent on the type of char
        %we have. We will simple not handle it for now... (as it is
        %currently only present in csense and dsense and whoever is doing it should
        %take care not to screw around when merging mets....
    end
    %now, having done this, we need to check whether we modified genes. If
    %so, we have to update the rules vector and the grRules vector, if
    %present
    if strcmp(type,'genes')
         if isfield(modelNew,'rules')
            for j = 1:numel(posToMerge)
                %Replace by new position.
                modelNew.rules = strrep(modelNew.rules,['x(' num2str(posToMerge(j)) ')'],['x(' num2str(posToKeep) ')']);
            end
         end
         %Update the grRules field (if necessary).
        if isfield(modelNew, 'grRules') && numel(unique(origNames)) > 1
            for j = 1:numel(positions)
                %Replace by new name.
                %First, replace all occurences, which only contain the
                %name.
                modelNew.grRules = regexprep(modelNew.grRules,['^' regexptranslate('escape',origNames{j}) '$'],modelNew.genes{posToKeep});
                %Then replace all occurences which are at the beginning or
                %end of a formula.
                modelNew.grRules = regexprep(modelNew.grRules,['^' regexptranslate('escape',origNames{j}) '([\) ]+)'],[ modelNew.genes{posToKeep} '$1']);
                modelNew.grRules = regexprep(modelNew.grRules,['([\( ]+)' regexptranslate('escape',origNames{j}) '$'],['$1' modelNew.genes{posToKeep} ]);
                %finally replace all  in the middle.
                modelNew.grRules = regexprep(modelNew.grRules,['([\( ]+)' regexptranslate('escape',origNames{j}) '([\) ]+)'],['$1' modelNew.genes{posToKeep} '$2']);
            end
        end
    end

end
%After merging remove the merged fields.

modelNew = removeFieldEntriesForType(modelNew,posToMerge,type,numel(modelNew.(type)));

end

function out = getSlice(A,idx,dim)
%Get a slice from a given matrix/vector
% USAGE:
%    out = getSlice(A,idx,dim)
%
% INPUTS:
%    A:               The matrix/vector to obtain a slice from
%    idx:             the indices in the given dimension to extract
%    dim:             The dimension to obtain the given indices from.
%
%
% OUTPUT:
%
%    out:             A slice from the matix/vector with the given indices
%
% Note:
%     Based on https://stackoverflow.com/questions/22537326/on-shape-agnostic-slicing-of-ndarrays
%
% .. Authors:
%                   - Thomas Pfau Sept 2017

slice = repmat({':'},1, ndims(A));
slice{dim} = idx;
out = A(slice{:});
end

function A = setSlice(A,idx,dim,B)
%Set a slice from a given matrix/vector
% USAGE:
%    A = setSlice(A,idx,dim,B)
%
% INPUTS:
%    A:               The matrix/vector to obtain a slice from
%    idx:             the indices in the given dimension to extract
%    dim:             The dimension to obtain the given indices from.
%    B:               The values to set at the given position
%
% OUTPUT:
%    A:             A slice from the matix/vector with the given indices
%
% Note:
%     Based on https://stackoverflow.com/questions/22537326/on-shape-agnostic-slicing-of-ndarrays
%
% .. Authors:
%                   - Thomas Pfau Sept 2017

slice = repmat({':'},1, ndims(A));
slice{dim} = idx;
%If B is a cell array of cell arrays
A(slice{:}) = B;
end

author: opencobra

repo: cobratoolbox


text:
function folder = layoutRoot()
%layoutRoot  Folder containing the GUI Layout Toolbox
%
%   folder = layoutRoot() returns the full path to the folder containing
%   the GUI Layout Toolbox.
%
%   Examples:
%   >> folder = layoutRoot()
%   folder = 'C:\tools\layouts2\layout'
%
%   See also: layoutVersion

%  Copyright 2009-2014 The MathWorks, Inc.
%  $Revision: 980 $ $Date: 2014-09-28 14:27:26 -0400 (Sun, 28 Sep 2014) $

folder = fileparts( mfilename( 'fullpath' ) );

end % layoutRoot

author: goGPS-Project


text:
function [ out ] = EMD1D( u, v )
%UNTITLED3 Summary of this function goes here
%   Detailed explanation goes here
N = length(u);
cumDif = zeros(N+1,1);

for i = 1:N
    cumDif(i+1) = cumDif(i) + u(i) - v(i);
end

cumDif = abs(cumDif);
out = sum(cumDif);

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function varargout = ft_read_vol(varargin)

% This function is a backward compatibility wrapper for existing MATLAB scripts
% that call a function that is not part of the FieldTrip toolbox any more.
%
% Please update your code to make it future-proof.

oldname = mfilename;
newname = 'ft_read_headmodel';

ft_warning('%s is only a backward compatibility wrapper, which will soon be removed. Please call %s instead.', upper(oldname), upper(newname));

[varargout{1:nargout}] = feval(newname, varargin{:});


text:
function [cij,flag] = makerandCIJdegreesfixed(in,out)
%MAKERANDCIJDEGREESFIXED        Synthetic directed random network
%
%   CIJ = makerandCIJdegreesfixed(N,K);
%
%   This function generates a directed random network with a specified
%   in-degree and out-degree sequence. The function returns a flag,
%   denoting whether the algorithm succeeded or failed.
%
%   Inputs:     in,     indegree vector
%               out,    outdegree vector
%
%   Output:     CIJ,    binary directed connectivity matrix
%               flag,   flag=1 if the algorithm succeeded; flag=0 otherwise
%
%
%   Notes:  Necessary conditions include:
%               length(in) = length(out) = n
%               sum(in) = sum(out) = k
%               in(i), out(i) < n-1
%               in(i) + out(j) < n+2
%               in(i) + out(i) < n
%
%           No connections are placed on the main diagonal
%
%
% Aviad Rubinstein, Indiana University 2005/2007

% intialize
n = length(in);
k = sum(in);
inInv = zeros(k,1);
outInv = inInv;
iIn = 1; iOut = 1;

for i = 1:n
    inInv(iIn:iIn+in(i) - 1) = i;
    outInv(iOut:iOut+out(i) - 1) = i;
    iIn = iIn+in(i);
    iOut = iOut+out(i);
end

cij = eye(n);
edges = [outInv(1:k)'; inInv(randperm(k))'];

% create cij, and check for double edges and self-connections
for i = 1:k
    if cij(edges(1,i),edges(2,i))
        warningCounter = 1;
        while (1)
            switchTo = ceil(k*rand);
            if ~(cij(edges(1,i),edges(2,switchTo)) || cij(edges(1,switchTo),edges(2,i)))
                cij(edges(1,i),edges(2,switchTo)) = 1;
                if switchTo < i
                    cij(edges(1,switchTo),edges(2,switchTo)) = 0;
                    cij(edges(1,switchTo),edges(2,i)) = 1;
                end
                temp = edges(2,i);
                edges(2,i) = edges(2,switchTo);
                edges(2,switchTo) = temp;
                break
            end
            warningCounter = warningCounter+1;
            % If there is a legitimate subtitution, it has a probability of 1/k of being done.
            % Thus it is highly unlikely that it will not be done after 2*k^2 attempts.
            % This is an indication that the given indegree / outdegree
            % vectors may not be possible.
            if warningCounter == 2*k^2
                flag = 0;  % no valid solution found
                return;
            end
        end
    else
        cij(edges(1,i),edges(2,i)) = 1;
    end
end

cij = cij - eye(n);

% a valid solution was found
flag = 1;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [p,n]=numSubplots(n)
% function [p,n]=numSubplots(n)
%
% Purpose
% Calculate how many rows and columns of sub-plots are needed to
% neatly display n subplots.
%
% Inputs
% n - the desired number of subplots.
%
% Outputs
% p - a vector length 2 defining the number of rows and number of
%     columns required to show n plots.
% [ n - the current number of subplots. This output is used only by
%       this function for a recursive call.]
%
%
%
% Example: neatly lay out 13 sub-plots
% >> p=numSubplots(13)
% p =
%     3   5
% for i=1:13; subplot(p(1),p(2),i), pcolor(rand(10)), end
%
%
% Rob Campbell - January 2010

while isprime(n) & n>4,
    n=n+1;
end

p=factor(n);

if length(p)==1
    p=[1,p];
    return
end

while length(p)>2
    if length(p)>=4
        p(1)=p(1)*p(end-1);
        p(2)=p(2)*p(end);
        p(end-1:end)=[];
    else
        p(1)=p(1)*p(2);
        p(2)=[];
    end
    p=sort(p);
end

%Reformat if the column/row ratio is too large: we want a roughly
%square design
while p(2)/p(1)>2.5
    N=n+1;
    [p,n]=numSubplots(N); %Recursive!
end

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function y = gpSubspaceOut(model,x)

% GPSUBSPACEOUT
%
% COPYRIGHT : Carl Henrik Ek, 2008

% GP

y = NaN.*ones(size(x,1),length(model.dim));
y(:,find(model.dim)) = gpOut(model,x);

return;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   MoveJ:
%   Caution. The global variable robot should be loaded
%   with a valid robot arm.
%   Example: MoveJ(pos_b_rec, 'v1000', 'z100', gripper, wobj0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function MoveJ(robtarget, speeddata, zonedata, gripper, Wobj)

global configuration robot

fprintf('\nCall to MoveJ');

%obtain current joint coordinates
q_current = robot.q;

%Ttool, transformation from the robot's end to the TCP
Ttool = transform_to_homogeneous(gripper(2:end));

%compute target joint coordinates, from the inverse kinematics
%select desired joint values by using conf data from robtarget
%select configuration
conf = get_conf_data(robtarget);
Ttotal = transform_to_homogeneous(robtarget);

%T = Ttotal*inv(Ttool);
T = Ttotal/Ttool;

%Obtain solutions from inversekinematic function
q_final=inversekinematic(robot, T, robot.q);

%select one of the solutions of the inverse kinematic problem
% solution according to configuration and language.
% The conf vector defines the quadrants of the joints [theta1 theta4 theta6 thetax]
% =[cf1 cf4 cf6 cfx]. Where the x angle is zero in this implementation.
q_final=select_configuration(robot, q_final, conf);
%(q_final*180/pi)'

speed = obtain_joint_speed(robot, speeddata);
%when fine is specified, the radius is zero, otherwise,
%the radius represents a sphere around the target where the
%movement is changed to the next target
radius = obtain_zone_data(zonedata);

%find time to perform movement
[actual_speed, tmax]=synchronize(q_current, q_final, speed, robot.accelmax);

%in this case, if the total time for the movement is lower than the
%configuration.delta_time time then 10 points are interpolated
if tmax <= configuration.delta_time
    t = [0:tmax/10:tmax]';
else
    %local time for the planning, normal case
    t = [0:configuration.delta_time:tmax]';
end

if radius==0
    %compute a smooth trajectory in q and qd
    %in this case, the final speed is zero
    [q, qd, qdd] = compute_joint_trajectory_indep(robot, q_current, q_final, robot.qd, zeros(1,robot.DOF), t);
else
    %in this case, the final speed is selected as the speed selected by the
    %variable speeddata
    [q, qd, qdd] = compute_joint_trajectory_indep(robot, q_current, q_final, robot.qd, actual_speed, t);
end
%finds the first joint coordinates that are inside a radius r of the target
%point
index = find_first_in_zone_data(robot, q, T, radius);

%the robot performs the movement until the index found. The coordinates, joint speed and acceleratin
%are stored and used in the planning of the next point
robot.q=q(:,index);
robot.qd=qd(:,index);
robot.qdd=qdd(:,index);

%store all the trajectory for plotting
%the joint trajectories, speeds and acceleration of susequent movements are
%store here
%robot.q_vector=[robot.q_vector q(:, 1:index)];
%robot.qd_vector=[robot.qd_vector qd(:, 1:index)];
%robot.qdd_vector=[robot.qdd_vector qdd(:, 1:index)];

robot.q_vector = q(:, 1:index);
robot.qd_vector = qd(:, 1:index);
robot.qdd_vector = qdd(:, 1:index);

%a global time for the planning is computed.
%in this way, the total trajectory of different movements can be plotted
if length(robot.time)==0
    tend = 0;
else
    tend = robot.time(end);
end
t = t + tend;
%store total time
robot.time=t(1:index)';

%Test whether there are joints outside mechanical limits
test_joint_limits(robot);

%Plot position, velocity and acceleration
%plot_joint_data(robot);
%Now, animate the robot in 3D
animate(robot, q);

%Plot the trajectory
pp = zeros(3,size(q,2));
for i=1:size(q,2),
    T=directkinematic(robot, q(:,i));
    %pp=[pp T(1:3,4)];
    pp(1:3,i)=T(1:3,4);
end
plot3(pp(1,:),pp(2,:),pp(3,:), 'k', 'LineWidth', 2)


text:
function [varargout] = rdir(rootdir,varargin)
% RDIR - Recursive directory listing
%
%  D = rdir(ROOT)
%  D = rdir(ROOT, TEST)
%  D = rdir(ROOT, TEST, RMPATH)
%  D = rdir(ROOT, TEST, 1)
%  D = rdir(ROOT, '', ...)
%  [D, P] = rdir(...)
%  rdir(...)
%
%
% *Inputs*
%
% * ROOT
%
% rdir(ROOT) lists the specified files.
% ROOT can be a pathname, filename, or can include both. One can use
% absolute and relative pathnames and wildcards (*). Wildcard can be placed
% anywhere and used many times like 'path*\*.m'
%
% One can also use a double wildcard (**) to match multiple directory
% levels. For example ROOT = 'path\**\*.m' will match all ".m" files in
% "path" and all subdirectories of "path".
%
% NOTE : ".svn" directories created by SubVersion (SVN) are excluded from
% the recursive listing.
%
% * TEST
%
% Optional test that can be performed on the returned files.
%
% TEST is a string indicating expression to be evaluated on selected field
% of rdir output.
% All fields (ie name, date, bytes, isdir and datenum) can be used.
%
% Tests are strings similar to what one would use in a "if" statement e.g.
%  'bytes>1024 & datenum>now-7'
%
% One can also use function like "regexp" or "strfind" with string fields
% like "name" and "date" e.g 'regexp(name, 'expr')'. In that case, tests
% that return a non empty value are considered as true.
%
% regexp(name, '(\.m$)|(\.mdl$)')
%
% Test can also be a function handle as used in arrayfun/cellfun, e.g.
% @(f)f.bytes>1024
%
% * RMPATH
%
% Optional path to remove from beginning of "name" field in returned
% output. Specified path must be common to all items found.
%
% If RMPATH = 1 or true, path to remove is part of ROOT before the first
% wildcard.
%
%
% *Outputs*
%
% * D
%
% D is a structure with the same fields as Matlab DIR output.
%
% The "name" field includes the relative path as well as the name to the
% file that was found. Path can be shorten or ommited when using 3rd
% argument RMPATH.
%
% * P
%
% Common path or RMPATH (if specified) for the file list returned in D.
%
% * Screen output
%
% If not output variable is specified then the output is sent to the
% screen.
%
%
% *Versions*
%
% * 1.0 - 2009, Gus Brown
% * 2.0 - 26/05/2011 Thomas Vanaret
%         No longer exclude all directories from a simple search (no *);
%         Fixing bug on returned path;
%         Exclude ".svn" directories;
%         Extended test possibilies;
%         Subfunctions created;
% * 2.1 - 14/07/2011 Thomas Vanaret
%         New argument allowing to remove common path from name;
%         Comments review;
% * 2.2 - 20/12/2011 Thomas Vanaret
%         Fixing bug on display with 0b files;
%         Specific display when no file match filter;
% * 2.3 - 19/01/2014 Thomas Vanaret
%         Adding improvements suggested by X. Mo :
%         - function handle as TEST input
%         - code optimisation (avoiding loop)
%         Fixing possible bug when using a wildcard at the beginning;
%         Common path as 2nd optionnal output;
%
%
% *Examples*
%
%   D = rdir('*.m');
%     for ii=1:length(D), disp(D(ii).name); end;
%
%   % to find all files in the current directory and sub directories
%   D = rdir('**\*')
%
%   % If no output is specified then the files are sent to
%   % the screen.
%   rdir('c:\program files\windows *\*.exe');
%   rdir('c:\program files\windows *\**\*.dll');
%
%   % Using the test function to find files modified today
%   rdir('c:\win*\*','datenum>floor(now)');
%   % Using the test function to find files of a certain size
%   rdir('c:\program files\win*\*.exe','bytes>1024 & bytes<1048576');
%   % Using the test function to find files modified in 2011
%   rdir('c:\win*\*','strfind(date, ''2011'')');
%
%   % Using the 3rd input to shorten output name
%   rdir([matlabroot, '\*.txt'], '', 'C:\Program Files\')
%   % Using the 3rd input to shorten output name
%   rdir([matlabroot, '\*.txt'], '', 1)
%
%
% See also DIR
%

%--------------------------------------------------------------------------
%% Input validation

% use the current directory if nothing is specified
if ~exist('rootdir','var'),
  rootdir = '*';
end

prepath = '';       % the path before the wild card
wildpath = '';      % the path wild card
postpath = rootdir; % the path after the wild card
I = find(rootdir==filesep,1,'last');

% Directory separator for current platform
if filesep == '\'
  % On PC, filesep is '\'
  anti_filesep = '/';
else
  % On UNIX system, filesep is '/'
  anti_filesep = '\';
end

if isempty(I) && ~isempty(strfind(rootdir, anti_filesep))
  error([mfilename, ':FileSep'],...
    'Use correct directory separator "%s".', filesep)
end

%--------------------------------------------------------------------------
%% Split rootdir
% split the file path around the wild card specifiers

if ~isempty(I),
  prepath = rootdir(1:I);
  postpath = rootdir(I+1:end);
  I = find(prepath=='*',1,'first');
  if ~isempty(I),
    postpath = [prepath(I:end) postpath];
    prepath = prepath(1:I-1);
    I = find(prepath==filesep,1,'last');
    if ~isempty(I),
      wildpath = prepath(I+1:end);
      prepath = prepath(1:I);
    end;
    I = find(postpath==filesep,1,'first');
    if ~isempty(I),
      wildpath = [wildpath postpath(1:I-1)];
      postpath = postpath(I:end);
    end;
  end;
end;

% disp([' "' prepath '" ~ "' wildpath '" ~ "' postpath '" ']);

%--------------------------------------------------------------------------
%% Recursive listing
% Search for matching files until all wildcards have been considered.

if isempty(wildpath)
  % If no directory wildcards then just get files and directories list

  D = dir([prepath postpath]);

  % Exclude ".", ".." and ".svn" directories from the list
  excl = isdotdir(D) | issvndir(D);
  D(excl) = [];

  if isdir([prepath postpath]);
    fullpath = [prepath postpath];
  else
    fullpath = prepath;
  end

  % Place directories on the top of the list
  is_dir = [D.isdir]';
  D = [D(is_dir); D(~is_dir)];

  % Add path before name
  for ii = 1:length(D)
    D(ii).name = fullfile(fullpath, D(ii).name);
  end

  % disp(sprintf('Scanning "%s"   %g files found',[prepath postpath],length(D)));

elseif strcmp(wildpath,'**')
  % A double wildcards directory means recurs down into sub directories

  % first look for files in the current directory (remove extra filesep)
  D = rdir([prepath postpath(2:end)]);

  % then look for sub directories
  D_sd = dir([prepath '*']);

  % Exclude ".", "..", ".svn" directories and files from the list
  excl = isdotdir(D_sd) | issvndir(D_sd) | ~([D_sd.isdir]');
  D_sd(excl) = [];

  % Process each sub directory found
  % Performance tweak: avoid growing array within loop (X. Mo)
  c_D = arrayfun(@(x) rdir([prepath x.name filesep wildpath postpath]),...
    D_sd, 'UniformOutput', false);

  D = [D; cell2mat( c_D ) ];

else
  % Process directory wild card looking for sub directories that match

  D_sd = dir([prepath wildpath]);

  % Exclude ".", "..", ".svn" directories and files from the list
  excl = isdotdir(D_sd) | issvndir(D_sd) | ~([D_sd.isdir]');
  D_sd(excl) = [];

  if ~isdir(prepath) || ( numel(D_sd)==1 && strcmp(D_sd.name, prepath))
    % Fix case like rdir('path*\...') where prepath is not a full directoty
    % name OR case were prepath match a unique directory.
    % Previous "dir" return then the matching directory name(s).
    % prepath is cleaned to use them.
    %
    % In else case, prepath is a valid path which must be kept.
    prepath = '';
  end

  % Process each directory found
  Dt = dir('');

  c_D = arrayfun(@(x) rdir([prepath x.name postpath]),...
    D_sd, 'UniformOutput', false);

  D = [Dt; cell2mat( c_D ) ];

end

%--------------------------------------------------------------------------
%% Apply filter
% If specified, apply the filter to refine the search.

nb_before_filt = length(D);
warning_msg = '';

if (nargin>=2 && ~isempty(varargin{1})),
  try
    if isa(varargin{1}, 'function_handle')
        test_tf = arrayfun(varargin{1}, D);
    else
        test_tf = evaluate(D, varargin{1});
    end

    D = D(test_tf);

  catch
    if isa(varargin{1}, 'function_handle')
      test_expr = func2str(varargin{1});
    else
      test_expr = varargin{1};
    end

    warning_msg = sprintf('Invalid TEST "%s" : %s', test_expr, lasterr);
  end
end

%--------------------------------------------------------------------------
%% Remove path
% If specified, remove given or common path from each returned path.

common_path = '';
if (nargin>=3 && ~isempty(varargin{2})),

  arg2 = varargin{2};
  if ischar(arg2)
    common_path = arg2;
  elseif (isnumeric(arg2) || islogical(arg2)) && arg2
    common_path = prepath;
  end

  rm_path = regexptranslate('escape', common_path);

  % Check that path is common to all
  start = regexp({D.name}', ['^', rm_path]);

  % Convert to a logical.
  is_common = not( cellfun(@isempty, start) );

  if all(is_common)
    for k = 1:length(D)
      D(k).name = regexprep(D(k).name, ['^', rm_path], '');
    end

  else
    common_path = '';
  end

  % 19/07/2012 : ajouter common_path en sortie optionnelle

end

%--------------------------------------------------------------------------
%% Display listing if no output variables are specified
% Screen display.

nout = nargout;

if nout == 0
  if isempty(D)
    if nb_before_filt == 0
      fprintf('%s not found.\n', rootdir)
    else
      fprintf('No item matching filter.\n')
    end
  else

    if ~isempty(common_path)
     fprintf('All in : %s\n', common_path)
    end

    pp = {'' 'k' 'M' 'G' 'T'};
    for ii = 1:length(D)
      if D(ii).isdir
        % Directory item : display name
        disp(sprintf(' %29s %-64s','',D(ii).name));
      else
        % File item : display size, modification date and name
        sz = D(ii).bytes;
        if sz > 0
          ss = min(4,floor(log2(sz)/10));
        else
          ss = 0;
        end
        disp(sprintf('%4.0f %1sb  %20s  %-64s ',...
          sz/1024^ss, pp{ss+1}, datestr(D(ii).datenum, 0), D(ii).name));
      end
    end
  end
elseif nout == 1
  % send list out
  varargout{1} = D;
else
  % send list and common path out
  varargout{1} = D;
  varargout{2} = common_path;
end;

if ~isempty(warning_msg)
  warning([mfilename, ':InvalidTest'],...
    warning_msg); % ap aff
end

%---------------------------- end of main function ------------------------

%% ------------------------------------------------------------------------
function tf = issvndir(d)
% True for ".svn" directories.
% d is a structure returned by "dir"
%

is_dir = [d.isdir]';

is_svn = strcmp({d.name}, '.svn')';
%is_svn = false; % uncomment to disable ".svn" filtering

tf = (is_dir & is_svn);

%---------------------------- end of subfunction --------------------------

%% ------------------------------------------------------------------------
function tf = isdotdir(d)
% True for "." and ".." directories.
% d is a structure returned by "dir"
%

is_dir = [d.isdir]';

is_dot = strcmp({d.name}, '.')';
is_dotdot = strcmp({d.name}, '..')';

tf = (is_dir & (is_dot | is_dotdot) );

%---------------------------- end of subfunction --------------------------

%% ------------------------------------------------------------------------
function tf = evaluate(d, expr)
% True for item where evaluated expression is correct or return a non empty
% cell.
% d is a structure returned by "dir"
%

% Get fields that can be used
name = {d.name}'; %#ok<NASGU>
date = {d.date}'; %#ok<NASGU>
datenum = [d.datenum]'; %#ok<NASGU>
bytes = [d.bytes]'; %#ok<NASGU>
isdir = [d.isdir]'; %#ok<NASGU>

tf = eval(expr); % low risk since done in a dedicated subfunction.

% Convert cell outputs returned by "strfind" or "regexp" filters to a
% logical.
if iscell(tf)
  tf = not( cellfun(@isempty, tf) );
end

%---------------------------- end of subfunction --------------------------

%---------------------------- END OF FUNCTION -----------------------------

author: peiyunh

repo: tiny

save_path: github-repos/MATLAB/peiyunh-tiny


text:
function RE = estimate_registration_error(fixedCoordinates,movingCoordinates,displacementField,spacing)
% ESTIMATE_REGISTRATION_ERROR Estimates the target registration error
%
% RE = estimate_registration_error(...
%   fixedCoordinates,movingCoordinates,displacementField,spacing)
%
% INPUT ARGUMENTS
% fixedCoordinates      - Coordinates in fixed image
% movingCoordinates     - Coordinates in moving image
% displacementField     - Displacement field, pointing from fixed to moving
% spacing               - Pixel spacing
%
% OPTIONAL INPUT ARGUMENTS
% N/A
%
% OUTPUT ARGUMENTS
% RE                    - Registration error

% Copyright (c) 2012 Daniel Forsberg
% danne.forsberg@outlook.com
%
% 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 this program.  If not, see <http://www.gnu.org/licenses/>.

dims = length(displacementField);
deformedFixedCoordinates = fixedCoordinates;
switch dims
    case 1
        xDisp = interp3(displacementField{1},fixedCoordinates(:,1));
        deformedFixedCoordinates(:,1) = fixedCoordinates(:,1) + xDisp;
        RE = deformedFixedCoordinates - movingCoordinates;
        RE = sqrt((RE(:,1)*spacing(1)).^2);
    case 2
        xDisp = interp3(displacementField{1},fixedCoordinates(:,1),fixedCoordinates(:,2));
        yDisp = interp3(displacementField{2},fixedCoordinates(:,1),fixedCoordinates(:,2));
        deformedFixedCoordinates(:,1) = fixedCoordinates(:,1) + xDisp;
        deformedFixedCoordinates(:,2) = fixedCoordinates(:,2) + yDisp;
        RE = deformedFixedCoordinates - movingCoordinates;
        RE = sqrt((RE(:,1).^spacing(1)).^2 + (RE(:,2)*spacing(2)).^2);
    case 3
        xDisp = interp3(displacementField{1},fixedCoordinates(:,1),fixedCoordinates(:,2),fixedCoordinates(:,3));
        yDisp = interp3(displacementField{2},fixedCoordinates(:,1),fixedCoordinates(:,2),fixedCoordinates(:,3));
        zDisp = interp3(displacementField{3},fixedCoordinates(:,1),fixedCoordinates(:,2),fixedCoordinates(:,3));
        deformedFixedCoordinates(:,1) = fixedCoordinates(:,1) + xDisp;
        deformedFixedCoordinates(:,2) = fixedCoordinates(:,2) + yDisp;
        deformedFixedCoordinates(:,3) = fixedCoordinates(:,3) + zDisp;
        RE = deformedFixedCoordinates - movingCoordinates;
        RE = sqrt((RE(:,1)*spacing(1)).^2 + (RE(:,2)*spacing(2)).^2 + (RE(:,3)*spacing(3)).^2);
    otherwise
        error('Not implemented for dimensions higher than three.')
end

author: fordanic

repo: image-registration


text:
%filename='test2d_micro';
filename = 'MicroQuad';
ptype = 'MICRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'circleInclusion';
cost={'chomog_alphabeta'};
weights=[1];
constraint = {'volumeConstraint'};
constraint_case = 'EQUALITY';
%incrementFactor = 1;
%designVariable = 'Density';
designVariable = 'LevelSet';
filterType = 'P1';
fracRadius = 0.51;
%optimizer = 'IPOPT';
optimizer = 'DualNestedInPrimal';
%optimizer = 'AlternatingPrimalDual';

%optimizerUnconstrained = 'PROJECTED GRADIENT';
%line_search_initiator = 'INCREASING LAST STEP';
%incrementFactor = 1.95;

optimizerUnconstrained = 'SLERP';

nsteps = 1;
Vfrac_final = 0.4;
Perimeter_target=1;
optimality_final = 1*1e-3;
constr_final =1e-12;

Vfrac_initial = 0.8;
optimality_initial = 1*1e-3;
constr_initial = 1e-12;

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 = true;
printing = true;
monitoring = true;
monitoring_interval = 1;
maxiter = 100;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ uv ] = coords2uv( coords, width, height )
%COORDS2UV Image coordinates (xy) to uv
%   Detailed explanation goes here
middleX = width/2+0.5;
middleY = height/2+0.5;
uv = [(coords(:,1)-middleX)./width*2*pi -(coords(:,2)-middleY)./height*pi];

end

author: zouchuhang

repo: LayoutNet


text:
function M = specialeuclideanfactory(n, k)
% Returns a manifold structure to optimize over the special Euclidean group
%
% function M = specialeuclideanfactory(n)
% function M = specialeuclideanfactory(n, k)
%
% The special Euclidean group (the manifold of rigid transformations):
% This is a product manifold of the rotations group SO(n) and the
% translation group R^n, copied k times.
%
% Points on the manifold are represented as structures X with two fields.
% X.R is a 3D array of size nxnxk such that each slice X.R(:, :, i)
% corresponds to a rotation matrix (orthogonal with determinant 1).
% X.t is a matrix of size nxk such that each column X.t(:, i) corresponds
% to a translation vector.
%
% Tangent vectors are represented as structures with the same fields. Note
% that rotational components of the tangent vectors are represented in the
% Lie algebra, i.e., each slice Xdot.R(:, :, i) is a skew-symmetric matrix.
% Use M.tangent2ambient(X, Xdot) to obtain a representation in the ambient
% space.
%
% This is a description of SE(n)^k with the induced metric from the
% embedding space (R^nxn)^k x (R^n)^k, i.e., this manifold is a Riemannian
% submanifold of the embedding Euclidean space with the usual inner
% product.
%
% By default, k = 1.
%
% This is a test geometry: it may not be the "appropriate" geometry to give
% to SE(n).
%
% See rotationsfactory and euclideanfactory for details.
%
% See also: rotationsfactory euclideanfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Sep. 23, 2014.
% Contributors:
% Change log:

    if ~exist('k', 'var') || isempty(k)
        k = 1;
    end

    elements = struct();
    elements.R = rotationsfactory(n, k);
    elements.t = euclideanfactory(n, k);

    M = productmanifold(elements);

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function [varargout] = plinproj(varargin)

% PLINPROJ projects a point onto a line or linepiece
%
% Use as
%   [proj, dist] = plinproj(l1, l2, r, flag)
% where l1 and l2 are the begin and endpoint of the linepiece, and r is
% the point that is projected onto the line
%
% the optional flag can be:
%   0 (default)  project the point anywhere on the complete line
%   1            project the point within or on the edge of the linepiece

% Copyright (C) 2002-2009, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% compile the missing mex file on the fly
% remember the original working directory
pwdir = pwd;

% determine the name and full path of this function
funname = mfilename('fullpath');
mexsrc  = [funname '.c'];
[mexdir, mexname] = fileparts(funname);

try
  % try to compile the mex file on the fly
  warning('trying to compile MEX file from %s', mexsrc);
  cd(mexdir);

  if ispc
    mex -I. -c geometry.c
    mex -I. -c plinproj.c ; mex plinproj.c plinproj.obj geometry.obj
  else
    mex -I. -c geometry.c
    mex -I. -c plinproj.c ; mex -o plinproj plinproj.o geometry.o
  end

  cd(pwdir);
  success = true;

catch
  % compilation failed
  disp(lasterr);
  error('could not locate MEX file for %s', mexname);
  cd(pwdir);
  success = false;
end

if success
  % execute the mex file that was just created
  funname   = mfilename;
  funhandle = str2func(funname);
  [varargout{1:nargout}] = funhandle(varargin{:});
end


text:
classdef ApoapsisAltitudeOptimizationVariable < AbstractOptimizationVariable
    %ApoapsisAltitudeOptimizationVariable Summary of this class goes here
    %   Detailed explanation goes here

    properties
        varObj(1,1) ApoapsisAltitudeTermCondition = ApoapsisAltitudeTermCondition(0);

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        useTf(1,1) = false;
    end

    methods
        function obj = ApoapsisAltitudeOptimizationVariable(varObj)
            obj.varObj = varObj;
            obj.varObj.optVar = obj;

            obj.id = rand();
        end

        function x = getXsForVariable(obj)
            x = [];

            if(obj.useTf)
                x = obj.varObj.apoalt;
            end
        end

        function [lb, ub] = getBndsForVariable(obj)
            lb = obj.lb(obj.useTf);
            ub = obj.ub(obj.useTf);
        end

        function [lb, ub] = getAllBndsForVariable(obj)
            lb = obj.lwrBnd;
            ub = obj.uprBnd;
        end

        function setBndsForVariable(obj, lb, ub)
            obj.lb = lb;
            obj.ub = ub;
        end

        function useTf = getUseTfForVariable(obj)
            useTf = obj.useTf;
        end

        function setUseTfForVariable(obj, useTf)
            obj.useTf = useTf;
        end

        function updateObjWithVarValue(obj, x)
            obj.varObj.apoalt = x;
        end

        function nameStrs = getStrNamesOfVars(obj, evtNum)
            nameStrs = {sprintf('Event %i Apoapsis Altitude Termination Condition', evtNum)};
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%function h = jim(varargin) - alias to "im" for julia compatibility
 function h = jim(varargin)
   out = im(varargin{:});
if nargout
   h = out;
end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [A, vA, vB, bb_rel] = crop_borders(A, bcol, padding)
%CROP_BORDERS Crop the borders of an image or stack of images
%
%   [B, vA, vB, bb_rel] = crop_borders(A, bcol, [padding])
%
%IN:
%   A - HxWxCxN stack of images.
%   bcol - Cx1 background colour vector.
%   padding - scalar indicating how much padding to have in relation to
%             the cropped-image-size (0<=padding<=1). Default: 0
%
%OUT:
%   B - JxKxCxN cropped stack of images.
%   vA     - coordinates in A that contain the cropped image
%   vB     - coordinates in B where the cropped version of A is placed
%   bb_rel - relative bounding box (used for eps-cropping)

% 06/03/15: Improved image cropping thanks to Oscar Hartogensis
% 08/06/15: Fixed issue #76: case of transparent figure bgcolor

    if nargin < 3
        padding = 0;
    end
    [h, w, c, n] = size(A);
    if isempty(bcol)  % case of transparent bgcolor
        bcol = A(ceil(end/2),1,:,1);
    end
    if isscalar(bcol)
        bcol = bcol(ones(c, 1));
    end

    % Crop margin from left
    bail = false;
    for l = 1:w
        for a = 1:c
            if ~all(col(A(:,l,a,:)) == bcol(a))
                bail = true;
                break;
            end
        end
        if bail
            break;
        end
    end

    % Crop margin from right
    bcol = A(ceil(end/2),w,:,1);
    bail = false;
    for r = w:-1:l
        for a = 1:c
            if ~all(col(A(:,r,a,:)) == bcol(a))
                bail = true;
                break;
            end
        end
        if bail
            break;
        end
    end

    % Crop margin from top
    bcol = A(1,ceil(end/2),:,1);
    bail = false;
    for t = 1:h
        for a = 1:c
            if ~all(col(A(t,:,a,:)) == bcol(a))
                bail = true;
                break;
            end
        end
        if bail
            break;
        end
    end

    % Crop margin from bottom
    bcol = A(h,ceil(end/2),:,1);
    bail = false;
    for b = h:-1:t
        for a = 1:c
            if ~all(col(A(b,:,a,:)) == bcol(a))
                bail = true;
                break;
            end
        end
        if bail
            break;
        end
    end

    % Crop the background, leaving one boundary pixel to avoid bleeding on resize
    %v = [max(t-padding, 1) min(b+padding, h) max(l-padding, 1) min(r+padding, w)];
    %A = A(v(1):v(2),v(3):v(4),:,:);
    if padding == 0  % no padding
        padding = 1;
    elseif abs(padding) < 1  % pad value is a relative fraction of image size
        padding = sign(padding)*round(mean([b-t r-l])*abs(padding)); % ADJUST PADDING
    else  % pad value is in units of 1/72" points
        padding = round(padding);  % fix cases of non-integer pad value
    end

    if padding > 0  % extra padding
        % Create an empty image, containing the background color, that has the
        % cropped image size plus the padded border
        B = repmat(bcol,(b-t)+1+padding*2,(r-l)+1+padding*2);
        % vA - coordinates in A that contain the cropped image
        vA = [t b l r];
        % vB - coordinates in B where the cropped version of A will be placed
        vB = [padding+1, (b-t)+1+padding, padding+1, (r-l)+1+padding];
        % Place the original image in the empty image
        B(vB(1):vB(2), vB(3):vB(4), :) = A(vA(1):vA(2), vA(3):vA(4), :);
        A = B;
    else  % extra cropping
        vA = [t-padding b+padding l-padding r+padding];
        A = A(vA(1):vA(2), vA(3):vA(4), :);
        vB = [NaN NaN NaN NaN];
    end

    % For EPS cropping, determine the relative BoundingBox - bb_rel
    bb_rel = [l-1 h-b-1 r+1 h-t+1]./[w h w h];
end

function A = col(A)
    A = A(:);
end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [Y,W,SetupStruc] = Process_FastICA_HO_Sawada(s,Transfer,SetupStruc)
K = SetupStruc.FastICA_HO_Sawada.K;
hop = SetupStruc.FastICA_HO_Sawada.hop;
win = hanning(K,'periodic');
win = win/sqrt(sum(win(1:hop:K).^2));
SetupStruc.FastICA_HO_Sawada.win = win;  % Preserve 'win' in 'SetupStruc'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = size(s,2);
for i = 1:N
    X(:,:,i) = fft(enframe(s(:,i),win,hop)');
end
frame_N = size(X,2);
K_m = K/2+1;
Num = size(Transfer,3);
Y = zeros((frame_N-1)*hop+K,Num);
Y_f = zeros(size(X,1),size(X,2),Num);
Y_P = zeros(frame_N,Num,K_m);
%%%%%%%%%%%%%%%%%%%%%%%%%% Obtain processing matrix 'W'
theta = 10^-6;
W = zeros(Num,N,K_m);
A = zeros(1001,K/2)-1; %%%% Show the decrease of non-linear correlation, ICA max iterations 1000
for i = 2:K_m
    X_f = permute(X(i,:,:),[3 2 1]);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% PCA and ICA processing
    [E,D] = PCA(X_f,1,Num);
    V = sqrt(D)\E';
    X_f = V*X_f;
    %%%%%%%% FastICA based on higher-order statistics
    W_f = eye(Num);
    pObj = inf;
    for i_iteration = 1:200
        Obj = 0;
        for i_n = 1:Num
            W_i = W_f(i_n,:);
            y_ = W_i*X_f;
            y_2 = real(y_).^(2)+imag(y_).^(2);
            %%%%%%%%%%%%%%%%%%%%% There are three different contrast functions.
%             %%%  G1(y) = sqrt(a1+y), g1(y) = 1/(2*sqrt(a1+y))
%             g = 0.5*(0.1+y_2).^(-1/2);
%             g_ = -0.25*(0.1+y_2).^(-3/2);
%             Obj = Obj+sum(sqrt(0.1+y_2))/frame_N;
            %%%  G2(y) = log(a2+y),  g2(y) = 1/(a2+y)
            g = (0.1+y_2).^(-1);
            g_ = -(0.1+y_2).^(-2);
            Obj = Obj+sum(log(0.1+y_2))/frame_N;
%             %%%  G3(y) = 0.5*y^2,    g3(y) = y     (Kurtosis)
%             g = y_2;
%             g_ = ones(1,frame_N);
%             Obj = Obj+0.5*sum(y_2.^2)/frame_N;
            %%%%%%%%%%%%%%%%%%%%%
%             W_i = sum(X_f.*repmat(conj(y_).*g,[Num,1]),2)/frame_N-sum(g+y_2.*g_,2)/frame_N*W_i';
            W_i = (X_f*(conj(y_).*g).'-(sum(g,2)+y_2*g_.')*W_i')/frame_N;
            W_i = W_i/norm(W_i);
            W_f(i_n,:) = W_i';
        end
        W_f = (W_f*W_f')^(-1/2)*W_f;
        dObj = pObj-Obj;
        pObj = Obj;
        A(i_iteration,i-1) = Obj;
        if(abs(dObj)/abs(Obj)<theta)
            break;
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     [Y_,W_ICA,A] = FDICA(X_f,eye(Num),A,i);  %%% 'A', 'i' record the decrease for observation
    Y_ = W_f*X_f;
    W(:,:,i) = W_f*V;
    Y_P(:,:,i) = Y_.';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Process the ambiguity of permutation and amplitude
P = Permu_Sawada(W,Y_P,SetupStruc,'all');  %%%% Options: 'DOA','cor', 'all'
for i = 2:K_m
    W(:,:,i) = P(:,:,i)*W(:,:,i);
    Y_ = permute(Y_P(:,:,i),[2 1 3]);
    Y_ = P(:,:,i)*Y_;
    Y_f(i,:,:) = Y_.';
    if(i~=K_m)
        Y_f(K+2-i,:,:) = Y_';
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Recover signals
if(K/hop==2)
    win = ones(K,1);
end
for i = 1:Num
    Y(:,i) = overlapadd(real(ifft(Y_f(:,:,i)))',win,hop);
end
return;

repo: Voice-Separation-and-Enhancement


text:
function sS=transpose(sS)
% transpose list of slipSystem

sS.b = sS.b.';
sS.n = sS.n.';
sS.CRSS = sS.CRSS.';


text:
function varargout = crossentropy(varargin)
% CROSSENTROPY
%
% y = CROSSENTROPY(x,y)
%
% Computes/declares cross entropy -sum(x.*log(y))
%
% See also ENTROPY, KULLBACKLEIBLER

switch class(varargin{1})

    case {'sdpvar','ndsdpvar'}

        varargin{1} = reshape(varargin{1},[],1);
        varargin{2} = reshape(varargin{2},[],1);

        if length(varargin{1})~=length(varargin{2})
            if length(varargin{1})==1
                varargin{1} = repmat(varargin{1},length(varargin{2}),1);
            elseif  length(varargin{2})==1
                varargin{2} = repmat(varargin{2},length(varargin{1}),1);
            else
                error('Dimension mismatch in crossentropy')
            end
        end

        varargout{1} = yalmip('define','crossentropy_internal',[varargin{1};varargin{2}]);

    otherwise
        error('SDPVAR/CROSSENTROPY called with CHAR argument?');
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [ AF ] = AutoCorrelationF( X,m )
AF=zeros(m,1);

for i=1:length(AF)
X1=X(1:end-i);
X1=X1-mean(X1);

X2=X(i+1:end);
X2=X2-mean(X2);
AF(i)=mean(X1.*X2)/sqrt((var(X1)*var(X2)));
end

end


text:
function [dat] = read_nex_data(filename, hdr, begsample, endsample, chanindx)

% READ_NEX_DATA for Plexon *.nex file
%
% Use as
%   [dat] = read_nex_data(filename, hdr, begsample, endsample, chanindx)
%
% See also READ_NEX_HEADER, READ_NEX_EVENT

% Copyright (C) 2007, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

try,
  % work with the original header, not the FieldTrip one
  hdr = hdr.orig;
catch
  % assume that we got the original header
end

numsmp = cell2mat({hdr.varheader.numsmp});
adindx = find(cell2mat({hdr.varheader.typ})==5);
smpfrq = hdr.varheader(adindx(1)).wfrequency;
sgn    = chanindx;
nsmp   = (endsample-begsample+1);
dat    = zeros(length(sgn), nsmp);

fid = fopen_or_error(filename, 'r', 'ieee-le');
for sgnlop=1:length(sgn)

  if hdr.varheader(sgn(sgnlop)).typ == 0
    % read a spike channel
    status = fseek(fid,hdr.varheader(sgn(sgnlop)).offset,'bof');

    % read the sample indices at which spikes occurred
    tim = fread(fid,hdr.varheader(sgn(sgnlop)).cnt,'int32');
    % downsample from 40kHz to the A/D sampling frequency
    tim = (tim ./ hdr.filheader.frequency) * smpfrq + 1;  % add one sample, since ts=0 corresponds to sample=1
    tim = round(tim);   % needed because of the edges in histc
    % select only samples between the desired begin and end
    % tim = tim(find(tim>=begsample & tim<=endsample));
    % convert sample indices into a continuous signal
    if ~isempty(tim)
      dum = histc(tim-begsample, 0:(nsmp-1), 1);
      dat(sgnlop,:) = dum(:)';
    end

  elseif hdr.varheader(sgn(sgnlop)).typ == 5
    % read an A/D channel
    status = fseek(fid,hdr.varheader(sgn(sgnlop)).offset,'bof');

    % this just reads the times of LFP starts
    tim = fread(fid,hdr.varheader(sgn(sgnlop)).cnt,'int32');
    % this just reads the indices of LFP starts
    ind = fread(fid,hdr.varheader(sgn(sgnlop)).cnt,'int32');
    if length(ind)>1
      ft_error('multiple A/D segments are not supported');
    end

    % convert from timestamps to samples, expressed in the sampling frequency of the AD channels
    tim = (tim ./ hdr.filheader.frequency) * smpfrq;
    tim = round(tim);

    ch_begsample = begsample - tim;
    ch_endsample = endsample - tim;

    if (ch_begsample<1)
      ft_error(sprintf('cannot read before the begin of the recorded data (channel %d)', sgn(sgnlop)));
    elseif (ch_endsample>hdr.varheader(sgn(sgnlop)).numsmp)
      ft_error(sprintf('cannot read beyond the end of the recorded data (channel %d)', sgn(sgnlop)));
    end

    % seek to the beginning of the interesting data, correct for the A/D card initialisation delay
    fseek(fid,(ch_begsample-1)*2, 'cof');
    % read the actual data for the whole channel
    dum = fread(fid,nsmp,'int16');
    % convert to mV
    dat(sgnlop,:) = dum(:)' * hdr.varheader(sgn(sgnlop)).adtomv;

  else
    % ft_warning('unsupported data format for channel %s', hdr.label{sgn(sgnlop)});
  end

end
status = fclose(fid);

author: fieldtrip

repo: fieldtrip


text:
function d = spline_pchip_set ( n, x, f )

%*****************************************************************************80
%
%% SPLINE_PCHIP_SET sets derivatives for a piecewise cubic Hermite interpolant.
%
%  Discussion:
%
%    This routine computes what would normally be called a Hermite
%    interpolant.  However, the user is only required to supply function
%    values, not derivative values as well.  This routine computes
%    "suitable" derivative values, so that the resulting Hermite interpolant
%    has desirable shape and monotonicity properties.
%
%    The interpolant will have an extremum at each point where
%    monotonicity switches direction.
%
%    The resulting piecewise cubic Hermite function may be evaluated
%    by SPLINE_PCHIP_VAL.
%
%    This routine was originally called "PCHIM".
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2005
%
%  Author:
%
%    Original FORTRAN77 version by Fred Fritsch.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Fred Fritsch, Ralph Carlson,
%    Monotone Piecewise Cubic Interpolation,
%    SIAM Journal on Numerical Analysis,
%    Volume 17, Number 2, April 1980, pages 238-246.
%
%    Fred Fritsch, J Butland,
%    A Method for Constructing Local Monotone Piecewise Cubic Interpolants,
%    LLNL Preprint UCRL-87559, April 1982.
%
%  Parameters:
%
%    Input, integer N, the number of data points.  N must be at least 2.
%
%    Input, real X(N), the strictly increasing independent
%    variable values.
%
%    Input, real F(N), dependent variable values to be interpolated.  This
%    routine is designed for monotonic data, but it will work for any F-array.
%    It will force extrema at points where monotonicity switches direction.
%
%    Output, real D(N), the derivative values at the
%    data points.  If the data are monotonic, these values will determine
%    a monotone cubic Hermite function.
%

%
%  Check the arguments.
%
  if ( n < 2 )
    ierr = -1;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SPLINE_PCHIP_SET - Fatal error!\n' );
    fprintf ( 1, '  Number of data points less than 2.\n' );
    error ( 'SPLINE_PCHIP_SET - Fatal error!' );
  end

  for i = 2 : n
    if ( x(i) <= x(i-1) )
      ierr = -3;
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SPLINE_PCHIP_SET - Fatal error!\n' );
      fprintf ( 1, '  X array not strictly increasing.\n' );
      error ( 'SPLINE_PCHIP_SET - Fatal error!' );
    end
  end

  ierr = 0;
  nless1 = n - 1;
  h1 = x(2) - x(1);
  del1 = ( f(2) - f(1) ) / h1;
  dsave = del1;
%
%  Special case N=2, use linear interpolation.
%
  if ( n == 2 )
    d(1) = del1;
    d(n) = del1;
    return
  end
%
%  Normal case, 3 <= N.
%
  h2 = x(3) - x(2);
  del2 = ( f(3) - f(2) ) / h2;
%
%  Set D(1) via non-centered three point formula, adjusted to be
%  shape preserving.
%
  hsum = h1 + h2;
  w1 = ( h1 + hsum ) / hsum;
  w2 = -h1 / hsum;
  d(1) = w1 * del1 + w2 * del2;

  if ( pchst ( d(1), del1 ) <= 0.0 )

    d(1) = 0.0
%
%  Need do this check only if monotonicity switches.
%
  elseif ( pchst ( del1, del2 ) < 0.0 )

     dmax = 3.0 * del1;

     if ( abs ( dmax ) < abs ( d(1) ) )
       d(1) = dmax;
     end

  end
%
%  Loop through interior points.
%
  for i = 2 : nless1

    if ( 2 < i )
      h1 = h2;
      h2 = x(i+1) - x(i);
      hsum = h1 + h2;
      del1 = del2;
      del2 = ( f(i+1) - f(i) ) / h2;
    end
%
%  Set D(I)=0 unless data are strictly monotonic.
%
    d(i) = 0.0;

    temp = pchst ( del1, del2 );

    if ( temp < 0.0 )

      ierr = ierr + 1;
      dsave = del2;
%
%  Count number of changes in direction of monotonicity.
%
    elseif ( temp == 0.0 )

      if ( del2 ~= 0.0D+00 )
        if ( pchst ( dsave, del2 ) < 0.0 )
          ierr = ierr + 1;
        end
        dsave = del2;
      end
%
%  Use Brodlie modification of Butland formula.
%
    else

      hsumt3 = 3.0 * hsum;
      w1 = ( hsum + h1 ) / hsumt3;
      w2 = ( hsum + h2 ) / hsumt3;
      dmax = max ( abs ( del1 ), abs ( del2 ) );
      dmin = min ( abs ( del1 ), abs ( del2 ) );
      drat1 = del1 / dmax;
      drat2 = del2 / dmax;
      d(i) = dmin / ( w1 * drat1 + w2 * drat2 );

    end

  end
%
%  Set D(N) via non-centered three point formula, adjusted to be
%  shape preserving.
%
  w1 = -h2 / hsum;
  w2 = ( h2 + hsum ) / hsum;
  d(n) = w1 * del1 + w2 * del2;

  if ( pchst ( d(n), del2 ) <= 0.0 )
    d(n) = 0.0;
  elseif ( pchst ( del1, del2 ) < 0.0 )
%
%  Need do this check only if monotonicity switches.
%
    dmax = 3.0 * del2;

    if ( abs ( dmax ) < abs ( d(n) ) )
      d(n) = dmax;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
width = 16;
height = 12;
fontsize = 9;
plot_folder = 'projects/seminar/output/';
% filename = ['data/eco/plugs/02/09/2012-07-13.mat'];
day = '2012-07-13';
household = 02;
dataset = 'eco';
plot_title = 'Laptop consumption - 2012-07-13';
plot_filename = 'laptop-2012-07-13';
appliance_id = getApplianceID('Laptop');
granularity = 1;
appliance_consumption = read_plug_data(dataset, ...
                                household, ...
                                appliance_id, ...
                                evaluation_days, ...
                                granularity);

idx_start = 1;
idx_stop = 86400;

% total_consumption = read_smartmeter_data(dataset, ...
%                                     num2str(household, '%02d'), ...
%                                     evaluation_days, ...
%                                     granularity, ...
%                                     'powerallphases');

%% plot
if ~exist(plot_folder, 'dir')
    mkdir(plot_folder);
end
fig = figure;
xvals = datenum(day)+(idx_start:idx_stop)/86400;
plot(xvals, appliance_consumption);
% dateformat: http://ch.mathworks.com/help/matlab/ref/datetick.html#btpnuk4-1
% dateFormat = 15;
dateFormat = 'HH:MM AM';
datetick('x',dateFormat);
% legend(['Smart meter', appliance_names_inferred], 'Location', 'NorthWest');
title(plot_title);
ylabel('Power [W]');
xlabel('Time of day');
fig = make_report_ready(fig, 'size', [width, height], 'fontsize', fontsize);
saveas(fig, [plot_folder, plot_filename, '.png'], 'png');
print('-depsc2', '-cmyk', '-r600', [plot_folder, plot_filename, '.eps']);
close(fig);

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
% Test for trigBary.m.
function pass = test_trigBary(pref)

if ( nargin < 1 )
    pref = chebfunpref();
end

% Generate a few random points in [-1 1] to use as test values.
seedRNG(3453);
xr = 2 * rand(1000, 1) - 1;

% Set an error tolerance.
tol = 1.0e-12;

% Check evaluation for a trigonometric polynomial.
p1 = @(x) cos(4*pi*x) - 2*sin(3*pi*x) + 3*sin(2*pi*x) - 2*cos(pi*x) + 1;

xk = trigpts(10, [-1,1]);
y = trigBary(xr, p1(xk), xk, [-1, 1]);
err = norm(p1(xr) - y, Inf);
pass(1) = err < tol;

xk = trigpts(1001, [-1,1]);
y = trigBary(xr, p1(xk), xk, [-1, 1]);
err = norm(p1(xr) - y, Inf);
pass(2) = err < tol;

% Test interpolation:
xk = trigpts(10, [-1,1]);
y = trigBary(xk, p1(xk), xk, [-1, 1]);
err = norm(p1(xk) - y, Inf);
pass(3) = err < tol;

xk = chebpts(8, [-1, 1], 1);
p2 = @(x) 2*sin(3*pi*x) + 3*sin(2*pi*x) - 2*cos(pi*x) + 1;
y = trigBary(xr, p2(xk), xk, [-1, 1]);
err = norm(p2(xr) - y, Inf);
pass(4) = err < tol;

%%
% Check evaluation for an array of two polynomials.
q = @(x) [.45 + sin(pi*x), .32 + sin(pi*x) + cos(2*pi*x)];
tol = 1e-10;
xk = -2 + 4*rand(8,1);
xr = -2 + 4*rand(100, 1);
difference = q(xr) - trigBary(xr, q(xk), xk, [-2, 2]);
err = norm(difference(:), Inf);
pass(5) = err < tol;

% Test interpolation:
difference = q(xk) - trigBary(xk, q(xk), xk, [-2, 2]);
err = norm(difference(:), Inf);
pass(6) = err < tol;

% Test trigBaryWeights():
n = 2500;
w = trigBaryWeights(sort(rand(n,1)));
pass(7) = length(w) == n;

% Check for #1743.
xk = [-1 (-1/3 - eps) (1/3 - eps) 1].';
fvals = [1 ; 1 ; 1 ; 1 + eps];
dom = [-pi pi];
y = trigBary(2, fvals, xk, dom);
pass(8) = ~isinf(y) && ~isnan(y) && (abs(y - 1) < 10*eps);

% Check for #1744.
w = trigBaryWeights(linspace(-1, 1, 4).');
w_ex = [0.38883657630910357 ; -1.0 ; 1.0 ; -0.38883657630910357];
err = norm(abs(w) - abs(w_ex), Inf);
tol = 10*eps;
pass(9) = err < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function varargout = chebcoeffs3(f)
%CHEBCOEFFS3   Trivariate Chebyshev coefficients
%   C = CHEBCOEFFS3(F) returns the tensor of trivariate coefficients such
%   that
%   C = sum_{i=0}^{m-1}  sum_{j=0}^{n-1} sum_{k=0}^{p-1}
%                                      C(i+1,j+1,k+1) T_i(x) T_j(y) T_k(z).
%
%   [CORE, C, R, T] = CHEBCOEFFS3(F) returns the same coefficients kept in
%   the Tucker form.
%
% See also CHEBCOEFFS, CHEBCOEFFS2 and COEFFS3

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

if ( isempty(f) )
    varargout = {[]};
    return
end

if ( iszero(f) )
    varargout = {0};
    return
end

% Get low rank representation of f:
[core, cols, rows, tubes] = tucker(f);

% Get the coeffs of the columns, rows and tubes:
cols_coeffs = chebcoeffs(cols);
rows_coeffs = chebcoeffs(rows);
tubes_coeffs = chebcoeffs(tubes);

if ( nargout <= 1 )
    % Return the tensor of coefficients:
    varargout = {chebfun3.txm(chebfun3.txm(chebfun3.txm(core, ...
        cols_coeffs, 1), rows_coeffs, 2), tubes_coeffs, 3)};
elseif ( nargout == 4 )
    varargout = {core, cols_coeffs, rows_coeffs, tubes_coeffs};
else
    % Two or three outputs are not allowed.
    error('CHEBFUN:CHEBFUN3:chebcoeffs3:outputs', ...
        'Incorrect number of outputs.');
end

end


text:
function [pass, u1, u2, info1, info2] = test_scalarODE_breakpoints(pref)
% A nonlinear CHEBOP test. This test tests a scalar ODE, where there are two
% breakpoints in the domain of the CHEBOP.  It solves the problem using chebcolloc1,
% chebcolloc2 and ultraS discretizations.
%
% Asgeir Birkisson, May 2014.

%% Setup
if ( nargin == 0 )
    pref = cheboppref;
end

dom = [0 1 2 pi];
pref.damping = 1;

N = chebop(@(x,u) diff(u,2) + sin(u-.2), dom);
N.lbc = @(u) u - 2;
N.rbc = @(u) u - 3;
rhs = 0;

%% Try different discretizations
% Start with chebcolloc2
pref.discretization = @chebcolloc2;
[u1, info1] = solvebvp(N, rhs, pref);

%% Change to ultraS
pref.discretization = @ultraS;
[u2, info2] = solvebvp(N, rhs, pref);

%% Change to chebcolloc1
pref.discretization = @chebcolloc1;
[u3, info3] = solvebvp(N, rhs, pref);

%% Did we pass?
% To pass, both residuals have to be small, but we should not expect u1 and u2
% to be identical!
tol = 5e4*pref.bvpTol;
err1 = normest(N(u1));
err2 = normest(N(u2));
err3 = normest(N(u3));
% TODO: This used to be 1*, 10* and 1*. Should try to restore once we tune the
% algorithms better.
pass(1) = err1 < tol;
pass(2) = err2 < tol;
pass(3) = err3 < tol;
pass(4) = ( (norm(u1 - u2) ~= 0) && (norm(u2 - u3) ~= 0) && ...
    (norm(u1 - u3) ~= 0));

end


text:
function [h,ax] = smooth(v,varargin)
% low level function for plotting functions on the sphere
%
% Syntax
%   smooth(v,values)
%
% Input
%
% Output
%
% Options
%
% See also
%

h = [];

% initialize spherical plot
opt = delete_option(varargin,{'lineStyle','lineColor','lineWidth','color'},1);
sP = newSphericalPlot(v,opt{:},'doNotDraw');

for j = 1:numel(sP)

  % ---------------- extract colors --------------------------

  % color given by the first argument?
  if ~isempty(varargin) && ~isempty(varargin{1}) ...
      && (islogical(varargin{1}) || isnumeric(varargin{1}))

    cdata = double(varargin{1});
    S2G = v;

  else % no color given -> do kernel density estimation

    sR = sP(j).sphericalRegion;
    if isfield(v.opt,'region'), sR = [sR,v.opt.region]; end %#ok<AGROW>
    S2G = plotS2Grid(sR);

    cdata = calcDensity(v(:),S2G,'halfwidth',5*degree,varargin{:});
    cdata = reshape(cdata,size(S2G));

  end

  sP(j).updateMinMax(cdata);

  % interpolate if no regular grid was given
  if ~isOption(S2G,'plot') || ~S2G.opt.plot

    if size(S2G,1) == 1 || size(S2G,2) == 1

      S2G = plotS2Grid(sP(j).sphericalRegion,varargin{:});
      cdata = interp(v,cdata,S2G,'cutOutside',varargin{:});

    elseif ~isa(sP(j).proj,'plainProjection')

      % close the gap between 0 and 2*pi
      varargin = set_option(varargin,'correctContour');

    end
  end

  % scale the data
  [cdata,colorRange] = scaleData(cdata,varargin{:});
  %if ~any(isnan(colorRange)), caxis(sP(j).ax,colorRange);end

  % ------------- compute contour lines ------------------------

  % number of contour lines
  contours = get_option(varargin,'contours',50);
  contours = get_option(varargin,{'contourf','contour'},contours,'double');

  % specify contourlines explicitely
  if length(contours) == 1
    if check_option(varargin,'log')
      contours = logspace(log10(colorRange(1)),log10(colorRange(2)),contours);
    else
      contours = linspace(colorRange(1),colorRange(2),contours);
    end
  end

  % ----------------- draw contours ------------------------------

  washold = ishold(sP(j).ax);
  hold(sP(j).ax,'on')

  % project data
  [x,y] = project(sP(j).proj,S2G,'noAntipodal');

  % extract non nan data
  data = reshape(cdata,size(x));

  % plot contours
  h = [h,betterContourf(sP(j).hgt,x,y,data,contours,varargin{:})]; %#ok<AGROW>

  if ~washold, hold(sP(j).ax,'off'); end

  % --------------- finalize the plot ---------------------------

  % adjust caxis according to colorRange
  if ~any(isnan(colorRange)), caxis(sP(j).ax,colorRange); end
  if check_option(varargin,'log'), set(sP(j).ax,'colorScale','log'); end

  % colormap
  if ~strcmpi(get_option(varargin,'fill'),'off')
    mtexColorMap(sP(j).ax,getMTEXpref('defaultColorMap'));
  end

  % bring grid in front
  sP(j).doGridInFront;

  sP(j).plotAnnotate(varargin{:})

end

% set styles
varargin = delete_option(varargin,'parent',1);
optiondraw(h,varargin{:});

if isappdata(sP(1).parent,'mtexFig')
  mtexFig = getappdata(sP(1).parent,'mtexFig');
  mtexFig.drawNow('figSize',getMTEXpref('figSize'),varargin{:});
end

if nargout == 0
  clear h;
else
  ax = [sP.ax];
end

end

% ------------------------------------------------------------
function h = betterContourf(ax,X,Y,data,contours,varargin)

if numel(unique(data)) == 1, data(1) = data(1) + 2*eps; end

% workauround for a MATLAB Bug
%if mean(X(:,1)) > mean(X(:,end))
%  X = fliplr(X);
%  Y = fliplr(Y);
%  data = flipdim(data,2);
%end

% contour correction
if check_option(varargin,'correctContour')
  X = [X;X(1,:)];
  Y = [Y;Y(1,:)];
  data = [data;data(1,:)];
end

if check_option(varargin,'pcolor')

  opt = {};
  if numel(data) >= 500
    opt = {'LineStyle','none','FaceColor','interp'};
  end

  % pcolor is actually surface
  h = surface(X,Y,zeros(size(X)),data,opt{:},'parent',ax);

else

  % extract style
  opt = extract_argoption([{'LineStyle','none','Fill','on'},varargin],{'LineStyle','Fill','LineColor'});

  if ~check_option(varargin,'contour')
    [CM,h] = contourf(X,Y,data,contours,opt{:},'parent',ax); %#ok<ASGLU>
  else
    [CM,h] = contour(X,Y,data,contours,opt{:},'parent',ax); %#ok<ASGLU>
  end

end

% do not display in the legend
set(get(get(h,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');

end


text:
function sphere_stereograph_test02 ( )

%*****************************************************************************80
%
%% SPHERE_STEREOGRAPH_TEST02 checks the generalized mapping.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 November 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_STEREOGRAPH_TEST02\n' );
  fprintf ( 1, '  SPHERE_STEREOGRAPH standard mapping from sphere to plane.\n' );
  fprintf ( 1, '  SPHERE_STEREOGRAPH2 generalized mapping:\n' );
  fprintf ( 1, '  (focus and center are arbitrary)\n' );
  fprintf ( 1, '  Check that these two functions can agree.\n' );

  m = 3;
  n = 100;

  focus = [ 0, 0, -1 ]';
  center = [ 0, 0, 0 ]';
%
%  Check #1.
%
  seed = 123456789;
  p1 = uniform_on_sphere01_map ( m, n, seed );

  q1 = sphere_stereograph ( p1 );

  q2 = sphere_stereograph2 ( p1, focus, center );

  dif = sqrt ( sum ( ( q1 - q2 ).^2, 1 ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Map points from sphere to plane.\n' );
  fprintf ( 1, '  Maximum difference for %d points was %f\n', n, max ( dif ) );
%
%  Check #2.
%
  q1 = zeros(m,n);
  q1(1:m-1,1:n) = rand ( m - 1, n );
  q1(m,1:n) = 1.0;

  p1 = sphere_stereograph_inverse ( q1 );

  p2 = sphere_stereograph2_inverse ( q1, focus, center );

  dif = sqrt ( sum ( ( p1 - p2 ).^2, 1 ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Map points from plane to sphere.\n' );
  fprintf ( 1, '  Maximum difference for %d points was %f\n', n, max ( dif ) );

  return
end