text:
classdef chebfun3
%CHEBFUN3   CHEBFUN3 class for representing functions on [a,b]x[c,d]x[e,g].
%   Class for approximating functions defined on finite cuboids. The
%   functions should be smooth.
%
%   CHEBFUN3(F) constructs a CHEBFUN3 object representing the function F on
%   [-1, 1] x [-1, 1] x [-1, 1].
%
%   The input F may be any of the following:
%   i)   A function handle like @(x,y,z) x.*y + cos(x).*z
%   ii)  A string like 'cos(x) + sin(x.*y.*z)'
%   iii) A discrete tensor corresponding to values of a function at points
%   generated by ndgrid.
%
%   F should be vectorized in the sense that it may be evaluated at a
%   tensor of points and returns a tensor in the output.
%
%   CHEBFUN3(F, 'eps', ep) specifies chebfun3eps to be ep.
%
%   CHEBFUN3(F, [A B C D E G]) specifies a cuboid [A B] x [C D] x [E G]
%   where the function F is defined. A, B, C, D, E and G must all be finite.
%
%   CHEBFUN3(F, [m n p]) returns a representation of a trivariate
%   polynomial of length (m, n, p), i.e., with degree m-1 in x, degree
%   n-1 in y and degree p-1 in z. The polynomial is compressed in low
%   trilinear rank form and the trilinear rank (rx, ry, rz) is still
%   determined adaptively.
%
%   CHEBFUN3(F, 'rank', [rx ry rz]) returns a CHEBFUN3 that is an
%   approximation to F and has the trilinear rank (rx, ry, rz).
%
%   CHEBFUN3(F, 'trig') or CHEBFUN3(F, 'periodic') constructs a trig-type
%   CHEBFUN3 object if the input function handle F is triply periodic on
%   [-1, 1] x [-1, 1] x [-1, 1]. In this case, factor quasimatrices of the
%   low rank format will be represented by Fourier instead of Chebyshev
%   technology.
%
%   If F is a discrete tensor, F = (f_{ijk}), the numbers f_{ijk} are used
%   as function values at tensor product Chebyshev points of the 2nd kind
%   generated by ndgrid.
%
%   CHEBFUN3(F, 'equi') can be used if F is a discrete tensor of values at
%   equispaced points in 3D.
%
%   CHEBFUN3(F, 'coeffs') returns a CHEBFUN3 in which the coefficients of
%   the Chebyshev expansion are the entries of the discrete tensor F.
%
%  Examples:
%
%   f = chebfun3(@(x,y,z) exp(-(x.^2 + y.^2 + z.^2))); isosurface(f)
%
%   f = chebfun3('sin(10*pi*x).^2 + sin(10*x.*y.*z)'); sum3(f)
%
%   cheb.gallery3
%
%   f1 = chebfun3('tanh(sqrt(3)*x)');   rank(f1), g1 = grad(f1); g1(0,0,0)
%   f2 = chebfun3('tanh(x+y+z)'); rank(f2), g2 = grad(f2); g2(0,0,0)
%
%   Chebfun3 is based on:
%
%   [1] B. Hashemi and L. N. Trefethen, Chebfun in three dimensions, SIAM
%   J. Sci. Comput., 39 (2017), C341-C363.
%
%   [2] S. Dolgov, D. Kressner, and C. Stroessner, Functional Tucker
%   approximation using Chebyshev interpolation, SIAM J. Sci. Comput.,
%   43 (2021), A2190-A2210.
%
%   By default since March 2023, chebfun3 objects are constructed as described
%   in [2].  The original constructor described in [1] can be invoked with
%   CHEBFUN3(F, 'classic').
%
% See also CHEBFUN3V and CHEBFUN3T.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CLASS PROPERTIES:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
properties
    % COLS: Mode-1 fibers, i.e., columns which are functions of x used in
    % Tucker representation.
    cols

    % ROWS: Mode-2 fibers, i.e. rows which are functions of y used in
    % Tucker representation.
    rows

    % TUBES: Mode-3 fibers, i.e. tubes which are functions of z used in
    % Tucker representation.
    tubes

    % CORE: discrete core tensor in Tucker representation
    core

    % DOMAIN: box of CHEBFUN3, default is [-1, 1] x [-1, 1] x [-1, 1].
    domain
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CLASS CONSTRUCTOR:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
methods
    function f = chebfun3(varargin)
        % The main CHEBFUN3 constructor!

        % Return an empty CHEBFUN3:
        if ( (nargin == 0) || isempty(varargin{1}) )
            return
        end

        % Call the constructor, all the work is done here:
        f = constructor(f, varargin{:});
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CLASS METHODS:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
methods (Access = public, Static = true)
    % Outer product of discrete tensors.
    varargout = outerProd(varargin);

    % Tensor x matrix.
    varargout = txm(varargin);

    % Unfold a discrete tensor to create a matrix.
    varargout = unfold(varargin);

    % Reshape a discrete matrix to get a tensor.
    varargout = fold(varargin);

    % Convert 3D values to 3D coefficients.
    coeffs3D = vals2coeffs(vals3D);

    % Convert 3D coefficients to 3D values.
    vals3D = coeffs2vals(coeffs3D);

    % Faster subscripts from linear index.
    varargout = myind2sub(varargin);

    % HOSVD of discrete tensors.
    varargout = discreteHOSVD(varargin);

end

methods (Access = public)
    % Retrieve and modify preferences for this class.
    varargout = subsref(f, index);

    % Get properties of a CHEBFUN3 object.
    out = get(f, idx);

    % Permute a CHEBFUN3.
    out = permute(f, varargin);

    % Evaluate a CHEBFUN3.
    out = feval(f, varargin);

    % Evaluate at vectors to get a tensor of values.
    out = fevalt(f, varargin);

    % Evaluate on a tensor product grid.
    varargout = chebpolyval3(varargin);

    % Sample a CHEBFUN3 on a tensor product grid
    varargout = sample(varargin);

    % Tucker rank of a CHEBFUN3 (i.e., size of the core in each direction).
    varargout = rank(f);

    % Length of a CHEBFUN3 (i.e., the number of Chebyshev or Fourier points
    % at each direction).
    varargout = length(f);

    % Size of a CHEBFUN3.
    varargout = size(f, varargin);

    % Minimum of a CHEBFUN3 along one dimension.
    varargout = min(varargin);

    % Maximum of a CHEBFUN3 along one dimension.
    varargout = max(varargin);

    % Minimum of a CHEBFUN3 along two dimensions.
    varargout = min2(varargin);

    % Maximum of a CHEBFUN3 along two dimensions.
    varargout = max2(varargin);

    % Global minimum of a CHEBFUN3.
    varargout = min3(f);

    % Global maximum of a CHEBFUN3.
    varargout = max3(f);

    % Global minimum and maximum of a CHEBFUN3.
    varargout = minandmax3(f);

    % Norm of a CHEBFUN3.
    varargout = norm(f, p);

    % Display a CHEBFUN3.
    varargout = disp(f, varargin);

    % Display a CHEBFUN3.
    varargout = display(varargin);

    % Simplify a CHEBFUN3.
    out = simplify(varargin);

    % Vertical scale of a CHEBFUN3.
    out = vscale(f);

    % Vertical concatenation of CHEBFUN3 objects.
    out = vertcat(varargin);

    % Just one common root of 3 CHEBFUN3 objects.
    varargout = root(f, g, h);

    % Number of degrees of freedom needed to represent a CHEBFUN3.
    out = ndf(f);

    % Get the low-rank representation (Tucker expansion) of a CHEBFUN3.
    varargout = tucker(f);

    % Definite integral of a CHEBFUN3 over the domain in one
    % direction. Output is a Chebfun2 object.
    out = sum(varargin);

    % Definite integral of a CHEBFUN3 over the domain in two directions.
    % Output is a Chebfun object.
    out = sum2(varargin);

    % Definite integral of a CHEBFUN3 over its domain.
    out = sum3(f);

    % Restrict the domain of a CHEBFUN3. Output might be a scalar,
    % a Chebfun, a Chebfun2 or a Chebfun3 object.
    out = restrict(f, varargin);

    % Volume of the domain of a CHEBFUN3.
    out = domainvolume(f);

    % Line integral of a CHEBFUN3 over a 3D parametric curve.
    out = integral(f, varargin);

    % Surface integral of a CHEBFUN3 over a parametric surface represented
    % as a CHEBFUN2V.
    out = integral2(f, varargin);

    % Integral of a CHEBFUN3 over a restricted cuboid.
    out = integral3(f, varargin);

    % Average or mean value of a CHEBFUN3 in one direction.
    out = mean(f, varargin);

    % Average or mean value of a CHEBFUN3 in two directions.
    out = mean2(f, varargin);

    % Average or mean value of a CHEBFUN3.
    out = mean3(f);

    % Standard deviation of a CHEBFUN3 along one variable.
    out = std(f, varargin);

    % Standard deviation of a CHEBFUN3 along two variables.
    out = std2(f, varargin);

    % Standard deviation of a CHEBFUN3.
    out = std3(f);

    % Squeeze a CHEBFUN3 into a CHEBFUN2 or a CHEBFUN.
    out = squeeze(f);

    % Create a scatter plot of the core tensor of a CHEBFUN3.
    varargout = coreplot(f, varargin);

    % Plot a CHEBFUN3.
    out = plot(f, varargin);

    % Plot slices of a CHEBFUN3.
    out = slice(f, varargin);

    % Scan plot of a CHEBFUN3.
    out = scan(f, varargin);

    % Isosurface plot of a CHEBFUN3.
    out = isosurface(f, varargin);

    % SURF for a CHEBFUN3 over its domain.
    varargout = surf(f, varargin);

    % plotcoeffs of a CHEBFUN3.
    varargout = plotcoeffs(f, varargin);

    % Tensor of coefficients of a CHEBFUN3.
    varargout = chebcoeffs3(f);

    % Tensor of coefficients of a CHEBFUN3 of specified size.
    varargout = coeffs3(f, varargin);
end

methods (Access = private, Static = true)
    % Classic constructor for CHEBFUN3 objects
    f = chebfun3classic(f, op, pref, dom, vectorize, fiberDim)

    % Constructor for CHEBFUN3 objects of tensor inputs.
    f = chebfun3double(f, op, dom, pref, isEqui)

    % Default constructor for CHEBFUN3 objects.
    f = chebfun3f(f, op, pref, dom, vectorize);
end

methods
    % Unary plus for a CHEBFUN3.
    out = uplus(f);

    % Plus for CHEBFUN3 objects.
    out = plus(f, g, tol);

    % Unary minus for a CHEBFUN3.
    out = uminus(f);

    % Subtraction of two CHEBFUN3 objects.
    out = minus(f, g);

    % Pointwise multiplication for CHEBFUN3 objects.
    out = times(f, g);

    % Pointwise multiplication for CHEBFUN3 objects.
    out = mtimes(f, g);

    % Pointwise power of a CHEBFUN3.
    out = power(varargin);

    % Pointwise right divide of CHEBFUN3 objects.
    out = rdivide(f, g);

    % Right divide for CHEBFUN3 objects.
    out = mrdivide(f, g);

    % Pointwise CHEBFUN3 left array divide.
    out = ldivide(f, g);

    % Left divide for CHEBFUN3 objects.
    out = mldivide(f, g);

    % Absolute value of a CHEBFUN3.
    out = abs(f);

    % Real part of a CHEBFUN3.
    out = real(f);

    % Imaginary part of a CHEBFUN3.
    out = imag(f);

    % Complex conjugate of a CHEBFUN3.
    out = conj(f);

    % Create f + i g for two CHEBFUN3 objects.
    out = complex(f, g);

    % Sine of a CHEBFUN3.
    out = sin(f);

    % Cosine of a CHEBFUN3.
    out = cos(f);

    % Tangent of a CHEBFUN3.
    out = tan(f);

    % Tangent of a CHEBFUN3 (in degrees).
    out = tand(f);

    % Hyperbolic tangent of a CHEBFUN3.
    out = tanh(f);

    % Exponential of a CHEBFUN3.
    out = exp(f);

    % Hyperbolic sine of a CHEBFUN3.
    out = sinh(f);

    % Hyperbolic cosine of a CHEBFUN3.
    out = cosh(f);

    % Compose command for CHEBFUN3 objects.
    out = compose(f, varargin);

    % Square root of a CHEBFUN3.
    out = sqrt(f, varargin);

    % Natural logarithm of a CHEBFUN3.
    out = log(f, varargin);

    % HOSVD of a CHEBFUN3.
    varargout = hosvd(f, varargin);

    % Test whether a CHEBFUN3 object is empty.
    out = isempty(f);

    % Test if a CHEBFUN3 is identically zero over its domain.
    varargout = iszero(f);

    % Real-valued CHEBFUN3 test.
    out = isreal(f);

    % Equality test for CHEBFUN3 objects.
    out = isequal(f, g);

    % Partial derivative of a CHEBFUN3 object.
    out = diff(f, varargin);

    % Partial derivative of a CHEBFUN3 object in the first variable.
    out = diffx(f, varargin);

    % Partial derivative of a CHEBFUN3 object in the second variable.
    out = diffy(f, varargin);

    % Partial derivative of a CHEBFUN3 object in the third variable.
    out = diffz(f, varargin);

    % Gradient of a CHEBFUN3.
    varargout = grad(f);

    % Gradient of a CHEBFUN3.
    varargout = gradient(f);

    % Normal vector of a CHEBFUN3.
    varargout = normal(f, varargin);

    % Determinant of Jacobian of three CHEBFUN3 objects.
    f = jacobian(f, g, h);

    % Laplacian of a CHEBFUN3.
    out = lap(f);

    % Laplacian of a CHEBFUN3.
    out = laplacian(f);

    % Biharmonic operator of a CHEBFUN3.
    out = biharm(f);

    % Biharmonic operator of a CHEBFUN3.
    out = biharmonic(f);

    % Indefinite integral of a CHEBFUN3 in one variable.
    out = cumsum(f, varargin);

    % Indefinite integral of a CHEBFUN3 in two variables.
    out = cumsum2(f, varargin);

    % Indefinite integral of a CHEBFUN3 in all variables.
    out = cumsum3(f);
end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% HIDDEN METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Hidden = true, Static = false )

        % Test if two CHEBFUN3 objects have the same domain.
        out = domainCheck(f, g);

    end
end


text:
function value = r4_tiny ( )

%*****************************************************************************80
%
%% R4_TINY returns the smallest positive R4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VALUE, a "tiny" value.
%
  value = realmin;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function isargcoord(varargin)
%ISARGCOORD throws an error if not a point in an euclidian coordinate system
%
%   Usage: isargcoord(args)
%
%   Input parameters:
%       args        - list of args
%
%   ISARGCOORD(args) tests if all given args are a point in an euclidian
%   coordinate system. This means they have to be a vector with a length of 1,2
%   or 3. Returns an error otherwise.
%
%   See also: isargvector

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% 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.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking for vector =============================================
for ii = 1:nargin
    if ~isnumeric(varargin{ii}) || ~isvector(varargin{ii}) || ...
        length(varargin{ii})>3
        error('%s need to be a vector with length <4.',inputname(ii));
    end
end

author: sfstoolbox

repo: sfs-matlab


text:
function y = cvx_constant( x )
y = reshape( full( x.basis_( 1, : ) ), x.size_ );
if cvx_use_sparse( y ), y = sparse( y ); end

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

author: yu-jiang

repo: radpbook

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


text:
function calcAllFusedTextures_batchHN(pathData,pathText,namePT,nameCT,nameROI,outcomes,featType,CTweight_mat,scale_mat,algo_cell,Ng_mat,nBatch,matlabPATH)
% -------------------------------------------------------------------------
% function calcAllFusedTextures_batchHN(pathData,pathText,namePT,nameCT,nameROI,outcomes,featType,CTweight_mat,scale_mat,Ng_mat,nBatch,matlabPATH)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes FUSED texture features for all patients, for all
% different combinations of the following texture extraction parameters:
% - CT weight: Weight given to MRI wavelet low-pass sub-bands in the
%               PET/CT fusion process.
% - Scale: Resolution at which the ROI is isotropically resampled.
% - Ng: Number of gray-levels in the quantization process.
%
% Different extraction parameters are passed as arrays or cells in the
% function in order to test all possible combinations. This function is
% used for FUSED scans specifically. See Ref. [1,2] and 'prepareVolume.m'
% for more details.
%
% Texture features are computed for all head and neck (HN) DICOM  imaging
% data downloaded from The Cancer Imaging Archive (TCIA) website at:
% <http://dx.doi.org/10.7937/K9/xxxxxxxxxxxxxxxxxxx, and first organized
% in a 'DATA' directory using the function readAllDICOM_HN.m.  Results are
% then saved in a folder 'TEXTURES' in the HN WORKSPACE.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% [2] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathData: Full path to the HN sData files directory.
%              --> Ex: '/myProject/WORKSPACE/DATA'
% 2. pathText: Full path to the HN non texture features directory.
%              --> Ex: '/myProject/WORKSPACE/FEATURES/TEXTURES'
% 3. namePT: Cell of strings of all PET sData files to read
%            --> Ex: {'HGJ_001_PT.PTscan.mat';'HGJ_022_PT.PTscan.mat'}
% 4. namePT: Cell of strings of all CT sData files to read
%            --> Ex: {'HGJ_001_CT.CTscan.mat';'HGJ_022_CT.CTscan.mat'}
% 5. nameROI: Cell of strings specifying the ROI names to analyze for the
%             patients defined by "namePT" and "nameCT"
%             --> Ex: {'GTV';'GTV-P'}
% 6. outcomes: Structure specifying the status (1 or 0) for different
%              outcomes in HN cancer. Contains: outcomes.Failure,
%              outcomes.Locoregional, outcomes.Distant. See ref.[1] for
%              more details.
% 7. featType: Either 'GTVp' for primary GTV, or 'GTVtot' for primaty GTV +
%              nodal GTVs
%              --> Ex: 'GTVp'
% 8. CTweight_mat: Array vector specifying the different CT weights to test
%                  --> Ex: [1/4,1/3,1/2,2/3,3/4]
% 9. scale_mat: Array vector specifying the different 'Scale' values to test.
%               --> Ex: [1,2,3,4,5]
% 10. Ng_mat: Array vector specifying the different 'Ng' values to test.
%            --> Ex: [8,16,32,64]
% 11. nBatch: Number of parallel batch.
%             --> Ex: 8
% 12.  matlabPATH: Full path to the MATLAB excutable on the system.
%      --> Ex: 'matlab' (symbolic link)
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: March 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
ind = strfind(namePT{1},'_'); cohortID = namePT{1}(1:ind-1);
cd(pathText), mkdir(['batchLog_',cohortID,'_',featType,'_FusText']), cd(['batchLog_',cohortID,'_',featType,'_FusText']), pathBatch = pwd;
time = 60; % Number of seconds to wait before checking if parallel computations are done
nameOutcomes = fieldnames(outcomes); nOutcomes = numel(nameOutcomes);
scans = {'PTCT'}; nScans = numel(scans);

% PRODUCE BATCH COMPUTATIONS
nPatient = numel(namePT); valid = ones(nPatient,1);
for i = 1:nPatient
    if isempty(nameROI{i})
        valid(i) = 0;
    end
end
indValid = find(valid); nPatient = numel(indValid);
if nPatient < nBatch
    nBatch = nPatient;
end
[patients] = batchPatients(nPatient,nBatch);
save('workspace','pathData','pathText','namePT','nameCT','nameROI','patients','indValid','featType','CTweight_mat','scale_mat','algo_cell','Ng_mat'), pause(5);
for i = 1:nBatch
    nameScript = ['batch',num2str(i),'_script.m'];
    fid = fopen(nameScript,'w');
    fprintf(fid,'load(''workspace'')\n');
    fprintf(fid,['calcAllFusedTextures_HN(pathData,pathText,namePT(indValid(patients{',num2str(i),'})),nameCT(indValid(patients{',num2str(i),'})),nameROI(indValid(patients{',num2str(i),'})),featType,CTweight_mat,scale_mat,algo_cell,Ng_mat)\n']);
    fprintf(fid,['system(''touch batch',num2str(i),'_end'');\n']);
    fprintf(fid,'clear all');
    fclose(fid);
    system([matlabPATH,' -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
end

% WAITING LOOP
waitBatch(pathBatch,time,nBatch)
delete('workspace.mat')

% GROUPING RESULTS FROM ALL BATCH
nPatient = numel(namePT);
for scan = 1:nScans
    cd(pathText)
    if exist(['HGJ_001_',scans{scan},'_',featType,'_text.mat'],'file')
        temp = load(['HGJ_001_',scans{scan},'_',featType,'_text']); temp = struct2cell(temp); temp = temp{1}; % In order to get the necessary 'nameType' and 'nameFeature' fields
    else
        temp = load(['HGJ_001_',scans{scan},'_','GTVp','_text']); temp = struct2cell(temp); temp = temp{1}; % In order to get the necessary 'nameType' and 'nameFeature' fields
    end
    nameType = fieldnames(temp.Experiment1); nameType(end) = []; nType = numel(nameType); % All texture types are the same
    text = cell(numel(CTweight_mat),numel(scale_mat),numel(algo_cell),numel(Ng_mat));
    tempText = cell(1,nPatient); % Cell used to load patient textures only once
    for p = 1:nPatient
        ind = strfind(namePT{p},'_'); namePatient = namePT{p}(1:ind(2)-1);
        if exist([namePatient,'_',scans{scan},'_',featType,'_text.mat'],'file')
            load([namePatient,'_',scans{scan},'_',featType,'_text.mat']) % Variable 'textures' is now in MATLAB workspace
        else
            load([namePatient,'_',scans{scan},'_','GTVp','_text.mat']) % Variable 'textures' is now in MATLAB workspace
        end
        tempText{p} = textures;
    end
    experiment = 0;
    for c = 1:numel(CTweight_mat)
        for s = 1:numel(scale_mat)
            for a = 1:numel(algo_cell)
                for n = 1:numel(Ng_mat)
                    text{c,s,a,n} = struct;
                    experiment = experiment + 1;
                    strExperiment = ['Experiment',num2str(experiment)];
                    for t = 1:nType
                        nameFeature = fieldnames(temp.(strExperiment).(nameType{t})); nFeature = numel(nameFeature);
                        for f = 1:nFeature
                            data = zeros(nPatient,1);
                            for p = 1:nPatient
                                data(p,1) = tempText{p}.(strExperiment).(nameType{t}).(nameFeature{f});
                            end
                            text{c,s,a,n}.(nameType{t}).(nameFeature{f}).Data = data;
                            for o = 1:nOutcomes
                                [text{c,s,a,n}.(nameType{t}).(nameFeature{f}).Spearman.(nameOutcomes{o}).rs,text{c,s,a,n}.(nameType{t}).(nameFeature{f}).Spearman.(nameOutcomes{o}).p] = corr(data,outcomes.(nameOutcomes{o}),'type','Spearman');
                            end
                        end
                    end
                end
            end
        end
    end
    cd .., save(['text_',cohortID,'_',scans{scan},'_',featType],'text')
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function err = fitFunction(var,funName,params,freeList,origVarargin)
%err = fitFunction(var,funName,params,freeList,origVarargin)
%
%Support function for 'fit.m'
%Written by G.M Boynton

%stick values of var into params

params = var2params(var,params,freeList);

%evaluate the function

evalStr = sprintf('err = %s(params',funName);
for i=1:length(origVarargin)
  evalStr= [evalStr,',origVarargin{',num2str(i),'}'];
end
evalStr = [evalStr,');'];
eval(evalStr);

author: Sable

repo: mcbench-benchmarks


text:
function test_approx_test12 ( )

%*****************************************************************************80
%
%% TEST_APPROX_TEST12 plots a Bernstein spline approximant for problem 7.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2012
%
%  Author:
%
%    John Burkardt
%
  approx_filename = 'test12_approx.txt';
  data_filename = 'test12_data.txt';
  nplot = 101;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_APPROX_TEST12\n' );
  fprintf ( 1, '  Plot a Bernstein approximant for problem 5.\n' );
  fprintf ( 1, '  Note that the Bernstein approximant requires equally\n' );
  fprintf ( 1, '  spaced data.\n' );

  prob = 5;
%
%  Get the data.
%
  data_num = p00_data_num ( prob );

  [ xdata, ydata ] = p00_dat ( prob, data_num );

  r8vec2_write ( data_filename, data_num, xdata, ydata )

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data values stored in "%s"\n', data_filename );
%
%  Evaluate the approximant function.
%
  a = xdata(1);
  b = xdata(data_num);

  xplot = zeros ( nplot, 1 );
  yplot = zeros ( nplot, 1 );

  for plot = 1 : nplot

    xval = ( ( nplot - plot     ) * a     ...
           + (         plot - 1 ) * b ) ...
           / ( nplot        - 1 );

    yval = bpab_approx ( data_num - 1, a, b, ydata, xval );

    xplot(plot) = xval;
    yplot(plot) = yval;

  end

  r8vec2_write ( approx_filename, nplot, xplot, yplot );

  fprintf ( 1, '  Approximant values stored in "%s"\n', approx_filename );

  return
end


text:
function [state,options] = psoplotbestf(options,state,flag)
% Plots the best, mean, and worst scores of particle swarm.

notinf = isfinite(state.Score) ;

if strcmp(flag,'init')
    set(gca,'NextPlot','add',...
        'XLabel',xlabel('Generation'),...
        'YLabel',ylabel('Score'))
%     line(state.Generation,max(state.Score),...
%         'Color','red',...
%         'Tag','Worst Scores',...
%         'Marker','.',...
%         'LineStyle','none')
    line(state.Generation,mean(state.Score(notinf)),...
        'Color','blue',...
        'Tag','Mean Scores',...
        'Marker','.',...
        'LineStyle','none')
    line(state.Generation,min(state.Score),...
        'Color','black',...
        'Tag','Best Scores',...
        'Marker','.',...
        'LineStyle','none')
elseif strcmp(flag,'done')
    legend({'Mean Score','Best Score'})
else
%     hworst = findobj(gca,'Tag','Worst Scores','Type','line') ;
    hmean = findobj(gca,'Tag','Mean Scores','Type','line') ;
    hbest = findobj(gca,'Tag','Best Scores','Type','line') ;
    x = [get(hmean,'XData'), state.Generation] ;
%     yworst = [get(hworst,'YData'), max(state.Score)] ;
    ymean = [get(hmean,'YData'), mean(state.Score(notinf))] ;
    ybest = [get(hbest,'YData'), min(state.Score)] ;
%     set(hworst,...
%         'XData',x,...
%         'YData',yworst)
    set(hmean,...
        'XData',x,...
        'YData',ymean)
    set(hbest,...
        'XData',x,...
        'YData',ybest)
    titletxt = sprintf('Best: %g Mean: %g',ybest(end),ymean(end)) ;
    title(titletxt)
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function mhd = read_itab_mhd(filename)

fid = fopen(filename, 'rb');

%  Name of structure
mhd.stname = fread(fid, [1 10], 'uint8=>char');           %  Header identifier  (VP_BIOMAG)
mhd.stver = fread(fid, [1 8], 'uint8=>char');             %  Header version
mhd.stendian = fread(fid, [1 4], 'uint8=>char');          %  LE (little endian) or BE (big endian) format

%  Subject's INFOs
mhd.first_name = fread(fid, [1 32], 'uint8=>char');       %  Subject's first name
mhd.last_name = fread(fid, [1 32], 'uint8=>char');        %  Subject's family name
mhd.id = fread(fid, [1 32], 'uint8=>char');               %  Subject's id
mhd.notes = fread(fid, [1 256], 'uint8=>char');           %  Notes on measurement

%  Other subj infos
mhd.subj_info.sex         = fread(fid, [1 1], 'uint8=>char');  %  Sex (M or F)
pad = fread(fid, [1 5], 'uint8');
mhd.subj_info.notes       = fread(fid, [1 256], 'uint8=>char'); %  Notes on subject
mhd.subj_info.height      = fread(fid, [1 1], 'float');  %  Height in cm
mhd.subj_info.weight      = fread(fid, [1 1], 'float');  %  Weight in kg
mhd.subj_info.birthday    = fread(fid, [1 1], 'int32'); %  Birtday (1-31)
mhd.subj_info.birthmonth  = fread(fid, [1 1], 'int32'); %  Birthmonth (1-12)
mhd.subj_info.birthyear   = fread(fid, [1 1], 'int32');  %  Birthyear (1900-2002)

%  Data acquisition INFOs
mhd.time = fread(fid, [1 12], 'uint8=>char');             %  time (ascii)
mhd.date = fread(fid, [1 16], 'uint8=>char');             %  date (ascii)

mhd.nchan     = fread(fid, [1 1], 'int32');          %  total number of channels
mhd.nelech    = fread(fid, [1 1], 'int32');          %  number of electric channels
mhd.nelerefch = fread(fid, [1 1], 'int32');          %  number of electric reference channels
mhd.nmagch    = fread(fid, [1 1], 'int32');          %  number of magnetic channels
mhd.nmagrefch = fread(fid, [1 1], 'int32');          %  number of magnetic reference channels
mhd.nauxch    = fread(fid, [1 1], 'int32');          %  number of auxiliary channels
mhd.nparamch  = fread(fid, [1 1], 'int32');          %  number of parameter channels
mhd.ndigitch  = fread(fid, [1 1], 'int32');          %  number of digit channels
mhd.nflagch   = fread(fid, [1 1], 'int32');          %  number of flag channels

mhd.data_type   = fread(fid, [1 1], 'int32');        %  0 - BE_SHORT      (HP-PA, big endian)
%  1 - BE_LONG       (HP-PA, big endian)
%  2 - BE_FLOAT      (HP-PA, big endian)
%  3 - LE_SHORT      (Intel, little endian)
%  4 - LE_LONG       (Intel, little endian)
%  5 - LE_FLOAT      (Intel, little endian)
%  6 - RTE_A_SHORT   (HP-A900, big endian)
%  7 - RTE_A_FLOAT   (HP-A900, big endian)
%  8 - ASCII

mhd.smpfq       = fread(fid, [1 1], 'single');        %  sampling frequency in Hz
mhd.hw_low_fr   = fread(fid, [1 1], 'single');        %  hw data acquisition low pass filter
mhd.hw_hig_fr   = fread(fid, [1 1], 'single');        %  hw data acquisition high pass filter
mhd.hw_comb     = fread(fid, [1 1], 'int32');        %  hw data acquisition 50 Hz filter  (1-TRUE)
mhd.sw_hig_tc   = fread(fid, [1 1], 'single');        %  sw data acquisition high pass time constant
mhd.compensation= fread(fid, [1 1], 'int32');        %  0 - no compensation    1 - compensation
mhd.ntpdata     = fread(fid, [1 1], 'int32');        %  total number of time points in data

mhd.no_segments    = fread(fid, [1 1], 'int32');     %  Number of segments described in the segment structure

%  INFOs on different segments
for i=1:5
  mhd.sgmt(i).start = fread(fid, [1 1], 'int32');       %  Starting time point from beginning of data
  mhd.sgmt(i).ntptot = fread(fid, [1 1], 'int32');      %  Total number of time points
  mhd.sgmt(i).type = fread(fid, [1 1], 'int32');
  mhd.sgmt(i).no_samples = fread(fid, [1 1], 'int32');
  mhd.sgmt(i).st_sample = fread(fid, [1 1], 'int32');
end

mhd.nsmpl          = fread(fid, [1 1], 'int32');     %  Overall number of samples

%  INFOs on different samples
for i=1:4096
  mhd.smpl(i).start = fread(fid, [1 1], 'int32');       %  Starting time point from beginning of data
  mhd.smpl(i).ntptot = fread(fid, [1 1], 'int32');      %  Total number of time points
  mhd.smpl(i).ntppre = fread(fid, [1 1], 'int32');      %  Number of points in pretrigger
  mhd.smpl(i).type = fread(fid, [1 1], 'int32');
  mhd.smpl(i).quality = fread(fid, [1 1], 'int32');
end

mhd.nrefchan        = fread(fid, [1 1], 'int32');    %  number of reference channels
mhd.ref_ch = fread(fid, [1 640], 'int32');          %  reference channel list
mhd.ntpref          = fread(fid, [1 1], 'int32');    %  total number of time points in reference

%  Header INFOs
mhd.raw_header_type   = fread(fid, [1 1], 'int32');  %  0  - Unknown header
%  2  - rawfile (A900)
%  3  - GE runfile header
%  31 - ATB runfile header version 1.0
%  41 - IFN runfile header version 1.0
%  51 - BMDSys runfile header version 1.0

mhd.header_type     = fread(fid, [1 1], 'int32');    %  0  - Unknown header
%  2  - rawfile (A900)
%  3  - GE runfile header
%  4  - old header
%  10 - 256ch normal header
%  11 - 256ch master header
%  20 - 640ch normal header
%  21 - 640ch master header
%  31 - ATB runfile header version 1.0
%  41 - IFN runfile header version 1.0
%  51 - BMDSys runfile header version 1.0

mhd.conf_file = fread(fid, [1 64], 'uint8=>char');        %  Filename used for data acquisition configuration

mhd.header_size      = fread(fid, [1 1], 'int32');   %  sizeof(header) at the time of file creation
mhd.start_reference  = fread(fid, [1 1], 'int32');   %  start reference
mhd.start_data       = fread(fid, [1 1], 'int32');   %  start data

mhd.rawfile          = fread(fid, [1 1], 'int32');   %  0 - not a rawfile      1 - rawfile
mhd.multiplexed_data = fread(fid, [1 1], 'int32');   %  0 - FALSE              1 - TRUE

mhd.isns        = fread(fid, [1 1], 'int32');     %  sensor code   1 - Single channel
%               28 - Original Rome 28 ch.
%               29 -  ..............
%               .. -  ..............
%               45 - Updated Rome 28 ch. (spring 2009)
%               .. -  ..............
%               .. -  ..............
%               55 - Original Chieti 55 ch. flat
%              153 - Original Chieti 153 ch. helmet
%              154 - Chieti 153 ch. helmet from Jan 2002

%  Channel's INFOs
for i=1:640
  mhd.ch(i).type = fread(fid, [1 1], 'uint8');
  pad = fread(fid, [1 3], 'uint8');
  %  type      0 - unknown
  %            1 - ele
  %            2 - mag
  %            4 - ele ref
  %            8 - mag ref
  %           16 - aux
  %           32 - param
  %           64 - digit
  %          128 - flag
  mhd.ch(i).number = fread(fid, [1 1], 'int32');         %  number
  mhd.ch(i).label = fixstr(fread(fid, [1 16], 'uint8=>char'));      %  label
  mhd.ch(i).flag = fread(fid, [1 1], 'uint8');
  pad = fread(fid, [1 3], 'uint8');
  %  on/off flag 0 - working channel
  %              1 - noisy channel
  %              2 - very noisy channel
  %              3 - broken channel
  mhd.ch(i).amvbit = fread(fid, [1 1], 'float');         %  calibration from LSB to mV
  mhd.ch(i).calib = fread(fid, [1 1], 'float');          %  calibration from mV to unit
  mhd.ch(i).unit = fread(fid, [1 6], 'uint8=>char');  %  unit label  (fT, uV, ...)
  pad = fread(fid, [1 2], 'uint8');

  %  Calibration to Tesla for (MEG, MEG REF) or Volt for (EEG, EEG REF)
  mhd.ch(i).unit = deblank(mhd.ch(i).unit);
  switch (mhd.ch(i).unit)
      case 'fT'
          mhd.ch(i).calib = mhd.ch(i).calib / 1.E-15;
      case 'pT'
          mhd.ch(i).calib = mhd.ch(i).calib / 1.E-12;
      case 'uV'
          mhd.ch(i).calib = mhd.ch(i).calib / 1.E-6;
      case 'mV'
          mhd.ch(i).calib = mhd.ch(i).calib / 1.E-3;
  end

  mhd.ch(i).ncoils = fread(fid, [1 1], 'int32');     %  number of coils building up one channel
  mhd.ch(i).wgt = fread(fid, [1 10], 'float');       %  weight of coils
  %  position and orientation of coils
  for j=1:10
    mhd.ch(i).position(j).r_s = fread(fid, [1 3], 'float');
    mhd.ch(i).position(j).u_s = fread(fid, [1 3], 'float');
  end
end

%  Sensor position INFOs

mhd.r_center = fread(fid, [1 3], 'float');          %  sensor position in convenient format
mhd.u_center = fread(fid, [1 3], 'float');
mhd.th_center      = fread(fid, [1 1], 'float');  %  sensor orientation as from markers fit
mhd.fi_center= fread(fid, [1 1], 'float');
mhd.rotation_angle= fread(fid, [1 1], 'float');
mhd.cosdir = fread(fid, [3 3], 'float');          %  for compatibility only
mhd.irefsys      = fread(fid, [1 1], 'int32');    %  reference system   0 - sensor reference system
%                     1 - Polhemus
%                     2 - head3
%                     3 - MEG

%  Marker positions for MRI integration
mhd.num_markers  = fread(fid, [1 1], 'int32');  %  Total number of markers
mhd.i_coil = fread(fid, [1 64], 'int32');       %  Markers to be used to find sensor position
mhd.marker = fread(fid, [3 64], 'float');       %  Position of all the markers  -  MODIFIED VP
mhd.best_chi     = fread(fid, [1 1], 'float');  %  Best chi_square value obtained in finding sensor position

mhd.cup_vertex_center  = fread(fid, [1 1], 'float');    %  dist anc sensor  cente vertex  (as entered
%  from keyboard)
mhd.cup_fi_center  = fread(fid, [1 1], 'float');        %  fi angle of sensor center
mhd.cup_rotation_angle  = fread(fid, [1 1], 'float');   %  rotation angle of sensor center axis

mhd.dist_a1_a2  = fread(fid, [1 1], 'float');           %  head informations
%  (used to find subject's head dimensions)
mhd.dist_inion_nasion = fread(fid, [1 1], 'float');
mhd.max_circ = fread(fid, [1 1], 'float');
mhd.nasion_vertex_inion = fread(fid, [1 1], 'float');

%  Data analysis INFOs
mhd.security  = fread(fid, [1 1], 'int32');          %  security flag
mhd.ave_alignement  = fread(fid, [1 1], 'int32');    %  average data alignement 0 - FALSE
%                          1 - TRUE

mhd.itri  = fread(fid, [1 1], 'int32');              %  trigger channel number
mhd.ntpch  = fread(fid, [1 1], 'int32');             %  no. of time points per channel
mhd.ntppre  = fread(fid, [1 1], 'int32');            %  no. of time points of pretrigger
mhd.navrg  = fread(fid, [1 1], 'int32');             %  no. of averages
mhd.nover  = fread(fid, [1 1], 'int32');             %  no. of discarded averages
mhd.nave_filt  = fread(fid, [1 1], 'int32');         %  no. of applied filters

%  Filters used before average
for i=1:15
  mhd.ave_filt(i).type = fread(fid, [1 1], 'int32');                 % type   0 - no filter
  %        1 - bandpass  - param[0]: highpass freq
  %                      - param[1]: lowpass freq
  %        2 - notch     - param[0]: notch freq 1
  %                        param[1]: notch freq 2
  %                        param[2]: notch freq 3
  %                        param[3]: notch freq 4
  %                        param[4]: span
  %        3 - artifact  - param[0]: True/False
  %        4 - adaptive  - param[0]: True/False
  %        5 - rectifier - param[0]: True/False
  %        6 - heart     - param[0]: True/False
  %        7 - evoked    - param[0]: True/False
  %        8 - derivate  - param[0]: True/False
  %        9 - polarity  - param[0]: True/False
  mhd.ave_filt(i).param = fread(fid, [1 5], 'float');             % up to 5 filter parameters
end

mhd.stdev  = fread(fid, [1 1], 'int32');             %  0 - not present
mhd.bas_start  = fread(fid, [1 1], 'int32');         %  starting data points for baseline
mhd.bas_end  = fread(fid, [1 1], 'int32');           %  ending data points for baseline

mhd.source_files  = fread(fid, [1 32], 'int32');  %  Progressive number of files (if more than one)

%  template INFOs
mhd.ichtpl = fread(fid, [1 1], 'int32');
mhd.ntptpl = fread(fid, [1 1], 'int32');
mhd.ifitpl = fread(fid, [1 1], 'int32');
mhd.corlim = fread(fid, [1 1], 'float');

%  Filters used before template
for i=1:15
  mhd.tpl_filt(i).type = fread(fid, [1 1], 'int32');                 % type   0 - no filter
  %        1 - bandpass  - param[0]: highpass freq
  %                      - param[1]: lowpass freq
  %        2 - notch     - param[0]: notch freq 1
  %                        param[1]: notch freq 2
  %                        param[2]: notch freq 3
  %                        param[3]: notch freq 4
  %                        param[4]: span
  %        3 - artifact  - param[0]: True/False
  %        4 - adaptive  - param[0]: True/False
  %        5 - rectifier - param[0]: True/False
  %        6 - heart     - param[0]: True/False
  %        7 - evoked    - param[0]: True/False
  %        8 - derivate  - param[0]: True/False
  %        9 - polarity  - param[0]: True/False
  mhd.tpl_filt(i).param = fread(fid, [1 5], 'float');             % up to 5 filter parameters
end

%  Just in case info
mhd.dummy = fread(fid, [1 64], 'int32');
% there seems to be more dummy data at the end...

fclose(fid);

function str = fixstr(str)
sel = find(str==0, 1, 'first');
if ~isempty(sel)
  str = str(1:sel-1);
end

author: brainstorm-tools


text:
function [residuals,W] = getresiduals (X,T,S,maxorder,order,orderoffset,...
    timelag,exptimelag,zeromean,W)
%
% Compute residuals - useful when using a global model to remove global trends
%
% INPUT
% X             time series
% T             length of series
%
%
% OUTPUT
% residuals
% W
%
% Author: Diego Vidaurre, OHBA, University of Oxford

if nargin<10, W = []; end

N = length(T);

if ~isempty(S) && any(S(:)==0)
    [~,order] = formorders(order,orderoffset,timelag,exptimelag);
    if isempty(W)
        [W,~,~,residuals] = mlmar(X,maxorder,T,double(S==0),order,orderoffset,timelag,...
        exptimelag,zeromean);
    else
        [~,~,~,residuals] = mlmar(X,maxorder,T,double(S==0),order,orderoffset,timelag,...
            exptimelag,zeromean,W);
    end
else
    W = [];
    residuals = zeros(sum(T)-length(T)*maxorder,size(X,2));
    for in=1:N
        t0 = sum(T(1:in-1));
        t = t0 - (in-1)*maxorder;
        residuals(t+1:t+T(in)-maxorder,:) = X(t0+maxorder+1:t0+T(in),:);
    end
end

author: OHBA-analysis

repo: HMM-MAR


text:
function a = triu(a,k)
%TRIU         Implements  triu(a,k)  for hessians
%
%   c = triu(a,k)
%
% functionality as Matlab function triu for matrices
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if nargin==1
    k = 0;
  end

  a.x = triu(a.x,k);
  index = ( triu( ones(size(a.x)) , k ) == 0 );
  a.dx(:,index) = 0;
  a.hx(:,index) = 0;

repo: OpenMAS


text:
function value = atan4 ( y, x )

%*****************************************************************************80
%
%% ATAN4 computes the inverse tangent of the ratio Y / X.
%
%  Discussion:
%
%    ATAN4 returns an angle whose tangent is ( Y / X ), a job which
%    the built in functions ATAN and ATAN2 already do.
%
%    However:
%
%    * ATAN4 always returns a positive angle, between 0 and 2 PI,
%      while ATAN and ATAN2 return angles in the interval [-PI/2,+PI/2]
%      and [-PI,+PI] respectively;
%
%    * ATAN4 accounts for the signs of X and Y, (as does ATAN2).  The ATAN
%     function by contrast always returns an angle in the first or fourth
%     quadrants.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real Y, X, two quantities which represent the tangent of
%    an angle.  If Y is not zero, then the tangent is (Y/X).
%
%    Output, real ATAN4, an angle between 0 and 2 * PI, whose tangent is
%    (Y/X), and which lies in the appropriate quadrant so that the signs
%    of its cosine and sine match those of X and Y.
%

%
%  Special cases:
%
  if ( x == 0.0E+00 )

    if ( 0.0E+00 < y )
      value = pi / 2.0E+00;
    elseif ( y < 0.0E+00 )
      value = 3.0E+00 * pi / 2.0E+00;
    elseif ( y == 0.0E+00 )
      value = 0.0E+00;
    end

  elseif ( y == 0.0E+00 )

    if ( 0.0E+00 < x )
      value = 0.0E+00;
    elseif ( x < 0.0E+00 )
      value = pi;
    end
%
%  We assume that ATAN2 is correct when both arguments are positive.
%
  else

    abs_y = abs ( y );
    abs_x = abs ( x );

    theta_0 = atan2 ( abs_y, abs_x );

    if ( 0.0E+00 < x & 0.0E+00 < y )
      value = theta_0;
    elseif ( x < 0.0E+00 & 0.0E+00 < y )
      value = pi - theta_0;
    elseif ( x < 0.0E+00 & y < 0.0E+00 )
      value = pi + theta_0;
    elseif ( 0.0E+00 < x & y < 0.0E+00 )
      value = 2.0E+00 * pi - theta_0;
    end

  end

  return
end


%  Script file: ode45_test1.m
%
%  Purpose:
%    This program solves a differential equation of the
%    form dy/dt + 2 * y = 0, with the initial condition
%    y(0) = 1.
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    03/15/07    S. J. Chapman        Original code
%
% Define variables:
%   odefun_handle -- Handle to function that defines the derivative
%   tspan         -- Duration to solve equation for
%   yo            -- Initial condition for equation
%   t             -- Array of solution times
%   y             -- Array of solution values

% Get a handle to the function that defines the
% derivative.
odefun_handle = @fun1;

% Solve the equation over the period 0 to 5 seconds
tspan = [0 5];

% Set the initial conditions
y0 = 1;

% Call the differential equation solver.
[t,y] = ode45(odefun_handle,tspan,y0);

% Plot the result
figure(1);
plot(t,y,'b-','LineWidth',2);
grid on;
title('\bfSolution of Differential Equation');
xlabel('\bfTime (s)');
ylabel('\bf\ity''');


text:
function test_suite = test_angle2Points
%TESTCLIPLINE  One-line description here, please.
%   output = testAngle2Points(input)
%
%   Example
%   testAngle2Points
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2009-04-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testHoriz(testCase) %#ok<*DEFNU>
% all points inside window, possibly touching edges

p1 = [0 0];
p2 = [10 0];
angle = angle2Points(p1, p2);
testCase.assertEqual(0, angle, 'AbsTol', .01);

angle = angle2Points(p2, p1);
testCase.assertEqual(pi, angle, 'AbsTol', .01);

function testVert(testCase)
% all points inside window, possibly touching edges

p1 = [0 0];
p2 = [0 10];
angle = angle2Points(p1, p2);
testCase.assertEqual(pi/2, angle, 'AbsTol', .01);

angle = angle2Points(p2, p1);
testCase.assertEqual(3*pi/2, angle, 'AbsTol', .01);

function testMultiMulti(testCase)
% all points inside window, possibly touching edges

p1 = [0 0;0 0;0 0;0 0];
p2 = [10 0;10 10;0 10;-10 10];
angle = angle2Points(p1, p2);
testCase.assertEqual(size(p1, 1), size(angle, 1), 'AbsTol', .01);

res = [0;pi/4;pi/2;3*pi/4];
testCase.assertEqual(res, angle, 'AbsTol', .01);

function testOneMulti(testCase)
% all points inside window, possibly touching edges

p1 = [0 0];
p2 = [10 0;10 10;0 10;-10 10];
angle = angle2Points(p1, p2);
testCase.assertEqual(size(p2, 1), size(angle, 1), 'AbsTol', .01);

res = [0;pi/4;pi/2;3*pi/4];
testCase.assertEqual(res, angle, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
classdef ObjectnessBING < handle
    %OBJECTNESSBING  The Binarized normed gradients algorithm for Objectness
    %
    % Implementation of BING for Objectness.
    %
    % ## Saliency API
    %
    % Many computer vision applications may benefit from understanding where
    % humans focus given a scene. Other than cognitively understanding the way
    % human perceive images and scenes, finding salient regions and objects in
    % the images helps various tasks such as speeding up object detection,
    % object recognition, object tracking and content-aware image editing.
    %
    % About the saliency, there is a rich literature but the development is
    % very fragmented. The principal purpose of this API is to give a unique
    % interface, a unique framework for use and plug sever saliency
    % algorithms, also with very different nature and methodology, but they
    % share the same purpose, organizing algorithms into three main
    % categories:
    %
    % * Static Saliency
    % * Motion Saliency
    % * Objectness
    %
    % Saliency UML diagram:
    %
    % ![image](https://docs.opencv.org/3.3.1/saliency.png)
    %
    % To see how API works, try tracker demo: `computeSaliency_demo.m`.
    %
    % ## Objectness Algorithms
    %
    % Objectness is usually represented as a value which reflects how likely
    % an image window covers an object of any category. Algorithms belonging
    % to this category, avoid making decisions early on, by proposing a small
    % number of category-independent proposals, that are expected to cover all
    % objects in an image. Being able to perceive objects before identifying
    % them is closely related to bottom up visual attention (saliency).
    %
    % Presently, the Binarized normed gradients algorithm [BING] has been
    % implemented.
    %
    % ## References
    % [BING]:
    % > Cheng, Ming-Ming, et al. "BING: Binarized normed gradients for
    % > objectness estimation at 300fps". In IEEE CVPR, 2014.
    %
    % See also: cv.StaticSaliencySpectralResidual,
    %  cv.MotionSaliencyBinWangApr2014
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % base for window size quantization. default 2
        Base
        % Size for non-maximal suppress. default 2
        NSS
        % As described in the paper: feature window size (W, W). default 8
        W
    end

    methods
        function this = ObjectnessBING()
            %OBJECTNESSBING  Constructor, creates a specialized saliency algorithm of this type
            %
            %     obj = cv.ObjectnessBING()
            %
            % See also: cv.ObjectnessBING.computeSaliency
            %
            this.id = ObjectnessBING_(0, 'new');
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.ObjectnessBING
            %
            if isempty(this.id), return; end
            ObjectnessBING_(this.id, 'delete');
        end

        function objectnessValues = getObjectnessValues(this)
            %GETOBJECTNESSVALUES  Return the list of the rectangles' objectness value
            %
            %     objectnessValues = obj.getObjectnessValues()
            %
            % ## Output
            % * __objectnessValues__ vector of floats in the same order as
            %   `objectnessBoundingBox` returned by the algorithm (in
            %   `computeSaliency` function). The bigger value these scores
            %   are, it is more likely to be an object window.
            %
            % See also: cv.ObjectnessBING.computeSaliency
            %
            objectnessValues = ObjectnessBING_(this.id, 'getobjectnessValues');
        end

        function setTrainingPath(this, trainingPath)
            %SETTRAININGPATH  This is a utility function that allows to set the correct path from which the algorithm will load the trained model
            %
            %     obj.setTrainingPath(trainingPath)
            %
            % ## Input
            % * __trainingPath__ trained model path.
            %
            % See also: cv.ObjectnessBING.setBBResDir
            %
            ObjectnessBING_(this.id, 'setTrainingPath', trainingPath);
        end

        function setBBResDir(this, resultsDir)
            %SETBBRESDIR  This is a utility function that allows to set an arbitrary path in which the algorithm will save the optional results
            %
            %     obj.setBBResDir(resultsDir)
            %
            % ## Input
            % * __resultsDir__ results' folder path.
            %
            % Writes on file the total number and the list of rectangles
            % returned by objectness, one for each row.
            %
            % See also: cv.ObjectnessBING.setTrainingPath
            %
            ObjectnessBING_(this.id, 'setBBResDir', resultsDir);
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.ObjectnessBING.empty, cv.ObjectnessBING.load
            %
            ObjectnessBING_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if detector object is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.ObjectnessBING.clear, cv.ObjectnessBING.load
            %
            b = ObjectnessBING_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.ObjectnessBING.load
            %
            ObjectnessBING_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.ObjectnessBING.save
            %
            ObjectnessBING_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.ObjectnessBING.save, cv.ObjectnessBING.load
            %
            name = ObjectnessBING_(this.id, 'getDefaultName');
        end
    end

    %% Saliency, Objectness
    methods
        function objectnessBoundingBox = computeSaliency(this, img)
            %COMPUTESALIENCY  Compute the saliency
            %
            %     objectnessBoundingBox = obj.computeSaliency(img)
            %
            % ## Input
            % * __img__ The input image, 8-bit color.
            %
            % ## Output
            % * __objectnessBoundingBox__ objectness Bounding Box vector.
            %   According to the result given by this specialized algorithm,
            %   the `objectnessBoundingBox` is a cell array of 4-element
            %   vectors. Each bounding box is represented by
            %   `[minX, minY, maxX, maxY]`.
            %
            % Performs all the operations and calls all internal functions
            % necessary for the accomplishment of the Binarized normed
            % gradients algorithm.
            %
            % See also: cv.ObjectnessBING.ObjectnessBING
            %
            objectnessBoundingBox = ObjectnessBING_(this.id, 'computeSaliency', img);
        end
    end

    %% Getters/Setters
    methods
        function value = get.Base(this)
            value = ObjectnessBING_(this.id, 'get', 'Base');
        end
        function set.Base(this, value)
            ObjectnessBING_(this.id, 'set', 'Base', value);
        end

        function value = get.NSS(this)
            value = ObjectnessBING_(this.id, 'get', 'NSS');
        end
        function set.NSS(this, value)
            ObjectnessBING_(this.id, 'set', 'NSS', value);
        end

        function value = get.W(this)
            value = ObjectnessBING_(this.id, 'get', 'W');
        end
        function set.W(this, value)
            ObjectnessBING_(this.id, 'set', 'W', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function ct3(im,w)
%function ct3(im,w)
%
%Allows simple browsing of a CT scan using the mousewheel.  Assumes a grid
%spacing ratio of 3:1 axial thickness:radial thickness for display purposes
%
%Usage: ct3(im) will open Figures 1,2, and 3 and use imshow() to display the
%       CT dataset contained in "im" in all 3 orthogonal planes (sagittal,
%       coronal, axial).  Scrolling the mousewheel will move to the
%       next/previous slice in the active figure window.  Clicking the mouse
%       will print the current slice number of all 3 figure windows in
%       the console.
%
%input:
%   -im: 3D array of grayscale values.  This function assumes that the
%        1st and 2nd dimensions of the array contain the radial (or
%        sagittal/coronal) data, and the 3rd dimension of "im" contains
%        axial data.  This function also assumes a voxel size ratio of
%        1:1:3.
%
%   -w : (optional) window to use for display.  Default=[-140 160]
%        1x2 numerical array that contains the range of values of interest
%        to be displayed using ct3
%        note - use Matlab built-in function "imcontrast" (Image Proc Toolbox)
%        or plot the image data using hist() to easily determine the proper
%        window for your dataset
%
%output:
%        none
%
%ex:
%     ct3(im)
%     ct3(im, [-140 160])

%reverse the z direction
im=im(:,:,end:-1:1);

f1=figure(1);
f2=figure(2);
f3=figure(3);

slice1=floor(size(im,1)/2);
slice2=floor(size(im,2)/2);
slice3=floor(size(im,3)/2);

mn=1;

mx1=size(im,1);
mx2=size(im,2);
mx3=size(im,3);

if nargin < 2
    w=[-140 160];
    %w=[766 961];
end

%% INITALIZE THE FIGURE WINDOWS

figure(3)
i3=imshow(im(:,:,slice3),w, 'parent',gca);
set(gca,'position',[0 0 1 1])
hold on
p311=plot([0 mx2/10],[slice1 slice1],'g', 'linewidth',2);
p312=plot([mx2*9/10 mx2],[slice1 slice1],'g', 'linewidth',2);
p321=plot([slice2 slice2],[0 mx1/10],'g', 'linewidth',2);
p322=plot([slice2 slice2],[mx1*9/10 mx1],'g', 'linewidth',2);
hold off

figure(1)
i1=imshow(permute(im(slice1,:,:),[3 2 1]),w,'parent',gca);
set(gca,'position',[0 0 1 1])
set(gca,'DataAspectRatio',[3 1 1])
%set(gcf, 'position',(get(f3,'position')))
hold on
p131=plot([0 mx2/10],[slice3 slice3],'g', 'linewidth',2);
p132=plot([mx2*9/10 mx2],[slice3 slice3],'g', 'linewidth',2);
p121=plot([slice2 slice2],[0 mx3/10],'g', 'linewidth',2);
p122=plot([slice2 slice2],[mx3*9/10 mx3],'g', 'linewidth',2);
hold off

figure(2)
i2=imshow(permute(im(:,slice2,:),[3 1 2]),w,'parent',gca);
set(gca,'position',[0 0 1 1])
hold on
p231=plot([0 mx1/10],[slice3 slice3],'g', 'linewidth',2);
p232=plot([mx1*9/10 mx1],[slice3 slice3],'g', 'linewidth',2);
p211=plot([slice1 slice1],[0 mx3/10],'g', 'linewidth',2);
p212=plot([slice1 slice1],[mx3*9/10 mx3],'g', 'linewidth',2);
set(gca,'DataAspectRatio',[3 1 1])
%set(gcf, 'position',(get(f3,'position')))
hold off

set(f1, 'WindowScrollWheelFcn', @wheel)
set(f1, 'WindowButtonDownFcn', @click)
set(f2, 'WindowScrollWheelFcn', @wheel)
set(f2, 'WindowButtonDownFcn', @click)
set(f3, 'WindowScrollWheelFcn', @wheel)
set(f3, 'WindowButtonDownFcn', @click)

%% CAPTURE MOUSEWHEEL EVENTS
    function wheel(src, evnt)
        cur=gcf;
        switch cur
            case 1
                if evnt.VerticalScrollCount > 0
                    slice1=slice1+1;
                else
                    slice1=slice1-1;
                end
                slice1=re_eval1(slice1);
            case 2
                if evnt.VerticalScrollCount > 0
                    slice2=slice2-1;
                else
                    slice2=slice2+1;
                end
                slice2=re_eval2(slice2);
            case 3
                if evnt.VerticalScrollCount > 0
                    slice3=slice3+1;
                else
                    slice3=slice3-1;
                end
                slice3=re_eval3(slice3);
            otherwise
                %do nothing
        end
    end

%% REDRAW FIGURES
    function slice1=re_eval1(slice1)
        if slice1>mx1
            slice1=mx1;
        elseif slice1<mn
            slice1=mn;
        else
            %imshow(flipud(imresize(squeeze(im(slice1,:,:)), size(squeeze(im(slice1,:,:))).*[1 3])'),w,'parent',gca)
            set(i1,'CData',permute(im(slice1,:,:),[3 2 1]))
        end
        do_lines
    end
    function slice2=re_eval2(slice2)
        if slice2>mx2
            slice2=mx2;
        elseif slice2<mn
            slice2=mn;
        else
            %imshow(flipud(imresize(squeeze(im(:,slice2,:)), size(squeeze(im(:,slice2,:))).*[1 3])'),w,'parent',gca)
            set(i2,'CData',permute(im(:,slice2,:),[3 1 2]))
        end
        do_lines
    end
    function slice3=re_eval3(slice3)
        if slice3>mx3
            slice3=mx3;
        elseif slice3<mn
            slice3=mn;
        else
            %imshow(im(:,:,slice3),w, 'parent',gca)
            set(i3,'CData',im(:,:,slice3))
        end
        do_lines
    end
%% CAPTURE MOUSE CLICKS
    function click(src,evnt)
        %account for reversal of z dimension
        disp(num2str([slice1 slice2 (mx3+1-slice3)]))

    end

%% REDRAW MARKER LINES
    function do_lines
        set(p131,'YData',[slice3 slice3]);
        set(p121,'XData',[slice2 slice2]);
        set(p132,'YData',[slice3 slice3]);
        set(p122,'XData',[slice2 slice2]);

        set(p231,'YData',[slice3 slice3]);
        set(p211,'XData',[slice1 slice1]);
        set(p232,'YData',[slice3 slice3]);
        set(p212,'XData',[slice1 slice1]);

        set(p311,'YData',[slice1 slice1]);
        set(p321,'XData',[slice2 slice2]);
        set(p312,'YData',[slice1 slice1]);
        set(p322,'XData',[slice2 slice2]);
    end

end


text:
function matching_channels = match_nirs_labels(labels)
% MATCHING_CHANNELS finds the labels of the signals associated to the
% same combination of receiver and transmitter and provides as output a
% binary matrix where 1 represents that channels in the corresponding
% row-column combination match.
%
% OUTPUT
%  matching_channels: nChan x nChan binary matrix; where 1 represents that
% channels in the corresponding row-column combination match, and 0 represents
% the cases where they do not match.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

num_channels=numel(labels);
matching_channels = ones(num_channels, num_channels);

for n = 1 : num_channels
    label_check1 = labels{n};
    label_check1_split = split(label_check1, ' [');
    label_check1_split = label_check1_split{1};
    for m = 1 : num_channels
        if(n == m)
            continue
        end
        label_check2 = labels{m};
        label_check2_split = split(label_check2, ' [');
        label_check2_split = label_check2_split{1};
        matching_channels(n,m) = strcmp(label_check1_split, label_check2_split);
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function y = nbBernPred(model, X)
% Prediction of naive Bayes classifier with independent Bernoulli.
% input:
%   model: trained model structure
%   X: d x n data matrix
% output:
%   y: 1 x n predicted class label
% Written by Mo Chen (sth4nth@gmail.com).
mu = model.mu;
w = model.w;
[~,y] = max(log(mu)'*X+log(1-mu)'*(1-X)+log(w(:)),[],1);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
% This function generate the log Mel filterbanks features of the training
% and evaluation data. Note that the filterbanks features generated here is
% slightly different from that by other toolkits, such as Kaldi.
% This recipe will always use the filterbank features generated in the same
% way as in this function, so we will be able to do joint training of
% feature extraction and acoustic modeling later. (e.g. you cannot use
% Kaldi fbank features in joint training unless you have Matlab version of
% the Kaldi feature extraction.

function GenFbankFeatures(dataset)
addpath('lib');

% load the wavelist

chime_root = ChoosePath4OS({'D:/Data/CHiME4', '/home/xiaoxiong/CHiME4'});   % you can set two paths, first for windows OS and second for Linux OS.
fbank_root = [chime_root '/fbank/plain'];

switch dataset
    case 'tr05_orig'
        wavlist = findFiles([chime_root '/audio/isolated/tr05_org']);
    case 'tr05'

    case 'dt05'

    case 'et05'

    otherwise
        fprintf('Unknown dataset: %s!\n', dataset);
end

[layer, para] = BuildFbankExtractionNet();

for i=1:length(wavlist)
    PrintProgress(i-1, length(wavlist), 100, 'Generate filterbank');
    words = ExtractWordsFromString_v2(wavlist{i}, '/');
    wav = audioread(wavlist{i});
    wav = StoreWavInt16(wav');
    Data(1).data{1} = wav;

    [fbank, layer2] = FeatureTree2(Data, para, layer);
    fbank = fbank{1}{1};

    fbank_file = [fbank_root '/' words{end-2} '/' words{end-1} '/' words{end}(1:end-4) '.fbank'];
    fbank_dir = fileparts(fbank_file);
    my_mkdir(fbank_dir);
    writeHTK(fbank_file, fbank', 'MFCC_0', 1);

end

author: singaxiong

repo: SignalGraph


text:
function y=nsct_fuse(I1,I2,nlevels)
%    NSCT
%    Input:
%    I1 - input image A
%    I2 - input image B
%    nlevels - number of directions in each decomposition level
%    Output:
%    y  - fused image
%
%    The code is edited by Yu Liu, 01-09-2014.

%nlevels = [2,3,3,4] ;
pfilter = 'pyrexc' ;
dfilter = 'vk' ;

I1=double(I1);
I2=double(I2);
coeffs_1 = nsctdec( I1, nlevels, dfilter, pfilter );
coeffs_2 = nsctdec( I2, nlevels, dfilter, pfilter );

[m,n]=size(I1);
coeffs=coeffs_2;
for i=2:numel(nlevels)+1

    if nlevels(i-1)==0
        E1=abs(coeffs_1{i});
        E2=abs(coeffs_2{i});
        %             map=E1>E2;
         um=3;
    A1 = ordfilt2(abs(es2(E1,floor(um/2))), um*um, ones(um));
  	A2 = ordfilt2(abs(es2(E2,floor(um/2))), um*um, ones(um));
    % second step
  	map= (conv2(double(A1 > A2), ones(um), 'valid')) > floor(um*um/2);
        coeffs{i}(map)=coeffs_1{i}(map);
    else
        for j=1:(2^nlevels(i-1))
            E1=abs(coeffs_1{i}{j});
            E2=abs(coeffs_2{i}{j});
%             map=E1>E2;
    um=3;
    A1 = ordfilt2(abs(es2(E1,floor(um/2))), um*um, ones(um));
  	A2 = ordfilt2(abs(es2(E2,floor(um/2))), um*um, ones(um));
    % second step
  	map= (conv2(double(A1 > A2), ones(um), 'valid')) > floor(um*um/2);
            coeffs{i}{j}(map)=coeffs_1{i}{j}(map);
        end
    end
end

% energy based method

coeffs{1}=(coeffs_1{1}+coeffs_2{1})/2;

y= nsctrec( coeffs, dfilter, pfilter ) ;

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [Model, Info] = linear_sparse_stepwise_reg(X,Y,Model,parm)
%  Estimate linear weight matrix for input-output mapping
%     regularization version
%
% Notice !!
%     Delay embedding is done in this module,
%     then input vector should be original input data without embedding
%
%   [Model, Info] = linear_sparse_stepwise_reg(X,Y,Model,parm)
%     Scalar iteration version
% --- Input
%  Y  : Output data ( N x T x Ntrial )
%  X  : Input data  ( M x (T + (Dtau-1)*Tau)  x Ntrial)
%  N  =  # of output
%  M  =  # of input (original input space dimension)
%  T  =  # of time sample
%
%  Estimate the following model
%    Y(t) = W * [X(:, t + (Dtau-1)*Tau ); ...; X(:, t)]
%
%  Model : Structure for estimated model
%  if Model is empty, initialization is done before training
%  if Model is previous training result, re-training us done
%
%  parm  : Structure for learning parameter
%  parm.Npre_train :  # of VB-update in initial training
%  parm.Ntrain :  # of training
%  parm.Nskip  :  skip # for print
%  parm.a_min  :  Min value for pruning small variance component
%  parm.Prune  :  = 1 : Prune small variance & irrelevant input dimension
%
%  parm.Tau       = Lag time
%  parm.Dtau      = Number of embedding dimension
%    Total input dimension in embedding space is (M * parm.Dtau)
% --- Output
%  Model : Structure for estimated model
%  Model.SY  :  Noise variance         ( 1 x 1 )
%  Model.W   :  Weight matrix          ( N x M*D ) , D = parm.Dtau
%  Model.A   :  Prior weight variance  ( 1 x M*D )
%  Model.ix_act : Active index for W after pruning
%
%  Info  : Structure for learning process history
%  Info.FE  = LP + H : Free energy
%  Info.LP  = Log likelihood
%  Info.H   = - Model entropy
%
% 2009-10-18 Made by M. Sato

MINVAL = 1.0e-15;

fprintf('linear_sparse_stepwise (regularization) start\n')

% Dimension
[N ,T  ,Ntrial]  = size(Y); % N = # of output
[M ,Tx ,Ntrialx] = size(X); % M = # of input without embedding

Tall = T * Ntrial;

if Ntrial~=Ntrialx, error('# of trial is different in X and Y'); end;

Nskip  = 100;   % skip steps for display info
a_min = 1e-14;	% Minimum value for weight pruning
Fdiff = 1e-12; % Threshold for convergence
Ncheck = 100;   % Minimum number of training iteration
Fstep  = 5;     % Free energy convergence check step
Prune  = 1;     % Prune mode
space_ARD=1;	% if space_ARD=1, ARD is done only for space dimension

if isfield(parm,'Nskip'), Nskip  = parm.Nskip; end;
if isfield(parm,'Fdiff'), Fdiff   = parm.Fdiff; end;
if isfield(parm,'a_min'), a_min   = parm.a_min ; end;
if isfield(parm,'Prune'), Prune = parm.Prune; end;
if isfield(parm,'Ncheck'), Ncheck = parm.Ncheck; end
if isfield(parm,'Fstep'),Fstep  = parm.Fstep ; end
if isfield(parm,'Fstep'),Fstep  = parm.Fstep ; end
if isfield(parm,'space_ARD'),space_ARD  = parm.space_ARD ; end

% # of embedding dimension
if isfield(parm,'Dtau')
	D    = parm.Dtau;
	tau  = parm.Tau;
else
	D    = 1;
	tau  = 1;
end
% Parameters
Ntrain = parm.Ntrain;

if Ntrain < 1, Info = []; return; end;

if isfield(parm,'Nupdate')
	Nupdate = parm.Nupdate;
else
	Nupdate = 1;
end

fprintf('--- Output Dimension  = %d\n',N)
fprintf('--- Input  Dimension  = %d\n',M)
fprintf('--- Embedding  Dimension  = %d\n',D)
fprintf('--- Number of trials  = %d\n',Ntrial)
fprintf('--- Number of training sample = %d\n',T)
fprintf('--- Total update iteration    = %d \n',Ntrain)

% original input dimension
Xdim  = M;
M = Xdim*D;
M_ALL = M;

%
% --- Initialization
%

% Input/Output variance
sx = mean((X(:) - mean(X(:))).^2);
sy = mean((Y(:) - mean(Y(:))).^2);

SY0 = mean(sy);
A0  = 1./mean(sx);

if isfield(parm,'Ta0') && parm.Ta0 > 0,
	Ta0 = parm.Ta0;
	a0  = parm.a0 * A0;
else
	Ta0 = 0;
	a0  = 1;
end

% Input covariance
XX  = sum(sum(X.^2,3),2)/(Tx*Ntrial);
XX  = repmat(XX', [1 D]);% 1 x M

if isfield(Model,'ix_act')
	SY  = mean(Model.SY);  % 1 x 1

	% Active index
	ix_act = Model.ix_act;
	A = zeros(1,M_ALL);
	W = zeros(N,M_ALL);

	A(ix_act)   = sum(Model.A,1);
	W(:,ix_act) = Model.W;

	if M_ALL ~= Xdim*D,
		fprintf('M_ALL=%d,Xdim=%d,D=%d\n',M_ALL ,Xdim ,D)
		error('M_ALL ~= Xdim*D')
	end

	% Active index in input space without embedding
	IX_dim = find( sum(reshape(A,[Xdim,D]),2) > 0);

	IX_act = repmat( (0:(D-1))* Xdim ,[length(IX_dim) 1]) ...
	       + repmat(IX_dim, [1 D]);
	IX_act = IX_act(:);

	W  = W(:,IX_act) ;  % N x (M*D)
	A  = A(IX_act) ;  % 1 x (M*D)
	X  = X(IX_dim,:,:);
	XX = XX(IX_act);

	M     = length(IX_act);
	Xdim  = length(IX_dim);
elseif isfield(Model,'A')
	A = Model.A;
	if size(A,1)==1 && size(A,2)==1
		A = repmat(A, [1,M_ALL]);
	end
	W = zeros(N,M_ALL);
	SY  = SY0;
	% Save original input for pruning
	IX_act  = 1:M_ALL;
	IX_dim  = 1:Xdim;
else
	A = repmat(A0, [1,M_ALL]);
	W = zeros(N,M_ALL);
	SY  = SY0;
	% Save original input for pruning
	IX_act  = 1:M_ALL;
	IX_dim  = 1:Xdim;
end

ix_act  = 1:M;
ix_dim  = 1:Xdim;
M_all   = M;

% Delay embedding index for W
Wid = [(0:(D-1))'* Xdim + 1 , (1:D)'* Xdim];

% Working variable
G_A = zeros(N,M);       % N x M

if size(A,1)==Xdim && size(A,2)==1
	A   = repmat(A',[1 D]);
elseif size(A,1)==1 && size(A,2)==Xdim
	A   = repmat(A,[1 D]);
elseif size(A,1)==1 && size(A,2)==1
	A   = repmat(A,[1 M_ALL]);
end

SW = SY./( Tall*XX + SY./A );

fprintf('a_min = %g\n', a_min)
fprintf('SY0 = %g\n', SY0)
fprintf('SY  = %g\n', SY)

% Free energy histry
FE  = zeros(Ntrain,1);
LP  = zeros(Ntrain,1);
H   = zeros(Ntrain,1);
Err = zeros(Ntrain,1);
Mhist = zeros(Ntrain,1);

% ARD hyper param. history
if isfield(parm,'Debug') && ~isempty(parm.Debug) && parm.Debug > 0
	Debug = 1;
	A_tmp = zeros(N*M_ALL, ceil(Ntrain/Nskip));
else
	Debug = 0;
end

% recover all component
A_all  = zeros(1,M_all);       % N x M

[dY] = error_delay_time(X,Y,W,tau);

k_save  = 0;

%%%%%% Learning Loop %%%%%%
for k=1:Ntrain
	% Ainv = alpha , A = 1/alpha
	Ainv    = SY./A;

	% E = (Y-W*X)^2/SY +  W^2/A
	%   = ( (Y-W*X)^2 +  W^2 * (SY/A) )/SY
    [W]  = weight_update_embed(X, dY, W, Tall*XX, Ainv, tau);
	[dY] = error_delay_time(X, Y, W, tau);

    dYY = sum(dY(:).^2)/(N*Tall);
	WW  = sum(W.^2,1);

    % Noise variance update
    SY  = dYY + sum(WW .* Ainv)/(N*Tall);
    % Prevent zero variance
    SY  = max( SY, MINVAL);

    % Log variance
    SWA     = max( SW .* Ainv , MINVAL);
    log_sw  = N*(sum( log(SWA) - SWA + 1 ));
    log_sy  = N*( log(SY) );
    log_a   = Ta0*(sum(log(Ainv) - a0.*Ainv + 1));

    % Free energy
    LP(k)  = - (0.5*Tall) * log_sy ;
    H(k)   = 0.5*( log_sw + log_a - sum(WW .*Ainv) );
    FE(k)  = LP(k) + H(k);
    Err(k) = sum(dYY)/(N*SY0);

    % Weight variance

    SW = SY./( Tall*XX + Ainv );

    % Hyper parameter for weight variance (ARD)
	if mod(k,Nupdate) == 0,
	    % A = Alpha^-1 * SY
		% Ainv = SY./A;

		% N*A = N * (Alpha * SY) =  (W.^2) + N * SW  ;
		% VB update rule (Stable)
		A  = (WW + N*SW + 2*Ta0*a0)/( N + 2*Ta0 );

		% regularization parameter
		Am = mean(A);
		A  = repmat(Am, 1,M);
		A = max(A,MINVAL);

	end
	% END of if mod(k,Nupdate) == 0

	Mhist(k) = M;
	SW  = SW /SY;

    if mod(k, Nskip)==0
        % Save history
		if Debug == 1
        	k_save = k_save + 1;
        	A_tmp(:,k_save) = A_all(:);
		end

        fprintf('Iter = %4d, err = %g, F = %g, SY= %g, H = %g\n', ...
               k, Err(k), FE(k), SY, - 2*H(k)/log(T));
    end

    % Convergence check
	if k > Ncheck,
		Fdif = (FE(k) - FE(k-Fstep))/(abs(FE(k))+eps);
	else
		Fdif = Fdiff + 1;
	end

	if (Fdiff > abs(Fdif)),
		fprintf('Converged : Free energy change = %g\n',Fdif)
		break;
	end;
end

ix_act = IX_act(ix_act);

% Active index
Model.ix_act = ix_act;
Model.M_all  = M_ALL ;

% Save output variable
Model.A    = A ;
Model.W    = W ;
Model.SY   = SY;
Model.method = 'linear_sparse_stepwise';
Model.mode   = 'scalar';
Model.sparse = 'sparse';

% Save history
Info.FE  = FE(1:k);
Info.LP  = LP(1:k);
Info.H   = H(1:k) ;
Info.Err = Err(1:k);
Info.M   = Mhist(1:k);

if exist('A_tmp','var')
	Info.A   = A_tmp(:,1:k_save) ;
end

% recover all component
%W_all   = zeros(N,M_all);
%SW_all  = zeros(N,M_all);
%
%A_all(:,ix_act)  = A;
%W_all(:,ix_act)  = W ;  % N x M
%SW_all(:,ix_act) = SW;  % N x M

%%% ---- Index transformation from old active_index to current active_index
function	jx = trans_index(ix,ix_old,M)

N = length(ix_old);
Itrans = zeros(M,1);
Itrans(ix_old) = 1:N;

jx = Itrans(ix);

author: KamitaniLab

repo: GenericObjectDecoding


text:
function edgeFaces = trimeshEdgeFaces(faces, varargin)
%TRIMESHEDGEFACES Compute index of faces adjacent to each edge of a triangular mesh.
%
%   EF = trimeshEdgeFaces(FACES)
%   EF = trimeshEdgeFaces(VERTICES, FACES)
%   EF = trimeshEdgeFaces(VERTICES, EDGES, FACES)
%   Compute index array of faces adjacent to each edge of a mesh.
%   FACES is a NF-by-3 array containing vertex indices of each face. The
%   result EF is a NE-by-2 array containing the indices of the two faces
%   incident to each edge. If an edge belongs to only one face, the other
%   face index is ZERO.
%
%   The list of edges (as array of source and target vertex indices) can be
%   obtained from the function 'meshEdges'.
%
%   Note: faces are listed in increasing order for each edge, and no
%   information is kept about relative orientation of edge and face.
%
%   Example
%     % compute incidence list of each edge of an octahedron. For example,
%     % first edge is incident to faces 1 and 5. Second edge is incident to
%     % faces 4 and 8, and so on.
%     [v, f] = createOctahedron;
%     ef = trimeshEdgeFaces(v, f)
%     ef =
%          1     5
%          4     8
%          4     1
%          5     8
%          2     6
%          1     2
%          6     5
%          3     7
%          2     3
%          7     6
%          3     4
%          7     8
%
%   See also
%   meshes3d, meshEdgeFaces, trimeshMeanBreadth, meshEdges

% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2015-08-19,    using Matlab 8.5.0.197613 (R2015a)
% Copyright 2015 INRA - Cepia Software Platform.

if nargin == 2
    faces = varargin{1};
elseif nargin == 3
    faces = varargin{2};
end

% compute vertex indices of each edge (in increasing index order)
edges = sort([faces(:,[1 2]) ; faces(:,[2 3]) ; faces(:,[3 1])], 2);

% create an array to keep indices of faces "creating" each edge
nFaces = size(faces, 1);
edgeFaceInds = repmat( (1:nFaces)', 3, 1);

% sort edges, keeping indices
[edges, ia, ib] = unique(edges, 'rows'); %#ok<ASGLU>
nEdges = size(edges, 1);

% allocate memory for result
edgeFaces = zeros(nEdges, 2);

% iterate over edges, to identify incident faces
for iEdge = 1:nEdges
    inds = find(ib == iEdge);
    edgeFaces(iEdge, 1:length(inds)) = edgeFaceInds(inds);
end


text:
function res = ne(a,b)
%NE           Implements  a ~= b  for slopes, compares only a.r and b.r
%

% written  12/06/98     S.M. Rump
% modified 09/28/01     S.M. Rump  matrices and multi-dimensional arrays
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 05/21/09     S.M. Rump  typo
%

  if ~isa(a,'slope')
    asize = size(a);
    bsize = b.size;
    res = ( a(:)~=b.r(:,1) );
  elseif ~isa(b,'slope')
    asize = a.size;
    bsize = size(b);
    res = ( a.r(:,1)~=b(:) );
  else
    asize = a.size;
    bsize = b.size;
    res = ( a.r(:,1)~=b.r(:,1) );
  end

  if ~isequal(asize,bsize) & ( prod(asize)~=1 ) & ( prod(bsize)~=1 )
    error('incompatible size for slope ~= ')
  end

  if prod(asize)==1
    res = reshape(res,bsize);
  else
    res = reshape(res,asize);
  end


text:
function varargout = imdown(varargin)
% VL_IMDOWN  Downsample an image by two
%   J = VL_IMDOWN(I) downsamples the image I by half by discarding
%   each other pixel.
%
%   VL_IMDOWN() accepts the following options:
%
%   Method:: Sample
%     'Sample' downsamples the image by discarding
%     pixels. 'Average' instead averages groups of 2x2 pixels.
%
%   See also: VL_IMUP(), VL_HELP().
[varargout{1:nargout}] = vl_imdown(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
%kill_spaces Remove all spaces from given string.

function r = kill_spaces(s)

r = s(setdiff(1:length(s),strfind(s,' ')));

author: strawlab

repo: MultiCamSelfCal


text:
%% FUNCTION Least_msmtfl_capL1
%   Multi-Stage Multi-Task Feature Learning
%
%% OBJECTIVE
%   min_W ||XW - Y||_F^2 + lambda*\sum_i \min{||W^i||_1,theta}
%
%   It's a nonconvex optimization problem, which can be relaxed into a Multi-Stage Convex optimization:
%   min_W ||XW - Y||_F^2 + \sum_i{gamma_i*||W^i||_1}
%
%% INPUT
%   X: {n * d} * t - input matrix
%   Y: {n * 1} * t - output matrix
%   lambda: regularized parameter (vector)
%   theta: theresholding paramter
%
%   (Optional)
%   opts.lFlag: estimate the upper bound of Lipschitz constant if nonzero, zero otherwise
%
%% OUTPUT
%   W: output weight
%   fun: function values
%
%% LICENSE
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as pubtolwlished by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou, Pinghua Gong and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on Dec 18, 2012.
%
%
%% RELATED PAPERS
%
% [1] Pinghua Gong, Jieping Ye, Changshui Zhang. Multi-Stage Multi-Task
%     Feature Learning. The 26th Annual Conference on Neural Information
%     Processing Systems (NIPS 2012), Lake Tahoe, Nevada, USA,
%     December 3-6, 2012.
%
%% RELATED FUNCTIONS
%  init_opts, combine_input (utils)

function [W,funcVal] = Least_msmtfl_capL1(X, Y, lambda, theta, opts)

if nargin <4
    error('\n Inputs: X, Y, and lambda1, and lambda2 should be specified!\n');
end
if nargin <5
    opts = [];
end

% initialize options.
opts=init_opts(opts);

% initial Lipschiz constant.
if isfield(opts, 'lFlag')
    lFlag = opts.lFlag;
else
    lFlag = false;
end

task_num = length(X);
[X, y, ~, samplesize] = combine_input(X, Y);
dimension = size(X, 2);

% initialize a starting point
if opts.init==2
    W0 = zeros(dimension, task_num);
elseif opts.init == 0
    W0 = randn(dimension, task_num);
else
    if isfield(opts,'W0')
        W0=opts.W0;
        if (nnz(size(W0)-[dimension, task_num]))
            error('\n Check the input .P0');
        end
    else
        W0=randn(dimension, task_num);
    end
end

% Set an array to save the objective value
funcVal = [];

W = W0;

[d,m] = size(W); % d: dimension, m: the number of tasks
X = diagonalize(X,samplesize);
XtX = X'*X; Xty = X'*y;

L1norm = max(sum(abs(X),1)); Linfnorm = max(sum(abs(X),2));

if lFlag
    % Upper bound for largest eigenvalue of Hessian matrix
    L = 2*min([L1norm*Linfnorm; size(X,1)*Linfnorm*Linfnorm; size(X,2)*L1norm*L1norm; size(X,1)*size(X,2)*max(abs(X(:)))]);
else
    % Lower bound for largest eigenvalue of Hessian matrix
    L = 2*max(L1norm*L1norm/size(X,1),Linfnorm*Linfnorm/size(X,2));
end
% Initial function value
funcVal = cat(1, funcVal, norm(X*W(:) - y)^2 + lambda*(sum(min(sum(abs(W),2),theta))));

tolw = 1e-5; % precision for inner iterations.

for iter = 1:opts.maxIter
    if iter == 1
        weight = lambda*ones(d,m);
    else
        weight = lambda*(repmat(sum(abs(W),2) < theta,1,m));
    end
    [W,~,iterw] = wLassomtl(X,y,XtX,Xty,weight,W,tolw,100,L,lFlag);

    if iterw == 1
        tolw = tolw/4;
    end

    funcVal = cat(1, funcVal, norm(X*W(:) - y)^2 + lambda*(sum(min(sum(abs(W),2),theta))));

    % stopping condition
    switch(opts.tFlag)
        case 0
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <= opts.tol)
                    break;
                end
            end
        case 1
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <=...
                        opts.tol* funcVal(end-1))
                    break;
                end
            end
        case 2
            if ( funcVal(end)<= opts.tol)
                break;
            end
        case 3
            if iter>=opts.maxIter
                break;
            end
    end

end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function p06_story ( )

%*****************************************************************************80
%
%% P06_STORY prints the "story" for problem p06.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carl DeBoor, John Rice,
%    Least-squares cubic spline approximation II - variable knots.
%    Technical Report CSD TR 21,
%    Purdue University, Lafayette, Indiana, 1968.
%
%    Carl DeBoor,
%    A Practical Guide to Splines,
%    Springer, 2001,
%    ISBN: 0387953663,
%    LC: QA1.A647.v27.
%
%  Parameters:
%
%    None
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data is due to deBoor and Rice.\n' );
  fprintf ( 1, '  The data represents a temperature dependent property of titanium.\n' );
  fprintf ( 1, '  The data has been used extensively as an example in spline\n' );
  fprintf ( 1, '  approximation with variably-spaced knots.\n' );
  fprintf ( 1, '  DeBoor considers two sets of knots:\n' );
  fprintf ( 1, '  (595,675,755,835,915,995,1075)\n' );
  fprintf ( 1, '  and\n' );
  fprintf ( 1, '  (595,725,850,910,975,1040,1075).\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [p] = gui_erp_open(p,command,parent)

% gui_erp_open - Load averaged ERP data into matlab workspace
%
% Usage: [p] = gui_erp_open(p,[command],[parent])
%
% p is a structure, generated by 'eeg_toolbox_defaults'
% command is either 'init' or 'load'
% parent is a handle to the gui that calls this gui, useful
% for updating the UserData field of the parent from this gui.
% The p structure is returned to the parent when the parent
% handle is given.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:56 $

% Licence:  GNU GPL, no express or implied warranties
% History:  01/2002, Darren.Weber_at_radiology.ucsf.edu
%           08/2002, Darren.Weber_at_radiology.ucsf.edu
%                    added EMSE erp handling
%                    added interpolation of zero point option
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('p','var'),
 [p] = eeg_toolbox_defaults;
elseif isempty(p),
 [p] = eeg_toolbox_defaults;
end

if ~exist('command','var'),
  command = 'init';
elseif isempty(command),
  command = 'init';
end

command = lower(command);

switch command,
  case 'init',
    if exist('parent','var'),
      ERPOpen = INIT(p,parent);
    else
      ERPOpen = INIT(p,'');
    end
  otherwise,
    ERPOpen = get(gcbf,'Userdata');
    set(ERPOpen.gui,'Pointer','watch');
end

switch command,

  case 'plot',

    ERPOpen.p = erp_open(ERPOpen.p,ERPOpen.gui);

    plotfig = figure('Name',ERPOpen.p.volt.file,...
      'NumberTitle','off',...
      'UserData',ERPOpen.p);
    movegui(plotfig,'center');

    plot(ERPOpen.p.volt.timeArray,ERPOpen.p.volt.data); axis tight;
    eeg_plot_metric;

    if isequal(get(ERPOpen.handles.Bhold,'Value'),1),
      parent = ERPOpen.gui;
    else
      if isfield(ERPOpen,'parent'),
        parent = ERPOpen.parent.gui;
      else
        parent = [];
      end
    end

    if isempty(parent), [Xpoint,Ypoint] = eeg_crosshair('init',ERPOpen.p);
    else                [Xpoint,Ypoint] = eeg_crosshair('init',ERPOpen.p,parent);
    end

  case 'return',

    ERPOpen.p = erp_open(ERPOpen.p);

  case 'save',

    fprintf('\nGUI_ERP_OPEN: Save As not implemented yet.\n');

  case 'cancel',

    GUI.parent = ERPOpen.parent;
    gui_updateparent(GUI);
    close gcbf;

  otherwise,

end

switch command,
  case {'init','cancel'},
  otherwise,
    set(ERPOpen.gui,'Pointer','arrow');
    set(ERPOpen.gui,'Userdata',ERPOpen);
    if isequal(get(ERPOpen.handles.Bhold,'Value'),1),
     [p] = gui_updateparent(ERPOpen,0);
    else
     [p] = gui_updateparent(ERPOpen);
      close gcbf;
    end
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ERPOpen] = INIT(p,parent)
% GUI General Parameters

GUIwidth  = 500;
GUIheight = 120;

version = '$Revision: 1.1 $';
name = sprintf('ERP File Open [v %s]\n',version(11:15));

GUI = figure('Name',name,'Tag','ERP_OPEN',...
  'NumberTitle','off',...
  'MenuBar','none','Position',[1 1 GUIwidth GUIheight]);
movegui(GUI,'center');

Font.FontName   = 'Helvetica';
Font.FontUnits  = 'Pixels';
Font.FontSize   = 12;
Font.FontWeight = 'normal';
Font.FontAngle  = 'normal';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Voltage Data Selection and Parameters

G.Title_data = uicontrol('Parent',GUI,'Style','text','Units','Normalized',Font, ...
  'Position',[.01 .75 .17 .2],...
  'String','Data Type:','HorizontalAlignment','left');

switch lower(p.volt.type),
  case 'ascii',   voltType = 1;
  case 'emse',    voltType = 2;
  case 'scan4x',  voltType = 3;
  case 'scan3x',  voltType = 4;
  case 'matlab',  voltType = 5;
  otherwise,      voltType = 1;
end
G.PvoltType = uicontrol('Tag','PvoltType','Parent',GUI,'Style','popupmenu',...
  'Units','Normalized',Font,  ...
  'Position',[.20 .75 .2 .2],...
  'String',{'ASCII' 'EMSE' 'Scan4x' 'Scan3x' 'Matlab'},'Value',voltType,...
  'Callback',strcat('ERPOpen = get(gcbf,''Userdata'');',...
  'ERPOpen.p.volt.type = popupstr(ERPOpen.handles.PvoltType);',...
  'set(gcbf,''Userdata'',ERPOpen); clear ERPOpen;'));

% ERP Parameters
G.Bascii = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
  'Position',[.45 .75 .2 .2],...
  'String','ERP Parameters','BusyAction','queue',...
  'BackgroundColor',[0.0 0.0 0.75],...
  'ForegroundColor',[1 1 1], 'HorizontalAlignment', 'center',...
  'Callback','tempgui = gui_eeg_ascii_parameters(gcbf); clear tempgui;');

% Interpolate Zero
G.BinterpZero = uicontrol('Parent',GUI,'Style','checkbox','Units','Normalized', Font, ...
  'Position',[.80 .75 .18 .2],'String','Interp 0','BusyAction','queue',...
  'TooltipString','Interpolate zero, mainly for ascii, scan3x, scan4x files.',...
  'Value',p.volt.interpZero,'HorizontalAlignment', 'center',...
  'Callback',strcat('ERPOpen = get(gcbf,''Userdata'');',...
  'ERPOpen.p.volt.interpZero = get(ERPOpen.handles.BinterpZero,''Value''); ',...
  'set(gcbf,''Userdata'',ERPOpen); clear ERPOpen;'));

G.Title_path = uicontrol('Parent',GUI,'Style','text','Units','Normalized',Font, ...
  'Position',[.01 .50 .17 .2],...
  'String','Path','HorizontalAlignment','left');
G.EvoltPath = uicontrol('Parent',GUI,'Style','edit','Units','Normalized',Font,  ...
  'Position',[.20 .50 .58 .2], 'String',p.volt.path,...
  'Callback',strcat('ERPOpen = get(gcbf,''Userdata'');',...
  'ERPOpen.p.volt.path = get(ERPOpen.handles.EvoltPath,''String'');',...
  'set(gcbf,''Userdata'',ERPOpen); clear ERPOpen;'));

G.Title_file = uicontrol('Parent',GUI,'Style','text','Units','Normalized',Font, ...
  'Position',[.01 .25 .17 .2],...
  'String','File','HorizontalAlignment','left');
G.EvoltFile = uicontrol('Parent',GUI,'Style','edit','Units','Normalized',Font,  ...
  'Position',[.20 .25 .58 .2], 'String',p.volt.file,...
  'Callback',strcat('ERPOpen = get(gcbf,''Userdata'');',...
  'ERPOpen.p.volt.file = get(ERPOpen.handles.EvoltFile,''String'');',...
  'set(gcbf,''Userdata'',ERPOpen); clear ERPOpen;'));

Font.FontWeight = 'bold';

% BROWSE: Look for the data
browsecommand = strcat('ERPOpen = get(gcbf,''Userdata'');',...
  'cd(ERPOpen.p.volt.path);',...
  '[file, path] = uigetfile(',...
  '{''*.dat;*.asc;*.txt'', ''ASCII Data Files (*.dat,*.asc,*.txt)'';', ...
  ' ''*.erp'', ''NeuroScan Average (*.erp)'';', ...
  ' ''*.erp'', ''EMSE Average (*.erp)'';', ...
  ' ''*.mat'', ''MATLAB Data Files (*.mat)'';', ...
  ' ''*.*'',   ''All Files (*.*)''},', ...
  '''Select ERP File'');',...
  'if ~isequal(path,0), ERPOpen.p.volt.path = path; end;',...
  'if ~isequal(file,0), ERPOpen.p.volt.file = file; end;',...
  'set(ERPOpen.handles.EvoltPath,''String'',ERPOpen.p.volt.path);',...
  'set(ERPOpen.handles.EvoltFile,''String'',ERPOpen.p.volt.file);',...
  'if ~isempty(findstr(file,''.dat'')) | ~isempty(findstr(file,''.txt'')) | ~isempty(findstr(file,''.asc'')), ',...
  '  set(ERPOpen.handles.PvoltType,''Value'',1); ',...
  '  ERPOpen.p.volt.type = popupstr(ERPOpen.handles.PvoltType);',...
  'elseif ~isempty(findstr(file,''.erp'')), ',...
  '  set(ERPOpen.handles.PvoltType,''Value'',3); ',...
  '  ERPOpen.p.volt.type = popupstr(ERPOpen.handles.PvoltType);',...
  'else, ',...
  '  set(ERPOpen.handles.PvoltType,''Value'',5); ',...
  '  ERPOpen.p.volt.type = popupstr(ERPOpen.handles.PvoltType);',...
  'end;',...
  'set(gcbf,''Userdata'',ERPOpen); clear ERPOpen file path;');
G.BvoltFile = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized',Font, ...
  'Position',[.01 .01 .17 .2], 'String','BROWSE',...
  'BackgroundColor',[0.8 0.8 0.0],...
  'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
  'Callback', browsecommand );

% PLOT: Load & plot the data!
G.Bplot = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
  'Position',[.20 .01 .18 .2],...
  'String','PLOT','BusyAction','queue',...
  'TooltipString','Plot the ERP data and return p struct.',...
  'BackgroundColor',[0.0 0.5 0.0],...
  'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
  'Callback',strcat('ERPOpen = get(gcbf,''Userdata'');',...
  'p = gui_erp_open(ERPOpen.p,''plot'');',...
  'clear ERPOpen;'));

% Save As
G.Bsave = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
  'Position',[.40 .01 .18 .2],'HorizontalAlignment', 'center',...
  'String','SAVE AS','TooltipString','ERP File Conversion Tool (not implemented yet)',...
  'BusyAction','queue',...
  'Visible','on',...
  'BackgroundColor',[0.0 0.0 0.75],...
  'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
  'Callback',strcat('ERPOpen = get(gcbf,''Userdata'');',...
  'p = gui_erp_open(ERPOpen.p,''save'');',...
  'clear ERPOpen;'));

% Quit, return file parameters
G.Breturn = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
  'Position',[.60 .01 .18 .2],...
  'String','RETURN','BusyAction','queue',...
  'TooltipString','Return p struct to workspace and parent GUI.',...
  'BackgroundColor',[0.75 0.0 0.0],...
  'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
  'Callback',strcat('ERPOpen = get(gcbf,''Userdata'');',...
  'p = gui_erp_open(ERPOpen.p,''return'');',...
  'clear ERPOpen;'));

% Cancel
G.Bcancel = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
  'Position',[.80 .01 .18 .2],...
  'String','CANCEL','BusyAction','queue',...
  'TooltipString','Close, do not return parameters.',...
  'BackgroundColor',[0.75 0.0 0.0],...
  'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
  'Callback',strcat('ERPOpen = get(gcbf,''Userdata'');',...
  'p = gui_erp_open(ERPOpen.p,''cancel'');',...
  'clear ERPOpen;'));

% Help
G.Bhelp = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
  'Position',[.80 .25 .18 .2],'String','Help','BusyAction','queue',...
  'BackgroundColor',[1 1 0],...
  'ForegroundColor',[0 0 0], 'HorizontalAlignment', 'center',...
  'Callback','doc erp_open;');

G.Bhold = uicontrol('Parent',GUI,'Style','checkbox','Units','Normalized', Font, ...
  'Position',[.80 .50 .18 .2],'String','Hold GUI','BusyAction','queue',...
  'TooltipString','ERP File Load GUI remains open after ''Plot'' or ''Return'' commands.',...
  'Value',p.hold,'HorizontalAlignment', 'center');

% Store userdata
if exist('parent','var'), ERPOpen.parent.gui = parent; end
ERPOpen.gui = GUI;
ERPOpen.handles = G;
ERPOpen.p = p;
set(GUI,'Userdata',ERPOpen);
set(GUI,'HandleVisibility','callback');

return

author: PatternRecognition

repo: OpenBMI


text:
function value = torus_area_3d ( r1, r2 )

%*****************************************************************************80
%
%% TORUS_AREA_3D returns the area of a torus in 3D.
%
%  Integration region:
%
%    Points (X,Y,Z) such that:
%
%      ( SQRT ( X*X + Y*Y ) - R1 )**2 + Z*Z = R2*R2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R1, R2, the two radii that define the torus.
%
%    Output, real TORUS_AREA_3D, the area of the torus.
%
  value = 4.0E+00 * pi * pi * r1 * r2;

  return
end


text:
function [OptExchRxns] = findOptExchRxns(model,Ex_Rxns,parFVA)
% USAGE:
%
%    [OptExchRxns] = findOptExchRxns(model, Ex_Rxns, parFVA)
%
% INPUTS:
%    model:         Metabolic model
%    Ex_Rxns:       Vector of exchange reactions for FVA
%    fastFVA:       use FastFVA (default=0)
%
% OUTPUTS:
%    OptExchRxns:
%
% .. Authors:
%       - Ines Thiele 2014
%       - Maike K. Aurich 27/05/15, remove hidden values

tol = -1e-6;  % Default tolerance

if ~exist('parFVA','var') || parFVA == 0

    matlabpool = 0;
    parFVA = 0; % default no parallelization of FVA
end

[a1(:,1),a1(:,2)]=fluxVariability(model,1,[],Ex_Rxns);
% find all uptakes that have ub=0;
if size(a1,1)==length(Ex_Rxns)
    OptSecr = Ex_Rxns;
    R = OptSecr(find(a1(:,2)<abs(tol)));
    R = [R; OptSecr(find(a1(:,1)>abs(tol)))];
    OptSecr(ismember(OptSecr, R))=[];% remove obligate uptake
    clear R
    OptUptake = Ex_Rxns;
    R = OptUptake(find(a1(:,2)<(tol)));
    R = [R; OptUptake(find(a1(:,1)>(tol)))];
    OptUptake(ismember(OptUptake, R))=[];% remove obligate uptake
    OptExchRxns = [OptUptake;OptSecr];
else
    OptExchRxns=[];
end

author: opencobra

repo: cobratoolbox