text:
function optimParam = tuneParam(LPProblem,contFunctName,timelimit,nrepeat,printLevel)
% Optimizes cplex parameters to make model resolution faster.
% Particularly interetsing for large-scale MILP models and repeated runs of
% optimisation.
% While, it won't optimize memory space nor model constraints for numerical
% infeasibilities, tuneParam will provide the optimal set of solver
% parameters for feasbile models. It requires IBM ILOG cplex (for now).
%
% USAGE:
%
%    optimalParameters = tuneParam(LPProblem,contFunctName,timelimit,nrepeat,printLevel);
%
% INPUT:
%         LPProblem:     MILP as COBRA LP problem structure
%         contFunctName: Parameters structure containing the name and value.
%                        A set of routine parameters will be added by the solver
%                        but won't be reported.
%         timelimit:     default is 10000 second
%         nrepeat:       number of row/column permutation of the original
%                        problem, reports robust results.
%                        sets the CPX_PARAM_TUNINGREPEAT parameter
%                        High values of nrepeat would require consequent
%                        memory and swap.
%         printLevel:    0/(1)/2/3
%
% OUTPUT:
%         optimParam: structure of optimal parameter values directly usable as
%                     contFunctName argument in solveCobraLP function
%
% .. Author: Marouen Ben Guebila 24/07/2017

if ~changeCobraSolver('ibm_cplex')
    error('This function requires IBM ILOG CPLEX');
end

if ~exist('printLevel','var')
    printLevel = 1;
end

if exist('timelimit','var')
    contFunctName.tune.timelimit = timelimit;
end
if exist('nrepeat','var')
    contFunctName.tune.repeat = nrepeat;
end
if exist('printLevel','var')
    contFunctName.tune.display = printLevel;
end

%read parameters
if isstruct(contFunctName)
    cpxControl=contFunctName;
else
    if ~isempty(contFunctName)
        %calls a user specified function to create a CPLEX control structure
        %specific to the users problem. A TEMPLATE for one such function is
        %CPLEXParamSet
        cpxControl=eval(contFunctName);
    else
        cpxControl=[];
    end
end

if ~isfield(LPProblem,'A')
    if ~isfield(LPProblem,'S')
            error('Equality constraint matrix must either be a field denoted A or S.')
    end
    LPProblem.A=LPProblem.S;
end

if ~isfield(LPProblem,'csense')
    nMet=size(LPProblem.A);
    if printLevel>0
        fprintf('%s\n','Assuming equality constraints, i.e. S*v=b');
    end
    %assuming equality constraints
    LPProblem.csense(1:nMet,1)='E';
end

if ~isfield(LPProblem,'osense')
    %assuming maximisation
    LPProblem.osense=-1;
    if printLevel>0
        fprintf('%s\n','Assuming maximisation of objective');
    end
end

if size(LPProblem.A,2)~=length(LPProblem.c)
    error('dimensions of A & c are inconsistent');
end

if size(LPProblem.A,2)~=length(LPProblem.lb) || size(LPProblem.A,2)~=length(LPProblem.ub)
    error('dimensions of A & bounds are inconsistent');
end

%get data
[c,x_L,x_U,b,csense,osense] = deal(LPProblem.c,LPProblem.lb,LPProblem.ub,...
    LPProblem.b,LPProblem.csense,LPProblem.osense);
%modify objective to correspond to osense
c=full(c*osense);

%cplex expects it dense
b=full(b);

%complex ibm ilog cplex interface
if ~isempty(csense)
    %set up constant vectors for CPLEX
    b_L(csense == 'E',1) = b(csense == 'E');
    b_U(csense == 'E',1) = b(csense == 'E');
    b_L(csense == 'G',1) = b(csense == 'G');
    b_U(csense == 'G',1) = Inf;
    b_L(csense == 'L',1) = -Inf;
    b_U(csense == 'L',1) = b(csense == 'L');
end

% Initialize the CPLEX object
try
    ILOGcplex = Cplex('fba');
catch ME
    error('CPLEX not installed or licence server not up')
end

ILOGcplex.Model.sense = 'minimize';

% Now populate the problem with the data
ILOGcplex.Model.obj   = c;
ILOGcplex.Model.lb    = x_L;
ILOGcplex.Model.ub    = x_U;
ILOGcplex.Model.A     = LPProblem.A;
ILOGcplex.Model.lhs   = b_L;
ILOGcplex.Model.rhs   = b_U;

%loop through parameters
ILOGcplex = setCplexParam(ILOGcplex, cpxControl, 1);
optimParam=cpxControl;

%Call parameter tuner
if ~ILOGcplex.tuneParam()
    fprintf('Optimal parameters found. \n');
    [paramList, paramPath] = getParamList(cpxControl, 1);
     for i=1:length(paramPath)
         try
             eval(['optimParam.' paramPath{i} '=ILOGcplex.Param.' paramPath{i} '.Cur;']);
         catch ME
             fprintf(['Parameter ' paramPath{i} ' was not found. \n']);
         end
     end
else
    fprintf('Optimisation failed. \n')
end

end


text:
% RESON_dyProg_mat.m
% Function to carry out dynamic programming componenent of SE_VQ
%
% Octave compatible
%
% Description
% Function to carry out dynamic programming method described in Ney (1989)
% and used previously in the ESPS GCI detection algorithm. The method
% considers target costs and transition costs which are accumulated in
% order to select the `cheapest' path, considering previous context
%
% Inputs
%  rel_amp         : [samples] [MxNcand] Relative residual amplitudes
%  GCI_N           : [samples] [MxNcand] matrix containing M by Ncandidate GCI locations
%  F0mean          : [Hz]      [1x1] Mean fundamental frequency
%  x               : [samples] [Nx1] Speech signal
%  fs              : [Hz]      [1x1] sampling frequency
%  trans_wgt       : [integer] [1x1] Weight of transition cost
%  relAmp_wgt      : [integer] [1x1] Weight of target cost
%
% Outputs
%  GCI_opt         : [samples] [Mx1] Optimal GCI path
%
% Example
%  GCI_opt =
%  RESON_dyProg_mat(GCI_relAmp,GCI_N,F0_mean,x,fs,trans_wgt,relAmp_wgt)
%
% References
%  [1] Kane, J., Gobl, C., (2013) `Evaluation of glottal closure instant
%       detection in a range of voice qualities', Speech Communication
%       55(2), pp. 295-314.
%
% Copyright (c) 2013 Trinity College Dublin
%
% License
%  This code is a part of the Voice Analysis Toolkit with the following
%  licence:
%  The software product (and any modifications thereof) is distributed under
%  a dual licence: an open source license for individual, non-commercial
%  purposes, and a commercial license. The opensource licence under which
%  the product is distributed is GNU GPL v2. For individual users, this
%  licence suits their use as these are not distributing proprietary
%  modifications, additions to, or derivatives of the product and don't
%  require legal protection of a commercial licence. For commercial users,
%  where open source does not meet their requirements, we offer commercial
%  licensing of the product. A commercial license permits customers to
%  modify, add or produce derivative products without the obligation of
%  making the subsequent code open source. For more information regarding
%  our commercial licence, please contact john.whelan@tcd.ie
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  John Kane kanejo@tcd.ie
%
% $Id <info set by the versioning system> $

function GCI_opt = RESON_dyProg_mat(GCI_relAmp,GCI_N,F0_mean,x,fs,trans_wgt,relAmp_wgt)

%% Initial settings
plots=0;
cost = (GCI_relAmp.*relAmp_wgt)';
ncands=size(GCI_N,1);
nframe=size(GCI_N,2);
GCI_N=GCI_N';
prev=zeros(nframe,ncands);      % traceback pointer
pulseLen = round(fs/F0_mean);
GCI_opt=zeros(1,nframe);

%% Do processing
for n=1:nframe

    if n>1
        costm=zeros(ncands,ncands);         % transition cost matrix: rows (previous), cols (current)

        for c=1:ncands
            % Transitions TO states in current frame
            start=GCI_N(n,c)-round(pulseLen/2);
            stop=GCI_N(n,c)+round(pulseLen/2);
            if stop > length(x)
                stop=length(x);
            end

            if start<1
                start=1;
            end
            pulse_cur = x(start:stop);

            for p=1:ncands
                % Transitions FROM states in previous frame
                start=GCI_N(n-1,p)-round(pulseLen/2);
                stop=GCI_N(n-1,p)+round(pulseLen/2);
                if start<1
                    start=1;
                end
                if stop > length(x)
                    stop=length(x);
                end
                if start<1
                    start=1;
                end

                pulse_prev = x(start:stop);

                if isempty(pulse_cur) ||isnan(pulse_cur(1)) || ...
                        isnan(pulse_prev(1))
                    costm(p,c)=0;
                else
                    if length(pulse_cur)~=length(pulse_prev)
                        cor_cur=0;
                    else
                        cor_cur=corrcoef(pulse_cur(:),pulse_prev(:));
			if length(cor_cur)==1
			   cor_cur=0;
			else cor_cur=cor_cur(2);
			end
                    end
                    costm(p,c) = (1-abs(cor_cur))*trans_wgt; % transition cost
                end
            end
        end

        costm=costm+repmat(cost(n-1,1:ncands)',1,ncands);  % add in cumulative costs
        [costi,previ]=min(costm,[],1);
        cost(n,1:ncands)=cost(n,1:ncands)+costi;
        prev(n,1:ncands)=previ;
    end

end

%% Do traceback
best=zeros(n,1);
[cbest,best(n)]=min(cost(n,1:ncands));
for i=n:-1:2
     best(i-1)=prev(i,best(i));
end

for n=1:nframe
    GCI_opt(n) = GCI_N(n,best(n));
end

%% Do plots
if plots
    GCI_norm=zeros(nframe,ncands);
    GCI_opt_norm=zeros(nframe,ncands);
    for n=1:nframe
        GCI_norm(n,:) = GCI_N(n,:)-GCI_N(n,1);
        GCI_opt_norm(n) = GCI_opt(n)-GCI_N(n,1);
    end
    subplot(211), plot(x), hold on, stem(GCI_N(:,1),ones(1,length(GCI_N(:,1)))*-.1,'r')
        stem(GCI_opt,ones(1,length(GCI_opt))*-.1,'k')
    subplot(212), plot(GCI_opt,GCI_norm,'rx'), ylim([-20 20]),
        hold on, plot(GCI_opt,GCI_opt_norm)
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
clear all; close all

%% case 1
fprintf('\n\n\nCase 1\n')

surf(peaks)
axis vis3d
h = camlight('left');
a = get(gca);
set(gca,'view',[30,60])
b = get(gca);
b.junk = 'test';
[match, er1, er2] = comp_struct(a,b,1);
match
er1
list_struct(er1)
er2
list_struct(er2)

fprintf('\n\n\n\n\n');

if menu('Run next test?','Yes','No') - 1; error('Stop run'); end

%% case 2

clear all; close all
fprintf('\n\n\nCase 2\n')

% first struct
a(1).a = 1;
a(2).a = 2;
a(1).b.a = 'test';
a(2).b.a = 'test';
a(1).b.b = cellstr('a');
a(2).b.b = cellstr('a');
a(1).b.c = 1==1;
a(2).b.c = 1==1;
a(1).c = single([1 2 3]);
a(2).c = single([1 2 3]);
a(1).d(1).a = int8(2);
a(1).d(2).a = uint8(2);
a(1).d(3).a = int16(2);
a(1).d(4).a = uint16(2);
a(1).d(5).a = int32(2);
a(1).d(6).a = uint32(2);
a(1).d(7).a = int64(2);
a(1).d(8).a = uint64(2);
a(2).d(1).a = int8(2);
a(2).d(2).a = uint8(2);
a(2).d(3).a = int16(2);
a(2).d(4).a = uint16(2);
a(2).d(5).a = int32(2);
a(2).d(6).a = uint32(2);
a(2).d(7).a = int64(2);
a(2).d(8).a = uint64(2);
a(1).e = @(x)sin(x);
a(2).e = @(x)sin(x);
a(1).f = 1;
a(1).j = [1.1 2.2 3.3];
a(1).k = @(x)sin(x);

% second struct
b = a;
b(2).a = 1;
b(2).b.a = 'x';
b(2).b.b = cellstr('b');
b(2).b.c = 1==2;
b(2).c = single([1 2 5]);
b(2).d(1).a = int8(4);
b(2).d(2).a = uint8(54);
b(2).d(3).a = int16(5);
b(2).d(4).a = uint16(65);
b(2).d(5).a = int32(7);
b(2).d(6).a = uint32(7);
b(2).d(7).a = int64(7);
b(2).d(8).a = uint64(7);
b(2).e = @(x)tan(x);
b(1).f = 1+1e-8;
b(1).j = single(a(1).j);

% non-match
a(2).h = 1;
b(1).g = 1;
b(2).i.a = 'r';
b(3) = b(2);
b(4) = b(2);
x = .5; a(2).k = @(x)sin(x);
x = .25; b(2).k = @(x)sin(x);
a(1).l.a = 1;
b(1).l.c = 2;

a(1).o = 't';
b(1).o = 5;

% wrong order fields
a(1).m.a = 1; a(1).m.b = 1;
a(1).n.a = 1; a(1).n.b = 1;
b(1).n.b = 1; b(1).n.a = 1;
b(1).m.b = 1; b(1).m.a = 1;

[match, er1, er2] = comp_struct(a,b,1,0,1e-3);

match
er1
er2

if menu('Run next test?','Yes','No') - 1; error('Stop run'); end

%% case 3
clear all; close all
fprintf('\n\n\nCase 3\n')

a.a = 1; b.a = 1;
a.b = 2; b.b = 1;
a.c.a = 1; a.c.b = 2; b.c = 1;
a.d = cellstr('test')'; b.d = cellstr('test')';
a.e = cellstr('test')'; b.e = cellstr('a')';
a.f = 'tt'; b.f = 'tt';
a.g = 'tt'; b.g = 't';
a.h.a = 1; b.h.b = 1;
a.i.a = 1;
a.j.a = @sin; b.j.a = @sin;
a.j.b = @sin; b.j.b = @cos;
x = .1; a.j.c = @(x) sin(x); x = .5; b.j.c = @(x) sin(x);
a(2) = a(1); b(2) = b(1); a(3) = a(1);

[match, er1, er2] = comp_struct(a,b);
match
list_struct(match)
er1
list_struct(er1)
er2
list_struct(er2)

if menu('Run next test?','Yes','No') - 1; error('Stop run'); end

%% case 4
clear all; close all
fprintf('\n\n\nCase 4\n')

str = char(97:122);

% build structure
for ii = 1:26
	stu1.(str(ii)) = rand(ceil(5*rand(1)),ceil(5*rand(1)));
end
stu2 = stu1;
for ii = 1:26
	if round(rand(1))
		if rand(1) > 0.3
			if round(rand)
				stu2.(str(ii)) = rand(size(stu2.(str(ii))));
			else
				stu2.(str(ii)) = stu2.(str(ii)) + 1e-8;
			end
		else
			stu2.(str(ii)) = rand(ceil(5*rand(1)),ceil(5*rand(1)));
		end
	end
end

[match, er1, er2] = comp_struct(stu1,stu2);
match
er1
list_struct(er1)
er2
list_struct(er2)

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function lm = gallager_ach(n, epsil, P)
% returns log M achievable by gallager's random coding.
%
% Note: In the case of uncostrained input gallager's bound can be written as
%
%  P_e <= M^s   V(s)^n
%
% Where V(s) = (   Sum_y [ Sum_x Q(x) W(y|x)^{1/1+s} ]^(1+s)   )
%
% Thus we could rewrite the bound as a lower bound on $M$ directly.
%
% For AWGN we can do the same in the case when we do not restrict codewords to a narrow shell
% around the power sphere. I did that and here is the result:
%
%   gallager_ach(3000, 1e-6, 1) = 1225
%   -betaq_low_v2(1-1e-6, 3000, 1) = 1287
%
%   new_gallager() = 1206 --- much worse than good gallager!
%
% Here's how to compute new_gallager:
% mu = chi2cdf(n,n);
% f = @(x) -( (1+x)./x * log2(mu) + n/2 .* log2( 1 + A^2 ./ (1+x)) + log2(epsil/2) ./ x )
% [x_best f_best] = fminbnd(f, 0, 1); new_gal = -f - 1;
%
% Note epsil/2 and also that we should subtract one bit.

R_up = cap_awgn(P);
R_down = 0;

% Check if for this n gallager works at all
Pe = gallager_pe(P, R_down, n);
if Pe > epsil;
    lm = 0;
    return;
end

precision = 1e-3 * R_up / n;

while 1;
	if ((R_up - R_down) < precision)
		break;
	end
	R = (R_up + R_down) / 2;

	Pe = gallager_pe(P, R, n);
	if( Pe < epsil)
		R_down = R;
	else
		R_up = R;
	end
end

lm = n * R_down;

function [pe Er] = gallager_pe(P, R, n, deltap_force)
%
% This function returns upper bound (achievability) on P_e using
% Gallager's random coding error exponents
%
% TODO: use expurgated exponent for low rates!
%

% conversion A->P. Old versions are all in terms of ``amplitude'' A.
A = sqrt(P);

if (nargin < 4)
	deltap_force = [];
end

if (R == 0)
	pe = 0;
	return;
end

Rcr = 1/2 * log2( 1/2 + A^2/4 + 1/2 * sqrt(1+A^4/4) );

% convert rate to nats
R = R / log2(exp(1));
Rcrn = Rcr / log2(exp(1));

if (R >= Rcrn)
	beta = exp(2*R);
	ro = A^2/(2*beta) * (1 + sqrt(1+ 4*beta/A^2/(beta-1) ) ) - 1;

	Er = A^2/(4*beta) * ...
		( (beta+1) - (beta-1)*sqrt(1 + 4*beta/A^2/(beta-1)) ) + ...
		1/2 * log(  beta - A^2*(beta-1)/2 * ( sqrt(1 + 4*beta/A^2/(beta-1)) - 1 ) );
else
	beta = 1/2 * (1 + A^2/2 + sqrt(1+A^4/4) );
	ro = 1;

	Er = 1 - beta + A^2/2 + 1/2 * log(beta - A^2/2) + 1/2*log(beta) - R;
end;

s = ro*A^2/(2 * (1+ro)^2 * beta);

%%
%% Now instead of maximizing over delta's we choose delta = 1/s  (see Gallager 7.4.39)
%% This is suboptimal, but let it be...
%%

if (max(size(deltap_force)) > 0)
	deltap = deltap_force * n;
else
	if (s == 0)
		deltap = n;
	else
		deltap = min(1/s, n);
	end;
end

%disp(sprintf('ro = %g, deltap = %g, s = %g', ro, deltap, s));

ccdfs = chi2cdf([n-deltap n], n);

mu = ccdfs(2) - ccdfs(1);
if (mu < 1e-10 * ccdfs(1))
	mu_new = quad(@(x)( chi2pdf(x, n) ), n - deltap, n);
	disp(sprintf( [	'--- gallager_pe: computing mu using quad()\n' ...
			'    `--> mu_new = %.3g, mu_old = %.3g (diff = %.2f%%)'], ...
			mu_new, mu, 200*abs(mu_new - mu)/(mu_new + mu)));
	mu = mu_new;
end;

multip = 2*exp(s*deltap)/mu;
multip_approx = ro*A^2*exp(1)*sqrt(4*pi*n)/(1+ro)^2 /beta;

%	disp(sprintf( [	'--- gallager_pe: multiplier comparsion\n' ...
%			'    `--> multip = %.3g, approx = %.3g (diff = %.2f%%)'], ...
%			multip, multip_approx, 100*abs(multip - multip_approx)/(multip)));

pe = multip*exp(-n*Er);

if(pe > 1) pe = 1; end;

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
%NCC Normalized cross correlation
%
% M = NCC(I1, I2) is the normalized cross-correlation between the
% two equally sized image patches I1 and I2.  The result M is a scalar in
% the interval -1 (non match) to 1 (perfect match) that indicates similarity.
%
% Notes::
% - A value of 1 indicates identical pixel patterns.
% - The NCC similarity measure is invariant to scale changes in image
%   intensity.
%
% See also ZNCC, SAD, SSD, ISIMILARITY.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function m = ncc(w1, w2)

	denom = sqrt( sum(sum(w1.^2))*sum(sum(w2.^2)) );

	if denom < 1e-10,
		m = 0;
	else
		m = sum(sum((w1.*w2))) / denom;
	end

author: petercorke

repo: machinevision-toolbox-matlab


text:
% ASTROTIK by Francesco Santilli
% R2BP (Restricted Two Bodies Problem)
% H2f converts hyperbolic anomaly to true anomaly for an hyperbolic orbit.
%
% Usage: f = H2f(H,e)
%
% where: H(k) = hyperbolic anomaly [rad]
%        e = eccentricity [-] (e>1)
%        f(k) = true anomaly [rad]

function f = H2f(H,e)

    if ~(nargin == 2)
        error('Wrong number of input arguments.')
    end

    check(H,1)
    check(e,0)

    if e <= 1
        error('e must be strictly major than 1.')
    end

    ee = sqrt((e+1)/(e-1));
    f = 2*atan( ee*tanh(H/2) );

    kk = H/(2*pi);
    k = round(kk) - ((kk-fix(kk)) == 0.5);
    f = f + k*(2*pi);

end

author: Sable

repo: mcbench-benchmarks


text:
function [MStitch, result] = RGBListMain_Process(file)

for i = 1 : length(file)-1
    if i == 1
        file1 = file{i};
        im1 = imread(file1);
        im1 = rgb2gray(im1);
        MStitch.im1 = double(im1);

        [Pheight, Pwidth] = size(im1);

        MStitch.Pwidth = Pwidth;
        MStitch.Pheight = Pheight;

        MStitch.W_min = round(0.60*Pwidth);
        MStitch.W_max = round(0.83*Pwidth);
        MStitch.H_min = round(0.98*Pheight);
        MStitch.minval = 255;

        im1 = imread(file1);
        MStitch.imrgb1 = double(im1);
        im1 = rgb2gray(im1);
        MStitch.im1 = double(im1);
    else
        MStitch.im1 = double(result1);
    end
    file2 = file{i+1};
    im2 = imread(file2);
    MStitch.imrgb2 = double(im2);
    im2 = rgb2gray(im2);
    im2 = double(im2);
    [W_box, H_box, bdown, MStitch] = Fun_Match(im2, MStitch);
    [MStitch, result1] = Fun_StitchRGB(im2, W_box, H_box, bdown, MStitch);
end
result = result1;

author: decouples


text:
function [] = outputReport(response, theTitle)

    fprintf('\nHED services report for %s\n', theTitle);
    fprintf('Error report:  [%s] %s\n', response.error_type, response.error_msg);

    %% Print out the results if available
    if ~isfield(response, 'results') || isempty(response.results)
        return
    end

    results = response.results;
    fprintf('[%s] status %s: %s\n', response.service, results.msg_category, results.msg);
    if isfield(results, 'schema_version')
       fprintf('HED version: %s\n', results.schema_version);
    end
    fprintf('\nReturn data for service %s [command: %s]:\n', ...
        response.service, results.command);
    data = results.data;
    if ~iscell(data)
        fprintf('%s\n', data);
    else
        for k = 1:length(data)
            if ~isempty(data{k})
                fprintf('[%d]: %s\n', k, data{k});
            end
        end
    end

    %% Output the spreadsheet if available
    if isfield(results, 'spreadsheet')
        fprintf('\n----Spreadsheet result----\n');
        fprintf(results.spreadsheet);
    end

    if isfield(results, 'definitions') && isstruct(results.definitions)
        fprintf('\n\n----------definitions---------\n');
        defNames = fieldnames(results.definitions);
        for k=1:length(defNames)
            name = defNames{k};
            fprintf('\n%s: %s\n', defNames{k}, results.definitions.(name));
        end
    end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function this = load(this, fileName, varargin)
% loads object from struct in .mat-file, different variable names possible
%
%
%   Y = MrCopyData()
%   Y.load(filName)
%   Y.load(filName, variableName)
%
%   Options for variable name
%   1) variableName is given in varargin - load if available
%   2a) no variableName is given, and only one variable stored in fileName,
%       load this one regardless of name
%   2b) no variableName is given, and several variables stored in fileName,
%       check whether 'objectAsStruct' XOR shortened  class name (without
%       'Mr', e.g. dimInfo) is available, load this
%   2c) no variableName is given, and several variables stored in fileName,
%       'objectAsStruct' and shortened class name available, nothing is
%       loaded
%   2d) no variableName is given, and several variables stored in fileName,
%       neither 'objectAsStruct' nor shortened class name available,
%       nothing is loaded
%
% This is a method of class MrCopyData.
%
% IN
%
% OUT
%
% EXAMPLE
%   load
%
%   See also MrCopyData MrCopyData.save MrCopyData.update_properties_from

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2018-08-14
% Copyright (C) 2018 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/>.

% get content of file
matObj = matfile(fileName);
details = whos(matObj);
nVars = numel(details);

% retrieve variable name if given
if nargin > 2
    objectName = varargin{1};
    hasVariableName = 1;
else
    hasVariableName = 0;
end

% search for possible matches if no variable name is given
if ~hasVariableName

    % if only one variable stored, use that one
    if nVars == 1
        objectName = details(1).name;
    else
        % check whether any of the following is available
        hasObjectAsStruct = strcmpi('objectAsStruct', {details(:).name});
        classNameThis = class(this);
        hasClassName = strcmpi(classNameThis(3:end), {details(:).name});
        if any(hasObjectAsStruct) && ~any(hasClassName)
            objectName = 'objectAsStruct';
            % has object has struct but not class Name
        elseif ~any(hasObjectAsStruct) && any(hasClassName)
            objectName = details(hasClassName).name;
            % has class name but not object as strcut
        else % has either both or none
            warning('uniqc:MrCopyData:load', ...
                ['Nothing has been loaded. Matching variables could not be found. \n', ...
                'Either no matching variable name or too many. ', ...
                'Valid variable names are ''objectAsStruct'' and the name of the class. \n', ...
                'You can also use an explicit variable name: dimInfo.load(fileName, varName)']);
            return;
        end
    end
end
% read data
objectAsStruct = matObj.(objectName);
% update properties
this.update_properties_from(objectAsStruct);

author: translationalneuromodeling

repo: tapas


text:
% icaproj() - project ICA component activations through the
%             associated weight matrices to reconstitute the
%             observed data using only the selected ICA components.
% Usage:
%   >> [icaprojdata] = icaproj(data,weights,compindex,[datameans],chansout);
%
% Inputs:
%   data        - data matrix (chans, frames*epochs)
%   weights     - unmixing weight matrix (e.g., weights*sphere from runica())
%   compindex   - vector of ICA component indices to project
%
% Optional inputs:
%   datamean    - Optional ICA row means (for each epoch) from runica()
%                 {default 0 -> distribute data offsets among the ICA components}
%   chansout    - Optional vector of channel indices to output {default: all}
%
%
% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 11-30-96
%
% See also: icavar(), runica()

% Copyright (C) 11-30-96 Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% 11-30-96 Scott Makeig CNL / Salk Institute, La Jolla as icaproject.m
% 12-24-96 added define for verbose -sm (V1.3)
% 2-11-97  use outer product math when only one component in compindex -sm
% 3-11-97  remove row means instead of grand mean -sm
% 3-19-97  used datamean argument instead of frames/baseframes -sm
% 4-03-97  changed name to icaproj() -sm
% 6-07-97  changed order of args to conform to runica -sm
% 6-10-97  fixed data mean handling -sm
% 6-23-97  trying pseudo-inverse for non-square weights -sm
% 7-23-97  distributed baseline offset (if any) among the activations -sm
% 10-31-97 removed errcode -sm
% 12-19-00 removed sphere, shifted order of args -sm
% 05-29-01 added chansout, made more efficient -sm
% 01-25-02 reformated help & license, added links -ad

function  [icaprojdata] = icaproj(data,weights,compindex,datamean,chansout)

verbose = 0;   % default no-verbose

if nargin<3   % need 3 args
    help icaproj
    return
end
if nargin<4
  datamean = 0;  % default
end
if isempty(datamean)
  datamean = 0;  % default
end

[chans,framestot] = size(data);
if nargin<5
  chansout = [];
end
if isempty(chansout) | chansout(1) == 0
  chansout = 1:chans;
end
if min(chansout)<1 | max(chansout)> chans
  fprintf('icaproj(): chansout variable out of 1:chans range.\n')
  return
end

[mchans,epochs] = size(datamean);
frames = floor(framestot/epochs);
if epochs*frames ~= framestot | frames < 1,
    fprintf(...
        'icaproj(): frames (%d) does not divide data length (%d)\n.',...
                frames,framestot);
    return
end

[ncomps,cols] = size(compindex);
if cols>1,
  if ncomps==1,    % if row vector,
      compindex = compindex';    % make col vector
      ncomps = cols;
  else
      fprintf('icaproj(): compindex must be a vector\n');
      return
  end
end
if ncomps>chans,
  fprintf('icaproj(): compindex must have <= %d entries\n',chans);
  return
end
for a=1:ncomps-1
  for b=a+1:ncomps
      if compindex(a,1)==compindex(b,1),
            fprintf('icaproj(): component index repeated in compindex\n.');
          return
      end
  end
end
for a=1:ncomps
    if compindex(a)>chans | compindex(a)<1
      fprintf('icaproj(): component index %d out of range!\n',compindex(a));
      return
      break
    end
end

if nargin<4
  datamean = 0; % default
end

if datamean ~= 0,
  %
  % Remove row means, e.g. those subtracted prior to ICA training by runica()
  %
  if verbose==1,
     fprintf('Removing data means of each channel and data epoch...\n');
  end
  for e=1:epochs
      data(:,(e-1)*frames+1:e*frames) = ...
         data(:,(e-1)*frames+1:e*frames) - datamean(:,e)*ones(1,frames);
  end;
end
if verbose == 1
  fprintf('Final input data range: %g to %g\n', ...
                min(min(data)),max(max(data)));
end

if size(weights,1) == size(weights,2)
  iweights    = inv(weights);             % inverse weight matrix
else
  iweights    = pinv(weights);            % pseudo-inverse weight matrix
end

activations = weights(compindex,:)*data;  % activation waveforms
                                          % at desired components
if ncomps==1,  % compute outer product only for single component projection
  if verbose==1,
    fprintf('icaproj(): Projecting data for ICA component %d\n',...
                     compindex(1));
  end
else % if ncomps > 1
  if verbose==1,
    fprintf('icaproj(): Projecting data for ICA components ');
    if ncomps<32
     for n=1:ncomps % for each included component
      fprintf('%d ',compindex(n));
     end
    else
      fprintf('specified.');
    end
    fprintf('\n');         % copy selected activations
  end
end

icaprojdata = iweights(chansout,compindex)*activations;
% reconstitute selected scalp data channels from selected ICA components

author: PatternRecognition

repo: OpenBMI


text:
function fem2d_pack_test17 ( )

%*****************************************************************************80
%
%% TEST17 demonstrates S_L2NORM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  element_num = 40;
  psi_num = element_num + 1;
  quad_num = 5;

  x_max = 1.0;
  x_min = 0.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST17\n' );
  fprintf ( 1, '  S_L2NORM computes the L2 norm of a scalar function\n' );
  fprintf ( 1, '  S(X) over a region (of any dimension), assuming:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  * that there is a set of finite element basis\n' );
  fprintf ( 1, '    functions PSI;\n' );
  fprintf ( 1, '  * that the region is decomposed into a number of\n' );
  fprintf ( 1, '    elements whose areas are known;\n' );
  fprintf ( 1, '  * that the integral is to be computed by a quadrature\n' );
  fprintf ( 1, '    rule applied in the same way to each element;\n' );
  fprintf ( 1, '  * that the value of the basis functions is given\n' );
  fprintf ( 1, '    at every quadrature node in every element;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Our example will have one spatial dimension.\n' );
  fprintf ( 1, '  XMIN = %f\n', x_min );
  fprintf ( 1, '  XMAX = %f\n', x_max );
  fprintf ( 1, '  The number of intervals will be    %d\n', element_num );
  fprintf ( 1, '  The number of basis functions is   %d\n', psi_num );
  fprintf ( 1, '  The number of quadrature points is %d\n', quad_num );
%
%  Set the nodes.
%
  for i = 1 : psi_num
    node_x(i) = ( ( psi_num - i     ) * x_min   ...
                + (           i - 1 ) * x_max ) ...
                / ( psi_num     - 1 );
  end
%
%  Set the element "areas".
%
  for i = 1 : element_num
    element_area(i) = node_x(i+1) - node_x(i);
  end
%
%  Set the quadrature weights and abscissas.
%
  quad_weight(1) = 0.236926885056189087514264040720;
  quad_weight(2) = 0.478628670499366468041291514836;
  quad_weight(3) = 0.568888888888888888888888888889;
  quad_weight(4) = 0.478628670499366468041291514836;
  quad_weight(5) = 0.236926885056189087514264040720;

  quad_weight(1:5) = quad_weight(1:5) / 2.0;

  quad_x(1) = - 0.906179845938663992797626878299;
  quad_x(2) = - 0.538469310105683091036314420700;
  quad_x(3) =   0.0;
  quad_x(4) =   0.538469310105683091036314420700;
  quad_x(5) =   0.906179845938663992797626878299;
%
%  Set the finite element coefficients of S.  In our formulation,
%  these finite element coefficients are simply the function values
%  at the nodes.
%
  for i = 1 : psi_num
    s_coef(i) = sin ( node_x(i) );
  end
%
%  For each basis function I,
%    for each element J,
%      for each quadrature point K,
%
%  ...evaluate the basis function in the element at the quadrature point.
%
  psi_quad(1:psi_num,1:element_num,1:quad_num) = 0.0;

  for psi = 1 : psi_num

    for element = 1 : element_num

      if ( element < psi - 1 )

        continue

      elseif ( element == psi - 1 )

        xl = node_x(element);
        xr = node_x(element+1);

        for quad = 1 : quad_num

          x = ( ( 1.0 - quad_x(quad) ) * xl   ...
              + ( 1.0 + quad_x(quad) ) * xr ) ...
              /   2.0;

          psi_quad(psi,element,quad) = ( x - xl ) / ( xr - xl );

        end

      elseif ( element == psi )

        xl = node_x(element);
        xr = node_x(element+1);

        for quad = 1 : quad_num

          x = ( ( 1.0 - quad_x(quad) ) * xl   ...
              + ( 1.0 + quad_x(quad) ) * xr ) ...
              /   2.0;

          psi_quad(psi,element,quad) = ( xr - x ) / ( xr - xl );

        end

      elseif ( psi < element )
        continue
      end

    end

  end
%
%  Ask S_L2NORM to compute the integral of L2 norm of S.
%
  l2norm = s_l2norm ( psi_num, element_num, quad_num, element_area, ...
    quad_weight, psi_quad, s_coef );

  fprintf ( 1, '\n');
  fprintf ( 1, '  The computed L2 norm is %f\n', l2norm );
%
%  The integral of (sin(x))**2 is x/2 - sin(2x)/4
%
  exact = sqrt ( ...
      ( 0.5 * x_max - 0.25 * sin ( 2.0 * x_max ) ) ...
    - ( 0.5 * x_min - 0.25 * sin ( 2.0 * x_min ) ) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The exact value is      %f\n', exact );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function c8mat_indicator_test ( )

%*****************************************************************************80
%
%% C8MAT_INDICATOR_TEST tests C8MAT_INDICATOR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2015
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8MAT_INDICATOR_TEST\n' );
  fprintf ( 1, '  C8MAT_INDICATOR returns the complex indicator matrix.\n' );

  a = c8mat_indicator ( m, n );

  c8mat_print ( m, n, a, '  The indicator matrix:' );

  return
end


text:
function slshowtext(T)
%SLSHOWTEXT Displays the text
%
% $ Syntax $
%   - slshowtext(T)
%
% $ Description $
%   - slshowtext(T) shows the text in cell array of lines.
%
% $ History $
%   - Created by Dahua Lin, on Aug 9th, 2006
%

if ~isempty(T)
    nlines = length(T);
    for i = 1 : nlines
        curline = T{i};
        if isempty(curline);
            disp(' ');
        else
            disp(curline);
        end
    end
end

author: lmthang


text:
function stuff = get(scnl, prop_name)
   %GET for the scnl object
   %  result = get(scnlobject, property), where PROPERTY is one of the
   %  following:
   %    'network', 'station', 'location', 'channel'
   %    'channeltag'
   %
   % If the results of a single SCNL are requested, then a string is returned.
   % Otherwise, a cell of values will be returned.

   lower_prop_name = lower(prop_name);
   if ~strcmp(lower_prop_name,prop_name)
      warning('SCNLOBJECT:get:propertyWarning',...
              'Use lowercase ''%s'' property name for consistency',lower(prop_name));
      prop_name = lower_prop_name;
   end

   switch prop_name

      case{'station','channel','network','location'}
         thesetags = [scnl.tag];
         stuff = {thesetags.(prop_name)};
      case{'nscl_string'}
         for n = numel(scnl): -1 : 1
            thistag = scnl(n).tag;
            stuff(n) = {[ thistag.network, '_', thistag.station, '_',...
               thistag.channel, '_', thistag.location ]};
         end
      case {'channeltag'}
         stuff = [scnl.tag];
         return
      otherwise
         error('SCNLOBJECT:get:UnrecognizedProperty',...
            'Unrecognized property name : %s',  prop_name);
   end

   %if a single scnl, then return the string representation, else return a
   %cell of strings.
   if numel(stuff) == 1
      stuff = stuff{1};
   else
      stuff = reshape(stuff,size(scnl));
   end;
end

author: geoscience-community-codes

repo: GISMO


text:
%
%   Calculate the scalar multiplication of a scalar t and vector u
%
%  Syntax:  >> v = VectorScale(s,u)
%
%   Input:  t --- (numeric) the scalar
%           u --- (matrix/string/cell) a general vector that is either
%                (i) an mxn matrix, or
%                (ii) a polynomial in character string, or
%                (iii) a 1xk cell array, each u{j} is either (i) or (ii)
%
%  Example:
% >> v = VectorScale(3,{[1 2 3; 4 5 6], '7+8*x*y+9*z^3'});
% >> v{:}
% ans =
%      3     6     9
%     12    15    18
% ans =
% 21 + 24*x*y + 27*z^3
%

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function MRSIStruct = setSpectralWidth(MRSIStruct, value)
    MRSIStruct.spectralWidth = value;
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function cn_leg_test ( n, expon )

%*****************************************************************************80
%
%% CN_LEG_TEST tests the rules for CN with Legendre weight on a monomial.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 March 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  N = %d\n', n );
  fprintf ( 1, '  EXPON = ' );
  for d = 1 : n
    fprintf ( 1, '  %2d', expon(d) );
  end
  fprintf ( 1, '\n' );
  d = sum ( expon(1:n) );
  fprintf ( 1, '  Degree = %d\n', d );
  fprintf ( 1, '\n' );

  exact = cn_leg_monomial_integral ( n, expon );

  p = 1;

  if ( d <= p )
    [ o, x, w ] = cn_leg_01_1 ( n );
    v = monomial_value ( n, o, x, expon );
    quad = w' * v';
    err = abs ( quad - exact );
    fprintf ( 1, '  CN_LEG_01_1:    %6d  %14.6g  %14.6e\n', o, quad, err );
  end

  p = 2;

  if ( d <= p )

    [ o, x, w ] = cn_leg_02_xiu ( n );
    v = monomial_value ( n, o, x, expon );
    quad = w' * v';
    err = abs ( quad - exact );
    fprintf ( 1, '  CN_LEG_02_XIU:  %6d  %14.6g  %14.6e\n', o, quad, err );

    gamma0 = 1.0;
    delta0 = 0.0;
    c1 = 1.0 / 3.0;
    volume_1d = 2.0;
    [ o, x, w ] = gw_02_xiu ( n, gamma0, delta0, c1, volume_1d );
    v = monomial_value ( n, o, x, expon );
    quad = w' * v';
    err = abs ( quad - exact );
    fprintf ( 1, '  GW_02_XIU:      %6d  %14.6g  %14.6e\n', o, quad, err );

  end

  p = 3;

  if ( d <= p )

    [ o, x, w ] = cn_leg_03_1 ( n );
    v = monomial_value ( n, o, x, expon );
    quad = w' * ( monomial_value ( n, o, x, expon ) )';
    err = abs ( quad - exact );
    fprintf ( 1, '  CN_LEG_03_1:    %6d  %14.6g  %14.6e\n', o, quad, err );

    [ o, x, w ] = cn_leg_03_xiu ( n );
    v = monomial_value ( n, o, x, expon );
    quad = w' * ( monomial_value ( n, o, x, expon ) )';
    err = abs ( quad - exact );
    fprintf ( 1, '  CN_LEG_03_XIU:  %6d  %14.6g  %14.6e\n', o, quad, err );

  end

  p = 5;

  if ( d <= p )

    if ( 4 <= n && n <= 6 )
      option = 1;
      [ o, x, w ] = cn_leg_05_1 ( n, option );
      v = monomial_value ( n, o, x, expon );
      quad = w' * ( monomial_value ( n, o, x, expon ) )';
      err = abs ( quad - exact );
      fprintf ( 1, '  CN_LEG_05_1(1): %6d  %14.6g  %14.6e\n', o, quad, err );
    end

    if ( 4 <= n && n <= 5 )
      option = 2;
      [ o, x, w ] = cn_leg_05_1 ( n, option );
      v = monomial_value ( n, o, x, expon );
      quad = w' * ( monomial_value ( n, o, x, expon ) )';
      err = abs ( quad - exact );
      fprintf ( 1, '  CN_LEG_05_1(2): %6d  %14.6g  %14.6e\n', o, quad, err );
    end

    if ( 2 <= n )
      [ o, x, w ] = cn_leg_05_2 ( n );
      v = monomial_value ( n, o, x, expon );
      quad = w' * ( monomial_value ( n, o, x, expon ) )';
      err = abs ( quad - exact );
      fprintf ( 1, '  CN_LEG_05_2:    %6d  %14.6g  %14.6e\n', o, quad, err );
    end

  end

  fprintf ( 1, '  EXACT:                  %14.6g\n', exact );

  return
end


text:
function [z_irm,z_ibm,t] = applyIdealMasks(xt,xi,nfft,hop,fs)
%APPLYIDEALMASKS Calculate and apply ideal masks via STFT
%
%   Z_IRM = IOSR.BSS.APPLYIDEALMASKS(XT,XI) calculates the ideal ratio mask
%   (IRM) and applies it to the mixture XT+XI, where XT is the target
%   signal and XI is the interference signal. The IRM is calculated via the
%   STFT using 1024-point windows with 512-point overlap. Z_IRM, XT, and XI
%   are vectors. If XT and XI are of different lengths then the shorter
%   signal is zero-padded in order to make them the same length; Z_IRM is
%   the same length as XT and XI.
%
%   Z_IRM = IOSR.BSS.APPLYIDEALMASKS(XT,XI,NFFT) uses NFFT-length segments
%   in the STFT.
%
%   Z_IRM = IOSR.BSS.APPLYIDEALMASKS(XT,XI,WINDOW) uses
%   LENGTH(WINDOW)-length segments in the STFT and applies WINDOW to each
%   segment.
%
%   Z_IRM = IOSR.BSS.APPLYIDEALMASKS(XT,XI,WINDOW,HOP) uses hop size HOP
%   for the STFT.
%
%   [Z_IRM,Z_IBM] = IOSR.BSS.APPLYIDEALMASKS(...) calculates the ideal
%   binary mask (IBM) and applies it to the mixture, returning the result
%   to Z_IBM.
%
%   [Z_IRM,Z_IBM,T] = IOSR.BSS.APPLYIDEALMASKS(XT,XI,WINDOW,HOP,FS) uses
%   sampling frequency FS to return the corresponding time T of each
%   element in Z_IRM and Z_IBM.
%
%   See also IOSR.DSP.STFT, IOSR.DSP.ISTFT, IOSR.BSS.IDEALMASKS,
%            IOSR.BSS.APPLYMASKS.

%   Copyright 2016 University of Surrey.

    %% check input

    % check signals
    assert(isvector(xt) && numel(xt)>1, 'iosr:applyIdealMasks:invalidXt', 'XT must be a vector')
    assert(isvector(xi) && numel(xi)>1, 'iosr:applyIdealMasks:invalidXi', 'XI must be a vector')

    % make equal length
    maxlength = max([length(xi) length(xt)]);
    xt = pad(xt,maxlength);
    xi = pad(xi,maxlength);

    % check nfft
    if nargin<3
        nfft = 1024;
    end

    % determine window
    if numel(nfft)>1
        win = nfft;
        assert(isvector(win), 'iosr:applyIdealMasks:invalidWin', 'WINDOW must be a vector')
        nfft = length(win);
    else
        assert(round(nfft)==nfft && nfft>0, 'iosr:applyIdealMasks:invalidNfft', 'NFFT must be a positive integer')
        win = hamming(nfft);
    end

    % check x length
    assert(length(xt)>=nfft, 'iosr:applyIdealMasks:invalidXt', 'XT must have at least NFFT samples')
    assert(length(xi)>=nfft, 'iosr:applyIdealMasks:invalidXi', 'XI must have at least NFFT samples')

    % determine hop
    if nargin<4
        hop = fix(nfft/2);
    else
        assert(isscalar(hop) & round(hop)==hop, 'iosr:applyIdealMasks:invalidHop', 'HOP must be an integer')
        assert(hop<=nfft && hop>0, 'iosr:applyIdealMasks:invalidHop', 'HOP must be less than or equal to NFFT, and greater than 0')
    end

    % determine fs
    if nargin<5
        fs = 1;
    else
        assert(isscalar(fs), 'iosr:applyIdealMasks:invalidFs', 'FS must be an scalar')
    end

    %% calculate outputs

    % STFTs of signals and mixture
    st = iosr.dsp.stft(xt,win,hop);
    si = iosr.dsp.stft(xi,win,hop);
    mix = iosr.dsp.stft(xt+xi,win,hop);

    % return ideal masks
    [irm,ibm] = iosr.bss.idealMasks(st,si);

    % apply IRM
    z_irm = iosr.bss.applyMask(mix,irm,nfft,hop,fs);
    z_irm = pad(z_irm,maxlength);

    % apply IBM
    if nargout>1
        z_ibm = iosr.bss.applyMask(mix,ibm,nfft,hop,fs);
        z_ibm = pad(z_ibm,maxlength);
    end

    % calculate t
    if nargout>2
        t = (0:length(z_irm)-1)./fs;
    end

end

function y = pad(x,dur)
%PAD Zero-pad a vector

    if length(x)<dur
        y = [x(:); zeros(dur-length(x),1)];
    else
        y = x(:);
    end

end

author: IoSR-Surrey

repo: MatlabToolbox


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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the target point positions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function targets = update_Target_Point_Positions(dt,current_time,targets)

%IDs = targets(:,1);                 % Stores Lag-Pt IDs in col vector
%xPts= targets(:,2);                 % Previous x-Values of x-Target Pts.
yPts= targets(:,3);                  % Previous y-Values of y-Target Pts.
%kStiffs = targets(:,4);             % Stores Target Stiffnesses
%N_target = length(targets(:,1));    % Gives total number of target pts!

%
% PUMPING PARAMETERS
%
freq = 2;     % pumping frequency
%G = 0.325;     % vertical gap between teeth
%occ = 0.95;   % occlusion total
a = 0.337;    % total vertical distance pump travels
N_top = 3276;  % # of points along top

%
% READ IN ORIGINAL yPT POSITIONS
%
yPTS = read_In_yPT_Positions('sawtooth.vertex');

%
% Create the pumping behavior
%
targets(1:N_top,3) = yPTS(1:N_top) - (a/2)*sin( 2*pi*freq*current_time );
%yPts(N_top+1:end) = yPts(N_top+1:end) + 0.95*(G/2)*sin( 2*pi*freq*current_time );

%
% Actually update target point positions
%
%targets(1:N_top,3) = yPts;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the # of vertex pts and all the vertex pts from the
%           .vertex file.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function PTS = read_In_yPT_Positions(struct_name)

filename = struct_name;  %Name of file to read in
fileID = fopen(filename);

% Read in the file, use 'CollectOutput' to gather all similar data together
% and 'CommentStyle' to to end and be able to skip lines in file.
C = textscan(fileID,'%f %f','CollectOutput',1);

fclose(fileID);     %Close the data file.

vertices = C{1};    %Stores all read in data in vertices (N+1,2) array

PTS = vertices(2:end,2);

author: nickabattista


text:
% This script evaluates the percentage of space time covered by
%alarms
%
re = [];

% Stefan Wiemer    4/95

report_this_filefun(mfilename('fullpath'));

abo = abo2;

for tre2 = min(abo(:,4)):0.1:max(abo(:,4)-0.1)
    abo = abo2;
    abo(:,5) = abo(:,5)* par1/365 + a(1,3);
    l = abo(:,4) >= tre2;
    abo = abo(l,:);
    l = abo(:,3) < tresh;
    abo = abo(l,:);
    hold on

    % space time volume covered by alarms
    if isempty(abo) == 1
        Va = 0;
    else
        Va = sum(pi*abo(:,3).^2)*iala;
    end

    % All space time
    [len, ncu] = size(cumuall);

    r = loc(3,:);
    %r = reshape(cumuall(len,:),length(gy),length(gx));
    %r=reshape(normlap2,length(yvect),length(xvect));
    l = r < tresh;
    V = sum(pi*r(l).^2*(teb-t0b));
    disp([' Zalarm = ' num2str(tre2)])
    disp([' =============================================='])
    disp([' Total space-time volume (R<Rmin):  ' num2str(V)])
    disp([' Space-time volume covered with alarms (R<Rmin):  ' num2str(Va)])
    disp([' Percent of total covered with alarms (R<Rmin):  ' num2str(Va/V*100) ' Percent' ])

    re = [re ; tre2 Va/V*100 ];
end   % for tre2

figure

matdraw
axis off

uicontrol('Units','normal',...
    'Position',[.0 .65 .08 .06],'String','Save ',...
     'Callback',{@calSave9, re(:,1), re(:,2)})

rect = [0.20,  0.10, 0.70, 0.60];
axes('position',rect)
hold on
pl = semilogy(re(:,1),re(:,2),'r');
set(pl,'LineWidth',1.5)
pl = semilogy(re(:,1),re(:,2),'ob');
set(pl,'LineWidth',1.5,'MarkerSize',10)
set(gca,'YScale','log')

set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.m,'FontWeight','bold',...
    'FontWeight','bold','LineWidth',1.5,...
    'Box','on')
grid

ylabel('Va/Vtotal in %')
xlabel('Zalarm ')
watchoff

author: CelsoReyes


text:
function varargout = DefaultArgs_ss(Args, DefArgs)
% auxillary function to replace argument check in the beginning and def. args assigment
% sets the absent or empty values of the Args (cell array, usually varargin)
% to their default values from the cell array DefArgs.
% Output should contain the actuall names of arguments that you use in the function

% e.g. : in function MyFunction(somearguments , varargin)
% calling [SampleRate, BinSize] = DefaultArgs(varargin, {20000, 20});
% will assign the defualt values to SampleRate and BinSize arguments if they
% are empty or absent in the varargin cell list
% (not passed to a function or passed empty)
if isempty(Args)
    Args ={[]};
end

% if iscell(Args) & isstr(Args{1}) & length(Args)==1
%     Args = Args{1};
% end

if ~iscell(DefArgs)
    DefArgs = {DefArgs};
end
nDefArgs = length(DefArgs);
nInArgs = length(Args);
%out = cell(nDefArgs,1);
if (nargout~=nDefArgs)
    error('number of defaults is different from assigned');
    %keyboard
end
for i=1:nDefArgs

    if (i>nInArgs | isempty(Args{i}))
        varargout(i) = {DefArgs{i}};
    else
        varargout(i) = {Args{i}};
    end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function mask=flatsegment(node,edge)
%
% mask=flatsegment(node,edge)
%
% decompose edge loops into flat segments alone arbitrary planes of the bounding box
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
% date: 2008/04/08
%
% this code is fragile: it can not handle curves with many co-linear
% nodes near the corner point
%
% input:
%    node:  x,y,z coordinates of each node of the mesh
%    edge:  input, a single vector separated by NaN, each segment
%           is a close-polygon consisted by node IDs
%
% output:
%    mask:  output, a cell, each element is a close-polygon
%           on x/y/z plane
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

idx=edge;
nn=length(idx);
val=zeros(nn,1);
for i=1:nn
  tet=mod(i:i+3,nn);
  tet(find(tet==0))=nn;
  val(i)=(abs(det([node(idx(tet),:),ones(4,1)]))>1e-5);
end

val(end+1:end+2)=val(1:2);
mask={};
oldend=1;
for i=1:nn
	if(val(i)==1&val(i+1)==1&val(i+2)==0)
            val(i+2)=2;
            mask{count}=idx(oldend:i+2);
            count=count+1;
            oldend=i+2;
        else
            mask{count}=[idx(oldend:end);mask{1}];
            break;
        end
end

repo: REST


text:
function [ fea, out ] = ex_planestress6( varargin )
%EX_PLANESTRESS6 Plane stress analysis of a pressure vessel.
%
%   [ FEA, OUT ] = EX_PLANESTRESS6( VARARGIN ) Model example for plain stress
%   approximation of a pressure vessel (annular cross section with symmetry).
%
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       E           scalar {207e9}         Modulus of elasticity
%       nu          scalar {0.27}          Poissons ratio
%       sfun        string {sflag1}        Shape function for displacements
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { ...
  'E',        207e9; ...
  'nu',       0.27; ...
  'sfun',     'sflag1'; ...
  'iplot',    1; ...
  'igrid',    1; ...
  'tol',      0.1; ...
  'fid',      1 };
[got,opt] = parseopt( cOptDef, varargin{:} );
fid       = opt.fid;

% Geometry and grid.
fea.sdim = { 'x' 'y' };   % Coordinate names.
fea.grid = ringgrid( 12, 216, 100e-3, 120e-3 );
fea.grid = delcells( fea.grid, selcells( fea.grid, '(x<=eps) | (y<=eps)') );
if( opt.igrid~=1 )
  fea.grid = quad2tri( fea.grid );
end
n_bdr = max(fea.grid.b(3,:));   % Number of boundaries.

% Problem definition.
fea = addphys( fea, @planestress );
fea.phys.pss.eqn.coef{1,end} = { opt.nu };
fea.phys.pss.eqn.coef{2,end} = { opt.E  };
fea.phys.pss.sfun            = { opt.sfun opt.sfun };

% Boundary conditions.
bctype = mat2cell( zeros(2,n_bdr), [1 1], ones(1,n_bdr) );
bctype{1,4} = 1;
bctype{2,3} = 1;
fea.phys.pss.bdr.coef{1,5} = bctype;

bccoef = mat2cell( zeros(2,n_bdr), [1 1], ones(1,n_bdr) );
bccoef{1,1} = '-nx*1e4';
bccoef{2,1} = '-ny*1e4';
fea.phys.pss.bdr.coef{1,end} = bccoef;

% Parse and solve problem.
fea       = parsephys( fea );
fea       = parseprob( fea );                          % Check and parse problem struct.
fea.sol.u = solvestat( fea, 'fid', fid, 'icub', 1+str2num(strrep(opt.sfun,'sflag','')) );   % Call to stationary solver.

% Postprocessing.
s_disp = fea.phys.pss.eqn.vars{2,end};
if( opt.iplot>0 )
  figure
  postplot( fea, 'surfexpr', s_disp )
  title( 'Total displacement' )
end

% Error checking.
s_sx     = fea.phys.pss.eqn.vars{5,end};
s_sy     = fea.phys.pss.eqn.vars{6,end};
s_sxy    = fea.phys.pss.eqn.vars{7,end};
s_sp1    = fea.phys.pss.eqn.vars{8,end};
s_sp3    = fea.phys.pss.eqn.vars{10,end};
s_ez     = fea.phys.pss.eqn.vars{13,end};
s_ep1    = fea.phys.pss.eqn.vars{15,end};
s_ep2    = fea.phys.pss.eqn.vars{16,end};
s_ep3    = fea.phys.pss.eqn.vars{17,end};
v_disp   = evalexpr( s_disp, [100e-3 120e-3-2*sqrt(eps);0 0]+sqrt(eps), fea )';
v_dref   = [2.809e-8 2.635e-8];
[v_sx(1),v_sx(2)] = minmaxsubd( s_sx, fea );
v_sxref  = [-10000 55454];
[v_sy(1),v_sy(2)] = minmaxsubd( s_sy, fea );
v_syref  = [-10000 55454];
[v_sxy(1),v_sxy(2)] = minmaxsubd( s_sxy, fea );
v_sxyref = [-32730 0];
[v_sp1(1),v_sp1(2)] = minmaxsubd( s_sp1, fea );
v_sp1ref = [4.5e4 55454];
[v_sp3(1),v_sp3(2)] = minmaxsubd( s_sp3, fea );
v_sp3ref = [-1e4 0];
[v_ez(1),v_ez(2)] = minmaxsubd( s_ez, fea );
v_ezref = [-5.929e-8 -5.929e-8];
[v_ep1(1),v_ep1(2)] = minmaxsubd( s_ep1, fea );
v_ep1ref = [2.196e-7 2.809e-7];
[v_ep2(1),v_ep2(2)] = minmaxsubd( s_ep2, fea );
v_ep2ref = [-5.929e-8 -5.929e-8];
[v_ep3(1),v_ep3(2)] = minmaxsubd( s_ep3, fea );
v_ep3ref = [-1.206e-7 -5.929e-8];
out.err  = [ abs([v_dref-v_disp])./v_dref ;
             abs([v_sxref-v_sx])./v_sxref ;
             abs([v_syref-v_sy])./v_syref ;
             abs([v_sxyref(1)-v_sxy(1)])./v_sxyref(1) 0 ;
             abs([v_sp1ref(2)-v_sp1(2)])./v_sp1ref(2) 0 ;
             abs([v_sp3ref(1)-v_sp3(1)])./v_sp3ref(1) 0 ;
             abs([v_ezref(1)-v_ez(1)])./v_ezref(1) 0 ;
             abs([v_ep1ref(1)-v_ep1(1)])./v_ep1ref(1) 0 ;
             abs([v_ep2ref(1)-v_ep2(1)])./v_ep2ref(1) 0 ;
             abs([v_ep3ref(1)-v_ep3(1)])./v_ep3ref(1) 0 ];
out.pass = all( out.err(:) <= opt.tol );

if( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
function keypointsAll = scoremaps2keypoints(expidx, image_set)

p = exp_params(expidx);

pad_orig = p.([image_set 'Pad']);
stride = p.stride;
inv_scale_factor = 1/p.scale_factor;
locref = p.locref;
scale_mul = sqrt(53);

pidxs = p.pidxs;
[~,parts] = util_get_parts24();

load(p.testGT,'annolist');

keypointsAll = [];

min_det_score = 1e-03;

save_dir = fullfile(p.exp_dir, image_set);

for imgidx = 1:length(annolist)
    fprintf('imgidx: %d/%d\n', imgidx, length(annolist));
    im_fn = annolist(imgidx).image.name;
    [~,im_name,~] = fileparts(im_fn);
    load(fullfile(p.unary_scoremap_dir, im_name), 'scoremaps');
    if locref
        load(fullfile(p.unary_scoremap_dir, [im_name '_locreg']), 'locreg_pred');
    end

    size2d = size(scoremaps(:,:,1));
    num_proposals = size(scoremaps, 1) * size(scoremaps, 2);
    num_channels = size(scoremaps, 3);
    unPosAll = zeros(num_proposals, 2);
    unProbAll = zeros(num_proposals, num_channels);
    if locref
        locationRefine = zeros(num_proposals, num_channels, 2);
    end

    for k = 1:num_proposals
        idx = k;
        [row, col] = ind2sub(size2d, idx);
        % transform heatmap to image coordinates
        im_y = (row-1)*stride;
        im_x = (col-1)*stride;
        unPosAll(k, :) = [pad_orig, pad_orig] + [im_x, im_y]*inv_scale_factor;
        unProbAll(k,:) = scoremaps(row,col,:);
        if locref
            locationRefine(k, :, :) = squeeze(locreg_pred(row, col, :, :))*scale_mul*inv_scale_factor;
        end
    end

    for i = 1:length(pidxs)
        pidx = pidxs(i);
        % part is a joint
        assert(parts(pidx+1).pos(1) == parts(pidx+1).pos(2));
        jidx = parts(pidx+1).pos(1);

        unPos = unPosAll;
        if locref
            unPos = unPos + squeeze(locationRefine(:,i,:));
        end
        unProb = unProbAll(:,i);

        I = unProb >= min_det_score;
        unProb = unProb(I,:);
        unPos = unPos(I,:);
        [unProb,I] = sort(unProb, 'descend');
        unPos = unPos(I,:);

        keypointsAll(imgidx).det{jidx+1} = [unPos unProb];
    end
end

mkdir_if_missing(save_dir);
save(fullfile(save_dir, 'keypointsAll'), 'keypointsAll');

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
%
%CA driver
%
%forest fire

clf
clear all

n=100;

Plightning = .000005;
Pgrowth = .01; %.01

z=zeros(n,n);
o=ones(n,n);
veg=z;
sum=z;

imh = image(cat(3,z,veg*.02,z));
set(imh, 'erasemode', 'none')
axis equal
axis tight

% burning -> empty
% green -> burning if one neigbor burning or with prob=f (lightning)
% empty -> green with prob=p (growth)
% veg = {empty=0 burning=1 green=2}
for i=1:3000
    %nearby fires?

     sum = (veg(1:n,[n 1:n-1])==1) + (veg(1:n,[2:n 1])==1) + ...
           (veg([n 1:n-1], 1:n)==1) + (veg([2:n 1],1:n)==1) ;

    veg = ...
         2*(veg==2) - ((veg==2) & (sum>0 | (rand(n,n)<Plightning))) + ...
         2*((veg==0) & rand(n,n)<Pgrowth) ;

    set(imh, 'cdata', cat(3,(veg==1),(veg==2),z) )
    drawnow
end

author: Eurus-Holmes


text:
function dn = Dmult(dn0,dn1)

% DMULT dual number multiplication
%
%   DN = DMULT(DN0,DN1) returns the dual number DN which is the dual number
%     multiplication of the dual numbers DN0 and DN1
%      - DN0 (resp. DN1) is a dual number (DN0 = a0 + eps*b0, eps^2 = 0).
%         It is a 2-vector or a 2*N array (column i represents dual number
%         i) where N is the number of dual numbers. DN0 and DN1 must have
%         the same size.
%      - DN is a dual number. It is a 2*N array (each column is the dual
%          multiplication of the corresponding columns in DN0 and DN1)
%
% See also DQINV

s0 = size(dn0);
s1 = size(dn1);
if s0 == [1 2], dn0 = dn0';s0 = size(dn0);end
if s1 == [1 2], dn1 = dn1';s1 = size(dn1);end

% wrong format
if s0(1) ~= 2 || s1(1) ~= 2
    error('DualQuaternion:Dmult:wrongsize',...
        '%d rows in the DN0 array and %d rows in the DN1 array. It should be 2 for both.',...
        s0(1),s1(1));
end

% sizes do not match
n1 = s0(2);
n2 = s1(2);
if n1 ~= n2
    error('DualQuaternion:Dmult:wrongsize',...
        '%d dual numbers in DN0 array and %d dual numbers in DN1 array.They should be equal.',...
        n1,n2);
end

dn = zeros(2,n1);
dn(1,:) = dn0(1,:).*dn1(1,:);
dn(2,:) = dn0(1,:).*dn1(2,:)+dn0(2,:).*dn1(1,:);

author: Sable


text:
function mask = coord2mask(M,N,vx,vy)
%COORD2MASK Generates a binary mask from given coordinates.
%   MASK = COORD2MASK(M,N,vx,yy) generates a binary mask, MASK, of size
%   M-by-N given a set of sequential coordinates [vx,yy] describing the
%   vertices of a polygon. The points inside and on the polygon are
%   labeled with 1s and all other points are labeled 0.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Create the mask by labeling as 1 all point inside or on the polygon
% defined by the vertices [vx,vy]. All other points in mask are labeled
% 0. This labeling is done automatically by MATLAB function poly2mask.
% NOTE: The syntax of this function is awkward. The order of vx (rows)
% and vy (columns) coordinates have to be reversed to match the book,
% but the M (number of rows) and N (number of columns) are input in the
% order used in the book.
mask = poly2mask(vy,vx,M,N);

author: dipum

repo: dipum-toolbox

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


text:
function value = r8_aint ( x )

%****************************************************************************80
%
%% R8_AINT truncates an R8 argument to an integer.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 October 2011
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the truncated version of X.
%
  if ( x < 0.0 )
    value = - floor ( abs ( x ) );
  else
    value =   floor ( abs ( x ) );
  end

  return
end


text:
function definput = arg_fwtcommon(definput)

%definput.flags.ansy = {'syn','ana'};
definput.keyvals.a = [];

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
  function strum_test
%|function strum_test
%| test strum object
%| Copyright 2006-1-19, Jeff Fessler, University of Michigan

dat.a = 'x';
dat.b = 'y';

imp10 = @(st) st.b;
imp11 = @(st, arg) [st.a arg];
imp12 = @(st, arg1, arg2) sum(arg1) + sum(arg2);

list = { ...
	'fun10', @strum_test_fun10, '()'; ...
	'fun11', @strum_test_fun11, '(arg)'; ...
	'fun12', @strum_test_fun12, '(arg1, arg2)'; ...
	'imp11', imp11, '(arg)'; ...
	'imp12', imp12, '(arg1, arg2)'; ...
};

% no longer attempting anything but exactly 1 output so comment out these:
%	'fun00', @strum_test_fun00, '()'; ...
%	'fun22', @strum_test_fun22, '(arg1, arg2)'; ...

for ii=1:2
	if ii==1
		st = strum(dat, list); % with comment
	else
		tmp = col(list(:,1:2)');
		st = strum(dat, tmp); % without comment
	end
	st.a;
	st.b;
	jf_equal(st.fun10(), dat.a)
	jf_equal(st.fun11('d'), [dat.a 'd'])
	jf_equal(st.fun12('d', 'e'), [dat.a 'd' 'e'])
	jf_equal(st.imp11('d'), imp11(dat, 'd'))
	jf_equal(st.imp12(1:2, 3:4), imp12(dat, 1:2, 3:4))

	% augment base strum with more data and methods
	new.c = {10, 20};
	new.d = 'z';
	if ii==1
		s2 = strum(new, {'imp10', imp10, 'imp10()'}, 'base', st);
	else
		s2 = strum(new, {'imp10', imp10}, 'base', st);
	end
	jf_equal(s2.d, 'z')
	jf_equal(s2.imp10(), st.b)

	% cell arguments
	new.c{1};
	s2.c{1};
	jf_equal(s2.c{1}, new.c{1})
	jf_equal(s2.c{2}, new.c{2})

%	s2.c{1:2} % fails due to matlab limitation
%	s2.c{:} % ""
	tmp = s2.c; jf_equal(new.c, tmp) % this is the workaround
end

if 0 % old attempt to see if two outputs would work (it does not)
	st = strum(dat, {'fun22', @strum_test_fun22});
	jf_equal(st.fun22('c', 'd'), 'xc')
%	[a b] = st.fun22('c', 'd'); % "Too many output arguments" says matlab

	try
		[a b] = st.fun22(3, 4)
		printm 'multiple outputs worked!'
	catch
		warn 'darn, matlab cannot handle multiple outputs, as matt said'
	end
end

function strum_test_fun00(ob) % 0 out, 0 in
printm 'ok'

function a = strum_test_fun10(ob) % 1 out, 0 in
a = ob.a;

function out = strum_test_fun11(ob, arg) % 1 out, 1 in
out = [ob.a arg];

function out = strum_test_fun12(ob, arg1, arg2) % 1 out, 2 in
out = [ob.a arg1 arg2];

function [a, b] = strum_test_fun22(ob, arg1, arg2) % 2 out, 2 in
a = [ob.a arg1];
b = [ob.b arg2];

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [mGal] = nmps22mGal(nmps2)
% Convert acceleration from nanometers per square-second to milligalileos
% Chad A. Greene 2012
mGal = nmps2*1e-4;


text:
% The COBRAToolbox: testFindRxnsFromMets.m
%
% Purpose:
%     - Tests the functionality of findRxnFromMets with all available
%     parameters
%
% Authors:
%     - Original file: Thomas Pfau Jan 2018
%

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testFindRxnsFromMets'));
cd(fileDir);

% load model
model = getDistributedModel('ecoli_core_model.mat');

% First, find all reactions involved with 3pg
singleTestMet = model.mets{3};
involvedReacs = {'Biomass_Ecoli_core_w_GAM';'PGK';'PGM'};
consuming = involvedReacs;
producing = involvedReacs(2:3);
reacs = findRxnsFromMets(model,singleTestMet);
assert(isempty(setxor(involvedReacs,reacs)));
reacs = findRxnsFromMets(model,singleTestMet,'consumersOnly',true);
assert(isempty(setxor(consuming,reacs)));
reacs = findRxnsFromMets(model,singleTestMet,'producersOnly',true);
assert(isempty(setxor(producing,reacs)));

% now, we change the bounds of PGK to not function in the forward direction
modelChanged = changeRxnBounds(model,'PGK',0,'l');
modelChanged = changeRxnBounds(modelChanged,'PGM',0,'u');
reacs = findRxnsFromMets(modelChanged,singleTestMet,'producersOnly',true);
assert(isempty(reacs));

% now, test the same with consumers and the PGM reaction, only biomass is
% left.
modelChanged = changeRxnBounds(model,'PGM',0,'l');
modelChanged = changeRxnBounds(modelChanged,'PGK',0,'u');
reacs = findRxnsFromMets(modelChanged,singleTestMet,'consumersOnly',true);
assert(isempty(setxor(involvedReacs(1),reacs)));

% lets test 2 metabolites (2pg and 3pg)
dualTestMet = model.mets(2:3);
involvedReacs = {'Biomass_Ecoli_core_w_GAM';'ENO';'PGM';'PGK'};
consuming = involvedReacs;
producing = involvedReacs(2:4);
reacs = findRxnsFromMets(model,dualTestMet);
assert(isempty(setxor(involvedReacs,reacs)));
reacs = findRxnsFromMets(model,dualTestMet,'consumersOnly',true);
assert(isempty(setxor(consuming,reacs)));
reacs = findRxnsFromMets(model,dualTestMet,'producersOnly',true);
assert(isempty(setxor(producing,reacs)));

% test printout.
% build comparison text
res = printRxnFormula(model,'PGM');
compText = res{1}; % this is sufficient;

% create diary
diaryFile = 'RxnsFromMetTest';
diary(diaryFile)
reacs = findRxnsFromMets(model,dualTestMet,'containsAll',true,'printFlag',1);
assert(isempty(setxor({'PGM'},reacs)));
diary off
text = importdata(diaryFile);
assert(~isempty(strfind(strrep(text,'\n',''),compText)));
%cleanup
delete(diaryFile);

% test implicit printOut is silent, if no output is created
diary(diaryFile)
[reacs,forms] = findRxnsFromMets(model,dualTestMet,'containsAll',true);
diary off;
text = importdata(diaryFile);
% ITs not found in the output
assert(isempty(strfind(text,compText)));
% but is equal to printRxnFormula
assert(isequal(forms,res))
%cleanup
delete(diaryFile);

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function ind = buscamin(x);
% It searches the first minimum of the modulus of x.

x = abs(x);
localmin = (x(2:end-1)<=x(1:end-2))...
         & (x(2:end-1)<=x(3:end));

ind = min(find(localmin));

author: marianux

repo: ecg-kit

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


text:
function sparse_grid_composite_test06 ( dim_num, level_max )

%*****************************************************************************80
%
%% SPARSE_GRID_COMPOSITE_TEST06 creates a sparse composite grid and writes it to a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the level.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_COMPOSITE_TEST06:\n' );
  fprintf ( 1, '  SPARSE_GRID_COMPOSITE makes a sparse composite grid.\n' );
  fprintf ( 1, '  Write the data to a set of quadrature files.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEVEL_MAX = %d\n', level_max );
  fprintf ( 1, '  Spatial dimension DIM_NUM = %d\n', dim_num );
%
%  Determine the number of points.
%
  point_num = sparse_grid_composite_size ( dim_num, level_max );

  r(1:dim_num,1) = -1.0;
  r(1:dim_num,2) = +1.0;
%
%  Compute the weights and points.
%
  [ w, x ] = sparse_grid_composite ( dim_num, level_max, point_num );
%
%  Write the data out.
%
  r_filename = sprintf ( 'composite_d%d_level%d_r.txt', dim_num, level_max );
  w_filename = sprintf ( 'composite_d%d_level%d_w.txt', dim_num, level_max );
  x_filename = sprintf ( 'composite_d%d_level%d_x.txt', dim_num, level_max );

  r8mat_write ( r_filename, dim_num, 2,         r );
  r8mat_write ( w_filename, 1,       point_num, w );
  r8mat_write ( x_filename, dim_num, point_num, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  R data written to "%s".\n', r_filename );
  fprintf ( 1, '  W data written to "%s".\n', w_filename );
  fprintf ( 1, '  X data written to "%s",\n', x_filename );

  return
end

author: johannesgerer


text:
function [R, t, X] = andreff(AA,BB)
% Solves the problem AX=XB
% using the formulation of
%
% On-line Hand-Eye Calibration.
% N. Andreff, R. Horaud, B. Espiau
%
% Mili Shah
% July 2014

[~,n] = size(AA); n = n/4;

A = zeros(12*n,12);
b = zeros(12*n,1);
for i = 1:n
    Ra = AA(1:3,4*i-3:4*i-1);
    Rb = BB(1:3,4*i-3:4*i-1);
    ta = AA(1:3,4*i);
    tb = BB(1:3,4*i);
    A(12*i-11:12*i-3,1:9) = eye(9) - kron(Rb,Ra);
    A(12*i-2:12*i,:) = [kron(tb',eye(3)) eye(3)-Ra];
    b(12*i-2:12*i) = ta;
end
x = pinv(A) * b;

X = reshape(x(1:9),3,3)';
X = sign(det(X))/abs(det(X))^(1/3)*X;

[u, ~, v] = svds(X, 3); X = u*v'; if det(X)<0, X = u*diag([1 1 -1])*v'; end
X = [X' x(10:12);[0 0 0 1]];
R = X(1 : 3, 1 : 3);
t = X(1 : 3, 4);

end


text:
function [infos, f_val, optgap, grad, gnorm, subgrad, subgnorm, smooth_grad, smooth_gnorm] = store_infos(problem, w, options, infos, epoch, grad_calc_count, elapsed_time)
% Function to store statistic information
%
% Inputs:
%       problem         function (cost/grad/hess)
%       w               solution
%       options         options
%       infos           struct to store statistic information
%       epoch           number of outer iteration
%       grad_calc_count number of calclations of gradients
%       elapsed_time    elapsed time from the begining
% Output:
%       infos           updated struct to store statistic information
%       f_val           cost function value
%       outgap          optimality gap
%       grad            gradient
%       gnorm           norm of gradient
%       subgrad         subgradient
%       subgnorm        norm of subgradient
%       smooth_grad     smoothed gradient
%       smooth_gnorm    norm of smoothed gradient
%
% This file is part of SGDLibrary.
%
% Created by H.Kasai on Sep. 25, 2017
% Modified by H.Kasai on Mar. 27, 2017
% Modified by H.Kasai on Oct. 30, 2020

    subgrad = [];
    subgnorm = [];

    if ~epoch

        infos.iter = epoch;
        infos.time = 0;
        infos.grad_calc_count = grad_calc_count;
        f_val = problem.cost(w);
        optgap = f_val - options.f_opt;
        % calculate norm of full gradient
        grad = problem.full_grad(w);
        gnorm = norm(grad);
        if ismethod(problem, 'full_subgrad')
            subgrad = problem.full_subgrad(w);
            subgnorm = norm(subgrad);
            infos.subgnorm = subgnorm;
        end

        if ismethod(problem, 'full_smooth_grad')
            smooth_grad = problem.full_smooth_grad(w);
            smooth_gnorm = norm(smooth_grad);
            infos.smooth_gnorm = smooth_gnorm;
        end

        infos.optgap = optgap;
        infos.best_optgap = optgap;
        infos.absoptgap = abs(optgap);
        infos.gnorm = gnorm;
        infos.cost = f_val;
        infos.best_cost = f_val;
        if ismethod(problem, 'reg')
            infos.reg = problem.reg(w);
        end
        if options.store_w
            infos.w = w;
        end
        if options.store_grad
            infos.grad = grad;
        end

    else

        infos.iter = [infos.iter epoch];
        infos.time = [infos.time elapsed_time];
        infos.grad_calc_count = [infos.grad_calc_count grad_calc_count];

        % calculate optimality gap
        f_val = problem.cost(w);
        optgap = f_val - options.f_opt;
        % calculate norm of full gradient
        grad = problem.full_grad(w);
        gnorm = norm(grad);
        if ismethod(problem, 'full_subgrad')
            subgrad = problem.full_subgrad(w);
            subgnorm = norm(subgrad);
            infos.subgnorm = [infos.subgnorm subgnorm];
        end

        if ismethod(problem, 'full_smooth_grad')
            smooth_grad = problem.full_smooth_grad(w);
            smooth_gnorm = norm(smooth_grad);
            infos.smooth_gnorm = [infos.smooth_gnorm smooth_gnorm];
        end

        infos.optgap = [infos.optgap optgap];
        infos.absoptgap = [infos.absoptgap abs(optgap)];
        if optgap < infos.best_optgap(end)
            infos.best_optgap = [infos.best_optgap optgap];
        else
            infos.best_optgap = [infos.best_optgap infos.best_optgap(end)];
        end
        infos.cost = [infos.cost f_val];
        if f_val < infos.best_cost(end)
            infos.best_cost = [infos.best_cost f_val];
        else
            infos.best_cost = [infos.best_cost infos.best_cost(end)];
        end
        infos.gnorm = [infos.gnorm gnorm];
        if ismethod(problem, 'reg')
            reg = problem.reg(w);
            infos.reg = [infos.reg reg];
        end
        if options.store_w
            infos.w = [infos.w w];
        end
        if options.store_grad
            infos.grad = [infos.grad grad];
        end

    end

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
%  Figure 7.55      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
% script for Fig. 7.55
%
clf;
dp=[1 1 0];
np=[1];
nc=conv([1 1],[8.32 0.8]);
dc=conv([1 4.08],[1 0.0196]);
num=conv(np,nc);
den=conv(dp,dc);
dcl=[0 0 num]+den; % closed-loop denominator
t=0:.1:5;
sys=tf(num,dcl);
y=step(sys,t);
plot(t,y)
grid;
xlabel('Time (sec)');
ylabel('y');
title('Fig.7.55 Step response for lag compensation design');
[a,b,c,d]=tf2ss(nc,dc);
% Samplig period
ts=0.1;
% Convert controller to discrete
[ad,bd,cd,dd]=c2dm(a,b,c,d,ts,'zoh');
[ncd,dcd]=ss2tf(ad,bd,cd,dd);
% response comparison of continuous and digital control.
sim('fig7_54')
figure(1)
plot(ycd(:,1),ycd(:,2))
hold on
plot(ycd(:,1),ycd(:,3),'g--')
title('Figure 7.53(a) Output Responses of Digital and Continuous Controllers')
xlabel('Time (sec)');
ylabel('y');
text(0.5,0.5,'continuous controller')
text(0.1,1.15,'digital controller')
nicegrid
hold off
%pause;
figure(2)
plot(ucd(:,1),ucd(:,2))
hold on
plot(ucd(:,1),ucd(:,3),'g--')
title('Figure 7.55(b)  Control Responses of Digital and Continuous Controllers')
xlabel('Time (sec)');
ylabel('u');
text(0.5,0.5,'continuous controller')
text(1.1,-1.5,'digital controller')
nicegrid
hold off

repo: mcbench-benchmarks


text:
function n = p22_n ( )

%*****************************************************************************80
%
%% P22_N returns the number of variables for problem 22.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 December 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N.  If N is positive, N represents the
%    only legal value for N for this problem.  If N is
%    negative, the absolute value of N is the least legal
%    value of N, but other values are allowable.
%
  n = 3;

  return
end


text:
classdef LevelSetWithSeveralHoles < LevelSetCreator

    properties (Access = private)
        nHoles
        rHoles
        phaseHoles
    end

    methods (Access = public)

        function obj = LevelSetWithSeveralHoles(cParams)
            obj.init(cParams);
            obj.compute(cParams);
        end

    end

    methods (Access = protected)

        function computeLevelSet(obj)
            obj.computeLevelSetValue();
        end

    end

    methods (Access = private)

        function computeLevelSetValue(obj)
            ls = ones(obj.lsSize);
            for idim = 1:obj.ndim
                coordV = obj.nodeCoord(:,idim);
                cosDim = obj.computeDirectionalCosinus(coordV,idim);
                ls = ls.*cosDim;
            end
            ls = ls + obj.rHoles-1;
            obj.levelSet = ls;
        end

        function init(obj,cParams)
            obj.nHoles     = cParams.nHoles;
            obj.rHoles     = cParams.rHoles;
            obj.phaseHoles = cParams.phaseHoles;
        end

        function cosDir  = computeDirectionalCosinus(obj,coord,dir)
            pos = coord;
            l = max(pos) - min(pos);
            n = obj.nHoles(dir);
            fase = obj.phaseHoles(dir);
            cosDir = cos((n + 1)*pos*pi/l + fase);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function ieFontInit(fig)
%Initialize the font size based on the preference setting
%
%   ieFontInit(fig)
%
% Font sizes are managed using the preference mechanism in Matlab. The font
% size is a field, fontDelta, in the mrVista preference structure. That
% information is retrieved here and the font size is changed based on this
% setting. This setting is stored across sessions.
%

vistaPref = getpref('mrVista');
if ~isempty(vistaPref)
    if checkfields(vistaPref,'fontDelta')
        fontDelta = getpref('vistaPref','fontDelta');
        ieFontChangeSize(fig,fontDelta);
    end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [irf_record,D_record,gamma_record,ETA_record,beta_gibbs,sigma_gibbs]...
    =IRFt6_Bayesian(betahat,IRFperiods,n,m,p,k,T,names,startdate,enddate,X,FEVDresperiods,Y,pref,IRFt,arvar,q,It,Bu,lambda1,lambda3,lambda4,strctident,favar)

%%Implementation of a Bayesian Proxy VAR based on the Codes published by Caldara and Herbst (2018)
%%Implementation by Ben Schumann

%% IV identification
% Load IV and make it comparable with the reduced form errors
[EPSIV,IVcut,~,sigmahatIV,sigma_hat,inv_sigma_hat,IV,txt,~,cut1,cut2,cut3,cut4]=...
    bear.loadIV(betahat,k,n,Y,X,T,p,names,startdate,enddate,strctident,pref);
% IV routine
[beta_draws,sigma_draws,IV_draws,C_draws]=...
    bear.irfIV_MH(EPSIV,IVcut,betahat,sigmahatIV,sigma_hat,inv_sigma_hat,IV,txt,cut1,cut2,cut3,cut4,names,It,Bu,n,arvar,lambda1,lambda3,lambda4,m,p,k,q,X,Y,T,startdate,enddate,pref,strctident,IRFperiods,IRFt);

%reorganize
%finaly tell bear how many draws we stored
It=size(IV_draws,2);
Bu=0; %set burn in to 0, as burn in is no longer requiered

%% check Restrictions
[irf_record,D_record,gamma_record,ETA_record,beta_gibbs,sigma_gibbs]...
        =bear.irfres(beta_draws,sigma_draws,C_draws,IV_draws,IRFperiods,n,m,p,k,T,Y,X,FEVDresperiods,strctident,pref,favar,IRFt,It,Bu);

author: european-central-bank

repo: BEAR-toolbox


text:
function DIC_3Djoined_results = join3DreconstructedPairs(DIC3DAllPairsResults)
%% function for joining multiple 3D-DIC pair results into one structure (plotting in step3)
%
% Inputs:
% * DIC3DAllPairsResults
%
% Outputs:
% * DIC_3Djoined_results
%

%%
nPairs=numel(DIC3DAllPairsResults);
nFrames=size(DIC3DAllPairsResults{1}.Points3D,1);

% create cells with all the variables for joining
Fi=cell(nPairs,1);
FCi=cell(nPairs,1);
Vi=cell(nPairs,nFrames);
corrCombi=cell(nPairs,nFrames);
FaceCorrCombi=cell(nPairs,nFrames);
DispVi=cell(nPairs,nFrames);
DispMi=cell(nPairs,nFrames);
VFci=cell(nPairs,nFrames);

for ip=1:nPairs
    DIC_3Dpair_struct=DIC3DAllPairsResults{ip};

    Fi{ip}=DIC_3Dpair_struct.Faces;
    FCi{ip}=DIC_3Dpair_struct.FaceColors;

    for it=1:nFrames
        Vi{ip,it}=DIC_3Dpair_struct.Points3D{it};
        corrCombi{ip,it}=DIC_3Dpair_struct.corrComb{it};
        FaceCorrCombi{ip,it}=DIC_3Dpair_struct.FaceCorrComb{it};
        DispVi{ip,it}=DIC_3Dpair_struct.Disp.DispVec{it};
        DispMi{ip,it}=DIC_3Dpair_struct.Disp.DispMgn{it};
        VFci{ip,it}=DIC_3Dpair_struct.FaceCentroids{it};
    end

end

% join variables that are time invariant
[Fj,~,~]=joinElementSets(Fi,Vi(:,1),[]);
[~,~,FCj]=joinElementSets(Fi,Vi(:,1),FCi);

for it=1:nFrames
    % join variables that are time variant and scalar

    [~,Vj{1,it},DispMj{1,it}]=joinElementSets(Fi,Vi(:,it),DispMi(:,it));
    [~,~,corrCombj{1,it}]=joinElementSets(Fi,Vi(:,it),corrCombi(:,it));
    [~,~,FaceCorrCombj{1,it}]=joinElementSets(Fi,Vi(:,it),FaceCorrCombi(:,it));

    % join variables that are time variant and vectors
    DispVj{1,it}=DispVi{1,it};
    VFcj{1,it}=VFci{1,it};

    for ip=2:nPairs
        DispVj{1,it}=[DispVi{1,it}; DispVi{ip,it}];
        VFcj{1,it}=[VFcj{1,it}; VFci{ip,it}];
    end
end

DIC_3Djoined_results=struct;
DIC_3Djoined_results.Faces=Fj;
DIC_3Djoined_results.FaceColors=FCj;
DIC_3Djoined_results.FaceCentroids=VFcj;
DIC_3Djoined_results.Points3D=Vj;
DIC_3Djoined_results.corrComb=corrCombj;
DIC_3Djoined_results.FaceCorrComb=FaceCorrCombj;

DIC_3Djoined_results.Disp.DispVec=DispVj;
DIC_3Djoined_results.Disp.DispMgn=DispMj;

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function test_issue1400

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_selectdata

ntrial = 2;

data = [];
data.label = {'1'};
for i=1:ntrial
  t = (1:1000)/1000 + i; % time vectors are unaligned to show the error
  s = sin(i*2*pi*t); % make a sine with 1, 2, 3, ... Hz
  data.time{i} = t;
  data.trial{i}(1,:) = s;  % insert it in the first channel
end

disp('Demonstrate averaging over time:')
disp('Averaging trials 1, 2')
cfg = [];
cfg.trials = [1,2];
cfg.avgoverrpt = 'yes';
cfg.keeprpt = 'no';
cfg.verbosity = 'off';
tmpdata = ft_selectdata(cfg, data);

% disp(['trialinfo(1, 1) = ', num2str(tmpdata.trialinfo(1,1))])
disp(['trial 1 time = ', regexprep(num2str(data.time{1}(1:5)),'\s+',', ')])
disp(['trial 2 time = ', regexprep(num2str(data.time{2}(1:5)),'\s+',', ')])
disp(['average time = ', regexprep(num2str(tmpdata.time{1}(1:5)),'\s+',', ')])

% add some trialinfo which is a mix of numbers and letters
disp('Demonstrate averaging over trialinfo:')
data.trialinfo = table([1;2],{'a'; 'b'});
disp(['trialinfo class is: ', class(data.trialinfo)])

% This would return an error in ft_selectdata, because it
% averaging the table prior to discarding the averaged trialinfo
disp('Averaging trials 1, 2')
cfg = [];
cfg.trials = [1,2];
cfg.avgoverrpt = 'yes';
cfg.keeprpt = 'no';
cfg.verbosity = 'off';
tmpdata = ft_selectdata(cfg, data);

author: fieldtrip

repo: fieldtrip


text:
% RGBTOGRAY.m
% -------------------------------------------------------------------
%
% Date:    27/04/2013
% Last modified: 27/04/2013
% -------------------------------------------------------------------

function gray = RGBTOGRAY(img)

    if size(img, 3) ~= 3,
        error('The input should be RGB');
    end

    if ~isinteger(img) && max(img(:)) > 1,
        gray = rgb2gray(uint8(img));
        gray = double(gray);
    else
        gray = rgb2gray(img);
    end

end

author: thfylsty


text:
function [model, time] = rbf_build(Xtr, Ytr, bf_type, bf_c, usePolyPart, verbose)
% RBFBUILD
% Builds a Radial Basis Function (RBF) interpolant using training data
%
% Call
%   [model, time] = rbf_build(Xtr, Ytr, bf_type, bf_c, usePolyPart, verbose)
%   [model, time] = rbf_build(Xtr, Ytr, bf_type, bf_c, usePolyPart)
%   [model, time] = rbf_build(Xtr, Ytr, bf_type, bf_c)
%   [model, time] = rbf_build(Xtr, Ytr, bf_type)
%   [model, time] = rbf_build(Xtr, Ytr)
%
% Input
% Xtr, Ytr    : Training data points (Xtr(i,:), Ytr(i)), i = 1,...,n
%               Note that the input variables must be scaled to e.g. [0,1]
%               or [-1,1] for better predictive performance.
% bf_type     : Type of the basis functions (default = 'MQ'):
%               'BH' = Biharmonic
%               'MQ' = Multiquadric
%               'IMQ' = Inverse Multiquadric
%               'TPS' = Thin plate spline
%               'G' = Gaussian
% bf_c        : Parameter c value (default = 1)
% usePolyPart : Use also the polynomial term P of the model Y = P + RBF
%               (default = 0, do not use)
% verbose     : Set to 0 for no verbose (default = 1)
%
% Output
% model     : RBF model - a struct with the following elements:
%    n      : Number of data points in the training data set
%    meanY  : Mean of Ytr
%    bf_type: Type of the basis functions
%    bf_c   : Parameter c value
%    poly   : Use also the polynomial term
%    coefs  : Coefficients of the model
% time      : Execution time
%
% Please give a reference to the software web page in any publication
% describing research performed using the software, e.g. like this:
% Jekabsons G. Radial Basis Function interpolation for Matlab, 2009,
% available at http://www.cs.rtu.lv/jekabsons/

% This source code is tested with Matlab version 7.1 (R14SP3).

% =========================================================================
% RBF interpolation
% Version: 1.1
% Date: August 12, 2009
% Author: Gints Jekabsons (gints.jekabsons@rtu.lv)
% URL: http://www.cs.rtu.lv/jekabsons/
%
% Copyright (C) 2009  Gints Jekabsons
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% =========================================================================

if nargin < 2
    error('Too few input arguments.');
else

    [n, d] = size(Xtr);
    [ny, dy] = size(Ytr);
    if (n < 2) || (d < 1) || (ny ~= n) || (dy ~= 1)
        error('Wrong training data sizes.');
    end

    if nargin < 3
        bf_type = 'MQ';
    end
    if nargin < 4
        bf_c = 1;
    end
    if nargin < 5
        usePolyPart = 1;
    end
    if nargin < 6
        verbose = 0;
    end

    tic;

    model.n = n;
    model.meanY = mean(Ytr);
    model.bf_type = bf_type;
    model.bf_c = bf_c;
    model.poly = usePolyPart;

    %calculate and transform distances between all the points in the training data
    dist = zeros(n, n);
    switch upper(model.bf_type)
        case 'BH'
            if verbose
                fprintf('Building RBF (biharmonic) model...\n');
            end
            for i = 1 : n
                %for j = i : n
                %    dist(i, j) = norm(Xtr(i,:) - Xtr(j,:));
                %end
                dist(i, i:n) = sqrt(sum((repmat(Xtr(i,:),n-i+1,1) - Xtr(i:n,:)).^2,2));
                dist(i+1:n, i) = dist(i, i+1:n);
            end
        case 'IMQ'
            if verbose
                fprintf('Building RBF (inverse multiquadric) model...\n');
            end
            for i = 1 : n
                %for j = i : n
                %    dist(i, j) = 1 / sqrt(sum((Xtr(i,:) - Xtr(j,:)).^2) + bf_c^2);
                %end
                dist(i, i:n) = 1 ./ sqrt(sum((repmat(Xtr(i,:),n-i+1,1) - Xtr(i:n,:)).^2,2) + bf_c^2);
                dist(i+1:n, i) = dist(i, i+1:n);
            end
        case 'TPS'
            if verbose
                fprintf('Building RBF (thin plate spline) model...\n');
            end
            for i = 1 : n
                %for j = i : n
                %    dist(i, j) = sum((Xtr(i,:) - Xtr(j,:)).^2);
                %    dist(i, j) = (dist(i, j) + bf_c^2) * log(sqrt(dist(i, j) + bf_c^2));
                %end
                dist(i, i:n) = sum((repmat(Xtr(i,:),n-i+1,1) - Xtr(i:n,:)).^2,2);
                dist(i, i:n) = (dist(i, i:n) + bf_c^2) .* log(sqrt(dist(i, i:n) + bf_c^2));
                dist(i+1:n, i) = dist(i, i+1:n);
            end
        case 'G'
            if verbose
                fprintf('Building RBF (Gaussian) model...\n');
            end
            for i = 1 : n
                %for j = i : n
                %    dist(i, j) = exp(-sum((Xtr(i,:) - Xtr(j,:)).^2) / (2*bf_c^2));
                %end
                dist(i, i:n) = exp(-sum((repmat(Xtr(i,:),n-i+1,1) - Xtr(i:n,:)).^2,2) / (2*bf_c^2));
                dist(i+1:n, i) = dist(i, i+1:n);
            end
        otherwise %MQ
            if verbose
                fprintf('Building RBF (multiquadric) model...\n');
            end
            for i = 1 : n
                %for j = i : n
                %    dist(i, j) = sqrt(sum((Xtr(i,:) - Xtr(j,:)).^2) + bf_c^2);
                %end
                dist(i, i:n) = sqrt(sum((repmat(Xtr(i,:),n-i+1,1) - Xtr(i:n,:)).^2,2) + bf_c^2);
                dist(i+1:n, i) = dist(i, i+1:n);
            end
    end

    %calculate coefs
    if model.poly == 0
        model.coefs = dist \ (Ytr - model.meanY);
    else
        A = [dist, ones(n,1), Xtr; [ones(n,1), Xtr]', zeros(d+1,d+1)];
        model.coefs  = A \ [Ytr; zeros(d+1,1)];
    end

    time = toc;

    if verbose
        fprintf('Execution time: %0.2f seconds\n', time);
    end

end
return

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function result = circle_cum ( func, xc, yc, radius, norder )

%*****************************************************************************80
%
%% CIRCLE_CUM approximates an integral on the circumference of a circle in 2D.
%
%  Integration region:
%
%    Points (X,Y) such that:
%
%      ( X - XC )^2 + ( Y - YC )^2 = RADIUS^2.
%
%  Discussion:
%
%    An NORDER point, (NORDER-1)-th degree formula is used,
%    Stroud number U2:M-1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2004
%
%  Reference:
%
%    Arthur H Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971.
%
%  Parameters:
%
%    Input, function VALUE = FUNC ( X, Y ), the name of the user supplied
%    function of two variables which is to be integrated.  The name should
%    be passed as a quoted string.
%
%    Input, real XC, YC, the coordinates of the center of
%    the circle.
%
%    Input, real RADIUS, the radius of the circle.
%
%    Input, integer NORDER, the number of points to use.
%
%    Output, real RESULT, the approximate integral of the function.
%
  quad = 0.0E+00;

  for i = 1 : norder
    angle = ( 2 * i ) * pi / norder;
    x = xc + radius * cos ( angle );
    y = yc + radius * sin ( angle );
    quad = quad + feval ( func, x, y );
  end

  quad = quad / norder;

  volume = pi * radius * radius;
  result = quad * volume;

  return
end

repo: jburkardt-m


text:
function [energyV,entropyV,sumAvgV,corrV,invDiffMomV,contrastV,...
    clustShadeV,clustProminV] = haralickTextureAllDir(scanArray3M, nL, offsetsM, flagv, hWait)
% [energyV,entropyV,sumAvgV,corrV,invDiffMomV,contrastV,...
%     clustShadeV,clustProminV] = haralickTextureAllDir(scanArray3M, nL, offsetsM, flagv, hWait)
%
% Haralick texture with cooccurrence matrix combined from all directions.
%
% APA, 05/12/2016

% Generate flags
if ~exist('flagv','var')
    flagv = ones(1,8);
elseif exist('flagv','var') && isempty(flagv)
    flagv = ones(1,8);
end

% % Flag to draw waitbar
% waitbarFlag = 0;
% if exist('hWait','var') && ishandle(hWait)
%     waitbarFlag = 1;
% end

% % Grid resolution
% deltaX = deltaXYZv(1);
% deltaY = deltaXYZv(2);
% deltaZ = deltaXYZv(3);
%
% % Get Block size to process
% slcWindow = floor(2*patchSizeV(3)/deltaZ);
% rowWindow = floor(2*patchSizeV(1)/deltaY);
% colWindow = floor(2*patchSizeV(2)/deltaX);
%
% % Make sure that the window is of odd size
% if mod(slcWindow,2) == 0
%     slcWindow = slcWindow + 1;
% end
% if mod(rowWindow,2) == 0
%     rowWindow = rowWindow + 1;
% end
% if mod(colWindow,2) == 0
%     colWindow = colWindow + 1;
% end

slcWindow = size(scanArray3M,3);
rowWindow = size(scanArray3M,1);
colWindow = size(scanArray3M,2);

% Build distance matrices
numColsPad = 1;
numRowsPad = 1;
numSlcsPad = 1;

% Get number of voxels per slice
[numRows, numCols, numSlices] = size(scanArray3M);
numVoxels = numRows*numCols;

minQ = min(scanArray3M(:));
maxQ = max(scanArray3M(:));
dq = (maxQ-minQ)/nL/4;
levels = linspace(minQ,maxQ,nL);
levels = linspace(dq,maxQ-dq,nL);
q1 = imquantize_cerr(scanArray3M,nL);
q2 = imquantize(scanArray3M, levels);
%levels = multithresh(scanArray3M, nL);
%q3 = imquantize(scanArray3M, levels);

% Pad doseArray2 so that sliding window works also for the edge voxels
%scanArrayTmp3M = padarray(scanArray3M,[numRowsPad numColsPad
%numSlcsPad],NaN,'both'); % aa commented
if exist('padarray.m','file')
    q = padarray(q1,[numRowsPad numColsPad numSlcsPad],NaN,'both');
else
    q = padarray_oct(q1,[numRowsPad numColsPad numSlcsPad],NaN,'both');
end

% Quantize the image
%nL = 16;
% if any(isnan(q)) % aa commented
%     nL = nL-1; % aa commented
% end % aa commented
%q = imquantize_cerr(scanArrayTmp3M,nL); % aa commented
%clear scanArrayTmp3M; % aa commented
qmax = max(q(:));
nanFlag = 0;
if any(isnan(q(:)))
    nanFlag = 1;
    q(isnan(q))=qmax+1;
end
qs=sort(unique(q));
lq=length(qs);

for k = 1:length(qs)
    q(q==qs(k)) = k;
end

q = uint16(q); % q is the quantized image
%numSlcsWithPadding = size(q,3);

% % Create indices for 2D blocks
% [m,n,~] = size(q);
% m = uint32(m);
% n = uint32(n);
% colWindow = uint32(colWindow);
% rowWindow = uint32(rowWindow);
% slcWindow = uint32(slcWindow);
%
% % Index calculation adapted from
% % http://stackoverflow.com/questions/25449279/efficient-implementation-of-im2col-and-col2im
%
% %// Start indices for each block
% start_ind = reshape(bsxfun(@plus,[1:m-rowWindow+1]',[0:n-colWindow]*m),[],1); %//'
%
% %// Row indices
% lin_row = permute(bsxfun(@plus,start_ind,[0:rowWindow-1])',[1 3 2]);  %//'
%
% %// Get linear indices based on row and col indices and get desired output
% % imTmpM = A(reshape(bsxfun(@plus,lin_row,[0:ncols-1]*m),nrows*ncols,[]));
% indM = reshape(bsxfun(@plus,lin_row,(0:colWindow-1)*m),rowWindow*colWindow,[]);

% Create indices list
%indCooccurM = reshape(1:nL*nL,nL,nL);

% Number of offsets
numOffsets = size(offsetsM,1);

% Indices of last level to filter out
nanIndV = false([lq*lq,1]);
nanIndV([lq:lq:lq*lq-lq, lq*lq-lq:lq*lq]) = true;

% Build levels vector for mu, sig
levRowV = repmat(1:lq,[1 lq]);
levColV = repmat(1:lq,[lq 1]);
levColV = levColV(:)';

% Build list of indices for px and contrast calculation
for n=0:lq-1
    % indices for p(x-y), contrast
    indCtrstV = false(lq*lq,1);
    indCtrst1V = 1:lq-n;
    indCtrst2V = 1+n:lq;
    indCtrstTmpV = indCtrst1V + (indCtrst2V-1)*lq;
    indCtrstTmpV = [indCtrstTmpV indCtrst2V + (indCtrst1V-1)*lq];
    indCtrstV(indCtrstTmpV) = 1;
    indCtrstV(nanIndV) = [];
    indCtrstC{n+1} = indCtrstV;

    % indices for p(x+y)
    indPxPlusYv = false(lq*lq,1);
    indPxPlusYv(levRowV + levColV == n+2) = 1;
    indPxPlusYv(nanIndV) = [];
    indPxPlusYc{n+1} = indPxPlusYv;

    % indices for px
    indPxV = false(lq*lq,1);
    indPxV(lq*n+1:lq*(n+1)) = true;
    indPxV(nanIndV) = [];
    indPxC{n+1} = indPxV;
end

% Filter NaN indices from row/col levels
levRowV(nanIndV) = [];
levColV(nanIndV) = [];

% Build linear indices column/row-wise for Symmetry
indRowV = zeros(1,lq*lq);
for i=1:lq
    indRowV((i-1)*lq+1:(i-1)*lq+lq) = i:lq:lq*lq;
end

% Initialize
energyV = [];
entropyV = [];
sumAvgV = [];
corrcorrV = [];
invDiffMomV = [];
contrastV = [];
clustShadeV = [];
clustProminV = [];

tic
% Initialize cooccurrence matrix (vectorized for speed)
cooccurV = zeros(lq*lq,1,'single');
for off = 1:numOffsets

    offset = offsetsM(off,:);
    slc1M = q(numRowsPad+(1:numRows),numColsPad+(1:numCols),...
        numSlcsPad+(1:numSlices));
    slc2M = circshift(q,offset);
    slc2M = slc2M(numRowsPad+(1:numRows),numColsPad+(1:numCols),numSlcsPad+(1:numSlices))...
        + (slc1M-1)*lq;
    cooccurV = cooccurV + accumarray(slc2M(:),1, [lq*lq,1]); % patch-wise cooccurance

end

cooccurV = cooccurV + cooccurV(indRowV,:); % for symmetry
cooccurV(nanIndV,:) = [];
cooccurV = cooccurV./sum(cooccurV);

% Calculate scalar texture for this offset
% Angular Second Moment (Energy)
if flagv(1)
    energyV = sum(cooccurV.^2);
end
% Entropy
if flagv(2)
    entropyV = -sum(cooccurV.*log2(cooccurV+1e-10));
end
% Contrast, inverse Difference Moment
contrastV = 0;
invDiffMomV = 0;
for n=0:nL-1
    % px
    px(n+1) = sum(cooccurV(indPxC{n+1},:),1);
    % p(x+y)
    pXplusY = cooccurV(indCtrstC{n+1},:);
    % p(x-y)
    pXminusY = cooccurV(indCtrstC{n+1},:);
    % contrast
    if flagv(6)
        contrastV = contrastV + ...
            sum(n^2*cooccurV(indCtrstC{n+1},:));
    end
    % inv diff moment
    if flagv(5)
        invDiffMomV = invDiffMomV + ...
            sum((1/(1+n^2))*cooccurV(indCtrstC{n+1},:));
    end
end
% weighted pixel average (mu), weighted pixel variance (sig)
mu = (1:nL) * px';
sig = ((1:nL)-mu).^2 * px';
%px = sum(cooccurV(indCooccurM),2);
muX = mean(px);
sigX = std(px);
%         for colNum = 1:numCalcVoxs
clstrV = (levRowV + levColV - 2*mu);
% Cluster Shade
if flagv(7)
    clustShadeV = clstrV.*clstrV.*clstrV * cooccurV(:,1);
end
% Cluster Prominence
%             clustProminTmpV(colNum) = (levRowV + levColV - 2*mu(colNum)).^4 ...
%                 * cooccurPatchM(:,colNum);
if flagv(8)
    clustProminV = clstrV.*clstrV.*clstrV.*clstrV ...
        * cooccurV(:,1);
end
% Sum Avg
if flagv(3)
    sumAvgV = (levRowV + levColV) * cooccurV(:,1);
end
% Correlation
if flagv(4)
    corrV = (levRowV .* levColV * cooccurV(:,1) - ...
        muX*muX) / sigX^2; %var(1); % check this
end

%     % Average texture from all directions
%     energyV = energyV / numOffsets;
%     entropyV = entropyV / numOffsets;
%     contrastV = contrastV / numOffsets;
%     invDiffMomV = invDiffMomV / numOffsets;
%     clustShadeV = clustShadeV / numOffsets;
%     clustProminV = clustProminV / numOffsets;

%     if flagv(1)
%         energyV = energyV / numOffsets;
%     end
%     if flagv(2)
%         entropyV = entropyV / numOffsets;
%     end
%     if flagv(3)
%         sumAvgV = sumAvgV / numOffsets;
%     end
%     if flagv(4)
%         corrV = corrV / numOffsets;
%     end
%
%     if flagv(6)
%         contrastV = contrastV / numOffsets;
%     end
%     if flagv(5)
%         invDiffMomV = invDiffMomV / numOffsets;
%     end
%     if flagv(7)
%         clustShadeV = clustShadeV / numOffsets;
%     end
%     if flagv(8)
%         clustProminV = clustProminV / numOffsets;
%     end

%     if waitbarFlag
%         set(hWait, 'Vertices', [[0 0 slcNum/numSlices slcNum/numSlices]' [0 1 1 0]']);
%         drawnow;
%     end

% end
toc

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR