text:
function im = visualize3Dpara(alpha, sincosTheta, tran, f, baseShape, edgeAdj, varargin)

% function im = visualize3Dpara(alpha, sincosTheta, tran, f, baseShape, varargin)
%
% Parameters:
%       alpha: nbasis x 1
%       sincosTheta: 6 x 1
%       tran: 2x1
%       f: 1
%       baseShape: 3 x np x nbasis
%       edgeAdj: nedge x 2
%
% Addition parameters:
%    h, w, lineWidth, addNode, circSize
%
% Example:
%   stickStruct = getStickFigure('class', 'chair');
%   im = visualize3Dpara(alpha, sincosTheta, tran, f, stickStruct.baseShape{1}, stickStruct.edgeAdj{1});

para.h = 320;
para.w = 240;
para.lineWidth = 6;
para.addNode = true;
para.circSize = 8;
para = propval(varargin, para);

theta = sctheta2theta(sincosTheta);
x = alpha2x_proj(tran,alpha,theta,f,baseShape, 'w', para.w, 'h',para.h);
im = renderImage(x,edgeAdj,'addNoise',false,'noiseLevel',0,'h',para.h,'w',para.w,'lineWidth',6,'parallel',false,...
    'addNode',para.addNode, 'circSize', para.circSize);

end

author: jiajunwu


text:
function h = orientedGaussianKernel(size, sigma, theta)
%ORIENTEDGAUSSIANKERNEL Oriented Gaussian kernel for directional filtering
%
%   H = orientedGaussianKernel(SIZE, SIGMA, THETA)
%   THETA is in degrees.
%
%   Example
%     orientedGaussianKernel
%
%   See also
%     imDirectionalFilter, orientedLaplacianKernel, imfilter

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2012-12-04,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2012 INRA - Cepia Software Platform.

% ensure size is a vector
if length(size) == 1
    size = [size size];
end

% ensure sigma is a vector
if length(sigma) == 1
    sigma = [sigma sigma];
end

% compute reference for each axis
nx = round((size(1) - 1) / 2);
ny = round((size(2) - 1) / 2);

% create 3D grid
[x, y] = meshgrid(-nx:nx, -ny:ny);

% precompute angles
cot = cosd(-theta);
sit = sind(-theta);

% Apply rotation by -theta
x2 = x * cot - y * sit;
y2 = x * sit + y * cot;

% compute gaussian function for each point of the grid
sx2 = 2 * sigma(1) ^ 2;
sy2 = 2 * sigma(2) ^ 2;
h = exp(-(x2.^2 / sx2)) .* exp(-(y2.^2 / sy2));

% normalize intensities
h = h / (sum(h(:)));

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [ y2, m2, d2, f2 ] = ymdf_inc_hebrew ( y1, m1, d1, f1, days )

%*****************************************************************************80
%
%% YMDF_INC_HEBREW increments a Hebrew YMDF date by DAYS days.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 April 2103
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1,  M1, D1, real F1,
%    the YMDF date.
%
%    Input, real DAYS, the number of days to advance the date.
%
%    Output, integer Y2, M2, D2, real F2,
%    the incremented YMDF date.
%

%
%  Copy the parameters.
%
  y2 = y1;
  m2 = m1;
  d2 = d1;
  f2 = f1 + days;
%
%  Check the parameters.
%
  [ y2, m2, d2, f2, ierror ] = ymdf_check_hebrew ( y2, m2, d2, f2 );

  return
end

author: johannesgerer

repo: jburkardt-m


function ground_truth = load_groundtruth(base_path, video)

%see if there's a suffix, specifying one of multiple targets, for
%example the dot and number in 'Jogging.1' or 'Jogging.2'.
if numel(video) >= 2 && video(end-1) == '.' && ~isnan(str2double(video(end))),
    suffix = video(end-1:end);  %remember the suffix
    video = video(1:end-2);  %remove it from the video name
else
    suffix = '';
end

%full path to the video's files
if base_path(end) ~= '/' && base_path(end) ~= '\',
    base_path(end+1) = '/';
end
video_path = [base_path video '/'];

%try to load ground truth from text file (Benchmark's format)
filename = [video_path 'groundtruth_rect' suffix '.txt'];
f = fopen(filename);
assert(f ~= -1, ['No initial position or ground truth to load ("' filename '").'])

%the format is [x, y, width, height]
try
    ground_truth = textscan(f, '%f,%f,%f,%f', 'ReturnOnError',false);
catch  %#ok, try different format (no commas)
    frewind(f);
    ground_truth = textscan(f, '%f %f %f %f');
end
ground_truth = cat(2, ground_truth{:});
fclose(f);

end


text:
% This is a small script to compute the Bethe approximation of a small
% factor graph. The factor graph is generated by taking the Bayes net of
% Lauritzen and Spiegelhalter (1988) and setting the following random
% variables as evidence
%
%   visit to Asia  = false
%   positive X-ray = true
%   smoking        = true
%
% For this particular set of evidence, the Bethe approximation is exact
% because the corresponding factor graph does not have any cycles.
%
% For more information, see: Lauritzen and Speigelhalter (1988). "Local
% computations with probabilities on graphical structures and their
% application to expert systems." Journal of the Royal Statistical
% Society, Series B, Vol. 50, No. 2, pp. 157-194.
%
% Copyright (C) 2007 Peter Carbonetto. All Rights Reserved.
% This code is published under the Eclipse Public License.
%
% Author: Peter Carbonetto
%         Dept. of Computer Science
%         University of British Columbia
%         May 19, 2007

verbose = true;

% The remaining random variables can either be true or false. They are:
%   1. tuberculosis
%   2. lung cancer
%   3. bronchitis
%   4. tuberculosis or lung cancer
%   5. dysponoea
nv = 5;             % The number of random variables.
nf = 6;             % The number of factors.
K  = 2*ones(nv,1);  % The number of possible discrete assignments to the
                    % random variable at each site.
C = cell(nf,1);     % The factor neighbourhoods.
f = cell(nf,1);     % The factors.

% Set up the factors.
% ------------------
% This is p(tuberculosis | visit to Asia = false).
C{1} = 1;
f{1} = [ 0.01 0.99 ]';

% This is p(lung cancer | smoking = true).
C{2} = 2;
f{2} = [ 0.1 0.9 ]';

% This is p(bronchitis | smoking = true).
C{3} = 3;
f{3} = [ 0.6 0.4 ]';

% This is p(tuberculosis or lung cancer | tuberculosis, lung cancer).
C{4}     = [4 1 2];
p        = zeros(2,2,2);
p(:,1,1) = [ 0.99 0.01 ]';
p(:,1,2) = [ 0.99 0.01 ]';
p(:,2,1) = [ 0.99 0.01 ]';
p(:,2,2) = [ 0.99 0.01 ]';
f{4}     = p;

% This is p(positive X-ray = true | tuberculosis or lung cancer).
C{5} = 4;
f{5} = [ 0.98 0.5 ]';

% This is p(dyspnoea | tuberculosis or lung cancer, bronchitis).
C{6}     = [5 4 3];
p        = zeros(2,2,2);
p(:,1,1) = [ 0.9 0.1 ]';
p(:,1,2) = [ 0.7 0.3 ]';
p(:,2,1) = [ 0.8 0.2 ]';
p(:,2,2) = [ 0.1 0.9 ]';
f{6}     = p;
clear p

% Set up the junction graph.
% -------------------------
% These are the large regions of the junction graph.
Rv{1} = [1 2 4];
Rf{1} = [1 2 4];
NR{1} = [1];

Rv{2} = [3 4 5];
Rf{2} = [3 5 6];
NR{2} = [1];

% These are the small regions of the junction graph.
Sv{1} = [4];
Sf{1} = [];
NS{1} = [1 2];

nr = length(Rv);  % The number of large regions.
ns = length(Sv);  % The number of small regions.

% Infer the region marginals.
% --------------------------
% Compute the junction graph approximation to the variational free
% energy. Since the junction graph is a tree in this case (trivially,
% because there are only two large regions), the junction graph
% approximation is exact.
qR = bopt(K,C,f,Rv,Rf,Sv,Sf,NS,verbose);

% Output the marginal probabilities.
% ---------------------------------
p = ndsum(qR{1},[2 3]);
fprintf('Pr(tuberculosis | evidence) = %0.2f \n', p(1));

p = ndsum(qR{1},[1 3]);
fprintf('Pr(lung cancer | evidence)  = %0.2f \n', p(1));

p = ndsum(qR{2},[2 3]);
fprintf('Pr(bronchitis | evidence)   = %0.2f \n', p(1));

p = ndsum(qR{2},[1 2]);
fprintf('Pr(dysponoea | evidence)    = %0.2f \n', p(1));

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function f=comp_inonsepdgtreal_quinqux(coef,g,a,M,do_timeinv)
%COMP_INONSEPDGTREAL_QUINQUX  Compute Inverse discrete Gabor transform
%   Usage:  f=inonsepdgt(c,g,a,M);
%
%   Input parameters:
%         c     : Array of coefficients.
%         g     : Window function.
%         a     : Length of time shift.
%         M     : Number of channels
%         do_timeinv : Do a time invariant phase ?
%   Output parameters:
%         f     : Signal.
%
%
%   This is a computational subroutine, do not call it directly.

%   AUTHOR : Nicki Holighaus and Peter L. Søndergaard
%   TESTING: TEST_NONSEPDGT
%   REFERENCE: OK

% Check input paramameters.

M2=size(coef,1);
N=size(coef,2);
W=size(coef,3);
L=N*a;

coef2=zeros(M,N,W,assert_classname(coef,g));

coef2(1:M2,:,:)=coef;
if rem(M,2)==0
    coef2(M2+1:M,1:2:N-1,:)=conj(coef(M2-1:-1:2,1:2:N-1,:));
    coef2(M2:M,2:2:N  ,:)  =conj(coef(M2-1:-1:1,2:2:N,:));
else
    coef2(M2+1:M,1:2:N-1,:)=conj(coef(M2:-1:2,1:2:N-1,:));
    coef2(M2+1:M,2:2:N  ,:)=conj(coef(M2-1:-1:1,2:2:N,:));
end;

coef=coef2;

lt=[1 2];
mwin=comp_nonsepwin2multi(g,a,M,lt,L);

% phase factor correction (backwards), for more information see
% analysis routine

E = exp(2*pi*i*a*kron(0:N/2-1,ones(1,2)).*...
        rem(kron(ones(1,N/2), 0:2-1),2)/M);

coef = bsxfun(@times,coef,E);

% simple algorithm: split into sublattices and add the result from eacg
% sublattice.
f=zeros(L,W,assert_classname(coef,g));
for ii=0:2-1
    % Extract sublattice
    sub=coef(:,ii+1:2:end,:);
    f=f+comp_idgt(sub,mwin(:,ii+1),2*a,[0 1],0,0);
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% Test file for @chebfun/times.m.

function pass = test_times(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebfunpref();
end

% Generate a few random points to use as test values.
seedRNG(6178);
x = 2 * rand(100, 1) - 1;

% Random numbers to use as arbitrary multiplicative constants.
alpha = -0.194758928283640 + 0.075474485412665i;

%% SCALAR-VALUED

% Check operation for empty inputs.
f = chebfun(@sin, [-1 1], pref);
pass(1) = isempty(f.*[]);
pass(2) = isempty(f.*chebfun());

% Turn on splitting, since we'll need it for the rest of the tests.
pref.splitting = 1;

% Test multiplication by scalars.
f1_op = @(x) sin(x).*abs(x - 0.1);
f1 = chebfun(f1_op, pref);
pass(3:4) = test_mult_function_by_scalar(f1, f1_op, alpha, x);

% Test multiplication of two chebfun objects.
g1_op = @(x) cos(x).*sign(x + 0.2);
g1 = chebfun(g1_op, pref);
pass(5:6) = test_mult_function_by_function(f1, f1_op, g1, g1_op, x);

%% ARRAY-VALUED

% Test operation for array-valued chebfuns.
f2_op = @(x) [sin(x).*abs(x - 0.1)  exp(x)];
f2 = chebfun(f2_op, pref);
pass(7:8) = test_mult_function_by_scalar(f2, f2_op, alpha, x);

g2_op = @(x) [cos(x).*sign(x + 0.2) tan(x)];
g2 = chebfun(g2_op, pref);
pass(9:10) = test_mult_function_by_function(f2, f2_op, g2, g2_op, x);

% Test operation for transposed chebfuns.
pass(11:12) = test_mult_function_by_scalar(f1.', @(x) f1_op(x).', alpha, x);
pass(13:14) = test_mult_function_by_function(f1.', @(x) f1_op(x).', ...
    g1.', @(x) g1_op(x).', x);

%% QUASIMATRICES
f2q = quasimatrix(f2_op, pref);
pass(15:16) = test_mult_function_by_scalar(f2q, f2_op, alpha, x);

% Quasi * array-cheb
pass(17:18) = test_mult_function_by_function(f2q, f2_op, g2, g2_op, x);

% Quasi * quasi
g2q = quasimatrix(g2_op, pref);
pass(19:20) = test_mult_function_by_function(f2q, f2_op, g2q, g2_op, x);

%% Check error conditions.

try
    h = f1.*'X';
    pass(21) = false;
catch ME
    pass(21) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN:times:unknown');
end

try
    h = f1.*f1.';
    pass(22) = false;
catch ME
    pass(22) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN:times:matdim');
end

try
    h = f1.*g2q.';
    pass(23) = false;
catch ME
    pass(23) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN:times:matdim');
end

%% Test on singular function:

dom = [-2 7];

% Generate a few random points to use as test values.
seedRNG(6178);
x = diff(dom) * rand(100, 1) + dom(1);

%% Case of a scalar and a function:
c = 3;
pow = -0.5;
op = @(x) (x - dom(2)).^pow.*sin(x);
op_exact = @(x) c*(x - dom(2)).^pow.*sin(x);
f = chebfun(op, dom, 'exps', [0 pow], 'splitting', 'on');
g = c.*f;
g_exact = chebfun(op_exact, dom, 'exps', [0 pow], 'splitting', 'on');

err = norm(feval(g, x) - feval(g_exact, x), inf);
pass(24) = ( err < 50*eps*norm(feval(g_exact, x), inf) );

%% Case of two functions: piecewise smooth chebfun - splitting on.
pow1 = -0.3;
pow2 = -0.5;
op1 = @(x) (x - dom(2)).^pow1.*sin(100*x);
op2 = @(x) (x - dom(2)).^pow2.*cos(300*x);
op_exact = @(x) (x - dom(2)).^(pow1+pow2).*sin(100*x).*cos(300*x);
f = chebfun(op1, dom, 'exps', [0 pow1], 'splitting', 'on');
g = chebfun(op2, dom, 'exps', [0 pow2] , 'splitting', 'on');
h = f.*g;
h_exact = chebfun(op_exact, dom, 'exps', [0 pow1+pow2], 'splitting', 'on');

err = norm(feval(h, x) - feval(h_exact, x), inf);
pass(25) = ( err < 1e4*eps*...
    norm(feval(h_exact, x), inf) );

%% Tests for function defined on unbounded domain:

% Functions on [-inf inf]:

% Set the domain:
dom = [-Inf Inf];
domCheck = [-1e2 1e2];

% Generate a few random points to use as test values:
x = diff(domCheck) * rand(100, 1) + domCheck(1);

opf = @(x) x.^2.*exp(-x.^2);
opg = @(x) (1-exp(-x.^2))./x;
oph = @(x) x.*exp(-x.^2).*(1-exp(-x.^2));
f = chebfun(opf, dom);
g = chebfun(opg, dom);
h = f.*g;
hVals = feval(h, x);
hExact = oph(x);
err = hVals - hExact;
pass(26) = norm(err, inf) < 10*eps*get(f,'vscale');

%% Test multiplication between a CHEBFUN and a TRIGFUN.

dom = [0 pi 2*pi];

% 1. One column case.
f = chebfun(@(x) x.^2, dom, pref);
g = chebfun(@(x) cos(x), [dom(1) dom(end)], 'periodic');
h1 = f.*g;
% We want the result to use the same tech as the one used by f.
pass(27) = strcmpi(func2str(get(h1.funs{1}.onefun, 'tech')), ...
                   func2str(get(f.funs{1}.onefun, 'tech')));
h2 = chebfun(@(x) (x.^2).*cos(x), dom, pref);
pass(28) = norm(h1-h2, inf) < 1e1*eps*get(h2,'vscale');

% 2. Quasimatrix case.
f = chebfun(@(x) [cos(x), sin(x)], [dom(1) dom(end)], 'periodic');
g = chebfun(@(x) [x, x.^3], dom, pref);
h1 = f.*g;
% We want the result to use the same tech as the one used by g.
pass(29) = strcmpi(func2str(get(h1(:,1).funs{1}.onefun, 'tech')), ...
                   func2str(get(g(:,1).funs{1}.onefun, 'tech')));
pass(30) = strcmpi(func2str(get(h1(:,2).funs{1}.onefun, 'tech')), ...
                   func2str(get(g(:,2).funs{1}.onefun, 'tech')));
h2 = chebfun(@(x) [x.*cos(x), x.^3.*sin(x)], dom, pref);
pass(31) = norm(h1-h2, inf) < 1e2*eps*get(h2,'vscale');

end

%% The tests

% Test the multiplication of a chebfun F, specified by F_OP, by a scalar ALPHA
% using a grid of points X in [-1  1] for testing samples.
function result = test_mult_function_by_scalar(f, f_op, alpha, x)
    g1 = f .* alpha;
    g2 = alpha .* f;
    result(1) = isequal(g1, g2);
    g_exact = @(x) f_op(x) .* alpha;
    err = feval(g1, x) - g_exact(x);
    result(2) = norm(err(:), inf) < 1e2*max(vscale(g1)*eps);

end

% Test the multiplication of two chebfuns F and G, specified by F_OP and G_OP,
% using a grid of points X in [-1  1] for testing samples.
function result = test_mult_function_by_function(f, f_op, g, g_op, x)
    h1 = f .* g;
    h2 = g .* f;
    result(1) = norm(h1 - h2) < 10*max(vscale(h1)*eps);
    h_exact = @(x) f_op(x) .* g_op(x);
    err = feval(h1, x) - h_exact(x);
    result(2) = norm(err(:), inf) < 1e2*max(vscale(h1)*eps);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


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

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

  if ( t <= 10.0 )
    jac(1,1) = 0.0;
  else
    jac(1,1) = -2.0;
  end

  return
end


%% Trend-Cycle tutorial: Predicting recessions
% Authors:   Filippo Ferroni and  Fabio Canova
% Date:     27/05/2020, revised  15/12/2020

close all; clc; clear all;

addpath ../../cmintools/
addpath ../../bvartools/

% predicting  recession and the start of  a recession in the  euro  area.
% use recession indicator constructed with example_2_dating.m; regressors
% labor  productivity, commodity prices, long term  rates  or spread.

% Euro area AWM DATABASE: Quarterly
    [a,b,~] = xlsread('awm19up18.xlsx');
    % names of variables
    varnames = b(1,2:end);

    % time convention: Q1 = .00 and Q4 =0.75
    time = 1970 : .25 : 2017.75;
    time_start = find(time==1970.50);
    time_start1 = find(time==1999.50);
    time_end   = find(time==2017.75);
    time_break   = find(time==2007.75);

    % The  CREDIT DATA: Quartely
    [e,f,~]=xlsread('ECB_Credit_gaps.xlsx');
    varnames2= f(1,2:end);

    % real GDP, real consumption, real  investment, GDP defl, HICP,
    % short and long term interest rate, commodity  prices,
    % labor  productivity, total  employment, credit to NFC to  GDP
    yy = [a(:,1) a(:,2) a(:,4) a(:,7) a(:,19) a(:,33) a(:,34) ...
        a(:,35) a(:,42) a(:,29) e(time_start-2:time_end,25)];
    % names: YER PCR ITR YED HICP STN LTN COMPR LPROD LNN

    % data transformations:
    % 1. the log of output, consumption,investment
    ddata(:,1:3) = log(yy(time_start+1 : time_end,1:3));
    % 2. the log/log difference of GDP deflator index and CPI
    %ddata(:,4:5) = log(yy(time_start+1  : time_end,4:5));
    ddata(:,4:5) = diff(log(yy(time_start : time_end,4:5)))*400;
    % 3. the level of short and long term interest rate
    ddata(:,6:7) = yy(time_start +1 : time_end,6:7);
    % 4. the log/log difference of commodity prices
    ddata(:,8) = log(yy(time_start+1  : time_end,8));
    %ddata(:,8) = diff(log(yy(time_start : time_end,8)))*400;
    % 5. the taking the log of labor  productivity and  employment
    ddata(:,9:10) = log(yy(time_start+1 : time_end,9:10));
    %ddata(:,9) = diff(log(yy(time_start+1 : time_end,9)));
    %ddata(:,10) = log(yy(time_start+1 : time_end,10));

    % 6. spread
    ddata(:,11)= ddata(:,7)-ddata(:,6);
    %  data for  credit  to  GDP starts  only  at  time_start1 (1999:25)
    ddata(:,12) = yy(time_start+1:time_end,11);
    endd=length(ddata);

 % cc=ones(time_end-time_start,1);  % constant
 % pick log labor productivity, log commodity prices, long term rate
 % z=ddata(:,[9 8 7]);
 % pick log labor productivity, log commodity prices, spread
 z=ddata(:, [9 8 11]);

load  Eurorec
% recind has  the  recession indicator  created  with  dating_exa.m
x=recind;

% use contemporanous  values
% [estimator,cov_Hessian,ME1,ME2,ME_std] = probit(x,z);

% use  one  period lagged  values
zz=zeros(size(z,1), size(z,2));
zz(2:length(z),:)=z(1:length(z)-1,:);
timeplot = time(time_start:time_end-1);

% Estimate Probit model and the marginal effects
[estimator,cov_Hessian,ME1,ME2,ME1_std] = probit(x,zz);

% prediction
pz=[ones(length(zz),1) zz];
px=normpdf(pz*estimator);
half=0.5*ones(length(z),1);

plot(timeplot,px,'r-.','Linewidth',2); hold on;
plot(timeplot,half,'b:','Linewidth',2); hold  on;
plot(timeplot,x,'k-','Linewidth',2); hold  off; axis  tight;
legend('predition','halfline', 'actual')
pause

close  all;

% do  estimation recursively: predicting  the  beginning  of  a  recession
% which is  the  same  as  predicting  a  peak
% recession dates
%rec1b = find(time==1974.00);
%rec1e = find(time==1975.25);
%[estimator1,cov_Hessian1,ME11,ME21,ME11_std] = probit(x(1:rec1b),zz(1:rec1b,:));

rec2b = find(time==1980.00);
rec2e = find(time==1984.25);
disp('Predicting 1980 recession')
[estimator2,cov_Hessian2,ME12,ME22,ME12_std] = probit(x(1:rec2b),zz(1:rec2b,:));
pz2=[ones(rec2b,1) zz(1:rec2b,:)];
px2=normpdf(pz2*estimator2);
half2=0.5*ones(rec2b,1);

plot(timeplot(1:rec2b),x(1:rec2b),'k-','Linewidth',2); hold  on;
plot(timeplot(1:rec2b),half2(1:rec2b),'b:','Linewidth',2); hold  on;
plot(timeplot(1:rec2b),px2(1:rec2b),'r-.','Linewidth',2); hold  off;axis tight;
legend('actual', 'halfline', 'predition')
pause

close all;

rec3b = find(time==1992.00);
rec3e = find(time==1993.75);
disp('Predicting  1992 recession')
[estimator3,cov_Hessian3,ME13,ME23,ME13_std] = probit(x(1:rec3b),zz(1:rec3b,:));
pz3=[ones(rec3b,1) zz(1:rec3b,:)];
px3=normpdf(pz3*estimator3);
half3=0.5*ones(rec3b,1);

plot(timeplot(1:rec3b),x(1:rec3b),'k-','Linewidth',2); hold  on;
plot(timeplot(1:rec3b),half3(1:rec3b),'b:','Linewidth',2); hold  on;
plot(timeplot(1:rec3b),px3(1:rec3b),'r-.','Linewidth',2); hold  off; axis  tight;
legend('actual', 'halfline','predition')
pause
close all;

rec4b = find(time==2001.00);
rec4e = find(time==2002.75);
disp('Predicting  2001 recession')
[estimator4,cov_Hessian4,ME14,ME24,ME14_std] = probit(x(1:rec4b),zz(1:rec4b,:));
pz4=[ones(rec4b,1) zz(1:rec4b,:)];
px4=normpdf(pz4*estimator4);
half4=0.5*ones(rec4b,1);

plot(timeplot(1:rec4b),x(1:rec4b),'k-','Linewidth',2); hold  on;
plot(timeplot(1:rec4b),half4(1:rec4b),'b:','Linewidth',2); hold  on;
plot(timeplot(1:rec4b),px4(1:rec4b),'r-.','Linewidth',2); hold  off; axis  tight;
legend('actual', 'halfline','predition')
pause
close all;

rec5b = find(time==2008.00);
rec5e = find(time==2009.50);
disp('Predicting 2008 recession')
[estimator5,cov_Hessian5,ME15,ME25,ME15_std] = probit(x(1:rec5b),zz(1:rec5b,:));
pz5=[ones(rec5b,1) zz(1:rec5b,:)];
px5=normpdf(pz5*estimator5);
half5=0.5*ones(rec5b,1);

plot(timeplot(1:rec5b),x(1:rec5b),'k-','Linewidth',2); hold  on;
plot(timeplot(1:rec5b),half5(1:rec5b),'b:','Linewidth',2); hold  on;
plot(timeplot(1:rec5b),px5(1:rec5b),'r-.','Linewidth',2); hold  off; axis tight;
legend('actual', 'halfline', 'predition')
pause
close all;

rec6b = find(time==2011.25);
rec6e = find(time==2013.00);
disp('Predicting  2011 recession')
[estimator6,cov_Hessian6,ME16,ME26,ME16_std] = probit(x(1:rec6b),zz(1:rec6b,:));
pz6=[ones(rec6b,1) zz(1:rec6b,:)];
px6=normpdf(pz6*estimator6);
half6=0.5*ones(rec6b,1);

plot(timeplot(1:rec6b),x(1:rec6b),'k-','Linewidth',2); hold  on;
plot(timeplot(1:rec6b),half6(1:rec6b),'b:','Linewidth',2); hold  on;
plot(timeplot(1:rec6b),px6(1:rec6b),'r-.','Linewidth',2); hold  off; axis  tight;
legend('actual', 'halfline', 'predition')

return


text:
classdef SettingsHomogenizedVarComputerFromInterpolation < ....
        SettingsHomogenizedVarComputer

    properties (Access = protected)
        defaultParamsName = 'paramsHomogenizedVarComputerFromInterpolation.json'
    end

    properties (Access = public)
        interpolation
        typeOfMaterial
        constitutiveProperties
        nelem
        ptype
        interpolationSettings
    end

    methods (Access = public)

        function obj = SettingsHomogenizedVarComputerFromInterpolation(varargin)
            if nargin == 1
                obj.loadParams(varargin{1})
            end
            obj.init();
        end

    end

    methods (Access = private)

        function init(obj)
            obj.initConstitutiveProperties();
            obj.initInterpolationSettings();
        end

        function initConstitutiveProperties(obj)
            s = obj.constitutiveProperties;
            obj.constitutiveProperties = SettingsConstitutiveProperties(s);
        end

        function initInterpolationSettings(obj)
            s.interpolation           = obj.interpolation;
            s.dim                     = obj.dim;
            s.typeOfMaterial          = obj.typeOfMaterial;
            s.constitutiveProperties  = obj.constitutiveProperties;
            obj.interpolationSettings = SettingsInterpolation(s);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [Feven,Fodd] = filterGabor2d( r, sig, lam, theta, omega, show )
% Creates an even/odd pair of 2D Gabor filters.
%
% Creates a pair of Gabor filters (one odd one even) at the specified
% orientation. For Thomas' ECCV98 filters, use sig=sqrt(2), lam=4.  Note
% that originally this function computed a quadratic masked with a
% Gaussian, and not a sin/cos masked with a Gaussian. Requires Matlab's
% 'Signal Processing Toolbox'.
%
% USAGE
%  [Feven,Fodd] = filterGabor2d( r, sig, lam, theta, [omega], [show] )
%
% INPUTS
%  r       - final mask will be 2r+1 x 2r+1
%  sig     - standard deviation of Gaussian mask
%  lam     - elongation of Gaussian mask
%  theta   - orientation (in degrees)
%  omega   - [1] wavlength of underlying sine (sould be >=1)
%  show    - [0] figure to use for optional display
%
% OUTPUTS
%  Feven   - even symmetric filter (-cosine masked with Gaussian)
%  Fodd    - odd symmetric filter (-sine masked with Gaussian)
%
% EXAMPLE
%  [Feven,Fodd] = filterGabor2d(15,sqrt(2),4,45,1,1);
%  [Feven,Fodd] = filterGabor2d(25,4,2,0,2,3);
%
% See also FILTERGABOR1D, FILTERGAUSS
%
% Piotr's Image&Video Toolbox      Version 2.12
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<5 || isempty(omega) ); omega=1; end
if( nargin<6 || isempty(show) ); show=0; end

% create even/odd Gabor filters
[x,y]=meshgrid(-r:r,-r:r);
mask = exp(-(y.^2)/(sig^2)-(x.^2)/(lam^2*sig^2));
Feven = -cos(2*pi*y/4/omega) .* mask;
%Feven = (4*(y.^2)/(sig^4)-2/(sig^2)).*mask; % original function was y.^2
Fodd = imag(hilbert(Feven));

% orient appropriately
Feven = imrotate(Feven,theta,'bil','crop');
Fodd = imrotate(Fodd,theta,'bil','crop');

% Set mean to 0 (should already be 0)
Feven=Feven-mean(Feven(:));
Fodd=Fodd-mean(Fodd(:));

% set L1norm to 0
Feven=Feven/norm(Feven(:),1);
Fodd=Fodd/norm(Fodd(:),1);

% display
if( show )
  filterVisualize( Feven, show );
  filterVisualize( Fodd, show+1 );
end

author: zouchuhang

repo: LayoutNet


text:
function [version, released] = manopt_version()
% Returns the version of the Manopt package you are running, as a vector.
%
% function [version, released] = manopt_version()
%
% version(1) is the primary version number.
% released is the date this version was released, in the same format as the
% date() function in Matlab.

    version = [7, 1, 0];
    released = '30-Sep-2022';

end

author: NicolasBoumal

repo: manopt


text:
function  [im_out, par] = WSC_Sigma_1AG(par)
im_out    =   par.nim;
% parameters for noisy image
[h,  w, ch]      =  size(im_out);
par.h = h;
par.w = w;
par.ch = ch;
par = SearchNeighborIndex( par );
for ite  =  1 : par.outerIter
    %     % iterative regularization
    im_out = im_out+par.delta*(par.nim - im_out);
    % image to patches and estimate local noise variance
    Y = Image2Patch( im_out, par );
    % estimation of noise variance
    dif = mean( mean( mean( (par.nim-im_out).^2 ) ) );
    par.sigma = sqrt( abs( par.nSig^2 - dif ) );
    % estimation of noise variance
    if mod(ite-1, par.innerIter)==0
        par.nlsp = par.nlsp - par.nlspgap;
        % searching  non-local patches
        blk_arr = Block_Matching( Y, par );
    end
    % Weighted Sparse Coding
    Y_hat = zeros(par.ps2ch, par.maxrc, 'single');
    W_hat = zeros(par.ps2ch, par.maxrc, 'single');
    for i = 1:par.lenrc
        index = blk_arr(:, i);
        nlY = Y( : , index );
        DC = mean(nlY, 2);
        nDCnlY = bsxfun(@minus, nlY, DC);
        % update D and S
        [D, S, ~] = svd( full(nDCnlY), 'econ' );
        S = diag(S);
        % update weight for sparse coding
        Wsc = repmat(bsxfun( @rdivide, par.lambda*par.sigma^2, S + eps ), [1 size(nDCnlY, 2)]);
        % update C by soft thresholding
        B = D' * nDCnlY;
        C = sign(B) .* max( abs(B) - Wsc, 0 );
        % update Y
        nDCnlYhat = D * C;
        % add back DC components
        nlYhat = bsxfun(@plus, nDCnlYhat, DC);
        % aggregation
        Y_hat(:, index) = Y_hat(:, index) + nlYhat;
        W_hat(:, index) = W_hat(:, index) + ones(par.ps2ch, par.nlsp);
    end
    % Reconstruction
    im_out = PGs2Image(Y_hat, W_hat, par);
    % calculate the PSNR and SSIM
    PSNR =   csnr( im_out*255, par.I*255, 0, 0 );
    SSIM      =  cal_ssim( im_out*255, par.I*255, 0, 0 );
    fprintf('Iter %d : PSNR = %2.4f, SSIM = %2.4f\n', ite, PSNR, SSIM);
    par.PSNR(ite, par.image) = PSNR;
    par.SSIM(ite, par.image) = SSIM;
end
return;

author: csjunxu


text:
function coords = rmGetCoords(vw, model)
% Get the coordinates of the voxels associated with a retinotopy model.
%
% coords = rmGetCoords(view, [model=selected model]);
%
% This is almost the same as rmGet(model, 'coords'), except it allows
% the user to enter the view (instead of grabbing the current view), so
% it should be more robust.
%
% ras, 12/2006.
if notDefined('vw'),  vw = getCurView;  end
if notDefined('model'), model = viewGet(vw, 'rmCurModel'); end
if isnumeric(model)
    allModels = viewGet(vw, 'rmModel');
    model = allModels{model};
end

if isfield(model, 'roi_coordinates') && ~isempty(model.roi_coordinates)
    val  = model.roi_coordinates;
else
    if isfield(vw, 'coords') && ~isempty(vw.coords)
        coords = vw.coords;   % gray / volume
    else        % inplane / flat
        sz = viewGet(vw,'Size');
        [X, Y, Z] = meshgrid(1:sz(2), 1:sz(1), 1:sz(3));
        coords = [X(:) Y(:) Z(:)]';
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function progressbar(n,N,w)
%PROGRESSBAR Display a progress bar.
%
%    progressbar(n,N,w);
%
% displays the progress of n out of N.
% n should start at 1.
% w is the width of the bar (default w=20).
%
%   Copyright (c) Gabriel Peyré 2006

if nargin<3
    w = 20;
end

% progress char
cprog = '.';
cprog1 = '*';
% begining char
cbeg = '[';
% ending char
cend = ']';

p = min( floor(n/N*(w+1)), w);

global pprev;
if isempty(pprev)
    pprev = -1;
end

if not(p==pprev)
    ps = repmat(cprog, [1 w]);
    ps(1:p) = cprog1;
    ps = [cbeg ps cend];
    if n>1
        % clear previous string
        fprintf( repmat('\b', [1 length(ps)]) );
    end
    fprintf(ps);
end
pprev = p;
if n==N
    fprintf('\n');
end

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% child script of canlab_glm_publish
% runs robust_results_batch in publishable manner
% runs in directory containing robust* directories

%% prep
z = '________________________________________________________________________________________________________';
assignin('base','z',z);
if ~exist('thresh','var'), thresh = [.001 .005 .05]; assignin('base','thresh',thresh); end
if ~exist('size','var'), size = [5 1 1]; assignin('base','size',size); end

if exist('EXPT','var')
    for i = 1:length(EXPT.SNPM.P)
        robregdirs{i} = fullfile(pwd,sprintf('robust%04d',i));
        connames{i} = EXPT.SNPM.connames(i,:);
    end
else
    d = filenames('robust[0-9][0-9][0-9][0-9]','absolute');
    if isempty(d), d{1} = pwd; end
    for i = 1:numel(d)
        robregdirs{i} = d{i};
        load(fullfile(robregdirs{i},'SETUP.mat'));
        load(fullfile(fileparts(SETUP.files(1,:)),'SPM.mat'));
        c = regexprep(SETUP.files(1,:),'.*con_0*([0-9]*)\.nii','$1');
        connames{i} = SPM.xCon(str2num(c)).name;
    end
end

%% write script
scriptfile = 'robfit_results.m';
scriptfileabs = fullfile(pwd,scriptfile);
fid = fopen(scriptfile,'w');

for i = 1:length(robregdirs)
    fprintf(fid,'%%%% contrast(%d) %s\n',i,connames{i});
    fprintf(fid,'cd(''%s'')\n',robregdirs{i});
    fprintf(fid,'fprintf(''%%s\\n%%s\\n%%s\\n%%s\\n%%s\\n%%s\\n'',z,z,''%s'',''%s'',z,z)\n',connames{i},robregdirs{i});

    fprintf(fid,'try\n');
    maskimg = fullfile(robregdirs{i}, 'rob_tmap_0001.img'); % any image in space with non-zero vals for all vox would do
    %fprintf(fid,'robust_results_batch(''thresh'', thresh, ''size'', size, ''prune'', ''mask'', ''%s'');\n',maskimg);
    fprintf(fid,'robust_results_batch;\n');

    fprintf(fid,'catch exc\n');
    fprintf(fid,'disp(getReport(exc,''extended''))\n');
    fprintf(fid,'end\n');
end
fprintf(fid,'close all\n');

fclose(fid);

%% publish script
outputdir = fullfile(pwd, 'Robust_Regression_Results_html');
mkdir(outputdir);

p = struct('useNewFigure', false, 'maxHeight', 1500, 'maxWidth', 1200, ...
    'outputDir', outputdir, 'showCode', false);

fout = publish(scriptfile, p);
fprintf('Created robfit results directory:\n%s\nHTML report: %s\n', outputdir, fout);

%% clean up
close all
delete(scriptfileabs)

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% function [lik,likv] = hmm_cl(X,T,K,Mu,Cov,P,Pi);
%
% Calculate Likelihood for Hidden Markov Model
%
% X - N x p data matrix
% T - length of each sequence (N must evenly divide by T, default T=N)
% K - number of states
% Mu - mean vectors
% Cov - output covariance matrix (full, tied across states)
% P - state transition matrix
% Pi - priors
%
% lik - log likelihood of X
% likv - vector of log likelihoods of each sequence
%
% If 0 or 1 output arguments requested, lik is returned. If 2 output
% arguments requested, [lik likv] is returned.
%
% Machine Learning Toolbox
% Version 1.0  01-Apr-96
% Copyright (c) by Zoubin Ghahramani
% http://mlg.eng.cam.ac.uk/zoubin/software.html
%
% ------------------------------------------------------------------------------
% The MIT License (MIT)
%
% Copyright (c) 1996, Zoubin Ghahramani
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.
% ------------------------------------------------------------------------------

function [lik,likv] = ZG_hmm_cl(X,T,K,Mu,Cov,P,Pi)

p = length(X(1,:));
N = length(X(:,1));
tiny = exp(-700);

if (rem(N,T) ~= 0)
  error('Data matrix length must be multiple of sequence length T');
end;
N = N/T;

alpha = zeros(T,K);
B = zeros(T,K);      % P( output | s_i)

k1 = (2*pi)^(-p/2);

Scale = zeros(T,1);

likv = zeros(1,N);

for n = 1:N

  B = zeros(T,K);
  iCov = inv(Cov);
  k2 = k1/sqrt(det(Cov));
  for i = 1:T
    for l = 1:K
      d = Mu(l,:)-X((n-1)*T+i,:);
      B(i,l) = k2*exp(-0.5*d*iCov*d');
    end;
  end;

  scale = zeros(T,1);
  alpha(1,:) = Pi(:)'.*B(1,:);
  scale(1) = sum(alpha(1,:));
  alpha(1,:) = alpha(1,:)/(scale(1)+tiny);
  for i = 2:T
    alpha(i,:) = (alpha(i-1,:)*P).*B(i,:);
    scale(i) = sum(alpha(i,:));
    alpha(i,:) = alpha(i,:)/(scale(i)+tiny);
  end;

  likv(n) = sum(log(scale+(scale == 0)*tiny));
  Scale = Scale+log(scale+(scale == 0)*tiny);
end;

lik = sum(Scale);

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function build_computeMinMarginalsBinaryMex
% build_computeMinMarginalsBinaryMex builds package computeMinMarginalsBinaryMex
%
% Anton Osokin,  12.04.2015

srcFiles = { 'computeMinMarginalsBinaryMex.cpp' };
allFiles = '';
for iFile = 1 : length(srcFiles)
    allFiles = [allFiles, ' ', srcFiles{iFile}];
end

cmdLine = ['mex ', allFiles, ' -output computeMinMarginalsBinaryMex -largeArrayDims '];
eval(cmdLine);

author: aosokin


text:
function op = proj_nuclear( q, sz )

%PROJ_NUCLEAR  Projection onto the set of matrices with nuclear norm less than or equal to q.
%    OP = PROJ_NUCLEAR( q ) returns a function that implements the
%    indicator for matrices with nuclear norm less than q.
%    Q is optional; if omitted, Q=1 is assumed. But if Q is supplied,
%    it must be a positive real scalar.
%    This function is like proj_psdUTrace.m but does not assume
%    that inputs are square Hermitian positive semidefinite matrices.
%
%    OP = PROJ_NUCLEAR( ..., SZ )
%    where SZ = [n1,n2], will assume the input has been vectorized
%    and thus will reshape it into a matrix of size n1 x n2. TFOCS
%    can handle matrix variables so this form of PROJ_NUCLEAR is not
%    always necessary, but it can be easier to find bugs if using
%    only vector variables.
%
% This version uses a dense svd decomposition; future versions
% of TFOCS will take advantage of low-rank and/or sparse structure.
% Dual function: prox_spectral.m
% See also prox_spectral, proj_psdUTrace.m, proj_simplex.m, prox_nuclear.m

% June 26, 2012

if nargin == 0,
	q = 1;
elseif ~isnumeric( q ) || ~isreal( q ) || numel( q ) ~= 1 || q <= 0,
	error( 'Argument must be positive.' );
end
if nargin < 2, sz = []; end
q = proj_simplex( q );
op = @(varargin)proj_nuclear_q( q, sz, varargin{:} );

function [ v, X ] = proj_nuclear_q( eproj, sz, X, t )

VECTORIZE = false;
% Input must be a matrix

% In proj_psdUTrace.m, we allow vector inputs because
%   we can reshape them (since we have square matrices).
%   For nuclear norm, inputs can be rectangular matrices,
%   so we do not know how to correctly reshape a vector in
%   this case.

% if size(X,1) ~= size(X,2)
%     n = sqrt(length(X));
%     X = reshape(X, n, n );
%     VECTORIZE   = true;
% end

% Update, Jan 2014, allowing vector variables, as long
%   as the 'sz' variable has been given
if ~isempty(sz)
    if size(X,2) > 1
        error('proj_nuclear: when using an explicit size value "sz", variable should be a vector');
    end
    X = reshape(X, sz(1), sz(2) );
    VECTORIZE   = true;
elseif size(X,2) == 1
    warning('proj_nuclear: appears variable is a vector not a matrix. Are you sure you want the SVD of a vector?');
end

v = 0;
if nargin > 3 && t > 0,
    [U,D,V] = svd(X,'econ');
    [dum,D] = eproj(diag(D),t);
    tt = D > 0;
    X  = U(:,tt)*diag(D(tt))*V(:,tt)';
    if VECTORIZE, X = X(:); end
elseif any(svd(X)<0),
    v = Inf;
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [node,elem,face]=s2m(v,f,keepratio,maxvol,method)
%
% [node,elem,face]=s2m(v,f,keepratio,maxvol)
%
% volumetric mesh generation from a closed surface, shortcut for surf2mesh
%
% author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
%
% inputs and outputs are similar to those defined in surf2mesh
%
% if method='cgalpoly', s2m will call cgals2m and keepratio should be a
% structure (as the 'opt' input in cgals2m)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

p0=min(v(:,1:3));
p1=max(v(:,1:3));
if(nargin>=5)
  if(strcmp(method,'cgalpoly'))
    [node,elem,face]=cgals2m(v,f,keepratio,maxvol);
    return;
  end
end
[node,elem,face]=surf2mesh(v,f,p0,p1,keepratio,maxvol,[],[]);

repo: REST


text:
function obj = read_from_file(obj)
% Reads data from image filenames into obj.dat
%
% Try obj = check_image_filenames(obj) first.
%
% This is automatically called if you create a new image_vector object with
% names but do not directly enter data. e.g., the commands below will load data:
%   - name = 'salientmap.nii';
%   - img = image_vector('image_names', name);
%

% disp('Reading image data into object .dat field.');

% if we have mask info, use that. Otherwise, create
if isempty(obj.volInfo)
    % disp('Creating mask info from first image and storing in .volInfo.');
    obj.volInfo = iimg_read_img(obj.fullpath(1, :), 2);
else
    % disp('Using mask and space-defining info already stored in .volInfo field.');
end

% Now extract the actual data from the mask
switch spm('Ver')

    case {'SPM8', 'SPM5','SPM12'}
        obj.dat = iimg_get_data(obj.volInfo, obj.fullpath, 'single', 'noexpand')';

    case {'SPM2', 'SPM99'}
        % legacy, for old SPM
        obj.dat = iimg_get_data(obj.volInfo, obj.fullpath, 'single')';

    otherwise
        error('Unknown version of SPM! Update code, check path, etc.');
end

end % function


text:
function pcut = addEvalVariableCuts(p)

pcut = p;
if ~isempty(p.evalMap)
    pcut = emptyNumericalModel;
    for i = 1:length(p.evalMap)
        y = p.evalVariables(i);
        x = p.evalMap{i}.variableIndex;
        xL = p.lb(x);
        xU = p.ub(x);

        % Generate a convex hull polytope
        if xL<xU
            % If there is no convex hull method available, we might
            % generate one if we are in a a region where convexity is known
            if isempty(p.evalMap{i}.properties.convexhull)
                p.evalMap{i}.properties.convexhull = createConvexHullMethod(p.evalMap{i},xL,xU);
                remove_auto_generated_convexhull = 1;
            else
                remove_auto_generated_convexhull = 0;
            end
            if ~isempty(p.evalMap{i}.properties.convexhull)
                % A convex hull generator function is available!
                % Might be able to reuse hull from last run node
                if isfield(p.evalMap{i},'oldhull') && isequal(p.evalMap{i}.oldhull.xL,xL) && isequal(p.evalMap{i}.oldhull.xU,xU)
                    [Ax,Ay,b,K] = getOldHull(p,i);
                else
                    [Ax,Ay,b,K,p] = updateHull(xL,xU,p,i);
                    if isempty(Ax)
                        % Operator bounder does not cover this interval so
                        % use the sample-based instead
                        [Ax,Ay,b,K,p] = convexhullSampled(xL,xU,p,i);
                    end
                end
                if remove_auto_generated_convexhull
                    p.evalMap{i}.properties.convexhull = [];
                end
            elseif ~(any(isinf(xL)) | any(isinf(xU)))
                [Ax,Ay,b,K] = convexhullSampled(xL,xU,p,i);
            else
                Ax = [];
                Ay = [];
                b = [];
                K = [];
            end
            if ~isempty(b)
                removeThese = find(any(isnan([Ax Ay b]),2) | any(isinf([Ax Ay b]),2));
                if any(removeThese)
                    Ax(removeThese,:) = [];
                    Ay(removeThese,:) = [];
                    b(removeThese) = [];
                    K.l = K.l - nnz(removeThese > K.f);
                    K.f = K.f - nnz(removeThese <= K.f);
                end
                F_structemp = zeros(size(b,1),length(p.c)+1);
                F_structemp(:,1+y) = -Ay;
                F_structemp(:,1+x) = -Ax;
                F_structemp(:,1) = b;
                localModel = createNumericalModel(F_structemp,K);
                pcut = mergeNumericalModels(pcut,localModel);
            end
        end
    end

    pcut = mergeNumericalModels(p,pcut);
end

function [Ax,Ay,b,K] = getOldHull(p,i)

Ax = p.evalMap{i}.oldhull.Ax;
Ay = p.evalMap{i}.oldhull.Ay;
b = p.evalMap{i}.oldhull.b;
K = p.evalMap{i}.oldhull.K;

function [Ax,Ay,b,K,p] = updateHull(xL,xU,p,i)

if strcmpi(p.evalMap{i}.fcn,'blackbox')
    [Ax,Ay,b,K]=feval(p.evalMap{i}.properties.convexhull,xL,xU);
else
    [Ax,Ay,b,K]=feval(p.evalMap{i}.properties.convexhull,xL,xU, p.evalMap{i}.arg{2:end-1});
end
if ~isempty(Ax)
    removeThese = find(any(isinf([Ax Ay b]),2) | any(isnan([Ax Ay b]),2));
    if any(removeThese)
        Ax(removeThese,:) = [];
        Ay(removeThese,:) = [];
        b(removeThese) = [];
        K.l = K.l - nnz(removeThese > K.f);
        K.f = K.f - nnz(removeThese <= K.f);
    end
end
p = saveOldHull(xL,xU,Ax,Ay,b,K,p,i);

function p = saveOldHull(xL,xU,Ax,Ay,b,K,p,i)
p.evalMap{i}.oldhull.xL = xL;
p.evalMap{i}.oldhull.xU = xU;
p.evalMap{i}.oldhull.Ax = Ax;
p.evalMap{i}.oldhull.Ay = Ay;
p.evalMap{i}.oldhull.b = b;
p.evalMap{i}.oldhull.K = K;

function [Ax,Ay,b,K,p] = convexhullSampled(xL,xU,p,i)

if length(xL)>1
    Ax = [];
    Ay = [];
    b = [];
    K = [];
    return
end
% sample function
z = linspace(xL,xU,100);

if isequal(p.evalMap{i}.fcn,'power_internal2')
    % Special code for automatically converting sigmonial
    % terms to be solvable with bmibnb
    fz = feval(p.evalMap{i}.fcn,z,p.evalMap{i}.arg{2});
else
    arg = p.evalMap{i}.arg;
    arg{1} = z;
    vectorized = 1;
    fz = real(feval(p.evalMap{i}.fcn,arg{1:end-1}));
    if length(fz) ~= length(z)
        % Operator is not vectorized?
        vectorized = 0;
        fz = [];
        for k = 1:length(z)
            arg{1} = z(k);
            fz = [fz real(feval(p.evalMap{i}.fcn,arg{1:end-1}))];
        end
    end
    if length(fz) ~= length(z)
        Ax = [];
        Ay = [];
        b = [];
        K = [];
        return
    end
    if isa(p.evalMap{i}.properties.stationary,'function_handle')
        [xS,fS] = feval(p.evalMap{i}.properties.stationary,xL,xU);
        if ~isempty(xS)
            % Add stationary points to the list of samples
            z = [z xS(:)'];
            fz = [fz fS(:)'];
        end
    end
    [minval,minpos] = min(fz);
    [maxval,maxpos] = max(fz);
    xtestmin = linspace(z(max([1 minpos-5])),z(min([100 minpos+5])),100);
    xtestmax = linspace(z(max([1 maxpos-5])),z(min([100 maxpos+5])),100);
    if vectorized
        arg{1} = xtestmin;
        fz1 = real(feval(p.evalMap{i}.fcn,arg{1:end-1}));
        arg{1} = xtestmax;
        fz2 = real(feval(p.evalMap{i}.fcn,arg{1:end-1}));
        z = [z(:);xtestmin(:);xtestmax(:)];
        fz = [fz(:);fz1(:);fz2(:)];
    else
        fz1 = [];
        for k = 1:length(xtestmin)
            arg{1} = xtestmin(k);
            fz1 = [fz1 real(feval(p.evalMap{i}.fcn,arg{1:end-1}))];
        end
        fz2 = [];
        for k = 1:length(xtestmax)
            arg{1} = xtestmax(k);
            fz2 = [fz2 real(feval(p.evalMap{i}.fcn,arg{1:end-1}))];
        end
    end
    [z,sorter] = sort(z);
    fz = fz(sorter);
    [z,ii,jj]=unique(z);
    fz = fz(ii);
end
if max(fz)-min(fz) < 1e-4
    Ax = [];
    Ay = [];
    b = [];
    K = [];
    return
end
[Ax,Ay,b,K] = convexhullFromSampled(z,fz,xL,xU);
removeThese = find(any(isnan([Ax Ay b]),2) | any(isinf([Ax Ay b]),2));
Ax(removeThese,:) = [];
Ay(removeThese,:) = [];
b(removeThese) = [];
K.l = K.l - length(removeThese);%[];
p = saveOldHull(xL,xU,Ax,Ay,b,K,p,i);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function hFig = bst_warp_prepare(ChannelFile, Options)
% BST_WARP_PREPARE:  Set up everything for a call to bst_warp (deformation of default anatomy)
%
% USAGE:  bst_warp_prepare(ChannelFile, Options)
%         bst_warp_prepare(ChannelFile)
%
% INPUTS:
%     - ChannelFile : Channel file that contains the HeadPoints for the deformation
%     - Options     : Structure of the options
%         |- tolerance    : Percentage of outliers head points, ignored in the calulation of the deformation.
%         |                 Set to more than 0 when you know your head points have some outliers.
%         |                 If not specified: asked to the user (default
%         |- isInterp     : If 0, do not do a full interpolation (default: 1)
%         |- isScaleOnly  : If 1, do not perform the full deformation but only a linear scaling in the three directions (default: 0)
%         |- isSurfaceOnly: If 1, do not warp/scale the MRI (default: 0)
%         |- isInteractive: If 0, do not ask anything to the user (default: 1)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2010-2022

%% ===== CHECK PARAMETERS =====
% Default options
defOptions = struct(...
    'tolerance',     [], ...
    'isInterp',      1, ...
    'isScaleOnly',   [], ...
    'isSurfaceOnly', 0, ...
    'isInteractive', 1);
% Copy default options
if (nargin < 2)
    Options = defOptions;
else
    Options = struct_copy_fields(Options, defOptions, 0);
end
% Additional non-interactive defaults
if ~Options.isInteractive && isempty(Options.isScaleOnly)
    Options.isScaleOnly = 0;
end
if ~Options.isInteractive && isempty(Options.tolerance)
    Options.tolerance = 2;
end
hFig = [];

%% ===== GET DEFAULT ANAT =====
% Start by unloading everything
bst_memory('UnloadAll', 'Forced');
% Progress bar
if Options.isInteractive
    bst_progress('start', 'Warp anatomy', 'Initialization...');
end
% Get default anatomy
sDefSubject = bst_get('Subject', 0);
% If no default anatomy defined: error
if isempty(sDefSubject.Anatomy) || isempty(sDefSubject.Surface) || isempty(sDefSubject.iScalp)
    local_error('Please define a default anatomy before doing this.', Options.isInteractive);
    return;
end
% Load fiducials from default MRI
srcMriFile = file_fullpath(sDefSubject.Anatomy(sDefSubject.iAnatomy).FileName);
srcMriMat = load(srcMriFile, 'SCS', 'Voxsize');
% Check whether fiducials are defined on MRI
if ~isfield(srcMriMat, 'SCS') || ~isfield(srcMriMat.SCS, 'NAS') || isempty(srcMriMat.SCS.NAS) || isempty(srcMriMat.SCS.LPA) || isempty(srcMriMat.SCS.RPA)
    local_error('Fiducials are not defined on default MRI.', Options.isInteractive);
    return;
end

%% ===== GET DIGITIZED HEAD POINTS =====
% Get head points
HeadPoints = channel_get_headpoints(ChannelFile, 1);
% Number of head points
if isempty(HeadPoints)
    nhp = 0;
    destPts = [];
else
    nhp = length(HeadPoints.Loc);
    destPts = HeadPoints.Loc';
end
% Not enough head points
if (nhp < 3)
    local_error('Not enough digitized head points to perform this operation.', Options.isInteractive);
    return;
% Few head points: Scale anatomy ?
elseif Options.isInteractive && isempty(Options.isScaleOnly)
    if (nhp < 70)
        msg = ['There are only ' num2str(nhp) ' digitized head points, this might be ' 10 ...
               'insufficient to deform the default anatomy precisely. ' 10 10 ...
               'You could instead perform a simple scaling of the anatomy.' 10 10];
    else
        msg = ['Do you want to perform a simple scaling (Scale)' 10 ...
               'or a deformation (Warp) of the default anatomy?' 10 10];
    end
    res = java_dialog('question', msg, 'Warp anatomy', [], {'Warp', 'Scale', 'Cancel'}, 'Scale');
    if isempty(res) || strcmpi(res, 'Cancel')
        bst_progress('stop');
        return;
    elseif strcmpi(res, 'Scale')
        Options.isScaleOnly = 1;
    else
        Options.isScaleOnly = 0;
    end
    if (nhp < 70)
        if ~java_dialog('confirm', ...
              ['There are only ' num2str(nhp) ' digitized head points, this might be ' 10 ...
               'insufficient to deform the default anatomy precisely. ' 10 10 ...
               'Proceed anyway?'], 'Warp default anatomy');
            bst_progress('stop');
            return;
        end
    end
end

%% ===== GET SUBJECT ANAT =====
% Get study for this channel file
sStudy = bst_get('ChannelFile', ChannelFile);
% Get subject for this channel file
[sSubject, iSubject] = bst_get('Subject', sStudy.BrainStormSubject, 1);
% If there is already an anatomy defined for this subject: warning
if (~isempty(sSubject.Anatomy) || ~isempty(sSubject.Surface))
    if Options.isInteractive
        isOk = java_dialog('confirm', ['Warning: there is already an anatomy defined for this subject.' 10 10 ...
                                       'This process will delete the existing anatomy, and replace it with ' 10 ...
                                       'a deformed version of the default anatomy.' 10 10 ...
                                       'Delete current subject anatomy ?' 10 10], ...
                                       'Warp default anatomy');
        if ~isOk
            bst_progress('stop');
            return
        end
    end
    % Delete previous surfaces
    if ~isempty(sSubject.Surface)
        SurfaceFiles = cellfun(@(c)file_fullpath(c), {sSubject.Surface.FileName}, 'UniformOutput', 0);
        file_delete(SurfaceFiles, 1);
    end
    if ~isempty(sSubject.Anatomy)
        MriFiles = cellfun(@(c)file_fullpath(c), {sSubject.Anatomy.FileName}, 'UniformOutput', 0);
        file_delete(MriFiles, 1);
    end
end
% Force subject to us default anatomy
s.UseDefaultAnat = 0;
s.Anatomy = [];
s.Surface = [];
% Update subject file
bst_save(file_fullpath(sSubject.FileName), s, 'v7', 1);
% Load default scalp
srcScalpFile = sDefSubject.Surface(sDefSubject.iScalp).FileName;
srcSurf = in_tess_bst(srcScalpFile);
% Display initial surface and head points
if Options.isInteractive
    view_headpoints(ChannelFile, srcScalpFile);
end

%% ===== PROJECT HEAD POINTS =====
% Compute center of mass of the head points, use as the center of projection
if (length(destPts) > 50)
    center = mean(destPts);
else
    center = [0 0 0];
end
% Project digitized head points on the scalp
[destPtsProj, dist] = project_on_surface(srcSurf, destPts, center);

%% ===== REMOVING OUTLIERS =====
if (length(destPts) > 15)
    % Ask the user the tolerance for outliers
    if Options.isInteractive && isempty(Options.tolerance)
        res = java_dialog('input', ['You can choose to ignore the digitized head points that are far ' 10 ...
                                    'away from the scalp in the calculation of the deformation field.' 10 10 ...
                                    'Percentage of head points to ignore [0-100]:'], ...
                                    'Warp default anatomy', [], '2');
        if isempty(res) || isnan(str2double(res))
            bst_progress('stop');
            return
        end
        Options.tolerance = str2double(res) / 100;
    end
    % Remove points that have the maximum distances
    if (Options.tolerance > 0)
        % Number of points to remove
        nRemove = ceil(Options.tolerance * length(destPts));
        disp(['WARP> Remove ' num2str(nRemove) ' points from the head points.']);
        % Sort points by distance to scalp
        [tmp__, iSort] = sort(dist, 1, 'descend');
        iRemove = iSort(1:nRemove);
        % Remove from list of destination points
        destPts(iRemove,:) = [];
        destPtsProj(iRemove,:) = [];
    end
end

%% ===== BUILD TRANSFORMATION =====
% Full deformation
if ~Options.isScaleOnly
    if Options.isInterp
        % Destination landmarks: Fit spherical harmonic to the digitized head points
        fvh = hsdig2fv(destPts, 5, 15/1000, 40*pi/180, 0);
        destPtsParam = fvh.vertices;
        % Source landmarks: Project remeshed digitized surface on scalp
        srcPtsParam = project_on_surface(srcSurf, destPtsParam, center);
    else
        % Destination landmarks: head points
        destPtsParam = destPts;
        % Source landmarks: head points projected on the scalp
        srcPtsParam = destPtsProj;
    end

% Scaling only
else
    % Downsample heavily the scalp surface
    [srcFacesParam, srcPtsParam] = reducepatch(srcSurf.Faces, srcSurf.Vertices, 300 / length(srcSurf.Vertices));
    % Destination surface: scale version of the downsampled scalp surface
    destPtsParam = srcPtsParam;
    % Scaling in each direction (x,y,z)
    for dim = 1:3
        % Compute scale factor
        fall = destPts(:,dim) ./ destPtsProj(:,dim);
        fall(isnan(fall) | (fall < .5) | (fall > 1.5)) = [];
        % Apply scale factor
        if ~isempty(fall)
            f = mean(fall);
            destPtsParam(:,dim) = f .* destPtsParam(:,dim);
        end
    end
end

%% ===== WARP =====
% Create list of files to process
SurfaceFiles     = {sDefSubject.Surface.FileName};
SurfaceFilesFull = cellfun(@(c)file_fullpath(c), SurfaceFiles, 'UniformOutput', 0);
MriFiles         = {sDefSubject.Anatomy([sDefSubject.iAnatomy, setdiff(1:length(sDefSubject.Anatomy), sDefSubject.iAnatomy)]).FileName};
MriFilesFull     = cellfun(@(c)file_fullpath(c), MriFiles, 'UniformOutput', 0);
OutputTag        = '_warped';
OutputDir        = bst_fileparts(file_fullpath(sSubject.FileName));
% Main call
bst_warp(destPtsParam, srcPtsParam, SurfaceFilesFull, MriFilesFull, OutputTag, OutputDir, Options.isSurfaceOnly);

%% ===== COPY ATLASES =====
% Get list of atlas files (scout files in the default anat folder)
atlasDir = bst_fileparts(file_fullpath(sDefSubject.FileName));
dirScout = dir(bst_fullfile(atlasDir, 'scout_*.mat'));
% Copy all the files file
for i = 1:length(dirScout)
    file_copy(bst_fullfile(atlasDir, dirScout(i).name), bst_fullfile(OutputDir, dirScout(i).name))
end

%% ===== UPDATE DATABASE =====
% Reload subject
db_reload_subjects(iSubject);
% Unload all the surfaces, close all the figures
bst_memory('UnloadAll', 'Forced');
% Get subject again
sSubject = bst_get('Subject', iSubject);
% Display warp head and cortex surfaces
if Options.isInteractive && ~isempty(sSubject.iScalp)
    ScalpFile = sSubject.Surface(sSubject.iScalp).FileName;
    hFig = view_surface(ScalpFile);
end
if Options.isInteractive && ~isempty(sSubject.iCortex)
    CortexFile = sSubject.Surface(sSubject.iCortex).FileName;
    hFig = view_surface(CortexFile);
end

% Close progress bar
if Options.isInteractive
    bst_progress('stop');
end

end

%% ========================================================================
%  ===== HELPER FUNCTIONS =================================================
%  ========================================================================

%% ===== PROJECT ON SURFACE =====
function [Q, dist] = project_on_surface(sSurf, P, center)
    % Center on (0,0,0)
    P = bst_bsxfun(@minus, P, center);
    sSurf.Vertices = bst_bsxfun(@minus, sSurf.Vertices, center);
    % Project points on surface
    Q = 0 .* P;
    for i = 1:length(P)
        proj = tess_ray_intersect(sSurf.Vertices, sSurf.Faces, [0 0 0], P(i,:))';
        if isempty(proj)
            Q(i,:) = P(i,:);
        elseif (size(proj,1) > 1)
            Q(i,:) = proj(1,:);
        else
            Q(i,:) = proj;
        end
    end
    % Compute the distance between each point and its projection (distance to the scalp surface)
    if (nargout >= 2)
        dist = sum(sqrt((P - Q).^2), 2);
    end
    % Restore center
    Q = bst_bsxfun(@plus, Q, center);
end

%% ===== LOCAL ERROR PROCESSING =====
function local_error(errMsg, isInteractive)
    if isInteractive
        bst_error(errMsg, 'Warp default anatomy', 0);
        bst_progress('stop');
    else
        bst_report('Error', 'process_warp', [], errMsg);
    end
end

author: brainstorm-tools


text:
function lattice_rule_test06 ( )

%*****************************************************************************80
%
%% LATTICE_RULE_TEST06 tests FIBONACCI_LATTICE_Q3;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ian Sloan, Stephen Joe,
%    Lattice Methods for Multiple Integration,
%    Oxford, 1994.
%
  dim_num = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LATTICE_RULE_TEST06\n' );
  fprintf ( 1, '  FIBONACCI_LATTICE_Q3 applies a Fibonacci lattice rule\n' );
  fprintf ( 1, '  to integrate a function over the unit square.\n' );
  fprintf ( 1, '  A nonlinear coordinate transformation is applied.\n' );
  fprintf ( 1, '  These Fibonacci rules are only available in 2D.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The spatial dimension DIM_NUM = %d\n', dim_num );

  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  exact = e_01_2d ( dim_num, a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '         K         M    EXACT       ESTIMATE    ERROR\n' );
  fprintf ( 1, '\n' );

  for k = 3 : 18

    quad = fibonacci_lattice_q3 ( k, @f_01_2d );

    error = abs ( exact - quad );
    m = fibonacci ( k );

    fprintf ( 1, '  %8d  %8d  %10.6f  %10.6f  %10.6e\n', k, m, exact, quad, error );

  end

  return
end


text:
addpath('..\mexopencv\');
dir='G:\data\20130808\109ND800_080813\';
prefix=sprintf('%sDSC_', dir);

img0 = takeImage_v002(prefix, 1023, 'jpg', 3);
mask=uint8(ones(size(img0,1),size(img0,2)));
mask(:,:)=1;
excluded_band=100;
mask(1:excluded_band,:)=0;
mask(end-excluded_band+1:end,:)=0;
mask(:,1:excluded_band)=0;
mask(:,end-excluded_band+1:end)=0;

 prevPts = cv.goodFeaturesToTrack(img0,'MaxCorners', 200, ...
 'QualityLevel', 0.02,'MinDistance', 10,'Mask', mask, 'BlockSize', 3, 'UseHarrisDetector', 0, 'K',0.04);
prevPts=cv.cornerSubPix(img0, prevPts);

% the argument order and vacancy does not matter in matlab in constrast to
% opencv C++
corners=zeros(2,size(prevPts,2));
    for jerk=1:size(prevPts,2)
        if(~isempty(prevPts{jerk}))
   corners(:,jerk)= prevPts{jerk}';
        end
    end

% increase lambda apparently increase number of corners
%larger gaussian window sz2 will slightly decrease the number of corners

figure(1); clf;
imshow(img0);
hold on;            %# Add subsequent plots to the image
%plot(corners(2,:),corners(1,:),'go');  %# NOTE: x_p and y_p are switched (see note below)!
plot(corners(1,:),corners(2,:),'ro','MarkerSize',5);

hold off;
[m,n,p]=size(img0);
winsz1=7;
sigma1=10;

imglast=img0(:,:,1);
for i=1024:1055
    im=takeImage_v002(prefix,i, 'jpg', 3);
            maxedge=max(size(im));
                while(maxedge>1000)
                    im=impyramid(im, 'reduce');% 8UC1
                    maxedge=round(maxedge/2);
                end
%     img= double((read(myVid,i-1)));
%     img2=double((read(myVid,i)));
%     for band=1:size(img,3)
%         img(:,:,band)=simgauss(img(:,:,band), winsz1, sigma1);
%         img2(:,:,band)=simgauss(img2(:,:,band), winsz1, sigma1);
%     end
    [prevPts, status, err] = cv.calcOpticalFlowPyrLK(imglast, im, prevPts, 'WinSize', [23, 23], ...
        'MaxLevel', 2,'Criteria', struct('type', 'Count+EPS', 'maxCount', 20, 'epsilon', 0.01));
    for jerk=1:length(status)
        if(status~=0)
        end
    end
    corners=zeros(2,size(prevPts,2));
    for jerk=1:size(prevPts,2)
        if(~isempty(prevPts{jerk}))
   corners(:,jerk)= prevPts{jerk}';
        end
    end

   % corners=SimpKLT(img,img2,corners,3);
    imshow(im);
    hold on
    plot(corners(1,:),corners(2,:),'ro','MarkerSize',5);
    hold off
    drawnow;
   imglast=im;

end

% close(writerObj);

repo: ekfmonoslam


text:
function [out] = interflow_9(S,p1,p2,p3,dt)
%interflow_9

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

% Flux function
% ------------------
% Description:  Non-linear interflow if storage exceeds a threshold
% Constraints:  f <= S-p2
%               S-p2 >= 0     prevents numerical issues with complex numbers
% @(Inputs):    p1   - time coefficient [d-1]
%               p2   - storage threshold for flow generation [mm]
%               p3   - exponential scaling parameter [-]
%               S    - current storage [mm]
%               dt   - time step size [d]

out = min(max((S-p2)/dt,0),(p1.*max(S-p2,0)).^p3);

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function shadedplot(c,scale,ord)

   % Private method. See ../plot for details.

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   % PREP PLOT
   figure('Color','w','Position',[50 50 850 1100]);
   box on; hold on;

   % GET TIME MAX AND MINS
   tmin =  999999;
   tmax = -999999;
   wstartrels = c.relativeStartTime(ord);
   wlengths = [c.traces(ord).nsamples];
   wFs = [c.traces(ord).samplerate];

   for i = 1:numel(ord)
      %wstartrel = 86400*(get(c.W(i),'START_MATLAB')-c.trig(i));	% relative start time (trigger is at zero)
      %tr = wstartrel + [ 0:get(c.W(i),'DATA_LENGTH')-1]'/get(c.W(i),'Fs');
      tr = wstartrels(i) + ( 0:wlengths(i)-1)'/wFs(i);

      % save min and max relative trace times
      if tr(1) < tmin
         tmin = tr(1);
      end;
      if tr(end) > tmax
         tmax = tr(end);
      end;
   end;

   % MAKE ALIGNED DATA MATRIX
   %p = 1/get(c.W(1),'Fs');         % assumes all sampling periods are thesame
   p = 1/wFs(1);         % assumes all sampling periods are the same

   tmin = round(tmin*wFs(i))*p;      % round to nearest sample
   tmax = round(tmax*wFs(i))*p;
   N = 1:length(ord);
   T =  tmin : p : tmax ;
   D = zeros(length(N),length(T));
   count = 0;
   absmax = max(abs(c.traces(ord)));
   absmax(absmax==0) = scale; %next line will be negated for zero-scale

   for n=1:numel(ord)
      c.traces(ord(n)) = c.traces(ord(n)) .* (scale ./ absmax(n));
   end;
   d = double(c.traces); %nan fill??
   % d = double(c.W(ord) .* (scale ./ absmax));
   for i = 1:numel(ord)
      count = count + 1;
      %d = get(c.W(ord(i)),'DATA');            %%%d = c.w(:,i);
      %if (max(abs(d)) ~= 0)
      %    d = scale * d/max(abs(d));		% apply a uniform amplitude scale;
      %end
      % DONE ABOVE: wstartrel = 86400*(get(c.W(i),'START_MATLAB')-c.trig(i));	% relative start time (trigger is at zero)
      [~,startindex] = min(abs(T-wstartrels(i)));
      D( count ,  startindex : (startindex+wlengths(i)-1) ) = d(:,i);
   end;
   imagesc(T,N,D);

   % ADJUST PLOT
   axis([tmin tmax 0.5 length(ord)+0.5]);
   set(gca,'YDir','reverse');
   n = length(ord);
   set(gca,'YTick',(1:round(n/50):n));
   yt = get(gca,'YTick');
   set(gca,'YTickLabel',datestr(c.trig(ord(yt))),'FontSize',6);
   xlabel('Relative Time,(s)','FontSize',8);

   % SET COLOR MAP
   cmap = [ 0 0 1;
      1 1 1
      1 0 0];
   cmap = interp1([-1 0 1],cmap, -1:.1:1, 'linear');
   colormap(cmap);

   maybeReplaceYticksWithStationNames(c,gca)

   %PRINT OUT FIGURE
   set(gcf, 'paperorientation', 'portrait');
   set(gcf, 'paperposition', [.25 .25 8 10.5] );
   % print(gcf, '-depsc2', 'FIG_alignwfm.ps')
   %!ps2pdf FIG_alignwfm.ps
   %!convert FIG_alignwfm.ps FIG_alignwfm.gif
   %!rm FIG_alignwfm.ps
end

author: geoscience-community-codes

repo: GISMO


text:
%%**********************************************************************
%% blktrace: compute <X1,Z1> + ... + <Xp,Zp>
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%**********************************************************************

  function trXZ = blktrace(blk,X,Z,parbarrier);

  if (nargin == 3)
     trXZ = 0;
     for p = 1:size(blk,1)
        pblk = blk(p,:);
        if strcmp(pblk{1},'s')
           if (length(pblk{2}) == 1)
              trXZ = trXZ + sum(sum(X{p}.*Z{p}));
           else
              xx = mexsvec(pblk,X{p},0);
              zz = mexsvec(pblk,Z{p});
              trXZ = trXZ + xx'*zz;
           end
        else
           trXZ = trXZ + sum(X{p}.*Z{p});
        end
     end
  elseif (nargin == 4)
     trXZ = 0;
     for p = 1:size(blk,1)
        pblk = blk(p,:);
        if (norm(parbarrier{p}) == 0)
           if strcmp(pblk{1},'s')
              if (length(pblk{2}) == 1)
                 trXZ = trXZ + sum(sum(X{p}.*Z{p}));
              else
                 xx = mexsvec(pblk,X{p},0);
                 zz = mexsvec(pblk,Z{p});
                 trXZ = trXZ + xx'*zz;
              end
           else
              trXZ = trXZ + sum(X{p}.*Z{p});
           end
        else
           idx = find(parbarrier{p} == 0);
           if ~isempty(idx)
              if strcmp(pblk{1},'s')
                 sumXZ = sum(X{p}.*Z{p});
                 ss = [0,cumsum(pblk{2})];
                 for k = 1:length(idx)
                    idxtmp = [ss(idx(k))+1:ss(idx(k)+1)];
                    trXZ = trXZ + sum(sumXZ(idxtmp));
                 end
              elseif strcmp(pblk{1},'q')
	         tmp = qops(pblk,X{p},Z{p},1);
                 trXZ = trXZ + sum(tmp(idx));
              elseif strcmp(pblk{1},'l')
                 trXZ = trXZ + sum(X{p}(idx).*Z{p}(idx));
              end
           end
        end
     end
  end
%%**********************************************************************


text:
function varargout = carromtable(X)
% Carrom table function
%
%   CARROMTABLE([x1, x2]) returns the value of the Carrom table
%   function at the specified points. [x1] and [x2] may be vectors.
%   The search domain is
%
%               -10 < x_i < 10
%
%   The four global minimum are found near the corners of the interval,
%   with
%
%               fmin = 24.1568155.

% Author: Rody P.S. Oldenhuis
% Delft University of Technology
% E-mail: oldenhuis@dds.nl
% Last edited 20/Jul/2009

    % if no input is given, return dimensions, bounds and minimum
    if (nargin == 0)
        varargout{1} = 2;  % # dims
        varargout{2} = [-10, -10]; % LB
        varargout{3} = [+10, +10]; % UB
        varargout{4} = [+9.646157266348881e+000, +9.646134286497169e+000
                        -9.646157266348881e+000, +9.646134286497169e+000
                        +9.646157266348881e+000, -9.646134286497169e+000
                        -9.646157266348881e+000, -9.646134286497169e+000]; % solution
        varargout{5} = -2.415681551650653e+001; % function value at solution

    % otherwise, output function value
    else

        % keep values in the serach interval
        X(X < -10) = inf;     X(X > 10) = inf;

        % split input vector X into x1, x2
        if size(X, 1) == 2
            x1 = X(1, :);        x2 = X(2, :);
        else
            x1 = X(:, 1);        x2 = X(:, 2);
        end

        % output function value
        varargout{1} = -((cos(x1).*cos(x2).*exp(abs(1 - sqrt(x1.^2 + x2.^2)/pi))).^2)/30;

    end

end

author: Sable

repo: mcbench-benchmarks


text:
%% Example using rbb3select

% This part from MATLAB(R)'s SCATTER3 example
[x,y,z] = sphere(16);
X = [x(:)*.5 x(:)*.75 x(:)];
Y = [y(:)*.5 y(:)*.75 y(:)];
Z = [z(:)*.5 z(:)*.75 z(:)];
S = repmat([1 .75 .5]*10,prod(size(x)),1);
C = repmat([1 2 3],prod(size(x)),1);
scatter3(X(:),Y(:),Z(:),S(:),C(:),'filled'), view(-60,60)

% Draw red circles around selected points
SS=rbb3select(X,Y,Z);
hold on
scatter3(X(find(SS)),Y(find(SS)),Z(find(SS)),'o','r')
hold off


text:
classdef MimExportEditedImage < MimGuiPlugin
    % MimExportEditedImage. Gui Plugin
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimExportEditedImage is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Export Edits'
        SelectedText = 'Export Edits'
        ToolTip = 'Exports the current edits to an external file'
        Category = 'Save / load edits'
        Visibility = 'Overlay'
        Mode = 'Edit'

        Icon = 'export_overlay.png'
        Location = 28

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.GetMode.ExportEdit();
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded && gui_app.ImagePanel.OverlayImage.ImageExists && ...
                isequal(gui_app.GetCurrentModeName, MimModes.EditMode) && gui_app.IsTabEnabled('Edit');
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [Simil2]=fuzzysimil2(A,B)

PA=sqrt((A(1)-A(2))^2+A(5)^2)+sqrt((A(3)-A(4))^2+A(5)^2)+(A(3)-A(2))+(A(4)-A(1));
PB=sqrt((B(1)-B(2))^2+B(5)^2)+sqrt((B(3)-B(4))^2+B(5)^2)+(B(3)-B(2))+(B(4)-B(1));
aA=1/2*(A(5)*(A(3)-A(2)+A(4)-A(1)));
aB=1/2*(B(5)*(B(3)-B(2)+B(4)-B(1)));
temp=(min(PA,PB)/max(PA,PB))*(min(aA,aB) + min(A(5),B(5)))/(max(aA,aB) + max(A(5),B(5)));

Simil2=(1-(sum(abs(A(1:4)-B(1:4))))/4)*temp;


text:
function c = mtimes(a,b)
	checkinputs(a,b);
	if numel(a)>1 && numel(b)>1
		error('Only elementwise multiplication supported.');
	end
	c = uint64times(uint64(a),uint64(b));
end


text:
function D = compute_deviator_tensor(T)

% compute_deviator_tensor - compute trace free tensor
%
%   D = compute_deviator_tensor(T);
%
%   D = T - trace(T)/2*Id
%
%   Copyright (c) Gabriel Peyre

n = size(T,1);

if size(T,3)==2 && size(T,4)==2
    t = (T(:,:,1)+T(:,:,4))/2;
    D = cat(3, t, zeros(n), zeros(n), t);
    D = T - reshape(D,[n n 2 2]);
elseif size(T,3)==3 && size(T,4)==1
    t = (T(:,:,1)+T(:,:,2))/2;
    D = cat(3, t, t, zeros(n));
    D = T - D;
else
    error('Wrong size');
end

author: gpeyre

repo: matlab-toolboxes


text:
% Test file for trigtech/vals2coeffs.m

function pass = test_vals2coeffs(varargin)

% Set a tolerance (pref.chebfuneps doesn't matter)
tol = 100*eps;

%%
% Test that a single value is converted correctly
vals = sqrt(2);
c = trigtech.vals2coeffs(vals);
pass(1) = (vals == c);

%%
% Simple data (odd case)
% Exact values
f_exact = @(x) 1 + cos(2*pi*x);
x_pts = trigtech.trigpts(5);
vals = f_exact(x_pts);
% Exact coefficients
cTrue = [0.5 0 1 0 0.5].';

%%
% Test real branch
c = trigtech.vals2coeffs(vals);
pass(2) = norm(c - cTrue, inf) < tol;

%%
% Test imaginary branch
c = trigtech.vals2coeffs(1i*vals);
pass(3) = norm(c - 1i*cTrue, inf) < tol;

%%
% Test general branch
c = trigtech.vals2coeffs((1+1i)*vals);
pass(4) = norm(c - (1+1i)*cTrue, inf) < tol;

%%
% Test for array input
c = trigtech.vals2coeffs([vals, -vals]);
pass(5) = norm(c(:,1) - cTrue, inf) < tol && ...
          norm(c(:,2) + cTrue, inf) < tol;

%%
% Simple data (even case)
x_pts = trigtech.trigpts(6);
% Exact values
f_exact = @(x) 1 + sin(2*pi*x) + cos(3*pi*x);
vals = f_exact(x_pts);
% Exact coefficients
cTrue = [1 0.5i 0 1 0 -0.5i].';

%%
% Test real branch
c = trigtech.vals2coeffs(vals);
pass(6) = norm(c - cTrue, inf) < tol;

%%
% Test imaginary branch
c = trigtech.vals2coeffs(1i*vals);
pass(7) = norm(c - 1i*cTrue, inf) < tol;

%%
% Test general branch
c = trigtech.vals2coeffs((1+1i)*vals);
pass(8) = norm(c - (1+1i)*cTrue, inf) < tol;

%%
% Test for array input
c = trigtech.vals2coeffs([vals, -vals]);
pass(9) = norm(c(:,1) - cTrue, inf) < tol && ...
          norm(c(:,2) + cTrue, inf) < tol;

%%
% Test for symmetry
x = trigpts(123);
vals = [cos(pi*x),sin(pi*x),cos(pi*x)+sin(pi*x)];
vals(1,2) = 0;
c = trigtech.vals2coeffs(vals);
pass(10) = norm(c(:,1) - flipud(c(:,1)), inf) == 0 && ...
           norm(c(:,2) + flipud(c(:,2)), inf) == 0 && ...
           norm(c(:,3) - flipud(c(:,3)), inf) > 0 && ...
           norm(c(:,3) + flipud(c(:,3)), inf) > 0;

end

author: chebfun

repo: chebfun


text:
function P = cpf_p(parameterization, step, z, V, lam, Vprv, lamprv, pv, pq)
%CPF_P Computes the value of the CPF parameterization function.
%   P = CPF_P(PARAMETERIZATION, STEP, Z, V, LAM, VPRV, LAMPRV, PV, PQ)
%
%   Computes the value of the parameterization function at the current
%   solution point.
%
%   Inputs:
%       PARAMETERIZATION : Value of cpf.parameterization option
%       STEP : continuation step size
%       Z : normalized tangent prediction vector from previous step
%       V : complex bus voltage vector at current solution
%       LAM : scalar lambda value at current solution
%       VPRV : complex bus voltage vector at previous solution
%       LAMPRV : scalar lambda value at previous solution
%       PV : vector of indices of PV buses
%       PQ : vector of indices of PQ buses
%
%   Outputs:
%       P : value of the parameterization function at the current point
%
%   See also CPF_PREDICTOR, CPF_CORRECTOR.

%   MATPOWER
%   Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC)
%   by Shrirang Abhyankar, Argonne National Laboratory
%   and Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% evaluate P(x0, lambda0)
if parameterization == 1        %% natural
    if lam >= lamprv
        P = lam - lamprv - step;
    else
        P = lamprv - lam - step;
    end
elseif parameterization == 2    %% arc length
    Va = angle(V);
    Vm = abs(V);
    Vaprv = angle(Vprv);
    Vmprv = abs(Vprv);
    P = sum(([Va([pv; pq]); Vm(pq); lam] - [Vaprv([pv; pq]); Vmprv(pq); lamprv]).^2) - step^2;
elseif parameterization == 3    %% pseudo arc length
    nb = length(V);
    Va = angle(V);
    Vm = abs(V);
    Vaprv = angle(Vprv);
    Vmprv = abs(Vprv);
    P = z([pv; pq; nb+pq; 2*nb+1])' * ...
        ( [Va([pv; pq]); Vm(pq); lam] - [Vaprv([pv; pq]); Vmprv(pq); lamprv] )...
        - step;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function ubj=saveubjson(rootname,obj,varargin)
%
% ubj=saveubjson(obj)
%    or
% ubj=saveubjson(rootname,obj,filename)
% ubj=saveubjson(rootname,obj,opt)
% ubj=saveubjson(rootname,obj,'param1',value1,'param2',value2,...)
%
% Convert a MATLAB object  (cell, struct, array, table, map, graphs ...)
% into a Universal Binary JSON (UBJSON, Draft-12) or a MessagePack binary stream
%
% author: Qianqian Fang (q.fang <at> neu.edu)
% initially created on 2013/08/17
%
% Format specifications:
%    Binary JData (BJData):https://github.com/NeuroJSON/bjdata
%    UBJSON:               https://github.com/ubjson/universal-binary-json
%    MessagePack:          https://github.com/msgpack/msgpack
%
% This function is the same as calling "savebj(..,'ubjson',1,'endian','B')"
% By default this function creates UBJSON-compliant output without the
% newly added uint16(u), uint32(m), uint64(M) and half-precision float (h)
% data types and use Big-Endian for all numerical values as in UBJSON
% Draft-12.
%
% This function by default still enables an optimized ND-array format for efficient
% array storage. To ensure the output compatible to UBJSON Draft-12, one should use
% "saveubjson(...,'NestArray',1)" or "savebj(...,'ubjson',1,'NestArray',1)"
%
% input:
%      rootname: the name of the root-object, when set to '', the root name
%           is ignored, however, when opt.ForceRootName is set to 1 (see below),
%           the MATLAB variable name will be used as the root name.
%      obj: a MATLAB object (array, cell, cell array, struct, struct array,
%           class instance)
%      filename: a string for the file name to save the output UBJSON data
%      opt: a struct for additional options, ignore to use default values.
%           opt can have the following fields (first in [.|.] is the default)
%
%           opt can be replaced by a list of ('param',value) pairs. The param
%           string is equivallent to a field in opt and is case sensitive.
%
%           Please type "help savebj" for details for all supported options.
%
% output:
%      ubj: a binary string in the UBJSON format (see http://ubjson.org)
%
% examples:
%      jsonmesh=struct('MeshVertex3',[0 0 0;1 0 0;0 1 0;1 1 0;0 0 1;1 0 1;0 1 1;1 1 1],...
%               'MeshTet4',[1 2 4 8;1 3 4 8;1 2 6 8;1 5 6 8;1 5 7 8;1 3 7 8],...
%               'MeshTri3',[1 2 4;1 2 6;1 3 4;1 3 7;1 5 6;1 5 7;...
%                          2 8 4;2 8 6;3 8 4;3 8 7;5 8 6;5 8 7],...
%               'MeshCreator','FangQ','MeshTitle','T6 Cube',...
%               'SpecialData',[nan, inf, -inf]);
%      saveubjson(jsonmesh)
%      saveubjson('',jsonmesh,'meshdata.ubj')
%      saveubjson('mesh1',jsonmesh,'FileName','meshdata.msgpk','MessagePack',1)
%      saveubjson('',jsonmesh,'KeepType',1)
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

if(nargin==1)
    ubj=savebj('',rootname,'ubjson',1,'endian','B');
elseif(length(varargin)==1 && ischar(varargin{1}))
    ubj=savebj(rootname,obj,'FileName',varargin{1},'ubjson',1,'endian','B');
else
    ubj=savebj(rootname,obj,varargin{:},'ubjson',1,'endian','B');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% istypeof
% Check the type of an object.

%%% Syntax
%  truth = istypeof(obj, typename)

%%% Description
% |istypeof(obj, typename)| returns |true| if |obj| is a scalar or array of the
% type described by a string |typename|. Exmaples of |typename| are |'int'|,
% |'real'|, |'complex'|, |'logical'| (for 'true' and 'false'), and the names of classes.
%
% To test if |obj| is a cell array whose elements are scalars or arrays of a
% specific type, a _cell type_ can be specified as |typename|.  A _cell_type_ is
% a conventional type name followed by 'cell', such as |'intcell'|,
% |'realcell'|, |'complexcell'|.
%
% For |obj == []|, the function returns |true| for any |typename| that is a
% non-cell type. For |obj == {}|, the function returns |true| for any |typename|
% that is a cell type.

%%% Example
%   % Test int.
%   none = []; s = 1; v = 1:10; m = ones(3, 3); type = 'int';
%   fprintf('Test 1: are all elements of the objects %s?  ', type);
%   if istypeof(none, type) && istypeof(s, type) && istypeof(v, type) && istypeof(m, type)
%       fprintf('Yes.\n');
%   else
%       fprintf('No.\n');
%   end
%
%   % Test real.
%   none = []; s = 2.4; v = linspace(1, 10, 17); m = [pi 1; exp(1) 0]; type = 'real';
%   fprintf('Test 2: are all elements of the objects %s?  ', type);
%   if istypeof(none, type) && istypeof(s, type) && istypeof(v, type) && istypeof(m, type)
%       fprintf('Yes.\n');
%   else
%       fprintf('No.\n');
%   end
%
%   % Test Axis.
%   none = []; s = Axis.x; v = Axis.elems; m = [Axis.x Axis.y; Axis.y Axis.z]; type = 'Axis';
%   fprintf('Test 3: are all elements of the objects %s?  ', type);
%   if istypeof(none, type) && istypeof(s, type) && istypeof(v, type) && istypeof(m, type)
%       fprintf('Yes.\n');
%   else
%       fprintf('No.\n');
%   end
%
%   % Test a cell type.
%   none = {}; c = {1, [2 3]; 6, [7 8; 9 10]}; type = 'intcell';
%   fprintf('Test 4: are all objects %s?  ', type);
%   if istypeof(none, type) && istypeof(c, type)
%       fprintf('Yes.\n');
%   else
%       fprintf('No.\n');
%   end

%%% See Also
% <istypesizeof.html istypesizeof>

function truth = istypeof(obj, typename)

chkarg(ischar(typename), '"typename" should be string.');
iscelltype = false;
if ~isempty(strfind(typename, 'cell'))  && isequal(typename(end-3:end), 'cell')  % typename ends with 'cell'
	typename = typename(1:end-4);
	iscelltype = true;
end

if isempty(obj)
	truth = ~xor(iscelltype, iscell(obj));  % true if "typename" is a cell type and obj == {}, or if "typename" is a non-cell type and obj == []
elseif iscelltype
	truth = iscell(obj);
	obj = obj(:);
	for i = 1:length(obj)
		truth = truth && istypeof(obj{i}, typename);
	end
elseif isequal(typename, 'int')
	truth = isnumeric(obj) && isint(obj);  % note isint('abc') == true
elseif isequal(typename, 'real')
	truth = isnumeric(obj) && isreal(obj);  % note isreal('abc') == true
elseif isequal(typename, 'complex')
	truth = isnumeric(obj);
elseif isequal(typename, 'arbitrary')
	truth = true;
else
	truth = isa(obj, typename);
end

function truth = isint(obj)

truth = isreal(obj) && all(mod(obj(:),1)==0);

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function test_suite = test_zinegbin

%   Run specific demo, save values and compare the results to the expected.
%   Works for both xUnit Test Framework package by Steve Eddins and for
%   the built-in Unit Testing Framework (as of Matlab version 2013b).
%
%   See also
%     TEST_ALL, DEMO_ZINEGBIN
%
% Copyright (c) 2011-2012 Ville Tolvanen
% Copyright (c) 2014 Tuomas Sivula

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  % Check which package to use
  if exist('initTestSuite', 'file')
    initTestSuite_path = which('initTestSuite');
  else
    initTestSuite_path = '';
  end
  if ~isempty(initTestSuite_path) && ...
     exist([initTestSuite_path(1:end-15) 'runtests'], 'file')
    % xUnit package
    initTestSuite;
  else
    % Built-in package
    % Use all functions except the @setup
    tests = localfunctions;
    tests = tests(~cellfun(@(x)strcmp(func2str(x),'setup'),tests));
    test_suite = functiontests(tests);
  end
end

% -------------
%     Tests
% -------------

function testRunDemo(testCase)
  % Run the correspondin demo and save the values. Note this test has to
  % be run at lest once before the other test may succeed.
  rundemo(getName(), ...
    {{'Ef' @(x)x(1:100)}, {'Varf', @(x)diag(x(1:100,1:100))}})
end

function testPredictionsZinegbin(testCase)
  verifyVarsEqual(testCase, getName(), {'Ef', 'Varf'}, ...
    'RelTolElement', 0.1, 'RelTolRange', 0.02)
end

% ------------------------
%     Helper functions
% ------------------------

function testCase = setup
  % Helper function to suply empty array into variable testCase as an
  % argument for each test function, if using xUnit package. Not to be
  % used with built-in test framework.
  testCase = [];
end

function name = getName
  % Helperfunction that returns the name of the demo, e.g. 'binomial1'.
  name = mfilename;
  name = name(6:end);
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
classdef medianFilter < EBSDFilter

  properties
    numNeighbours % number of neigbours to consider (default 1)
  end

  methods

    function F = medianFilter(varargin)
      %

      F.numNeighbours = get_option(varargin,'neighbours',1);

      addlistener(F,'isHex','PostSet',@check);

      function check(varargin)
        if F.isHex
          warning(['Hexagonal grids are not yet fully supportet for the medianFilter. ' ...
            'It might give reasonable results anyway']);
        end
      end

    end

    function ori = smooth(F,ori,quality)

      ori(quality==0) = nan;

      % this projects to the fundamental region around the mean
      [~,ori] = mean(ori);

      % make verything quaternion
      q = quaternion(ori);

      % some shortcuts
      nn = F.numNeighbours;
      dn = 1+2*nn;

      % the mean distance from every candiate to all others
      meanDist = zeros([size(q),2*F.numNeighbours+1,2*F.numNeighbours+1]);

      % make q a bit larger
      q = [quaternion.nan(F.numNeighbours,size(q,2)+2*F.numNeighbours);...
        [quaternion.nan(size(q,1),F.numNeighbours),...
        q,quaternion.nan(size(q,1),F.numNeighbours)];...
        quaternion.nan(F.numNeighbours,size(q,2)+2*F.numNeighbours)];

      % compute for any candiate the mean distance to all other points
      % the first two loops are for the candidate
      for i1 = 1:dn
        for j1 = 1:dn

          % the candidate
          qq = q(i1+(0:end-dn),j1+(0:end-dn));
          count = zeros(size(qq));

          % compute the distance from the candidate to all other candidates
          for i2 = 1:dn
            for j2 = 1:dn

              omega = angle(qq,q(i2+(0:end-dn),j2+(0:end-dn)));
              [meanDist(:,:,i1,j1),count] = nanplus(meanDist(:,:,i1,j1),omega,count);

            end
          end

          meanDist(:,:,i1,j1) = meanDist(:,:,i1,j1) ./ count;

        end
      end

      % find median
      meanDist = reshape(meanDist,[size(qq),(2*F.numNeighbours+1)^2,]);
      [mm,id] = min(meanDist,[],3);

      [i,j] = ind2sub(size(qq),1:length(qq));
      [ii,jj] = ind2sub([2*F.numNeighbours+1 2*F.numNeighbours+1],id);

      % in regions where everything is nan take simply the center point
      % we may later weaken this to allow inpainting
      ii(isnan(mm)) = nn+1;
      jj(isnan(mm)) = nn+1;

      % compute the final indece to the median
      ind = sub2ind(size(q),i(:)+ii(:)-1,j(:)+jj(:)-1);

      % switch to median
      ori(1:length(ori)) = q(ind);

    end
  end
end

author: mtex-toolbox

repo: mtex

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


text:
classdef DNSGAII < ALGORITHM
 % <multi> <real/integer/label/binary/permutation> <constrained/none> <dynamic>
 % Dynamic NSGA-II
 % type ---   1 --- 1. Mutation based reinitialization 2. Random reinitialization
 % zeta --- 0.2 --- Ratio of reinitialized solutions

%------------------------------- Reference --------------------------------
% K. Deb, U. Bhaskara Rao N., and S. Karthik, Dynamic multi-objective
% optimization and decision-making using modified NSGA-II: A case study on
% hydro-thermal power scheduling, Proceedings of the International
% Conference on Evolutionary Multi-Criterion Optimization, 2007, 803-817.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [type,zeta] = Algorithm.ParameterSet(1,0.2);
            % Reset the number of saved populations (only for dynamic optimization)
            Algorithm.save = sign(Algorithm.save)*inf;

            %% Generate random population
            Population = Problem.Initialization();
            [~,FrontNo,CrowdDis] = EnvironmentalSelection(Population,Problem.N);
            % Archive for storing all populations before each change
            AllPop = [];

            %% Optimization
            while Algorithm.NotTerminated(Population)
                if Changed(Problem,Population)
                    % Save the population before the change
                    AllPop = [AllPop,Population];
                    % React to the change
                    [Population,FrontNo,CrowdDis] = Reinitialization(Problem,Population,type,zeta);
                end
                MatingPool = TournamentSelection(2,Problem.N,FrontNo,-CrowdDis);
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                [Population,FrontNo,CrowdDis] = EnvironmentalSelection([Population,Offspring],Problem.N);
                if Problem.FE >= Problem.maxFE
                    % Return all populations
                    Population = [AllPop,Population];
                    [~,rank]   = sort(Population.adds(zeros(length(Population),1)));
                    Population = Population(rank);
                end
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function p=paramsetch(d,q,m,c,t)
%PARAMSETCH update and check parameter values p=(d,q,m,c,t)
% Usage: (1) function x=func(y,q)
%            d=struct('a',1,'b',2,'c',3); % default parameters
%            p=paramsetch(d,q); % update selected parameters
%
%        (2) function x=func(y,q)
%            d=struct('a',1,'b',2,'c',3); % default parameters
%            c={'p.a>0 && p.a<5','p.b>0'};
%            p=paramsetch(d,q,'E',c); % check parameter ranges
%
%        (3) t={'a','description of parameter a';'c','and of parameter c'}
%            p=paramsetch(d,q,'l',c,t); % list values with optional descritions
%
%  Inputs:
%       d  default parameter structure
%       q  new parameter values either a struct or alternatively matrix with
%          each row a different variable in the same order as the fields of d
%       m  mode string: any combination of the following
%           'a' ignore additional fields in q
%           'A' additional fields in q constitute an error
%           'e' print errors but don't exit
%           'E' print errors and exit
%           'l' list fields and their values (default if no output)
%       c  cell array with parameter checking conditions e.g. 'p.a>3' (use p for structure name)
%       t  cell array with descriptive text for each field in a new row. Either in
%          the form t(:,*)={'field' 'description'} or a single column of
%          descriptions in the same order as the fields of d
%
% Outputs:
%       p  output parameter structure
%

%      Copyright (C) Mike Brookes 2017
%      Version: $Id: paramsetch.m 9302 2017-01-18 16:19:20Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p=d;   % initialize to the default values
numerr=0; % initialize error count
dn=fieldnames(d);
ndn=length(dn);
dup=zeros(length(dn),1); % update flag for existing fields
% sort out input arguments
if nargin<5
    t={'' ''}; % define an empty description array
    if nargin<4
        c=cell(0); % define an empty check condition array
    end
    if nargin<3
        m='';
        if nargin>=2 && numel(q)>0
        end
    end
end
% now update the selected fields
if nargin>1 && numel(q)>0
    if isstruct(q)
        qn=fieldnames(q);
        addnew=~any(lower(m)=='a');
        adderr=any(m=='A');
        for i=1:length(qn)
            fi=qn{i};
            old=isfield(p,fi);
            if addnew || old
                p.(fi)=q.(fi);
            end
            if old
                dup(find(strcmp(fi,dn),1))=1; % indicate field is updated
            end
            if adderr && ~old
                fprintf(2,'%s is an unknown parameter field\n',fi);
                numerr=numerr+1;
            end
        end

    else
        nq=min(size(q,1),ndn); % number of fields to update
        dup(1:nq)=1; % indicate which fields are updated
        for i=1:nq
            p.(dn{i})=q(i,:);
        end
        if size(q,1)>nq && any(lower(m)=='e')
            fprintf(2,'More than %d parameters specified\n',nq);
            numerr=numerr+1;
        end
    end
end
% Apply parameter checks
if any(lower(m)=='e') && numel(c)>0
    for i=1:numel(c)
        if any(~eval(c{i}))
            numerr=numerr+1;
            fprintf(2,'Parameter check failed: %s\n',c{i});
        end
    end
end
% print out a list of the parameters if requested
if ~nargout || any(m=='l')
    pn=fieldnames(p);
    nf=length(pn);
    st=size(t);
    for i=1:nf
        fi=pn{i};
        vi=p.(fi);
        if i>ndn
            cat='+';
        else
            cat='-'+('*'-'-')*dup(i);
        end
        if st(2)>1
            jti=find(strcmp(fi,t(:,1)),1);
            if ~isempty(jti)
                jti=t{jti,2}; % description string
            end
        elseif i<=st(1)
            jti=t{i,1}; % description string
        else
            jti=[];
        end
        if isnumeric(vi) && length(vi)==numel(vi) && isreal(vi) % can print on one line
            fit=fi;
            if size(vi,1)>1
                fit=[fi ''''];
            end
            fprintf('%3d%c %s =',i,cat,fit);
            fprintf(' %g',vi);
            if isempty(jti)
                fprintf('\n');
            else
                fprintf(' = %s\n',jti);
            end
        else
            fprintf('%3d%c %s =',i,cat,fi);
            if isempty(jti)
                fprintf('\n');
            else
                fprintf(' %s =\n',jti);
            end
            disp(vi);
        end
    end
end
if numerr>0 && any(m=='E')
    error('%d error%c in parameter specification',numerr,(' '+(numerr>1)*('s'-' ')));
end

author: singaxiong

repo: SignalGraph


text:
function [T]=writeNodeFile_tetGen(inputStruct)

%%

dispStartTitleGibbonCode('Writing NODE file');

%% PARSE INPUT STRUCTURE

V=inputStruct.Nodes;
nodeFileName=inputStruct.modelName;

%Force extension to be .node
[pathstr,name,~] = fileparts(nodeFileName);
nodeFileName=fullfile(pathstr,[name,'.node']);

%%

V_id=1:1:size(V,1);
V_field=[V_id(:) V];
V_char=sprintf('%i %0.16e %0.16e %0.16e \n',V_field');
V_cell = regexp(V_char, '\n', 'split');

if numel(V_cell)>1
    V_cell=V_cell(1:end-1);
end

T={'#PART 1 - Node list';'#num nodes, num dimensions, num attributes, num boundary markers'};
numNodes=size(V,1);
numDims=3;
numAtr=0;
boundMarker=0;
doubleList=[numNodes numDims numAtr boundMarker];
charList=sprintf('%i %i %i %i',doubleList');
T(end+1)={charList};
T(end+1)={'#Node ID, x, y, z,attribute,boundary marker'};
T(end+1:end+numel(V_cell))=V_cell;

%% SAVING TXT FILE

cell2txtfile(nodeFileName,T,0,0);
dispDoneGibbonCode;

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON