text:
function [bndinfo, pbim, gconf, bndinfo_all] =  processIm2Occlusion(im, varargin)

%% Set parameters

%% Read image

if max(size(im))>640
  fprintf('Warning, this image is pretty big...\n');
  %im = imresize(im, 640/max(size(im)), 'bilinear');
end

%% Get occlusion info
[bndinfo, pbim, gconf, bndinfo_all] = im2boundariesTopLevel(im);
gconf = single(gconf);
pbim = single(pbim);

%save(outname, 'bndinfo', 'pbim', 'gconf', 'bndinfo_all');

author: Cloud-CV

repo: object-proposals


text:
function calpak_test018 ( )

%*****************************************************************************80
%
%% CALPAK_TEST018 tests JED_TO_YMDF_ENGLISH and YMDF_TO_JED_ENGLISH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST018\n' );
  fprintf ( 1, '  For the English calendar,\n' );
  fprintf ( 1, '  JED_TO_YMDF_ENGLISH: JED -> YMDF.\n' );
  fprintf ( 1, '  YMDF_TO_JED_ENGLISH: YMDF -> JED.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  JED (in)    YMDF               JED (out)\n' );
  fprintf ( 1, '\n' );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed1 = jed_test ( i );

    if ( jed1 < 0.0 )
      break
    end

    [ y2, m2, d2, f2 ] = jed_to_ymdf_english ( jed1 );

    s2 = ymdf_to_s_english ( y2, m2, d2, f2 );

    jed3 = ymdf_to_jed_english ( y2, m2, d2, f2 );

    fprintf ( 1, '  %11.2f  %20s  %11.2f\n', jed1, s2, jed3 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function spm_defaults
% Set the defaults which are used by SPM
%__________________________________________________________________________
%
% If you want to customise some defaults for your installation, do not
% modify this file directly, but create a file named spm_my_defaults.m
% instead, accessible from MATLAB search path; e.g., it can be saved in
% MATLAB Startup Folder: userhome/Documents/MATLAB.
%
% Example: create the following file to change the image file extension:
% ----------- file /home/karl/Documents/MATLAB/spm_my_defaults.m -----------
% global defaults
% defaults.images.format = 'img';
%--------------------------------------------------------------------------
%
% spm_defaults should not be called directly in any script or function
% (apart from SPM internals).
% To load the defaults, use spm('Defaults',modality).
% To get/set the defaults, use spm_get_defaults.
%
%                 ** This file should not be edited **
%__________________________________________________________________________
% Copyright (C) 1994-2014 Wellcome Trust Centre for Neuroimaging

% SPM
% $Id: spm_defaults.m 6255 2014-11-05 09:55:02Z christophe $

global defaults

% Command Line Usage default
%==========================================================================
defaults.cmdline  = 0;

% User Interface defaults
%==========================================================================
defaults.ui.monitor = NaN;
defaults.ui.colour  = [0.58 0.77 0.57];
defaults.ui.fs      = 14;  % unused
defaults.ui.print   = 'ps';
defaults.renderer   = 'zbuffer';

% File format specific
%==========================================================================
% Note that defaults.analyze.flip is no longer used.  Specifying the
% left/right handedness of the voxel indices is now done entirely by
% spm_flip_analyze_images.m
defaults.images.format  = 'nii'; % options: 'img', 'nii'
defaults.mat.format     = '-v6'; % options: '-mat', '-v6', '-v7.0', '-v7.3'

% Toolboxes defaults
%==========================================================================
defaults.tbx.dir = { fullfile(spm('Dir'),'toolbox') };

% DICOM Import defaults
%==========================================================================
defaults.dicom.root = 'flat'; % Folder hierarchy

% fMRI design defaults
%==========================================================================
defaults.stats.fmri.t   = 16;
defaults.stats.fmri.t0  = 8;
defaults.stats.fmri.hpf = 128;
defaults.stats.fmri.cvi = 'AR(1)';

% Mask defaults
%==========================================================================
defaults.mask.thresh    = 0.8;

% Stats defaults
%==========================================================================
defaults.stats.maxmem      = 2^26;
defaults.stats.maxres      = 64;
defaults.stats.resmem      = false;
defaults.stats.fmri.ufp    = 0.001;  % Upper tail F-probability
defaults.stats.pet.ufp     = 0.05;
defaults.stats.eeg.ufp     = 0.05;
defaults.stats.topoFDR     = 1;
defaults.stats.rft.nonstat = 0;
defaults.stats.results.volume.distmin =  8;
defaults.stats.results.volume.nbmax   =  3;
defaults.stats.results.svc.distmin    =  4;
defaults.stats.results.svc.nbmax      = 16;
defaults.stats.results.mipmat         = {fullfile(spm('dir'),'MIP.mat')};

% Filename prefix defaults
%==========================================================================
defaults.slicetiming.prefix     = 'a';
defaults.realign.write.prefix   = 'r';
defaults.coreg.write.prefix     = 'r';
defaults.unwarp.write.prefix    = 'u';
defaults.normalise.write.prefix = 'w';
defaults.smooth.prefix          = 's';
defaults.imcalc.prefix          = 'i';

% Realignment defaults
%==========================================================================
defaults.realign.estimate.quality = 0.9;
defaults.realign.estimate.interp  = 2;
defaults.realign.estimate.wrap    = [0 0 0];
defaults.realign.estimate.sep     = 4;
defaults.realign.estimate.fwhm    = 5;
defaults.realign.estimate.rtm     = 1;
defaults.realign.write.mask       = 1;
defaults.realign.write.interp     = 4;
defaults.realign.write.wrap       = [0 0 0];
defaults.realign.write.which      = [2 1];

% Unwarp defaults
%==========================================================================
defaults.unwarp.estimate.rtm      = 0;
defaults.unwarp.estimate.fwhm     = 4;
defaults.unwarp.estimate.basfcn   = [12 12];
defaults.unwarp.estimate.regorder = 1;
defaults.unwarp.estimate.regwgt   = 1e5;
defaults.unwarp.estimate.foe      = [4 5];
defaults.unwarp.estimate.soe      = [];
defaults.unwarp.estimate.rem      = 1;
defaults.unwarp.estimate.jm       = 0;
defaults.unwarp.estimate.noi      = 5;
defaults.unwarp.estimate.expround = 'Average';
%
% Unwarp uses defaults.realign.write defaults for writing.
%

% Coregistration defaults
%==========================================================================
defaults.coreg.estimate.cost_fun = 'nmi';
defaults.coreg.estimate.sep      = [4 2];
defaults.coreg.estimate.tol      = [0.02 0.02 0.02 0.001 0.001 0.001 0.01 0.01 0.01 0.001 0.001 0.001];
defaults.coreg.estimate.fwhm     = [7 7];
defaults.coreg.write.interp      = 4;
defaults.coreg.write.wrap        = [0 0 0];
defaults.coreg.write.mask        = 0;

% Spatial Normalisation defaults
%==========================================================================
defaults.normalise.write.preserve   = 0;
defaults.normalise.write.bb         = [[-78 -112 -70];[78 76 85]];
defaults.normalise.write.vox        = [2 2 2];
defaults.normalise.write.interp     = 4;

% Old Spatial Normalisation defaults
%==========================================================================
defaults.old.normalise.estimate.smosrc  = 8;
defaults.old.normalise.estimate.smoref  = 0;
defaults.old.normalise.estimate.regtype = 'mni';
defaults.old.normalise.estimate.weight  = '';
defaults.old.normalise.estimate.cutoff  = 25;
defaults.old.normalise.estimate.nits    = 16;
defaults.old.normalise.estimate.reg     = 1;
defaults.old.normalise.write.preserve   = 0;
defaults.old.normalise.write.bb         = [[-78 -112 -70];[78 76 85]];
defaults.old.normalise.write.vox        = [2 2 2];
defaults.old.normalise.write.interp     = 1;
defaults.old.normalise.write.wrap       = [0 0 0];
defaults.old.normalise.write.prefix     = 'w';

% VBM Preprocessing defaults
%==========================================================================
defaults.old.preproc.tpm     = cellstr(char(...
    fullfile(spm('Dir'),'toolbox','OldSeg','grey.nii'),...
    fullfile(spm('Dir'),'toolbox','OldSeg','white.nii'),...
    fullfile(spm('Dir'),'toolbox','OldSeg','csf.nii'))); % Prior probability maps
defaults.old.preproc.ngaus          = [2 2 2 4]'; % Gaussians per class
defaults.old.preproc.warpreg        = 1;          % Warping Regularisation
defaults.old.preproc.warpco         = 25;         % Warp Frequency Cutoff
defaults.old.preproc.biasreg        = 0.0001;     % Bias regularisation
defaults.old.preproc.biasfwhm       = 60;         % Bias FWHM
defaults.old.preproc.regtype        = 'mni';      % Affine Regularisation
defaults.old.preproc.fudge          = 5;          % Fudge factor, can not be set in batch GUI
defaults.old.preproc.samp           = 3;          % Sampling distance
defaults.old.preproc.output.GM      = [0 0 1];
defaults.old.preproc.output.WM      = [0 0 1];
defaults.old.preproc.output.CSF     = [0 0 0];
defaults.old.preproc.output.biascor = 1;
defaults.old.preproc.output.cleanup = 0;

% Smooth defaults
%==========================================================================
defaults.smooth.fwhm = [8 8 8];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-Prevent users from making direct calls to spm_defaults
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent runOnce
try
    if ~isdeployed && isempty(runOnce)
        d = dbstack;
        if isempty(intersect({'spm','spm_get_defaults'},{d.name}))
            fprintf(['Direct calls to spm_defauts are deprecated.\n' ...
                'Please use spm(''Defaults'',modality) ' ...
                'or spm_get_defaults instead.\n']);
            runOnce = 1;
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-Execute user-specified defaults files
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def = defaults;
user_defaults = 'spm_my_defaults.m';
if exist(user_defaults,'file')
    if isdeployed && exist(fullfile(spm('Dir'),user_defaults),'file')
        user_defaults_file = cellstr(fullfile(spm('Dir'),user_defaults));
    else
        user_defaults_file = which(user_defaults,'-ALL');
    end
    for i=1:numel(user_defaults_file)
        try
            spm('run', user_defaults_file{i});
        catch
            lr = lasterror;
            warning(lr.message);
        end
    end
    if spm_check_version('matlab','8.0') >= 0, my_isequaln = @isequaln;
    else my_isequaln = @isequalwithequalnans; end
    if ~my_isequaln(def,defaults)
        fprintf('Defaults settings have been modified by file(s):\n');
        for i=1:numel(user_defaults_file)
            fprintf('  %s\n',user_defaults_file{i});
        end
        fn0 = fieldnames(def);
        mf = fn0(~cellfun(@(x) my_isequaln(def.(x),defaults.(x)),fn0));
        if ~isempty(mf)
            fprintf('Modified fields: ');
            for i=1:numel(mf)
                fprintf('%s ',mf{i});
            end
            fprintf('\n');
        end
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function X = QLiftRec2(C, S, filtername)
%-----------------------------------------------------------------------------
% QLiftRec2
% Multilevel 2-D reconstruction by inverting the lifting scheme and using
% quincunx grids
%
% Syntax: X = QLiftRec2(C, S, filtername)
%
% QLiftRec2 performs the reconstruction of a two-dimensional signal (matrix,
% image) X by inverting the lifting scheme using prediction and update filters
% that are indicated by filtername.
% The reconstruction involves the vector of coefficients C and the bookkeeping
% matrix S. The structure and dimensions of C and S are supposed to be
% consistent with how they are produced by QLiftDec2.
% QLiftRec2 is the inverse function of QLiftDec2.
% filtername is supposed to be the same string as chosen for QLiftDec2.
%
% Output X is a two-dimensional signal (image).
%
% Calls for: QLiftRec2Nevill, QLiftRec2MaxMin, QLiftRec2MinMin, QLiftRec2MaxMax
%
% See also: QLiftDec2
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: February 17, 2003.
% (c) 1999-2003 Stichting CWI, Amsterdam
%-----------------------------------------------------------------------------
% Note: future argument list of this function might be extended with an
% argument that points to another level than 1, nargin could be checked
% for this.
if isempty(C)
  error(' QLiftRec2 - empty decomposition ');
else
  if isempty(S)
    error(' QLiftRec2 - empty bookkeeping ');
  else
    if strncmpi(filtername,'Neville',7)
       X = QLiftRec2Nevill(C,S,filtername);
    elseif strncmpi(filtername,'MaxMin',6)
       X = QLiftRec2MaxMin(C,S);
    elseif strncmpi(filtername,'MinMin',6)
       X = QLiftRec2MinMin(C,S);
    elseif strncmpi(filtername,'MaxMax',6)
       X = QLiftRec2MaxMax(C,S);
%   elseif strncmpi(filtername,'some',4)
%      stencilP = something;
%      centerP =  something;
%      stencilU = something;
%      centerU =  something;
%      X = QLiftRec2Custom(C, S, stencilP, centerP, stencilU, centerU);
    else
       error([' QLiftRec2 - unknown filter ' filtername]);
    end
  end
end
%
%-----------------------------------------------------------------------------

author: Sable

repo: mcbench-benchmarks


text:
function [sos,g] = driving_function_imp_nfchoa_ls(N,R,r,conf)
%DRIVING_FUNCTION_IMP_NFCHOA_LS second-order section representation for a
%line source in NFC-HOA
%
%   Usage: sos = driving_function_imp_nfchoa_ls(N,R,r,conf)
%
%   Input parameters:
%       N       - order of spherical Hankel function
%       R       - radius of secondary source array / m
%       r       - distance of line source from array center / m
%       conf    - configuration struct (see SFS_config)
%
%   Output parameters:
%       sos     - second-order section representation
%       g       - scalar gain factor
%
%   See also: sound_field_imp, sound_field_imp_nfchoa,
%       driving_function_imp_nfchoa

%*****************************************************************************
% 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 of input  parameters ==================================
nargmin = 4;
nargmax = 4;
narginchk(nargmin,nargmax);
isargpositivescalar(N,R,r);
isargstruct(conf);

%% ===== Configuration ==================================================
% Speed of sound
dimension = conf.dimension;
driving_functions = conf.driving_functions;

%% ===== Computation =====================================================
% Find spherical Hankel function zeros
z = sphbesselh_zeros(N);

% Get the delay and weighting factors
if strcmp('2D',dimension)

    % === 2-Dimensional ==================================================

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        to_be_implemented;
    otherwise
        error(['%s: %s, this type of driving function is not implemented', ...
            'for a 2D line source.'],upper(mfilename),driving_functions);
    end

elseif strcmp('2.5D',dimension)

    % === 2.5-Dimensional ================================================

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        to_be_implemented;
    otherwise
        error(['%s: %s, this type of driving function is not implemented', ...
            'for a 2.5D line source.'],upper(mfilename),driving_functions);
    end

elseif strcmp('3D',dimension)

    % === 3-Dimensional ==================================================

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        to_be_implemented;
    otherwise
        error(['%s: %s, this type of driving function is not implemented', ...
            'for a 3D line source.'],upper(mfilename),driving_functions);
    end

else
    error('%s: the dimension %s is unknown.',upper(mfilename),dimension);
end

author: sfstoolbox

repo: sfs-matlab


text:
function traj = dirTrans_Euler(config)

n = config.grid.nTrajPts;

% Create the initial guess:
guess.time = linspace(config.guess.time(1),config.guess.time(end),n);
guess.state = interp1(...
    config.guess.time', config.guess.state', guess.time')';
guess.control = interp1(...
    config.guess.time', config.guess.control', guess.time')';

% Create a list of all linear constraints, then add function handles:
[problem, pack] = buildConstraints(guess,config);
problem.objective = @(z)( costFunctionWrapper(z,pack) );
problem.nonlcon = @(z)( nonLinCon(z,pack, config) );

% Solve using fmincon:
[zSoln,fSoln,exitFlag] = fmincon(problem);

% Post-processing:
[t,x,u] = unPackDecVar(zSoln,pack);

traj.time = linspace(t(1),t(2),n);
traj.state = x;
traj.control = u;
traj.objVal = fSoln;
traj.exitFlag = exitFlag;
traj.interp.state = @(tt)( interp1(traj.time',traj.state',tt')' );
traj.interp.control = @(tt)( interp1(traj.time',traj.control',tt')' );

end

function [C, Ceq] = nonLinCon(z,pack,config)
%
% This function enforces the dynamics of the cart-pole system
%

% Based on the Trapazoid Method for discretization, as defined in Bett's book,
% chapter 4.

n = pack.nState(2);
[t,x,u] = unPackDecVar(z,pack);
dt = (t(2)-t(1))/(n-1);

% Evaluate the dynamics at each collocation point
dx = cartPoleDynamics(x,[u; zeros(size(u))],config.dyn);

% Trapazoid rule:
idxLow = 1:(n-1);
idxUpp = 2:n;
intStateTrap = 0.5*dt*(dx(:,idxLow) + dx(:,idxUpp));
intStateCol = x(:,idxUpp)-x(:,idxLow);

% Defect constraint:
defect = intStateTrap - intStateCol;
%
% % user-defined boundary constraints:
% [bndIneq, bndEq] = boundaryConstraint(t,x(:,1),x(:,end),config.userData);

C = [];%bndIneq;
Ceq = [reshape(defect,numel(defect),1)];% bndEq];

end

function cost = costFunctionWrapper(z,pack)

[t,x,u] = unPackDecVar(z,pack);

% Trapazoid rule to integrate cost function:
tt = linspace(t(1),t(end),pack.nState(2));  %Time vector
dc = costFunction(tt,x,u);
nTime = size(x,2);
dt = (t(2)-t(1))/(nTime-1);
w = ones(nTime,1); w([1,end]) = 0.5;
cost = dt*dc*w;

end

author: MatthewPeterKelly

repo: dscTutorials


text:
clear Models MM GA

% User-Specified Parameters
% ===============================================================

% Conditions and stimuli
% ---------------------------------------------------------------
GA.conditions = [1 2 3 4 5];
	% a vector of integers 1:# of event-related conditions (e.g., [1 2 3 4])
	% to add "rest" intervals, add an extra condition and set its contrast weight
	% to zero in all contrasts.
	% e.g., condition 3 in this design could be passive rest.

GA.freqConditions = [.20 .20 .20 .20 .20];
	% vector of frequencies of each trial type; should sum to 1

GA.scanLength = 360;
	% how long your run is, in seconds.
    % This and the ISI determine how many stimuli are in each design vector

GA.ISI = 2;
	% how long between stimulus presentations? (you can also include "rest" presentations)
	% also the time resolution of stimulus condition function (list of stimuli)
	% designs will be constructed in time units of this resolution

GA.TR = 2;
    % the TR (sampling resolution) of your experiment; time for volume acquisition

GA.epochdur = 2;        % build epochs instead of events; enter dur in sec (default is events, 0 sec)

% Hox optimization parameters
% A hox gene is a master gene.  These numbers control the stimulus
% parameters, which can be optimized within the GA rather than
% pre-specified.
% Hox elements in this GA go in the stimlist as the first several elements.
% hox sequence currently codes, in number order:
% ISI TR cuelen cuerest stimlen stimrest resplen resprest (in s)
% the last 6 parameters are for use with doepochs.
% They code epoch lengths for 6 periods within each trial.
% TR does not work right now
% use zeros to use default rather than variable hox parameters
% ---------------------------------------------------------------
GA.numhox = 0;
GA.hoxrange = [];
    % rows are each hox gene's allowable range, cols index hox elements
	% use this if you want to determine the ISI, etc. using the GA
	% also may require significant user input/program modification at this stage.

% Genetic algorithm parameters
% ---------------------------------------------------------------
nmodels = 1;
	% how many runs do you want to optimize?
	% GA runs one separate optimization for each model.

GA.cbalColinPowerWeights = [0 1 0 .3];	% 1 = cbal, 2 = eff, 3 = hrf shape, 4 = freq
	% first element: counterbalancing of stimuli
	% second: contrast detection efficiency
	% third: hrf shape estimation efficiency
	% maintenance of input frequencies for each trial type
	% fitness scores for each measure are multiplied by these values
	% before collapsing to a single, final fitness measure for each design.
	% does not have to sum to 1.

GA.numGenerations = 10;
	% how many iterations of the GA to run.
GA.sizeGenerations = 40;
	% how many designs to test  for each generation?  Population size.
GA.maxTime = 30;
	% max time to run in s, or Inf for infinite time
	% The GA stops when either numGenerations or maxTime is reached.

GA.alph = 2.1;
	% "selection pressure": higher is more extreme selection; always pick the best 50% of designs
	%1 is no selection pressure at all, or random selection of designs for recombination
	% selection pressure is like in evolution; refers to whether the best designs are selected to
	% continue on to the next generation.
	% an intermediate value is best, to keep the population heterogeneity high.

GA.plotFlag = 1;
	% plot results after each run.
	% Recommended to leave this off if nmodels > 1, and then plot your final results later.

% Filtering, counterbalancing, and design tolerance
% ---------------------------------------------------------------
% Empty brackets indicate the option is not to be used.

GA.lowerLimit = [];
GA.HPlength = [120];
	% high-pass filter length, in s, for analysis; [] for no HP filter
	% Used to filter out noise at lower frequencies than design during data analysis
	% The cutoff you should use depends on how much power in your design is below the cutoff -
	% you want all the power in your design to be above the cutoff.

GA.LPsmooth = [];
	% low-pass filter type for analysis; [] for no LP filter
	% choices are 'hrf' or []

GA.maxOrder = 1;
    % order of counterbalancing to use
    % 1 = each trial type follows each other with equal probabilities, adjusted for base frequencies
    % 2 = one-back and two-back counterbalancing.  3 = 1+2+3 back, etc.

% Hard constraints
% ---------------------------------------------------------------
% Here you can specify tolerances for designs
% If designs do not meet these criteria, they will receive -infinity fitness scores
% Empty brackets indicate the option is not to be used.

GA.NumStimthresh = [];      % maximum number of repeats of any one event type in a row
GA.maxCbalDevthresh = []; 	% maximum acceptable counterbalancing deviation (actual freq - expected freq)
GA.maxFreqDevthresh = .2;   % maximum acceptable deviation from input frequencies (GA.freqConditions)

% Contrast setup
% ---------------------------------------------------------------
% Here you specify contrasts across conditions, for use with the contrast estimation fitness measure
% Each contrasts should be a row in the matrix GA.contrasts
%
GA.contrasts = [1 0 0 0 0; ...
                0 1 0 0 0; ...
                0 0 1 0 0; ...
                0 0 0 1 0];
    % there should be one column per condition in your design, not including the intercept.
    %if trans2switch or trans2block = 1, double the number of columns.
    % when using epoch design, three elements per condition for epochs (3 epochs)!
GA.contrastweights = [1 1 1 1];
    % Weighting function for contrasts (rows of GA.contrasts)
    % Contrast efficiencies will be multiplied by these weights before computing overall design fitness
    % If no contrasts are specified, this vector can specify weights for predictors (columns)

% Autocorrelation and special options
% ---------------------------------------------------------------
AutocorrelationFileName = 'myscannerxc';
    % This is the name of a mat file without the .mat extension
    % In the mat file, there should be a variable called myscannerxc
    % which is a row vector containing the autocorrelation fucntion you wish to use.
    % This function is used as the intrinsic noise autocorrelation estimate.
    % Included with the GA scripts are these functions:
    % myscannerxc   (U of M 3T, 7 subjects)
    % hiautocorr    (1/f)
    % hiautocorr3   (another 1/f)
    % noautocorr    (identity; no autocorrelation)

GA.restlength = [];
    % if inserting probe or instruction periods, modeled with periodic rests in the simulation,
    % this is the length of the rest/probe/instr periods to insert, in units of the ISI
GA.restevery = [];
    % This is how many regular events you want to have (in ISIs) between rests.
    % Leave restlength and restevery blank to avoid using these options.
GA.trans2switch = 0;
    % 1 or 0.  This option works if the conditions field is 2 elements long.
    % This option creates new conditions 3 and 4 that occur whenever 1 is followed by 2 or 2 by 1
    % So this option creates trial history-dependent predictors.  1 and 2 are repeats, 3 and 4 are switches
    % Useful when studying habituation, switching between items, etc.
GA.trans2block = 0;
    % 1 or 0.  This option doubles the input length of conditions - so [1 2] is transformed to [1 2 3 4]
    % This option creates alternating ABAB blocks of your input event types (e.g., 1 2) with the new
    % event types (e.g., 3 4).  This is useful when simulating mixed block/event related designs.
    % Trial types 1-2, for example, can be two events within block A, and 3-4 can be the same (or different)
    % events within block B.  The alternation frequency (in ISIs) is specified by GA.restevery.
    % If you use this option, the contrasts field must contain twice as many columns as you have input conditions.
    % for example, if conditions = [1 2 3 4], contrasts = [1 1 1 1 -1 -1 -1 -1] specifies block A vs block B.
    % contrasts = [1 1 1 1 -1 -1 -1 -1; 1 1 -1 -1 1 1 -1 -1] specifies A vs B and events [1 2 5 6] - [3 4 7 8].
GA.dofirst = 0;
    % 1 or 0.  This creates a separate trial type for the first trial following rests in GA.restlength.
    % It only works if you're using trans2switch and trans2block together; unreliable under other circumstances.
GA.nonlinthreshold = [2];
    % If the value here is x, predictor heights are clipped (thresholded) at x times the unit HRF height.
    % If empty, no thresholding is performed.
    % This is useful for modeling nonlinearities in the BOLD signal, which is important if your design uses
    % short ISIs (generally, 2 s or less).  The clipping acts as a 'saturation' factor.
    % I like to use 2, because data in our lab suggests the the nonlinearities in the HRF can be
    % roughly approximated by a clipping function with this parameter value.

% ---------------------------------------------------------------
%
% * E N D   U S E R   I N P U T
%
% ---------------------------------------------------------------
eval(['load ' AutocorrelationFileName]);
GA.xc = myscannerxc;

% =============================================
% * vary by parameter - set this in this script
% * if running multiple models, nmodels > 1
% ---------------------------------------------

varyparam = 0;
fieldname = 'alph';		% or 'freqConditions(end), etc.
incrementby = 0.3;
incrementevery = 5;		% every n models
contrastsofinterest = 1;% contrasts or predictors, if no cons specified

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

if varyparam
	eval(['paramvalue = GA.' fieldname ';'])
	disp(' ');disp('* *********************************************************************')
	disp('*  You have selected to vary a parameter across models in ga_gui_script')
	disp(['*  GA.' fieldname ' starting at ' num2str(paramvalue)])
	disp('* *********************************************************************')
end

for nm = 1:nmodels

   eval(['diary model' num2str(nm) '.log'])
   disp(['Starting Model ' num2str(nm)])
    % freqConditions

   M = optimizeGA(GA);

   Models(:,nm) = M.stimlist;
   MM{nm} = M;
   save GAworkspace

   if varyparam
   	if isfield(M,'consebeta'),
		Fit(1,nm) = mean(M.consebeta(contrastsofinterest));
	else
		Fit(1,nm) = mean(M.sebeta(contrastsofinterest));
	end
   	eval(['Fit(2,nm) = GA.' fieldname ';'])
   	FC(nm,:) = GA.freqConditions;
   	ALPH(nm,:) = GA.alph;
	str = ['if mod(nm,' num2str(incrementevery) ') == 0,GA.' fieldname ' = GA.' fieldname '+' num2str(incrementby) ';,end']
   	eval(str);
	eval(['paramvalue = GA.' fieldname ';'])
	disp('* *********************************************************************')
	disp('*  Varying parameter option on.')
	disp(['*  GA.' fieldname ' is now ' num2str(paramvalue)])
	disp('* *********************************************************************')
   end

   eval(['save model' num2str(nm) '_' M.date ' M'])
   save GAworkspace
   diary off
   fprintf('\n')
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function m = grammar2simple(gm)
% m = grammar2simple(gm)
%
% Convert a simple grammar model (with some strong assumptions on
% its structure) into a simple format for use with the cascade code.
% This function works with models produced by voc-release5.
% This structure roughly matches the model format from voc-release3.
%
% gm  grammar model returned by project_model.m

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2009-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

%% copy scalars
m.sbin      = gm.sbin;
m.thresh    = gm.thresh;
m.maxsize   = gm.maxsize;
m.minsize   = gm.minsize;
m.interval  = gm.interval;
m.numblocks = gm.numblocks;
m.class     = gm.class;
m.numcomponents ...
            = length(gm.rules{gm.start});
m.pca_coeff = gm.pca_coeff;
m.year      = gm.year;
m.note      = gm.note;
m.features  = gm.features;
if isfield(gm, 'bboxpred')
  m.bboxpred  = gm.bboxpred;
end

% Copy root filters
%% model.rootfilter{i} =
%     size: [6 13]
%     w: [6x13x31 double]
%     blocklabel: 2

% assume the first numcomponents filters are the root filters
m.rootfilters = cell(1, m.numcomponents);
for i = 1:m.numcomponents
  bl = gm.filters(i).blocklabel;
  m.rootfilters{i}.size = gm.filters(i).size;
  if gm.filters(i).flip
    m.rootfilters{i}.w    = single(gm.blocks(bl).w_flipped);
    m.rootfilters{i}.wpca = single(gm.blocks(bl).w_pca_flipped);
  else
    m.rootfilters{i}.w    = single(gm.blocks(bl).w);
    m.rootfilters{i}.wpca = single(gm.blocks(bl).w_pca);
  end
  m.rootfilters{i}.blocklabel = gm.filters(i).blocklabel;
end

% Copy component offsets
%% model.offsets{i} =
%     w: -4.2323
%     blocklabel: 1

m.offsets = cell(1, m.numcomponents);
for i = 1:length(gm.rules{gm.start})
  m.offsets{i}   = gm.rules{gm.start}(i).offset;
  m.offsets{i}.w = model_get_block(gm, gm.rules{gm.start}(i).offset) ...
                   * gm.features.bias;
end

m.loc_w = cell(1, m.numcomponents);
for i = 1:length(gm.rules{gm.start})
  m.loc{i}   = gm.rules{gm.start}(i).loc;
  m.loc{i}.w = model_get_block(gm, gm.rules{gm.start}(i).loc);
end

% Set up component structure
%% model.components{i} =
% x   rootindex: 1
% x   offsetindex: 1
% x   parts: {[1x1 struct]  [1x1 struct]  [1x1 struct]  [1x1 struct]  [1x1 struct]  [1x1 struct]}

numparts = 0;
m.components = cell(1, m.numcomponents);
for i = 1:m.numcomponents
  m.components{i}.rootindex = i;
  m.components{i}.offsetindex = i;
  n = length(gm.rules{gm.start}(i).rhs)-1;
  m.components{i}.parts = cell(1, n);
  numparts = numparts + n;
end

% Copy part filters, deformation models, and set up
% components{i}.parts{j}
%% model.partfilters{i} =
%    w: [9x6x31 double]
%    blocklabel: 5
%% model.defs{i} =
%    anchor: [1 4]
%    w: [0.0148 -4.1908e-04 0.0106 -0.0034]
%    blocklabel: 6
%% model.components{i}.parts{j} =
%    partindex: 1
%    defindex: 1

m.partfilters = cell(1, numparts);
m.defs = cell(1, numparts);
n = 0;
for i = 1:m.numcomponents
  m.components{i}.parts = cell(1, numparts/m.numcomponents);
  % assume first symbol is the root filter
  for j = 2:length(gm.rules{gm.start}(i).rhs)
    defsym = gm.rules{gm.start}(i).rhs(j);
    partsym = gm.rules{defsym}.rhs;
    filterid = gm.symbols(partsym).filter;
    bl = gm.filters(filterid).blocklabel;
    n = n + 1;
    if gm.filters(filterid).flip
      m.partfilters{n}.w    = single(gm.blocks(bl).w_flipped);
      m.partfilters{n}.wpca = single(gm.blocks(bl).w_pca_flipped);
    else
      m.partfilters{n}.w    = single(gm.blocks(bl).w);
      m.partfilters{n}.wpca = single(gm.blocks(bl).w_pca);
    end
    m.partfilters{n}.blocklabel = gm.filters(filterid).blocklabel;
    m.defs{n}.w = model_get_block(gm, gm.rules{defsym}.def);
    m.defs{n}.blocklabel = gm.rules{defsym}.def.blocklabel;
    m.defs{n}.anchor = gm.rules{gm.start}(i).anchor{j}(1:2);
    m.components{i}.parts{j-1}.partindex = n;
    m.components{i}.parts{j-1}.defindex = n;
  end
end

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function out = isMisorientation(o)
% check whether o is a misorientation

out = isa(o.SS,'crystalSymmetry') && isa(o.CS,'crystalSymmetry');

author: mtex-toolbox

repo: mtex

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


text:
function opt = opt_proplistToStruct(varargin)
% OPT_PROPLISTTOSTRUCT - Make options struct from parameter/value list
%
%Synopsis:
%  OPT= opt_proplistToStruct(<'Param1', VALUE1, 'Param2', VALUE2, ...>)
%  OPT= opt_proplistToStruct(OPT_IN, <'Param1', VALUE1, 'Param2', VALUE2, ...>)
%
%Arguments:
%  OPT_IN:      STRUCT of optional properties
%  'Param1', VALUE1, ...: Property/value list
%
%Returns:
%  OPT:  STRUCT with (new) fields created from the property/value list
%
%Description:
%  Generates a struct OPT (or updates a given struct OPT_IN) with field
%  'Param1' set to value VALUE1, field 'Param2' set to value VALUE2, and
%  so forth (property/value list).
%
%See also opt_setDefaults opt_structToProplist

% 06-2012 Benjamin Blankertz

opt= [];
if nargin==0,
  return;
end

if isstruct(varargin{1}),
  % First input argument is already a structure: Start with that, write
  % the additional fields
  opt= varargin{1};
  iListOffset= 1;
elseif isempty(varargin{1}),
  % First input argument is empty, which can happen under some
  % circumstances. This is not an opt, but also not part of the list
  iListOffset = 1;
else
  % First argument is not a structure: Assume this is the start of the
  % parameter/value list
  iListOffset = 0;
end

nFields= (nargin-iListOffset)/2;
if nFields~=round(nFields),
  error('Invalid parameter/value list');
end

for ff= 1:nFields,
  fld = varargin{iListOffset+2*ff-1};
  if ~ischar(fld),
    error('Invalid parameter/value list');
  end
  opt.(fld)= varargin{iListOffset+2*ff};
end

author: PatternRecognition

repo: OpenBMI


text:
clc;

%% compile
% mex kdtree_build.cpp
% mex kdtree_k_nearest_neighbors.cpp
% disp('compiled.');

%% create data and execute query
rand('seed',1);
p = rand( 30, 2 ); % input data
q = [.5,.5]; % query data
tree = kdtree_build( p );
idxs = kdtree_k_nearest_neighbors(tree,q,10);

%% visualize
close all;
xlim( [0 1] );
ylim( [0 1] );
hold on; axis equal; axis off;
plot( p(:,1), p(:,2), '.b');
plot(q(1), q(2),'.r');
plot(p(idxs,1), p(idxs,2),'or');
legend('database', 'query', 'query results');

%% text lables
for i=1:numel(idxs)
    text(p(idxs(i),1), p(idxs(i),2),sprintf(' %d',i));
end

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function a = r8cvv_iset ( mn, a, m, roff, i, j, aij )

%*****************************************************************************80
%
%% R8CVV_ISET sets item J from row I in an R8CVV.
%
%  Discussion:
%
%    An R8CVV is a "vector of vectors" of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 November 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer MN, the size of the cell array.
%
%    Input, real A(MN), the cell array.
%
%    Input, integer M, the number of rows in the array.
%
%    Input, integer ROFF(M+1), the row offsets.
%
%    Input, integer I, the row of the item.
%    1 <= I <= M.
%
%    Input, integer J, the column of the item.
%    1 <= J <= NR(I).
%
%    Input, real AIJ, the new value of item A(I,J).
%
%    Output, real A(MN), the updated cell array.
%
  k = roff(i) + j;
  a(k) = aij;

  return
end


text:
function filter_spec = ext2filter_spec(exts)
%EXT2FILTER_SPEC Generates a filter specification from a list of file extensions.
% Usage:
%   filter_spec = ext2filter_spec(ext) % single extension
%   filter_spec = ext2filter_spec(exts)
%
% See also: uigetfile, uigetdir, uibrowse

if ischar(exts)
    exts = {exts};
end

if ~iscellstr(exts)
    error('Expected a string or cell array of strings.')
end

% Keep only characters afte the period
exts = regexprep(exts, '.*\.', '');

% Create filter specification string
filter_spec = ['*.' strjoin(exts, ';*.')];

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function []=repop_testcases(testType)
%
% This file contains lots of test-cases to test the performance of the repop
% files vs. the matlab built-ins.
%
% N.B. there appears to be a bug in MATLAB when comparing mixed
% complex/real + double/single values in a max/min
%
% Copyright 2006-     by Jason D.R. Farquhar (jdrf@zepler.org)
% Permission is granted for anyone to copy, use, or modify this
% software and accompanying documents for any uncommercial
% purposes, provided this copyright notice is retained, and note is
% made of any changes that have been made. This software and
% documents are distributed without any warranty, express or
% implied
if ( nargin<1 || isempty(testType) ) testType={'acc','timing'}; end;

if ( ~isempty(strmatch('acc',testType)) )
fprintf('-------------------  Accuracy tests -------------------\n');

X=complex(randn(10,100),randn(10,100));
Y=complex(randn(size(X)),randn(size(X)));
fprintf('\n****************\n Double Real X, Double Real Y\n******************\n')
accuracyTests(real(X),real(Y),'dRdR')
fprintf('\n****************\n Double Complex X, Double Real Y\n******************\n')
accuracyTests(X,real(Y),'dCdR')
fprintf('\n****************\n Double Real X, Double Complex Y\n******************\n')
accuracyTests(real(X),Y,'dRdC')
fprintf('\n****************\n Double Complex X, Double Complex Y\n******************\n')
accuracyTests(X,Y,'dCdC')

fprintf('\n****************\n Double Real X, Single Real Y\n******************\n')
accuracyTests(real((X)),real(single(Y)),'dRsR')
fprintf('\n****************\n Double Complex X, Single Real Y\n******************\n')
accuracyTests((X),real(single(Y)),'dCsR')
fprintf('\n****************\n Double Real X, Single Complex Y\n******************\n')
accuracyTests((real(X)),single(Y),'dRsC')
fprintf('\n****************\n Double Complex X, Single Complex Y\n******************\n')
accuracyTests((X),single(Y),'dCsC')

fprintf('\n****************\n Single Real X, Double Real Y\n******************\n')
accuracyTests(real(single(X)),real((Y)),'sRdR')
fprintf('\n****************\n Single Complex X, Double Real Y\n******************\n')
accuracyTests(single(X),real((Y)),'sCdR')
fprintf('\n****************\n Single Real X, Double  Complex Y\n******************\n')
accuracyTests(single(real(X)),(Y),'sRdC')
fprintf('\n****************\n Single Complex X,Double Complex Y\n******************\n')
accuracyTests(single(X),(Y),'sCdC')

fprintf('\n****************\n Single Real X, Single Real Y\n******************\n')
accuracyTests(real(single(X)),real(single(Y)),'sRsR')
fprintf('\n****************\n Single Complex X, Single Real Y\n******************\n')
accuracyTests(single(X),real(single(Y)),'sCsR')
fprintf('\n****************\n Single Real X, Single Complex Y\n******************\n')
accuracyTests(single(real(X)),single(Y),'sRsC')
fprintf('\n****************\n Single Complex X, Single Complex Y\n******************\n')
accuracyTests(single(X),single(Y),'sCsC')

fprintf('All tests passed\n');
end

if ( ~isempty(strmatch('timing',testType)) )
fprintf('-------------------  Timing tests -------------------\n');

X=complex(randn(100,1000),randn(100,1000));
Y=complex(randn(size(X)),randn(size(X)));
timingTests(real(X),real(Y),'[100x1000] RR');
timingTests(X,Y,'[100x1000] CC');

end

return;

function []=accuracyTests(X,Y,str)
% PLUS
unitTest([str ' Matx + col Vec'],X,'+',Y(:,1),repop(X,'+',Y(:,1)),X+repmat(Y(:,1),[1,size(X,2)]));
unitTest([str ' Matx + row Vec'],X,'+',Y(1,:),repop(X,'+',Y(1,:)),X+repmat(Y(1,:),[size(X,1),1]));
unitTest([str ' Matx + Matx(:,1:2)'],X,'+',Y(:,1:2),repop(X,'+',Y(:,1:2),'m'),X+repmat(Y(:,1:2),[1,size(X,2)/2]));
unitTest([str ' Matx + Matx'],X,'+',Y(:,:),repop(X,'+',Y(:,:)),X+repmat(Y(:,:),[1,1]));
% TIMES
unitTest([str ' Matx * col Vec'],X,'*',Y(:,1),repop(X,'*',Y(:,1)),X.*repmat(Y(:,1),[1,size(X,2)]));
unitTest([str ' Matx * row Vec'],X,'*',Y(1,:),repop(X,'*',Y(1,:)),X.*repmat(Y(1,:),[size(X,1),1]));
unitTest([str ' Matx * Matx(:,1:2)'],X,'*',Y(:,1:2),repop(X,'*',Y(:,1:2),'m'),X.*repmat(Y(:,1:2),[1,size(X,2)/2]));
unitTest([str ' Matx * Matx'],X,'*',Y(:,:),repop(X,'*',Y(:,:)),X.*repmat(Y(:,:),[1,1]));

% other operations
unitTest([str ' Matx - row vec'],X,'-',Y(:,1),repop(X,'-',Y(:,1)),X-repmat(Y(:,1),[1,size(X,2)]));
unitTest([str ' Matx / row vec'],X,'/',Y(:,1),repop(X,'/',Y(:,1)),X./repmat(Y(:,1),[1,size(X,2)]));
unitTest([str ' Matx \ row vec'],X,'\',Y(:,1),repop(X,'\',Y(:,1)),X.\repmat(Y(:,1),[1,size(X,2)]));
unitTest([str ' Matx ^ row vec'],X,'^',Y(:,1),repop(X,'^',Y(:,1)),X.^repmat(Y(:,1),[1,size(X,2)]),1e-5);
unitTest([str ' Matx == row vec'],X,'==',Y(:,1),repop(X,'==',Y(:,1)),X==repmat(Y(:,1),[1,size(X,2)]));
unitTest([str ' Matx ~= row vec'],X,'~=',Y(:,1),repop(X,'~=',Y(:,1)),X~=repmat(Y(:,1),[1,size(X,2)]));
% N.B. repop tests with complex inputs use the magnitude of the input!
tX=X; tY=Y; if( ~isreal(X) | ~isreal(Y) ) tX=abs(X); tY=abs(Y); end;
unitTest([str ' Matx < row vec'],X,'<',Y(:,1),repop(X,'<',Y(:,1)),tX<repmat(tY(:,1),[1,size(X,2)]));
unitTest([str ' Matx > row vec'],X,'>',Y(:,1),repop(X,'>',Y(:,1)),tX>repmat(tY(:,1),[1,size(X,2)]));
unitTest([str ' Matx <= row vec'],X,'<=',Y(:,1),repop(X,'<=',Y(:,1)),tX<=repmat(tY(:,1),[1,size(X,2)]));
unitTest([str ' Matx >= row vec'],X,'>=',Y(:,1),repop(X,'>=',Y(:,1)),tX>=repmat(tY(:,1),[1,size(X,2)]));
%unitTest([str ' min Matx, row vec'],repop('min',X,Y(:,1)),min(X,repmat(Y(:,1),[1, size(X,2)])));
%unitTest([str ' max Matx, row vec'],repop('max',X,Y(:,1)),max(X,repmat(Y(:,1),[1,size(X,2)])));
%return;
%function []=inplaceaccuracyTests(X,Y,str)

return;

% Inplace operations test
% PLUS
% N.B. need the Z(1)=Z(1); to force to make a "deep" copy, i.e. not just
% equal pointers
Z=X;unitTest([str ' Matx + col Vec (inplace)'],Z,'+',Y(:,1),repop(Z,'+',Y(:,1),'i'),X+repmat(Y(:,1),[1,size(X,2)]));
Z=X;unitTest([str ' Matx + row Vec (inplace)'],Z,'+',Y(1,:),repop(Z,'+',Y(1,:),'i'),X+repmat(Y(1,:),[size(X,1),1]));
Z=X;unitTest([str ' Matx + Matx(:,1:2) (inplace)'],Z,'+',Y(:,1:2),repop(Z,'+',Y(:,1:2),'mi'),X+repmat(Y(:,1:2),[1,size(X,2)/2]));
Z=X;unitTest([str ' Matx + Matx (inplace)'],Z,'+',Y(:,:),repop(Z,'+',Y(:,:),'i'),X+repmat(Y(:,:),[1,1]));
% TIMES
Z=X;unitTest([str ' Matx * col Vec (inplace)'],Z,'*',Y(:,1),repop(Z,'*',Y(:,1),'i'),X.*repmat(Y(:,1),[1,size(X,2)]));
Z=X;unitTest([str ' Matx * row Vec (inplace)'],Z,'*',Y(1,:),repop(Z,'*',Y(1,:),'i'),X.*repmat(Y(1,:),[size(X,1),1]));
Z=X;unitTest([str ' Matx * Matx(:,1:2) (inplace)'],Z,'*',Y(:,1:2),repop(Z,'*',Y(:,1:2),'mi'),X.*repmat(Y(:,1:2),[1,size(X,2)/2]));
Z=X;unitTest([str ' Matx * Matx (inplace)'],Z,'*',Y(:,:),repop(Z,'*',Y(:,:),'i'),X.*repmat(Y(:,:),[1,1]));
% other operations
Z=X;unitTest([str ' Matx - row vec (inplace)'],Z,'-',Y(:,1),repop(Z,'-',Y(:,1),'i'),X-repmat(Y(:,1),[1,size(X,2)]));
Z=X;unitTest([str ' Matx \ row vec (inplace)'],Z,'\',Y(:,1),repop(Z,'\',Y(:,1),'i'),X.\repmat(Y(:,1),[1,size(X,2)]));
Z=X;unitTest([str ' Matx / row vec (inplace)'],Z,'/',Y(:,1),repop(Z,'/',Y(:,1),'i'),X./repmat(Y(:,1),[1,size(X,2)]));
Z=X;unitTest([str ' Matx ^ row vec (inplace)'],Z,'^',Y(:,1),repop(Z,'^',Y(:,1),'i'),X.^repmat(Y(:,1),[1,size(X,2)]),1e-5);
%unitTest([str ' min Matx, row vec (inplace)'],repop('min',X,Y(:,1),'i'),min(X,repmat(Y(:,1),[1, size(X,2)])));
%unitTest([str ' max Matx, row vec (inplace)'],repop('max',X,Y(:,1),'i'),max(X,repmat(Y(:,1),[1,size(X,2)])));

function []=timingTests(X,Y,str)
%PLUS
fprintf('%s Matx + Scalar\n',str);
tic; for i=1:1000; Z=repop(X,'+',10); end;
fprintf('%30s %gs\n','repop',toc/1000);
tic; Z=X;Z(1)=Z(1);for i=1:1000; Z=repop(Z,'+',10,'i'); end;
fprintf('%30s %gs\n','repop (inplace)',toc/1000);
tic; for i=1:1000; T=X+10;end;
fprintf('%30s %gs\n','MATLAB',toc/1000);

fprintf('%s Matx + col vec\n',str);
tic, for i=1:1000; Z=repop(X,'+',Y(:,1)); end;
fprintf('%30s %gs\n','repop',toc/1000); % = .05  / .01
tic, Z=X;Z(1)=Z(1); for i=1:1000; Z=repop(Z,'+',Y(:,1),'i'); end;
fprintf('%30s %gs\n','repop (inplace)',toc/1000); % = .05  / .01
tic, for i=1:1000; Z=X+repmat(Y(:,1),1,size(X,2));end;
fprintf('%30s %gs\n','MATLAB',toc/1000); % = .05  / .01

fprintf('%s Matx + row vec\n',str);
tic, for i=1:1000; Z=repop(X,'+',Y(1,:)); end;
fprintf('%30s %gs\n','repop',toc/1000); % = .05  / .01
tic, Z=X;Z(1)=Z(1); for i=1:1000; Z=repop(Z,'+',Y(1,:),'i'); end;
fprintf('%30s %gs\n','repop (inplace)',toc/1000); % = .05  / .01
tic, for i=1:1000; Z=X+repmat(Y(1,:),size(X,1),1);end;
fprintf('%30s %gs\n','MATLAB',toc/1000); % = .05  / .01

fprintf('%s Matx + Matx(:,1:2)\n',str);
tic; for i=1:1000; Z=repop(X,'+',Y(:,1:2),'m'); end;
fprintf('%30s %gs\n','repop',toc/1000); % = .05  / .01
tic, Z=X;Z(1)=Z(1); for i=1:1000; Z=repop(Z,'+',Y(:,1:2),'im'); end;
fprintf('%30s %gs\n','repop (inplace)',toc/1000); % = .05  / .01
tic; for i=1:1000; T=X+repmat(Y(:,1:2),[1,size(X,2)/2]);end;
fprintf('%30s %gs\n','MATLAB',toc/1000); % = .05  / .01

%TIMES
fprintf('%s Matx * col Vec\n',str);
tic; for i=1:1000; Z=repop(X,'*',Y(:,1)); end;
fprintf('%30s %gs\n','repop',toc/1000);
tic; Z=X;Z(1)=Z(1); for i=1:1000; Z=repop(Z,'*',Y(:,1),'i'); end;
fprintf('%30s %gs\n','repop (inplace)',toc/1000);
tic; for i=1:1000; Z=X.*repmat(Y(:,1),1,size(X,2));end;
fprintf('%30s %gs\n','MATLAB',toc/1000);
tic; for i=1:1000; Z=spdiags(Y(:,1),0,size(X,1),size(X,1))*X;end;
fprintf('%30s %gs\n','MATLAB (spdiags)',toc/1000);

fprintf('%s Matx * row Vec\n',str);
tic; for i=1:1000; Z=repop(X,'*',Y(1,:)); end;
fprintf('%30s %gs\n','repop',toc/1000);
tic; Z=X;Z(1)=Z(1); for i=1:1000; Z=repop(Z,'*',Y(1,:),'i'); end;
fprintf('%30s %gs\n','repop (inplace)',toc/1000);
tic; for i=1:1000; Z=X.*repmat(Y(:,1),1,size(X,2));end;
fprintf('%30s %gs\n','MATLAB',toc/1000);
tic; for i=1:1000; Z=X*spdiags(Y(1,:)',0,size(X,2),size(X,2));end;
fprintf('%30s %gs\n','MATLAB (spdiags)',toc/1000);

fprintf('%s Matx * Matx(:,1:2)\n',str);
tic; for i=1:1000; Z=repop(X,'*',Y); end;
fprintf('%30s %gs\n','repop',toc/1000);
tic; Z=X;Z(1)=Z(1); for i=1:1000; Z=repop(Z,'*',Y,'i'); end;
fprintf('%30s %gs\n','repop (inplace)',toc/1000);
tic; for i=1:1000; Z=X.*repmat(Y(:,1:2),[1,size(X,2)/2]);end;
fprintf('%30s %gs\n','MATLAB',toc/1000);

return;

% simple function to check the accuracy of a test and report the result
function [testRes,trueRes,diff]=unitTest(testStr,A,op,B,testRes,trueRes,tol)
global LOGFILE;
if ( ~isempty(LOGFILE) ) % write tests and result to disc
   writeMxInfo(LOGFILE,A);
   fprintf(LOGFILE,'%s\n',op);
   writeMxInfo(LOGFILE,B);
   fprintf(LOGFILE,'=\n');
   writeMxInfo(LOGFILE,trueRes);
   fprintf(LOGFILE,'\n');
end
if ( nargin < 7 )
   if ( isa(trueRes,'double') ) tol=1e-11;
   elseif ( isa(trueRes,'single') ) tol=1e-5;
   elseif ( isa(trueRes,'integer') )
      warning('Integer inputs!'); tol=1;
   elseif ( isa(trueRes,'logical') ) tol=0;
   end
end;
diff=abs(testRes-trueRes)./max(1,abs(testRes+trueRes));
fprintf('%45s = %0.3g ',testStr,max(diff(:)));
if ( max(diff(:)) > tol )
   if ( exist('mimage') )
      mimage(squeeze(testRes),squeeze(trueRes),squeeze(diff))
   end
   warning([testStr ': failed!']);
   fprintf('Type return to continue\n'); keyboard;
else
   fprintf('Passed \n');
end


text:
function timeString = BF_TheTime(tsec,formatLong)
% BF_TheTime    Converts a duration (seconds) to a human-interpretable string
%
% e.g., converts to minutes or hours or days as appropriate) output is something
% like '25.5 minutes' or '3.2 days' -- always displays to one decimal place.
%
% This code is useful for displaying user feedback on tic/toc statements.
%
%---INPUTS:
% tsec, the duration in seconds
% formatLong, (i) 0: display short units (like 's' instead of 'seconds')
%                    [default]
%             (ii) 1: display long units (like 'seconds' instead of 's')
%
%---OUTPUT:
% timeString, an interpretable text version of the input time.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

if nargin < 2 || isempty(formatLong)
    formatLong = false; % Set to 1 to use a longer format for the unit
end

if tsec < 1E-3
    if formatLong
        timeString = '< 1 milliseconds';
    else
        timeString = '< 1ms';
    end
elseif tsec < 1 % less than a second, display in integer number of milliseconds
    if formatLong
        timeString = sprintf('%.0f milliseconds',tsec*1000);
    else
        timeString = sprintf('%.0fms',tsec*1000);
    end
elseif tsec <= 60 % less than a minute, display in seconds
    if formatLong
        timeString = sprintf('%.1f seconds',tsec);
    else
    	timeString = sprintf('%.1fs',tsec);
    end
elseif tsec <= 60*60 % less than an hour, display in minutes
    if formatLong
        timeString = sprintf('%.1f minutes',tsec/60);
    else
        timeString = sprintf('%.1fmin',tsec/60);
    end
elseif tsec <= 60*24*60 % less than a day, display in hours
    if formatLong
        timeString = sprintf('%.1f hours',tsec/60/60);
    else
        timeString = sprintf('%.1fh',tsec/60/60);
    end
% elseif tsec<=60*24*7*60 % less than a week, display in days
else % display in days
	timeString = sprintf('%.1f days',tsec/60/60/24);
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
classdef CEC2008_F5 < PROBLEM
% <single> <real> <large/none> <expensive/none>
% Shifted Griewank's function

%------------------------------- Reference --------------------------------
% K. Tang, X. Yao, P. N. Suganthan, C. MacNish, Y.-P. Chen, C.-M. Chen, and
% Z. Yang, Benchmark functions for the CEC'2008 special session and
% competition on large scale global optimization, Nature Inspired
% Computation and Applications Laboratory, USTC, China, 2007.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;  % Optimal decision vector
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2008.mat'),'Data');
            obj.O = Data{5};
            obj.M = 1;
            if isempty(obj.D); obj.D = 100; end
            obj.D = min(obj.D,length(obj.O));
            obj.lower    = zeros(1,obj.D) - 600;
            obj.upper    = zeros(1,obj.D) + 600;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopObj = 1/4000*sum(Z.^2,2) - prod(cos(Z./sqrt(repmat(1:size(Z,2),size(Z,1),1))),2) + 1;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ft] = au2ft(au)
% Convert length from astronomical units to feet.
% Chad A. Greene 2012
ft = au*490806662372;


text:
function res = imMeanFilter(img, se, varargin)
%IMMEANFILTER Compute mean value in the neighboorhood of each pixel
%
%   RES = imMeanFilter(IMG, SE)
%   Computes the mean filter of image IMG, using structuring element SE.
%   The goal of this function is to provide the same interface as for
%   other image filters (imopen, imerode ...), and to allow the use of
%   mean filter with user-defined structuring element.
%   This function can be used for directional filtering.
%
%
%   RES = imMeanFilter(IMG, SE, PADOPT)
%   also specify padding option. PADOPT can be one of X (numeric value),
%   'symmetric', 'replicate', 'circular' (see imfilter for details).
%   Default is 'replicate'.
%
%   See Also:
%     imBoxFilter, imMedianFilter, imDirectionalFilter, imfilter
%

%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 16/02/2004.
%

%   HISTORY
%   2004-02-17 add support for 'strel' objects.
%   2004-02-20 add PADOPT option, and documentation.
%   2011-11-05 rename from imMean to imMeanFilter

% transform STREL object into single array
if strcmp(class(se), 'strel') %#ok<STISA>
    se = getnhood(se);
end

% get Padopt option
padopt = 'replicate';
if ~isempty(varargin)
    padopt = varargin{1};
end

% perform filtering
res = imfilter(img, se ./ sum(se(:)), padopt);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function a = jordan_inverse ( n, alpha )

%*****************************************************************************80
%
%% JORDAN_INVERSE returns the inverse of the JORDAN matrix.
%
%  Formula:
%
%    if ( I <= J )
%      A(I,J) =  -1 * (-1/ALPHA)^(J+1-I)
%    else
%      A(I,J) = 0
%
%  Example:
%
%    ALPHA = 2, N = 4
%
%    1/2 -1/4  1/8 -1/16
%    0    1/2 -1/4  1/8
%    0    0    1/2 -1/4
%    0    0    0    1/2
%
%  Properties:
%
%    A is upper triangular.
%
%    A is Toeplitz: constant along diagonals.
%
%    A is generally not symmetric: A' /= A.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    The inverse of A is the Jordan block matrix, whose diagonal
%    entries are ALPHA, whose first superdiagonal entries are 1,
%    with all other entries zero.
%
%    det ( A ) = (1/ALPHA)^N.
%
%    LAMBDA(1:N) = 1 / ALPHA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real ALPHA, the eigenvalue of the Jordan matrix.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  if ( alpha == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'JORDAN_INVERSE - Fatal error!\n' );
    fprintf ( 1, '  The input parameter ALPHA was 0.\n' );
    error ( 'JORDAN_INVERSE - Fatal error!' );
  end

  for i = 1 : n
    for j = 1 : n

      if ( i <= j )
        a(i,j) =  - ( - 1.0 / alpha )^(j+1-i);
      else
        a(i,j) = 0.0;
      end

    end
  end

  return
end


text:
function lambda = bab_eigenvalues ( n, alpha, beta )

%*****************************************************************************80
%
%% BAB_EIGENVALUES returns the eigenvalues of the BAB matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real ALPHA, BETA, the parameters.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  for i = 1 : n
    angle = i * pi / ( n + 1 );
    lambda(i,1) = alpha + 2.0 * beta * cos ( angle );
  end

  return
end


text:
% MKLJAC  Estimate the Jacobian of a function via Numerical Differences
%
% mklJac uses the Intel Math Kernel Library (MKL) djacobi function to
% estimate the Jacobian of a function using central differences.
%
%   jac = mklJac(fun,x) uses the supplied function handle fun and the
%   current state vector x to estimate the gradient of the function.
%
%   jac = mklJac(fun,x,nrow) specifies the number of rows in the vector
%   returned from fun, used to determine the return size of the Jacobian.
%   If the number of rows is not specified, then a dummy function call is
%   used in order to determine the number of rows.
%
%   jac = mklJac(fun,x,nrow,tol) specifies the tolerance of the numerical
%   difference algorithm of all variables. By default tol is 1e-6.
%
%   [jac,status] = mklJac(fun,x,nrow,tol) also returns 1 if the algorithm
%   was successful, or 0 if it failed.
%
%
%   Copyright (C) 2011 Jonathan Currie (I2C2)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function y=gauss_res2(p)
global grab;
global xpix;
global ypix;
global psf_w02;

yfit=p(3)*exp(-2*((xpix-p(1)).*(xpix-p(1))+(ypix-p(2)).*(ypix-p(2)))/psf_w02);
ydev=yfit-double(grab);
ydev2=ydev.*ydev;
y=sum(sum(ydev2));
%y=100*(p(2)-p(1)^2)^2+(1-p(1))^2;
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [errorcode, rows, columns] = rndcheck(nargs,nparms,arg1,arg2,arg3,arg4,arg5)
    %RNDCHECK error checks the argument list for the random number generators.

    %	B.A. Jones  1-22-93
    %	Copyright (c) 1993 by The MathWorks, Inc.
    %	$Revision: 1399 $  $Date: 2006-08-11 11:19:27 +0200 (Fr, 11 Aug 2006) $

    report_this_filefun(mfilename('fullpath'));

    sizeinfo = nargs - nparms;
    errorcode = 0;

    if nparms == 3
        [r1, c1] = size(arg1);
        [r2, c2] = size(arg2);
        [r3, c3] = size(arg3);
    end

    if nparms == 2
        [r1, c1] = size(arg1);
        [r2, c2] = size(arg2);
    end

    if sizeinfo == 0
        if nparms == 1
            [rows, columns] = size(arg1);
        end

        if nparms == 2
            scalararg1 = isscalar(arg1);
            scalararg2 = isscalar(arg2);
            if ~scalararg1  &&  ~scalararg2
                if r1 ~= r2 | c1 ~= c2
                    errorcode = 1;
                    return;
                end
            end
            if ~scalararg1
                [rows, columns] = size(arg1);
            elseif ~scalararg2
                [rows, columns] = size(arg2);
            else
                [rows, columns] = size(arg1);
            end
        end

        if nparms == 3
            scalararg1 = isscalar(arg1);
            scalararg2 = isscalar(arg2);
            scalararg3 = isscalar(arg3);

            if ~scalararg1  &&  ~scalararg2
                if r1 ~= r2 | c1 ~= c2
                    errorcode = 1;
                    return;
                end
            end

            if ~scalararg1  &&  ~scalararg3
                if r1 ~= r3 | c1 ~= c3
                    errorcode = 1;
                    return;
                end
            end

            if ~scalararg3  &&  ~scalararg2
                if r3 ~= r2 | c3 ~= c2
                    errorcode = 1;
                    return;
                end
            end
            if ~scalararg1
                [rows, columns] = size(arg1);
            elseif ~scalararg2
                [rows, columns] = size(arg2);
            else
                [rows, columns] = size(arg3);
            end
        end
    end

    if sizeinfo == 1
        scalararg1 = isscalar(arg1);
        if nparms == 1
            if prod(size(arg2)) ~= 2
                errorcode = 2;
                return;
            end
           if  ~scalararg1  && arg2 ~= size(arg1)
                errorcode = 3;
                return;
            end
            if (arg2(1) < 0 | arg2(2) < 0 | arg2(1) ~= round(arg2(1)) | arg2(2) ~= round(arg2(2)))
                errorcode = 4;
                return;
            end
            rows    = arg2(1);
            columns = arg2(2);
        end

        if nparms == 2
            if prod(size(arg3)) ~= 2
                errorcode = 2;
                return;
            end
            scalararg2 = isscalar(arg2);
            if ~scalararg1  &&  ~scalararg2
                if r1 ~= r2 | c1 ~= c2
                    errorcode = 1;
                    return;
                end
            end
            if (arg3(1) < 0 | arg3(2) < 0 | arg3(1) ~= round(arg3(1)) | arg3(2) ~= round(arg3(2)))
                errorcode = 4;
                return;
            end
            if ~scalararg1
                if any(arg3 ~= size(arg1))
                    errorcode = 3;
                    return;
                end
                [rows, columns] = size(arg1);
            elseif ~scalararg2
                if any(arg3 ~= size(arg2))
                    errorcode = 3;
                    return;
                end
                [rows, columns] = size(arg2);
            else
                rows    = arg3(1);
                columns = arg3(2);
            end
        end

        if nparms == 3
            if prod(size(arg4)) ~= 2
                errorcode = 2;
                return;
            end
            scalararg1 = isscalar(arg1);
            scalararg2 = isscalar(arg2);
            scalararg3 = isscalar(arg3);

            if (arg4(1) < 0 | arg4(2) < 0 | arg4(1) ~= round(arg4(1)) | arg4(2) ~= round(arg4(2)))
                errorcode = 4;
                return;
            end

            if ~scalararg1  &&  ~scalararg2
                if r1 ~= r2 | c1 ~= c2
                    errorcode = 1;
                    return;
                end
            end

            if ~scalararg1  &&  ~scalararg3
                if r1 ~= r3 | c1 ~= c3
                    errorcode = 1;
                    return;
                end
            end

            if ~scalararg3  &&  ~scalararg2
                if r3 ~= r2 | c3 ~= c2
                    errorcode = 1;
                    return;
                end
            end
            if ~scalararg1
                if any(arg4 ~= size(arg1))
                    errorcode = 3;
                    return;
                end
                [rows, columns] = size(arg1);
            elseif ~scalararg2
                if any(arg4 ~= size(arg2))
                    errorcode = 3;
                    return;
                end
                [rows, columns] = size(arg2);
            elseif ~scalararg3
                if any(arg4 ~= size(arg3))
                    errorcode = 3;
                    return;
                end
                [rows, columns] = size(arg3);
            else
                rows    = arg4(1);
                columns = arg4(2);
            end
        end
    end

    if sizeinfo == 2
        if nparms == 1
            scalararg1 = isscalar(arg1);
            if ~scalararg1
                [rows, columns] = size(arg1);
                if rows ~= arg2 | columns ~= arg3
                    errorcode = 3;
                    return;
                end
            end
            if (arg2 < 0 | arg3 < 0 | arg2 ~= round(arg2) | arg3 ~= round(arg3))
                errorcode = 4;
                return;
            end
            rows = arg2;
            columns = arg3;
        end

        if nparms == 2
            scalararg1 = isscalar(arg1);
            scalararg2 = isscalar(arg2);
            if ~scalararg1  &&  ~scalararg2
                if r1 ~= r2 | c1 ~= c2
                    errorcode = 1;
                    return;
                end
            end
            if ~scalararg1
                [rows, columns] = size(arg1);
                if rows ~= arg3 | columns ~= arg4
                    errorcode = 3;
                    return;
                end
            elseif ~scalararg2
                [rows, columns] = size(arg2);
                if rows ~= arg3 | columns ~= arg4
                    errorcode = 3;
                    return;
                end
            else
                if (arg3 < 0 | arg4 < 0 | arg3 ~= round(arg3) | arg4 ~= round(arg4))
                    errorcode = 4;
                    return;
                end
                rows = arg3;
                columns = arg4;
            end
        end

        if nparms == 3
            scalararg1 = isscalar(arg1);
            scalararg2 = isscalar(arg2);
            scalararg3 = isscalar(arg3);

            if ~scalararg1  &&  ~scalararg2
                if r1 ~= r2 | c1 ~= c2
                    errorcode = 1;
                    return;
                end
            end

            if ~scalararg1  &&  ~scalararg3
                if r1 ~= r3 | c1 ~= c3
                    errorcode = 1;
                    return;
                end
            end

            if ~scalararg3  &&  ~scalararg2
                if r3 ~= r2 | c3 ~= c2
                    errorcode = 1;
                    return;
                end
            end

            if ~scalararg1
                [rows, columns] = size(arg1);
                if rows ~= arg4 | columns ~= arg5
                    errorcode = 3;
                    return;
                end
            elseif ~scalararg2
                [rows, columns] = size(arg2);
                if rows ~= arg4 | columns ~= arg5
                    errorcode = 3;
                    return;
                end
            elseif ~scalararg3
                [rows, columns] = size(arg3);
                if rows ~= arg4 | columns ~= arg5
                    errorcode = 3;
                    return;
                end
            else
                if (arg4 < 0 | arg5 < 0 | arg4 ~= round(arg4) | arg5 ~= round(arg5))
                    errorcode = 4;
                    return;
                end
                rows    = arg4;
                columns = arg5;
            end
        end
    end

author: CelsoReyes


text:
% DESCRIPTION:
%       subscript to create the absorption variables
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 26th November 2010
%       last update - 11th February 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

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

% define the lossy derivative operators and proportionality coefficients
if strcmp(equation_of_state, 'absorbing')

    % make sure the operators are positive and real
    medium.alpha_coeff = abs(real(medium.alpha_coeff));
    medium.alpha_power = abs(real(medium.alpha_power));

    % convert the absorption coefficient to nepers.(rad/s)^-y.m^-1
    medium.alpha_coeff = db2neper(medium.alpha_coeff, medium.alpha_power);

    % compute the absorbing fractional Laplacian operator and coefficient
    if ~(isfield(medium, 'alpha_mode') && strcmp(medium.alpha_mode, 'no_absorption'))
        absorb_nabla1 = (kgrid.k).^(medium.alpha_power-2);
        absorb_nabla1(isinf(absorb_nabla1)) = 0;
        absorb_nabla1 = ifftshift(absorb_nabla1);
        absorb_tau = -2*medium.alpha_coeff.*medium.sound_speed.^(medium.alpha_power - 1);
    else
        absorb_nabla1 = 0;
        absorb_tau = 0;
    end

    % compute the dispersive fractional Laplacian operator and coefficient
    if ~(isfield(medium, 'alpha_mode') && strcmp(medium.alpha_mode, 'no_dispersion'))
        absorb_nabla2 = (kgrid.k).^(medium.alpha_power-1);
        absorb_nabla2(isinf(absorb_nabla2)) = 0;
        absorb_nabla2 = ifftshift(absorb_nabla2);
        absorb_eta = 2*medium.alpha_coeff.*medium.sound_speed.^(medium.alpha_power)*tan(pi*medium.alpha_power/2);
    else
        absorb_nabla2 = 0;
        absorb_eta = 0;
    end

    % pre-filter the absorption parameters if alpha_filter is defined (this
    % is used for time-reversal photoacoustic image reconstruction
    % with absorption compensation)
    if isfield(medium, 'alpha_filter');

        % update command line status
        disp('  filtering absorption variables...');

        % frequency shift the absorption parameters
        absorb_nabla1 = fftshift(absorb_nabla1);
        absorb_nabla2 = fftshift(absorb_nabla2);

        % apply the filter
        absorb_nabla1 = absorb_nabla1.*medium.alpha_filter;
        absorb_nabla2 = absorb_nabla2.*medium.alpha_filter;

        % shift the parameters back
        absorb_nabla1 = ifftshift(absorb_nabla1);
        absorb_nabla2 = ifftshift(absorb_nabla2);

    end

    % modify the sign of the absorption operators if alpha_sign is defined
    % (this is used for time-reversal photoacoustic image reconstruction
    % with absorption compensation)
    if isfield(medium, 'alpha_sign')
       if numel(medium.alpha_sign) == 2
           % if two parameters are given, apply the first to the absorption
           % parameter and the second to the disperion parameters
           absorb_tau = sign(medium.alpha_sign(1))*absorb_tau;
           absorb_eta = sign(medium.alpha_sign(2))*absorb_eta;
       else
           error('medium.alpha_sign must be given as a 2 element array controlling absorption and dispersion, respectively.');
       end
    end
end

repo: Aero-matlab


text:
function plot_time(catalogObject)
    %CATALOG.PLOT_TIME Plot magnitude and depth against time
    %   catalogObject.plot_time()

    % Glenn Thompson 2014/06/01

    symsize = get_symsize(catalogObject);
    timerange = catalogObject.gettimerange();
    xlims = [floor(timerange(1)) ceil(timerange(2))];

    % time-depth
    if all(isnan(catalogObject.depth))
        warning('No depth data to plot');
    else
        figure;
        set(gcf,'Color', [1 1 1]);
        subplot(2,1,1);
        scatter(catalogObject.otime, catalogObject.depth, symsize);
        set(gca, 'XLim', xlims);
        datetick('x');
        xlabel('Date');
        ylabel('Depth (km)');
        set(gca, 'YDir', 'reverse');
        grid on;

        % time-mag
        subplot(2,1,2);
    end
    if all(isnan(catalogObject.mag))
        warning('No magnitude data to plot');
    else
        scatter(catalogObject.otime, catalogObject.mag, symsize);
        %stem(catalogObject.otime, catalogObject.mag);
        set(gca, 'XLim', xlims);
        datetick('x');
        xlabel('Date');
        ylabel('Magnitude');
        grid on;
    end
end

author: geoscience-community-codes

repo: GISMO


text:
function varargout = isreal(varargin)
%ISREAL   Real-valued SPHEREFUN test.
%   ISREAL(F) returns logical true if F does not have an imaginary part and
%   false otherwise.
%
%   ~ISREAL(F) detects SPHEREFUN object that have an imaginary part even if it i
%   all zero.
%
%   Since only real-valued SPHEREFUNS are presently supported, this
%   function always returns a logical true.

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

[varargout{1:nargout}] = isreal@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [Model, Info] = linear_sparse_stepwise(X,Y,Model,parm)
%  Estimate linear weight matrix for input-output mapping
%     Automatic Relevance Prior for each input dimension
%     is imposed to get sparse weight matrix
%
% Notice !!
%     Delay embedding is done in this module,
%     then input vector should be original input data without embedding
%
%   [Model, Info] = linear_sparse_stepwise(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 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,'Npre_train')
	Npre_train = parm.Npre_train;
else
%	Npre_train = Ntrain;
	if Tall >= 2*M*D
		Npre_train = 0;
	elseif Tall >= M*D
		Npre_train = fix(Ntrain/2);
	else
		Npre_train = Ntrain;
	end
end

if Npre_train > Ntrain, Npre_train = Ntrain; 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',Tall)
fprintf('--- Total update iteration    = %d (%d)\n',Ntrain,Npre_train)

% 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 variance
XX  = sum(sum(X.^2,3),2)/(Tx*Ntrial);
XX  = repmat(XX', [1 D]);% 1 x Xdim*D

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);
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]);
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 = 1./A;
	% E = ( (Y-W*X)^2 +  W^2 * Ainv )/SY
    [W]  = weight_update_embed(X, dY, W, Tall*XX, Ainv, tau);
	[dY] = error_delay_time(X, Y, W, tau);

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

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

    % Weight variance
    SW   = 1./( Tall*XX + Ainv );

    % 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 per data
%    LP(k)  = - (0.5) * log_sy ;
    LP(k)  = - (0.5) * (log_sy + N*sum(SW.*XX));
    H(k)   = 0.5*( log_sw + log_a );
%    H(k)   = 0.5*( log_sw + log_a - sum(WW .*Ainv) );
    FE(k)  = LP(k) + H(k)/Tall;
    Err(k) = sum(dYY)/(N*SY0);

	% E = ( (Y-W*X)^2 +  W^2 * Ainv )/SY
%    SW = SY./( Tall*XX + Ainv );

    % Hyper parameter for weight variance (ARD)
	if mod(k,Nupdate) == 0,
	    % SW  = 1./( Tall*XX + Ainv );
	    % G_A = 1 - SW./A;
	    %     = (Tall*XX + Ainv - 1/A) * SW
%	    G_A = Tall.*XX.*SW;
	    G_A = Tall.*XX.*SW;
	    G_A = max((G_A), MINVAL);

		% N*A =  (W.^2)/SY + N * SW  ;
		if k <= Npre_train,
			% VB update rule (Stable)
			%  ARD for each weight
			A  = (WW./SY + N*SW + 2*Ta0*a0)/( N + 2*Ta0 );
		else
			% Accelerated update rule (Unstable)
		    A  = sqrt( A.* (WW./SY)./(G_A * N));
		    %A  = (WW + 2*Ta0*a0)./(G_A * N + 2*Ta0);
		end

		if space_ARD==1
			Am = mean(reshape(A,[M/D,D]),2);
			A  = repmat(Am', 1,D);
		end

	    % Prune small variance
	    if Prune == 1
		    % Find active input dimension
		    ix_act_old = ix_act;
		    ix_dim_old = ix_dim;

		    % recover all component
			A_all  = zeros(1,M_all);       % 1 x M
		    %  A_all(:,ix_act)  = A;
		    A_all(ix_act) = WW/max(WW);

		    % find active input dimension (absolute index)
		    ix_dim = find( sum(reshape(A_all,[Xdim,D]),2) > a_min );
		    ix_act = repmat(ix_dim, [1 D]) ...
		           + repmat([0:D-1]*Xdim, [length(ix_dim) 1]);
		    ix_act = ix_act(:);

		    Mnew   = length(ix_act);  		% # of effective input

		    if Mnew < M,
			    % convert to relative index
			    jx_act = trans_index(ix_act,ix_act_old,M_all);
			    jx_dim = trans_index(ix_dim,ix_dim_old,Xdim);

			    M   = Mnew;
			    A   = A(jx_act) ;  % N x M
			    W   = W(:,jx_act) ;  % N x M
			    SW  = SW(jx_act);  % N x M

			    X	 = X(jx_dim,:,:);  	 	% M x T
			    XX	 = XX(jx_act);  	 	% M x T
			end
	    end
	    % END of if Prune == 1

		A = max(A,MINVAL);

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

	Mhist(k) = M;

    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, M = %4d, err = %g, F = %g, SY= %g, H = %g\n', ...
               k, length(ix_act), Err(k), FE(k), SY, - 2*H(k)/log(Tall));
    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.SW   = SW; % = 1./(Tall*XX + Ainv)

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:
%PREX_CONFMAT PRTools example on confusion matrix, scatterplot and gridsize
%
% Prtools example code to show the use of confusion matrix,
% scatterplot and gridsize.

help prex_confmat

delfigs
randstate = randreset;
echo on
                % Load 8-class 2D problem
  randn('state',1);
  rand('state',1);
  a = gendatm;
                % Compute the Nearest Mean Classifier
  w = nmc(a);
                % Scatterplot
  figure;
  gridsize(30);
  scatterd(a,'legend');
                % Plot the classifier
  plotc(w);
  title([getname(a) ', Gridsize 30']);
                % Set higher gridsize
  gridsize(100);
  figure;
  scatterd(a,'legend');
  plotc(w);
  title([getname(a) ', Gridsize 100']);

                % Classify training set
  d = a*w;
                % Look at the confusion matrix and compare it to the scatterplot
  confmat(d);

echo off
randreset(randstate);
showfigs
c = num2str(gridsize);
disp(' ')
disp('   Classifier plots are inaccurate for small gridsizes. The standard');
disp('   value of 30 is chosen because of the speed, but it is too low to');
disp('   ensure good plots. Other gridsizes may be set by gridsize(n).')
disp('   Compare the two figures and appreciate the difference.')

author: marianux

repo: ecg-kit

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


text:
function [OfConstraint, OfAll] = determineFluxValuesOnBoundary(model, solution)
% This function determines the number of reactions in the flux distributions that are on
% the boundaries
%
% [OfConstraint, OfAll] = determineFluxValuesOnBoundary(model, solution)
%
% INPUT
% model         Model structure
% solution      Solution structure
%
% OUTPUT
% OfConstraint  Fraction of flux values that are on the lower and upper bounds
%               of all constrained reactions (assuming a minimum infinity
%               of -1,000,000 and a  maximum infinity of 1,000,000
% OfAll         Fraction of flux values that are on the lower and upper bounds
%               of all reactions in the model
%
% Ines Thiele, December 2018
%
%
%
minInf = -1000000;
maxInf = 1000000;
% find reactions that have flux values on upper bound
I = find(abs(solution.full)>1e-6);
% find all non-zero and non-inf bounds
J = find(abs(model.ub)~=0);
Ji = find(abs(model.ub)<=abs(minInf));
IJ = intersect(I,J);
IJ = intersect(IJ,Ji);
Usedub = model.ub(IJ);
UsedF = solution.full(IJ);
UsedR = model.rxns(IJ);
OnUB= length(find(abs(Usedub-UsedF)<1e-5));

% find reactions that have flux values on lower bound
I = find(abs(solution.full)>1e-6);
% find all non-zero and non-inf bounds
J = find(abs(model.lb)~=0);
Ji = find(abs(model.lb)<=abs(maxInf));
IJ = intersect(I,J);
IJ = intersect(IJ,Ji);
Usedlb = model.lb(IJ);
UsedF = solution.full(IJ);
UsedR = model.rxns(IJ);
OnLB= length(find(abs(Usedlb-UsedF)<1e-5));

% constraint reactions
minConstraints = length(intersect(find(model.lb>minInf),find(model.lb)));
maxConstraints =length(intersect(find(model.ub<maxInf),find(model.ub)));
PercentageConstraintRxns_model = (minConstraints + maxConstraints)*100/length(model.ub);
%percentage of all reactions
OfAll = (OnLB+OnUB)/length(model.rxns);
OfConstraint = (OnLB+OnUB)/(minConstraints + maxConstraints);

author: opencobra

repo: cobratoolbox


text:
function dag = sample_dag(P)
% SAMPLE_DAG Create a random directed acyclic graph with edge probabilities P(i,j)
% dag = sample_dag(P)
%
% This uses rejection sampling to reject graphs with directed cycles.

done = 0;
directed = 1;
iter = 1;
while ~done
  dag = binornd(1, P); % each edge is an indep Bernoulli (0/1) random variable
  dag = setdiag(dag, 0);
  done = acyclic(dag, directed);
  iter = iter + 1
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


function printMetrics(metrics, metricsInfo, dispHeader,dispMetrics,padChar)
% print metrics
%
% ...
%

% default names
if nargin==1
    metricsInfo.names.long = {'Recall','Precision','False Alarm Rate', ...
        'GT Tracks','Mostly Tracked','Partially Tracked','Mostly Lost', ...
        'False Positives', 'False Negatives', 'ID Switches', 'Fragmentations', ...
        'MOTA','MOTP', 'MOTA Log'};

    metricsInfo.names.short = {'Rcll','Prcn','FAR', ...
        'GT','MT','PT','ML', ...
        'FP', 'FN', 'IDs', 'FM', ...
        'MOTA','MOTP', 'MOTAL'};

    metricsInfo.widths.long = [6 9 16 9 14 17 11 15 15 11 14 5 5 8];
    metricsInfo.widths.short = [5 5 5 3 3 3 3 5 5 4 4 5 5 5];

    metricsInfo.format.long = {'.1f','.1f','.2f', ...
        'i','i','i','i', ...
        'i','i','i','i', ...
        '.1f','.1f','.1f'};

    metricsInfo.format.short=metricsInfo.format.long;
end

namesToDisplay=metricsInfo.names.long;
widthsToDisplay=metricsInfo.widths.long;
formatToDisplay=metricsInfo.format.long;

namesToDisplay=metricsInfo.names.short;
widthsToDisplay=metricsInfo.widths.short;
formatToDisplay=metricsInfo.format.short;

if nargin<3, dispHeader=1; end
if nargin<4
    dispMetrics=1:length(metrics);
end
if nargin<5
    padChar={' ',' ','|',' ',' ',' ','|',' ',' ',' ','| ',' ',' ',' '};
end

if dispHeader
    for m=dispMetrics
        printString=sprintf('fprintf(''%%%is%s'',char(namesToDisplay(m)))',widthsToDisplay(m),char(padChar(m)));
        eval(printString)
    end
    fprintf('\n');
end

for m=dispMetrics
    printString=sprintf('fprintf(''%%%i%s%s'',metrics(m))',widthsToDisplay(m),char(formatToDisplay(m)),char(padChar(m)));
    eval(printString)
end

% if standard, new line
if nargin<4
    fprintf('\n');
end


text:
classdef DOC8 < PROBLEM
% <multi> <real> <constrained>
% Benchmark MOP with constraints in both decision and objective spaces

%------------------------------- Reference --------------------------------
% Z. Liu and Y. Wang, Handling constrained multiobjective optimization
% problems with constraints in both the decision and objective spaces. IEEE
% Transactions on Evolutionary Computation, 2019, 23(5): 870-884.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 3;
            obj.D = 10;
            obj.lower    = [0 0 500 1000 5000 100 100 100 100 100];
            obj.upper    = [1 1 1000 2000 6000 500 500 500 500 500];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values and constraint violations
        function Population = Evaluation(obj,varargin)
            X = varargin{1};
            X = max(min(X,repmat(obj.upper,size(X,1),1)),repmat(obj.lower,size(X,1),1));
            g_temp = X(:, 3) + X(:, 4) + X(:, 5);
            g = g_temp-7049.2480205286 +1;
            PopObj(:,1) = (X(:,1).*X(:,2)).*g;
            PopObj(:,2) = (X(:,1).*(1 - X(:,2))).*g;
            PopObj(:,3) = (1-X(:,1)).*g;
            % Constraints in objective space
            c(:,1) = max( - (PopObj(:,3) - 0.4).*(PopObj(:,3) - 0.6), 0);
            % Constraints in decision space
            c(:,2) = -1 + 0.0025 * (X(:, 6) + X(:, 8));
            c(:,3) = -1 + 0.0025 * (X(:, 7) + X(:, 9) - X(:, 6));
            c(:,4) = -1 + 0.01 * (X(:, 10) - X(:, 7));
            c(:,5) = -X(:, 3).* X(:, 8) + 833.33252 * X(:, 6) + 100 * X(:, 3) - 83333.333;
            c(:,6) = -X(:, 4).* X(:, 9) + 1250 * X(:, 7) + X(:, 4).* X(:, 6) - 1250 * X(:, 6);
            c(:,7) = -X(:, 5).* X(:, 10) + 1250000 + X(:, 5).* X(:, 7) - 2500 * X(:, 7);
            Population  = SOLUTION(X,PopObj,c,varargin{2:end});
            obj.FE      = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
             R = UniformPoint(N,3);
             R(0.4<R(:,3)&R(:,3)<0.6,:) = [];
        end
        %% Generate the feasible region
        function R = GetPF(obj)
            a = linspace(0,1,20)';
            x = a*a';
            y = a*(1-a');
            z = (1-a)*ones(size(a'));
            z(0.4<z&z<0.6) = nan;
            R = {x,y,z};
        end
    end
end


text:
function hRD = plugin_rd_tx
% Reference design definition

% Call the common reference design definition function
hRD = AnalogDevices.adrv9364z7020.common.plugin_rd('ccbob_lvds','Tx');
AnalogDevices.adrv9364z7020.common.add_io(hRD, 'Tx');

author: analogdevicesinc


text:
function symb_pvec = sdisplay(pvec,symbolicname)
%SDISPLAY Symbolic display of SDPVAR expression
%
% Note that the symbolic display only work if all
% involved variables are explicitely defined as
% scalar variables.
%
% Variables that not are defined as scalars
% will be given the name ryv(i). ryv means
% recovered YALMIP variables, i indicates the
% index in YALMIP (i.e. the result from getvariables)
%
% If you want to change the generic name ryv, just
% pass a second string argument
%
% EXAMPLES
%  sdpvar x y
%  sdisplay(x^2+y^2)
%    ans =
%       'x^2+y^2'
%
%  t = sdpvar(2,1);
%  sdisplay(x^2+y^2+t'*t)
%    ans =
%      'x^2+y^2+ryv(5)^2+ryv(6)^2'

% Author Johan L�fberg
% $Id: sym.m,v 1.1 2005-02-22 16:50:11 johanl Exp $
allnames = {};
for pi = 1:size(pvec,1)
    for pj = 1:size(pvec,2)
        Y.type = '()';
        Y.subs = [{pi} {pj}];
        p = subsref(pvec,Y);

        if isnumeric(p)
            symb_p = num2str(p);
        else
            LinearVariables = depends(p);
            x = recover(LinearVariables);
            exponent_p = full(exponents(p,x));
            names = cell(length(x),1);
            for i = 1:length(names)
                names{i} = ['x' num2str(LinearVariables(i))];
                allnames{end+1} = names{i};
            end

            symb_p = '';
            if all(exponent_p(1,:)==0)
                symb_p = num2str(full(getbasematrix(p,0)));
                exponent_p = exponent_p(2:end,:);
            end

            for i = 1:size(exponent_p,1)
                coeff = full(getbasematrixwithoutcheck(p,i));
                switch coeff
                    case 1
                        coeff='+';
                    case -1
                        coeff = '-';
                    otherwise
                        if isreal(coeff)
                        if coeff >0
                            coeff = ['+' num2str2(coeff)];
                        else
                            coeff=[num2str2(coeff)];
                        end
                        else
                            coeff = ['+' '(' num2str2(coeff) ')' ];
                        end
                end
                symb_p = [symb_p coeff symbmonom(names,exponent_p(i,:))];
            end
            if symb_p(1)=='+'
                symb_p = symb_p(2:end);
            end
        end
        symb_pvec{pi,pj} = symb_p;
    end
end
allnames = unique(allnames);
for i = 1:length(allnames)
    evalin('caller',['syms ' allnames{i}]);
end

S = '';
for pi = 1:size(pvec,1)
    ss = '';
    for pj = 1:size(pvec,2)
        ss = [ss ' ' symb_pvec{pi,pj} ','];
    end
    S = [S ss ';'];
end
S = ['[' S ']']   ;
symb_pvec = evalin('caller',S);

function s = symbmonom(names,monom)
s = '';
for j = 1:length(monom)
    if abs( monom(j))>0
        s = [s names{j}];
        if monom(j)~=1
            s = [s '^' num2str(monom(j))];
        end
        s =[s '*'];
    end

end
if isequal(s(end),'*')
    s = s(1:end-1);
end

function s = num2str2(x)
s = num2str(full(x));
if isequal(s,'1')
    s = '';
end
if isequal(s,'-1')
    s = '-';
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function  [Itt,para]=color_transfer(Ic,It,ind)

[m1,n1,c]=size(Ic);
X=reshape(Ic,[m1*n1,c]);

[m2,n2,c]=size(It);
Y=reshape(It,[m2*n2,c]);

if nargin <3
   ind=[1:(m2*n2)]';
end

mu_x=mean(X);

S_x=(X-repmat(mu_x,[m1*n1,1]))'*(X-repmat(mu_x,[m1*n1,1]))/(m1*n1);
[Ux,Dx,~]=svd(S_x);
mu_y=mean(Y(ind,:));
S_y=(Y(ind,:)-repmat(mu_y,[size(ind,1),1]))'*(Y(ind,:)-repmat(mu_y,[size(ind,1),1]))/(size(ind,1));
[Uy,Dy,~]=svd(S_y);

A=Ux*diag(diag(Dx).^(0.5))*Ux'*Uy*diag(diag(Dy).^(-0.5))*Uy';
b=mu_x'-A*mu_y';

Z=A*Y'+repmat(b,[1,m2*n2]);
Z=Z';
para.A=A;
para.b=b;
%mu_z=mean(Z);
%S_z=(Z-repmat(mu_z,[m2*n2,1]))'*(Z-repmat(mu_z,[m2*n2,1]))/(m2*n2);
Itt=reshape(Z,[m2,n2,c]);
end

author: csjcai

repo: RealSR

save_path: github-repos/MATLAB/csjcai-RealSR


text:
function E = foe_energy_a(x, im, basis, W)
%FOE_ENERGY_A Compute value of FoE energy function for a specific a
%
%   E = foe_energy_a(x, im, W)
%
% Computes the enegery value of the Field of Experts model defined by
% parameters W and a at point x. The energy is returned in E.
%
% This function is only used to check the gradient FOE_ENERGY_GRAD_A.
%
%
% (C) Laurens van der Maaten, 2009
% Delft University of Technology

    % Compute energy
    a = x;
    E = foe_energy(im, basis, W, a);

author: lbasek

repo: image-denoising-benchmark