text:
function mesh = ft_defacemesh(cfg, mesh)

% FT_DEFACEMESH allows you to de-identify a scalp surface mesh by erasing specific
% regions, such as the face and ears. The interactive graphical user interface allows
% you to position a box over the anatomical data inside which all vertices will be
% removed. You might have to call this function multiple times when both face and
% ears need to be removed. Following defacing, you should check the result with
% FT_PLOT_MESH.
%
% Use as
%   mesh = ft_defacevolume(cfg, mesh)
%
% The configuration can contain the following options
%   cfg.translate  = initial position of the center of the box (default = [0 0 0])
%   cfg.scale      = initial size of the box along each dimension (default is automatic)
%   cfg.rotate     = initial rotation of the box (default = [0 0 0])
%   cfg.selection  = which vertices to keep, can be 'inside' or 'outside' (default = 'outside')
%
% See also FT_ANONYMIZEDATA, FT_DEFACEVCOLUME, FT_ANALYSISPIPELINE, FT_PLOT_MESH

% Copyright (C) 2015-2022, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar    mesh
ft_preamble provenance mesh

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% the actual work is done by FT_DEFACEVOLUME
tmpcfg = cfg;
tmpcfg.showcallinfo = 'no';
tmpcfg.method = 'interactive';
mesh = ft_defacevolume(tmpcfg, mesh);
% remember the output rotate, scale and translate
rotate    = mesh.cfg.rotate;
scale     = mesh.cfg.scale;
translate = mesh.cfg.translate;
% restore provenance information and put back cfg.callinfo
tmpcallinfo = cfg.showcallinfo;
[cfg, mesh] = rollback_provenance(cfg, mesh);
% store these in the output configuration
cfg.showcallinfo = tmpcallinfo;
cfg.rotate       = rotate;
cfg.scale        = scale;
cfg.translate    = translate;

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous mesh
ft_postamble provenance mesh
ft_postamble history mesh
ft_postamble savevar mesh

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% THE QUATERNION DIFFERENCE
function [dq] = qDifference(q1,q2)
% This function gets the quaternion that will rotate from the
% q1 orientation to the q2 orientation [q1->q2].

assert(size(q1,1) == 4,'The initial quaternion (q1) must be a 4x1 column vector')
assert(size(q1,1) == 4,'The terminal quaternion (q2) must be a 4x1 column vector')

% Get the quaternion inverse
[q1_inv] = qInverse(q1);
% Calculate difference
dq = qMultiply(q1_inv,q2);
% Re-normalise
dq = unit(dq);
end

repo: OpenMAS


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

% parameters
dataset = 'thun';
granularity = 1;
appliance_name = 'Laptop';
households = [2];
dates = ['2012-07-01'];
start_date = '2012-07-01';
end_date = '2012-07-01';
consumptionArray = zeros((24*60*60)/granularity,5);

% obtain consumption data
appliance_id =  getApplianceID(appliance_name);
houses = intersect(findHouseholds(appliance_id), households);
numOfDaysPerHousehold = zeros(length(houses), 1);
counter = 0;
for house = houses
    consumption = read_plug_data(dataset, house, appliance_id, dates, granularity );
    counter = counter + 1;
    consumptionArray(1:length(consumption),counter) = consumption;
end

% plot consumption data
f = figure();
startDate = datenum(strcat(start_date, ' 00:00:00'));
endDate = datenum(strcat(end_date, ' 23:59:59'));
timeTics = linspace(startDate,endDate,86400/granularity);
plot(timeTics, consumptionArray);
datetickzoom('x', 'HH:MM');
ylabel('Power [W]');
xlabel('Time of day');
%legend('Household 1', 'Household 2', 'Household 4', 'Household 5', 'Household 6');
grid on;
%title(sprintf('Power Consumption of %s', appliance_name));
axis tight;

author: beckel

repo: nilm-eval

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


text:
function test_ft_denoise_pca

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

fs = 500;
nchan = 32;
start_time = -1; %s
end_time = 2.5; %s
nsamples = (end_time - start_time) * fs + 1;

data = [];
data.time{1} = linspace(start_time, end_time, nsamples);
data.trial{1} = randn(nchan,nsamples);
data.label = cellstr(num2str((1:nchan).'));

cfg = [];
cfg.refchannel = 'all';
cfg.channel = 'all';
dataout = ft_denoise_pca(cfg,data);


text:
function hrf = glm_hrf(params, tSeries, stim)
%
% hrf = glm_hrf(params, <tSeries>, <stim>);
%
% Returns the hemodynamic impulse response function to use for GLM
% analyses, given the current event analysis settings. The hrf will
% always be returned in units of MR frames.
%
% 'tSeries' and 'stim' are optional input arguments, which are only
% needed if params.glmHRF==1: compute the HRF from the data and
% selected conditions in params.snrConds.
%
% ras, 01/2007.
verbose = prefsVerboseCheck;

if ischar(params.glmHRF)
    hrfPath = fullfile(hrfDir, params.glmHRF);
    load(hrfPath, 'hrf', 'timeWindow', 'tr');
%     if verbose, fprintf('Loading HRF from file %s\n', hrfPath); end
    return
end

% for the pre-defined (Boynton, SPM, Dale&Buckner) HRFs, we'll need
% to get the time window sampled in units of MR frames, excluding pre-
% onset periods:
if ismember(params.glmHRF, [2 3 4])
    tr = params.framePeriod;
    maxT = max(params.timeWindow);          % in secs
    t = 0:tr:maxT;                          % also in secs
end

switch params.glmHRF
    case 1,
        % estimate hrf as mean time course to all stim
        if isempty(tSeries)
            errmsg = ['Need tSeries to estimate mean trial response ' ...
                '(glmHRF option 1)'];
            error(errmsg);
        end

%         if verbose, disp('Estimating mean trial response for GLM'); end
        tSeries = squeeze(mean(tSeries,2)); % mean across voxels
        anal = er_chopTSeries2(tSeries(:),stim,params);
        postStim = find(anal.timeWindow>=0);
        hrf = anal.meanTcs(postStim,anal.condNums>0);
        hrf = hrf(:,params.snrConds);
        if size(hrf,2) > 1
            hrf = nanmeanDims(hrf,2);
        end

    case 2,  % boynton gamma function
        % create sub-params struct
%         if verbose, disp('Boynton & Heeger gamma HRF'); end
        [p params] = glm_getHrfParams(params);
        n = params.glmHRF_params(1);
        tau = params.glmHRF_params(2);
        delay = params.glmHRF_params(3);
        hrf = boyntonHIRF(t, n, tau, delay);

    case 3,  % spm difference-of-gammas
%         if verbose, disp('SPM difference of gamma HRF'); end
        [p params] = glm_getHrfParams(params);
        hrf = spm_hrf(tr, p);

    case 4, % dale & buckner HRF
%         if verbose, disp('Dale and Buckner 2000 HRF: t^2*exp(-t)'); end
        [p params] = glm_getHrfParams(params);
        hrf = fmri_hemodyn(t, p(1), p(2));

    case 5, % saved HRF from file
        % choose from file
%         if verbose, disp('HRF selected from file \n'); end
        [f p] = myUiGetFile(hrfDir, '*.mat', 'Select a saved HRF File...');
        if f==0 % user canceled
            disp('User canceled -- setting HRF to Boynton Gamma')
            params.glmHRF = 2;
        else
            params.glmHRF = f(1:end-4);
        end
        hrf = glm_hrf(params);
        % TO DO: deal w/ situations in which the HRF was
        % saved using diff't TR, time window

    otherwise, error('Invalid specification for params.glmHRF.')
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [err,Eq] = Plane_Equation (Triplets, verbose)
%
%   [err,Eq] = Plane_Equation (Triplets, [verbose])
%
%Purpose:
%   Determine the equation of the plane passing through three points
%
%
%Input Parameters:
%   Triplets is an Nx1 vector of strucutres. Each structure defines a plane
%    .XYZ : is a 3x3 matrix containing the XYZ of each of the three points
%           Each row is a point. ie: XYZ = [0 0 0; 1 0 0; 1 1 1]; is for the
%           three points (0 0 0), (1 0 0) and (1 1 1).
%           If the three points are colinear, Eq = [0 0 0 0]
%
%   verbose (0/1), default is 1
%
%Output Parameters:
%   err : 0 No Problem
%       : 1 Mucho Problems
%
%   Eq is a Nx4 matrix containing the equation of the plane containing each
%       triplet in  Triplets. The plane passing by triplet i is speicifed in
%       Eq(i,:) the plane would be Eq(i,1)x + Eq(i,2)y + Eq(i,3)z + Eq(i,4) = 0
%
%More Info :
%
%   see also ShowPlane
%   try
%       Triplets(1).XYZ = [0 0 0; 1 0 0; 1 1 0];
%       Triplets(2).XYZ = [0 0 0; 1 0 0; 1 1 1];
%       Triplets(3).XYZ = [0 5 0; 1 5 0; 1 1 1];
%
%       [err,Eq] = Plane_Equation (Triplets);
%       [err,PatchHandles] = ShowPlane (Eq); view(3)
%
%     Author : Ziad Saad
%     Date : Thu Oct 22 16:09:56 CDT 1998

%Define the function name for easy referencing
FuncName = 'Plane_Equation';

%initailize return variables
err = 1;

if (nargin == 1), 	verbose = 1;	end

if (is_row(Triplets) == -1),	err = ErrEval(FuncName,'Err_Triplets must be an Nx1 vector');	return;	end

Triplets = Triplets(:);
Nplanes = size(Triplets,1);
Eq = zeros(Nplanes,4); %allocate

for (i=1:1:Nplanes),
	%Form the equation of the plane
	x1 = Triplets(i).XYZ(1,1); y1 = Triplets(i).XYZ(1,2); z1 = Triplets(i).XYZ(1,3);
	x2 = Triplets(i).XYZ(2,1); y2 = Triplets(i).XYZ(2,2); z2 = Triplets(i).XYZ(2,3);
	x3 = Triplets(i).XYZ(3,1); y3 = Triplets(i).XYZ(3,2); z3 = Triplets(i).XYZ(3,3);
	Eq(i,1) = y1.*(z2-z3) + y2.*(z3-z1) + y3.*(z1-z2);
	Eq(i,2) = z1.*(x2-x3) + z2.*(x3-x1) + z3.*(x1-x2);
	Eq(i,3) = x1.*(y2-y3) + x2.*(y3-y1) + x3.*(y1-y2);
	Eq(i,4) = -x1.*(y2.*z3 - y3.*z2) - x2.*(y3.*z1 - y1.*z3) - x3.*(y1.*z2 - y2.*z1);

	if (verbose),
		if (~rem (i,200)),
			fprintf (1,'[%g/%g]\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b',i,Nplanes);
		end

	end
end

err = 0;
return;


text:
function isRenamed = db_rename_protocol(newName)
% DB_RENAME_PROTOCOL: Rename current protocol

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

global GlobalData;
isRenamed = 0;

% Get protocol to rename
iProtocol = bst_get('iProtocol');
ProtocolInfo = bst_get('ProtocolInfo');
if isempty(iProtocol) || (iProtocol < 1)
    bst_error('No protocol selected.', 'Rename protocol', 0);
    return;
end

% Get new protocol name
if (nargin < 1) || isempty(newName)
    newName = java_dialog('input', 'Enter new protocol name:', 'Rename protocol', [], ProtocolInfo.Comment);
end
if isempty(newName) || isequal(newName, ProtocolInfo.Comment)
    return;
elseif ~isequal(newName, file_standardize(newName))
    bst_error('Invalid protocol name.', 'Rename protocol', 0);
    return;
elseif any(strcmpi(newName, {GlobalData.DataBase.ProtocolInfo.Comment}))
    bst_error(['A protocol named "' newName '" already exists in your database.'], 'Rename protocol', 0);
    return;
end

% New protocol folder
[oldPathAnat, anatName] = bst_fileparts(ProtocolInfo.SUBJECTS);
[oldPathData, dataName] = bst_fileparts(ProtocolInfo.STUDIES);
[basePath, oldName] = bst_fileparts(oldPathData);
newPath = bst_fullfile(basePath, newName);

% Check that anat and data folders are in the same folder
if ~isequal(oldPathData, oldPathAnat)
    bst_error('Only protocols with anat and data folders in the same folder can be renamed.', 'Rename protocol', 0);
    return;
end
% Check if folder already exists
if file_exist(newPath)
    bst_error(['This folder already exists:' 10 newPath], 'Rename protocol', 0);
    return;
end

% Rename protocol folder
isMoved = file_move(oldPathData, newPath);
if ~isMoved
    bst_error(['Could not rename protocol folder.' 10 'Source: ' oldPathData 10 'Destination: ' newPath], 'Rename protocol', 0);
    return;
end

% Change protocol folders
ProtocolInfo.Comment = newName;
ProtocolInfo.SUBJECTS = bst_fullfile(newPath, anatName);
ProtocolInfo.STUDIES = bst_fullfile(newPath, dataName);
db_edit_protocol('edit', ProtocolInfo, iProtocol);

isRenamed = 1;

% Refresh protocols list
gui_brainstorm('UpdateProtocolsList');

author: brainstorm-tools


text:
function [V,openVariables] = branchVolume(p,openVariables)

if nargin == 1
    d = p.ub(p.branch_variables)-p.lb(p.branch_variables);
    openVariables = find(d~=0);
    d = d(openVariables);
    if any(d<0)
        % This is infeasible!
        V = 0;
    else
        V = geomean(d);
    end
else
    if ~p.feasible
        V = 0;
    else
        d = p.ub(p.branch_variables)-p.lb(p.branch_variables);
        d = d(openVariables);
        if any(d<0)
            % This has been propagated to infeasibility
            V = 0;
        else
            V = geomean(d);
        end
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function ns2de_gnuplot ( header, n, x, y, u, v, s )

%*****************************************************************************80
%
%% NS2DE_GNUPLOT writes the Navier-Stokes velocity field to files for GNUPLOT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string HEADER, a header to be used to name the files.
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), Y(N), the coordinates of the evaluation points.
%
%    Input, real U(N), V(N), the velocity components.
%
%    Input, real S, a scale factor for the velocity vectors.
%

%
%  Write the data file.
%
  data_filename = strcat ( header, '_data.txt' );

  data_unit = fopen ( data_filename, 'wt' );

  for i = 1 : n
    fprintf ( data_unit, '  %10.4f  %10.4f  %10.4f  %10.4f\n', ...
      x(i), y(i), s * u(i), s * v(i) );
  end

  fclose ( data_unit );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data written to "%s".\n', data_filename );
%
%  Write the command file.
%
  command_filename = strcat ( header, '_commands.txt' );
  plot_filename = strcat ( header, '.png' );

  command_unit = fopen ( command_filename, 'wt' );

  fprintf ( command_unit, '#  %s\n', command_filename );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, 'set term png\n' );
  fprintf ( command_unit, 'set output "%s"\n', plot_filename );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, '#  Add titles and labels.\n' );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, 'set xlabel "<--- X --->"\n' );
  fprintf ( command_unit, 'set ylabel "<--- Y --->"\n' );
  fprintf ( command_unit, 'set title "Navier-Stokes velocity field"\n' );
  fprintf ( command_unit, 'unset key\n' );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, '#  Add grid lines.\n' );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, 'set grid\n' );
  fprintf ( command_unit, 'set size ratio -1\n' );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, '#  Timestamp the plot.\n' );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, 'set timestamp\n' );
  fprintf ( command_unit, 'plot "%s" using 1:2:3:4 with vectors \\\n', ...
    data_filename );
  fprintf ( command_unit, '  head filled lt 2 linecolor rgb "blue"\n' );
  fprintf ( command_unit, 'quit\n' );

  fclose ( command_unit );

  fprintf ( 1, '  Commands written to "%s".\n', command_filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [irf_record_allt,favar]=tvbvarirf2(beta_gibbs,D_record,It,Bu,IRFperiods,n,m,p,k,T,favar)

% Favar preliminaries
if favar.FAVAR==1
    signresX_index=[];
    npltX=favar.npltX;
    if npltX>0
        favar_FAVAR=1;
        plotX_index=favar.plotX_index;
    else
        favar_FAVAR=0;
        plotX_index=[];
    end
    %relevant loadings
    relLindex=[signresX_index;plotX_index];
    % reshape the L gibbs draws
    Lgibbs=reshape(favar.L_gibbs,size(favar.L,1),size(favar.L,2),It-Bu);
    Lgibbs=Lgibbs(relLindex,:,:);
else
    npltX=0;
    favar_FAVAR=0;
    Lgibbs=NaN;
end

% create the cell aray that will store the values from the simulations
irf_record_allt=cell(n,n);
favar.IRF.favar_irf_record_allt=cell(npltX,n);

% loop over sample periods
for tt=1:T
    % loop over iterations
    for jj=1:It-Bu
        % draw beta
        beta=beta_gibbs{tt}(:,jj);
        D=reshape(D_record(:,It-Bu),n,n);
        if favar_FAVAR==1
            Lg=squeeze(Lgibbs(:,:,It-Bu)); % this is different for IRFt4, we need to record the index there and adjust the case
        end
        [~,ortirfmatrix]=bear.irfsim(beta,D,n,m,p,k,IRFperiods);

        % if we have FAVAR restrictions we scale the ortirfmatrix from the previous step
            if favar_FAVAR==1
                favar_ortirfmatrix=[];
                % scale with loading
                for uu=1:npltX %over variables in X that we choose to plot
                    for mm=1:IRFperiods
                        for ll=1:n % over shocks
                            favar_ortirfmatrix(uu,ll,mm)=Lg(uu,:)*ortirfmatrix(:,ll,mm);
                        end
                    end
                end
            end

            % save
            for kk=1:n
                for ll=1:n
                    for mm=1:IRFperiods
                        irf_record_allt{kk,ll}(jj,mm,tt)=ortirfmatrix(kk,ll,mm);
                    end
                end
            end

            if favar_FAVAR==1
                for uu=1:npltX % loop over variables
                    for ll=1:n % loop over shocks
                        for mm=1:IRFperiods
                        favar.IRF.favar_irf_record_allt{uu,ll}(jj,mm,tt)=favar_ortirfmatrix(uu,ll,mm);
                        end
                    end
                end
            end
     end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function blas2_test01 ( )

%*****************************************************************************80
%
%% BLAS2_TEST01 tests DGEMV.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 February 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BLAS2_TEST01\n' );
  fprintf ( 1, '  For a general matrix A,\n' );
  fprintf ( 1, '  DGEMV computes y := alpha * A * x + beta * y\n' );
  fprintf ( 1, '  or             y := alpha * A'' * x + beta * y.\n' );
%
%  y = alpha * A * x + beta * y
%
  trans = 'N';
  m = 5;
  n = 4;
  alpha = 2.0;
  lda = m;
  a = r8mat_test ( trans, lda, m, n );
  x = zeros ( n, 1 );
  for i = 1 : n
    x(i) = i;
  end
  incx = 1;
  beta = 3.0;
  y = zeros ( m, 1 );
  for i = 1 : m
    y(i) = 10 * i;
  end
  incy = 1;

  r8mat_print ( m, n, a, '  Matrix A:' );
  r8vec_print ( n, x, '  Vector X:' );
  r8vec_print ( m, y, '  Vector Y:' );

  y = dgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy );

  r8vec_print ( m, y, '  Result Y = alpha * A  * x + beta * y' );
%
%  y = alpha * A' * x + beta * y
%
  trans = 'T';
  m = 5;
  n = 4;
  alpha = 2.0;
  lda = m;
  a = r8mat_test ( trans, lda, n, m );
  x = zeros ( m, 1 );
  for i = 1 : m
    x(i) = i;
  end
  incx = 1;
  beta = 3.0;
  y = zeros ( n, 1 );
  for i = 1 : n
    y(i) = 10 * i;
  end
  incy = 1;

  r8mat_print ( m, n, a, '  Matrix A:' );
  r8vec_print ( m, x, '  Vector X:' );
  r8vec_print ( n, y, '  Vector Y:' );

  y = dgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy );

  r8vec_print ( n, y, '  Result Y = alpha * A'' * x + beta * y' );

  return
end


text:
A=randn(5000,500);
A=A'*A;
b=ones(500,1);
x0=b;
tol=1e-4;
itermax=0.5*length(b);

tic
for ii = 1:20
x1 = mexConjGrad(A,b,x0,tol,itermax);
end
t=toc;
fprintf('mex-file time: %fs\n',t);

tic
for ii = 1:20
x2 = pcg(A,b);
end
t=toc;
fprintf('Matlab time: %fs\n',t);
sum((x1(:)-x2(:)).^2)

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function g = sigmoid(z)

% Computes sigmoid of z

g = 1.0 ./ (1.0 + exp(-z));

author: Sable

repo: mcbench-benchmarks


text:
classdef S2VectorFieldTri < S2VectorField
% a class represeneting a function on the sphere

  properties
    tri       % S2Triangulation
    values = vector3d  % function values
  end

  properties (Dependent = true)
    vertices
    antipodal
  end

  methods

    function sVF = S2VectorFieldTri(nodes,values)
      % initialize a spherical vector field

      if nargin == 0, return; end

      if isa(nodes,'function_handle')
        n = equispacedS2Grid('resolution',1.5*degree);
        values = nodes(n);
        nodes = n;
      end

      if isa(nodes,'S2Triangulation')
        sVF.tri = nodes;
      else
        sVF.tri = S2Triangulation(nodes);
      end

      sVF.values = values;

    end

    function v = get.vertices(S2F)
      v = S2F.tri.vertices;
    end

    function v = get.antipodal(S2F)
      v = S2F.tri.antipodal;
    end

    function S2F = set.vertices(S2F,v)
      if ~isempty(S2F.values), S2F.values = S2F.eval(v); end
      S2F.tri.vertices = v;
      S2F.tri.update;
    end

  end

end

author: mtex-toolbox

repo: mtex

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


text:
function varargout = subsref(varargin)
%SUBSREF (overloaded)

% Stupid first slice call (supported by MATLAB)
% x = sdpvar(2);x(1,:,:)
Y = varargin{2};
if length(Y)==1
    if  length(Y.subs) > 2 && isequal(Y.type,'()') && ~( isa(Y(1).subs{1},'constraint') || isa(Y(1).subs{1},'lmi') )
        i = 3;
        ok = 1;
        while ok && (i <= length(Y.subs))
            ok = ok && (isequal(Y.subs{i},1) || isequal(Y.subs{i},':'));
            i = i + 1;
        end
        if ok
            Y.subs = {Y.subs{1:2}};
        else
            error('??? Index exceeds matrix dimensions.');
        end
    end
end

X = varargin{1};

try
    switch Y(1).type
        case '()'
            if  isa(Y(1).subs{1},'constraint') || isa(Y(1).subs{1},'lmi')
                z = Y(1).subs{1};
                for i = 1:length(Y(1).subs)
                    z = [z, Y(1).subs{i}];
                end
                % z == 0 => b+Ax == 0 => x = -A\b
                z = sdpvar(z);
                B = getbase(z);
                A = B(:,2:end);
                b = B(:,1);
                x0 = -A\b;
                x = recover(getvariables(z));
                varargout{1} = replace(X,x,x0);
                return
            end
            % Check for simple cases to speed things up (yes, ugly but we all want speed don't we!)
            switch size(Y(1).subs,2)
                case 1
                    y = subsref1d(X,Y(1).subs{1},Y);
                case 2
                    y = subsref2d(X,Y.subs{1},Y(1).subs{2},Y);
                otherwise
                    if all( [Y(1).subs{3:end}]==1)
                        y = subsref2d(X,Y.subs{1},Y(1).subs{2},Y);
                    else
                        error('Indexation error.');
                    end
            end
        case '{}'
            varargout{nargout} = [];

            % it could be the case that we have an extended variable
            % This is a bit tricky, so we do the best we can; assume that
            % we want to replace the internal argument wih the new
            % expression
            OldArgument = recover(depends(X));
            vars = getvariables(X);
            mpt_solution = 1;
            if all(ismembc(vars,yalmip('extvariables')))
                for i = 1:length(X)
                    nonlinearModel = yalmip('extstruct',vars);
                    if isequal(nonlinearModel{1}.fcn,'pwa_yalmip') | isequal(nonlinearModel{1}.fcn,'pwq_yalmip')
                    else
                        mpt_solution = 0;
                    end
                end
                if mpt_solution
                    assign(nonlinearModel{1}.arg{2},Y(1).subs{:});
                    XX = value(X);
                    varargout{1} = value(X);
                    return
                end
            end
            vars = getvariables(X);
            if (length(vars) == 1) & ismembc(vars,yalmip('extvariables'))
                nonlinearModel = yalmip('extstruct',vars);
                OldArgument = [];
                for i = 1:length(nonlinearModel.arg)
                    if isa(nonlinearModel.arg{i},'sdpvar')
                        OldArgument = [OldArgument;  nonlinearModel.arg{i}];
                    end
                end
                if isnumeric([Y.subs{:}])
                    assign(reshape(OldArgument,[],1),reshape([Y(1).subs{:}],[],1));
                    varargout{1} = value(X);
                    return
                end
            end
            y = replace(X,OldArgument,[Y(1).subs{:}]);
            if isnumeric(y)
                varargout{1} = y;
                return
            end

        case '.'
            switch Y(1).subs
                case {'minimize','maximize'}
                    options = [];
                    constraints = [];
                    objective = varargin{1};
                    opsargs = {};
                    if length(Y)==2
                        if isequal(Y(2).type,'()')
                            for i = 1:length(Y(2).subs)
                                switch class(Y(2).subs{i})
                                    case {'lmi','constraint'}
                                        constraints = [constraints, Y(2).subs{i}];
                                    case 'struct'
                                        options = Y(2).subs{i};
                                    case {'double','char','gem','sgem'}
                                        opsargs{end+1} = Y(2).subs{i};
                                    otherwise
                                        error('Argument to minimize should be constraints or options');
                                end
                            end
                        else
                            error(['What do you mean with ' Y(2).type '?']);
                        end
                    end
                    if length(opsargs)>0
                        if isempty(options)
                           options = sdpsettings(opsargs{:});
                        else
                            options = sdpsettings(options,opsargs{:});
                        end
                    end
                    if isequal(Y(1).subs,'minimize')
                        sol = solvesdp(constraints,objective,options);
                    else
                        sol = solvesdp(constraints,-objective,options);
                    end
                    varargout{1} = varargin{1};
                    varargout{2} = sol;
                    return
                case 'derivative'
                    try
                        m = model(varargin{1});
                        varargout{1} = m{1}.derivative;
                    catch
                        varargout{1} = 1;
                    end
                    return
                otherwise
                    error(['Indexation  ''' Y.type Y.subs ''' not supported']) ;
            end
        otherwise
            error(['Indexation with ''' Y.type ''' not supported']) ;
    end
catch
    error(lasterr)
end
if isempty(y.lmi_variables)
    y = full(reshape(y.basis(:,1),y.dim(1),y.dim(2)));
else
    % Reset info about conic terms
    y.conicinfo = [0 0];
    y.extra.createTime = definecreationtime;
end
varargout{1} = y;

function X = subsref1d(X,ind1,Y)

% Get old and new size
n = X.dim(1);
m = X.dim(2);

% Convert to linear indecicies
if islogical(ind1)
    ind1 = double(find(ind1));
elseif ischar(ind1)
    X.dim(1) = n*m;
    X.dim(2) = 1;
    return;
elseif ~isnumeric(ind1)
    X = milpsubsref(X,Y);
    return
end

% Detect X(scalar)
if length(ind1) == 1 & ind1 <= n*m

    Z = X.basis(ind1,:);
    nnew = 1;
    mnew = 1;

else

    % What would the size be for a double
    dummy = reshape(X.basis(:,1),n,m);
    dummy = dummy(ind1);
    nnew = size(dummy,1);
    mnew = size(dummy,2);
    [nx,mx] = size(X.basis);

    if length(ind1) > 1
        try
            % row-based subsref for sparse objects can be very slow and
            % take a lot of memory, transposing a big sparse array as in
            %   Z = X.basis.';
            %   Z = Z(:,ind1);
            %   Z = Z.';
            % can also be quite slow for very large matrices, so we use
            % some custom code
            [ix,jx,sx] = find(X.basis);
            if length(ix) == length(unique(ix))
                % We never have two elements on the same row of X.basis
                [isNnz,loc] = ismember(ind1(:),ix);
                sel = loc(isNnz);
                ix = find(isNnz);
                jx = jx(sel);
                sx = sx(sel);
                Z = sparse(ix,jx,sx,numel(ind1),mx);
            else
                % We can have seveal elements on a given row of X.basis
                if numel(ind1) == length(unique(ind1))
                    % Every row of X.basis is selected at most once
                    [keep,loc] = ismember(ix,ind1(:));
                    ix = loc(keep);
                    jx = jx(keep);
                    sx = sx(keep);
                    Z = sparse(ix,jx,sx,numel(ind1),mx);
                else
                    % Complicated case, we use more generica code
                    Z = sparse(1:numel(ind1),ind1,1,numel(ind1),size(X,1))*X.basis;
                end
            end

            if false
                % For checking purpose
                Z0 = X.basis.';
                Z0 = Z0(:,ind1);
                Z0 = Z0.';
                assert(isequal(Z, Z0));
            end
        catch
            Z = X.basis(ind1,:);
        end
    else
        Z = X.basis(ind1,:);
    end
end

% Find non-zero basematrices
nzZ = find(any(Z(:,2:end),1));
if ~isempty(nzZ)
    X.dim(1) = nnew;
    X.dim(2) = mnew;
    X.lmi_variables =  X.lmi_variables(nzZ);
    X.basis = Z(:,[1 1+nzZ]);
else
    bas = reshape(X.basis(:,1),n,m);
    X.dim(1) = nnew;
    X.dim(2) = mnew;
    X.lmi_variables = [];
    X.basis = reshape(bas(ind1),nnew*mnew,1);
end

function X = subsref2d(X,ind1,ind2,Y)

if isnumeric(ind1)
elseif ischar(ind1)
    ind1 = 1:X.dim(1);
elseif islogical(ind1)
    ind1 = double(find(ind1));
elseif ~isnumeric(ind1)
    X = milpsubsref(X,Y);
    return
end
if isnumeric(ind2)
elseif ischar(ind2)
    ind2 = 1:X.dim(2);
elseif islogical(ind2)
    ind2 = double(find(ind2));
elseif  ~isnumeric(ind2)
    X = milpsubsref(X,Y);
    return
end

n = X.dim(1);
m = X.dim(2);
lind2 = length(ind2);
lind1 = length(ind1);
if lind2 == 1
    ind1_ext = ind1(:);
else
    ind1_ext = kron(ones(lind2,1),ind1(:));
end
if lind1 == 1
    ind2_ext = ind2(:);
else
    ind2_ext = kron(ind2(:),ones(lind1,1));
end

if any(ind1 > n) || any(ind2 > m)
    error('Index exceeds matrix dimensions.');
end

if lind1==1 && lind2==1
    if isequal(X.conicinfo,[-1 0])
        X.basis = [0 1];
        X.lmi_variables = X.lmi_variables(1)+ind1+(ind2-1)*n-1;
        X.dim = [1 1];
        X.conicinfo = [0 0];
        return
    end
end

if prod(size(ind1_ext))==0 | prod(size(ind2_ext))==0
    linear_index = [];
else
    % Speed-up for some bizarre code with loads of indexing of vector
    if m==1 & ind2_ext==1
        linear_index = ind1_ext;
    elseif length(ind2_ext)==1 && length(ind1_ext)==1
        linear_index = ind1_ext + (ind2_ext-1)*n;
    else
        linear_index = sub2ind([n m],ind1_ext,ind2_ext);
    end
end
nnew = length(ind1);
mnew = length(ind2);

% Put all matrices in vectors and extract sub matrix
Z = X.basis(linear_index,:);
% Find non-zero basematrices
%nzZ = find(any(Z(:,2:end),1));
nzZ = find(any(Z,1))-1;
if numel(nzZ)>0
    if nzZ(1)==0
    nzZ = nzZ(2:end);
    end
end
if ~isempty(nzZ)
    X.dim(1) = nnew;
    X.dim(2) = mnew;
    X.lmi_variables =  X.lmi_variables(nzZ);
    X.basis = Z(:,[1 1+nzZ]);
else
    bas = reshape(X.basis(:,1),n,m);
    X.dim(1) = nnew;
    X.dim(2) = mnew;
    X.lmi_variables = [];
    X.basis = reshape(bas(linear_index),nnew*mnew,1);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% CVX: Top-level commands to control CVX.
%   cvx_begin         - Starts a new CVX specification.
%   cvx_begin_set     - Starts a new cvx CVX specification.
%   cvx_cd            - Change current working directory to a CVX subdirectory.
%   cvx_clear         - Clears all active CVX data.
%   cvx_end           - Completes a cvx specification.
%   cvx_end_set       - Completes a cvx set specification.
%   cvx_gp_precision  - Changes the precision of CVX's GP approximation.
%   cvx_pause         - Pauses the processing of CVX models.
%   cvx_power_warning - Controls the CVX warning message for x.^p expressions.
%   cvx_precision     - Controls CVX solver precision.
%   cvx_profile       - CVX-specific profiler control.
%   cvx_quiet         - CVX output control.
%   cvx_solver        - CVX solver selection.
%   cvx_where         - Returns the location of the CVX system.

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [C,fval]=optimPID(G,ctype,idx)
% OPTIMPID  Optimal PID tuning based on integral performance criteria
%
% [C,fval]=optimPID(G,ctype,idx) returns the optimal PID paraters based on
% specified controller type and performance criterion.
%
% Inputs:
%               G: The plant model as an LTI object
%           ctype: Controler type (1 = P, 2* = PI, 3 = PID)
%             idx: Performance criterion
%                  1  - ISE
%                  2  - IAE
%                  3  - ITSE
%                  4* - ITAE
% Outputs:
%               C: Controller transfer function as an LTI object
%            fval: optimal performance criterion
% Example:
%{
G=tf(1,[1 6 11 6 0]);
C1=optimPID(G,3,1);   % PID-Control, ISE index
C2=optimPID(G,3,2);   % PID-Control, IAE index
C3=optimPID(G,3,3);   % PID-Control, ITSE index
C4=optimPID(G,3,4);   % PID-Control, ITAE index
K=znpidtuning(G,3);   % Ziegler-Nichols stability margin tuning
t=0:0.1:30;
y1=step(feedback(C1*G,1),t);
y2=step(feedback(C2*G,1),t);
y3=step(feedback(C3*G,1),t);
y4=step(feedback(C4*G,1),t);
y=step(feedback(G*(K.kc*(1+tf(1,[K.ti 0])+tf([K.td 0],1))),1),t);
plot(t,y1,t,y2,t,y3,t,y4,t,y,'--','Linewidth',2)
legend('ISE','IAE','ITSE','ITAE','Z-N')
grid
%}
% By Yi Cao at Cranfield University on 8th Feb 2008
%

% Check inputs and outputs
error(nargchk(1,3,nargin));
error(nargoutchk(0,3,nargout));
assert(isa(G,'lti'),'G must be an LTI object.');

% default setting
if nargin<3
    idx=4;
end
if nargin<2
    ctype=2;
end
% Initial parameters using stability based tuning
[Gm,Pm,Wcg]=margin(G);
pu=2*pi/Wcg;
ku=Gm;
x=ku/2;
den=1;
if ctype==2
    x=ku/2.2*[1 1.2/pu];
    den=[1 0];
elseif ctype==3
    x=ku*2/pu/1.7*[pu/8 1 2/pu];
    den=[1 0];
end
% closed-loop response of initial tuning to find dt and tend
[y,t]=step(feedback(tf(x,den)*G,1));
% reduce dt by half for possible improvement in response speed
dt=(t(2)-t(1))/2;
% exptend tend twice to ensure closed-loop stability
t=0:dt:t(end)*2;
% redefine cost function to facilitate optimization
cost = @(x) iecost(x,G,den,t,dt,idx);
opt=optimset('display','off','TolX',1e-9,'TolFun',1e-9,'LargeScale','off');
flag=0;
while ~flag % if flag=0 restart optimization from current solution
    [x,fval,flag]=fminunc(cost,x,opt);
end
% the transfer function of optimal PID controller
C=tf(x,den);

function J=iecost(x,G,den,t,dt,idx)
% control error of step response
e=1-step(feedback(G*tf(x,den),1),t);
% performance calculation
switch idx
    case 1  % ISE
        J=e'*e*dt;
    case 2  % IAE
        J=sum(abs(e)*dt);
    case 3  % ITSE
        J=(t.*e'*dt)*e;
    case 4  % ITAE
        J=sum(t'.*abs(e)*dt);
end


text:
% set the dimensions of the two leg links

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

% create the leg links based on DH parameters
%                    theta   d     a  alpha
links(1) = Link([    0       0    0   pi/2 ], 'standard');
links(2) = Link([    0       0    L1   0   ], 'standard');
links(3) = Link([    0       0   -L2   0   ], 'standard');

% now create a robot to represent a single leg
leg = SerialLink(links, 'name', 'leg', 'offset', [pi/2   0  -pi/2]);

% define the key parameters of the gait trajectory, walking in the
% x-direction
xf = 5; xb = -xf;   % forward and backward limits for foot on ground
y = 5;              % distance of foot from body along y-axis
zu = 2; zd = 5;     % height of foot when up and down
% define the rectangular path taken by the foot
segments = [xf y zd; xb y zd; xb y zu; xf y zu] * 0.01;

% build the gait. the points are:
%   1 start of walking stroke
%   2 end of walking stroke
%   3 end of foot raise
%   4 foot raised and forward
%
% The segments times are :
%   1->2  3s
%   2->3  0.5s
%   3->4  1s
%   4->1  0.5ss
%
% A total of 4s, of which 3s is walking and 1s is reset.  At 0.01s sample
% time this is exactly 400 steps long.
%
% We use a finite acceleration time to get a nice smooth path, which means
% that the foot never actually goes through any of these points.  This
% makes setting the initial robot pose and velocity difficult.
%
% Intead we create a longer cyclic path: 1, 2, 3, 4, 1, 2, 3, 4. The
% first 1->2 segment includes the initial ramp up, and the final 3->4
% has the slow down.  However the middle 2->3->4->1 is smooth cyclic
% motion so we "cut it out" and use it.
segments = [segments; segments];
tseg = [3 0.25 0.5 0.25]';
tseg = [1; tseg; tseg];
x = mstraj(segments, [], tseg, segments(1,:), 0.01, 0.1);

% pull out the cycle
xcycle = x(100:500,:);
qcycle = leg.ikine( transl(xcycle), [], [1 1 1 0 0 0], 'pinv' );

% dimensions of the robot's rectangular body, width and height, the legs
% are at each corner.
W = 0.1; L = 0.2;

% a bit of optimization.  We use a lot of plotting options to
% make the animation fast: turn off annotations like wrist axes, ground
% shadow, joint axes, no smooth shading.  Rather than parse the switches
% each cycle we pre-digest them here into a plotopt struct.
% plotopt = leg.plot({'noraise', 'nobase', 'noshadow', ...
%     'nowrist', 'nojaxes'});
% plotopt = leg.plot({'noraise', 'norender', 'nobase', 'noshadow', ...
%     'nowrist', 'nojaxes', 'ortho'});

plotopt = {'noraise', 'nobase', 'noshadow', 'nowrist', 'nojaxes', 'delay', 0};

% create 4 leg robots.  Each is a clone of the leg robot we built above,
% has a unique name, and a base transform to represent it's position
% on the body of the walking robot.
legs(1) = SerialLink(leg, 'name', 'leg1');
legs(2) = SerialLink(leg, 'name', 'leg2', 'base', transl(-L, 0, 0));
legs(3) = SerialLink(leg, 'name', 'leg3', 'base', transl(-L, -W, 0)*trotz(pi));
legs(4) = SerialLink(leg, 'name', 'leg4', 'base', transl(0, -W, 0)*trotz(pi));

% create a fixed size axis for the robot, and set z positive downward
clf; axis([-0.3 0.1 -0.2 0.2 -0.15 0.05]); set(gca,'Zdir', 'reverse')
hold on
% draw the robot's body
patch([0 -L -L 0], [0 0 -W -W], [0 0 0 0], ...
    'FaceColor', 'r', 'FaceAlpha', 0.5)
% instantiate each robot in the axes
for i=1:4
    legs(i).plot(qcycle(1,:), plotopt{:});
end
hold off

% walk!
k = 1;
%A = Animate('walking');
%while 1
for i=1:500
    legs(1).animate( gait(qcycle, k, 0,   0));
    legs(2).animate( gait(qcycle, k, 100, 0));
    legs(3).animate( gait(qcycle, k, 200, 1));
    legs(4).animate( gait(qcycle, k, 300, 1));
    drawnow
    k = k+1;
    %A.add();
end

author: Allopart

repo: rbpf-gmapping


text:
function x = uniform_01_cdf_inv ( cdf )

%*****************************************************************************80
%
%% UNIFORM_01_CDF_INV inverts the Uniform 01 CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%    0.0 <= CDF <= 1.0.
%
%    Output, real X, the corresponding argument.
%
  if ( cdf < 0.0 | 1.0 < cdf )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'UNIFORM_01_CDF_INV - Fatal error!\n' );
    fprintf ( 1, '  CDF < 0 or 1 < CDF.\n' );
    error ( 'UNIFORM_01_CDF_INV - Fatal error!' );
  end

  x = cdf;

  return
end


text:
%% DEMO_stent_hexahedral_sweeping
% Below is a demonstration for:
%
% * Creating a hexahedral mesh for a vascular stent by sweeping allong a
% curve and copying over the segments.

%% Keywords
% * Sweeping, sweepLoft
% * Hexahedral mesh
% * stent, vascular
% * Exporting Abaqus, .inp

%%
clear; close all; clc;

%%
% plot settings

fontSize=25;
markerSize=10;
lineWidth=1;

%% Contol parameters

stentRadius=3; %The outer radius of the stent
numPeriodsWave=10; %The number of periods to use for a sinusoidal modulation
numStepsPeriod=100; %Number of sweeping steps allong a single period segment for sweeping
waveAmplitude=0.9; %Amplitude of the sinusoidal modulation
stentSectionHeight=0.1; %Height of the stent wire
stentSectionWidth=0.1; %Width of the stent wire
numStepsCircumference=(numPeriodsWave*numStepsPeriod)+1; %Number of sweeping steps across curve
overSampleFactorCurve=10; %Oversample factor curve
numSegments=8; %Number of stent segments to stack
segmentAxialOffset=waveAmplitude*2; %Axial offset for stacking stents

%% Build stent section
% The rectangular stent wire section is created here.

V_section=[-stentSectionWidth/2  stentSectionHeight/2 0; ...
            stentSectionWidth/2  stentSectionHeight/2 0; ...
            stentSectionWidth/2 -stentSectionHeight/2 0; ...
           -stentSectionWidth/2 -stentSectionHeight/2 0; ...
            ];

%%
% Visualize stent section
cFigure; hold on;
title('Stent section','fontSize',fontSize);
plotV(V_section,'b.-','lineWidth',lineWidth,'MarkerSize',markerSize);
view(2); axis tight; axis equal; grid on; box on;
set(gca,'fontSize',fontSize);
drawnow;

%% Create guide curve
% The sweepLoft (see |HELP_sweepLoft|) is created here. First and angle
% based parameterization is created. Next this curve is evenly sample
% across the curve length (see |HELP_evenlySampleCurve|).
t=linspace(0,2*pi,numStepsCircumference*overSampleFactorCurve); %Angles
t=t(1:end-1); %Remove last point so it is not closed for resampling
x=stentRadius.*sin(t); %x-coordinates
y=stentRadius.*cos(t); %y-coordinates
z=waveAmplitude.*sin(numPeriodsWave*t); %z-coordinates
V_guide_curve=[x(:) y(:) z(:)]; %Collected curve nodes
[V_guide_curve] = evenlySampleCurve(V_guide_curve,numStepsCircumference-1,'pchip',1); %Resample curve evenly
V_guide_curve(end+1,:)=V_guide_curve(1,:); %Append start to end so it is a closed loop

%%
% Visualize guide curve
cFigure; hold on;
title('Stent guide curve','fontSize',fontSize);
plotV(V_guide_curve,'k.-','lineWidth',lineWidth,'MarkerSize',markerSize);
axisGeom;
drawnow;

%% Position stent section at the start and end of the guide curve
% Next the section is translated and rotated so it is placed at the start
% of the guide curve such that the curve normal points allong the curve.

% Create rotation matrix
n3=vecnormalize(V_guide_curve(2,:)-V_guide_curve(1,:)); %Out of section normal direction z ish direction
[~,indMin]=min(dot(n3(ones(1,2),:),[1 0 0; 0 1 0],2)); %Get index most appropriate initial other axis
switch indMin
    case 1
        n1=[1 0 0]; %Initialized x direction
        n2=vecnormalize(cross(n3,n1)); %y ish direction
        n1=vecnormalize(cross(n2,n3)); %Proper x ish direction
        R=[n1; n2; n3]; %Rotation matrix
    case 2
        n2=[0 1 0]; %Initialized y direction
        n1=vecnormalize(cross(n2,n3)); %x ish direction
        n2=vecnormalize(cross(n3,n1)); %Proper y ish direction
        R=[n1; n2; n3]; %Rotation matrix
end

p1=V_guide_curve(1,:); %The start node
V_section=V_section*R; %Rotate the section
V_section=V_section+p1(ones(size(V_section,1),1),:); % Translate coordinate to start

%%
% Visualize guide curve
cFigure; hold on;
title('Stent section positioned on guide curve','fontSize',fontSize);
plotV(V_guide_curve,'k-','lineWidth',1);
plotV(V_section,'k.-','lineWidth',lineWidth,'MarkerSize',markerSize);
quiverVec(p1,n1,1,'r');
quiverVec(p1,n2,1,'g');
quiverVec(p1,n3,1,'b');
axisGeom;
drawnow;

%% Sweeping section allong curve
% Normally |sweepLoft| produces patch data as an output (e.g. faces and
% vertices). However these outputs are supressed here and the coordinate
% mesh output is instead used to create a hexahedral mesh. See also |HELP_sweepLoft|

numTwist=0; %Number of additional twists of loft feature around guide curve
numStepsSweep=numStepsCircumference; %Number of steps for loft feature from sketch 1 to sketch 2
plotOn=0; %Turn on plotting to view lofting behaviour
[~,~,~,S]=sweepLoft(V_section,V_section,n3,n3,V_guide_curve,numStepsSweep,numTwist,plotOn);

%% Construct hexahedral mesh

X=S.X'; Y=S.Y'; Z=S.Z'; %Coordinate matrices
V=[X(:) Y(:) Z(:)]; %Create node list

F=reshape((1:1:size(V,1)),4,size(V,1)/4)'; %All top and bottom faces
E=[F(2:end,:) F(1:end-1,:)]; %The hexahedral elements
[E,V]=mergeVertices(E,V); %Merge nodes (start and end are not shared yet)
C=hexVol(E,V); %Get hexahedral element volumes

[F,CF]=element2patch(E,C); %Create face data for plotting

%%
% Visualize hexahedral mesh

cFigure; hold on;
title('Stent hexahedral mesh','fontSize',fontSize);
plotV(V_guide_curve,'k-','lineWidth',3);
gpatch(F,V,CF,'k',1);
% patchNormPlot(F,V); %Check normal directions
colormap gjet; colorbar;
axisGeom;
camlight headlight;
drawnow;

%% Refine mesh
% The swept mesh can be refined through slitting. The splitting can be
% homogeneous or only in a particular direction (see HELP_subHex|)
% Split method explanation:
% 1: Overall splitting in all directions
% 2: Split allong curve direction
% 3: Split axially
% 4: Splint inward

splitMethod=3;
nRefine=1;
[E,V,Cs]=subHex(E,V,nRefine,splitMethod);

splitMethod=4;
nRefine=1;
[E,V,Css]=subHex(E,V,nRefine,splitMethod);
C=hexVol(E,V); %Get hexahedral element volumes

Cs=Cs(Css); %Colors for original element indices (and sweeping steps)

[F,CF]=element2patch(E,C); %Create face data

%%
% Visualize hexahedral mesh
cFigure; hold on;
title('Stent hexahedral mesh','fontSize',fontSize);
plotV(V_guide_curve,'k-','lineWidth',3);
gpatch(F,V,CF,'k',1);
% patchNormPlot(F,V);
colormap gjet; colorbar;
axisGeom;
camlight headlight;
drawnow;

%% Create additional segments
E_cell=repmat({E},1,numSegments);
V_cell=repmat({V},1,numSegments);
for q=1:1:numSegments
    V_cell{q}(:,3)=V_cell{q}(:,3)+((q-1)*segmentAxialOffset);
end
[ET,VT,CT]=joinElementSets(E_cell,V_cell);

[FT,CTF]=element2patch(ET,CT); %Create face data for plotting

%% Create face labelling to construct node sets

[indBoundary]=tesBoundary(FT,VT);
faceMarker=ones(size(ET,1),1)*(1:6); %The 6 face colors for the hexahedral faces
faceMarker=faceMarker(:); %Force as a column
FTb=FT(indBoundary,:); %Select the boundary faces (which will exclude tops (1) and bottoms (2))
faceBoundaryMarker=faceMarker(indBoundary,:)-2; %Get boundary colors and subtract 2 so they are 1-4

%%
% Visualize face labels

cFigure; hold on;
title('Stent hexahedral mesh','fontSize',fontSize);
gpatch(FTb,VT,faceBoundaryMarker,'none',1);
% patchNormPlot(F,V); %Check normal directions
colormap gjet; colorbar;
axisGeom;
camlight headlight;
drawnow;

%%
% Visualize composed segments and animate sweep

hf=cFigure;
gtitle('Stent hexahedral mesh',fontSize);
subplot(1,2,1); hold on;
plotV(V_guide_curve,'k-','lineWidth',3);
hp=gpatch(F,V,'gw','k',1);
plotV(V_guide_curve,'k-','lineWidth',3);
axisGeom;
camlight headlight;
axis off; axis manual;

subplot(1,2,2); hold on;
gpatch(FT,VT,'gw','none',1);
axisGeom;
camlight headlight;
lighting gouraud
axis off;
drawnow;

%%

[F,CFs]=element2patch(E,Cs); %Create face data for plotting

nSteps=40; %Number of animation steps
animStruct.Time=round(linspace(1,max(CFs(:)),nSteps));%Create the time vector

for q=1:1:nSteps
    F_now=F(ismember(CFs,1:animStruct.Time(q)),:);
    %Set entries in animation structure
    animStruct.Handles{q}=hp; %Handles of objects to animate
    animStruct.Props{q}={'Faces'}; %Properties of objects to animate
    animStruct.Set{q}={F_now}; %Property values for to set in order to animate
end

anim8(hf,animStruct);

%% Create node set for inner surface nodes

FTb_inner=FTb(faceBoundaryMarker==2,:); %Inner faces
indicesNodesInner=unique(FTb_inner(:)); %Inner nodes
indicesElementsInner=find(sum(ismember(ET,indicesNodesInner),2)==4); %Inner elements

%%
% Visualize face labels

[F_inner_elements]=element2patch(ET(indicesElementsInner,:)); %Create face data for plotting
C_side=ones(numel(indicesElementsInner),1)*(1:1:6);
C_side=C_side(:);

cFigure; hold on;
title('Stent hexahedral mesh','fontSize',fontSize);
gpatch(FTb,VT,'kw','none',0.2);
% gpatch(FTb_inner,VT,'gw','k',1);
gpatch(F_inner_elements,VT,C_side,'none',1);
% plotV(VT(indicesNodesInner,:),'g.','MarkerSize',25);
colormap(gjet(6)); icolorbar;
axisGeom;
camlight headlight;
drawnow;

%% Setup structure to define an Abaqus inp file

%%--> Heading
abaqus_spec.Heading.COMMENT{1}='Job name: ABAQUS inp file creation demo';
abaqus_spec.Heading.COMMENT{2}='Generated by: GIBBON';

%%--> Preprint
abaqus_spec.Preprint.ATTR.echo='NO';
abaqus_spec.Preprint.ATTR.model='NO';
abaqus_spec.Preprint.ATTR.history='NO';
abaqus_spec.Preprint.ATTR.contact='NO';

%--> Part

% Node
nodeIds=(1:1:size(VT,1))';
abaqus_spec.Part.COMMENT='This section defines the part geometry in terms of nodes and elements';
abaqus_spec.Part.ATTR.name='Stent';
abaqus_spec.Part.Node={nodeIds,VT};

% Element
elementIds=(1:1:size(ET,1))';
abaqus_spec.Part.Element{1}.ATTR.type='C3D8';%'C3D8R';
abaqus_spec.Part.Element{1}.VAL={elementIds,ET};

% Element sets
abaqus_spec.Part.Elset{1}.ATTR.elset='Set-1';
abaqus_spec.Part.Elset{1}.VAL=elementIds;

surfaceElementSetName='elementSetSideSurface';
abaqus_spec.Part.Elset{2}.ATTR.elset=surfaceElementSetName;
abaqus_spec.Part.Elset{2}.ATTR.internal=''; %Remains hidden uppon import
abaqus_spec.Part.Elset{2}.VAL=indicesElementsInner(:);

% Surfaces
sidePick=5;
abaqus_spec.Part.Surface{1}.ATTR.type='ELEMENT';
abaqus_spec.Part.Surface{1}.ATTR.name=[surfaceElementSetName,'_side',num2str(sidePick)];
abaqus_spec.Part.Surface{1}.VAL={surfaceElementSetName,['S',num2str(sidePick)]};

% Sections
abaqus_spec.Part.Solid_section.ATTR.elset='Set-1';
abaqus_spec.Part.Solid_section.ATTR.material='Elastic';

%--> Assembly
abaqus_spec.Assembly.ATTR.name='Assembly-1';
abaqus_spec.Assembly.Instance.ATTR.name='Stent-assembly';
abaqus_spec.Assembly.Instance.ATTR.part='Stent';

abaqus_spec.Assembly.Nset{1}.ATTR.nset='Set-1';
abaqus_spec.Assembly.Nset{1}.ATTR.instance=abaqus_spec.Assembly.Instance.ATTR.name;
abaqus_spec.Assembly.Nset{1}.VAL=indicesNodesInner(:);

%%

%%--> Material
abaqus_spec.Material.ATTR.name='Elastic';
abaqus_spec.Material.Elastic=[0.5 0.4];

%%--> Step
abaqus_spec.Step.ATTR.name='Step-1';
abaqus_spec.Step.ATTR.nlgeom='YES';
abaqus_spec.Step.Static=[0.1 1 1e-5 0.1];

% % Boundary
% abaqus_spec.Step.Boundary{1}.VAL={'Set-1',[1,1]};
% abaqus_spec.Step.Boundary{2}.VAL={'Set-2',[2,2]};
% abaqus_spec.Step.Boundary{3}.VAL={'Set-3',[3,3]};
% abaqus_spec.Step.Boundary{4}.VAL={'Set-4',[3,3],-0.1};

%Output
abaqus_spec.Step.Restart.ATTR.write='';
abaqus_spec.Step.Restart.ATTR.frequency=0;

abaqus_spec.Step.Output{1}.ATTR.field='';
abaqus_spec.Step.Output{1}.ATTR.variable='PRESELECT';
abaqus_spec.Step.Output{2}.ATTR.history='';
abaqus_spec.Step.Output{2}.ATTR.variable='PRESELECT';
% abaqus_spec.Step.Node_print.ATTR.nset='all';
% abaqus_spec.Step.Node_print.ATTR.frequency = 1;
% abaqus_spec.Step.Node_print.VAL='COORD';
% abaqus_spec.Step.El_print.VAL='S';

%% Creating the INP file
% You can use |abaqusStruct2inp| to write the structure data to a file.

%Create file name for INP file
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');
fileName=fullfile(savePath,'tempModel.inp');
[~,fileNamePart,~]=fileparts(fileName);

% Export INP file
abaqusStruct2inp(abaqus_spec,fileName);

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function varargout = graph2Contours(nodes, edges) %#ok<INUSL>
%GRAPH2CONTOURS Convert a graph to a set of contour curves.
%
%   CONTOURS = GRAPH2CONTOURS(NODES, EDGES)
%   NODES, EDGES is a graph representation (type "help graphs" for details)
%   The algorithm assume every node has degree 2, and the set of edges
%   forms only closed loops. The result is a list of indices arrays, each
%   array containing consecutive point indices of a contour.
%
%   To transform contours into drawable curves, please use :
%   CURVES{i} = NODES(CONTOURS{i}, :);
%
%
%   NOTE : contours are not oriented. To manage contour orientation, edges
%   also need to be oriented. So we must precise generation of edges.

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2004-08-05
% Copyright 2004-2022 INRA - TPV URPOI - BIA IMASTE

curves = {};
c = 0;

while size(edges,1)>0
	% find first point of the curve
	n0 = edges(1,1);
    curve = n0;

    % second point of the curve
	n = edges(1,2);
	e = 1;

	while true
        % add current point to the curve
		curve = [curve n];         %#ok<AGROW>

        % remove current edge from the list
        edges = edges((1:size(edges,1))~=e,:);

		% find index of edge containing reference to current node
		e = find(edges(:,1)==n | edges(:,2)==n);
        e = e(1);

		% get index of next current node
        % (this is the other node of the current edge)
		if edges(e,1)==n
            n = edges(e,2);
		else
            n = edges(e,1);
		end

        % if node is same as start node, loop is closed, and we stop
        % node iteration.
        if n==n0
            break;
        end
	end

    % remove the last edge of the curve from edge list.
    edges = edges((1:size(edges,1))~=e,:);

    % add the current curve to the list, and start a new curve
    c = c+1;
    curves{c} = curve; %#ok<AGROW>
end

if nargout == 1
    varargout = {curves};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [ r, seed ] = r8mat_uniform_01 ( m, n, seed )

%*****************************************************************************80
%
%% R8MAT_UNIFORM_01 returns a unit pseudorandom R8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in the array.
%
%    Input, integer SEED, the integer "seed" used to generate
%    the output random number.
%
%    Output, real R(M,N), an array of random values between 0 and 1.
%
%    Output, integer SEED, the updated seed.  This would
%    normally be used as the input seed on the next call.
%
  i4_huge = 2147483647;
  r = zeros ( m, n );

  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_UNIFORM_01 - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R8MAT_UNIFORM_01 - Fatal error!' );
  end

  for j = 1 : n
    for i = 1 : m

      seed = floor ( seed );

      seed = mod ( seed, i4_huge );

      if ( seed < 0 )
        seed = seed + i4_huge;
      end

      k = floor ( seed / 127773 );

      seed = 16807 * ( seed - k * 127773 ) - k * 2836;

      if ( seed < 0 )
        seed = seed + i4_huge;
      end

      r(i,j) = seed * 4.656612875E-10;

    end
  end

  return
end


text:
function m = nanmean(x,dim)
%NANMEAN Mean value, ignoring NaNs.
%   M = NANMEAN(X) returns the sample mean of X, treating NaNs as missing
%   values.  For vector input, M is the mean value of the non-NaN elements
%   in X.  For matrix input, M is a row vector containing the mean value of
%   non-NaN elements in each column.  For N-D arrays, NANMEAN operates
%   along the first non-singleton dimension.
%
%   NANMEAN(X,DIM) takes the mean along dimension DIM of X.
%
%   See also MEAN, NANMEDIAN, NANSTD, NANVAR, NANMIN, NANMAX, NANSUM.

%   Copyright 1993-2004 The MathWorks, Inc.
%   $Revision: 1.1.8.1 $  $Date: 2010/03/16 00:15:50 $

% Find NaNs and set them to zero
nans = isnan(x);
x(nans) = 0;

if nargin == 1 % let sum deal with figuring out which dimension to use
    % Count up non-NaNs.
    n = sum(~nans);
    n(n==0) = NaN; % prevent divideByZero warnings
    % Sum up non-NaNs, and divide by the number of non-NaNs.
    m = sum(x) ./ n;
else
    % Count up non-NaNs.
    n = sum(~nans,dim);
    n(n==0) = NaN; % prevent divideByZero warnings
    % Sum up non-NaNs, and divide by the number of non-NaNs.
    m = sum(x,dim) ./ n;
end

author: marianux

repo: ecg-kit

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


text:
function X=stfth(x,wl)
disp('Qverlooping of window is 50%');
disp('1 Rectangular Window, 2 Hamming window,  3 Hanning window');
window=input('Enter your choice');
L=length(x);
if L<wl
    z=wl-L;
    x=[x,zeros(1,z)];
end
switch window  % Window functions are (all =0 outside the interval 0<=n<=wl)
    case 1
        win=ones(1,wl); % Rectangular Window
    case 2
        win=hamming(wl)'; % Hanning window w[n]=0.5-0.5*cos(2*pi*n/M), for 0<=n<=M
    case 3
        win=hann(wl)';  % Hamming window w[n]=0.54-.46*cos(2*pi*n/M) for 0<=n<=M
    otherwise
        win=ones(1,window_length);
        disp('Not a right option, By default rectangular window is taken.');
end

L=length(x);
hop=ceil(wl/2);  % Hoop size of window
if hop<1
    hop=wl;
end
i=1;str=1; len=wl; X=[];
while (len<=L || i<2)
    if i==1
    if  len>L   % If window size excceds the L of signal for 1st time
        z=len-L;
        x=[x,zeros(1,z)]; % padding zeros
        i=1+1;
    end
    x1=x(str:len);
    X=[X;fft(x1.*win)];  % Matrix mul
    str=str+hop; len=str+wl-1; % to make window overlapping
        end
    end
figure,subplot(2,1,1)
imagesc(abs(X)); title('Result by imagesc function')

subplot(2,1,2)
surf(abs(X)); title('Result by surf function')


text:
function[x]=cellreal(x)
%CELLREAL  Real part of each element in a cell array.
%
%   XR=CELLREAL(X) where X is a cell array of N arrays, is equivalent to
%
%      XR{1}=REAL(X{1}), XR{2}=REAL(X{2}), ..., XR{N}=REAL(X{N})
%
%   thus returning the real part of each element in the cell array.
%
%   See also JCELL.
%
%   Usage: xr=cellreal(x);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2009--2019 J.M. Lilly --- type 'help jlab_license' for details

if ~iscell(x)
    error('X must be a cell array.')
end
for i=1:length(x)
    x{i}=real(x{i});
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function bbci= bbci_save(bbci, data)
%BBCI_SAVE - Save BBCI classifier
%
%Synopsis:
%  BBCI= bbci_save(BBCI, <DATA>)
%
%To get a description on the structures 'bbci', type
%help bbci_calibrate_structures

% 09-2011 Benjamin Blankertz

bbci= bbci_save_setDefaults(bbci);

if nargin<2,
  data= [];
end
BS= bbci.calibrate.save;

file= fullfile(BS.folder, BS.file);
num= 1;
while ~BS.overwrite && exist([file '.mat'], 'file'),
  num= num + 1;
  file= fullfile(BS.folder, [BS.file, sprintf('%02d',num)]);
end

save(file, '-STRUCT', 'bbci');
msg= sprintf('Classifier saved to <%s.mat>', file);

if ~isempty(data),
  if ~BS.raw_data,
    data= rmfield(data, {'cnt','mrk','mnt'});
  end
  if isequal(BS.data, 'separately'),
    save([file '_data'], '-STRUCT', 'data');
    msg= [msg ' and data separately to *_data.mat'];
  else
    save(file, 'data', '-APPEND');
  end
end

if BS.figures,
  if ~isfield(data, 'figure_handles'),
    data.figure_handles= sort(findobj('Type','figure'));
  end
  fig_folder= strcat(file, '_figures');
  if ~exist(fig_folder, 'dir'),
    mkdir(fig_folder);
  end
  for ff= data.figure_handles(:)',
    figure(ff);
    fig_name= strrep(get(ff,'Name'), ' ', '_');
    if ispc,
      fig_name= strrep(fig_name, '.', '_');
    end
    filename= fullfile(fig_folder, sprintf('Fig-%02d_%s', ff, fig_name));
    util_printFigure(filename, BS.figures_spec{:});
  end
  msg= [msg ' and figures in subfolder _figures'];
end
fprintf('%s.\n', msg);

author: bbci


text:
function p = multipdf(x,theta)
%MULTIPDF Multinomial probability density function.
%   p = multipdf(x,theta) returns the probabilities of
%   vector x, under the multinomial distribution
%   with parameter vector theta.
%
%   Author: David Ross

%--------------------------------------------------------
% Check the arguments.
%--------------------------------------------------------
error(nargchk(2,2,nargin));

% make sure theta is a vector
if ndims(theta) > 2 | all(size(theta) > 1)
    error('theta must be a vector');
end

% make sure x is of the appropriate size
if ndims(x) > 2 | any(size(x) ~= size(theta))
    error('columns of X must have same length as theta');
end

%--------------------------------------------------------
% Main...
%--------------------------------------------------------
p = prod(theta .^ x);
p = p .* factorial(sum(x)) ./ prod(factorial_v(x));

%--------------------------------------------------------
% Function factorial_v(x): computes the factorial function
% on each element of x
%--------------------------------------------------------
function r = factorial_v(x)

if size(x,2) == 1
    x = x';
end

r = [];
for y = x
    r = [r factorial(y)];
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text: save  data1  a  b  c

author: Eurus-Holmes


text:
% edge_response_1d
% show how edge response is nonlinear for edge-preserving regularization

if ~isvar('ytrue'), printm 'setup'
	nx = 2^5;
	ny = nx;
	ny = 20*nx;
	steps = [1 2 4 8];
	nstep = numel(steps);
	xtrue = zeros(nx, ny, 'single');
	ny_per = ny / nstep;
	for ii=1:nstep
		iy{ii} = [1:ny_per] + (ii-1)*ny_per;
		xtrue(nx/2+1:end, iy{ii}) = steps(ii); % for-step
%		xtrue(round(ii/(nstep+1)*nx), iy{ii}) = steps(ii); % impulse
	end

	psf = [1 2 1]';
	mask = true(nx,ny);
	A = Gblur(mask, 'psf', psf);
	ytrue = A * xtrue;

	if 1
		im('notick', xtrue, ' ')
%		ir_savefig edge_response_1d_xtrue
	prompt
	end

	if 1
		im('notick', ytrue, ' ')
%		ir_savefig edge_response_1d_ytrue
	prompt
	end
end

if ~isvar('R'), printm 'R'
	l2b = -0;
	bet = 2^l2b; % essentially product of beta * sigma^2
%	f.pot = {'quad'};
	f.pot = {'hyper3', 1.0};

	% todo: use periodic boundary conditions?
% todo: compare noise with step
	R = Reg1(mask, 'beta', bet, 'type_penal', 'mat', ...
		'offsets', 1, 'pot_arg', f.pot); % 1d regularization
	if 1
		qpwls_psf(A, R, 1, mask);
	prompt
	end
end

rng(0)
%sig = 0.0;
sig = 0.1;
yi = ytrue + sig * randn(size(ytrue));

if ~isvar('xh'), printm 'xh'
	xh = pwls_pcg1(xtrue(mask), A, 1, yi(:), R, 'niter', 30);
	xh = embed(xh, mask);

	if 1
		im('notick', xh, ' ')
%		cbar, colormap(jet)
%		ir_savefig edge_response_1d_xh
	prompt
	end
end

if 1 % plots
	tmp = zeros(nx, nstep);
	stds = zeros(nx, nstep);
	leg = cell(nstep,1);
	for ii=1:nstep
		tmp(:,ii) = mean(xh(:, iy{ii}),2) / steps(ii); % normalize
		stds(:,ii) = std(xh(:, iy{ii}),0,2);
		leg{ii} = num2str(steps(ii));
	end
	if 0 % for-step
		plot(-nx/2:nx/2-1, tmp, '.-')
		axis([0*nx/2+[-1 1]*5 -0.1 1.1])
	else
		plot(tmp, '.-')
		axis([0 nx -0.1 1.1])
		plot(stds, '.-')
	end
	ir_legend(leg)
%	ytick([0 1])
	xlabel 'horizontal location'
	ylabel 'normalized profile'
%	ir_savefig eps_c edge_response_1d_p1
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%{
    Tests the Dantzig Selector

    min_x ||x||_1
s.t.
    || D*A'*(A*x - b) || <= delta

The solvers solve a regularized version, using
    ||x||_1 + mu/2*||x-x_0||_2^2

see also test_sDantzig.m

This demo shows three formulations of the Dantzig selector
Instead of calling a pre-built solver, we show how to call
tfocs_SCD directly.

%}

% Before running this, please add the TFOCS base directory to your path

% Try to load the problem from disk
mu = 0;
fileName = fullfile('reference_solutions','dantzig_problem1_smoothed_noisy');
randn('state',34324);
rand('state',34324);
N = 1024;
M = round(N/2);
K = round(M/5);
A = randn(M,N);
if exist([fileName,'.mat'],'file')
    load(fileName);
    fprintf('Loaded problem from %s\n', fileName );
else
    disp('Please run test_sDantzig.m to setup the file');
end

[M,N]           = size(A);
K               = nnz(x_original);
norm_x_ref      = norm(x_ref);
norm_x_orig     = norm(x_original);
er_ref          = @(x) norm(x-x_ref)/norm_x_ref;
er_signal       = @(x) norm(x-x_original)/norm_x_orig;
resid           = @(x) norm(A*x-b)/norm(b);  % change if b is noisy

fprintf('\tA is %d x %d, original signal has %d nonzeros\n', M, N, K );
fprintf('\tl1-norm solution and original signal differ by %.2e (mu = %.2e)\n', ...
    norm(x_ref - x_original)/norm(x_original),mu );

%% Call the TFOCS solver
er              = er_ref;  % error with reference solution (from IPM)
opts = [];
opts.restart    = 500;
opts.errFcn     = { @(f,dual,primal) er(primal), ...
                    @(f,dual,primal) obj_ref - f  };
opts.maxIts     = 1000;
opts.printEvery = 100;
z0  = [];   % we don't have a good guess for the dual

%% Method 1: use the epigraph of the l_infinity norm as the cone
% Note: this is equivalent to calling:
% [ x, out, optsOut ] = solver_sDantzig( {A,D}, b, delta, mu, x0, z0, opts );

DAtb = D.*(A'*b);
DD = @(x) D.*(x);
objectiveF = prox_l1;
affineF     = {linop_matrix(diag(D)*(A'*A)), -DAtb };
dualproxF  = prox_l1( delta );
[ x, out, opts ] = tfocs_SCD( objectiveF, affineF, dualproxF, mu, x0, z0, opts );
x1 = x;
out1 = out;

fprintf('Solution has %d nonzeros.  Error vs. IPM solution is %.2e\n',...
    nnz(x), er(x) );

% Check that we are within allowable bounds
if out.err(end,1) < 1e-3
    disp('Everything is working');
else
    error('Failed the test');
end
%% Method 2: use the LP formulation
% Instead of the constraint ||Ax-b||_infty <= delta
% (where "A" is really DA'A),
% think of it as
% -(Ax-b) + delta >= 0
%  (Ax-b) + delta >= 0

objectiveF = prox_l1;
affineF     = {linop_matrix(-diag(D)*(A'*A)),  DAtb + delta;...
               linop_matrix( diag(D)*(A'*A)), -DAtb + delta; };

dualproxF  = { proj_Rplus; proj_Rplus };
[ x, out, opts ] = tfocs_SCD( objectiveF, affineF, dualproxF, mu, x0, z0, opts );
x2 = x;
out2 = out;

% Check that we are within allowable bounds
if out.err(end,1) < 1e-3
    disp('Everything is working');
else
    error('Failed the test');
end
%% Method 3: put objective into constraint
% This is the trick we do with solver_sDantzig_W, to deal with ||Wx||_1
% Instead of minimizing ||x||_1, we minimize t,
%   with the constraint that ||x||_1 <= t
% This version has some scaling considerations -- if we
% are careful about scaling, the dual problem will be solved much faster.

% Note: we still have to deal with the original constraints.  We
% can use either method 1 or method 2 above.  Here, we'll use
% method 1.

objectiveF  = [];    % tfocs_SCD recognizes this special objective
normA2      = norm( diag(D)*A'*A )^2;
scale       = 1/sqrt(normA2);
affineF     = {linop_matrix(diag(D)*(A'*A)), -DAtb; linop_scale(1/scale), 0 };
dualproxF   = { prox_l1(delta); proj_linf(scale) };
[ x, out, opts ] = tfocs_SCD( objectiveF, affineF, dualproxF, mu, x0, z0, opts );
x3 = x;
out3 = out;

% Check that we are within allowable bounds
if out.err(end,1) < 1e-2
    disp('Everything is working');
else
    error('Failed the test');
end
%% plot
figure;
semilogy( out1.err(:,1) );
hold all
semilogy( out2.err(:,1) );
semilogy( out3.err(:,1) );
legend('Method 1 (epigraph cone)','Method 2 (LP)','Method 3 (W=I)' );

%% close all plots
close all

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

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [A uTrue uEst] = buildAlignmentMatrix( zTrue, zEst )
% Given true and estimated label sequence, construct align matrix A
%   where A( ke, kt ) = fraction of labels ke aligned
%                          from ESTIMATED label ke to TRUE label kt
% Note that each ROW of A will sum to 1

uTrue = unique( zTrue );
uEst  = unique( zEst  );

for ue = 1:length( uEst )
    ts = zEst == uEst(ue);
    A(ue,:) = histc( zTrue(ts),  uTrue );
end
A = bsxfun( @rdivide, A, sum(A,2)  );

author: michaelchughes

repo: NPBayesHMM


text:
function planC = contourSUV(structNum,percent,planC)
%function contourSUV(structNum,percent,planC)
%
%This function creates structure at percent% SUV level for structNum
%
%APA,12/15/2006
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global stateS
if ~exist('planC','var')
    global planC
end
indexS = planC{end};

scanNum                             = getStructureAssociatedScan(structNum,planC);
[rasterSegments, planC, isError]    = getRasterSegments(structNum,planC);
if isempty(rasterSegments)
    warning('Could not create conotour.')
    return
end
[mask3M, uniqueSlices]              = rasterToMask(rasterSegments, scanNum, planC);
scanArray3M                         = double(getScanArray(planC{indexS.scan}(scanNum)));
%%%% reverse data
%maxscan=max(scanArray3M(:));
%scanArray3M =maxscan-scanArray3M;
SUVvals3M                           = mask3M.*scanArray3M(:,:,uniqueSlices);
maxSUVinStruct                      = max(SUVvals3M(:));
cutoff                              = percent/100*maxSUVinStruct;
[xVals, yVals, zVals]               = getScanXYZVals(planC{indexS.scan}(scanNum));
newStructNum                        = length(planC{indexS.structures}) + 1;

newStructS = newCERRStructure(scanNum, planC);
for slcNum = 1:length(uniqueSlices)
    C = contourc(xVals, yVals, SUVvals3M(:,:,slcNum),[cutoff cutoff]);
    indC = getSegIndices(C);
    if ~isempty(indC)
        for seg = 1:length(indC)
            points = [C(:,indC{seg})' zVals(uniqueSlices(slcNum))*ones(length(C(1,indC{seg})),1)];
            newStructS.contour(uniqueSlices(slcNum)).segments(seg).points = points;
        end
    else
        newStructS.contour(uniqueSlices(slcNum)).segments.points = [];
    end
end

for l = max(uniqueSlices)+1 : length(planC{indexS.scan}(scanNum).scanInfo)
    newStructS.contour(l).segments.points = [];
end

newStructS.structureName    = [planC{indexS.structures}(structNum).structureName 'SUVMax',num2str(percent)];

planC{indexS.structures} = dissimilarInsert(planC{indexS.structures}, newStructS, newStructNum);
planC = getRasterSegs(planC, newStructNum);
planC = updateStructureMatrices(planC, newStructNum, uniqueSlices);

if ~isempty(stateS) && isfield(stateS,'handle') && isfield(stateS.handle,'CERRSliceViewer') && ishandle(stateS.handle.CERRSliceViewer)
    stateS.structsChanged = 1;
    % Refresh View
    CERRRefresh
end

return;

function indC = getSegIndices(C)
% function getSegIndices(C)
%
%This function returns the indices for each segment of input contour C.
%C is output from in-built "contourc" function
%
%APA, 12/15/2006

start = 1;
counter = 1;
indC = [];
while start < length(C(2,:))
    numPts = C(2,start);
    indC{counter} = [(start+1):(start+numPts) start+1];
    start = start + numPts + 1;
    counter = counter + 1;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function varargout = HoughObject(varargin)
% HOUGHOBJECT M-file for HoughObject.fig
%      HOUGHOBJECT, by itself, creates a new HOUGHOBJECT or raises the existing
%      singleton*.
%
%      H = HOUGHOBJECT returns the handle to a new HOUGHOBJECT or the handle to
%      the existing singleton*.
%
%      HOUGHOBJECT('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in HOUGHOBJECT.M with the given input arguments.
%
%      HOUGHOBJECT('Property','Value',...) creates a new HOUGHOBJECT or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before HoughObject_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to HoughObject_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help HoughObject

% Last Modified by GUIDE v2.5 07-Dec-2005 20:15:07

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @HoughObject_OpeningFcn, ...
                   'gui_OutputFcn',  @HoughObject_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before HoughObject is made visible.
function HoughObject_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to HoughObject (see VARARGIN)

% Choose default command line output for HoughObject
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes HoughObject wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = HoughObject_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pbGenerate.
function pbGenerate_Callback(hObject, eventdata, handles)
% hObject    handle to pbGenerate (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
load ObjectTemplate
S1 = zeros(180,180);
Str = {'Sq','Sc','St'};

a = 1; b = 3;
x = round(a + (b-a) * rand(1));
S2 = eval(Str{x});

a = 1; b = 90;
x = round(a + (b-a) * rand(1));
S3 = imrotate(S2,x);

S = S1;
sz = size(S3);

a = 1; b = min(sz);
x = round(a + (b-a) * rand(1));

S(x:sz(1)+x-1,x:sz(2)+x-1)=S3;
S = im2bw(S);
axes(handles.axes1);

imshow(S);
title('Original Image');
handles.S = S;
guidata(hObject, handles);

% --- Executes on button press in pnHough.
function pnHough_Callback(hObject, eventdata, handles)
% hObject    handle to pnHough (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
S = handles.S;
[H, theta,rho]=hough(S);
axes(handles.axes2);
imshow(H,[],'xdata',theta,'ydata',rho);
xlabel('\theta'),ylabel('\rho')
axis on, axis normal;
title('Hough Matrix');

clear data;
for cnt = 1:max(max(H))
    data(cnt) = sum(sum(H == cnt));
end
axes(handles.axes3);
%data(data==0)=NaN;
plot(data,'--x');
xlabel('Hough Matrix Intensity'),ylabel('Counts')
handles.data = data;
guidata(hObject, handles);

% --- Executes on button press in pbDetect.
function pbDetect_Callback(hObject, eventdata, handles)
% hObject    handle to pbDetect (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
data = handles.data;
[maxval,maxind] = max(data);
medval = median(data);

[p]=polyfit(1:maxind-5,data(1:maxind-5),2);

if maxval<3*medval
    set(handles.txtResult,'string','Triangle');
elseif  p(3)>100
    set(handles.txtResult,'string','Square');
else
    set(handles.txtResult,'string','Round');
end

% --- Executes on button press in pbInfo.
function pbInfo_Callback(hObject, eventdata, handles)
% hObject    handle to pbInfo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
web('http://basic-eng.blogspot.com', '-browser');


text:
function ncorr_util_colormap(handle_fig)
% This function adjusts the colormap correctly. This was taken from some
% matlab code.
%
% Inputs -----------------------------------------------------------------%
%   handle_fig - handle; handle of the figure;
%
% Outputs ----------------------------------------------------------------%
%   none;

    m = size(get(handle_fig,'colormap'),1);
    n = ceil(m/4);
    u = [(1:1:n)/n ones(1,n-1) (n:-1:1)/n]';
    g = ceil(n/2) - (mod(m,4)==1) + (1:length(u))';
    r = g + n;
    b = g - n;
    g(g>m) = [];
    r(r>m) = [];
    b(b<1) = [];
    J = zeros(m,3);
    J(r,1) = u(1:length(r));
    J(g,2) = u(1:length(g));
    J(b,3) = u(end-length(b)+1:end);

    if verLessThan('matlab','9.1')
        set(handle_fig,'Colormap',J);
    else
        % Get all axes and set their colormap. This is a fix for 2016b+
        % (9.1+) which now sets colormaps on a per-axis basis for imshow()
        % instead of for the entire figure.
        handle_axes = findobj(handle_fig,'type','axes');
        for i = 1:length(handle_axes)
            colormap(handle_axes(i),J);
        end
    end
end

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
% The OxyGEN proyect by Profoty.xyz
% Together against the Covid-19
% V4.0 rev.17/03/2020
%
% Website/blog: oxygen.protofy.xyz
% Contact: oxygen@protofy.xyz

%GEOMETRICAL PARAMETERS

%Length from the bearing to the hinge
l0 = 27.3;

%Length of the vertical wall that supports the hinge
h0 = 6.5;

%Lenght of the vertical support of the bearing
hb = 5.0;

%Bearing radius
br = 1.1;

%Array of different max and min positions of the bearing
h1a = [6.5, 8.5, 10.5, 12.5];
h2a = [20.5, 20.5, 20.5, 20.5];

%Number of cycles per cam turn (1 to 3).
nc = 1;

%Minimum radius of the camshaft
rmin = 5;

%BREATHING CYCLE PARAMETERS

%Duration of the inhale cycle / duration of the whole cycle
lambda1 = 0.500;
lambda2 = 0.060;

%Soft transition between inhale and exhale cycles
dpsi21 = 0.12;
dpsi12 = 0.15;

%Adjust parameters of the inhale curve
ga1 = pi;
gb1 = 1.2;
ff1 = 20;
sp1 = 1;

%Adjust parametes of the exhale curve
ga2 = 9;
gb2 = .5;
ff2 = 20;

%GENERATION OF THE CURVES AND THE CAMSHAFT

%Generation of the geometry

l = sqrt(l0^2+hb^2);

for i = 1:numel(h1a);

    ymin(i) = h1a(i);
    ymax(i) = h2a(i);

    alphamin(i) = acos((h0 - h1a(i)) / l);
    alphamax(i) = acos((h0 - h2a(i)) / l);

    xmin(i) = l*sin(alphamin(i));
    xmax(i) = l*sin(alphamax(i));

    d(i) = sqrt((xmin(i)-xmax(i))^2 + (ymin(i)-ymax(i))^2);

    alphatan(i) = (alphamin(i) + alphamax(i)) / 2;

    xtan(i) = l*sin(alphatan(i));
    ytan(i) = h0 - l*cos(alphatan(i));

    xsup(i) = xtan(i) + (d(i)/2)*cos(alphatan(i));
    xinf(i) = xtan(i) - (d(i)/2)*cos(alphatan(i));

    ysup(i) = ytan(i) + (d(i)/2)*sin(alphatan(i));
    yinf(i) = ytan(i) - (d(i)/2)*sin(alphatan(i));

    xcam(i) = xtan(i) + (d(i)/2 + rmin + br)*cos(alphatan(i));
    ycam(i) = ytan(i) + (d(i)/2 + rmin + br)*sin(alphatan(i));

end;

%Time increment
dt = 0.01;

%time coordinate during the whole cycle
theta = 0:dt:2*pi;

%Generation of the soft transition between inhale and exhale curves
psi1 = [];
psi2 = [];

for i = 1:numel(theta);

    if theta(i) < (lambda2-dpsi21/2)*2*pi;
        psi1(i) = 0;

    elseif theta(i) < (lambda2+dpsi21/2)*2*pi;
        psi1(i) = 0.5 + 0.5*cos((theta(i)-(lambda2+dpsi21/2)*2*pi)/(2*dpsi21));

    else theta(i) < (lambda1-dpsi12/2)*2*pi;
        psi1(i) = 1;

    end;

    psi2(i) = 1 - psi1(i);

end;

%Generation of the complete breathing cycle rho(theta)

rho = [];
rho1 = [];
rho2 = [];
rho2next = [];

rhomin = 1000;
rhomax = 0;

for i = 1:numel(theta);

    %Inhale curve
    rho1(i) = ff1*normpdf(sp1*theta(i), ga1, gb1);
    rho1next(i) = ff1*normpdf(sp1*theta(i)+sp1*2*pi, ga1, gb1);

    %Exhale curve
    rho2(i) = ff2*gampdf(theta(i), ga2, gb2);
    rho2next(i) = ff2*gampdf(theta(i)+2*pi, ga2, gb2);

    rho(i) = rho1(i);

    %Capturing min and max in order to generate the normalized curve
    if rho(i) > rhomax
        rhomax = rho(i);
    end;

    if rho(i) < rhomin
        rhomin = rho(i);
    end;

end;

%Generation of a normalised curve and camshaft

rhonorm = [];
rhocam = [];

for i = 1:numel(theta);

    rhonorm(i) = (rho(i)-rhomin)/(rhomax-rhomin);

    for n = 1:numel(d)

        rhocam(n,i) = rmin + rhonorm(i)*d(n);

    end;

end;

%Generation of the first derivate of the camshaft geometry to analize and
%validate the design

drho = [];
drhonorm = [];
drhocam = [];

a = rmin;
b = 1/dt;

for i = 1 : numel(rho)-1;

    drho(i) = b*(rho(i+1)-rho(i));
    drhonorm(i) = b*(rhonorm(i+1)-rhonorm(i));
    drhocam(i) = b*(rhocam(i+1)-rhocam(i)) + a;

end;

drho(numel(rho)) = b*(rho(1)-rho(numel(rho)));
drhonorm(numel(rhonorm)) = b*(rhonorm(1)-rhonorm(numel(rhonorm)));
drhocam(numel(rhocam)) = b*(rhocam(1)-rhocam(numel(rhocam))) + a;

%X and Y coordinates during two cycles (for 2-cycle camshaft plot)

theta2 = [theta/2, (2*pi+theta)/2];

rhocam2 = [rhocam, rhocam];
drhocam2 = [drho, drho];

%X and Y coordinates during three cycles (for 3-cycle camshaft plot)

theta3 = [theta/3, (2*pi+theta)/3, (4*pi+theta)/3];

rhocam3 = [rhocam, rhocam, rhocam];
drhocam3 = [drho, drho, drho];

%GENERATION OF THE PLOTS

%Trying to print it out in real size (FAIL)
%set(gcf,'PaperUnits','centimeters');
%set(gcf,'PaperSize',[42 29.7]);

fpos = figure('Name', 'Dimensions', 'Units', 'centimeters', 'NumberTitle', 'off');
set(gcf,'PaperUnits','centimeters', 'PaperSize', [42/2 27.9/2]);
set(gcf, 'units', 'centimeters', 'position', [0, 0, 28, 20]);
fpos = gcf;
hold on

xlim([-5 35]);
ylim([0 30]);

for i = 1:numel(xmin)

    plot([0, xmin(i)], [h0, ymin(i)], '-x')
    plot([0, xmax(i)], [h0, ymax(i)], '-x')

    plot([0, xtan(i)], [h0, ytan(i)], '--xb')

    plot([xsup(i), xinf(i)], [ysup(i), yinf(i)], '--xr')
    scatter(xcam(1), ycam(1))

end;

hold off

%Plot of the breathing cycle interpolation by curves

fcycle = figure('Name', 'Breath Cycle curves', 'Units', 'centimeters', 'NumberTitle', 'off');

hold on
plot(theta, psi1, '-k')
plot(theta, psi2, '-k')
plot(theta, rho1, '-g')
plot(theta, rho1next, '-g')
plot(theta, rho2, '-g')
plot(theta, rho2next, '-g')
%plot(theta, drho, '-r')
plot(theta, rho, '-b')
hold off

set(gcf,'PaperUnits','centimeters', 'PaperSize', [42 27.9]);
set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
fcam3 = gcf;

%Plot of the normalized breathing cycle and first derivate

fnorm = figure('Name', 'Normalized Breath Cycle', 'Units', 'centimeters', 'NumberTitle', 'off');

hold on
%plot(theta, drhonorm, '-r')
plot(theta, rhonorm, '-b')
hold off

set(gcf,'PaperUnits','centimeters', 'PaperSize', [42/2 27.9/2]);
set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
fcam3 = gcf;

%Plot of a 1-cycle camshaft

for n = 1:numel(d)

    fcam1 = figure('Name', '1-Cycle Camshaft', 'Units', 'centimeters', 'NumberTitle', 'off');
    %hold on

    %polarplot(theta, drhocam, '-r')
    polar(theta, rhocam(n,:), '-b')

    %rlim([0 25]);
    hold off

    set(gcf,'PaperUnits','centimeters', 'PaperSize', [42 27.9]);
    set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
    fcam3 = gcf;

end

% %Plot of a 2-cycle camshaft
%
% fcam2 = figure('Name', '2-Cycle', 'Units', 'centimeters', 'NumberTitle', 'off');
%
% %hold on;
% %polar(theta, drhocam2, '-r')
% polar(theta2, rhocam2, '-b')
% %hold off;
%
% set(gcf,'PaperUnits','centimeters', 'PaperSize', [42 27.9]);
% set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
% fcam3 = gcf;
%
%
% %Plot of a 3-cycle camshaft
%
% fcam3 = figure('Name', '3-Cycle', 'Units', 'centimeters', 'NumberTitle', 'off');
%
% %hold on;
% %polar(theta,drhocam3)
% polar(theta3, rhocam3)
% %hold off;
%
% set(gcf,'PaperUnits','centimeters', 'PaperSize', [42 27.9]);
% set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
% fcam3 = gcf;

% The OxyGEN proyect by Profoty.xyz
% Together against the Covid-19
% V4.0 rev.17/03/2020
%
% Website/blog: oxygen.protofy.xyz
% Contact: oxygen@protofy.xyz

author: ProtofyTeam

repo: OxyGEN

save_path: github-repos/MATLAB/ProtofyTeam-OxyGEN


text:
% visualize_rbm - Visualize
% Copyright (C) 2011 KyungHyun Cho, Tapani Raiko, Alexander Ilin
%
% 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.
%
function [f] = visualize_rbm(f, R, v0, vf, h0, hf, W_grad, vbias_grad, hbias_grad);

    n_visible = floor(sqrt(size(R.W,1))).^2;

    set(0, 'CurrentFigure', f);

    subplot(3,3,1);
    visualize(R.W(1:n_visible,:), 1);
    title ('Weights');

    subplot(3,3,2);
    visualize(R.vbias(1:n_visible), 1);

    if R.adaptive_lrate.use == 1
        subplot(3,3,3);
        semilogy(R.signals.lrates);
        title('Learning rates');
    end

    subplot(3,3,4);
    visualize(v0(:,1:n_visible)', 1);
    title ('Data particles');

    subplot(3,3,5);
    visualize(vf(:,1:n_visible)', 1);
    title ('Fantasy particles');

    subplot(3,3,6);
    plot(R.signals.recon_errors);
    title ('Reconstruction errors');

    subplot(3,3,7);
    visualize(W_grad(1:n_visible,:), 1);
    title ('Gradient (Weights)');

    subplot(3,3,8);
    hist(mean(h0, 1));
    xlim([0 1]);
    title ('Average activation (data)');

    subplot(3,3,9);
    hist(mean(hf, 1));
    xlim([0 1]);
    title ('Average activation (model)');

author: kyunghyuncho

repo: deepmat


text:
function sclMask3M = perturbImageVolume(mask3M,scale)
% function sclMask3M = perturbImageVolume(mask3M,scale)
%
% Scales mask3M by amount scale.
%
% APA, 2/25/2019

sclMask3M = zeros(size(mask3M),'like',mask3M);
sizV = size(mask3M);
iCtr = round(sizV(1)/2);
jCtr = round(sizV(2)/2);
% scale = normrnd(0,2,1);
for slc = 1:size(mask3M,3)
%     sizV = size(mask3M(:,:,slc));
%     [iV,jV] = find3d(mask3M(:,:,slc));
%     maskM = mask3M(min(iV):max(iV),min(jV):max(jV),slc);
    maskM = mask3M(:,:,slc);
%     di = max(iV) - min(iV);
%     dj = max(jV) - min(jV);
%     iCtr = ceil((max(iV) + min(iV))/2);
%     jCtr = ceil((max(jV) + min(jV))/2);
%     minRow = -di*pctJitter/100;
%     maxRow = di*pctJitter/100;
%     minCol = -dj*pctJitter/100;
%     maxCol = dj*pctJitter/100;
%     numRows = minRow + (maxRow-minRow)*rand(1);
%     numCols = minCol + (maxCol-minCol)*rand(1);
%     minAng = -pctJitter/100*180;
%     maxAng = pctJitter/100*180;
%     angl = minAng + (maxAng-minAng)*rand(1);
%     scl = (100-pctJitter)/100 + 2*pctJitter/100*rand(1);
%     maskM = imtranslate(maskM,[numRows,numCols],'nearest','FillValues',0);
%     maskM = imrotate(maskM,angl,'bilinear');
    maskM = imresize(maskM, scale, 'nearest');
    newSiz = size(maskM);
    iStart = 1;
    jStart = 1;
    iEnd = 0;
    jEnd = 0;
    iMin = iCtr - ceil(newSiz(1)/2);
    if iMin < 0
        iStart = 1-iMin;
        iMin = 1;
    else
        iMin = iMin + 1;
    end
    iMax = iCtr + floor(newSiz(1)/2);
    if iMax > sizV(1)
        iEnd = iMax - sizV(1);
        iMax = sizV(1);
    end
    jMin = jCtr - ceil(newSiz(2)/2);
    if jMin < 0
        jStart = 1-jMin;
        jMin = 1;
    else
        jMin = jMin + 1;
    end
    jMax = jCtr + floor(newSiz(2)/2);
    if jMax > sizV(2)
        jEnd = jMax - sizV(2);
        jMax = sizV(2);
    end
    %mask3M(:,:,slc) = 0;
    sclMask3M(iMin:iMax,jMin:jMax,slc) = maskM(iStart:end-iEnd,jStart:end-jEnd);
end


text:
function testLandmarks

FV.vertices= double(reshape(BFM.shapeMU,3,length(BFM.shapeMU)/3)');
FV.faces = BFM.faces;

figure;
subplot(1,3,1); plotLandmarks(FV.vertices*roty(-30),faces,idx);
subplot(1,3,2); plotLandmarks(FV.vertices,faces,idx)
subplot(1,3,3); plotLandmarks(FV.vertices*roty(30),faces,idx)
function plotLandmarks(vertices,faces,idx)
patch('Vertices',vertices,'Faces',faces,'FaceColor', [1 1 1], 'EdgeColor', 'none', 'FaceLighting', 'phong'); axis equal; axis off; light;
plot3(vertices(idx,1),vertices(idx,2),vertices(idx,3),'.');
text(vertices(idx,1),vertices(idx,2),vertices(idx,3)*2,cellstr(num2str((1:length(idx))')));

author: anilbas

repo: BFMLandmarks


text:
function v=add_rotation(v,r,f,n)
% add global rotation
w=2*pi*f;
wv=[0;0;w];
for nc=1:n
    rx=r(1,nc,1);
    ry=r(1,nc,2);
    rv=[rx;ry;0];
    vv=cross(wv,rv);
    v(1,nc,1)=v(1,nc,1)+vv(1);
    v(1,nc,2)=v(1,nc,2)+vv(2);

end

author: Sable

repo: mcbench-benchmarks


text:
function test01 ( prefix )

%*****************************************************************************80
%
%% TEST01 tests XY_TO_DIST_TABLE.
%
%  Discussion:
%
%    Get the XY coordinates of a set of cities, and compute
%    the city-to-city distance table.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 September 009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string PREFIX, the common file prefix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01'
  fprintf ( 1, '  POINT_TO_DIST_TABLE computes a distance table from a'
  fprintf ( 1, '  list of point locations.'

  main_filename = sprintf ( '%s_main.txt', prefix );
  point_filename = sprintf ( '%s_points.txt', prefix );
  dist_table_filename = sprintf ( '%s_dist_table.txt', prefix );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The main filename is "%s".\n', main_filename );
  fprintf ( 1, '  The point filename is "%s".\n', point_filename );
  fprintf ( 1, '  The distance table filename will be "%s".\n', dist_table_filename );

  [ dim_num, point_num ] = r8mat_header_read ( point_filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The spatial dimension is ', dim_num
  fprintf ( 1, '  The number of points is  ', point_num

  point = r8mat_data_read ( point_file, dim_num, point_num );

  r8mat_transpose_print ( dim_num, point_num, point, '  The points:' );

  dist_table = point_to_dist_table ( dim_num, point_num, point );

  dist_table(1:n,1:n) = round ( dist_table(1:n,1:n) );

  r8mat_print_some ( point_num, point_num, dist_table, 1, 1, 5, 5, ...
    '  Initial 5x5 distance subtable:' );

  r8mat_write ( dist_table_filename, point_num, point_num, dist_table );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% $Id: publish_lagextraction_plugin.m 4 2009-08-15 21:10:35Z gramfort $
% $LastChangedBy: gramfort $
% $LastChangedDate: 2009-08-15 17:10:35 -0400 (Sam, 15 aoû 2009) $
% $Revision: 4 $

clear
opts.stylesheet = 'doc/agramfort.xsl';
opts.outputDir = 'doc';
opts.maxWidth = 600;
opts.catchError = false;
opts.evalCode = true;
% opts.evalCode = false;

% publish('lagextraction_demo',opts);
publish('lagextraction_tutorial',opts);
zip('../lagextraction_eeglab_plugin.zip','../lagextraction');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [relativearea, name] = objectareas(D)

k = 0;
for i = 1:length(D);
    if isfield(D(i).annotation, 'object')
        Nobjects = length(D(i).annotation.object);

        nrows = D(i).annotation.imagesize.nrows;
        ncols = D(i).annotation.imagesize.ncols;

        for n = 1:Nobjects
            [X,Y] = getLMpolygon(D(i).annotation.object(n).polygon);

            area = polyarea(X,Y); % ignores intersections

            k = k+1;
            relativearea(k) = area/(nrows*ncols);
            name{k} = D(i).annotation.object(n).name;
        end
    end
end

author: CSAILVision

repo: LabelMeToolbox


text:
load -ascii mushrooms.dat
mushroomsD=mushrooms';
clear mushrooms

[N, m] = size(mushroomsD);
class = 1

%rand('state',0); randn('state',0);
%abalone = abalone(:,randperm(m));

for node = 1:N
  UNI=setdiff(unique(mushroomsD(node,:)),-9999);
  for val = 1:length(UNI)
    [I,J]=find(mushroomsD(node,:)==UNI(val));
    mushroomsD(node,J)=val;
  end
end

ns = max(mushroomsD');
seul=find(ns==1);
mushroomsD=mushroomsD(setdiff(1:N,seul),:);
[N, m] = size(mushroomsD)
ns = max(mushroomsD')

Napp = ceil(m*2/3);
Ntest = m-Napp

app  = mushroomsD(:,1:Napp);size(app)
test = mushroomsD(:,Napp+1:end);size(test)

unique(app(class,:))
unique(test(class,:))

clear mushroomsD seul UNI node I J