text:
function varargout = process_pac_analysis( varargin )
% PROCESS_PAC_ANALYSIS: Further analysis of tpac maps.
%
% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Soheila Samiee, 2014-2017
%   - 2.0: SS. Aug. 2017
%                - Imported in public brainstorm rep
%   - 2.1: SS. Jul. 2018
%                - Bug fix in average over sources (mean)
%
eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Basic Analysis of tPAC maps';
    sProcess.FileTag     = '';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = {'Frequency','Time-resolved Phase-Amplitude Coupling'};
    sProcess.Index       = 1019;
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'timefreq'};
    sProcess.OutputTypes = {'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    sProcess.isSeparator = 0;

    % === ANALYSIS TO BE DONE
    sProcess.options.label.Comment = 'Analysis:';
    sProcess.options.label.Type    = 'label';
    sProcess.options.analyze_type.Comment = {'Mean (Over sources)', ...
        'Median (Over sources)','Z-score on time (If no negative time, on total recording)', ...
        'Mean (Over time)'};
    sProcess.options.analyze_type.Type    = 'radio';
    sProcess.options.analyze_type.Value   = 1;

        % === Using phase
    sProcess.options.usePhase.Comment = 'Use phase in averaging (mean)';
    sProcess.options.usePhase.Type    = 'checkbox';
    sProcess.options.usePhase.Value   = 0;

end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInput) %#ok<DEFNU>
    OutputFiles = {};
    isMean   = 0;
    isMedian = 0;
    isZscore = 0;
    isTimeMean = 0;

    usePhase = sProcess.options.usePhase.Value;

    % Get options
    if sProcess.options.analyze_type.Value ==1
        isMean  = 1;
        tag = '| mean';
    elseif sProcess.options.analyze_type.Value ==2
        isMedian  = 1;
        tag = '| median';
    elseif sProcess.options.analyze_type.Value ==3
        isZscore = 1;
        tag = '| zscore';
    elseif sProcess.options.analyze_type.Value ==4
        isTimeMean = 1;
        tag = '| TimeMean';
    end

    % Load TF file
    tpacMat = in_bst_timefreq(sInput(1).FileName, 0);
    % Error
    if isempty(tpacMat)
        bst_report('Error', sProcess, sInput, Messages);
        return;
    end

    % Apply the appropriate function
    tpac_avg = tpacMat.sPAC.DynamicPAC;
    if usePhase
        tpac_avg_phase = tpacMat.sPAC.DynamicPhase;
    end

    if isZscore
        iBaseline = find(tpacMat.Time<0);
        if isempty(iBaseline)
            iBaseline = 1:length(tpacMat.Time);
        end
        tpac_avg = process_zscore('Compute', tpac_avg, iBaseline);

    elseif isMean
        if usePhase
            tmp = mean(tpac_avg.*exp(1i*tpac_avg_phase),1);
            tpac_avg = abs(tmp);
            tpac_avg_phase = angle(tmp);
        else
            tpac_avg = mean(tpac_avg,1);
        end

        if length(sInput)>1
            N = length(sInput);
            for iFile = 2:N
                TimefreqMat2 = in_bst_timefreq(sInput(iFile).FileName, 0);
                if usePhase
                    tmp = mean(TimefreqMat2.sPAC.DynamicPAC.*exp(1i*TimefreqMat2.sPAC.DynamicPhase),1);
                    TimefreqMat2.sPAC.DynamicPAC  = abs(tmp);
                    TimefreqMat2.sPAC.DynamicPhase  = angle(tmp);
                else
                    TimefreqMat2.sPAC.DynamicPAC = mean(TimefreqMat2.sPAC.DynamicPAC,1);
                end
                TimefreqMat2.TF = TimefreqMat2.sPAC.DynamicPAC;
                %Saving the files
                TimefreqMat2.Comment = [TimefreqMat2.Comment, ' ', tag];
                % Output filename: add file tag
                FileTag = strtrim(strrep(tag, '|', ''));
                pathName = file_fullpath(sInput(iFile).FileName);
                OutputFile = strrep(pathName, '.mat', ['_' FileTag '.mat']);
                OutputFile = file_unique(OutputFile);
                % Save file
                bst_save(OutputFile, TimefreqMat2, 'v6');
                % Add file to database structure
                db_add_data(sInput(iFile).iStudy, OutputFile, TimefreqMat2);
                OutputFiles{end + 1} = OutputFile;
            end
        end

    elseif isMedian
        tpac_avg = median(tpac_avg,1);
        if length(sInput)>1
            N = length(sInput);
            for iFile = 2:N
                TimefreqMat2 = in_bst_timefreq(sInput(iFile).FileName, 0);
                TimefreqMat2.sPAC.DynamicPAC = median(TimefreqMat2.sPAC.DynamicPAC,1);
                TimefreqMat2.TF = median(TimefreqMat2.sPAC.DynamicPAC,1);
                %Saving the files
                TimefreqMat2.Comment = [TimefreqMat2.Comment, ' ', tag];
                % Output filename: add file tag
                FileTag = strtrim(strrep(tag, '|', ''));
                pathName = file_fullpath(sInput(iFile).FileName);
                OutputFile = strrep(pathName, '.mat', ['_' FileTag '.mat']);
                OutputFile = file_unique(OutputFile);
                % Save file
                bst_save(OutputFile, TimefreqMat2, 'v6');
                % Add file to database structure
                db_add_data(sInput(iFile).iStudy, OutputFile, TimefreqMat2);
                OutputFiles{end + 1} = OutputFile;
            end
        end

    elseif isTimeMean
        if length(sInput)==1
            tpac_avg = repmat(mean(tpac_avg,2),[1,size(tpac_avg,2),1]);
            [PACmax,tmp] = max(abs(tpac_avg),[],1);
            tpacMat.TF = squeeze(PACmax)';
        elseif length(sInput)>1
            N = length(sInput);
            for iFile = 2:N
                tPACMat2 = in_bst_timefreq(sInput(iFile).FileName, 0);
%                 TimefreqMat2.sPAC.DynamicPAC = median(TimefreqMat2.sPAC.DynamicPAC,1);
%                 TimefreqMat2.TF = median(TimefreqMat2.sPAC.DynamicPAC,1);
                tpac_avg = repmat(mean(tPACMat2.sPAC.DynamicPAC,2),[1,size(tPACMat2.sPAC.DynamicPAC,2),1]);
                [PACmax,tmp] = max(abs(tpac_avg),[],1);
                tPACMat2.TF = tpac_avg;%squeeze(PACmax)';
                tPACMat2.sPAC.DynamicPAC = tpac_avg;

                %Saving the files
                tPACMat2.Comment = [tPACMat2.Comment, ' ', tag];
                % Output filename: add file tag
                FileTag = strtrim(strrep(tag, '|', ''));
                pathName = file_fullpath(sInput(iFile).FileName);
                OutputFile = strrep(pathName, '.mat', ['_' FileTag '.mat']);
                OutputFile = file_unique(OutputFile);
                % Save file
                bst_save(OutputFile, tPACMat2, 'v6');
                % Add file to database structure
                db_add_data(sInput(iFile).iStudy, OutputFile, tPACMat2);
                OutputFiles{end + 1} = OutputFile;
            end
        end
        tpacMat.sPAC.DynamicPAC = tpac_avg;
        tpacMat.TF = tpac_avg;
        if usePhase
            tpacMat.sPAC.DynamicPhase = tpac_avg_phase;
        end
    end

    tpacMat.TF = tpac_avg;
    tpacMat.sPAC.DynamicPAC = tpac_avg;
    if usePhase
        tpacMat.sPAC.DynamicPhase = tpac_avg_phase;
    end

    % === SAVING THE DATA IN BRAINSTORM ===
    % Getting the study
    [sOutputStudy, iOutputStudy] = bst_process('GetOutputStudy', sProcess, sInput(1));
    % Comment
    tpacMat.Comment = [tpacMat.Comment, ' ', tag];
    % Output filename: add file tag
    FileTag = strtrim(strrep(tag, '|', ''));
    pathName = file_fullpath(sInput(1).FileName);
    % Preparing the output file
    OutputFile = strrep(pathName, '.mat', ['_' FileTag '.mat']);
    OutputFile = file_unique(OutputFile);
    % Save on disk
    bst_save(OutputFile, tpacMat, 'v6');
    % Register in database
    db_add_data(iOutputStudy, OutputFile, tpacMat);
    OutputFiles{end + 1} = OutputFile;
end


text:
%
% Visualize a MotorProgram
%
% additional arguments 'motor' plots the motor trace,
% and the arguments 'image' and 'pmap' can select
% betwen the image and the probability map
%
%  Other parameters
%  'size_start' : where 0 if we don't want to plot the start points
%  'lw' : line width for strokes
function vizMP(M,varargin)

    assert(isa(M,'MotorProgram'));

    n = length(varargin);
    use_pimg = any(strcmp('pmap',varargin)) || any(strcmp('pimg',varargin));
    use_img = any(strcmp('image',varargin));
    assert(~(use_pimg && use_img));
    use_motor = any(strcmp('motor',varargin));

    % what is the size of the start marker?
    fstart = find(strcmp('size_start',varargin));
    size_start = [];
    if ~isempty(fstart)
        size_start = varargin{fstart+1};
    end
    assert(isnumeric(size_start));

    % what is the line width?
    flw = find(strcmp('lw',varargin));
    lw = [];
    if ~isempty(flw)
       lw = varargin{flw+1};
    end

    if use_motor
        if use_pimg
            plot_motor_to_image(M.pimg,M.motor_warped,size_start,lw);
        else
            plot_motor_to_image(M.I,M.motor_warped,size_start,lw);
        end
    else
        if use_pimg
            plot_image_only(M.pimg);
        else
            plot_image_only(M.I);
        end
    end

end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
%
% Type-level resampling using MCMC
%
% Sample a chain using ps.mcmc.nsamp_type_chain iterations
%
% Store a subset (ps.mcmc.nsamp_type_store) of these samples linearly
% spaced along the chain
%
% Return
%   samples_store: [nsamp_store x 1 cell]
%
function samples_store = RunMCMCType(M,lib)

    ps = defaultps;

    % fill-in type-level shape variables
    Q = M.copy();
    if isempty(Q.S{1}.shapes_type)
       for sid=1:M.ns
          Q.S{sid}.shapes_type = Q.S{sid}.shapes_token;
       end
    end

    % fill-in type-level eval variable
    ns = Q.ns;
    for sid=1:ns
       R = Q.S{sid}.R;
       if strcmp(R.type,'mid') && isempty(Q.S{sid}.R.eval_spot_type)
           Q.S{sid}.R.eval_spot_type = Q.S{sid}.R.eval_spot_token;
       end
    end

    % run mcmc
    samples_chain = mcmc_all(Q,lib,ps.mcmc.nsamp_type_chain,'type');

    % pick well-spaced set
    int = ps.mcmc.nsamp_type_chain ./ ps.mcmc.nsamp_type_store;
    indx = round(linspace(int,ps.mcmc.nsamp_type_chain,ps.mcmc.nsamp_type_store));
    samples_store = samples_chain(indx);

end


text:
function [ cvx_ver2, cvx_bld ] = cvx_version
cvx_ver = 1.21;
cvx_bld = '795';
cvx_bdate = '2010-05-31 08:26:19';
cvx_dbld = '795';
cvx_ddate = '2010-05-31 08:26:19';
if nargout == 0,
   fprintf( '\n' );
   fprintf( 'CVX version %g\n', cvx_ver );
   fprintf( '    Code: build %s, %s\n', cvx_bld, cvx_bdate );
   fprintf( '    Documentation: build %s, %s\n', cvx_dbld, cvx_ddate );
   if exist( 'OCTAVE_VERSION', 'var' ),
       fprintf( 'GNU Octave %s on %s\n', version, computer );
       fprintf( 'NOTE: Sorry, Octave support is not yet functional.\n' );
   else
       verd = ver('MATLAB');
       fprintf( 'MATLAB version %s %s on %s\n', verd.Version, verd.Release, computer );
   end
   fprintf( '\n' );
else
    cvx_ver2 = cvx_ver;
end

repo: REST


text:
% STD_FILECHECK - Check if ERSP or SPEC file contain specific parameters.
%                  This file must contain a Matlab structure with a field
%                  named 'parameter'. The content of this field will be
%                  compared to the 'params' input. If they are identical
%                  the output flag will indicate that recomputing this
%                  file is not necessary. If they are different, the user
%                  is queried ('guion' option) to see if he wishes to use
%                  the new parameters and recompute the file (not done in
%                  this function) or if he wishes to use the parameters
%                  of the file on disk.
%
% Usage:
% >> [ recompflag params ] = std_filecheck(filename, params, mode, ignorefields);
%
% Inputs:
%   filename     - [string] file containing a given measure (ERSP data for
%                  instance).
%   params       - [cell array or structure] cell array of parameters or
%                  structure. This is compared to the 'parameters' field
%                  in the data file.
%   mode         - ['guion'|'usedisk'|'recompute'] 'guion' query the user
%                  if the disk and input parameters are different. The
%                  outcome may be either 'usedisk' or 'recompute'. See
%                  recompflag output for more information.
%   ignorefields - [cell array] list fields to ignore
%
% Outputs:
%   recompflag - ['same'|'different'|'recompute'|'usedisk'|'cancel'] 'same'
%                (resp. 'different') indicates that the parameters in the
%                data file are identical (resp. different). 'recompute'
%                indicate that the measure should be recomputed and the
%                file has been erased. 'usedisk' indicate that the user
%                wishes (from the GUI) to use the version on disk.
%   params     - [structure] final parameter. This is usually identical
%                to the 'params' input unless the user choose to use
%                parameters from the file on disk. These are then
%                copied to this output structure.
%
% Authors: Arnaud Delorme, 2006-

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

function [ res, params2 ] = std_filecheck(filename, params2, guiflag, ignorefields);

    if nargin < 2
        help std_filecheck;
        return;
    end
    if nargin < 3
        guiflag = 'guion';
    end
    if nargin < 4
        ignorefields = {};
    end

    if ~exist( filename ), res = guiflag; return; end

    params1 = load('-mat', filename, 'parameters');
    params1 = finputcheck( params1.parameters, { 'tmp' 'real' [] NaN}, '', 'ignore'); % allow to tackle duplicate fields
    params1 = rmfield(params1, 'tmp');
    if iscell(params2), params2 = struct(params2{:}); end

    % test if the fields are different
    % --------------------------------
    params1 = orderfields(params1);
    params2 = orderfields(params2);
    fields1 = fieldnames( params1 ); fields1 = setdiff_bc( fields1, ignorefields);
    fields2 = fieldnames( params2 ); fields2 = setdiff_bc( fields2, ignorefields);
    allfields = union_bc(fields1, fields2);

    % make fields the same
    % --------------------
    res = 'same';
    if ~isequal( fields1, fields2 ),
        for ind = 1:length(allfields)
            if strcmpi(allfields{ind}, 'plotitc'), adsfads; end
            if ~isfield( params1, allfields{ind})
                if ischar(getfield(params2, allfields{ind}))
                     params1 = setfield(params1, allfields{ind}, '');
                else params1 = setfield(params1, allfields{ind}, []);
                end
                res = 'different';
            end
            if ~isfield( params2, allfields{ind})
                if ischar(getfield(params1, allfields{ind}))
                     params2 = setfield(params2, allfields{ind}, '');
                else params2 = setfield(params2, allfields{ind}, []);
                end
                res = 'different';
            end
        end
    end

    % data type
    % ---------
    if ~isempty(strmatch('cycles', allfields)), strcom = 'ERSP';
    else                                        strcom = 'SPECTRAL';
    end

    % compare fields
    % --------------
    txt = {};
    for ind = 1:length(allfields)
        if ~strcmp(allfields{ind},'baseline')
            val1 = getfield(params1, allfields{ind});
            val2 = getfield(params2, allfields{ind});
            val1str = fastif(isempty(val1), 'not set', vararg2str({val1(1:min(3, length(val1)))}));
            val2str = fastif(isempty(val2), 'not set', vararg2str({val2(1:min(3, length(val2)))}));

            tmptxt = sprintf('        ''%s'' is %s in the file (vs. %s)', allfields{ind}, val1str, val2str);
            if length(val1) ~= length(val2)
                res        = 'different';
                txt{end+1} = tmptxt;
            elseif ~isequal(val1, val2)
                if ~isnan(val1) && ~isnan(val2)
                    res        = 'different';
                    txt{end+1} = tmptxt;
                end
            end
        end
    end

    % build gui or return
    % -------------------
    if strcmpi(guiflag, 'usedisk'),
        if isempty(txt)
            params2 = params1;
            res     = 'usedisk';
            disp(['Using file on disk: ' filename ]);
            return;
        else
            strvcat(txt{:})
            error([ 'Two ' strcom ' files had different parameters and the ' strcom ' function' 10 ...
                          'cannot handle that. We suggest that you delete all files. See command line details' ]);
        end
    elseif strcmpi(guiflag, 'recompute'),
        res     = 'recompute';
        disp(['Deleting and recomputing file: ' filename ]);
        return;
    elseif strcmpi(res, 'same') && ( strcmpi(guiflag, 'guion') || strcmpi(guiflag, 'same') )
        disp(['Using file on disk: ' filename ]);
        return;
    end

    set_yes = [ 'set(findobj(''parent'', gcbf, ''tag'', ''ersp_no''), ''value'', 0);'];
    set_no  = [ 'set(findobj(''parent'', gcbf, ''tag'', ''ersp_yes''), ''value'', 0);' ];
    textgui1 = strvcat( [  upper(strcom)  ' info exists in file: ' filename ], ...
                'However, as detailed below, it does not fit with the requested values:');
    textgui2 = strvcat(txt{:});

    %textgui2 = ['wavelet cycles - [' num2str(params1.cycles(1)) ' ' num2str(params2.cycles(2)) ...
    %                '] instead of [' num2str(params1.cycles(1)) ' ' num2str(params2.cycles(2)) ...
    %            '], padratio - ' num2str(params1.padratio) ' instead of ' num2str(param2.padratio) ...
    %            ', and bootstrap significance - ' num2str(params1.alpha) ' instead of ' num2str(params2.alpha) ];

    uilist = { {'style' 'text' 'string' textgui1 } ...
               {'style' 'text' 'string' textgui2 } {} ...
               {'style' 'text' 'string' ['Would you like to recompute ' upper(strcom) ' and overwrite those values?' ]} ...
               {'style' 'checkbox' 'tag' 'ersp_yes' 'string' 'Yes, recompute' 'value' 1 'Callback' set_yes }  ...
               {'style' 'checkbox' 'tag' 'ersp_no' 'string' 'No, use data file parameters (and existing ERSP info)' 'value' 0 'Callback' set_no } };

    ersp_ans = inputgui('geometry', {[1] [1] [1] [1] [0.5 1] }, 'geomvert', [2 max(length(txt),1)*0.7 1 1 1 1], 'uilist', uilist, ...
                        'helpcom', '', 'title', ['Recalculate ' upper(strcom) ' parameters -- part of std_ersp()']);

    if isempty(ersp_ans), res = 'cancel'; return; end
    if find(celltomat(ersp_ans))  == 2 % use existing ERSP info from this dataset
        disp(['Using file on disk: ' filename ]);
        params2 = params1;
        res              = 'usedisk';
    else % Over write data in dataset
        disp(['Deleting and recomputing file: ' filename ]);
        res              = 'recompute';
        %delete(filename);
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [ o, x, w ] = epn_glg_00_1 ( n, alpha )

%*****************************************************************************80
%
%% EPN_GLG_00_1 implements the "midpoint rule" for region EPN_GLG.
%
%  Discussion:
%
%    The rule has order O = 1.
%
%    The rule has precision P = 0.
%
%    EPN_GLG is the N-dimensional positive space [0,+oo)^N with generalized
%    Laguerre weight function:
%
%      w(alpha;x) = product ( 1 <= i <= n ) x(i)^alpha exp ( - x(i) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Input, real ALPHA, the exponent of X in the weight function.
%    -1.0 < ALPHA.
%
%    Input, integer O, the order.
%
%    Output, real X(N,O), the abscissas.
%
%    Output, real W(O), the weights.
%
  if ( alpha <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'EPN_GLG_00_1 - Fatal error!\n' );
    fprintf ( 1, '  ALPHA <= -1.0\n' );
    error ( 'EPN_GLG_00_1 - Fatal error!' );
  end

  expon = 0;
  volume = ep1_glg_monomial_integral ( expon, alpha );
  volume = volume ^ n;

  o = 1;

  x = zeros ( n, o );
  w = zeros ( o, 1 );

  k = 0;
%
%  1 point.
%
  k = k + 1;
  x(1:n,k) = 1.0;
  w(k) = volume;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% ------------------------------------------------------
% SwarmOps - Heuristic optimization for Matlab
% Copyright (C) 2003-2010 Magnus Erik Hvass Pedersen.
% Please see the file license.txt for license details.
% SwarmOps on the internet: http://www.Hvass-Labs.org/
% ------------------------------------------------------

% Create data-struct for Griewank problem.
% Parameters:
%     dim; the dimensionality of the search-space, e.g. 10.
%     maxEvaluations; the maximum number of fitness evaluations
%                     to perform in optimization.
% Returns:
%     data; the data-struct.
function data = griewankdata(dim, maxEvaluations)
    data = benchmarkdata(dim, 0.001, maxEvaluations, 300, 600, -600, 600);

    % Array used by griewank function to save time.
    data.steps = sqrt([1:dim]);
end

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

author: Sable

repo: mcbench-benchmarks


text:
function setforsave(h,linewidthvalue, fontsize)
% setforsave(h,linewidthvalue)
% Sets paramaters of the graph for saving. (linewidt, fontsize, minorticks)
if ~exist('fontsize','var')
    fontsize = 20;
end
if ~exist('linewidthvalue','var')
    linewidthvalue = 2;
end
if ~exist('h', 'var')
    h = gcf;
elseif isempty(h)
    h=gcf;
end

gh=get(h);
for ii=1:length(gh.Children)
    ga{ii} = get(gh.Children(ii));
end
figure (h)
setfontsizefigure(12)
for ii=1:length(gh.Children)
set(ga{ii}.Children, 'linewidth',linewidthvalue)

end
set(gh.Children,'XMinorTick','on', 'YMinorTick','on')
setfontsizefigure(fontsize)
grid on

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef GraphTrainLog
    properties
        nCost = 1;
        cost;
        cost_cv;
        cost_cv_hour;

        costMinibatch;
        costBlock;

        nHourSeen = 0;
        actual_LR = [];
    end

    methods
        function obj = GraphTrainLog(sigGraph)
            obj.nCost = length(sigGraph.costWeight);
            obj.cost = GraphCost(obj.nCost, 0);
            obj.cost_cv = GraphCost(obj.nCost, 0);
            obj.cost_cv_hour = GraphCost(obj.nCost, 0);
            obj.costMinibatch = GraphCost(obj.nCost, 0);
            obj.costBlock = GraphCost(obj.nCost, 0);
        end

        function obj = initialize(obj)
            obj.actual_LR = [];
        end

        function obj = initializeBlockCost(obj)
            obj.costBlock = GraphCost(obj.nCost, 0);
        end

        function obj = initializeMiniBatchCost(obj,nCost, nBatch, useGPU)
            obj.costMinibatch = GraphCost(nCost, nBatch);
        end

        function obj = accumulateMiniBatchCost(obj, batch_i, cost_func)
            obj.costMinibatch = cost_func.copyCost(obj.costMinibatch, batch_i);
        end

        function obj = accumulateBlockCost(obj)
            obj.costBlock = obj.costMinibatch.appendCost(obj.costBlock);
        end

        function obj = accumulateItrCost(obj)
            obj.cost = obj.costBlock.appendMeanCost(obj.cost);
        end

        function obj = accumulateTrainDataAmount(obj, minibatch, data)
            for i=1:length(minibatch)
                [D,D2,T,N] = size(minibatch{i});
                if data.streams(i).frameRate==0
                    frameRate = 1;  % frameRate=0 means one data point per sequence. assume the sequence is 1s long.
                else
                    frameRate = data.streams(i).frameRate;
                end
                duration(i) = T*N/frameRate;
            end
            obj.nHourSeen = obj.nHourSeen + max(duration) / 3600;
        end

        function SaveModel(obj, itr, para, sigGraph)
            if mod(itr,para.saveModelEveryXIter)==0
                modelfile = [para.output sprintf('.itr%d.LR%s', itr, FormatFloat4Name(learning_rate))];
                if obj.cost_pure_cv(end)>0.1
                    modelfile = sprintf('%s.CV%2.3f.mat', modelfile, cost_pure_cv(end));
                else
                    modelfile = sprintf('%s.CV%s.mat', modelfile, FormatFloat4Name(cost_pure_cv(end)));
                end
                LOG = obj;
                SaveDNN(sigGraph, para, LOG, modelfile);
            end
        end

        function DisplayCostMinibatch(obj, itr, blk_i, nBlock, batch_i, nMiniBatch, displayInterval, displayTag)
            if itr<=0
                ItrStr = 'CV cost';
            else
                ItrStr = sprintf('Cost at itr %i', itr);
            end

            if mod(batch_i,displayInterval)==0 || (batch_i==nMiniBatch && displayInterval>batch_i)
                costMB = obj.costMinibatch;
                meanTotalCost = mean(costMB.totalCost(max(1,batch_i-displayInterval+1) : batch_i));
                if isnan(meanTotalCost)
                    pause(0.1);
                end
                fprintf('%s, blk %d/%d, MB %d/%d = %.4g / %.4g', ...
                    ItrStr, blk_i, nBlock, batch_i, nMiniBatch, ...
                    meanTotalCost, ...
                    mean(costMB.taskCost(max(1,batch_i-displayInterval+1) : batch_i)));
                if obj.nCost>1
                    fprintf(' - Subcosts: ');
                    for i = 1:obj.nCost
                        fprintf('%d) %f; ', i, mean(costMB.subCost(i,max(1,batch_i-displayInterval+1) : batch_i)));
                    end
                end
                fprintf(' - %s - %s\n', datestr(now), displayTag);
                pause(0.01);
            end
        end

        function DisplayCostIteration(obj, itr)
            if itr<=0
                ItrStr = 'CV cost';
            else
                ItrStr = sprintf('Cost at itr %i', itr);
            end

            fprintf('%s = %f / %f\t - %s\n', ...
                ItrStr, obj.cost.totalCost(end), sum(obj.cost.taskCost(:,end)), datestr(now));
            if length(obj.cost.subCost(:,end))>1
                fprintf('        ----subcosts = %f\n', obj.cost.subCost(:,end));
            end
            if sum(obj.cost.subAcc(:,end)>=0)
                fprintf('    ----(Sub)task accuracy = %2.2f%%\n', 100*obj.cost.subCost(:,end));
            end
        end
    end

    methods

    end

end

author: singaxiong

repo: SignalGraph


text:
classdef ShFunc_Compliance < ShFunWithElasticPdes

    properties (Access = private)
        compliance
        fieldToPrint
        adjointProblem
        gradientGauss
    end

    methods (Access = public)

        function obj = ShFunc_Compliance(cParams)
            cParams.filterParams.quadratureOrder = 'LINEAR';
            obj.init(cParams);
            fileName = cParams.femSettings.fileName;
            obj.createEquilibriumProblem(fileName);
            obj.createOrientationUpdater();
        end

        function fP = addPrintableVariables(obj)
            phy = obj.getPdesVariablesToPrint();
            fP{1}.value = phy{1};
            fP{2}.value = obj.compliance/obj.value0;
            fP{3}.value = obj.designVariable.alpha;
            fP{4}.value = abs(obj.designVariable.alpha);
            fP{5}.value = permute(obj.gradientGauss,[3 1 2]);
            fP = obj.addHomogVariables(fP);
        end

        function v = getVariablesToPlot(obj)
            v{1} = obj.value*obj.value0;
        end

        function t = getTitlesToPlot(obj)
            t{1} = 'Compliance non scaled';
        end

        function fP = createPrintVariables(obj)
            types = {'Elasticity','ScalarGauss','VectorGauss'...
                        'VectorGauss','VectorGauss'};
            names = {'Primal','ComplianceGauss','AlphaGauss',...
                        'AlphaAbsGauss','GradientGauss'};
            fP = obj.obtainPrintVariables(types,names);
            fP = obj.addHomogPrintVariablesNames(fP);
        end

        function [fun, funNames] = getFunsToPlot(obj)
            mesh = obj.designVariable.mesh;
            phy = obj.physicalProblem;
            strain = phy.strainFun{1}; % !!!
            stress = phy.stressFun{1}; % !!!
            displ  = phy.uFun{1}; % !!!
            compl  = obj.compliance/obj.value0;

            quad = Quadrature.set(mesh.type);
            quad.computeQuadrature('LINEAR');

            aa.mesh       = mesh;
            aa.quadrature = quad;
            aa.fValues    = permute(compl, [3 2 1]);
            complFun = FGaussDiscontinuousFunction(aa);

            bb.mesh    = mesh;
            bb.fValues = obj.designVariable.alpha';
            alphaFun = P0Function(bb);

            fun      = {complFun, strain, stress, displ};
            funNames = {'compliance', 'strain', 'stress', 'u'};

            cc.mesh     = mesh;
            cc.filename = 'shfunc_compliance';
            cc.fun      = fun;
            cc.funNames = funNames;
%             pvPst = ParaviewPostprocessor(cc);
%             pvPst.print();
%             fp = FunctionPrinter(cc);
%             fp.print();
        end
    end

    methods (Access = protected)

        function solveState(obj)
            obj.physicalProblem.setC(obj.homogenizedVariablesComputer.C) % (:,:,7200,4); cmat
%             obj.physicalProblem.computeVariables();
            obj.physicalProblem.solve();
        end

        function solveAdjoint(obj)
            obj.adjointProblem = obj.physicalProblem;
        end

        function computeFunctionValue(obj)
            phy = obj.physicalProblem;
            dvolum = phy.getDvolume()';
            stress = phy.variables.stress;
            strain = phy.variables.strain;
            ngaus  = size(strain,1);
            nelem  = size(strain,3);

            c = zeros(nelem,ngaus);
            for igaus = 1:ngaus
                stressG = squeeze(stress(igaus,:,:));
                strainG  = squeeze(strain(igaus,:,:));
                e = stressG.*strainG;
                c(:,igaus) = c(:,igaus) + sum(e)';
            end
            obj.compliance = c;
            int = c.*dvolum;
            obj.value = sum(int(:));
        end

        function computeGradientValue(obj)
            obj.computeGradientInGauss();
            obj.gradient = obj.gradientGauss;
        end

        function computeGradientInGauss(obj)
            phy = obj.physicalProblem;
            ep    = phy.variables.strain;
            ngaus  = size(ep,1);
            nstre  = size(ep,2);
            nelem  = size(ep,3);
            g = zeros(nelem,ngaus,obj.nVariables);
            for igaus = 1:ngaus
                for istre = 1:nstre
                    for jstre = 1:nstre
                        eu_i = squeeze(ep(igaus,istre,:));
                        ep_j = squeeze(ep(igaus,jstre,:));
                        for ivar = 1:obj.nVariables
                            dCij = squeeze(obj.homogenizedVariablesComputer.dC(istre,jstre,ivar,:,igaus));
                            g(:,igaus,ivar) = g(:,igaus,ivar) + (-eu_i.*dCij.*ep_j);
                        end
                    end
                end
            end
            obj.gradientGauss = g;
        end

        function f = getPdesVariablesToPrint(obj)
            f{1} = obj.getPdeVariableToPrint(obj.physicalProblem);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text: %SUM Datafile overload

author: marianux

repo: ecg-kit

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


text:
function [X,T] = readOff(filename)

fid = fopen(filename,'r');
if( fid==-1 )
    error('Cannot open the file.');
    return;
end

str = fgets(fid);   % -1 if eof

if strcmp(str(1:4), 'COFF')
    [X,T,~] = readCoff(filename,4); % assume 4 color channels
    return;
end

if ~strcmp(str(1:3), 'OFF')
    error('The file is not a valid OFF one.');
end

str = fgets(fid);
sizes = sscanf(str, '%d %d', 2);
while length(sizes) ~= 2
    str = fgets(fid);
    sizes = sscanf(str, '%d %d', 2);
end
nv = sizes(1);
nf = sizes(2);

% Read vertices
[X,cnt] = fscanf(fid,'%lf %lf %lf\n', [3,nv]);
if cnt~=3*nv
    warning('Problem in reading vertices.');
end
X = X';

[T,cnt] = fscanf(fid,'3 %ld %ld %ld\n', [3,inf]);
T = T'+1;

fclose(fid);

author: gpeyre


text:
function title = p04_title ( )

%*****************************************************************************80
%
%% P04_TITLE returns the title of problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'The jump data, 0 0 0 0 0 1/2 1 1 1 1 1';

  return
end


text:
function T = real(T)

T.M = real(T.M);

author: mtex-toolbox

repo: mtex

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


text:
function f_x = ParFor10(in1)
%PARFOR10
%    F_X = PARFOR10(IN1)

%    This function was generated by the Symbolic Math Toolbox version 8.1.
%    16-Jul-2019 11:57:29

u = in1(:,1);
uxx = in1(:,5);
f_x = (uxx.*-1.57974309709607)./(u-1.4728218492819);

author: dynamicslab

repo: SINDy-PI


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   [alpha, beta, gamma] = rot2euler(R, convention)
%   Returns the Euler angles alpha, beta and gamma that yield a given matrix
%   R. The convention specifies the order and axes of rotations.
%   Currently, only the XYZ convention is supported.
%
%   In particular, the function computes the angles alpha, beta and gamma that
%   allow to compute R as.
%
%   R = Rot(alpha,'x')*Rot(beta,'y')*Rot(gamma,'z')
%
%   Author: Arturo Gil. Universidad Miguel Hernandez de Elche. email:
%   arturo.gil@umh.es date:   11/11/2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function [sol1, sol2] = rot2euler(R, convention)

if convention=='XYZ'
    [sol1, sol2]=conventionXYZ(R);
elseif convention=='ZYX'
    [sol1, sol2]=conventionZYX(R);
else
    'Unknown convention. Only XYZ and ZYX are supported'
end

function [sol1, sol2] = conventionXYZ(R)
%'R(1,3)=sen(beta)=1??'
if abs(R(1,3)) == 1
    % degenerate case in which sen(beta)=+-1 and cos(beta)=0
    alpha1 = 0; % arbitrarily set alpha to zero
    alpha2 = pi; % arbitrarily set alpha to pi
    beta1 = asin(R(1,3));
    beta2 = pi-beta1;
    if sin(beta1) > 0
        gamma1 = atan2(R(2,1), -R(3,1));
        gamma2 = atan2(R(2,1), -R(3,1)) - alpha2;
    else
        gamma1 = atan2(R(2,1), R(3,1));
        gamma2 = atan2(R(2,1), R(3,1))+alpha2;
    end
else
    beta1 = asin(R(1,3));
    beta2 = pi-beta1;

    % standard way to compute alpha beta and gamma
    alpha1 = -atan2(R(2,3)/cos(beta1), R(3,3)/cos(beta1));
    alpha2 = -atan2(R(2,3)/cos(beta2), R(3,3)/cos(beta2));
    gamma1 = -atan2(R(1,2)/cos(beta1), R(1,1)/cos(beta1));
    gamma2 = -atan2(R(1,2)/cos(beta2), R(1,1)/cos(beta2));
end
% Normalize all angles to -pi, pi
sol1 = [alpha1, beta1, gamma1];
sol2 = [alpha2, beta2, gamma2];
sol1 = normalize_angles(sol1);
sol2 = normalize_angles(sol2);

function [sol1, sol2] = conventionZYX(R)
%R(3,1)=sin(beta)=+-1?? --> degenerate case
% degenerate case in which sen(beta)=+-1 and cos(beta)=0
if abs(R(3,1)) == 1
    alpha1 = 0; % arbitrarily set alpha to zero
    alpha2 = pi;%arbitrarily set to pi
    beta1 = asin(-R(3,1));
    beta2 = pi-beta1;
    if sin(beta1) > 0
        gamma1 = atan2(R(1,2), R(2,2));
        gamma2 = atan2(R(1,2), R(2,2)) + alpha2;
    else
        gamma1 = atan2(-R(1,2), R(2,2));
        gamma2 = atan2(-R(1,2), R(2,2)) - alpha2;
    end
else
    % standard way to compute alpha beta and gamma outside of the gimbal
    % lock
    beta1 = asin(-R(3,1));
    beta2 = pi-beta1;

    alpha1 = atan2(R(2,1)/cos(beta1), R(1,1)/cos(beta1));
    alpha2 = atan2(R(2,1)/cos(beta2), R(1,1)/cos(beta2));
    gamma1 = atan2(R(3,2)/cos(beta1), R(3,3)/cos(beta1));
    gamma2 = atan2(R(3,2)/cos(beta2), R(3,3)/cos(beta2));
end
% Normalize all angles to -pi, pi
sol1 = [alpha1, beta1, gamma1];
sol2 = [alpha2, beta2, gamma2];
sol1 = normalize_angles(sol1);
sol2 = normalize_angles(sol2);

function sol_norm = normalize_angles(sol)
sol_norm = [0 0 0];
for i=1:3
    sol_norm(i) = atan2(sin(sol(i)), cos(sol(i)));
end

repo: ARTE


text:
% UMD_SRC Creates a source for the UMD Texture dataset
%
% Usage
%    src = UMD_SRC(directory)
%
% Input
%    directory: The directory containing the UMD Texture dataset.
%
% Output
%    src: The UMD source.
%

function src = umd_src(directory)
	if (nargin<1)
		directory = '/Users/laurentsifre/TooBigForDropbox/Databases/umd';
	end
	src = create_src(directory, @uiuc_extract_objects_fun);
end

function [objects, classes] = uiuc_extract_objects_fun(file)
	objects.u1 = [1, 1];
	objects.u2 = [960, 1280];
	path_str = fileparts(file);
	path_parts = regexp(path_str, filesep, 'split');
	classes = {path_parts{end}};
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function matrix = genMatCurlVIFE3D(pde,mesh,fem,meshI,femI,option)
%% Generate global matrices and load vector of FEM for 3D ellipic eq
%     -div(A grad u)  = f,    x\in \Omega
% INPUTS:
% pde --- given data function from equation, e.g.
%         pde.A --- diffusion coefficient
%         pde.f --- right hand side function
%         pde.gD --- Dirichlet boundary value function
%         pde.one --- constant function 1.
% mesh --- mesh structure.
% fem --- global degree of freedom of FEM
%
% OUTPUTS:
% matrix.S --- stiffness matrix (w/o boundary condition)
% matrix.A --- final FEM matrix (after boundary condition)
% matrix.rhsF --- load vector (w/o boundary condition)
% matrix.f --- final RHS matrix (after boundary condition)

% Last Modified: 08/07/2020 by Xu Zhang

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g2ldofNint = femI.g2ldofNint;
node = meshI.node;
gdof = femI.gdof;
NEdof = size(gdof,1);

% %% 1. Stiffness Matrix
feEvalBas1 = @EvalNed1Bas3D;
feEvalBas2 = @EvalNed1Bas3D;

dof1 = 6; dof2 = 6; nloc = dof1*dof2;
ntID = find(mesh.tLoc > 0); ntN = length(ntID);

%% 1. Matrix on noninterface elements
MatInd = 1;
AN = fem.area(ntID);
gxN = fem.gx(ntID,:); gyN = fem.gy(ntID,:); gzN = fem.gz(ntID,:); gw = fem.gw;
XN = zeros(nloc*ntN, 1);

coefN = feval(pde.A,gxN,gyN,gzN);
Ibasx = cell(dof1,1); Ibasy = cell(dof1,1); Ibasz = cell(dof1,1);
Jbasx = cell(dof2,1); Jbasy = cell(dof2,1); Jbasz = cell(dof2,1);
t_e_orit = fem.t_e_orit(ntID,:);

for i = 1:dof1
    Ibasx{i} = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, MatInd, 1).*t_e_orit(:,i);
    Ibasy{i} = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, MatInd, 2).*t_e_orit(:,i);
    Ibasz{i} = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, MatInd, 3).*t_e_orit(:,i);
end
for j = 1:dof2
    Jbasx{j} = feEvalBas2(fem.bas, ntID, gxN, gyN, gzN, j, MatInd, 1).*t_e_orit(:,j);
    Jbasy{j} = feEvalBas2(fem.bas, ntID, gxN, gyN, gzN, j, MatInd, 2).*t_e_orit(:,j);
    Jbasz{j} = feEvalBas2(fem.bas, ntID, gxN, gyN, gzN, j, MatInd, 3).*t_e_orit(:,j);
end

IN = reshape(repmat(g2ldofNint(:,1:6),6,1),nloc*ntN,1);
JN = repmat(reshape(g2ldofNint(:,1:6),dof2*ntN,1),6,1);
ind = 0;
for i = 1:dof1
    for j = 1:dof2
        XN(ind+1:ind+ntN) = AN.*(sum(((Ibasx{i}.*(coefN.*Jbasx{j})).*gw'),2) + ...
            sum(((Ibasy{i}.*(coefN.*Jbasy{j})).*gw'),2) + ...
            sum(((Ibasz{i}.*(coefN.*Jbasz{j})).*gw'),2));
        ind = ind + ntN;
    end
end
ID = find(XN~=0);
SN = sparse(IN(ID),JN(ID),XN(ID),NEdof,NEdof);

%%%%%%

MatInd = 0;
AN = fem.area(ntID);
gxN = fem.gx(ntID,:); gyN = fem.gy(ntID,:); gzN = fem.gz(ntID,:); gw = fem.gw;
XN = zeros(nloc*ntN, 1);

coefN = feval(pde.B,gxN,gyN,gzN);
Ibasx = cell(dof1,1); Ibasy = cell(dof1,1); Ibasz = cell(dof1,1);
Jbasx = cell(dof2,1); Jbasy = cell(dof2,1); Jbasz = cell(dof2,1);

for i = 1:dof1
    Ibasx{i} = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, MatInd, 1).*t_e_orit(:,i);
    Ibasy{i} = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, MatInd, 2).*t_e_orit(:,i);
    Ibasz{i} = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, MatInd, 3).*t_e_orit(:,i);
end
for j = 1:dof2
    Jbasx{j} = feEvalBas2(fem.bas, ntID, gxN, gyN, gzN, j, MatInd, 1).*t_e_orit(:,j);
    Jbasy{j} = feEvalBas2(fem.bas, ntID, gxN, gyN, gzN, j, MatInd, 2).*t_e_orit(:,j);
    Jbasz{j} = feEvalBas2(fem.bas, ntID, gxN, gyN, gzN, j, MatInd, 3).*t_e_orit(:,j);
end

IN = reshape(repmat(g2ldofNint(:,1:6),6,1),nloc*ntN,1);
JN = repmat(reshape(g2ldofNint(:,1:6),dof2*ntN,1),6,1);
ind = 0;
for i = 1:dof1
    for j = 1:dof2
        XN(ind+1:ind+ntN) = AN.*(sum(((Ibasx{i}.*(coefN.*Jbasx{j})).*gw'),2) + ...
            sum(((Ibasy{i}.*(coefN.*Jbasy{j})).*gw'),2) + ...
            sum(((Ibasz{i}.*(coefN.*Jbasz{j})).*gw'),2));
        ind = ind + ntN;
    end
end
ID = find(XN~=0);
MN = sparse(IN(ID),JN(ID),XN(ID),NEdof,NEdof);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dof = 6;  nloc = dof;
X = zeros(nloc*ntN, 1);

fN1 = feval(pde.f1,gxN,gyN,gzN);
fN2 = feval(pde.f2,gxN,gyN,gzN);
fN3 = feval(pde.f3,gxN,gyN,gzN);
ind = 0;
I = reshape(g2ldofNint(:,1:6),nloc*ntN,1);
for i = 1:dof
    ibas1 = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, 0, 1);
    ibas2 = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, 0, 2);
    ibas3 = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, 0, 3);
    X(ind+1:ind+ntN) = AN.*sum((ibas1.*fN1+ibas2.*fN2+ibas3.*fN3).*gw',2).*t_e_orit(:,i);
    ind = ind + ntN;
end
rhsN = sparse(I,1,X,NEdof,1);

if isfield(option,'HM')
    if option.HM == 1
        Atotal = SN - pde.kappa*MN + femI.KI - pde.kappa*femI.SI;
    end
else
    Atotal = SN + MN + femI.KI + femI.SI;
end

rhs = rhsN + femI.b;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% 3. Dirichlet Boundary Conditions

% non-interface edges
vSign = meshI.vSign;
eSign = sum(vSign(femI.gdof),2)/2;
eIdm = find(eSign==-1/2); eIdp = find(eSign==1/2);
[gew,gex,gey,gez] = gaussPedge(node(gdof(:,1),:),node(gdof(:,2),:),2);
% tu1 = sum(feval(pde.exactu1,gex,gey,gez).*gew,2);
% tu2 = sum(feval(pde.exactu2,gex,gey,gez).*gew,2);
% tu3 = sum(feval(pde.exactu3,gex,gey,gez).*gew,2);
tu1 = sum(feval(pde.exactu1,gex,gey,gez).*gew,2);
tu2 = sum(feval(pde.exactu2,gex,gey,gez).*gew,2);
tu3 = sum(feval(pde.exactu3,gex,gey,gez).*gew,2);
tu1(eIdm) = sum(feval(pde.um1,gex(eIdm,:),gey(eIdm,:),gez(eIdm,:)).*gew(eIdm,:),2);
tu1(eIdp) = sum(feval(pde.up1,gex(eIdp,:),gey(eIdp,:),gez(eIdp,:)).*gew(eIdp,:),2);
tu2(eIdm) = sum(feval(pde.um2,gex(eIdm,:),gey(eIdm,:),gez(eIdm,:)).*gew(eIdm,:),2);
tu2(eIdp) = sum(feval(pde.up2,gex(eIdp,:),gey(eIdp,:),gez(eIdp,:)).*gew(eIdp,:),2);
tu3(eIdm) = sum(feval(pde.um3,gex(eIdm,:),gey(eIdm,:),gez(eIdm,:)).*gew(eIdm,:),2);
tu3(eIdp) = sum(feval(pde.up3,gex(eIdp,:),gey(eIdp,:),gez(eIdp,:)).*gew(eIdp,:),2);
tgt = node(gdof(:,2),:) - node(gdof(:,1),:);
tgt = tgt./sum(tgt.^2,2).^(1/2);
tu = tu1.*tgt(:,1) + tu2.*tgt(:,2) + tu3.*tgt(:,3);

bcind = fem.bcind;
[bc,mapper] = boundaryEdge3D(node,gdof,bcind);
bdidx = zeros(NEdof,1);
isBdEdge = true(NEdof,1);
isBdEdge(mapper) = false;
bdidx(isBdEdge) = 1;
Tbd = spdiags(bdidx,0,NEdof,NEdof);
T = spdiags(1-bdidx,0,NEdof,NEdof);
A = T*Atotal*T + Tbd;

ub = tu;
ub(mapper) = 0;
rhsB = Atotal*ub;
f = rhs - rhsB;
f(isBdEdge) = tu(isBdEdge);

%% Outputs
matrix = struct('A', A, 'f', f,'tu',tu,'isBdEdge',isBdEdge);

repo: ifem


text:
function y = mytestNEW(B,X)

y = B - X*X;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% Newton下山法 当初值选择不当时 调整下山因子lamda
clear;
format long;
tol = 1e-5;
N = 100;
x0 = 0.6;
lamda = 1;
f = @(x) x^3 - x - 1;  %f(x)表达式
df = @(x) 3*x^2 - 1;
fprintf('f(x)的初值: %d\n', abs(f(x0)));
for k = 1 : N
    x1 = x0 - f(x0)/ df(x0);
     while abs(f(x1)) > abs(f(x0))
         lamda = lamda / 2;   % 下山因子减半
         fprintf('下山因子减半后的值: %d\n', lamda);
         x1 = x0 - lamda * f(x0) / df(x0);
     end
    fprintf('本次迭代f(x)的值: %d\n', abs(f(x1)));
    if abs(x1 - x0) < tol
        fprintf('迭代次数: %d\n', k);
        fprintf('方程的正根: %10.8f\n', x1);
        break;
    end
    x0 = x1;
end
if k == N
    fprintf('迭代方法失败\n');
end

repo: NumericalAnalysis


text:
function [x cutvalue cutvalue_upperbound Y] = maxcut(L, r)
% Algorithm to (try to) compute a maximum cut of a graph, via SDP approach.
%
% function x = maxcut(L)
% function [x cutvalue cutvalue_upperbound Y] = maxcut(L, r)
%
% L is the Laplacian matrix describing the graph to cut. The Laplacian of a
% graph is L = D - A, where D is the diagonal degree matrix (D(i, i) is the
% sum of the weights of the edges adjacent to node i) and A is the
% symmetric adjacency matrix of the graph (A(i, j) = A(j, i) is the weight
% of the edge joining nodes i and j). If L is sparse, this will be
% exploited.
%
% If the graph has n nodes, then L is nxn and the output x is a vector of
% length n such that x(i) is +1 or -1. This partitions the nodes of the
% graph in two classes, in an attempt to maximize the sum of the weights of
% the edges that go from one class to the other (MAX CUT problem).
%
% cutvalue is the sum of the weights of the edges 'cut' by the partition x.
%
% If the algorithm reached the global optimum of the underlying SDP
% problem, then it produces an upperbound on the maximum cut value. This
% value is returned in cutvalue_upperbound if it is found. Otherwise, that
% output is set to NaN.
%
% If r is specified (by default, r = n), the algorithm will stop at rank r.
% This may prevent the algorithm from reaching a globally optimal solution
% for the underlying SDP problem (but can greatly help in keeping the
% execution time under control). If a global optimum of the SDP is reached
% before rank r, the algorithm will stop of course.
%
% Y is a matrix of size nxp, with p <= r, such that X = Y*Y' is the best
% solution found for the underlying SDP problem. If cutvalue_upperbound is
% not NaN, then Y*Y' is optimal for the SDP and cutvalue_upperbound is its
% cut value.
%
% By Goemans and Williamson 1995, it is known that if the optimal value of
% the SDP is reached, then the returned cut, in expectation, is at most at
% a fraction 0.878 of the optimal cut. (This is not exactly valid because
% we do not use random projection here; sign(Y*randn(size(Y, 2), 1)) will
% give a cut that respects this statement -- it's usually worse though).
%
% The algorithm is essentially that of:
% Journee, Bach, Absil and Sepulchre, 2010
% Low-rank optimization on the code of positive semidefinite matrices.
%
% It is itself based on the famous SDP relaxation of MAX CUT:
% Goemans and Williamson, 1995
% Improved approximation algorithms for maximum cut and satisfiability
% problems using semidefinite programming.

% This file is part of Manopt and is copyrighted. See the license file.
%
% Main author: Nicolas Boumal, July 18, 2013
% Contributors:
%
% Change log:
%

    % If no inputs are provided, generate a random Laplacian.
    % This is for illustration purposes only.
    if ~exist('L', 'var') || isempty(L)
        n = 20;
        A = triu(randn(n) <= .4, 1);
        A = A+A';
        D = diag(sum(A, 2));
        L = D-A;
    end

    n = size(L, 1);
    assert(size(L, 2) == n, 'L must be square.');

    if ~exist('r', 'var') || isempty(r) || r > n
        r = n;
    end

    % We will let the rank increase. Each rank value will generate a cut.
    % We have to go up in the rank to eventually find a certificate of SDP
    % optimality. This in turn will give us an upperbound on the MAX CUT
    % value and assure us that we're doing well, according to Goemans and
    % Williamson's argument. In practice though, the good cuts often come
    % up for low rank values, so we better keep track of the best one.
    best_x = ones(n, 1);
    best_cutvalue = 0;
    cutvalue_upperbound = NaN;

    time = [];
    cost = [];

    for rr = 2 : r

        manifold = elliptopefactory(n, rr);

        if rr == 2

            % At first, for rank 2, generate a random point.
            Y0 = manifold.rand();

        else

            % To increase the rank, we could just add a column of zeros to
            % the Y matrix. Unfortunately, this lands us in a saddle point.
            % To escape from the saddle, we may compute an eigenvector of
            % Sy associated to a negative eigenvalue: that will yield a
            % (second order) descent direction Z. See Journee et al ; Sy is
            % linked to dual certificates for the SDP.
            Y0 = [Y zeros(n, 1)];
            LY0 = L*Y0;
            Dy = spdiags(sum(LY0.*Y0, 2), 0, n, n);
            Sy = (Dy - L)/4;
            % Find the smallest (the "most negative") eigenvalue of Sy.
            [v, s] = eigs(Sy, 1, 'SA');
            % If there is no negative eigenvalue for Sy, than we are not at
            % a saddle point: we're actually done!
            if s >= -1e-8
                % We can stop here: we found the global optimum of the SDP,
                % and hence the reached cost is a valid upper bound on the
                % maximum cut value.
                cutvalue_upperbound = max(-[info.cost]);
                break;
            end

            % This is our escape direction.
            Z = manifold.proj(Y0, [zeros(n, rr-1) v]);

            % % These instructions can be uncommented to see what the cost
            % % function looks like at a saddle point. But will require the
            % % problem structure which is not defined here: see the helper
            % % function.
            % plotprofile(problem, Y0, Z, linspace(-1, 1, 101));
            % drawnow; pause;

            % Now make a step in the Z direction to escape from the saddle.
            % It is not obvious that it is ok to do a unit step ... perhaps
            % need to be cautious here with the stepsize. It's not too
            % critical though: the important point is to leave the saddle
            % point. But it's nice to guarantee monotone decrease of the
            % cost, and we can't do that with a constant step (at least,
            % not without a proper argument to back it up).
            stepsize = 1;
            Y0 = manifold.retr(Y0, Z, stepsize);

        end

        % Use the Riemannian optimization based algorithm lower in this
        % file to reach a critical point (typically a local optimizer) of
        % the max cut cost with fixed rank, starting from Y0.
        [Y info] = maxcut_fixedrank(L, Y0);

        % Some info logging.
        thistime = [info.time];
        if ~isempty(time)
            thistime = time(end) + thistime;
        end
        time = [time thistime]; %#ok<AGROW>
        cost = [cost [info.cost]]; %#ok<AGROW>

        % Time to turn the matrix Y into a cut.
        % We can either do the random rounding as follows:
        % x = sign(Y*randn(rr, 1));
        % or extract the "PCA direction" of the points in Y and cut
        % orthogonally to that direction, as follows:
        [u, ~, ~] = svds(Y, 1);
        x = sign(u);

        cutvalue = (x'*L*x)/4;
        if cutvalue > best_cutvalue
            best_x = x;
            best_cutvalue = cutvalue;
        end

    end

    x = best_x;
    cutvalue = best_cutvalue;

    plot(time, -cost, '.-');
    xlabel('Time [s]');
    ylabel('Relaxed cut value');
    title('The relaxed cut value is an upper bound on the optimal cut value.');

end

function [Y info] = maxcut_fixedrank(L, Y)
% Try to solve the (fixed) rank r relaxed max cut program, based on the
% Laplacian of the graph L and an initial guess Y. L is nxn and Y is nxr.

    [n r] = size(Y);
    assert(all(size(L) == n));

    % The fixed rank elliptope geometry describes symmetric, positive
    % semidefinite matrices of size n with rank r and all diagonal entries
    % are 1.
    manifold = elliptopefactory(n, r);

    % % If you want to compare the performance of the elliptope geometry
    % % against the (conceptually simpler) oblique manifold geometry,
    % % uncomment this line.
    % manifold = obliquefactory(r, n, true);

    problem.M = manifold;

    % % For rapid prototyping, these lines suffice to describe the cost
    % % function and its gradient and Hessian (here expressed using the
    % % Euclidean gradient and Hessian).
    % problem.cost = @(Y)  -trace(Y'*L*Y)/4;
    % problem.egrad = @(Y) -(L*Y)/2;
    % problem.ehess = @(Y, U) -(L*U)/2;

    % Instead of the prototyping version, the functions below describe the
    % cost, gradient and Hessian using the caching system (the store
    % structure). This alows to execute exactly the required number of
    % multiplications with the matrix L. These multiplications are counted
    % using the Lproducts_counter and registered for each iteration in the
    % info structure outputted by solvers, via the statsfun function.
    % Notice that we do not use the store structure to count: this does not
    % behave well in general and is not advised.

    Lproducts_counter = 0;

    % For every visited point Y, we will need L*Y. This function makes sure
    % the quantity L*Y is available, but only computes it if it wasn't
    % already computed.
    function store = prepare(Y, store)
        if ~isfield(store, 'LY')
            store.LY = L*Y;
            Lproducts_counter = Lproducts_counter + 1;
        end
    end

    problem.cost = @cost;
    function [f store] = cost(Y, store)
        store = prepare(Y, store);
        LY = store.LY;
        f = -(Y(:)'*LY(:))/4; % = -trace(Y'*LY)/4;
    end

    problem.grad = @grad;
    function [g store] = grad(Y, store)
        store = prepare(Y, store);
        LY = store.LY;
        g = manifold.egrad2rgrad(Y, -LY/2);
    end

    problem.hess = @hess;
    function [h store] = hess(Y, U, store)
        store = prepare(Y, store);
        LY = store.LY;
        LU = L*U;
        Lproducts_counter = Lproducts_counter + 1;
        h = manifold.ehess2rhess(Y, -LY/2, -LU/2, U);
    end

    % statsfun is called exactly once after each iteration (including after
    % the evaluation of the cost at the initial guess). We then register
    % the value of the Lproducts counter (which counts how many product
    % were needed since the last iteration), and reset it to zero.
    options.statsfun = @statsfun;
    function stats = statsfun(problem, Y, stats, store) %#ok
        stats.Lproducts = Lproducts_counter;
        Lproducts_counter = 0;
    end

    % % Diagnostics tools: to make sure the gradient and Hessian are
    % % correct during the prototyping stage.
    % checkgradient(problem); pause;
    % checkhessian(problem); pause;

    % % To investigate the effect of the rotational invariance when using
    % % the oblique or the elliptope geometry, or to study the saddle point
    % % issue mentioned above, it is sometimes interesting to look at the
    % % spectrum of the Hessian. For large dimensions, this is slow!
    % stairs(sort(hessianspectrum(problem, Y)));
    % drawnow; pause;

    % % When facing a saddle point issue as described in the master
    % % function, and when no sure mechanism exists to find an escape
    % % direction, it may be helpful to set useRand to true and raise
    % % miniter to more than 1, when using trustregions. This will tell the
    % % solver to not stop before at least miniter iterations were
    % % accomplished (thus disregarding the zero gradient at the saddle
    % % point) and to use random search directions to kick start the inner
    % % solve (tCG) step. It is not as efficient as finding a sure escape
    % % direction, but sometimes it's the best we have.
    % options.useRand = true;
    % options.miniter = 5;

    options.verbosity = 2;
    Lproducts_counter = 0;
    [Y Ycost info] = trustregions(problem, Y, options); %#ok

    % fprintf('Products with L: %d\n', sum([info.Lproducts]));

end

author: andrewssobral

repo: lrslibrary


text:
Ts_INS = 1/150;
fprintf('Sample time set: Ts_INS = %f seconds\n',Ts_INS);


text:
function strMacro = macro_method()
% MACRO_METHOD: Script to insert at the beginning of all the brainstorm class functions

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

strMacro = [...
    'if (nargin >= 1), ' ...
        '[varargout{1:nargout}] = feval(varargin{1}, varargin{2:end});', ...
    'end'];

author: brainstorm-tools


text:
%% Real Time Data Stream Plotting Example
function real_time_data_stream_plotting
%%
% This example demonstrates how to automatically read a set number of data bytes as and
% when they are available. This MATLAB(R) script also generates a real time plot of streaming
% data collected from the TCPIP server.
%
% The script may be updated to use any instrument/device/TCPIP server
% to collect real time data. You may need to update the IP address and
% port.
%
% To generate a report of this entire script, you may use the PUBLISH
% command at the MATLAB(R) command line as follows:
%
% publish(real_time_data_plot);

% Author: Ankit Desai
% Copyright 2010 - The MathWorks, Inc.

%% Create the interface object
% Create a TCPIP object listening to port 19 (Character Generator).
%
% *Note* : To enable character generator service at port 19 on a Windows platform, enable:
%
%  Control Panel > Add Remove Programs > Add/Remove Windows Component > Networking Services
%
interfaceObject = tcpip('localhost',19);

%%
% Setup a figure window and define a callback function for close operation
figureHandle = figure('NumberTitle','off',...
    'Name','Live Data Stream Plot',...
    'Color',[0 0 0],...
    'CloseRequestFcn',{@localCloseFigure,interfaceObject});

%%
% Setup the axes
axesHandle = axes('Parent',figureHandle,...
    'YGrid','on',...
    'YColor',[0.9725 0.9725 0.9725],...
    'XGrid','on',...
    'XColor',[0.9725 0.9725 0.9725],...
    'Color',[0 0 0]);

xlabel(axesHandle,'Number of Samples');
ylabel(axesHandle,'Value');

%%
% Initialize the plot and hold the settings on
hold on;
plotHandle = plot(axesHandle,0,'-y','LineWidth',1);

%% Setup interface object to read chunks of data
% Set the number of bytes to read at a time
bytesToRead = 500;

%%
% Define a callback function to be executed when desired number of bytes
% are available in the input buffer
interfaceObject.BytesAvailableFcn = {@localReadAndPlot,plotHandle,bytesToRead};
interfaceObject.BytesAvailableFcnMode = 'byte';
interfaceObject.BytesAvailableFcnCount = bytesToRead;

%%
% Open the interface object
fopen(interfaceObject);
pause(3);
snapnow;
%% Implement the bytes available callback
function localReadAndPlot(interfaceObject,~,figureHandle,bytesToRead)

%%
% Read the desired number of data bytes
data = fread(interfaceObject,bytesToRead);

%%
% Update the plot
set(figureHandle,'Ydata',data);

%% Implement the close figure callback
function localCloseFigure(figureHandle,~,interfaceObject)

%%
% Clean up the interface object
fclose(interfaceObject);
delete(interfaceObject);
clear interfaceObject;

%%
% Close the figure window
delete(figureHandle);


text:
% The COBRAToolbox: testSbmlTestModelToMat.m
%
% Purpose:
%     - tests the batch conversion of SBML models to .mat files.
%
% Authors:
%     - Original file: Jacek Wachowiak
%     - Updated to check models for similarity: Thomas Pfau - Sept 2017
%

global CBTDIR

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testSbmlTestModelToMat.m'));
cd(fileDir);

%Create Temporary Folders
SBMLFolder = tempname;
mkdir(SBMLFolder);

MATFolder = tempname;
mkdir(MATFolder);

%copy all xml files from the models to the temp folder
modeldir = getDistributedModelFolder('iIT341.xml');

modelfiles = dir(modeldir);

modelNames = cell(0);
models = cell(0);
%We will use up to 3 model files for the test.
for i = 1:size(modelfiles)
     if numel(models) >= 3
        break;
    end
    if ~isempty(regexp(modelfiles(i).name,'.*\.xml$'))
        copyfile([modeldir filesep, modelfiles(i).name],SBMLFolder);
        models{end+1} = readCbModel([modeldir filesep modelfiles(i).name]);
        modelNames{end+1} = modelfiles(i).name;
    end
end
sbmlTestModelToMat(SBMLFolder,MATFolder)
matModels = cell(0);
for i = 1:size(models,2)
    %We load them, as otherwise the ID would change to the mat file and
    %they would no longer be equivalent.
    load([MATFolder filesep strrep(modelNames{i},'.xml','.mat')]);
    matModels{end+1} = model;

end

%Test that the models are the same
assert(numel(matModels) == numel(models));
for i = 1:numel(matModels)
    assert(isSameCobraModel(matModels{i},models{i}));
end

%clear all the xml files from the temporary dir.
delete([SBMLFolder filesep '*.xml'])
delete([MATFolder filesep '*.mat'])

%create an invalid xml file (and test with the defaults)
tempfolder = tempname;
mkdir(tempfolder);
cd(tempfolder);
defaultname = 'm_model_collection';
mkdir(defaultname);

invalidFile = fopen([defaultname filesep 'Invalid.xml'],'w');
fprintf(invalidFile,'Not a vlid XML\n');
fclose(invalidFile);
%Test with defaults, and nothing works
sbmlTestModelToMat()
defaultfoldercontent = dir(defaultname);

% the invalid file, '.' and '..'
assert(size(defaultfoldercontent,1) == 3);

cd(fileDir);

%Clean up the folders. (In a try/catch for windows)
%Even if this does not work, those folders are just temporary folders.
try
   rmdir(MATFolder,'s');
   rmdir(SBMLFolder,'s');
   rmdir([tempfolder filesep defaultname],'s');
   isdeleted = true;
catch

end

cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function [q1New,q2New,dq1New,dq2New] = autoGen_heelStrike(q1,q2,dq1,dq2,m,I,d,l)
%AUTOGEN_HEELSTRIKE
%    [Q1NEW,Q2NEW,DQ1NEW,DQ2NEW] = AUTOGEN_HEELSTRIKE(Q1,Q2,DQ1,DQ2,M,I,D,L)

%    This function was generated by the Symbolic Math Toolbox version 6.2.
%    07-Sep-2015 19:06:42

q1New = q2;
if nargout > 1
    q2New = q1;
end
if nargout > 2
    t2 = d.^2;
    t3 = m.^2;
    t4 = l.^2;
    t5 = I.^2;
    t6 = t2.^2;
    t7 = q1-q2;
    t8 = cos(t7);
    t9 = t3.*t6;
    t10 = t2.*t3.*t4.*(3.0./2.0);
    t11 = I.*m.*t2.*2.0;
    t12 = I.*m.*t4.*2.0;
    t13 = q1.*2.0;
    t18 = q2.*2.0;
    t14 = t13-t18;
    t15 = cos(t14);
    t16 = t5+t9+t10+t11+t12-d.*l.*t2.*t3.*2.0-t2.*t3.*t4.*t15.*(1.0./2.0)-I.*d.*l.*m.*2.0;
    t17 = 1.0./t16;
    dq1New = t17.*(dq2.*t5+dq2.*t3.*t6+I.*dq2.*m.*t2.*2.0-I.*d.*dq2.*l.*m+I.*dq1.*m.*t4.*t8.*2.0-d.*dq2.*l.*t2.*t3+dq1.*t2.*t3.*t4.*t8-I.*d.*dq1.*l.*m.*t8-d.*dq1.*l.*t2.*t3.*t8);
end
if nargout > 3
    dq2New = t17.*(dq1.*t5+dq1.*t3.*t6+dq1.*t2.*t3.*t4.*3.0+I.*dq1.*m.*t2.*2.0+I.*dq1.*m.*t4.*2.0-I.*d.*dq1.*l.*m.*3.0-d.*dq1.*l.*t2.*t3.*3.0-d.*dq1.*l.*t3.*t4-dq2.*t2.*t3.*t4.*t8-dq1.*t2.*t3.*t4.*t15+I.*d.*dq2.*l.*m.*t8+d.*dq2.*l.*t2.*t3.*t8+d.*dq1.*l.*t3.*t4.*t15);
end

author: MatthewPeterKelly

repo: OptimTraj


text:
function x = r8vec_cheby1space ( n, a, b )

%*****************************************************************************80
%
%% R8VEC_CHEBY1SPACE creates a vector of Type 1 Chebyshev spaced values in [A,B].
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real A, B, the first and last entries.
%
%    Output, real X(N,1), a vector of Chebyshev spaced data.
%
  x = zeros ( n, 1 );

  if ( n == 1 )

    x(1) = ( a + b ) / 2.0;

  else

    for i = 1 : n

      theta = ( 2 * ( n - i ) + 1 ) * pi / ( 2 * n );

      c = cos ( theta );

      if ( mod ( n, 2 ) == 1 )
        if ( 2 * i - 1 == n )
          c = 0.0;
        end
      end

      x(i) = ( ( 1.0 - c ) * a  ...
             + ( 1.0 + c ) * b ) ...
             /   2.0;

    end

  end

  return
end


text:
function [locate_voi_min_x,locate_voi_max_x,locate_voi_min_y,locate_voi_max_y,locate_voi_min_z,locate_voi_max_z] = dicomrt_voiboundaries(dose_xmesh,dose_ymesh,dose_zmesh,VOI,voi2use,PatientPosition)
% dicomrt_voiboundaries(dose_xmesh,dose_ymesh,dose_zmesh,VOI,voi2use,PatientPosition)
%
% Returns the voxel index of the smallest and largest coordinate of the selected VOI in the specified Z slice.
%
% slice is the slice number where the mins and maxs will be calculated.
%      Slice number is relative to voiref.
% voiref is the voi of reference
% dose_xmesh,dose_ymesh, are coordinates of the center of the dose-pixel
% VOI is a cell array which contain the patients VOIs as read by dicomrt_loadvoi
% voi2use is a vector pointing to the number of VOIs to be used
%
% Example:
%
% [xmin,xmax,ymin,ymax]=dicomrt_voiboundariesZ(12,1,xmesh,xmesh,VOI,5)
%
% returns in xmin,xmax,ymin,ymax the pixel numbers which correspond to the edges of
% structure "5" (e.g. target) in slice number 12 with reference to structure number "1".
%
% See also dicomrt_voiboundaries, dicomrt_surfdose
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check input
[VOI_temp]=dicomrt_checkinput(VOI);
VOI=dicomrt_varfilter(VOI_temp);

nslices=size(VOI{voi2use,2},1);

voiref=1;

for kk=1:nslices
    [locate_voi_min_1,locate_voi_max_1,locate_voi_min_2,locate_voi_max_2]=dicomrt_voiboundaries_single(kk,dose_xmesh,dose_ymesh,VOI_temp,voi2use,PatientPosition);
    if kk==1
        locate_voi_min_x=locate_voi_min_1;
        locate_voi_max_x=locate_voi_max_1;
        locate_voi_min_y=locate_voi_min_2;
        locate_voi_max_y=locate_voi_max_2;
    else
        if locate_voi_min_1 < locate_voi_min_x
            locate_voi_min_x=locate_voi_min_1;
        end
        if locate_voi_max_1 > locate_voi_max_x
            locate_voi_max_x=locate_voi_max_1;
        end
        if locate_voi_min_2 < locate_voi_min_y
            locate_voi_min_y=locate_voi_min_2;
        end
        if locate_voi_max_2 > locate_voi_max_y
            locate_voi_max_y=locate_voi_max_2;
        end
    end
end

% z mesh is always sorted
min_z_voi=VOI_temp{2,1}{voi2use,2}{1}(1,3);
max_z_voi=VOI_temp{2,1}{voi2use,2}{end}(1,3);

locate_voi_min_z=dicomrt_findpointVECT(dose_zmesh,min_z_voi);
locate_voi_max_z=dicomrt_findpointVECT(dose_zmesh,max_z_voi);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function G = curl(F)
%CURL Curl of a BALLFUNV.
%   G = CURL(F) returns the BALLFUNV of the curl of F.
%
% See also DIV.

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

if isempty( F )
    G = ballfunv();
    return
end

% Extract the components of a BALLFUNV:
Fc = F.comp;

% Formula for 3D curl
G = [ diff(Fc{3},2) - diff(Fc{2},3); ...
      diff(Fc{1},3) - diff(Fc{3},1); ...
      diff(Fc{2},1) - diff(Fc{1},2) ];
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [error_x, error_y, fex, fey, ae] = stokespost_q1q1_p(q1q1sol,jmpx,jmpy,els,xy,ev)
%stokespost_q1q1_p  computes Poisson error estimator for Q1-Q1
%   [err_x, err_y, fex, fey, ae] = stokespost_q1p0_p(q1q1sol,jmpx,jmpy,els,xy,ev);
%   input
%          q1q1sol        Q1-Q1 flow solution
%          jmpx, jmpy     component elementwise edge stress jumps
%          els            elementwise edge lengths
%          xy             vertex coordinate vector
%          ev             element mapping matrix
%   output
%          err_x, err_y   component of velocity elementwise error estimate
%          fex, fey       component elementwise rhs vectors
%          ae             elementwise Poisson problem matrices
%
%   IFISS function: DJS; 8 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
      fprintf('computing local error estimator... ')
      x=xy(:,1); y=xy(:,2);nvtx=length(x);
      nel=length(ev(:,1));
      error_x=zeros(nel,1);  error_y=zeros(nel,1);
      psol=q1q1sol(2*nvtx+1:end);
%
% set up 3x3 Gauss points
      gpt=sqrt(0.6);
      s(1) = -gpt; t(1) = -gpt; wt(1)=25/81;
      s(2) =  gpt; t(2) = -gpt; wt(2)=25/81;
      s(3) =  gpt; t(3) =  gpt; wt(3)=25/81;
      s(4) = -gpt; t(4) =  gpt; wt(4)=25/81;
      s(5) =  0.0; t(5) = -gpt; wt(5)=40/81;
      s(6) =  gpt; t(6) =  0.0; wt(6)=40/81;
      s(7) =  0.0; t(7) =  gpt; wt(7)=40/81;
      s(8) = -gpt; t(8) =  0.0; wt(8)=40/81;
      s(9) =  0.0; t(9) =  0.0; wt(9)=64/81;
%
% inner loop over elements
        for ivtx = 1:4
        xl_v(:,ivtx) = x(ev(:,ivtx));
        yl_v(:,ivtx) = y(ev(:,ivtx));
        psol_v(:,ivtx) = psol(ev(:,ivtx));
		end

        ae = zeros(nel,5,5); elerrx=zeros(5,nel);  elerry=zeros(5,nel);
        fex = zeros(nel,5); fey = zeros(nel,5);
% loop over Gauss points
         for igpt = 1:9
         sigpt=s(igpt);
         tigpt=t(igpt);
         wght=wt(igpt);
% evaluate derivatives etc
         [jac_v,invjac_v,phi_v,dphidx_v,dphidy_v] = deriv(sigpt,tigpt,xl_v,yl_v);
         [psi_v,dpsidx_v,dpsidy_v] = qderiv(sigpt,tigpt,xl_v,yl_v);

            for j = 1:5
               for i = 1:5
               ae(:,i,j) = ae(:,i,j)+wght*dpsidx_v(:,i+4).*dpsidx_v(:,j+4).*invjac_v(:);
               ae(:,i,j) = ae(:,i,j)+wght*dpsidy_v(:,i+4).*dpsidy_v(:,j+4).*invjac_v(:);
               end
               for ss=1:4
               fex(:,j) = fex(:,j)-wght*dphidx_v(:,ss).*psol_v(:,ss).*psi_v(:,j+4);
               fey(:,j) = fey(:,j)-wght*dphidy_v(:,ss).*psol_v(:,ss).*psi_v(:,j+4);
               end
            end
% end of Gauss point loop
         end
%
% include edge jumps (evaluated at the midpoint)
         for ee = 1:4
         fex(:,ee) = fex(:,ee) - jmpx(:,ee) .* els(:,ee)*(1/3);
         fey(:,ee) = fey(:,ee) - jmpy(:,ee) .* els(:,ee)*(1/3);
         end
%
% solve for local estimate
         for ielem = 1:nel
		 elerrx(:,ielem) = squeeze(ae(ielem,1:5,1:5))\(fex(ielem,1:5)');
		 elerry(:,ielem) = squeeze(ae(ielem,1:5,1:5))\(fey(ielem,1:5)');
	     end
%%
         for ivtx=1:5,
	     error_x(:) = error_x(:) + fex(:,ivtx) .* elerrx(ivtx,:)';
		 error_y(:) = error_y(:) + fey(:,ivtx) .* elerry(ivtx,:)';
         end
%%
         fprintf('done.\n')
		 return

author: johannesgerer

repo: jburkardt-m


text:
function V = scn_write_plane(filenames_or_V, dat, wh_slice, varargin)
% Images can be 3D or 4D (4D is possible with SPM2+ using ,xx indexing, or SPM5+).
%
% :Usage:
% ::
%
%     V = scn_write_plane(filenames_or_V, dat, wh_slice, [exampleV])
%
% :Inputs:
%
%   **ofilenames_or_V:**
%        is string matrix or cell array of names, or structures created with spm_vol(names)
%
%   **dat:**
%        is 3-D array of data, vox x vox x slices, with data on 3rd dim being
%        written to separate image files
%
%   **wh_slice:**
%        slice number (voxel space)
%
%   **exampleV:**
%        is optional, but required if filenames are passed in.
%
% Write a plane, given filenames or spm_vol structures and 4-D data with
% the slice to write.
%
% SPM2/5 compatible, and creates image if necessary from file name and
% example V structure.
%
% ::
%
%    P = char({'one.img', 'two.img', 'three.img'})
%    dat = randn(64, 64, 3);
%    Vout = scn_write_plane(P, dat, wh_slice, V)
%    spm_image('init', Vout(1).fname);
%
%
% Notes on SPM5 and why we need to do what we do the way we do it:
% Write data...in SPM5, by accessing file_array object in
% V.private directly (with spm_write_plane. spm_write_vol
% first creates NIFTI obj/file array and then uses spm_write_plane.)
%
% The file_array object points to data in
% the actual file, so when values are assigned to the array
% object, they are written directly in the file.
% These values depend on the offset and slope values (spm
% scaling factors!) in V.private.dat as well, so care must be
% taken to assign data to a file_array object with the correct
% scaling factors.  This is why it is better to load the
% structure one wants to write to with spm_vol first, so that
% the name in V.private.dat.fname will be correct, and
% V.private.dat.scl_slope and ...inter will be correct as well.
% Vout = spm_vol(V(i));  % loads correct .private info from V.fname
%
% in SPM5, this simply assigns data in Vout.private.dat
% in SPM2, it does something different, but should be
% compatible, since spm_vol was used above...
% ::
%
%    spm_write_plane(Vout, fsl(:, :, i), slicei);
%    if isfield(V, 'dt'), Vout.dt = V.dt; end      % SPM5 only
%    if isfield(V, 'n'), Vout.n = V.n;  end        % SPM5 only
%        Vout = spm_create_vol(Vout);
%    else
%        Vout = spm_vol(V(i).fname);
%    end
%    spm_write_plane(Vout, fsl(:, :, i), slicei);
%
% ..
%    Tor Wager, March 2008
% ..

% ..
%    Make sure all images exist, and create if not
%    Return V structure of spm_vol mapped image structures
% ..

if ~isstruct(filenames_or_V)
    % Assume we have a string matrix of file names
    % Make sure they exist...
    % create volumes, if necessary

    % Make sure char array, from cell...
    if iscell(filenames_or_V), filenames_or_V = char(filenames_or_V{:}); end

    for i = 1:size(filenames_or_V, 1)
        currentfilename = deblank(filenames_or_V(i, :));

        if ~(exist(currentfilename, 'file'))

            fprintf('Creating Image: %3.0f ', i);

            if isempty(varargin) && ~exist('exampleV', 'var')
                error('Must pass in example V struct!');
            elseif ~exist('exampleV', 'var')
                exampleV = varargin{1};
            end

            % Create empty float image with no spm scaling (pinfo intercept is 0,
            % slope is 1)
            % --------------------------------------------------------------
            %       V.fname - the filename of the image.
            %       V.dim   - the x, y and z dimensions of the volume
            %       V.dt    - A 1x2 array.  First element is datatype (see spm_type).
            %                 The second is 1 or 0 depending on the
            %                 endian-ness.
            %                   * 1 is big-endian, 0 is little-endian
            %       V.mat   - a 4x4 affine transformation matrix mapping from
            %                 voxel coordinates to real world coordinates.
            %       V.pinfo - plane info for each plane of the volume.
            %              V.pinfo(1,:) - scale for each plane
            %              V.pinfo(2,:) - offset for each plane
            %                 The true voxel intensities of the jth image are given
            %                 by: val*V.pinfo(1,j) + V.pinfo(2,j)
            %              V.pinfo(3,:) - offset into image (in bytes).
            %                 If the size of pinfo is 3x1, then the volume is assumed
            %                 to be contiguous and each plane has the same scalefactor
            %                 and offset.

            Vout = struct('fname', currentfilename, 'mat', exampleV.mat, 'dim', exampleV.dim, 'pinfo', [1 0 0]');
            %if isfield(exampleV, 'dt'), Vout.dt = exampleV.dt; end      % SPM5 only; data type
            if isfield(exampleV, 'n'), Vout.n = exampleV.n;  end        % SPM5 only

            % Make sure datatype is float / float32
            switch(spm('Ver'))
                case 'SPM2'
                    Vout.dim(4) = spm_type('float');

                case 'SPM5'
                    Vout.dt = [spm_type('float32') 0];
                    if isfield(exampleV, 'dt'), Vout.dt(2) = exampleV.dt(2); end      % use endian-ness of example.

                case 'SPM8'
                    Vout.dt = [spm_type('float32') 0];
                    if isfield(exampleV, 'dt'), Vout.dt(2) = exampleV.dt(2); end      % use endian-ness of example.

                otherwise
                    error('Unknown SPM version.');
            end

            try
                V(i) = spm_create_vol(Vout);
                create_empty(V(i));

            catch
                try
                    % You cannot create only some volumes in a list, if
                    % others exist...work-around...
                    Vtmp = spm_create_vol(Vout);
                    Vtmp = spm_vol(Vtmp);
                    create_empty(Vtmp);
                    V(i) = spm_vol(Vtmp.fname);
                catch

                    disp('Cannot create vol.  This could occur for several reasons:')
                    disp('1) you may not have write permission.')
                    disp('2) Other reasons...?');
                    error('Quitting.');
                end
            end

        else
            % it exists; get it
            V(i) = spm_vol(currentfilename);
        end
    end

else
    % fprintf('Checking spm_vol structures ');
    % Make sure spm_vol structure (SPM5 only)
    % alter data type if needed
    switch(lower(spm('Ver')))
        case 'spm2'
            V = filenames_or_V; % Skip because hard to check.
        case {'spm5', 'spm8', 'spm12'}
            V = spm_vol(filenames_or_V);
%         otherwise
%             error('Fix or verify image writing for SPMxx. Not implemented yet.');
        otherwise
            error('Unknown SPM version.');
    end
end

str = sprintf('Writing slice data: %04d', 0);
fprintf('%s', str)

% V is a structure array, one structure per filename
% For each file to write to...write a plane

for i = 1:length(V)
    if mod(i, 10) == 0, fprintf('\b\b\b\b%04d', i); end

    spm_write_plane(V(i), dat(:, :, i), wh_slice);
end

erase_string(str)
%fprintf('\n');

end

function create_empty(V)
% Write empty data
switch(spm('Ver'))
    case 'SPM2'
        spm_write_vol(V, zeros(V.dim(1:3)));
    case {'SPM5', 'SPM8'}
        % assign directly
        for i = 1:length(V), V(i).private.dat(:, :, :) = 0; end
    otherwise
        error('Unknown SPM version.')
end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function dcData=dc_wavelet(dcfile)
[m,n]=size(dcfile);   % 4096x80
nw=32*32;
nbcol=size(colormap(gray),1);

for j=1:n
    X=double(reshape(dcfile(:,j),64,64));
    [cA,cH,cV,cD]=dwt2(X,'haar');
    cod_cH1=wcodemat(cH,nbcol);
    cod_cV1=wcodemat(cV,nbcol);
    cod_edge=cod_cH1+cod_cV1;
    dcData(:,j)=reshape(cod_edge,nw,1);
end


text:
function sys = probability(c)
% PROBABILITY Create basis for chance constraint
%
% EXAMPLE:
% The following example computes the largest value t such that the
% probability that a zero mean unit variance of a Gaussian variable is
% larger than t, is larger than 0.9
%
%  w = sdpvar(1,1);
%  t = sdpvar(1);
%  Model = [probability(a >= t) >= 0.9,uncertain(a,'normal',0,eye(1))];
%  solvesdp(derandomize(Model),-t)

if isa(c,'constraint') | isa(c,'lmi')
    if ~is(c,'elementwise')
        error('Probability constraints only applicable to single elementwise constraints');
    else
        sys.Constraint = c;
        sys = class(sys,'probability');
    end
else
    error('Probability constraints only applicable to single elementwise constraints');
end


text:
function p = predictOneVsAll(all_theta, X)
%PREDICT Predict the label for a trained one-vs-all classifier. The labels
%are in the range 1..K, where K = size(all_theta, 1).
%  p = PREDICTONEVSALL(all_theta, X) will return a vector of predictions
%  for each example in the matrix X. Note that X contains the examples in
%  rows. all_theta is a matrix where the i-th row is a trained logistic
%  regression theta vector for the i-th class. You should set p to a vector
%  of values from 1..K (e.g., p = [1; 3; 1; 2] predicts classes 1, 3, 1, 2
%  for 4 examples)

m = size(X, 1);
num_labels = size(all_theta, 1);

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

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters (one-vs-all).
%               You should set p to a vector of predictions (from 1 to
%               num_labels).
%
% Hint: This code can be done all vectorized using the max function.
%       In particular, the max function can also return the index of the
%       max element, for more information see 'help max'. If your examples
%       are in rows, then, you can use max(A, [], 2) to obtain the max
%       for each row.
%
% temp = (X*all_theta');
% p(:) = max(temp,[],2);

h_theta = sigmoid(X * all_theta');
p = max(h_theta, [], 2);

% temp is of dimension 5000x10

% Dimention of X: 5000x401
% Dimension of all_theta: 10x401
% Dimension of X*all_theta': 5000x10

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

end

author: anirudhjayaraman

repo: Machine-Learning