text:
function [mu,sig,alpha,beta] = fit_eeg_distribution(X,min_clean_fraction,max_dropout_fraction,quants,step_sizes,beta)
% Estimate the mean and standard deviation of clean EEG from contaminated data.
% [Mu,Sigma,Alpha,Beta] = fit_eeg_distribution(X,MinCleanFraction,MaxDropoutFraction,FitQuantiles,StepSizes,ShapeRange)
%
% This function estimates the mean and standard deviation of clean EEG from a sample of amplitude
% values (that have preferably been computed over short windows) that may include a large fraction
% of contaminated samples. The clean EEG is assumed to represent a generalized Gaussian component in
% a mixture with near-arbitrary artifact components. By default, at least 25% (MinCleanFraction) of
% the data must be clean EEG, and the rest can be contaminated. No more than 10%
% (MaxDropoutFraction) of the data is allowed to come from contaminations that cause lower-than-EEG
% amplitudes (e.g., sensor unplugged). There are no restrictions on artifacts causing
% larger-than-EEG amplitudes, i.e., virtually anything is handled (with the exception of a very
% unlikely type of distribution that combines with the clean EEG samples into a larger symmetric
% generalized Gaussian peak and thereby "fools" the estimator). The default parameters should be
% fine for a wide range of settings but may be adapted to accomodate special circumstances.
%
% The method works by fitting a truncated generalized Gaussian whose parameters are constrained by
% MinCleanFraction, MaxDropoutFraction, FitQuantiles, and ShapeRange. The alpha and beta parameters
% of the gen. Gaussian are also returned. The fit is performed by a grid search that always finds a
% close-to-optimal solution if the above assumptions are fulfilled.
%
% In:
%   X : vector of amplitude values of EEG, possible containing artifacts
%       (coming from single samples or windowed averages)
%
%   MinCleanFraction : Minimum fraction of values in X that needs to be clean
%                      (default: 0.25)
%
%   MaxDropoutFraction : Maximum fraction of values in X that can be subject to
%                        signal dropouts (e.g., sensor unplugged) (default: 0.1)
%
%   FitQuantiles : Quantile range [lower,upper] of the truncated generalized Gaussian distribution
%                  that shall be fit to the EEG contents (default: [0.022 0.6])
%
%   StepSizes : Step size of the grid search; the first value is the stepping of the lower bound
%               (which essentially steps over any dropout samples), and the second value
%               is the stepping over possible scales (i.e., clean-data quantiles)
%               (default: [0.01 0.01])
%
%   ShapeRange : Range that the clean EEG distribution's shape parameter beta may take (default:
%                1.7:0.15:3.5)
%
% Out:
%   Mu : estimated mean of the clean EEG distribution
%
%   Sigma : estimated standard deviation of the clean EEG distribution
%
%   Alpha : estimated scale parameter of the generalized Gaussian clean EEG distribution (optional)
%
%   Beta : estimated shape parameter of the generalized Gaussian clean EEG distribution (optional)
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2013-08-15

% Copyright (C) Christian Kothe, SCCN, 2013, christiankothe@gmail.com
%
% This program is free software; you can redistribute it and/or modify it under the terms of the GNU
% General Public License as published by the Free Software Foundation; either version 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

% assign defaults
if ~exist('min_clean_fraction','var') || isempty(min_clean_fraction)
    min_clean_fraction = 0.25; end
if ~exist('max_dropout_fraction','var') || isempty(max_dropout_fraction)
    max_dropout_fraction = 0.1; end
if ~exist('quants','var') || isempty(quants)
    quants = [0.022 0.6]; end
if ~exist('step_sizes','var') || isempty(step_sizes)
    step_sizes = [0.01 0.01]; end
if ~exist('beta','var') || isempty(beta)
    beta = 1.7:0.15:3.5; end

% sanity checks
if ~isvector(quants) || numel(quants) > 2
    error('Fit quantiles needs to be a 2-element vector (support for matrices deprecated).'); end
if any(quants(:)<0) || any(quants(:)>1)
    error('Unreasonable fit quantiles.'); end
if any(step_sizes<0.0001) || any(step_sizes>0.1)
    error('Unreasonable step sizes.'); end
if any(beta>=7) || any(beta<=1)
    error('Unreasonable shape range.'); end

% sort data so we can access quantiles directly
X = double(sort(X(:)));
n = length(X);

% calc z bounds for the truncated standard generalized Gaussian pdf and pdf rescaler
for b=1:length(beta)
    zbounds{b} = sign(quants-1/2).*gammaincinv(sign(quants-1/2).*(2*quants-1),1/beta(b)).^(1/beta(b)); %#ok<*AGROW>
    rescale(b) = beta(b)/(2*gamma(1/beta(b)));
end

% determine the quantile-dependent limits for the grid search
lower_min = min(quants);                    % we can generally skip the tail below the lower quantile
max_width = diff(quants);                   % maximum width is the fit interval if all data is clean
min_width = min_clean_fraction*max_width;   % minimum width of the fit interval, as fraction of data

% get matrix of shifted data ranges
X = X(bsxfun(@plus,(1:round(n*max_width))',round(n*(lower_min:step_sizes(1):lower_min+max_dropout_fraction))));
X1 = X(1,:); X = bsxfun(@minus,X,X1);

opt_val = Inf;
% for each interval width...
for m = round(n*(max_width:-step_sizes(2):min_width))
    % scale and bin the data in the intervals
    nbins = round(3*log2(1+m/2));
    H = bsxfun(@times,X(1:m,:),nbins./X(m,:));
    logq = log(histc(H,[0:nbins-1,Inf]) + 0.01);

    % for each shape value...
    for b=1:length(beta)
        bounds = zbounds{b};
        % evaluate truncated generalized Gaussian pdf at bin centers
        x = bounds(1)+(0.5:(nbins-0.5))/nbins*diff(bounds);
        p = exp(-abs(x).^beta(b))*rescale(b); p=p'/sum(p);

        % calc KL divergences
        kl = sum(bsxfun(@times,p,bsxfun(@minus,log(p),logq(1:end-1,:)))) + log(m);

        % update optimal parameters
        [min_val,idx] = min(kl);
        if min_val < opt_val
            opt_val = min_val;
            opt_beta = beta(b);
            opt_bounds = bounds;
            opt_lu = [X1(idx) X1(idx)+X(m,idx)];
        end
    end
end

% recover distribution parameters at optimum
alpha = (opt_lu(2)-opt_lu(1))/diff(opt_bounds);
mu = opt_lu(1)-opt_bounds(1)*alpha;
beta = opt_beta;

% calculate the distribution's standard deviation from alpha and beta
sig = sqrt((alpha^2)*gamma(3/beta)/gamma(1/beta));

repo: REST


text:
function A=polyarea_signed(V)

% function A=polyarea_signed(V)
%-------------------------------------------------------------------------
%
%
%
% Background (https://demonstrations.wolfram.com/SignedAreaOfAPolygon/):
% The formula for the area of a simple polygon can be elegantly derived
% using Green's theorem and extended to moments of the region. S. F.
% Bockman, "Generalizing the Formula for Areas of Polygons to Moments,"
% Amer. Math. Monthly, 96(2), 1989 pp. 131-132.
%
%-------------------------------------------------------------------------
%%

E=[(1:size(V,1))' [(2:size(V,1))'; 1]]; %Edge array
X=V(:,1); %X coordinates
Y=V(:,2); %Y coordinates
XE=X(E);  %X coordinates of edge points
YE=Y(E); %Y coordinates of edge points
A=sum(0.5*((XE(:,1).*YE(:,2))-(XE(:,2).*YE(:,1)))); %Signed area

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function setLastLineStyle(h, line_style)

% setLastLineStyle get all the lines handler and change them to a <linestyle> value
%
% SINTAX:
%   setLastLineStyle(<h>, line_style)
%   setLastLineStyle(line_style)
%
% EXAMPLE:
%   setLastLineStyle(gcf, ':');
%   setLastLineStyle(':');
%
% INPUT:
%   h          = handler to the figure to modify        <optional argument>
%   line_style = requested line style for the line
%
% DEFAULT VALUES:
%   h       = gcf
%

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Andrea Gatti
%  Contributors:     ...
%  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
    line_style = h;
    h = gcf;
end
hline = findobj(h, 'type', 'line');
set(hline(1), 'LineStyle', line_style);

author: goGPS-Project


text:
function c = truncCoeffs(f, N)
%TRUNCCOEFFS   Trigonometric least square coefficients of a CLASSICFUN.
%   C = TRUNCCOEFFS(F, N) returns the 'middle' 2*N+1 trigonometric
%   Fourier coefficients of F.

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

% Call TRUNCCOEFFS() of the .ONEFUN:
c = truncCoeffs(f.onefun, N);
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% SYNTAX:
% data_dod = hmrR_MotionCorrectRLOESS(data_dod, span, turnon)
%
% UI NAME:
% Motion_Correct_RLOESS
%
% DESCRIPTION:
%
% INPUTS:
% data_dod: SNIRF data structure containing delta_OD
% span:
% turnon:   Optional argument to enable/disable this function in a processing stream chain
%
% OUTPUTS:
% data_dod: SNIRF data structure containing delta_OD after motion correction,
%           same size as dod (Channels that are not in the active ml remain unchanged)
%
% USAGE OPTIONS:
% Motion_Correct_RLOESS: dod = hmrR_MotionCorrectRLOESS(dod, span, turnon)
%
% PARAMETERS:
% span: 0.02
% turnon: 1
%
% PREREQUISITES:
% Intensity_to_Delta_OD: dod = hmrR_Intensity2OD( intensity )
%
% LOG:
%
function data_dod = hmrR_MotionCorrectRLOESS(data_dod, span, turnon)

% span = 0.02 (default)
if span<0
    return
end

% Meryem Yucel, Oct, 2017
% Added turn on/off option Meryem Nov 2017
if ~exist('turnon','var')
    turnon = 1;
end
if turnon==0
    return;
end

for iBlk=1:length(data_dod)
    dod = data_dod(iBlk).GetDataTimeSeries();
    t   = data_dod(iBlk).GetTime();
    for i=1:size(dod,2)
        dod(:,i) = smooth(t, dod(:,i), span, 'rloess');
    end
    data_dod(iBlk).SetDataTimeSeries(dod);
end

author: BUNPC


text:
clear all
B0 = 1; Rtriple = 2; Damnio = 3;
B1 = 4; Ramnio = 5; Dabort = 6;
B2 = 7; U = 8;

N = 8;
dag = zeros(N,N);
dag(B0, [Rtriple B1 Ramnio]) = 1;
dag(Rtriple, [Damnio Dabort]) = 1;
dag(Damnio, [B1 Ramnio]) = 1;
dag(B1, B2) = 1;
dag(Ramnio, [Dabort U]) = 1;
dag(Dabort, B2) = 1;
dag(B2, U) = 1;

ns = zeros(1,N);
ns(B0) = 2;
ns(B1) = 3;
ns(B2) = 4;
ns(Rtriple) = 2;
ns(Ramnio) = 3;
ns(Damnio) = 2;
ns(Dabort) = 2;
ns(U) = 1;

limid = mk_limid(dag, ns, 'chance', [B0 B1 B2], ...
		 'decision', [Damnio Dabort], 'utility', [U]);

% states of nature
healthy = 1; downs = 2; miscarry = 3; aborted = 4;
% test results
pos = 1; neg = 2; unk = 3;
% actions
yes = 1; no = 2;

% Prior probability baby has downs syndrome
tbl = zeros(2,1);
p = 1/1000; % from www.downs-syndrome.org.uk figure
p = 24/10000; % www-personal.umich.edu/~bobwolfe/560/review/Downs.pdf (for women agen 35-40)
tbl(healthy) = 1-p;
tbl(downs) = p;
limid.CPD{B0} = tabular_CPD(limid, B0, tbl);

% Reliability of triple screen test
% Unreliable sensor
% B0 -> Rtriple
tbl = zeros(2,2); % Rtriple = pos, neg
p = 0.5;  % high false positive rate (guess)
tbl(healthy, :) = [p 1-p];
p = 0.6; % low detection rate (march of dimes figure)
tbl(downs, :) = [p 1-p];
limid.CPD{Rtriple} = tabular_CPD(limid, Rtriple, tbl);

limid.CPD{Damnio} = tabular_decision_node(limid, Damnio);

% Effect of amnio on baby  B0,Damnio -> B1
 % 1/200 risk of miscarry
p = 1/200; % (march of dimes figure)
tbl = zeros(2, 2, 3); % B1 = healthy, downs, miscarry
tbl(healthy, no, :) =  [1     0     0];
tbl(downs, no, :) =    [0     1     0];
tbl(healthy, yes, :) = [1-p     0   p];
tbl(downs, yes, :) =   [0     1-p   p];
limid.CPD{B1} = tabular_CPD(limid, B1, tbl);

% Reliability of amnio  B0, Damnio -> Ramnio
% Perfect sensor
tbl = zeros(2,2,3); % Ramnio = pos, neg, unk
tbl(:, no, :) =        repmat([0 0 1], 2 ,1);
tbl(healthy, yes, :) = [0 1 0];
tbl(downs, yes, :) =   [1 0 0];
limid.CPD{Ramnio} = tabular_CPD(limid, Ramnio, tbl);

limid.CPD{Dabort} = tabular_decision_node(limid, Dabort);

% Effect of abortion on baby  B1, Dabort -> B2
tbl = zeros(3, 2, 4); % B2 = healthy, downs, miscarry, aborted
tbl(:, yes, :) =       repmat([0 0 0 1], 3, 1);
tbl(healthy, no, :) =  [1 0 0 0];
tbl(downs, no, :) =    [0 1 0 0];
tbl(miscarry, no, :) = [0 0 1 0];
limid.CPD{B2} = tabular_CPD(limid, B2, tbl);

% Utility U(Ramnio, B2)
tbl = zeros(3, 4);
tbl(:, healthy) = 5000;
tbl(:, downs) = -50000;
tbl(:, miscarry) = -1000;
tbl(:, aborted) = -1000;

if 0
%tbl(unk, miscarry) = 0; % this case is impossible
tbl(pos, miscarry) = -1;
tbl(neg, miscarry) = -1000;
if 1
  tbl(unk, aborted) = -100;
  tbl(pos, aborted) = -1;
  tbl(neg, aborted) = -500;
else % pro-life utility fn
  tbl(unk, aborted) = -500000;
  tbl(pos, aborted) = -500000;
  tbl(neg, aborted) = -500000;
end
end

limid.CPD{U} = tabular_utility_node(limid, U,  tbl);

engine = jtree_limid_inf_engine(limid);
[strategy, MEU] = solve_limid(engine);

% Rtriple U(Damnio=1=yes)  U(Damnio=2=no)
% 1=pos    0               1
% 2=neg    0               1
dispcpt(strategy{Damnio})
if isequal(strategy{Damnio}(1,:), strategy{Damnio}(2,:))
  % Rtriple result irrelevant
  doAmnio = argmax(strategy{Damnio}(1,:))
else
  doAmnio = 1;
end

% Rtriple Ramnio U(Dabort=yes=1) U(Dabort=no=2)
% 1=pos   1=pos  1               0
% 2=neg   1=pos  1               0
% 1=pos   2=neg  0               1
% 2=neg   2=neg  0               1
% 1=pos   3=unk  0               1
% 2=neg   3=unk  0               1
dispcpt(strategy{Dabort})

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [ Wopt ] =update_beam_v2( H,K,M,grt,Pt,beta,omega )
%     Hc=H*W_old;
%     He=abs(Hc).^2;
    %% update W
    A=zeros(M,M);
    for k0=1:K
        tmp=abs(beta(k0))^2*H(k0,:)'*H(k0,:);
        A=A+tmp;
    end
    %% init
    lambda_min=0;
    lambda_max=real(sum(sum(A)));
    flag=0;
    while(1)
        [Wn]=downbeam_lambda(A,H,K,M,grt,beta,lambda_max);
        power=power_W( Wn );
        if power>Pt
            lambda_min=lambda_max;
            lambda_max=lambda_max*2;
        elseif power==Pt
            flag=1;
            break
        else
            break
        end
    end
    if flag~=1
        rho=Pt/power;
        Wopt=Wn.*sqrt(rho);
        [ ~,~,f0 ] = update_SINR( H,Wopt,K,omega );
        while(1)
%             [ lambda ] = update_lambda( lambda_min,lambda_max,A,M );
%             lambda=(lambda_min+lambda_max)/2;
            [ lambda,lambda_min ] = update_lambda_v2( lambda_min,lambda_max,A,M );
            [Wn]=downbeam_lambda(A,H,K,M,grt,beta,lambda);
            power=power_W( Wn );
            rho=Pt/power;
            Wn=Wn.*sqrt(rho);
            [ ~,~,f1 ] = update_SINR( H,Wn,K,omega );
            %%
            if power>Pt
                lambda_min=lambda;
            else
                lambda_max=lambda;
            end
            if f1>f0
                Wopt=Wn;
            end
            f0=f1;
            if abs(lambda_max-lambda_min)<1e-3 && abs(power-Pt)<1e-5 && abs(f1-f0)<1e-5
                break
            elseif abs(lambda_max-lambda_min)<1e-7 && lambda_min==0
                break
            end
        end
    end
end

author: guohuayan


text:
function chopDoseOutsideStruct(doseNum,structureNum)
%function chopDoseOutsideStruct(doseNum,structureNum)
%
%This function creates a new dose which is same as doseNum, but which lies
%within structureNum.
%
%Usage: planC = chopDoseOutsideStruct(1,6)
%
%APA, 04/01/2010
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

global stateS planC
indexS = planC{end};

%Get associated scan number
scanNum = getStructureAssociatedScan(structureNum);
assocScanUID = planC{indexS.dose}(doseNum).assocScanUID;
scanNumDose = getAssociatedScan(assocScanUID);
tmDose = getTransM('dose',doseNum,planC);
if isempty(tmDose)
    tmDose = eye(4);
end
tmScan = getTransM('scan',scanNum,planC);
if isempty(tmScan)
    tmScan = eye(4);
end
if isempty(scanNumDose) && ~isequal(tmDose,tmScan)
    error('This function currently supports dose and structure with same transformation matrix')
end

%Get dose grid
[xDoseVals, yDoseVals, zDoseVals] = getDoseXYZVals(planC{indexS.dose}(doseNum));

%Get Uniformized structure and grid
[xUnifVals, yUnifVals, zUnifVals] = getUniformScanXYZVals(planC{indexS.scan}(scanNum));
structureMask3M = getUniformStr(structureNum);

[xDoseValsM, yDoseValsM, zDoseValsM] = meshgrid(xDoseVals, yDoseVals, zDoseVals);
xDoseValsV = xDoseValsM(:);
yDoseValsV = yDoseValsM(:);
zDoseValsV = zDoseValsM(:);

%Interpolate structure on to dose grid
structureOnDoseV = interp3(xUnifVals, yUnifVals, zUnifVals, single(structureMask3M), xDoseValsV, yDoseValsV, zDoseValsV, 'nearest',0);

structureOnDoseM = reshape(structureOnDoseV,[length(yDoseVals), length(xDoseVals), length(zDoseVals)]);

planC{indexS.dose}(end+1) = planC{indexS.dose}(doseNum);
planC{indexS.dose}(end).doseArray = planC{indexS.dose}(end).doseArray .* structureOnDoseM;
planC{indexS.dose}(end).doseUID = createUID('dose');

stateS.doseSet = length(planC{indexS.dose});

sliceCallBack('refresh');

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [x cost info] = conjugategradient(problem, x, options)
% Conjugate gradient minimization algorithm for Manopt.
%
% function [x cost info] = conjugategradient(problem)
% function [x cost info] = conjugategradient(problem, x0)
% function [x cost info] = conjugategradient(problem, x0, options)
% function [x cost info] = conjugategradient(problem, [], options)
%
% Apply the conjugate gradient minimization algorithm to the problem
% defined in the problem structure, starting at x0 if it is provided
% (otherwise, at a random point on the manifold). To specify options whilst
% not specifying an initial guess, give x0 as [] (the empty matrix).
%
% In most of the examples bundled with the toolbox (see link below), the
% solver can be replaced by the present one if need be.
%
% The outputs x and cost are the best reached point on the manifold and its
% cost. The struct-array info contains information about the iterations:
%   iter : the iteration number (0 for the initial guess)
%   cost : cost value
%   time : elapsed time in seconds
%   gradnorm : Riemannian norm of the gradient
%   stepsize : norm of the last tangent vector retracted
%   beta : value of the beta parameter (see options.beta_type)
%   linesearch : information logged by options.linesearch
%   And possibly additional information logged by options.statsfun.
% For example, type [info.gradnorm] to obtain a vector of the successive
% gradient norms reached.
%
% The options structure is used to overwrite the default values. All
% options have a default value and are hence optional. To force an option
% value, pass an options structure with a field options.optionname, where
% optionname is one of the following and the default value is indicated
% between parentheses:
%
%   tolgradnorm (1e-6)
%       The algorithm terminates if the norm of the gradient drops below this.
%   maxiter (1000)
%       The algorithm terminates if maxiter iterations have been executed.
%   maxtime (Inf)
%       The algorithm terminates if maxtime seconds elapsed.
%   minstepsize (1e-10)
%       The algorithm terminates if the linesearch returns a displacement
%       vector (to be retracted) smaller in norm than this value.
%   beta_type ('H-S')
%       Conjugate gradient beta rule used to construct the new search
%       direction, based on a linear combination of the previous search
%       direction and the new (preconditioned) gradient. Possible values
%       for this parameter are:
%           'S-D', 'steep' for beta = 0 (preconditioned steepest descent)
%           'F-R' for Fletcher-Reeves's rule
%           'P-R' for Polak-Ribiere's modified rule
%           'H-S' for Hestenes-Stiefel's modified rule
%           'H-Z' for Hager-Zhang's modified rule
%       See Hager and Zhang 2006, "A survey of nonlinear conjugate gradient
%       methods" for a description of these rules in the Euclidean case and
%       for an explanation of how to adapt them to the preconditioned case.
%       The adaption to the Riemannian case is straightforward: see in code
%       for details. Modified rules take the max between 0 and the computed
%       beta value, which provides automatic restart, except for H-Z which
%       uses a different modification.
%   orth_value (Inf)
%       Following Powell's restart strategy (Math. prog. 1977), restart CG
%       (that is, make a -preconditioned- gradient step) if two successive
%       -preconditioned- gradients are "too" parallel. See for example
%       Hager and Zhang 2006, "A survey of nonlinear conjugate gradient
%       methods", page 12. An infinite value disables this strategy. See in
%       code formula for the specific criterion used.
%   linesearch (@linesearch_adaptive)
%       Function handle to a line search function. The options structure is
%       passed to the line search too, so you can pass it parameters. See
%       each line search's documentation for info. Another available line
%       search in manopt is @linesearch, in /manopt/linesearch/linesearch.m
%   statsfun (none)
%       Function handle to a function that will be called after each
%       iteration to provide the opportunity to log additional statistics.
%       They will be returned in the info struct. See the generic Manopt
%       documentation about solvers for further information.
%   stopfun (none)
%       Function handle to a function that will be called at each iteration
%       to provide the opportunity to specify additional stopping criteria.
%       See the generic Manopt documentation about solvers for further
%       information.
%   verbosity (3)
%       Integer number used to tune the amount of output the algorithm
%       generates during execution (mostly as text in the command window).
%       The higher, the more output. 0 means silent.
%   storedepth (2)
%       Maximum number of different points x of the manifold for which a
%       store structure will be kept in memory in the storedb. If the
%       caching features of Manopt are not used, this is irrelevant. For
%       the CG algorithm, a store depth of 2 should always be sufficient.
%
%
% See also: steepestdescent trustregions manopt/solvers/linesearch manopt/examples

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, Dec. 30, 2012.
% Contributors: Nicolas Boumal
% Change log:
%
%   March 14, 2013, NB:
%       Added preconditioner support : see Section 8 in
%       https://www.math.lsu.edu/~hozhang/papers/cgsurvey.pdf
%
%   Sept. 13, 2013, NB:
%       Now logging beta parameter too.
%
%	Nov. 7, 2013, NB:
%       The search direction is not normalized before it is passed to the
%       linesearch anymore. This way, it is up to the designers of the
%       linesearch to decide whether they want to use the norm of the
%       search direction in their algorithm or not. There are reasons
%       against it, but practical evidence that it may help too, so we
%       allow it. The default linesearch_adaptive used does exploit the
%       norm information. The base linesearch does not. You may select it
%       by setting options.linesearch = @linesearch;
%
%	Nov. 29, 2013, NB:
%       Documentation improved: options are now explicitly described.
%       Removed the Daniel rule for beta: it was not appropriate for
%       preconditioned CG and I could not find a proper reference for it.

% Verify that the problem description is sufficient for the solver.
if ~canGetCost(problem)
    warning('manopt:getCost', ...
        'No cost provided. The algorithm will likely abort.');
end
if ~canGetGradient(problem)
    warning('manopt:getGradient', ...
        'No gradient provided. The algorithm will likely abort.');
end

% Set local defaults here
localdefaults.minstepsize = 1e-10;
localdefaults.maxiter = 1000;
localdefaults.tolgradnorm = 1e-6;
localdefaults.linesearch = @linesearch_adaptive;
localdefaults.storedepth = 2;
% Changed by NB : H-S has the "auto restart" property.
% See Hager-Zhang 2005/2006 survey about CG methods.
% Well, the auto restart comes from the 'max(0, ...)', not so much from the
% reason stated in Hager-Zhang I believe. P-R also has auto restart.
localdefaults.beta_type = 'H-S';
localdefaults.orth_value = Inf; % by BM as suggested in Nocedal and Wright

% Merge global and local defaults, then merge w/ user options, if any.
localdefaults = mergeOptions(getGlobalDefaults(), localdefaults);
if ~exist('options', 'var') || isempty(options)
    options = struct();
end
options = mergeOptions(localdefaults, options);

% for convenience
inner = problem.M.inner;
lincomb = problem.M.lincomb;

% Create a store database
storedb = struct();

timetic = tic();

% If no initial point x is given by the user, generate one at random.
if ~exist('x', 'var') || isempty(x)
    x = problem.M.rand();
end

% Compute objective-related quantities for x
[cost grad storedb] = getCostGrad(problem, x, storedb);
gradnorm = problem.M.norm(x, grad);
[Pgrad storedb] = getPrecon(problem, x, grad, storedb);
gradPgrad = inner(x, grad, Pgrad);

% Iteration counter (at any point, iter is the number of fully executed
% iterations so far)
iter = 0;

% Save stats in a struct array info and preallocate,
% see http://people.csail.mit.edu/jskelly/blog/?x=entry:entry091030-033941
stats = savestats();
info(1) = stats;
info(min(10000, options.maxiter+1)).iter = [];

% Initial linesearch memory
lsmem = [];

if options.verbosity >= 2
    fprintf(' iter\t    cost val\t grad. norm\n');
end

% Compute a first descent direction (not normalized)
desc_dir = lincomb(x, -1, Pgrad);

% Start iterating until stopping criterion triggers
while true

    % Display iteration information
    if options.verbosity >= 2
        fprintf('%5d\t%+.4e\t%.4e\n', iter, cost, gradnorm);
    end

    % Start timing this iteration
    timetic = tic();

    % Run standard stopping criterion checks
    [stop reason] = stoppingcriterion(problem, x, options, info, iter+1);

    % Run specific stopping criterion check
    if ~stop && abs(stats.stepsize) < options.minstepsize
        stop = true;
        reason = 'Last stepsize smaller than minimum allowed.';
    end

    if stop
        if options.verbosity >= 1
            fprintf([reason '\n']);
        end
        break;
    end

    % The line search algorithms require the directional derivative of the
    % cost at the current point x along the search direction.
    df0 = inner(x, grad, desc_dir);

    % If we didn't get a descent direction: restart, i.e., switch to the
    % negative gradient. Equivalent to resetting the CG direction to a
    % steepest descent step, which discards the past information.
    if df0 >= 0

        % Or we switch to the negative gradient direction.
        if options.verbosity >= 3
            fprintf(['Conjugate gradient info: got an ascent direction '...
                     '(df0 = %2e), reset to the (preconditioned) '...
                     'steepest descent direction.\n'], df0);
        end
        % Reset to negative gradient: this discards the CG memory.
        desc_dir = lincomb(x, -1, Pgrad);
        df0 = -gradPgrad;

    end

    % Execute line search
    [stepsize newx storedb lsmem lsstats] = options.linesearch(...
                 problem, x, desc_dir, cost, df0, options, storedb, lsmem);

    % Compute the new cost-related quantities for x
    [newcost newgrad storedb] = getCostGrad(problem, newx, storedb);
    newgradnorm = problem.M.norm(newx, newgrad);
    [Pnewgrad storedb] = getPrecon(problem, x, newgrad, storedb);
    newgradPnewgrad = inner(newx, newgrad, Pnewgrad);

    % Apply the CG scheme to compute the next search direction.
    %
    % This paper https://www.math.lsu.edu/~hozhang/papers/cgsurvey.pdf
	% by Hager and Zhang lists many known beta rules. The rules defined
    % here can be found in that paper (or are provided with additional
    % references), adapted to the Riemannian setting.
	%
    if strcmpi(options.beta_type, 'steep') || ...
       strcmpi(options.beta_type, 'S-D')              % Gradient Descent

        beta = 0;
        desc_dir = lincomb(x, -1, Pnewgrad);

    else

        oldgrad = problem.M.transp(x, newx, grad);
        orth_grads = inner(newx, oldgrad, Pnewgrad)/newgradPnewgrad;

        % Powell's restart strategy (see page 12 of Hager and Zhang's
        % survey on conjugate gradient methods, for example)
        if abs(orth_grads) >= options.orth_value,
            beta = 0;
            desc_dir = lincomb(x, -1, Pnewgrad);

        else % Compute the CG modification

            desc_dir = problem.M.transp(x, newx, desc_dir);

            if strcmp(options.beta_type, 'F-R')  % Fletcher-Reeves
                beta = newgradPnewgrad / gradPgrad;

            elseif strcmp(options.beta_type, 'P-R')  % Polak-Ribiere+
                % vector grad(new) - transported grad(current)
                diff = lincomb(newx, 1, newgrad, -1, oldgrad);
                ip_diff = inner(newx, Pnewgrad, diff);
                beta = ip_diff/gradPgrad;
                beta = max(0, beta);

            elseif strcmp(options.beta_type, 'H-S')  % Hestenes-Stiefel+
                diff = lincomb(newx, 1, newgrad, -1, oldgrad);
                ip_diff = inner(newx, Pnewgrad, diff);
                beta = ip_diff / inner(newx, diff, desc_dir);
                beta = max(0, beta);

            elseif strcmp(options.beta_type, 'H-Z') % Hager-Zhang+
                diff = lincomb(newx, 1, newgrad, -1, oldgrad);
                Poldgrad = problem.M.transp(x, newx, Pgrad);
                Pdiff = lincomb(newx, 1, Pnewgrad, -1, Poldgrad);
                deno = inner(newx, diff, desc_dir);
                numo = inner(newx, diff, Pnewgrad);
                numo = numo - 2*inner(newx, diff, Pdiff)*...
                                       inner(newx, desc_dir, newgrad)/deno;
                beta = numo/deno;

                % Robustness (see Hager-Zhang paper mentioned above)
                desc_dir_norm = problem.M.norm(newx, desc_dir);
                eta_HZ = -1/(desc_dir_norm * min(0.01, gradnorm));
                beta = max(beta,  eta_HZ);

            else
                error(['Unknown options.beta_type. ' ...
                       'Should be steep, S-D, F-R, P-R, H-S or H-Z.']);
            end
            desc_dir = lincomb(newx, -1, Pnewgrad, beta, desc_dir);
        end

    end

    % Make sure we don't use too much memory for the store database.
    storedb = purgeStoredb(storedb, options.storedepth);

    % Update iterate info
    x = newx;
    cost = newcost;
    grad = newgrad;
    Pgrad = Pnewgrad;
    gradnorm = newgradnorm;
    gradPgrad = newgradPnewgrad;

    % iter is the number of iterations we have accomplished.
    iter = iter + 1;

    % Log statistics for freshly executed iteration
    stats = savestats();
    info(iter+1) = stats; %#ok<AGROW>

end

info = info(1:iter+1);

if options.verbosity >= 1
    fprintf('Total time is %f [s] (excludes statsfun)\n', info(end).time);
end

% Routine in charge of collecting the current iteration stats
    function stats = savestats()
        stats.iter = iter;
        stats.cost = cost;
        stats.gradnorm = gradnorm;
        if iter == 0
            stats.stepsize = nan;
            stats.time = toc(timetic);
            stats.linesearch = [];
            stats.beta = 0;
        else
            stats.stepsize = stepsize;
            stats.time = info(iter).time + toc(timetic);
            stats.linesearch = lsstats;
            stats.beta = beta;
        end
        stats = applyStatsfun(problem, x, storedb, options, stats);
    end

end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
% Introduction to DS Method - Series LCR Circuit
% File:  dsintro.m
% 12/13/04
clc;clear;
format short g
% Unit suffixes
uF=1e-6;mH=1e-3;KHz=1e3;
%
% Assign component and source values
%
L1=1*mH;C2=0.253303*uF;R3=5;I1=1;E2=1;Ein=1;
%
% Assign circuit parameters.
%
N=2; % Order of circuit; i.e., number of L's and/or C's.
M=1; % Number of independent inputs.
U=2; % Number of unknown nodes in converted circuit.
Y=2; % Output node.
%
% Form A1 & B2
%
%  V1 V2 eL1 iC2  column order of A1
A1=[0 0 0 1;
   0 1 0 0;
   1 -1 0 0;
   1 0 1 0];
%
%  I1 E2 Ein  column order of B2
B2=[I1 0 0;
   I1*R3 0 0;
   0 E2 0;
   0 0 Ein];
%
P=diag([L1 C2]); % Order L1 C2 corresponds to order of I1 and E2 in B2.
%
% As noted previously in the Word file dsintro.doc, if the arrays A1, B2, and P are
% set up correctly as shown here, the coding from here on never changes from
% circuit to circuit, and is hence guaranteed to be correct.
% The remaining code, starting with V=A1\B2, will always be the same,
% i.e., it is a template or "boilerplate".
%
% Uncomment the following two lines to echo the screen display to a text file qbout.txt
%
%fid=fopen('c:\M_files\qbout.txt','w'); % Use local directory
%diary c:\M_files\qbout.txt; % Use local directory
%
% Solve for superposed dc node voltages, dc voltages across the 1A sources, (replacing
% the L's) and dc currents through the 1V sources (replacing the C's).
%
V=A1\B2
%
% H is the last N rows extracted from V.
%
H=V(U+1:U+N,1:N+M)
%
% Get AB containing A and B arrays
%
AB=P\H
%
% Extract A and B from AB
%
A=AB(1:N,1:N)
B=AB(1:N,N+1:N+M)
%
% Extract D and E from V
%
D=V(Y:Y,1:N)
E=V(Y:Y,N+1:N+M)
%
% Optional:  Get Eigenvalues (poles of circuit transfer function)
%
L=eig(A)/(2*pi) % complex conjugate poles
fo=abs(L(1)) % magnitude of one pole
%
% Dc analysis
%
X=-A\B
Vdc=D*X+E
%
% AC analysis
%
BF=5*KHz; % BF = Beginning Frequency in Hz
LF=15*KHz; % LF = Last Frequency in Hz
NP=101; % NP = Number of points
I=eye(N); % Nth order identity matrix
F=linspace(BF,LF,NP); % Vector of frequency points
%
for i=1:Lit
   s=2*pi*F(i)*j; % j = sqrt(-1)
   cv=D*((s*I-A)\B)+E; % complex value of output
   Vo(i)=abs(cv); % Output magnitude in Volts at node Y = 2.
end
%
% plot Vo
%
h=plot(F/KHz,Vo,'k');
grid on;
set(h,'LineWidth',2);
title('Output at node Y=2');
xlabel('Freq (kHz)');ylabel('Volts');
XT=linspace(BF/KHz,LF/KHz,11);
set(gca,'xtick',XT); % force X-axis tick marks
figure(1);
% If opened above, close the opened file qbout.txt
%diary off
%status=fclose(fid);

author: Sable

repo: mcbench-benchmarks


text:
function [A] = tapas_uniqc_spm_matrix(P, order)
% Return an affine transformation matrix
% copy of spm_matrix from Version 6906 (SPM12) 20-Oct-16
% FORMAT [A] = spm_matrix(P [,order])
% P(1)  - x translation
% P(2)  - y translation
% P(3)  - z translation
% P(4)  - x rotation about - {pitch} (radians)
% P(5)  - y rotation about - {roll}  (radians)
% P(6)  - z rotation about - {yaw}   (radians)
% P(7)  - x scaling
% P(8)  - y scaling
% P(9)  - z scaling
% P(10) - x affine
% P(11) - y affine
% P(12) - z affine
%
% order - application order of transformations [Default: 'T*R*Z*S']
%
% A     - affine transformation matrix
%__________________________________________________________________________
%
% spm_matrix returns a matrix defining an orthogonal linear (translation,
% rotation, scaling or affine) transformation given a vector of
% parameters (P).  By default, the transformations are applied in the
% following order (i.e., the opposite to which they are specified):
%
% 1) shear
% 2) scale (zoom)
% 3) rotation - yaw, roll & pitch
% 4) translation
%
% This order can be changed by calling spm_matrix with a string as a
% second argument. This string may contain any valid MATLAB expression
% that returns a 4x4 matrix after evaluation. The special characters 'S',
% 'Z', 'R', 'T' can be used to reference the transformations 1)-4)
% above. The default order is 'T*R*Z*S', as described above.
%
% SPM uses a PRE-multiplication format i.e. Y = A*X where X and Y are 4 x n
% matrices of n coordinates.
%__________________________________________________________________________
%
% See also: tapas_uniqc_spm_imatrix.m
%__________________________________________________________________________
% Copyright (C) 1994-2011 Wellcome Trust Centre for Neuroimaging

% Karl Friston

%-Special case: translation only
%--------------------------------------------------------------------------
if numel(P) == 3
    A = eye(4);
    A(1:3,4) = P(:);
    return;
end

%-Pad P with 'null' parameters
%--------------------------------------------------------------------------
q  = [0 0 0 0 0 0 1 1 1 0 0 0];
P  = [P q((length(P) + 1):12)];

%-Translation / Rotation / Scale / Shear
%--------------------------------------------------------------------------
T  =   [1   0   0   P(1);
        0   1   0   P(2);
        0   0   1   P(3);
        0   0   0   1];

R1  =  [1   0           0           0;
        0   cos(P(4))   sin(P(4))   0;
        0  -sin(P(4))   cos(P(4))   0;
        0   0           0           1];

R2  =  [cos(P(5))   0   sin(P(5))   0;
        0           1   0           0;
       -sin(P(5))   0   cos(P(5))   0;
        0           0   0           1];

R3  =  [cos(P(6))   sin(P(6))   0   0;
       -sin(P(6))   cos(P(6))   0   0;
        0           0           1   0;
        0           0           0   1];

R   = R1*R2*R3;

Z   =  [P(7)   0       0       0;
        0      P(8)    0       0;
        0      0       P(9)    0;
        0      0       0       1];

S   =  [1      P(10)   P(11)   0;
        0      1       P(12)   0;
        0      0       1       0;
        0      0       0       1];

%-Affine transformation matrix
%--------------------------------------------------------------------------
if nargin < 2
    A = T*R*Z*S;
else
    A = eval(sprintf('%s;', order));
    if ~isnumeric(A) || ~isequal(size(A),[4 4])
        error('tapas:uniqc:SpmMatrix:InvalidOrderExpression', ...
            'Invalid order expression ''%s''.', order);
    end
end

author: translationalneuromodeling

repo: tapas


text:
classdef TestFisheyeCalibrate
    %TestFisheyeCalibrate

    methods (Static)
        function test_1
            %TODO
        end

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

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function inside = checkInside(sR,v,varargin)
% check for points to be inside the spherical region
%
% Syntax
%
%   inside = checkInside(sR,v)
%   inside = checkInside(sR,v,'noAntipodal')
%
% Input
%  sR - @sphericalRegion
%  v  - @vector3d
%
% Output
%  isInside - logical
%
% Options
%  noAntipodal - ignore antipodal symmetry
%

% for antipodal symmetry check v and -v
if (sR.antipodal || check_option(varargin,'antipodal')) && ...
    ~check_option(varargin,'noAntipodal')
  sR.antipodal = false;
  inside = checkInside(sR,v) | checkInside(sR,-v);
  return
end

if check_option(varargin,'noAntipodal')
  v.antipodal = false;
end

% verify all conditions are satisfies
inside = true(size(v));
v = normalize(vector3d(v));
for i = 1:length(sR.N)
  inside = inside & (dot(subSet(sR.N,i),v) >= sR.alpha(i)-1e-4);
end

end

author: mtex-toolbox

repo: mtex

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


text:
% Batch script for sparse prediction
clear all

% --- training method ID
method_id = 8;
% 1: linear_sparse_space
% 2: linear_sparse_space_vec
% 3: linear_sparse_stepwise
% 4: linear_sparse_stepwise_vec
% 5: linear_sparse_stepwise_aa
% 6: linear_sparse_seq
% 7: linear_stepwise_reg_bayes

% --- Basic Learning Parameters
parm.Ntrain = 1000; % # of total training iteration
parm.Nskip  = 100;	% skip steps for display info

% --- Time delay embedding parameter
parm.Tau   = 1 ;   % Lag time steps
parm.Dtau  = 3 ;   % Number of embedding dimension
parm.Tpred = 0 ;   % Prediction time step : y(t+Tpred) = W * x(t)

% --- Normalization parameter
parm.data_norm  = 1;	% Normalize input and output

% File name of old training result for retraining
old_file = [];

% File name for test data
datafile  = ['./test/test.mat'];
modelfile = ['./test/model'];
newdata   = 0; % = 1: make new data

% --- Make or Load training & test data
if newdata == 0 && exist(datafile,'file')
	load(datafile, ...
		'xdata','ydata','xtest','ytest')
else
	% Training & test data setting
	Tdata  = 1000 ; % # of training data
	Ttest  = 1000 ; % # of test data

	parm.Ydim = 4;
	parm.Xdim = 200;   % Input dim
	parm.Meff = 30;    % Effective Input
	parm.Sdim = 10; % number of sinusoidal input
	parm.Tmin = 50; % minimum period
	parm.Tmax = 200;% maximum period
	parm.SY   = [0.3; 0.1; 0.5; 1.0];% output moise variance
	parm.Ntrial = 1;

	parm.sy = parm.SY;
	[xdata,ydata,Wout,Weff,Wrdn,Win,Wirr] = ...
		make_train_data(parm, Tdata);

	parm.sy = 0.0;
	[xtest,ytest] = ...
		make_train_data(parm, Ttest, Wout,Wrdn,Win,Wirr);

	if ~isempty(datafile)
		save(datafile, ...
			'xdata','ydata','xtest','ytest','Wout','Weff','parm')
	end
end

% Normalize input data
[X,nparm] = normalize_data(xdata, parm.data_norm);
parm.xmean = nparm.xmean;
parm.xnorm = nparm.xnorm;

% Normalize output data
[Y,nparm] = normalize_data(ydata, parm.data_norm);
parm.ymean = nparm.xmean;
parm.ynorm = nparm.xnorm;

% Time alignment for prediction using embedding input
[tx,ty] = pred_time_index(X,parm);
X = X(:,tx,:);
Y = Y(:,ty,:);

% --- Initialization of Model Parameters
if ~isempty(old_file)
	% Start from old result
	load([old_file], 'Model')
else
	Model = [];
end

%profile_on = 1;
%profile_start(profile_on);

%
% --- Sparse estimation
%
switch	method_id
case	1
	[Model, Info] = linear_sparse_space(X, Y, Model, parm);
case	2
	[Model, Info] = linear_sparse_space_vec(X, Y, Model, parm);
case	3
	[Model, Info] = linear_sparse_stepwise(X, Y, Model, parm);
case	4
	[Model, Info] = linear_sparse_stepwise_vec(X, Y, Model, parm);
case	5
	[Model, Info] = linear_sparse_stepwise_aa(X, Y, Model, parm);
case	6
	[Model, Info] = linear_sparse_seq(X, Y, Model, parm);
case	7
	[Model, Info] = linear_stepwise_reg_bayes(X, Y, Model, parm);
case	8
	[Model, Info] = linear_fit(X, Y, Model, parm);
case	10
	% --- Delay embedding is done before training
	Xd = delay_embed(X, parm.Tau, parm.Dtau);
	[Model, Info] = linear_map_sparse_cov(Xd, Y, Model, parm);
end

%
% --- Estimate prediction error for test data
%

% Time alignment for prediction using embedding input
[tx,ty] = pred_time_index(xtest,parm);
xtest = xtest(:,tx,:);
ytest = ytest(:,ty,:);

% Use normalization constant calculated by training data
xtest = normalize_data(xtest, parm.data_norm, parm);

% --- Prediction for test data
ypred = predict_output(xtest, Model, parm);
err   = sum((ytest(:)-ypred(:)).^2)/sum(ytest(:).^2)

if ~isempty(modelfile)
	fsave = [modelfile sprintf('_id%d.mat',method_id)];
	save(fsave, 'Model', 'Info', 'parm');
end

plot_predict

%profile_end(profile_on)
return

author: KamitaniLab

repo: GenericObjectDecoding


text:
% test code for tt_qlaplacex_dd()
%
% September 22, 2010
% Vladimir Kazeev
% vladimir.kazeev@gmail.com
% INM RAS
% Moscow, Russia
%
% Look for details in the Preprint No. 75, 2010 of
% Max-Planck Institute for Mathematics in the Sciences
% Vladimir A. Kazeev and Boris N. Khoromskij
% On explicit QTT representation of Laplace operator and its inverse
% http://www.mis.mpg.de/publications/preprints/2010/prepr2010-75.html

% d is the only parameter
d=[3,4,5];
%

D=size(d,2);
tt=tt_qlaplacex_dd(d);
tt=tt_mat_to_vec(tt);

full=nd_to_full(tt);

disp('inv computation... ');
Li=cell(D,1);
for k=1 : D
	L=2*eye(2^d(k));
	for i=1 : 2^d(k)-1
		L(i,i+1)=-1;
		L(i+1,i)=-1;
	end
	Li{k}=inv(L);
end
disp('OK');

disp('kron computation... ');
Z=Li{1};
for k=2 : D
	Z=kron(Z,Li{k});
end
disp('OK');

err=norm(full-Z,'fro');
fprintf('fro err = %e\n', err);

author: oseledets

repo: TT-Toolbox


text:
function dy = diffxy(x,y,varargin)
% DIFFXY - accurate numerical derivative/differentiation of Y w.r.t X.
%
%   DY = DIFFXY(X,Y) returns the derivative of Y with respect to X using a
%        pseudo second-order accurate method. DY has the same size as Y.
%   DY = DIFFXY(X,Y,DIM) returns the derivative along the DIM-th dimension
%        of Y. The default is differentiation along the first
%        non-singleton dimension of Y.
%   DY = DIFFXY(X,Y,DIM,N) returns the N-th derivative of Y w.r.t. X.
%        The default is 1.
%
%   Y may be an array of any dimension.
%   X can be any of the following:
%       - array X with size(X) equal to size(Y)
%       - vector X with length(X) equal to size(Y,DIM)
%       - scalar X denotes the spacing increment
%   DIM and N are both integers, with 1<=DIM<=ndims(Y)
%
%   DIFFXY has been developed especially to handle unequally spaced data,
%   and features accurate treatment for end-points.
%
%   Example:
%   % Data with equal spacing
%     x = linspace(-1,2,20);
%     y = exp(x);
%
%     dy = diffxy(x,y);
%     dy2 = diffxy(x,dy);  % Or, could use >> dy2 = diffxy(x,y,[],2);
%     figure('Color','white')
%     plot(x,(y-dy)./y,'b*',x,(y-dy2)./y,'b^')
%
%     Dy = gradient(y)./gradient(x);
%     Dy2 = gradient(Dy)./gradient(x);
%     hold on
%     plot(x,(y-Dy)./y,'r*',x,(y-Dy2)./y,'r^')
%     title('Relative error in derivative approximation')
%     legend('diffxy: dy/dx','diffxy: d^2y/dx^2',...
%            'gradient: dy/dx','gradient: d^2y/dx^2')
%
%   Example:
%   % Data with unequal spacing.
%     x = 3*sort(rand(20,1))-1;
%     % Run the example above from y = exp(x)
%
%   See also DIFF, GRADIENT
%        and DERIVATIVE on the File Exchange

% for Matlab (should work for most versions)
% version 1.0 (Nov 2010)
% (c) Darren Rowland
% email: darrenjrowland@hotmail.com
%
% Keywords: derivative, differentiation

[h,dy,N,perm] = parse_inputs(x,y,varargin);
if isempty(dy)
    return
end
n = size(h,1);
i1 = 1:n-1;
i2 = 2:n;

for iter = 1:N
    v = diff(dy)./h;
    if n>1
        dy(i2,:) = (h(i1,:).*v(i2,:)+h(i2,:).*v(i1,:))./(h(i1,:)+h(i2,:));
        dy(1,:) = 2*v(1,:) - dy(2,:);
        dy(n+1,:) = 2*v(n,:) - dy(n,:);
    else
        dy(1,:) = v(1,:);
        dy(n+1,:) = dy(1,:);
    end
end

% Un-permute the derivative array to match y
dy = ipermute(dy,perm);

%%% Begin local functions %%%
function [h,dy,N,perm] = parse_inputs(x,y,v)

numvarargs = length(v);
if numvarargs > 2
    error('diffxy:TooManyInputs', ...
        'requires at most 2 optional inputs');
end

h = [];
N = [];
perm = [];

% derivative along first non-singleton dimension by default
dim = find(size(y)>1);
% Return if dim is empty
if isempty(dim)
    dy = [];
    return
end
dim = dim(1);

% Set defaults for optional arguments
optargs = {dim 1};
newVals = ~cellfun('isempty', v);
optargs(newVals) = v(newVals);
[dim, N] = optargs{:};

% Error check on inputs
if dim<1 || dim>ndims(y) || dim~=fix(dim) || ~isreal(dim)
    error('diffxy:InvalidOptionalArg',...
        'dim must be specified as a non-negative integer')
end
if N~=fix(N) || ~isreal(N)
    error('diffxy:InvalidOptionalArg',...
        'N must be an integer')
end

% permutation which will bring the target dimension to the front
perm = 1:length(size(y));
perm(dim) = [];
perm = [dim perm];
dy = permute(y,perm);

if length(x)==1  % Scalar expansion to match size of diff(dy,[],1)
    sizeh = size(dy);
    sizeh(1) = sizeh(1) - 1;
    h = repmat(x,sizeh);
elseif ndims(x)==2 && any(size(x)==1) % Vector x expansion
    if length(x)~=size(dy,1)
        error('diffxy:MismatchedXandY',...
            'length of vector x must match size(y,dim)')
    end
    x = x(:);
    sizeh = size(dy);
    sizeh(1) = 1;
    h = repmat(diff(x),sizeh);
else
    if size(y) ~= size(x)
        error('diffxy:MismatchedXandY',...
            'mismatched sizes of arrays x and y');
    end
    % Permute x as for y, then diff
    h = diff(permute(x,perm),[],1);
end


text:
function Nodes = traverseGrid(Nodes, Edges, verbose)
%% Traverse the color grid graph to find the neighbors of each Node

%% License
% ACADEMIC OR NON-PROFIT ORGANIZATION NONCOMMERCIAL RESEARCH USE ONLY
% Copyright (c) 2018 Bingyao Huang
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:

% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.

% If you publish results obtained using this software, please cite our paper.

% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
% SOFTWARE.

%%

if(nargin < 3)
    verbose = 0;
end

%% create the ajacency matrix and graph from Edges and Nodes
numNodes = numel(Nodes);
numEdges = numel(Edges);

A = zeros(numNodes, numNodes); % for both horizontal and vertical
% hA = zeros(numNodes, numNodes); % for horizontal only

for i = 1:numEdges
    curEdge = Edges(i);

    if (numel(curEdge.nodes) == 2)% ignore end point and end edge
        edgeLength = numel(curEdge.PixelIdxList);
        A(curEdge.nodes(1), curEdge.nodes(2)) = edgeLength;
        A(curEdge.nodes(2), curEdge.nodes(1)) = edgeLength;
%         if(curEdge.isH)
%             hA(curEdge.nodes(1), curEdge.nodes(2)) = edgeLength;
%             hA(curEdge.nodes(2), curEdge.nodes(1)) = edgeLength;
%         end
    end
end

% vA = A - hA; % for vertical only

% create a graph object of ajacency matrix
% then we can use all the functions in http://www.mathworks.com/help/matlab/ref/graph-object.html
G = graph(A);

%% Traverse
[Nodes.neighbor] = deal([]);
[Nodes.N] = deal([-1]);
[Nodes.S] = deal([-1]);
[Nodes.W] = deal([-1]);
[Nodes.E] = deal([-1]);
[Nodes.error] = deal([]);

for i =1:numNodes
    myRow = Nodes(i).Centroid(2);
    myCol = Nodes(i).Centroid(1);
    neighborNodes = neighbors(G, i);
    Nodes(i).neighbor = neighborNodes;

    numNeighbors = numel(neighborNodes);
    count = 0;
    for j = 1:numNeighbors
        curNode = neighborNodes(j);
        row = Nodes(curNode).Centroid(2);
        col = Nodes(curNode).Centroid(1);

        rowDiff = abs(row - myRow);
        colDiff = abs(col - myCol);
        flag = rowDiff >= colDiff;

        if(myRow >= row && flag)
            if(Nodes(i).N < 0)
                Nodes(i).N = curNode;
                count= count+1;
            end
        elseif(myRow <= row && flag)
            if(Nodes(i).S < 0)
                Nodes(i).S = curNode;
                count= count+1;
            end
        elseif(myCol <= col)
            if(Nodes(i).E < 0)
                Nodes(i).E = curNode;
                count= count+1;
            end
        elseif(myCol >= col)
            if(Nodes(i).W < 0)
                Nodes(i).W = curNode;
                count= count+1;
            end
        end
    end

    if(count ~= numNeighbors)
        Nodes(i).error = 1;
    end
end

%% double check if we have unmatched neighbors
for i =1:numNodes
    N = Nodes(i).N;
    S = Nodes(i).S;
    W = Nodes(i).W;
    E = Nodes(i).E;

    if(N > 0)
        if(Nodes(N).S ~= i)
            Nodes(N).error = i;
            Nodes(i).error = N;
        end
    end

    if(S > 0)
        if(Nodes(S).N ~= i)
            Nodes(S).error = i;
            Nodes(i).error = S;
        end
    end

    if(W > 0)
        if(Nodes(W).E ~= i)
            Nodes(W).error = i;
            Nodes(i).error = W;
        end
    end

    if(E > 0)
        if(Nodes(E).W ~= i)
            Nodes(E).error = i;
            Nodes(i).error = E;
        end
    end
end

errorNodes = [Nodes.error];
if(verbose)
    disp(['Get ', num2str(nnz(errorNodes)), ' error nodes during traverse'])
end
end

author: BingyaoHuang

repo: single-shot-pro-cam-calib


text:
function [fits] = tapas_sem_generate_fits(data, samples, model, time)
%% Generate the fits of the data.
%
% Input
%
% Output
%

% aponteeduardo@gmail.com
% copyright (C) 2019
%

CONG = 0;
INCONG = 1;

n = 3;

n = n + 1;
if nargin < n
    t = data.y.t;
    maxt = max(t);
    mint = min(t);

    nv = 300;

    time = linspace(mint * 0.9, maxt * 1.1, nv)';

end

nv = numel(time);

conds = unique(data.u.tt);
nconds = numel(conds);

fits = struct('pro', cell(nconds, 1), 'anti', [], 't', []);

for i = 1:nconds
    fits(i).t = time;

    % Make sure the dimension are correct
    samples = reshape(samples, 1, numel(samples));

    % Iterate over arrays because the interface is cluncky for single times
    for a = [CONG, INCONG]
        llh = zeros(nv, 1);
        for j = 1:nv
            tdata = struct(...
                'y', struct('t', time(j), 'a', a), ...
                'u', struct('tt', conds(i)));
            tllh = model.llh(tdata, samples);
            % Average the log likelihood
            tllh = mean(exp(tllh));
            llh(j) = tllh;
        end
        switch a
            case CONG
                fits(i).pro = llh;
            case INCONG
                fits(i).anti = llh;
        end
    end
end

end


text:
function this = add_dims(this, iDim, varargin)
% Adds additional dimensions to existing dimInfo
%
%   Y = MrDimInfo()
%   Y.add_dims(iDims, varargin)
%
% This is a method of class MrDimInfo.
%
% IN
%   iDim        (vector of) dimension indices to be added (e.g. 3 for 3rd
%               dimension)
%                   or
%               (cell of) strings of dimension names
%   varargin    PropertyName/Value pairs to set parameters of the new
%               dimension
%
% OUT
%
% EXAMPLE
%   add_dims('coil', 'units', 'nil', 'samplingPoints', {1:8})
%
%   See also MrDimInfo set_dims

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2016-06-18
% Copyright (C) 2016 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% nothing to do here...
if isempty(iDim)
    return
end

%% To update a dimension (set_dims), it has to have a label first
% This can be given explicitly or using the default dimLabels of MrDimInfo

% later permutation at specified index, setup here
doPermute = isnumeric(iDim);
if doPermute
    iDimNumeric = iDim;
end

% direct adding of a labeled dimension, e.g., add_dims('coil', ...)
isStringiDimInput = ischar(iDim) || (iscell(iDim) && ischar(iDim{1}));
if isStringiDimInput
    additionalDimLabels = cellstr(iDim);
else
    % add_dims via their index, e.g. add_dims([4:6}, propName/Value...)
    if ~iscell(iDim)
        iDim = num2cell(iDim);
    end

    if isempty(iDim{1})
        % empty set of dimLabels to add, therefore return
        return
    else

        additionalDimLabels = cellfun(@(x) this.get_default_dim_labels(x), iDim, ...
            'UniformOutput', false);
    end
end

% empty defaults here will be overwritten by set_dims
nDimsOld = this.nDims;
nDimsAdditional = numel(additionalDimLabels);
this.dimLabels = [this.dimLabels additionalDimLabels];
this.samplingPoints(nDimsOld+(1:nDimsAdditional)) = {[]};
this.resolutions(nDimsOld+(1:nDimsAdditional)) = 1;
this.samplingWidths(nDimsOld+(1:nDimsAdditional)) = NaN;
this.units(nDimsOld+(1:nDimsAdditional)) = {''};
% set_dims also needed to add samplingPoints (e.g. via Nsamples/resolution)
if nargin > 2
    this.set_dims(additionalDimLabels, varargin{:});
end

% permute here to corresponding index
if doPermute
    iDim = iDimNumeric;
    nDimsNew = nDimsOld + nDimsAdditional;

    % current order of dimensions after appending the dimensions at the end
    dimOrderNew = [setdiff(1:nDimsNew, iDim), iDim];

    % desired permutation to have right dimensions at their place
    invDimOrderNew(dimOrderNew) = 1:length(dimOrderNew); % inverse permutation
    this.permute(invDimOrderNew);
end


text:
function [y, dzdg, dzdb] = vl_nngnorm(x, g, b, varargin)
%VL_NNGNORM CNN group normalization.
%   Y = VL_NNGNORM(X,G,B) applies group normalization
%   to the input X with shape HxWxCxN. Group normalization is defined as:
%
%      Y(i,j,k,t) = G(k',t) * X_HAT(i,j,k,t) + B(k',t)
%
%   where
%      k' = group_idx(k,C,G), where N_G is the number of groups and
%        group_idx(k,C,G) := floor(k / (C/N_G)).
%      X_HAT(i,j,k,t) = (X_HAT(i,j,k,t) - mu(k',t)) / sigma(k',t)
%      mu(k',t) = mean_ijk'' X(i,j,k'',t),
%      sigma2(k',t) = mean_ijk'' (X(i,j,k'',t) - mu(k'',t))^2,
%      sigma(k',t) = sqrt(sigma2(k) + EPSILON)
%        where k'' takes values such that group_idx(k'',C,G) == group_idx(k,C,G)
%
%   VL_NNGNORM(..., 'option', value, ...) takes the following option:
%
%   `numGroups`:: 32
%    The number of groups used to split the channels when computing
%    normalization statistics.
%
%   `epsilon`:: 1e-4
%    A parameter to add stability to the normalization operation.
%
%   Notes: GroupNorm is introduced in the paper:
%      `Group Normalization, Yuxin Wu, Kaiming He,
%      arXiv preprint arXiv:1803.08494 (2018)
%
% Copyright (C) 2018 Samuel Albanie
% All rights reserved.

  opts.numGroups = 32 ;
  opts.epsilon = 1e-4 ;
  [opts,dzdy] = vl_argparsepos(opts, varargin) ;

  bsize = size(x, 4) ;
  expectedSz = [1 1 size(x,3) 1] ;
  sg = size(g) ; sb = size(b) ;
  assert(all(expectedSz(1:numel(sg)) == sg), 'GAINS have unexpected size') ;
  assert(all(expectedSz(1:numel(sb)) == sb), 'BIASES have unexpected size') ;

  szX = size(x) ; % store original shape

  % compute statistics per group for current minibatch and normalize
  x = reshape(x, size(x,1), size(x,2), [], opts.numGroups, bsize) ;

  mu = groupAvg(x) ;
  sigma2 = groupAvg(bsxfun(@minus, x, mu).^ 2) ;
  sigma = sqrt(sigma2 + opts.epsilon) ;
  x_hat = bsxfun(@rdivide, bsxfun(@minus, x, mu), sigma) ;

  if isempty(dzdy)
    x_hat_ = reshape(x_hat, szX) ;
    y = bsxfun(@times, g, x_hat_) ; % apply gain
    y = bsxfun(@plus, y, b) ; % add bias
  else
    dzdy = dzdy{1} ;
    dzdb = chanSum(dzdy) ;
    x_hat_ = reshape(x_hat, szX) ; dzdg = chanSum(x_hat_ .* dzdy) ;
    dzdy = reshape(dzdy, size(x,1), size(x,2), [], opts.numGroups, bsize) ;

    g_ = reshape(g, 1, 1, size(dzdy, 3), []) ;
    dzdx_hat = bsxfun(@times, dzdy, g_) ;
    t1 = bsxfun(@minus, x, mu) ;
    m = prod([size(x,1) size(x,2) size(x,3)]) ;
    dzdsigma = groupSum((-1/2) * dzdx_hat .* bsxfun(@rdivide, t1, sigma.^3)) ;

    dzdmu = groupSum(bsxfun(@rdivide, dzdx_hat, -sigma)) + ...
                bsxfun(@times, dzdsigma, -2 * groupAvg(t1)) ;

    t4 = bsxfun(@rdivide, dzdx_hat, sigma) + ...
         bsxfun(@times, dzdsigma,  (2 / m) * t1) ;
    dzdx = bsxfun(@plus, t4, dzdmu * (1/m)) ;
    y = reshape(dzdx, szX) ;
  end

% ----------------------------------------
function avg = groupAvg(x)
% ----------------------------------------
  avg = mean(mean(mean(x, 1), 2), 3) ;

% -----------------------
function res = groupSum(x)
% -----------------------
  res = sum(sum(sum(x, 1), 2), 3) ;

% -----------------------
function res = chanSum(x)
% -----------------------
  res = sum(sum(sum(x, 1), 2), 4) ;

repo: MFT


text:
%IM_FEATURES Fixed mapping computing image features by the Matlab image toolbox
%
%		F = IM_FEATURES(A,GRAY,FEATURES)
%
% INPUT
%   A        Dataset with binary object images dataset
%   GRAY     Gray-valued images (matched with A, optional)
%   FEATURES Features to be computed
%
% OUTPUT
%   F        Dataset with computed features
%
% This function is a *replacement* for the function IM_MEASURE that is
% based on the DipLib measure.m. This implementation is using the Matlab
% function REGIONPROPS from the image toolbox. Use HELP REGIONPROPS to
% find out which features are exactly supported.
% In each image of the measurement set GRAY the features given in FEATURES
% are measured. In A a segmented version of GRAY has to be supplied.
% When no GRAY is supplied, the binary images in A are used.
%
% Use FEATURES = 'all' for computing all features.
%
% SEE ALSO
% DATASETS, DATAFILES, REGIONPROPS

% Copyright: D.M.J.Tax, D.M.J.Tax@prtools.org.
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function b = im_features(a,gray,features)

if nargin < 3  features = []; end
if nargin < 2  gray = []; end

if nargin < 1 | isempty(a)
   b = prmapping(mfilename,'fixed',{gray,features});
   b = setname(b,'Image features');
elseif isa(a,'prdataset')
   if (nargin < 3 | isempty(features)) & ~isdataset(gray)
      features = gray;
      gray = a;
   end
   if ~strcmp(class(a),class(gray))
      error('Binary and gray images should be both datasets or datafiles')
   end
   fsize = getfeatsize(a);
   if any(getfeatsize(gray) ~= fsize)
      error('Image structures of binary and gray images should be identical')
   end
   if length(fsize) == 2, fsize = [fsize 1]; end
   if size(a,1) ~= size(gray,1)
      error('Same number of binary and gray images expected')
   end
   if isempty(features)
      features = 'all';
   end
   % first generate the correct feature labels, and a structure PROPS
   % that can be used in REGIONPROPS of Matlab:
   im = data2im(a(1,:));
   [props,featlab] = getlabfromregionprops(im,data2im(gray(1,:)),features);
   % RD 27Nov13, compute dat size before
   ftrial = getfeatfromregionprops(im,data2im(gray(1,:)),props);
   out = zeros(size(a,1),numel(ftrial));
   % Run over all images, and compute the features:
   for n = 1:size(a,1)
     out(n,:) = getfeatfromregionprops(data2im(a,n),data2im(gray,n),props);
   end
   % store in a dataset:
   %b = prdataset(out,getlabels(a));
   b = prdataset(out);
   b = setlablist(b,getlablist(a));
   b = setnlab(b,getnlab(a));
   b = setfeatlab(b,featlab);
   b = setprior(b,getprior(a,0));
   b = setname(b,getname(a));

elseif isa(a,'double')
   if (nargin < 3 | isempty(features)) & ~isdataset(gray)
      features = gray;
      gray = a;
   end
   [props,featlab] = getlabfromregionprops(a,gray,features);
   b = getfeatfromregionprops(a,gray,props);
else
   error('Wrong input')
end

return

% [props,featlab] = getlabfromregionprops(im,gray,features)
%
% Generate a structure PROPS and a cell-array FEATLAB that contains the
% features that are extracted from the images. PROPS is used in the
% Matlab function REGIONPROPS.
function [props,featlab] = getlabfromregionprops(im,gray,features)

% compute the properties:
%keyboard
thr = min(im(:)) + (max(im(:))-min(im(:)))/2;
r = regionprops((im>thr),gray,features);
props = fieldnames(r);
% we do it complicated like this, because bwprops could have been 'all',
% and by this we now get all possible (black and white) properties

% remove the things that we are not interested in:
% RD 27Nov13, PixelValues added
rmfeatures = {'PixelList', 'SubarrayIdx', 'ConvexHull','ConvexImage', ...
   'Image', 'FilledImage', 'PixelIdxList', 'Extrema', 'PixelValues'};
props = setdiff(props,rmfeatures);

% generate the feature labels
featlab = {};
nr = 0;
for i=1:length(props)
   nrf = size(getfield(r,props{i}),2);
   if nrf>1 % there are more measurements per single property
      for j=1:nrf
         nr = nr+1;
         featlab{nr} = [props{i},num2str(j)];
      end
   else % one measurement per property
      nr = nr+1;
      featlab{nr} = props{i};
   end
end
featlab = strvcat(featlab);

return

%   f = getfeatfromregionprops(im,gray,props)
%
% Compute the features F from the images IM. Which features are computed
% is defined in PROPS.
function f = getfeatfromregionprops(im,gray,props)

thr = min(im(:)) + (max(im(:))-min(im(:)))/2;
r = regionprops((im>thr),gray,props);
f = [];
for i=1:length(props)
   f = [f getfield(r,props{i})];
end

author: marianux

repo: ecg-kit

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


text:
% This file is part of the project NILM-Eval (https://github.com/beckel/nilm-eval).
% Licence: GPL 2.0 (http://www.gnu.org/licenses/gpl-2.0.html)
% Copyright: ETH Zurich, 2014
% Author: Romano Cicchetti

% Returns a threshold that defines whether
% the specified appliance is in state 'on' or 'off'
function [threshold] = getThresholdDiffOnOff(applianceID)

threshold_vector = [
        15;  % fridge
        15;  % freezer
        500; % microwave
        500; % dishwasher
        15;  % entertainment
        500; % kettle
        500; % stove
        15;  % coffee machine
        500; % washing machine
        300; % dryer
        15;  % lamp
        15;  % PC
        15;  % laptop
        15;  % TV
        15;  % Stereo
        5;   % Tablet
        5;   % Router
        5;   % Illuminated fountain
    ]; % stereo

     threshold = threshold_vector(applianceID,1);
end

author: beckel

repo: nilm-eval

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


text:
function [oz] = yd32oz(yd3)
% Convert volume from cubic yards to US liquid ounces.
% Chad Greene 2012
oz = yd3*25852.675324;


text:
function F = complements(C1,C2)
%COMPLEMENTS Defines complementary constraints
%
%   F = COMPLEMENTS(C1,C2)

F = complements(lmi(C1),lmi(C2));

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function wimg = warpimg(img, p, sz)
% function wimg = warpimg(img, p, sz)
%
%    img(h,w)
%    p(6,n) : mat format
%    sz(th,tw)
%

%% Copyright (C) 2005 Jongwoo Lim and David Ross.
%% All rights reserved.

if (nargin < 3)
    sz = size(img);
end
if (size(p,1) == 1)
    p = p(:);
end
w = sz(2);  h = sz(1);  n = size(p,2);
%[x,y] = meshgrid(1:w, 1:h);
[x,y] = meshgrid([1:w]-w/2+0.5, [1:h]-h/2);
pos = reshape(cat(2, ones(h*w,1),x(:),y(:)) ...
              * [p(1,:) p(2,:); p(3:4,:) p(5:6,:)], [h,w,n,2]);
cn=size(img,3);
wimg=zeros([sz cn]);
for i=1:cn
    wimg(:,:,i) = squeeze(interp2(img(:,:,i), pos(:,:,:,1), pos(:,:,:,2)));
end
wimg(find(isnan(wimg))) = 0;

repo: Context-Aware-CF-Tracking


text:
function r = mtimes(a,b)
%MTIMES       Taylor multiplication  a * b   (same as  a .* b )
%

% written  05/21/09     S.M. Rump
%

  r = a.*b;

repo: OpenMAS


text:
function x = p15_start ( option, nvar )

%*****************************************************************************80
%
%% P15_START returns a starting point for problem 15.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Output, real X(NVAR), the starting point.
%
  x = zeros ( nvar, 1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [stimset,stim_full] = StimulusKeyPreprocessing(frame_turn,i_fish,isplotting)
% [stimset,stim_full] = StimulusKeyPreprocessing(frame_turn,i_fish,'isplotting');

if exist('isplotting','var'),
    isPlotting = true; % plot parsing process for manual inspection!
else
    isPlotting = false;
end
% Manual inputs needed for this function:
% - protocol sequence and naming, as in 'block_raw' and 'stimset' initialization
% - substitution array 'M_range', to substitute original stimulus code to a
% tighter range that is better defined for plotting (both for plotting values and stimulus bar)

%% Load stimulus

st = frame_turn(:,17)'; % from raw file

% cap abnormally high values in raw values
temp = sort(st,'descend');
thr = temp(round(length(st)/100)); % arb. thr % used /10 up till 1/27/16 for all fish except Fish12
st(st>thr) = thr+1; % assign cap value to be thr+1

xv = 0:thr+1; % x-vector for histogram bins

% Manual inspection
% figure;
% hist(st,xv);

% (continue by default) find bins with significantly high counts, i.e.
% whole range/set of raw stimulus keys, store in 'M_range_raw'
% use to manually define standardized 'M_range'
counts = hist(st,xv);
thr_counts = round(length(st)/500);
M_range_raw = xv(counts>thr_counts);
if isPlotting,
    disp(['M_range_raw = ' num2str(M_range_raw)]); % use this to adjust the substitution array, M_range, manually after this inspection!%%%%%%%%%%%%%%%%%%%%%%
end

%% Raw Stimlus Key ----------------------------------------------------------------------------- MANUAL INPUT!

% Manually input protocol sequence
% stored in 'block', based on the raw stimulus code
%% Manual: params for each fish

% write out protocol sequence for each block; write stimulus sets within block in arrays within cell
if i_fish == 8,
    nBlocks = 3;
    block_raw = cell(nBlocks,1); % number of blocks
    block_raw{1} = {[2,3,4,5],[12,13,12,14,12,15],99,};
    block_raw{2} = {[2,100,2,3,4,100,4,5],[2,100,2,3,4,100,4,5],[12,100,12,13,12,100,12,14,12,100,12,15],99};
    block_raw{3} = {[2,3,4,5],[12,13,12,14,12,15],99};

    stimset = [];
    stimset(1).name = 'PT';
    stimset(1).ij = [1,1; 3,1;];
    stimset(2).name = 'OMR';
    stimset(2).ij = [1,2; 3,2;];
    stimset(3).name = 'Spontaneous';
    stimset(3).ij = [1,3; 2,4;3,3];
    stimset(4).name = 'PT-shock';
    stimset(4).ij = [2,1; 2,2;];
    stimset(5).name = 'OMT-shock';
    stimset(5).ij = [2,3];

    % M_range_raw = [2,3,4,5,12,13,14,15,99,100]; % for Fish 8
    M_range =       [3,1,3,2, 3,10,11,12, 0,16]; % standardized

elseif i_fish >=9 && i_fish <=11,
    nBlocks = 3;
    block_raw = cell(nBlocks,1); % number of blocks
    block_raw{1} = {[2,3,4,5],[12,13,12,14,12,15],99,[23,22],[30,31,30,32],[42,43,44,45]};
    block_raw{2} = {[2,100,2,3,4,100,4,5],99,[12,100,12,13,12,100,12,14,12,100,12,15],99,[23,22],[30,31,30,32],[42,100,42,43,44,100,44,45]};
    block_raw{3} = {[2,3,4,5],[12,13,12,14,12,15],99,[23,22],[30,31,30,32],[42,43,44,45],[2,3,4,5],[12,13,12,14,12,15]};

    stimset = [];
    stimset(1).name = 'PT';
    stimset(1).ij = [1,1; 3,1; 3,7];
    stimset(2).name = 'OMR';
    stimset(2).ij = [1,2; 3,2; 3,8];
    stimset(3).name = 'Spontaneous';
    stimset(3).ij = [1,3; 2,4; 3,3];
    stimset(4).name = 'Dot/prey';
    stimset(4).ij = [1,4; 2,5; 3,4];
    stimset(5).name = 'Looming';
    stimset(5).ij = [1,5; 2,6; 3,5];
    stimset(6).name = 'Red|Blue';
    stimset(6).ij = [1,6; 3,6];
    stimset(7).name = 'PT-shock';
    stimset(7).ij = [2,1];
    stimset(8).name = 'OMR-shock';
    stimset(8).ij = [2,3];
    stimset(9).name = 'RB-shock';
    stimset(9).ij = [2,7];

    % M_range_raw = [2,3,4,5,12,13,14,15,22,23,30,31,32,42,43,44,45,99,100] % for Fish 10
    M_range =       [3,1,3,2, 3,10,11,12,13, 0, 3,14,15,23,21,23,22, 0,16]; % standardized

elseif i_fish >=12 && i_fish <=14,
    nBlocks = 2;
    block_raw = cell(nBlocks,1); % number of blocks
    block_raw{1} = {[2,3,4,5],[12,13,12,14,12,15,12,16],[0,12],99,[30,31,30,32]};
    block_raw{2} = {[2,3,4,5],[12,13,12,14,12,15,12,16],[0,12],99,[30,31,30,32]};

    stimset = [];
    stimset(1).name = 'PT';
    stimset(1).ij = [1,1; 2,1];
    stimset(2).name = 'OMR';
    stimset(2).ij = [1,2; 2,2];
    stimset(3).name = 'DF';
    stimset(3).ij = [1,3; 2,3];
    stimset(4).name = 'Spontaneous';
    stimset(4).ij = [1,4; 2,4];
    stimset(5).name = 'Looming';
    stimset(5).ij = [1,5; 2,5];

    % M_range_raw = [0,2,3,4,5,12,13,14,15,16,30,31,32,99] % for Fish 12
    M_range =       [0,3,1,3,2, 3,10, 9,11,12, 3,14,15, 4]; % standardized

% elseif i_fish == 12 || i_fish == 13 || i_fish == 14,
%     nBlocks = 1;
%     block_raw = cell(nBlocks,1); % number of blocks
%
%     block_raw{1} = {[0],[3],[4]};
%
%     stimset = [];
%     stimset(1).name = 'pre-para';
%     stimset(1).ij = [1,1];
%     stimset(2).name = 'add para';
%     stimset(2).ij = [1,2];
%     stimset(3).name = 'post-para';
%     stimset(3).ij = [1,3];
%
%     % M_range_raw = [1,3,4] % for Fish 10
%     M_range =       [0,3,4]; % standardized
elseif ismember(i_fish,[15,17,18]),
%     M_range_raw = [2   3   4   5  12  13  14  15  16  17  18  99]
    nBlocks = 2;
    block_raw = cell(nBlocks,1); % number of blocks
    block_raw{1} = {[2,3,4,5],[12,13,12,14,12,15,12,16],[17,18],99,[30,31,30,32]};
    block_raw{2} = {[2,3,4,5],[12,13,12,14,12,15,12,16],[17,18],99,[30,31,30,32]};

    stimset = [];
    stimset(1).name = 'PT';
    stimset(1).ij = [1,1; 2,1];
    stimset(2).name = 'OMR';
    stimset(2).ij = [1,2; 2,2];
    stimset(3).name = 'DF';
    stimset(3).ij = [1,3; 2,3];
    stimset(4).name = 'Spontaneous';
    stimset(4).ij = [1,4; 2,4];
    stimset(5).name = 'Looming';
    stimset(5).ij = [1,5; 2,5];

    % M_range_raw = [2,3,4,5,12,13,14,15,16,17,18,30,31,32,99] % for Fish 12
    M_range =       [3,1,3,2, 3,10, 9,11,12, 3, 0, 3,14,15, 4]; % standardized

% elseif i_fish==16,
%     %     M_range_raw = [2   3   4   5  12  13  14  15  16  17  18  99]
%     nBlocks = 1;
%     block_raw = cell(nBlocks,1); % number of blocks
%     block_raw{1} = {31,[2,3,4,5],[12,13,12,14,12,15,12,16],31};
%
%     stimset = [];
%     stimset(1).name = 'NA';
%     stimset(1).ij = [1,1; 2,1];
%     stimset(2).name = 'PT';
%     stimset(2).ij = [1,1; 2,1];
%     stimset(3).name = 'OMR';
%     stimset(3).ij = [1,2; 2,2];
%     stimset(4).name = 'NA';
%     stimset(4).ij = [1,3; 2,3];
%
%
%     % M_range_raw = [2,3,4,5,12,13,14,15,16,31] % for Fish 12
%     M_range =       [3,1,3,2, 3,10, 9,11,12,0]; % standardized

else % inspect manually to set these manual params
    % set M_range here after first inspection of M_range_raw:
    M_range =       [3,1,3,2, 9,10,11,12, 0,13, 3,14,15, 4,16]; % standardized

    M_replace = zeros(1,max(M_range_raw)+1); % starting from zero
    for i = 0:max(M_range_raw),
        ix = find(i==M_range_raw);
        if ~isempty(ix),
            M_replace(i+1) = M_range(ix);
        end
    end

    stim_full = zeros(size(st));
    for i = 1:length(M_range_raw),
        stim_full(st==M_range_raw(i)) = M_range(i);
    end

    figure;
    plot(stim_full);

    stimset = [];
    return;
end
%% Raw stimulus key:
% PT:Red/Dark
% 2,4: Red|Red
% 3: Dark|Red
% 5: Red|Dark
%
% PT:Red/Blue
% 42,44: Red|Red
% 43: Blue|Red
% 45: Red|Blue

% OMR
% 12: Red
% 13: OMR F
% 14: OMR R
% 15: OMR L
% 16: OMR?

% Spont
% 99: same as black
%
% Dot(prey)
% 22: Dot
% 23: No dot
%
% Looming
% 30: Red
% 31: Blob R
% 32: Blob L

%% Standardized stimulus key:
% (Define normalized code, e.g. same stimulus under different raw stim keys can be unified)

% 0 = all black; 1 = phototaxis R; 2 = phototaxis L; 3 = all white; 4 = all gray;
% 5 = gray/black; 6 = white/gray; 7 = black/gray; 8 = gray/white.
% 9 = backward grating
% 10 = forward grating (very slow, more for calibration)
% 11 = rightward grating
% 12 = leftward grating
% 13 = Dot
% 14 = looming L %Blob R
% 15 = looming R %Blob L

% 16 = electric shock (spike)

% 21,22,23 = 'red|blue R','blue|red L','red|red'

% % 0 = all black; 1 = black/white; 2 = white/black; 3 = all white; 4 = all gray;
% % 5 = gray/black; 6 = white/gray; 7 = black/gray; 8 = gray/white.
% % 9 = backward grating... %%OMR baseline = not moving?? grey??
% % 10 = forward grating (very slow, more for calibration)
% % 11 = rightward grating
% % 12 = leftward grating
% % 13 = Dot
% % 14 = Blob R
% % 15 = Blob L
%
% % 16 = electric shock (spike)
%
% % 21,22,23 = 'red|blue','blue|red','red|red'

%% correct transient frames
% short method: st_rounded = interp1(M_range,M_range,st,'nearest');
% but this doesn't work if a frame is different from both the one
% before and the one after. Use loop below to force a manual 'round'
% to the closer one of the 2 neighbors.
temp = diff([0,st]);
switches = find(temp);
trans = switches(diff(switches)==1);

for i = 1:length(trans),
    j = trans(i);
    if j>1,
        st(j) = st(j-1) + round((st(j)-st(j-1))/abs(st(j+1)-st(j-1))); % round to nearest of st(j-1) and st(j+1)
    end
end

% (loop until all transient ones are corrected)
temp = diff([0,st]);
switches = find(temp);
trans = switches(diff(switches)==1);
count = 0;
while ~isempty(trans),
    for i = 1:length(trans),
        j = trans(i);
        st(j) = st(j-1) + round((st(j)-st(j-1))/abs(st(j+1)-st(j-1))); % round to nearest of st(j-1) and st(j+1)
    end
    temp = diff([0,st]);
    switches = find(temp);
    trans = switches(diff(switches)==1);
    count = count+1;
    if count>5,
        disp('count>5??');
        break;
    end
end

%% Standardize stimulus keys for both 'stim' and 'block', based on 'M_range_raw'->'M_range'
% (tricky for cell arrays) custom method:
% make a replacement matrix 'M_replace' to substitude raw stimulus key in 'block_raw'
%  to normalized key, stored in new cell array 'block'
M_replace = zeros(1,max(M_range_raw)+1); % starting from zero
for i = 0:max(M_range_raw),
    ix = find(i==M_range_raw);
    if ~isempty(ix),
        M_replace(i+1) = M_range(ix);
    end
end

stim_full = zeros(size(st));
block = cell(nBlocks,1);
for i = 1:length(M_range_raw),
    stim_full(st==M_range_raw(i)) = M_range(i);
    for j = 1:nBlocks,
        block{j} = cellfun(@(x) M_replace(x+1),block_raw{j},'UniformOutput',0);
    end
end
stim_full_1 = horzcat(stim_full,-1);% set a virtual next frame to a invalid number for segmentation later

%% Stimulus set segmentation, based on 'block' (manual input of protocol sequence)
% reduce 'stim' to a sequence of keys without consecutive repeats
ix_singles = [1,find(diff(stim_full_1))+1]; % index array to map back to raw indices (frame-number)
singles = stim_full_1(ix_singles); % sequence of keys without consecutive repeats.
% '_sg' below stands for 'singles'. Manipulations below in both raw and singles, in parallel.

jump = 1; % searching for 'jumps' in 'singles'
% Based on the known protocol sequence for each set, look for first
% mismatch in 'singles', i.e. position where the next set begins, stored
% as 'block_change(_sg)' and 'set_start(_sg)'.
block_change = cell(nBlocks,1);
block_change_sg = cell(nBlocks,1);
flag_isbreakagain = false;
for I = 1:nBlocks,
    nSets = length(block{I});
    block_change{I} = zeros(1,nSets);
    block_change_sg{I} = zeros(1,nSets);
    for J = 1:nSets,
        jump_last = jump;
        pattern = block{I}{J}; % protocol sequence for this set

        % this is sort of awkward. Instead of padding the 'already-found' sets with zeros,
        % one could also search the cropped array and adjusted the index.
        % Now 'comparison' is aligned with stim, works.
        textile = repmat(pattern,1,ceil(length(singles)/length(pattern))); % tile, for use below

        template = zeros(size(singles));
        ixs_fill = jump_last:length(singles);
        template(ixs_fill) = textile(1:length(ixs_fill)); % crop the 'textile' to fill the space after jump_last

        singles_temp = singles;
        if jump_last>1,
            singles_temp(1:jump_last-1) = 0;
        end
        comparison = (singles_temp-template);
        % find the position where the next set begins
        jump = find((comparison),1,'first');

        %% contingency plan: if start of the new set is different from
        % expected protocol ('block')
        count_stuck = 0;
        count_stuck2 = 0;
        while jump-jump_last<length(pattern),
            count_stuck = count_stuck+1;
            if count_stuck>length(pattern),
                count_stuck = 0;
                jump_last = jump_last+1;
                count_stuck2 = count_stuck2+1;
                disp(['count_stuck2=' num2str(count_stuck2)]);
                if count_stuck2>10, % arb. thres
                    disp('count_stuck2>10');
                    break;
                end
            end

            pattern_ = circshift(pattern,-count_stuck,2);
            textile = repmat(pattern_,1,ceil(length(singles)/length(pattern_)));

            template = zeros(size(singles));
            ixs_fill = jump_last:length(singles);
            template(ixs_fill) = textile(1:length(ixs_fill));

            singles_temp = singles;
            if jump_last>1,
                singles_temp(1:jump_last-1) = 0;
            end
            comparison = (singles_temp-template);

            jump = find((comparison),1,'first');

        end
        %% save
        if isempty(jump), % reached end of experiment!
            flag_isbreakagain = true;
            break;
        end
        block_change_sg{I}(J) = jump;
        block_change{I}(J) = ix_singles(jump);
    end
    if flag_isbreakagain,
        break;
    end
end

% eliminate empty sets
for I = 1:nBlocks,
    nSets = length(block_change{I});
    if nSets == 0,
        block_change{I} = [];
        block_change_sg{I} = [];
        block{I} = [];
        nBlocks = I-1;
    else
        for J = 1:nSets,
            if block_change{I}(J) == 0,
                block_change{I}(J:end) = [];
                block_change_sg{I}(J:end) = [];
                if length(block{I})>= J,
                    block{I}(J:end) = [];
                end
                break;
            end
        end
    end
end

% block_change is basically start of next set. shift by one set to get 'set_start'.
set_start = cell(nBlocks,1);
set_start_sg = cell(nBlocks,1);
set_stop = cell(nBlocks,1);
set_stop_sg = cell(nBlocks,1);
for I = 1:nBlocks,
    nSets = length(block{I});
    set_start{I} = zeros(1,nSets);
    set_start_sg{I} = zeros(1,nSets);
    if I == 1,
        set_start{I}(1) = 1;
        set_start_sg{I}(1) = 1;
    else
        set_start{I}(1) = block_change{I-1}(end);
        set_start_sg{I}(1) = block_change_sg{I-1}(end);
    end
    set_start{I}(2:end) = block_change{I}(1:end-1);
    set_start_sg{I}(2:end) = block_change_sg{I}(1:end-1);
    set_stop{I} = block_change{I}-1;
    set_stop_sg{I} = block_change_sg{I}-1;
end

%% Visualize current segmentation of blocks
if isPlotting,
    figure;
    hold on;
    plot(stim_full_1);
    for I = 1:nBlocks,
        for J = 1:length(block_change{I}),
            x = block_change{I}(J);
            plot([x,x],[-1,max(M_range)+1],'r--')
            ylim([-1,max(M_range)+1]);
        end
    end
end
%% Find shift-corrections for each stim-set, used for averaging later

%% Find periods for each set-stype, and organize info by type into 'stimset'.
for i_ss = 1:length(stimset),
    % find empty sets (for this particular experiment) and delete ij's in 'stimset'
    for i_SetRep = 1:size(stimset(i_ss).ij,1),
        I = stimset(i_ss).ij(i_SetRep,1);
        J = stimset(i_ss).ij(i_SetRep,2);
        if I>length(block) || J>length(block{I}),
            stimset(i_ss).ij(i_SetRep:end,:) = [];
            break;
        end
    end
    nSetReps = size(stimset(i_ss).ij,1);
    if nSetReps == 0,
       stimset(i_ss) = [];
       break;
    end

    %% get period from first set of block
    I = stimset(i_ss).ij(1,1);
    J = stimset(i_ss).ij(1,2);
    pattern = block{I}{J}; % protocol sequence (singles) in normalized code
    stimset(i_ss).pattern = pattern;

    % find periodicity
    if length(pattern)==1, % exception, e.g. 'Spontaneous' set, single stimlus key held for full duration of set
        if J < length(block{I}),
            next_start = set_start{I}(J+1);
        elseif J == length(block{I}),
            if I < nBlocks,
                next_start = set_start{I+1}(1);
            else
                next_start = length(stim_full_1);
            end
        end
        stimset(i_ss).period = next_start - set_start{I}(J);
    else % find periodicity based on stimlus key repetition
        ix_start = set_start_sg{I}(J);
        segment = singles(ix_start:ix_start+length(pattern)-1);

        for i_segment = 1:length(segment),
            if length(find(segment==segment(i_segment)))==1, % use a stimulus that appears only once per period
                % find period
                start_s = set_start_sg{I}(J) + i_segment -1;
                stop_s = start_s + find(singles(1,start_s+1:end)==segment(i_segment),1,'first');
                stimset(i_ss).period = ix_singles(stop_s+1)-ix_singles(start_s+1); % shift by 1!! in case first stim is not intact
                break;
            end
        end
    end

    %%
    %     stimset(i_ss).period % assigned above

    stimset(i_ss).rawstarts = zeros(1,nSetReps);
    stimset(i_ss).rawstops = zeros(1,nSetReps);
    stimset(i_ss).starts = zeros(1,nSetReps);
    stimset(i_ss).stops = zeros(1,nSetReps);
    stimset(i_ss).nReps = zeros(1,nSetReps);
    for i_SetRep = 1:nSetReps,
        I = stimset(i_ss).ij(i_SetRep,1);
        J = stimset(i_ss).ij(i_SetRep,2);
        stimset(i_ss).rawstarts(i_SetRep) = set_start{I}(J);
        stimset(i_ss).rawstops(i_SetRep) = block_change{I}(J)-1;

        if length(pattern)==1, % exception, e.g. 'Spontaneous' set, single stimlus key held for full duration of set
            stimset(i_ss).starts(i_SetRep) = stimset(i_ss).rawstarts(i_SetRep);
            % length of following occurences need to match period
            % (period primitively extracted from first occurence)
%             i_start = stimset(i_ss).rawstarts(i_SetRep);
            i_stop = stimset(i_ss).rawstops(i_SetRep);

            stimset(i_ss).nReps(i_SetRep) = 1;
            ix = stimset(i_ss).starts(i_SetRep) + stimset(i_ss).period - 1;
            if ix < i_stop,
                stimset(i_ss).stops(i_SetRep) = ix;
            else
                stimset(i_ss).stops(i_SetRep) = i_stop;
            end

        else
            % circshift to find intact periods
            pattern = stimset(i_ss).pattern;

            ix_start = set_start_sg{I}(J);
            ix_stop = set_stop_sg{I}(J);

            % find 'shift' (shift = 1 means starting from 2nd in 'singles')
            segment = singles(ix_start:ix_start+length(pattern)-1);
            shift = 0;
            while ~isequal(segment,pattern),
                shift = shift+1;
                segment = singles(ix_start+shift:ix_start+length(pattern)-1+shift);
            end
            stimset(i_ss).starts(i_SetRep) = ix_singles(ix_start+shift);
            % round to period within this set
            nReps = floor((ix_stop-ix_start-shift+1)/length(pattern));
            stimset(i_ss).nReps(i_SetRep) = nReps;

            ix = stimset(i_ss).starts(i_SetRep) + nReps*stimset(i_ss).period - 1;
            if ix>stimset(i_ss).rawstops(i_SetRep),
                nReps = nReps-1;
            end
            stimset(i_ss).stops(i_SetRep) = stimset(i_ss).starts(i_SetRep) + nReps*stimset(i_ss).period - 1;
        end
    end
end

% discard very first period because of potential artifacts at beginning of
% entire experiment
if stimset(1).nReps(1)>1,
    stimset(1).nReps(1) = stimset(1).nReps(1) - 1;
    stimset(1).starts(1) = stimset(1).starts(1) + stimset(1).period;
end

end

author: xiuyechen

repo: FishExplorer


text:
%% newton raphson example
% Find the Darcy friction factor for pipe flow using the Colebrook
% equation.
fprintf('\n**************************************************\n')
fprintf('NON-LINEAR SYSTEM OF EQUATIONS - PIPE FLOW EXAMPLE\n')
fprintf('**************************************************\n')
%% References:
% * http://en.wikipedia.org/wiki/Darcy_friction_factor_formulae
% * http://en.wikipedia.org/wiki/Darcy%E2%80%93Weisbach_equation
% * http://en.wikipedia.org/wiki/Moody_chart
% * http://www.mathworks.com/matlabcentral/fileexchange/35710-iapwsif97-functional-form-with-no-slip
%% inputs:
p = 0.68; % [MPa] water pressure (100 psi)
dp = -0.068*1e6; % [Pa] pipe pressure drop (10 psi)
T = 323; % [K] water temperature
D = 0.10; % [m] pipe hydraulic diameter
L = 100; % [m] pipe length
roughness = 0.00015; % [m] cast iron pipe roughness
rho = 1./IAPWS_IF97('v_pT',p,T); % [kg/m^3] water density (988.1 kg/m^3)
mu = IAPWS_IF97('mu_pT',p,T); % [Pa*s] water viscosity (5.4790e-04 Pa*s)
Re = @(u) rho*u*D/mu; % Reynolds number
%% governing equations
% Use Colebrook and Darcy-Weisbach equation to solve for pipe flow.

% friction factor (Colebrook eqn.)
residual_friction = @(u, f) 1/sqrt(f) + 2*log10(roughness/3.7/D + 2.51/Re(u)/sqrt(f));
% pressure drop (Darcy-Weisbach eqn.)
residual_pressdrop = @(u, f) rho*u^2*f*L/2/D + dp;
% residuals
fun = @(x) [residual_friction(x(1),x(2)), residual_pressdrop(x(1),x(2))];
%% solve
x0 = [1,0.01]; % initial guess
fprintf('\ninitial guess: u = %g[m/s], f = %g\n',x0) % display initial guess
options = optimset('TolX',1e-12); % set TolX
[x, resnorm, f, exitflag, output, jacob] = newtonraphson(fun, x0, options);
fprintf('\nexitflag: %d, %s\n',exitflag, output.message) % display output message
%% results
fprintf('\nOutputs:\n')
properties = {'Pressure','Pressure-drop','Temperature','Diameter','Length', ...
    'roughness','density','viscosity','Reynolds-number','speed','friction'};
units = {'[Pa]','[Pa]','[C]','[cm]','[m]','[mm]','[kg/m^3]','[Pa*s]','[]','[m/s]','[]'};
values = {p*1e6,dp,T-273.15,D*100,L,roughness*1000,rho,mu,Re(x(1)),x(1),x(2)};
fprintf('%15s %10s %10s\n','Property','Unit','Value')
results = [properties; units; values];
fprintf('%15s %10s %10.4g\n',results{:})
%% comparison
% solve using Haaland
Ntest = 10;
u0 = linspace(x(1)*0.1, x(1)*10, Ntest); % [m/s]
Re0 = Re(u0);
f0 = (1./(-1.8*log10((roughness/D/3.7)^1.11 + 6.9./Re0))).^2;
u0 = sqrt(-dp/rho./(f0*L/2/D));
% plot
plot(u0, f0, '-', u0, x(2)*ones(1,Ntest), '--', x(1)*ones(1,Ntest), f0, '--')
grid
title('Pipe flow solution using Haaland equation.')
xlabel('water speed, u [m/s]'),ylabel('friction factor, f')
legend('f_{Haaland}',['f = ',num2str(x(2))], ['u = ',num2str(x(1)),' [m/s]'])
%% LSQ Curve Fitting
fprintf('\n**********************************************\n')
fprintf('LEAST SQUARES CURVE FITTING WITH NEWTONRAPHSON\n')
fprintf('**********************************************\n')
% independent variables
[x,y] = meshgrid(0:10,0:10);
% bivariate distribution
bivar = @(x1,x2,sig,u1,u2)1/2/pi/sig^2*exp(-1/2*(((x1-u1).^2+(x2-u2).^2)/sig));
sigma = 3; ux = 4; uy = 5; % std dev, x & y means
z = bivar(x,y,sigma,ux,uy); % dist
% plot
figure,contour(x,y,z),hold('all')
title('lsq curve-fitting bivariate pdf with newtonraphson')
xlabel('x-coord'),ylabel('y-coord') % axes titles
grid,colorbar % show colorbar and grid
z_meas = z + (2*rand(11)-1)/1e4; % measured data
% fitting function
lsqfitfun = @(c)z_meas-bivar(x,y,c(1),c(2),c(3));
% fit coefficients to fun
c0 = [1,2,3]; % initial guess
fprintf('\ninitial guess: sigma = %g, ux = %g, uy = %g\n',c0) % display initial guess
options = optimset('TolX',1e-12); % set TolX
[c, ~, ~, exitflag, output] = newtonraphson(lsqfitfun, c0, options);
fprintf('\nexitflag: %d, %s\n',exitflag, output.message) % display output message
fprintf('\ncurve-fit coefficients: sigma=%g, ux=%g, uy=%g\n',c)
lines = plot(repmat(c(2),1,11),0:10,'--',0:10,repmat(c(3),1,11),'--', ...
    c(2),c(3),'o');
set(lines,'LineWidth',2);
legend('bivariate distribution','u_x','u_y','center')


text:
function [d,a] = training(a,d)

  a.dat=d;
  d=set_x(d,feval('get_distance',a,d));

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function desc = calcLabHist(rgb_im, seg, numRegion)
    if ~isa(rgb_im,'uint8'),
        rgb_im = im2uint8(rgb_im);
    end

    cform = makecform('srgb2lab');
    im = applycform(rgb_im,cform);

    binNum = 50;
    binVal = 0:256/(binNum):256;
    desc = zeros([numRegion binNum*3]);

    cnt = 0;
    ind={};
    for iReg=1:numRegion
        ind{iReg} = seg(:)==iReg;
    end

    for ch=1:3
        for bin=1:binNum
            cnt = cnt + 1;
            I = im(:,:,ch);
            I = ( (I>=binVal(bin)) & (I<binVal(bin+1)) );
            for iReg=1:numRegion
                desc(iReg, cnt) = sum(I(ind{iReg}));
            end
        end
    end

    tmp = sum(desc, 2);
    desc = (desc ./ repmat(tmp(:), [1 size(desc,2)]))*3;
end

author: kittenish

repo: Image-Shadow-Detection-and-Removal


text:
function x = r8bb_sl ( n1, n2, ml, mu, a_lu, pivot, b )

%*****************************************************************************80
%
%% R8BB_SL solves a R8BB system factored by R8BB_FA.
%
%  Discussion:
%
%    Note that in C++ and FORTRAN, we can look at A as an abstract
%    vector, but then look at parts of A as storing a two dimensional
%    array.  MATLAB assigns an inherent dimensionality to a data object,
%    and gets very unhappy when you try to manipulate the data yourself.
%    This means that the MATLAB implementation of this routine requires
%    the use of temporary 2D arrays.
%
%    The R8BB storage format is for a border banded matrix.  Such a
%    matrix has the logical form:
%
%      A1 | A2
%      ---+---
%      A3 | A4
%
%    with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
%    are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
%    respectively.
%
%    A should be defined as a vector.  The user must then store
%    the entries of the four blocks of the matrix into the vector A.
%    Each block is stored by columns.
%
%    A1, the banded portion of the matrix, is stored in
%    the first (2*ML+MU+1)*N1 entries of A, using standard LINPACK
%    general band format.  The reason for the factor of 2 in front of
%    ML is to allocate space that may be required if pivoting occurs.
%
%    The following formulas should be used to determine how to store
%    the entry corresponding to row I and column J in the original matrix:
%
%    Entries of A1:
%
%      1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
%
%      Store the I, J entry into location
%      (I-J+ML+MU+1)+(J-1)*(2*ML+MU+1).
%
%    Entries of A2:
%
%      1 <= I <= N1, N1+1 <= J <= N1+N2.
%
%      Store the I, J entry into location
%      (2*ML+MU+1)*N1+(J-N1-1)*N1+I.
%
%    Entries of A3:
%
%      N1+1 <= I <= N1+N2, 1 <= J <= N1.
%
%      Store the I, J entry into location
%      (2*ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
%
%    Entries of A4:
%
%      N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
%
%      Store the I, J entry into location
%      (2*ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
%      (same formula used for A3).
%
%    The linear system A * x = b is decomposable into the block system:
%
%      ( A1 A2 ) * (X1) = (B1)
%      ( A3 A4 )   (X2)   (B2)
%
%    All the arguments except B are input quantities only, which are
%    not changed by the routine.  They should have exactly the same values
%    they had on exit from R8BB_FA.
%
%    If more than one right hand side is to be solved, with the same matrix,
%    R8BB_SL should be called repeatedly.  However, R8BB_FA only needs to be
%    called once to create the factorization.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N1, N2, the order of the banded and dense blocks.
%    N1 and N2 must be nonnegative, and at least one must be positive.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%    ML and MU must be nonnegative and no greater than N1-1.
%
%    Input, real A_LU( (2*ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the factor information
%    computed by R8BB_FA.
%
%    Input, integer PIVOT(N1+N2), the pivoting information from R8BB_FA.
%
%    Input, real B(N1+N2), the right hand side of the linear system.
%
%    Output, real X(N1+N2), the solution.
%
  nband = (2*ml+mu+1)*n1;
%
%  Set B1 := inverse(A1) * B1.
%  Copy the banded matrix out of A_LU and into A1_LU.
%
  if ( 0 < n1 )

    a1_lu(1:2*ml+mu+1,1:n1) = r8vec_to_r8gb ( n1, n1, ml, mu, a_lu(1:nband) );

    job = 0;

    x(1:n1) = r8gb_sl ( n1, ml, mu, a1_lu, pivot, b, job );

  end
%
%  Modify the right hand side of the second linear subsystem.
%  Set B2 := B2 - A3*B1.
%
  for i = 1 : n2
    for j = 1 : n1
      ij = nband + n1*n2 + (j-1)*n2 + i;
      b(n1+i) = b(n1+i) - a_lu(ij) * x(j);
    end
  end
%
%  Set B2 := inverse(A4) * B2.
%  Copy the dense matrix out of A_LU and into A4_LU.
%
  if ( 0 < n2 )

    a4_lu(1:n2,1:n2) = r8vec_to_r8ge ( ...
      n2, n2, a_lu(nband+2*n1*n2+1:nband+2*n1*n2+n2*n2) );

    job = 0;
    x(n1+1:n1+n2) = r8ge_sl ( n2, a4_lu, pivot(n1+1:n1+n2), b(n1+1:n1+n2), job );

  end
%
%  Modify the first subsolution.
%  Set B1 := B1 + A2*B2.
%
  for i = 1 : n1
    for j = 1 : n2
      ij = nband + (j-1)*n1 + i;
      x(i) = x(i) + a_lu(ij) * x(n1+j);
    end
  end

  return
end