text:
%% Copyright (C) 2014, 2016, 2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym jacobian (@var{f})
%% @defmethodx @@sym jacobian (@var{f}, @var{x})
%% Symbolic Jacobian of symbolic expression.
%%
%% The Jacobian of a scalar expression is:
%% @example
%% @group
%% syms f(x, y, z)
%% jacobian(f)
%%   @result{} (sym 1×3 matrix)
%%       ⎡∂               ∂               ∂             ⎤
%%       ⎢──(f(x, y, z))  ──(f(x, y, z))  ──(f(x, y, z))⎥
%%       ⎣∂x              ∂y              ∂z            ⎦
%% @end group
%% @end example
%%
%% @var{x} can be a scalar, vector or cell list.  If omitted,
%% it is determined using @code{symvar}.
%%
%% Example:
%% @example
%% @group
%% f = sin(x*y);
%% jacobian(f)
%%   @result{} (sym) [y⋅cos(x⋅y)  x⋅cos(x⋅y)]  (1×2 matrix)
%%
%% jacobian(f, [x y z])
%%   @result{} (sym) [y⋅cos(x⋅y)  x⋅cos(x⋅y)  0]  (1×3 matrix)
%% @end group
%% @end example
%%
%% For vector input, the output is a matrix:
%% @example
%% @group
%% syms f(x,y,z) g(x,y,z)
%% jacobian([f; g])
%%   @result{} (sym 2×3 matrix)
%%       ⎡∂               ∂               ∂             ⎤
%%       ⎢──(f(x, y, z))  ──(f(x, y, z))  ──(f(x, y, z))⎥
%%       ⎢∂x              ∂y              ∂z            ⎥
%%       ⎢                                              ⎥
%%       ⎢∂               ∂               ∂             ⎥
%%       ⎢──(g(x, y, z))  ──(g(x, y, z))  ──(g(x, y, z))⎥
%%       ⎣∂x              ∂y              ∂z            ⎦
%% @end group
%% @end example
%%
%% Example:
%% @example
%% @group
%% jacobian([2*x + 3*z; 3*y^2 - cos(x)])
%%   @result{} (sym 2×3 matrix)
%%       ⎡  2      0   3⎤
%%       ⎢              ⎥
%%       ⎣sin(x)  6⋅y  0⎦
%% @end group
%% @end example
%% @seealso{@@sym/divergence, @@sym/gradient, @@sym/curl, @@sym/laplacian,
%%          @@sym/hessian}
%% @end defmethod

function g = jacobian(f, x)

  assert (isvector(f), 'jacobian: defined only for vectors expressions')

  if (nargin == 1)
    x = symvar(f);
    if (isempty(x))
      x = sym('x');
    end
  elseif (nargin == 2)
    % no-op
  else
    print_usage ();
  end

  if (~iscell(x) && isscalar(x))
    x = {x};
  end

  cmd = { '(f, x) = _ins'
          'if not f.is_Matrix:'
          '    f = Matrix([f])'
          'G = f.jacobian(x)'
          'return G,' };

  g = pycall_sympy__ (cmd, sym(f), x);

end

%!error jacobian (sym(1), 2, 3)
%!error <only for vectors> jacobian ([sym(1) 2; sym(3) 4])

%!shared x,y,z
%! syms x y z

%!test
%! % 1D
%! f = x^2;
%! assert (isequal (jacobian(f), diff(f,x)))
%! assert (isequal (jacobian(f,{x}), diff(f,x)))
%! assert (isequal (jacobian(f,x), diff(f,x)))

%!test
%! % const
%! f = sym(1);
%! g = sym(0);
%! assert (isequal (jacobian(f), g))
%! assert (isequal (jacobian(f,x), g))

%!test
%! % double const
%! f = 1;
%! g = sym(0);
%! assert (isequal (jacobian(f,x), g))

%!test
%! % diag
%! f = [x y^2];
%! g = [sym(1) 0; 0 2*y];
%! assert (isequal (jacobian(f), g))
%! assert (isequal (jacobian(f, [x y]), g))
%! assert (isequal (jacobian(f, {x y}), g))

%!test
%! % anti-diag
%! f = [y^2 x];
%! g = [0 2*y; sym(1) 0];
%! assert (isequal (jacobian(f), g))
%! assert (isequal (jacobian(f, {x y}), g))

%!test
%! % shape
%! f = [x y^2];
%! assert (isequal (size(jacobian(f, {x y z})), [2 3]))
%! assert (isequal (size(jacobian(f, [x y z])), [2 3]))
%! assert (isequal (size(jacobian(f, [x; y; z])), [2 3]))
%! assert (isequal (size(jacobian(f.', {x y z})), [2 3]))

%!test
%! % scalar f
%! f = x*y;
%! assert (isequal (size(jacobian(f, {x y})), [1 2]))
%! g = gradient(f, {x y});
%! assert (isequal (jacobian(f, {x y}), g.'))

%!test
%! % vect f wrt 1 var
%! f = [x x^2];
%! assert (isequal (size(jacobian(f, x)), [2 1]))
%! f = f.';  % same shape output
%! assert (isequal (size(jacobian(f, x)), [2 1]))

repo: octsympy


text:
function calcAllAertsFeatures_HN(pathData,pathText,namePT,nameCT,nameROI,featType)
% -------------------------------------------------------------------------
% function calcAllSeparateTextures_HN(pathData,pathText,namePT,nameCT,nameROI,featType,scale_mat,Ng_mat)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes SEPARATE texture features for all patients, for
% all different combinations of the following texture extraction parameters:
% - Scale: Resolution at which the ROI is isotropically resampled.
% - Ng: Number of gray-levels in the quantization process.
%
% Different extraction parameters are passed as arrays or cells in the
% function in order to test all possible combinations. This function is
% used for SEPARATE scans specifically. See Ref. [1,2] and 'prepareVolume.m'
% for more details.
%
% Texture features are computed for all head and neck (HN) DICOM imaging
% data downloaded from The Cancer Imaging Archive (TCIA) website at:
% <http://dx.doi.org/xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx>, and first organized
% in a 'DATA' directory using the function readAllDICOM_HN.m.  Results are
% then saved in a folder 'TEXTURES' in the HN WORKSPACE.
% -------------------------------------------------------------------------
% REFERENCES:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% [2] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathData: Full path to the HN sData files directory.
%              --> Ex: '/myProject/WORKSPACE/DATA'
% 2. pathNonText: Full path to the HN non texture features directory.
%              --> Ex: '/myProject/WORKSPACE/FEATURES/NON_TEXTURES'
% 3. namePT: Cell of strings of all PET sData files to read
%            --> Ex: {'HGJ_001_PT.PTscan.mat';'HGJ_022_PT.PTscan.mat'}
% 4. namePT: Cell of strings of all CT sData files to read
%            --> Ex: {'HGJ_001_CT.CTscan.mat';'HGJ_022_CT.CTscan.mat'}
% 5. nameROI: Cell of strings specifying the ROI names to analyze for the
%             patients defined by "namePT" and "nameCT"
%             --> Ex: {'GTV';'GTV-P'}
% 6. featType: Either 'GTVp' for primary GTV, or 'GTVtot' for primary GTV +
%              nodal GTVs
% 7. scale_mat: Array vector specifying the different 'Scale' values to test.
%               --> Ex: [1,2,3,4,5]
% 8. Ng_mat: Array vector specifying the different 'Ng' values to test.
%            --> Ex: [8,16,32,64]
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: March 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath=pwd;

% INITIALIZATION
scans = {'PT','CT','CTorig'}; nScans = numel(scans); names = {namePT,nameCT,nameCT};
scanTypes = {'PET','CT','CT'};
nPatient = numel(namePT);

% COMPUTATION
fprintf('\n')
for i = 1:nPatient
    for j = 1:nScans
        cd(pathData)
        tStart = tic;
        if strcmp(scans{j},'CTorig')
            fprintf(['\n*********************** COMPUTING AERTS SIGNATURE: %s, ORIGINAL ***********************'],names{j}{i});
        else
            fprintf(['\n*********************** COMPUTING AERTS SIGNATURE: %s, ***********************'],names{j}{i});
        end
        scanType = scanTypes{j};
        load(names{j}{i}) % Variable 'sData' now in MATLAB Workspace;
        spatialRef = sData{2}.scan.volume.spatialRef;
        mask = zeros(size(sData{2}.scan.volume.data));
        [contourVect] = findContours_HN(sData,{nameROI{i}},scanType);
        for c = 1:numel(contourVect)
            temp = zeros(size(sData{2}.scan.volume.data));
            [~,ROImask,bb] = getROI(sData,contourVect(c));
            temp(bb(1,1):bb(1,2),bb(2,1):bb(2,2),bb(3,1):bb(3,2)) = ROImask;
            mask = mask + temp;
        end
        mask(mask > 0) = 1;
        [bb] = computeBoundingBox(mask); a = 0; b = 0; c = 0;
        sizeROI = [bb(1,2)-bb(1,1)+1,bb(2,2)-bb(2,1)+1,bb(3,2)-bb(3,1)+1];
        if mod(sizeROI(1),2)
            a = a + 1;
        end
        if mod(sizeROI(2),2)
            b = b + 1;
        end
        if mod(sizeROI(3),2)
            c = c + 1;
        end
        mask = mask(bb(1,1):bb(1,2)+a,bb(2,1):bb(2,2)+b,bb(3,1):bb(3,2)+c);
        vol = sData{2}.scan.volume.data(bb(1,1):bb(1,2)+a,bb(2,1):bb(2,2)+b,bb(3,1):bb(3,2)+c);
        if strcmp(scanType,'PET')
            [vol] = computeSUVmap(vol,sData{3}(1));
        end
        if strcmp(scans{j},'CTorig')
            [signature] = getAertsSign(vol,mask,spatialRef,scanType); % Calculating the radiomic signature as originially defined by Aerts et al.
            suffix = '_orig';
        else
            [signature] = getAertsSignMod(vol,mask,spatialRef,scanType); % Improved version (dimensionless compactness, 1 GLRLM matrix, GLN normalization)
            suffix = '';
        end
        ind = strfind(names{j}{i},'.'); namePatientScan = names{j}{i}(1:ind(1)-1);
        cd(pathText), save([namePatientScan,'_',featType,'_sign',suffix],'signature')
        time = toc(tStart);
        fprintf('TOTAL TIME: %.2f seconds\n',time)
    end
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
% IM = mkR(SIZE, EXPT, ORIGIN)
%
% Compute a matrix of dimension SIZE (a [Y X] 2-vector, or a scalar)
% containing samples of a radial ramp function, raised to power EXPT
% (default = 1), with given ORIGIN (default = (size+1)/2, [1 1] =
% upper left).  All but the first argument are optional.

% Eero Simoncelli, 6/96.

function [res] = mkR(sz, expt, origin)

sz = sz(:);
if (size(sz,1) == 1)
  sz = [sz,sz];
end

% -----------------------------------------------------------------
% OPTIONAL args:

if (exist('expt') ~= 1)
  expt = 1;
end

if (exist('origin') ~= 1)
  origin = (sz+1)/2;
end

% -----------------------------------------------------------------

[xramp,yramp] = meshgrid( [1:sz(2)]-origin(2), [1:sz(1)]-origin(1) );

res = (xramp.^2 + yramp.^2).^(expt/2);

repo: SelfExSR


text:
filename='Cantileverbeam_Tetrahedra_Linear_Structured_Fine';%Cantilever_tetrahedra';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
%cost = {'compliance','perimeter'};
%weights = [1 0.1];
cost = {'compliance'};
weights = [1];
constraint = {'volumeConstraint'};
%optimizer = 'DualNestedInPrimal';
optimizer = 'AlternatingPrimalDual';

optimizerUnconstrained = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.2;
optimality_final =1e-3;
constr_final = 1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-1;
constr_initial = 1e-1;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

plotting = false;
printing = false;
monitoring = false;
monitoring_interval = 10;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [a, b] = isintent(this,intent)
% Correspondance between fieldnames and NIfTI intent codes
% FORMAT ind = isintent(this,intent)
% this    -  GIfTI object
% intent  -  fieldnames
% a       -  indices of found intent(s)
% b       -  indices of dataarrays of found intent(s)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: isintent.m 6345 2015-02-20 12:25:50Z guillaume $

a = [];
b = [];
if ischar(intent), intent = cellstr(intent); end
for i=1:length(this(1).data)
    switch this(1).data{i}.attributes.Intent(14:end)
        case 'POINTSET'
            [tf, loc] = ismember('vertices',intent);
            if tf
                a(end+1) = loc;
                b(end+1) = i;
            end
            [tf, loc] = ismember('mat',intent);
            if tf
                a(end+1) = loc;
                b(end+1) = i;
            end
        case 'TRIANGLE'
            [tf, loc] = ismember('faces',intent);
            if tf
                a(end+1) = loc;
                b(end+1) = i;
            end
        case 'VECTOR'
            [tf, loc] = ismember('normals',intent);
            if tf
                a(end+1) = loc;
                b(end+1) = i;
            end
        case 'NODE_INDEX'
            [tf, loc] = ismember('indices',intent);
            if tf
                a(end+1) = loc;
                b(end+1) = i;
            end
        case cdata
            [tf, loc] = ismember('cdata',intent);
            if tf
                a(end+1) = loc;
                b(end+1) = i;
            end
            if strcmp(this(1).data{i}.attributes.Intent(14:end),'LABEL')
                [tf, loc] = ismember('labels',intent);
                if tf
                    a(end+1) = loc;
                    b(end+1) = i;
                end
            end
        otherwise
            fprintf('Intent %s is ignored.\n',this.data{i}.attributes.Intent);
    end
end
%[d,i] = unique(a);
%if length(d) < length(a)
%    warning('Several fields match intent type. Using first.');
%    a = a(i);
%    b = b(i);
%end

function c = cdata

c = {
'NONE'
'CORREL'
'TTEST'
'FTEST'
'ZSCORE'
'CHISQ'
'BETA'
'BINOM'
'GAMMA'
'POISSON'
'NORMAL'
'FTEST_NONC'
'CHISQ_NONC'
'LOGISTIC'
'LAPLACE'
'UNIFORM'
'TTEST_NONC'
'WEIBULL'
'CHI'
'INVGAUSS'
'EXTVAL'
'PVAL'
'LOGPVAL'
'LOG10PVAL'
'ESTIMATE'
'LABEL'
'NEURONAMES'
'GENMATRIX'
'SYMMATRIX'
'DISPVECT'
'QUATERNION'
'DIMLESS'
'TIME_SERIES'
'RGB_VECTOR'
'RGBA_VECTOR'
'SHAPE'
'CONNECTIVITY_DENSE'
'CONNECTIVITY_DENSE_TIME'
'CONNECTIVITY_PARCELLATED'
'CONNECTIVITY_PARCELLATED_TIME'
'CONNECTIVITY_CONNECTIVITY_TRAJECTORY'
};

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Y = num2strs( X, varargin )
% Applies num2str to each element of an array X.
%
% USAGE
%  Y = num2strs( X, [varargin] )
%
% INPUTS
%  X           - array of number to convert to strings
%  varargin    - [] additional input to num2str
%
% OUTPUTS
%  Y           - cell array of strings
%
% EXAMPLE
%  Y = num2strs( [1.3 2.6; 3 11] )
%
% See also NUM2STR
%
% Piotr's Computer Vision Matlab Toolbox      Version 1.5
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

Y = cell(size(X));
for i=1:numel(X)
  Y{i} = num2str( X(i), varargin{:} );
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function DEM_spatial_deconvolution
% FORMAT DEM_spatial_deconvolution
%--------------------------------------------------------------------------
% This (toy) demonstration routine illustrates spatiotemporal
% deconvolution of regional responses from imaging time-series.  The
% generative model assumes the data are generated by a small number of
% anatomical parcels that are smoothly displaced. The resulting data are
% then convolved spatially with a smoothly varying spatial kernel. The smooth
% displacement and dispersion are modelled in the usual way using discrete
% cosine basis set.  The model operates on reduced data features, using
% the eigenvariates of the original time-series - this supplements the
% implicit deconvolution with eigen-de-noising.  The ensuing estimates are
% anatomically informed because the generative model stars with a parcellation
% scheme.
%__________________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEM_spatial_deconvolution.m 6975 2016-12-18 20:27:00Z karl $

% parcellation scheme (U) (make an atlas on the fly)
%==========================================================================
rng('default')

ATLAS = spm_conv(randn(64,16),6,6);
AP    = ATLAS;
a     = spm_vec(ATLAS);
[q,x] = hist(a,4);
A     = zeros(size(a));
for i = 1:length(x)
    A(:,i) = a < x(i);
    j      = find(A(:,i));
    a(j)   = Inf;
    AP(j)  = i;
end

U.u   = A;
U.dt  = size(ATLAS);

% Specify the generative model and generate some noisy data
%==========================================================================

% parameters and priors
%--------------------------------------------------------------------------
nD    = numel(U.dt);       % number of dimensions
nB    = 2;                 % number of basis functions
nP    = nB^nD;             % number of parameters
nY    = size(U.u,2);       % number of parcels
nT    = 32;

pE.D  = 0;
pE.S  = zeros(nP,nD);
pC.D  = pE.D + 1;
pC.S  = pE.S + 1;

P.D   = 0;                % true dispersion (log precision)
P.S   = randn(nP,nD)/2;    % true dislocation
G     = randn(nY,nT);      % true time series
G     = spm_conv(G,0,2);

% gradient functions for speed (not implemented here)
%--------------------------------------------------------------------------
M.FS   = @(y,M)M.U'*y;
M.IS   = @spm_gen_image;
M.G    = @(g,M)g';

% create data
%--------------------------------------------------------------------------
y      = spm_gen_image(P,M,U)*M.G(G)';
y      = y + randn(size(y))/128;

% dimension reduction
%--------------------------------------------------------------------------
[u,s,v] = spm_svd(y);
try
    M.U = u(:,1:(nY*2));
    M.V = v(:,1:(nY*2));
catch
    M.U = u;
    M.V = v;
end

% complete data and model specification
%==========================================================================
Y.y    = y*M.V;

M.pE   = pE;
M.pC   = pC;
M.gE   = zeros(nY,size(Y.y,2));
M.gC   = speye(spm_length(M.gE));
M.hE   = 4;
M.hC   = 1/16;
M.Nmax = 32;

% invert (deconvolve)
%==========================================================================
[Ep,Eg] = spm_nlsi_N(M,U,Y);
EG      = Eg*M.V';

% reconstructed data (with no distortions or dispersion)
%--------------------------------------------------------------------------
EP    = Ep;
EP.D  = 3;
EP.S  = EP.S*0;
YQ    = spm_gen_image(EP,M,U)*M.G(EG)';

PS    = P;
PS.D  = 2;
YS    = spm_gen_image(PS,M,U)*M.G(G)';
PS.S  = PS.S*0;
YP    = spm_gen_image(PS,M,U)*M.G(G)';

% Graphics
%==========================================================================
spm_figure('GetWin','Figure 1'); clf

subplot(3,4,1), imagesc(spm_unvec(YP(:,1),ATLAS))
title('True signal','FontSize',16), axis image
subplot(3,4,2), imagesc(spm_unvec(YS(:,1),ATLAS))
title('Displaced','FontSize',16), axis image
subplot(3,4,3), imagesc(spm_unvec(y(:,1),ATLAS))
title('BOLD data','FontSize',16), axis image
subplot(3,4,4), imagesc(spm_unvec(YQ(:,1),ATLAS))
title('Deconvolved','FontSize',16), axis image

subplot(3,1,2),
plot(1:nT,EG,'-.'), hold on
plot(1:nT,G),            hold off
title('True and reconstructed time courses','FontSize',16)
xlabel('Time'), ylabel('regional signal')

subplot(3,2,5), plot(G,EG,'.')
title('Deconvolution','FontSize',16)
xlabel('True'), ylabel('predicted')

subplot(3,2,6), plot(G,(A*diag(1./sum(A)))'*y,'.')
title('Standard ROI averages','FontSize',16)
xlabel('True'), ylabel('predicted')

return

function y = spm_gen_image(P,M,U)
% generative model of image
% FORMAT y = spm_gen_image(P,M,U)
%__________________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging

% convolve
%--------------------------------------------------------------------------
y  = K(P,M,U)*U.u;

function K = K(P,M,U)
% convolution kernel
% FORMAT K = K(P)
%__________________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging

nV  = U.dt;        % number of voxels
nD  = numel(nV);   % number of dimensions
nP  = size(P.S,1); % number of parameters
nB  = nP^(1/nD);   % number of basis functions

B     = 1;
for i = 1:nD
    b = spm_dctmtx(nV(i),nB)*sqrt(nV(i));
    B = kron(B,b);
end

% dispersion and displacement
%--------------------------------------------------------------------------
D     = exp(P.D)/2;
S     = B*P.S;

% kernel
%--------------------------------------------------------------------------
nK     = size(B,1);
K      = sparse(nK,nK);
[I,J]  = ind2sub(nV,1:nK);
X(:,1) = I;
X(:,2) = J;

k     = 0;
for i = 1:nD
 d    = (X(:,i) + S(:,i))*ones(1,nK) - ones(nK,1)*X(:,i)';
 k    = k + D*(d.^2);
end

% cosine aproximation to Gaussian
%--------------------------------------------------------------------------
i     = find(abs(k) < pi);
K(i)  = (1 + cos(k(i)/pi))/2;
K     = diag(1./sum(K,2))*K;

author: spm


text:
% This subroutine assigns Parameter values for the grid
% which is used to calculate a Max Z value map.
% dx, dy is the grid spacing in degrees.
% For each one of the grid points, Ni events are counted.
%

report_this_filefun(mfilename('fullpath'));

% initial values
%
dx = 1.00;
dy = 1.00 ;
dz = 10.00 ;
ni = 100;

%
% make the interface
%
figure_w_normalized_uicontrolunits(...
    'Name','Grid Input Parameter',...
    'NumberTitle','off', ...
    'MenuBar','none', ...
    'NextPlot','new', ...
    'units','points',...
    'Visible','off', ...
    'Position',[ wex+200 wey-200 450 250]);
axis off

% creates a dialog box to input grid parameters
%
freq_field=uicontrol('Style','edit',...
    'Position',[.60 .50 .22 .10],...
    'Units','normalized','String',num2str(ni),...
    'Callback','ni=str2double(get(freq_field,''String'')); set(freq_field,''String'',num2str(ni));');

freq_field2=uicontrol('Style','edit',...
    'Position',[.60 .40 .22 .10],...
    'Units','normalized','String',num2str(dx),...
    'Callback','dx=str2double(get(freq_field2,''String'')); set(freq_field2,''String'',num2str(dx));');

freq_field3=uicontrol('Style','edit',...
    'Position',[.60 .30 .22 .10],...
    'Units','normalized','String',num2str(dy),...
    'Callback','dy=str2double(get(freq_field3,''String'')); set(freq_field3,''String'',num2str(dy));');

freq_field6=uicontrol('Style','edit',...
    'Position',[.60 .20 .22 .10],...
    'Units','normalized','String',num2str(par1),...
    'Callback','par1=str2double(get(freq_field6,''String'')); set(freq_field6,''String'',num2str(par1));');

freq_field7=uicontrol('Style','edit',...
    'Position',[.60 .60 .22 .10],...
    'Units','normalized','String',num2str(dz),...
    'Callback','dz=str2double(get(freq_field7,''String'')); set(freq_field7,''String'',num2str(dz));');

uicontrol('Units','normal','Position',[.1 .90 .15 .12],'String','LoadGrid', 'Callback','close;loadgrid')

close_button=uicontrol('Style','Pushbutton',...
    'Position',[.60 .05 .15 .12 ],...
    'Units','normalized','Callback','close;done','String','Cancel');

go_button1=uicontrol('Style','Pushbutton',...
    'Position',[.20 .05 .15 .12 ],...
    'Units','normalized',...
    'Callback',' gomakegr',...
    'String','ZmapGrid');

go_button2=uicontrol('Style','Pushbutton',...
    'Position',[.40 .05 .15 .12 ],...
    'Units','normalized',...
    'Callback','close;think; genascum',...
    'String','GenasGrid');

txt3 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0.30 0.94 0 ],...
    'Rotation',0 ,...
    'FontSize',ZmapGlobal.Data.fontsz.l ,...
    'FontWeight','bold',...
    'String',' Grid Parameter ');

txt5 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0. 0.42 0 ],...
    'Rotation',0 ,...
    'FontSize',ZmapGlobal.Data.fontsz.m ,...
    'FontWeight','bold',...
    'String','Spacing in x (dx) in deg:');

txt6 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0. 0.32 0 ],...
    'Rotation',0 ,...
    'FontSize',ZmapGlobal.Data.fontsz.m ,...
    'FontWeight','bold',...
    'String','Spacing in y (dy) in deg:');

text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0. 0.20 0 ],...
    'Rotation',0 ,...
    'FontSize',ZmapGlobal.Data.fontsz.m ,...
    'FontWeight','bold',...
    'String','Time steps in days:');

txt1 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0. 0.53 0 ],...
    'Rotation',0 ,...
    'FontSize',ZmapGlobal.Data.fontsz.m,...
    'FontWeight','bold',...
    'String','Number of Events (Ni):');

text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0. 0.63 0 ],...
    'Rotation',0 ,...
    'FontSize',ZmapGlobal.Data.fontsz.m,...
    'FontWeight','bold',...
    'String','Depth increment in (km) :');

set(gcf,'visible','on');
watchoff

author: CelsoReyes


text:
function missingStereo = checkForMissingStereo(model, trainingModel)
% USAGE:
%
%    missingStereo = checkForMissingStereo(model, trainingModel)
%
% INPUTS:
%    model:    structure with fields:
%
%                * .mets
%                * .inchi.standard
%                * .inchi.standardWithStereo
%
%    trainingModel:     structure with fields:
%
%                *.inchi.standard
%                *.inchi.standardWithStereo
%
% OUTPUTS:
%    missingStereo:

nistStdBool = false(length(trainingModel.inchi.standard));
for n = 1:length(trainingModel.inchi.standard)
   if ~isempty(trainingModel.inchi.standard{n})
       if ~any(nistStdBool(:,n))
           nistStdBool(n,n) = true;
           nistStdBool(n,strcmp(trainingModel.inchi.standard{n},trainingModel.inchi.standard)) = true;
       end
   end
end

nistStdStereoBool = false(length(trainingModel.inchi.standardWithStereo));
for n = 1:length(trainingModel.inchi.standardWithStereo)
   if ~isempty(trainingModel.inchi.standardWithStereo{n})
       if ~any(nistStdStereoBool(:,n))
           nistStdStereoBool(n,n) = true;
           nistStdStereoBool(n,strcmp(trainingModel.inchi.standardWithStereo{n},trainingModel.inchi.standardWithStereo)) = true;
       end
   end
end

nistStdStereoBool(diag(true(size(nistStdStereoBool,1),1))) = false;
nistBool = nistStdBool & ~nistStdStereoBool;
nistBool = nistBool(sum(nistBool,2)>1,:);

nistStdInchi = trainingModel.inchi.standard(any(nistBool));
modelStdInchi = model.inchi.standard;
modelStdInchi(cellfun(@isempty,modelStdInchi)) = {'N/A'};

modelStdInchiStereo = model.inchi.standardWithStereo(ismember(modelStdInchi,nistStdInchi));
modelMets = model.mets(ismember(modelStdInchi,nistStdInchi));
modelMets = regexprep(modelMets,'\[\w\]','');
[modelMets,crossi] = unique(modelMets);
modelStdInchiStereo = modelStdInchiStereo(crossi);
missingStereo = modelMets(cellfun('isempty',regexp(modelStdInchiStereo,'[tms]')));

end

author: opencobra

repo: cobratoolbox


text:
% Outputs:  (1) r: steady state interest rate
%           (2) w: steady state wage
%           (3) K: capital stock
%           (4) A: matrix for computing value function
%           (5) u: utility over grid
%           (5) c: consumption over grid
%           (6) V: value function over grid
%			(7) g: distribution
%           (8) dV_Upwind: derivative of value function by upwind scheme
%           (9) dVf: derivative of value function by forward difference
%           (10) dVb: derivative of value function by backward difference
%           (11) If: indicator for forward drift in savings
%           (12) Ib: indicator for backward drift in savings
%           (13) I0: indicator for no drift in savings

function [r,w,K,A,u,c,V,g,dV_Upwind,dVf,dVb,If,Ib,I0] = compute_steady_state()

%----------------------------------------------------------------
% Housekeeping
%----------------------------------------------------------------

% Declare global variables
global ggamma rrho ddelta aalpha ssigmaTFP rrhoTFP z lla mmu ttau I amin amax a da aa ...
	zz Aswitch rmin rmax r0 maxit crit Delta Ir crit_S zAvg A

% Initialze variables for iteration
dVf = zeros(I,2);
dvB = zeros(I,2);
c = zeros(I,2);
KS = zeros(Ir,1);
r = r0;
KD = (((aalpha) / (r + ddelta)) ^ (1 / (1 - aalpha))) * zAvg;
w = (1 - aalpha) * (KD ^ aalpha) * ((zAvg) ^ (-aalpha));
v0(:,1) = (w*mmu*(1-z(1)) + r.*a).^(1-ggamma)/(1-ggamma)/rrho;
v0(:,2) = (w*(1-ttau)*z(2) + r.*a).^(1-ggamma)/(1-ggamma)/rrho;

%----------------------------------------------------------------
% Iterate to find steady state interest rate
%----------------------------------------------------------------

for ir=1:Ir

    r_r(ir)=r;
    rmin_r(ir)=rmin;
    rmax_r(ir)=rmax;

    KD(ir,1) = (((aalpha) / (r + ddelta)) ^ (1 / (1 - aalpha))) * zAvg;
	w = (1 - aalpha) * (KD(ir) ^ aalpha) * ((zAvg) ^ (-aalpha));

    if ir>1

        v0 = V_r(:,:,ir-1);

    end

    v = v0;

    %%%%
    % Solve for value function given r
    %%%%

    for n=1:maxit

        V = v;
        V_n(:,:,n)=V;

        % Compute forward difference
        dVf(1:I-1,:) = (V(2:I,:)-V(1:I-1,:))/da;
        dVf(I,:) = (w*((1 - ttau) * z + mmu * (1 - z)) + r.*amax).^(-ggamma); %will never be used, but impose state constraint a<=amax just in case

        % Compute backward difference
        dVb(2:I,:) = (V(2:I,:)-V(1:I-1,:))/da;
        dVb(1,:) = (w*((1 - ttau) * z + mmu * (1 - z)) + r.*amin).^(-ggamma); %state constraint boundary condition

        % Compute consumption and savings with forward difference
        cf = dVf.^(-1/ggamma);
        ssf = w*((1 - ttau) * zz + mmu * (1 - zz)) + r.*aa - cf;

        % Compute consumption and savings with backward difference
        cb = dVb.^(-1/ggamma);
        ssb = w*((1 - ttau) * zz + mmu * (1 - zz)) + r.*aa - cb;

        % Compute consumption and derivative of value function for no drift
        c0 = w*((1 - ttau) * zz + mmu * (1 - zz)) + r.*aa;
        dV0 = c0.^(-ggamma);

        % Compute upwind differences
        If = ssf > 0;       %positive drift --> forward difference
        Ib = ssb < 0;       %negative drift --> backward difference
        I0 = (1-If-Ib);     %no drift
        dV_Upwind = dVf.*If + dVb.*Ib + dV0.*I0;
        c = dV_Upwind.^(-1/ggamma);
        u = c.^(1-ggamma)/(1-ggamma);
        savingsSS = w*((1 - ttau) * zz + mmu * (1 - zz)) + r.*aa - c;

        % Construct matrix for updating implicit scheme
        X = -min(ssb,0)/da;
        Y = -max(ssf,0)/da + min(ssb,0)/da;
        Z = max(ssf,0)/da;

        A1=spdiags(Y(:,1),0,I,I)+spdiags(X(2:I,1),-1,I,I)+spdiags([0;Z(1:I-1,1)],1,I,I);
        A2=spdiags(Y(:,2),0,I,I)+spdiags(X(2:I,2),-1,I,I)+spdiags([0;Z(1:I-1,2)],1,I,I);
        A = [A1,sparse(I,I);sparse(I,I),A2] + Aswitch;

        B = (1/Delta + rrho)*speye(2*I) - A;

        u_stacked = [u(:,1);u(:,2)];
        V_stacked = [V(:,1);V(:,2)];
        b = u_stacked + V_stacked/Delta;

        % Solve system of equations for updating implicit scheme
        V_stacked = B\b;

        V = [V_stacked(1:I),V_stacked(I+1:2*I)];

        % Update value function and check convergence
        Vchange = V - v;
        v = V;

        dist(n) = max(max(abs(Vchange)));
        if dist(n)<crit

            %disp('Value Function Converged, Iteration = ')
            %disp(n)
            break

        end

    end

    %%%%
    % Solve for stationary distribution
    %%%%

    % Preallocate matrices for solving linear system
    AT = A';
    b = zeros(2*I,1);

    % Normalization so pdf integrates to 1
    i_fix = 1;
    b(i_fix)=.1;
    row = [zeros(1,i_fix-1),1,zeros(1,2*I-i_fix)];
    AT(i_fix,:) = row;

    %Solve linear system for distribution
    gg = AT\b;
    g_sum = gg'*ones(2*I,1)*da;
    gg = gg./g_sum;
    g = [gg(1:I),gg(I+1:2*I)];

    % Compute objects from this iteration
    g_r(:,:,ir) = g;
    adot(:,:,ir) = w * ((1 - ttau) * zz + mmu * (1 - zz))+ r.*aa - c;
    V_r(:,:,ir) = V;

    KS(ir,1) = g(:,1)'*a*da + g(:,2)'*a*da;
    S(ir,1) = KS(ir,1) - KD(ir,1);

    % Update interest rate
    if S(ir)>crit_S

        %disp('Excess Supply')
        rmax = r;
        r = 0.5*(r+rmin);

    elseif S(ir)<-crit_S;

        %disp('Excess Demand')
        rmin = r;
        r = 0.5*(r+rmax);

    elseif abs(S(ir))<crit_S;

        display('Steady State Found, Interest rate =')
        disp(r)
        break

    end

end

% Save steady state aggregate capital stock
K = KS(ir,1);
% investment
savings = w*((1 - ttau) * zz + mmu * (1 - zz)) + r.*aa - c;
invest = sum((reshape(savings,I*2,1) .* gg)'*da);

author: gregkaplan

repo: phact

save_path: github-repos/MATLAB/gregkaplan-phact


text:
%% Copyright (C) 2014, 2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @defun mat_rccross_access (@var{A}, @var{r}, @var{c})
%% Private helper routine for symbolic array access.
%%
%% Access entries of @var{A} that are the cross product of vectors
%% @var{r} and @var{c}.  @var{r} and @var{c} could be strings.
%% Namely @code{':'}.
%%
%% @var{r} and @var{c} could contain duplicates.  This is one
%% reason why this code doesn't easily replace
%% @code{mat_rclist_access}.
%%
%% @var{z} return is small matrix with @code{len(@var{c})} columns
%% and @code{len(@var{r})} rows.
%%
%% @end defun

function z = mat_rccross_access(A, r, c)

  if (ischar(r) && ischar(c) && strcmp(r, ':') && strcmp(c, ':'))
    z = A;
    return
  end

  %if both expressible as py slices...
  % FIXME: Could optimize these cases

  [n, m] = size(A);

  if (isnumeric(r) && isempty(r))
    % no-op
  elseif (isnumeric(r) && isvector(r))
    assert(all(r >= 1) && all(r <= n), 'index out of range')
  elseif (strcmp(r, ':'))
    r = 1:n;
  elseif (islogical(r) && isvector(r) && (length(r) == n))
    I = r;
    r = 1:n;  r = r(I);
  else
    error('unknown 2d-indexing in rows')
  end

  if (isnumeric(c) && isempty(c))
    % no-op
  elseif (isnumeric(c) && isvector(c))
    assert(all(c >= 1) && all(c <= m), 'index out of range')
  elseif (strcmp(c,':'))
    c = 1:m;
  elseif (islogical(c) && isvector(c) && (length(c) == m))
    J = c;
    c = 1:m;  c = c(J);
  else
    error('unknown 2d-indexing in columns')
  end

  cmd = { '(A, rr, cc) = _ins'
          'if not A.is_Matrix:'
          '    A = sp.Matrix([A])'
          'M = [[0]*len(cc) for i in range(len(rr))]'
          'for i in range(0, len(rr)):'
          '    for j in range(0, len(cc)):'
          '        M[i][j] = A[rr[i], cc[j]]'
          '# M = sp.Matrix(M)'
          '# not quite b/c of preserving 3x0 and 0x3 for example'
          'M = sp.Matrix(len(rr), len(cc), [c for r in M for c in r])'
          'return M,' };

  rr = num2cell(int32(r-1));
  cc = num2cell(int32(c-1));
  z = pycall_sympy__ (cmd, A, rr, cc);

  % FIXME: here's some code could be used for slices
  if (1==0)
  cmd = { 'A = _ins[0]'  ...
          'r = slice(_ins[1],_ins[2])'  ...
          'c = slice(_ins[3],_ins[4])'  ...
          'M = A[r,c]'  ...
          'return M,' };
  z = pycall_sympy__ (cmd, A, r1-1, r2, c1-1, c2);
  end
end


text:
function v=edge_weight_vector(As,A,varargin)
% EDGE_WEIGHT_VECTOR Returns input for the edge_weight option for a weighted
% matrix
%
% Given a structural and weighted matrix pair (As,A), this function returns
% input for the edge_weight option.  The structural matrix As can have
% arbitrary non-zero values, but the non-zero structure of A must be a
% subset of the non-zero structure of As.  In terms of graphs, think about
% it as: the weights come from A, but the edges come from As.
%
% Example:
%   n = 8; u = 1; v = 2;
%   E = [1:n 2:n 1; 2:n 1 1:n]';
%   w = [1 zeros(1,n-1) 1 zeros(1,n-1)]';
%   A = sparse(E(:,1), E(:,2), w, n, n); % create weighted sparse matrix
%   As = sparse(E(:,1), E(:,2), true, n, n); % create structural sparse matrix
%   [d pred] = shortest_paths(As,u,struct('edge_weight',edge_weight_vector(As,A)));
%   d(v)

% David Gleich
% Copyright, Stanford University, 2008

%% History
%  2008-09-24: Initial coding
%%
[trans check] = get_matlab_bgl_options(varargin{:});
if check
    check_matlab_bgl(A,struct());
    if any(size(As)~=size(A)), error('matlab_bgl:edge_weight_vector', ...
            'The structral and weight matrix must be the same size'); end
end
n = size(A,1);
if trans, [j i] = find(As');
else [i j] = find(As);
end
Ai = sparse(i,j,1:nnz(As),n,n);
inds = nonzeros(Ai.*(A~=0));
if length(inds) ~= nnz(A), error('matlab_bgl:edge_weight_vector', ...
    'The matrix A cannot have additional non-zeros not in As.'); end
v = zeros(nnz(As),1);
v(inds) = nonzeros(A);
if check && ~isequal(A,sparse(j,i,v,n,n)')
    error('matlab_bgl:edge_weight_vector','error in weight function');
end

%% Todo
% Check that is actually works with the is trans option, I'm not convinced
% the test coverage works.

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%% Example of using KST class for interfacing with KUKA iiwa robots

% Using the gamepad for controlling the end-effector of the KUKA
% iiwa 7 R 800 robot.

% 1- Left joystick of the game pad is used to control the X-Y motion of the end effector.
%   a- Move the left joystick right and left to move the end-effector along the X axes.
%   b- Move the left joystick up or down to move the end-effector along the Y axes.
%   c- Combination is used to move the end-effector along X and Y
%   simultaneously
% -------------------------------------
% 2- Right joystick of the game pad is used to control the orientation of the end effector around the X and the Y axes.
%   a- Move the right joystick right and left to rotate the end-effector around the X axes.
%   b- Move the right joystick up or down to rotate the end-effector around the Y axes.
%   c- Combination is used to rotate the end-effector around X and Y
%   simultaneously
% -------------------------------------

% Warning: keep away from joint limits, and singularities they are not
% yet taken into consideratio.
% The Toolbox is a work in progress, other functionalities are going to be
% added with time.

% This example works with Sunrise application version KST_1.7  and higher.
% Copyright: Mohammad SAFEEA, 25th of June 2018

close all;clear;clc;
warning('off')
% Start the joystick, make sure that the joystick is connected
instrreset;
ID=1;
joy=vrjoystick(ID);
%% Add path of KST class to work space
cDir = pwd;
cDir=getTheKSTDirectory(cDir);
addpath(cDir);

%% Instantiate the KST object
ip='172.31.1.147'; % The IP of the controller
arg1=KST.LBR7R800; % choose the robot iiwa7R800 or iiwa14R820
arg2=KST.Medien_Flansch_elektrisch; % choose the type of flange
Tef_flange=eye(4); % transofrm matrix of EEF with respect to flange
iiwa=KST(ip,arg1,arg2,Tef_flange); % create the object

%% Start a connection with the server
flag=iiwa.net_establishConnection();
if flag==0
	fprintf('Can not connect to KST \n');
    fprintf('Program terminated \n');
    return;
end

% Initial configuration
jPos={0, pi / 180 * 30, 0, -pi / 180 * 60, 0,...
                        pi / 180 * 90, 0};
% Move robot to initial configuration
relVel=0.15;
iiwa.movePTPJointSpace(jPos, relVel); % move to initial configuration

% Start direct servo in joint space
iiwa.realTime_startDirectServoJoints();

% IK solver parameters
numberOfIterations=10;
lambda=0.1;
TefTool=eye(4);

% Initial configuration vector
qin=zeros(7,1);
for i=1:7
    qin(i)=jPos{i};
end

% Ttransformation matrix at initial configurarion
[Tt,j]=iiwa.directKinematics(qin);

% max angular velocity, at which the end-effector can rotate.
w=5*pi/180; % rad/sec
% max linear velocity, at which the end-effector can move.
v=0.05; % m/sec
% Joint space control
firstExecution=0;

pause(0.1);

% filtering buffer, value
c=0.9;
joyStatus=read(joy);
filter=zeros(size(joyStatus));

% Control loop
while true
    joyStatus=read(joy);
    % Remove the bias in the input analog signal
    analogPrecission=1/20;
    for i=1:4
        if abs( joyStatus(i))<analogPrecission
            joyStatus(i)=0;
        end
    end
    filter=c*filter+(1-c)*joyStatus;
    joyStatus=filter;

    % About the variable (joyStatus)
    % joyStatus: is a 4x1 vector, the first and the second elements of this
    % vector correspond to the analog values of the left analog-stick
    % psoition. the third and the fourth elements of this
    % vector correspond to the analog values of the right analog-stick
    % psoition.
    %
    % All foour elemets of the analog signal of the vector  (joyStatus) are used for
    % controlling the end-effector.
    %

    % Calculate the elapsed time between updates
    if(firstExecution==0)
        firstExecution=1;
        a=datevec(now);
        timeNow=a(6)+a(5)*60+a(4)*60*60; % calculate time at this instant
        dt=0; % elapsed is zero at first excution
        time0=timeNow;
    else
        a=datevec(now);
        timeNow=a(6)+a(5)*60+a(4)*60*60; % calculate time at this instant
        dt=timeNow-time0; % calculate elapsed time interval
        time0=timeNow;
    end
 % Construct motion control command
     n=[joyStatus(4);
         joyStatus(3);
         0];
     k=[0;0;1];
     w_control=w*cross(n,k);
     v_control=v*[joyStatus(2);
         joyStatus(1);
         0];
% Calculate target transform using command, and current transform
     Tt=updateTransform(Tt,w_control,v_control,dt);

     for i=1:7
         qin(i)=jPos{i};
     end

      qt=kukaDLSSolver( qin, Tt, TefTool,numberOfIterations,lambda );

      if(size(qt)==[7,1])
          errorFlag=false;
          for i=1:7
              if abs(qt(i)-qin(i))>pi*0.5/180
                  errorFlag=true;
              end
          end
      % check joints limits
     qt=iiwa.gen_InverseKinematics( qin, Tt,numberOfIterations,lambda );
     flag=jointLimitReached(qt);
     % if the motion causes joint limit violation, stay in the current configuration
     if flag==true;
         qt=qin;
     end

          if errorFlag==false
              for i=1:7
                  jPos{i}=qt(i);
              end
          end
      end
      iiwa.sendJointsPositions(jPos);

     % A condition to break the loop, when x button of the gamepad is
     % pressed
     b=button(joy,1);
     if(b==1)
         break;
     end

end
% turn off the server
% close connection
iiwa.net_turnOffServer();

repo: KST-Kuka-Sunrise-Toolbox


text:
function output = BalancePrep(meta_meta)
% Balances the reactions, used in `ReactionChecks`
%
% USAGE:
%
%    output = BalancePrep(meta_meta)
%
% INPUT:
%    meta_meta:    matrix from `ReconstructionTool`
%
% OUTPUT:
%    output:   `leftside`, `rightside`, `charge_l` and `charge_r` in `varargout`
%
% .. Author: - Stefan G. Thorleifsson 2010
%
% .. rBioNet is published under GNU GENERAL PUBLIC LICENSE 3.0+
% .. Thorleifsson, S. G., Thiele, I., rBioNet: A COBRA toolbox extension for
% .. reconstructing high-quality biochemical networks, Bioinformatics, Accepted.
% .. rbionet@systemsbiology.is

output = [];
if isempty(meta_meta)
    return
end

S = size(meta_meta);

for i = 1:S(1);
    if ismember('',meta_meta(i,6))
        msgbox(['Reaction cannot be created with metabolites that are'...
            'missing the charged formula. ' meta_meta(i,1)  ],...
            'Missing charge formula','warn');
        return
    end
end
%---------- Creating formula ---------------- begin

%initialize variables
leftside = []; rightside = [];
charge_left = []; charge_right = [];

for i = 1:S(1)

    %creating metabolite [compartment] ready
    comp_str = meta_meta{i,4};
    comp_cnt = regexpi(comp_str,'\(');
    comp = comp_str(comp_cnt+1);

    if meta_meta{i,3} == 1 % No numbers
        newmetab = [meta_meta{i,1} '[' comp ']'];
    else
        newmetab = [num2str(meta_meta{i,3}) ' ' meta_meta{i,1} '[' comp ']'];
    end

    % Dividing left from right.
    if strmatch(meta_meta(i,5),'Substrate','exact')
        charge_left = [charge_left str2double(meta_meta{i,7})*meta_meta{i,3}]; %adding up the charge

        if isempty(leftside)
            leftside = newmetab;
        else
            leftside = [leftside ' + ' newmetab];
        end
    elseif strmatch(meta_meta(i,5),'Product','exact')
        charge_right = [charge_right str2double(meta_meta{i,7})*meta_meta{i,3}];%adding up the charge
        if isempty(rightside)
            rightside = newmetab;
        else
            rightside = [rightside ' + ' newmetab];
        end
    else
        msgbox('Something is wrong with Substrate and Product.',...
            'Reaction balance.','error');
        return
    end
end

charge_l = sum(charge_left);
charge_r = sum(charge_right);

output = {leftside,rightside,charge_l,charge_r};


text:
function tests = test_ft_connectivity_mim

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_connectivity_mim

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  func = localfunctions;
  for i=1:numel(func)
    fprintf('evaluating %s\n', func2str(func{i}));
    feval(func{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function test_chan_chan_freq(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nchan   = 6;
nfreq   = 5;
dimord  = 'chan_chan_freq';
input   = zeros(nchan, nchan, nfreq);

for freq=1:nfreq
  fdat = 10 * randn(nchan, 1) + 1i * 10 * randn(nchan, 1);
  input(:,:,freq) = fdat * ctranspose(fdat); % compute the cross-spectral density
end

result = {};
result{end+1} = ft_connectivity_cancorr(input, 'dimord', dimord, 'indices', [1 1 2 2 3 3]); % e.g. three pair of planar channels
result{end+1} = ft_connectivity_cancorr(input, 'dimord', dimord, 'indices', [1 1 1 2 2 2]); % e.g. two dipole moments (xyz)

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequaln(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end


text:
% std_comppol() - inverse component polarity in a component cluster
%
% Usage: [compout pol] = std_comppol(compin);
%
% Inputs:
%    compin  - component scalp maps, one per column.
%
% Outputs:
%    compout - component scalp maps some of them with inverted
%              polarities, one per column.
%    pol     - logical vector of component with inverted
%              polarities (same length as the number of rows in
%              compin)
%
% Author: Arnaud Delorme & Hilit Serby, SCCN, INC, UCSD, 2004

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

function [compin, pol] = std_comppol(compin);

if nargin < 1
    help std_comppol;
    return;
end;

% remove the NaN
% --------------
for index = 1:size(compin,2)
    compin(isnan(compin(:,index)),:) =[];
end;

% run several iterations
% ----------------------
pol     = ones(1,size(compin,2));
for repeat=1:3
    compave = mean(compin,2);
    for index = 1:size(compin,2)

        % remove diagonal and put 0 and 1
        % -------------------------------
        if ~all(compin(:,index) == 0)
             r = corrcoef(compave, compin(:,index) );
        else r = zeros(2,2);
        end;

        % invert component polarities
        % ---------------------------
        if r(2) < 0
            compin(:,index) = -compin(:,index);
            pol(index)      = -pol(index);
        end;
    end;
end;

author: PatternRecognition

repo: OpenBMI


text:
classdef BatchNorm < dagnn.ElementWise
    properties
        numChannels
        epsilon = 1e-3
        opts = {'NoCuDNN'} % ours seems slightly faster
    end

    properties (Transient)
        moments
    end

    methods
        function outputs = forward(obj, inputs, params)
                outputs{1} = vl_nnbnorm(inputs{1}, params{1}, params{2}, ...
                    'moments', params{3}, ...
                    'epsilon', obj.epsilon, ...
                    obj.opts{:}) ;
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
            [derInputs{1}, derParams{1}, derParams{2}, derParams{3}] = ...
                vl_nnbnorm(inputs{1}, params{1}, params{2}, derOutputs{1}, ...
                'epsilon', obj.epsilon, ...
                'moments', params{3}, ...
                obj.opts{:}) ;
            obj.moments = [] ;
            % multiply the moments update by the number of images in the batch
            % this is required to make the update additive for subbatches
            % and will eventually be normalized away
            derParams{3} = derParams{3} * size(inputs{1},4) ;
        end

        % ---------------------------------------------------------------------
        function obj = BatchNorm(varargin)
            obj.load(varargin{:}) ;
        end

        function params = initParams(obj)
            params{1} = ones(obj.numChannels,1,'single') ;
            params{2} = zeros(obj.numChannels,1,'single') ;
            params{3} = zeros(obj.numChannels,2,'single') ;
        end

        function attach(obj, net, index)
            attach@dagnn.ElementWise(obj, net, index) ;
            p = net.getParamIndex(net.layers(index).params{3}) ;
            net.params(p).trainMethod = 'average' ;
            net.params(p).learningRate = 0.1 ;
        end
    end
end

author: layumi

repo: Image-Text-Embedding


text:
function coeffsRet=polyIntMultiDim(coeffs,intDim,intConst)
%%POLYINTMULTIDIM Compute the indefinite integral of a multivariate
%                 polynomial with respect to a particular dimension given a
%                 hypermatrix of its coefficients (including cross terms).
%
%INPUTS: coeffs A hypermatrix of the coefficients for the multivariate
%               polynomial. These are arranged such that
%               coeffs(a1,a2,a3...an) corresponds to the coefficient of an
%               x1^(a1-1)*x2^(a2-1)*x3^(a3-1)...xn^(an-1) term.  Thus, the
%               number of indices coeffs takes is equal to the
%               dimensionality of x (not counting singleton dimensions at
%               the end of coeffs). Note that this ordering is the reverse
%               that used in the 1D polyval function that is built into
%               Matlab. The number of elements for each index in coeffs is
%               the maximum order of that dimension +1.
%        intDim The dimension of x (index of coeffs) with respect to which
%               the integral is to be taken.
%      intConst The integration constant. If this parameter is omitted or
%               an empty matrix is passed, a default value of 0 is used.
%
%OUTPUTS: coeffsRet The indefinite integral of the coeffs hypermatrix.
%
%This function is just implemented using the basic rules of polynomial
%integration.
%
%As an example, consider the multivariate polynomial:
%14+3*x1^2-18*x2+12*x1*x3-3*x2*x3
%The integrals with respect to the first, second and third variables are
% coeffs=zeros(3,2,2);
% coeffs(0+1,0+1,0+1)=14;
% coeffs(2+1,0+1,0+1)=3;
% coeffs(0+1,1+1,0+1)=-18;
% coeffs(1+1,0+1,1+1)=12;
% coeffs(0+1,1+1,1+1)=-3;
% intCoeffs1=polyIntMultiDim(coeffs,1);
% intCoeffs2=polyIntMultiDim(coeffs,2);
% intCoeffs3=polyIntMultiDim(coeffs,3);
%One will see that intCoeffs1 corresponds to
%14*x1+x1^3-18*x1*x2+6*x1^2*x3-3*x1*x2*x3
%intCoeffs2 corresponds to
%14*x2*+3*x1^2*x2-9*x2^2+12*x1*x2*x3-(3/2)*x2^2*x3
%and intCoeffs3 corresponds to.
%14*x3+3*x1^2*x3-18*x2*x3+6*x1*x3^2-(3/2)*x2*x3^2
%
%Note that the fucntion can be used with variables that are not present in
%the original polynomial. In the above example,
% intCoeffs4=polyIntMultiDim(coeffs,4);
%provides coefficients corresponding to
%14*x4+3*x1^2*x4-18*x2*x4+12*x1*x3*x4-3*x2*x3*x4
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(intConst))
   intConst=0;
end

dimSizeList=size(coeffs);
numIdx=length(dimSizeList);
totalEls=prod(dimSizeList);

%If the dimension of integration is a singleton dimension that has been
%suppressed from the end of dimSizeList, then we must add in all of the
%missing singleton dimensions.
dimRetSizeList=[dimSizeList,ones(1,intDim-numIdx)];
%Integration will increase the maximum degree of the term over which
%integration is performed by one, meaning that the matrix has to be made
%larger.
dimRetSizeList(intDim)=dimRetSizeList(intDim)+1;

%Allocate space for the return coefficients.
coeffsRet=zeros(dimRetSizeList);

%We will now loop through all of the values in the original coefficient
%set, adjusting them and inserting them into the proper location in the new
%modified coefficient set, after multiplication by the appropriate
%constant.

if(intDim>numIdx)
%This is used when integrating over a trailing singleton dimension.
    indicesFull=ones(numIdx,1);
    indicesFull(intDim)=2;
end
for curEl=1:totalEls
    curCoeff=coeffs(curEl);
    indices=index2NDim(dimSizeList,curEl);

    if(intDim<=numIdx)
        constVal=1/indices(intDim);
        indices(intDim)=indices(intDim)+1;
        coeffsRet(nDim2Index(dimRetSizeList,indices))=constVal*curCoeff;
    else%If we are integration over a trailing singleton dimension.
        indicesFull(1:numIdx)=indices;
        coeffsRet(nDim2Index(dimRetSizeList,indicesFull))=curCoeff;
    end
end

%Add in the integration constant. It is the zeroth-order term.
coeffsRet(1)=intConst;

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function d = ComputeFlyAngleSubtended(obj,n,fly1)

nflies = obj.nfliespermovie(n);
flyidx1 = obj.getFlyIdx(n,fly1);

% initialize
x_mm1 = obj.GetPerFrameData('x_mm',n,fly1);
y_mm1 = obj.GetPerFrameData('y_mm',n,fly1);
a_mm1 = obj.GetPerFrameData('a_mm',n,fly1);
b_mm1 = obj.GetPerFrameData('b_mm',n,fly1);
theta_mm1 = obj.GetPerFrameData('theta_mm',n,fly1);
nframes1 = obj.nframes(flyidx1);
firstframe1 = obj.firstframes(flyidx1);
endframe1 = obj.endframes(flyidx1);
d = nan(nflies,nframes1);

for fly2 = 1:nflies,
  if fly2 == fly1, continue; end

  flyidx2 = obj.getFlyIdx(n,fly2);
  firstframe2 = obj.firstframes(flyidx2);
  endframe2 = obj.endframes(flyidx2);

  % get start and end frames of overlap
  t0 = max(firstframe1,firstframe2);
  t1 = min(endframe1,endframe2);

  % no overlap
  if t1 < t0, continue; end

  % indices for these frames
  offi = firstframe1-1;
  offj = firstframe2-1;

  x_mm2 = obj.GetPerFrameData('x_mm',n,fly2);
  y_mm2 = obj.GetPerFrameData('y_mm',n,fly2);
  theta_mm2 = obj.GetPerFrameData('theta_mm',n,fly2);
  a_mm2 = obj.GetPerFrameData('a_mm',n,fly2);
  b_mm2 = obj.GetPerFrameData('b_mm',n,fly2);

  % distance from fly1's nose to fly2
  for t = t0:t1,
    i = t - offi;
    j = t - offj;

    d(fly2,i) = anglesubtended(...
      x_mm1(i),y_mm1(i),a_mm1(i),b_mm1(i),theta_mm1(i),...
      x_mm2(j),y_mm2(j),a_mm2(j),b_mm2(j),theta_mm2(j),...
      obj.fov);

  end

end

author: kristinbranson

repo: JAABA


text:
function [ state ] =movePTP_ConditionalTorque_ArcYZ_AC(t_Kuka,theta,c,VEL,joints_indices,max_torque,min_torque)
%% This function is used for moving the endeffector on an arc in the YZ plane, for the KUKA iiwa 7 R 800.
% this motion is interruptible when one or more of the joints-torques exceed the predefined limits

%% Syntax:
% [ state ] =movePTP_ConditionalTorque_ArcYZ_AC(t_Kuka,theta,c,VEL,joints_indices,max_torque,min_torque)

%% About:
% This function is used to move the end-effector on an arc in the YZ plane,

%% Arreguments:
% t_Kuka: is the TCP/IP connection
% theta: is the arc angle, in radians
% c: the YZ coordinates of the center of the circle, it is 1x2 vector.
% VEL : is a double, defines the motion velocity mm/sec.
% joints_indices: is a vector of the indices of the joints where the
% torques limits are to be imposed, the joints are indexed starting from one.
% max_torque: is a vector of the maximum torque limits for the joints
% specified in the (joints_indices) vector.
% min_torque: is a vecotr of the minimum torque limits for the joints
% specified in the (joints_indices) vector.

%% Return value:
% state: a number equals to one if the motion is completed sccessfully or a
% zero if the motion was interrupted due to external contact with the
% robot.
% if there is an error, the return value is minus one

% Copy right, Mohammad SAFEEA, 9th of April 2018

if((size(VEL,1)==1)&&(size(VEL,2)==1))
else
    disp('Error, velocity is a double and shall not be an array');
    state=-1;
    return;
end

k=[1;0;0];
pos=getEEFPos( t_Kuka );
c=colVec(c);
c1=[pos{1};c(1);c(2)];
[ state ] =movePTP_ConditionalTorque_Arc_AC(t_Kuka,theta,c1,k,VEL,joints_indices,max_torque,min_torque);

end

function [ y ] = colVec( v)
% Convert a vector to a column vector:
    if(size(v,2)==1)
        y=v;
    else
        y=v';
    end
end


text:
% CORRECT_MC - compute an upper limit for the number of independent
%                time-frequency estimate in a given time-frequency image.
%                This number can be used to correct for multiple comparisons.
%
% Usage:
%   [ncorrect array] = correct_mc( EEG, cycles, maxfreq, timesout);
%
% Inputs:
%    EEG       - EEGLAB structure
%    cycles    - [float] same as the cycle input to TIMEF. Default is [3 0.5].
%    freqrange - [float] minimum and maximum frequency. Default is [2 50] Hz.
%    timesout  - [integer] array of number of time points to test.
%
% Output:
%    ncorrect - number of independent tf estimate in the time-freq image
%    array    - array of size (freqs x timesout) containing pvalues.
%
% Method details:
%
% Dividing by the total number of time-frequency estimate in the 2-D
% time-frequency image decomposition would be too conservative since
% spectral estimates of neighboring time-frequency points are highly
% correlated. One must thus estimate the number of independent
% time-frequency points in the TF image. Here, I used geometrical wavelets
% which are optimal in terms of image compression, so neighboring
% frequencies can be assume to carry independent spectral estimates.
% We thus had time-frequency decompositions at only X frequencies (e.g. 120,
% 60, 30, 15, 7.5, 3.25, 1.625 Hz). For each frequency, I then found
% the minimum number of time points for which there was a significant
% correlation of the spectral estimates between neighboring time points
% (for each frequency and number of time point, I computed the correlation
% from 0 to 1 for all data channel to obtain an a probability distribution
% of correlation; we then fitted this distribution using a 4th order curve
% (Ramberg, J. S., E. J. Dudewicz, et al. (1979). "A probability
% distribution and its uses in fitting data." Technometrics 21(2)) and
% assessed the probability of significance for the value 0 (no correlation)
% to be within the distribution of estimated correlation). For instance,
% using 28 time points at 120 Hz, there was no significant (p>0.05 taking
% into account Bonferoni correction for multiple comparisons) correlation
% between neighboring time-frequency power estimate, but there was a
% significant correlation using 32 time points instead of 28 (p<0.05).
% Applying the same approach for the X geometrical frequencies and summing
% the minimum number of time points for observing a significant correlation
% at all frequencies, ones obtain in general a number below 200 (with the
% defaults above and 3-second data epochs) independent estimates. In all
% the time-frequency plots, one has to used a significance mask at p<0.00025
% (0.05/200). An alternative method for correcting for multiple comparisons
% is presented in Nichols & Holmes, Human Brain Mapping, 2001.
%
% Author: Arnaud Delorme, SCCN, Jan 17, 2004

% Copyright (C) 2004 Arnaud Delorme, SCCN, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [ncorrect, pval] = correct_mc( EEG, cycles, freqrange, timesout);

    if nargin < 1
        help correct_mc;
        return;
    end
    if nargin < 2
        cycles  = [3 0.5];
    end
    if nargin < 3
        freqrange = [2 50];
    end
    if nargin < 4
        % possible number of time outputs
        % -------------------------------
        timesout = [5 6 7 8 9 10 12 14 16 18 20 24 28 32 36 40];
    end
    nfreqs = ceil(log2(freqrange(2)));

    % scan times
    % ----------
    for ti = 1:length(timesout)
        clear tmpf

        % scan data channels
        % ------------------
        for index = 1:EEG.nbchan

            clf; [ersp,itc,powbase,times,freqs,erspboot,itcboot] = newtimef(EEG.data(index,:),EEG.pnts, ...
                             [EEG.xmin EEG.xmax]*1000,EEG.srate, cycles, 'timesout', timesout(ti), ...
                             'freqscale', 'log', 'nfreqs', nfreqs, 'freqrange', freqrange, 'plotitc', 'off', 'plotersp', 'off');

            % compute correlation
            % -------------------
            for fi = 1:length(freqs)
                tmp      = corrcoef(ersp(fi,1:end-1), ersp(fi,2:end));
                tmpf(index,fi) = tmp(2,1);
            end

        end

        % fit curve and determine if the result is significant
        % ----------------------------------------------------
        for fi = 1:length(freqs)
            pval(fi, ti) = rsfit(tmpf(:,fi)', 0);
            if pval(fi,ti) > 0.9999, pval(fi,ti) = NaN; end
        end
    end

    % find minimum number of points for each frequency
    % ------------------------------------------------
    ncorrect = 0;
    threshold = 0.05 / prod(size(pval));
    for fi = 1:size(pval,1)
        ti = 1;
        while ti <= size(pval,2)
            if pval(fi,ti) < threshold
                ncorrect = ncorrect +  timesout(ti);
                ti = size(pval,2)+1;
            end
            ti = ti+1;
        end
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function sst = sort_sst(sst)
%SORT_SST: Sort SST (Event Start/Stop Times) based on start times
%
%USAGE: sst = sort_sst(sst)
%
%INPUTS: sst --> nx2 Start/Stop Time (datenum) array (unsorted)
%
%OUTPUTS: N --> nx2 Start/Stop Time (datenum) array (sorted)
%
% See also ADD_SST, CHK_T, COMPARE_SST, DELETE_SST, EXTRACT_SST, IS_SST,
%          ISTEQUAL, MERGE_SST, SEARCH_SST, SORT_SST, NAN2SST, SSD2SST,
%          SST2NAN, SST2SSD, SST2VAL, SST2WFA, WFA2SST
%
% Author: Dane Ketner, Alaska Volcano Observatory
% $Date$
% $Revision$

[S N] = sort(sst(:,1));
sst = [S; sst(N,2)];

author: geoscience-community-codes

repo: GISMO


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

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

  if ( mod ( floor ( t ), 2 ) == 0 )
    jac(1,1) = - 1.5;
  else
    jac(1,1) = - 0.5;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% DIRECT JACOBIAN DEMO

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE 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.
%
% ARTE 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 ARTE.  If not, see <http://www.gnu.odrg/licenses/>.

close all;

fprintf('\nThe demo shows how to compute the end effectors speed as a function of the joint speeds')

robot =  load_robot('example','scara');

T=1; %seconds. Tiempo que dura el movimiento

q1=0:0.01:pi/2;
q2=0:0.01:pi/2;

q_v=pi/2/T; %rad/s
V=zeros(3, length(q1));
for i=1:length(q1),
   v = compute_end_velocity(robot, [q1(i) -q2(i) 0 0], [q_v q_v 0 0]);
   V(:,i) = v(1:3);
end

figure, hold
plot(V(1, :), 'r')
plot(V(2, :), 'g')
plot(V(3, :), 'b')
legend('V_x', 'V_y', 'V_z')
title('End effector speed (m/s)')
xlabel('time (s)')

figure, hold
plot(q1, 'r')
plot(q2, 'g')
legend('q_1', 'q_2')
title('Joint positions (rad)')
xlabel('time (s)')

repo: ARTE


text:
function DrawAllJoints(j)
global uLINK
radius    = 0.02;
len       = 0.06;
joint_col = 0;

if j ~= 0
    if ~isempty(uLINK(j).vertex)
        vert = uLINK(j).R * uLINK(j).vertex;
        for k = 1:3
            vert(k,:) = vert(k,:) + uLINK(j).p(k); % adding x,y,z to all vertex
        end
        DrawPolygon(vert, uLINK(j).face,0);
    end

    hold on

    i = uLINK(j).mother;
    if i ~= 0
        Connect3D(uLINK(i).p,uLINK(j).p,'k',2);
    end
    DrawCylinder(uLINK(j).p, uLINK(j).R * uLINK(j).a, radius,len, joint_col);

    DrawAllJoints(uLINK(j).child);
    DrawAllJoints(uLINK(j).sister);
end

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
function s_mm=rmNeighborsCompare(vw,model)
%  rmNeighborsCompare - compute pRF size in mm weighted by variance
%  explained
%
% 2009/02: SD & BH wrote it.

if ~exist('vw','var') || isempty(vw), error('Need view struct'); end
if ~exist('model','var') || isempty(model), error('Need rm model file'); end

% get gray connection structure
gNodes = viewGet(vw,'nodes');
gEdges = viewGet(vw,'edges');
coords = viewGet(vw,'coords');

% load model params
x = rmGet(model,'x');
y = rmGet(model,'y');
s = rmGet(model,'sigma');
ve = rmGet(model,'varexp');

% output
s_mm = zeros(size(s));

nGnodes=length(gNodes);

tic;
fprintf('[%s]:Computing...',mfilename);
for t=1:nGnodes % for each gNode...
    % Find its edges (the nodes of the things that it's connected to...)
    thisOffset=gNodes(5,t);
    thisNumEdges=gNodes(4,t);
    theseEdges=gEdges(thisOffset:(thisOffset-1+thisNumEdges)); %thisoffset-1 or 0?

    % variance explained for the neighbors
    ven = ve(theseEdges);
    ven = ven./sum(ven);

    % compute cortical distance from neighbors
    cdist = coords(:,theseEdges);
    cdist = cdist - (coords(:,t)*ones(1,size(cdist,2)));
    cdist = sum(sqrt(sum(cdist.^2)).*ven);    % distance

    % compute distance from neighbors in visual field
    vfdist = [x(theseEdges); y(theseEdges)];
    vfdist = vfdist - ([x(t); y(t)]*ones(1,size(vfdist,2)));
    vfdist = sum(sqrt(sum(vfdist.^2)).*ven);

    % compute cortical pRF size
    s_mm(t) = s(t) * (cdist./vfdist);
end

% some are not finite when vfdist == 0, we interpolate these values
% we do this only for spurious voxels, large patches will be set to global
% mean
ii = find(~isfinite(s_mm));
for n=1:5
    for t=ii
        % Find its edges (the nodes of the things that it's connected to...)
        thisOffset=gNodes(5,t);
        thisNumEdges=gNodes(4,t);
        theseEdges=gEdges(thisOffset:(thisOffset-1+thisNumEdges)); %thisoffset-1 or 0?

        % lookup neighboring values with data
        nb = s_mm(theseEdges);
        nb = nb(isfinite(nb));

        if ~isempty(nb)
            s_mm(t) = mean(nb);
        end
    end
    ii = find(~isfinite(s_mm));
end

% any left we set to the global mean
ii = ~isfinite(s_mm);
s_mm(ii) = mean(s_mm(~ii));

fprintf('Done[%dsecs].\n',round(toc));

% rmGet(model,'s_mm');

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef TimeUnitEnum < matlab.mixin.SetGet
    %TimeUnitEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        Seconds('Seconds')
        Minutes('Minutes')
        Hours('Hours')
        Days('Days')
        Years('Years')
    end

    properties
        name(1,:) char
    end

    methods
        function obj = TimeUnitEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function [listBoxStr, m] = getListboxStr()
            m = enumeration('TimeUnitEnum');
            listBoxStr = {m.name};
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function analyzePhase()
addpath('phase_feature_extraction');
fs = 16000;
if 1
    clean_root = 'F:\Data\ReverbChallenge\REVERB_WSJCAM0_et\data\cln_test\secondary_microphone\si_et_2\c3e';
    reverb_root = 'F:\Data\ReverbChallenge\REVERB_WSJCAM0_et\data\far_test\secondary_microphone\si_et_2\c3e';
    dereverb_root = 'F:\Data\ReverbChallenge\WavOutDNN_bigDNN19fr\REVERBWSJCAM0\SimData_et_for_1ch_far_room3_A';
    uttID = 'c3ec020a';

    wav_clean = audioread([clean_root '\' uttID '.wav']);
    wav_noisy = audioread([reverb_root '\' uttID '_ch1.wav']);
%     [grp_phase_clean, cep_clean] = modified_group_delay_feature([clean_root '\' uttID '.wav'], 0.4, 0.9, 12);
%     [grp_phase_noisy, cep_noisy] = modified_group_delay_feature([reverb_root '\' uttID '_ch1.wav'], 0.4, 0.9, 12);
else
    wav_clean = wavread('E:\Workspace2\Data\AURORA4\WAV\test_clean_wv1\440_16k\440c020a.wv1');
    wav_noisy = wavread('E:\Workspace2\Data\AURORA4\WAV\test_train_wv1\440_16k\440c020a.wv1');
end

frame_shift = 100/fs;
frame_length = 0.02;
[~,FT_clean] = wav2abs(wav_clean, fs, frame_shift, frame_length);
[~,FT_noisy] = wav2abs(wav_noisy, fs, frame_shift, frame_length);
FT_clean = FT_clean(1:257,:);
FT_noisy = FT_noisy(1:257,:);

mag_clean = abs(FT_clean);
mag_noisy = abs(FT_noisy);
phase_clean = angle(FT_clean);
phase_noisy = angle(FT_noisy);

[instan_freq_clean] = comp_instan_freq(phase_clean);
[instan_freq_noisy] = comp_instan_freq(phase_noisy);

[BPD_clean, IF_clean] = comp_BPD(phase_clean, 512, frame_shift*fs);
[BPD_noisy, IF_clean] = comp_BPD(phase_noisy, 512, frame_shift*fs);

[group_delay_clean] = comp_group_delay(phase_clean);
[group_delay_noisy] = comp_group_delay(phase_noisy);

[mgd_clean, log_mgd_clean] = modified_group_delay_feature(wav_clean, fs, 1, 1);
[mgd_noisy, log_mgd_noisy] = modified_group_delay_feature(wav_noisy, fs, 1, 1);

subplot(4,2,1); imagesc(log(mag_clean));
subplot(4,2,3); imagesc(log_mgd_clean);
subplot(4,2,5); imagesc(instan_freq_clean);
subplot(4,2,7); imagesc(BPD_clean);

subplot(4,2,2); imagesc(log(mag_noisy));
subplot(4,2,4); imagesc(log_mgd_noisy);
subplot(4,2,6); imagesc(instan_freq_noisy);
subplot(4,2,8); imagesc(BPD_noisy);

end

author: singaxiong

repo: SignalGraph


text:
function value = s_scrabble_points ( s )

%*****************************************************************************80
%
%% S_SCRABBLE_POINTS returns the Scrabble point value of a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string.
%
%    Output, integer S_SCRABBLE_POINTS, the point value of
%    the string.
%
  s_length = length ( s );

  value = 0;
  for i = 1 : s_length
    value = value + ch_scrabble_points ( s(i:i) );
  end

  return
end


text:
function [Pro,Res] = interpolationAMGn(A,isC)
%% INTERPOLATIONAMGS construct prolongation and restriction matrices
%
% [Pro,Res] = interpolatoinAMGs(A,isC) construct prolongation and
% restriction matrices use standard matrix-dependent interpolation.

% In the input, A is a SPD matrix and isC is a logical array to indicate
% nodes in coarse matrix. In the output Pro and Res are prolongation and
% restriction matrices satisfying Res = Pro'.
%
% Example
%   load lakemesh
%   A = assemblematrix(node,elem);
%   [isC,As] = coarsenAMGc(A);
%   [Ac,Pro,Res] = interpolatoinAMGs(As,isC);
%
% See also: coarsenAMGc, amg
%
% New interpolation added by Xiaozhe Hu. 12/10/2011.

N = size(A,1);

%% Index map between coarse grid and fine grid
allNode = (1:N)';
fineNode = allNode(~isC);
Nf = length(fineNode);
Nc = N-length(fineNode);
coarseNode = (1:Nc)';    % coarse node index
coarse2fine = find(isC); % coarse node index in the fine grid
% fine2coarse(isC) = coarseNode;
ip = coarse2fine;  % coarse node index in the fine grid
jp = coarseNode;   % coarse node index
sp = ones(Nc,1);   % identity matrix for coarse nodes in the fine grid

%% Construct prolongation and restriction operator
Afc = A(fineNode, coarse2fine);
Dsum = spdiags(A(fineNode,fineNode),0);
alpha = (sum(A(fineNode,:),2) - Dsum)./sum(Afc,2);
Dsum = spdiags(((-1).*alpha)./Dsum, 0, Nf, Nf);
[ti,tj,tw] = find(Dsum*Afc);
ip = [ip; fineNode(ti)];  % fine node index
jp = [jp; tj];            % coarse node index
sp = [sp; tw];            % weight
Pro = sparse(ip,jp,sp,N,Nc);
Res = Pro';

repo: ifem


text:
function im = fftc(d)
% Function performs a centered fft
im = fftshift(fft(fftshift(d)));


text:
%% SVMSGD Interactive Classification
% Train with SVMSGD algorithm.
% The classifier can handle linearly separable 2-class dataset.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/train_svmsgd.cpp>
%

function varargout = train_svmsgd_demo_gui()
    % data
    sz = [594 841];  % [height width]
    samples = [];    % Set of train samples. Contains points on image
    responses = [];  % Set of responses for train samples

    % create the UI
    h = buildGUI();
    onHelp([],[]);   % display instructions
    if nargout > 0, varargout{1} = h; end

    function h = buildGUI()
        %BUILDGUI  Creates the UI

        h = struct();
        h.fig = figure('Name','Train SVMSGD', 'Menubar','none', ...
            'Position',[200 200 sz(2) sz(1)]);
        if ~mexopencv.isOctave()
            %HACK: not implemented in Octave
            movegui(h.fig, 'center');
        end
        h.ax = axes('Parent',h.fig, 'Units','normalized', 'Position',[0 0 1 1]);
        img = zeros([sz 3], 'uint8');
        if ~mexopencv.isOctave()
            h.img = imshow(img, 'Parent',h.ax);
        else
            %HACK: https://savannah.gnu.org/bugs/index.php?45473
            axes(h.ax);
            h.img = imshow(img);
        end

        % register mouse button handlers and change cursor
        set(h.fig, 'Pointer','cross', 'WindowKeyPressFcn',@onType, ...
            'WindowButtonDownFcn',@onMouseDown);
    end

    function onHelp(~,~)
        %ONHELP  Display usage help dialog

        hd = helpdlg({
            'Left-click the mouse to add a positive sample.'
            'Right-click the mouse to add a negative sample.'
            'Hot keys:'
            '  h: usage dialog'
            '  r: reset'
            '  s: save current output image as PNG image'
            '  e: export current data as MAT-file'
            '  q: quit the program'
        }, 'Interactive SVMSGD demo');
        set(hd, 'WindowStyle','modal');
    end

    function onReset(~,~)
        %ONRESET  Restart from scratch

        % reset data
        samples = [];
        responses = [];

        % update plot
        img = zeros([sz 3], 'uint8');
        set(h.img, 'CData',img);
        drawnow;
    end

    function onType(~,e)
        %ONTYPE  Event handler for key press on figure

        % handle keys
        switch e.Key
            case 'r'
                onReset([],[]);
            case 'h'
                onHelp([],[]);
            case 's'
                fname = [tempname() '.png'];
                imwrite(img, fname);
                fprintf('Output saved as "%s"\n', fname);
            case 'e'
                uisave({'samples', 'responses'}, 'data.mat');
            case {'q', 'escape'}
                close(h.fig);
        end
    end

    function onMouseDown(~,~)
        %ONMOUSEDOWN  Event handler for mouse down on figure

        % get current location of mouse pointer
        p = get(h.ax, 'CurrentPoint');
        p = round(p(1,1:2));

        % add point to train set with corresponding positive/negative class
        samples(end+1,:) = p;
        if strcmp(get(h.fig,'SelectionType'), 'normal')
            responses(end+1) = +1;
        else
            responses(end+1) = -1;
        end

        % process (train model and draw results on image)
        [weights, shift] = doTrain(samples, responses(:));
        pts = doFindPointsForLine(sz, weights, shift);
        img = doRedraw(sz, samples, responses, pts);

        % update plot
        set(h.img, 'CData',img);
        drawnow;
    end
end

function [weights, shift] = doTrain(samples, responses)
    %DOTRAIN  Train with SVMSGD algorithm
    %
    %     [weights, shift] = doTrain(samples, responses)
    %
    % ## Input
    % * __samples__, __responses__ train set
    %
    % ## Output
    % * __weights__, __shift__ vector of decision function of SVMSGD algorithm
    %

    weights = [];
    shift = [];
    if numel(unique(responses)) < 2
        % ensure we have at least one point from each class
        return;
    end

    model = cv.SVMSGD();
    model.train(samples, responses);

    if model.isTrained()
        weights = model.getWeights();
        shift = model.getShift();

        display(model)
        fprintf('%f*x + %f*y + %f = 0\n', weights(1), weights(2), shift);
    end
end

function pts = doFindPointsForLine(sz, weights, shift)
    %DOFINDPOINTSFORLINE  Find two points for drawing decision function line (w*x = 0)

    pts = [];
    if isempty(weights)
        return;
    end

    % axis-aligned border segments
    segments = {
        [sz(2) 0; sz(2) sz(1)];  % right
        [0 sz(1); sz(2) sz(1)];  % top
        [0 0; sz(2) 0];          % bottom
        [0 0; 0 sz(1)]           % left
    };

    % test intersection against each segment until we collect two points
    for i=1:numel(segments)
        pt = doFindCrossPointWithBorders(weights, shift, segments{i});
        if ~isempty(pt)
            pts(end+1,:) = pt;
            if size(pts,1) >= 2
                return;
            end
        end
    end
end

function pt = doFindCrossPointWithBorders(weights, shift, seg)
    %DOFINDCROSSPOINTWITHBORDERS  Find intersection of line (w*x = 0) and segment
    %
    % (y = HEIGHT, 0 <= x <= WIDTH) or (x = WIDTH, 0 <= y <= HEIGHT)
    %
    % decision function line equation: w(1)*x + w(2)*y + s = 0
    % border equations either x=c or y=c
    %

    xmn = min(seg(:,1));
    xmx = max(seg(:,1));
    ymn = min(seg(:,2));
    ymx = max(seg(:,2));
    pt = [];
    if xmn == xmx && weights(2) ~= 0
        % intersect with vertical border
        x = xmn;
        y = floor(-(weights(1) * x + shift) / weights(2));
        if ymn <= y && y <= ymx
            pt = [x y];
        end
    elseif ymn == ymx && weights(1) ~= 0
        % intersect with horizontal border
        y = ymn;
        x = floor(-(weights(2) * y + shift) / weights(1));
        if xmn <= x && x <= xmx
            pt = [x y];
        end
    end
end

function img = doRedraw(sz, samples, responses, pts)
    %DOREDRAW  Redraw point set and decision function line (w*x = 0)

    img = zeros([sz 3], 'uint8');

    if ~isempty(samples)
        img = cv.circle(img, samples(responses==+1,:), 6, ...
            'Color',[255 0 0], 'Thickness','Filled');
        img = cv.circle(img, samples(responses==-1,:), 6, ...
            'Color',[0 0 255], 'Thickness','Filled');
    end

    if ~isempty(pts)
        img = cv.line(img, pts(1,:), pts(2,:), 'Color',[0 255 0]);
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
classdef PoleFigure < dynProp & dynOption
%
% The class *PoleFigure* is used to store experimetnal pole figure
% intensitied, i.e., XRD, synchrotron or neuron data. It provides several
% <PoleFigureCorrection.html data correction methods> as well as the
% <PoleFigure2ODF.html reconstruction of an orientation density function
% (ODF)>. Importing pole figure data is explained in <PoleFigureImport.html
% this section>.
%
% Input
%  h     - crystal directions (@vector3d | @Miller)
%  r     - specimen directions (@S2Grid)
%  intensities - diffraction counts (double)
%  CS,SS - crystal, specimen @symmetry
%
% Options
%  superposition - weights for superposed crystal directions
%  background    - background intensities
%
% Class Properties
%  allH - cell of @Miller
%  allR - cell of @vector3d
%  allI - cell of diffraction intensities
%  c    - structure coefficients
%  SS   - specimen symmetry
%
% Dependent Class Properties
%  CS      - @crystalSymmetry
%  h       - @Miller direction of single pole figure
%  r       - specimen directions
%  intensities - diffraction intensities
%  antipodal -
%
% See also
% ImportPoleFigureData loadPoleFigure loadPoleFigure_generic
% This section describes the class *PoleFigure* and gives an overview of
% the functionality MTEX offers to analyze pole figure data.

  properties
    allH = {}           % crystal directions
    allR = {}           % specimen directions
    allI = {}           % intensities
    c = {}              % structure coefficients for superposed pole figures
    SS = specimenSymmetry       % specimen symmetry
  end

  properties (Dependent = true)
    CS                  % crystal symmetry
    h                   % crystal direction of single pole figure
    r                   % specimen directions
    intensities         % diffraction intensities
    antipodal
  end

  methods

    function pf = PoleFigure(h,r,intensities,varargin)
      % constructor

      if nargin == 0, return;end

      pf.allH = ensurecell(h);
      pf.allR = ensurecell(r);
      if numel(pf.allR) == 1, pf.allR = repmat(pf.allR,size(pf.allH));end
      if ~check_option(varargin,'complete'), pf.allR{1}.antipodal = true;end
      pf.allI = ensurecell(intensities);

      pf.c = ensurecell(get_option(varargin,'superposition',...
        cellfun(@(x) ones(1,length(x)),pf.allH,'uniformoutput',false)));

      % normalize structure coefficients
      %pf.c = cellfun(@(x) x./sum(x),pf.c,'uniformOutput',false);

      % extract symmetries
      pf.CS = getClass(varargin,'crystalSymmetry',pf.CS);
      pf.SS = getClass(varargin,'specimenSymmetry',pf.SS);

    end

    function pf = set.CS(pf,CS)

      for i = 1:length(pf.allH)
        pf.allH{i}.CS = CS;
      end
    end

    function CS = get.CS(pf)

      CS = pf.allH{1}.CS;

    end

    function h = get.h(pf)
      h = pf.allH;
      h = horzcat(h{:});
    end

    function r = get.r(pf)
      try
        r = [pf.allR{:}];
      catch
        for i = 1:numel(pf.allR)
          pf.allR{i} = pf.allR{i}(:);
        end
        r = vertcat(pf.allR{:});
      end
    end

    function i = get.intensities(pf)
      try
        i = [pf.allI{:}];
      catch
        for i = 1:numel(pf.allI)
          pf.allI{i} = pf.allI{i}(:);
        end
        i = vertcat(pf.allI{:});
      end
    end

    function pf = set.intensities(pf,i)

      if numel(i) == 1
        for ipf = 1:numel(pf.allI)

          pf.allI{ipf} = i*ones(size(pf.allI{ipf}));

        end
      else
        cs = cumsum([0,cellfun('prodofsize',pf.allI)]);

        for ipf = 1:numel(pf.allI)
          pf.allI{ipf} = reshape(i(cs(ipf)+1:cs(ipf+1)),size(pf.allI{ipf}));
        end
      end
    end

    function out = get.antipodal(pf)
      out = pf.allR{1}.antipodal;
    end

    function pf = set.antipodal(pf,value)
      for i = 1:pf.numPF
        pf.allR{i}.antipodal = value;
      end
    end

    function varargout = size(pf,varargin)
      [varargout{1:nargout}] = size(pf.r,varargin{:});
    end

    function varargout = length(pf,ip)
      if nargin == 2
        if isempty(ip)
          varargout{1} = cellfun(@length,pf.allR);
        else
          [varargout{1:nargout}] = length(pf.allR{ip});
        end
      else
        [varargout{1:nargout}] = length(pf.r);
      end
    end

    function n = numPF(pf)
      n = numel(pf.allH);
    end

    function e = end(pf,i,n)
      % overloaded end function

      if n==1
        e = numel(pf.r);
      else
        e = size(pf.r,i);
      end
    end

  end

  methods (Static = true)
    [pf,interface,options] = load(fname,varargin)
  end

end

author: mtex-toolbox

repo: mtex

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


text:
function [X_poly] = polyFeatures(X, p)
%POLYFEATURES Maps X (1D vector) into the p-th power
%   [X_poly] = POLYFEATURES(X, p) takes a data matrix X (size m x 1) and
%   maps each example into its polynomial features where
%   X_poly(i, :) = [X(i) X(i).^2 X(i).^3 ...  X(i).^p];
%

% You need to return the following variables correctly.
X_poly = zeros(numel(X), p);

% ====================== YOUR CODE HERE ======================
% Instructions: Given a vector X, return a matrix X_poly where the p-th
%               column of X contains the values of X to the p-th power.
%
%

% Think of some vectorized implementation
for i = 1:numel(X)
    for j = 1:p
        X_poly(i, j) = X(i)^j;
    end
end

% =========================================================================

end

author: SaveTheRbtz

repo: ml-class


text:
function c = band(a,p,q)
%BAND         Extract band from matrix a, lower bandwidth p, upper bandwidth q
%   if parameter q is omitted, q:=p
%
%   c = band(a,p,q)
%

% written  10/16/98     S.M. Rump
% modified 11/30/98     S.M. Rump  improved speed
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if nargin<3
    q = p;
  end

  c = a;

  if a.complex
    c.mid = tril(triu(a.mid,-p),q);
    c.rad = tril(triu(a.rad,-p),q);
  else
    c.inf = tril(triu(a.inf,-p),q);
    c.sup = tril(triu(a.sup,-p),q);
  end

repo: OpenMAS


text:
%% hip_implant
% Below is a demonstration of the features of the |hip_implant| function

%%
clear; close all; clc;

%% Syntax
% |[F,V]=hip_implant;|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<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 y = logical(x)
% Allows char input
if ischar(x)
    x = double(x);
end
y = builtin('logical', x);
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL