text:
function [y,L] = mdwt(x,h,L);
%    [y,L] = mdwt(x,h,L);
%
%    Function computes the discrete wavelet transform y for a 1D or 2D input
%    signal x using the scaling filter h.
%
%    Input:
%	x : finite length 1D or 2D signal (implicitly periodized)
%       h : scaling filter
%       L : number of levels. In the case of a 1D signal, length(x) must be
%           divisible by 2^L; in the case of a 2D signal, the row and the
%           column dimension must be divisible by 2^L. If no argument is
%           specified, a full DWT is returned for maximal possible L.
%
%    Output:
%       y : the wavelet transform of the signal
%           (see example to understand the coefficients)
%       L : number of decomposition levels
%
%    1D Example:
%       x = makesig('LinChirp',8);
%       h = daubcqf(4,'min');
%       L = 2;
%       [y,L] = mdwt(x,h,L)
%
%    1D Example's  output and explanation:
%
%       y = [1.1097 0.8767 0.8204 -0.5201 -0.0339 0.1001 0.2201 -0.1401]
%       L = 2
%
%    The coefficients in output y are arranged as follows
%
%       y(1) and y(2) : Scaling coefficients (lowest frequency)
%       y(3) and y(4) : Band pass wavelet coefficients
%       y(5) to y(8)  : Finest scale wavelet coefficients (highest frequency)
%
%    2D Example:
%
%       load test_image
%       h = daubcqf(4,'min');
%       L = 1;
%       [y,L] = mdwt(test_image,h,L);
%
%    2D Example's  output and explanation:
%
%       The coefficients in y are arranged as follows.
%
%              .------------------.
%              |         |        |
%              |    4    |   2    |
%              |         |        |
%              |   L,L   |   H,L  |
%              |         |        |
%              --------------------
%              |         |        |
%              |    3    |   1    |
%              |         |        |
%              |   L,H   |  H,H   |
%              |         |        |
%              `------------------'
%
%       where
%            1 : High pass vertically and high pass horizontally
%            2 : Low pass vertically and high pass horizontally
%            3 : High pass vertically and low  pass horizontally
%            4 : Low pass vertically and Low pass horizontally
%                (scaling coefficients)
%
%
%
%
%    See also: midwt, mrdwt, mirdwt
%

%File Name: mdwt.m
%Last Modification Date: 08/07/95	15:13:25
%Current Version: mdwt.m	2.4
%File Creation Date: Wed Oct 19 10:51:58 1994
%Author: Markus Lang  <lang@jazz.rice.edu>
%
%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
%Created by Markus Lang, Department of ECE, Rice University.
%
%This software is distributed and licensed to you on a non-exclusive
%basis, free-of-charge. Redistribution and use in source and binary forms,
%with or without modification, are permitted provided that the following
%conditions are met:
%
%1. Redistribution of source code must retain the above copyright notice,
%   this list of conditions and the following disclaimer.
%2. Redistribution in binary form must reproduce the above copyright notice,
%   this list of conditions and the following disclaimer in the
%   documentation and/or other materials provided with the distribution.
%3. All advertising materials mentioning features or use of this software
%   must display the following acknowledgment: This product includes
%   software developed by Rice University, Houston, Texas and its contributors.
%4. Neither the name of the University nor the names of its contributors
%   may be used to endorse or promote products derived from this software
%   without specific prior written permission.
%
%THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
%AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
%BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
%FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
%OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
%EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
%PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
%OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
%WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
%OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
%USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
%For information on commercial licenses, contact Rice University's Office of
%Technology Transfer at techtran@rice.edu or (713) 348-6173
%
%Change History:
%
%Modification #1
%Mon Aug  7 11:42:11 CDT 1995
%Rebecca Hindman <hindman@ece.rice.edu>
%Added L to function line so that it can be displayed as an output
%
%Change History:
%
%Modification #1
%Thu Mar  2 13:07:11 CDT 2000
%Ramesh Neelamani<neelsh@ece.rice.edu>
%Revamped the help file
%

author: thomaskuestner


text:
function [ mX ] = ProjectPosSemiDefinite( mY, numIterations, stopTol )
% ----------------------------------------------------------------------------------------------- %
% [ mX ] = ProjectPosSemiDefinite( mY, numIterations, stopTol )
%   Solves \arg \min_{X} 0.5 || X - Y ||, s.t. X \in \mathcal{S}_{+}^{n} by
%   Iterative Projection onto the set of Symmetric Matrices and PSD
%   Matrices.
% Input:
%   - mY            -   Input Matrix.
%                       Input model matrix.
%                       Structure: Matrix (m x n).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - numIterations -   Number of Iterations.
%                       Sets the number of iterations for the algorithm to
%                       run.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range {1, 2, ...}.
%   - stopTol       -   Stopping Condition Tolerance.
%                       Sets the stopping threshold for the L Inf (Maximum
%                       Absolute Value) of the change between 2 iterations
%                       of the algorithm.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range [0, inf).
% Output:
%   - mX            -   Solution Matrix.
%                       A PSD Symmetric matrix which is the closest in
%                       Frobenius Norm sense to 'mY'.
%                       Structure: Matrix (m x n).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% References
%   1.  A
% Remarks:
%   1.  Since both the Symmetric Matrices Set and the PSD Matrices Set are
%       sub spaces the iterative (POCS) projection onto each is the
%       orthogonal projection onto the intersection. See https://math.stackexchange.com/questions/1492095.
% TODO:
%   1.  C
% Release Notes:
%   -   1.0.000     12/04/2020  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

vectorMode = OFF;

if(size(mY, 2) == 1)
    numRows     = sqrt(size(mY, 1));
    vectorMode  = ON;
    mX          = reshape(mY, numRows, numRows);
else
    numRows = size(mY, 1);
    mX      = mY;
end

mQ      = zeros(numRows, numRows);
mD      = zeros(numRows, numRows);
mXPrev  = mX;

for ii = 1:numIterations
    % Projection onto Symmetric Matrices Set
    mX(:) = (mX.' + mX) / 2;

    % Projection onto PSD Matrices Set
    [mQ(:), mD(:)] = eig(mX);
    mD(:) = max(mD, 0);
    mX(:) = mQ * mD * mQ.';

    if(max(abs(mX(:) - mXPrev(:))) <= stopTol)
        break;
    end

    mXPrev(:) = mX;

end

if(vectorMode == ON)
    mX = mX(:);
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function b = sprepmat(a,m,n)
% function B = sprepmat(A,M,N)
%
% DESCRIPTION
%   B is a sparse matrix consisting of an MxN tiling of the
%   sparse matrix A.
%   [repmat fails on scalar, sparse matrices in Matlab 6.5:
%     >>repmat(sparse(5),2,3)
%    This yields an error.]
%
% INPUTS
%   A: sparse matrix
%   M,N: Copies of A in the row and column directions
%
% OUTPUTS
%   B: sparse matrix
%
% SYNTAX
%   B = sprepmat(A,M,N);

% This file is part of SOSTOOLS - Sum of Squares Toolbox ver 3.00.
%
% Copyright (C)2002, 2004, 2013  A. Papachristodoulou (1), J. Anderson (1),
%                                G. Valmorbida (1), S. Prajna (2),
%                                P. Seiler (3), P. A. Parrilo (4)
% (1) Department of Engineering Science, University of Oxford, Oxford, U.K.
% (2) Control and Dynamical Systems - California Institute of Technology,
%     Pasadena, CA 91125, USA.
% (3) Aerospace and Engineering Mechanics Department, University of
%     Minnesota, Minneapolis, MN 55455-0153, USA.
% (4) Laboratory for Information and Decision Systems, M.I.T.,
%     Massachusetts, MA 02139-4307
%
% Send bug reports and feedback to: sostools@cds.caltech.edu
%
% 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, see <http://www.gnu.org/licenses/>.

% 1/30/2003: PJS  Initial Coding

[nra,nca] = size(a);
ridx = (1:nra)';
ridx = ridx(:,ones(1,m));
cidx = (1:nca)';
cidx = cidx(:,ones(1,n));
b = a(ridx,cidx);

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [movie_frame] = show_s2_partition(N,varargin)
%SHOW_S2_PARTITION 3D illustration of an EQ partition of S^2
%
%Syntax
% [movie_frame] = show_s2_partition(N,options);
%
%Description
% SHOW_S2_PARTITION(N) uses a 3d plot to illustrate the partition of
% the unit sphere S^2 into N regions.
%
% MOVIE_FRAME = SHOW_S2_PARTITION(N) sets MOVIE_FRAME to be an array of
% movie frames for use with MOVIE. The movie frames will contain the region by
% region build-up of the illustration.
%
% SHOW_S2_PARTITION(N,'offset','extra') uses experimental extra offsets.
% For more detail on partition options, see HELP PARTITION_OPTIONS.
%
% SHOW_S2_PARTITION(N,options) also recognizes a number of illustration
% options, which are specified as name, value pairs.
% Any number of pairs can be used, in any order.
%
% The following illustration options are used.
%
% SHOW_S2_PARTITION(N,'fontsize',size)
% Font size used in titles (numeric, default 16).
%
% SHOW_S2_PARTITION(N,'title','show')
% SHOW_S2_PARTITION(N,'title','hide')
% Show or hide title (default 'show').
%
% SHOW_S2_PARTITION(N,'points','show')
% SHOW_S2_PARTITION(N,'points','hide')
% Show or hide center points (default 'show').
%
% SHOW_S2_PARTITION(N,'sphere','show')
% SHOW_S2_PARTITION(N,'sphere','hide')
% Show or hide the unit sphere S^2 (default 'show').
%
% For more detail on illustration options, see HELP ILLUSTRATION_OPTIONS.
%
%Examples
% > show_s2_partition(10)
% > frames=show_s2_partition(9,'offset','extra')
% frames =
% 1x10 struct array with fields:
%     cdata
%     colormap
% > show_s2_partition(99,'points','hide')
%
%See also
% MOVIE, PARTITION_OPTIONS, ILLUSTRATION_OPTIONS, PROJECT_S2_PARTITION

% Copyright 2004-2005 Paul Leopardi for the University of New South Wales.
% $Revision 1.10 $ $Date 2005-06-01 $
% Function changed name from s2x to polar2cart
% Documentation files renamed
% $Revision 1.00 $ $Date 2005-02-13 $
%
% For licensing, see COPYING.
% For references, see AUTHORS.
% For revision history, see CHANGELOG.

pdefault.extra_offset =  false;
popt = partition_options(pdefault, varargin{:});

gdefault.fontsize = 16;
gdefault.show_title  = true;
gdefault.show_points = true;
gdefault.show_sphere = true;
gopt = illustration_options(gdefault, varargin{:});

dim = 2;

surf_jet;

if gopt.show_title
    if gopt.show_points
        pointstr = ', showing the center point of each region';
    else
        pointstr = '';
    end
    titlestr = sprintf(...
        '\nRecursive zonal equal area partition of {S^2} \n into %d regions%s.',...
        N,pointstr);
    title(titlestr,'FontWeight','bold','FontUnits','normalized',...
        'FontSize',gopt.fontsize/512);
end

frame_no = 1;
if nargout > 0
    movie_frame(frame_no) = getframe(gcf);
    frame_no = frame_no + 1;
end

if gopt.show_sphere
    show_s2_sphere;
    hold on
    if nargout > 0
        movie_frame(frame_no) = getframe(gcf);
        frame_no = frame_no + 1;
    end
end

R = eq_regions(dim,N,popt.extra_offset);
top_colat = 0;
for i = N:-1:2
    if top_colat ~= R(2,1,i)
        top_colat = R(2,1,i);
        pause(0);
    end
    show_s2_region(R(:,:,i),N);
    if nargout > 0
        movie_frame(frame_no) = getframe(gcf);
        frame_no = frame_no + 1;
    end
end

if gopt.show_points
    x = eq_point_set(dim,N,popt.extra_offset);
    show_r3_point_set(x,'sphere','hide','title','hide');
    hold on
    if nargout > 0
        movie_frame(frame_no) = getframe(gcf);
        frame_no = frame_no + 1;
    end
end

hold off
%
% end function

function show_s2_region(region,N)
%SHOW_S2_REGION Illustrate a region of S^2
%
%Syntax
% show_s2_region(region,N);
%
%Description
% SHOW_S2_REGION(REGION,N) uses 3D surface plots to illustrate a region of S^2.
% The region is given as a 2 x 2 matrix in spherical polar coordinates

tol = eps*2^5;

dim = size(region,1);
t = region(:,1);
b = region(:,2);

if abs(b(1)) < tol
    b(1) = 2*pi;
end
pseudo = 0;
if abs(t(1)) < tol && abs(b(1)-2*pi) < tol
    pseudo = 1;
end
n = 21;
delta = 1/(n-1);
h = 0:delta:1;
t_to_b = zeros(dim,n);
b_to_t = t_to_b;
r = sqrt(1/N)/12;
for k = 1:dim
    if ~pseudo || k < 2
        L = 1:dim;
        j(L) = mod(k+L,dim)+1;
        t_to_b(j(1),:) = t(j(1))+(b(j(1))-t(j(1)))*h;
        t_to_b(j(2),:) = t(j(2))*ones(1,n);
        t_to_b_x = polar2cart(t_to_b);
        [X,Y,Z] = fatcurve(t_to_b_x,r);
        surface(X,Y,Z,-ones(size(Z)),...
       'FaceColor','interp','FaceLighting','phong','EdgeColor','none')
        axis equal
        hold on
    end
end
grid off
axis off
%
% end function

repo: gspbox


text:
function c=mapcolor(d,map,clim)
% Usage: c=mapcolor(d,[map],[clim])
% Given an array "d", mapcolor uses the colormap "map" to assign a color to
% each entry in "d", returning the results as RGB triplets in "c". Provide
% a two-element vector "clim" to impose limits on the colormap. Basically
% this does exactly what Matlab's built-in colormaps do, but does not
% affect the "Colormap" property of the current figure--so you can use
% multiple colormaps on the same figure. See also mapcolorbar.m.
%
%   Example:
%      load clown
%      mylim=[min(X(:)) max(X(:))];
%      Xhot=mapcolor(X,hot,mylim);
%      subplot(211);
%      imagesc(Xhot);
%      mapcolorbar(hot,mylim);
%      Xjet=mapcolor(X,jet,mylim);
%      subplot(212);
%      imagesc(Xjet);
%      mapcolorbar(jet,mylim);

% Written 29 April 2010 by Douglas H. Kelley, dhk [at] dougandneely.com.

mapdefault=jet(256); % specifying a size prevents a figure from popping up

if nargin<1
    error(['Usage: c = ' mfilename '(d,[map],[clim])'])
end
if ~exist('map','var') || isempty(map)
    map=mapdefault;
elseif ischar(map)
    warning(['MATLAB:' mfilename ':mapIsString'], ...
        ['Map given as a string; attempting to evaluate ''' map '''.'])
    map=eval(map);
end
if ~exist('clim','var') || isempty(clim)
    clim=[min(d(:)) max(d(:))];
end

N=numel(d);
NN=size(d);
mapsize=size(map,1);
binnum=floor( mapsize/diff(clim)*(d-clim(1))+1 ); % A linear mapping...
binnum(isnan(binnum))=1; % ...setting NaNs to the mininum of the map, ...
binnum(binnum<1)=1; % ...then correcting for saturation at bottom...
binnum(binnum>mapsize)=mapsize; % ...and top.
if isvector(d)
    c=nan*ones(N,3);
    for ii=1:N
        c(ii,:)=map(binnum(ii),:);
    end
elseif ndims(d)==2
    c=nan*ones([size(d) 3]);
    [i,j]=ind2sub(NN,[1:N]);
    for ii=1:N
        c(i(ii),j(ii),:)=map(binnum(ii),:);
    end
else
    error(['Sorry, mapcolor does not support arrays of dimension '...
        'greater than two.'])
end

author: Sable

repo: mcbench-benchmarks


text:
function compDDSampling( aid )
%COMPDDSAMPLING Summary of this function goes here
%   Detailed explanation goes here
global config;
bufname = config.bufname;

% load([bufname config.objHypoFile3D],'CAN_POOL');
A = load([bufname config.roomModelFile],'orientation','surfacelabel_ori','colorName_rot');
omap = A.orientation;
gc = A.surfacelabel_ori;
cn = A.colorName_rot;
clear A;

% bufName = sprintf([bufname '/' config.objTestImgFeature], aid);
bufName = [bufname '/' config.objTestImgFeature];
% gndBufName = sprintf([bufname '/' config.objGndImgFeature], aid);
gndBufName = [bufname '/' config.objGndImgFeature];
load([bufname config.vpEstimationFile]);
load(config.annotationfile);
H2G_R = ANNO_ALL(aid).ANNO3D.R * ANNO_ALL(aid).ANNO3D.Rc / R;
clear ANNO_ALL;

load(config.groundtruthfile);
gnd = ALL_GNDS_ROT{aid,1};
GND_COMP = getHypInfo({gnd});
clear ALL_GNDS_ROT
load(config.roomtransformfile);
VALID_TRANS_GNDS = TRANS_GNDS(config.anno_valid);
clear TRANS_GNDS

load([bufname config.objHypoFile3D]);
ALL_POOL = CAN_POOL;
clear CAN_POOL

rot_file_name = [bufname config.IMGLIST(aid).name config.rotateImageSuffix];
rotImg = imread(rot_file_name);

%% sample
for rid = 1:length(ALL_POOL)
    fprintf('Room ID: %d\n', rid);
    if config.USEPREVIOUS
        if exist(sprintf([bufname '/' config.sceneFeatFile], rid), 'file') ...
                && exist(sprintf([bufname '/' config.objFeatFile], rid), 'file')
            continue;
        end
    end

    CAN_POOL = ALL_POOL(rid);
    if isempty(CAN_POOL{1})
        objfea = [];
        roomfea = [];
        anglesid = [];
        ALL_HYPS = [];
        if config.UPDATE_TO_DISK
            save(sprintf([bufname '/' config.objFeatFile], rid),'objfea','roomfea','anglesid');
            save(sprintf([bufname '/' config.sceneFeatFile], rid), 'ALL_HYPS');
        end
        fprintf('Empty pool, done.\n');
        continue;
    end

    [objfea, roomfea, anglesid] = compObjHypsFeature( rotImg, omap, gc, cn, CAN_POOL{1}, config.typenum, bufName );
    if config.UPDATE_TO_DISK
        save(sprintf([bufname '/' config.objFeatFile], rid),'objfea','roomfea','anglesid');
    end

    CAN_POOL{1}.sel_hyps.objfea = objfea;
    CAN_POOL{1}.room.roomfea = roomfea;
    CAN_POOL{1}.sel_hyps.anglesid = anglesid;
    ALL_HYPS = globalSampling( CAN_POOL{1}.sel_hyps, config.samplenum, config.typenum );
    [ sceneImgFea ] = compSceneHypsFeatureA( ALL_HYPS, CAN_POOL{1}.sel_hyps, CAN_POOL{1}.room, omap, gc, cn );
    for i = 1:config.samplenum
        ALL_HYPS(i).sceneImgFea = sceneImgFea(i,:);
    end

    ALL_SCENE = packupScene(ALL_HYPS, CAN_POOL{1});
    [ MINSCORE, MINTRANS] = compRoomMatchScore( ALL_SCENE, VALID_TRANS_GNDS );
    for i = 1:length(ALL_SCENE)
        ALL_HYPS(i).MINSCORE = MINSCORE(i,:);
        ALL_HYPS(i).MINTRANS = MINTRANS(i,:);
    end

    if strcmp(config.lossfunc, 'cost')
        for i = 1:length(ALL_SCENE)
            if isempty(ALL_SCENE{i})
                ALL_HYPS(i).COST = 4.5;
                continue;
            end
            ALL_HYPS(i).COST = roomLossFunction3D( gnd, ALL_SCENE{i}, H2G_R );
        end
    elseif strcmp(config.lossfunc, 'align')
        for i = 1:length(ALL_SCENE)
            if isempty(ALL_SCENE{i})
                ALL_HYPS(i).COST = 4.5;
                continue;
            end
            HYP_COMP = getHypInfo(ALL_SCENE{i});
            ALL_HYPS(i).COST = roomAlignmentMex(HYP_COMP, GND_COMP, 3, 0);
        end
    end
    if config.UPDATE_TO_DISK
        save(sprintf([bufname '/' config.sceneFeatFile], rid), 'ALL_HYPS');
    end
end

% sample near ground truth, only for training
if config.SAMPLE_CLOSE_GND
    [GOOD_HYPS, BESTRHID] = gndtruthSampling( ALL_POOL, gnd, config.samplenum );
    load(sprintf([bufname '/' config.objFeatFile], BESTRHID));
    ALL_POOL{BESTRHID}.sel_hyps.objfea = objfea;
    ALL_POOL{BESTRHID}.room.roomfea = roomfea;
    ALL_POOL{BESTRHID}.sel_hyps.anglesid = anglesid;
    [ sceneImgFea ] = compSceneHypsFeatureA( GOOD_HYPS, ALL_POOL{BESTRHID}.sel_hyps, ALL_POOL{BESTRHID}.room, omap, gc, cn );
    for i = 1:config.samplenum
        GOOD_HYPS(i).sceneImgFea = sceneImgFea(i,:);
    end

    ALL_SCENE = packupScene( GOOD_HYPS, ALL_POOL{BESTRHID});
    [ MINSCORE, MINTRANS] = compRoomMatchScore( ALL_SCENE, VALID_TRANS_GNDS );
    for j = 1:length(ALL_SCENE)
        GOOD_HYPS(j).MINSCORE = MINSCORE(j,:);
        GOOD_HYPS(j).MINTRANS = MINTRANS(j,:);
    end

    if strcmp(config.lossfunc, 'cost')
        for i = 1:length(ALL_SCENE)
            if isempty(ALL_SCENE{i})
                GOOD_HYPS(i).COST = 4.5;
                continue;
            end
            GOOD_HYPS(i).COST = roomLossFunction3D( gnd, ALL_SCENE{i}, H2G_R );
        end
    elseif strcmp(config.lossfunc, 'align')
        for i = 1:length(ALL_SCENE)
            if isempty(ALL_SCENE{i})
                GOOD_HYPS(i).COST = 4.5;
                continue;
            end
            HYP_COMP = getHypInfo(ALL_SCENE{i});
            GOOD_HYPS(i).COST = roomAlignmentMex(HYP_COMP, GND_COMP, 3, 0);
        end
    end

    % compute ground truth
    gnd_rotImg = imread([bufname '/' config.IMGLIST(aid).name config.gndImageSuffix]);
    gnd_omap = rotatePanorama( omap, H2G_R);
    gnd_gc = rotatePanorama( gc, H2G_R);
    gnd_cn = rotatePanorama( cn, H2G_R);

    [ GND_HYPS, ~ ] = groundTruthHypothesisNewFunc( gnd, gnd_omap, gnd_gc, gnd_cn, gnd_rotImg, config.typenum, gndBufName);
    [ MINSCORE, MINTRANS] = compRoomMatchScore( {gnd}, VALID_TRANS_GNDS );
    GND_HYPS.MINSCORE = MINSCORE;
    GND_HYPS.MINTRANS = MINTRANS;
    GND_HYPS.COST = 0;
    if config.UPDATE_TO_DISK
        save([bufname config.goodSceneFeatFile], 'GOOD_HYPS', 'GND_HYPS', 'BESTRHID');
    end
end

end

author: zouchuhang

repo: LayoutNet


text:
function [g1, g2] = lfmXrbfKernGradient(lfmKern, rbfKern, t1, t2, covGrad, meanVector)

% LFMXRBFKERNGRADIENT Compute gradient between the LFM and RBF kernels.
% FORMAT
% DESC computes the gradient of an objective function with respect
% to cross kernel terms between LFM and RBF kernels for
% the multiple output kernel.
% ARG lfmKern : the kernel structure associated with the LFM
% kernel.
% ARG rbfKern : the kernel structure associated with the RBF
% kernel.
% ARG t : inputs for which kernel is to be computed.
% RETURN g1 : gradient of objective function with respect to kernel
% parameters of LFM kernel.
% RETURN g2 : gradient of objective function with respect to kernel
% parameters of RBF kernel.
%
% FORMAT
% DESC computes the gradient of an objective function with respect
% to cross kernel terms between LFM and RBF kernels for
% the multiple output kernel.
% ARG lfmKern : the kernel structure associated with the LFM
% kernel.
% ARG rbfKern : the kernel structure associated with the RBF
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% RETURN g1 : gradient of objective function with respect to kernel
% parameters of LFM kernel.
% RETURN g2 : gradient of objective function with respect to kernel
% parameters of RBF kernel.
%
% FORMAT
% DESC computes the gradient of an objective function with respect
% to cross kernel terms between LFM and RBF kernels for
% the multiple output kernel.
% ARG lfmKern : the kernel structure associated with the LFM
% kernel.
% ARG rbfKern : the kernel structure associated with the RBF
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG meanVec : precomputed factor that is used for the switching dynamical
% latent force model.
% RETURN g1 : gradient of objective function with respect to kernel
% parameters of LFM kernel.
% RETURN g2 : gradient of objective function with respect to kernel
% parameters of RBF kernel.
%
% SEEALSO : multiKernParamInit, multiKernCompute, lfmKernParamInit, rbfKernParamInit
%
% COPYRIGHT : David Luengo, 2007, 2008
%
% MODIFICATIONS : Neil D. Lawrence, 2007
%
% MODIFICATIONS : Mauricio A. Alvarez, 2008, 2010

% KERN

subComponent = false; % This is just a flag that indicates if this kernel is part of a bigger kernel (SDLFM)

if nargin == 4
    covGrad = t2;
    t2 = t1;
elseif nargin == 6
    subComponent = true;
    if numel(meanVector)>1
        if size(meanVector,1) == 1,
            if size(meanVector, 2)~=size(covGrad, 2)
                error('The dimensions of meanVector don''t correspond to the dimensions of covGrad')
            end
        else
            if size((meanVector'), 2)~=size(covGrad,2)
                error('The dimensions of meanVector don''t correspond to the dimensions of covGrad')
            end
        end
    else
        if numel(t1)==1 && numel(t2)>1
            % matGrad will be row vector and so should be covGrad
            dimcovGrad = length(covGrad);
            covGrad = reshape(covGrad, [1 dimcovGrad]);
        elseif numel(t1)>1 && numel(t2)==1
            % matGrad will be column vector and sp should be covGrad
            dimcovGrad = length(covGrad);
            covGrad = reshape(covGrad, [dimcovGrad 1]);
        end
    end
end

if size(t1, 2) > 1 || size(t2, 2) > 1
    error('Input can only have one column');
end
if lfmKern.inverseWidth ~= rbfKern.inverseWidth
    error('Kernels cannot be cross combined if they have different inverse widths.')
end

if lfmKern.isNormalised ~= rbfKern.isNormalised
  error('Kernels cannot be cross combined if they have different normalization settings.')
end

m = lfmKern.mass;
D = lfmKern.spring;
C = lfmKern.damper;
S = lfmKern.sensitivity;

alpha = C/(2*m);
omega = sqrt(D/m-alpha^2);

sigma2 = 2/lfmKern.inverseWidth;
sigma = sqrt(sigma2);

if isreal(omega)
    gamma = alpha + j*omega;
    ComputeUpsilon1 = lfmComputeUpsilonMatrix(gamma,sigma2,t1, t2);
    if lfmKern.isNormalised
        K0 = S/(2*sqrt(2)*m*omega);
    else
        K0 = sigma*sqrt(pi)*S/(2*m*omega);
    end
else
    gamma1 = alpha + j*omega;
    gamma2 = alpha - j*omega;
    ComputeUpsilon1 = lfmComputeUpsilonMatrix(gamma2,sigma2,t1, t2);
    ComputeUpsilon2 = lfmComputeUpsilonMatrix(gamma1,sigma2,t1, t2);
    if lfmKern.isNormalised
        K0 = (S/(j*4*sqrt(2)*m*omega));
    else
        K0 = sigma*sqrt(pi)*S/(j*4*m*omega);
    end
end

g1 = zeros(1,5);

% Gradient with respect to m, C and D

for ind = 1:3 % Parameter (m, D or C)
    switch ind
        case 1  % Gradient wrt m
            gradThetaM = 1;
            gradThetaAlpha = -C/(2*(m^2));
            gradThetaOmega = (C^2-2*m*D)/(2*(m^2)*sqrt(4*m*D-C^2));
        case 2  % Gradient wrt D
            gradThetaM = 0;
            gradThetaAlpha = 0;
            gradThetaOmega = 1/sqrt(4*m*D-C^2);
        case 3  % Gradient wrt C
            gradThetaM = 0;
            gradThetaAlpha = 1/(2*m);
            gradThetaOmega = -C/(2*m*sqrt(4*m*D-C^2));
    end

    % Gradient evaluation

    if isreal(omega)
        gamma = alpha + j*omega;
        gradThetaGamma = gradThetaAlpha + j*gradThetaOmega;
        matGrad = -K0*imag(lfmGradientUpsilonMatrix(gamma,sigma2,t1, t2)*gradThetaGamma ...
            - (gradThetaM/m + gradThetaOmega/omega) ...
            * ComputeUpsilon1);
    else
        gamma1 = alpha + j*omega;
        gamma2 = alpha - j*omega;
        gradThetaGamma1 = gradThetaAlpha + j*gradThetaOmega;
        gradThetaGamma2 = gradThetaAlpha - j*gradThetaOmega;
        matGrad = K0*(lfmGradientUpsilonMatrix(gamma2,sigma2, t1, t2)*gradThetaGamma2 ...
            -  lfmGradientUpsilonMatrix(gamma1,sigma2, t1, t2)*gradThetaGamma1 ...
            - (gradThetaM/lfmKern.mass + gradThetaOmega/omega) ...
            * (ComputeUpsilon1 - ComputeUpsilon2));

    end
    if subComponent
        if size(meanVector,1) ==1,
            matGrad = matGrad*meanVector;
        else
            matGrad = (meanVector*matGrad).';
        end
    end
    g1(ind) = sum(sum(matGrad.*covGrad));
end

% Gradient with respect to sigma

if isreal(omega)
    gamma = alpha + j*omega;
    if lfmKern.isNormalised
        matGrad = -K0*imag(lfmGradientSigmaUpsilonMatrix(gamma,sigma2,t1, t2));
    else
        matGrad = -(sqrt(pi)*S/(2*m*omega)) ...
            * imag(ComputeUpsilon1 ...
            + sigma*lfmGradientSigmaUpsilonMatrix(gamma,sigma2,t1, t2));
    end
else
    gamma1 = alpha + j*omega;
    gamma2 = alpha - j*omega;
    if lfmKern.isNormalised
        matGrad = K0*(lfmGradientSigmaUpsilonMatrix(gamma2,sigma2,t1,t2) ...
            - lfmGradientSigmaUpsilonMatrix(gamma1,sigma2,t1,t2));
    else
        matGrad = (sqrt(pi)*S/(j*4*m*omega)) ...
            *(ComputeUpsilon1 - ComputeUpsilon2 ...
            + sigma*(lfmGradientSigmaUpsilonMatrix(gamma2,sigma2,t1,t2) ...
            - lfmGradientSigmaUpsilonMatrix(gamma1,sigma2,t1,t2)));
    end
end;

if subComponent
    if size(meanVector,1) ==1,
        matGrad = matGrad*meanVector;
    else
        matGrad = (meanVector*matGrad).';
    end
end

g1(4) = sum(sum(matGrad.*covGrad))*(-(sigma^3)/4); % temporarly introduced by MA
g2(1) = g1(4);

% Gradient with respect to S

if isreal(omega)
    if lfmKern.isNormalised
        matGrad = -(1/(2*sqrt(2)*m*omega))* imag(ComputeUpsilon1);
    else
        matGrad = -(sqrt(pi)*sigma/(2*m*omega))* imag(ComputeUpsilon1);
    end
else
    if lfmKern.isNormalised
        matGrad = (1/(j*4*sqrt(2)*m*omega))*(ComputeUpsilon1 - ComputeUpsilon2);
    else
        matGrad = (sqrt(pi)*sigma/(j*4*m*omega))*(ComputeUpsilon1 - ComputeUpsilon2);
    end
end;

if subComponent
    if size(meanVector,1) ==1,
        matGrad = matGrad*meanVector;
    else
        matGrad = (meanVector*matGrad).';
    end
end

g1(5) = sum(sum(matGrad.*covGrad));
g1 = real(g1);
% Gradient with respect to the "variance" of the RBF
g2(1) = 0; % Otherwise is counted twice, temporarly changed by MA
g2(2) = 0;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function S = cgf_fitS_sc2(A,X, sparsity, noise_var, beta, epsilon, sigma, tol, disp_ocbsol, disp_patnum, disp_stats, Sinit)
% cgf_fitS -- fit internal vars S to the data X using fast congugate gradient
%   Usage
%     S = cgf_fitS(A,X,noise_var,beta,sigma,
%                  [tol, disp_ocbsol, disp_patnum, disp_stats])
%   Inputs
%      A             basis functions
%      X             data vectors
%      noise_var     variance of the noise (|x-As|^2)
%      beta          steepness term for prior
%      sigma         scaling term for prior
%      tol           solution tolerance (default 0.001)
%      disp_ocbsol   display info from the fitting process
%      disp_patnum   display the pattern number
%      disp_stats    display summary statistics for the fit
%   Outputs
%      S             the estimated coefficients

maxiter=100;

[L,M] = size(A);
N = size(X,2);

if ~exist('tol','var');		tol = 0.001;			end
if ~exist('disp_ocbsol','var');	disp_ocbsol = 0;		end
if ~exist('disp_patnum','var');	disp_patnum = 1;		end
if ~exist('disp_stats','var');	disp_stats = 1;			end
if ~exist('maxiter','var');		maxiter = 8;			end
if ~exist('reduction','var');	reduction = 8;			end

% XXX: we don't use initialization for "log" sparsity function because of local optima
if ~exist('Sinit','var') %|| strcmp(sparsity, 'log') || strcmp(sparsity, 'huberL1') || strcmp(sparsity, 'epsL1')
	Sinit=A'*X;
	normA2=sum(A.*A)';
	for i=1:N
        Sinit(:,i)=Sinit(:,i)./normA2;
	end
    initiated = 0;
else
    initiated = 1;
end

if ~strcmp(sparsity, 'log') && ~strcmp(sparsity, 'huberL1') && ~strcmp(sparsity, ...
                                                      'epsL1')
	error('sparsity function is not properly specified!\n');
end

lambda=1/noise_var;

if strcmp(sparsity, 'huberL1') || strcmp(sparsity, 'epsL1')
	if ~exist('epsilon','var') || isempty(epsilon) || epsilon==0
		error('epsilon was not set properly!\n')
	end
end

S = zeros(M,N);
tic
if ~initiated
    if strcmp(sparsity, 'log')
        [S niters nf ng] = cgf_sc2(A,X,Sinit,0,lambda,beta,sigma,tol,maxiter, disp_ocbsol,disp_patnum);
    elseif strcmp(sparsity, 'huberL1')
        [S niters nf ng] = cgf_sc2(A,X,Sinit,1,lambda,beta,sigma,tol,maxiter, disp_ocbsol,disp_patnum, epsilon);
    elseif strcmp(sparsity, 'epsL1')
        [S niters nf ng] = cgf_sc2(A,X,Sinit,2,lambda,beta,sigma,tol,maxiter, disp_ocbsol,disp_patnum, epsilon);
    end
else
    if strcmp(sparsity, 'log')
        [S niters nf ng] = cgf_sc2(A,X,Sinit,0,lambda,beta,sigma,tol,maxiter, disp_ocbsol,disp_patnum);
    elseif strcmp(sparsity, 'huberL1')
        [S niters nf ng] = cgf_sc2(A,X,Sinit,1,lambda,beta,sigma,tol,maxiter, disp_ocbsol,disp_patnum, epsilon);
    elseif strcmp(sparsity, 'epsL1')
        [S niters nf ng] = cgf_sc2(A,X,Sinit,2,lambda,beta,sigma,tol,maxiter, disp_ocbsol,disp_patnum, epsilon);
    end
%      for i=1:size(X,2)
%          [aa,bb] = sort(abs(Sinit(:,i)));
%          bb = flipud(bb);
%          active = bb(1:M/reduction);
%          if strcmp(sparsity, 'log')
%              [S2 niters nf ng] = cgf_sc2(A(:,active),X(:,i),Sinit(:,i),0,lambda,beta,sigma,tol,maxiter, disp_ocbsol,disp_patnum);
%          elseif strcmp(sparsity, 'huberL1')
%              [S2 niters nf ng] = cgf_sc2(A(:,active),X(:,i),Sinit(:,i),1,lambda,beta,sigma,tol,maxiter, disp_ocbsol,disp_patnum, epsilon);
%          elseif strcmp(sparsity, 'epsL1')
%              [S2 niters nf ng] = cgf_sc2(A(:,active),X(:,i),Sinit(:,i),2,lambda,beta,sigma,tol,maxiter, disp_ocbsol,disp_patnum, epsilon);
%          end
%          S(active,i) = S2;
%      end
%      fprintf('%d',reduction);
end
t = toc;

if (disp_stats)
  fprintf(' aits=%6.2f af=%6.2f ag=%6.2f  at=%7.4f\n', ...
      niters/N, nf/N, ng/N, t/N);
end

author: thomas-koehler

repo: SupER


text:
function [ fx, fy ] = f01_f1 ( n, x, y )

%*****************************************************************************80
%
%% F01_F1 returns first derivatives of function 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N,1), Y(N,1), the evalution points.
%
%    Output, real FX(N,1), FY(N,1), the derivative values.
%
  t1(1:n,1) = exp ( - ( ( 9.0 * x(1:n,1) - 2.0 ).^2 ...
                    + ( 9.0 * y(1:n,1) - 2.0 ).^2 ) / 4.0 );
  t2(1:n,1) = exp ( - ( ( 9.0 * x(1:n,1) + 1.0 ).^2 ) / 49.0 ...
                    - ( 9.0 * y(1:n,1) + 1.0 ) / 10.0 );
  t3(1:n,1) = exp ( - ( ( 9.0 * x(1:n,1) - 7.0 ).^2 ...
                    + ( 9.0 * y(1:n,1) - 3.0 ).^2 ) / 4.0 );
  t4(1:n,1) = exp ( -   ( 9.0 * x(1:n,1) - 4.0 ).^2 ...
                    - ( 9.0 * y(1:n,1) - 7.0 ).^2 );

  fx(1:n,1) = ...
    - 3.375           * ( 9.0 * x(1:n,1) - 2.0 ) * t1 ...
    - ( 27.0 / 98.0 ) * ( 9.0 * x(1:n,1) + 1.0 ) * t2 ...
    - 2.25            * ( 9.0 * x(1:n,1) - 7.0 ) * t3 ...
    + 3.6             * ( 9.0 * x(1:n,1) - 4.0 ) * t4;

  fy(1:n,1) = ...
    - 3.375 * ( 9.0 * y(1:n,1) - 2.0 ) * t1 ...
    - 0.675                          * t2 ...
    - 2.25  * ( 9.0 * y(1:n,1) - 3.0 ) * t3 ...
    + 3.6   * ( 9.0 * y(1:n,1) - 7.0 ) * t4;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function point_num = gp_se_size ( dim_num, level_max )

%*****************************************************************************80
%
%% GP_SE_SIZE: Gauss Patterson Slow Exponential Growth.
%
%  Discussion:
%
%    The Gauss Patterson Slow Exponential Growth family assumes that, for the
%    underlying 1D rules, a precision of 2*L+1 is needed at level L.  Therefore,
%    the lowest possible order Gauss-Patterson rule is chosen that will achieve
%    that precision.  This retains a combination of the advantages of
%    nestedness and high accuracy.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Output, integer POINT_NUM, the total number of unique
%    points in the grids.
%

%
%  Special case.
%
  if ( level_max < 0 )
    point_num = 0;
    return
  end

  if ( level_max == 0 )
    point_num = 1;
    return
  end
%
%  Count the points in the 1D rule.
%
  order_1d = zeros ( level_max+1, 1 );

  order_1d(0+1) = 1;
  for level = 1 : level_max
    p = 5;
    o = 3;
    while ( p < 2 * level + 1 )
      p = 2 * p + 1;
      o = 2 * o + 1;
    end
    order_1d(level+1) = o;
  end
%
%  Count the new points in the 1D rule.
%
  new_1d = zeros(level_max+1,1);

  new_1d(0+1) = 1;
  for level = 1 : level_max
    new_1d(level+1) = order_1d(level+1) - order_1d(level-1+1);
  end
%
%  Count the number of points by counting the number of new points
%  associated with each level vector.
%
  level_1d = zeros ( dim_num, 1 );

  point_num = 0;

  for level = 0 : level_max

    level_1d = [];
    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );

      point_num = point_num + prod ( new_1d(level_1d(1:dim_num)+1) );

      if ( ~more )
        break
      end

    end

  end

  return
end


text:
function [out] = area_1(p1,p2,S,Smin,Smax,varargin)
%area_1

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Auxiliary function that calculates a variable contributing area.
% Constraints:  A <= 1
% @(Inputs):    p1   - linear scaling parameter [-]
%               p2   - exponential scaling parameter [-]
%               S    - current storage [mm]
%               Smin - minimum contributing storage [mm]
%               Smax - maximum contributing storage [mm]
%               varargin(1) - smoothing variable r (default 0.01)
%               varargin(2) - smoothing variable e (default 5.00)

if size(varargin,2) == 0
    out = min(1,p1.*(max(0,S-Smin)./(Smax-Smin)).^p2).*...
            (1-smoothThreshold_storage_logistic(S,Smin));                       % default smoothing
elseif size(varargin,2) == 1
    out = min(1,p1.*(max(0,S-Smin)./(Smax-Smin)).^p2).*...
            (1-smoothThreshold_storage_logistic(S,Smin,varargin(1)));           % user-specified smoothing
elseif size(varargin,2) == 2
     out = min(1,p1.*(max(0,S-Smin)./(Smax-Smin)).^p2).*...
            (1-smoothThreshold_storage_logistic(S,Smin,varargin(1),varargin(2))); % user-specified smoothing
end

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function u = log10(a)
%LOG10        slope base 10 logarithm  log10(a)
%

% written  12/06/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  u = a;

  u.r = log10(a.r);
  u.s = slopeconvexconcave('log10','1./(log(intval(10))*(%))',a,0);

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function dep = cfg_dep(varargin)

% This is the configuration dependency class
%
% Data structure
% ==============
% Description fields
%    * sname        - display name of dependency source
%    * src_exbranch - subsref/subsasgn struct referencing the dependency
%                     source exbranch
%    * src_output   - subsref/subsasgn struct referencing the dependency
%                     source output item
%    * tname        - display name of dependency target
%    * tgt_exbranch - subsref/subsasgn struct referencing the dependency
%                     target exbranch in the config tree
%    * tgt_input    - subsref/subsasgn struct referencing the dependency
%                     target item in the config tree
%    * tgt_spec     - an cfg_findspec that can be used to restrict matches
%                     of this dependency to certain input items - this will
%                     be checked in subsasgn checks for cfg_items. Defaults
%                     to {} - match all cfg_items.
%    * jtsubs       - subsref/subsasgn struct referencing the dependency
%                     target item in the job tree (this is currently not
%                     used and may be removed in future)
%
% Public Methods
% ==============
%
% Public internal Methods
% =======================
%    * subsasgn
%    * subsref
%    * display
%    * disp
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_dep.m 4863 2012-08-27 08:09:23Z volkmar $

rev = '$Rev: 4863 $'; %#ok

dep = class(struct('tname','Target', ...
                   'tgt_exbranch', struct('type',{},'subs',{}), ...
                   'tgt_input', struct('type',{},'subs',{}), ...
                   'tgt_spec', {{}}, ...
                   'jtsubs', struct('type',{},'subs',{}), ...
                   'sname','Source', ...
                   'src_exbranch', struct('type',{},'subs',{}), ...
                   'src_output', struct('type',{},'subs',{})), ...
            'cfg_dep');
switch nargin
    case 0
        return;
    case 1
        if isa(varargin{1},'cfg_dep')
            dep = varargin{1};
        else
            dep.sname                          = varargin{1};
        end;
    case 2
        dep.sname                              = varargin{1};
        dep.src_exbranch(1:numel(varargin{2})) = varargin{2};
    case 3
        dep.sname                              = varargin{1};
        dep.src_exbranch(1:numel(varargin{2})) = varargin{2};
        dep.src_output(1:numel(varargin{3}))   = varargin{3};
    case 4
        dep.sname                              = varargin{1};
        dep.src_exbranch(1:numel(varargin{2})) = varargin{2};
        dep.src_output(1:numel(varargin{3}))   = varargin{3};
        dep.tgt_spec                           = varargin{4};
    case 5
        dep.sname                              = varargin{1};
        dep.src_exbranch(1:numel(varargin{2})) = varargin{2};
        dep.src_output(1:numel(varargin{3}))   = varargin{3};
        dep.tgt_spec                           = varargin{4};
        dep.jtsubs(1:numel(varargin{5}))       = varargin{5};
    case 6
        dep.sname                              = varargin{1};
        dep.src_exbranch(1:numel(varargin{2})) = varargin{2};
        dep.src_output(1:numel(varargin{3}))   = varargin{3};
        dep.tgt_spec                           = varargin{4};
        dep.jtsubs(1:numel(varargin{5}))       = varargin{5};
        dep.tname                              = varargin{6};
    case 7
        dep.sname                              = varargin{1};
        dep.src_exbranch(1:numel(varargin{2})) = varargin{2};
        dep.src_output(1:numel(varargin{3}))   = varargin{3};
        dep.tgt_spec                           = varargin{4};
        dep.jtsubs(1:numel(varargin{5}))       = varargin{5};
        dep.tname                              = varargin{6};
        dep.tgt_exbranch(1:numel(varargin{7})) = varargin{7};
    case 8
        dep.sname                              = varargin{1};
        dep.src_exbranch(1:numel(varargin{2})) = varargin{2};
        dep.src_output(1:numel(varargin{3}))   = varargin{3};
        dep.tgt_spec                           = varargin{4};
        dep.jtsubs(1:numel(varargin{5}))       = varargin{5};
        dep.tname                              = varargin{6};
        dep.tgt_exbranch(1:numel(varargin{7})) = varargin{7};
        dep.tgt_input(1:numel(varargin{8}))    = varargin{8};
    otherwise
        cfg_message('matlabbatch:constructor:nargin', 'Wrong number of arguments.');
end;

author: spm


text:
% Most Software Machine Data File

%% scanimage.SI (ScanImage)

% Global microscope properties
objectiveResolution = 27.5;     % Resolution of the objective in microns/degree of scan angle

% Data file location

% Custom Scripts
startUpScript = 'ScanImageStartScript';     % Name of script that is executed in workspace 'base' after scanimage initializes
shutDownScript = '';     % Name of script that is executed in workspace 'base' after scanimage exits

fieldCurvatureZs = [];     % Field curvature for mesoscope
fieldCurvatureRxs = [];     % Field curvature for mesoscope
fieldCurvatureRys = [];     % Field curvature for mesoscope
useJsonHeaderFormat = false;     % Use JSON format for TIFF file header

fieldCurvatureTip = 0;
fieldCurvatureTilt = 0;

%% scanimage.components.Motors (SI Motors)
% SI Stage/Motor Component.
motorXYZ = {'ScopeHolder' 'ScopeHolder' 'ScopeHolder'};     % Defines the motor for ScanImage axes X Y Z.
motorAxisXYZ = [2 1 3];     % Defines the motor axis used for Scanimage axes X Y Z.
scaleXYZ = [1 -1 1];     % Defines scaling factors for axes.
backlashCompensation = [0 0 0];     % Backlash compensation in um (positive or negative)

%% scanimage.components.scan2d.RggScan (MINI2P_001)

acquisitionDeviceId = 'vDAQ0';     % RDI Device ID
acquisitionEngineIdx = 1;

resonantScanner = 'MEMS2000Hz';     % Name of the resonant scanner
xGalvo = '';     % Name of the x galvo scanner
yGalvo = 'MEMS2000Hz_slowaxis';     % Name of the y galvo scanner
beams = {'920nm AOM'};     % beam device names
fastZs = {'TLens'};     % fastZ device names
shutters = {'Laser' 'MEMS' 'PMT' 'LED'};     % shutter device names

channelsInvert = [true true false false];     % Logical: Specifies if the input signal is inverted (i.e., more negative for increased light signal)
keepResonantScannerOn = false;     % Always keep resonant scanner on to avoid drift and settling time issues

externalSampleClock = false;     % Logical: use external sample clock connected to the CLK IN terminal of the FlexRIO digitizer module
externalSampleClockRate = 8e+07;     % [Hz]: nominal frequency of the external sample clock connected to the CLK IN terminal (e.g. 80e6); actual rate is measured on FPGA
externalSampleClockMultiplier = 1;     % Multiplier to apply to external sample clock

extendedRggFov = 0;     % If true and x galvo is present, addressable FOV is combination of resonant FOV and x galvo FOV.

% Advanced/Optional
PeriodClockDebounceTime = 1e-07;     % [s] time the period clock has to be stable before a change is registered
TriggerDebounceTime = 5e-07;     % [s] time acquisition, stop and next trigger to be stable before a change is registered
reverseLineRead = 0;     % flips the image in the resonant scan axis
defaultFlybackTimePerFrame = 0.001;     % [s] default time to allow galvos to fly back after one frame is complete. overridden by cfg file
defaultFlytoTimePerScanfield = 0.001;     % [s] time to allow galvos to fly from one scanfield to the next. overridden by cfg file

% Aux Trigger Recording, Photon Counting, and I2C are mutually exclusive

% Aux Trigger Recording
auxTriggersTimeDebounce = 1e-07;     % [s] time after an edge where subsequent edges are ignored
auxTriggerLinesInvert = [false false false false];     % [logical] 1x4 vector specifying polarity of aux trigger inputs
auxTrigger1In = '';     % Digital input lines for aux trigger 1
auxTrigger2In = '';     % Digital input lines for aux trigger 2
auxTrigger3In = '';     % Digital input lines for aux trigger 3
auxTrigger4In = '';     % Digital input lines for aux trigger 4

% Signal Conditioning
disableMaskDivide = [false false false false];     % disable averaging of samples into pixels; instead accumulate samples

% I2C
i2cEnable = false;
i2cSdaPort = '';
i2cSclPort = '';
i2cAddress = 0;     % [byte] I2C address of the FPGA
i2cDebounce = 1e-07;     % [s] time the I2C signal has to be stable high before a change is registered
i2cStoreAsChar = false;     % if false, the I2C packet bytes are stored as a uint8 array. if true, the I2C packet bytes are stored as a string. Note: a Null byte in the packet terminates the string
i2cSendAck = true;     % When enabled FPGA confirms each packet with an ACK bit by actively pulling down the SDA line

% Laser Trigger
LaserTriggerPort = '';     % Digital input where laser trigger is connected.

% Trigger Outputs
frameClockOut = '/vDAQ0/D3.7';     % Output line for the frame clock
lineClockOut = '';     % Output line for the line clock
beamModifiedLineClockOut = '';     % Output line for beam clock
volumeTriggerOut = '';     % Output line for the volume clock

% Calibration data
scannerToRefTransform = [1 0 0;0 1 0;0 0 1];
virtualChannelSettings = [];
LaserTriggerDebounceTicks = 1;

enableHostPixelCorrection = false;
hostPixelCorrectionMultiplier = 500;

%% dabs.mirrorcle.ResonantAxis (MEMS2000Hz)
hAOZoom = '/vDAQ0/AO1';     % Analog output channel to be used to control resonant axis of the mirror e.g. '/vDAQ0/AO0'
hDOSync = '/vDAQ0/D1.7';     % String identifying the digital output for the sync signal e.g. '/vDAQ0/D0.1
hDOFilterX = '/vDAQ0/D3.2';     % String identifying the digital output for the X filter clock e.g. '/vDAQ0/D0.2
hDOFilterY = '/vDAQ0/D3.1';     % String identifying the digital output for the Y filter clock e.g. '/vDAQ0/D0.3

inputVoltageRange_Vpp = 8;     % Max input voltage range of the controller
angularRange_deg = 16;     % Max angular range of the device

% Default scan settings
syncPhase_deg = 100;
rampTime_s = 0;
xFilterClockFreq_Hz = 300000;
yFilterClockFreq_Hz = 150000;
xFilterClockEnable = 1;
yFilterClockEnable = 1;
nominalFrequency_Hz = 2000;

% Calibration Settings
amplitudeToLinePhaseMap = [4 -2.912e-06;4.5 -3.144e-06;5 -3.056e-06;5.333 -3.48e-06;5.714 -2.56e-06;6 -3.432e-06;6.154 -3.48e-06;6.667 -3.352e-06;7.273 -5.7912e-05;7.5 -3.416e-06;7.692 -3.504e-06;7.727 -3.416e-06;7.72727 0;8 -3.512e-06;8.182 6.3208e-05;8.333 -3.472e-06;8.5 -3.728e-06;8.636 -3.48e-06;9 -3.576e-06;9.091 -3.48e-06;9.5 -3.576e-06;10 0.00011808;11 -3.072e-06];     % translates an amplitude (degrees) to a line phase (seconds)
amplitudeLUT = zeros(0,2);     % translates a nominal amplitude (degrees) to an output amplitude (degrees)

%% dabs.generic.DigitalShutter (Laser)
DOControl = '/vDAQ0/D0.0';     % control terminal  e.g. '/vDAQ0/DIO0'
invertOutput = false;     % invert output drive signal to shutter
openTime_s = 0.1;     % settling time for shutter in seconds

%% dabs.generic.DigitalShutter (MEMS)
DOControl = '/vDAQ0/D0.1';     % control terminal  e.g. '/vDAQ0/DIO0'
invertOutput = false;     % invert output drive signal to shutter
openTime_s = 0.1;     % settling time for shutter in seconds

%% dabs.generic.DigitalShutter (PMT)
DOControl = '/vDAQ0/D0.2';     % control terminal  e.g. '/vDAQ0/DIO0'
invertOutput = false;     % invert output drive signal to shutter
openTime_s = 0.1;     % settling time for shutter in seconds

%% dabs.generic.DigitalShutter (LED)
DOControl = '/vDAQ0/D0.3';     % control terminal  e.g. '/vDAQ0/DIO0'
invertOutput = false;     % invert output drive signal to shutter
openTime_s = 0.1;     % settling time for shutter in seconds

%% dabs.generic.BeamModulatorFastAnalog (920nm AOM)
AOControl = '/vDAQ0/AO2';     % control terminal  e.g. '/vDAQ0/AO0'
AIFeedback = '';     % feedback terminal e.g. '/vDAQ0/AI0'

outputRange_V = [0 2];     % Control output range in Volts
feedbackUsesRejectedLight = false;     % Indicates if photodiode is in rejected path of beams modulator.
calibrationOpenShutters = {};     % List of shutters to open during the calibration. (e.g. {'Shutter1' 'Shutter2'}

powerFractionLimit = 1;     % Maximum allowed power fraction (between 0 and 1)

% Calibration data for the laser
powerFraction2ModulationVoltLut = [0 0;0.0157 0.2;0.0541 0.4;0.1165 0.6;0.2009 0.8;0.3047 1;0.4247 1.2;0.5583 1.4;0.7005 1.6;0.8493 1.8;1 2];
powerFraction2PowerWattLut = [0 0;1 164.6];
powerFraction2FeedbackVoltLut = zeros(0,2);
feedbackOffset_V = 0;

% Advanced Settings. Note: these settings are unused for vDAQ based systems
modifiedLineClockIn = '';     % Terminal to which external beam trigger is connected. Leave empty for automatic routing via PXI/RTSI bus
frameClockIn = '';     % Terminal to which external frame clock is connected. Leave empty for automatic routing via PXI/RTSI bus
referenceClockIn = '';     % Terminal to which external reference clock is connected. Leave empty for automatic routing via PXI/RTSI bus
referenceClockRate = 1e+07;     % if referenceClockIn is used, referenceClockRate defines the rate of the reference clock in Hz. Default: 10e6Hz

%% dabs.legacy.motor.LegacyMotor (ScopeHolder)
% Motor used for X/Y/Z motion, including stacks.

controllerType = 'PI C884/C863';     % If supplied, one of {'sutter.mp285', 'sutter.mpc200', 'thorlabs.mcm3000', 'thorlabs.mcm5000', 'scientifica', 'pi.e665', 'pi.e816', 'npoint.lc40x', 'bruker.MAMC'}.
comPort = '';     % Integer identifying COM port for controller, if using serial communication
customArgs = {};     % Additional arguments to stage controller. Some controller require a valid stageType be specified
invertDim = '$$$''---''';     % string with one character for each dimension specifying if the dimension should be inverted. '+' for normal, '-' for inverted
positionDeviceUnits = [-0.001 -0.001 -0.001];     % 1xN array specifying, in meters, raw units in which motor controller reports position. If unspecified, default positionDeviceUnits for stage/controller type presumed.
velocitySlow = [];     % Velocity to use for moves smaller than motorFastMotionThreshold value. If unspecified, default value used for controller. Specified in units appropriate to controller type.
velocityFast = [];     % Velocity to use for moves larger than motorFastMotionThreshold value. If unspecified, default value used for controller. Specified in units appropriate to controller type.
moveCompleteDelay = 0.1;     % Delay from when stage controller reports move is complete until move is actually considered complete. Allows settling time for motor
moveTimeout = '';     % Default: 2s. Fixed time to wait for motor to complete movement before throwing a timeout error
moveTimeoutFactor = '';     % (s/um) Time to add to timeout duration based on distance of motor move command

%% dabs.generic.FastZPureAnalog (TLens)
AOControl = '/vDAQ0/AO3';     % control terminal  e.g. '/vDAQ0/AO0'
AIFeedback = '';     % feedback terminal e.g. '/vDAQ0/AI0'
FrameClockIn = '';     % frame clock input terminal e.g. '/Dev1/PFI0'

parkPositionUm = 0;     % park position in micron
travelRangeUm = [-250 10];     % travel range in micron

voltsPerUm = -0.032;     % volts per micron
voltsOffset = 0;     % volts that sets actuator to zero position

% Calibration Data for uTLens
positionLUT = [-240 -267.5;-230 -250.625;-220 -233.75;-210 -217.5;-200 -205;-190 -191.25;-180 -180;-170 -168.958;-160 -159.167;-150 -150;-140 -139.375;-130 -130;-120 -121.25;-110 -113.542;-100 -104.583;-90 -96.25;-80 -86.875;-70 -80;-60 -71.925;-50 -65;-40 -56.2;-30 -46.2;-20 -34.95;-10 -21.2;0 0;10 25];     % Position LUT
feedbackVoltLUT = zeros(0,2);     % [Nx2] lut translating feedback volts into position volts

%% dabs.generic.GalvoPureAnalog (MEMS2000Hz_slowaxis)
AOControl = '/vDAQ0/AO0';     % control terminal  e.g. '/vDAQ0/AO0'
AOOffset = '';     % control terminal  e.g. '/vDAQ0/AO0'
AIFeedback = '';     % feedback terminal e.g. '/vDAQ0/AI0'

angularRange = 6;     % total angular range in optical degrees (e.g. for a galvo with -20..+20 optical degrees, enter 40)
voltsPerOpticalDegrees = 1.3;     % volts per optical degrees for the control signal
parkPosition = 0;     % park position in optical degrees
slewRateLimit = Inf;     % Slew rate limit of the analog output in Volts per second

% Calibration settings
feedbackVoltLUT = zeros(0,2);     % [Nx2] lut translating feedback volts into position volts
offsetVoltScaling = 1;     % scalar factor for offset volts

voltsOffset = -0.3;

author: kavli-ntnu


text:
%DISP_HHS2  display Hilbert-Huang spectrum
%
% DISP_HHS2(im,t,inf,sub,color)
% displays in a new figure the spectrum contained in matrix "im"
% (amplitudes in dB).
%
% inputs:  - im: image matrix (e.g., output of "toimage")
%          - t (optional): time instants (e.g., output of "toimage")
%          - inf (optional): -dynamic range in dB (wrt max)
%            default: inf = -20
%          - fs: sampling frequency
%          - sub: subset ratio
%          - color: 0 = grayscale ; 1 = color (default)
%
% use:  disp_hhs(im) ; disp_hhs(im,t) ; disp_hhs(im,inf)
%       disp_hhs(im,t,inf) ; disp_hhs(im,inf,fs) ; disp_hhs(im,[],fs)
%       disp_hhs(im,t,[],fs) ; disp_hhs(im,t,inf,fs)
%
%
% See also
%  emd, hhspectrum, toimage
%
% Modification of G. Rilling code
% gabriel.rilling@ens-lyon.fr

function disp_hhs2(varargin)

error(nargchk(1,6,nargin));
fs = 0;
inf = -20;
im = varargin{1};
t = 1:size(im,2);
sub = 1;
color = 1;
switch nargin
  case 1
    %raf
  case 2
    if isscalar(varargin{2})
      inf = varargin{2};
    else
      t = varargin{2};
    end
  case 3
    if isvector(varargin{2})
      t = varargin{2};
      inf = varargin{3};
    else
      inf = varargin{2};
      fs = varargin{3};
    end
  case 4
    t = varargin{2};
    inf = varargin{3};
    fs = varargin{4};
  case 5
    if isscalar(varargin{5})
      t = varargin{2};
      inf = varargin{3};
      fs = varargin{4};
      sub = varargin{5};
    end
  case 6
    if isscalar(varargin{6})
      t = varargin{2};
      inf = varargin{3};
      fs = varargin{4};
      sub = varargin{5};
      color = varargin{6};
    end
end

if isempty(inf)
  inf = -20;
end

if inf > 0
  inf = -inf;
elseif inf == 0
  error('inf must be nonzero')
end

if color~=0
    color = 1;
end

M=max(max(im));

warning off
im = 10*log10(im/M);
warning on

figure
subplot(6,1,[2:6]);   %AJOUT
if fs == 0
  imagesc(t,[0,pi/sub],im(1:floor(end/sub),:),[inf,0]);
  ylabel('frequency (rad/s)')
else
  imagesc(t,[0,0.5*fs],im,[inf,0]);
  ylabel('frequency')
end
if color == 0
    colormap(1.-gray);
end
set(gca,'YDir','normal')
xlabel('time')


text:
function J = computeCostMulti(X, y, theta)
%COMPUTECOSTMULTI Compute cost for linear regression with multiple variables
%   J = COMPUTECOSTMULTI(X, y, theta) computes the cost of using theta as the
%   parameter for linear regression to fit the data points in X and y

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta
%               You should set J to the cost.

y_pred = X * theta; % X.shape=(m,3) theta.shape=(3,1)
J = 1/(2*m) * sum((y_pred - y).^2);

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

end

author: imLogM


text:
function bvhPlayFile(fileName)

% BVHPLAYFILE Play motion capture data from a bvh file.
% FORMAT
% DESC plays motion capture data by reading in a bvh file from disk.
% ARG fileName : the name of the file to read in, in bvh format.
%
% COPYRIGHT : Neil D. Lawrence, 2005
%
% SEEALSO : acclaimPlayFile, bvhReadFile, skelPlayData

% MOCAP

[skel, channels, frameLength] = bvhReadFile(fileName);
skelPlayData(skel, channels, frameLength);


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open 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.
%
% dslam_open 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 dslam_open. If not, see <http://www.gnu.org/licenses/>.

function plotAccuracy(accuracy_measurements, num_robots, vo_end_time, ...
    orb_ate)

times = cellfun(@(x) x.time, accuracy_measurements);
num_meas = numel(accuracy_measurements);

robot_accuracies = cell(num_robots, 1);
mg_sizes = cell(num_robots, 1);
og_sizes = cell(num_robots, 1);
for robot_i = 1:num_robots
    robot_accuracies{robot_i} = zeros(num_meas, 1);
    og_sizes{robot_i} = zeros(num_meas, 1);
    mg_sizes{robot_i} = zeros(num_meas, 1);
    for meas_i = 1:num_meas
        meas = accuracy_measurements{meas_i};
        for mg_i = 1:numel(meas.matched_groups)
            if (any(meas.matched_groups(mg_i).members == robot_i))
                mg_sizes{robot_i}(meas_i) = ...
                    meas.matched_groups(mg_i).num_frames;
                break;
            end
        end
        for og_i = 1:numel(meas.optimized_groups)
            if (any(meas.optimized_groups(og_i).members == robot_i))
                robot_accuracies{robot_i}(meas_i) = ...
                    meas.optimized_groups(og_i).ATE;
                og_sizes{robot_i}(meas_i) = ...
                    meas.optimized_groups(og_i).num_frames;
                break;
            end
        end
    end
end
yyaxis left;
ate = plot(times, robot_accuracies{1}, '-');
hold on;
for i = 2:num_robots
    plot(times, robot_accuracies{i}, '-');
end
%orb = plot(xlim, [orb_ate orb_ate], 'b--');
hold off;
ylabel('ATE [m]');
xlabel('time [s]');
grid on;

yyaxis right;
opted = plot(times, og_sizes{1}, '-');
hold on;
for i = 2:num_robots
    plot(times, og_sizes{i}, '-');
end
%for i = 1:num_robots
%    matched = plot(times, mg_sizes{i}, '.');
%end
et = plot([vo_end_time vo_end_time], ylim, 'k--');
hold off;
ylabel('frame count in CC');
title('accuracy development');
set(gca, 'XLim',[min(times) max(times)]);
%legend([ate, opted, matched, et, orb], 'ATE', 'optimized', 'matched', ...
%    'VO end time', 'ORB SLAM','Location', 'NorthWest');
legend(et, 'VO end time', 'Location', 'NorthWest');

end

author: uzh-rpg


text:
function w = cvx( v, b )

%CVX   The CVX disciplined convex programming system.
%   CVX is a modeling framework for building, constructing, and solving
%   disciplined convex programs. CVX has online help for many functions and
%   operations, divided into several subsection:
%      cvx/commands   - Top-level commands to control CVX
%      cvx/keywords   - Keywords for declaring variables and objectives
%      cvx/builtins   - Built-in operators and functions supported in CVX models
%      cvx/functions  - Additional functions added by CVX
%      cvx/sets       - Definitions of common convex sets
%      cvx/structures - Matrix structure definitions
%   CVX also provides an extensive user guide in PDF format, which is found
%   in its top directory. This directory can be found by typing
%      cvx_where
%   at the command prompt.

global cvx___
switch nargin,
    case 2,
        w = cvx___.obj;
        switch numel( v ),
            case 2,
                w.size_ = v;
            case 1,
                w.size_ = [ v, 1 ];
            case 0,
                % w.size_ = [ 0, 1 ];
            otherwise,
                if v(end) == 1,
                    v = v(1:find(v>1,1,'last'));
                    v(1,end+1:2) = 1;
                end
                w.size_ = v;
        end
        if issparse( b )
            w.basis_ = b;
        elseif isempty( b )
            w.basis_ = sparse( 1, prod( v ) );
        else
            w.basis_ = sparse( b );
        end
    case 1,
        if isnumeric( v ),
            w = cvx___.obj;
            w.size_ = size( v );
            if issparse( v ),
                w.basis_ = reshape( v, 1, prod(s) );
            else
                w.basis_ = sparse( v(:)' );
            end
        elseif isa( v, 'cvx' ),
            w = v;
            return
        else
            w = class( v, 'cvx' );
        end
    case 0,
        w = cvx___.obj;
end
w.id_ = cvx___.id + 1;
cvx___.id = w.id_;

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.


text:
randomseed(0);
x = randgamma(repmat(3,1,3));
randomseed(0);
x2 = randgamma(repmat(3,1,3));
assert(all(x == x2));
seed = randomseed;
x = randgamma(repmat(3,1,3));
randomseed([4 5 6]);
randomseed(seed);
x2 = randgamma(repmat(3,1,3));
assert(all(x == x2));
fprintf('randomseed test passed.\n');

author: andrewssobral

repo: lrslibrary


text:
function [r_j] = calcResidual(p_f_G, camStates, measurements)
%CALCRESIDUAL Calculates the residual for a feature position

% measurements is 2 x M_j
% camStates is a cell array of the camState structs for the states
%   included in measurements
    r_j = NaN(2*size(camStates,2), 1);
    for i = 1:size(camStates,2)

        C_CG = quatToRotMat(camStates{i}.q_CG);
        p_f_C = C_CG * (p_f_G - camStates{i}.p_C_G);

        zhat_i_j = p_f_C(1:2)/p_f_C(3);

        iStart = 2*(i-1)+1;
        iEnd = 2*i;
        r_j(iStart:iEnd) = measurements(:,i) - zhat_i_j;
    end

end

author: utiasSTARS

repo: msckf-swf-comparison


text:
function z = rsum(fit)
%
% function to extract log-likelihood and degrees-of-freedom
% from locfit fit.
%
% order of returned vector:   df0 df1 llk.
%

fp = fit.fit_points;
gf = fp.kappa;

z = gf([2 3 1]);

return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function smean = strongMean(data, thr_perc, thr_perc_global, n_sigma)
% Returns the mean per column with a certaint percentile of values in data
% The code uses only the data filtered by the requested percentile to estimate
% an std column by column
% With this estimation only the data within n_sigma range are used for the
% robust mean estimation
%
% SYNTAX:
%   smean = strongMean(data, thr_perc, thr_perc_global, n_sigma)
%
% INPUT:
%   data                matrix of values
%   thr_perc            percentile requested [0 1] per column
%   thr_perc_global     percentile requested [0 1] global
%   n_sigma             number of sigmas
%
% OUTPUT:
%   smean   strong mean per column

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Andrea Gatti
%  Contributors:     Andrea Gatti ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------
    if nargin < 2
        thr_perc = 1;
    end

    if nargin < 4
        n_sigma = 2;
    end

    % remove the median
    sensor = data;
    smean = zeros(1, size(data, 2));
    for c = 1 : size(data, 2)
        if any(nan2zero(sensor(:, c)))
            sensor(:, c) = sensor(:, c) - median(sensor(:, c), 'omitnan');
        end
    end

    if nargin >= 3 && ~isnan(thr_perc_global) && (thr_perc_global < 1)
        thr = perc(serialize(noNaN(abs(sensor))), thr_perc_global);
        id_ko = abs(sensor) > thr;
    else
        id_ko = false(size(sensor));
    end

    for c = 1 : size(data, 2)
        if any(nan2zero(sensor(:, c)))
            if ~isempty(thr_perc) && thr_perc < 1
                id_ok = abs(sensor(:, c)) < perc(noNaN(abs(sensor(:, c))), thr_perc) & ~id_ko(:, c);
            else
                id_ok = ~id_ko(:, c);
            end

            id_ok = abs(sensor(:, c)) < n_sigma * std(sensor(id_ok, c), 'omitnan');
            if any(id_ok)
                smean(c) = mean(data(id_ok, c), 'omitnan');
            else
                smean(c) = median(data(:, c), 'omitnan');
            end
        end
    end
end

author: goGPS-Project


text:
function a = wilson_eigen_right ( )

%*****************************************************************************80
%
%% WILSON_EIGEN_RIGHT returns the right eigenvectors of the WILSON matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A(4,4), the right eigenvector matrix.
%

%
%  Note that the matrix entries are listed by row.
%
  a = [ ...
   0.380262074390714, ...
   0.396305561186082, ...
   0.093305039089285, ...
   0.830443752841578; ...
   0.528567849528642, ...
   0.614861280394151, ...
  -0.301652326903523, ...
  -0.501565058582058; ...
   0.551954849631663, ...
  -0.271601039711768, ...
   0.760318430013036, ...
  -0.208553600252039; ...
   0.520924780743657, ...
  -0.625396181050490, ...
  -0.567640668325261, ...
   0.123697458332363 ];

  return
end


text:
function p = norm_density(x, mu, sigma)
% Compute density of normal distribution function
%
% file:      	norm_density.m, (c) Matthew Roughan, Tue Jul 21 2009
% directory:   /home/mroughan/src/matlab/NUMERICAL_ROUTINES/
% created: 	Tue Jul 21 2009
% author:  	Matthew Roughan
% email:   	matthew.roughan@adelaide.edu.au
%
% Calculate normal density in 1D.
%
% INPUTS:
%   x      = a vector of points at which to calculate the normal CDF
%   mu     = mean of the normal
%   sigma  = std dev. of the normal
%
%
p = exp( -(x-mu).^2 / (2*sigma^2) ) / (sigma * sqrt(2*pi) );


text:
function [yr] = us2yr(us)
% Convert time from microseconds to Julian years (365.25 days).
% Chad Greene 2012
yr = us*3.168808781403e-14 ;


text:
% a very basic UTE-like sequence, without ramp-sampling, ramp-RF and other
% tricks yet. Achieves TE in the range of 300-400 us

% set system limits
sys = mr.opts('MaxGrad', 28, 'GradUnit', 'mT/m', ...
    'MaxSlew', 100, 'SlewUnit', 'T/m/s', 'rfRingdownTime', 20e-6, ...
    'rfDeadTime', 100e-6, 'adcDeadTime', 10e-6);

seq=mr.Sequence(sys);           % Create a new sequence object
fov=250e-3; Nx=256;             % Define FOV and resolution
alpha=10;                       % flip angle
sliceThickness=3e-3;            % slice
TR=10e-3;                       % TR
Nr=128;                         % number of radial spokes
delta= 2* pi / Nr;              % angular increment; try golden angle pi*(3-5^0.5) or 0.5 of it
ro_duration=2.56e-3;            % read-out time: controls RO bandwidth and T2-blurring
ro_os=2;                        % oversampling
ro_asymmetry=1;                 % 0: fully symmetric 1: half-echo
minRF_to_ADC_time=50e-6;        % the parameter wich defines TE (together with the RO asymmetyry)

% more in-depth parameters
rfSpoilingInc=117;              % RF spoiling increment

% Create alpha-degree slice selection pulse and gradient
[rf, gz, gzReph] = mr.makeSincPulse(alpha*pi/180,'Duration',1e-3,...
    'SliceThickness',sliceThickness,'apodization',0.5,'timeBwProduct',2,...
    'centerpos',1,'system',sys);

% Align RO assymmetry to ADC samples
Nxo=round(ro_os*Nx);
ro_asymmetry = round(ro_asymmetry*Nxo/2)/Nxo*2; % check whether we need to use 2Nx or so...
% Define other gradients and ADC events
deltak=1/fov/(1+ro_asymmetry);
ro_area=Nx*deltak;
gx = mr.makeTrapezoid('x','FlatArea',ro_area,'FlatTime',ro_duration,'system',sys);
adc = mr.makeAdc(Nxo,'Duration',gx.flatTime,'Delay',gx.riseTime,'system',sys);
gxPre = mr.makeTrapezoid('x','Area',-(gx.area-ro_area)/2 -gx.amplitude*adc.dwell/2 - ro_area/2*(1-ro_asymmetry),'system',sys);

% gradient spoiling
gxSpoil=mr.makeTrapezoid('x','Area',0.2*Nx*deltak,'system',sys);

% Calculate timing
%ceil((TE - mr.calcDuration(gxPre) - gz.fallTime - gz.flatTime/2 ...
%    - mr.calcDuration(gx)/2)/seq.gradRasterTime)*seq.gradRasterTime;
% calculate actual achieved TE
TE = gz.fallTime + mr.calcDuration(gxPre,gzReph)+gx.riseTime + adc.dwell*Nxo/2*(1-ro_asymmetry);
delayTR=ceil((TR - mr.calcDuration(gxPre,gzReph) - mr.calcDuration(gz) ...
    - mr.calcDuration(gx))/seq.gradRasterTime)*seq.gradRasterTime;
assert(all(delayTR>=mr.calcDuration(gxSpoil)));

fprintf('TE= %d us\n', round(TE*1e6));

if mr.calcDuration(gzReph) > mr.calcDuration(gxPre)
    gxPre.delay=mr.calcDuration(gzReph) - mr.calcDuration(gxPre);
end

rf_phase=0;
rf_inc=0;

for i=1:Nr
    for c=1:2
        rf.phaseOffset=rf_phase/180*pi;
        adc.phaseOffset=rf_phase/180*pi;
        rf_inc=mod(rf_inc+rfSpoilingInc, 360.0);
        rf_phase=mod(rf_phase+rf_inc, 360.0);
        % UTE: alternate GZ
        gz.amplitude=-gz.amplitude;
        gzReph.amplitude=-gzReph.amplitude;
        %
        seq.addBlock(rf,gz);
        phi=delta*(i-1);
        gpc=gxPre;   gps=gxPre;   gpc.amplitude=gxPre.amplitude*cos(phi);   gps.amplitude=gxPre.amplitude*sin(phi);   gps.channel='y';
        grc=gx;      grs=gx;      grc.amplitude=gx.amplitude*cos(phi);      grs.amplitude=gx.amplitude*sin(phi);      grs.channel='y';
        gsc=gxSpoil; gss=gxSpoil; gsc.amplitude=gxSpoil.amplitude*cos(phi); gss.amplitude=gxSpoil.amplitude*sin(phi); gss.channel='y';
        seq.addBlock(gpc,gps,gzReph);
        seq.addBlock(grc,grs,adc);
        seq.addBlock(gsc,gss,mr.makeDelay(delayTR));
    end
end

%% check whether the timing of the sequence is correct
[ok, error_report]=seq.checkTiming;

if (ok)
    fprintf('Timing check passed successfully\n');
else
    fprintf('Timing check failed! Error listing follows:\n');
    fprintf([error_report{:}]);
    fprintf('\n');
end

%%
seq.plot();

%% plot gradients to check for gaps and optimality of the timing
gw=seq.waveforms_and_times();
figure; plot(gw{1}(1,:),gw{1}(2,:),gw{2}(1,:),gw{2}(2,:),gw{3}(1,:),gw{3}(2,:)); % plot the entire gradient shape

%% k-space trajectory calculation
[ktraj_adc, t_adc, ktraj, t_ktraj, t_excitation, t_refocusing] = seq.calculateKspacePP();

% plot k-spaces
figure; plot(t_ktraj, ktraj'); % plot the entire k-space trajectory
hold; plot(t_adc,ktraj_adc(1,:),'.'); % and sampling points on the kx-axis
figure; plot(ktraj(1,:),ktraj(2,:),'b'); % a 2D plot
axis('equal'); % enforce aspect ratio for the correct trajectory display
hold;plot(ktraj_adc(1,:),ktraj_adc(2,:),'r.'); % plot the sampling points

%
seq.setDefinition('FOV', [fov fov sliceThickness]);
seq.setDefinition('Name', 'ute');

seq.write('ute.seq');       % Write to pulseq file

%seq.install('siemens');

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
% runpca() -  perform principal component analysis (PCA) using singular value
%             decomposition (SVD) using Matlab svd() or svds()
%                        >> inv(eigvec)*data = pc;
% Usage:
%    >> [pc,eigvec,sv] = runpca(data);
%    >> [pc,eigvec,sv] = runpca(data,num,norm)
%
% Inputs:
%   data   - input data matrix (rows are variables, columns observations)
%   num    - number of principal comps to return  {def|0|[] -> rows in data}
%   norm   - 1/0 = do/don't normalize the eigvec's to be equivariant
%                                                {def|0 -> no normalization}
% Outputs:
%   pc     - the principal components, i.e.        >> inv(eigvec)*data = pc;
%   eigvec - the inverse weight matrix (=eigenvectors). >> data = eigvec*pc;
%   sv     - the singular values (=eigenvalues)
%
% Author: Colin Humphries, CNL / Salk Institute, 1997
%
% See also: runica()

% Copyright (C) Colin Humphries, CNL / Salk Institute, Aug, 1997
%
% 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% 01/31/00 renamed runpca() and improved usage message -sm
% 01-25-02 reformated help & license, added links -ad

function [pc,M,S] = runpca(data,N,norm)

BIG_N = 50; % for efficiency, switch to sdvs() when BIG_N<=N or N==rows

if nargin < 1
  help runpca
  return
end

rows = size(data,1);

% remove the mean
for i = 1:rows
  data(i,:) = data(i,:) - mean(data(i,:));
end

if nargin < 3
  norm = 0;
elseif isempty(norm)
  norm = 0;
end

if nargin < 2
  N = 0;
end
if isempty(N)
  N = 0;
end

if N == 0  | N == rows
  N = rows;
  [U,S,V] = svd(data',0);   % performa SVD
  if norm == 0
    pc = U';
    M = (S*V')';
  else % norm
    pc = (U*S)';
    M = V;
  end
else
  if N > size(data,1)
    error('N must be <= the number of rows in data.')
  end
  %if N <= BIG_N | N == rows
  %[U,S,V] = svd(data',0);
  %else
  [U,S,V] = svds(data',N);
  %end
  if norm == 0
    pc = U';
    M = (S*V')';
  else % norm
    pc = (U*S)';
    M = V;
  end
  %if N > BIG_N & N < rows
  %pc = pc(1:N,:);
  %M = M(:,1:N);
  %end
end
%S = diag(S(1:N,1:N));

author: PatternRecognition

repo: OpenBMI


text:
function jac = p32_jac ( neqn, t, y )

%*****************************************************************************80
%
%% P32_JAC evaluates the jacobian for problem p32.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the jacobian.
%
%    Output, real JAC(NEQN,NEQN), the jacobian matrix.
%
  jac = zeros ( neqn, neqn );

  beta = p32_param ( 'GET', 'BETA', [] );
  rho = p32_param ( 'GET', 'RHO', [] );
  sigma = p32_param ( 'GET', 'SIGMA', [] );

  jac(1,1) = - sigma;
  jac(1,2) =   sigma;
  jac(1,3) =   0.0;

  jac(2,1) =   rho - y(3);
  jac(2,2) = - 1.0;
  jac(2,3) = - y(1);

  jac(3,1) =   y(2);
  jac(3,2) =   y(1);
  jac(3,3) =   - beta;

  return
end


text:
report_this_filefun(mfilename('fullpath'));

update(mainmap())
[x, y] = ginput(1) ;   %use you cursor to fix the location of interest

di =  sqrt(((a.Longitude-x)*cos(pi/180*y)*111).^2 + ((a.Latitude-y)*111).^2) ;

[s,is] = sort(di);
di = sort(di);
newt2 = a(is(:,1),:) ;

bv2 = [];
bv3 = [] ;
me = [];
def = {'150'};
ni2 = inputdlg('Number of events in each window?','Input',1,def);
l = ni2{:};
ni = str2double(l);

think

for i = 1:ni/10:length(newt2)-ni
    % [bv magco,  stan] =  bvalcalc(newt2(i:i+ni,:));
    [bv magco stan ] =  bvalca2(newt2(i:i+ni,:));

    bv2 = [bv2 ; magco di(i)];
end

% Find out of figure already exists
%
[existFlag,figNumber]=figure_exists('Mc with time',1);
newdepWindowFlag=~existFlag;
bdep= figNumber;

% Set up the window

if newdepWindowFlag
    Mcfig = figure_w_normalized_uicontrolunits( ...
        'Name','Mc with time',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'NextPlot','replace', ...
        'backingstore','on',...
        'Visible','on');

    uicontrol('Units','normal',...
        'Position',[.0 .85 .08 .06],'String','Info ',...
         'Callback','infoz(1)');
    matdraw
end

hold on
figure_w_normalized_uicontrolunits(Mcfig)
hold on
delete(gca)
delete(gca)
axis off

rect = [0.15 0.30 0.7 0.45];
axes('position',rect)
pl = plot(bv2(:,2),bv2(:,1),'^r');
set(pl,'LineWidth',1.5,'MarkerSize',10,...
    'MarkerFaceColor','y','MarkerEdgeColor','r')
hold on
pl = plot(bv2(:,2),bv2(:,1),'b')
set(pl,'LineWidth',1.0)

grid
set(gca,'Color',color_bg)
set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
%
ylabel('Mc')
%set(gca,'Xlim',[t0b teb]);

xlabel('Distance')
tist = [  name ' - b(t), ni = ' num2str(ni) ];
title(tist)

nt = [];
con = 0;

ms = round(bv2(:,1)*10)/10;
for  m = max(bv2(:,1)):-0.1:min(bv2(:,1))

    % find comp level and times.
    i = max(find(abs(m-ms) < 0.01));
    if isempty(i) == 0
        con = con+1;
        nt = [nt ; m bv2(i,2)];
        if con > 1 &  nt(con,2) < nt(con-1,2) ; nt(con,:) = []; con = con-1; end
    end
end

nt(1,2) = min(newt2.Date)
i = max(find((ms-min(ms) > 0.01)));
% nt(con,2) = bv2(i+1,2);

%hold on
%plot(nt(:,2),nt(:,1),'h')
%hold on
%plot(nt(:,2),nt(:,1),'r','LineWidth',2);

done

author: CelsoReyes


text:
function [] = autoWrite(S,fs,method)
% 'S':Len * Num
if(~exist('method','var'))
    method = [];
end
num = size(S,2);
S = S/max(max(abs(S)));
for i = 1:num
    eval(strcat('audiowrite(''S_',method,num2str(i),'.wav''',',S(:,',num2str(i),'),','fs);'));
end

return

repo: Voice-Separation-and-Enhancement


text:
function [kcal] = cal2kcal(cal)
% Convert energy or work from calories to kilocalories.
% Note: these calories are different from the capitalized Calories found on
% American cereal boxes.  I suppose we'd be afraid to eat a candy bar with
% 250,000 calories in it.  Because then we'd feel like fatties.
% Chad A. Greene 2012
kcal = cal/1000;

author: Sable


text:
% [INPUT]
% data = A float t-by-n matrix (-Inf,Inf) representing the model input.
% normalize = A boolean that indicates whether to normalize the model input (optional, default=true).
%
% [OUTPUT]
% coefficients = A float n-by-n matrix (-Inf,Inf) representing the PCA coefficients.
% scores = A float t-by-n matrix (-Inf,Inf) representing the PCA scores.
% explained = A column vector of floats [0,100] of length n representing the percentage of total variance explained by each PCA component.

function [coefficients,scores,explained] = pca_shorthand(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('data',@(x)validateattributes(x,{'double'},{'real' '2d' 'nonempty'}));
        ip.addOptional('normalize',true,@(x)validateattributes(x,{'logical'},{'scalar'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    data = validate_input(ipr.data);
    normalize = ipr.normalize;

    nargoutchk(3,3);

    [coefficients,scores,explained] = pca_shorthand_internal(data,normalize);

end

function [coefficients,scores,explained] = pca_shorthand_internal(data,normalize)

    if (normalize)
        for i = 1:size(data,2)
            c = data(:,i);

            m = mean(c,'omitnan');

            s = std(c,'omitnan');
            s(s == 0) = 1;

            data(:,i) = (c - m) ./ s;
        end
    end

    [coefficients,scores,~,~,explained] = pca(data,'Economy',false);

end

function data = validate_input(data)

    nan_indices = any(isnan(data),1);
    data(:,nan_indices) = [];

    [t,n] = size(data);

    if ((t < 5) || (n < 2))
        error('The value of ''data'' is invalid. Expected input to be a matrix with a minimum size of 5x2, after the exclusion of time series containing NaN values.');
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function[]=yoffset(i1,i2)
%YOFFSET  Offsets lines in the y-direction after plotting.
%
%   YOFFSET allows data to be manipulated after it has been plotted, by
%   operating on the data stored in the figure itself.  This may be used
%   only with 2-D line plots (not symbol plots, contour plots, etc.).
%
%   YOFFSET offsets all lines in the current axes a specified amount.
%
%   YOFFSET(N) or YOFFSET N, where N is a *real* number, offsets each
%   column of the Y-data by an amount N from the previous column.  N must
%   be a number, not a variable whose value is a number.
%
%   YOFFSET(N) or YOFFSET N, where N is an *imaginary* number, offsets
%   each column of the Y-data by an amount imag(N)*DY from the previous
%   column, where DY is the Y-axis length of the original (unoffset) plot.
%
%   YOFFSET(M,N) or YOFFSET M N, applies the same offsets to each of M
%   groups of lines. For example, for complex-valued data Z use UVPLOT(Z)
%   followed by YOFFSET 2 N.
%
%   YOFFSET with no arguments returns the data to its original (unoffset)
%   orientation [as do YOFFSET(0) and YOFFSET(0i)].
%
%   YOFFSET also allow multiple axes to be manipulated simultaneously.
%
%   YOFFSET LOCK and YOFFSET UNLOCK lock and unlock all Y-data in the
%   current figure. When LOCK is on, calls to YOFFSET are applied to each
%   subplot individually.
%
%   See also YOFFSET, LINESTYLE.
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2000--2015 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(i1,'--t')
    return
end

bcontinue=1;
if nargin==1;
  N=10000;  %Very large number of lines
end
if nargin==2
  N=i1;
  i1=i2;
  if ischar(N)
    N=str2double(N);
  end
end

if nargin>0
   if ischar(i1)
      if isempty(str2double(i1))
           i1=deblank(i1);
           i1=fliplr(deblank(fliplr(i1)));
           if strcmpi(i1,'lock')
               if isappdata(gca,'lastyoffset')
                    lastoff=getappdata(gca,'lastyoffset');
                    setappdata(gcf,'yoffsetlock',1)
                    yoffsetloop(0,N);
                    yoffsetloop(lastoff,N);
               else
                    setappdata(gcf,'yoffsetlock',1)
               end
               bcontinue=0;
           end
           if strcmpi(i1,'unlock')
               setappdata(gcf,'yoffsetlock',0)
               bcontinue=0;
           end
      end
   end
end

if bcontinue
   if nargin>0
      delta=i1;
   else
      delta=0;
   end

   h=linehandles;
   if plotmodified(h)
   %If the plot has been changed, the old "lastyoffset" may no
   %longer be correct.  If so, we set it to zero
       setappdata(gca,'lastyoffset',0)
   end
   yoffsetloop(delta,N);
end

function[]=yoffsetloop(delta,N)

h=gca;
if isappdata(gcf,'yoffsetlock')
   if getappdata(gcf,'yoffsetlock')
      h=axeshandles(gcf);
   end
end

for i=1:length(h)
    axes(h(i))
    if isappdata(gca,'lastyoffset')
       lastdelta=getappdata(gca,'lastyoffset');
       yoffsetapply(-1*lastdelta,N);
       setappdata(gca,'lastyoffset',0)
    end
    yoffsetapply(delta,N);
end

axes(h(1))

function[]=yoffsetapply(delta,N)
h=linehandles;

if ischar(delta),delta=str2double(delta);end
yy=get(h,{'ydata'});

if length(delta)==1 && ~isreal(delta)
   ax=axis;
   dy=ax(4)-ax(3);
   delta=dy.*imag(delta)/100;
end

if length(delta)==1
   for i=length(yy):-1:1
       yy{i}=yy{i}+delta*(mod(length(yy)-i,round(length(yy)/N)));
   end
elseif length(delta)==length(h)
   for i=length(yy):-1:1
       yy{i}=yy{i}+delta(i);
   end
else
   error('Number of line handles does not equal offset array length.')
end

set(h,{'ydata'},yy);
setappdata(gca,'lastyoffset',delta)
setappdata(gca,'lasthandles',h)

function[b]=plotmodified(h)
b=0;
if isappdata(gca,'lasthandles')
  b=1;
  h2=getappdata(gca,'lasthandles');
  if length(h2)==length(h)
      if all(h2==h)
          b=0;
      end
  end
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function imo = cnn_get_batch(images, varargin)
% Modified from CNN_IMAGENET_GET_BATCH
%
% - added `pad` option
% - deals with images of types other than jpeg

opts.imageSize = [227, 227] ;
opts.border = [29, 29] ;
opts.pad = 0;  % [TOP BOTTOM LEFT RIGHT]
opts.keepAspect = true ;
opts.numAugments = 1 ;
opts.transformation = 'none' ;
opts.averageImage = [] ;
opts.rgbVariance = zeros(0,3,'single') ;
opts.interpolation = 'bilinear' ;
opts.numThreads = 1 ;
opts.prefetch = false ;
opts = vl_argparse(opts, varargin);

% if only one value is given, apply the same amount of padding to all borders
if numel(opts.pad)==1, opts.pad = repmat(opts.pad,[1 4]); end
if numel(opts.border)==1, opts.border = repmat(opts.border,[1 2]); end

% fetch is true if images is a list of filenames (instead of
% a cell array of images)
fetch = numel(images) >= 1 && ischar(images{1}) ;

% isjpg is true if all images to fetch are of jpeg format
isjpg = fetch && strcmpi(images{1}(end-3:end),'.jpg');

if opts.prefetch
  if isjpg, vl_imreadjpeg(images, 'numThreads', opts.numThreads, 'prefetch'); end
  imo = [] ;
  return ;
end

if fetch
  if isjpg,
    im = vl_imreadjpeg(images,'numThreads', opts.numThreads) ;
  else
    im = cell(size(images));
  end
else
  im = images ;
end

tfs = [] ;
switch opts.transformation
  case 'none'
    tfs = [
      .5 ;
      .5 ;
       0 ] ;
  case 'f5'
    tfs = [...
      .5 0 0 1 1 .5 0 0 1 1 ;
      .5 0 1 0 1 .5 0 1 0 1 ;
       0 0 0 0 0  1 1 1 1 1] ;
  case 'f25'
    [tx,ty] = meshgrid(linspace(0,1,5)) ;
    tfs = [tx(:)' ; ty(:)' ; zeros(1,numel(tx))] ;
    tfs_ = tfs ;
    tfs_(3,:) = 1 ;
    tfs = [tfs,tfs_] ;
  case 'stretch'
  otherwise
    error('Uknown transformations %s', opts.transformation) ;
end
[~,transformations] = sort(rand(size(tfs,2), numel(images)), 1) ;

if ~isempty(opts.rgbVariance) && isempty(opts.averageImage)
  opts.averageImage = zeros(1,1,3) ;
end
if numel(opts.averageImage) == 3
  opts.averageImage = reshape(opts.averageImage, 1,1,3) ;
end

imo = zeros(opts.imageSize(1), opts.imageSize(2), 3, ...
            numel(images)*opts.numAugments, 'single') ;

si = 1 ;
for i=1:numel(images)

  % acquire image
  if isempty(im{i})
    imt = imread(images{i}) ;
    imt = single(imt) ; % faster than im2single (and multiplies by 255)
  else
    imt = im{i} ;
  end
  if size(imt,3) == 1
    imt = cat(3, imt, imt, imt) ;
  end

  % pad
  if ~isempty(opts.pad) && any(opts.pad>0),
    imtt = imt;
    imt = 255*ones(size(imtt,1)+sum(opts.pad(1:2)), ...
      size(imtt,2)+sum(opts.pad(3:4)), 3, 'like', imtt);
    imt(opts.pad(1)+(1:size(imtt,1)), opts.pad(3)+(1:size(imtt,2)),:) = imtt;
  end

  % resize
  w = size(imt,2) ;
  h = size(imt,1) ;
  factor = [(opts.imageSize(1)+opts.border(1))/h ...
            (opts.imageSize(2)+opts.border(2))/w];

  if opts.keepAspect
    factor = max(factor) ;
  end
  if any(abs(factor - 1) > 0.0001)
    imt = imresize(imt, ...
                   'scale', factor, ...
                   'method', opts.interpolation) ;
  end

  % crop & flip
  w = size(imt,2) ;
  h = size(imt,1) ;
  for ai = 1:opts.numAugments
    switch opts.transformation
      case 'stretch'
        sz = round(min(opts.imageSize(1:2)' .* (1-0.1+0.2*rand(2,1)), [h;w])) ;
        dx = randi(w - sz(2) + 1, 1) ;
        dy = randi(h - sz(1) + 1, 1) ;
        flip = rand > 0.5 ;
      otherwise
        tf = tfs(:, transformations(mod(ai-1, numel(transformations)) + 1)) ;
        sz = opts.imageSize(1:2) ;
        dx = floor((w - sz(2)) * tf(2)) + 1 ;
        dy = floor((h - sz(1)) * tf(1)) + 1 ;
        flip = tf(3) ;
    end
    sx = round(linspace(dx, sz(2)+dx-1, opts.imageSize(2))) ;
    sy = round(linspace(dy, sz(1)+dy-1, opts.imageSize(1))) ;
    if flip, sx = fliplr(sx) ; end

    if ~isempty(opts.averageImage)
      offset = opts.averageImage ;
      if ~isempty(opts.rgbVariance)
        offset = bsxfun(@plus, offset, reshape(opts.rgbVariance * randn(3,1), 1,1,3)) ;
      end
      imo(:,:,:,si) = bsxfun(@minus, imt(sy,sx,:), offset) ;
    else
      imo(:,:,:,si) = imt(sy,sx,:) ;
    end
    si = si + 1 ;
  end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [UX, DX, VX] = fadi( A, B, M, N, p, q )
%FADI   Factored alternating direction implicit method.
% X = FADI( A, B, M, N, p, q ) solves the Sylvester equation
%
%            A*X - X*B = M*N.'
%
% using the FADI method with ADI shifts p and q. The righthand side must be
% given in low-rank form, i.e., M and N where rhs = M*N.'.

%% Reference:
%
% [1] Benner, Peter, Ren-Cang Li, and Ninoslav Truhar.
% "On the ADI method for Sylvester equations." J. of Comp. and App. Math.
% 233.4 (2009): 1035-1045.

[m, rho] = size( M );
n = size(N, 1);
RankOfSolution = rho * numel(p);
UX = zeros(m, RankOfSolution);
VX = zeros(n, RankOfSolution);
DX = diag( kron(q-p, ones(1,rho)) );
Im = speye(m);
In = speye(n);
UX(:, 1:rho) = (A+p(1)*Im)\M;
VX(:, 1:rho) = (B+q(1)*In)\N;
for j = 1:numel(p)-1
    UX(:, j*rho+(1:rho)) = (A+q(j)*Im)*((A+p(j+1)*Im)\UX(:,(j-1)*rho+(1:rho)));
    VX(:, j*rho+(1:rho)) = (B+p(j)*In)*((B+q(j+1)*In)\VX(:,(j-1)*rho+(1:rho)));
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function imdb = dtd_get_database(dtdDir, varargin)
opts.seed = 1 ;
opts.joint = 0;
opts = vl_argparse(opts, varargin) ;

rng(opts.seed) ;

imdb.imageDir = fullfile(dtdDir, 'images') ;
imdb.maskDir = fullfile(dtdDir, 'mask') ;

cats = dir(imdb.imageDir) ;
cats = cats([cats.isdir] & ~ismember({cats.name}, {'.','..'})) ;
imdb.classes.name = {cats.name} ;
imdb.images.id = [] ;
imdb.sets = {'train', 'val', 'test'} ;

for c=1:numel(cats)
  ims = dir(fullfile(imdb.imageDir, imdb.classes.name{c}, '*.jpg'));
  imdb.images.name{c} = cellfun(@(S) fullfile(imdb.classes.name{c}, S), ...
    {ims.name}, 'Uniform', 0);
  imdb.images.label{c} = c * ones(1,numel(ims)) ;
  if numel(ims) ~= 120, error('DTD folder appears inconsistent') ; end
end
imdb.images.name = horzcat(imdb.images.name{:}) ;
imdb.images.label = horzcat(imdb.images.label{:}) ;
imdb.images.id = 1:numel(imdb.images.name) ;

for s = 1:3
  list = textread(fullfile(dtdDir, 'labels', ...
    sprintf('%s%d.txt', imdb.sets{s}, opts.seed)),'%s') ;
  ok = ismember(imdb.images.name, list) ;
  imdb.images.set(find(ok)) = s ;
end
if any(~ismember(imdb.images.set, 1:3)), error('DTD appears inconsistent') ; end

imdb.segments = imdb.images ;
imdb.segments.imageId = imdb.images.id ;
imdb.segments.mask = strrep(imdb.images.name, 'image', 'mask') ;

imdb.meta.classes = imdb.classes.name ;
imdb.meta.inUse = true(1, numel(imdb.meta.classes)) ;
imdb.segments.difficult = false(1, numel(imdb.segments.id)) ;

if (opts.joint)
    annoPath = fullfile(dtdDir, 'labels', 'labels_joint_anno.txt');
    imdb.images.label = -1 * ones(numel(imdb.meta.classes), numel(imdb.images.name));
    fid = fopen(annoPath);
    if (fid > 0)
        lines = textscan(fid, '%s', 'Delimiter', '\n');
        lines = lines{1};
        for ii = 1 : numel(lines)
            parts = strsplit(lines{ii}, ' ');
            parts = parts(1: end - 1);
            imName = parts{1};
            parts = parts{2 : end};
            keyAttr = imName(1 : find(imName == '/') - 1);
            imId = strcmp(imdb.images.name, imName);
            imdb.images.label(ismember(imdb.meta.classes, parts), imId) = 1;
            imdb.images.label(ismember(imdb.meta.classes, keyAttr), imId) = 1;
        end
    end
    imdb.segments.label = imdb.images.label;
end

author: mcimpoi

repo: deep-fbanks

save_path: github-repos/MATLAB/mcimpoi-deep-fbanks


text:
function ViewClusters_Callback(hObject, eventdata, handles)
[ClusteringData,clustAssign] = CreateClusteringData(handles, 'forClustering', false);

[~, clusterName, rejected, finished, clustAssign] = clusteringGUI(clustAssign, ClusteringData);

% Save the clusters
if finished == 1
    saveChoice =  questdlg('Update files with new clusters?','Save clusters','Yes','No','No');
    switch saveChoice
        case 'Yes'
            [~, ~, clustAssign] = unique(clustAssign);
            UpdateCluster(ClusteringData, clustAssign, clusterName, rejected)
        case 'No'
            return
    end
end

end

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak