text:
function planC = translateDose(doseNum,xyzT,doseName,planC)
%
% function planC = translateDose(structNum,xyzT,structName,planC)
%
% This function creates a new dose matrix by translating doseNum by amounts
% xT, yT and zT in x,y,z directions respectively. Note that xyzT = [xT yT zT].
% doseName must be a string to name the new dose matrix.
%
% Example:
%   planC = translateDose(1,[-1 0 0],'moveDose1');
%
% Emiliano Spezi, 11 Sep 2013
%
% See also GETUNIFORMSTR MASKTOCERRSTRUCTURE
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

if ~exist('planC')
    global planC
end

indexS = planC{end};
planC{indexS.dose}(end+1) = planC{indexS.dose}(doseNum);
planC{indexS.dose}(end).fractionGroupID = doseName;
planC{indexS.dose}(end).doseUID = createUID('dose');

planC{indexS.dose}(end).coord1OFFirstPoint = planC{indexS.dose}(doseNum).coord1OFFirstPoint + xyzT(1);
planC{indexS.dose}(end).coord2OFFirstPoint = planC{indexS.dose}(doseNum).coord2OFFirstPoint + xyzT(2);
planC{indexS.dose}(end).zValues = [planC{indexS.dose}(doseNum).zValues] + xyzT(3);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%============================================================================
% Copyright (C) 2014, Heikki Hyyti
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the "Software"),
% to deal in the Software without restriction, including without limitation
% the rights to use, copy, modify, merge, publish, distribute, sublicense,
% and/or sell copies of the Software, and to permit persons to whom the
% Software is furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
% FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
% DEALINGS IN THE SOFTWARE.
%============================================================================

function [vals, lags] = rootMeanSquaredErrors(data, ref, n_lags)
    lags = -n_lags:1:n_lags;
    vals = zeros(size(lags));
    base_start = 1 + n_lags;
    base_stop = length(data) - n_lags;
    for i = 1:length(lags)
        start_idx = base_start + lags(i);
        stop_idx = base_stop + lags(i);
        vals(i) = sqrt(mean((data(start_idx:stop_idx) - ref(base_start:base_stop)).^2));
    end
end

author: hhyyti

repo: dcm-imu

save_path: github-repos/MATLAB/hhyyti-dcm-imu


text:
function [trial_indx, subject_indx, datavar_matrix, datavar_cell] = select_trials_and_subjects(DAT, conditioning_text_var_name, text_pattern, varargin)
% Filter an event-level variable:
% - Return a logical vector of which events and subjects have values on a
% string variable (conditioning_text_var_name) that match a string pattern
% (text_pattern)
%
% - If the name of a second variable is entered as an optional argument,
% extract data from that variable, selecting only events that match the
% filter. [var name to extract data] in the command string below is the
% optional argument.
%
% conditioning_text_var_name is the name of an Event_Level variable in your dataset
% text_pattern is a string pattern to match with regexp
%
% [trial_indx, subject_indx, datavar_matrix, datavar_cell] = select_trials_and_subjects(DAT, conditioning_text_var, text_pattern, [var name to extract data])
%
% datavar_matrix is matrix of subjects x trials with extracted data
% datavar_cell is cell array of length subjects containing column vectors with extracted data

[datavar_matrix, datavar_cell] = deal([]);

[~, block] = get_var(DAT, conditioning_text_var_name);

% This function turns a string into a vector of which trials match that string.
subj_indx = @(subjcell) ~cellfun(@isempty, regexp(subjcell, text_pattern, 'match'));

% This function tests for all empty or nans in a cell.
nanfun = @(block) all(cellfun(@(x) all(isempty(x) | isnan(x)), block{1}));
%nanfun(block(1))

% This function tests whether any trials match the string
anyfun = @(x) any(subj_indx(x{1}));

% This function operates on one cell, for use with cellfun wrapper
onesubjfun = @(x) subj_indx(x{1});

% subj_indx(block{5})
% onesubjfun(block(5))

% doesn't work for unknown reasons - actually, can't handle nans.
% indic = cellfun(onesubjfun, block, 'UniformOutput', false)
% indic = cellfun(@(x) any(x), block)

clear trial_indx subject_indx
for i = 1:length(block)

    if nanfun(block(i))
        trial_indx{i} = false(size(block{i}));
        subject_indx(i, 1) = false;
    else
        trial_indx{i} = subj_indx(block{i});
        subject_indx(i, 1) = anyfun(block(i));

    end

end

% retrieve target data variable, removing NaNs and concatenating
if length(varargin) > 0

    var_name_to_extract = varargin{1};

    [~, datavar_cell, ~, ~, ~, textflag] = get_var(DAT, var_name_to_extract);

    for i = 1:length(trial_indx)

        datavar_cell{i} = datavar_cell{i}(trial_indx{i});

    end

    % pad to longest and concatenate into matrix

    %[~, textflag] = get_varlevel(DAT, var_name_to_extract);

    d = pad_cells_with_nan(datavar_cell, textflag);

    datavar_matrix = cat(2, d{:})';

end

end % function

function d = pad_cells_with_nan(d, textflag)

slen = max(cellfun(@length, d)); % max length for any subject

refvec = ones(slen, 1);

for i = 1:length(d)

    if textflag % text var
        reftext = cell(1, abs(length(d{i}) - length(refvec)))';
        [reftext{:}] = deal('');

        d{i} = [d{i}; reftext];

    else % numeric
        d{i} = padwithnan(d{i}, refvec, 1);
    end

end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% test for steerable transform

n = 512;
name = 'turbulence';
name = 'barb';
name = 'disk';
name = 'lena';
M = load_image(name,n);

k = 4;
options.nb_orientations = k;
J = 4;
Jmax = log2(n)-1;
Jmin = Jmax-J+1;
MW = perform_steerable_transform(M,Jmin,options);
M2 = perform_steerable_transform(MW,Jmin,options);

save_image = 1;
rep = ['results/steerable/'];

% reconstruction error
disp(['--> Reconstruction error: ' num2str(psnr(M,M2)) 'dB']);

if ~exist(rep)
    mkdir(rep);
end

% display and save the images
m = 0;
clf;
warning off;
for j=1:J
    for s=1:k
        m = m+1;
        A = MW{m+1};
        % A = rescale_wavelet_coefs(M,eta);
        str = ['j' num2str(j) '-s' num2str(s)];
        imageplot(A, str, k,J,m)
        if save_image
            imwrite(rescale(MW{m+1}), [rep name '-' str '.png'], 'png');
%            imwrite(rescale(MW{m+1}), [rep name '-' str '.jpg'], 'jpg');
        end
    end
end
colormap gray(256);

if save_image
    imwrite(rescale(M), [rep name '-original.png'], 'png');
    imwrite(rescale(MW{1}), [rep name '-high.png'], 'png');
    imwrite(rescale(MW{end}), [rep name '-low.png'], 'png');
end
warning on;

author: gpeyre

repo: matlab-toolboxes


text:
%CHERNOFFM Suboptimal discrimination linear mapping (Chernoff mapping)
%
%	W = CHERNOFFM(A,N,R)
%
% INPUT
%	A   Dataset
%	N   Number of dimensions to map to, N < C, where C is the number of classes
%	    (default: min(C,K)-1, where K is the number of features in A)
%	R   Regularization variable, 0 <= r <= 1, default is r = 0, for r = 1 the
%	    Chernoff mapping is (should be) equal to the Fisher mapping
%
% OUTPUT
% 	W   Chernoff mapping
%
% DESCRIPTION
% Finds a mapping of the labeled dataset A onto an N-dimensional linear
% subspace such that it maximizes the heteroscedastic Chernoff criterion
% (also called the Chernoff mapping).
%
% REFERENCE
% M. Loog and R.P.W. Duin, Linear Dimensionality Reduction via a Heteroscedastic
% Extension of LDA: The Chernoff Criterion, IEEE Transactions on pattern
% analysis and machine intelligence, vol. PAMI-26, no. 6, 2004, 732-739.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, FISHERM, NLFISHERM, KLM, PCA

% Copyright: M. Loog, marco@isi.uu.nl
% Image Sciences Institute, University Medical Center Utrecht
% P.O. Box 85500, 3508 GA Utrecht, The Netherlands

% $Id: chernoffm.m,v 1.3 2010/02/08 15:31:48 duin Exp $

function W = chernoffm(a,n,r)

		if (nargin < 3)
        r = 0;
	end

	if (nargin < 2)
		prwarning (4, 'number of dimensions to map to not specified, assuming min(C,K)-1');
		n = [];
	end

	% If no arguments are specified, return an untrained mapping.

	if (nargin < 1) | (isempty(a))
		W = prmapping('chernoffm',{n,r});
		W = setname(W,'Chernoff mapping');
		return
	end

	isvaldset(a,2,2); % at least 2 objects per class, 2 classes

	% If N is not given, set it.

	[m,k,c] = getsize(a); if (isempty(n)), n = min(k,c)-1; end
	if (n >= m)
		error('The dataset is too small for the requested dimensionality')
	end

	% To simplify computations, the within scatter W is set to the identity matrix
    % therefore A is centered and sphered by KLMS.

	v = klms(a); a = a*v;
	k = size(v,2); % dimensionalities may have changed for small training sets

	% Now determine a solution to the Chernoff criterion

	[U,G] = meancov(a);

	CHERNOFF = zeros(k);
    p = getprior(a);

    for i = 1:c-1 % loop over all class pairs
        for j = i+1:c
            p_i = p(i)/(p(i)+p(j)); p_j = p(j)/(p(i)+p(j));
            G_i = (1-r)*G(:,:,i) + r*eye(k); G_j = (1-r)*G(:,:,j) + r*eye(k);
            G_ij = p_i*G_i + p_j*G_j;
            m_ij = sqrtm(real(prinv(G_ij)))*(U(i,:) - U(j,:))';
            CHERNOFF = CHERNOFF + p(i) * p(j) * ...
                ( m_ij * m_ij' ...
                + 1/(p_i*p_j) * ( logm(G_ij) ...
                - p_i * logm(G_i) - p_j * logm(G_j) ) );
        end
    end

	[F,V] = preig(CHERNOFF);
	[dummy,I] = sort(-diag(V));
	I = I(1:n);
	rot = F(:,I);
	off = -mean(a*F(:,I));
	preW = affine(rot,off,a);
	W = v*preW;					% Construct final mapping.
	W = setname(W,'Chernoff mapping');

return

author: marianux

repo: ecg-kit

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


text:
function s = ySegment(x,y1,y2,z)

% YSEGMENT   Horizontal segment in the Y direction
%   YSEGMENT(X,Y1,Y2,Z) makes an horizontal segment at the plane
%   location (X,Z) from widths Y1 to Y2.

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

s = [x;y1;z;x;y2;z];

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [mwcn,wpixc, globmax]=analyzeWconfPSF(w,peval,p)
% [mwcn,wpixc, globmax]=analyzeWconfPSF(w,peval,p)
% Finds local maxima in the results w (nx*ny,ncomp) convoluted with PSF
% (generated from parameters p.lambda, p.NA, p.ri,p.pixelsize)
% wpixc - w convolved with estimated PSF.
% mwcn = image of local maxima (non zero pixels) with values of hte convoluted image, nornalised to the maximum of each column.
% globmax = maximum of each column of W convolved with PSF. - PSF and W are L1 normalised - teh maximum can be used as a "quality measure."
%
% To get number of local maxima with relative strength < .5: sum(mwcn>.5)
ncomp=size(w,2);
o = kSimPSF( {'lambdaEm',p.lambda;'Pi4Em',0;'relEmInt',1;'relEmPhase',0;'na',p.NA;'ri',p.ri;'sX',peval.nx;'sY',peval.ny;'sZ',1;'scaleX',p.pixelsize;'scaleY',p.pixelsize;'twophoton',0;'confocal',0;'nonorm',0;'pinhole',1;'Pi4Ex',0;'relExInt',1});
psf=double(o); % it is normalized to 1
autoconvmax=max(max(conv2(psf,psf,'same'))); % maximum of the "autoconvolution" - use to normalise the convolution with w.
wpix=reshape(w,peval.nx, peval.ny, ncomp); % each frame normalized to 1
wpixc=(convstack(wpix,psf,'same'));

wc=reshape(wpixc,peval.nx*peval.ny,ncomp);
mpix=maximastack(wpixc); % binary image of local maxima locations in each frame
m=reshape(mpix,peval.nx*peval.ny,ncomp);
mwc=m.*wc; % pixels indicates locations and value of hte local maxima
globmax=max(mwc)/autoconvmax; % The value of the global maximum for each column of W normalised with respect to the 'autoconvolution' of the PSF.
mwcn=normcMax(mwc); % normalised to the maximum of each column.

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function m1 = placematrix(m1,m2,pos)

% function m1 = placematrix(m1,m2,pos)
%
% <m1> is a 2D matrix, with potentially some extra stuff in the third dimension.
% <m2> is a 2D matrix, with potentially some extra stuff in the third dimension.
%   note that size(m1,3) must equal size(m2,3).
% <pos> (optional) is [R C] with a position.  R and C can be any integers.
%   special case is [] which means to center <m2> with respect to <m1>.
%   if exact centering can't be achieved, we shift down and right.
%   default: [].
%
% place <m2> in <m1> positioned with first element at <pos>.
% if any part of <m2> lies outside of <m1>, it just gets ignored.
%
% example:
% isequal(placematrix([1 2 3; 4 5 6; 7 8 9],[10 10; 10 10],[0 0]),[10 2 3; 4 5 6; 7 8 9])

%% SEE ALSO padarray.m ?
%% see also assignsubvolume2.m???

% input
if ~exist('pos','var') || isempty(pos)
  pos = [];
end

m1r = size(m1,1);
m1c = size(m1,2);
m2r = size(m2,1);
m2c = size(m2,2);

if isempty(pos)
  pos = [1+ceil((m1r-m2r)/2) 1+ceil((m1c-m2c)/2)];
end

badup = max(0,1-pos(1));               % how many bad pixels to the up
badleft = max(0,1-pos(2));             % how many bad pixels to the left
baddown = max(0,(pos(1)+m2r-1)-m1r);   % how many bad pixels to the bottom
badright = max(0,(pos(2)+m2c-1)-m1c);  % how many bad pixels to the right

m1(pos(1)+badup:pos(1)+m2r-1-baddown,pos(2)+badleft:pos(2)+m2c-1-badright,:) = ...
  m2(1+badup:end-baddown,1+badleft:end-badright,:);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
Wd = function (matvec)
%take a 12 param matvec, typically WARPDRIVE_MATVEC_INV_000000 and
%change it to a WARP_DATA
%Used to test transformations needed to make AFNI use 3dWarpDrive's output
%as a TLRC transform

if (nargin == 0) matvec = [];

if (isempty(matvec)),
   WARNING: Using preset matvecs!\n\n');
   %test case, matvecs from DemoSubj_SurfVol_Alnd_Exp_at2+tlrc
   if (0),
      matvec = [  0.9709206    -0.04199445     0.01165925     -0.7014965
               0.07297896    0.9860933    -0.08915191     -0.2648731
               0.005342953    -0.02306779  0.9538299      -6.124281
               ];
   else
      %that would be WARPDRIVE_MATVEC_INV_000000
      %the proper one to use for WARP_DATA, the one that takes you
      %from +orig to +tlrc
      matvec = [  1.026725     0.04352641   -0.008482005      0.6798272
                  -0.07667368  1.013075     0.09562658      0.8001939
                  -0.007605589   0.02425676 1.050765       6.436271
               ];
   end
else
   if (length(matvec) ~= 12),
      fprintf(2,'Error: matvec must be 12 elements long\n');
   end
end

Mfor = [matvec; 0 0 0 1]; %the forward transform
Mbac = inv(Mfor); %the inverse transform

%begin transformation
mfor = Mfor(1:3, 1:3);
mbac = Mbac(1:3, 1:3);
bvec = -Mfor(1:3,4); %Don't ask about the minus sign, see README.attributes:
                     %Under WARP_DATA, The forward transformation is  [x_map] = [mfor] [x_in]  - [bvec];
                     %Under Registration Attributes: [xyz_out] = [mat] ([xyz_in]-[xyz_cen]) + [vec] + [xyz_cen]
svec = -Mbac(1:3,4);
bot = [-80 ; -80 ; -65]; %good enough for TT box
top = [80 ; 110 ; 85];

Wd = [reshape(mfor',1,9) reshape(mbac',1,9) reshape(bvec,1,3) reshape(svec,1,3) [bot'] [top']];

fprintf(1,'Wd = \n\n');
fprintf(1,'       %g   %g   %g   %g   %g   %g\n',Wd)
fprintf(1,'\n\n');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% calculate AUC

for j=1:15
[~,~,~,AUC(j)] = perfcurve(ans_auc(:,j),prop(:,j),0);
end

mean_AUC = mean(AUC);
disp(mean_AUC)

author: PatternRecognition

repo: OpenBMI


text:
function [desc, f] = at_cnnfeat2vlfeat(x)

[u,v] = meshgrid(1:size(x,2),1:size(x,1));
f = [u(:)'; v(:)'];
desc1 = reshape(shiftdim(x(:)),size(x,1)*size(x,2),[])';

desc = yael_vecs_normalize(desc1);

author: HajimeTaira


% sin_x.m: This m-file calculates and plots the
% function sin(x) for 0 <= x <= 6.
x = 0:0.1:6;
y = sin(x);
plot(x,y);


text:
function h = rdivide(f, g, pref)
%./   Pointwise CHEBFUN right divide.
%   F./G returns a CHEBFUN that represents the function F(x)/G(x).
%   If F and G are array-valued column (row) CHEBFUNs, they must have the same
%   number of columns (rows).
%
% See also MRDIVIDE, TIMES.

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

% Trivial empty case:
if ( isempty(f) || isempty(g) )
    h = chebfun();
    return
end

% Trivial zero numerator case:
if ( isnumeric(f) && ~any(f) )
    h = 0*g;
    return
end

if ( nargin < 3 )
    pref = chebfunpref();
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CHEBFUN ./ CHEBFUN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ( isa(f,'chebfun') && isa(g, 'chebfun') )

    % Check the number of columns match:
    if ( numColumns(f) ~= numColumns(g) )
        if ( numColumns(f) == 1 )
            h = g;
            for k = 1:numColumns(g)
                h(k) = rdivide(f, g(k));
            end
        elseif ( numColumns(g) == 1 )
            h = f;
            for k = 1:numColumns(f)
                h(k) = rdivide(f(k), g);
            end
        else
            error('CHEBFUN:CHEBFUN:rdivide:quasi', ...
                'Chebfun quasimatrix dimensions must agree.')
        end
        return
    end

    if ( numel(f) == 1 && numel(g) == 1 )
        % CHEBFUN case :

        % If one of the two CHEBFUNs uses a PERIODICTECH reprensetation,
        % cast it to a NONPERIODICTECH.
        if ( ~isPeriodicTech(f.funs{1}) && isPeriodicTech(g.funs{1}) )
            g = chebfun(g, g.domain, 'tech', get(f.funs{1}, 'tech'));
        elseif ( isPeriodicTech(f.funs{1}) && ~isPeriodicTech(g.funs{1}) )
            f = chebfun(f, f.domain, 'tech', get(g.funs{1}, 'tech'));
        end

        % Array-valued CHEBFUN case:
        h = columnRdivide(f, g, pref);
    else
        % QUASIMATRIX case:

        % Convert to a cell array:
        f = mat2cell(f);
        g = mat2cell(g);
        % Loop over the columns:
        for k = numel(f):-1:1
            h(k) = columnRdivide(f{k}, g{k}, pref);
        end
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CHEBFUN ./ constant %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ( isa(f, 'chebfun') )

    numColsF = numColumns(f);
    numelF = numel(f);
    numelG = numel(g);

    % Different cases:
    if ( numColsF == 1 )
        % e.g., x./[1 2 3]
        h = columnRdivide(f, g, pref);
    elseif ( numelG == 1 )
        % e.g., [x sin(x)]./2
        for k = numelF:-1:1
            h(k) = columnRdivide(f(k), g, pref);
        end
    elseif ( numelG == numColsF )
        % e.g., [x sin(x) exp(x)]./[1 2 3]
        if ( numel(f) == 1 )
            h = columnRdivide(f, g, pref);
        else
            f = mat2cell(f);
            for k = numColsF:-1:1
                h(k) = columnRdivide(f{k}, g(k), pref);
            end
        end
    else
        error('CHEBFUN:CHEBFUN:rdivide:dim', ...
            'Chebfun quasimatrix dimensions must agree.');
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%% constant ./ CHEBFUN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else

    numColsG = numColumns(g);
    numelF = numel(f);
    numelG = numel(g);

    % Different cases:
    if ( numColsG == 1 )
        % e.g., [1 2 3]./(2+x)
        for k = numelF:-1:1
            h(k) = columnRdivide(f(k), g, pref);
        end
        h = quasi2cheb(h);
    elseif ( numelF == 1 )
        % e.g., 2./(2+[x x abs(x)])
        g = mat2cell(g);
        for k = numColsG:-1:1
            h(k) = columnRdivide(f, g{k}, pref);
        end
        try h = quasi2cheb(h); catch, end
    elseif ( numelF == numColsG )
        % e.g., [1 2]./(2+[x sin(x)])
        g = mat2cell(g);
        for k = numColsG:-1:1
            h(k) = columnRdivide(f(k), g{k}, pref);
        end
        try h = quasi2cheb(h); catch, end
    else
        error('CHEBFUN:CHEBFUN:rdivide:dim', ...
            'Chebfun quasimatrix dimensions must agree.');
    end

end

end

function h = columnRdivide(f, g, pref)

% If g is numeric then call TIMES():
if ( isnumeric(g) )
    if ( any(g == 0) )
        % note g could be a vector, with some zero components
        % TODO:  Return identically Inf/NaN CHEBFUN instead?
        error('CHEBFUN:CHEBFUN:rdivide:columnRdivide:divisionByZero', ...
            'Division by zero.')
    end
    h = f.*(1./g);
    return
end

% Check for zero FUNs:
for k = 1:numel(g.funs)
    if ( iszero(g.funs{k}) )
        % TODO:  Return CHEBFUN with identically Inf/NaN FUN instead?
        error('CHEBFUN:CHEBFUN:rdivide:columnRdivide:divisionByZeroChebfun', ...
            'Division by CHEBFUN with identically zero FUN.');
    end
end

% Add breaks at the roots of g:
g = addBreaksAtRoots(g, pref);

if ( isa(f, 'chebfun') )

    % Check that the domains are the same:
    if ( ~domainCheck(f, g) )
        error('CHEBFUN:CHEBFUN:rdivide:columnRdivide:domain', ...
            'Inconsistent domains.');
    end

    % Check that the orientation is the same:
    if ( xor(f.isTransposed, g.isTransposed) )
        error('CHEBFUN:CHEBFUN:rdivide:columnRdivide:dim', ...
            'Matrix dimension do not agree (transposed)');
    end

    % Introduce matching breakpoints in f and g:
    [f, g] = overlap(f, g);

    % Copy g to h in preparation for output:
    h = g;

    % Loop over the FUNS:
    for k = 1:numel(g.funs)
        h.funs{k} = rdivide(f.funs{k}, g.funs{k});
    end

    % Divide the pointValues:
    h.pointValues = f.pointValues./g.pointValues;

    % Avoid NaN or Inf from divide by zero where possible:
    hpv = h.pointValues; fpv = f.pointValues; gpv = g.pointValues;
    idx = (isnan(hpv) & ~(isnan(fpv) | isnan(gpv))) | ...
        (isinf(hpv) & (abs(fpv) < 100*chebfuneps*vscale(f)));
    if ( any(idx) )
        pvavg = (feval(h, h.domain, 'left') + feval(h, h.domain, 'right'))/2;
        h.pointValues(idx) = pvavg(idx);
    end

else

    % Copy g to h in preparation for output:
    h = g;

    % Loop over the FUNS:
    for k = 1:numel(g.funs)
        h.funs{k} = rdivide(f, g.funs{k});
    end

    % Divide the pointValues:
    h.pointValues = f./g.pointValues;

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function answer = strstart (str1, str2)
    case_sensitive = true;
    answer = strstarti (str1, str2, case_sensitive);
end

%!test
%! myassert(strstart('haha', 'haha/private'));
%! myassert(~strstart('Haha', 'haha/private'));
%! myassert(~strstart('hahablah', 'haha/private2'));
%! myassert(~strstart('packed:rcond:complexNotSupported', ...
%!     'packed:chol:complexNotSupported'));
%! out = strstart('haha', {'haha/private', 'blah/private', ''});
%! myassert(out, [true, false, false])

author: Sable

repo: mcbench-benchmarks


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

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

% Last Modified by GUIDE v2.5 04-Jun-2013 16:23:04

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

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

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

% Choose default command line output for JAABAInitOpen
handles.output = hObject;
handles.sel = struct('val',[],'edit',false);

figureJLabel = varargin{1};
% Change a few things so they still work well on Mac
adjustColorsIfMac(hObject);

centerOnParentFigure(hObject,figureJLabel);

% Update handles structure
guidata(hObject, handles);

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

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

% Get default command line output from handles structure
varargout{1} = handles.sel;
delete(handles.figure1);

% --- Executes on button press in pushbutton_new.
function pushbutton_new_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_new (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.sel.val = 'New';
guidata(hObject,handles);
figure1_CloseRequestFcn(handles.figure1,[],handles);

% --- Executes on button press in pushbutton_open.
function pushbutton_open_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_open (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.sel.val = 'Open';
guidata(hObject,handles);
figure1_CloseRequestFcn(handles.figure1,[],handles);

% --- Executes on button press in pushbutton_opengt.
function pushbutton_opengt_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_opengt (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.sel.val = 'OpenGT';
guidata(hObject,handles);
figure1_CloseRequestFcn(handles.figure1,[],handles);

% --- Executes on button press in checkbox_edit.
function checkbox_edit_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox_edit

handles.sel.edit = get(hObject,'Value');
guidata(hObject,handles);

% --- Executes when user attempts to close figure1.
function figure1_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: delete(hObject) closes the figure
% delete(hObject);
uiresume(hObject);

author: kristinbranson

repo: JAABA


text:
function [imgs, interpolated] = scaleup_Zeyde(conf, imgs)

% Super-Resolution Iteration
for j = 1:conf.level
    fprintf('Scale-Up Zeyde et al. #%d', j);
    midres = resize(imgs, conf.upsample_factor, conf.interpolate_kernel);
    interpolated = resize(imgs, conf.scale, conf.interpolate_kernel);

    for i = 1:numel(midres)
        features = collect(conf, {midres{i}}, conf.upsample_factor, conf.filters);
        features = double(features);
        % Encode features using OMP algorithm

        coeffs = omp(double(conf.dict_lores), conf.V_pca' * features, [], 3);

        % Reconstruct using patches' dictionary
        patches = conf.dict_hires * full(coeffs);

        % Add low frequencies to each reconstructed patch
        patches = patches + collect(conf, {interpolated{i}}, conf.scale, {});

        % Combine all patches into one image
        img_size = size(imgs{i}) * conf.scale;
        grid = sampling_grid(img_size, ...
            conf.window, conf.overlap, conf.border, conf.scale);
        result = overlap_add(patches, img_size, grid);
        imgs{i} = result; % for the next iteration
        fprintf('.');
    end
end
fprintf('\n');

author: thomas-koehler

repo: SupER


text:
function [ mO ] = ApplySideWindowFiltering( mI, boxRadius, numIterations )
% ----------------------------------------------------------------------------------------------- %
% [ mO ] = ApplySideWindowFiltering( mI, boxRadius, numIterations )
%   Applying Image Edge Preserving Filter using the Side Window Box
%   Filtering algorithm.
% Input:
%   - mI            -   Input Image.
%                       Structure: Image Matrix (Single Channel)
%                       Type: 'Single' / 'Double'.
%                       Range: [0, 1].
%   - boxRadius     -   Box Radius.
%                       Sets the radius of the box filter.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: {1, 2, 3, ...}.
%   - numIterations -   Number of Iterations.
%                       Sets the number of iterations to apply the filter.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: {1, 2, 3, ...}.
% Output:
%   - mO            -   Output Image.
%                       Structure: Image Matrix (Single Channel)
%                       Type: 'Single' / 'Double'.
%                       Range: [0, 1].
% References
%   1.  Side Window Filtering (https://arxiv.org/abs/1905.07177).
% Remarks:
%   1.  Can be much faster if written In Place.
% TODO:
%   1.  U.
% Release Notes:
%   -   1.0.000     24/04/2021  Royi Avital     RoyiAvital@yahoo.com
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

dataType    = class(mI);
boxLength   = (2 * boxRadius) + 1;

% Atoms of the Filter Bank (Since any of the filters is separable)
vK  = ones(boxLength, 1, dataType) / boxLength;
vKU = ones(boxLength, 1, dataType) / (boxRadius + 1); %<! Up
vKD = ones(boxLength, 1, dataType) / (boxRadius + 1); %<! Down

vKU((boxRadius + 2):boxLength)  = 0;
vKD(1:boxRadius)                = 0;

mO = padarray(mI, [boxRadius, boxRadius], 'both', 'replicate');
mF = zeros(size(mO, 1), size(mO, 2), 8, dataType);
mM = zeros(size(mO), dataType); %<! Minimum Index

for kk = 1:numIterations
    % Written for clarity, not performance
    mF(:, :, 1) = conv2(vK, vKU, mO, 'same'); %<! Left Box
    mF(:, :, 2) = conv2(vK, vKD, mO, 'same'); %<! Right Vox
    mF(:, :, 3) = conv2(vKU, vK, mO, 'same'); %<! Up Box
    mF(:, :, 4) = conv2(vKD, vK, mO, 'same'); %<! Down Box
    mF(:, :, 5) = conv2(vKU, vKU, mO, 'same'); %<! NW Box
    mF(:, :, 6) = conv2(vKU, vKD, mO, 'same'); %<! NE Box
    mF(:, :, 7) = conv2(vKD, vKU, mO, 'same'); %<! SW Box
    mF(:, :, 8) = conv2(vKD, vKD, mO, 'same'); %<! SE Box

    [~, mM(:, :)]   = min(abs(mF - mO), [], 3); %<! Index which minimizes
    for jj = 1:size(mO, 2)
        for ii = 1:size(mO, 1)
            mO(ii, jj) = mF(ii, jj, mM(ii, jj));
        end
    end

end

% Removing boundary
mO = mO((boxRadius + 1):(end - boxRadius), (boxRadius + 1):(end - boxRadius));
mO = min(max(mO, 0), 1);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function scaledImg = scaleImg(img,range)
% scale image to range
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

if(nargin < 2)
    range = [0 1];
end

if(length(range) == 1)
    range = [0 range];
end

scaledImg = ((img - min(img(:))) * (range(2)-range(1)))./(max(img(:)) - min(img(:)));

end

author: thomaskuestner


text:
function [rot,zoom,frustum,origin] = mrmGetRotation(id)
% Get rotation and zoom
%
%   [rot,zoom,frustum] = mrmGetRotation(id)
%

p.actor=0; p.get_all=1;
[id,stat,r] = mrMesh('localhost', id, 'get', p);
zoom = diag(chol(r.rotation'*r.rotation))';
rotMat = r.rotation/diag(zoom);
frustum = r.frustum;
origin = r.origin;
% Note- there may be slight rounding errors allowing the inputs to
% asin/atan go outside of the range (-1,1). May want to clip those.
rot(2) = asin(rotMat(1,3));
if (abs(rot(2))-pi/2).^2 < 1e-9,
    rot(1) = 0;
    rot(3) = atan2(-rotMat(2,1), -rotMat(3,1)/rotMat(1,3));
else
    c      = cos(rot(2));
    rot(1) = atan2(rotMat(2,3)/c, rotMat(3,3)/c);
    rot(3) = atan2(rotMat(1,2)/c, rotMat(1,1)/c);
end
rot(1) = -rot(1); % flipped OpenGL Y-axis.
%fprintf('rot=[%0.6f %0.6f %0.6f];\nzoom=[%0.3f %0.3f %0.3f];\nfrustum=[%0.6f %0.6f %0.6f %0.6f];\n',rot,zoom,r.frustum);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function planC = quality_assure_planC(fileName, planC, forceSaveFlag)
%function planC = quality_assure_planC(fileName, planC, forceSaveFlag)
%
% This function quality assures planC.
%
%APA, 04/07/2011

if ~exist('planC','var')
    global planC
end
if ~exist('forceSaveFlag','var')
    forceSaveFlag = 0;
end

global stateS
indexS = planC{end};

% Quality Assurance
bug_found = 0;

% Detect and Fix incorrect rasterSegments
if length(planC{indexS.header})>0 && isfield(planC{indexS.header}(1), 'lastSavedInVer')
    lastSavedInVer = planC{indexS.header}(1).lastSavedInVer;
else
    lastSavedInVer = '';
end

if ~isempty(str2double(lastSavedInVer)) && (isempty(lastSavedInVer) || str2double(lastSavedInVer) < 4.1) && ~isempty(planC{indexS.structures})
    planC = reRasterAndUniformize(planC);
    bug_found = 1;
end

%Update imageType for CT scans where required
for nScan = 1:length(planC{indexS.scan})
    if strcmpi(planC{indexS.scan}(nScan).scanInfo(1).imageType,'CT')
        [planC{indexS.scan}(nScan).scanInfo(:).imageType] = deal('CT SCAN');
    end
end

% %Check for mesh representation and load meshes into memory
% currDir = cd;
% meshDir = fileparts(which('libMeshContour.dll'));
% cd(meshDir)
% for strNum = 1:length(planC{indexS.structures})
%     if isfield(planC{indexS.structures}(strNum),'meshRep') && ~isempty(planC{indexS.structures}(strNum).meshRep) && planC{indexS.structures}(strNum).meshRep
%         try
%             calllib('libMeshContour','loadSurface',planC{indexS.structures}(strNum).strUID,planC{indexS.structures}(strNum).meshS)
%         catch
%             planC{indexS.structures}(strNum).meshRep    = 0;
%             planC{indexS.structures}(strNum).meshS      = [];
%         end
%     end
% end
% cd(currDir)

if ~isfield(stateS,'optS')
    stateS.optS = opts4Exe([getCERRPath,'CERROptions.json']);
end

%Check color assignment for displaying structures
[assocScanV,relStrNumV] = getStructureAssociatedScan(1:length(planC{indexS.structures}),planC);
for scanNum = 1:length(planC{indexS.scan})
    scanIndV = find(assocScanV==scanNum);
    for i = 1:length(scanIndV)
        strNum = scanIndV(i);
        colorNum = relStrNumV(strNum);
        if isempty(planC{indexS.structures}(strNum).structureColor)
            color = stateS.optS.colorOrder( mod(colorNum-1, size(stateS.optS.colorOrder,1))+1,:);
            planC{indexS.structures}(strNum).structureColor = color;
        end
    end
    % Check for scanType field and populate it with Series description
    if isempty(planC{indexS.scan}(scanNum).scanType) && ...
            isfield(planC{indexS.scan}(scanNum).scanInfo(1),'DICOMHeaders') && ...
            isfield(planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders,'SeriesDescription')
        planC{indexS.scan}(scanNum).scanType = planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders.SeriesDescription;
    end
    % Save SOPInstanceUID and SOPClassUID within scanInfo
    if ~isfield(planC{indexS.scan}(scanNum).scanInfo(1),'sopInstanceUID') || ...
            isempty(planC{indexS.scan}(scanNum).scanInfo(1).sopInstanceUID) || ...
            isempty(planC{indexS.scan}(scanNum).scanInfo(1).patientID) || ...
            isempty(planC{indexS.scan}(scanNum).scanInfo(1).frameOfReferenceUID)
        if isfield(planC{indexS.scan}(scanNum).scanInfo(1),'DICOMHeaders') ...
                && ~isempty(planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders)
            for slcNum = 1:length(planC{indexS.scan}(scanNum).scanInfo)
                planC{indexS.scan}(scanNum).scanInfo(slcNum).sopInstanceUID = ...
                    planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders.SOPInstanceUID;
                planC{indexS.scan}(scanNum).scanInfo(slcNum).sopClassUID = ...
                    planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders.SOPClassUID;
                planC{indexS.scan}(scanNum).scanInfo(slcNum).seriesInstanceUID = ...
                    planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders.SeriesInstanceUID;
                %                 planC{indexS.scan}(scanNum).scanInfo(slcNum).patientBirthDate = ...
                %                     planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders.PatientBirthDate;
                if isfield(planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders,'FrameofReferenceUID')
                    planC{indexS.scan}(scanNum).scanInfo(slcNum).frameOfReferenceUID = ...
                        planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders.FrameofReferenceUID;
                end
                if isfield(planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders,'PatientID')
                    planC{indexS.scan}(scanNum).scanInfo(slcNum).patientID = ...
                        planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders.PatientID;
                end
                if isfield(planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders,'StudyDate')
                    planC{indexS.scan}(scanNum).scanInfo(slcNum).studyDate = ...
                        planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders.StudyDate;
                end
                if isfield(planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders,'StudyTime')
                    planC{indexS.scan}(scanNum).scanInfo(slcNum).studyTime = ...
                        planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders.StudyTime;
                end

            end
        end
    end
    % Save StudyDate and StudyTime within scanInfo
    if ~isfield(planC{indexS.scan}(scanNum).scanInfo(1),'StudyDate') || ...
            isempty(planC{indexS.scan}(scanNum).scanInfo(1).StudyDate)
        if isfield(planC{indexS.scan}(scanNum).scanInfo(1),'DICOMHeaders') ...
                && ~isempty(planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders)
            for slcNum = 1:length(planC{indexS.scan}(scanNum).scanInfo)
                if isfield(planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders,'StudyDate')
                    planC{indexS.scan}(scanNum).scanInfo(slcNum).studyDate = ...
                        planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders.StudyDate;
                end
                if isfield(planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders,'StudyTime')
                    planC{indexS.scan}(scanNum).scanInfo(slcNum).studyTime = ...
                        planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders.StudyTime;
                end
            end
        end
    end
    % Save image orientation to scanInfo
    if ~isfield(planC{indexS.scan}(scanNum).scanInfo(1),'imageOrientationPatient') || ...
            isempty(planC{indexS.scan}(scanNum).scanInfo(1).imageOrientationPatient)
        if isfield(planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders,'ImageOrientationPatient')
            imgOriV = planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders.ImageOrientationPatient;
            [planC{indexS.scan}(scanNum).scanInfo(:).imageOrientationPatient] = deal(imgOriV);
        end
    end
    if ~isfield(planC{indexS.scan}(scanNum).scanInfo(1),'imagePositionPatient') || ...
            isempty(planC{indexS.scan}(scanNum).scanInfo(1).imagePositionPatient)
        if isfield(planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders,'ImagePositionPatient')
            for slc = 1:length(planC{indexS.scan}(scanNum).scanInfo)
                planC{indexS.scan}(scanNum).scanInfo(slc).imagePositionPatient = ...
                    planC{indexS.scan}(scanNum).scanInfo(slc).DICOMHeaders.ImagePositionPatient;
            end
        end
    end
end

%Check dose-grid
for doseNum = 1:length(planC{indexS.dose})
    if length(planC{indexS.dose}(doseNum).zValues) > 1
        if planC{indexS.dose}(doseNum).zValues(2) - planC{indexS.dose}(doseNum).zValues(1) < 0
            planC{indexS.dose}(doseNum).zValues = flipud(planC{indexS.dose}(doseNum).zValues);
            %planC{indexS.dose}(doseNum).doseArray = flipdim(planC{indexS.dose}(doseNum).doseArray,3);
            planC{indexS.dose}(doseNum).doseArray = flip(planC{indexS.dose}(doseNum).doseArray,3);
        end
    end
    % Save image orientation and position to dose field
    if ~isfield(planC{indexS.dose}(doseNum),'imageOrientationPatient') || ...
            isempty(planC{indexS.dose}(doseNum).imageOrientationPatient)
        if isfield(planC{indexS.dose}(doseNum).DICOMHeaders,'ImageOrientationPatient')
            imgOriV = planC{indexS.dose}(doseNum).DICOMHeaders.ImageOrientationPatient;
            planC{indexS.dose}(doseNum).imageOrientationPatient = imgOriV;
        end
    end
    if ~isfield(planC{indexS.dose}(doseNum),'imagePositionPatient') || ...
            isempty(planC{indexS.dose}(doseNum).imagePositionPatient)
        if isfield(planC{indexS.dose}(doseNum).DICOMHeaders,'ImagePositionPatient')
            imgPosV = planC{indexS.dose}(doseNum).DICOMHeaders.ImagePositionPatient;
            planC{indexS.dose}(doseNum).imagePositionPatient = imgPosV;
        end
    end
end

%Check whether uniformized data is in cellArray format.
if ~isempty(planC{indexS.structureArray}) && iscell(planC{indexS.structureArray}(1).indicesArray)
    planC = setUniformizedData(planC,planC{indexS.CERROptions});
    indexS = planC{end};
end

if length(planC{indexS.structureArrayMore}) ~= length(planC{indexS.structureArray})
    for saNum = 1:length(planC{indexS.structureArray})
        if saNum == 1
            planC{indexS.structureArrayMore} = struct('indicesArray', {[]},...
                'bitsArray', {[]},...
                'assocScanUID',{planC{indexS.structureArray}(saNum).assocScanUID},...
                'structureSetUID', {planC{indexS.structureArray}(saNum).structureSetUID});

        else
            planC{indexS.structureArrayMore}(saNum) = struct('indicesArray', {[]},...
                'bitsArray', {[]},...
                'assocScanUID',{planC{indexS.structureArray}(saNum).assocScanUID},...
                'structureSetUID', {planC{indexS.structureArray}(saNum).structureSetUID});
        end
    end
end

% Get CERR version of last save
if length(planC{indexS.header})>0 &&...
        isfield(planC{indexS.header},'lastSavedInVer')...
        && ~isempty(planC{indexS.header}.lastSavedInVer)
    CERRImportVersion = planC{indexS.header}.lastSavedInVer;
elseif length(planC{indexS.header})>0 && ...
        isfield(planC{indexS.header},'CERRImportVersion')...
        && ~isempty(planC{indexS.header}.CERRImportVersion)
    CERRImportVersion = planC{indexS.header}.CERRImportVersion;
else
    CERRImportVersion = '0';
end

%Check DSH Points for old CERR versions
if str2num(CERRImportVersion(1)) < 4
    bug_found = 1;
    for structNum = 1:length(planC{indexS.structures})
        if ~isempty(planC{indexS.structures}(structNum).DSHPoints)
            planC = getDSHPoints(planC, stateS.optS, structNum);
        end
    end
end

% Fix RTOG orientation for HFP scans
if str2num(strtok(CERRImportVersion, ',')) < 5.2
    for scanNum = 1:length(planC{indexS.scan})
        if isfield(planC{indexS.scan}(scanNum).scanInfo(1),'DICOMHeaders') ...
                && ~isempty(planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders) ...
                && isfield(planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders,'PatientPosition')
            %pPos = planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders.PatientPosition;
            imgOriV = planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders.ImageOrientationPatient;
        else
            %pPos = '';
            imgOriV = [];
        end
        if  ~isempty(imgOriV) && max(abs((imgOriV(:) - [-1 0 0 0 -1 0]'))) < 1e-3 % Position: HFP
            planC = flipAlongX(scanNum, planC);
            bug_found = 1;
        end
    end
end

% Add radiopharma params to scanInfo
for scanNum = 1:length(planC{indexS.scan})
    if isfield(planC{indexS.scan}(scanNum).scanInfo(1),'DICOMHeaders') ...
            && ~isempty(planC{indexS.scan}(scanNum).scanInfo(1).DICOMHeaders) ...
            && any(strcmpi(planC{indexS.scan}(scanNum).scanInfo(1).imageType,{'PT','PET'})) ...
            && isempty(planC{indexS.scan}(scanNum).scanInfo(1).decayCorrection)
        for slcNum = 1:size(planC{indexS.scan}(scanNum).scanArray,3)
            dicomhd = planC{indexS.scan}(scanNum).scanInfo(slcNum).DICOMHeaders;
            ptweight = [];
            if isfield(dicomhd,'PatientWeight')
                ptweight = dicomhd.PatientWeight;
            elseif isfield(dicomhd,'PatientsWeight')
                ptweight = dicomhd.PatientsWeight;
            end
            if isempty(dicomhd.RadiopharmaceuticalInformationSequence)
                continue;
            end
            injectionTime = ...
                dicomhd.RadiopharmaceuticalInformationSequence.Item_1.RadiopharmaceuticalStartTime;
            halfLife = ...
                dicomhd.RadiopharmaceuticalInformationSequence.Item_1.RadionuclideHalfLife;
            injectedDose = ...
                dicomhd.RadiopharmaceuticalInformationSequence.Item_1.RadionuclideTotalDose;
            seriesTime = dicomhd.SeriesTime;
            decayCorrection = dicomhd.DecayCorrection;
            correctedImage = dicomhd.CorrectedImage;
            imageUnits = dicomhd.Units;
            petActivityConcentrationScaleFactor = [];
            if isfield(dicomhd,'ActivityConcentrationScaleFactor')
                petActivityConcentrationScaleFactor = ...
                    dicomhd.ActivityConcentrationScaleFactor;
            end
            patientSize = [];
            if isfield(dicomhd,'PatientSize')
                patientSize = dicomhd.PatientSize;
            end

            planC{indexS.scan}(scanNum).scanInfo(slcNum).patientWeight = ...
                ptweight;
            planC{indexS.scan}(scanNum).scanInfo(slcNum).patientSize = ...
                patientSize;
            planC{indexS.scan}(scanNum).scanInfo(slcNum).petActivityConcentrationScaleFactor = ...
                petActivityConcentrationScaleFactor;
            planC{indexS.scan}(scanNum).scanInfo(slcNum).acquisitionTime = ...
                dicomhd.AcquisitionTime;
            planC{indexS.scan}(scanNum).scanInfo(slcNum).injectionTime = ...
                injectionTime;
            planC{indexS.scan}(scanNum).scanInfo(slcNum).halfLife = ...
                halfLife;
            planC{indexS.scan}(scanNum).scanInfo(slcNum).injectedDose = ...
                injectedDose;
            planC{indexS.scan}(scanNum).scanInfo(slcNum).seriesTime = ...
                seriesTime;
            planC{indexS.scan}(scanNum).scanInfo(slcNum).decayCorrection = ...
                decayCorrection;
            planC{indexS.scan}(scanNum).scanInfo(slcNum).correctedImage = ...
                correctedImage;
            planC{indexS.scan}(scanNum).scanInfo(slcNum).imageUnits = ...
                imageUnits;

        end
    end
end

% Check for GSPS and make it empty if no objects are present
if length(planC{indexS.GSPS}) == 1 && isempty(planC{indexS.GSPS}.SOPInstanceUID)
    planC{indexS.GSPS}(1) = [];
end

% Check for IM UIDs
if ~isempty(planC{indexS.IM}) && isfield(planC{indexS.IM}(1),'IMUID') && isempty(planC{indexS.IM}(1).IMUID)
    planC = createIMuids(planC);
end

% Check visible flag for structures. Set it to 1 if it is an empty string
for i = 1:length(planC{indexS.structures})
    if isempty(planC{indexS.structures}(i).visible)
        planC{indexS.structures}(i).visible = 1;
    end
end

% Overwrite the existing CERR file if a bug is found and fixed
if forceSaveFlag == 1 || (~isempty(stateS) && isfield(stateS.optS,'overwrite_CERR_File') && stateS.optS.overwrite_CERR_File == 1 && bug_found)
    try
        if exist('fileName','var')
            % do nothing
        elseif isfield(stateS,'CERRFile')
            fileName = stateS.CERRFile;
        end
        planC = save_planC(planC, stateS.optS, 'passed', fileName);
    catch
        disp('Could not overwrite the exisitng file. Please save manually')
    end
end

author: cerr


text:
classdef IBEA < ALGORITHM
% <multi/many> <real/integer/label/binary/permutation>
% Indicator-based evolutionary algorithm
% kappa --- 0.05 --- Fitness scaling factor

%------------------------------- Reference --------------------------------
% E. Zitzler and S. Kunzli, Indicator-based selection in multiobjective
% search, Proceedings of the International Conference on Parallel Problem
% Solving from Nature, 2004, 832-842.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            kappa = Algorithm.ParameterSet(0.05);

            %% Generate random population
            Population = Problem.Initialization();

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = TournamentSelection(2,Problem.N,-CalFitness(Population.objs,kappa));
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                Population = EnvironmentalSelection([Population,Offspring],Problem.N,kappa);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the resting lengths and/or spring stiffness
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function springs = update_Springs(dt,current_time,xLag,yLag,springs)

% dt: time-step
% current_time: current time of simulation
% xLag: Vector of all xLagrangian Pts
% springs: col 1: starting spring pt (by lag. discretization)
%          col 2: ending spring pt. (by lag. discretization)
%          col 3: spring stiffness
%          col 4: spring resting lengths

s_1 = springs(:,1);    % MASTER Node Indices
s_2 = springs(:,2);    % SLAVE Node Indices
k_Vec = springs(:,3);  % Springs Stiffnesses Vector
RL_Vec= springs(:,4);  % Resting Lengths Vectors

N = length(xLag);                   % Gives total number of Lagrangian pts!
N_springs = length( springs(:,1) ); % Gives total number of springs!

freq = 10;  % Frequency for pumping
d = 1.0;    % Diameter of Heart Tube

for i=(N-2)+10:(N-2)+20            % Loops over desired springs!

    RL_Vec(i) = d - abs( 0.9*d*sin( 2*pi*freq*current_time ) );

end

springs(:,4) = RL_Vec; % Update the springs_info

author: nickabattista


text:
function SimCurveResults = SIRFSE_SimCurve(Fit, Prot, FitOpt, interpolate )
if ~exist('interpolate','var'), interpolate=1; end

FitOpt.R1map = 0;
FitOpt.R1reqR1f = 0;

xFit = [Fit.F, Fit.kr, Fit.R1f, Fit.R1r, -1*Fit.Sf, Fit.Sr, Fit.M0f];
Fit.table = xFit';

if (all(Prot.td == Prot.td(1))) && interpolate
    Prot.ti = sort(Prot.ti);
    Fit.ti = (0.5*Prot.ti(1):0.001:1.5*Prot.ti(end))';
    Fit.td = Prot.td(1)*ones(length(Fit.ti),1);
    Fit.curve = SIRFSE_fun(xFit,[Fit.ti, Fit.td], FitOpt);
else
    [Prot.ti, order] = sort(Prot.ti);
    Fit.ti = Prot.ti;
    Fit.td = Prot.td(order);
    Fit.curve(order,1) = SIRFSE_fun(xFit,[Fit.ti, Fit.td], FitOpt);
end

SimCurveResults = Fit;


text:
function report = generateIssueReport(EEG)
% Generates an issue report for an EEG structure that has been robustly referenced
aveCorrThreshold = 0.91;
medCorrThreshold = 0.95;
if nargin < 1
    error('generateIssueReport:NotEnoughArgs', ...
        'generateIssueReport requires an EEG structure argument');
end
report = '';

if ~isstruct(EEG)
    report = 'Input not an EEG structure';
    return;
end

nDetect = getFieldIfExists(EEG, {'etc', 'noiseDetection'});
if isempty(nDetect)
    report = 'Input EEG has not been robustly referenced';
    return;
end

status = getFieldIfExists(nDetect, {'errors', 'status'});
if ~isempty(status) && ~strcmpi(status, 'good');
    report = ['EEG referencing has following errors: [' ...
               getStructureString(nDetect.errors) ']'];
    return;
end

badChans = getFieldIfExists(nDetect, ...
    {'reference', 'noisyStatistics', 'noisyChannels' 'all'});
if ~isempty(badChans)
    report = [report ...
        sprintf('Noisy channels not interpolated after referencing: %s\n', ...
        getListString(badChans))];
end

evaluationChans = length(getFieldIfExists(nDetect, ...
    {'reference', 'noisyStatistics', 'evaluationChannels'}));
if 0.25*evaluationChans < length(badChans)
     report = [report sprintf('Data set has %d of %d EEG channels bad\n', ...
               length(badChans), evaluationChannels)];
end

%% Median
theCorr = getFieldIfExists(nDetect, ...
         {'reference', 'noisyStatistics', 'maximumCorrelations'});
if ~isempty(theCorr)
    meanCorr = mean(theCorr(:));
    medianCorr = median(theCorr(:));
    if meanCorr> aveCorrThreshold || medianCorr > medCorrThreshold
      report = [report ...
          sprintf('Max win correlation [median=%g, mean=%g]\n', ...
                medianCorr, meanCorr)];
    end
end

devRef = getFieldIfExists(nDetect, ...
         {'reference', 'noisyStatistics', 'channelDeviations'});
devOrig = getFieldIfExists(nDetect, ...
         {'reference', 'noisyStatisticsOriginal', 'channelDeviations'});
if ~isempty(devRef) && ~isempty(devOrig)
   devRef = median(devRef(:));
   devOrig = median(devOrig(:));
   if devOrig < devRef
   report = [report  ...
        sprintf('Referencing did not improve amplitude [ref=%g, orig=%g]\n', ...
                devRev, devOrig)];
   end
end

author: VisLab

repo: EEG-Clean-Tools


text:
function c = mvn_vbcost(x,L,mu,L_mu,L_K,logdet_K)
warning('This function is deprecated')

% c = mvnvbcost(x,Covx,mu,Covmu,pCholCov,pLogDetCov)
%
% Calculates c = -KL(q||p) = <log p(X)> - <log q(X)>
% where the expectation <.> is taken over q(X). This is used for
% calculating the lower bound of the marginal loglikelihood in VB
% models.
%
% p(X) = N(M,K)
% q(X) = N(x,Covx)
%
% Rest of the parameters are defined as:
% q(M) = N(mu,Covmu)
% <inv K> = inv(L_K * L_K')
% <log det K> = logdet_K
%
% If Covx==[], then the term <log q(X)> is not calculated.
% This is useful when X is, e.g., observations.

c = 0;
d = length(x);

% Cost from q-posterior
if ~isempty(L)
  c = mnorm_entropy(L);
else
  L = 0;
end

% Use Cholesky of the prior covariance
z = solve_tril(L_K, x-mu);

% Cost from prior

% Below: (x-mu)' * Cov^(-1) * (x-mu) + trace(Cov^(-1)*(Covx+Covmu))
V = linsolve_chol(L_K, L*L'+L_mu*L_mu')
err2 = z'*z + trace(V, V); % TODO: Optimize!
% $$$ err2 = z'*z + trace(solve_triu(pCholCov',solve_tril(pCholCov,Covx+Covmu)));
c = c - 0.5*logdet_K - 0.5*err2 - 0.5*d*log(2*pi);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function g = coeff_var_dct(I)

temp1=dct2(I);
temp2=temp1(:);
temp3=temp2(2:end);

%g=kurtosis(temp3);
g=coeff_var_gen_gauss(temp3);

repo: sr-metric


text:
classdef FWModes
    enumeration
        FAT32_ZIP, DISK_IMAGE
    end
end

author: analogdevicesinc


text:
function [outputFile] = reportPDF(model, microbeID, biomassReaction, inputDataFolder, reportFolder, ncbiID)
% Runs all test functions and writes a report in a .tex file. Writes a PDF
% using pdf2latex.
% Requires a LaTex installation, e.g. MiKTex (https://miktex.org/download)
% and pdftex (https://ctan.org/pkg/pdftex).
%
% INPUT
% model             COBRA model structure
% microbeID         Microbe ID in carbon source data file
% biomassReaction   Model biomass reaction
% inputDataFolder   Folder with experimental data and database files
%                   to load
% reportFolder      Path to folder where report documents should be written
%                   (e.g., 'C:\Reports')
%
% OPTIONAL INPUT
% ncbiID            Organism NCBI ID
%
% OUTPUT
% outputFile        Name of the output file with report
%
% Stefania Magnusdottir, Nov 2017
% Almut Heinken, Sep 2018-adapted nomenclature

% load database
database=loadVMHDatabase;

% switch directory to create the LaTex file
currentDir = pwd;
cd(reportFolder)

% Open file to write report (using LaTeX)
% fileNameTex = [reportFolder filesep datestr(datetime('now'), 'yyyymmdd'), '_', microbeID, '.tex'];
fileNameTex = [reportFolder filesep  microbeID, '.tex'];
fid = fopen(fileNameTex, 'w');
cd(currentDir);

% Print preamble to file
fprintf(fid, '\\documentclass[11pt,a4paper]{article}\n');
fprintf(fid, '\\title{%s}\n', strrep(microbeID, '_', '\_'));
fprintf(fid, '\\date{\\today}\n');
fprintf(fid, '\\begin{document}\n');
fprintf(fid, '\\maketitle\n');
fprintf(fid, '\n');

% Taxonomic lineage
fprintf(fid, '\\section{Taxonomic lineage}\n');
if ~isempty(ncbiID)
    taxonomy = parseNCBItaxonomy(ncbiID);
    taxFields = fieldnames(taxonomy);
    if ~isempty(taxFields)
        fprintf(fid, '\\begin{tabular}{rl}\n');
        for i = 1:length(taxFields)
            fprintf(fid, '%s: & %s \\\\\n', taxFields{i}, strrep(taxonomy.(taxFields{i}), '_', '\_'));
        end
        fprintf(fid, '\\end{tabular}\n');
    end
    fprintf(fid, '\n');
end

% Reconstruction properties
fprintf(fid, '\\section{Reconstruction properties}\n');
fprintf(fid, '\\begin{tabular}{rl}\n');
fprintf(fid, 'Reactions: & %d \\\\\n', length(model.rxns));
fprintf(fid, 'Metabolites: & %d \\\\\n', length(model.mets));
fprintf(fid, 'Genes: & %d \\\\\n', length(model.genes));
fprintf(fid, '\\end{tabular}\n');
fprintf(fid, '\n');

% Mass and charge balance
fprintf(fid, '\\section{Mass and charge balance}\n');
[massImbalancedRxns, chargeImbalancedRxns, ~, metsMissingFormulas] = ...
    testModelMassChargeBalance(model, true, biomassReaction);
if size(massImbalancedRxns, 1) > 1
    fprintf(fid, '\\subsection{Mass imbalanced reactions}\n');
    fprintf(fid, '\\begin{tabular}{lc}\n');
    fprintf(fid, 'Reaction & Element imbalances\\\\\n');
    fprintf(fid, '\\hline\\\\\n');
    for i = 2:size(massImbalancedRxns, 1)
        fprintf(fid, '%s & %s \\\\\n', strrep(massImbalancedRxns{i, 1}, '_', '\_'), ...
                massImbalancedRxns{i, 3});
    end
    fprintf(fid, '\\end{tabular}\n');
    fprintf(fid, '\n');
end
if size(chargeImbalancedRxns, 1) > 1
    fprintf(fid, '\\subsection{Charge imbalanced reactions}\n');
    fprintf(fid, '\\begin{tabular}{lc}\n');
    fprintf(fid, 'Reaction & Charge imbalance\\\\\n');
    fprintf(fid, '\\hline\\\\\n');
    for i = 2:size(chargeImbalancedRxns, 1)
        fprintf(fid, '%s & %s \\\\\n', strrep(chargeImbalancedRxns{i, 1}, '_', '\_'), ...
                chargeImbalancedRxns{i, 3});
    end
    fprintf(fid, '\\end{tabular}\n');
    fprintf(fid, '\n');
end
if ~isempty(metsMissingFormulas)
    fprintf(fid, '\\subsection{Metabolites missing formulas}\n');
    fprintf(fid, '\\begin{tabular}{ll}\n');
    fprintf(fid, 'Metabolite & Name \\\\\n');
    fprintf(fid, '\\hline\\\\\n');
    for i = 1:size(metsMissingFormulas, 1)
        fprintf(fid, '%s & %s\\\\\n', strrep(metsMissingFormulas{i, 1}, '_', '\_'), ...
                strrep(metsMissingFormulas{i, 2}, '_', '\_'));
    end
    fprintf(fid, '\\end{tabular}\n');
    fprintf(fid, '\n');
end

% Leak test
fprintf(fid, '\\section{Leak test}\n');
leakingMets = testLeakingMetabolites(model);
if ~isempty(leakingMets)
    fprintf(fid, '\\begin{tabular}{l}\n');
    for i = 1:length(leakingMets)
        fprintf(fid, '%s\\\\\n', strrep(leakingMets{i, 1}, '_', '\_'));
    end
    fprintf(fid, '\\end{tabular}\n');
    fprintf(fid, '\n');
else
    fprintf(fid, 'No leaking metabolites.\n');
end

% ATP production
fprintf(fid, '\\section{ATP production}\n');
[atpFluxAerobic, atpFluxAnaerobic] = testATP(model);
fprintf(fid, '\\begin{tabular}{ll}\n');
fprintf(fid, 'ATP demand (aerobic): & %4.4f mmol/gDW/h\\\\\n', atpFluxAerobic);
fprintf(fid, 'ATP demand (anaerobic): & %4.4f mmol/gDW/h\\\\\n', atpFluxAnaerobic);
fprintf(fid, '\\end{tabular}\n');
fprintf(fid, '\n');

% Growth
fprintf(fid, '\\section{Growth}\n');
[AerobicGrowth, AnaerobicGrowth] = testGrowth(model, biomassReaction);
fprintf(fid, '\\begin{tabular}{ll}\n');
fprintf(fid, 'Unlimited growth (aerobic): & %4.4f mmol/gDW/h\\\\\n', AerobicGrowth(1, 1));
fprintf(fid, 'Unlimited growth (anaerobic): & %4.4f mmol/gDW/h\\\\\n', AnaerobicGrowth(1, 1));
fprintf(fid, 'C omplex medium growth (aerobic): & %4.4f mmol/gDW/h\\\\\n', AerobicGrowth(1, 2));
fprintf(fid, 'Complex medium growth (anaerobic): & %4.4f mmol/gDW/h\\\\\n', AnaerobicGrowth(1, 2));
fprintf(fid, '\\end{tabular}\n');
fprintf(fid, '\n');

% Carbon sources
fprintf(fid, '\\section{Carbon sources}\n');
[TruePositives, FalseNegatives] = testCarbonSources(model, microbeID, biomassReaction, database, inputDataFolder);
if ~isempty(FalseNegatives) || ~isempty(TruePositives)
    fprintf(fid, '\\begin{tabular}{ll}\n');
    fprintf(fid, '\\textit{In vitro} carbon source & Taken up by model\\\\\n');
    fprintf(fid, '\\hline\\\\\n');
    if ~isempty(FalseNegatives)
        for i = 1:length(FalseNegatives)
            fprintf(fid, '%s & %s\\\\\n', FalseNegatives{i}, 'FALSE');
        end
    end
    if ~isempty(TruePositives)
        for i = 1:length(TruePositives)
            fprintf(fid, '%s & %s\\\\\n', TruePositives{i}, 'TRUE');
        end
    end
    fprintf(fid, '\\end{tabular}\n');
    fprintf(fid, '\n');
else
    fprintf(fid, 'No carbon sources reported for organism.\n');
end

% Fermentation products
fprintf(fid, '\\section{Fermentation products}\n');
[TruePositives, FalseNegatives] = testFermentationProducts(model, microbeID, biomassReaction, database, inputDataFolder);
if ~isempty(FalseNegatives) || ~isempty(TruePositives)
    fprintf(fid, '\\begin{tabular}{ll}\n');
    fprintf(fid, '\\textit{In vitro} fermentation product & Secreted by model\\\\\n');
    fprintf(fid, '\\hline\\\\\n');
    if ~isempty(FalseNegatives)
        for i = 1:length(FalseNegatives)
            fprintf(fid, '%s & %s\\\\\n', FalseNegatives{i}, 'FALSE');
        end
    end
    if ~isempty(TruePositives)
        for i = 1:length(TruePositives)
            fprintf(fid, '%s & %s\\\\\n', TruePositives{i}, 'TRUE');
        end
    end
    fprintf(fid, '\\end{tabular}\n');
    fprintf(fid, '\n');
else
    fprintf(fid, 'No fermentation products reported for organism.\n');
end

% Uptake of known consumed metabolites
fprintf(fid, '\\section{Known consumed metabolites}\n');
[TruePositives, FalseNegatives] = testMetaboliteUptake(model, microbeID, biomassReaction, database, inputDataFolder);
if ~isempty(FalseNegatives) || ~isempty(TruePositives)
    fprintf(fid, '\\begin{tabular}{ll}\n');
    fprintf(fid, 'Consumed metabolites in model\\\\\n');
    fprintf(fid, '\\hline\\\\\n');
    if ~isempty(FalseNegatives)
        for i = 1:length(FalseNegatives)
            fprintf(fid, '%s & %s\\\\\n', FalseNegatives{i}, 'FALSE');
        end
    end
    if ~isempty(TruePositives)
        for i = 1:length(TruePositives)
            fprintf(fid, '%s & %s\\\\\n', TruePositives{i}, 'TRUE');
        end
    end
    fprintf(fid, '\\end{tabular}\n');
    fprintf(fid, '\n');
else
    fprintf(fid, 'No known consumed metabolites reported for organism.\n');
end

% Production of known secretion products
fprintf(fid, '\\section{Known secretion products}\n');
[TruePositives, FalseNegatives] = testSecretionProducts(model, microbeID, biomassReaction, database, inputDataFolder);
if ~isempty(FalseNegatives) || ~isempty(TruePositives)
    fprintf(fid, '\\begin{tabular}{ll}\n');
    fprintf(fid, 'Secretion products in model\\\\\n');
    fprintf(fid, '\\hline\\\\\n');
    if ~isempty(FalseNegatives)
        for i = 1:length(FalseNegatives)
            fprintf(fid, '%s & %s\\\\\n', FalseNegatives{i}, 'FALSE');
        end
    end
    if ~isempty(TruePositives)
        for i = 1:length(TruePositives)
            fprintf(fid, '%s & %s\\\\\n', TruePositives{i}, 'TRUE');
        end
    end
    fprintf(fid, '\\end{tabular}\n');
    fprintf(fid, '\n');
else
    fprintf(fid, 'No known secretion products reported for organism.\n');
end

% Growth on defined media
[growsOnDefinedMedium,constrainedModel,growthOnKnownCarbonSources] = testGrowthOnDefinedMedia(model, microbeID, biomassReaction, inputDataFolder);
fprintf(fid, '\\section{Growth on defined medium}\n');
if ~strcmp(growsOnDefinedMedium,'NA')
    fprintf(fid, '\\begin{tabular}{ll}\n');
    if growsOnDefinedMedium==1
        fprintf(fid, 'The model grows on defined medium with at least one carbon source.\n');
    elseif growsOnDefinedMedium==0
        fprintf(fid, 'The model cannot on defined medium reported for the organism.\n');
    end
        fprintf(fid, '\\end{tabular}\n');
    fprintf(fid, '\n');
%     fprintf(fid, '\\subsection{Growth rates on defined medium}\n');
%     fprintf(fid, '\\begin{tabular}{ll}\n');
%     fprintf(fid, 'Carbon source & Aerobic & Anaerobic\\\\\n');
%     fprintf(fid, '\\hline\\\\\n');
%     for i=1:size(growthOnKnownCarbonSources,1)
%         fprintf(fid, '%s & %s\\\\\n', growthOnKnownCarbonSources{i,1}, str2double(growthOnKnownCarbonSources{i,2}), str2double(growthOnKnownCarbonSources{i,3}));
%     end
%     fprintf(fid, '\\end{tabular}\n');
%     fprintf(fid, '\n');
else
    fprintf(fid, 'No growth media reported for organism.\n');
end

% Bile acid biosynthesis
fprintf(fid, '\\section{Bile acid biosynthesis}\n');
[TruePositives, FalseNegatives] = testBileAcidBiosynthesis(model, microbeID, biomassReaction, database);
if ~isempty(FalseNegatives) || ~isempty(TruePositives)
    fprintf(fid, '\\begin{tabular}{ll}\n');
    fprintf(fid, 'Bile acid biosynthesis in model\\\\\n');
    fprintf(fid, '\\hline\\\\\n');
    if ~isempty(FalseNegatives)
        for i = 1:length(FalseNegatives)
            fprintf(fid, '%s & %s\\\\\n', FalseNegatives{i}, 'FALSE');
        end
    end
    if ~isempty(TruePositives)
        for i = 1:length(TruePositives)
            fprintf(fid, '%s & %s\\\\\n', TruePositives{i}, 'TRUE');
        end
    end
    fprintf(fid, '\\end{tabular}\n');
    fprintf(fid, '\n');
else
    fprintf(fid, 'No bile acid biosynthesis reported for organism.\n');
end

% Close document
fprintf(fid, '\\end{document}\n');
fclose(fid);

% write PDF
cd(reportFolder)
system(['pdflatex ' fileNameTex]);
cd(currentDir);

outputFile = strcat(microbeID, '.pdf');
end

author: opencobra

repo: cobratoolbox


text:
function [ w, xy ] = triangle_unit_o07 ( )

%*****************************************************************************80
%
%% TRIANGLE_UNIT_O07 returns a 7 point quadrature rule for the unit triangle.
%
%  Discussion:
%
%    This rule is precise for monomials through degree 5.
%
%    The integration region is:
%
%      0 <= X
%      0 <= Y
%      X + Y <= 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carlos Felippa,
%    A compendium of FEM integration formulas for symbolic work,
%    Engineering Computation,
%    Volume 21, Number 8, 2004, pages 867-890.
%
%  Parameters:
%
%    Output, real W(7), the weights.
%
%    Output, real XY(2,7), the abscissas.
%
  w(1:7,1) = [ ...
    0.12593918054482715260, ...
    0.12593918054482715260, ...
    0.12593918054482715260, ...
    0.13239415278850618074, ...
    0.13239415278850618074, ...
    0.13239415278850618074, ...
    0.22500000000000000000 ];

  xy(1:2,1:7) = [ ...
    0.79742698535308732240,  0.10128650732345633880; ...
    0.10128650732345633880,  0.79742698535308732240; ...
    0.10128650732345633880,  0.10128650732345633880; ...
    0.059715871789769820459,  0.47014206410511508977; ...
    0.47014206410511508977,  0.059715871789769820459; ...
    0.47014206410511508977,  0.47014206410511508977; ...
    0.33333333333333333333,  0.33333333333333333333 ]';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function A = bezout(p,q)
%BEZOUT       Bezout matrix of two univariate polynomials p and q
%
%   A = bezout(p,q)
%
%Both polynomials must be univariate with the same independent variable.
%The (symmetric) Bezout matrix is singular iff p and q have a root in common.
%For p or q being interval polynomials, A is an interval matrix.
%The call
%
%   A = bezout(p)
%
%is the same as   A = bezout(p,p') .
%The definition follows Fiedler, Special matrices and their applications.
%

% written  02/03/04     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if nargin==1
    q = p';
  end

  if ( size(p.v)>1 ) | ( size(q.v)>1 )
    error('both polynomials must be univariate')
  end
  if p.v~=q.v
    error('both polynomials must depend on the same variable')
  end

  np = p.e;
  nq = q.e;
  n = max(np,nq);

  A = toeplitz([p.c(np+1) zeros(1,n-1)],[fliplr(p.c) zeros(1,2*n-np-1)]);
  B = toeplitz([q.c(nq+1) zeros(1,n-1)],[fliplr(q.c) zeros(1,2*n-nq-1)]);
  A = flipud( A(:,n+1:end)*B(:,1:n) - B(:,n+1:end)*A(:,1:n) );

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function [rdiff, Z, pval, stats] = correl_compare_indep_inputr(r1, r2, N1, N2, varargin)
% Compare two Pearson's correlation values collected from independent samples
% :Based on: www.stat-help.com/
% :Usage:
% ::
%     [rdiff, Z, pval, stats] = correl_compare_indep_inputr(r1, r2, N1, N2, varargin)
%
% :Inputs:
%
%   **r1:**
%        Correlation coefficient for group 1
%
%   **r2:**
%        Correlation coefficient for group 2
%
%   **N1:**
%        Sample size for group 1
%
%   **N2:**
%        Sample size for group 2
%
% :Outputs:
%
%   **rdiff:**
%        Difference in correlation coefficient matrices (r1 - r2)
%
%   **Z:**
%        Z score for difference in correlation coefficients.
%
%   **pval:**
%        p value for two-tailed z-test
%
%   **stats:**
%        Structure of variables including r1, r2, N1, N2, rdiff, zdiff, Z,
%        pval, myalpha, sig
%
% :Examples:
% ::
%
%    %generate two random input matrices (20 subjects x 5 variables)
%    r1 = .77;
%    r2 = .33;
%    N1 = 100;
%    N2 = 100;
%
%   [rdiff, Z, pval, stats] = correl_compare_indep_inputr(r1, r2, N1, N2);
%
%
% :See also:
%   - correl_compare_indep_inputr, correl_compare_permute*
%
%
% ..
%    Tor Wager, March 2010
% ..

    if nargin == 0
        disp('Using sample values for r1 r2 N1 N2')
        N1 = 100;
        N2 = 100;
        r1 = .33;
        r2 = .77;
    end

    rdiff = r1 - r2;

    sediff = sqrt(1/(N1-3) + 1/(N2-3));  % Standard error of difference between two standard normals, with df = N - 3 for correlation coeff
    zdiff = fisherz(r1) - fisherz(r2);  % difference in Fisher's z values

    Z = zdiff ./ sediff;   % Z-score (inferential stat)

    pval = 2 * (1 - normcdf(abs(Z))); % p-value

    myalpha = .05;  % ROI
    sig = pval < myalpha;

    stats = struct('r1', r1, 'r2', r2, 'N1', N1, 'N2', N2, 'diff_r_values', rdiff, ...
        'fishers_z_diff', zdiff, 'Z_stat', Z, 'p', pval, 'alpha', myalpha, 'sig', sig);

end

%out = correl_compare_dep(y1,y2, 'alpha',myalpha)

author: canlab

repo: CanlabCore


text:
 function y = ir_imfilter_many(x, psf, varargin)
%function y = ir_imfilter_many(x, psf, varargin)
%|
%| as of 2015-04-12, octave's imfilter can handle only 2D input
%| whereas matlab can handle 3D inputs
%| this is a kludge to deal with that.
%|
%| 2015-04-12, Jeff Fessler, University of Michigan

if nargin < 2, ir_usage, end

if ~ir_is_octave || ndims(x) == ndims(psf)
	y = imfilter(x, psf, varargin{:}); % matlab or usual sizes
return
end

% hereafter for octave with unequal sizes

if ndims(x) == 3 && ndims(psf) == 2
	tmp = imfilter(x(:,:,1), psf, varargin{:});
	nz = size(x,3);
	y = zeros([size(tmp) nz]);
	y(:,:,1) = tmp;
	for iz = 2:nz
		y(:,:,iz) = imfilter(x(:,:,iz), psf, varargin{:});
	end
return
end

pr size(x)
pr size(psf)
fail('imfilter in octave does not support such sizes, nor does my kludge')

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% see also retreating_stall_envelope.m

profile_cd = .008; %blade Cd0
ClimbEfficiency = .9;
Kind = 1.05;
max_tip_M = .92;


text:
function p = predict(theta, X)
%PREDICT Predict whether the label is 0 or 1 using learned logistic
%regression parameters theta
%   p = PREDICT(theta, X) computes the predictions for X using a
%   threshold at 0.5 (i.e., if sigmoid(theta'*x) >= 0.5, predict 1)

m = size(X, 1); % Number of training examples

% You need to return the following variables correctly
p = zeros(m, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters.
%               You should set p to a vector of 0's and 1's
%

for i=1:m,
	if sigmoid(X(i,:)*theta) >= 0.5,
		p(i) = 1;
	else
		p(i) = 0;
	end;
end;

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

end

author: UtkarshPathrabe


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function fmpart
%FMPART	Unit test for the function FMPAR.

%	O. Lemoine - February 1996.

N=256;
t=1:N;

% NARGIN=2
P1=[0.5 -0.007 2.7*10^(-5)];
[x,iflaw]=fmpar(N,P1);

if any(abs(iflaw-(P1(1)+P1(2).*t+P1(3).*t.^2)')>sqrt(eps))~=0,
  error('fmpar test 1 failed');
end

% NARGIN=4
P1=[10 0.4];
P2=[120 0.05];
P3=[246 0.39];
[x,iflaw]=fmpar(N,P1,P2,P3);
if abs(iflaw(P1(1))-P1(2))>sqrt(eps) | abs(iflaw(P2(1))-P2(2))>sqrt(eps) | ...
   abs(iflaw(P3(1))-P3(2))>sqrt(eps),
  error('fmpar test 2 failed');
end

P1=[20 0.13];
P2=[145 0.45];
P3=[226 0.19];
[x,iflaw]=fmpar(N,P1,P2,P3);
if abs(iflaw(P1(1))-P1(2))>sqrt(eps) | abs(iflaw(P2(1))-P2(2))>sqrt(eps) | ...
   abs(iflaw(P3(1))-P3(2))>sqrt(eps),
  error('fmpar test 3 failed');
end

N=251;
t=1:N;

% NARGIN=2
P1=[0.5 -0.007 2.7*10^(-5)];
[x,iflaw]=fmpar(N,P1);

if any(abs(iflaw-(P1(1)+P1(2).*t+P1(3).*t.^2)')>sqrt(eps))~=0,
  error('fmpar test 4 failed');
end

% NARGIN=4
P1=[10 0.4];
P2=[120 0.05];
P3=[246 0.39];
[x,iflaw]=fmpar(N,P1,P2,P3);
if abs(iflaw(P1(1))-P1(2))>sqrt(eps) | abs(iflaw(P2(1))-P2(2))>sqrt(eps) | ...
   abs(iflaw(P3(1))-P3(2))>sqrt(eps),
  error('fmpar test 5 failed');
end

P1=[20 0.13];
P2=[145 0.45];
P3=[226 0.19];
[x,iflaw]=fmpar(N,P1,P2,P3);
if abs(iflaw(P1(1))-P1(2))>sqrt(eps) | abs(iflaw(P2(1))-P2(2))>sqrt(eps) | ...
   abs(iflaw(P3(1))-P3(2))>sqrt(eps),
  error('fmpar test 6 failed');
end

author: HeLiangHIT


text:
function [lags,use_ascend,dc,sc,E] = gc_lags(points,options)

% GC_LAGS       Run lag extraction using a binary graph cut
%
%   Run lag extraction on a reordered raster plot (2D Image) using a binary graph cut
%
%   SYNTAX
%       [LAGS,USE_ASCEND,DC,SC] = GC_LAGS(POINTS,OPTIONS)
%
%   DC is the data cost
%   SC is the smoothness cost
%

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

if nargin<2
    options.null = 0;
end

if ~isfield(options, 'alpha')
    options.alpha = 0.01;
end
alpha = options.alpha;

if ~isfield(options, 'verbose')
    options.verbose = true;
end
verbose = options.verbose;

if ~isfield(options, 'use_lcurve')
    options.use_lcurve = false;
end
use_lcurve = options.use_lcurve;

if ~isfield(options, 'maxit') % maximum number of iterations for alpha search in l-curve
    options.maxit = 40;
end
maxit = options.maxit;

if ~isfield(options, 'dc_tol') % tolerance on data cost for alpha search in l-curve
    options.dc_tol = 1e-25;
end
dc_tol = options.dc_tol;

if isfield(options, 'null')
    options = rmfield(options,'null');
end

%%%%%%%% END OPTIONS %%%%%%%%%

points = - points; % min cut through maximum of points values
points = points - min(points(:)); % only positive capacities for the maxflow

npoints = size(points,1);

if options.disp_log
    disp(['Nb trials : ',num2str(npoints)]);
    disp(['Nb time samples : ',num2str(size(points,2))]);
    disp(['Lambda : ',num2str(alpha)]);
end

if use_lcurve
    disp('Computing best alpha with L-curve');
    verbose = false;
    sc = Inf;
    ht = java.util.Hashtable;

    % Add alpha = 0 to ht
    [lags0,use_ascend,dc0,sc0] = gc_lags_aux(0.0);

    % Find the biggest alpha
    while sc > 0
        fprintf('.');
        [lags,use_ascend,dc,sc] = gc_lags_aux(alpha);
        ht.put(alpha,[dc sc]);
        alpha = 3*alpha;
        % alpha = 2.5*alpha;
        % alpha = 5*alpha;
        % alpha = 2*alpha;
        % alpha = 1.5*alpha;
    end
    fprintf('\n Biggest alpha : %f\n',alpha);
    % Try to find best alpha iteratively
    goon = true;
    niter = 1;
    while goon
        [alphas,dcs,scs] = ht_to_array();

        % Make sure dcs and scs values are unique
        [tmp,dcs_idx] = unique(dcs,'first');
        [tmp,scs_idx] = unique(scs,'first');
        good_idx = sort(intersect(dcs_idx,scs_idx));
        dcs = dcs(good_idx);
        scs = scs(good_idx);
        alphas = alphas(good_idx);
        clear good_idx

        scs_intermediate = diff(dcs+alphas.*scs) ./ diff(alphas);
        dcs_intermediate = (dcs(1:end-1)+alphas(1:end-1).*scs(1:end-1)) - alphas(1:end-1) .* scs_intermediate;
        dcs_interp = linspace(min(dcs),max(dcs),500);
        scs_interp = interp1([dcs;dcs_intermediate],[scs;scs_intermediate],dcs_interp,'linear');

        l1 = 1./alphas(1:end-1);
        l2 = 1./alphas(2:end);
        angle_intermediate = acos((1+l1.*l2)./(sqrt(1+l1.^2).*sqrt(1+l2.^2))); % angle between neighboring segments

        angle_intermediate(1) = 0;
        angle_intermediate(end) = 0;
        midx = localmax(angle_intermediate);

        [tmp,best_midx] = max(angle_intermediate);

        best_dcs = dcs_intermediate(best_midx);
        best_scs = scs_intermediate(best_midx);

        best_dcs_all = best_dcs;
        best_scs_all = best_scs;

        % get alpha by computing harmonic mean of neighboring alphas (better than arithmetic mean)
        new_alphas = alphas(midx) .* alphas(midx+1) ./ (alphas(midx)+alphas(midx+1)) * 2;

        % Sort alphas using order of angle differences
        [tmp,I] = sort(angle_intermediate(midx));
        new_alphas = new_alphas(I);

        precision = 0;

        needs_update = false;
        for ii = 1:length(new_alphas)
            alpha = new_alphas(ii); % take in between alpha

            if isempty(ht.get(alpha)) % Check if solution with this alpha has already between computed
                disp(['    Computing solution for alpha : ',num2str(alpha)]);
                [lags,use_ascend,dc,sc] = gc_lags_aux(alpha);
                ht.put(alpha,[dc sc]);
                needs_update = true; % Say if at least one new alpha is added
            else
                sc_dc = ht.get(alpha);
                dc = sc_dc(1);
                sc = sc_dc(2);
            end
            precision = max(precision,min(abs(dcs-dc)));
        end

        if ~needs_update
            disp('Stopping search (all alphas already computed)');
            goon = false;
        end

        disp([' Error : ',num2str(precision)]);
        if precision < dc_tol | niter >= maxit
            goon = false;
        end
        niter = niter+1;

        if 1 % set best final alpha with interpolation
            if ~goon
                [alpha,best_dcs,best_scs] = get_best_alpha();
            end
        end

        if 1 % show l-curve
            smart_figure('Lag extraction L-curve'); clf
            title('Lag extraction L-curve')
            if goon
                color = {'g','b'};
                linewidth = 1;
            else
                color = {'r','r'};
                linewidth = 3;
            end
            USE_LOG = false;
            % USE_LOG = true;
            if USE_LOG % plot log-log or not
                loglog(dcs_interp(dcs_interp>0),scs_interp(dcs_interp>0),color{2});
                hold on
                loglog(dcs(dcs>0),scs(dcs>0),[color{2},'x']);
            else
                plot(dcs_interp,scs_interp,color{2},'LineWidth',linewidth);
                hold on
                plot(dcs,scs,[color{2},'x']);
            end

            for k=1:length(dcs)
                if alphas(k) > 0
                    flow = dcs(k)+alphas(k)*scs(k);
                    line([0 , flow],[flow/alphas(k), 0],'Color','k');
                end
            end

            axis tight
            line([min(dcs_interp)/100,best_dcs],[best_scs,best_scs], ...
                 'LineStyle','--','Color',color{1},'LineWidth',linewidth);
            line([best_dcs,best_dcs],[min(scs_interp)/100,best_scs], ...
                 'LineStyle','--','Color',color{1},'LineWidth',linewidth);
            hold off;
            xlabel('Data term')
            ylabel('Smoothness term')

            if 1 % show max curvature points
                hold on
                if USE_LOG
                    loglog(best_dcs_all,best_scs_all,'bo');
                else
                    plot(best_dcs_all,best_scs_all,'bo');
                end
                hold off
            end
        end
    end
    disp(['L-curve done with final alpha : ',num2str(alpha)]);
    disp(['         tol   : ',num2str(precision)]);
    disp(['         niter : ',num2str(niter)]);
else
    [lags,use_ascend,dc,sc] = gc_lags_aux(alpha);
end

E = dc+alpha*sc;

if options.disp_log
    disp('E = DC + alpha * SC');
    disp(sprintf('%f = %f + %f * %f',dc+alpha*sc,dc,alpha,sc));
end

% ****************************************************************************** %
function [best_alpha,best_dcs,best_scs] = get_best_alpha()

[local_alphas,local_dcs,local_scs] = ht_to_array();
offset = 1; % don't consider extremities of the curve
local_dcs = local_dcs(1+offset:end-offset);
local_scs = local_scs(1+offset:end-offset);
local_alphas = local_alphas(1+offset:end-offset);

local_dcs_interp = linspace(min(local_dcs),max(local_dcs),100);
local_scs_interp = interp_cubic_herm(local_dcs,local_scs,-1./local_alphas,local_dcs_interp);
local_curvature = curvature_2D_diff(local_dcs_interp,local_scs_interp);
local_curvature = abs(local_curvature); % use absolute value
local_curvature(1:5) = NaN;
local_curvature(end-5:end) = NaN;

% [tmp,idx] = max(local_curvature);
[tmp,idx] = min(local_curvature);
best_dcs = local_dcs_interp(idx);
best_scs = local_scs_interp(idx);

[tmp,local_alpha_idx] = min(abs(local_dcs-best_dcs));

best_alpha = local_alphas(local_alpha_idx);

if 0
    smart_figure('Lag extraction L-curve'); hold on
    plot(local_dcs_interp,local_scs_interp,'g')
end

if 0
    smart_figure('SC vs alpha');
    plot(local_alphas,local_scs)
    xlabel('alpha')
    ylabel('SC')
    hold on
    plot(best_alpha,best_scs,'ro');
    hold off
end

if verbose
    smart_figure('Curve interpolation');
    plot(local_dcs_interp,local_scs_interp);
    hold on
    plot(best_dcs,best_scs,'ro');
    hold off
    title('Curve interpolation');
end

end % end get_best_alpha function

% ****************************************************************************** %
function [alphas,dcs,scs] = ht_to_array()
    alphas = zeros(ht.size(),1);
    dcs = zeros(ht.size(),1);
    scs = zeros(ht.size(),1);
    keys = ht.keys();
    cnt = 1;
    while keys.hasNext()
       alpha = keys.next();
       dc_sc = ht.get(alpha);
       dc = dc_sc(1);
       sc = dc_sc(2);
       alphas(cnt) = alpha;
       dcs(cnt) = dc;
       scs(cnt) = sc;
       cnt = cnt+1;
    end;
    [tmp,I] = sort(alphas);
    alphas = alphas(I);
    dcs = dcs(I);
    scs = scs(I);

    dcs = dcs - dc0; % remove offset
end

% ****************************************************************************** %
function idx = localmax(x)
% Find local maxima in an array

idx = find( diff( sign( diff([0; x(:); 0]) ) ) < 0 );

end % end function

% ****************************************************************************** %

function [lags,use_ascend,dc,sc] = gc_lags_aux(alpha)

if options.disp_log
    disp('Testing ascend')
end

order = [npoints:-1:1];
[lags_ascend,flow_ascend,labels_ascend] = gc_aux_mex(points(order,:),alpha);
lags_ascend = lags_ascend(order);
labels_ascend = labels_ascend(order,:);

% Compute manually resulting data cost and smoothness cost
%
% lag_ind = sub2ind(size(points(order,:)),[1:size(points,1)]',lags_ascend(:));
% dc = sum(points(lag_ind))
% sc = sum(abs(diff(lags_ascend)))
% flow_ascend
% fl = dc + alpha * sc

if options.disp_log
    disp('Testing descend')
end
[lags_descend,flow_descend,labels_descend] = gc_aux_mex(points,alpha);

if flow_ascend < flow_descend
    lags = lags_ascend;
    use_ascend = true;
    flow = flow_ascend;
else
    lags = lags_descend;
    use_ascend = false;
    flow = flow_descend;
end

sc = sum(abs(diff(lags))); % Smoothness cost
dc = flow - alpha * sc; % Data cost

end %  function

end

author: fieldtrip

repo: fieldtrip


text:
function [mriCoord] = cs_mni2mri(MRI,mniCoord)
% CS_MNI2MRI: Transform MNI point coordinates (in mm) to MRI coordinate system (in mm)
%
% DEPRECATED: Replace with cs_convert()

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

disp('BST> WARNING: Deprecated function cs_mni2mri(), use cs_convert() instead.');

mriCoord = cs_convert(MRI, 'mni', 'mri', mniCoord' / 1000)' * 1000;

author: brainstorm-tools


text:
function b = bell ( n )

%*****************************************************************************80
%
%% BELL returns the Bell numbers from 0 to N.
%
%  Discussion:
%
%    The Bell number B(N) is the number of restricted growth functions
%    on N.
%
%    Note that the Stirling numbers of the second kind, S^m_n, count the
%    number of partitions of N objects into M classes, and so it is
%    true that
%
%      B(N) = S^1_N + S^2_N + ... + S^N_N.
%
%  Definition:
%
%    The Bell number B(N) is defined as the number of partitions (of
%    any size) of a set of N distinguishable objects.
%
%    A partition of a set is a division of the objects of the set into
%    subsets.
%
%  Examples:
%
%    There are 15 partitions of a set of 4 objects:
%
%      (1234), (123)(4), (124)(3), (12)(34), (12)(3)(4),
%      (134)(2), (13)(24), (13)(2)(4), (14)(23), (1)(234),
%      (1)(23)(4), (14)(2)(3), (1)(24)(3), (1)(2)(34), (1)(2)(3)(4)
%
%    and so B(4) = 15.
%
%  First values:
%
%     N         B(N)
%     0           1
%     1           1
%     2           2
%     3           5
%     4          15
%     5          52
%     6         203
%     7         877
%     8        4140
%     9       21147
%    10      115975
%
%  Recursion:
%
%    B(I) = sum ( 1 <= J <= I ) Binomial ( I-1, J-1 ) * B(I-J)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of Bell numbers desired.
%
%    Output, integer B(1:N+1), the Bell numbers from 0 to N.
%
  b(1) = 1;

  for i = 1 : n
    b(i+1) = 0;
    for j = 1 : i
      b(i+1) = b(i+1) + i4_choose ( i-1, j-1 ) * b(i-j+1);
    end
  end

  return
end


text:
function cullArriveDepartCombosPerCB(mainGUIHandle)
%cullArriveDepartCombosPerCB Summary of this function goes here
%   Detailed explanation goes here

    %Set up CB/Arrive/Depart combo boxes
    [cbBodyNames, arriveDepartNames] = genCBArriveDepartComboBoxStrs(mainGUIHandle);

    %Strings for CB combo boxes
    cbBodyNames = cap1stLetter(cbBodyNames);
    cbComboString = char(cbBodyNames);
    hCentralBodyCombo = findobj(mainGUIHandle,'Tag','centralBodyCombo');

    try
        set(hCentralBodyCombo,'String',cbComboString);
    catch
        hCentralBodyCombo.Items = cbBodyNames;
    end

    %Strings for Depart/Arrival combo boxes
    arriveDepartNames = cap1stLetter(arriveDepartNames);
    departArriveComboString = char(arriveDepartNames);

    if(size(departArriveComboString,1)>5)
        dBValue = 3;
        aBValue = 4;
    elseif(size(departArriveComboString,1)>=2)
        dBValue = 1;
        aBValue = 2;
    else
        dBValue = 1;
        aBValue = 1;
    end

    hDepartBodyCombo = findobj(mainGUIHandle,'Tag','departBodyCombo');
    try
        set(hDepartBodyCombo,'Value',dBValue);
        set(hDepartBodyCombo,'String',departArriveComboString);
    catch
        hDepartBodyCombo.Items = arriveDepartNames;
        if(numel(arriveDepartNames) > 0)
            hDepartBodyCombo.Value = arriveDepartNames{dBValue};
        end
    end

    hArriveBodyCombo = findobj(mainGUIHandle,'Tag','arrivalBodyCombo');
    try
        set(hArriveBodyCombo,'Value',aBValue);
        set(hArriveBodyCombo,'String',departArriveComboString);
    catch
        hArriveBodyCombo.Items = arriveDepartNames;
        if(numel(arriveDepartNames) > 0)
            hArriveBodyCombo.Value = arriveDepartNames{aBValue};
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%  Figure 7.80      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
%% fig7_80.m
%% Heat Exchanger
clf;
[t]=sim('Fig7_78cl');
plot(t,u);
grid;
title('Fig. 7.80: Control effort for heat echanger');
xlabel('Time (sec)');
ylabel('Control, u');

author: Sable


text:
function pols = ortho2eva0 ( mmax, z )

%*****************************************************************************80
%
%% ORTHO2EVA0 evaluates the orthonormal polynomials on the triangle.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    28 June 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer MMAX, the maximum order to which the polynomials are
%    to be evaluated.
%
%    Input, real Z(2), the coordinates of the evaluation point.
%
%    Output, real POLS((mmax+1)*(mmax+2)/2), the orthogonal
%    polynomials evaluated at the point Z.
%
  zero = 0.0;
  sqrt2 = sqrt ( 2.0 );
  sqrt3 = sqrt ( 3.0 );
  r11 = -1.0 / 3.0;
  r12 = -1.0 / sqrt3;
  r21 = - 1.0 / 3.0;
  r22 = 2.0 / sqrt3;

  a = z(1);
  b = z(2);
%
%  Map the reference triangle to the right
%  triangle with the vertices (-1,-1), (1,-1), (-1,1)
%
  x = r11 + r12 * b + a;
  y = r21 + r22 * b;
%
%  Evaluate the Koornwinder's polynomials via the three term recursion.
%
  par1 = ( 2.0 * x + 1.0 + y ) / 2.0;
  par2 = ( 1.0 - y ) / 2.0;
  f1 = klegeypols ( par1, par2, mmax );

  f2 = zeros(mmax+1,mmax+1);
  for m = 0 : mmax
    par1 = 2 * m + 1;
    f2(1:mmax+1,m+1) = kjacopols ( y, par1, zero, mmax - m );
  end

  npols = ( ( mmax + 1 ) * ( mmax + 2 ) ) / 2;
  pols = zeros(npols,1);

  kk = 0;
  for m = 0 : mmax
    for n = 0 : m
      kk = kk + 1;
%
%  Evaluate the polynomial (m-n, n)
%
      pols(kk) = f1(m-n+1) * f2(n+1,m-n+1);
%
%  Normalize.
%
      scale = sqrt ...
      ( ...
        ( ...
          ( 1 + ( m - n ) + n ) * ...
          ( 1 + ( m - n ) + ( m - n ) ) ...
        ) / sqrt3 ...
      );

      pols(kk) = pols(kk) * scale;

    end

  end

  return
end