text:
classdef TestDrawLineMatches
    %TestDrawLineMatches

    properties (Constant)
        img1 = fullfile(mexopencv.root(),'test','left01.jpg');
        img2 = fullfile(mexopencv.root(),'test','right01.jpg');
    end

    methods (Static)
        function test_1
            im1 = imread(TestDrawLineMatches.img1);
            im2 = imread(TestDrawLineMatches.img2);
            obj = cv.BinaryDescriptor();
            [klines1, feat1] = obj.detectAndCompute(im1);
            [klines2, feat2] = obj.detectAndCompute(im2);
            matcher = cv.BinaryDescriptorMatcher();
            m = matcher.match(feat1, feat2);

            im1 = cv.cvtColor(im1, 'GRAY2RGB');
            im2 = cv.cvtColor(im2, 'GRAY2RGB');
            out = cv.drawLineMatches(im1, klines1, im2, klines2, m);
            validateattributes(out, {class(im1)}, {'size',...
                [max(size(im1,1),size(im2,1)), size(im1,2)+size(im2,2), 3]});
        end

        function test_2
            im1 = imread(TestDrawLineMatches.img1);
            im2 = imread(TestDrawLineMatches.img2);
            obj = cv.BinaryDescriptor();
            [klines1, feat1] = obj.detectAndCompute(im1);
            [klines2, feat2] = obj.detectAndCompute(im2);
            matcher = cv.BinaryDescriptorMatcher();
            m = matcher.match(feat1, feat2);

            [~,ord] = sort([m.distance]);
            mask = false(size(ord));
            mask(ord(1:min(20,end))) = true;

            out = cv.drawLineMatches(im1, klines1, im2, klines2, m, ...
                'MatchesMask',mask, ...
                'OutImage',repmat(cat(2,im1,im2), [1 1 3]), ...
                'MatchColor',[255 0 0], 'SingleLineColor',[0 255 0], ...
                'NotDrawSingleLines',true);
            validateattributes(out, {class(im1)}, {'size',...
                [max(size(im1,1),size(im2,1)), size(im1,2)+size(im2,2), 3]});
        end

        function test_error_argnum
            try
                cv.drawLineMatches();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end


text:
function out = issing(f)
%ISSING   Test if a FUN object is built upon SINGFUN.
%   out = ISSING(F) returns logical true if its ONEFUN is a SINGFUN and false
%   otherwise.

% [TODO]: Should the name of this function go camel case? So far, all is*
% function are named without camel case.

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

out = isa(f.onefun, 'singfun');

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function makedir(dirname)

% create directory dirname

[path,name] = fileparts(dirname);

if ~exist(path,'dir')
    makedir(path);
end

if isunix
    if unix(['mkdir ' dirname])~=0
        error(['Cannot create directory ' dirname])
    end
else
    mkdir(path, name);
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [f] = spm_SHC_fx(x,v,P,varargin)
% equations of motion for Lotka-Volterra dynamics
% FORMAT [f] = spm_SHC_fx(x,v,P)
%
% x   - hidden states
% v   - exogenous inputs
% P.f - lateral connectivity
% P.k - rate [default 1]
%
% returns f = dx/dt = P.f*S(x) - x/8 + 1;
%              S(x) = 1./(1 + exp(-x))
%
% where C determines the order of unstable fixed points visited in the
% stable heteroclinic channel.
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_SHC_fx.m 3265 2009-07-10 14:02:22Z karl $

% intialise
%==========================================================================

% SHC states
%--------------------------------------------------------------------------
f.x  = P.f*spm_phi(x.x) - x.x/32 + 1;
f.c  = (spm_phi(x.x) - x.c)/64;

f    = spm_vec(f)*16;

author: spm


text:
function [X_norm, mu, sigma] = featureNormalize(X)
%FEATURENORMALIZE Normalizes the features in X
%   FEATURENORMALIZE(X) returns a normalized version of X where
%   the mean value of each feature is 0 and the standard deviation
%   is 1. This is often a good preprocessing step to do when
%   working with learning algorithms.

% You need to set these values correctly
X_norm = X;
mu = zeros(1, size(X, 2));
sigma = zeros(1, size(X, 2));

% ====================== YOUR CODE HERE ======================
% Instructions: First, for each feature dimension, compute the mean
%               of the feature and subtract it from the dataset,
%               storing the mean value in mu. Next, compute the
%               standard deviation of each feature and divide
%               each feature by it's standard deviation, storing
%               the standard deviation in sigma.
%
%               Note that X is a matrix where each column is a
%               feature and each row is an example. You need
%               to perform the normalization separately for
%               each feature.
%
% Hint: You might find the 'mean' and 'std' functions useful.
%

mu = mean(X);
sigma = std(X);
X_norm = (X - mu) ./ sigma;
#X_norm = (X - min(X)) ./ (max(X) - min(X))

% ============================================================

end

repo: machine-learning


text:
function out = NL_TSTL_dimensions(y,numBins,embedParams)
% NL_TSTL_dimensions Box counting, information, and correlation dimension of a time series.
%
% Computes the box counting, information, and correlation dimension of a
% time-delay embedded time series using the TSTOOL code 'dimensions'.
% This function contains extensive code for estimating the best scaling range to
% estimate the dimension using a penalized regression procedure.
%
% cf. TSTOOL, http://www.physik3.gwdg.de/tstool/
%
%---INPUTS:
% y, column vector of time series data
% numBins, maximum number of partitions per axis
% embedParams, embedding parameters to feed BF_Embed() for embedding the
%              signal in the form {tau,m}
%
%---OUTPUTS:
% A range of statistics are returned about how each dimension estimate changes
% with m, the scaling range in r, and the embedding dimension at which the best
% fit is obtained.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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/>.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%% Preliminaries, check inputs
% ------------------------------------------------------------------------------
doPlot = false; % plot outputs to screen

% (1) Maximum number of bins, numBins
if nargin < 2 || isempty(numBins)
    numBins = 50; % 50 points
    fprintf(1,'Using a default of 50 bins per axis\n');
end

% (2) Set embedding parameters to defaults
if nargin < 3 || isempty(embedParams)
    embedParams = {'ac','fnnmar'};
    fprintf(1,'Using default time-delay embedding parameters: autocorrelation and cao');
else
    if length(embedParams) ~= 2
        error('Embedding parameters are incorrectly formatted -- need {tau,m}')
    end
end

% ------------------------------------------------------------------------------
%% Embed the signal
% ------------------------------------------------------------------------------
% Convert to embedded signal object for TSTOOL
s = BF_Embed(y,embedParams{1},embedParams{2},1);

if ~isa(s,'signal') && isnan(s); % embedding failed
    error('Time-delay embedding for TSTOOL failed')
end

if size(data(s),2) < 3 % embedded with dimension < 3
    % note the 'true' predicted embedding dimension
    mopt = size(data(s),2);
    % embed with dimension m = 3
    s = BF_Embed(y,embedParams{1},3,1);
    fprintf(1,'Re-embedded with embedding dimension 3\n');
else
	mopt = size(data(s),2);
end

% ------------------------------------------------------------------------------
%% Run the TSTOOL function:
% ------------------------------------------------------------------------------
% This looks for the dimensions file in the tstoolbox/@signal/dimensions directory
if ~exist(fullfile('tstoolbox','@signal','dimensions'),'file')
    error('Cannot find the code ''dimensions'' from the TSTOOL package. Is it installed and in the Matlab path?');
end
try
    [bc, ~, co] = dimensions(s,numBins);
catch me
    error('Error running TSTOOL code dimensions: %s',me.message);
end

% We now have the scaling of the boxcounting dimension, D0, the information
% dimension D1, and the correlation dimension D2.

% It seems like there's not extra information in the in dimension, D1, by these
% estimates -- focus on the bc, D0 and correlation, D2.
% Can switch this here:
compute_in = 0;

% ------------------------------------------------------------------------------
%% Convert output to vectors
% ------------------------------------------------------------------------------
% calculations for each dimension up to the maximum:
% Seems to be in units of log_2 -- log2, so actually doesn't span a very wide
% range of length scales... -- although I think the maximum length is at 1,
% so I think it might be referring to fractions of the 'attractor size'...

% (1) Boxcounting dimension (BC)
bc_logN = data(bc); % this is log(N(r)) -- number within radius (look for this to scale linearly)
bc_logr = spacing(bc); % this is log(r) -- length scale
bc_logNlogr = bc_logN./(ones(size(bc_logN,2),1)*bc_logr)'; % look for this to be constant

% (2) Information dimension (IN)
if compute_in
    in_logl = data(in); % I think this is log(l(r)), look for this to scale linearly
    in_logr = spacing(in);
    in_logllogr = in_logl./(ones(size(in_logl,2),1)*in_logr)'; % look for this to be constant
end

% (3) Correlation dimension (CO)
co_logC = data(co); % look for this to scale linearly
co_logr = spacing(co);
co_logClogr = co_logC./(ones(size(co_logC,2),1)*co_logr)'; % look for this to be constant

if doPlot
    plot(bc_logr,bc_logNlogr,'o-')
    plot(bc_logr,bc_logN,'o-')
    input('BC')
    if compute_in
        plot(in_logr,in_logllogr,'o-')
        plot(in_logr,in_logl,'o-')
        input('IN')
    end
    plot(co_logr,co_logClogr,'o-')
    plot(co_logr,co_logC,'o-')
    input('CO')
end

% *** We now have to look for scaling regimes in each of these dimensions

% ------------------------------------------------------------------------------
%% Basic statistics on curves
% ------------------------------------------------------------------------------
out = struct;

%-------------------------------------------------------------------------------
%% How do curves change with m?
%-------------------------------------------------------------------------------
% Use SUB_mch

% Box counting dimension:
out = SUB_mch(bc_logr,bc_logN,'bc',out);

% Information dimension:
if compute_in
    out = SUB_mch(bc_logr,bc_logN,'in',out);
end

% Correlation dimension:
out = SUB_mch(co_logr,co_logC,'co',out);

% ------------------------------------------------------------------------------
%% What is the scaling range in r?
% ------------------------------------------------------------------------------
% ... and how good is the fit over this range?
% Use SUB_ScalingRange

% Box counting dimension, m = 1
out = SUB_ScalingRange(bc_logr,bc_logN(:,1),'scr_bc_m1',out);

% Box counting dimension m = 2
out = SUB_ScalingRange(bc_logr,bc_logN(:,2),'scr_bc_m2',out);

% Box counting dimension m = 3
out = SUB_ScalingRange(bc_logr,bc_logN(:,3),'scr_bc_m3',out);

% Box counting dimension m = chosen/given
out = SUB_ScalingRange(bc_logr,bc_logN(:,mopt),'scr_bc_mopt',out);

if compute_in
    % Information dimension, m = 1
    out = SUB_ScalingRange(in_logr,in_logl(:,1),'scr_in_m1',out);

    % Information dimension m = 2
    out = SUB_ScalingRange(in_logr,in_logl(:,2),'scr_in_m2',out);

    % Information dimension m = 3
    out = SUB_ScalingRange(in_logr,in_logl(:,3),'scr_in_m3',out);

    % Information dimension m = chosen/given
    out = SUB_ScalingRange(in_logr,in_logl(:,mopt),'scr_in_mopt',out);
end

% Correlation dimension, m = 1
out = SUB_ScalingRange(co_logr,co_logC(:,1),'scr_co_m1',out);

% Correlation dimension m = 2
out = SUB_ScalingRange(co_logr,co_logC(:,2),'scr_co_m2',out);

% Correlation dimension m = 3
out = SUB_ScalingRange(co_logr,co_logC(:,3),'scr_co_m3',out);

% Correlation dimension m = chosen/given
out = SUB_ScalingRange(co_logr,co_logC(:,mopt),'scr_co_mopt',out);

% ------------------------------------------------------------------------------
%% What m gives best fit?
% ------------------------------------------------------------------------------
% Use SUB_bestm

% Box counting dimension
out = SUB_bestm(bc_logr,bc_logN,'bc',out);

if compute_in
    % Information dimension
    out = SUB_bestm(in_logr,in_logl,'in',out);
end

% Correlation dimension
out = SUB_bestm(co_logr,co_logC,'co',out);

% ------------------------------------------------------------------------------
% ------------------------------------------------------------------------------
function out = SUB_mch(logr,logN,prefix,out)
    % looks at how changes with m. Since m will in general be different for each
    % different time series (i.e., if choosing an automatic method for
    % determining the embedding parameters), we have that m is at least
    % 3 here so that we can do statistics on at least these ones...

    % (i) on average the raw means at each m up to m = 3
    out.([prefix,'_meanm1']) = mean(logN(:,1));
    out.([prefix,'_meanm2']) = mean(logN(:,2));
    out.([prefix,'_meanm3']) = mean(logN(:,3));
    out.([prefix,'_meanmmax']) = mean(logN(:,end));

    % (ii) raw minimum at each m up to m = 3
    out.([prefix,'_minm1']) = min(logN(:,1));
    out.([prefix,'_minm2']) = min(logN(:,2));
    out.([prefix,'_minm3']) = min(logN(:,3));
    out.([prefix,'_minmmax']) = min(logN(:,end));

    % (iii) range at each m up to m = 3
    out.([prefix,'_range1']) = range(logN(:,1));
    out.([prefix,'_range2']) = range(logN(:,2));
    out.([prefix,'_range3']) = range(logN(:,3));
    out.([prefix,'_rangemmax']) = range(logN(:,end));

    % (iv) increments with m
    out.([prefix,'_mindiff']) = mean([min(logN(:,2))-min(logN(:,1)),min(logN(:,3))-min(logN(:,2))]);
    out.([prefix,'_meandiff']) = mean([mean(logN(:,2))-mean(logN(:,1)),mean(logN(:,3))-mean(logN(:,2))]);

    % (v) slopes and goodness of fit across whole r range
    [out.([prefix,'_lfitm1']), out.([prefix,'_lfitb1']), out.([prefix,'_lfitmeansqdev1'])] = subsublinfit(logr,logN(:,1)');
    [out.([prefix,'_lfitm2']), out.([prefix,'_lfitb2']), out.([prefix,'_lfitmeansqdev2'])] = subsublinfit(logr,logN(:,2)');
    [out.([prefix,'_lfitm3']), out.([prefix,'_lfitb3']), out.([prefix,'_lfitmeansqdev3'])] = subsublinfit(logr,logN(:,3)');
    [out.([prefix,'_lfitmmax']), out.([prefix,'_lfitbmax']), out.([prefix,'_lfitmeansqdevmax'])] = subsublinfit(logr,logN(:,end)');

    function [m, b, meansqdev] = subsublinfit(x,y)
        p1 = polyfit(x,y,1);
        pfit = p1(1)*x + p1(2);
        res = y-pfit;
        m = p1(1); % gradient
        b = p1(2); % intercept
        meansqdev = mean(res.^2);
    end
end

% ------------------------------------------------------------------------------
function out = SUB_ScalingRange(logr,logN,prefix,out)
    % determines the scaling range in r for some m
    % we remove points from either extreme in r until minimize some
    % error measure
    % two dimensional optimization: over starting point and ending
    % point.

    l = length(logr);
    stptr = 1:floor(l/2)-1; % must be in the first half (not necessarily, but for here)
    endptr = ceil(l/2)+1:l; % must be in second half (not necessarily, but for here)
    mybad = zeros(length(stptr),length(endptr));
    for i = 1:length(stptr)
        for j = 1:length(endptr)
            mybad(i,j) = lfitbadness(logr(stptr(i):endptr(j)),logN(stptr(i):endptr(j))');
        end
    end
    [a,b] = find(mybad == min(min(mybad))); % this defines the 'best' scaling range
%         plot(logr,logN,'o-b'); hold on; plot(logr(stptr(a):endptr(b)),logN(stptr(a):endptr(b)),'o-r');
%         hold off
%         disp(['keep from ' num2str(stptr(a)) ' to ' num2str(endptr(b))])

    out.([prefix,'_logrmin']) = logr(stptr(a)); % minimum of scaling range
    out.([prefix,'_logrmax']) = logr(endptr(b)); % maximum of scaling range
    out.([prefix,'_logrrange']) = logr(endptr(b)) - logr(stptr(a)); % range of scaling... range
    out.([prefix,'_pgone']) = (stptr(a)-1 + l - endptr(b))/length(logr); % number of points removed in process
                                               				  % of choosing the optimum scaling range

	% Do the optimum fit again
	x = logr(stptr(a):endptr(b));
	y = logN(stptr(a):endptr(b))';
    p = polyfit(x,y,1);
    pfit = p(1)*x + p(2);
    res = pfit - y;
    out.([prefix,'_meanabsres']) = mean(abs(res));
	out.([prefix,'_meansqres']) = mean(res.^2);
    out.([prefix,'_scaling_exp']) = p(1);
	out.([prefix,'_scaling_int']) = p(2);
	out.([prefix,'_minbad']) = min(min(mybad));

    function badness = lfitbadness(x,y)
        gamma = 0.02; % reguralization parameter gamma selected empirically, could be tweaked in future work
        p = polyfit(x,y,1);
        pfit = p(1)*x+p(2);
        res = pfit - y;
        badness = mean(abs(res)) - gamma*length(x); % want to still maximize length(x)
    end
end

% ------------------------------------------------------------------------------
function out = SUB_bestm(logr,logNN,prefix,out)
    % logNN is a matrix... logN is a vector for a given m
	% determines the scaling range in r for some m
    % we remove points from either extreme in r until minimize some
    % error measure
    % two dimensional optimization: over starting point and ending
    % point.

	store_scalingexps = zeros(size(logNN,2),1);
	store_meansqres = zeros(size(logNN,2),1);
	for k = 1:size(logNN,2);
		logN = logNN(:,k); % take this element

        l = length(logr);
        stptr = 1:floor(l/2)-1; % must be in the first half (not necessarily, but for here)
        endptr = ceil(l/2)+1:l; % must be in second half (not necessarily, but for here)
        mybad = zeros(length(stptr),length(endptr));
        for i = 1:length(stptr)
            for j = 1:length(endptr)
                mybad(i,j) = lfitbadness(logr(stptr(i):endptr(j)),logN(stptr(i):endptr(j))');
            end
        end
        [a, b] = find(mybad == min(min(mybad))); % this defines the 'best' scaling range

		% Do the optimum fit again
		x = logr(stptr(a):endptr(b));
		y = logN(stptr(a):endptr(b))';
        p = polyfit(x,y,1);
        pfit = p(1)*x+p(2);
        res = pfit - y;
		% subout.meanabsres = mean(abs(res));
		% subout.meansqres = mean(res.^2);
		% subout.scaling_exp = p(1);
		% subout.scaling_int = p(2);
		% subout.minbad = min(min(mybad));

		store_scalingexps(k) = p(1);
		store_meansqres(k) = mean(res.^2);
	end

	out.([prefix,'_minscalingexp']) = min(store_scalingexps);
	out.([prefix,'_meanscalingexp']) = mean(store_scalingexps);
	out.([prefix,'_maxscalingexp']) = max(store_scalingexps);
	out.([prefix,'_mbestfit']) = find(store_meansqres == min(store_meansqres),1,'first');

    function badness = lfitbadness(x,y)
        gamma = 0.02; % reguralization parameter gamma selected empirically, could be tweaked in future work
        p = polyfit(x,y,1);
        pfit = p(1)*x + p(2);
        res = pfit - y;
        badness = mean(abs(res)) - gamma*length(x); % want to still maximize length(x)
    end
end

end


text:
function areasmooth = SmoothAreaOutliers(area,filterorder,maxfreq,maxerr)

areasmooth = LowPassFilterArea(area,filterorder,maxfreq);
isoutlier = isinf(area) | isnan(area) | abs(areasmooth - area) > maxerr;
[starts,ends] = get_interval_ends(isoutlier);
ends = ends - 1;
areasmooth = area;
nframes = numel(area);
for i = 1:numel(starts),
  if starts(i) == 1 && ends(i) == nframes,
    break;
  elseif starts(i) == 1,
    areasmooth(starts(i):ends(i)) = area(ends(i)+1);
  elseif ends(i) == nframes,
    areasmooth(starts(i):ends(i)) = area(starts(i)-1);
  else
    areasmooth(starts(i):ends(i)) = (area(starts(i)-1)+area(ends(i)+1))/2;
  end
end

author: kristinbranson

repo: JAABA


text:
% BETALPR - Beta Distribution - Log Probability Ratio
% Copyright (c) 1998, Harvard University. Full copyright in the file Copyright
%
%   [ lpr ] = betalpr(p1,p2,alpha,beta)
%
% returns the log of the p(p1) / p(p2) when both
% are distributed Beta(alpha,beta).
%
% See also: METROP, *LPR

function [ lpr ] = betalpr(p1,p2,alpha,beta)
lpr = (alpha-1) * log(p1/p2) + (beta-1) * log((1-p1)/(1-p2)) ;

author: Sable

repo: mcbench-benchmarks


text:
function plot_laplace( x,params,hAx,plot_num,fontsize )
% plot the laplace distribution with parameter "u" and "b"
%
% the distribution is given by:
%
%        p(x) = 1/(2*b)*exp(-abs(x-u)/b)
%
% format:   plot_laplace( x,params,hAx,plot_num,fontsize )
%
% input:    x         - X axis, for the plot
%           params    - the distribution parameter, RMS error and VAR or CRB
%           hAx       - where to plot the distribution curve
%           plot_num  - since the curve is added with a text to the axes,
%                       this parameter specifies where the text should be displayed
%                       and what color to choose for the curve
%           fontsize  - size of the font of the text, default 9
%
%
% example:  plot_laplace( x,fit_ML_laplace(data),hAx,3 )
%

% init graphic parameters
switch plot_num
case 1, cl = [1 0 0];
case 2, cl = [0 1 0];
case 3, cl = [1 0 1];
case 4, cl = [0 1 1];
case 5, cl = [0.5 0.5 0];
end
if ~exist('fontsize')
    fontsize = 9;
end
p       = plot_num*0.15 + 0.3;
ylimit  = ylim(hAx);
xlimit  = xlim(hAx);
fnc_txt = '\it{1/2\bfb}\rm \bf{\cdot e}^{-\mid\bf{x-\mu}\mid/\bf{b}}\rm';
if isfield( params,'VAR' )
    txt     = sprintf( '\\fontsize{%d}\\bfLaplace PDF\\rm with %s:  %s\n\\mu = %g    b = %g\nVAR(b) = %1.3g\nRMS err = %1.3g\n',...
        fontsize,params.type,fnc_txt,params.u,params.b,params.VAR_b,params.RMS );
else
    txt     = sprintf( '\\fontsize{%d}\\bfLaplace PDF\\rm with %s:  %s\n\\mu = %g    b = %g\nCRB(b) = %1.3g\nRMS err = %1.3g\n',...
        fontsize,params.type,fnc_txt,params.u,params.b,params.CRB_b,params.RMS );
end

% calculate distribution
u       = params.u;
b       = params.b;
y       = 1/(2*b)*exp(-abs(x-u)/b);

% plot and write the text
line( 'parent',hAx,'xdata',x,'ydata',y,'linewidth',2,'color',cl );
hTxt    = text( 1.2*mean(x),ylimit(2)*p,txt,'parent',hAx );
ext     = get( hTxt,'Extent' );
line( ext(1) - xlimit/15 ,(ext(2)+ext(4)/2)*[1 1],'linewidth',3,'color',cl );
% ext(3)  = ext(3)+xlimit(2)/15;
% rectangle( 'position',ext );

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function Y = khatri_rao(A, B)
%KHATRI_RAO Performs Khatri-Rao product (column-wise Kronecker product).
%Syntax:
%   Y = KHATRI_RAO(A, B)

[r1, c1] = size(A);
[r2, c2] = size(B);
if c1 ~= c2
    error('Column size mismatch');
end

Y = zeros(r1 * r2, c1);
for i = 1:c1
    Y(:,i) = kron(A(:,i), B(:,i));
end

end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
function genes = drosTargets(tf),

% DROSTARGETS returns a list of hand selected enchanced targets for given tf
%
% Usage:
%   genes = drosTargets(tf)
% where tf is one of 'bin', 'mef2', 'tin', 'twi'
%
% COPYRIGHT : Antti Honkela, 2007

% 'FBgn0029082', 'FBgn0036459',
tingenes1 = {'FBgn0003117', 'FBgn0033459', 'FBgn0051368'};
tingenenames = {'pannier', 'CG12744', 'CG31368'};

tingenes2 = {'FBgn0003117', 'FBgn0015904', 'FBgn0030054', ...
	     'FBgn0027611', 'FBgn0029082'};

mefgenes1 = {'FBgn0000412', ...  % ???
	     'FBgn0000416', ...
	     'FBgn0000497', ...
	     'FBgn0001624', ...  % Twist bound nearby, probably for another gene
	     'FBgn0004401', ...
	     'FBgn0010109'};     % Lots of other activity nearby, probably for other genes

mefgenes2 = {'FBgn0031717', 'FBgn0035158', 'FBgn0036282', 'FBgn0004028', ... % !!!
	     'FBgn0032919', ... % This or CG1099
	     'FBgn0033688'};

twigenes1 = {'FBgn0003062', 'FBgn0003520', 'FBgn0014033', 'FBgn0025776', ...
	     'FBgn0028999', 'FBgn0029146'};
%twigenes1 = {'FBgn0003062', 'FBgn0003313', 'FBgn0004646', 'FBgn0003520', 'FBgn0014033', 'FBgn0020391'};

twigenes2 = {'FBgn0030749', ... % Mef2 and bin activity nearby
	     'FBgn0040108', ... % !!!
	     'FBgn0003062', ... % Mef2 nearby, probably for other genes
	     'FBgn0003308', 'FBgn0003607', 'FBgn0004635'}; % Good targets

bingenes1 = {'FBgn0004880', 'FBgn0010225', 'FBgn0015011', 'FBgn0026562', ...
	     'FBgn0033661', 'FBgn0035091', 'FBgn0038881', 'FBgn0039704'};

bingenes2 = {'FBgn0022073', 'FBgn0037240', 'FBgn0039268', 'FBgn0001197', ...
	     'FBgn0003890', 'FBgn0015011'};

bapgenes1 = {'FBgn0031349'};

switch tf,
 case 'twi',
  %genes = twigenes1([1 2 3 4 6]);
  genes = twigenes2;
 case 'mef2',
  %genes = mefgenes1([2 3 5]);
  %genes = mefgenes2([1 2 3 5 6]);
  genes = mefgenes2;
 case 'tin',
  %genes = tingenes1;
  genes = tingenes2;
 case 'bin',
  %genes = bingenes1(1:6);
  genes = bingenes2;
 %case 'bin1',
 % genes = bingenes1;
 otherwise,
  error('Unknown TF');
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function plot(catalogObject, varargin)
    %CATALOG.PLOT Plot hypocenters in map view and cross-sections
    %
    %
    %   Optional name/value pairs:
    %     'nsigma' - controls how zoomed-in the axes are (default
    %     5)

    % Glenn Thompson 2014/06/01
    if all(isnan(catalogObject.lat))
        warning('No hypocenter data to plot');
        return
    end
    p = inputParser;
    p.addParamValue('nsigma', '5', @isstr);
    p.parse(varargin{:});
    nsigma = p.Results.nsigma;

    % change region
    region = get_region(catalogObject, nsigma);

    % Compute Marker Size
    symsize = get_symsize(catalogObject);

    figure;
    set(gcf,'Color', [1 1 1]);

    % lon-lat plot
    axes('position',[0.05 0.45 0.5 0.5]);
    scatter(catalogObject.lon, catalogObject.lat, symsize);
    grid on;
    %set(gca, 'XLim', [region(1) region(2)]);
    %set(gca, 'YLim', [region(3) region(4)]);
    xlabel('Longitude');

    % depth-longitude
    axes('position',[0.05 0.05 0.5 0.35]);
    scatter(catalogObject.lon, catalogObject.depth, symsize);
    ylabel('Depth (km)');
    xlabel('Longitude');
    grid on;
    set(gca, 'YDir', 'reverse');
    %set(gca, 'XLim', [region(1) region(2)]);

    % depth-lat
    axes('position',[0.6 0.45 0.35 0.5]);
    scatter(catalogObject.depth, catalogObject.lat, symsize);
    xlabel('Depth (km)');
    %set(gca, 'XDir', 'reverse');
    ylabel('Latitude');
    grid on;
    %set(gca, 'YLim', [region(3) region(4)]);

%             % world map
%             load plotEarthquakes/plates.mat
%             coast = load('coast');
%             figure
%             worldmap world
%             setm(gca,'mlabelparallel',-90,'mlabellocation',90)
%             plotm(coast.lat,coast.long,'Color','k')
%             plotm(lat,lon,'LineWidth',2)
%
%             %% Find the first plate
%             % Look for the first NaN and stop there.
%             ind = find(isnan(lat),1,'first')
%             plotm(lat(1:ind),lon(1:ind),'Color','red','Linewidth',3)

end

author: geoscience-community-codes

repo: GISMO


text:
function a = downshift_inverse ( n )

%*****************************************************************************80
%
%% DOWNSHIFT_INVERSE returns the inverse of the DOWNSHIFT matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real A(N,N), the inverse.
%
  a = upshift ( n );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function nsub = subset_enum ( n )

%*****************************************************************************80
%
%% SUBSET_ENUM enumerates the subsets of a set with N elements.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of elements in the set.
%    N must be at least 0.
%
%    Output, integer NSUB, the number of distinct elements.
%
  nsub = 2^n;

  return
end


text:
function CompareClassifierWithLabels(clfile,ncomparisons,explist)
% Compares the classifier's predictions with a classifier trained with the current labels.

if nargin <2
  ncomparisons = 3;
end

if nargin<3,
  explist = {};
end

Q = load(clfile);
configfile = Q.configfilename;

data = JLabelData(configfile,'openmovie',false);

if isempty(explist),
  data.SetClassifierFileName(clfile);
  allOrig = getCurrentScores(data);

else
  data.SetClassifierFileNameWoExp(clfile);
  for ndx = 1:numel(explist)
    [success,msg] = data.AddExpDir(explist{ndx});
    if ~success,
      warning(msg);
      continue;
    end
  end

  scores = [];
  for endx = 1:data.nexps
    sfn = data.GetFile('scores',endx);
    Q = load(sfn);
    allScores = Q.allScores;

    for ndx = 1:numel(allScores.scores);
      ts = allScores.tStart(ndx);
      te = allScores.tEnd(ndx);
      scores = [scores allScores.scores{ndx}(ts:te)];
    end
  end
  allOrig = scores;

end

allNew = {};
for ndx = 1:ncomparisons
  data.Train();
  allNew{ndx} = getCurrentScores(data);
end

origdiff = [];
newdiff = [];
origdiffpos = [];
newdiffpos = [];
origdiffneg = [];
newdiffneg = [];
sn = data.windowdata.scoreNorm;
bins = linspace(-sn,sn,20);
origmat = zeros(2,2,ncomparisons);
newmat = zeros(2,2,ncomparisons);
for ndx = 1:ncomparisons
  t = histc(allOrig - allNew{ndx},[-inf bins inf]);
  t(1) = []; t(end) = [];
  origdiff(end+1,:) = t;

  origmat(1,1,ndx) = nnz(allOrig>0 & allNew{ndx}>0);
  origmat(1,2,ndx) = nnz(allOrig>0 & allNew{ndx}<0);
  origmat(2,1,ndx) = nnz(allOrig<0 & allNew{ndx}>0);
  origmat(2,2,ndx) = nnz(allOrig<0 & allNew{ndx}<0);

  popul = setdiff(1:ncomparisons,ndx);
  if isempty(popul),
    ndx1 = 1;
  else
    ndx1 = randsample(setdiff(1:ncomparisons,ndx),1);
  end
  t = histc(allNew{ndx1} - allNew{ndx},[-inf bins inf]);
  t(1) = []; t(end) = [];
  newdiff(end+1,:) = t;

  newmat(1,1,ndx) = nnz(allNew{ndx1}>0 & allNew{ndx}>0);
  newmat(1,2,ndx) = nnz(allNew{ndx1}>0 & allNew{ndx}<0);
  newmat(2,1,ndx) = nnz(allNew{ndx1}<0 & allNew{ndx}>0);
  newmat(2,2,ndx) = nnz(allNew{ndx1}<0 & allNew{ndx}<0);

  pos = allOrig>0;
  t = histc(allOrig(pos) - allNew{ndx}(pos),[-inf bins inf]);
  t(1) = []; t(end) = [];
  origdiffpos(end+1,:) = t;

  t = histc(allNew{ndx1}(pos) - allNew{ndx}(pos),[-inf bins inf]);
  t(1) = []; t(end) = [];
  newdiffpos(end+1,:) = t;

  pos = allOrig<0;
  t = histc(allOrig(pos) - allNew{ndx}(pos),[-inf bins inf]);
  t(1) = []; t(end) = [];
  origdiffneg(end+1,:) = t;

  t = histc(allNew{ndx1}(pos) - allNew{ndx}(pos),[-inf bins inf]);
  t(1) = []; t(end) = [];
  newdiffneg(end+1,:) = t;

end

origmat = round(mean(origmat,3));
newmat = round(mean(newmat,3));

figure;
subplot(2,2,[1 2]); hold on;
for ndx = 1:ncomparisons
  plot(bins/sn,origdiff(ndx,:),'r');
  plot(bins/sn,newdiff(ndx,:),'b');
end

numPosOrig = nnz(allOrig>0);
numNegOrig = nnz(allOrig<0);
numPosNew = 0; numNegNew = 0;
for ndx = 1:numel(allNew)
  numPosNew = numPosNew + nnz(allNew{ndx}>0);
  numNegNew = numNegNew + nnz(allNew{ndx}<0);

end
numPosNew = numPosNew/numel(allNew);
numNegNew = numNegNew/numel(allNew);

title({'Red: Orig - new, Blue: new - new'});
msg = {...
 sprintf('Frames predicted by the orig, pos :%d neg:%d',numPosOrig,numNegOrig),...
 sprintf('Frames predicted by the new (avg), pos :%d neg:%d',round(numPosNew),round(numNegNew)),...
 sprintf('Orig Vs New'),...
 sprintf('           New Pos          New Neg     '),...
 sprintf('Orig Pos     %d          %d     ',origmat(1,1),origmat(1,2)),...
 sprintf('Orig Neg     %d          %d     ',origmat(2,1),origmat(2,2)),...
 sprintf('New Vs New'),...
 sprintf('           New Pos          New Neg     '),...
 sprintf('New Pos     %d          %d     ',newmat(1,1),newmat(1,2)),...
 sprintf('New Neg     %d          %d     ',newmat(2,1),newmat(2,2)),...
};
warndlg(msg);

subplot(2,2,3); hold on;
for ndx = 1:ncomparisons
  plot(bins/sn,origdiffpos(ndx,:),'r');
  plot(bins/sn,newdiffpos(ndx,:),'b');
end
title('On frames predicted positive by orig');

subplot(2,2,4); hold on;
for ndx = 1:ncomparisons
  plot(bins/sn,origdiffneg(ndx,:),'r');
  plot(bins/sn,newdiffneg(ndx,:),'b');
end
title('On frames predicted negative by orig');

function scores = getCurrentScores(data)
scores = [];
for expi = 1:data.nexps
  allScores = data.PredictWholeMovie(expi);
  for ndx = 1:numel(allScores.scores);
    ts = allScores.tStart(ndx);
    te = allScores.tEnd(ndx);
    scores = [scores allScores.scores{ndx}(ts:te)];
  end
end


text:
function mn = i4cvv_size ( m, nr )

%*****************************************************************************80
%
%% I4CVV_SIZE determines the size of an I4CVV.
%
%  Discussion:
%
%    An I4CVV is a "vector of vectors" of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in the array.
%
%    Input, integer NR(M), the size of each row.
%
%    Output, integer MN, the size of the cell array.
%
  mn = sum ( nr(1:m) );

  return
end


text:
function outputarray = linarray(start,spacing,numpoints)
% outputarray = linarray(start,spacing,numpoints)
%
% The function linarray generates an row vector based on a start value, a
% spacing, and a number of points.  For example: linarray(1,.1,5) returns:
% [1.0000    1.1000    1.2000    1.3000    1.4000]
%
% 3/2/11 (c) James F. Mack

outputarray = [start:spacing:start+(numpoints-1)*spacing];


text:
function polygons = splitPolygons(polygon)
%SPLITPOLYGONS Convert a NaN separated polygon list to a cell array of polygons.
%
%   POLYGONS = splitPolygons(POLYGON);
%   POLYGON is a N-by-2 array of points, possibly with pairs of NaN values.
%   The functions separates each component separated by NaN values, and
%   returns a cell array of polygons.
%
%   See also
%   polygons2d
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2007-10-12, using Matlab 7.4.0.287 (R2007a)
% Copyright 2007-2022 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas

if iscell(polygon)
    % case of a cell array
    polygons = polygon;

elseif sum(isnan(polygon(:))) == 0
    % single polygon -> no break
    polygons = {polygon};

else
    % find indices of NaN couples
    inds = find(sum(isnan(polygon), 2) > 0);

    % number of polygons
    N = length(inds) + 1;
    polygons = cell(N, 1);

    % iterate over NaN-separated regions to create new polygon
    inds = [0; inds; size(polygon, 1)+1];
    for i = 1:N
        polygons{i} = polygon((inds(i)+1):(inds(i+1)-1), :);
    end
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [ pr, pq ] = plane_normal_basis_3d ( pp, normal )

%*****************************************************************************80
%
%% PLANE_NORMAL_BASIS_3D finds two perpendicular vectors in a plane in 3D.
%
%  Discussion:
%
%    The normal form of a plane in 3D is:
%
%      PP is a point on the plane,
%      N is a normal vector to the plane.
%
%    The two vectors to be computed, PQ and PR, can be regarded as
%    the basis of a Cartesian coordinate system for points in the plane.
%    Any point in the plane can be described in terms of the "origin"
%    point PP plus a weighted sum of the two vectors PQ and PR:
%
%      P = PP + a * PQ + b * PR.
%
%    The vectors PQ and PR have unit length, and are perpendicular to N
%    and to each other.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real PP(3), a point on the plane.  (Actually,
%    we never need to know these values to do the calculation!)
%
%    Input, real NORMAL(3), a normal vector N to the plane.  The
%    vector must not have zero length, but it is not necessary for N
%    to have unit length.
%
%    Output, real PR(3), a vector of unit length,
%    perpendicular to the vector N and the vector PQ.
%
%    Output, real PQ(3), a vector of unit length,
%    perpendicular to the vector N and the vector PR.
%
  dim_num = 3;
%
%  Compute the length of NORMAL.
%
  normal_norm = norm ( normal );

  if ( normal_norm == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PLANE_NORMAL_BASIS_3D - Fatal error!\n' );
    fprintf ( 1, '  The normal vector is 0.\n' );
    error ( 'PLANE_NORMAL_BASIS_3D - Fatal error!' );
  end
%
%  Find a vector PQ that is normal to NORMAL and has unit length.
%
  pq = r8vec_any_normal ( dim_num, normal );
%
%  Now just take the cross product NORMAL x PQ to get the PR vector.
%
  pr = r8vec_cross_product_3d ( normal, pq );

  pr_norm = norm ( pr );

  pr(1:dim_num) = pr(1:dim_num) / pr_norm;

  return
end


text:
function y = nigpdf(x, alpha, beta, mu, delta)
%NIGPDF Probability density function (pdf) for Normal-Inverse-Gaussian distribution.
%   Y = NIGPDF(X, ALPHA, BETA, MU, DELTA) returns the pdf of the
%   Normal-Inverse-Gaussian distribution with the parameter BETA which
%   determines the skewness, shape parameter ALPHA, location parameter MU
%   and scale parameter DELTA, evaluated at the values in X.
%
%   ALPHA, BETA, MU, DELTA must be scalar values.
%   ALPHA and DELTA must be positive values.
%   ALPHA > |BETA| must hold.
%
%   Default values for ALPHA = 1, BETA = 0; MU = 0; DELTA = 1.
%
%   See also NIGCDF, NIGINV, NIGRND, NIGSTATS.
%
%   References:
%   [1] Prause, K. (1999). The Generalized Hyperbolic Model

% -------------------------------------------------------------------------
%
% Allianz, Group Risk Controlling
% Risk Methodology
% Koeniginstr. 28
% D-80802 Muenchen
% Germany
% Internet:    www.allianz.de
% email:       ralf.werner@allianz.de
%
% Implementation Date:  2006 - 05 - 01
% Author:               Dr. Ralf Werner
%
% -------------------------------------------------------------------------

%% Default input values
if nargin < 2
    alpha = 1;
end
if nargin < 3
    beta = 0;
end
if nargin < 4
    mu = 0;
end
if nargin < 5
    delta = 1;
end

%% Constraints for the parameters
if alpha <= 0
    error('ALPHA must be positive.');
end
if delta <= 0
    error('DELTA must be positive.');
end
if abs(beta)>=alpha
    error('|BETA| must be smaller than ALPHA');
end

% transform input into column vector
[nx, mx] = size(x);
x = reshape(x, nx * mx, 1);

% calculate leading factor
q = (delta * alpha / pi) * exp(delta * sqrt(alpha^2 - beta^2));
xbar = x - mu;
z = sqrt(delta^2 + xbar .* xbar);

% modified Bessel function of the third kind
K = besselk(1, (alpha * z));

y = q ./ z .* K .* exp(beta * (x - mu));

% transform input vector back to input format
y = reshape(y, nx, mx);
y(isinf(x)) = 0;


text:
function point = projPointOnLine3d(point, line)
%PROJPOINTONLINE3D Project a 3D point orthogonally onto a 3D line.
%
%   PT2 = projPointOnLine3d(PT, LINE).
%   Computes the (orthogonal) projection of 3D point PT onto the 3D line
%   LINE.
%
%   Function works also for multiple points and lines. In this case, it
%   returns multiple points.
%   Point PT1 is a N-by-3 array, and LINE is a N-by-6 array.
%   Result PT2 is a N-by-3 array, containing coordinates of orthogonal
%   projections of PT1 onto lines LINE.
%
%
%   See also:
%   projPointOnLine, distancePointLine3d
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 2012-08-23.
%

%   HISTORY

% direction vector of the line
vx = line(:, 4);
vy = line(:, 5);
vz = line(:, 6);

% difference of point with line origin
dx = point(:,1) - line(:,1);
dy = point(:,2) - line(:,2);
dz = point(:,3) - line(:,3);

% Position of projection on line, using dot product
delta = vx .* vx + vy .* vy + vz .* vz;
tp = (dx .* vx + dy .* vy + dz .* vz) ./ delta;

% convert position on line to cartesian coordinates
point = [line(:,1) + tp .* vx, line(:,2) + tp .* vy, line(:,3) + tp .* vz];

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ fn ] = getBatchWrapper_ucf101_imgs(opts, numThreads, trainopts)
fn = @(imdb,batch, moreopts) getBatch(imdb,batch,opts,numThreads, trainopts, moreopts) ;

end

% -------------------------------------------------------------------------
function [im,labels] = getBatch(imdb, batch, opts, numThreads, trainopts, moreopts)
% -------------------------------------------------------------------------
opts.nFrames = imdb.images.nFrames(batch);
images = cell(1,numel(opts.nFrames));
labels = imdb.images.label(batch) ;
for k = 1:numel(opts.nFrames)
  images{k}  = strcat([imdb.images.name{batch(k)}(1:end-4) filesep], images{k} ) ;
end
if ~isempty(moreopts) || nargin < 6
  for f = fieldnames(moreopts)'
    f = char(f) ;
    trainopts.(f) = moreopts.(f);
  end
end

opts.nFramesPerVid = trainopts.nFramesPerVid;
opts.numAugments = trainopts.numAugments;
opts.frameSample = trainopts.frameSample;
opts.augmentation =  trainopts.augmentation;
if isfield(trainopts, 'temporalStride')
  opts.temporalStride =  trainopts.temporalStride(randi(numel(trainopts.temporalStride))) ; % shuffle;
end
if isfield(trainopts, 'keepFramesDim')
  keepFramesDim = trainopts.keepFramesDim;
else
  keepFramesDim = false;
end
if isfield(trainopts, 'cheapResize')
  opts.cheapResize =  trainopts.cheapResize;
end
if isfield(trainopts, 'prefetch')
    opts.prefetch = trainopts.prefetch;
end
if isfield(trainopts, 'frameList')
    opts.frameList = trainopts.frameList;
end

im = cnn_ucf101_get_frame_batch(images, opts, ...
                            'numThreads', numThreads, ...
                            'augmentation', opts.augmentation , ...
                          'imageDir', imdb.imageDir) ;

if iscell(im), return; end

labels = imdb.images.label(batch) ;

if nargout == 1 % to work with dag code
  if ndims(im) > 4 && ~keepFramesDim
    sz = size(im);
    nFrames = sz(5:end);
    im = permute(im, [1 2 3 5 4]);
    im = reshape(im, sz(1), sz(2), sz(3), []);
  end
  if size(im,3) > 1
    im = {'input', im(:,:,end-2:end,:,:), 'label', labels, 'input_flow', im(:,:,1:end-3,:,:)} ;
  else
    im = {'input', im, 'label', labels} ;
  end
end

end

author: feichtenhofer

repo: st-resnet


text:
function rx = rxRecenterRx(rx,recenter3D);
%
% rx = rxRecenterRx(rx,[recenter3D]);
%
% Recenter a prescription using the mouse,
% on the Rx Figure.
%
% The clicked-on location will be the center
% of the prescription in the two dimensions
% being displayed in the Rx Figure's image, but
% the third dimension will be ignored. E.g.,
% if you're looking at the prescription on axials,
% it will recenter along the sagittal and coronal,
% but not axial, dimensions. This can be overridden
% if the recenter3D argument is entered as 1 [default
% is 0, ignore 3rd dimension].
%
%
% ras 03/05.
% ras 01/06: think I fixed why it wasn't always working -- the mapping
% of necessary rotations for the xform matrix (which goes:
%   coronal rotation, axial rotation, sagittal rotation
% ), and the orientation button order (which goes:
%   axial, coronal, sagittal
%   and which corresponds to the mrVista conventions) is weird. This
% seems to work, as long as you set 3D recenter to be the default.
% setting it to 0 is nice, if the rx is orthogonal to where you're viewing
% the volume, but weird if it's oblique w.r.t. the slice axis.
if ieNotDefined('rx')
    cfig = findobj('Tag','rxControlFig');
    rx = get(cfig,'UserData');
end

if ieNotDefined('recenter3D')
    recenter3D = 1;
end

if ~ishandle(rx.ui.rxFig)
    % exit quietly
    return
end

%%%%%params
slice = get(rx.ui.volSlice.sliderHandle,'Value');
ori = findSelectedButton(rx.ui.volOri);

pt = get(gca,'CurrentPoint');
x = pt(1,1); y = pt(1,2);

%%%%%figure out offset from current location
% center of current rx (vol coordinate system)
cen = rx2vol(rx,[rx.rxDims/2]')';

% desired location (vol coordinate system)
switch ori
    case 1, tgt = [y slice x]; % axial
    case 2, tgt = [slice y x]; % coronal
    case 3, tgt = [x y slice]; % sagittal
    otherwise, return; % not yet implemented
end

offset = tgt - cen;

if recenter3D ~= 1
    % ignore dimension orthogonal to view
    order = [2 1 3]; % weird mapping of orientation buttons -> trans axes
    offset(order(ori)) = 0;
end

%%%%%update existing rx
[trans rot scale skew] = affineDecompose(rx.xform);
trans = trans + offset;
newXform = affineBuild(trans, rot, scale, skew);
rx = rxSetXform(rx, newXform);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [limitV,valV,ntcpLim,tcpLim,toolTipC] = getLimits(critS,type,strNum,doseNum,scale)

global planC;
indexS = planC{end};
structListC = {planC{indexS.structures}.structureName};

% Get DVH
[doseBinsV,volHistV] = getDVH(strNum,doseNum,planC);
doseBinsV = doseBinsV.*scale;

%Get val, limit
fieldsC = fieldnames(critS.(type));
valV = zeros(1,numel(fieldsC));
passV = zeros(1,numel(fieldsC));
nFail = 0;
ntcpLim = [];
tcpLim = [];
for m = 1:numel(fieldsC)
    lim = critS.(type).(fieldsC{m}).limit;
    limitV(m) = lim(1); %TEMP SKIP RANGE
    if strcmpi(fieldsC{m},'ntcp')
        ntcpLim = [limitV(m),m];
    end
    if strcmpi(fieldsC{m},'tcp')
        tcpLim = [limitV(m),m];
    end
    if isfield(critS.(type).(fieldsC{m}),'function')
        fn = critS.(type).(fieldsC{m}).function;
        if isfield(critS.(type).(fieldsC{m}),'additionalInputs')
            additionalInputs = critS.(type).(fieldsC{m}).additionalInputs;
            if numel(additionalInputs)==1  %%%TEMP FIX%%%%%%
             valV(m) = feval(fn,doseBinsV,volHistV,additionalInputs);
            else
                valV(m) = feval(fn,doseBinsV,volHistV,additionalInputs(1),additionalInputs(2));
            end
        else
            valV(m) = feval(fn,doseBinsV,volHistV);
        end
    end
    passV(m) = valV(m)<=limitV(m);
    toolTipC{m} = [structListC{strNum},' ',fieldsC{m},': ',num2str(valV(m))];
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [v,H,var,n,vsat,idx]=robust_spp(v0,H0,var0,nv,vsat,idx,flag)

[nr,nl]=size(H0);
ave_distance=zeros(nv,1); exc=zeros(nv,1); n=0; %#ok
v=zeros(nr,1); H=zeros(nr,nl); var=zeros(nr,nr);

if flag==0
    MAX_DISTANCE=3; %for pseudorange
else
    MAX_DISTANCE=0.2; %for doppler
end

for i=1:nv
    vi=abs(v0(i));
    if vi<MAX_DISTANCE,continue;end
    v_other=abs(v0);v_other(i)=[];v_other(nv:end,:)=[];

    if vi>3*sum(v_other)/(nv-1)
        exc(i)=1;
        vsat(idx(i))=0;
    end
end

for i=1:nv
    if exc(i)==1,continue;end
    v(n+1)=v0(i);
    H(n+1,:)=H0(i,:);
    var(n+1,n+1)=var0(i,i);
    n=n+1;
end

return

repo: GINav


text:
function [mi3] = cc2mi3(cc)
% Convert volume from cubic centimeters to cubic miles.
% Chad Greene 2012
mi3 = cc*2.3991275858e-16;


text:
function varargout = struct2single(varargin)

% This function is a backward compatibility wrapper. It allows existing
% MATLAB scripts that do not use the new FieldTrip ft_xxx function naming
% scheme to work with recent versions of the FieldTrip toolbox.
%
% Please look in ft_xxx for the help of the function that you are looking
% for, where xxx is the name of the function that you were looking for.

% Copyright (C) 2009, Robert Oostenveld

prefix    = 'ft_';
funname   = mfilename;
warning([upper(mfilename) ' is only a compatibility wrapper, which will soon be removed. Please instead call ' upper(prefix) upper(funname) '.']);
funhandle = str2func([prefix funname]);
[varargout{1:nargout}] = funhandle(varargin{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


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

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

% Last Modified by GUIDE v2.5 18-Feb-2011 17:46:26

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

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

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

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

% Update handles structure
guidata(hObject, handles);

date1 = date;
set(handles.date1, 'String', date);
initialize_gui(hObject, handles, false);

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

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

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

% --- Executes during object creation, after setting all properties.
function numbera_CreateFcn(hObject, eventdata, handles)
% hObject    handle to numbera (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function numbera_Callback(hObject, eventdata, handles)
% hObject    handle to numbera (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of numbera as text
%        str2double(get(hObject,'String')) returns contents of numbera as a double

numbera = str2double(get(hObject, 'String'));
if isnan(numbera)
    set(hObject, 'String', 0);
    errordlg('Input must be a number','Error');
end

% Save the new density value
handles.metricdata.numbera = numbera;
guidata(hObject,handles)

% --- Executes during object creation, after setting all properties.
function numberb_CreateFcn(hObject, eventdata, handles)
% hObject    handle to numberb (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function numberb_Callback(hObject, eventdata, handles)
% hObject    handle to numberb (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of numberb as text
%        str2double(get(hObject,'String')) returns contents of numberb as a double

numberb = str2double(get(hObject, 'String'));
if isnan(numberb)
    set(hObject, 'String', 0);
    errordlg('Input must be a number','Error');
end

% Save the new density value
handles.metricdata.numberb = numberb;
guidata(hObject,handles)

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

result = handles.metricdata.numbera + handles.metricdata.numberb;
set(handles.result, 'String', result)

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

result = handles.metricdata.numbera - handles.metricdata.numberb;
set(handles.result, 'String', result)

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

result = handles.metricdata.numbera / handles.metricdata.numberb;
set(handles.result, 'String', result)

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

result = handles.metricdata.numbera * handles.metricdata.numberb;
set(handles.result, 'String', result)

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

result = (log(handles.metricdata.numbera));
set(handles.result, 'String', result)

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

result = (1 / handles.metricdata.numbera);
set(handles.result, 'String', result)

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

result = (handles.metricdata.numbera)^3;
set(handles.result, 'String', result)

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

result = (handles.metricdata.numbera)^(handles.metricdata.numberb);
set(handles.result, 'String', result)

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

result = (sqrt(handles.metricdata.numbera));
set(handles.result, 'String', result)

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

result = log10(handles.metricdata.numbera);
set(handles.result, 'String', result)

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

result = (handles.metricdata.numbera)^(1 /handles.metricdata.numberb);
set(handles.result, 'String', result)

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

N = 1;
for n = 1:1:handles.metricdata.numbera
    N = N.*n;
end
result = N;
set(handles.result, 'String', result)

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

result = (handles.metricdata.numbera)^(2);
set(handles.result, 'String', result)

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

initialize_gui(gcbf, handles, true);

function initialize_gui(fig_handle, handles, isreset)
% If the metricdata field is present and the reset flag is false, it means
% we are we are just re-initializing a GUI by calling it from the cmd line
% while it is up. So, bail out as we dont want to reset the data.
if isfield(handles, 'metricdata') && ~isreset
    return;
end

handles.metricdata.numbera = 0;
handles.metricdata.numberb  = 0;

set(handles.numbera, 'String', handles.metricdata.numbera);
set(handles.numberb,  'String', handles.metricdata.numberb);
set(handles.result, 'String', 0);
set(handles.date1, 'String', date);

% Update handles structure
guidata(handles.figure1, handles);

function result_Callback(hObject, eventdata, handles)
% hObject    handle to result (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of result as text
%        str2double(get(hObject,'String')) returns contents of result as a double

% --- Executes during object creation, after setting all properties.
function result_CreateFcn(hObject, eventdata, handles)
% hObject    handle to result (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


text:
function [] = writefort19( fort19dat, finame )
%
%

if ( nargin == 1 )
  finputname = 'fort.19_1' ;
else
  finputname = finame ;
end

fid = fopen(finputname,'w') ;

fprintf( fid, '%d\n', fort19dat.ETIMINC ) ;

valpernode = size( fort19dat.Val,2 ) ;
str = [];
for ll = 1: valpernode
    str = [str '%.8f '] ;
end
str = [str '\n' ] ;

fprintf( fid, str, fort19dat.Val' ) ;

fclose(fid) ;
%EOF
end

author: CHLNDDEV


text:
function [x,info] = opti_zmumps(A,b)
%OPTI_ZMUMPS Solve a system of complex sparse linear equations using MUMPS
%
%   x = opti_zmumps(A,b) solves the system of complex sparse linear equations Ax = b.
%
%   THIS IS A WRAPPER FOR ZMUMPS USING THE MEX INTERFACE
%   See supplied License

%   Copyright (C) 2013 Jonathan Currie (IPL)

t = tic;

%Initialize MUMPS
id = struct('SYM',0,'JOB',-1,'ICNTL',zeros(1,40)-9998,'CNTL',zeros(1,15)-9998,...
            'PERM_IN',-9999,'COLSCA',-9999,'ROWSCA',-9999,'RHS',-9999,'INFOG',zeros(1,40)-9998,...
            'RINFOG',zeros(1,40)-9998,'VAR_SCHUR',-9999,'SCHUR',-9999,'INST',-9999,'SOL',-9999,...
            'REDRHS',-9999,'PIVNUL_LIST',-9999,'MAPPING',-9999,'SYM_PERM',-9999,'UNS_PERM',-9999,'TYPE',0);

id = zmumps(id);

id.Type = 2; %arithmetic type
id.JOB = 6; %analyze + factorize + solve

%Setup and Solve
id.RHS = b;
if(~issparse(A))
    optiwarn('opti:sparse','The A matrix should be sparse, correcting: [sparse(A)]');
    A = sparse(A);
end
id = zmumps(id,A);
x = id.SOL;

%Destroy mumps instance
id.JOB = -2;
id = zmumps(id); %#ok<NASGU>

info.Iterations = [];
info.Time = toc(t);
info.Algorithm = 'MUMPS: Complex Sequential Build';
info.Status = [];


text:
clearvars;
close all;
clc;
rng default;

mf = spx.graphics.Figures;

mf.new_figure();
[signal, fs] = spx.dsp.dtmf({'4', '5', '0', '7'});
time = (0:(numel(signal) - 1)) / fs;
plot(1e3*time, signal);
xlabel('Time (ms)');
ylabel('Amplitude');
grid on;
saveas(gcf, 'images/dtmf_4507.png');

mf.new_figure();
envelope_signal = envelope(signal, 80,'rms');
plot(1e3*time, envelope_signal);
xlabel('Time (ms)');
ylabel('Amplitude');
grid on;
saveas(gcf, 'images/dtmf_4507_envelope.png');

mf.new_figure();
pulsewidth(envelope_signal,fs)
saveas(gcf, 'images/dtmf_4507_pulses.png');

mf.new_figure();
periodogram(signal, [], [], fs);
saveas(gcf, 'images/dtmf_4507_periodogram.png');

% find the power spectrum and corresponding frequencies
[pxx,f]=periodogram(signal,[],[],fs);

n_freqs = 5;
% Let's find the most important frequencies
% 'NPeaks', n_freqs
[peak_values, peak_freqs] = findpeaks(pxx, f, 'SortStr','descend', 'MinPeakHeight', max(pxx) / 10);
peak_freqs = round(peak_freqs');
% Let's also mark them on a plot
mf.new_figure();
findpeaks(pxx, f, 'SortStr','descend','MinPeakHeight', max(pxx) / 10);

mf.new_figure();
spectrogram(signal, [], [], [], fs, 'yaxis');
% restrict the y-axis between 500Hz to 1500 Hz.
ylim([0.5 1.5]);
saveas(gcf, 'images/dtmf_4507_spectrogram.png');

mf.new_figure();
window_length = floor(fs * 50 / 1000);
overlap_length = floor(window_length / 2);
n_fft = 2^nextpow2(window_length);
spectrogram(signal,hamming(window_length),overlap_length,n_fft, fs, 'yaxis');
ylim([0.5 1.5]);
saveas(gcf, 'images/dtmf_4507_spectrogram_50ms.png');

mf.new_figure();
overlap_length = floor(0.8 * window_length );
spectrogram(signal,hamming(window_length),overlap_length,n_fft, fs, 'yaxis', 'MinThreshold', -50);
ylim([0.5 1.5]);
saveas(gcf, 'images/dtmf_4507_spectrogram_50ms_40ms_50db.png');

mf.new_figure();
spectrogram(signal,hamming(window_length),overlap_length,n_fft, fs, 'yaxis', 'MinThreshold', -50, 'reassigned');
ylim([0.5 1.5]);
saveas(gcf, 'images/dtmf_4507_spectrogram_50ms_40ms_50db_reassigned.png');

author: indigits

repo: sparse-plex


text:
function S = sample_mc_endstate(startprob, trans, endprob)
% SAMPLE_MC_ENDSTATE Generate a random sequence from a Markov chain until enter the endstate.
% seq = sample_mc(startprob, trans, endprob)

% add an end state
Q = size(trans,1);
transprob = zeros(Q,Q+1);
end_state = Q+1;
for i=1:Q
  for j=1:Q
    transprob(i,j) = (1-endprob(i)) * trans(i,j);
  end
  transprob(i,end_state) = endprob(i);
  %assert(approxeq(sum(transprob(i,:)), 1))
end

S = [];
S(1) = sample_discrete(startprob);
t = 1;
p = endprob(S(t));
stop = (S(1) == end_state);
while ~stop
  S(t+1) = sample_discrete(transprob(S(t),:));
  stop = (S(t+1) == end_state);
  t = t + 1;
end
S = S(1:end-1); % don't include end state


text:
function t = basicinfo(fid,seekset,t)
% T = BASICINFO(FID,SEEKSET,T);
%   FID = valid filepointer to a sqdfile
%   SEEKSET = starting point of file read
%   T = SQDHANDLE OBJECT/STRUCTURE
% Gets the basic info from the file pointed to by fid and returns
% a SQDHANDLE OBJECT if second argument is passed as sqdhandle or
% if no arguemnt is passed or a structure is passsed, it returns a
% structure
if nargin<1
    error('First argument must be valid file-pointer to a sqd-file');
elseif nargin<2
    seekset = -1;
end;

% Get offset of basic information
fseek( fid, 16, seekset );
basic_offset = fread(fid,1,'long');
newline = sprintf('\n');
% Read basic information
fseek( fid, basic_offset, seekset );
tmp.Version	= fread(fid,1,'int');
tmp.Revision	= fread(fid,1,'int');
tmp.SystemID	= fread(fid,1,'int');
tmp.SystemName	= fread(fid,128,'uchar');
tmp.SystemName	= setstr(tmp.SystemName(1:...
    min(find(tmp.SystemName==0)))');
if ~isempty(tmp.SystemName)
    newindx = findstr(tmp.SystemName,newline);
    tmp.SystemName(newindx(end):end) = '';
end;
tmp.ModelName	= fread(fid,128,'uchar');
tmp.ModelName	= setstr(tmp.ModelName(1:...
    min(find(tmp.ModelName==0)))');
if ~isempty(tmp.ModelName)
    newindx = findstr(tmp.ModelName,newline);
    if size(newindx) > 0, tmp.ModelName(newindx(end):end) = ''; end;
end;
tmp.ChannelCount	= fread(fid,1,'int');
tmp.Comment	= fread(fid,256,'uchar');
tmp.Comment	= setstr(tmp.Comment');
if nargin>=2
    taglist = fieldnames(tmp);
    if isa(t,'sqdhandle')
        for i =1:length(taglist)
            t = set(t,taglist{i},getfield(tmp,taglist{i}));
        end;
    elseif istruct(t)
        for i =1:length(taglist)
            t = setfield(t,taglist{i},getfield(tmp,taglist{i}));
        end;
    end;
else
    t = tmp;
end;


text:
function Out = LimitFigSize(FigH, Prop, Ext)
    % LimitFigSize - Set and get minimum or maximum figure size
    % Out = LimitFigSize(FigH, Prop, Ext)
    % INPUT:
    %   FigH: Handle of a Matlab figure.
    %         Optional, default: current figure.
    %   Prop: String, command:
    %         'min': Set minimum extent.
    %         'max': Set maximum extent.
    %         'get': Reply a struct with fields MinSize and MaxSize. If not limited,
    %                the empty matrix [] is replied.
    %         'clear': Clear the min and max limits.
    %         Optional, default: 'min'
    %   Ext:  [width, height] vector to determing the extent of the inner position
    %         in pixels. Optional, default: Current extent.
    %
    % OUTPUT:
    %   Out:  Output for the 'get' command.
    %
    % EXAMPLES:
    %   FigH = figure;
    %   LimitFigSize(FigH, 'min', [200, 200])
    %   LimitFigSize('max', [400, 400])        % Current figure as default
    %   LimitFigSize('min')                    % Current figure, current extent
    %   Limit = LimitFigSize(FigH, 'get')
    %
    % NOTE: This function is based on undocumented functions, such that it might
    % fail in future versions.
    %
    % Tested: Matlab 7.7, 7.8, 7.13, 8.6, WinXP/32, Win7/64
    %         Not compatible to Matlab 6.5.
    % Author: Jan Simon, Heidelberg, (C) 2012-2015 matlab.2010(a)n(MINUS)simon.de
    %
    % See also MaximizeFig, AlphaFig, ClipFig, TopFig, ToMonitorFig, WindowAPI.
    % $JRev: R-i V:008 Sum:6o/8OiL6gm/Z Date:20-Dec-2015 17:20:36 $
    % $License: BSD (use/copy/change/redistribute on own risk, mention the author) $
    % $UnitTest: uTest_LimitFigSize $
    % $File: Tools\GLGui\LimitFigSize.m $
    % History:
    % 001: 07-Oct-2012 15:02, First version.
    % 005: 02-Nov-2014 02:31, Consider HG2 client.
    %      BUGFIX: Replies empty matrix for 'get' if size has not been limited.
    % Initialize: ==================================================================
    % Global Interface: ------------------------------------------------------------
    % Consider dependencies to Matlab versions:
    MatlabVer = [100, 1] * sscanf(version, '%d.', 2);
    hasHG2    = (MatlabVer >= 804);
    hasHG1    = (MatlabVer >= 706);
    % hasHG2 = isa(FigH, 'matlab.ui.Figure');
    % hasHG2 = ~isnumeric(FigH);
    if ~usejava('jvm') || MatlabVer < 700
        % Should this be an error, a warning or a message??
        warning(['JSimon:', mfilename, ':NeedJava'], ...
            '*** %s: Require java and Matlab >= R14.', mfilename);
        return;
    end
    % Initial values: --------------------------------------------------------------
    % Program Interface: -----------------------------------------------------------
    % Parse the inputs:
    switch nargin
        case 0
            FigH = gcf;
            Prop = 'min';
            Ext  = [];
        case 1
            if numel(FigH) == 1 && ishandle(FigH)  % LimitFigSize(FigH)
                Prop = 'min';
                Ext  = [];
            elseif ischar(FigH)                    % LimitFigSize('string')
                Prop = FigH;
                Ext  = [];
                FigH = GetCurrentFigure(hasHG2);
            else                                   % LimitFigSize([W, H])
                Prop = 'min';
                Ext  = FigH;
                FigH = GetCurrentFigure(hasHG2);
            end
        case 2
            if ischar(FigH)                        % LimitFigSize('string', [W, H])
                Ext  = Prop;
                Prop = FigH;
                FigH = GetCurrentFigure(hasHG2);
            elseif numel(FigH) == 1 && ishandle(FigH)
                if ischar(Prop)                     % LimitFigSize(FigH, 'string')
                    Ext  = [];
                else                                % LimitFigSize(FigH, [W, H])
                    Ext  = Prop;
                    Prop = 'min';
                end
            else
                error(['JSimon:', mfilename, ':BadHandle'], ...
                    '*** %s: 1st input is not a figure handle or string.', mfilename);
            end
        case 3
            % Nothing to do
        otherwise
    end
    % Check if input is a figure handle:
    if numel(FigH) ~= 1 || ~ishandle(FigH) || ...     % Short-circuit
            ~strcmpi(get(FigH, 'Type'), 'figure')
        error(['JSimon:', mfilename, ':BadHandle'], ...
            '*** %s: 1st input is not a figure handle.', mfilename);
    end
    % User Interface: --------------------------------------------------------------
    % Do the work: =================================================================
    % Get the Java frame:
    if hasHG2
        % Suppress warnings:
        bakWarn = warning('off','MATLAB:HandleGraphics:ObsoletedProperty:JavaFrame');
        jFrame  = get(handle(FigH), 'JavaFrame');  % Fails in Matlab < 7.0
        warning(bakWarn);
    else
        jFrame = get(handle(FigH), 'JavaFrame');  % Fails in Matlab < 7.0
    end
    try
        if hasHG2      % HG2 since R2014b:
            jClient = jFrame.fHG2Client;
        elseif hasHG1  % R2008a up to 2014a:
            jClient = jFrame.fHG1Client;
        else           % R14 (7.0) to R2011a:
            jClient = jFrame.fFigureClient;
        end
    catch ME          % Getting the client is not documented and subject to changes:
        error(['JSimon:', mfilename, ':NoJClient'], ...
            'Cannot get Java Figure Client: %s', ME.message);
    end
    % The com.mathworks.hg.peer.FigureFrameProxy is empty, if no DRAWNOW allowed to
    % update the figure since its creation:
    jWindow = jClient.getWindow;
    if isempty(jWindow)
        drawnow;
        pause(0.2);
        jWindow = jClient.getWindow;
    end
    % Apply, clear or get the limits:
    if strncmpi(Prop, 'min', 3)
        Ext = GetExtent(FigH, Ext);
        jWindow.setMinimumSize(java.awt.Dimension(Ext(1), Ext(2)));
    elseif strncmpi(Prop, 'max', 3)
        Ext = GetExtent(FigH, Ext);
        jWindow.setMaximumSize(java.awt.Dimension(Ext(1), Ext(2)));
    elseif strcmpi(Prop, 'clear')
        jWindow.setMinimumSize([]);
        jWindow.setMaximumSize([]);
    elseif strcmpi(Prop, 'get')
        % Get difference to outer limits:
        bakUnits = get(FigH, 'Units');
        set(FigH, 'Units', 'pixels');
        DiffPos = get(FigH, 'OuterPosition') - get(FigH, 'Position');
        set(FigH, 'Units', bakUnits);
        jDim = jWindow.getMinimumSize;
        if isempty(jDim)  % Size has not been limited before:
            Out.MinSize = [];
        else
            Out.MinSize = [jDim.getWidth, jDim.getHeight] - DiffPos(3:4);
        end

        jDim = jWindow.getMaximumSize;
        if isempty(jDim)
            Out.MaxSize = [];
        else
            Out.MaxSize = [jDim.getWidth, jDim.getHeight] - DiffPos(3:4);
        end
    end
    % return;
    % ******************************************************************************
function FigH = GetCurrentFigure(hasHG2)
    % A DRAWNOW is required when the figure is created to update the Java handle.
    % See GCF.
    if hasHG2
        FigH = get(groot,'CurrentFigure');
    else
        FigH = get(0, 'CurrentFigure');
    end
    if isempty(FigH)
        FigH = figure;
        drawnow;
    end
    % return;
    % ******************************************************************************
function Ext = GetExtent(FigH, Ext)
    % Get or convert extent to outer position.
    bakUnits = get(FigH, 'Units');
    set(FigH, 'Units', 'pixels');
    OuterPos = get(FigH, 'OuterPosition');
    if isempty(Ext)
        Ext = OuterPos(3:4);
    elseif isnumeric(Ext) && numel(Ext) == 2 && all(Ext >= 0)
        InnerPos = get(FigH, 'Position');
        Ext      = [Ext(1) + OuterPos(3) - InnerPos(3), ...
            Ext(2) + OuterPos(4) - InnerPos(4)];
    else
        error(['JSimon:', mfilename, ':BadSize'], ...
            '*** %s: Extent must be a positive [1 x 2] vector.', mfilename);
    end
    set(FigH, 'Units', bakUnits);
    % return;


text:
function kSum=KroneckerSum(A,B)
%KRONECKERSUM Find the Kronecker sum of an nXn matrix A and an mXm matrix
%             B. This is defined as kron(A,eye(m))+kron(eye(n),B)
%
%INPUTS: A An nXn matrix.
%        B An mXm matrix.
%
%OUTPUTS: kSum The (n*m)X(n*m) Kronecker sum of A and B.
%
%The Kronecker sum arises in numerous matrix identities. For example, it is
%used  in [1] for expressions for eigenvectors.
%
%REFERENCES:
%[1] J. W. Brewer, "Kronecker products and matrix calculus in systems
%    theory," IEEE Transactions on Circuits and Systems, vol. CAS-25, no.
%    9, pp. 772-781, Sep. 1978.
%
%December 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n=size(A,1);
m=size(B,1);

kSum=kron(A,eye(m))+kron(eye(n),B);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function out = false( in )
%FALSE for initialization of nested cell arrays
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

out = cellfun(@(x) false(x), in.data{2,1}, 'UniformOutput', false);
out = reconFlatCellMatrix(out, in.data{1,1});
out = TRAFO(out,in.meta);

end

author: thomaskuestner


text:
function X_hat = prox_l1(X, tau)
    X_hat = max(abs(X) - tau, 0) .* sign(X);
end

author: YimianDai


text:
% dbm_get_hidden_raw
% Copyright (C) 2011 KyungHyun Cho, Tapani Raiko, Alexander Ilin
%
%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 should have received a copy of the GNU General Public License
%along with this program; if not, write to the Free Software
%Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
%
function [h_mf] = dbm_get_hidden_raw(x0, binary, layers, ...
    W, biases, sigmas, ...
    max_iter, tol, reg, ...
    do_centering, centers)

layers = layers;
n_layers = length(layers);

h_mf = cell(n_layers, 1);
h_mf{1} = x0;

for l = 2:n_layers
    if l == n_layers
        mult = 1;
    else
        mult = 2;
    end
    if l == 2
        if binary == 1
            h_mf{l} = sigmoid(bsxfun(@plus, mult * h_mf{l-1} * W{l-1}, biases{l}'));
        else
            h_mf{l} = sigmoid(bsxfun(@plus, mult * bsxfun(@rdivide, h_mf{l-1}, sigmas.^2') * W{l-1}, biases{l}'));
        end
    else
        h_mf{l} = sigmoid(bsxfun(@plus, mult * h_mf{l-1} * W{l-1}, biases{l}'));
    end
end

h_mf_prev = h_mf;

for iter = 1:max_iter
    diff_err = 0;
    for oddeven = [0 1]
        for l = 2:n_layers
            if mod(l, 2) == oddeven
                continue;
            end
            h_mf{l} = h_mf{l} * 0;
            if do_centering
                if l > 1
                    if l == 2
                        if binary == 1
                            h_mf{l} = h_mf{l} + bsxfun(@minus, h_mf{l-1}, centers{l-1}') * W{l-1};
                        else
                            h_mf{l} = h_mf{l} + bsxfun(@rdivide, h_mf{l-1}, sigmas.^2') * W{l-1};
                        end
                    else
                        h_mf{l} = h_mf{l} + bsxfun(@minus, h_mf{l-1}, centers{l-1}') * W{l-1};
                    end
                end

                if l < n_layers
                    h_mf{l} = h_mf{l} + bsxfun(@minus, h_mf{l+1}, centers{l+1}') * W{l}';
                end
            else
                if l > 1
                    if l == 2
                        if binary == 1
                            h_mf{l} = h_mf{l} + h_mf{l-1} * W{l-1};
                        else
                            h_mf{l} = h_mf{l} + bsxfun(@rdivide, h_mf{l-1}, sigmas.^2') * W{l-1};
                        end
                    else
                        h_mf{l} = h_mf{l} + h_mf{l-1} * W{l-1};
                    end
                end

                if l < n_layers
                    h_mf{l} = h_mf{l} + h_mf{l+1} * W{l}';
                end
            end

            h_mf{l} = sigmoid(bsxfun(@plus, h_mf{l}, biases{l}'));

            if reg > 0
                h_mf{l} = max(h_mf{l} - reg, 0);
            end

            diff_err = diff_err + sum(sum((h_mf_prev{l} - h_mf{l}).^2));
        end
    end

    %fprintf(2, '%d\n', diff_err);

    if diff_err < tol
        break;
    end

    h_mf_prev = h_mf;
end

clear h_mf_prev;

author: kyunghyuncho

repo: deepmat


text:
function cycle_floyd_test04 ( )

%*****************************************************************************80
%
%% CYCLE_FLOYD_TEST04 tests CYCLE_FLOYD for F4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CYCLE_FLOYD_TEST04\n' );
  fprintf ( 1, '  Test CYCLE_FLOYD for F4().\n' );
  fprintf ( 1, '  f4(i) = mod ( 31421 * i + 6927, 65536 ).\n' );

  x0 = 1;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Starting argument X0 = %d\n', x0 );

  [ lam, mu ] = cycle_floyd ( @f4, x0 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reported cycle length is %d\n', lam );
  fprintf ( 1, '  Expected value is 65536\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reported distance to first cycle element is %d\n', mu );
  fprintf ( 1, '  Expected value is 0\n' );

  return
end


text:
% tests for denoising using tv and wavelets

% 1D or 2D test
test = 1;

%% load data
if test==1
    name = 'piece-polynomial';
    n = 1024;
    x0 = load_signal(name, n);
    sigma = .03;
else
    name = 'barb';
    n = 256;
    x0 = load_image(name);
    x0 = crop(x0,n);
    sigma = .05;
end

eta = .05;
x0 = rescale(x0,eta,1-eta);

randn('state',123456);
x = x0 + sigma*randn(size(x0));

%% wavelet domain denoising
%options for wavelets
options.ti = 1;
if test==1
    options.wavelet_vm = 4;
    options.wavelet_type = 'haar';
else
    options.wavelet_vm = 3;
    options.wavelet_type = 'biorthogonal';
    options.decomp_type = 'quad';
end
Jmin = 3;

niter = 10;
Tlist = sigma*linspace(.5,2.5,niter);
err = [];
if test==1
    xw = perform_wavelet_transform(x, Jmin, +1, options);
else
    xw = perform_atrou_transform(x, Jmin, options);
end
for i=1:niter
    xwT = perform_thresholding(xw,Tlist(i), 'soft');
    if test==1
        x1 = perform_wavelet_transform(xwT, Jmin, -1, options);
    else
        x1 = perform_atrou_transform(xwT, Jmin, options);
    end
    err(i) = norm(x0-x1,'fro');
    if i>1 && err(i)<min(err(1:i-1))
        xwav = x1;
    elseif i==1
        xwav = x1;
    end
end

%% TV denoising
options.niter = 5000;

options.etgt = norm(x-x0);
options.etgt = [];

if test==1
    options.lambda_max = .4;
    options.lambda_min = .1;
else
    options.lambda_max = .1;
    options.lambda_min = .0;
end
options.x0 = x0;
[xtv,err] = perform_tv_denoising(x,options);

pwav = psnr(x0,xwav,1);
ptv = psnr(x0,xtv,1);
pnoisy = psnr(x0,x1,1);
disp(['Wav=' num2str(pwav) 'dB, TV=' num2str(ptv) 'dB.' ]);

rep = 'results/tv-denoising/';
if not(exist(rep))
    mkdir(rep);
end

% save in a file the PSNR
filename = [rep 'results.txt'];
fid = fopen(filename, 'a');
fprintf( fid, '%s - noisy=%fdB - wav(%s)=%fdB - tv=%fdB\n', name, pnoisy, options.wavelet_type, pwav, ptv );
fclose(fid);

% display
if test==1
    fs = 20;
    clf;
    subplot(3,1,1);
    plot(x, 'k'); axis([1 n 0 1]);
    set(gca, 'FontSize', fs);
    subplot(3,1,2);
    plot(xwav, 'k'); axis([1 n 0 1]);
    set(gca, 'FontSize', fs);
    subplot(3,1,3);
    plot(xtv, 'k'); axis([1 n 0 1]);
    set(gca, 'FontSize', fs);
    saveas(gcf, [rep name '-tv-denoising.png'], 'png');
    saveas(gcf, [rep name '-tv-denoising.eps'], 'eps');
else
    imageplot({x0 x xwav xtv}, ...
            {'Original' ['Noisy=' num2str(pnoisy)] ['Wav=' num2str(pwav)]  ['TV=' num2str(ptv)] }, 2,2);
    saveas(gcf, [rep name '-tv-denoising.png'], 'png');
    warning off;
    imwrite(clamp(x0), [rep name '-original.png'], 'png' );
    imwrite(clamp(x), [rep name '-noisy.png'], 'png' );
    imwrite(clamp(xwav), [rep name '-wav.png'], 'png' );
    imwrite(clamp(xtv), [rep name '-tv.png'], 'png' );
    warning on;
end

author: gpeyre

repo: matlab-toolboxes


text:
% l1_softth - soft threshold function for L1 regularization
%
% Copyright(c) 2009 Ryota Tomioka
% This software is distributed under the MIT license. See license.txt

function [vv,ss]=l1_softth(vv,lambda,info)

n = size(vv,1);

Ip=find(vv>lambda);
In=find(vv<-lambda);

vv=sparse([Ip;In],1,[vv(Ip)-lambda;vv(In)+lambda],n,1);

ss=abs(vv);

repo: REST


text:
function [d] = lbfgsProd(g,S,Y,YS,lbfgs_start,lbfgs_end,Hdiag)
% BFGS Search Direction
%
% This function returns the (L-BFGS) approximate inverse Hessian,
% multiplied by the negative gradient

% Set up indexing
[nVars,maxCorrections] = size(S);
if lbfgs_start == 1
	ind = 1:lbfgs_end;
	nCor = lbfgs_end-lbfgs_start+1;
else
	ind = [lbfgs_start:maxCorrections 1:lbfgs_end];
	nCor = maxCorrections;
end
al = zeros(nCor,1);
be = zeros(nCor,1);

d = -g;
for j = 1:length(ind)
	i = ind(end-j+1);
	al(i) = (S(:,i)'*d)/YS(i);
%     al(i) = sum(S(:, i) .*d(:)) / YS(i);
	d = d-al(i)*Y(:,i);
end

% Multiply by Initial Hessian
d = Hdiag*d;

for i = ind
	be(i) = (Y(:,i)'*d)/YS(i);
%     be(i) = sum(Y(:, i) .* d(:)) / YS(i);
	d = d + S(:,i)*(al(i)-be(i));
%     if i == 1
%         fprintf("albe %.15f %.15f %.15f %.15f\n", al(i)-be(i), al(i), be(i), S(i));
%     end
end

repo: pTVreg


text:
fileID = fopen('rgd_cropped320/rgb_0000_annotationsCropped320.txt','r');
sizeA = [2 inf];
A = fscanf(fileID, '%f %f', sizeA)
imread('rgd_cropped320/rgd_0000Cropped320.png');

figure;
imshow(ans);
hold on;
[row col] = size(A);

for idx = 1:4:col
    x = A(1, [idx:idx+3 idx]);
    y = A(2, [idx:idx+3 idx]);
    plot(x,y);
    hold on
end
fclose(fileID);

author: ivalab


text:
%% Anaylsis_MI
clear all; clc; close all;
%% initialization
DATADIR = 'WHERE\IS\DATA';
%% MI
MIDATA = 'EEG_MI.mat';
STRUCTINFO = {'EEG_MI_train', 'EEG_MI_test'};
SESSIONS = {'session1', 'session2'};
TOTAL_SUBJECTS = 54;
FS=100;

%% PERFORMANCE PARAMETERS
params = { 'task',{'mi_off','mi_on'}; ...
    'channel_index', [8:11 13:15 18:21 33:41]; ...
    'band', [8 30]; ...
    'time_interval', [1000 3500]; ...
    'CSPFilter', 2; ...
    };

% for MI_cv
Niteration = 10;
% for CSSP
tau=[0.01:0.01:0.15]*1000;
% for FBCSP
filterbank = [4 8;8 12;12 16;16 20;20 24;24 28;28 32;32 36;36 40];
NUMfeat = 4;
% for BSSFO
bssfo_param = {'init_band', [4 40]; ...
    'numBands', 30; ...
    'numIteration', 10; ...
    'mu_band', [7 15]; ...
    'beta_band', [14 30]; ...
    };
%% validation
for sessNum = 1:length(SESSIONS)
    session = SESSIONS{sessNum};
    fprintf('\n%s validation\n',session);
    for subNum = 1:TOTAL_SUBJECTS
        subject = sprintf('s%d',subNum);
        fprintf('LOAD %s ...\n',subject);

        data = importdata(fullfile(DATADIR,session,subject,MIDATA));

        CNT{1} = prep_resample(data.(STRUCTINFO{1}), FS,{'Nr', 0});
        CNT{2} = prep_resample(data.(STRUCTINFO{2}), FS,{'Nr', 0});

        ACC.MI_cv(subNum,sessNum) = mi_performance(CNT,params,Niteration);
        ACC.MI_off2on(subNum,sessNum) = mi_performance_off2on(CNT,params);
        ACC.MI_CSSP(subNum,sessNum) = cssp_off2on(CNT,params,tau);
        ACC.MI_FBCSP(subNum,sessNum) = fbcsp_off2on(CNT,params,filterbank,NUMfeat);
        ACC.MI_BSSFO(subNum,sessNum) = bssfo_off2on(CNT,params,bssfo_param);

        fprintf('CSP_crsval = %f\n',ACC.MI_cv(subNum,sessNum));
        fprintf('CSP = %f\n',ACC.MI_off2on(subNum,sessNum));
        fprintf('CSSP = %f\n',ACC.MI_CSSP(subNum,sessNum));
        fprintf('FBCSP = %f\n',ACC.MI_FBCSP(subNum,sessNum));
        fprintf('BSSFO = %f\n',ACC.MI_BSSFO(subNum,sessNum));
        clear CNT
    end
end

author: PatternRecognition

repo: OpenBMI


text:
function [u,s,U_s,U_pk,U_pd,U_idp]  = projIdpPntIntoPinHole(Sf, Spk, Spd, idp)

% PROJIDPPNTINTOPINHOLE Project Idp pnt into pinhole.
%    [U,S] = PROJIDPPNTINTOPINHOLE(RF, SF, SPK, SPD, L) projects 3D Inverse
%    Depth points into a pin-hole camera, providing also the non-measurable
%    depth. The input parameters are:
%       SF : pin-hole sensor frame
%       SPK: pin-hole intrinsic parameters [u0 v0 au av]'
%       SPD: radial distortion parameters [K2 K4 K6 ...]'
%       L  : 3D inverse depth point [x y z pitch yaw rho]'
%    The output parameters are:
%       U  : 2D pixel [u v]'
%       S  : non-measurable depth
%
%    The function accepts an idp points matrix L = [L1 ... Ln] as input.
%    In this case, it returns a pixels matrix U = [U1 ... Un] and a depths
%    row-vector S = [S1 ... Sn].
%
%    [U,S,U_S,U_K,U_D,U_L] = ... gives also the jacobians of the
%    observation U wrt all input parameters. Note that this only works for
%    single points.
%
%    See also PINHOLE, PY2VEC, PROJIDPPNTINTOPINHOLEONROB.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

p0 = idp(1:3,:);   % origin
py = idp(4:5,:);   % pitch and roll
r  = idp(6,:);     % inverse depth

t = Sf.t; % frame things
q = Sf.q;

if nargout <= 2  % No Jacobians requested

    m = py2vec(py);    % vector from anchor
    v = m - (t-p0).*r;  % vector from sensor
    w = Rtp(q,v);
    u = pinHole(w, Spk, Spd); % pixel

    e = w./r;  % euclidean
    s = e(3,:);

else            % Jacobians requested

    if size(idp,2) == 1

        % function calls
        [m,M_py] = py2vec(py);    % vector from anchor
        v        = m - (t-p0)*r;  % vector from sensor
        V_m      = 1;
        V_p0     = r;  %  r*eye(3);
        V_t      = -r; % -r*eye(3);
        V_r      = p0-t;
        [w, W_q, W_v] = Rtp(q,v);
        [u, ~, U_w, U_pk, U_pd] = pinHole(w, Spk, Spd); % pixel

        e = w/r;   % euclidean
        s = e(3);  % depth of point

        % chain rule
        U_v  = U_w*W_v;
        U_p0 = U_v*V_p0;
        U_py = U_v*V_m*M_py;
        U_r  = U_v*V_r;

        U_t  = U_v*V_t;
        U_q  = U_w*W_q;

        U_idp = [U_p0 U_py U_r];
        U_s   = [U_t U_q];

    else
        error('??? Jacobians not available for multiple IDP points.')

    end

end

return

%% jac
syms x y z a b c d real
syms au av u0 v0 real
syms d1 d2 d3 real
syms x0 y0 z0 p y r real

Sf.x = [x;y;z;a;b;c;d];
Sf   = updateFrame(Sf);
Spk  = [u0;v0;au;av];
Spd  = [d1;d2;d3];
idp  = [x0;y0;z0;p;y;r];

[u,s,U_s,U_pk,U_pd,U_idp]  = projIdpPntIntoPinHole(Sf, Spk, Spd, idp);
u,s  = projIdpPntIntoPinHole(Sf, Spk, Spd, idp);

simplify(U_pk  - jacobian(u,Spk))
simplify(U_pd  - jacobian(u,Spd))
% simplify(U_idp - jacobian(u,idp))

% ========== 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:
%% housekeeping
close all
clc
clear

%% add a path

addpath('routines')

%% rise the model

m=rise('lwz_switch_model','solve_linear',true);

%% get the parameters

[p,priors]=create_parameters_switch(true);

m=set(m,'parameters',p);

%% bring in the data

data=create_data();

data=pages2struct(data);

vnames=fieldnames(data);

figure('name','Observables')

for ii=1:numel(vnames)

    v=vnames{ii};

    subplot(3,2,ii)

    plot(data.(v))

    title(v)

end

%% Normalization

myNormalization={
    'sig_Eps_phi_hetero_1<=sig_Eps_phi_hetero_2'
    };

%% estimate the model
clc

%do_estimate=false;
%
%if do_estimate

    [ms,filtration]=estimate(m,'data',data,...
        'estim_priors',priors,...
        'kf_presample',3,...
        'kf_init_variance',10,...
        'estim_nonlinear_restrictions',myNormalization);

    % save the mode
    %--------------

    switch_mode=get(ms,'mode');

%    save switch_estimates ms switch_mode filtration
%
%else
%
%    load switch_estimates
%
%end

%% irfs

clc

myirfs=irf(ms);

quick_irfs(ms,myirfs)

%% probability plots

plot_probabilities(ms)

%% probability plots

plot_data_against_probabilities(ms)

%% replicate figure 9
close all

probs=filtration.smoothed_state_probabilities;
% I substract a constant since I don't know the correct level of the data.
% Maybe Dan can help with this?
lp=cumsum(data.DLogQl)-0.55;

plotyy(lp,probs.hetero_2,'linewidth',2)

%% (Approximate) Historical decomposition

h=historical_decomposition(ms);

obslist=get(ms,'obs_list');

dropped={'init'};

figure('name','Observables')

for ii=1:numel(obslist)

    v=obslist{ii};

    subplot(3,2,ii)

    plot_decomp(drop(h.(v),dropped))

    title(v)

    if ii==1

        l=legend(h.(v).varnames-dropped);

    end

end

set(l,'interpreter','none')

%% forecast

%% conditional forecast

%% counterfactual

%% posterior sampling
clc

do_sampling=false;

if do_sampling

    [objective,lb,ub,mu,SIG]=pull_objective(ms,...
        'solve_check_stability',false,'fix_point_TolFun',1e-6);

    sampling_options=struct('MaxTime',3600*24*7,...
        'N',10000,...
        'nchain',2,...
        'thin',5,...
        'adapt_covariance',true);

    results=mh_sampler(objective,lb,ub,sampling_options,mu);

    save switch_sampling results objective

end

do_mdd=true;

if do_mdd

    load switch_sampling

    % Marginal data density
    algorithms={'mhm','swz','mueller','bridge','is','ris','cj','laplace','laplace_mcmc'};
    options=struct('log_post_kern',objective,'L',10000);
    nalgos=numel(algorithms);
    log_mdd=cell(nalgos,3);
    log_mdd(:,1)=algorithms(:);
    tictoc=nan(nalgos,1);

    pop=results{1}.pop;
    for jj=2:numel(results)
        pop=[pop;results{jj}.pop];
    end

    parfor ii=1:nalgos
        opt=options;
        opt.algorithm=algorithms{ii};
        tic
        log_mdd{ii,2} = mcmc_mdd(pop,lb,ub,opt);
        tictoc(ii)=toc;
    end
    log_mdd(:,3)=num2cell(tictoc(:));

    save switch_sampling log_mdd -append

end

% the benchmark DSGE model (2344.0), the DSGE with two volatility regimes
% (2354.1), and the DSGE model with three volatility regimes (2353.2)
%%

pnames={ms.estimation.priors.name};

obj=mcmc(results,pnames)%,start_from,trimming

%% Trace plots

chain_id=1;

plotfunc=@(name)traceplot(obj,name,chain_id);

hfig=utils.plot.multiple(plotfunc,pnames,'trace plots',3,3)

%% Density plots

chain_id=1;

plotfunc=@(name)densplot(obj,name,chain_id,250);

hfig=utils.plot.multiple(plotfunc,pnames,'density plots',3,3)

%% Mean plots

chain_id=[];

plotfunc=@(name)meanplot(obj,name,chain_id);

hfig=utils.plot.multiple(plotfunc,pnames,'Mean plots',3,3);

%% Autocorrelation plots

chain_id=[];

order=10;

plotfunc=@(name)autocorrplot(obj,name,chain_id,order);

hfig=utils.plot.multiple(plotfunc,pnames,'Autocorrelation plots',3,3)

%% Gelman-Rubin plots

plotfunc=@(name)psrf_plot(obj,name);

hfig=utils.plot.multiple(plotfunc,pnames,'Gelman-Rubin (PSRF) plots',3,3)

%% Scatter plots

np=numel(pnames);

scatnames=cell(1,np^2);

iter=0;
for ii=1:np
    for jj=ii+1:np
        iter=iter+1;
        scatnames{iter}={pnames{ii},pnames{jj}};
    end
end
scatnames=scatnames(1:iter);

chain_id=[];

plotfunc=@(name)scatterplot(obj,name{1},name{2},chain_id);%,varargin

hfig=utils.plot.multiple(plotfunc,scatnames,'Scatter plots',7,7)

%% plot priors

[pdata,hdl]=plot_priors(ms,pnames)

%% plot posteriors

plot_posteriors(ms,results.pop,pnames)

%% priors and posteriors

plot_priors_and_posteriors(ms,results.pop,pnames)

%% mode curvature

mode_curvature(ms)

author: jmaih