text:
function handles = vis_plotModelValidation(whitestats,PCstats,stabilitystats,varargin)
% Visualize the results of model validation as computed by
% pop_est_validateMVAR()
%
% Inputs:
%
%       whitestats:     Cell array of structs each containing results of
%                       tests for whiteness of residuals for a single
%                       dataset/condition
%                       (as output by est_checkMVARWhiteness())
%                       If not available, pass in an empty matrix to skip
%                       plotting this result.
%       PCstats:        Cell array of structs each containing results of
%                       tests for consistency of the model for a single
%                       dataset/condition
%                       (as output by est_checkMVARConsistency())
%                       If not available, pass in an empty matrix to skip
%                       plotting this result.
%       stabilitystats: Cell array of structs each containing results of
%                       tests for stability of the model for a single
%                       dataset/condition
%                       (as output by est_checkMVARStability())
%                       If not available, pass in an empty matrix to skip
%                       plotting this result.
%
%
% Optional:
%
%       whitenessCriteria:   Cell array of whiteness criteria
%       checkWhiteness:      (boolean) plot whiteness results
%       checkConsistency:    (boolean) plot consistency results
%       checkStability:      (boolean) plot stability results
%       windowTimes:        Times (sec) of windows for x-axis labeling.
%                           This MUST be the same length as number of
%                           time windows (e.g.) whitestats{1}.winStartIdx
%       conditions:         Cell array of condition labels
%                           e.g. conditions = {ALLEEG.condition}
%
% Outputs:
%
%       Figures will be rendered showing model validation results for each
%       dataset. The handles to the figures are returned.
%
% See Also: pop_est_validateMVAR(),est_checkMVARWhiteness(),
%           est_checkMVARStability(), est_checkMVARConsistency()
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual. Chapter 6.
%   Available at: http://www.sccn.ucsd.edu/wiki/Sift
%
% Author: Tim Mullen, 2010, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% if ~isempty(varargin) && isstruct(varargin{1})
%     varargin = [hlp_struct2varargin(varargin{1}) varargin];
% end

whitenessCriteria = {};

% find the names of all whiteness tests stored in whitestats
if ~isempty(whitestats{1})
    fn = fieldnames(whitestats{1});
    for i=1:length(fn)
        if isfield(whitestats{1}.(fn{i}),'fullname')
            whitenessCriteria = [whitenessCriteria fn{i}];
        end
    end
end

num_conds = max([length(whitestats),length(PCstats),length(stabilitystats)]);

g = finputcheck(varargin, ...
    {'whitenessCriteria'   'cell'  whitenessCriteria   whitenessCriteria; ...
    'checkWhiteness',     'boolean'   []          true; ...
    'checkConsistency'    'boolean'   []          true; ...
    'checkStability'      'boolean'   []          true; ...
    'windowTimes'         'real'      []          [];   ...
    'conditions'          'cell'      {}          {}; ...
    },'mode','ignore','quiet');

if isempty(whitestats{1}),         g.checkWhiteness    = false;    end
if isempty(PCstats{1}),            g.checkConsistency  = false;    end
if isempty(stabilitystats{1}),     g.checkStability    = false;    end
if isempty(g.conditions),          g.conditions        = cell(1,num_conds); end

numrows = sum([g.checkWhiteness g.checkConsistency g.checkStability]);
numcols = 1;

for cond = 1:num_conds

    if isempty(g.conditions{cond})
        g.conditions{cond} = sprintf('Condition %d',cond);
    end

    % plot results
    handles(cond) = figure('Name',sprintf('%s - Model Validation Results',g.conditions{cond}));
    curplot=1;

    if g.checkWhiteness
        % Plot results of residual whiteness checks
        if ~iscell(whitestats), whitestats = {whitestats}; end

        ax=subplot(numrows,numcols,curplot);
        for i = 1:length(g.whitenessCriteria)
            wcstr = lower(hlp_variableize(g.whitenessCriteria{i}));
            wc = whitestats{cond}.(wcstr);
            pvals(i,:) = wc.pval;
            if size(pvals,2)>1
                lgnd{i} = sprintf('%s (%d/%d white)',wc.fullname, sum(wc.w),length(wc.w));
            else
                lgnd{i} = sprintf('%s (%swhite)',wc.fullname, fastif(wc.w,'','not '));
            end
        end

        if size(pvals,2)>1
            % more than one window -- make lineplot
            if ~isempty(g.windowTimes)
                abscissa = g.windowTimes;
            else
                abscissa = 1:length(whitestats{cond}.winStartIdx);
            end
            plot(ax,abscissa,pvals','Marker','.');
            xlabel(ax,fastif(isempty(g.windowTimes),'Window number','Time (sec)'));

            legend(ax,lgnd);
            set(ax, ...
                'Xlim',[abscissa(1)-abs(diff(abscissa(1:2))) abscissa(end)+abs(diff(abscissa(1:2)))], ...
                'Ylim',[max(0,min(pvals(:))-0.5), min(1,max(pvals(:))+0.5)]);
            axcpystr = 'legend(''''show'''');';
        else
            % single window -- create bar plots
            h = bar(ax,pvals);
            ch = get(h,'Children');

            set(ax,'xticklabel',lgnd);
            colors = [[1 0 0];[0 0 1];[0 1 0];[0 0 0];[1 0 1];[0 1 1]];
            colors = colors(1:length(pvals),:);
            set(ch,'FaceVertexCData',colors);
            %             set(gca,'xlim',[0 length(length(g.whitenessCriteria))+1],'Ylim',[max(0,min(pvals(:))-0.5), min(1,max(pvals(:))+0.5)]);
            axcpystr = '';
        end

        hl=hline(whitestats{cond}.alpha,'b','1-p_{port}',[1.01 -0.01],ax);
        set(hl,'linestyle','--','linewidth',2);
        ylabel({'Whiteness Significance'});

        curplot=curplot+1;

        if ismember_bc('acf',lower(g.whitenessCriteria))
            hl=hline(1-whitestats{cond}.alpha,'r','1-P_{acf}',[1.01 -0.01],ax);
            set(hl,'linestyle','--','linewidth',2);
        end

        axcopy(ax,axcpystr);
    end

    if g.checkConsistency
        % Plot results of consistency checks

        if ~iscell(PCstats), PCstats = {PCstats}; end

        ax=subplot(numrows,numcols,curplot);
        if length(PCstats{cond}.PC)>1
            % more than one window -- make lineplot
            meanpc = mean(PCstats{cond}.PC);
            stdpc  = std(PCstats{cond}.PC);
            if ~isempty(g.windowTimes)
                abscissa = g.windowTimes;
            else
                abscissa = 1:1:length(PCstats{cond}.winStartIdx);
            end
            plot(ax,abscissa,PCstats{cond}.PC,'Marker','.');
            xlabel(ax,fastif(isempty(g.windowTimes),'Window number','Time (sec)'));
            set(ax,...
            'Xlim',[abscissa(1)-abs(diff(abscissa(1:2))) abscissa(end)+abs(diff(abscissa(1:2)))], ...
            'Ylim',[min(PCstats{cond}.PC)-50 min(max(PCstats{cond}.PC)+50,100)]);
            % make legend
            text(0.98,0.9,sprintf('Mean PC: %0.2f+/-%0.3g%%',meanpc,stdpc), ...
                    'units','normalized','horizontalalignment','right', ...
                    'edgecolor','k','backgroundcolor','w','parent',ax);
            % make a small histogram on right side of plot
            %             axpos = get(ax,'Position');
            %             axhist = axesRelative(ax, 'Position',[1.01 0 0.1 1], 'Units','Normalized');   %axes('Position',[axpos(1)+axpos(3)+0.01 axpos(2) 0.05 axpos(4)]);
            %             hist(axhist,PCstats{cond}.PC,10);
            %             vline(mean(PCstats{cond}.PC),':r');
            %             set(axhist,'View',[90 90]);
            %             set(axhist,'xdir','rev');
        else
            % single window -- make barplot
            bar(ax,PCstats{cond}.PC);
            % make legend
            text(0.98,0.9,sprintf('(%0.2f%% Consistent)',PCstats{cond}.PC), ...
                'units','normalized','horizontalalignment','right', ...
                'edgecolor','k','backgroundcolor','w','parent',ax);
%             legend(sprintf('(%0.2f%% Consistent)',PCstats{cond}.PC));
            set(ax,'Ylim',[min(0,PCstats{cond}.PC)-1 max(PCstats{cond}.PC,100)+1]);
            set(ax,'xtick',[]);
        end

        ylabel('Percent Consistency');
        axcopy(ax);
        curplot = curplot+1;
    end

    if g.checkStability
        % Plot results of stability checks

        if ~iscell(stabilitystats), stabilitystats = {stabilitystats}; end

        % plot stability results
        ax=subplot(numrows,numcols,curplot);
        if length(stabilitystats{cond}.stability)>1
            % more than one window -- make lineplot
            %boxplot(real(lambda)');
            maxlambda = max(real(stabilitystats{cond}.lambda),[],2);
            if ~isempty(g.windowTimes)
                abscissa = g.windowTimes;
            else
                abscissa = 1:length(stabilitystats{cond}.winStartIdx);
            end
            plot(ax,abscissa,maxlambda,'Marker','.');
            xlabel(ax,fastif(isempty(g.windowTimes),'Window number','Time (sec)'));

            set(ax,'Xlim',[abscissa(1)-abs(diff(abscissa(1:2))) abscissa(end)+abs(diff(abscissa(1:2)))], ...
                   'Ylim',[1.2*min(maxlambda(:)) max(0.01,1.3*max(maxlambda(:)))]);
        else
            % single window -- make barplot
            maxlambda = max(real(stabilitystats{cond}.lambda(:)));
            bar(ax,maxlambda);
            set(ax,'xtick',[]);
        end

        %     set(gca,'ylim',[max(0,0.7*min(abs(lambda(:)))) max(1.3,1.3*max(abs(lambda(:))))]);
        %     axis auto
        hl=hline(0,[],[],[],ax);
        set(hl,'linestyle','--','linewidth',2);
        ylabel(ax,{'Stability Index','(should be < 0)'});
        numstable = sum(stabilitystats{cond}.stability);
        % make legend
        text(0.98,0.9,sprintf('(%d/%d stable)',numstable,length(stabilitystats{cond}.stability)), ...
                'units','normalized','horizontalalignment','right', ...
                'edgecolor','k','backgroundcolor','w','parent',ax);
%         legend(ax,sprintf('(%d/%d stable)',numstable,length(stabilitystats{cond}.stability)));

        axcopy(ax);
    end
end

try
    icadefs; set(handles, 'color', BACKCOLOR);
catch
end;

repo: REST


text:
function [ vX ] = ProjectL2Ball( vY, ballRadius )
% ----------------------------------------------------------------------------------------------- %
% [ vX ] = ProjectL2Ball( vY, ballRadius, stopThr )
%   Solving the Orthogonal Projection Problem of the input vector onto the
%   L2 Ball.
% Input:
%   - vY            -   Input Vector.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - ballRadius    -   Ball Radius.
%                       Sets the Radius of the L2 Ball. For Unit L2 Ball
%                       set to 1.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: (0, inf).
% Output:
%   - vX            -   Output Vector.
%                       The projection of the Input Vector onto the L2
%                       Ball.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% References
%   1.  h
% Remarks:
%   1.  a
% TODO:
%   1.  U.
% Release Notes:
%   -   1.0.000     29/06/2017  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

vX = min((ballRadius / norm(vY, 2)), 1) * vY;

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function impred =  runPatchVDSR(net, imlow, gpu, rf)
v = ceil(size(imlow, 1)/2);
h = ceil(size(imlow, 2)/2);

imTL = imlow(1:v+rf,   1:h+rf);
imBL = imlow(v-rf+1:end, 1:h+rf);
imTR = imlow(1:v+rf,   h-rf+1:end);
imBR = imlow(v-rf+1:end, h-rf+1:end);

if gpu, imTL = gpuArray(imTL); end;
impredTL = runVDSR(net, imTL, gpu);
impredTL = impredTL(1:v, 1:h);

if gpu, imBL = gpuArray(imBL); end;
impredBL = runVDSR(net, imBL, gpu);
impredBL = impredBL(rf+1:end, 1:h);

if gpu, imTR = gpuArray(imTR); end;
impredTR = runVDSR(net, imTR, gpu);
impredTR = impredTR(1:v, rf+1:end);

if gpu, imBR = gpuArray(imBR); end;
impredBR = runVDSR(net, imBR, gpu);
impredBR = impredBR(rf+1:end, rf+1:end);

impredL = cat(1, impredTL, impredBL);
impredR = cat(1, impredTR, impredBR);
impred = cat(2, impredL, impredR);

author: thomas-koehler

repo: SupER


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function outNames = SpharmMatAlignment(confs, objs, method)

numSbj = length(objs);
if ~exist(confs.OutDirectory,'dir')
    mkdir(confs.OutDirectory);
end
if (~exist([confs.OutDirectory '/Logs'],'dir'))
    mkdir([confs.OutDirectory '/Logs']);
end
outNames = {};

%class(confs) - struct
%confs.vars -list of variable
%class(objs) - cell

h = waitbar(0,'Please wait...');
for i = 1:numSbj
    file = objs{i};
    [path, name, ext] = fileparts(file);

    switch method
        case 'AligSHREC'
            diary(fullfile([confs.OutDirectory '/Logs'],[date, '_', num2str(round(cputime)), '_' name '_AlignSHREC.log']));
            % Call AligSHREC method
            if isempty(confs.Template)
                disp('SHREC needs a template object');
                return;
            end
            [vertices, sph_verts, faces, fvec, outNames{end+1}] = align_CPS_SHREC(file, confs);
            clear('vertices', 'sph_verts', 'faces', 'fvec');
        case 'AligFOE'
            diary(fullfile([confs.OutDirectory '/Logs'],[date, '_', num2str(round(cputime)), '_' name '_AlignFOE.log']));
            % Call AligFOE method
            [vertices, sph_verts, faces, fvec, outNames{end+1}] = align_FOE(file, confs);
            clear('vertices', 'sph_verts', 'faces', 'fvec');
        case 'AligLandmark'
            diary(fullfile([confs.OutDirectory '/Logs'],[date, '_', num2str(round(cputime)), '_' name '_AlignLDMK.log']));
            % Call AligLandmark method
            disp('Not Implmented');
    end
    diary('off');
    waitbar(i/numSbj)
end
close(h);

return

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function printv(v)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2020-2025, by Kai Chen, All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=size(v,1);
fprintf('v     =\n');
for i=1:n
    fprintf('%11.6f',v(i));
    fprintf('\n');
end
fprintf('nv=%d\n',n);

return

repo: GINav


text:
function boxes = masksToBoxes( segs, props )
    r = cell(length(segs));
    for i = 1:length(segs)
        r{i} = maskToBox( segs{i}, props{i} );
    end
    boxes = cat(1,r);
end

author: Cloud-CV

repo: object-proposals


text:
% SIZE - size of memory mapped underlying array
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, Nov. 2008

% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [s varargout] = size(obj,dim)

    if obj.transposed
        if length(obj.dimensions) ~= 2 && length(obj.dimensions) ~= 3
            error('Transposed array must be 2 or 3 dims');
        end
        if length(obj.dimensions) == 2 tmpdimensions = [obj.dimensions(2) obj.dimensions(1)];
        else                           tmpdimensions = [obj.dimensions(3) obj.dimensions(1:end-1)];
        end
    else
        tmpdimensions = obj.dimensions;
    end

    s = tmpdimensions;

    if nargin > 1
        if dim >length(s)
            s = 1;
        else
            s = s(dim);
        end
    else
        if nargout > 1
            s = [s 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
            alls = s;
            s = s(1);
            for index = 1:max(nargout,1)-1
                varargout{index} = alls(index+1);
            end
        end
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [pdf]=gsp_erdos_renyi_density_normalized(y,N,p)
%GSP_ERDOS_RENYI_DENSITY_NORMALIZED  The asymptotic density function of the normalized graph Laplacian
%eigenalues of an Erdos-Renyi random graph
%   Usage:  pdf=gsp_erdos_renyi_density_normalized(y,N,p);
%
%   Input parameters:
%         y               : A vector of values at which the density will be evaluated.
%         N               : Number of vertices in the graph.
%         p               : Edge probability.
%   Output parameters:
%         pdf             : The values of the density function at the points y.
%
%   'gsp_erdos_renyi_density_normalized(y,N,p)' computes the asymptotic density
%   function (for large $N$) of the normalized graph Laplacian
%   eigenvalues of an Erdos-Renyi random graph with $N$ vertices and edge
%   probability $p$. This is a semicircular density function.
%

%   AUTHOR : David I Shuman.
%   TESTING:
%   REFERENCE:

pdf=sqrt((p*N)/(1-p))*(1/(2*pi))*real((y>1-2*sqrt((1-p)/(p*N))).*(y<1+2*sqrt((1-p)/(p*N))).*sqrt(4-p*N/(1-p)*(1-y).^2));

repo: gspbox


text:
%SOLVELP  Solve given (non-integer) linear programming problem using the Simplex Algorithm
%
%     [z, res] = cv.solveLP(Func, Constr)
%     [...] = cv.solveLP(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __Func__ This row-vector corresponds to `c` in the LP problem formulation
%   (see below). It should contain 32- or 64-bit floating-point numbers. As a
%   convenience, column-vector may be also submitted, in the latter case it is
%   understood to correspond to `c'`.
% * __Constr__ m-by-n+1 matrix, whose rightmost column corresponds to `b` in
%   formulation above and the remaining to `A`. It should contain 32 or 64-bit
%   floating-point numbers.
%
% ## Output
% * __z__ The solution will be returned here as a column-vector, it
%   corresponds to `x` in the formulation above. It will contain 64-bit
%   floating-point numbers.
% * __res__ Return code. One of:
%   * __Unbounded__ problem is unbounded (target function can achieve
%     arbitrary high values)
%   * __Unfeasible__ problem is unfeasible (there are no points that satisfy
%     all the constraints imposed)
%   * __Single__ there is only one maximum for target function
%   * __Multi__ there are multiple maxima for target function, the arbitrary
%     one is returned
%
% What we mean here by "linear programming problem" (or LP problem, for short)
% can be formulated as:
%
%     Maximize c*x
%       subject to  A*x <= b
%                   x >= 0
%
% Where `c` is fixed 1-by-n row-vector, `A` is fixed m-by-n matrix, `b` is
% fixed m-by-1 column vector and `x` is an arbitrary n-by-1 column vector,
% which satisfies the constraints.
%
% Simplex algorithm is one of many algorithms that are designed to handle this
% sort of problems efficiently. Although it is not optimal in theoretical
% sense (there exist algorithms that can solve any problem written as above in
% polynomial type, while simplex method degenerates to exponential time for
% some special cases), it is well-studied, easy to implement and is shown to
% work well for real-life purposes.
%
% The particular implementation is taken almost verbatim from:
%
% > Introduction to Algorithms, 3rd edition
% > by T. H. Cormen, C. E. Leiserson, R. L. Rivest and Clifford Stein.
%
% In particular, the [Bland's rule](https://en.wikipedia.org/wiki/Bland%27s_rule)
% is used to prevent cycling.
%
% See also: linprog
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [inH2O] = atm2inH2O(atm)
% Convert pressure from atmospheres to inches of water column.
% Chad A. Greene of chadagreene.com fame
inH2O = atm*406.782;

author: Sable

repo: mcbench-benchmarks


text:
%NU_SVR Support Vector Classifier: NU algorithm
%
% 	[W,J,C] = NU_SVR(A,TYPE,PAR,C,SVR_TYPE,NU_EPS,MC,PD)
%
% INPUT
%   A	    Dataset
%   TYPE  Type of the kernel (optional; default: 'p')
%   PAR   Kernel parameter (optional; default: 1)
%   C     Regularization parameter (0 < C < 1): expected fraction of SV
%         (optional; default: 0.25)
%   SVR_TYPE  This type can be 'nu' or 'epsilon'
%   NU_EPS The corresponding value for NU or epsilon
%   MC    Do or do not data mean-centering (optional; default: 1 (to do))
%   PD    Do or do not the check of the positive definiteness (optional;
%         default: 1 (to do))
%
% OUTPUT
%   W     Mapping: Support Vector Classifier
%   J     Object identifiers of support objects
%   C     Equivalent C regularization parameter of SVM-C algorithm
%
% DESCRIPTION
% Optimizes a support vector classifier for the dataset A by
% quadratic programming. The classifier can be of one of the types
% as defined by PROXM. Default is linear (TYPE = 'p', PAR = 1). In J
% the identifiers of the support objects in A are returned.
%
% C belogs to the interval (0,1). C close to 1 allows for more class
% overlap.  Default C = 0.25.
%
% C is bounded from above by NU_MAX = (1 - ABS(Lp-Lm)/(Lp+Lm)), where
% Lp (Lm) is the number of positive (negative) samples. If NU > NU_MAX
% is supplied to the routine it will be changed to the NU_MAX.
%
% If C is less than some NU_MIN which depends on the overlap between
% classes algorithm will typically take long time to converge (if at
% all).  So, it is advisable to set NU larger than expected overlap.
%
% Output is rescaled in a such manner as if it were returned by SVC with
% the parameter C.
%
%
% SEE ALSO (PRTools Guide)
% NU_SVRO, SVO, SVC, MAPPINGS, DATASETS, PROXM

% Copyright: S.Verzakov, s.verzakov@ewi.tudelft.nl
% Based on SVC.M by D.M.J. Tax, D. de Ridder, R.P.W. Duin
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: nu_svr.m,v 1.2 2009/01/31 18:43:11 duin Exp $

function [W, J, epsilon_or_nu] = nu_svcr(a,type,par,C,svr_type,nu_or_epsilon,mc,pd)
if nargin < 2 | ~isa(type,'prmapping')
	if nargin < 8
		pd = 1;
	end
	if nargin < 7
		mc = 1;
	end
	if nargin < 6
		nu_or_epsilon = [];
	end
	if nargin < 5 | isempty(svr_type)
		svr_type = 'epsilon';
	end

	switch svr_type
	case 'nu'
		if isempty(nu_or_epsilon)
			prwarning(3,'nu is not specified, assuming 0.25.');
			nu_or_epsilon = 0.25;
		end
	%nu = nu_or_epsilon;
	case {'eps', 'epsilon'}
		svr_type = 'epsilon';
		if isempty(nu_or_epsilon)
			prwarning(3,'epsilon is not specified, assuming 1e-2.');
			nu_or_epsilon = 1e-2;
		end
	%epsilon = nu_or_epsilon;
	end

	if nargin < 4 | isempty(C)
		prwarning(3,'C set to 1\n');
		C = 1;
	end

	if nargin < 3 | isempty(par)
		par = 1;
		prwarning(3,'Kernel parameter par set to 1\n');
	end
	if nargin < 2 | isempty(type)
		type = 'p';
		prwarning(3,'Polynomial kernel type is used\n');
	end
	if nargin < 1 | isempty(a)
		W = prmapping(mfilename,{type,par,C,svr_type,nu_or_epsilon,mc,pd});
		W = setname(W,['Support Vector Regression (' svr_type ' algorithm)']);
		return;
	end

	islabtype(a,'targets');
	[m,k] = getsize(a);
	y = gettargets(a);
	% The 1-dim SVR
	if size(y,2) == 1   % 1-dim regression
		uy = mean(y);
		y  = y - uy;
		if mc
			u  = mean(a);
			a  = a - ones(m,1)*u;
		else
			u  = [];
		end

		K = a*proxm(a,type,par);
		% Perform the optimization:
		[v,J,epsilon_or_nu] = nu_svro(+K,y,C,svr_type,nu_or_epsilon,pd);
		% Store the results:
		v(end) = v(end)+uy;
		W = prmapping(mfilename,'trained',{u,a(J,:),v,type,par},getlablist(a),k,1);
		W = setname(W,['Support Vector Regression (' svr_type ' algorithm)']);
		%W = setcost(W,a);
		J = getident(a,J);
		%J = a.ident(J);

	else
		error('multivariate SVR is not supported');
	end

else % execution
	w = +type;
	m = size(a,1);

	% The first parameter w{1} stores the mean of the dataset. When it
	% is supplied, remove it from the dataset to improve the numerical
	% precision. Then compute the kernel matrix using proxm.

	if isempty(w{1})
		d = a*proxm(w{2},w{4},w{5});
	else
		d = (a-ones(m,1)*w{1})*proxm(w{2},w{4},w{5});
	end

	% When Data is mapped by the kernel, now we just have a linear
	% regression  w*x+b:
	d = [d ones(m,1)] * w{3};
	W = setdat(a,d,type);
end

return;

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [hyperplane_normal, hyperplane_offset] = linortfitn(data)
% LINORTFITN  Fit a line to data by ORTHOGONAL least-squares.
%    [N,C] = LINORTFITN(DATA) finds the coefficients of a hyperplane (in
%    Hessian normal form) that best fits the data in an ORTHOGONAL
%    least-squares sense.  Consider the hyperplane
%       H = {x | dot(N,x) + C == 0},
%    and the minimum (Euclidean) distance between this hyperplane and each
%    datapoint DATA(i,:) -- LINORTFITN finds N and C such that the sum of
%    squared distances is minimized.

[M,N] = size(data);
if M <= N,
    error('linortfitn:DegenerateProblem',...
        'There are fewer datapoints than dimensions: the data is perfectly fit by a hyperplane.');
end

[U,S,V] = svd(data - repmat(mean(data),M,1), 0);
hyperplane_normal = V(:,end);
hyperplane_offset = - mean(data * hyperplane_normal);


text:
% This is the main script runner for training, it collects the
% training samples, followed by model training

clear

% define the root name of database
root = '../data_preparation/prepared_data/';

% which scales we're doing
scales = [0.25, 0.35, 0.5];

% the data generation parameters
sigma = 1;
num_samples = 5e5;
ratio_neg = 5;
norm = 1;
normalisation_size = 19;

patch_types = {'reg', 'grad'};

% Should a frontal view be created
frontalView = [1];

% The other views to be used
profileViewInds = [2];
upDownViewInds = [];

version = 'wild';

% the naming of
wild_loc = 'wild_';

for s=scales
    Train_all(root, frontalView, profileViewInds, upDownViewInds, s,...
        sigma, version, norm, 'ratio_neg', ratio_neg,...
        'num_samples', num_samples, 'data_loc', wild_loc,...
        'patch_types', patch_types, 'normalisation_size', normalisation_size);
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function [mn cr]=finminexceptTabu(tmp,Tabu)
%% Written by Muhammet Balcilar, France
% all rights reverved
nodes=tmp(:,1:end-1);
cost=tmp(:,end);
for i=1:size(Tabu,1)
    tt=sum(abs(nodes-Tabu(i,:)),2);
    I=find(tt==0);
    if length(I)>0
        cost(I(1))=inf;
    end
end
[mn cr]=min(cost);

author: balcilar

repo: Multi-Robot-Path-Planning-on-Graphs


text:
%kksinc 'Output = sin(Input)/Input'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros ksinc.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'Input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kksinc(i, {'i','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% ksinc - Output = sin(Input)/Input
%
%  DESCRIPTION
% The "Sinc" operator calculates the sinc function (sin(x)/x) of each data
% point in the input data object, \fBInput".
%
% Executing "Sinc" runs the program \fIkarith1\fP with the -sinc flag.
%
%  "\fBData Type"" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "\fBMap Data"" 5
% .cI $DATAMANIP/repos/shared/man/sections/map_1input
%
%  "\fBValidity Mask"" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "\fBExplicit Location and Time Data"" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "\fBFailure Modes""
% .cI $DATAMANIP/repos/shared/man/sections/fail_1input
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::karith1
%
%  RESTRICTIONS
% Operations on complex data are not supported at this time.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kksinc(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kksinc(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output'};
maxval={0,0};
minval={0,0};
istoggle=[0,0];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'karith1"  -sinc'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [] = testUnpackTheta()

visibleDim = 5;  % width/height of input
visibleSize = visibleDim^2;

hiddenDimL1 = 3;  % width/height of hidden layer 1
hiddenSizeL1 = hiddenDimL1^2;
hiddenViewDimL1 = 3;  % width/height of local receptive area for hidden layer 1
hiddenViewSizeL1 = hiddenViewDimL1^2;  % local receptive area size for hidden layer 1
hiddenViewStepL1 = 1;  % step for receptive area moving across input
assert (((hiddenDimL1-1) * hiddenViewStepL1 + hiddenViewDimL1) == visibleDim);

hiddenDimL2 = 2;  % width/height of hidden layer 2
hiddenSizeL2 = hiddenDimL2^2;

hiddenViewDimL2 = 2;  % width/height of local receptive area for hidden layer 2
hiddenViewSizeL2 = hiddenViewDimL2^2;  % local receptive area size for hidden layer 2
hiddenViewStepL2 = 1;  % step for receptive area moving across hidden layer 1
assert (((hiddenDimL2-1) * hiddenViewStepL2 + hiddenViewDimL2) == hiddenDimL1);

[W1, b1] = initializeOneLayerParams(hiddenSizeL1, hiddenViewSizeL1);
[W2, b2] = initializeOneLayerParams(hiddenSizeL2, hiddenViewSizeL2);
[W3, b3] = initializeOneLayerParams(1, hiddenSizeL2);

theta = [W1(:);b1(:);W2(:);b2(:);W3(:);b3(:)];
[tW1, tb1, tW2, tb2, tW3, tb3] = unpackTheta(theta, hiddenSizeL1, hiddenViewSizeL1, hiddenSizeL2, hiddenViewSizeL2);

assert(all(W1(:) == tW1(:)));
assert(all(b1(:) == tb1(:)));
assert(all(W2(:) == tW2(:)));
assert(all(b2(:) == tb2(:)));
assert(all(W3(:) == tW3(:)));
assert(all(b3(:) == tb3(:)));

% Real sizes

visibleDim = 96;  % width/height of input
visibleSize = visibleDim^2;

hiddenDimL1 = 18;  % width/height of hidden layer 1
hiddenSizeL1 = hiddenDimL1^2;
hiddenViewDimL1 = 11;  % width/height of local receptive area for hidden layer 1
hiddenViewSizeL1 = hiddenViewDimL1^2;  % local receptive area size for hidden layer 1
hiddenViewStepL1 = 5;  % step for receptive area moving across input
assert(all(((hiddenDimL1-1) * hiddenViewStepL1 + hiddenViewDimL1) == visibleDim));

hiddenDimL2 = 8;  % width/height of hidden layer 2
hiddenSizeL2 = hiddenDimL2^2;

hiddenViewDimL2 = 11;  % width/height of local receptive area for hidden layer 2
hiddenViewSizeL2 = hiddenViewDimL2^2;  % local receptive area size for hidden layer 2
hiddenViewStepL2 = 1;  % step for receptive area moving across hidden layer 1
assert(all(((hiddenDimL2-1) * hiddenViewStepL2 + hiddenViewDimL2) == hiddenDimL1));

[W1, b1] = initializeOneLayerParams(hiddenSizeL1, hiddenViewSizeL1);
[W2, b2] = initializeOneLayerParams(hiddenSizeL2, hiddenViewSizeL2);
[W3, b3] = initializeOneLayerParams(1, hiddenSizeL2);

theta = [W1(:);b1(:);W2(:);b2(:);W3(:);b3(:)];
[tW1, tb1, tW2, tb2, tW3, tb3] = unpackTheta(theta, hiddenSizeL1, hiddenViewSizeL1, hiddenSizeL2, hiddenViewSizeL2);

assert(all(W1(:) == tW1(:)));
assert(all(b1(:) == tb1(:)));
assert(all(W2(:) == tW2(:)));
assert(all(b2(:) == tb2(:)));
assert(all(W3(:) == tW3(:)));
assert(all(b3(:) == tb3(:)));

end

author: zellyn


text:
function rs = db(s, dbmin)

%tstoolbox/@signal/db
%   Syntax:
%     * db(s, dbmin)
%
%   Compute decibel values of signal relative to a reference value that is
%   determined by the signal's yunit values below dbmin are set to dbmin.
%   If dbmin is ommited it is set to -120.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt
narginchk(1,2);
if nargin < 2, dbmin = -120; end

yu = yunit(s);

ref = dbref(yu);
scf = dbscale(yu);

c = db(s.core, ref, scf, dbmin); 	% call real working routine for parent core object
rs = signal(c, s);				% special constructor calling syntax for working routines

rs = setyunit(rs, unit(['dB' label(yu)]));
rs = addhistory(rs,  ['Calculated decibel values (stretch=' num2str(scf) ',ref=' num2str(ref) ')'] );
rs = addcommandlines(rs, 's = db(s', dbmin);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function cmap = autumnCmap(numGrays,numColors)
%
% cmap = autumnCmap(numGrays,numColors)
%
% Makes colormap array with:
%   gray scale - 1:numGrays
%   autumn colors - numGrays+1:numGrays+numColors
%
% djh 1/98
if ~exist('numGrays','var')
  numGrays=128;
end
if ~exist('numColors','var')
  numColors=96;
end

cmap = zeros(numGrays+numColors,3);
cmap(1:numGrays+numColors,:) = [gray(numGrays); autumn(numColors)];

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [kli, set_kli]= KLIndex(ar_est,ma_est,varx_est,ar,ma,varx,n_obs)

%KLIndex Kullback-Leibler index.
%   kli = KLIndex(a_est,b_est,varx_est,a,b,varx,nobs) is the exact Kullbach-Leibler
%   Index between two time series models for nobs observations.
%   The estimated proces a_est, b_est,var_est and the true process a,b,varx are assumed
%   to be normally distributed.
%
%   Possible second output argument: set of KLI's of AR models of increasing autoregressive
%   model order:
%   [kli set_kli] = KLIndex(a_est,b_est,varx_est,a,b,varx,nobs)

%S. de Waele, Februari 2003.

rc = ar_ma2rc(ar,ma,n_obs);
ar_order = length(rc)-1;
rc_est = ar_ma2rc(ar_est,ma_est,n_obs);
ar_est_order = length(rc_est)-1;

%Kullback-Leiber index of the estimated process
%with respect to the true process: kli_fef.
varp_est = varx_est*[1 cumprod(1-rc_est(2:end).^2)];

%Prediction errors (vareta)
set_ar_est = rc2arset(rc_est,0:ar_est_order);
cov = varx*arma2cor(ar,ma,ar_est_order);
covs = [cov(end:-1:2) cov];
for L = 0:ar_est_order,
    par_est = set_ar_est{1+L};
    sh = ar_est_order+L+1;
    vareta(L+1) = convol(par_est,convol(covs,par_est(end:-1:1)),sh,sh);
end

L = 0:ar_est_order;
set_kli_fef = cumsum(log(varp_est)) + cumsum(vareta./varp_est) ...
    + (n_obs-L-1).*(log(varp_est) + vareta./varp_est);

set_kli = set_kli_fef;
kli = set_kli(end);


text:
% 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 qp = mtimes(q1, q2)
%Quaternion.mtimes Multiply a quaternion object
%
% Q1*Q2   is a quaternion formed by Hamilton product of two quaternions.
% Q*V     is the vector V rotated by the quaternion Q
% Q*S     is the element-wise multiplication of quaternion elements by by the scalar S

    if isa(q1, 'Quaternion') & isa(q2, 'Quaternion')
    %QQMUL  Multiply unit-quaternion by unit-quaternion
    %
    %   QQ = qqmul(Q1, Q2)
    %
    %   Return a product of unit-quaternions.
    %
    %   See also: TR2Q

        % decompose into scalar and vector components
        s1 = q1.s;  v1 = q1.v;
        s2 = q2.s;  v2 = q2.v;

        % form the product
        qp = Quaternion([s1*s2-v1*v2' s1*v2+s2*v1+cross(v1,v2)]);

    elseif isa(q1, 'Quaternion') & isa(q2, 'double')

    %QVMUL  Multiply vector by unit-quaternion
    %
    %   VT = qvmul(Q, V)
    %
    %   Rotate the vector V by the unit-quaternion Q.
    %
    %   See also: QQMUL, QINV

        if length(q2) == 3
			qp = q1 * Quaternion([0 q2(:)']) * inv(q1);
            qp = qp.v(:);
        elseif length(q2) == 1
            qp = Quaternion(double(q1)*q2);
        else
            error('quaternion-vector product: must be a 3-vector or scalar');
        end

    elseif isa(q2, 'Quaternion') & isa(q1, 'double')
        if length(q1) == 3
            qp = q2 * Quaternion([0 q1(:)']) * inv(q2);
            qp = qp.v;
        elseif length(q1) == 1
            qp = Quaternion(double(q2)*q1);
        else
            error('quaternion-vector product: must be a 3-vector or scalar');
        end
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function ShowImageCell(imageCell, n, m, figurename, imageNames)
% ShowImageCell(imageCell, n, m, figurename, imageNames)
%
% Generate a figure with thumbnails of the images in the imageCell.
%
% imageCell:            Cell array with images which can be displayed
%                       with imshow.
% n:                    number of thumbnail rows per figure.
% m:                    number of thumbnail columns per figure.
% figurename:           Name of the figures (optional).
%
%     Jasper Uijlings - 2013

totImages = length(imageCell);
numFigures = ceil(totImages / (n * m));

if nargin < 4
    figurename = 'untitled';
end

if nargin < 5
    imageNames = cell(length(imageCell));
end

if ~iscell(imageNames)
    imageNamesC = cell(length(imageNames));
    for i=1:length(imageNames)
        imageNamesC{i} = sprintf('%g', imageNames(i));
    end
    imageNames = imageNamesC;
end

idx = 1;
screenSize = get(0, 'ScreenSize');

for i=1:numFigures
    if ispc
        figure('Position', [1, 1, screenSize(3), screenSize(4)], 'Name', figurename);
%         figure('Position', [1 49 1920 946] , 'Name', figurename);
    else
        figure('Position', [1, 1, screenSize(3)/2, screenSize(4)], 'Name', figurename);
    end
    clf;
    for j = 1:n * m
        if(idx <= totImages)
            subplot(n, m, j);
            imshow(imageCell{idx});
            xlabel(imageNames{idx});
            idx = idx + 1;
        end
    end
end

author: zouchuhang

repo: LayoutNet


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  battistn@tcnj.edu
% Date Created: October 8th, 2018
% Institution: TCNJ
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs / non-invariant beams *)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (battistn@tcnj.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: provides data from previous .vtk information for restarting a
%           simulation that has ended because of power failure, etc.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [current_time,cter,ctsave,U,V,xLag,yLag,xLag_P,yLag_P] = pass_Back_Data_For_Restart(dt,ctsave,print_dump,path_to_data);

%
% Automated
%
cter = ctsave * print_dump; % Total # of time-steps thus far up to and included last data point saved
current_time = cter * dt;   % Current time in simulation when last time-step was saved

%
% Read in LAGRANGIAN data for last and second to last timepoint for U,V,mVelocity,F_Poro, xLag,yLag,xLag_P,yLag_P
%
[xLag,yLag,xLag_P,yLag_P] = please_Give_Saved_Lag_Point_Info(ctsave,path_to_data);

%
% Read in EULERIAN data for last timepoint for U,V
%
[U,V] = please_Give_Saved_Eulerian_Info(ctsave,path_to_data);

%
% Update for next iteration (*pretends data from last timepoint was just saved*)
%
ctsave = ctsave + 1;                % Update for next ctsave number (always increases by 1 after data is saved)
current_time = current_time + dt;   % Update for next time-step
cter = cter + 1;                    % Update for next time-step

cd ../ % Go back into Example Directory

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: passes back EULERIAN data from last saved timepoint
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [U,V] = please_Give_Saved_Eulerian_Info(ctsave,path_to_data)

    % Points to desired data viz_IB2d data file for CURRENT time-step
    if ctsave<10
       numSim = ['000', num2str(ctsave) ];
    elseif ctsave<100
       numSim = ['00', num2str(ctsave) ];
    elseif ctsave<1000
       numSim = ['0', num2str(ctsave)];
    else
       numSim = num2str(ctsave);
    end

    % Imports (x,y) grid values and ALL EULERIAN DATA %
    %                      DEFINITIONS
    %          x: x-grid                y: y-grid
    %       Omega: vorticity           P: pressure
    %    uMag: mag. of velocity
    %    uX: mag. of x-Velocity   uY: mag. of y-Velocity
    %    U: x-directed velocity   V: y-directed velocity
    %    Fx: x-directed Force     Fy: y-directed Force
    %
    %  Note: U(j,i): j-corresponds to y-index, i to the x-index
    %
    %
    Eulerian_Flags(1) = 0;   % OMEGA
    Eulerian_Flags(2) = 0;   % PRESSURE
    Eulerian_Flags(3) = 0;   % uMAG
    Eulerian_Flags(4) = 0;   % uX (mag. x-component of velocity)
    Eulerian_Flags(5) = 0;   % uY (mag. x-component of velocity)
    Eulerian_Flags(6) = 1;   % uVEC (vector components of velocity: U,V)
    Eulerian_Flags(7) = 0;   % Fx (x-component of force )
    Eulerian_Flags(8) = 0;   % Fy (y-component of force)
    %
    [~,~,~,~,~,~,~,U,V,~,~] = import_Eulerian_Data(path_to_data,numSim,Eulerian_Flags);

    % TURNS OUT UNNECESSARY
    %U = U'; % For IB2d convention from .vtk
    %V = V'; % For IB2d convention from .vtk

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: passes back LAGRANGIAN data from last two saved timepoints
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag,xLag_P,yLag_P] = please_Give_Saved_Lag_Point_Info(ctsave,path_to_data)

    iP = ctsave - 1; % Previous time-point

    % Points to desired data viz_IB2d data file for CURRENT time-step
    if ctsave<10
       numSim = ['000', num2str(ctsave) ];
    elseif ctsave<100
       numSim = ['00', num2str(ctsave) ];
    elseif ctsave<1000
       numSim = ['0', num2str(ctsave)];
    else
       numSim = num2str(ctsave);
    end

    % Points to desired data viz_IB2d data file for PREVIOUS time-step data
    if iP<10
       numSim_Prev = ['000', num2str(iP) ];
    elseif iP<100
       numSim_Prev = ['00', num2str(iP) ];
    elseif iP<1000
       numSim_Prev = ['0', num2str(iP)];
    else
       numSim_Prev = num2str(iP);
    end

    % Imports immersed boundary positions %
    [xLag,yLag] = give_Lag_Positions(path_to_data,numSim);
    [xLag_P,yLag_P] = give_Lag_Positions(path_to_data,numSim_Prev);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: gives (x,y) positions of the immersed boundary at a single step
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Lag_Positions(path,numSim)

[xLag,yLag] = read_Lagrangian_Data_From_vtk(path,numSim);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in (x,y) positions of the immersed boundary from .vtk
%           format
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = read_Lagrangian_Data_From_vtk(path,simNums)

cd(path);

filename = ['lagsPts.' num2str(simNums) '.vtk'];  % desired lagsPts.xxxx.vtk file

fileID = fopen(filename);
if ( fileID== -1 )
    error('\nCANNOT OPEN THE FILE!');
end

str = fgets(fileID); %-1 if eof
if ~strcmp( str(3:5),'vtk');
    error('\nNot in proper VTK format');
end

% read in the header info %
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID); % Up to white space in previous lagsPts.X files (pre June 2021)

% Check whether VTK file has time info. This is a VTK file with time,
%       need to read the next 3 lines to have read in appropriate
%       number of header lines.
if ~strcmp( str(1), 'P')
	str = fgets(fileID);
    str = fgets(fileID);
    str = fgets(fileID);
end

% stores # of Lagrangian Pts. as stated in .vtk file
numLagPts = sscanf(str,'%*s %f %*s',1);

% read in the vertices %
[mat,count] = fscanf(fileID,'%f %f %f',3*numLagPts);
if count ~= 3*numLagPts
   error('\nProblem reading in Lagrangian Pts.');
end

mat = reshape(mat, 3, count/3); % Reshape vector -> matrix (every 3 entries in vector make into matrix row)
vertices = mat';              % Store vertices in new matrix

fclose(fileID);               % Closes the data file.

xLag = vertices(:,1);         % x-Lagrangian Pts.
yLag = vertices(:,2);         % y-Lagrangian Pts.

%---------------------------------------------------------------
% NEEDS TO BE COMMENTED OUT!
% (when compared to other file in data analysis toolbox)
%---------------------------------------------------------------
%cd ..;                        % Change directory back to ../viz_IB2d/ directory

clear mat str filename fileID count;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: imports all Eulerian Data at a single step
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x,y,Omega,P,uMag,uX,uY,U,V,Fx,Fy] = import_Eulerian_Data(path,numSim,Eulerian_Flags)

    %
    % EULERIAN FLAGS FOR WHAT GETS SPIT OUT %
    %
    %Eulerian_Flags(1):   OMEGA
    %Eulerian_Flags(2):   PRESSURE
    %Eulerian_Flags(3):   uMAG
    %Eulerian_Flags(4):   uX (mag. x-component of velocity)
    %Eulerian_Flags(5):   uY (mag. x-component of velocity)
    %Eulerian_Flags(6):   uVEC (x,y-components of velocity: U,V)
    %Eulerian_Flags(7):   Fx (x-component of force )
    %Eulerian_Flags(8):   Fy (y-component of force)
    %

%analysis_path = pwd

% read in Vorticity %
if Eulerian_Flags(1)
    strChoice = 'Omega'; first = 1;
    [Omega,x,y] = read_Eulerian_Data_From_vtk(path,numSim,strChoice,first);
else
    Omega=[];
end

% read in Pressure %
if Eulerian_Flags(2)
    strChoice = 'P'; first = 1;
    [P,x,y] = read_Eulerian_Data_From_vtk(path,numSim,strChoice,first);
else
    P=[];
end

% read in Velocity Magnitude %
if Eulerian_Flags(3)
    strChoice = 'uMag'; first = 1;
    [uMag,x,y] = read_Eulerian_Data_From_vtk(path,numSim,strChoice,first);
else
    uMag=[];
end

% read in x-directed Velocity Magnitude %
if Eulerian_Flags(4)
    strChoice = 'uX'; first = 1;
    [uX,x,y] = read_Eulerian_Data_From_vtk(path,numSim,strChoice,first);
else
    uX=[];
end

% read in y-directed Velocity Magnitude %
if Eulerian_Flags(5)
    strChoice = 'uY'; first = 1;
    [uY,x,y] = read_Eulerian_Data_From_vtk(path,numSim,strChoice,first);
else
    uY=[];
end

% read in x-directed Forces %
if Eulerian_Flags(7)
    strChoice = 'Fx'; first = 1;
    [Fx,x,y] = read_Eulerian_Data_From_vtk(path,numSim,strChoice,first);
else
    Fx=[];
end

% read in y-directed Forces %
if Eulerian_Flags(8)
    strChoice = 'Fy'; first = 1;
    [Fy,x,y] = read_Eulerian_Data_From_vtk(path,numSim,strChoice,first);
else
    Fy=[];
end

% read in Velocity Field %
if Eulerian_Flags(6)
    [U,V] = read_Eulerian_Velocity_Field_vtk(path,numSim);
else
    U=[];
    V=[];
end

% Default for x,y values
if max(Eulerian_Flags([1:5,7,8]))==0
    x=[];
    y=[];
end

%cd(analysis_path);

clear analysis_path;

clear strChoice first;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the velocity data field from .vtk format
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%s

function [U,V] = read_Eulerian_Velocity_Field_vtk(path,simNums)

analysis_path = pwd; % Store path to analysis folder!

cd(path);

filename = ['u.' num2str(simNums) '.vtk'];  % desired EULERIAN-DATA.xxxx.vtk file

fileID = fopen(filename);
if ( fileID== -1 )
    error('\nCANNOT OPEN THE FILE!');
end

str = fgets(fileID); %-1 if eof
if ~strcmp( str(3:5),'vtk')
    error('\nNot in proper VTK format');
end

% read in the header info %
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID);

% Check whether VTK file has time info. This is a VTK file with time,
%       need to read the next 3 lines to have read in appropriate
%       number of header lines.
if ~strcmp( str(1:3), 'DIM')
	str = fgets(fileID);
    str = fgets(fileID);
    str = fgets(fileID);
end

% Store grid info
Nx = sscanf(str,'%*s %f %*f %*s',1);
Ny = sscanf(str,'%*s %*f %f %*s',1);

% bypass lines in header %
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID);

% get formatting for reading in data from .vtk in fscanf %
strVec = '%f';
for i=2:3*Nx
    strVec = [strVec ' %f'];
end

% read in the vertices %
[e_Data,count] = fscanf(fileID,strVec,3*Nx*Ny);
if count ~= 3*Nx*Ny
   error('Problem reading in Eulerian Data.');
end

% reshape the matrix into desired data type %
e_Data = reshape(e_Data, 3, count/3); % Reshape (3*Nx*Nx,1) vector to (Nx*Nx,3) matrix
e_Data = e_Data';                     % Store vertices in new matrix

U = e_Data(:,1);       % Store U data
V = e_Data(:,2);       % Store V data

U = reshape(U,Nx,Ny)';  % Reshape (Nx*Nx,1) matrix to  (Nx,Nx)
V = reshape(V,Nx,Ny)';  % Reshape (Nx*Nx,1) matrix to  (Nx,Nx)

fclose(fileID);         % Closes the data file.

cd ..;                  % Change directory back to ../viz_IB2d/ directory

cd(analysis_path);      % Change directory back to Data Analysis Folder

clear filename fileID str strVec count analysis_path e_Data Nx;

author: nickabattista


text:
classdef TestPCTSignatures
    %TestPCTSignatures

    methods (Static)
        function test_1
            img = imread(fullfile(mexopencv.root(),'test','books_left.jpg'));

            pct = cv.PCTSignatures();
            sig = pct.computeSignature(img);
            validateattributes(sig, {'numeric'}, {'2d', 'size',[NaN 8]});

            sigs = pct.computeSignatures({img, img});
            validateattributes(sigs, {'cell'}, {'vector', 'numel',2});

            out = cv.PCTSignatures.drawSignature(img, sig);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_2
            initPoints = cv.PCTSignatures.generateInitPoints(2000, 'Uniform');
            pct = cv.PCTSignatures(initPoints, 400);
            pct.GrayscaleBits = 4;
            assert(isequal(pct.GrayscaleBits, 4));
        end
    end

end


text:
function [prob,sol,fmin] = lp_prob(varargin)
%LP_PROB  Return an OPTI LP
%
%   prob = lp_prob(no) return a pre-built optiprob of a saved LP.
%
%   [prob,sol,fmin] = lp_prob(no) returns the optimum solution and function
%   eval at the optimum
%
%   no = lp_prob() returns the number of problems available for testing.

%   (C) 2011 Jonathan Currie (IPL)

%Check if just returning no problems
if(nargin < 1)
    prob = 11; sol = []; fmin = [];
    return;
else
    no = varargin{1};
end

%Big switch yard
switch(no)
    case 1
        f = -[-1, 2]';
        A = [2, 1;-4, 4];
        b = [5, 5]';
        e = -[1, 1];
        prob = optiprob('f',f,'mix',A,b,e,'name','TestLP1');
        sol = [1.25;2.5];
        fmin = -3.75;

    case 2
        f = -[50, 100];
        A = [10, 5;4, 10; 1, 1.5];
        b = [2500, 2000, 450]';
        e = [-1, -1, -1];
        prob = optiprob('f',f,'mix',A,b,e,'name','TestLP2');
        sol = [187.5;125];
        fmin = -21875;

    case 3
        f = [40, 36];
        A = [5, 3];
        b = 45;
        e = 1;
        ub = [8, 10];
        prob = optiprob('f',f,'mix',A,b,e,'ub',ub,'name','TestLP3');
        sol = [8;5/3];
        fmin = 380;

    case 4
        f = [3, -7, -12];
        A = [-3, 6, 8;6, -3, 7;-6, 3, 3];
        b = [12, 8, 5]';
        e = [-1, -1, -1];
        prob = optiprob('f',f,'mix',A,b,e,'name','TestLP4');
        sol = [-0.0666666666;0.23333333333;1.3];
        fmin = -17.43333333333333;

    case 5
        f = [2, 3, 7, 7];
        A = [1, 1, -2, -5;-1, 2, 1, 4];
        b = [2, -3]';
        e = [1, 1];
        lb = zeros(4,1);
        ub = [30 100 20 1]';
        prob = optiprob('f',f,'mix',A,b,e,'bounds',lb,ub,'name','TestLP5');
        sol = [2;0;0;0];
        fmin = 4;

    case 6
        f = [1, 2, 3, 7, 8, 8];
        A = [5, -3, 2, -3, -1, 2; -1, 0, 2, 1, 3, -3;1, 2, -1, 0, 5, -1];
        b = [-5, -1, 3]';
        e = [1, 1, 1];
        lb = zeros(6,1);
        ub = 10*ones(6,1);
        prob = optiprob('f',f,'mix',A,b,e,'bounds',lb,ub,'name','TestLP6');
        sol = [0;1.5;0;0;0;0];
        fmin = 3;

    case 7
        n = 40;
        t = (0:n-1)';
        y = 3.5 -.2*t;
        b = y + 0.5*ones(size(y));
        m = [ones(n,1),t(:)];
        A = [m,-m,eye(n)];
        f = [sum(m),sum(-m),2*ones(1,n)];
        e = ones(n,1);
        lb = zeros(n+4,1);
        ub = [10, 10, 10, 10, 5*ones(1,n)];
        prob = optiprob('f',f,'mix',A,b,e,'bounds',lb,ub,'name','TestLP7');
        sol = [4;0;0;0.2;zeros(40,1)];
        fmin = 4.0000000000000044;

    case 8
        f = -[8, 15];
        A = [10, 21;2, 1];
        b = [156, 22]';
        e = [-1, -1];
        prob = optiprob('f',f,'mix',A,b,e,'name','TestLP8');
        sol = [9.5625;2.875];
        fmin = -119.625;

    case 9
        f = -[3, 13];
        A = [2, 9;11, -8];
        b = [40, 82]';
        e = [-1, -1];
        prob = optiprob('f',f,'mix',A,b,e,'name','TestLP9');
        sol = [9.2;2.4];
        fmin = -58.8;

    case 10
        f = -[592, 381, 273, 55, 48, 37, 23];
        A = [3534, 2356, 1767, 589, 528, 451, 304];
        b = 119567;
        e = -1;
        lb = zeros(7,1);
        ub = [100 50 33 20 77 44 20]';
        prob = optiprob('f',f,'mix',A,b,e,'bounds',lb,ub,'name','TestLP10');
        sol = [33.8333333333333;zeros(6,1)];
        fmin = -20029.33333333333;

    case 11
        f = [-1 -1 -3 -2 -2]';
        A = [-1 -1 1 1 0;
             1 0 1 -3 0];
        b = [30;30];
        ub = [40;1;inf;inf;1];
        prob = optiprob('f',f,'ineq',A,b,'ub',ub,'name','TestLP11');
        sol = [40;1;50.75;20.25;1];
        fmin = -235.7500;

    otherwise
        error('Problem not available or not implemented yet');
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% VL_TWISTER Random number generator
%    VL_TWISTER() is essentially equivalent to MATLAB native RAND()
%    when using the Twister random number generator. VL_TWISTER(),
%    VL_TWISTER(M,N,P,...)  and VL_TWISTER([M N P ...]) are equivalent
%    to RAND(), RAND(M,N,P,...)  and RAND([M N P ...]) respectively.
%
%    The state of the random generator can be seeded by
%    VL_TWISTER('STATE', X), where X is a DOUBLE scalar (this is
%    equivalent to RAND('TWISTER', X)). The state can be read by
%    VL_TWISTER('STATE') (equivalent to RAND('TWISTER')) and set by
%    VL_TWISTER('STATE', STATE) (equivalent to RAND('TWISTER',
%    STATE)). Here STATE is a vector of 625 elements of class
%    UINT32. Finally VL_TWISTER('STATE',KEY) seeds the generator by a
%    vector of DOUBLE of length not greater than 624.
%
%    VL_TWISTER() is slightly faster than RAND(). Moreover it can be
%    used to control the state of the random number generator used by
%    all VLFEAT functions.
%
%    See also: VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function f_x = ParFor10(in1)
%PARFOR10
%    F_X = PARFOR10(IN1)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    31-Mar-2020 17:03:12

r = in1(:,1);
s = in1(:,13);
u = in1(:,19);
ut = in1(:,20);
ux = in1(:,21);
uxx = in1(:,22);
uy = in1(:,23);
uyy = in1(:,24);
z = in1(:,7);
t2 = u.^2;
f_x = r.*-4.845016406907234+s.*6.96612441301113+ut.*5.457984690783633e-1-ux.*1.530748504418966e-2+uxx.*2.135644249305187-uy.*3.570878181824355e-2+uyy.*5.02823932467436-z.*1.323504040227272e2+r.*u.*6.191151455437648-s.*u.*9.105578580027213-t2.*u.*1.077295023277402e2+u.*ut.*5.65024831920482e-2+u.*ux.*1.830355084143775e-2-u.*uxx.*4.132096419489244+u.*uy.*4.442896435307375e-2-u.*uyy.*6.867807677976089+u.*z.*1.83605552916415e2+4.028421339788474e1;

author: dynamicslab

repo: SINDy-PI


text:
% Normalizes all of the field components so that the mode has
% unity power (Poynting vector integrated over the cross section.)

function [ex,ey,ez,hx,hy,hz] = normalize(dx,dy,EX,EY,EZ,HX,HY,HZ)

Z0 = 119.9169832*pi; % vacuum impedance

[nx,ny] = size(EX);

if (length(dx) ~= nx),
  dx = dx*ones(nx,1);
end
if (length(dy) ~= ny),
  dy = dy*ones(1,ny);
end

ii = zeros(nx+1,ny+1);
ii(:) = (1:(nx+1)*(ny+1));

i1 = ii(1:nx,2:ny+1);
i2 = ii(1:nx,1:ny);
i3 = ii(2:nx+1,1:ny);
i4 = ii(2:nx+1,2:ny+1);

HXp = (HX(i1) + HX(i2)+ HX(i3) + HX(i4))/4;
HYp = (HY(i1) + HY(i2)+ HY(i3) + HY(i4))/4;

SZ = Z0*(conj(EX).*HYp - conj(EY).*HXp + EX.*conj(HYp) - EY.*conj(HXp))/4;
dA = dx*dy;
N = sqrt(sum(SZ(:).*dA(:)));
ex = Z0*EX/N;
ey = Z0*EY/N;
ez = Z0*EZ/N;
hx = HX/N;
hy = HY/N;
hz = HZ/N;


text:
function order_num = tetrahedron_ncc_order_num ( rule )

%*****************************************************************************80
%
%% TETRAHEDRON_NCC_ORDER_NUM returns the order of an NCC rule for the tetrahedron.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Peter Silvester,
%    Symmetric Quadrature Formulae for Simplexes,
%    Mathematics of Computation,
%    Volume 24, Number 109, January 1970, pages 95-100.
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule.
%
%    Output, integer ORDER_NUM, the order (number of points) of the rule.
%
  suborder_num = tetrahedron_ncc_suborder_num ( rule );

  suborder(1:suborder_num) = tetrahedron_ncc_suborder ( rule, suborder_num );

  order_num = sum ( suborder(1:suborder_num) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Machine Learning Online Class
%  Exercise 6 | Support Vector Machines
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  exercise. You will need to complete the following functions:
%
%     gaussianKernel.m
%     dataset3Params.m
%     processEmail.m
%     emailFeatures.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% =============== Part 1: Loading and Visualizing Data ================
%  We start the exercise by first loading and visualizing the dataset.
%  The following code will load the dataset into your environment and plot
%  the data.
%

fprintf('Loading and Visualizing Data ...\n')

% Load from ex6data1:
% You will have X, y in your environment
load('ex6data1.mat');

% Plot training data
plotData(X, y);

%fprintf('Program paused. Press enter to continue.\n');
%pause;

%% ==================== Part 2: Training Linear SVM ====================
%  The following code will train a linear SVM on the dataset and plot the
%  decision boundary learned.
%

% Load from ex6data1:
% You will have X, y in your environment
load('ex6data1.mat');

fprintf('\nTraining Linear SVM ...\n')

% You should try to change the C value below and see how the decision
% boundary varies (e.g., try C = 1000)
C = 0.1;
model = svmTrain(X, y, C, @linearKernel, 1e-3, 20);
visualizeBoundaryLinear(X, y, model);

%fprintf('Program paused. Press enter to continue.\n');
%pause;

%% =============== Part 3: Implementing Gaussian Kernel ===============
%  You will now implement the Gaussian kernel to use
%  with the SVM. You should complete the code in gaussianKernel.m
%
fprintf('\nEvaluating the Gaussian Kernel ...\n')

x1 = [1 2 1]; x2 = [0 4 -1]; sigma = 2;
sim = gaussianKernel(x1, x2, sigma);

fprintf(['Gaussian Kernel between x1 = [1; 2; 1], x2 = [0; 4; -1], sigma = 0.5 :' ...
         '\n\t%f\n(this value should be about 0.324652)\n'], sim);

%fprintf('Program paused. Press enter to continue.\n');
%pause;

%% =============== Part 4: Visualizing Dataset 2 ================
%  The following code will load the next dataset into your environment and
%  plot the data.
%

fprintf('Loading and Visualizing Data ...\n')

% Load from ex6data2:
% You will have X, y in your environment
load('ex6data2.mat');

% Plot training data
plotData(X, y);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ========== Part 5: Training SVM with RBF Kernel (Dataset 2) ==========
%  After you have implemented the kernel, we can now use it to train the
%  SVM classifier.
%
fprintf('\nTraining SVM with RBF Kernel (this may take 1 to 2 minutes) ...\n');

% Load from ex6data2:
% You will have X, y in your environment
load('ex6data2.mat');

% SVM Parameters
C = 1; sigma = 0.1;

% We set the tolerance and max_passes lower here so that the code will run
% faster. However, in practice, you will want to run the training to
% convergence.
model = svmTrain(X, y, C, @(x1, x2) gaussianKernel(x1, x2, sigma));
visualizeBoundary(X, y, model);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =============== Part 6: Visualizing Dataset 3 ================
%  The following code will load the next dataset into your environment and
%  plot the data.
%

fprintf('Loading and Visualizing Data ...\n')

% Load from ex6data3:
% You will have X, y in your environment
load('ex6data3.mat');

% Plot training data
plotData(X, y);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ========== Part 7: Training SVM with RBF Kernel (Dataset 3) ==========

%  This is a different dataset that you can use to experiment with. Try
%  different values of C and sigma here.
%

% Load from ex6data3:
% You will have X, y in your environment
load('ex6data3.mat');

% Try different SVM Parameters here
[C, sigma] = dataset3Params(X, y, Xval, yval);

% Train the SVM
model = svmTrain(X, y, C, @(x1, x2) gaussianKernel(x1, x2, sigma));
visualizeBoundary(X, y, model);

fprintf('Program paused. Press enter to continue.\n');
pause;

author: worldveil

repo: coursera-ml


text:
function writehtk(file,d,fp,tc)
%WRITEHTK write data in HTK format []=(FILE,D,FP,TC)
%
% Inputs:
%    FILE = name of file to write (no default extension)
%       D = data to write: one row per frame
%      FP = frame period in seconds
%      TC = type code = the sum of a data type and (optionally) one or more of the listed modifiers
%             0  WAVEFORM     Acoustic waveform
%             1  LPC          Linear prediction coefficients
%             2  LPREFC       LPC Reflection coefficients:  -lpcar2rf([1 LPC]);LPREFC(1)=[];
%             3  LPCEPSTRA    LPC Cepstral coefficients
%             4  LPDELCEP     LPC cepstral+delta coefficients (obsolete)
%             5  IREFC        LPC Reflection coefficients (16 bit fixed point)
%             6  MFCC         Mel frequency cepstral coefficients
%             7  FBANK        Log Fliter bank energies
%             8  MELSPEC      linear Mel-scaled spectrum
%             9  USER         User defined features
%            10  DISCRETE     Vector quantised codebook
%            11  PLP          Perceptual Linear prediction
%            12  ANON
%            64  _E  Includes energy terms                  hd(1)
%           128  _N  Suppress absolute energy               hd(2)
%           256  _D  Include delta coefs                    hd(3)
%           512  _A  Include acceleration coefs             hd(4)
%          1024  _C  Compressed                             hd(5)
%          2048  _Z  Zero mean static coefs                 hd(6)
%          4096  _K  CRC checksum (not implemented yet)     hd(7) (ignored)
%          8192  _0  Include 0'th cepstral coef             hd(8)
%         16384  _V  Attach VQ index                        hd(9)
%         32768  _T  Attach delta-delta-delta index         hd(10)

%   Thanks to Scott Otterson for fixing a bug in writing ultra-long frames.

%      Copyright (C) Mike Brookes 2005
%      Version: $Id: writehtk.m,v 1.7 2006/11/16 08:49:29 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 2 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 can obtain a copy of the GNU General Public License from
%   ftp://prep.ai.mit.edu/pub/gnu/COPYING-2.0 or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fid=fopen(file,'w','b');
if fid < 0
    error(sprintf('Cannot write to file %s',file));
end
tc=bitset(tc,13,0);                 % silently ignore a checksum request

[nf,nv]=size(d);
nhb=10;                             % number of suffix codes
ndt=6;                              % number of bits for base type
hb=floor(tc*pow2(-(ndt+nhb):-ndt));
hd=hb(nhb+1:-1:2)-2*hb(nhb:-1:1);   % extract bits from type code
dt=tc-pow2(hb(end),ndt);            % low six bits of tc represent data type

if ~dt & (size(d,1)==1)             % if waveform is a row vector
    d=d(:);                         % ... convert it to a column vector
    [nf,nv]=size(d);
end

if hd(5)                            % if compressed
    dx=max(d,[],1);
    dn=min(d,[],1);
    a=ones(1,nv);                   % default compression factors for cols with max=min
    b=dx;
    mk=dx>dn;
    a(mk)=65534./(dx(mk)-dn(mk));   % calculate compression factors for each column
    b(mk)=0.5*(dx(mk)+dn(mk)).*a(mk);
    d=d.*repmat(a,nf,1)-repmat(b,nf,1); % compress the data
    nf=nf+4;                        % adjust frame count to include compression factors
end
fwrite(fid,nf,'long');              % write frame count
fwrite(fid,round(fp*1.E7),'long');  % write frame period (in 100 ns units)
if any(dt==[0,5,10]) | hd(5)        % write data as shorts
    if dt==5                        % IREFC has fixed scale factor
        d=d*32767;
        if hd(5)
            error('Cannot use compression with IREFC format');
        end
    end
    nby=nv*2;
    if nby<=32767
        fwrite(fid,nby,'short');        % write byte count
        fwrite(fid,tc,'short');         % write type code
        if hd(5)
            fwrite(fid,a,'float');        % write compression factors
            fwrite(fid,b,'float');
        end
        fwrite(fid,d.','short');        % write data array
    end
else
    nby=nv*4;
    if nby<=32767
        fwrite(fid,nby,'short');        % write byte count
        fwrite(fid,tc,'short');         % write type code
        fwrite(fid,d.','float');        % write data array
    end
end
fclose(fid);
if nby>32767
    delete(file);                       % remove file if byte count is rubbish
    error(sprintf('byte count of frame is %d which exceeds 32767 (is data transposed?)',nby));
end

author: singaxiong

repo: SignalGraph


text:
function output = callcsdp(interfacedata);

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
x0      = interfacedata.x0;
pars   = interfacedata.options.csdp;
ub      = interfacedata.ub;
lb      = interfacedata.lb;

% Bounded variables converted to constraints
if ~isempty(ub)
    [F_struc,K] = addStructureBounds(F_struc,K,ub,lb);
end

pars.printlevel=options.verbose;

if options.savedebug
    save csdpdebug F_struc c K pars
end

if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end
solvertime = tic;
if options.verbose==0 % to fix display bug reported from user
    evalc('[x_s,y_s,z_s,info]=csdp(-F_struc(:,2:end),-c,F_struc(:,1),K,pars);');
else
    [x_s,y_s,z_s,info]=csdp(-F_struc(:,2:end),-full(c),F_struc(:,1),K,pars);
end
solvertime = toc(solvertime);

% We solve dual problem with CSDP
Dual = x_s;
x    = y_s;

switch info
case 0
    problem = 0;
case 1
    problem = 2;
case 2
    problem = 1;
case {3,5,6,7,8,9}
    problem = 4;
case 4
    problem = 3;
otherwise
    problem = 9;
end
infostr = yalmiperror(problem,interfacedata.solver.tag);

% Save ALL data sent to solver
if options.savesolverinput
    solverinput.A = -F_struc(:,2:end);
    solverinput.c = F_struc(:,1);
    solverinput.b = -c;
    solverinput.K = K;
    solverinput.pars = pars;
else
    solverinput = [];
end

% Save ALL data from the solution?
if options.savesolveroutput
    solveroutput.x = x_s;
    solveroutput.y = y_s;
    solveroutput.info = info;
else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(x,Dual,z_s,problem,infostr,solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [A, a, B, b] = adjacencyAndCoords(Lmk,Frm,Fac)

nframes  = sum([Frm.used]);
nlmks    = sum([Lmk.used]);

% Motion
A = sparse(nframes,nframes);
a = zeros(nframes,3);

% Measurements
B = sparse(nframes+nlmks,nframes+nlmks);
b = zeros(nframes+nlmks,3);

for fac = find([Fac.used])

    frames = Fac(fac).frames;
    lmk = Fac(fac).lmk;

    switch Fac(fac).type
        case'motion'
            A(frames(1),frames(2))   = 1;
            a(frames(1),1:3) = Frm(frames(1)).state.x(1:3)';
            a(frames(2),1:3) = Frm(frames(2)).state.x(1:3)';
        case 'absolute'
        case 'measurement'
            B(frames,lmk+nframes) = 1;
            b(frames,1:3) = Frm(frames).state.x(1:3)';
            b(lmk+nframes,1:3) = Lmk(lmk).state.x(1:3)';
    end

end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [f, g] = gradFuncWrapper(x, func, grad, varargin)

% GRADFUNCWRAPPER Wrapper function to enable use of Carl Rasmussen's minimze function.

% OPTIMI

func = str2func(func);
grad = str2func(grad);
f = func(x', varargin{:});

if nargout > 1
  g = grad(x', varargin{:})';
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
function  write_boxes_to_file( stats, filename )

% Create the columns for different recall levels
str_header  = 'ncands\t';
str_pattern = '%d\t';
for ii=1:length(stats.overlap_levels)
    str_header = [str_header '\trec_at_' num2str(stats.overlap_levels(ii))]; %#ok<AGROW>
    str_pattern = [str_pattern '\t%f']; %#ok<AGROW>
end
str_header  = [str_header  '\n'];
str_pattern = [str_pattern '\n'];

fid = fopen(filename,'w');
fprintf(fid, str_header);
fprintf(fid, str_pattern, [stats.mean_n_masks; stats.rec_at_overlap]);
fclose(fid);

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
function depth_first_search(A,u,dfs_visitor,varargin)
% DEPTH_FIRST_SEARCH Fully wrap the Boost depth_first_search call
% including the dfs_visitor.
%
% depth_first_search(A,u,dfs_visitor) performs a depth first traversal
% of A starting from vertex u.  For each event defined by the dfs_visitor
% structure below, the visitor is called with the either the name of the
% vertex (u), or the edge index and it's source and target (ei,u,v).
%
% See http://www.boost.org/libs/graph/doc/DFSVisitor.html for a description
% of the events.
%
% dfs_visitor is a struct with the following optional fields
%    vis.initialize_vertex(u)
%    vis.start_vertex(u)
%    vis.discover_vertex(u)
%    vis.examine_edge(ei,u,v)
%    vis.tree_edge(ei,u,v)
%    vis.back_edge(ei,u,v)
%    vis.forward_or_cross_edge(ei,u,v)
%    vis.finish_vertex(u)
% Each dfs_visitor parameter should be a function pointer, which returns 0
% if the dfs should stop.  (If the function does not return anything, the
% dfs continues.)
%
% This method works on directed graphs.
% The runtime is O(V+E), excluding the complexity of the visitor
% operations.
%
% Realistically, this function must be used with the
% pass-by-reference/in-place modification library.
%
% ... = depth_first_search(A,u,vis,...) takes a set of
% key-value pairs or an options structure.  See set_matlab_bgl_options
% for the standard options.
%   options.full: compute the full dfs instead of the dfs of
%      the current component (see Note 1) [{0} | 1]
%
% Note 1: When computing the full dfs, the vertex u is ignored, vertex 1 is
% always used as the starting vertex.
%
% Note: this function does not depend upon the non-zero values of A, but
% only uses the non-zero structure of A.
%
% Example:
%   This example finds the distance to a single point and stops the search.
%   function dist_uv(A,u,v)
%     vstar = v;
%     dmap = ipdouble(zeros(size(A,1),1));
%     function stop=on_tree_edge(ei,u,v)
%       dmap(v) = dmap(u)+1;
%       stop = (v ~= vstar);
%     end
%     depth_first_search(A,u,struct('tree_edge',@on_tree_edge));
%   end
%
% See also DFS

% David Gleich
% Copyright, Stanford University, 2007-2008

%% History
%  2006-05-21: Initial version
%  2006-05-31: Added full2sparse check
%  2007-07-24: Fixed example
%  2008-10-07: Changed options parsing
%%

[trans check full2sparse] = get_matlab_bgl_options(varargin{:});
if full2sparse && ~issparse(A), A = sparse(A); end

if check
    % no additional input checks
    check_matlab_bgl(A,struct());
end

if trans, A = A'; end

% parse the optional parameters
full = 0;
if ~isempty(varargin)
    optionsu = merge_options(struct(),varargin{:});
    if (isfield(optionsu,'full'))
        full = optionsu.full;
    end
end

if full
    % 202 is the call for dfs with full searches
    bfs_dfs_vis_mex(A,u,dfs_visitor,202);
else
    % 201 is the call for dfs with partial searches
    bfs_dfs_vis_mex(A,u,dfs_visitor,201);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function y = cvx_isconvex( x, full )
error( nargchk( 1, 2, nargin ) );
y = cvx_vexity( x );
if nargin < 2,
    y = nnz( y < 0 ) == 0;
else
    y = cvx_reshape( y >= 0, x.size_ );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function [r, c] = partial_corr_coef(S, i, j, Y)
% PARTIAL_CORR_COEF Compute a partial correlation coefficient
% [r, c] = partial_corr_coef(S, i, j, Y)
%
% S is the covariance (or correlation) matrix for X, Y, Z
% where X=[i j], Y is conditioned on, and Z is marginalized out.
% Let S2 = Cov[X | Y] be the partial covariance matrix.
% Then c = S2(i,j) and r = c / sqrt( S2(i,i) * S2(j,j) )
%

% Example: Anderson (1984) p129
% S = [1.0 0.8 -0.4;
%     0.8 1.0 -0.56;
%     -0.4 -0.56 1.0];
% r(1,3 | 2) = 0.0966
%
% Example: Van de Geer (1971) p111
%S = [1     0.453 0.322;
%     0.453 1.0   0.596;
%     0.322 0.596 1];
% r(2,3 | 1) = 0.533

X = [i j];
i2 = 1; % find_equiv_posns(i, X);
j2 = 2; % find_equiv_posns(j, X);
S2 = S(X,X) - S(X,Y)*inv(S(Y,Y))*S(Y,X);
c = S2(i2,j2);
r = c / sqrt(S2(i2,i2) * S2(j2,j2));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip