text:
%ELLIPSES_DEMO  One-line description here, please.
%
%   output = ellipses_demo(input)
%
%   Example
%   ellipses_demo
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2019-08-09,    using Matlab 9.6.0.1072779 (R2019a)
% Copyright 2019 INRA - Cepia Software Platform.

load fisherIris;
figure; hold on; set(gca, 'fontsize', 14);
colors = {'b', 'g', 'm'};
hi = zeros(1,3);
for i = 1:3
  pts = meas((1:50)+(i-1)*50, 3:4);
  hi(i) = drawPoint(pts, 'marker', 'o', 'color', colors{i}, 'markerfacecolor', colors{i});
  drawEllipse(equivalentEllipse(pts), 'color', colors{i}, 'linewidth', 2);
end
legend(hi, species([1 51 101]), 'Location', 'NorthWest');

print(gcf, 'equivalentEllipse_demo.png');

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [varargout, peval] = separcomp(dpix, peval, winit_pix, hinit)
% separcomp(dpix, p_script, savethis, winit, hinit)
% separate components
% V ~ WH
% V -> N_pix x N_t
% W -> N_pix x N_comp - xth pixel of the ith components
% H -> N_copm x N_t - contribution of the i-th component in the time t

[peval.nx, peval.ny, peval.nt] = size(dpix);
peval.numpix = peval.nx*peval.ny;
peval.meandata = mean(dpix(:));

dvec = reshape(dpix,peval.numpix, peval.nt);
dpix_dip = dip_image(dpix);

% background subtraction with empirical values:
peval.bg_clip = 'no';
peval.bg_fs_var = 5;
peval.bg_perc = 20;
peval.bg_ob_dist = 8;
if ~isfield(peval, 'bg')
    [out_nobg, peval.bg, bg_im]=backgroundoffset(dpix_dip, peval.bg_clip, peval.bg_fs_var, peval.bg_perc, peval.bg_ob_dist);
end

if ~isfield(peval, 'ncomp')
    peval.ncomp=estimate_ncpomp(dvec);
end

[winit, hinit] = initwh(winit_pix, hinit, peval); %initialization of w and h

if ~isfield(peval, 'w_fixvec') peval.w_fixvec=[]; end
if ~isfield(peval, 'h_fixvec') peval.h_fixvec=[]; end

if isempty(peval.w_fixvec)
    peval.w_fixvec = peval.ncomp + 1; %fixing background component
elseif ~(peval.w_fixvec(end) == peval.ncomp + 1)
    peval.w_fixvec = [peval.w_fixvec, peval.ncomp + 1];
end

if isempty(peval.h_fixvec)
    peval.h_fixvec = peval.ncomp + 1; %fixing background component
elseif ~(peval.h_fixvec(end) == peval.ncomp + 1)
    peval.h_fixvec = [peval.h_fixvec, peval.ncomp + 1];
end

verbose = 1;
% if strcmp (peval.method, 'nmf_classic') %classical nmf updates
%     [w,h,peval, dtrace, htrace]=nmf_classic(dvec,winit,hinit,peval,verbose);
% elseif strcmp (peval.method, 'nmf_conjgrad_penalty') %classical nmf updates
%     [w,h,peval, dtrace, htrace]=nmf_conjgrad_penalty(dvec,winit,hinit,peval,verbose);
% end
eval (['[w,h,peval, dtrace, htrace]=' peval.method '(dvec,winit,hinit,peval,verbose);']);

varargout = struct('w',w,'h',h, 'dtrace', dtrace, 'htrace', htrace);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


% Take in a vector of complex samples representing individual QPSK constellation points with the constellation rotated
% such that the points are ideally at 1+i, -1+i, -1-i, and 1-i.
%
% This function uses hard decision, so it's not what you want to use in low SNR environments
%
% The constellation mapping is:
%     1+i == 0b00
%     1-i == 0b01
%    -1+i == 0b10
%    -1-i == 0b11
%
% Which comes from https://github.com/ttsou/openphy/blob/master/src/lte/qam.c#L35
%
% @param data_carriers Row or column vector of complex samples
% @return quantized_bits Vector of 1/0 values that make up the bits demapped from the provided sample vector
function [quantized_bits] = quantize_qpsk(data_carriers)
    assert(iscolumn(data_carriers) || isrow(data_carriers), "Data carriers must be row/column vector");

    quantized_bits = zeros(length(data_carriers), 1);

    % Track where in the `quantized_bits` vector the new bits should be placed
    bits_offset = 1;

    % Walk through each complex sample in the input vector
    for sample_idx = 1:length(data_carriers)
        sample = data_carriers(sample_idx);

        % Determine bit mapping based on the quadrant that the sample is located
        if (real(sample) > 0 && imag(sample) > 0)
            bits = [0, 0];
        elseif (real(sample) > 0 && imag(sample) < 0)
            bits = [0, 1];
        elseif (real(sample) < 0 && imag(sample) > 0)
            bits = [1, 0];
        elseif (real(sample) < 0 && imag(sample) < 0)
            bits = [1, 1];
        else
            bits = [0, 0];
        end

        % Save off the quatized bits and move the counter ahead by 2
        quantized_bits(bits_offset:bits_offset+1) = bits;
        bits_offset = bits_offset + 2;
    end
end


text:
% test for tensor voting
% Use the tensor structure direction field as an input.

sigma = 0.5;
mu = 3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a test image
disp('--> Building test image.');
Radius = 0.9;

n = 75;
x = -1:2/(n-1):1;
y = linspace(0,1,n);
[X,Y] = meshgrid(x,x);

M = double( abs( X.^2+Y.^2-Radius^2 )<0.03 );
M = double( X.^2+Y.^2<Radius^2 );

% add some noise in the middle of the image
sigma = 0.8;    % level of the noise
r = 0.05;        % width of the gap
I = find(abs(X-0.5)<r);
M = M + sigma*randn(n).*(abs(X-0.5)<r);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% structure tensor computation
disp('--> Compute tensor structure.');
h = compute_gaussian_filter([11 11],0.6/n,[n n]);
Mh = perform_convolution(M,h);
ST = compute_rigidity_tensor(Mh);
h = compute_gaussian_filter([11 11],0.3/n,[n n]);
for i=1:4   % convole each entry
    ST(:,:,i) = perform_convolution( ST(:,:,i) ,h);
end
[e1,e2,l1,l2] = perform_tensor_decomp(ST);   % e1 is a regularized gradient

% As the main direction is the direction orthogonal to the discontinuity
% we must swap the eigenvectors.
T = perform_tensor_recomp(e2,e1,l1,l2);

clf;
subplot(2,2,1);
plot_vf(e2,M);
axis off;
colormap gray(256);
title('Original Image and Field');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% an example of field
A = compute_stick_tf([1,0.3],n);
[e1,e2,l1,l2] = perform_tensor_decomp(A);

subplot(2,2,2);
plot_vf(e1,l1);
axis off;
title('Example of Stick Field');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% peform voting
disp('--> Performing tensor voting.');
sigma = 0.2;
c = 0.1*sigma;
T = perform_voting(T,sigma,c);
disp('--> Performing tensor decomposition.');
[e1,e2,l1,l2] = perform_tensor_decomp(T,'abs');

% re-orient vf
e1 = perform_vf_reorientation(e1);
e2 = perform_vf_reorientation(e2);

% plot e1
subplot(2,2,3);
plot_vf(e1,l1);
axis off;
title('New Field with Eigenvalue');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% curve extraction via zeros crossing of A=ds/dn where s=l1-l2, n=e2,
% orthogonal to e1
s = l1-l2;
g = compute_grad(s);
A = prod_vf_vf(g,e2);

% extract only the 20 longest curves
options.min_length = 20;
c_list = perform_curve_extraction( A,0,10 );

% plot curves
subplot(2,2,4);
plot_curve(c_list,A);
axis off;
title('Extracted Curves');
saveas(gcf, 'test_structure', 'png');

author: gpeyre

repo: matlab-toolboxes


text:
function plotBoundingBox(img, bb, score, colors)
% Visualices and image and a set of bounding boxes
%
%   img = image
%   bb = each row is [xmin ymin xmax ymax]

Nbb = size(bb,1);
if nargin==2
    score = 3*ones(Nbb,1);
end

if ~isempty(img)
    image(img); axis('off'); axis('equal')
    if size(img,3)==1; colormap(gray(256)); end
    hold on
end

if nargin<4
    colors = 'rrgbcmy';
end
for i = 1:Nbb
    plot([bb(i,1) bb(i,3) bb(i,3) bb(i,1) bb(i,1)], [bb(i,2) bb(i,2) bb(i,4) bb(i,4) bb(i,2)], 'k', 'linewidth', score(i))
    plot([bb(i,1) bb(i,3) bb(i,3) bb(i,1) bb(i,1)], [bb(i,2) bb(i,2) bb(i,4) bb(i,4) bb(i,2)], colors(mod(i,6)+1), 'linewidth', score(i)-2)
end

%title(sprintf('There are %d boxes', Nbb))

author: CSAILVision

repo: LabelMeToolbox


text:
% BONMIN  Solve a MINLP using BONMIN
%
% THIS IS A LOW LEVEL FUNCTION - USE opti_bonmin() INSTEAD!
%
% bonmin uses the Basic Open source Nonlinear Mixed Integer library.
%
%   [x,output] = bonmin(x0,funcs,opts)
%
%   Input arguments:
%       x0 - initial solution guess
%       funcs - functions structure (see below)
%       opts - options structure (see below)
%
%   Return arguments:
%       x - solution vector
%       output - return information structure
%
%   Funcs Fields:
%       objective - objective function handle
%       gradient - objective gradient function handle
%       constraints - constraints function handle (optional)
%       jacobian - constraints Jacobian function handle (sparse)
%       jacobianstructure - constraints Jacobian sparsity structure (sparse)
%       hessian - objective Hessian of the Lagrangian (sparse,optional)
%       hessianstructure - objective Hessian sparsity structure (sparse)
%
%   Opts Fields:
%       lb - decision variable lower bounds
%       ub - decision variable upper bounds
%       cl - constraints lower bounds
%       cu - constraints upper bounds
%       ipopt - relaxed solver settings structure (see ipoptset.m)
%       bonmin - bonmin settings (see bonminset.m)
%       var_type - integer variable indicies (-1 binary, 0 cont, 1 integer)
%       var_lin - decision variable linearity (0 nonlinear, 1 linear)
%       cons_lin - constraint linearity (0 nonlinear, 1 linear)
%       display - log level (0 off, 1-5 increasing)
%
%   Return Status:
%       0 - proven optimal
%      -1 - exceeded iterations
%       1 - proven infeasible
%       2 - feasible
%       3 - unbounded or infeasible
%       4 - no solution known
%       other - bonmin error / other
%
%
%   Based almost entirely on the MEX interface to IPOPT by Dr. Peter Carbonetto.
%
%   Copyright (C) 2011 Jonathan Currie (IPL)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function value = p13_exact ( dim_num )

%*****************************************************************************80
%
%% P13_EXACT returns the exact integral for problem 13.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Output, real VALUE, the exact value of the integral.
%
  value = 0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function  makeVideoFrame(Fig,filename,FigOpt,ExpOpt)

% MAKEVIDEOFRAME  Export figure to image.
%   MAKEVIDEOFRAME(FIG,FNAME,FIGOPT,EXPOPT) exports the constents of figure
%   FIG to an image with name filename, using FIGOPT and EXPOPT as options.

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

dir = [ExpOpt.root 'figures/' ExpOpt.site '/' ExpOpt.lmkTypes '/images/'];
if ~isdir(dir)
    mkdir(dir);
end
filepath = [dir filename];
imgFrame(Fig.fig,filepath);

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

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

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function points = intersectLineCylinder(line, cylinder, varargin)
%INTERSECTLINECYLINDER Compute intersection points between a line and a cylinder.
%
%   POINTS = intersectLineCylinder(LINE, CYLINDER)
%   Returns intersection points between a line and a cylinder.
%
%   Input parameters:
%   LINE     = [x0 y0 z0  dx dy dz]
%   CYLINDER = [x1 y1 z1 x2 y2 z2 R]
%
%   Output:
%   POINTS   = [x1 y1 z1 ; x2 y2 z2]
%
%   POINTS = intersectLineCylinder(LINE, CYLINDER, 'checkBounds', B)
%   Where B is a boolean (TRUE by default), check if the points are within
%   the bounds defined by the two extreme points. If B is false, the
%   cylinder is considered to be infinite.
%
%   Example
%     % Compute intersection between simple vertical cylinder and line
%     line = [60 60 60 1 2 3];
%     cylinder = [20 50 50 80 50 50 30];
%     points = intersectLineCylinder(line, cylinder);
%     % Display the different shapes
%     figure;
%     drawCylinder(cylinder);
%     hold on; light;
%     axis([0 100 0 100 0 100]);
%     drawLine3d(line);
%     drawPoint3d(points, 'ko');
%
%
%     % Compute intersections when one of the points is outside the
%     % cylinder
%     line = [80 60 60 1 2 3];
%     cylinder = [20 50 50 80 50 50 30];
%     intersectLineCylinder(line, cylinder)
%     ans =
%           67.8690   35.7380   23.6069
%
%
%   See also
%   lines3d, intersectLinePlane, drawCylinder, cylinderSurfaceArea
%
%   References
%   See the link:
%   http://www.gamedev.net/community/forums/topic.asp?topic_id=467789
%

% ------
% Author: David Legland, from a file written by Daniel Trauth (RWTH)
% E-mail: david.legland@inra.fr
% Created: 2007-01-27
% Copyright 2007-2022

%% Parse input arguments

% default arguments
checkBounds = true;
% type of cylinder, one of {'closed', 'open', 'infinite'}
type = 'closed';

% parse inputs
while length(varargin)>1
    var = varargin{1};
    if strcmpi(var, 'checkbounds')
        checkBounds = varargin{2};
    elseif strcmpi(var, 'type')
        type = varargin{2};
    else
        error(['Unkown argument: ' var]);
    end
    varargin(1:2) = [];
end

%% Parse cylinder parameters

% Starting point of the line
l0 = line(1:3);

% Direction vector of the line
dl = line(4:6);

% position of cylinder extremities
c1 = cylinder(1:3);
c2 = cylinder(4:6);

% Direction vector of the cylinder
dc = c2 - c1;

% Radius of the cylinder
r = cylinder(7);

%% Resolution of a quadratic equation to find the increment

% normalisation coefficient corresponding to direction of vector
coef = dc / dot(dc, dc);

% Substitution of parameters
e = dl - dot(dl,dc) * coef;
f = (l0-c1) - dot(l0-c1, dc) * coef;

% Coefficients of 2-nd order equation
A = dot(e, e);
B = 2 * dot(e,f);
C = dot(f,f) - r^2;

% compute discriminant
delta = B^2 - 4*A*C;

% check existence of solution(s)
if delta < 0
    points = zeros(0, 3);
    return;
end

% extract roots
pos1 = (-B + sqrt(delta)) / (2*A);
pos2 = (-B - sqrt(delta)) / (2*A);
posList = [pos1;pos2];

%% Estimation of point positions

% process the smallest position
pos1 = min(posList);

% Point on the line: l0 + x*dl = p
point1 = l0 + pos1 * dl;

% process the greatest position
pos2 = max(posList);

% Point on the line: l0 + x*dl = p
point2 = l0 + pos2 * dl;

% Format result
points = [point1 ; point2];

%% Check if points are located between bounds

% if checkBounds option is not set, we can simply skip the rest
if ~checkBounds || strncmpi(type, 'infinite', 1)
    return;
end

% compute cylinder axis
axis = [c1 dc];

% compute position on axis
ts = linePosition3d(points, axis);

% check bounds for open cylinder
% (keep only intersection points whose projection is between the two
% cylinder extremities)
if strncmpi(type, 'open', 1)
    ind = ts>=0 & ts<=1;
    points = points(ind, :);
    return;
end

% which intersection fall before and after bounds
ind1 = find(ts < 0);
ind2 = find(ts > 1);

% case of both intersection on the same side -> no intersection
if length(ind1) == 2 || length(ind2) == 2
    points = zeros(0, 3);
    return;
end

% Process the remaining case of closed cylinder
% -> compute eventual intersection(s) with end faces
if ~isempty(ind1)
    plane = createPlane(c1, dc);
    points(ind1, :) = intersectLinePlane(line, plane);
end
if ~isempty(ind2)
    plane = createPlane(c2, dc);
    points(ind2, :) = intersectLinePlane(line, plane);
end


text:
function [same,sloc] = setset2(iset,jset)
%SETSET2 a (much) faster variant of ISMEMBER for edge lists.
%   [IN] = SETSET2(ISET,JSET) returns an I-by-1 array IN,
%   with IN(K) = TRUE if ISET(K,:) is present in JSET. This
%   routine is essentially an optimised ISMEMBER variant de-
%   signed for processing lists of edge indexing. ISET is an
%   I-by-2 array of "query" edges, JSET is a J-by-2 array of
%   edges to test against.
%
%   See also ISMEMBER

%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 29/01/2017

%---------------------------------------------- basic checks
    if ( ~isnumeric(iset) || ~isnumeric(jset) )
        error('setset2:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end

%---------------------------------------------- basic checks
    if (ndims(iset) ~= +2 || ndims(jset) ~= +2)
        error('setset2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end
    if (size(iset,2)~= +2 || size(jset,2)~= +2)
        error('setset2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

%---------------------------------------------- set v1 <= v2
    iset = sort(iset,2) ;
    jset = sort(jset,2) ;

%-- this is the slow, but easy-to-undertsand version of what
%-- is happening here...

  % if (nargout == +1)
  % same = ismember(iset,jset,'rows') ;
  % else
  % [same,sloc] = ...
  %        ismember(iset,jset,'rows') ;

%-- as above, the 'ROWS' based call to ISMEMBER can be sped
%-- up by casting the edge lists (i.e. pairs of UINT32 valu-
%-- es) to DOUBLE, and performing the sorted queries on vec-
%-- tor inputs!
    if (nargout == +1)
    same       = ismember( ...
        iset*[2^31;1], jset*[2^31;1]) ;
    else
   [same,sloc] = ismember( ...
        iset*[2^31;1], jset*[2^31;1]) ;
    end

end

author: CHLNDDEV


text:
function [theta] = tapas_sem_sample_gaussian_uniform_priors(ptheta)
%% Sample from a Gaussian prior.
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

np = size(ptheta.jm, 2);

sample_pars = logical(sum(ptheta.jm, 2));

if size(ptheta.pm, 2) ~= np
    pm = diag(ptheta.pm);
else
    pm = ptheta.pm;
end

lt = chol(pm);
theta = ptheta.mu +  lt \ ptheta.jm * randn(np, 1);

% Beta parameters.
bdist = zeros(size(ptheta.pm));
bdist(ptheta.bdist) = 1;
bdist = logical(bdist .* sample_pars);
vals = betarnd(abs(ptheta.mu), abs(ptheta.pm));
vals = ptheta.jm * ptheta.sm' * vals;
theta(bdist) = log(vals(bdist) ./ ( 1 - vals(bdist)));

theta(~sample_pars) = ptheta.p0(~sample_pars);

end

author: translationalneuromodeling

repo: tapas


text:
%% Loading file paths
clear
close all
addpath('src/utility/')
addpath('src/')
addpath('../mesh2tri/')

[file,path] = uigetfile({'*.csv'},...
    'Please select the .csv file... containing SDF of an object', ...
    '../', 'MultiSelect','off');
if isequal(file, 0)
    error('No file selected.');
end

file_path = [path, file];
sdf = csvread(file_path)';
voxelGrid.size = ones(1, 3) * sdf(1);
voxelGrid.range = sdf(2 : 7);
sdf = sdf(8 : end);

voxelGrid.x = linspace(voxelGrid.range(1), voxelGrid.range(2), voxelGrid.size(1));
voxelGrid.y = linspace(voxelGrid.range(3), voxelGrid.range(4), voxelGrid.size(2));
voxelGrid.z = linspace(voxelGrid.range(5), voxelGrid.range(6), voxelGrid.size(3));
[x, y, z] = ndgrid(voxelGrid.x, voxelGrid.y, voxelGrid.z);
voxelGrid.points = reshape(cat(4, x, y, z), [], 3)';
voxelGrid.interval = (voxelGrid.range(2) - voxelGrid.range(1)) / (voxelGrid.size(1) -1);
voxelGrid.truncation = 1.2 * voxelGrid.interval; %1.2
voxelGrid.disp_range = [-inf, voxelGrid.truncation];
voxelGrid.visualizeArclength = 0.01 * sqrt(voxelGrid.range(2) - voxelGrid.range(1));
clearvars x y z

sdf = min(max(sdf, -voxelGrid.truncation), voxelGrid.truncation);

%% marching-primitives
tic
[x] = MPS(sdf, voxelGrid, 'minArea', 5);
toc

%% triangularization and compression
[mesh_original] = meshSuperquadrics(x, 'Arclength', voxelGrid.visualizeArclength);
% compression
mesh = reducepatch(mesh_original.f, mesh_original.v, 0.1); %0.05-2
stl = triangulation(mesh.faces, mesh.vertices);

ifsave = true;
[pathname,name,ext] = fileparts([path, file]);
if ifsave
    x_save = single(x);
    save(fullfile(pathname, [name,'_sq.mat']),'x_save')
    stlwrite(stl,fullfile(pathname, [name,'_sq.stl']), 'binary')
end

%% visualize

close all
view_vector = [151, -40];
light_vector = [190,10];
camera_roll = 50;
color = [145,163,176] ./255;

% rearrange sdf to 3D array for region connection checking
sdf3d_region = reshape(sdf, voxelGrid.size(1), voxelGrid.size(2), voxelGrid.size(3));
% rearrange sdf to 3D array for visualization
sdf3d = permute(sdf3d_region, [2, 1, 3]);

[mesh_gt.f, mesh_gt.v] = plyread(...
    fullfile(pathname, [name,'_watertight.ply']),'tri');

figure(1)
trisurf(mesh_gt.f,mesh_gt.v(:,1),mesh_gt.v(:,2),mesh_gt.v(:,3), ...
    'FaceColor', color, 'FaceAlpha', 1, 'EdgeColor', 'none')
axis equal
view(view_vector)
camroll(camera_roll)
light
lightangle(light_vector(1), light_vector(2))
material dull
axis(voxelGrid.range)
grid off
axis off
title('ground truth mesh from marching cubes')

figure(2)
trisurf(mesh.faces,mesh.vertices(:,1),mesh.vertices(:,2),mesh.vertices(:,3), ...
    'FaceColor', color, 'FaceAlpha', 1, 'EdgeColor', 'none')
axis equal
view(view_vector)
camroll(camera_roll)
light
lightangle(light_vector(1), light_vector(2))
material dull
axis(voxelGrid.range)
grid off
axis off
title('superquadrics representation from marching primitives')

figure(3)
trisurf(mesh_gt.f,mesh_gt.v(:,1),mesh_gt.v(:,2),mesh_gt.v(:,3), ...
    'FaceColor', 'g', 'FaceAlpha', 0.5, 'EdgeColor', 'none')
hold on
trisurf(mesh.faces,mesh.vertices(:,1),mesh.vertices(:,2),mesh.vertices(:,3), ...
    'FaceColor', color, 'FaceAlpha', 1, 'EdgeColor', 'none')
view(view_vector)
camroll(camera_roll)
axis equal
light
lightangle(light_vector(1), light_vector(2))
material dull
axis(voxelGrid.range)
grid off
axis off
hold off
title('overlapping recovered representation with the ground truth')

author: ChirikjianLab

repo: Marching-Primitives


text:
function response = Forgetful(n,p)

persistent ncount trustcount roundstore
if n==1; ncount=0; trustcount=0; roundstore=0; end
if ncount==1 && (p==1 || p==5); trustcount=trustcount+1; end

if n>5
    if isequal( roundstore(end-3:end), [5 0 5 0] );
        trustcount = 2;
        ncount = 5;
    end
end

if ( p==1 || p==5 || n==20 ) && ncount<5 && n~=1
    response = 'defect';
else
    response = 'cooperate';
    if trustcount>2; response = 'defect'; end
    ncount = 0;
end

roundstore(n) = p;

ncount = ncount + 1;

author: Sable

repo: mcbench-benchmarks


text:
function varargout = size(F, dim)
%SIZE   Size of a SPHEREFUNV object
%   D = SIZE(F) returns a three-element vector D = [K, M, N]. If F is a column
%   SPHEREFUNV object then K is the number of components in F, N and M are INF.
%   If F is a row vector then K and M are INF and N is the number of components
%   of F.
%
%   [K, M, N] = SIZE(F) returns the dimensions of F as separate output
%   variables.
%
%   D = SIZE(F, DIM) returns the dimensions specified by the dimension DIM.
%

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

% Empty check:
if ( isempty(F) )
    varargout = {[], []};
    return
end

nF = 3;
tr = F.isTransposed;

% Manually work out the size:
if ( ~tr )
    K = nF;
    M = inf;
    N = inf;
else
    N = nF;
    K = inf;
    M = inf;
end

% Default to vector:
if ( nargin == 1 )
    dim = 0;
end

% Manually work out what should be displayed:
if ( dim == 1 )
    varargout = { K };
elseif ( dim == 2 )
    varargout = { M };
elseif ( dim == 3 )
    varargout = { N };
elseif ( ( dim == 0 ) && ( nargin == 1 ) )
    if ( nargout <= 1 )
        varargout = { [K, M, N] };
    else
        varargout = { K, M, N };
    end
else
    error('SPHEREFUN:SPHEREFUNV:size:dim', 'Unrecognised dimension.');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% OPTI MEX Install Script

% This file enables me to quickly rebuild all solvers - it is not intended
% to be called by the user!

% My build platform:
% - Windows 7 x64
% - Visual Studio 2015
% - Intel Compiler XE (C++ & FORTRAN)
% - Intel Math Kernel Library

% YOU MUST BE IN THE BASE DIRECTORY OF THE OPTI TOOLBOX!

addpath([cd '/Solvers/Source'])
addpath([cd '/Utilities/Source'])

% %% -- Intel C++ Solvers --
% clc
% fprintf(2,'Please Specify Intel C++ as your Compiler with VS2015 Linker (Option 3, Don''t Look for Installed Compilers)...\n\n');
% mex -setup
% %% CLP
% opti_CLP_Install

%% -- VS2017 Solvers --
clc
fprintf(2,'Please Specify Visual Studio 2017 as your Compiler...\n\n');
mex -setup
%% CLP
opti_CLP_Install

%% BONMIN
opti_BONMIN_Install

%% CBC
opti_CBC_Install

%% CSDP
opti_CSDP_Install

%% DSDP
opti_DSDP_Install

%% FILTERSD
opti_FILTERSD_Install

%% GLPK
opti_GLPK_Install

%% GSL
opti_GSL_Install

%% IPOPT
opti_IPOPT_Install

%% L-BFGS-B
opti_LBFGSB_Install

%% LEVMAR
opti_LEVMAR_Install

%% LP_SOLVE
opti_LPSOLVE_Install

%% M1QN3
opti_M1QN3_Install

%% MINPACK
opti_MINPACK_Install

%% MKL TR NLS
opti_MKLTRNLS_Install

%% MUMPS
opti_MUMPS_Install

%% NL2SOL
opti_NL2SOL_Install

%% NLOPT
opti_NLOPT_Install

%% NOMAD
opti_NOMAD_Install

%% OOQP
opti_OOQP_Install

%% PSWARM
opti_PSWARM_Install

%% SCIP
opti_SCIP_Install

%% -- Utilities --
%% MKL JAC
opti_MKLJAC_Install

%% CoinUtils
opti_COINUTILS_Install

%% AMPL
opti_ASL_Install

%% RMATHLIB
opti_RMathlib_Install

%% Remove Paths
rmpath([cd '/Solvers/Source'])
rmpath([cd '/Utilities/Source'])


text:
function pot = rescale_pot(pot, s)
% RESCALE_POT Add a constant to the mpot scale factor.
% pot = rescale_pot(pot, s)

pot.p = pot.p*s;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function X = matrandcong(m,n,gamma)
%MATRANDCONG Create a random matrix with a fixed congruence.
%
%   X = MATRANDCONG(M,N,GAMMA) creates a matrix X of size M x N such
%   that each column of X has norm 1 and any two columns of X have an inner
%   product equal to GAMMA.
%
%   Based on code from Evrim Acar and the paper G. Tomasi and R. Bro, A
%   comparison of algorithms for fitting the PARAFAC model, Computational
%   Statistics & Data Analysis, 50: 1700-1734, 2006.
%
%   See also MATRANDORTH, MATRANDNORM, CREATE_PROBLEM, CREATE_GUESS.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

CG = gamma * ones(n,n) + (1-gamma) * eye(n);
CGR = chol(CG);
X = randn(m,n);
[Q,~] = qr(X,0);
X = Q * CGR;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
clear all;

% load previously generated fig file
figFile = 'multiple.fig';
open(figFile)

% change properties
opt.XLabel = 'Time, t (ms)'; % xlabel
opt.YLabel = 'Voltage, V (V)'; %ylabel
opt.YTick = [-10, 0, 10]; %[tick1, tick2, .. ]
opt.XLim = [0, 80]; % [min, max]
opt.YLim = [-11, 11]; % [min, max]

opt.Colors = [ % three colors for three data set
    1,      0,       0;
    0.25,   0.25,    0.25;
    0,      0,       1;
    ];

opt.LineWidth = [2, 2, 2]; % three line widths
opt.LineStyle = {'-', '-', '-'}; % three line styles
opt.Markers = {'o', '', 's'};
opt.MarkerSpacing = [15, 15, 15];
opt.Legend = {'\theta = 0^o', '\theta = 45^o', '\theta = 90^o'}; % legends

% Save? comment the following line if you do not want to save
opt.FileName = 'plotMarkers.png';

% create the plot
setPlotProp(opt);

author: masumhabib

repo: PlotPub

save_path: github-repos/MATLAB/masumhabib-PlotPub


text:
function [p, v, q] = ch_nav_equ_local_tan(p, v, q ,acc, gyr, dt, gN)
%  惯导解算更新，当地直角坐标系，不考虑地球自转
% p          位置 XYZ 单位 m
% v          速度 XYZ 单位 m/s
% q         Qb2n姿态,四元数表示
% acc      比力， 加速度计测量值 单位  (m/s^2),
% gyr      角速度 (rad/s)]
% dt        dt (s) 积分间隔如 0.01s
% gn       当地重力向量

old_v = v;

sf = acc;

%  姿态结算
q = ch_att_upt(q, gyr, dt);

% 速度解算
sf = ch_qmulv(q, sf);
sf = sf + gN;
v = old_v + dt *sf;

% 位置解算
p = p + (old_v + v) *dt/2;

end

%
%
% function x = ch_nav_equ_local_tan(x ,u, dt, gN)
%
% persistent a_old;
% if isempty(a_old)
%    a_old= u(1:3);
% end
%
% old_v = x(4:6);
%
% a_new =u(1:3);
% %sf = sf + 0.5*cross(u(4:6)*dt, sf);
%
% %  姿态结算
% gyr = u(4:6);
% q_old = x(7:10);
% x(7:10) = ch_att_upt(x(7:10), gyr, dt);
% q_new = x(7:10);
%
% % 速度解算
%
% x(4:6) = old_v + ((ch_qmulv(q_new, a_new) + ch_qmulv(q_old, a_old) )/2 + gN) *dt;
%
% % 位置解算
% x(1:3) = x(1:3) + (old_v + x(4:6)) *dt/2;
% a_old = a_new;
% end
%
%

author: yandld


text:
function B=unimodalcrossproducts(XtX,XtY,Bold)

%UNIMODALCROSSPRODUCTS
% Solves the problem min|Y-XB'| subject to the columns of
% B are unimodal and nonnegative. The algorithm is iterative and
% only one iteration is given, hence the solution is only improving
% the current estimate
%
% I/O B=unimodalcrossproducts(XtX,XtY,Bold)
% Modified from unimodal.m to handle crossproducts in input 1999
% Reference
% Bro and Sidiropoulos, "Journal of Chemometrics", 1998, 12, 223-247.

% Copyright (C) 1995-2006  Rasmus Bro & Claus Andersson
% Copenhagen University, DK-1958 Frederiksberg, Denmark, rb@life.ku.dk
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA  02110-1301, USA.

B=Bold;
F=size(B,2);
for f=1:F
   xty = XtY(f,:)-XtX(f,[1:f-1 f+1:F])*B(:,[1:f-1 f+1:F])';
   beta=pinv(XtX(f,f))*xty;
   B(:,f)=ulsr(beta',1);
end

function [b,All,MaxML]=ulsr(x,NonNeg);

% ------INPUT------
%
% x          is the vector to be approximated
% NonNeg     If NonNeg is one, nonnegativity is imposed
%
%
%
% ------OUTPUT-----
%
% b 	     is the best ULSR vector
% All 	     is containing in its i'th column the ULSRFIX solution for mode
% 	     location at the i'th element. The ULSR solution given in All
%            is found disregarding the i'th element and hence NOT optimal
% MaxML      is the optimal (leftmost) mode location (i.e. position of maximum)
%
% ___________________________________________________________
%
%
%               Copyright 1997
%
% Nikos Sidiroupolos
% University of Maryland
% Maryland, US
%
%       &
%
% Rasmus Bro
% Royal Veterinary & Agricultural University
% Denmark
%
%
% ___________________________________________________________

% This file uses MONREG.M

x=x(:);
I=length(x);
xmin=min(x);
if xmin<0
  x=x-xmin;
end

% THE SUBSEQUENT
% CALCULATES BEST BY TWO MONOTONIC REGRESSIONS

% B1(1:i,i) contains the monontonic increasing regr. on x(1:i)
[b1,out,B1]=monreg(x);

% BI is the opposite of B1. Hence BI(i:I,i) holds the monotonic
% decreasing regression on x(i:I)
[bI,out,BI]=monreg(flipud(x));
BI=flipud(fliplr(BI));

% Together B1 and BI can be concatenated to give the solution to
% problem ULSR for any modloc position AS long as we do not pay
% attention to the element of x at this position

All=zeros(I,I+2);
All(1:I,3:I+2)=B1;
All(1:I,1:I)=All(1:I,1:I)+BI;
All=All(:,2:I+1);
Allmin=All;
Allmax=All;
% All(:,i) holds the ULSR solution for modloc = i, disregarding x(i),

iii=find(x>=max(All)');
b=All(:,iii(1));
b(iii(1))=x(iii(1));
Bestfit=sum((b-x).^2);
MaxML=iii(1);
for ii=2:length(iii)
  this=All(:,iii(ii));
  this(iii(ii))=x(iii(ii));
  thisfit=sum((this-x).^2);
  if thisfit<Bestfit
    b=this;
    Bestfit=thisfit;
    MaxML=iii(ii);
  end
end

if xmin<0
  b=b+xmin;
end

% Impose nonnegativity
if NonNeg==1
  if any(b<0)
    id=find(b<0);
    % Note that changing the negative values to zero does not affect the
    % solution with respect to nonnegative parameters and position of the
    % maximum.
    b(id)=zeros(size(id))+0;
  end
end

function [b,B,AllBs]=monreg(x);

% Monotonic regression according
% to J. B. Kruskal 64
%
% b     = min|x-b| subject to monotonic increase
% B     = b, but condensed
% AllBs = All monotonic regressions, i.e. AllBs(1:i,i) is the
%         monotonic regression of x(1:i)
%
%
% Copyright 1997
%
% Rasmus Bro
% Royal Veterinary & Agricultural University
% Denmark
% rb@kvl.dk
%

I=length(x);
if size(x,2)==2
   B=x;
else
   B=[x(:) ones(I,1)];
end

   AllBs=zeros(I,I);
   AllBs(1,1)=x(1);
   i=1;
   while i<size(B,1)
      if B(i,1)>B(min(I,i+1),1)
          summ=B(i,2)+B(i+1,2);
          B=[B(1:i-1,:);[(B(i,1)*B(i,2)+B(i+1,1)*B(i+1,2))/(summ) summ];B(i+2:size(B,1),:)];
          OK=1;
          while OK
             if B(i,1)<B(max(1,i-1),1)
                summ=B(i,2)+B(i-1,2);
                B=[B(1:i-2,:);[(B(i,1)*B(i,2)+B(i-1,1)*B(i-1,2))/(summ) summ];B(i+1:size(B,1),:)];
                i=max(1,i-1);
             else
                OK=0;
             end
          end
          bInterim=[];
          for i2=1:i
             bInterim=[bInterim;zeros(B(i2,2),1)+B(i2,1)];
          end
          No=sum(B(1:i,2));
          AllBs(1:No,No)=bInterim;
      else
          i=i+1;
          bInterim=[];
          for i2=1:i
             bInterim=[bInterim;zeros(B(i2,2),1)+B(i2,1)];
          end
          No=sum(B(1:i,2));
          AllBs(1:No,No)=bInterim;
      end
  end

  b=[];
  for i=1:size(B,1)
    b=[b;zeros(B(i,2),1)+B(i,1)];
 end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function qmr_test ( )

%*****************************************************************************80
%
%% QMR_TEST tests QMR.
%
%  Modified:
%
%    27 March 2006
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QMR_TEST:\n' );
  fprintf ( 1, '  QMR uses the Quasi-Minimum Residual \n' );
  fprintf ( 1, '  iterative method to approximate the solution \n' );
  fprintf ( 1, '  of a linear system A * x = b.\n' );

  n = 10;

  A = zeros ( n, n );

  for i = 1 : n
    A(i,i) = 2.0;
  end
  for i = 1 : n-1
    A(i,i+1) = -1;
  end
  for i = 2 : n
    A(i-1,i) = -1;
  end
  x = [ 1 : n ]';
  b = A * x;
  x = ones ( n, 1 );

  M = 0;
  max_it = 10;
  tol = 0.0001;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  For this example, the order of the system is N = %d\n', n );
  fprintf ( 1, '  The matrix A is the simple tridiagonal -1, 2, -1.\n' );
  fprintf ( 1, '  The correct solution is x = [ 1, 2, ..., n].\n' );
  fprintf ( 1, '  The right hand side b is determined by computing A * x.\n' );
  fprintf ( 1, '  The exact x is then replaced by a vector of all 1''s for\n' );
  fprintf ( 1, '  use as a starting guess.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Other parameters are set as follows:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The maximum number of steps is %d.\n', max_it );
  fprintf ( 1, '  The error tolerance is %f\n', tol );

  [ x, error_norm, iter, flag ] = qmr ( A, x, b, M, max_it, tol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The QMR routine has returned with FLAG = %d\n', flag );
  if ( flag == 0 )
    fprintf ( 1, '  This indicates that the iteration has converged.\n' );
  elseif ( flag == 1 )
    fprintf ( 1, '  This indicates that the iteration has NOT converged.\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of iterations taken was %d\n', iter );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The L2 norm of the error per iteration:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : iter
    fprintf ( 1, '  %4d  %f\n', i, error_norm(i) );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The computed solution vector X:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %4d  %f\n', i, x(i) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'QMR_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end


text:
function varargout = panel_bem(varargin)
% PANEL_BEM: Options for the construction of BEM layers (GUI).
%
% USAGE:  bstPanelNew = panel_bem('CreatePanel')
%                   s = panel_bem('GetPanelContents')

% @=============================================================================
% 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, 2011

eval(macro_method);
end

%% ===== CREATE PANEL =====
function [bstPanelNew, panelName] = CreatePanel()  %#ok<DEFNU>
    panelName = 'BemOptions';
    % Java initializations
    import java.awt.*;
    import javax.swing.*;
    % Create main main panel
    jPanelNew = gui_river();
    % Default OpenMeeg options
    BemOptions.nvert = [1922, 1922, 1922];
    BemOptions.thickness = [7 4 3];  % outer-skin, inner-outer, brain-inner
    VertList = {'162', '273', '362', '482', '642', '812', '1082', '1442', '1922', '2432', '2562', '3242', '4322', '5762', '7682', '7292', '9722', '10242', '12962'};

    % ===== BEM LAYERS =====
    jPanelVert = gui_river([4,4], [3,15,10,10], 'Number of vertices per layer');
        gui_component('label', jPanelVert, '', 'Best results are obtained with:
outer skull=inner skull

', [], [], [], []);
        % Scalp
        gui_component('label', jPanelVert, 'br', 'Scalp: ', [], [], [], []);
        jComboVertScalp = gui_component('combobox', jPanelVert, 'tab', [], {VertList}, [], [], []);
        jComboVertScalp.setSelectedIndex(8);
        % Outer skull
        gui_component('label', jPanelVert, 'br', 'Outer skull: ', [], [], [], []);
        jComboVertOuter = gui_component('combobox', jPanelVert, 'tab', [], {VertList}, [], [], []);
        jComboVertOuter.setSelectedIndex(8);
        % Inner skull
        gui_component('label', jPanelVert, 'br', 'Inner skull: ', [], [], [], []);
        jComboVertInner = gui_component('combobox', jPanelVert, 'tab', [], {VertList}, [], [], []);
        jComboVertInner.setSelectedIndex(8);
    jPanelNew.add('br hfill', jPanelVert);

    % ===== BEM THICKNESS =====
    jPanelThick = gui_river([4,4], [0,15,12,10], 'Thickness of the layers');
%         gui_component('label', jPanelThick, '', 'Relative values, scaled for each subject
by the distance between brain and skin.

', [], [], [], []);
%         % Scalp-Outer
%         gui_component('label', jPanelThick, 'br', 'Scalp / outer skull: ', [], [], [], []);
%         jTextThickScalp = gui_component('texttime', jPanelThick, 'tab', num2str(BemOptions.thickness(1), '%d'), [], [], [], []);
%         % Outer-Inner
%         gui_component('label', jPanelThick, 'br', 'Outer skull / inner skull: ', [], [], [], []);
%         jTextThickOuter = gui_component('texttime', jPanelThick, 'tab', num2str(BemOptions.thickness(2), '%d'), [], [], [], []);
%         % Inner-Brain
%         gui_component('label', jPanelThick, 'br', 'Inner skull / brain: ', [], [], [], []);
%         jTextThickInner = gui_component('texttime', jPanelThick, 'tab', num2str(BemOptions.thickness(3), '%d'), [], [], [], []);
        gui_component('label', jPanelThick, 'br', 'Skull (mm): ', [], [], [], []);
        jTextThickOuter = gui_component('texttime', jPanelThick, 'tab', num2str(BemOptions.thickness(2), '%d'), [], [], [], []);
    jPanelNew.add('br hfill', jPanelThick);

    % ===== VALIDATION BUTTONS =====
    gui_component('button', jPanelNew, 'br right', 'Cancel', [], [], @ButtonCancel_Callback, []);
    gui_component('button', jPanelNew, [], 'OK', [], [], @ButtonOk_Callback, []);

    % ===== PANEL CREATION =====
    % Return a mutex to wait for panel close
    bst_mutex('create', panelName);
    % Controls list
    ctrl = struct('BemOptions',       BemOptions, ...
                  'jComboVertScalp',  jComboVertScalp, ...
                  'jComboVertOuter',  jComboVertOuter, ...
                  'jComboVertInner',  jComboVertInner, ...
                  'jTextThickOuter',  jTextThickOuter);
%                   'jTextThickScalp',  jTextThickScalp, ...
%                   'jTextThickOuter',  jTextThickOuter, ...
%                   'jTextThickInner',  jTextThickInner);
    % Create the BstPanel object that is returned by the function
    % => constructor BstPanel(jHandle, panelName, sControls)
    bstPanelNew = BstPanel(panelName, jPanelNew, ctrl);
    % Update panel
    UpdatePanel();

%% =================================================================================
%  === INTERNAL CALLBACKS ==========================================================
%  =================================================================================
%% ===== CANCEL BUTTON =====
    function ButtonCancel_Callback(hObject, event)
        % Close panel without saving (release mutex automatically)
        gui_hide(panelName);
    end

%% ===== OK BUTTON =====
    function ButtonOk_Callback(varargin)
        % Release mutex and keep the panel opened
        bst_mutex('release', panelName);
    end

%% ===== UPDATE PANEL =====
    function UpdatePanel(varargin)

    end
end

%% =================================================================================
%  === EXTERNAL CALLBACKS ==========================================================
%  =================================================================================
%% ===== GET PANEL CONTENTS =====
function s = GetPanelContents() %#ok<DEFNU>
    % Get panel controls
    ctrl = bst_get('PanelControls', 'BemOptions');
    % Bem layers
    nvert1 = str2double(char(ctrl.jComboVertScalp.getSelectedItem()));
    nvert2 = str2double(char(ctrl.jComboVertOuter.getSelectedItem()));
    nvert3 = str2double(char(ctrl.jComboVertInner.getSelectedItem()));
    %thick1 = str2double(char(ctrl.jTextThickScalp.getText()));
    thick2 = str2double(char(ctrl.jTextThickOuter.getText()));
    %thick3 = str2double(char(ctrl.jTextThickInner.getText()));
    % Check for errors
    if isnan(nvert1) || isnan(nvert2) || isnan(nvert3) || isnan(thick2)
        error('Invalid values.');
    end
    s.nvert = [nvert1 nvert2 nvert3];
    s.thickness = [ctrl.BemOptions.thickness(1) thick2 ctrl.BemOptions.thickness(3)];
end

author: brainstorm-tools


text:
function [pb,theta] = pbTG(im,radius,norient)
% function [pb,theta] = pbTG(im,radius,norient)
%
% Compute probability of boundary using TG.
%
% David R. Martin <dmartin@eecs.berkeley.edu>
% April 2003

if nargin<2, radius=0.02; end
if nargin<3, norient=8; end

% beta from logistic fits (trainTG.m)
if radius==0.02, % 64 textons
  beta = [ -4.7151584e+00  1.2222425e+00 ];
  fstd = [  1.0000000e+00  1.9171689e-01 ];
  beta = beta ./ fstd;
else
  error(sprintf('no parameters for radius=%g\n',radius));
end

% get gradients
[tg,gtheta] = detTG(im,radius,norient);

% compute oriented pb
[h,w,unused] = size(im);
pball = zeros(h,w,norient);
for i = 1:norient,
  t = tg(:,:,i); t = t(:);
  x = [ones(size(t)) t];
  pbi = 1 ./ (1 + (exp(-x*beta')));
  pball(:,:,i) = reshape(pbi,[h w]);
end

% nonmax suppression and max over orientations
[unused,maxo] = max(pball,[],3);
pb = zeros(h,w);
theta = zeros(h,w);
r = 2.5;
for i = 1:norient,
  mask = (maxo == i);
  a = fitparab(pball(:,:,i),r,r,gtheta(i));
  pbi = nonmax(max(0,a),gtheta(i));
  pb = max(pb,pbi.*mask);
  theta = theta.*~mask + gtheta(i).*mask;
end
pb = max(0,min(1,pb));

% mask out 1-pixel border where nonmax suppression fails
pb(1,:) = 0;
pb(end,:) = 0;
pb(:,1) = 0;
pb(:,end) = 0;

author: Cloud-CV

repo: object-proposals


text:
function sqi = ssqi(signal)
%sSQI Skewness SQI
%
% Returns the skewness of a signal
%
% Reference:
% Li, Q., Mark, R. G., Clifford, G. D., & Li. (2008). Robust heart rate
% estimation from multiple asynchronous noisy sources using signal quality
% indices and a Kalman filter. Physiol. Meas., 29(1), 15–32.
% http://doi.org/10.1088/0967-3334/29/1/002
%
% Input:
%   signal:         single channel (F)ECG [1xN double]
%
% Output:
%   sqi:            resulting sSQI for segment
%
% Fetal Extraction Toolbox, version 1.0, February 2014
% Released under the GNU General Public License
%
% Copyright (C) 2014 Fernando Andreotti
% Dresden University of Technology, Institute of Biomedical Engineering
% fernando.andreotti@mailbox.tu-dresden.de
%
% Last updated : 09-03-2014
%
%
% 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.

sqi=skewness(signal);

end

author: fernandoandreotti


text:
function [lambda, r2] = duration_amplitude(self, law, min_amplitude, mag_zone)
    %DURATION_AMPLITUDE Use the duration-amplitde
    %Compute the fraction of a data series above each
    %amplitude level, plot the duration vs amplitude data, and then
    %allow user input to fit a regression line.
    %   rsamObject.duration_amplitude(law, min_amplitude)
    %
    %   Inputs:
    %       law = 'exponential' or 'power'
    %       min_amplitude = (Optional) the smallest amplitude to use on
    %       the x-axis. Otherwise will be 10^10 times smallest than the
    %       largest amplitude.
    %
    %   Outputs:
    %       (None) A graph is plotted, the user clicks two points, and
    %       from that slope, the characteristic amplitude is computed
    %       and shown on the screen.

    y = self.data;
    n = length(y);
    a = abs(y);
    max_amplitude = max(a);
    if ~exist('min_amplitude', 'var')
        min_amplitude = max([min(a) max(a)*1e-10]);
    end

    % Method 1
    index=0;
    x = min_amplitude;
    while x < max_amplitude,
        i = find(a>x);
        f = length(i);
        index = index+1;
        frequency(index) = f;
        threshold(index) = x;
        x = x * 1.2;
    end
    clear x y a  f  n  min_amplitude max_amplitude index ;

    %             % Method 2
    %             threshold = [0.0 logspace(min_amplitude, max_amplitude, 50)];
    %             nsamples=[];
    %             for d = 1:length(threshold)
    %                 i = find(a > threshold(d));
    %                 frequency(d) = length(i)/length(y);
    %             end
    %             clear d, nsamples, i, y

    %% PLOT_DURATION_AMPLITUDE
    % plot graph, user select two points, compute
    % characteristic from slope.
    % Use different method depending on whether it is a
    % power law or exponential.

    % define x and y
    switch law
        case {'exponential'}
            x = threshold;
            xlabelstr = 'RMS Displacement(nm)';
        case {'power'}
            x = log10(threshold);
            xlabelstr = 'log10(RMS Displacement(nm))';
        otherwise
            error('law unknown')
    end
    y=log10(frequency);

    % plot duration-amplitude data as circles
    figure
    plot(x,y,'o');
    xlabel(xlabelstr);
    ylabel('log10(Cumulative Minutes)');
    hold on;
    %set(gca,'XLim',[xmin xmax]);

    lambda=0;
    r2=0;

    % check if we have pre-set the magnitude range, effectively
    % our x1 and x2 click points with ginput
    if exist('mag_zone','var') % no user select
        switch law
            case {'power'}
                x1=min(mag_zone);
                x2=max(mag_zone);
            case {'exponential'}
                x1=10^min(mag_zone);
                x2=10^max(mag_zone);

            otherwise
                error('law unknown')
        end
        if x1<min(x)
            x2=min(x);
        end
        y1=interp1(x,y,x1);
        if x2>max(x)
            x2=max(x);
        end
        y2=interp1(x,y,x2);

        % draw a dotted line to show where user selected
        %plot([x1 x2], [y1 y2], '-.');

        % select requested data range and do a least squares fit
        ii = find(x >= x1 & x <= x2);
        wx = x(ii);
        wy = y(ii);
        [p,S]=polyfit(wx,wy,1);
        yfit = polyval(p,wx);
        thiscorr = corrcoef(wy, yfit)
        %try
        if numel(thiscorr)>1
            r2 = thiscorr(1,2);

            % compute lambda
            switch law
                case {'exponential'}
                    lambda = -p(1)/log10(exp(1));
                case {'power'}
                    lambda = -p(1);
                otherwise
                    error('law unknown')
            end

            disp(sprintf('characteristic D_R_S=%.2f cm^2, R^2=%.2f',lambda,r2));

            % draw the fitted line
            xf = [min(wx) max(wx)];
            yf = xf * p(1) + p(2);
            plot(xf, yf,'-');

            %ylabel('log10(t/t0)');
            %xlabel(sprintf('D_R_S (%s) (cm^2)',measure));

            % Add legend
            yrange=get(gca,'YLim');
            xlim = get(gca,'XLim');
            xmax=max(xlim);

            xpos = xmax*0.65;
            ypos = (yrange(2)-yrange(1))*0.8;
            r2str=sprintf('%.2f',r2);
            lambdastr=sprintf('%.2f',lambda);
            if strcmp(law,'exponential')
                tstr = [' \lambda=',lambdastr,' R^2=',r2str];
            else
                tstr = [' \gamma=',lambdastr,' R^2=',r2str];
            end

            text(xpos, ypos, tstr, ...
                'FontName','Helvetica','FontSize',[14],'FontWeight','bold');
        else
            lambda=NaN;
            r2=NaN;
        end

    else

        % user select a range of data
        disp('Left-click Select lowest X, any other mouse button to ignore this station')
        [x1, y1, button1]=ginput(1);
        if button1==1
            disp('Left-click Select highest X, any other mouse button to ignore this station')
            [x2, y2, button2]=ginput(1);
            if button2==1
                if x2>x1
                   % draw a dotted line to show where user selected
                    plot([x1 x2], [y1 y2], '-.');

                    % select requested data range and do a least squares fit
                    ii = find(x >= x1 & x <= x2);
                    wx = x(ii);
                    wy = y(ii);
                    [p,S]=polyfit(wx,wy,1);
                    yfit = polyval(p,wx);
                    thiscorr = corrcoef(wy, yfit)

                    r2 = thiscorr(1,2);

                    % compute lambda
                    switch law
                        case {'exponential'}
                            lambda = -p(1)/log10(exp(1));
                        case {'power'}
                            lambda = -p(1);
                        otherwise
                            error('law unknown')
                    end

                    disp(sprintf('characteristic D_R_S=%.2f cm^2, R^2=%.2f',lambda,r2));

                    % draw the fitted line
                    xf = [min(wx) max(wx)];
                    yf = xf * p(1) + p(2);
                    plot(xf, yf,'-');

                    %ylabel('log10(t/t0)');
                    %xlabel(sprintf('D_R_S (%s) (cm^2)',measure));

                    % Add legend
                    yrange=get(gca,'YLim');
                    xlim = get(gca,'XLim');
                    xmax=max(xlim);

                    xpos = xmax*0.65;
                    ypos = (yrange(2)-yrange(1))*0.8;
                    r2str=sprintf('%.2f',r2);
                    lambdastr=sprintf('%.2f',lambda);
                    if strcmp(law,'exponential')
                        tstr = [self.sta,' \lambda=',lambdastr,' R^2=',r2str];
                    else
                        tstr = [self.sta,' \gamma=',lambdastr,' R^2=',r2str];
                    end

                    text(xpos, ypos, tstr, ...
                        'FontName','Helvetica','FontSize',[14],'FontWeight','bold');

                end
            end
        end
    end
end

author: geoscience-community-codes

repo: GISMO


text:
function R = Rroll(phi)

c = cos(phi); s = sin(phi);
R = [1 0 0; 0 c -s; 0 s c];

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
function cosmo_disp(x,varargin)
% display the input as a string representation
%
% cosmo_disp(x,opt)
%
% Inputs:
%   x              any type of data element (can be a dataset struct)
%   opt            Optional struct with fields
%     .threshold   If the number of values in an array along a dimension
%                  exceeds threshold, then an array is showed in summary
%                  style along that dimension. Default: 5
%     .edgeitems   When an array is shown in summary style, edgeitems sets
%                  the number of items at the beginning and end of the
%                  array to be shown (separated by '...' in rows and by ':'
%                  in columns).
%                  Default: 3
%     .precision   Numeric precision, indicating number of decimals after
%                  the floating point
%                  Default: 3
%     .strlen      Maximal string lenght, if a string is longer the
%                  beginning and end are shown separated by ' ... '.
%                  Default: 20
%     .depth       Maximum recursion depth
%                  Default: 6
%
% Side effect:     Calling this function caused the representation of x
%                  to be displayed.
%
%
% Examples:
%     % display a complicated data structure
%     x=struct();
%     x.a_cell={[],{'cell within cell',[1 2; 3 4]}};
%     x.small_matrix=[10 11 12; 13 14 15];
%     x.big_matrix=reshape(1:200,10,20);
%     x.huge=2^40;
%     x.tiny=2^-40;
%     x.a_string='hello world';
%     x.a_struct.another_struct.name='me';
%     x.a_struct.another_struct.func=@abs;
%     cosmo_disp(x);
%     %|| .a_cell
%     %||   { [  ]  { 'cell within cell'  [ 1         2
%     %||                                   3         4 ] } }
%     %|| .small_matrix
%     %||   [ 10        11        12
%     %||     13        14        15 ]
%     %|| .big_matrix
%     %||   [  1        11        21  ...  171       181       191
%     %||      2        12        22  ...  172       182       192
%     %||      3        13        23  ...  173       183       193
%     %||      :         :         :        :         :         :
%     %||      8        18        28  ...  178       188       198
%     %||      9        19        29  ...  179       189       199
%     %||     10        20        30  ...  180       190       200 ]@10x20
%     %|| .huge
%     %||   [ 1.1e+12 ]
%     %|| .tiny
%     %||   [ 9.09e-13 ]
%     %|| .a_string
%     %||   'hello world'
%     %|| .a_struct
%     %||   .another_struct
%     %||     .name
%     %||       'me'
%     %||     .func
%     %||       @abs
%     %
%     cosmo_disp(x.a_cell)
%     %|| { [  ]  { 'cell within cell'  [ 1         2
%     %||                                 3         4 ] } }
%     cosmo_disp(x.a_cell{2}{2})
%     %|| [ 1         2
%     %||   3         4 ]
%
%     % make a cell in a cell in a cell in a cell ...
%     m={{{{{{{{{{{'hello'}}}}}}}}}}};
%     cosmo_disp(m)
%     %|| { { { { { { <cell> } } } } } }
%     cosmo_disp(m,'depth',8)
%     %|| { { { { { { { { <cell> } } } } } } } }
%     cosmo_disp(m,'depth',Inf)
%     %|| { { { { { { { { { { { 'hello' } } } } } } } } } } }
%
%     % illustrate 'threshold' and 'edgeitems' arguments
%     cosmo_disp(num2cell('a':'k'))
%     %|| { 'a'  'b'  'c' ... 'i'  'j'  'k'   }@1x11
%     cosmo_disp(num2cell('a':'k'),'threshold',Inf)
%     %|| { 'a'  'b'  'c'  'd'  'e'  'f'  'g'  'h'  'i'  'j'  'k' }
%     cosmo_disp(num2cell('a':'k'),'edgeitems',2)
%     %|| { 'a'  'b' ... 'j'  'k'   }@1x11
%
%     % illustrate 'precision' argument
%     cosmo_disp(pi*[1 2],'precision',1);
%     %|| [ 3       6 ]
%     cosmo_disp(pi*[1 2],'precision',3);
%     %|| [ 3.14      6.28 ]
%     cosmo_disp(pi*[1 2],'precision',5);
%     %|| [ 3.1416      6.2832 ]
%     cosmo_disp(pi*[1 2],'precision',7);
%     %|| [ 3.141593      6.283185 ]
%
%     % illustrate n-dimensional arrays
%     x=zeros([2 2 1 2 3]);
%     x(:)=2*(1:numel(x));
%     cosmo_disp(x)
%     %|| <double>@2x2x1x2x3
%     %||    (:,:,1,1,1) =  [ 2         6
%     %||                     4         8 ]
%     %||    (:,:,1,2,1) =  [ 10        14
%     %||                     12        16 ]
%     %||    (:,:,1,1,2) =  [ 18        22
%     %||                     20        24 ]
%     %||    (:,:,1,2,2) =  [ 26        30
%     %||                     28        32 ]
%     %||    (:,:,1,1,3) =  [ 34        38
%     %||                     36        40 ]
%     %||    (:,:,1,2,3) =  [ 42        46
%     %||                     44        48 ]
%     cosmo_disp(reshape(char(65:72),[2 2 2]))
%     %|| <char>@2x2x2
%     %||    (:,:,1) = 'AC
%     %||               BD'
%     %||    (:,:,2) = 'EG
%     %||               FH'
%     cosmo_disp(zeros([2 3 5 7 0 2]))
%     %|| <double>@2x3x5x7x0x2 (empty)
%
%     % illustrate non-singleton structs
%     x=struct('x',{1 2; 3 4});
%     cosmo_disp(x);
%     %|| <struct>@2x2
%     %||    (1,1).x
%     %||           [ 1 ]
%     %||    (2,1).x
%     %||           [ 3 ]
%     %||    (1,2).x
%     %||           [ 2 ]
%     %||    (2,2).x
%     %||           [ 4 ]
%     x3=cat(3,x,x,x);
%     cosmo_disp(x3);
%     %|| <struct>@2x2x3
%     %||    (1,1,1).x
%     %||             [ 1 ]
%     %||    (2,1,1).x
%     %||             [ 3 ]
%     %||    (1,2,1).x
%     %||             [ 2 ]
%     %||      :        :
%     %||    (2,1,3).x
%     %||             [ 3 ]
%     %||    (1,2,3).x
%     %||             [ 2 ]
%     %||    (2,2,3).x
%     %||             [ 4 ]
%
% Notes:
%   - Unlike the builtin 'disp' function, this function shows the contents
%     of the input using recursion. For example if a cell contains a
%     struct, then the contents of that struct is shown as well
%   - A use case is displaying dataset structs
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    defaults.threshold=5;    % max #items before triggering summary style
    defaults.edgeitems=3;    % #items at edges in summary style
    defaults.precision=3;    % show floats with 3 decimals
    defaults.strlen=20;      % insert '...' with strings more than 20 chars
    defaults.depth=6;        % maximal depth
    defaults.show_size=false;% whether to show size of matrices

    opt=cosmo_structjoin(defaults,varargin);

    % get string representation of x
    s=disp_helper(x, opt);

    % print string representation of x
    disp(s);

function s=disp_helper(x, opt)
    % general helper function to get a string representation. Unlike the
    % main function this function returns a string, which makes it suitable
    % for recursion
    depth=opt.depth;
    if depth<=0
        s=any2summary_str(x,opt);
        return
    end

    opt.depth=depth-1;

    if ~has_size(x)
        s=any2summary_str(x,opt);
    else
        s=nd_any2str(x,opt);
    end

function s=any2summary_str(x,unused)
    if has_size(x)
        sz=size(x);
    else
        sz=[1 1];
    end
    s=surround_with(true,'<',class(x),'>',sz);

function tf=has_size(x)
    % helper, because some classes have no 'size'
    try
        size(x);
        tf=true;
    catch
        tf=false;
    end

function y=nd_any2str(x,opt)
    if isstruct(x)
        if numel(x)<=1
            y=struct2str(x,opt);
        else
            y=multi_any2string(x,0,opt);
        end
    elseif numel(size(x))==2
        if iscell(x)
            y=cell2str(x,opt);
        elseif isnumeric(x) || islogical(x)
            y=matrix2str(x,opt);
        elseif ischar(x)
            y=string2str(x,opt);
        elseif isa(x, 'function_handle')
            y=function_handle2str(x,opt);
        else
            y=any2summary_str(x,opt);
        end
    else
        y=multi_any2string(x,2,opt);
    end

function y=multi_any2string(x,ndim_post,opt)
    sz=size(x);
    sz_rest=sz((ndim_post+1):end);
    ndim_rest=numel(sz_rest);
    n_rest=prod(sz_rest);

    parts=cell(ndim_rest,1);
    for k=1:ndim_rest
        parts{k}=num2cell(1:sz_rest(k));
    end

    p=cosmo_cartprod(parts);
    if ndim_post>0
        xflat=reshape(x,[sz(1:ndim_post) n_rest]);
    else
        xflat=x(:);
    end

    [pre,post]=get_mx_idxs(xflat,opt.edgeitems,opt.threshold,ndim_post+1);

    header=any2summary_str(x,opt);
    s_pre=nd_any2str_helper(xflat,p, pre,opt);

    if isempty(post)
        s_post={'',''};
        s_dots={'',''};
    else
        s_post=nd_any2str_helper(xflat,p,post,opt);
        s_dots=cell(1,2);
        for k=1:2
            szs=cellfun(@(x)size(x,2),s_post(:,k));
            pos=round(max(szs)/2);
            s_dots{k}=[spaces(1,pos) ':'];
        end
    end
    s_all=cat(1,s_pre,s_dots,s_post);

    y=strcat_({header;strcat_(s_all)});

function s=nd_any2str_helper(xflat, p, idxs, opt)
    n_rest=numel(idxs);
    s=cell(n_rest,2);
    sz=size(xflat);
    npre=numel(sz)-1;
    for k=1:n_rest
        idx=idxs(k);
        switch npre
            case 1
                % for struct
                v=xflat(idx);
                idx_prefix='';
                idx_postfix='';
            case 2
                % anything else
                v=xflat(:,:,idx);
                idx_prefix=':,:,';
                idx_postfix=' = ';
            otherwise
                assert(false);
        end

        v_str=disp_helper(v,opt);
        idx_str=sprintf(',%d',p(idx,:));
        s{k,1}=sprintf('   (%s%s)%s',idx_prefix,idx_str(2:end),idx_postfix);
        s{k,2}=v_str;
    end

function y=strcat_(xs)
    if isempty(xs)
        y='';
        return
    end

    % all elements in xs are char
    [nr,nc]=size(xs);
    ys=cell(1,nc);

    % sizes for each element
    width_per_col=max_element_size(xs,2);
    height_per_row=max_element_size(xs,1);
    for k=1:nc
        xcol=cell(nr,1);
        width=width_per_col(k);
        row_pos=0;
        for j=1:nr
            height=height_per_row(j);
            if height==0
                continue;
            end

            x=xs{j,k};
            if ~ischar(x) && isempty(x)
                x='';
            end

            sx=size(x);
            to_add=[height width]-sx;

            % pad with spaces
            row_pos=row_pos+1;
            xcol{row_pos}=[[x spaces(sx(1),to_add(2))];...
                        spaces(to_add(1), width)];
        end
        ys{k}=char(xcol{1:row_pos});
    end
    y=[ys{:}];

function m=max_element_size(x,dim)
    % faster than cellfun
    n=numel(x);
    sizes=zeros(size(x));
    for k=1:n
        sizes(k)=size(x{k},dim);
    end
    m=max(sizes,[],3-dim);

function y=spaces(nx,ny)
    % faster than repmat(' ',nx,ny)
    if nx>0 && ny>0
        y(nx,ny)=' ';
        y(:)=' ';
    else
        if nx<0
            nx=0;
        end
        if ny<0
            ny=0;
        end
        y=reshape('',nx,ny);
    end

function y=struct2str(x,opt)
    if numel(x)==0
        show_size=opt.show_size;
        y=[surround_with(show_size,'', 'struct', '', size(x)) ' (empty)'];
        return;
    end

    assert(numel(x)==1)
    fns=fieldnames(x);
    n=numel(fns);

    if n==0
        show_size=opt.show_size;
        y=[surround_with(show_size,'', 'struct', '', size(x)) ' (empty)'];
    else
        r=cell(n*2,1);
        for k=1:n
            fn=fns{k};
            r{k*2-1}=['.' fn];
            d=disp_helper(x.(fn),opt);
            r{k*2}=[spaces(size(d,1),2) d];
        end
        y=strcat_(r);
    end

function s=function_handle2str(x,opt)
    s_with_quotes=string2str(func2str(x),opt);
    s=['@' s_with_quotes(2:(end-1))];

function s=string2str(x, opt)
    if ~ischar(x), error('expected a char'); end
    [nrows,ncols]=size(x);

    if ncols>opt.strlen
        infix=' ... ';
        h=floor((opt.strlen-numel(infix))/2);
        x=strcat_({x(:,1:h), infix ,x(:,ncols+((1-h):0))});
    end
    quote='''';
    pre=quote;
    post=[spaces(nrows-1,1);quote];
    s=strcat_({pre,x,post});

function s=cell2str(x, opt)
    % display a cell

    edgeitems=opt.edgeitems;
    threshold=opt.threshold;

    % get indices of rows and columns to show
    [r_pre, r_post]=get_mx_idxs(x, edgeitems, threshold, 1);
    [c_pre, c_post]=get_mx_idxs(x, edgeitems, threshold, 2);

    part_idxs={{r_pre, r_post}, {c_pre, c_post}};

    nrows=numel([r_pre r_post])+~isempty(r_post);
    ncols=numel([c_pre c_post])+~isempty(c_post);

    sinfix=cell(nrows,ncols*2+1);
    for k=1:(ncols-1)
        sinfix{1,k*2+2}='  ';
    end

    cpos=1;
    for cpart=1:2
        col_idxs=part_idxs{2}{cpart};
        nc=numel(col_idxs);

        rpos=0;
        for rpart=1:2
            row_idxs=part_idxs{1}{rpart};

            nr=numel(row_idxs);
            if nr==0
                continue
            end
            for ci=1:nc
                col_idx=col_idxs(ci);
                trgc=cpos+ci*2;
                for ri=1:nr
                    row_idx=row_idxs(ri);
                    sinfix{rpos+ri,trgc}=disp_helper(x{row_idx,...
                                                             col_idx},opt);
                    if cpart==2 && ci==1 && nc>0
                        sinfix{rpos+ri,cpos+ci*2-1}=' ... ';
                    end
                end

                if rpart==2
                    max_length=max(cellfun(@numel,sinfix(:,trgc)));
                    pre_spaces=spaces(1,floor(max_length/2-1));
                    sinfix{rpos,cpos+ci*2}=[pre_spaces ':'];
                end
            end
            rpos=rpos+nr+1;
        end
        cpos=cpos+nc*2;
    end

    show_size=opt.show_size || ~isempty(r_post) || ~isempty(c_post);
    s=surround_with(show_size,'{ ', strcat_(sinfix), ' }', size(x));

function pre_infix_post=surround_with(show_size, pre, infix, post, matrix_sz)
    % surround infix by pre and post, doing
    n=prod(matrix_sz);
    if show_size && n~=1
        size_str=sprintf('x%d',matrix_sz);
        size_str(1)='@';
        if n==0
            size_str=[size_str ' (empty)'];
        end
    else
        size_str='';
    end
    post=strcat_({spaces(size(infix,1)-1,1); [post size_str]});
    pre_infix_post=strcat_({pre, infix, post});

function s=matrix2str(x,opt)
    if isempty(x)
        show_size=opt.show_size;
        s=surround_with(show_size,'[','  ',']',size(x));
        return
    end

    % display a matrix
    edgeitems=opt.edgeitems;
    threshold=opt.threshold;
    precision=opt.precision;

    % get indices of rows and columns to show
    [r_pre, r_post]=get_mx_idxs(x, edgeitems, threshold, 1);
    [c_pre, c_post]=get_mx_idxs(x, edgeitems, threshold, 2);

    % data to be shown
    y=x([r_pre r_post],[c_pre c_post]);

    % convert to string
    s=num2str(y,precision);

    % number of characters in first and second dimension
    [nc_row,nc_col]=size(s);

    % see where each column is a space; that's a potential split point
    sp_col=sum(s==' ',1)==nc_row;

    % col_index has value k for characters in the k-th column, else zero
    col_index=zeros(1,nc_col);
    col_count=1;
    in_num=true;
    for k=1:nc_col
        if in_num
            if sp_col(k)
                col_count=col_count+1;
                in_num=false;

            else
                col_index(k)=col_count;
            end
        elseif ~sp_col(k)
            in_num=true;
            col_index(k)=col_count;
        end
    end

    % deal with rows
    row_blocks=cell(3,1);
    if isempty(r_post)
        row_blocks{1,1}=s;
    else
        % insert ':' for each column
        line=spaces(1,nc_col);
        for k=1:max(col_index)
            idxs=find(col_index==k);
            median_pos=round(mean(idxs));
            line(median_pos)=':';
        end
        row_blocks{1}=s(1:edgeitems,:);
        row_blocks{2}=line;
        row_blocks{3}=s(edgeitems+(1:edgeitems),:);
    end

    % deal with columns
    row_and_col_blocks=cell(3,3);
    for row=1:3
        if isempty(c_post)
            row_and_col_blocks{row}=row_blocks{row};
        else
            % insert ' ... ' halfway each row
            pre_end=find(col_index==edgeitems,1,'last')+1;
            post_start=find(col_index==(edgeitems+1),1,'first')-1;

            r=row_blocks{row,1};
            if isempty(r)
                continue;
            end
            row_and_col_blocks{row,1}=r(:,1:pre_end);
            if row~=2
                row_and_col_blocks{row,2}=repmat(' ... ',size(r,1),1);
            end
            row_and_col_blocks{row,3}=r(:,post_start:end);
        end
    end

    show_size=opt.show_size || ~isempty(r_post) || ~isempty(c_post);
    s=surround_with(show_size,'[ ',strcat_(row_and_col_blocks),' ]',...
                                                                size(x));

function [pre,post]=get_mx_idxs(x, edgeitems, threshold, dim)
    % returns the first and last indices for showing an array along
    % dimension dim. If size(x,dim)<2*edgeitems, then pre has all the
    % indices, otherwise pre and post have the first and last edgeitems
    % indices, respectively
    n=size(x,dim);

    if n>max(threshold,2*edgeitems) % properly deal with Inf values
        pre=1:edgeitems;
        post=n-edgeitems+(1:edgeitems);
    else
        pre=1:n;
        post=[];
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [elem,HB] = uniformcoarsenquadred(elem)
%% UNIFORMCOARSEQUADNRED uniform coarsening of red refinement
%
% [elem,HB] = uniformcoarsenred(elem) remove grid points added by uniform
% refinement. See the illustration below:
%
% 4  - 7 -  3
% | t4 | t3 |
% 8 -  9 -  6
% | t1 | t2 |
% 1 -  5 -  2
%
% It is mainly used to get multilevel decomposition in multigrid methods.
% The input matrix elem stands for the fine mesh and the output one for the
% coarse mesh. The HB records the hierarchical structure of added points
% going from the coarse to the fine mesh such that HB(:,2:3) are two parent
% nodes of HB(:,1).
%
%   See also: uniformcoarsen, coarsen, bisect, uniformcoarsen3, mg
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

HB = [];
NT = size(elem,1);
if mod(NT,4)==0
    NTc = NT/4; % number of triangles in the coarse grid
else
%     display('Not from red refinement');
    return
end

%% Find points
t1 = 1:NTc; t2 = t1+NTc; t3 = t2+NTc; t4 = t3+NTc;
if any(elem(t1,2)~=elem(t2,1)) || any(elem(t1,3)~=elem(t3,1)) || ...
   any(elem(t4,1)~=elem(t1,4)) || any(elem(t4,2)~=elem(t3,1))
%     display('Not from red refinement');
    return
end
p1 = elem(t1,1);
p2 = elem(t2,2);
p3 = elem(t3,3);
p4 = elem(t4,4);
p5 = elem(t1,2);
p6 = elem(t2,3);
p7 = elem(t3,4);
p8 = elem(t4,1);
p9 = elem(t1,3);

%% Remove quad
elem(t1,:) = [p1 p2 p3 p4];
elem = elem(t1,:);

%% Record HB
HB(p9,:) = [p9 p1 p3];
HB(p5,:) = [p5 p1 p2];
HB(p6,:) = [p6 p2 p3];
HB(p7,:) = [p7 p4 p3];
HB(p8,:) = [p8 p1 p4];
Nc = max(elem(:));
HB = HB(Nc+1:end,:);

repo: ifem


text:
% [INPUT]
% va = A vector of floats [0,Inf) of length n representing the market values of assets.
% vap = Input argument representing the distributional parameters of assets whose type depends on the chosen option pricing model:
%   - for Black-Scholes-Merton, a float [0,Inf) representing the annualized volatility of assets;
%   - for Gram-Charlier, a row vector of floats (-Inf,Inf) of length 3 whose values represent respectively the annualized volatility, skewness and excess kurtosis of assets.
% cds = A vector of floats [0,Inf) of length n representing the credit default swap spreads.
% db = A float or a vector of floats [0,Inf) of length n representing the default barrier.
% r = A float or a vector of floats (-Inf,Inf) of length n representing the annualized risk-free interest rate.
% t = A float or a vector of floats (0,Inf) of length n representing the time to maturity of default barrier.
%
% [OUTPUT]
% el = A column vector of floats [0,Inf) of length n representing the expected losses.
% cl = A column vector of floats [0,Inf) of length n representing the contingent liabilities.
% a = A column vector of floats [0,1] of length n representing the contingent alphas.

function [el,cl,a] = contingent_claims_analysis(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('va',@(x)validateattributes(x,{'double'},{'real' 'finite' 'nonnegative' '2d' 'nonempty'}));
        ip.addRequired('vap',@(x)validateattributes(x,{'double'},{'real' 'finite' 'vector' 'nonempty'}));
        ip.addRequired('cds',@(x)validateattributes(x,{'double'},{'real' 'nonnegative' 'vector' 'nonempty'}));
        ip.addRequired('db',@(x)validateattributes(x,{'double'},{'real' 'finite' 'nonnegative' 'vector' 'nonempty'}));
        ip.addRequired('r',@(x)validateattributes(x,{'double'},{'real' 'finite' 'vector' 'nonempty'}));
        ip.addRequired('t',@(x)validateattributes(x,{'double'},{'real' 'finite' '>' 0 'vector' 'nonempty'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    [va,vap,cds,db,r,t] = validate_input(ipr.va,ipr.vap,ipr.cds,ipr.db,ipr.r,ipr.t);

    nargoutchk(3,3);

    [el,cl,a] = contingent_claims_analysis_internal(va,vap,cds,db,r,t);

end

function [el,cl,a] = contingent_claims_analysis_internal(va,vap,cds,db,r,t)

    s = vap(1);
    st = s * sqrt(t);

    dbd = db .* exp(-r .* t);

    d1 = (log(va ./ db) + ((r + (0.5 * s^2)) .* t)) ./ st;
    d2 = d1 - st;

    put_price = (dbd .* normcdf(-d2)) - (va .* normcdf(-d1));

    if (numel(vap) == 3)
        g = vap(2);
        k = vap(3);

        t1 = (g / 6) .* ((2 * s) - d1);
        t2 = (k / 24) .* (1 - d1.^2 + (3 .* d1 .* s) - (3 * s^2));

        put_price = put_price - (va .* normcdf(d1) .* s .* (t1 - t2));
    end

    put_price = max(0,put_price);

    rd = dbd - put_price;

    cds_put_price = dbd .* (1 - exp(-cds .* max(0.5,((db ./ rd) - 1)) .* t));
    cds_put_price = min(cds_put_price,put_price);

    a = max(0,min(1 - (cds_put_price ./ put_price),1));
    a(~isreal(a)) = 0;

    el = put_price;
    cl = el .* a;

end

function [va,vap,cds,db,r,t] = validate_input(va,vap,cds,db,r,t)

    va = va(:);
    va_len = numel(va);

    if (va_len < 5)
        error('The value of ''va'' is invalid. Expected input to be a vector containing at least 5 elements.');
    end

    vap_len = numel(vap);

    if ((vap_len ~= 1) && (vap_len ~= 3))
        error('The value of ''vap'' is invalid. Expected input to be a vector containing either 1 or 3 elements.');
    end

    if (vap(1) < 0)
        error('The value of ''vap'' is invalid. Expected input first element to be greater than or equal to 0.');
    end

    cds = cds(:);

    if (numel(va) ~= va_len)
        error(['The value of ''cds'' is invalid. Expected input to be a vector of length ' num2str(va_len) ' elements.']);
    end

    if (all(cds >= 1))
        cds = cds ./ 10000;
    end

    data = {db(:) r(:) t(:)};

    l = unique(cellfun(@numel,data));
    l_scalar = (l == 1);

    if (any(l_scalar))
        if (any(l(~l_scalar) ~= va_len))
            error(['The number of elements of ''db'', ''r'' and ''t'' must be either 1 or equal to ' num2str(va_len) '.']);
        end
    else
        if (any(l ~= va_len))
            error(['The number of elements of ''db'', ''r'' and ''t'' must be either 1 or equal to ' num2str(va_len) '.']);
        end
    end

    for i = 1:numel(data)
        data_i = data{i};

        if (numel(data_i) == 1)
            data{i} = repmat(data_i,va_len,1);
        end
    end

    [db,r,t] = deal(data{:});

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function [A,B,C,G] = sopls(X,Y,nhidden,reg,verbose)
% SOPLS  Sparse orthogonalized partial least squares where the elastic
% net is implemented using either a native implementation or the glmnet package
%
% The native implementation supports arbitrary alpha (ridge penalty) and
% lambda (lasso penalty) values and supports a coupling matrix for alpha to
% implement smoothing. It requires a fixed value for lambda.
%
% The glmnet implementation is much faster but uses alpha to mix between
% ridge and lasso while lambda determines that amount of regularization. It
% also does not support a coupling matrix. It learns the optimal lambda
% using inner cross-validation.
%
% X: ninput x nsamples input data matrix
% Y: noutput x nsamples output data matrix
% nhidden: number of components [1]
% reg: employed regression method
%
% A: noutput x nhidden weight matrix
% B: ninput x nhidden weight matrix
% C: 1 x nhidden bias vector
% G: output parameters for debugging purposes

% Parse inputs

if nargin < 3, nhidden = 1; end
if nargin < 4, error('please specify regularizer'); end
if nargin < 5, verbose = false; end

ninput = size(X,1);
noutput = size(Y,1);
nsamples = size(X,2);

if nhidden > 1,

    G = cell(1,nhidden);

    % Run nhidden times sequentially

    R = Y;
    A = zeros(noutput,nhidden);
    B = zeros(ninput,nhidden);
    C = zeros(1,nhidden);
    for i=1:nhidden,

      if verbose
        fprintf('learning hidden variable %d of %d\n',i,nhidden);
      end

      [A(:,i),B(:,i),C(i),G{i}] = sopls(X,R,1,reg);
      if i < nhidden,    % deflate
        Z = B(:,i)'*X + C(i); % hidden activations
        R = R - A(:,i)*Z; % Y activations after deflation
      end

    end
else

    % Run for a single hidden unit

    B = zeros(ninput,1);
    C = 0;
    iter = 0;
    maxiter = 100;
    tol = nhidden*noutput*(1e-10);

    % Initialize A to first principal component of Y

    optseig.disp = 0;
    if nsamples < noutput,
        [d1,d2] = eigs(Y'*Y,[],1,'LM',optseig);
        A = Y*d1;
        A = A/sqrt(A'*A);
    else
        [A,d2] = eigs(Y*Y',[],1,'LM',optseig);
    end

    Aold = A;
    while iter < maxiter,

        Z = A'*Y;    % reconstruct Z given A from output Y

        try

          f = reg.train(X',Z');

          % all regularizers should use this convention
          B = f.model.weights;
          C = f.model.bias;

        catch

          % this should not happen
          warning(lasterr);
          B = zeros(ninput,1);
          C = 0;

        end

        Z = B'*X + C;   % reconstruct Z given B and C

        % Find optimal A under constraint A'*A = 1

        Syz = Y*Z'/nsamples;
        denom = sqrt(Syz'*Syz);
        if denom,
            A = Syz/denom;
        else
            A = Aold;
        end
        if any(isnan(A(:))),    % check...
            error('nans!!!\n');
        end

        if sumsqr(A - Aold) < tol,
          iter = maxiter;
        else
          iter = iter + 1;
          Aold = A;
        end
    end

    G = f;

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function value = r8_sign_match ( r1, r2 )

%*****************************************************************************80
%
%% R8_SIGN_MATCH is TRUE if two R8's are of the same sign.
%
%  Discussion:
%
%    This test could be coded numerically as
%
%      if ( 0 <= r1 * r2 ) then ...
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R1, R2, the values to check.
%
%    Output, logical VALUE, is TRUE if ( R1 <= 0 and R2 <= 0 )
%    or ( 0 <= R1 and 0 <= R2 ).
%
  value = ( r1 <= 0.0 && r2 <= 0.0 ) || ( 0.0 <= r1 && 0.0 <= r2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ret = write_rgbe(img, filename, hdr_info)
%
%       ret = write_rgbe(img, filename, hdr_info)
%
%       This function writes an image using RGBE encoding.
%	Exposure and gamma, in the hdr_info, have to be applied before
%	by the user.
%
%
%        Input:
%           -img: the image to write on the hard disk
%           -filename: the name of the image to write
%           -hdr_info: RGBE format extra datum such as: RLE comression, exposure, gamma, etc.
%
%        Output:
%           -ret: a boolean value, it is set to 1 if the write succeed, 0 otherwise
%
%     Copyright (C) 2011-14  Francesco Banterle
%
%     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/>.
%

if(~exist('hdr_info', 'var'))
    hdr_info = struct('exposure', 1.0, 'gamma', 1.0, 'bRLE', 1);
end

if(~exist('hdr_info.bRLE','var'))
    bRLE = 1;
else
    bRLE = hdr_info.bRLE;
end

if(~exist('hdr_info.exposure', 'var'))
    exposure = 1.0;
else
    exposure = hdr_info.exposure;
end

if(~exist('hdr_info.gamma', 'var'))
    gamma = 1.0;
else
    gamma = hdr_info.gamma;
end

ret = 0;

[n, m, c] = size(img);

if(c ~= 3)
    error('RGBE encoding requires a three color channels image!');
end

fid = fopen(filename, 'W');
fprintf(fid,'#?RADIANCE\n');
fprintf(fid,'#Generated by the HDR Toolbox by Francesco Banterle\n');
fprintf(fid,'FORMAT=32-bit_rle_rgbe\n');
fprintf(fid,'EXPOSURE= %f\n\n', exposure);
fprintf(fid,'-Y %d +X %d\n', n, m);

%convert from float to RGBE
RGBEbuffer = uint8(float2RGBE(img));

if((n < 8) || (n > 32767)) %RLE encoding is not allowed in these cases
    bRLE = 0;
end

%MIN_RLE_RUN = 4;

if(bRLE)
    buffer_line_start = uint8([2, 2, bitshift(m, -8), bitand(m, 255)]);

    for i=1:n
        fwrite(fid, buffer_line_start, 'uint8');
        for j=1:4
            buffer_line = reshape(RGBEbuffer(i,:,j), m, 1);
            write_rgbe_line(buffer_line, fid);
        end
    end
else
    %data reshape
    tot = n * m * 4;
    data = zeros(tot, 1, 'uint8');
    for i=1:4
        C = i:4:(tot);
        data(C) = reshape(RGBEbuffer(:,:,i)', m * n, 1);
    end

    fwrite(fid, data, 'uint8');

end

fclose(fid);

ret = 1;

end

author: banterle


text:
function Solution = GhiaSolution
%GhiaSolution returns a benchmark solution for lid-driven cavity flow
%   The data is taken from:
%     Ghia, U. K. N. G., Kirti N. Ghia, and C. T. Shin.
%     "High-Re solutions for incompressible flow using the Navier-Stokes
%     equations and a multigrid method."
%     Journal of computational physics 48.3 (1982): 387-411.

Solution.Re=[100 400 1000];

Solution.y = [1.0000,...
              0.9766,...
              0.9688,...
              0.9609,...
              0.9531,...
              0.8516,...
              0.7344,...
              0.6172,...
              0.5000,...
              0.4531,...
              0.2813,...
              0.1719,...
              0.1016,...
              0.0703,...
              0.0625,...
              0.0547,...
              0.0000];
% Re=100
Solution.u{1}=   [1.0000,...
                  0.84123,...
                  0.78871,...
                  0.73722,...
                  0.68717,...
                  0.23151,...
                  0.00332,...
                 -0.13641,...
                 -0.20581,...
                 -0.21090,...
                 -0.15662,...
                 -0.10150,...
                 -0.06434,...
                 -0.04775,...
                 -0.04192,...
                 -0.03717,...
                  0.00000];
% Re=400
Solution.u{2}=   [1.00000,...
                  0.75837,...
                  0.68439,...
                  0.61756,...
                  0.55892,...
                  0.29093,...
                  0.16256,...
                  0.02135,...
                 -0.11477,...
                 -0.17119,...
                 -0.32726,...
                 -0.24299,...
                 -0.14612,...
                 -0.10338,...
                 -0.09266,...
                 -0.08186,...
                  0.00000];
%Re=1000
Solution.u{3}=   [1.0000,...
                  0.65928,...
                  0.57492,...
                  0.51117,...
                  0.46604,...
                  0.33304,...
                  0.18719,...
                  0.05702,...
                 -0.06080,...
                 -0.10648,...
                 -0.27805,...
                 -0.38289,...
                 -0.29730,...
                 -0.22220,...
                 -0.20196,...
                 -0.18109,...
                  0.00000];

Solution.x = [1.0000,...
              0.9688,...
              0.9609,...
              0.9531,...
              0.9453,...
              0.9063,...
              0.8594,...
              0.8047,...
              0.5000,...
              0.2344,...
              0.2266,...
              0.1563,...
              0.0938,...
              0.0781,...
              0.0703,...
              0.0625,...
              0.0000];
%Re=100
Solution.v{1}=[0.00000,...
              -0.05906,...
              -0.07391,...
              -0.08864,...
              -0.10313,...
              -0.16914,...
              -0.22445,...
              -0.24533,...
               0.05454,...
               0.17527,...
               0.17507,...
               0.16077,...
               0.12317,...
               0.10890,...
               0.10091,...
               0.09233,...
               0.00000];
% Re=400
Solution.v{2}=[0.00000,...
              -0.12146,...
              -0.15663,...
              -0.19254,...
              -0.22847,...
              -0.23827,...
              -0.44993,...
              -0.38598,...
               0.05186,...
               0.30174,...
               0.30203,...
               0.28124,...
               0.22965,...
               0.20920,...
               0.19713,...
               0.18360,...
               0.00000];
%Re=1000
Solution.v{3}=[0.00000,...
              -0.21388,...
              -0.27669,...
              -0.33714,...
              -0.39188,...
              -0.51550,...
              -0.42665,...
              -0.31966,...
               0.02526,...
               0.32235,...
               0.33075,...
               0.37095,...
               0.32627,...
               0.30353,...
               0.29012,...
               0.27485,...
               0.00000];

end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function GenerateFunctionMFile(ask)
%GENERATEFUNCTIONMFILE Generates an m-file suitable for a function,
% with header text.
%
% GENERATEFUNCTIONMFILE() generates an m-file suitable for a function,
% with header text, saves it in your current directory, and opens it in
% the editor.
%
% GENERATEFUNCTIONMFILE(ASK) works as above, but prompts for the name of
% the function.

% $Author: rcotton $  $Date: 2010/04/12 10:36:31 $ $Revision: 1.2 $
% Copyright: Health and Safety Laboratory 2010

persistent fileNameCount

if nargin > 0 && ask
   fnName = input('Enter the name of your function >> ', 's');
   if ~isvarname(fnName)
      error('GenerateFunctionMFile:InvalidFunctionName', ...
         [fnName ' is not a valid function name']);
   end
else
   if isempty(fileNameCount)
      fileNameCount = 1;
   else
      fileNameCount = fileNameCount + 1;
   end
   while exist(['Untitled' num2str(fileNameCount)], 'file') >= 2
      fileNameCount = fileNameCount + 1;
   end
   fnName = ['Untitled' num2str(fileNameCount)];
end

fileName = fullfile(cd, [fnName '.m']);

if exist(fnName, 'file')
   error('GenerateFunctionMFile:FileAlreadyExists', ...
      [fnName '.m already exists in the current directory']);
end

WriteFile(fnName, fileName)
edit(fileName);

end

function WriteFile(fnName, fileName)
% Engine for writing the file.

now = clock;
nl = char(10);

   function PrintCommentLine(fid)
   fprintf(fid, '%% %s', nl);
   end

   function PrintBlankLine(fid)
   fprintf(fid, '%s', nl);
   end

fid = fopen0(fileName, 'w');
c = onCleanup(@() fclose(fid));

fprintf(fid, 'function [outputArgs] = %s(inputArgs)%s', ...
   fnName, nl);
fprintf(fid, '%%%s Summary of this function goes here%s', ...
   upper(fnName), nl);
PrintCommentLine(fid);
fprintf(fid, '%% [OUTPUTARGS] = %s(INPUTARGS) Explain usage here%s', ...
   upper(fnName), nl);
PrintCommentLine(fid);
fprintf(fid, '%% Examples: %s', nl);
PrintCommentLine(fid);
fprintf(fid, '%% Provide sample usage code here%s', nl);
PrintCommentLine(fid);
fprintf(fid, '%% See also: List related files here%s', nl);
fprintf(fid, nl);
% Need to keep $ separate in next line, so CVS doesn't replace values
% Also, this line will not work with GIT, since it messes with the SHA-1
fprintf(fid, '%% %sAuthor: %s $\t%sDate: %s $\t%sRevision: 0.1 $%s', ...
   '$', GetUserLoginName(), '$', ...
   datestr(now, 'yyyy/mm/dd HH:MM:SS'), '$', nl);
fprintf(fid, '%% Copyright: Health and Safety Laboratory %s%s', ...
   datestr(now, 'yyyy'), nl);
PrintBlankLine(fid);
PrintBlankLine(fid);
PrintBlankLine(fid);
fprintf(fid, 'end%s', nl);
end

author: Sable


text:
function [ feature_image ] = get_fhog(im, fparam, gparam )
%extract fhog features using piotrs toolbox. Currently takes no parameters
%except hog-cell-size
if ~isfield(fparam, 'nOrients')
    fparam.nOrients = 9;
end
[im_height, im_width,~, num_images] = size(im);
feature_image = zeros(floor(im_height/gparam.cell_size), floor(im_width/gparam.cell_size), fparam.nDim, num_images, 'single');
for k = 1:num_images
    hog_image = fhog(single(im(:,:,:,k)), gparam.cell_size, fparam.nOrients);
    %the last dimension is all 0 so we can discard it
    feature_image(:,:,:,k) = hog_image(:,:,1:end-1);
end
end

repo: AutoTrack


text:
function [onOff] = convertTfToOnOff(tf)
    if(tf)
        onOff = 'on';
    else
        onOff = 'off';
    end
end

author: Arrowstar

repo: ksptot


text:
function [y] = favar_transx(x,tcode)
%% Transform  x
% Stock & Watson (2016)
%   Return Series with same dimension and corresponding dates
%   Missing values where not calculated
%   -- Tcodes:
%            1 Level
%            2 First Difference
%            3 Second Difference
%            4 Log-Level
%            5 Log-First-Difference
%            6 Log-Second-Difference
% 		     7 First difference of percent change: (x(t)/x(t-1)-1)-(x(t-1)/x(t-2)-1), FRED-MD data (McCracken & Ng)
%

small=1.0e-06;

n=size(x,1);
y=NaN*zeros(n,1);

 if tcode==1
  y=x;

 elseif tcode==2
  y(2:n)=x(2:n)-x(1:n-1);

 elseif tcode==3
  y(3:n)=x(3:n)-2*x(2:n-1)+x(1:n-2);

 elseif tcode==4
  if min(x) < small
   y=NaN;
  else
   x=log(x);
   y=x;
  end

 elseif tcode==5
  if min(x) < small
   y = NaN;
  else
   x=log(x);
   y(2:n)=x(2:n)-x(1:n-1);
  end

 elseif tcode==6
  if min(x) < small
   y = NaN;
  else
   x=log(x);
   y(3:n)=x(3:n)-2*x(2:n-1)+x(1:n-2);
  end

  elseif tcode==7
  if min(x) < small
   y = NaN;
  else
   y1(2:n)=(x(2:n)-x(1:n-1))./x(1:n-1);
   y(3:n)=y1(3:n)-y1(2:n-1);
  end
 else
  y = NaN;
 end

author: european-central-bank

repo: BEAR-toolbox


text:
function [optimState,t_func] = initdesign_vbmc(optimState,Ns,funwrapper,t_func,options)
%INITDESIGN_VBMC Initial sample design (provided or random box).

x0 = optimState.Cache.X_orig;
[N0,D] = size(x0);

if N0 <= Ns
    Xs = x0;
    ys = optimState.Cache.y_orig;
    if N0 < Ns
        switch lower(options.InitDesign)
            case 'plausible'
                % Uniform random samples in the plausible box (in transformed space)
                Xrnd = bsxfun(@plus,bsxfun(@times,rand(Ns-N0,D),optimState.PUB-optimState.PLB),optimState.PLB);
            case 'narrow'
                xstart = warpvars_vbmc(x0(1,:),'dir',optimState.trinfo);
                Xrnd = bsxfun(@plus,bsxfun(@times,rand(Ns-N0,D)-0.5,0.1*(optimState.PUB-optimState.PLB)),xstart);
                Xrnd = bsxfun(@min,bsxfun(@max,Xrnd,optimState.PLB),optimState.PUB);
            otherwise
                error('Unknown initial design for VBMC.');
        end
        Xrnd = warpvars_vbmc(Xrnd,'inv',optimState.trinfo);  % Convert back to original space
        Xs = [Xs; Xrnd];
        ys = [ys; NaN(Ns-N0,1)];
    end
    idx_remove = true(N0,1);

elseif N0 > Ns
    % Cluster starting points
    kmeans_options = struct('Display','off','Method',2,'Preprocessing','whiten');
    idx = fastkmeans(x0,Ns,kmeans_options);

    % From each cluster, take points with higher density in original space
    Xs = NaN(Ns,D); ys = NaN(Ns,1); idx_remove = false(N0,1);
    for iK = 1:Ns
        idxK = find(idx == iK);
        xx = optimState.Cache.X_orig(idxK,:);
        yy = optimState.Cache.y_orig(idxK);
        [~,idx_y] = max(yy);
        Xs(iK,:) = xx(idx_y,:);
        ys(iK) = yy(idx_y);
        idx_remove(idxK(idx_y)) = true;
    end
end
% Remove points from starting cache
optimState.Cache.X_orig(idx_remove,:) = [];
optimState.Cache.y_orig(idx_remove) = [];

Xs = warpvars_vbmc(Xs,'d',optimState.trinfo);

for is = 1:Ns
    timer_func = tic;
    if isnan(ys(is))    % Function value is not available
        [~,optimState] = funlogger_vbmc(funwrapper,Xs(is,:),optimState,'iter');
    else
        [~,optimState] = funlogger_vbmc(funwrapper,Xs(is,:),optimState,'add',ys(is));
    end
    t_func = t_func + toc(timer_func);
end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function [P, Q] = invMyForm(A, B, C)
    % return inv(myForm(A, B, C))
    % -----------------------------------------------
    % Author: Tiep Vu, thv102@psu.edu, 4/8/2016
    %         (http://www.personal.psu.edu/thv102/)
    % -----------------------------------------------
    if nargin == 0
        d = 10;
        C = 10;
        A = rand(d, d);
        B = rand(d, d);
        M = myForm(A, B, C);
    end
    %%
    invA = inv(A);
    invB = inv(B);
    % Q = myForm(invA, -invA*inv(invB + C*invA)*invA, C);
    P = invA;
    Q = -invA*inv(invB + C*invA)*invA;
    % toc;
    if nargin == 0
        M1 = inv(M);
        M2 = myForm(P, Q, C);
        norm(M1 - M2)
        pause
    end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function env  = mrvGetEvironment()
% function env  = mrvGetEvironment()
%
% Copyright Stanford team, mrVista, 2011
%
% FP, 6/30/2011
%
% See also mrvValidate.m and mrvValidateAll.m

% save the matlab version being used
env.matlabVer = ver;

% save the architecture
env.computer = computer;

% vistasoft revision number
env.vistaDir = mrvRootPath;
tmp = gitInfo(vistaRootPath);
env.gitOrigin = tmp.origin;
env.gitchecksum = tmp.checksum;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef homochoricPlot < axisAnglePlot

  methods

    function oP = homochoricPlot(varargin)
      % create a 3d Euler angle plot

      oP = oP@axisAnglePlot(varargin{:});

     end

     function [x,y,z] = project(oP,ori,varargin)

      if ~check_option(varargin,'noBoundaryCheck')
        switch oP.fRMode
          case 'project2FundamentalRegion'
            ori = project2FundamentalRegion(ori);
          case 'restrict2FundamentalRegion'
            ori(~oP.oR.checkInside) = NaN;
        end
      end

      [x,y,z] = double(homochoric(ori));

    end

    function ori = iproject(oP,x,y,z,varargin)
      ori = orientation.id;
    end

    function ori = makeGrid(oP,varargin)

      [ori,S2G,omega] = makeGrid@axisAnglePlot(oP,varargin{:});

      [oP.plotGrid.x,oP.plotGrid.y,oP.plotGrid.z] = ...
        double( S2G .* (3./4 * (omega - sin(omega))).^(1/3));

    end
  end
end

author: mtex-toolbox

repo: mtex

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


text:
function symObj = symmetrise(obj)
% symmetrise embedding
%
% Syntax
%
%   symE = symmetrise(e)
%
% Input
%  e - @embedding
%
% Output
%  symE - @embedding invariant with respect to all symmetry elements
%

symObj = mean( rotate_outer(obj, obj.CS.properGroup.rot),1);
symObj = reshape(symObj, size(obj));

end


text:
% The COBRAToolbox: testConvertHypergraphToBipartiteGraph.m
%
% Purpose:
%     - testconvertHypergraphToBipartiteGraph tests the convertHypergraphToBipartiteGraph
%     function and its different methods
%
% Note:
%      - test can be extended to test the performance of ConvertHypergraphToBipartiteGraph
%        with B2 = hypergraph2bipartgraph(S~=0);
%
% Author:
%     - original file: Marouen BEN GUEBILA - 10/02/2017
%     - integration of test to CI: Laurent Heirendt - February 2017

global CBTDIR

% save the current path
currentDir = pwd;

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

for flag = 1:2
    if flag == 1
        S = [-1  0  1 ;
              1 -1  0 ;
              0  1 -1 ];
    elseif flag == 2
        S = [-1; 1; 1; 0];
    end

    [A, B1] = convertHypergraphToBipartiteGraph(S);

    fprintf('\nB1\n');
    disp(full(B1));
    fprintf('\nA\n');
    disp(full(A));

    %Compute the strongly connected components of a graph
    [sci, sizes] = scomponents(A);

    if flag == 1
        assert(isequal(sci, [1; 1; 1; 1; 1; 1]));
        assert(sizes == 6)
    elseif flag == 2
        assert(isequal(sci, [1; 1; 1; 2; 1;]));
        assert(isequal(sizes, [4; 1]))
    end
end

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

for printLevel = 0:1
    [A, B] = convertHypergraphToBipartiteGraph(model.S, printLevel);

    % compare test data and results
    assert(isequal(A, testA))
    assert(isequal(B, testB))
end

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox