text:
% plotfresp() - Plot a filter's impulse, step, magnitude, and phase response
%
% Usage:
%   >> plotfresp(b, a, nfft, fs, dir);
%
% Inputs:
%   b     - vector numerator coefficients
%
% Optional inputs:
%   a     - scalar or vector denominator coefficients (IIR support is
%           experimental!) {default 1}
%   nfft  - scalar number of points {default 512}
%   fs    - scalar sampling frequency {default 1}
%   dir   - string filter direction {default 'onepass'}
%
% Author: Andreas Widmann, University of Leipzig, 2005
%
% See also:
%   firws

%123456789012345678901234567890123456789012345678901234567890123456789012

% Copyright (C) 2005-2014 Andreas Widmann, University of Leipzig, widmann@uni-leipzig.de
%
% 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
%
% $Id$

function plotfresp(b, a, nfft, fs, dir)

if nargin < 5 || isempty(dir)
    dir = 'onepass';
end
if nargin < 4  || isempty(fs)
    fs = 1;
end
if nargin < 3 || isempty(nfft)
    nfft = 512;
end
if nargin < 2 || isempty(a)
    a = 1;
end
if nargin < 1
    ft_error('Not enough input arguments.');
end

% FIR?
if isscalar(a) && a == 1
    isFIR = true;
else
    isFIR = false;
end

% Linear phase FIR
if isFIR && all(b(:)' == fliplr(b(:)')) % TODO: antisymmetric
    isLinPhaseFir = true;
else
    isLinPhaseFir = false;
end

% Twopass/zerophase?
if strncmp('twopass', dir, 7)
    isTwopass = true;
    isZerophase = true;
elseif strcmp('onepass-zerophase', dir) || strcmp('onepass-reverse-zerophase', dir)
    if ~isLinPhaseFir
        ft_error('Onepass-zerophase filtering is only allowed for linear-phase FIR filters.')
    end
    isTwopass = false;
    isZerophase = true;
else
    isTwopass = false;
    isZerophase = false;
end

% Impulse response
if isFIR
    impresp = b(:)';
else
    if ~exist('impz', 'file')
        ft_warning('Plotting IIR filter responses requires signal processing toolbox.')
        return
    end
    impresp = impz(b, a)';
end

% Twopass
if isTwopass
    impresp = conv(impresp, fliplr(impresp));
end
n = length(impresp);

% Zerophase
if isZerophase
    groupdelay = (n - 1) / 2;
    x = -groupdelay:groupdelay;
else
    x = 0:n - 1;
end

nfft = max([2^ceil(log2(n)) nfft]); % Do not truncate impulse response
f = linspace(0, fs / 2, nfft / 2 + 1);
z = fft(impresp, nfft);
z = z(1:nfft / 2 + 1);

% Find open figure window
H = findobj('Tag', 'plotfiltresp', 'type', 'figure');
if ~isempty(H)
    figure(H);
else
    H = figure;
    set(H, 'Tag', 'plotfiltresp');
    posArray = get(H, 'Position');
    posArray(3) = posArray(4) * 1.6;
    set(H, 'Position', posArray);
end

% Formatting
titlePropArray = {'Fontweight', 'bold'};
axisPropArray = {'NextPlot', 'add', 'XGrid', 'on', 'YGrid', 'on', 'Box', 'on'};

% Impulse resonse
ax(1) = subplot(2, 3, 1, axisPropArray{:});
stem(x, impresp, 'fill')
title('Impulse response', titlePropArray{:});
ylabel('Amplitude');

% Step response
ax(4) = subplot(2, 3, 4, axisPropArray{:});
stem(x, cumsum(impresp), 'fill');
title('Step response', titlePropArray{:});
ylimArray = ylim;
if ylimArray(2) < -ylimArray(1) + 1;
    ylimArray(2) = -ylimArray(1) + 1;
    ylim(ylimArray);
end
xMin = []; xMax = [];
childrenArray = get(ax(4), 'Children');
for iChild =1:length(childrenArray)
    xData = get(childrenArray(iChild), 'XData');
    xMin = min([xMin min(xData)]);
    xMax = max([xMax max(xData)]);
end
set(ax([1 4]), 'XLim', [xMin xMax]);
ylabel('Amplitude');

% Magnitude response
ax(2) = subplot(2, 3, 2, axisPropArray{:});
plot(f, abs(z));
title('Magnitude response', titlePropArray{:});
ylabel('Magnitude (linear)');

ax(5) = subplot(2, 3, 5, axisPropArray{:});
plot(f, 20 * log10(abs(z)));
title('Magnitude response', titlePropArray{:});
ylimArray = ylim;
if ylimArray(1) < -200
    ylimArray(1) = -200;
    ylim(ylimArray);
end
ylabel('Magnitude (dB)');

% Phase response
ax(3) = subplot(2, 3, 3, axisPropArray{:});
phaseresp = unwrap(angle(z));
if isZerophase % Correct delay for zero-phase FIR filter?
    delay = -f / fs * groupdelay * 2 * pi;
    phaseresp = phaseresp - delay;
    phaseresp = mod(round(phaseresp / pi), 2) * pi; % Avoid rounding errors; linear-phase FIR only!
end
plot(f, phaseresp);
title('Phase response', titlePropArray{:});
ylabel('Phase (rad)');

% Formatting
xlabelArray = get(ax(1:5), 'XLabel');
if fs == 1
    set([xlabelArray{[2 3 5]}], 'String', 'Normalized frequency (2 \pi rad / sample)');
else
    set([xlabelArray{[2 3 5]}], 'String', 'Frequency (Hz)');
end
set([xlabelArray{[1 4]}], 'String', 'n (samples)');
set(ax([2 3 5]), 'XLim', [0 fs / 2]);
set(ax(1:5), 'ColorOrder', circshift(get(ax(1), 'ColorOrder'), -1));

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [newMap] = changeMetColor(map, metList, color)
% Change color of every metabolite from a list of Names
%
% USAGE:
%
%   [newMap] = changeMetColor(map, metList, color);
%
% INPUTS:
%   map:            file from CellDesigner parsed to MATLAB format
%   metList:        List of metabolites names
%
% OPTIONAL INPUT:
%   color:          New color of metabolites from list(default: RED)
%
% OUTPUT:
%   newMap:         Matlab structure of map with reaction modifications
%
% .. Authors:
%       - A.Danielsdottir 17/07/2017 LCSB. Belval. Luxembourg
%       - N.Sompairac - Institut Curie, Paris, 17/07/2017.

    if nargin < 3
        color = 'RED';
    end

    newMap = map;

    colors = createColorsMap();
    % Index for specName is the same as for corresponding specID
    specID = newMap.specID(ismember(newMap.specName, metList));
    index = find(ismember(newMap.molID, specID));

    % Change color
    for i = index'
        newMap.molColor{i} = colors(color);
    end

end

author: opencobra

repo: cobratoolbox


text:
clear
% define the root name of database
root = '../data_preparation/prepared_data/';

% which scales we're doing
sigma = 1;
num_samples = 5e5;

scales = [0.25,0.35,0.5];
frontalView = 1;

profileViewInds = [];

version = 'cofw';
ratio_neg = 5;
norm = 1;

data_loc = 'cofw_';
rng(0);

similarities = {};
sparsity = 1;
sparsity_types = [];

lambda_a = 100;
lambda_b = 1000;
lambda_th = 1;
num_layers = 7;

for s=scales

    Train_all_patch_experts(root, frontalView, profileViewInds,...
        s, sigma, version, 'ratio_neg', ratio_neg,...
        'num_samples', num_samples, 'data_loc', data_loc,...
        'normalisation_size', 19, 'similarity_types', similarities, 'sparsity', sparsity,...
        'sparsity_types', sparsity_types, 'lambda_a', lambda_a, 'lambda_b', lambda_b, 'lambda_th', lambda_th, 'num_layers', num_layers);
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function [X_den,P]=denoise_TV_One(Xobs,lambda,l,u,P_init,pars)
%This function implements the FISTA method for TV denoising problems.
%
% INPUT
% Xobs ..............................an observed noisy image.
% lambda ........................ parameter
% pars.................................parameters structure
% pars.MAXITER ..................... maximum number of iterations
%                                                      (Default=100)
% pars.epsilon ..................... tolerance for relative error used in
%                                                       the stopping criteria (Default=1e-4)
% pars.print ..........................  1 if a report on the iterations is
%                                                       given, 0 if the  report is silenced
% pars.tv .................................. type of total variation
%                                                      penatly.  'iso' for isotropic (default)
%                                                      and 'l1' for nonisotropic
%
% OUTPUT
% X_den ........................... The solution of the problem
%                                            min{||X-Xobs||^2+2*lambda*TV(X)}

%%% Written by Chen Chen and Junzhou Huang at UT Arlington
%%% April. 13, 2012

%%% Related paper: Chen Chen and Junzhou Huang, "Compressive Sensing MRI
%%% with Wavelet Tree Sparsity"

%Define the Projection onto the box
if((l==-Inf)&&(u==Inf))
    project=@(x)x;
elseif (isfinite(l)&&(u==Inf))
    project=@(x)(((l<x).*x)+(l*(x<=l)));
elseif (isfinite(u)&&(l==-Inf))
     project=@(x)(((x<u).*x)+((x>=u)*u));
elseif ((isfinite(u)&&isfinite(l))&&(l<u))
    project=@(x)(((l<x)&(x<u)).*x)+((x>=u)*u)+(l*(x<=l));
else
    error('lower and upper bound l,u should satisfy l<u');
end

% Assigning parameres according to pars and/or default values
flag=exist('pars', 'var');
if (flag&&isfield(pars,'MAXITER'))
    MAXITER=pars.MAXITER;
else
    MAXITER=100;
end
% if (flag&&isfield(pars,'epsilon'))
%     epsilon=pars.epsilon;
% else
%     epsilon=1e-4;
% end
% if(flag&&isfield(pars,'print'))
%     prnt=pars.print;
% else
%     prnt=1;
% end
if(flag&&isfield(pars,'tv'))
    tv=pars.tv;
else
    tv='iso';
end

[m,n]=size(Xobs);
% clear P; clear R;
if(isempty(P_init))
    P{1}=zeros(m-1,n);    P{2}=zeros(m,n-1);
    R{1}=zeros(m-1,n);    R{2}=zeros(m,n-1);
else
    P{1}=P_init{1};    P{2}=P_init{2};
    R{1}=P_init{1};    R{2}=P_init{2};
end
tk=1;tkp1=1;count=0;i=0;

D=zeros(m,n);%fval=inf;fun_all=[];
while((i<MAXITER)&&(count<5))
%    fold=fval;
    i=i+1;
%     Dold=D;
    Pold=P;
    tk=tkp1;
    D=project(Xobs-lambda*Lforward(R, m, n));
    Q=Ltrans(D, m, n);
    %%%%%%%%%%
    % Taking a step towards minus of the gradient
    P{1}=R{1}+1/(8*lambda)*Q{1};
    P{2}=R{2}+1/(8*lambda)*Q{2};

    %%%%%%%%%%
    % Peforming the projection step
    switch tv
        case 'iso'
            A=[P{1}.^2;zeros(1,n)]+[P{2}.^2,zeros(m,1)];
            A=sqrt(max(A,1));
            P{1}=P{1}./A(1:m-1,:); P{2}=P{2}./A(:,1:n-1);
        case 'l1'
            P{1}=P{1}./(max(abs(P{1}),1));
            P{2}=P{2}./(max(abs(P{2}),1));
        otherwise
            error('unknown type of total variation. should be iso or l1');
    end

    %%%%%%%%%%
    %Updating R and t
    tkp1=(1+sqrt(1+4*tk^2))/2;

    R{1}=P{1}+(tk-1)/(tkp1)*(P{1}-Pold{1});
    R{2}=P{2}+(tk-1)/tkp1*(P{2}-Pold{2});

%     re=norm(D-Dold,'fro')/norm(D,'fro');
%     if (re<epsilon)
%         count=count+1;
%     else
%         count=0;
%     end
    C=Xobs-lambda*Lforward(P, m, n);
    PC=project(C);
%     fval=-norm(C-PC,'fro')^2+norm(C,'fro')^2;
%     fun_all=[fun_all;fval];
end
X_den=D;iter=i;

function X=Lforward(P, m, n)

%       [m2,n2]=size(P{1});
%       [m1,n1]=size(P{2});
%
%       if (n2~=n1+1)
%           error('dimensions are not consistent')
%       end
%       if(m1~=m2+1)
%           error('dimensions are not consistent')
%       end
%
%       m=m2+1;
%       n=n2;

      X=zeros(m,n);
      X(1:m-1,:)=P{1};
      X(:,1:n-1)=X(:,1:n-1)+P{2};
      X(2:m,:)=X(2:m,:)-P{1};
      X(:,2:n)=X(:,2:n)-P{2};
   end

   function P=Ltrans(X, m, n)
%       [m,n]=size(X);
      P{1}=X(1:m-1,:)-X(2:m,:);
      P{2}=X(:,1:n-1)-X(:,2:n);
   end
end

author: thomaskuestner


text:
function y = gt(X,Y)
%GT (overloaded)

if isa(X,'blkvar')
    X = sdpvar(X);
end

if isa(Y,'blkvar')
    Y = sdpvar(Y);
end

warning('YALMIP:strict','Strict inequalities are not supported. A non-strict has been added instead'')');

try
    y = constraint(X,'>=',Y);
catch
    error(lasterr)
end

author: shamilmamedov


text:
function l = IsLastBefore(timestamps,list,varargin)

%IsLastBefore - Identify last item before each of a list of timestamps.
%
% For each element i of a list of test timestamps, find in a list of reference timestamps
% the last element j smaller than i. Returns a list of logical indices (the elements of the
% test list which are smaller than all the elements of the reference timestamps are ignored).
%
%  USAGE
%
%    logical = IsLastBefore(timestamps,list,<options>)
%
%    timestamps     a list of reference timestamps
%    list           a list of test timestamps
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'strict'      'on' for strict (<) comparisons, 'off' otherwise (<=)
%                   (default = 'off')
%    =========================================================================
%
%  SEE
%
%    See also IsFirstAfter.

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% 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.

strict = 'off';

if nargin < 2 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help IsLastBefore'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
  if ~ischar(varargin{i}),
    error(['Parameter ' num2str(i+1) ' is not a property (type ''help IsLastBefore'' for details).']);
  end
  switch(lower(varargin{i})),
    case 'strict',
		strict = lower(varargin{i+1});
		if ~isstring_FMAT(strict,'on','off'),
			error('Incorrect value for property ''strict'' (type ''help IsLastBefore'' for details).');
		end
    otherwise,
      error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help IsLastBefore'' for details).']);
  end
end

% Make sure 'timestamps' looks like [t1;t2;...;tn] and 'list' like [l1 l2 ... lm]
if min(size(timestamps)) ~= 1 | min(size(list)) ~= 1,
	error('Incorrect sizes (must be vectors)');
end
if size(timestamps,1) == 1,
	timestamps = timestamps';
end
if size(list,2) == 1,
	list = list';
end

% Now compare [t1 t1 ... t1;t2 t2 ... t2;...;tn tn ... tn] with [l1 l2 ... lm;l1 l2 ... lm;...;l1 l2 ... lm]
if strcmp(strict,'on'),
	where = repmat(timestamps,1,length(list)) < repmat(list,length(timestamps),1);
else
	where = repmat(timestamps,1,length(list)) <= repmat(list,length(timestamps),1);
end
% This yielded something like [1 1 1;1 1 1;...;0 1 1;...;0 0 1;...;0 0 0;...;0 0 0]
% where the transition from 1 to 0 in each column indicates the index of the last element
% in 'timestamps' inferior to the corresponding element in 'list'. Now, find these transitions.
% (we append ones before and zeros after 'where' to handle the special cases where the first
% element in 'timestamps' inferior to the element of 'list' is the last element of 'timestamps'
% or does not exist).
where = diff([ones(1,size(where,2));where;zeros(1,size(where,2))]);
% Convert the resulting logical matrix into a list of logical indices
l = sum(where,2) == -1;
l = l(2:length(timestamps)+1);
l = logical(l);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function net = fcnInitializeModel8sGeneric(labelCount, net)
%FCNINITIALIZEMODEL8S Initialize the FCN-16S model from FCN-32

%% Remove the last layer
net.removeLayer('deconv16') ;

%% Add the first deconv layer
filters = single(bilinear_u(4, 1, labelCount)) ;
net.addLayer('deconv2bis', ...
  dagnn.ConvTranspose('size', size(filters), ...
                      'upsample', 2, ...
                      'crop', 1, ...
                      'hasBias', false), ...
             'x41', 'x42', 'deconv2bisf') ;
f = net.getParamIndex('deconv2bisf') ;
net.params(f).value = filters ;
net.params(f).learningRate = 0 ;
net.params(f).weightDecay = 1 ;

%% Add a convolutional layer that take as input the pool3 layer
net.addLayer('skip3', ...
     dagnn.Conv('size', [1 1 256 labelCount]), ...
     'x17', 'x43', {'skip3f','skip3b'});

f = net.getParamIndex('skip3f') ;
net.params(f).value = zeros(1, 1, 256, labelCount, 'single') ;
net.params(f).learningRate = 0.01 ;
net.params(f).weightDecay = 1 ;

f = net.getParamIndex('skip3b') ;
net.params(f).value = zeros(1, 1, labelCount, 'single') ;
net.params(f).learningRate = 2 ;
net.params(f).weightDecay = 1 ;

%% Add the sumwise layer
net.addLayer('sum2', dagnn.Sum(), {'x42', 'x43'}, 'x44') ;

%% Add deconvolutional layer implementing bilinear interpolation
filters = single(bilinear_u(16, labelCount, labelCount)) ;
net.addLayer('deconv8', ...
  dagnn.ConvTranspose('size', size(filters), ...
                      'upsample', 8, ...
                      'crop', 4, ...
                      'numGroups', labelCount, ...
                      'hasBias', false), ...
             'x44', 'prediction', 'deconv8f') ;

f = net.getParamIndex('deconv8f') ;
net.params(f).value = filters ;
net.params(f).learningRate = 0 ;
net.params(f).weightDecay = 1 ;

% Make the output of the bilinear interpolator is not discared for
% visualization purposes
net.vars(net.getVarIndex('prediction')).precious = 1 ;

% empirical test
if 0
  figure(100) ; clf ;
  n = numel(net.vars) ;
  for i=1:n
    vl_tightsubplot(n,i) ;
    showRF(net, 'input', net.vars(i).name) ;
    title(sprintf('%s', net.vars(i).name)) ;
    drawnow ;
  end
end

author: nightrome

repo: matconvnet-calvin


text:
function k=FibonacciNumInv(n,choice)
%FIBONACCINUMINV Given a real number >=1, determine the k such that
%                n==FibonacciNum(k). If n given is not a Fibonacci number,
%                then the behaviour is determined by the choice option.
%
%INPUTS: n A scalar or matrix of values for which the argument of
%          FibonacciNum is desired, n>=1.
%   choice An optional parameter that determines what is returned if n
%          is not an exact Fibonacci number. If this parameter is omitted,
%          then the default is one (the next lower value).
%          0 means return the closest value.
%          1 means return the next lower value.
%          2 means return the next higher value.
%
%OUTPUTS: k The argument of FibonacciNum that gives n, or an appropriately
%           close value selected by the choice option. k will always be
%           >=2.
%
%A non-recursive expression for Fibonacci numbers is given in [1] and can
%be inverted. However, due to finite prevision issues, for n in the mid
%70's the value returned can be off. Thus, testing needs to be done for all
%inputs in the choice option to make sure that the proper value is chosen.
%
%REFERENCES:
%[1] Chandra, Pravin and Weisstein, Eric W. "Fibonacci Number." From
%    MathWorld--A Wolfram Web Resource.
%    http://mathworld.wolfram.com/FibonacciNumber.html
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(nargin<2||isempty(choice))
        choice=1;
    end

    %The golden ratio.
    phi=(1+sqrt(5))/2;
    k=fix(log(n*sqrt(5))/log(phi)+(1/2));

    %Due to finite precision errors, even if n is a Fibonacci number, k
    %might be the correct value, it might be 1 less or it might be 1 too
    %much. We will adjust k to deal with the situation. Note that from
    %above, k>=2, so there are no issues with k-1 being <1.
    switch(choice)
        case 0%Return the closest value
            valsLower=FibonacciNum(k-1);
            valsCur=FibonacciNum(k);
            valsUpper=FibonacciNum(k+1);

            diffLower=abs(n-valsLower);
            diffCur=abs(n-valsCur);
            diffUpper=abs(n-valsUpper);

            sel=(diffLower<diffCur)&(diffLower<diffUpper);
            k(sel)=k(sel)-1;
            sel=(valsUpper<diffCur)&(diffUpper<diffLower);
            k(sel)=k(sel)+1;
        case 1%Return the next lower value.
            %We have to do a check, because finite precision error
            %above might make k be too large.
            valsCur=FibonacciNum(k);

            sel=valsCur>n;
            k(sel)=k(sel)-1;
        case 2%Return the next higher value.
            valsCur=FibonacciNum(k);
            sel=valsCur<n;
            k(sel)=k(sel)+1;
        otherwise
            error('Unknown choice value specified.')
    end
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 fixedStatus=fixFlatMesh1File(inFileName,outFileName,grayFile)
% function fixedStatus=fixFlatMesh1File(inFileName,outFileName)
% mrFlatMesh v1.0 had a nasty bug : the 3dGlocs of layer 1 gnodes were
% wrong by -1 voxel in all directions (i.e. a node at [1,2,3] was
% given a 3D location of [0,1,2]
% This has been fixed in v1.5 upwards but old flat.mat files need to be fixed.
% This routine reads in a flat.mat file and tries to fix it by adding 1 to all the layer 1 gLocs3ds
% This is not simple because the layer information has been lost by this point. However,
% infoStr.grayFile should contain the name of the gray file that was used to
% make the flat.mat file. This gray file can be used to identify l1 gNodes.
% If it doesn't exist, the user is prompted to browse for a gray file.
% You can enter this routine with no parameters in which case you'll be asked
% to browse for the input and output files as well
%
% Output files contain an extra field in the infoStr: gloc3d_fixed (==indices of fixed nodes if the fix was successful, 0 otherwise)

if (~exist('inFileName','var'))
    [inFileName,inPathName]=uigetfile('','Locate a flat file');
    inFileName=[inPathName,inFileName];
end

flatFile=load(inFileName);

if (~isfield(flatFile,'gLocs3d'))
    error('No gLocs3d in this file. Can''t go on');

end

gLocs3d=flatFile.gLocs3d;

if (~isfield(flatFile,'infoStr'))
    disp('Warning - no infoStr field found.');
    thisInfo='';
else
    thisInfo=flatFile.infoStr;
end

% Check version information on this file to see if it needs to be fixed.
if (isfield(thisInfo,'Version'))
    if (thisInfo.Version>1.5)
        disp('This file is generated by mrFlatMesh > 1.5. It does not need to be fixed.');
        return;
    end
else
    disp('No version info in file. It may need fixing');
end

% Okay, try to load in the gray file

grayFileName=thisInfo.grayFile;
disp('Attempting to load in the gray file that this flat.mat was made from..');
[nodes, edges, vSize] = readGrayGraph(grayFileName);
if (~nodes & ~edges & ~vSize) % This happens when the file could not be read
    disp(['Failed to read the gray file: ',grayFileName]);
    disp('This may be because you are running on a different platform. Or the gray file may have moved.');
    disp('Try to locate it now....');
    [grayFileName,grayPathName]=uigetfile('','Locate a gray file');
    grayFileName=[grayPathName,grayFileName];
    if (~grayFileName)
        disp('Exiting');
        return;
    else % if a gray file was selected
        % try to load this in..
        [nodes, edges, vSize] = readGrayGraph(grayFileName);
        if (~nodes & ~edges & ~vSize) % This happens when the file could not be read
            error ('Could not read this gray file');
        end % endif bad gray file read
    end % endif gray file selected with uigetfile
end % endif couldn't load the gray file from the infoStr

% We now have nodes, edges, vSize
% nodes(6,:) give the layer...
l1NodeIndices=find(nodes(6,:)==1);
l1NodeCoords=nodes(1:3,l1NodeIndices);
% We now have the coordinates of the l1 nodes in the full gray. We need to know
% which of these are also in gLocs3d

% Note that this isn't as simple as it sounds: we can subtract [1 1 1] from nodes and then
% do an intersection with gLocs3d but we're also going to catch some l2,3,4 nodes as well
% The only thing  that saves us is that we know that gLocs3d was created by concatenating
% l1,l2,l3,l4 so we know that all the l1 nodes should be in a contiguous block at the
% start. We might get a few errors at the end of the block though.
%
% First add 1 to all of gLocs3d
gLocsInc3d=gLocs3d+1;

% Now use sub2ind to prepare for the intersect
maxDim=max([gLocsInc3d(:);l1NodeCoords(:)]);
gLocsIncInd=sub2ind([maxDim,maxDim,maxDim],gLocsInc3d(:,1),gLocsInc3d(:,2),gLocsInc3d(:,3));
l1NodeInd=sub2ind([maxDim,maxDim,maxDim],l1NodeCoords(1,:),l1NodeCoords(2,:),l1NodeCoords(3,:));

% Find candidate l1nodes
[candidateL1Nodes,indAB,indBA]=intersect(gLocsIncInd,l1NodeInd);

% This returns the nodes that have the same locations in the gray file and the flat gLocs3d
% indAB is an index vector such that candidateL1Nodes=gLocsIncInd(indAB)
% indBA is an index such that candiadateL1Nodes=l1NodeInd(indBA)
% So: indAB should contain a contiguous run from 1:(numL1GNodes) where numL1GNodes is the number of l1GNodes
% it may also contain extra entries where it intersected with other random l2,l3,l4 nodes

% We only want the initial, contiguous indices.
indAB=indAB(:); % Make sure
indABPlus1=shift(indAB,[1,0]); % Circular shift the entries up by 1
diffAB=indAB-indABPlus1; % This will be 1 for the first set of contiguous entries
endofBlock=find(diffAB~=1);
endofBLock=endofBlock(1);

% We're almost done :)
gLocs3d(1:endOfBlock,:)=gLocs3d(1:endOfBlock)+1;
% TOTO Mark field to say that we've fixed it.
% Save out the new flat file

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [g2vMap, distSq] = mrmMapGrayToVertices(grayNodes, vertexCoords, mmPerVox, distThresh)
%
% g2vMap = mrmMapGrayToVertices(grayNodes,vertexCoords, mmPerVox)
%
% Finds a map between all gray nodes and the mesh vertices
% To see the coordinate of a mesh node nearest to a gray node, we
% can use
%
%      initVertices(1:3,g2vMap(idx))
%
%
% HISTORY:
%  2004.03.26 ARW wade@ski.org. Based on RFD's routine mrmVerticesToGray
if notDefined('mmPerVox'), error('Voxel size (mm per vox) is required.'); end;

if notDefined('distThresh')
    if (ispref('VISTA','defaultSurfaceWMMapDist'))
        distThresh = getpref('VISTA','defaultSurfaceWMMapDist');
	else
		if prefsVerboseCheck>=1,
	       disp('No dist threhold preference set. Setting distance threshold to 2 by default');
		end
		distThresh=2;
    end

end

% The gray coordinates are in voxels in the vAnatomy file.  This scales
% them into real physical (mm) coordinates.  And transposes them.
grayCoords = grayNodes([1,2,3], :);
grayCoords = [grayCoords(1,:).*mmPerVox(1); ...
        grayCoords(2,:).*mmPerVox(2); ...
        grayCoords(3,:).*mmPerVox(3) ]';

% Transposes these mesh coordinates, which were already built in real
% physical coordinates.
% Major comments needed here.
vertexCoords = vertexCoords' + 1;

[g2vMap,distSq]=nearpoints(double(grayCoords'),vertexCoords');
badPoints=find(distSq>(distThresh.^2));
g2vMap(badPoints)=0;

g2vMap = int32(g2vMap);

return;


text:
function newParameters = load(filepath)
% load   Load GPT-2
%
%   parameters = load(filepath) will load a GPT-2 model from the directory
%   specified by filepath.

% Load the parameters
s = load(filepath);

% First, assign the hyperparameters and then remove them from the loaded
% structure.
newParameters = struct;
newParameters.Hyperparameters = s.hyperparameters;
s = rmfield(s, 'hyperparameters');

% Next, assign the weights
newParameters.Weights = struct;

originalWeightNames = fieldnames(s);
newWeightNames = erase(originalWeightNames, 'model_');

% Permute all of the weights to match the MATLAB format
for i = 1:numel(originalWeightNames)
    if startsWith(newWeightNames{i},'h')
        nameSplit = split(newWeightNames{i}, '_');
        newParameters.Weights.(char(nameSplit(1))).(char(join(nameSplit(2:end),'_'))) = ...
            s.(originalWeightNames{i});
    else
        newParameters.Weights.(char(newWeightNames{i})) = s.(originalWeightNames{i});
    end
end

newParameters.Weights = iStructRecurseFun(@dlarray, newParameters.Weights);

end

function s = iStructRecurseFun(F, s, varargin)
if ~isstruct(s)
    s = F(s, varargin{:});
else
    fn = fieldnames(s);
    for i = 1:numel(fn)
        others = cellfun(@(o)o.(fn{i}), varargin, 'UniformOutput', false);
        s.(fn{i}) = iStructRecurseFun(F, s.(fn{i}), others{:});
    end
end
end

author: matlab-deep-learning

repo: transformer-models


text:
% Quadratic discrimination (separating ellipsoid)
% Section 8.6.2, Boyd & Vandenberghe "Convex Optimization"
% Original by Lieven Vandenberghe
% Adapted for CVX by Joelle Skaf - 10/16/05
% (a figure is generated)
%
% The goal is to find an ellipsoid that contains all the points
% x_1,...,x_N but none of the points y_1,...,y_M. The equation of the
% ellipsoidal surface is: z'*P*z + q'*z + r =0
% P, q and r can be obtained by solving the SDP feasibility problem:
%           minimize    0
%               s.t.    x_i'*P*x_i + q'*x_i + r >=  1   for i = 1,...,N
%                       y_i'*P*y_i + q'*y_i + r <= -1   for i = 1,...,M
%                       P <= -I

% data generation
n = 2;
rand('state',0);
randn('state',0);
N=50;
X = randn(2,N);  X = X*diag(0.99*rand(1,N)./sqrt(sum(X.^2)));
Y = randn(2,N);  Y = Y*diag((1.02+rand(1,N))./sqrt(sum(Y.^2)));
T = [1 -1; 2 1];  X = T*X;  Y = T*Y;

% Solution via CVX
fprintf(1,'Find the optimal ellipsoid that seperates the 2 classes...');

cvx_begin sdp
    variable P(n,n) symmetric
    variables q(n) r(1)
    P <= -eye(n); %#ok
    sum((X'*P).*X',2) + X'*q + r >= +1; %#ok
    sum((Y'*P).*Y',2) + Y'*q + r <= -1; %#ok
cvx_end

fprintf(1,'Done! \n');

% Displaying results
r = -r; P = -P; q = -q;
Pq = P \ q;
c = 0.25*q'*Pq - r;
xc = -0.5*Pq;
nopts = 1000;
angles = linspace(0,2*pi,nopts);
ell = sqrtm(P/c)\[cos(angles); sin(angles)] + repmat(xc,1,nopts);
graph=plot(X(1,:),X(2,:),'o', Y(1,:), Y(2,:),'o', ell(1,:), ell(2,:),'-');
set(graph(2),'MarkerFaceColor',[0 0.5 0]);
set(gca,'XTick',[]); set(gca,'YTick',[]);
title('Quadratic discrimination');
% print -deps ellips.eps

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function       X = Image2Patch( im_out, par )
% record the non-local patch set and the index of each patch in
% of seed patches in image
im_out         =  single(im_out);
X          =  zeros(par.ps2, par.maxrc, 'double');
k    =  0;
for l = 1:par.ch
    for i = 1:par.ps
        for j = 1:par.ps
            k    =  k+1;
            blk  = im_out(i:end-par.ps+i,j:end-par.ps+j, l);
            X(k,:) = blk(:)';
        end
    end
end

author: csjunxu


text:
%wfBrickIO Wifi interface between MATLAB and the brick
%
% Methods::
%
%  wfBrickIO    Constructor, initialises and opens the wifi connection
%  delete       Destructor, closes the wifi connection
%
%  open         Open a wifi connection to the brick
%  close        Close the wifi connection to the brick
%  read         Read data from the brick through wifi
%  write        Write data to the brick through wifi
%
% Example::
%           wfbrick = wfBrickIO(1,'192.168.1.104',5555,'0016533dbaf5')
%
% Notes::
% - Java string casting is inconsistent in MATLAB across different
% operating systems. The purpose of OutputStreamSend is to ensure that the
% correct string sequence is sent to the brick, no matter what operating
% system MATLAB is running under.
% - The OutputStreamSend class file will need to be present in the java
% path before a wfBrickIO object is made.
% - This class works across all 3 operating systems.

classdef wfBrickIO < BrickIO
    properties
        % connection handle
        handle
        % debug input
        debug = 0;
        % socket input strem
        inputStream
        % socket output stream
        outputStream
        % brick IP address
        addr = '192.168.1.104';
        % brick tcp port (default is 5555)
        port = 5555;
        % brick serial number
        serialNum = '0016533dbaf5';
    end

    methods

        function brickIO = wfBrickIO(debug,addr,port,serialNum)
            %wfBrickIO.wfBrickIO Create a wfBrickIO object
            %
            % wfbrick = wfBrickIO(debug,addr,port,serialNum) is an object
            % which initialises and opens a wifi connection between MATLAB
            % and the brick using java functions. The following string is
            % sent to the brick
            % ['GET /target?sn=' brickIO.serialNum ' VMTP1.0' char(10) 'Protocol: EV3']
            % which initialises the conncetion. The brick returns
            % ['Accept:EV340']
            %
            % Notes::
            % - debug is a flag specifying output printing (0 or 1).
            % - addr is the IP address of the brick
            % - port is the TCP port of the brick (default is 5555)
            % - serialNum is the serial number of the brick. The serial
            % number can be found in the menu on the EV3 brick or obtained
            % through the emitted UDP packet on port 3015.

            if nargin > 0
                brickIO.debug = debug;
                brickIO.addr = addr;
                brickIO.port = port;
                brickIO.serialNum = serialNum;
            end
            if brickIO.debug > 0
                fprintf('wfBrickIO init\n');
            end
            % import the required java libraries
            import java.io.*;
            import java.net.*;
            % add the java path to the class path
            javaaddpath .
            % open the brick IO connection
            brickIO.handle = Socket(brickIO.addr, brickIO.port);
            % set the input stream
            brickIO.inputStream = DataInputStream(brickIO.handle.getInputStream);
            % set the output stream
            brickIO.outputStream = OutputStreamSend();
            % connection message to the brick
            wmsg = ['GET /target?sn=' brickIO.serialNum ' VMTP1.0' char(10) 'Protocol: EV3'];
            if brickIO.debug > 0
                fprintf('Sent: [ %s ]\n',char(wmsg));
            end
            % send the connection message
            brickIO.write(wmsg);
            % receive the reply
            rmsg = brickIO.read;
            if brickIO.debug > 0
                fprintf('Returned: [ %s ]\n',char(rmsg));
            end
        end

        function delete(brickIO)
            %wfBrickIO.delete Delete the wfBrickIO object
            %
            % delete(brickIO) closes the wifi connection handle

            if brickIO.debug > 0
                fprintf('wfBrickIO delete\n');
            end
            % delete the wf handle
            brickIO.close
        end

        function handle = open(brickIO)
            %wfBrickIO.open Open the wfBrickIO object
            %
            % handle = wfBrickIO.open() returns a handle to the wifi
            % socket connection using the initialised IP address and TCP
            % port values from the wfBrickIO constructor.

            if brickIO.debug > 0
                fprintf('wfBrickIO open\n');
            end
            % open the socket handle
            handle = Socket(brickIO.addr, brickIO.port);
        end

        function close(brickIO)
            %wfBrickIO.close Close the wfBrickIO object
            %
            % wfBrickIO.close() closes the wifi socket connection

            if brickIO.debug > 0
                fprintf('wfBrickIO close\n');
            end
            % close the close handle
            brickIO.handle.close();
        end

        function rmsg = read(brickIO)
            %wfBrickIO.read Read data from the wfBrickIO object
            %
            % rmsg = wfBrickIO.read() reads data from the brick through wifi
            % and returns the data in uint8 format.
            %
            % Notes::
            % - This function is blocking with no time out in the current
            % implementation.

            if brickIO.debug > 0
                fprintf('wfBrickIO read\n');
            end
            % block until bytes have been received
            while (brickIO.inputStream.available == 0)

            end
            % get the number of bytes to be read from the input stream
            nLength = brickIO.inputStream.available;
            % read the bytes from the input stream
            rmsg = zeros(1,nLength);
            for ii=1:nLength
                rmsg(ii) =  brickIO.inputStream.readByte;
            end
            % convert from double to int8
            rmsg = cast(rmsg,'int8');
            % change from int8 to uint8 (cannot do this in one step)
            rmsg = typecast(rmsg,'uint8');
        end

        function write(brickIO,wmsg)
            %wfBrickIO.write Write data to the wfBrickIO object
            %
            % wfBrickIO.write(wmsg) writes data to the brick through wifi.
            %
            % Notes::
            % - wmsg is the data to be written to the brick via wifi in
            % uint8 format.

            if brickIO.debug > 0
                fprintf('wfBrickIO write\n');
            end
            % add the message to the output stream
            brickIO.outputStream.addtoBufferN(char(wmsg),length(wmsg));
            % send the output stream data
            brickIO.outputStream.send(brickIO.handle.getOutputStream());
            % clear the output stream
            brickIO.outputStream.clear();
        end
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function flux = EulerRoe2D(nx, ny, QM, QP, gamma)

% function flux = EulerRoe2D(nx, ny, QM, QP, gamma)
% Purpose: compute surface fluxes for Euler's equations using the approximate
%          Riemann solver based on Roe averages

% Rotate "-" trace momentum to face normal-tangent coordinates
rhouM = QM(:,:,2); rhovM = QM(:,:,3);
QM(:,:,2) = nx.*rhouM + ny.*rhovM; QM(:,:,3) =-ny.*rhouM + nx.*rhovM;

% Rotate "+" trace momentum to face normal-tangent coordinates
rhouP = QP(:,:,2); rhovP = QP(:,:,3);
QP(:,:,2) = nx.*rhouP + ny.*rhovP; QP(:,:,3) =-ny.*rhouP + nx.*rhovP;

% Compute fluxes and primitive variables in rotated coordinates
[fxQM,fyQM,rhoM,uM,vM,pM] = EulerFluxes2D(QM, gamma);
[fxQP,fyQP,rhoP,uP,vP,pP] = EulerFluxes2D(QP, gamma);
EnerM = QM(:,:,4); EnerP = QP(:,:,4);

Nfields = 4;
HM   = (EnerM+pM)./rhoM; cM = sqrt(gamma*pM./rhoM);
HP   = (EnerP+pP)./rhoP; cP = sqrt(gamma*pP./rhoP);

% Compute Roe average variables
rhoMs = sqrt(rhoM); rhoPs = sqrt(rhoP);

rho = rhoMs.*rhoPs;
u   = (rhoMs.*uM + rhoPs.*uP)./(rhoMs + rhoPs);
v   = (rhoMs.*vM + rhoPs.*vP)./(rhoMs + rhoPs);
H   = (rhoMs.*HM + rhoPs.*HP)./(rhoMs + rhoPs);

c2  = (gamma-1)*(H - 0.5*(u.^2 + v.^2)); c = sqrt(c2);

% Compute estimate of waves speeds
SL = min(uM-cM, u-c); SR = max(uP+cP, u+c);

% Compute HLL flux
t1 = (min(SR,0)-min(0,SL))./(SR-SL);
t2 = 1-t1;
t3 = (SR.*abs(SL)-SL.*abs(SR))./(2*(SR-SL));

for n=1:4
  fx(:,:,n) = t1.*fxQP(:,:,n) + t2.*fxQM(:,:,n) - t3.*(QP(:,:,n)-QM(:,:,n));
end

% rotate flux back into Cartesian coordinates
flux(:,:,1) = fx(:,:,1);
flux(:,:,2) = nx.*fx(:,:,2) - ny.*fx(:,:,3);
flux(:,:,3) = ny.*fx(:,:,2) + nx.*fx(:,:,3);
flux(:,:,4) = fx(:,:,4);;
return;

repo: Aero-matlab


text:
function theta = moment_estimation( M, default_beta, default_dev )

% Parameter estimates from moment equations
% Note that parameter estimation is bad for small data sets: The
% variances can be negative and the beta's are numerically too large.
%
% Syntax:
%   E = moment_estimation( M, n_dev[, def_b, max_v] )
%
% Input:
%   M       : Matrix with moments from compute_moments.
%
%   def_b   : If beta is estimated to be negative or too large it is
%             replaced by DEF_B
%			  Default: 1
%
%   def_std : If sigma and kappa are estimated to be negative they are
%             replaced by STD_VAL.
%             Default: 0.5
%
%
% Output:
%   theta : L-by-5-by-D matrix where L is the number of levels in the
%           wavelet transform, D is the number of directions and the
%           (l+1)'th row is estimates of
%           mu_l, sigma_l, alpha_l, beta_l, kappa_l
%
% See also: COMPUTE_MOMENTS

% Default values
if ~exist( 'default_beta', 'var' )
    default_beta = 1;
end

if ~exist( 'std_val', 'var' )
    default_dev = 0.5;
end

% Preallocate output
L = size(M, 1) - 1;
theta = zeros( L+1, 5, size(M,3) );

% Indices for beta's, alpha's and kappa's
idx = 1:L;

% sigma's
sigma2 = log(M(:,2,:)/3) - 2*log(M(:,1,:));
theta(:, 2, :) = sqrt( sigma2 );

% mu's
theta(:, 1, :) = log(M(:,1,:)) - sigma2/2;

% TODO: Is there a good way to test if the beta estimate is too large

% beta's
theta(1+idx, 4, :) = ( log(M(1+idx,3,:)) - 2*log(M(1+idx,1,:)) ) ./ ...
    ( log(M(idx,4,:)) - log(M(1+idx,1,:)) - log(M(idx,1,:)) );

% Restrict large beta's
beta_idx = false( size(theta) );
beta_idx(1+idx, 4, :) = true;

positive_idx = real(theta) > 2;
negative_idx = real(theta) < 0;

theta( beta_idx & (positive_idx | negative_idx) ) = default_beta;

% alpha's
theta(1+idx, 3, :) = theta(1+idx, 1, :) - theta(1+idx, 4, :).*theta(idx, 1, :);

% kappa's
kappa2 = sigma2(1+idx, :, :) - theta(1+idx, 4, :).^2 .* sigma2(idx, :, :);
theta(1+idx, 5, :) = sqrt( kappa2 );

% ----------------------------------------------------------------------
%
% If any of the variances are estimated to be negative, we estimate
% parameters in the homogeneous GLG model

ind = find( imag(theta) );
[~, ~, D] = ind2sub( size(theta), ind );
D = unique( D );

if isempty(D)
    return
end

mu    = theta(1,1,:);
sigma = theta(1,2,:);
alpha = theta(2,3,:);
beta  = theta(2,4,:);
kappa = theta(2,5,:);

% Default values
sigma( imag(sigma) ~= 0 ) = default_dev;
kappa( imag(kappa) ~= 0 ) = default_dev;
alpha( imag(alpha) ~= 0 ) = -default_beta;
beta( imag(beta) ~= 0 )   = default_beta;

homogeneous_theta = homogeneous_to_full_GLG( L+1, mu, sigma, alpha, beta, kappa );
theta(:,:,D) = homogeneous_theta(:,:,D);

end

author: Sable

repo: mcbench-benchmarks


text:
%CMDL_HAT Create a CollisionModel object of the human HAT
%
% Returns a CollisionModel object of the human head and torso segment
% (also includes neck). The model may be given for any shoulder frame
% transformation (synonymous with the base transform of a HAL object).
% The head is a sphere, the neck a cylinder, the shoulders an ellipsoid
% and the torso an elliptical cylinder.
%
% Copyright (C) Bryan Moutrie, 2013-2014
% Licensed under the GNU Lesser General Public License
% see full file for full statement
%
% Syntax:
%  (1) hat = cmdl_hat(Tg)
%  (2) hat = cmdl_hat()
%
%  (2) is as per (1), using HAL's default value of Tg
%
% Outputs:
%  hat : CollisionModel object of the human HAT
%
% Inputs:
%  Tg : Transformation frame of the shoulder
%
% See also anthroData cmdl_arm HAL

% LICENSE STATEMENT:
%
% This file is part of pHRIWARE.
%
% pHRIWARE 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.
%
% pHRIWARE 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 Lesser General Public
% License along with pHRIWARE.  If not, see <http://www.gnu.org/licenses/>.

function hat = cmdl_hat(Tg)

if ~nargin, Tg = trotx(pi/2); end

skin = pHRIWARE('skin');
shirt = pHRIWARE('orange');

% Anthropometric Data
[refnote, pb, sb, pn, sn, pg, sg, pt, st] = ...
    anthroData('pb', 'rb', 'pn', 'sn', 'pg', 'sg', 'pt', 'st');

% Head
gTb = [eye(3), pb; 0, 0, 0, 1];
Tb = Tg * gTb;
Head = Sphere(Tb, sb, 'FaceColor', skin, 'EdgeColor', 'none');

% Neck
gTn = [rotx(-pi/2), pn; 0, 0, 0, 1];
Tn = Tg * gTn;
Neck = Cylinder(Tn, sn, 'FaceColor', skin, 'EdgeColor', 'none');

% Shoulders
gT = [eye(3), pg; 0, 0, 0, 1];
T = Tg * gT;
Shoulders = Ellipsoid(T, sg, 'FaceColor', shirt, 'EdgeColor', 'none');

% Torso
gTt = [rotx(pi/2), pt; 0, 0, 0, 1];
Tt = Tg * gTt;
Torso = Cylinder(Tt, st, 'FaceColor', shirt, 'EdgeColor', 'none');

hat = CollisionModel(refnote, Head, Neck, Shoulders, Torso);

end

author: Allopart

repo: rbpf-gmapping


text:
function f2=f2(x)
Bound=[-10 10];

if nargin==0
    f2 = Bound;
else
    f2=sum(abs(x))+prod(abs(x));
end

author: sfvsfv

repo: Mathematical-modeling


text:
function sys = flatten(sys)

% Go from an internal format which is hierarchical and performs better
% when adding many constraint objects.
if length(sys.clauses)>0
    if isa(sys.clauses{1},'cell')
        sys.clauses = [sys.clauses{:}];
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function test_bug2558

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

nsubj = 23;
nchan = 274;
ntime = 1560;

gavg_erf7.trial = randn(nsubj,nchan,ntime);
%gavg_erf7.individual = randn(nsubj,nchan,ntime);
gavg_erf7.time = (1:ntime)/1000;
for i=1:nchan
  gavg_erf7.label{i} = num2str(i);
end
gavg_erf7.dimord = 'subj_chan_time';
%gavg_erf7.avg = randn(nchan, ntime);
gavg_erf7.cfg = [];

gavg_erf8.trial = randn(nsubj,nchan,ntime);
%gavg_erf8.individual = randn(nsubj,nchan,ntime);
gavg_erf8.time = (1:ntime)/1000;
for i=1:nchan
  gavg_erf8.label{i} = num2str(i);
end
gavg_erf8.dimord = 'subj_chan_time';
%gavg_erf8.avg = randn(nchan, ntime);
gavg_erf8.cfg = [];

cfg = [];
cfg.channel = 'all'; % modification w.r.t. original bug report
cfg.latency = [0.4 0.7];
cfg.avgovertime = 'no';
cfg.avgoverchan = 'no';
%cfg.avgoverfreq = 'yes'; % <--- is this necessary?
cfg.parameter = 'trial';
cfg.method = 'montecarlo';
%cfg.correctm = 'cluster';
cfg.statistic = 'ft_statfun_depsamplesT';
cfg.alpha = 0.05;
%cfg.clusteralpha = 0.05;
%cfg.neighbours  = []; % modification w.r.t. original bug report
cfg.correctm    = 'no';
cfg.correcttail = 'prob';
cfg.numrandomization = 1000;

cfg.design(1,1:2*nsubj) =  [ones(1,nsubj) 2*ones(1,nsubj)];
cfg.design(2,1:2*nsubj) = [1:nsubj 1:nsubj];
cfg.ivar                = 1; % the 1st row in cfg.design contains the independent variable
cfg.uvar                = 2; % the 2nd row in cfg.design contains the subject numbe

stat = ft_timelockstatistics(cfg, gavg_erf7,gavg_erf8);

author: fieldtrip


text:
function [ g ] = gsp_design_simoncelli(G, param)
%GSP_DESIGN_SIMONCELLI Create a Simoncelli filterbank
%   Usage: g = gsp_design_simoncelli( G );
%          g = gsp_design_simoncelli( G, param );
%
%   Inputs parameters:
%       G       : Graph structure or lmax
%       param   : Structure of optional parameters
%
%   Outputs parameters:
%       g       : filterbank
%
%   This function creates a parseval filterbank of $2$ filters. The low-pass
%   filter is defined by a function $f_l(x)$:
%
%   ..              /  1                                   if x <= a
%   ..    f_l(x) = |   cos( pi/2 * log(x/a) / log(2) )     if a < x <= 2a
%   ..              \  0                                   if x > 2a
%
%   .. math:: f_{l}=\begin{cases} 1 & \mbox{if }x\leq a\\ \cos\left(\frac{\pi}{2}\frac{\log\left(\frac{x}{2}\right)}{\log(2)}\right) & \mbox{if }a<x\leq2a\\ 0 & \mbox{if }x>2a \end{cases}
%
%   The high pass filter is adaptated to obtain a tight frame.
%
%   This function will compute the maximum eigenvalue of the laplacian. To
%   be more efficient, you can precompute it using::
%
%       G = gsp_estimate_lmax(G);
%
%   Example:::
%
%         G = gsp_sensor(100);
%         G = gsp_estimate_lmax(G);
%         g = gsp_design_simoncelli(G);
%         gsp_plot_filter(G,g);
%         [A,B] = gsp_filterbank_bounds(G,g)
%
%   *param* is an optional structure containing the following fields
%
%   * *param.verbose*: verbosity level. 0 no log - 1 display warnings.
%     (default 1)
%   * *param.a*: see equations above for this parameter. Note that the
%     spectrum is scaled between 0 and 2 (default 2/3).
%

% Author: Nathanael Perraudin, David Shuman
% Date  : 21 June 2014
% Testing: test_filter

if nargin < 2
    param = struct;
end

if ~isfield(param,'verbose'), param.verbose = 1; end
if ~isfield(param,'a'), param.a = 2/3; end

if isstruct(G)
    if ~isfield(G,'lmax')
        if param.verbose
            fprintf('GSP_DESIGN_SIMONCELLI has to compute lmax \n')
        end
        G = gsp_estimate_lmax(G);
    end
   lmax = G.lmax;
else
   lmax = G;
end

a = param.a;

g = cell(2,1);
g{1} = @(x) simoncelli(x*(2/lmax),a);
g{2} = @(x) real(sqrt(1-(simoncelli(x*(2/lmax),a)).^2));

end

function y = simoncelli(val,a)

y = zeros(size(val));

l1 = a;
l2 = 2*a;

r1ind = val >= 0     &    val < l1;
r2ind = val >= l1    &    val < l2;
r3ind = val >= l2;

y(r1ind) = 1;
y(r2ind) = cos(pi/2*log(val(r2ind)/a)/log(2));
y(r3ind) = 0;

end

repo: gspbox


text:
function out = chebcoeffs(f, varargin)
%CHEBCOEFFS   Chebyshev polynomial coefficients of a CLASSICFUN.
%   CHEBCOEFFS(F) returns the Chebyshev coefficients of F.ONEFUN.
%
% See also LEGPOLY, TRIGCOEFFS.

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

% Call CHEBCOEFFS() of the .ONEFUN:
out = chebcoeffs(f.onefun, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function const_pts = RectQAM_const(M)
%RectQAM_const Rectangular QAM Constellation points with Gray mapping.
%   C = rectQAM_const(M) returns the M-ary rectangular QAM constellation
%   points for M = 2^k, where k=2,3,...,10.
%
%   For more information, see
%   [1] Cho, K., and Yoon, D., "On the general BER expression of one- and
%        two-dimensional amplitude modulations", IEEE Trans. Commun.,
%        Vol. 50, Number 7, pp. 1074-1080, 2002.
%
%   See also PAM_Gray_Code, QAMMOD.

%   Written by Idin Motedayen-Aval
%   Applications Engineer
%   The MathWorks, Inc.
%   zq=[4 2 5 -15 -1 -3 24 -57 45 -12 19 -12 15 -8 3 -7 8 -69 53 12 -2];
%   char(filter(1,[1,-1],[105 zq])), clear zq

% This function constructs the rectangular QAM constellations by
% doing two PAM modulations: one for the in-phase and one for the
% quadrature dimension (see [1] for more details).

% M-PAM bit/symbol ordering
[bo so] = PAM_Gray_Code;

k = log2(M);
I = 2^ceil(k/2);
Q = 2^floor(k/2);

% Produce PAM constellation points for in-phase and quadrature
I_t = -(I-1):2:I-1;
Q_t = Q-1:-2:-(Q-1);

% Re-order based on the PAM Gray mapping
I_pts = zeros(size(I_t)); % pre-allocate
Q_pts = zeros(size(Q_t)); % pre-allocate
for i=1:length(I_t)
    I_pts(so{log2(I)}(i)+1) = I_t(i);
end
for i=1:length(Q_t)
    Q_pts(so{log2(Q)}(i)+1) = Q_t(i);
end

% Map the symbols 0:M-1 to constellation points
QAM_I=zeros(1,M); % pre-allocate
QAM_Q=zeros(1,M);
for i=0:M-1
    MSBs = floor(i/I);
    LSBs = i-MSBs*I;
    QAM_I(i+1) = I_pts(LSBs+1);
    QAM_Q(i+1) = Q_pts(MSBs+1);
end

const_pts=complex(QAM_I,QAM_Q);


text:
%FBA
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 3/24/2011
classdef FBA
    %printing
    methods (Static)
        function run(sim, fileName)
            import edu.stanford.covert.cell.sim.analysis.FBA;
            import edu.stanford.covert.cell.sim.util.PlotUtil;
            import edu.stanford.covert.cell.sim.util.PrintUtil;

            %excel file
            [content, colLabels, indentation] = FBA.printNetworkReduction(sim);
            if nargin == 1
                PrintUtil.printToStdIO(content, colLabels, struct('indentation', indentation));
            else
                PrintUtil.printToFile(content, colLabels, [fileName '.xls'], 'NetworkReduction', struct('indentation', indentation));
            end

            %plots
            if nargin == 1
                FBA.plotNetworkReduction(sim, PlotUtil.newAxesHandle());
            else
                [axesHandle, figHandle] = PlotUtil.newAxesHandle();

                cla(axesHandle);
                FBA.plotNetworkReduction(sim, axesHandle);
                saveas(figHandle, [fileName '-NetworkReduction.pdf']);

                close(figHandle);
            end
        end

        function [content, colLabels, indentation] = printNetworkReduction(sim)
            p = sim.process('Metabolism');

            [subCmpRemoved, rxnRemoved] = p.formulateFBA([], [], true);

            content = cell(0, 6);
            colLabels = {'ID', 'Compartment', 'Reason', 'Iteration-1', 'Iteration-2'};
            reasons = {
                'Disconnected'
                'Substrate cannot be balanced: only appears in single reaction'
                'Substrate cannot be balanced: move in same direction in all reactions in which it appears'
                'Reaction has zero flux in all solutions of S*v = 0: reaction has no non-zero component in null space of S';
                };

            subIDs = p.substrateWholeCellModelIDs;
            cmpIDs = p.compartment.wholeCellModelIDs(p.substrateMetaboliteCompartmentIndexs(1, :));
            content = [content;{
                0 'Substrates' [] [] [] []}];
            for i = 1:size(subCmpRemoved, 1)
                if subCmpRemoved(i, 1) == 0
                    continue;
                end
                [sIdx, cIdx] = ind2sub([size(p.reactionStoichiometryMatrix, 1) size(p.reactionStoichiometryMatrix, 3)], i);
                content = [content;{
                    1 subIDs{sIdx} cmpIDs{cIdx} reasons{subCmpRemoved(i, 3)} subCmpRemoved(i, 1) subCmpRemoved(i, 2)}]; %#ok<AGROW>
            end

            rxnIDs = p.reactionWholeCellModelIDs;
            content = [content;{
                0 'Reactions' [] [] [] []}];
            for i = 1:size(rxnRemoved, 1)
                if rxnRemoved(i, 3) == 0
                    continue;
                end
                content = [content;{
                    1 rxnIDs{i} [] reasons{rxnRemoved(i, 3)} rxnRemoved(i, 1) rxnRemoved(i, 2)}]; %#ok<AGROW>
            end

            %format output
            indentation = cell2mat(content(:, 1));
            content = content(:, 2:end);
        end

        function plotNetworkReduction(sim, axesHandle)
            p = sim.process('Metabolism');

            [subCmpRemoved, rxnRemoved] = p.formulateFBA([], [], true);
            subCmpRemoved(subCmpRemoved(:, 1) == 0, 1) = max(subCmpRemoved(:, 1)) + 1;
            rxnRemoved(rxnRemoved(:, 1) == 0, 1) = max(rxnRemoved(:, 1)) + 1;

            levels = sortrows(unique([subCmpRemoved; rxnRemoved], 'rows'), [1 2 3]);

            [~, subCmpIdxs] = sortrows(subCmpRemoved, 1:3);
            [~, rxnIdxs] = sortrows(rxnRemoved, 1:3);

            [~, subCmpLevels] = ismember(subCmpRemoved, levels, 'rows');
            [~, rxnLevels] = ismember(rxnRemoved, levels, 'rows');

            hold on;

            x = [
                0.5
                numel(rxnIdxs)+0.5
                numel(rxnIdxs)+0.5
                0.5
                ];
            y = [
                0.5
                0.5
                numel(subCmpIdxs)+0.5
                numel(subCmpIdxs)+0.5
                ];
            patch(x, y, 1, 'FaceColor', [1 1 1], 'Parent', axesHandle, 'FaceAlpha', 1);
            patch(x, y, 1, 'CDataMapping', 'direct', 'Parent', axesHandle, 'FaceAlpha', 0.25, 'EdgeAlpha', 1);
            for i = 2:size(levels, 1)
                j = find(subCmpLevels(subCmpIdxs) < i, 1, 'last');
                k = find(rxnLevels(rxnIdxs) < i, 1, 'last');
                if isempty(j), j = 0; end;
                if isempty(k), k = 0; end;

                x = [
                    k+0.5
                    numel(rxnIdxs)+0.5
                    numel(rxnIdxs)+0.5
                    k+0.5
                    ];
                y = [
                    j+0.5
                    j+0.5
                    numel(subCmpIdxs)+0.5
                    numel(subCmpIdxs)+0.5
                    ];
                patch(x, y, 1, 'FaceColor', [1 1 1], 'Parent', axesHandle, 'FaceAlpha', 1);
                patch(x, y, i, 'CDataMapping', 'direct', 'Parent', axesHandle, 'FaceAlpha', 0.25, 'EdgeAlpha', 1);
            end

            rxnSMat = reshape(permute(p.reactionStoichiometryMatrix, [2 1 3]), ...
                size(p.reactionStoichiometryMatrix, 2), [])';

            [y, x] = find(...
                ((rxnSMat(subCmpIdxs, rxnIdxs) < 0 & repmat(p.reactionBounds(rxnIdxs, 2)' > 0, numel(subCmpIdxs), 1)) | ...
                (rxnSMat(subCmpIdxs, rxnIdxs) > 0 & repmat(p.reactionBounds(rxnIdxs, 1)' < 0, numel(subCmpIdxs), 1))) & ...
                ~((rxnSMat(subCmpIdxs, rxnIdxs) > 0 & repmat(p.reactionBounds(rxnIdxs, 2)' > 0, numel(subCmpIdxs), 1)) | ...
                (rxnSMat(subCmpIdxs, rxnIdxs) < 0 & repmat(p.reactionBounds(rxnIdxs, 1)' < 0, numel(subCmpIdxs), 1))));
            h1 = plot(axesHandle, x, y, 'r.', 'MarkerSize', 6);

            [y, x] = find(...
                ~((rxnSMat(subCmpIdxs, rxnIdxs) < 0 & repmat(p.reactionBounds(rxnIdxs, 2)' > 0, numel(subCmpIdxs), 1)) | ...
                (rxnSMat(subCmpIdxs, rxnIdxs) > 0 & repmat(p.reactionBounds(rxnIdxs, 1)' < 0, numel(subCmpIdxs), 1))) &  ...
                ((rxnSMat(subCmpIdxs, rxnIdxs) > 0 & repmat(p.reactionBounds(rxnIdxs, 2)' > 0, numel(subCmpIdxs), 1)) | ...
                (rxnSMat(subCmpIdxs, rxnIdxs) < 0 & repmat(p.reactionBounds(rxnIdxs, 1)' < 0, numel(subCmpIdxs), 1))));
            h2 = plot(axesHandle, x, y, 'g.', 'MarkerSize', 6);

            [y, x] = find(...
                ((rxnSMat(subCmpIdxs, rxnIdxs) < 0 & repmat(p.reactionBounds(rxnIdxs, 2)' > 0, numel(subCmpIdxs), 1)) | ...
                (rxnSMat(subCmpIdxs, rxnIdxs) > 0 & repmat(p.reactionBounds(rxnIdxs, 1)' < 0, numel(subCmpIdxs), 1))) &  ...
                ((rxnSMat(subCmpIdxs, rxnIdxs) > 0 & repmat(p.reactionBounds(rxnIdxs, 2)' > 0, numel(subCmpIdxs), 1)) | ...
                (rxnSMat(subCmpIdxs, rxnIdxs) < 0 & repmat(p.reactionBounds(rxnIdxs, 1)' < 0, numel(subCmpIdxs), 1))));
            h3 = plot(axesHandle, x, y, 'b.', 'MarkerSize', 6);

            legend([h1 h2 h3], {'Drained'; 'Produced'; 'Bidirectional'}, 'Location', 'NorthEastOutside');

            colormap(jet);

            xlim([0.5 numel(rxnIdxs)+0.5])
            ylim([0.5 numel(subCmpIdxs)+0.5]);

            box('on');
            xlabel('Reactions', 'FontSize', 12);
            ylabel('Substrates', 'FontSize', 12);
            set(axesHandle, 'YDir', 'reverse');
            set(axesHandle, 'XAxisLocation', 'top');
            axis('square')
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function result=gcaleventor(userName, password, title, content, location, smsMe, emailMe, DelayInMin)
% author: Ofir Bibi
% This script uses the GData interface for the Google Calendar.
% Using this script you can create a new event in your calendar.
% userName - your google username (including gmail.com domain)
% password - your google password
% title - The title of the event
% content - The content of the event
% location - the location of the event
% smsMe - true if you want an sms reminder
% emailMe - true if you want an email reminder
% DelayInMin - gap in minutes from now to the time of the event.
% returns 0 if the event was created, negative number for errors.

import java.io.*;
import java.net.*;
import java.lang.*;
import java.text.DateFormat;
import java.text.SimpleDateFormat;
import java.util.Date;

MSMIN=60000;
MAXITER=10;
result=0;

eventpre=['<entry xmlns=''http://www.w3.org/2005/Atom''' ...
' xmlns:gd=''http://schemas.google.com/g/2005''>' ...
'<category scheme=''http://schemas.google.com/g/2005#kind''' ...
' term=''http://schemas.google.com/g/2005#event''></category>' ...
'<title type=''text''>' title '</title>' ...
'<content type=''text''>' content '</content>' ...
'<gd:transparency value=''http://schemas.google.com/g/2005#event.opaque''>' ...
'</gd:transparency><gd:eventStatus' ...
' value=''http://schemas.google.com/g/2005#event.confirmed''>' ...
'</gd:eventStatus><gd:where valueString=''' location '''></gd:where>' ...
'<gd:when startTime='''];
eventmid=''' endTime=''';
eventend='''>';
if (smsMe)
   eventend=[eventend '<gd:reminder method=''sms'' minutes=''0''/>'];
end
if (emailMe)
    eventend=[eventend '<gd:reminder method=''email'' minutes=''0''/>'];
end
eventend=[eventend '</gd:when></entry>'];
posturlString = ['http://www.google.com/calendar/feeds/' userName '/private/full'];

[authorized, aToken]=connectNAuthorize(userName, password);
if (~authorized)
    result=-1;
    return;
end
success=false;
try
    % try to create new event
    safeguard=0;
    while (~success && safeguard<MAXITER)
        safeguard=safeguard+1;
        url = URL( posturlString );
        con = url.openConnection();
        con.setInstanceFollowRedirects(false);
        % set up url connection to get retrieve information back
        con.setRequestMethod( 'POST' );
        con.setDoOutput( true );
        con.setDoInput( true );
        con.setRequestProperty('content-type','application/atom+xml;charset=UTF-8');
        % stuff the Authorization request header
        con.setRequestProperty('Authorization',String('GoogleLogin ').concat(aToken));
        now = Date();
        now.setTime(now.getTime()+MSMIN*(1+DelayInMin));
        df = SimpleDateFormat('yyyy-MM-dd HH:mm:ss');
        time=df.format(now);
        time=time.replace(' ', 'T');
        event=[eventpre char(time) eventmid char(time) eventend];
        ps = PrintStream(con.getOutputStream());
        ps.print(event);
        ps.close();
        % pull the information back from the URL
        is = con.getInputStream();
        if (con.getResponseCode()==302)
            posturlString = con.getHeaderField('Location');
            con.disconnect();
            continue;
        end
        con.disconnect();
        success=true;
    end
catch ME
    disp(ME.message);
    disp(['from line: ' num2str(ME.stack.line)]);
    result = -2;
    return;
end
if (~success)
    result = -3;
end

end

function [authorized, aToken]=connectNAuthorize(userName, password)
import java.io.*;
import java.net.*;
import java.lang.*;
loginurlString = 'https://www.google.com/accounts/ClientLogin';
params = {'Email', 'Passwd', 'source', 'service'};
vals = {userName,password,'OPB-MatlabCal-1','cl'};
authorized = true;
try
    % open url connection
    url = URL( loginurlString );
    con = url.openConnection();
    % set up url connection to get retrieve information back
    con.setRequestMethod( 'POST' );
    con.setDoOutput( true );
    con.setRequestProperty('content-type','application/x-www-form-urlencoded');
    ps = PrintStream(con.getOutputStream());
    for i=1:length(params)
        if (i > 0)
            ps.print('&');
        end
        param = URLEncoder.encode(params{i}, 'UTF-8');
        value = URLEncoder.encode(vals{i}, 'UTF-8');
        ps.print(param.concat('=').concat(value));
    end
    ps.close();
    % pull the information back from the URL
    is = con.getInputStream();
    buf = java.lang.StringBuffer();
    c = int32(is.read());
    while( c ~= -1 )
        buf.append( char(c) );
        c = int32(is.read());
    end
    con.disconnect();
    % get the Authorization Token
    aToken=java.lang.String('auth').concat(buf.substring(buf.lastIndexOf('Auth')+4,buf.length()-1));
catch ME
    disp(['from line: ' num2str(ME.stack.line)]);
    error(ME.message);
end

end


text:
function [Xs, delta, C, B, names] = get_session_X(obj, s)
% Get design matrix (predictors) for one session of fmri_model object, using
% basis functions defined in the object and onsets for one session (s).
%
% :Usage:
% ::
%
%     [Xs, delta, C, B, names] = get_session_X(obj, session number)
%

% ..
%    Define sessions and number of conditions
% ..

nsess = length(obj.Sess);

if s > nsess, error('Session %3.0f does not exist', s); end

TR = obj.xY.RT;

nconds = length(obj.Sess(s).U);

[ons, name] = deal(cell(1, nconds));

[ons{:}] = deal(obj.Sess(s).U(:).ons);
[name{:}] = deal(obj.Sess(s).U(:).name);

% make sure onsets are in TRs, not secs
switch obj.xBF(1).UNITS
    case 'secs'
        % onsets are in sec, convert
        for i = 1:nconds
            ons{i} = ons{i} ./ TR;
        end

    case {'tr', 'TR', 'trs'}
        % ok, do nothing
end

% ----------------------------------------------
% Predictors
% ----------------------------------------------

delta = onsets2delta(ons, obj.nscan(s));

% Of-interest part of design matrix
% time res is defined as TR / 16, so build and downsample by 16 to TR
% Allow for different basis sets for each condition.

Xs = cell(1, nconds);

for i = 1:nconds

    bf = obj.xBF(i).bf;

    Xs{i} = getPredictors(delta(:, i), bf, 16);
end

Xs = cat(2, Xs{:});

% covariate part of design matrix
C = [];
if ~isempty(obj.Sess(s).C) && isfield(obj.Sess(s).C, 'C')
    C = obj.Sess(s).C.C;
end

% baseline
% ----------------------------------------------

B = ones(obj.nscan(s), 1);

% ----------------------------------------------
% modulators
% ----------------------------------------------

Xs_pm = cell(1, nconds);

for i = 1:nconds
    is_pm = isfield(obj.Sess(s).U(i), 'P') && ~isempty(obj.Sess(s).U(i).P) && isfield(obj.Sess(s).U(i).P, 'P') && ~isempty(obj.Sess(s).U(i).P.P);

    if is_pm

        pm_vals = {obj.Sess(s).U(i).P.P};

        model = onsets2parametric_mod_X(ons(i), pm_vals, obj.nscan(s), obj.xBF(i).bf, 16);

        Xs_pm{i} = model;

    end

end

Xs_pm = cat(2, Xs_pm{:});

% add to Xs
Xs = [Xs Xs_pm];

% ----------------------------------------------
% get names for each BF
% Allow variable number of basis fcns for each condition
% ----------------------------------------------

names = {};

for i = 1:nconds

    nbf = size(obj.xBF(i).bf, 2);

    for j = 1:nbf
        myname = [name{i} ' BF' num2str(j)];
        myname = replaceblanks(myname);

        names{end+1} = myname;
    end
end

all_pmnames = {};

for i = 1:nconds
    is_pm = isfield(obj.Sess(s).U(i), 'P') && ~isempty(obj.Sess(s).U(i).P) && isfield(obj.Sess(s).U(i).P, 'P') && ~isempty(obj.Sess(s).U(i).P.P);

    if is_pm
        % Get param mod names
        % ------------------------------------------------
        pmnames = cell(1, nbf);

        for j = 1:nbf
            myname = [obj.Sess(s).U(i).P.name ' BF' num2str(j)];
            myname = replaceblanks(myname);

            pmnames{1, j} = myname;
        end

        all_pmnames = cat(2, all_pmnames, pmnames{:});

    end
end

names = names(:)';

names = [names all_pmnames];

end

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

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

function myname = replaceblanks(myname)

myname(myname == ' ') = '-';
% this works in 2010b but not a...
%tmp = diff(num2str(myname)) == 0;

tmp = diff(double(myname)) == 0;
myname([false tmp] & myname == '-') = [];

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function neighborsCell = identifyNeighborsAuto(inputImages, inputSignals, varargin)
	% This code automatically sorts through to find all obj neighbors within a certain distance of the target (boundary to boundary). The output is a cell array with vectors of the neighbor indices to each obj.
	% Biafra Ahanonu
	% started: 2013.11.01
	% based on code by laurie burns, started: sept 2010.
	% inputs
		% inputImages - [x y nCells] matrices containing each set of filters
		% inputSignals - [nFilters frames] matrices containing each set of filter traces
	% options
		% _
	% outputs
		% _
	% changelog
		% 2013.11.01 refactored so accepts ICA filters and traces as input, also changed display to it matches ICAchooser, e.g. it applies a cellmap along with the identified and new cells
		% 2017.06.29 - change to using centroids since faster than imdilate and more accurate.
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	options.plottingOn = 0;
	options.overlapradius = 10;
	% overlap radius in pixels (so convert properly)
	options.overlapDistance = 10;
	% 'centroid', 'imdilate'
	options.neighborMethod = 'centroid';
	%
	options.waitbarOn = 1;
	% Input pre-computed x,y coordinates for objects in images
	options.xCoords = [];
	options.yCoords = [];

	options.inputImagesThres = [];;

	%
	options.startCellNo = 1;
	%
	options.cropSizeLength = 20;

	% get options
	options = getOptions(options,varargin);
	% unpack options into current workspace
	fn=fieldnames(options);
	for i=1:length(fn)
		eval([fn{i} '=options.' fn{i} ';']);
	end
	%========================

	% cell vector info
	cellvec(1) = 1;
	cellvec(2) = size(inputImages,3);

	%% look for overlap of IC and dilated IC
	neighborsCell=cell(cellvec(2),1);
	% reshape to (x,y,z) indexing for compatibility
	% inputImages = permute(inputImages, [2 3 1]);

	reverseStr = '';
	nSignals = cellvec(2);
	if options.plottingOn==1
		figure(222)
	end

	ignoreDistanceReplace = 1e7;

	% [xCoords yCoords] = findCentroid(signalImages,'thresholdValue',0.8,'imageThreshold',0.3);
	% coordsGlobal = [xCoords; yCoords]';
	% distanceMatrix = squareform(pdist(coordsGlobal));
	switch options.neighborMethod
		case 'centroid'
			% body
			% find all cells below distance threshold
			distanceThreshold = options.overlapDistance; % 1px = 1um, 20 um per Jesse's procedures
			nCells = size(inputImages,3);

			% get the centroids and other info for movie
			if isempty(options.xCoords)
				% [xCoords, yCoords] = findCentroid(inputImagesThres,'waitbarOn',options.waitbarOn,'runImageThreshold',0);
				[thresholdedImages boundaryIndices] = thresholdImages(inputImages,'binary',1,'getBoundaryIndex',1,'threshold',0.35,'imageFilter','');
				[xCoords yCoords] = findCentroid(thresholdedImages,'runImageThreshold',0);
			else
				xCoords = options.xCoords(:)';
				yCoords = options.yCoords(:)';
			end

			coordsGlobal = [xCoords; yCoords]';
			% size(coordsGlobal)
			% nCells
			distanceMatrix = squareform(pdist(coordsGlobal));
			% distanceMatrix = distanceMatrix.*~diag(1e5*ones(1,nCells));
			distanceMatrix(logical(eye(size(distanceMatrix)))) = ignoreDistanceReplace;

			distanceMatrixThres = distanceMatrix<distanceThreshold;
			neighborsCell = {};
			reverseStr = '';
			for cellNo = 1:nCells
				neighborIdx = find(distanceMatrixThres(cellNo,:));
				neighborsCell{cellNo,1} = neighborIdx;
				% if mod(cellNo,50)==1
				%     fprintf('up to cell number %d of %d \n',cellNo,nCells)
				% end
				if (cellNo==1||mod(cellNo,50)==0||cellNo==nCells)&options.waitbarOn==1
					reverseStr = cmdWaitbar(cellNo,nCells,reverseStr,'inputStr','Finding cell neighbors');
				end
			end
			fprintf('Done assigning neighbor IDs...\n')
		case 'imdilate'
			% do nothing
		otherwise
			for c = 1:nSignals
				se = strel('disk',overlapradius,0);
				% threshold ICs
				inputImages = thresholdImages(inputImages);
				% get a dilated version of the thresholded image
				thisCellDilateCopy = repmat(imdilate(squeeze(inputImages(:,:,c)),se),[1 1 size(inputImages,3)]);
				% thisCellDilateCopy = permute(thisCellDilateCopy, [3 1 2]);
				% matrix multiple, any overlap will be labeled
				res = inputImages.*thisCellDilateCopy;
				res = squeeze(sum(sum(res,2),1));
				% all cells above threshold are neighbors
				res = find(res>1);
				neighborsCell{c,1} = setdiff(res,c);
				% reduce waitbar access
				if c==1||mod(c,7)==0|c==nSignals
					reverseStr = cmdWaitbar(c,nSignals,reverseStr,'inputStr','identifying neighboring cells');
				end
				if options.plottingOn==1
					keyIn = get(gcf,'CurrentCharacter');
					if strcmp(keyIn,'f')%user wants to exit
						set(gcf,'currentch','3');drawnow;
						keyIn = get(gcf,'CurrentCharacter');
						break
					end
				end
			end
	end
	if plottingOn
		viewNeighborsAuto(inputImages, inputSignals, neighborsCell);
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch)
%MAKEYBUS   Builds the bus admittance matrix and branch admittance matrices.
%   [YBUS, YF, YT] = MAKEYBUS(MPC)
%   [YBUS, YF, YT] = MAKEYBUS(BASEMVA, BUS, BRANCH)
%
%   Returns the full bus admittance matrix (i.e. for all buses) and the
%   matrices YF and YT which, when multiplied by a complex voltage vector,
%   yield the vector currents injected into each line from the "from" and
%   "to" buses respectively of each line. Does appropriate conversions to p.u.
%   Inputs can be a MATPOWER case struct or individual BASEMVA, BUS and
%   BRANCH values. Bus numbers must be consecutive beginning at 1
%   (i.e. internal ordering).
%
%   See also MAKEJAC, MAKESBUS, EXT2INT.

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

%% extract from MPC if necessary
if nargin < 3
    mpc     = baseMVA;
    baseMVA = mpc.baseMVA;
    bus     = mpc.bus;
    branch  = mpc.branch;
end

%% constants
nb = size(bus, 1);          %% number of buses
nl = size(branch, 1);       %% number of lines

%% define named indices into bus, branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;

%% check that bus numbers are equal to indices to bus (one set of bus numbers)
if any(bus(:, BUS_I) ~= (1:nb)')
    error('makeYbus: buses must be numbered consecutively in bus matrix; use ext2int() to convert to internal ordering')
end

%% for each branch, compute the elements of the branch admittance matrix where
%%
%%      | If |   | Yff  Yft |   | Vf |
%%      |    | = |          | * |    |
%%      | It |   | Ytf  Ytt |   | Vt |
%%
stat = branch(:, BR_STATUS);                    %% ones at in-service branches
Ys = stat ./ (branch(:, BR_R) + 1j * branch(:, BR_X));  %% series admittance
Bc = stat .* branch(:, BR_B);                           %% line charging susceptance
tap = ones(nl, 1);                              %% default tap ratio = 1
i = find(branch(:, TAP));                       %% indices of non-zero tap ratios
tap(i) = branch(i, TAP);                        %% assign non-zero tap ratios
tap = tap .* exp(1j*pi/180 * branch(:, SHIFT)); %% add phase shifters
Ytt = Ys + 1j*Bc/2;
Yff = Ytt ./ (tap .* conj(tap));
Yft = - Ys ./ conj(tap);
Ytf = - Ys ./ tap;

%% compute shunt admittance
%% if Psh is the real power consumed by the shunt at V = 1.0 p.u.
%% and Qsh is the reactive power injected by the shunt at V = 1.0 p.u.
%% then Psh - j Qsh = V * conj(Ysh * V) = conj(Ysh) = Gs - j Bs,
%% i.e. Ysh = Psh + j Qsh, so ...
Ysh = (bus(:, GS) + 1j * bus(:, BS)) / baseMVA; %% vector of shunt admittances

%% bus indices
f = branch(:, F_BUS);                           %% list of "from" buses
t = branch(:, T_BUS);                           %% list of "to" buses

%% for best performance, choose method based on MATLAB vs Octave and size
if nb < 300 || have_feature('octave')   %% small case OR running on Octave
    %% build Yf and Yt such that Yf * V is the vector of complex branch currents injected
    %% at each branch's "from" bus, and Yt is the same for the "to" bus end
    i = [1:nl 1:nl]';                           %% double set of row indices
    Yf = sparse(i, [f; t], [Yff; Yft], nl, nb);
    Yt = sparse(i, [f; t], [Ytf; Ytt], nl, nb);

    %% build Ybus
    Ybus = sparse([f;f;t;t], [f;t;f;t], [Yff;Yft;Ytf;Ytt], nb, nb) + ... %% branch admittances
            sparse(1:nb, 1:nb, Ysh, nb, nb);        %% shunt admittance
else                                %% large case running on MATLAB
    %% build connection matrices
    Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);      %% connection matrix for line & from buses
    Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);      %% connection matrix for line & to buses

    %% build Yf and Yt such that Yf * V is the vector of complex branch currents injected
    %% at each branch's "from" bus, and Yt is the same for the "to" bus end
    Yf = sparse(1:nl, 1:nl, Yff, nl, nl) * Cf + sparse(1:nl, 1:nl, Yft, nl, nl) * Ct;
    Yt = sparse(1:nl, 1:nl, Ytf, nl, nl) * Cf + sparse(1:nl, 1:nl, Ytt, nl, nl) * Ct;

    %% build Ybus
    Ybus = Cf' * Yf + Ct' * Yt + ...            %% branch admittances
            sparse(1:nb, 1:nb, Ysh, nb, nb);    %% shunt admittance
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [nfm_Sharp_lunwrap, mask_sharp] = backgroundRemovalSharp(phase_lunwrap, mask_pad, filterMode)
%BACKGROUNDREMOVALSHARP Background phase removal using SHARP (Sophisticated harmonic artifact reduction for phase data)
%   Code refractored from Berkin Bilgic's scripts: "script_Laplacian_unwrap_Sharp_Fast_TV_gre3D.m"
%   and "script_Laplacian_unwrap_Sharp_Fast_TV_gre3D.m"
%   Original source: https://martinos.org/~berkin/software.html
%
%   phase_lunwrap: unwrapped phase volume
%   mask_pad: zero-padded mask volume
%   filterMode: 'once' or 'iterative'
%
%   Original reference:
%   Bilgic et al. (2014), Fast quantitative susceptibility mapping with
%   L1-regularization and automatic parameter selection. Magn. Reson. Med.,
%   72: 1444-1459. doi:10.1002/mrm.25029
%
%   ... which references:
%
%   Li W, Wu B, Liu C. Quantitative susceptibility mapping of
%   human brain reflects spatial variation in tissue composition.
%   NeuroImage 2011; 55(4): 1645?1656.
%

    if nargin < 4
        filterMode = 'once';
    end
    switch filterMode
        case 'once'
            [nfm_Sharp_lunwrap, mask_sharp] = sharp_once(phase_lunwrap, mask_pad);
        case 'iterative'
            [nfm_Sharp_lunwrap, mask_sharp] = sharp_iterative(phase_lunwrap, mask_pad);
    end

end

function [nfm_Sharp_lunwrap, mask_sharp] = sharp_once(phase_lunwrap, mask_pad)

    N = size(mask_pad);

    ksize = [9, 9, 9];                % Sharp kernel size
    threshold = .05;                  % truncation level

    % calculate del kernel and its inverse
    del_sharp = calc_del_kernel(ksize, N);

    delsharp_inv = zeros(size(del_sharp));
    delsharp_inv( abs(del_sharp) > threshold ) = 1 ./ del_sharp( abs(del_sharp) > threshold );

    % erode mask to remove convolution artifacts
    mask_sharp = erode_mask(mask_pad, ksize);

    % apply Sharp to Laplacian wrapped phase
    phase_del = ifftn(fftn(phase_lunwrap) .* del_sharp);
    Phase_Del = phase_del .* mask_sharp;

    phase_Sharp_lunwrap = real( ifftn(fftn(Phase_Del) .* delsharp_inv) .* mask_sharp );

    nfm_Sharp_lunwrap = phase_Sharp_lunwrap;
end

function [nfm_Sharp_lunwrap, mask_sharp] = sharp_iterative(phase_lunwrap, mask_pad)

    N = size(mask_pad);

    threshold = .05;                     % truncation level

    Kernel_Sizes = 9:-2:3;

    % initiate volumes
    Phase_Del = zeros(N);
    mask_prev = zeros(N);

    for k = 1:length(Kernel_Sizes)

        disp(['Kernel size: ', num2str(Kernel_Sizes(k))])

        Kernel_Size = Kernel_Sizes(k);
        ksize = [Kernel_Size, Kernel_Size, Kernel_Size];                % Sharp kernel size

        % calculate del kernel and its inverse
        del_sharp = calc_del_kernel(ksize, N);

        if k == 1
            delsharp_inv = zeros(size(del_sharp));
            delsharp_inv( abs(del_sharp) > threshold ) = 1 ./ del_sharp( abs(del_sharp) > threshold );
        end

        % erode mask to remove convolution artifacts
        mask_sharp = erode_mask(mask_pad, ksize);

        % apply Sharp to Laplacian unwrapped phase
        phase_del = ifftn(fftn(phase_lunwrap) .* del_sharp);
        Phase_Del = Phase_Del + phase_del .* (mask_sharp - mask_prev);

        mask_prev = mask_sharp;

    end

    phase_Sharp_lunwrap = real( ifftn(fftn(Phase_Del) .* delsharp_inv) .* mask_sharp );
    nfm_Sharp_lunwrap = phase_Sharp_lunwrap;

end

function del_sharp = calc_del_kernel(ksize, N)

    khsize = (ksize-1)/2;
    [a,b,c] = meshgrid(-khsize(2):khsize(2), -khsize(1):khsize(1), -khsize(3):khsize(3));

    kernel = (a.^2 / khsize(1)^2 + b.^2 / khsize(2)^2 + c.^2 / khsize(3)^2 ) <= 1;
    kernel = -kernel / sum(kernel(:));
    kernel(khsize(1)+1,khsize(2)+1,khsize(3)+1) = 1 + kernel(khsize(1)+1,khsize(2)+1,khsize(3)+1);

    Kernel = zeros(N);
    Kernel( 1+N(1)/2 - khsize(1) : 1+N(1)/2 + khsize(1), 1+N(2)/2 - khsize(2) : 1+N(2)/2 + khsize(2), 1+N(3)/2 - khsize(3) : 1+N(3)/2 + khsize(3) ) = -kernel;

    del_sharp = fftn(fftshift(Kernel));

end

function mask_sharp = erode_mask(mask_pad, ksize)

    erode_size = ksize + 1;

    mask_sharp = imerode(mask_pad, strel('line', erode_size(1), 0));
    mask_sharp = imerode(mask_sharp, strel('line', erode_size(2), 90));
    mask_sharp = permute(mask_sharp, [1,3,2]);
    mask_sharp = imerode(mask_sharp, strel('line', erode_size(3), 0));
    mask_sharp = permute(mask_sharp, [1,3,2]);

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab