text:
% ------------------------------------------------------------------- %
% Matlab files listed in Appendix B of the following book by          %
% Xin-She Yang, Mathematical Modelling for Earth Sciences,            %
%               Dunedin Academic Presss, Edinburgh, UK, (2008).       %
% http://www.dunedinacademicpress.co.uk/ViewReviews.php?id=100        %
% ------------------------------------------------------------------- %

% ------------------------------------------------------------------- %
% Gauss Numerical Integration of I=\int_a^b f(x) dx
% by 7-point Gaussian Quadrature
% Programmed by Xin-She Yang (Cambridge University) @ 2008
% Usage: Gauss_quad(function_str,a,b)
%   E.g. Gauss_quad('(sin(x)./x).^2',0,pi);
% ------------------------------------------------------------------- %

function [I]=Gauss_quad(fstr,a,b)
format long;
if nargin<3,
 disp('Usage:Gauss_quad(integrand,a,b)');
 disp('E.g., Gauss_quad(''exp(-x^2)*2/sqrt(pi)'',-1,1);');
end
% Default function and values if no input
if nargin<1,
   help Gauss_quad.m;
   fstr='exp(-x^2)*2/sqrt(pi)';
   a=-1.0; b=1.0;
end
% Converting the input integrand, fstr, to a function f(x)
f=inline(vectorize(fstr));
% Seven-point integration scheme so zeta_1 to zeta_7
zeta=[-0.9491079123; -0.7415311855; -0.4058451513; 0.0;
       0.4058451513; 0.7415311855; 0.9491079123];
% Weighting coefficients
w=[0.1294849661; 0.2797053914; 0.3818300505; 0.4179591836;
   0.3818300505; 0.2797053914; 0.1294849661];
% Index for the seven points
Index=1:7;
Iunit=sum(w(Index).*f((b-a).*(zeta(Index)+1)/2+a));
I=(b-a)/2*Iunit;

disp(' '); disp(strcat('The integral is = ',num2str(I)));
% ------------------------end of program ---------------------------


text:
%% QUAD_FSU uses the fsuClusterMatlab command to run the QUAD computation.
%
%  Discussion:
%
%    The arguments to fsuClusterMatlab mean:
%
%    [] allows us to specify an output directory;
%    [] allows us to specify switches to the MOAB queue handler, such as more time.
%    'm" indicates this is an "MPI-like" job;
%    'w' indicates that the MATLAB session should wait for completion;
%    4 is the number of processors requested;
%    @quad_fun indicates the function we want to execute;
%    {} would contain input arguments if there were any.
%
%    The function "quad_fun" must correspond to a MATLAB M-file "quad_fun.m"
%    and that file must be in MATLAB's path.
%
%    The easiest way to ensure that the file is in MATLAB's path is to
%    create a subdirectory called "matlab" immediately below your home directory,
%    and place the function file there (as well as any other files needed by the job).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 March 2010
%
%  Author:
%
%    John Burkardt
%
  fsuClusterMatlab ( [], [], 'm', 'w', 4, @quad_fun, { 100000 } )

author: johannesgerer

repo: jburkardt-m


text:
function isloaded = loadGlobalNamespace()
   % loadGlobalNamespace handles the global variables and routines used within
   % wavefrom.
   %

   %global mep2dep dep2mep

   global WAVEFORM_HISTORY
   WAVEFORM_HISTORY = true; %keep track of waveform modifications

   %time converter functions
   % matlab date (# days since 0 JAN 0000)
   %mep2dep = inline ('(matepoch - 719529) * 24 * 3600','matepoch');
   % antelope epoch( # secs since 1 JAN 1970 )
   %dep2mep = inline ('(dep / 86400 + 719529)', 'dep');

   %default period for envelopes
   global DEFAULT_PERIOD
   DEFAULT_PERIOD = 10; %seconds

   global SECONDS_IN_DAY
   SECONDS_IN_DAY = 24 * 60 * 60;

   isloaded = true;
end

author: geoscience-community-codes

repo: GISMO


text:
function w = apply_filter(w, PARAMS)
    if ~exist('PARAMS','var')
        PARAMS.filterObj = filterobject('h',0.1,2);
    end
    for c=1:numel(w)
        try
            w(c) = filtfilt(PARAMS.filterObj, w(c));
        end
    end
end


text:
function Cube = CreateNewGame(hAxes,Cube)
% to create new game
n = 20; % number of random rotations
Axe = unidrnd(3,1,n); % axes
Side = unidrnd(2,1,n)*2-3; % sides
Direction = unidrnd(2,1,n)*2-3; % directions
for k=1:n, % random rotations
  Cube = RotateLayer(hAxes,Cube,Axe(k),Side(k),Direction(k),0);
end
PlotCube(hAxes,Cube);
return

author: Sable

repo: mcbench-benchmarks


text:
function my_stokes_GUI(nr,ntheta)
% close all;
% DEFINE function
r = linspace(0,1,100);
theta = linspace(0,2*pi,100);

[R,T]=meshgrid(r,theta);
x = R.*cos(T);
y = R.*sin(T);

% figure;
plot_curls(nr,ntheta,3);
camlight
caxis([-3 3])
zlim([-3 3])
grid on;
hold on;

% slider
set(gcf,'position',[410 150 560 626]);
figsize= get(gcf,'position');
SliderH = uicontrol('style','slider','position',[figsize(3)*0.95, figsize(4)*0.1 20 300],...
    'min', -3, 'max', 3,'SliderStep',[0.001, 0.01],'Value',3);

btn_on = uicontrol('style','pushbutton','String','LIGHT ON',...
    'position',[20 45 70 20],...
    'callback',@light_button_press);

btn_off = uicontrol('style','pushbutton','String','LIGHT OFF',...
    'position',[20 20 70 20],...
    'callback',@light_off_button_press);

addlistener(SliderH, 'Value','PostSet',@callbackfn);

    function callbackfn(source, eventdata)
        [az,el]=view;
        num = get(eventdata.AffectedObject, 'Value');
        cla
        plot_curls(nr,ntheta,num)
        view([az,el]);
        zlim([-3 3])

        caxis([-3 3])
        camlight
        grid on;
    end

    function light_off_button_press(source,event)
        delete(findall(gcf,'Type','light'))

    end

    function light_button_press(source,event)
        camlight;
    end
end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function res = normF2(A)
	%% ================== File info ==========================
	% Author		: Tiep Vu (http://www.personal.psu.edu/thv102/)
	% Time created	: Tue Jan 26 22:33:47 2016
	% Last modified	: Tue Jan 26 22:33:48 2016
	% Description	: square of the Frobenius Norm of A
	%% ================== end File info ==========================
	res = norm(A, 'fro')^2;
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function [metaboliteList, stoichCoeffList, revFlag] = parseRxnFormula(formula)
% Parses reaction formula into a list of metabolites and a
% list of `S` coefficients
%
% USAGE:
%
%    [metaboliteList, stoichCoeffList, revFlag] = parseRxnFormula(formula)
%
% INPUT:
%    formula:           Reaction formula, may contain symbols '+', '->', '<=>' in
%                       addition to stoichiometric coefficients and metabolite names
%                       examples:
%                       '0.01 cdpdag-SC[m] + 0.01 pg-SC[m]  -> 0.01 clpn-SC[m] + cmp[m] + h[m]' (irreversible reaction)
%                       'cit[c] + icit[x]  <=> cit[x] + icit[c] ' (reversible reaction)
%                       If no stoichiometric coefficient is provided, it is assumed
%                       to be = 1.
%                       Reaction formula should be a string, not a cell array
%
% OUTPUTS:
%    metaboliteList:     Cell array with metabolite names
%    stoichCoeffList:    List of `S` coefficients
%    revFlag:            Indicates whether the reaction is reversible (true) or not (false)
%
% EXAMPLE:
%
%    formula = '0.01 cdpdag-SC[m] + 0.01 pg-SC[m]  -> 0.01 clpn-SC[m] + cmp[m] + h[m]'
%
%    [metaboliteList, stoichCoeffList, revFlag] = parseRxnFormula(formula)
%
%    %metaboliteList = 'cdpdag-SC[m]'    'pg-SC[m]'    'clpn-SC[m]'    'cmp[m]'    'h[m]'
%    %stoichCoeffList = -0.01 -0.01 0.01 1 1
%    %revFlag = false
%
% .. Authors:
%       - Markus Herrgard 6/1/07
%       - Richard Que 1/25/10 Modified to handle '-->' and '<==>' as arrows
%         as well as reactionsformatted as '[compartment] : A --> C'.
%       - IT May 2012 Modified to handle '=>'

tokens = splitString(formula);

stoichCoeffList = [];
metaboliteList = {};
revFlag = true;

% Marks the start of a new stoichiometry + metabolite block
newMetFlag = true;
% Designates products vs reactants
productFlag = false;
compartment = '';
for i = 1:length(tokens)
    t = tokens{i};
    if strcmp(t(1),'[')
        %set compartment
        compartment = t;
    elseif strcmp(t,':')
        %Do nothing
    elseif strcmp(t,'+')
        % Do nothing
        newMetFlag = true;
    elseif strcmp(t,'->') || strcmp(t,'-->') || strcmp(t,'=>')
        % Irreversible
        revFlag = false;
        productFlag = true;
        newMetFlag = true;
    elseif strcmp(t,'<=>') || strcmp(t,'<==>')
        % Reversible
        revFlag = true;
        productFlag = true;
        newMetFlag = true;
    else
        sCoeff = str2double(t);
        if (~isnan(sCoeff))
            % Stoich coefficient
            if ~productFlag
                sCoeff = -sCoeff;
            end
            stoichCoeffList(end+1) = sCoeff;
            newMetFlag = false;
        else
            % Metabolite name
            metaboliteList{end+1} = strcat(t,compartment);
            if newMetFlag
                if ~productFlag
                    stoichCoeffList(end+1) = -1;
                else
                    stoichCoeffList(end+1) = 1;
                end
                newMetFlag = true;
            end
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
% This M-file is to change the settings of the LMS sensor in MATLAB to do
% the following:
% 1) Change the baudrate to 500k using RS422 communication setup via Fastcom 422/2-PCI-335 board
% 2) Permanently set this baudrate as data transfer rate for the sensor (To memorize even the power is turned off)

%% Create a serial object for MATLAB to talk with the board
s=serial('COM3','BaudRate',9600);
fopen(s)
%% At this time, switch on the Sick Sensor connected to the PC either using the RS-232 or RS-422 connection
% Read the initial response message from the LMS (30sec after the unit is turned on when the LED turns GREEN)
Initial_Response = fread(s,29,'uint8');
%% Status Request Telegram
% This telegram can be sent at any time to verify if you are able to
% communicate with the sensor.
% NOTE: A response telegram from the sensor if it starts with
% 06hex (6 decimal) indicates success
% 15hex (21 decimal) indicates failure in the sense that the sensor is unable to send any response for the telegram
Status_Chk = uint8([2 0 1 0 49 21 18]);
fwrite(s,Status_Chk);
Status_Ack = fread(s,161,'uint8');
%% Change to 500K BaudRate (high-speed mode) Telegram
Bd500k = uint8([2 0 2 0 32 72 88 8]);
fwrite(s,Bd500k);
Bd500k_Ack = fread(s,10,'uint8');
% change the board rate on the PC as well before any further communications
s.BaudRate = 9600; % the board assumes 9600 would actually correspond to highest possible BaudRate for the corresponding clock rate
%% close the serial object
fclose(s);

%% At this time, change the board settings from the PC by right clicking on
% My Computer->Properties->Hardware->Device Manager->Ports and change the clock rate to 8 MHz and check the "enable custom baudrate" option

%%  Reopen the serial object
fopen(s);
%% Change to Installation mode (password is hex SICK_LMS)
Install_Mode = uint8([2 0 10 0 32 0 83 73 67 75 95 76 77 83 190 197]);
fwrite(s,Install_Mode);
InstallMode_Ack = fread(s);
%% Set the current BaudRate as permanent
Perm_BdRate = uint8([2 0 2 0 102 1 157 78]);
fwrite(s,Perm_BdRate);
Perm_BdRate_Ack = fread(s);
%% Close the serial port object in MATLAB
fclose(s);


text:
classdef prtBrvVbMembershipModel < prtBrvMembershipModel

    methods (Abstract)

        y = conjugateVariationalAverageLogLikelihood(self, x)
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function create(obj,wrt)
% Create a NIFTI-1 file
% FORMAT create(obj)
% Write out the header information for the nifti object
%
% FORMAT create(obj,wrt)
% Also write out an empty image volume if wrt==1
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

%
% $Id: create.m 7147 2017-08-03 14:07:01Z spm $

for i=1:numel(obj)

    o = obj(i);

    if ~isa(o.dat,'file_array'), error('Data must be a file_array.'); end

    fname = o.dat.fname;
    if isempty(fname), error('No filename to write to.'); end

    %-Write NIFTI header
    sts = write_hdr_raw(fname, o.hdr, o.dat.dtype(end-1)=='B');
    if ~sts, error('Unable to write header for "%s".',fname); end

    %-Write extra information
    write_extras(fname,o.extras);

    %-Create an empty image file if necessary
    if nargin>1 && any(wrt==1)
        [pth,nam] = fileparts(fname);
        if any(strcmp(o.hdr.magic(1:3),{'n+1','n+2'}))
            ext = '.nii';
        else
            ext = '.img';
        end
        o.dat.fname = fullfile(pth,[nam ext]);

        initialise(o.dat);
    end

end

author: spm


text:
function test_fileformat_asa

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_read_sens ft_read_headmodel ft_read_headshape ft_read_mri read_asa read_asa_dip	read_asa_mri read_asa_vol read_asa_bnd read_asa_elc read_asa_msr

elcfile  = dccnpath('/home/common/matlab/fieldtrip/data/test/original/electrodes/asa/standard_primed.elc');
volfile  = dccnpath('/home/common/matlab/fieldtrip/data/test/original/headmodel/asa/standard.vol');
skinfile = dccnpath('/home/common/matlab/fieldtrip/data/test/original/headshape/asa/standard_skin_14038.vol');
mrifile  = dccnpath('/home/common/matlab/fieldtrip/data/test/original/mri/asa/standard.mri');

elec = ft_read_sens(elcfile);
vol  = ft_read_headmodel(volfile);
skin = ft_read_headshape(skinfile);
mri  = ft_read_mri(mrifile);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% The COBRAToolbox: testSampleCbModel.m
%
% Purpose:
%     - tests the sampleCbModel function using the E. coli Core Model
%

% Check Requirements
% quad minos and dqqMinos cannot be used due to parallel processing
% for some reason, matlab fails on this system if it runs in parallel.
solverPkgs = prepareTest('needsUnix',true, 'needsLP',true, 'excludeSolvers',{'mosek', 'dqqMinos','quadMinos','matlab','pdco'},'requiredToolboxes', {'statistics_toolbox'}); %TODO: Check, whether UNIX is really still required for the test.

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testSampleCbModel'));
cd(fileDir);
% define the samplers
samplers = {'CHRR','CHRR_EXP','ACHR'}; %'MFE'

% create a parallel pool (if possible)
try
    minWorkers = 2;
    myCluster = parcluster(parallel.defaultClusterProfile);

    if myCluster.NumWorkers >= minWorkers
        poolobj = gcp('nocreate');  % if no pool, do not create new one.
        if isempty(poolobj)
            parpool(minWorkers);  % launch minWorkers workers
        end
    end
catch
    %No pool
end

fprintf('   Testing sampleCbModel ... \n' );
for k = 1:length(solverPkgs.LP)
    fprintf(' -- Running testSampleCbModel using the solver interface: %s ... ', solverPkgs.LP{k});

    % set the solver
    solverOK = changeCobraSolver(solverPkgs.LP{k}, 'LP', 0);

    % Load model
    model = getDistributedModel('ecoli_core_model.mat');

    for i = 1:length(samplers)

        samplerName = samplers{i};

        switch samplerName
            case 'ACHR'
                fprintf('\nTesting the artificial centering hit-and-run (ACHR) sampler\n.');

                options.nFiles = 4;
                options.nStepsPerPoint = 5;
                options.nPointsReturned = 20;
                options.nPointsPerFile = 5;
                options.nFilesSkipped = 0;
                options.nWarmupPoints = 200;

                [modelSampling, samples, volume] = sampleCbModel(model, 'EcoliModelSamples', 'ACHR', options);

                % check if sample files created
                assert(exist('EcoliModelSamples_1.mat', 'file') == 2 && exist('EcoliModelSamples_2.mat', 'file') == 2 && exist('EcoliModelSamples_3.mat', 'file') == 2 && exist('EcoliModelSamples_4.mat', 'file') == 2)
                removedRxns = find(~ismember(model.rxns, modelSampling.rxns));
                assert(all(removedRxns == [26; 27; 29; 34; 45; 47; 52; 63]))

            case 'CHRR'
                fprintf('\nTesting the coordinate hit-and-run with rounding (CHRR) sampler\n.');

                options.nStepsPerPoint = 1;
                options.nPointsReturned = 10;
                options.toRound = 1;

                [modelSampling, samples, volume] = sampleCbModel(model, 'EcoliModelSamples', 'CHRR', options);

                assert(norm(samples) > 0)

            case 'CHRR_EXP'
                fprintf('\nTesting the coordinate hit-and-run with rounding (CHRR) sampler, with exponential target distribution.\n.');

                options.nStepsPerPoint = 1;
                options.nPointsReturned = 10;
                options.toRound = 1;
                numRxns = length(model.c);
                options.lambda = 0*model.c + 1;

                [modelSampling, samples, volume] = sampleCbModel(model, 'EcoliModelSamples', 'CHRR_EXP', options);

                assert(norm(samples) > 0)
        end
    end
end

% print a line for success of loop i
fprintf(' Done.\n');

% change the directory
cd(currentDir)


text:
% pop_rejmenu() - Main menu for rejecting trials in an EEG dataset
%
% Usage: >> pop_rejmenu(INEEG, typerej);
%
% Inputs:
%   INEEG      - input dataset
%   typerej    - data to reject on (0 = component activations;
%                1 = raw electrode data). {Default: 1 = reject on raw data}
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: eeglab(), pop_eegplot(), pop_eegthresh, pop_rejtrend()
% pop_rejkurt(), pop_jointprob(), pop_rejspec()

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function pop_rejmenu( EEG, icacomp );

if icacomp == 0
	if isempty( EEG.icasphere )
		disp('Error: you must first run ICA on the data'); return;
	end;
end;

if icacomp == 1 	rejtitle = 'Reject trials using data statistics - pop_rejmenu()'; tagmenu = 'rejtrialraw';
else            	rejtitle = 'Reject trials using component activity statistics - pop_rejmenu()'; tagmenu = 'rejtrialica';
end;

if ~isempty( findobj('tag', tagmenu))
	error('cannot open two identical windows; close the first one first');
end;

figure('visible', 'off', 'numbertitle', 'off', 'name', rejtitle, 'tag', tagmenu);

% definition of callbacks
% -----------------------
checkstatus = [ 'rejstatus = get( findobj(''parent'', gcbf, ''tag'', ''rejstatus''), ''value'');' ...
                'if rejstatus == 3,' ...
                '    EEG.reject.disprej = {};' ...
				'    if get( findobj(''parent'', gcbf, ''tag'', ''IManual''), ''value''), EEG.reject.disprej{1} = ''manual''; end;' ...
                '    if get( findobj(''parent'', gcbf, ''tag'', ''IThresh''), ''value''), EEG.reject.disprej{2} = ''thresh''; end;' ...
                '    if get( findobj(''parent'', gcbf, ''tag'', ''IConst''), ''value''), EEG.reject.disprej{3} = ''const''; end;' ...
                '    if get( findobj(''parent'', gcbf, ''tag'', ''IEnt''), ''value''), EEG.reject.disprej{4} = ''jp''; end;' ...
                '    if get( findobj(''parent'', gcbf, ''tag'', ''IKurt''), ''value''), EEG.reject.disprej{5} = ''kurt''; end;' ...
                '    if get( findobj(''parent'', gcbf, ''tag'', ''IFreq''), ''value''), EEG.reject.disprej{6} = ''freq''; end;' ...
                'end;' ...
                'rejstatus = rejstatus-1;' ]; % from 1-3 range, go to 0-2 range

if icacomp, ICAPREFIX = '';
else        ICAPREFIX = 'ica';
end;
% tmp_comall is used when returning from eegplot
tmp_comall =      [ 'set(findobj(''''parent'''', findobj(''''tag'''', ''''' tagmenu ...
					'''''), ''''tag'''', ''''mantrial''''), ''''string'''', num2str(sum(EEG.reject.' ICAPREFIX 'rejmanual)));' ...
                    'set(findobj(''''parent'''', findobj(''''tag'''', ''''' tagmenu ...
					'''''), ''''tag'''', ''''threshtrial''''), ''''string'''', num2str(sum(EEG.reject.' ICAPREFIX 'rejthresh)));' ...
                    'set(findobj(''''parent'''', findobj(''''tag'''', ''''' tagmenu ...
					'''''), ''''tag'''', ''''freqtrial''''), ''''string'''', num2str(sum(EEG.reject.' ICAPREFIX 'rejfreq)));' ...
                    'set(findobj(''''parent'''', findobj(''''tag'''', ''''' tagmenu ...
					'''''), ''''tag'''', ''''consttrial''''), ''''string'''', num2str(sum(EEG.reject.' ICAPREFIX 'rejconst)));' ...
                    'set(findobj(''''parent'''', findobj(''''tag'''', ''''' tagmenu ...
					'''''), ''''tag'''', ''''enttrial''''), ''''string'''', num2str(sum(EEG.reject.' ICAPREFIX 'rejjp)));' ...
                    'set(findobj(''''parent'''', findobj(''''tag'''', ''''' tagmenu ...
					'''''), ''''tag'''', ''''kurttrial''''), ''''string'''', num2str(sum(EEG.reject.' ICAPREFIX 'rejkurt)));' ];

cb_manual = [ 	checkstatus ...
				'pop_eegplot( EEG, ' int2str( icacomp ) ', rejstatus, 0,''' tmp_comall ''');' ...
				'clear rejstatus;' ];

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

cb_compthresh = [ '  posthresh = get( findobj(''parent'', gcbf, ''tag'', ''threshpos''), ''string'' );' ...
                  '  negthresh = get( findobj(''parent'', gcbf, ''tag'', ''threshneg''), ''string'' );' ...
                  '  startime  = get( findobj(''parent'', gcbf, ''tag'', ''threshstart''), ''string'' );' ...
                  '  endtime   = get( findobj(''parent'', gcbf, ''tag'', ''threshend''), ''string'' );' ...
                  '  elecrange = get( findobj(''parent'', gcbf, ''tag'', ''threshelec''), ''string'' );' ...
                  checkstatus ...
				  '[EEG Itmp LASTCOM] = pop_eegthresh( EEG,' int2str(icacomp) ...
				  ', elecrange, negthresh, posthresh, str2num(startime)/1000, str2num(endtime)/1000, rejstatus, 0,''' tmp_comall ''');' ...
				  'EEG = eegh(LASTCOM, EEG);' ...
				  'clear com Itmp elecrange posthresh negthresh startime endtime rejstatus;' ];

%                  '  set(findobj(''parent'', gcbf, ''tag'', ''threshtrial''), ''string'', num2str(EEG.trials - length(Itmp)));' ...
%				  'eegh(LASTCOM);' ...
%				  'clear Itmp elecrange posthresh negthresh startime endtime rejstatus;' ];

% -----------------------------------------------------
cb_compfreq   = [ '  posthresh = get( findobj(''parent'', gcbf, ''tag'', ''freqpos''), ''string'' );' ...
                  '  negthresh = get( findobj(''parent'', gcbf, ''tag'', ''freqneg''), ''string'' );' ...
                  '  startfreq = get( findobj(''parent'', gcbf, ''tag'', ''freqstart''), ''string'' );' ...
                  '  endfreq   = get( findobj(''parent'', gcbf, ''tag'', ''freqend''), ''string'' );' ...
                  '  elecrange = get( findobj(''parent'', gcbf, ''tag'', ''freqelec''), ''string'' );' ...
                  checkstatus ...
				  '[EEG Itmp LASTCOM] = pop_rejspec( EEG,' int2str(icacomp) ...
				  ', elecrange, negthresh, posthresh, startfreq, endfreq, rejstatus, 0,''' tmp_comall ''');' ...
				  'EEG = eegh(LASTCOM, EEG);' ...
				  'clear Itmp elecrange posthresh negthresh startfreq endfreq rejstatus;' ];

% -----------------------------------------------------
cb_compconstrej = [ '  minslope   = get( findobj(''parent'', gcbf, ''tag'', ''constpnts''), ''string'' );' ...
					'  minstd    = get( findobj(''parent'', gcbf, ''tag'', ''conststd''), ''string'' );' ...
					'  elecrange = get( findobj(''parent'', gcbf, ''tag'', ''constelec''), ''string'' );' ...
                    checkstatus ...
				    '[rej LASTCOM] = pop_rejtrend(EEG,' int2str(icacomp) ', elecrange, ''' ...
					int2str(EEG.pnts) ''', minslope, minstd, rejstatus, 0,''' tmp_comall ''');' ...
				    'EEG = eegh(LASTCOM, EEG);' ...
				    'clear rej elecrange minslope minstd rejstatus;' ];

% -----------------------------------------------------
cb_compenthead = [ '  locthresh  = get( findobj(''parent'', gcbf, ''tag'', ''entloc''), ''string'' );', ...
				   '  globthresh = get( findobj(''parent'', gcbf, ''tag'', ''entglob''), ''string'' );', ...
				   '  elecrange  = get( findobj(''parent'', gcbf, ''tag'', ''entelec''), ''string'' );' ];
cb_compenttail = [ '  set( findobj(''parent'', gcbf, ''tag'', ''entloc''), ''string'', num2str(locthresh) );', ...
				   '  set( findobj(''parent'', gcbf, ''tag'', ''entglob''), ''string'', num2str(globthresh) );', ...
				   '  set( findobj(''parent'', gcbf, ''tag'', ''enttrial''), ''string'', num2str(nrej) );' ...
				   'EEG = eegh(LASTCOM, EEG);' ...
				   'clear nrej elecrange locthresh globthresh rejstatus;' ];
cb_compentplot =  [ cb_compenthead ...
				    '[EEG locthresh globthresh nrej LASTCOM] = pop_jointprob( EEG, ' int2str(icacomp) ...
					', elecrange, locthresh, globthresh, 0, 0,''' tmp_comall ''');', ...
				    cb_compenttail ];
cb_compentcalc =  [ cb_compenthead ...
				    '[EEG locthresh globthresh nrej LASTCOM] = pop_jointprob( EEG, ' int2str(icacomp) ...
					', [ str2num(elecrange) ], str2num(locthresh), str2num(globthresh), 0, 0);', ...
 				    cb_compenttail ];
cb_compenteeg  =  [ cb_compenthead ...
                    checkstatus ...
				    '[EEG locthresh globthresh nrej LASTCOM] = pop_jointprob( EEG, ' int2str(icacomp) ...
					', elecrange, locthresh, globthresh, rejstatus, 0, 1,''' tmp_comall ''');', ...
 				    cb_compenttail ];

% -----------------------------------------------------
cb_compkurthead =[ '  locthresh  = get( findobj(''parent'', gcbf, ''tag'', ''kurtloc''), ''string'' );', ...
				   '  globthresh = get( findobj(''parent'', gcbf, ''tag'', ''kurtglob''), ''string'' );', ...
				   '  elecrange  = get( findobj(''parent'', gcbf, ''tag'', ''kurtelec''), ''string'' );' ];
cb_compkurttail =[ '  set( findobj(''parent'', gcbf, ''tag'', ''kurtloc''), ''string'', num2str(locthresh) );', ...
				   '  set( findobj(''parent'', gcbf, ''tag'', ''kurtglob''), ''string'', num2str(globthresh) );', ...
				   '  set( findobj(''parent'', gcbf, ''tag'', ''kurttrial''), ''string'', num2str(nrej) );' ...
				   'EEG = eegh(LASTCOM, EEG);' ...
				   'clear nrej elecrange locthresh globthresh rejstatus;' ];
cb_compkurtplot = [ cb_compkurthead ...
				    '[EEG locthresh globthresh nrej LASTCOM] = pop_rejkurt( EEG, ' int2str(icacomp) ...
					', elecrange, locthresh, globthresh, 0, 0,''' tmp_comall ''');', ...
				    cb_compkurttail ];
cb_compkurtcalc = [ cb_compkurthead ...
				    '[EEG locthresh globthresh nrej LASTCOM] = pop_rejkurt( EEG, ' int2str(icacomp) ...
					', [ str2num(elecrange) ], str2num(locthresh), str2num(globthresh), 0, 0);', ...
 				    cb_compkurttail ];
cb_compkurteeg  = [ cb_compkurthead ...
                    checkstatus ...
				    '[EEG locthresh globthresh nrej LASTCOM] = pop_rejkurt( EEG, ' int2str(icacomp) ...
					', elecrange, locthresh, globthresh, rejstatus, 0, 1,''' tmp_comall ''');', ...
 				    cb_compkurttail ];

% -----------------------------------------------------
cb_reject =      [ 'set( findobj(''parent'', gcbf, ''tag'', ''rejstatus''), ''value'', 3);' ... % force status to 3
				    checkstatus ...
				   '[EEG LASTCOM] = eeg_rejsuperpose(EEG,' int2str(icacomp) ',1,1,1,1,1,1,1); EEG = eegh(LASTCOM, EEG);' ...
                   'if isempty(find(EEG.reject.rejglobal)),' ...
                   '   warndlg2(strvcat(''No epoch selected...'', ''When using thresholding, click update'',''marks in the EEG plotting window''));' ...
				   'else,' ...
                   '   [EEG LASTCOM] = pop_rejepoch( EEG, EEG.reject.rejglobal, 1);' ...
				   '   if ~isempty(LASTCOM), ' ...
				   '      EEG = eegh(LASTCOM, EEG); [ALLEEG EEG CURRENTSET LASTCOM] = pop_newset(ALLEEG, EEG, CURRENTSET); eegh(LASTCOM);' ...
				   '   end; eeglab redraw; close(gcbf);' ...
                   'end;' ];

cb_clear =       [ 'close gcbf; EEG = rmfield( EEG, ''reject''); EEG.reject.rejmanual = [];' ...
				   'EEG=eeg_checkset(EEG); pop_rejmenu(' inputname(1) ',' int2str(icacomp) ');' ];

cb_close =       [ 'close gcbf;' ...
				   'disp(''Marks stored in dataset'');' ...
				   '[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);' ...
				   'eegh(''[ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);'');'];

lisboxoptions = { 'string', [ 'Show only the new trials marked for rejection by the measure selected above|' ...
                  'Show previous and new trials marked for rejection by the measure selected above|' ...
                  'Show all trials marked for rejection by the measure selected above or checked below'], 'tag', 'rejstatus', 'value', 1, 'callback', ...
                  [ 'if get(gcbo, ''value'') == 3,' ...
                    '   set(findobj(''parent'', gcbf, ''style'', ''checkbox''), ''enable'', ''on'');' ...
                    'else' ...
                    '   set(findobj(''parent'', gcbf, ''style'', ''checkbox''), ''enable'', ''off'');' ...
                    'end;' ] };

	chanliststr = [fastif(icacomp,'Electrode(s)','Component(s)') ];
	chanlistval = fastif(icacomp, [ '1:' int2str(EEG.nbchan) ], [ '1:' int2str(size(EEG.icaweights,1)) ]);

	% assess previous rejections
	% --------------------------
	sizeman = 0;
	sizethresh = 0;
	sizetrend = 0;
	sizejp = 0;
	sizekurt = 0;
	sizespec = 0;
	if icacomp == 1
		if ~isempty(EEG.reject.rejmanual), sizeman   = length(find(EEG.reject.rejmanual)); end;
		if ~isempty(EEG.reject.rejconst),  sizetrend = length(find(EEG.reject.rejconst)); end;
		if ~isempty(EEG.reject.rejjp),     sizejp    = length(find(EEG.reject.rejjp)); end;
		if ~isempty(EEG.reject.rejkurt),   sizekurt  = length(find(EEG.reject.rejkurt)); end;
		if ~isempty(EEG.reject.rejfreq),   sizespec  = length(find(EEG.reject.rejfreq)); end;
	else
		if ~isempty(EEG.reject.icarejmanual), sizeman   = length(find(EEG.reject.icarejmanual)); end;
		if ~isempty(EEG.reject.icarejconst),  sizetrend = length(find(EEG.reject.icarejconst)); end;
		if ~isempty(EEG.reject.icarejjp),     sizejp    = length(find(EEG.reject.icarejjp)); end;
		if ~isempty(EEG.reject.icarejkurt),   sizekurt  = length(find(EEG.reject.icarejkurt)); end;
		if ~isempty(EEG.reject.icarejfreq),   sizespec  = length(find(EEG.reject.icarejfreq)); end;
	end;

    stdl = [0.25 1.2 0.8 1.2 0.8]; % standard line
	titl = [0.9 0.18 1.55]; % title line
    geometry = { 	[0.883 0.195 .2 .45 .4 .4] ...
					[1] titl stdl stdl stdl stdl ...
					[1] titl stdl stdl stdl  ...
					[1] titl stdl stdl stdl  ...
					[1] titl stdl stdl stdl  ...
					[1] titl stdl stdl stdl stdl ...
					[1] [1] [1] [1 1 1] [1 1 1]  ...
					[1] [1 1 1]};

	listui = {{'Style', 'text', 'string', 'Mark trials by appearance', 'fontweight', 'bold' }, ...
    {  'Style', 'pushbutton', 'string', '', 'tag', 'butmanual', ...
	   'callback', [ 'tmpcolor = uisetcolor(EEG.reject.rejmanualcol); if length(tmpcolor) ~= 1,' ...
	   'EEG.reject.rejmanualcol=tmpcolor; set(gcbo, ''backgroundcolor'', tmpcolor); end; clear tmpcolor;'] },...
	{ } { 'Style', 'pushbutton', 'string', fastif(icacomp,'Scroll Data','Scroll Acts.'), 'callback', cb_manual }, ...
	{ 'Style', 'text', 'string', 'Marked trials' }, ...
	{ 'Style', 'text', 'string', int2str(sizeman), 'tag', 'mantrial' }, ...
	{ }, ...
	... % ---------------------------------------------------------------------------
	{ 'Style', 'text', 'string', 'Find abnormal values', 'fontweight', 'bold' }, ...
    {  'Style', 'pushbutton', 'string', '', 'tag', 'butthresh',  ...
	   'callback', [ 'tmpcolor = uisetcolor(EEG.reject.rejthreshcol); if length(tmpcolor) ~= 1,' ...
	   'EEG.reject.rejthreshcol=tmpcolor; set(gcbo, ''backgroundcolor'', tmpcolor); end; clear tmpcolor;'] }, { },...
	...
	{ }, { 'Style', 'text', 'string', ['Upper limit(s) ' fastif(icacomp,'(uV)', '(std. dev.)')] }, ...
	{ 'Style', 'edit', 'string', '25', 'tag', 'threshpos' }, ...
	{ 'Style', 'text', 'string', ['Lower limit(s) ' fastif(icacomp,'(uV)', '(std. dev.)')] }, ...
	{ 'Style', 'edit', 'string', '-25', 'tag', 'threshneg' }, ...
	...
	{ }, { 'Style', 'text', 'string', 'Start time(s) (ms)' }, ...
	{ 'Style', 'edit', 'string', int2str(EEG.xmin*1000), 'tag', 'threshstart' }, ...
	{ 'Style', 'text', 'string', 'Ending time(s) (ms)' }, ...
	{ 'Style', 'edit', 'string', int2str(EEG.xmax*1000), 'tag', 'threshend' }, ...
	...
	{ }, { 'Style', 'text', 'string', chanliststr }, ...
	{ 'Style', 'edit', 'string', chanlistval, 'tag', 'threshelec' }, ...
	{ 'Style', 'text', 'string', 'Currently marked trials' }, ...
	{ 'Style', 'text', 'string', int2str(sizethresh), 'tag', 'threshtrial' }, ...
	...
	{ }, { 'Style', 'pushbutton', 'string', 'Calc / Plot', 'callback', cb_compthresh }, ...
	{ }, { },{ 'Style', 'pushbutton', 'string', 'HELP', 'callback', 'pophelp(''pop_eegthresh'');' }, ...
	... % ---------------------------------------------------------------------------
	{ }, { 'Style', 'text', 'string', 'Find abnormal trends', 'fontweight', 'bold' }, ...
    {  'Style', 'pushbutton', 'string', '', 'tag', 'buttrend', ...
	   'callback', [ 'tmpcolor = uisetcolor(EEG.reject.rejconstcol); if length(tmpcolor) ~= 1,' ...
	   'EEG.reject.rejconstcol=tmpcolor; set(gcbo, ''backgroundcolor'', tmpcolor); end; clear tmpcolor;'] }, { },...
	...
	{ }, { 'Style', 'text', 'string', ['Max slope ' fastif(icacomp, ...
           '(uV/epoch)', ...
           '(std. dev./epoch)') ] }, ...
	{ 'Style', 'edit', 'string', '50', 'tag', 'constpnts' }, ...
	{ 'Style', 'text', 'string', 'R-squared limit (0 to 1)' }, ...
	{ 'Style', 'edit', 'string', '0.3', 'tag', 'conststd'  }, ...
    ...
	{ }, { 'Style', 'text', 'string', chanliststr }, ...
	{ 'Style', 'edit', 'string', chanlistval, 'tag', 'constelec'  }, ...
	{ 'Style', 'text', 'string', 'Currently marked trials' }, ...
	{ 'Style', 'text', 'string', int2str(sizetrend), 'tag', 'consttrial'  }, ...
	...
	{ }, { 'Style', 'pushbutton', 'string', 'Calc / Plot', 'callback', cb_compconstrej }, ...
	{ }, { }, { 'Style', 'pushbutton', 'string', 'HELP', 'callback', 'pophelp(''pop_rejtrend'');' }, ...
	...  % ---------------------------------------------------------------------------
	{ }, { 'Style', 'text', 'string', 'Find improbable data', 'fontweight', 'bold' }, ...
    {  'Style', 'pushbutton', 'string', '',  'tag', 'butjp', ...
	   'callback', [ 'tmpcolor = uisetcolor(EEG.reject.rejjpcol); if length(tmpcolor) ~= 1,' ...
	   'EEG.reject.rejjpcol=tmpcolor; set(gcbo, ''backgroundcolor'', tmpcolor); end; clear tmpcolor;'] }, { },...
	...
	{ }, { 'Style', 'text', 'string', fastif(icacomp, 'Single-channel limit (std. dev.)', 'Single-comp. limit (std. dev.)') }, ...
	{ 'Style', 'edit', 'string', fastif(icacomp, '5', '20'), 'tag', 'entloc' }, ...
	{ 'Style', 'text', 'string', fastif(icacomp, 'All channels limit (std. dev.)', 'All comp. limit (std. dev.)') }, ...
	{ 'Style', 'edit', 'string', '5', 'tag', 'entglob' }, ...
    ...
	{ }, { 'Style', 'text', 'string', chanliststr }, ...
	{ 'Style', 'edit', 'string', chanlistval, 'tag', 'entelec' }, ...
	{ 'Style', 'text', 'string', 'Currently marked trials' }, ...
	{ 'Style', 'text', 'string', int2str(sizejp), 'tag', 'enttrial' }, ...
	...
	{ }, { 'Style', 'pushbutton', 'string', 'Calculate', 'callback', cb_compentcalc }, ...
	{ 'Style', 'pushbutton', 'string', 'Scroll Data', 'callback', cb_compenteeg }, ...
	{ 'Style', 'pushbutton', 'string', 'PLOT', 'callback', cb_compentplot }, ...
	{ 'Style', 'pushbutton', 'string', 'HELP', 'callback', 'pophelp(''pop_jointprob'');' }, ...
	...  % ---------------------------------------------------------------------------
	{ }, { 'Style', 'text', 'string', 'Find abnormal distributions', 'fontweight', 'bold' }, ...
    {  'Style', 'pushbutton', 'string', '', 'tag', 'butkurt', ...
	   'callback', [ 'tmpcolor = uisetcolor(EEG.reject.rejkurtcol); if length(tmpcolor) ~= 1,' ...
	   'EEG.reject.rejkurtcol=tmpcolor; set(gcbo, ''backgroundcolor'', tmpcolor); end; clear tmpcolor;'] }, { },...
	...
	{ }, { 'Style', 'text', 'string', fastif(icacomp, 'Single-channel limit (std. dev.)', 'Single-comp. limit (std. dev.)') }, ...
	{ 'Style', 'edit', 'string', fastif(icacomp, '5', '20'), 'tag', 'kurtloc' }, ...
	{ 'Style', 'text', 'string', fastif(icacomp, 'All channels limit (std. dev.)', 'All comp. limit (std. dev.)') }, ...
	{ 'Style', 'edit', 'string', '5', 'tag', 'kurtglob' }, ...
    ...
	{ }, { 'Style', 'text', 'string', chanliststr }, ...
	{ 'Style', 'edit', 'string', chanlistval, 'tag', 'kurtelec' }, ...
	{ 'Style', 'text', 'string', 'Currently marked trials' }, ...
	{ 'Style', 'text', 'string', int2str(sizekurt), 'tag', 'kurttrial' }, ...
	...
	{ }, { 'Style', 'pushbutton', 'string', 'Calculate', 'callback', cb_compkurtcalc }, ...
	{ 'Style', 'pushbutton', 'string', 'Scroll Data', 'callback', cb_compkurteeg }, ...
	{ 'Style', 'pushbutton', 'string', 'PLOT', 'callback', cb_compkurtplot }, ...
	{ 'Style', 'pushbutton', 'string', 'HELP', 'callback', 'pophelp(''pop_rejkurt'');' }, ...
	... % ---------------------------------------------------------------------------
	{ }, { 'Style', 'text', 'string', 'Find abnormal spectra (slow)', 'fontweight', 'bold' }, ...
    {  'Style', 'pushbutton', 'string', '', 'tag', 'butspec', ...
	   'callback', [ 'tmpcolor = uisetcolor(EEG.reject.rejfreqcol); if length(tmpcolor) ~= 1,' ...
	   'EEG.reject.rejfreqcol=tmpcolor; set(gcbo, ''backgroundcolor'', tmpcolor); end; clear tmpcolor;'] }, { },...
	...
	{ }, { 'Style', 'text', 'string', 'Upper limit(s) (dB)' }, ...
	{ 'Style', 'edit', 'string', '25', 'tag', 'freqpos' }, ...
	{ 'Style', 'text', 'string', 'Lower limit(s) (dB)' }, ...
	{ 'Style', 'edit', 'string', '-25', 'tag', 'freqneg' }, ...
	...
	{ }, { 'Style', 'text', 'string', 'Low frequency(s) (Hz)' }, ...
	{ 'Style', 'edit', 'string', '0', 'tag', 'freqstart' }, ...
	{ 'Style', 'text', 'string', 'High frequency(s) (Hz)' }, ...
	{ 'Style', 'edit', 'string', '50', 'tag', 'freqend' }, ...
	...
	{ }, { 'Style', 'text', 'string', chanliststr }, ...
	{ 'Style', 'edit', 'string', chanlistval, 'tag', 'freqelec' }, ...
	{ 'Style', 'text', 'string', 'Currently marked trials' }, ...
	{ 'Style', 'text', 'string', int2str(sizespec), 'tag', 'freqtrial' }, ...
	...
	{ }, { 'Style', 'pushbutton', 'string', 'Calc / Plot', 'callback', cb_compfreq }, ...
	{ }, { },{ 'Style', 'pushbutton', 'string', 'HELP', 'callback', 'pophelp(''pop_rejspec'');' }, ...
    ...
    {}, ...  % ---------------------------------------------------------------------------
	...
	{ 'Style', 'text', 'string', 'Plotting options', 'fontweight', 'bold' }, ...
	{ 'style', 'listbox', lisboxoptions{:}, 'value', 3 }, ...
	{ 'style', 'checkbox', 'String', 'Abnormal appearance', 'tag', 'IManual', 'value', 1, ...
      'callback', 'set(gcbo, ''value'', 1); disp(''This checkbox must be set so that you may mark/unmark data epochs'');'}, ...
	{ 'style', 'checkbox', 'String', 'Abnormal values', 'tag', 'IThresh', 'value', 1}, ...
	{ 'style', 'checkbox', 'String', 'Abnormal trends', 'tag', 'IConst', 'value', 1}, ...
	{ 'style', 'checkbox', 'String', 'Improbable epochs', 'tag', 'IEnt', 'value', 1}, ...
	{ 'style', 'checkbox', 'String', 'Abnormal distributions', 'tag', 'IKurt', 'value', 1}, ...
	{ 'style', 'checkbox', 'String', 'Abnormal spectra', 'tag', 'IFreq', 'value', 1}, ...
	...
	{ }, ...
 	{ 'Style', 'pushbutton', 'string', 'CLOSE (KEEP MARKS)', 'callback', cb_close }, ...
 	{ 'Style', 'pushbutton', 'string', 'CLEAR ALL MARKS', 'callback', cb_clear  }, ...
	{ 'Style', 'pushbutton', 'string', 'REJECT MARKED TRIALS', 'callback', cb_reject }};

	allh = supergui( gcf, geometry, [], listui{:});

%	{ 'style', 'checkbox', 'String', ['Include ' fastif(icacomp, 'ica data', 'raw data')], 'tag', 'IOthertype', 'value', 1}, { }, ...
set(gcf, 'userdata', { allh });
set(findobj('parent', gcf', 'tag', 'butmanual'), 'backgroundcolor', EEG.reject.rejmanualcol);
set(findobj('parent', gcf', 'tag', 'butthresh'), 'backgroundcolor', EEG.reject.rejthreshcol);
set(findobj('parent', gcf', 'tag', 'buttrend'),  'backgroundcolor', EEG.reject.rejconstcol);
set(findobj('parent', gcf', 'tag', 'butjp'),     'backgroundcolor', EEG.reject.rejjpcol);
set(findobj('parent', gcf', 'tag', 'butkurt'),   'backgroundcolor', EEG.reject.rejkurtcol);
set(findobj('parent', gcf', 'tag', 'butspec'),   'backgroundcolor', EEG.reject.rejfreqcol);
set( findobj('parent', gcf, 'tag', 'rejstatus'), 'style', 'popup');

author: PatternRecognition

repo: OpenBMI


text:
function gf = gradspars(Hkt_r, varargin)
% function gf = gradddivH(Hkt_r, varargin)
% gf is 1 x k*t row vector
% Hkt_r = reshape(Hkt,1,k*t) -> row vector
% Vxt = varargin{1};  %data
% Wxk = varargin{2};  %W matrix
% Wxk_fix = varargin{3}; %fixed part of the Wxk matrix (e.g. background) ->rows
% Hkt_fix = varargin{4}; %fixed part (lines) of the H matrix (e.g. background)

alphaH=1; %for now....
Vxt = varargin{1};  %data
Wxk_tmp = varargin{2};  %W matrix
Wxk_fix = varargin{3}; %fixed part of the Wxk matrix (e.g. background) ->rows
Hkt_fix = varargin{4}; %fixed part (lines) of the H matrix (e.g. background)
peval = varargin{5}; %parameters

if ~isfield(peval, 'w_lambda') peval.w_lambda=0; end

t=size(Vxt,2);
k=length(peval.h_dovec);

Hkt_tmp = exp(reshape(Hkt_r,k,t));

Wxk = zeros(peval.numpix, peval.ncomp+1);
Hkt = zeros(peval.ncomp+1, peval.nt);

Wxk(:,peval.w_dovec)=Wxk_tmp;
Hkt(peval.h_dovec,:)=Hkt_tmp;

Wxk(:,peval.w_fixvec)=Wxk_fix;
Hkt(peval.h_fixvec,:)=Hkt_fix;

deltasum=sum(sum(Wxk_tmp,1))-k;
if and(~isempty(Wxk_tmp), abs(deltasum)>10^-6)
    error('Wxk is not correctly normalized! (sum(Wxk_tmp,1)<>1)\n sum(Wxk_tmp,1)=%f',deltasum)
end

% ap = peval.alphapen;
% sumH_t = sum(Hkt(peval.h_dovec,:),2);
sumH_t = sum(Hkt,2);
sumH_t_sq = sum(sumH_t.^2);
sumH = sum(sumH_t);
nh = length(sumH_t);

sparsness = (sqrt(nh)-sumH/(sqrt(sumH_t_sq)))/(sqrt(nh)-1);
% gradspars=(1/(sqrt(nh)-1))*(sumH/(sumH_t_sq)^1.5*sumH_t - 1/(sqrt(sumH_t_sq)));
gradspars= (1/(sqrt(nh)-1)*(sumH_t*(sumH)-sumH_t_sq)/(sumH_t_sq^1.5));
gradspars_kt = repmat(gradspars, 1, t);

% fkt = ddivergence(Vxt,Wxk*Hkt);
% gfkt = (1-Wxk'*(Vxt./(Wxk*Hkt)))*alphaH.*Hkt; %d/dh(d-divergence)
% gfkt = (1-Wxk'*(Vxt./(Wxk*Hkt))).*Hkt; %d/dh(d-divergence)
%(fg)'=(f'g-fg')/g^2;
% gfkt_all=(gfkt.*sparsness-fkt.*gradspars_kt)./(sparsness).^2;
% gfkt_all = gradspars_kt.*Hkt(peval.h_dovec, :); %because of exp...
gfkt_all = gradspars_kt.*Hkt; %because of exp...
% one is tehre because Wxt is normalized: sum(Wxt,1)=1
gf=reshape(gfkt_all(peval.h_dovec,:),1,k*t); %making row vector

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%  Save NIFTI dataset. Support both *.nii and *.hdr/*.img file extension.
%  If file extension is not provided, *.hdr/*.img will be used as default.
%
%  Usage: save_nii(nii, filename, [old_RGB])
%
%  nii.hdr - struct with NIFTI header fields (from load_nii.m or make_nii.m)
%
%  nii.img - 3D (or 4D) matrix of NIFTI data.
%
%  filename - NIFTI file name.
%
%  old_RGB    - an optional boolean variable to handle special RGB data
%       sequence [R1 R2 ... G1 G2 ... B1 B2 ...] that is used only by
%       AnalyzeDirect (Analyze Software). Since both NIfTI and Analyze
%       file format use RGB triple [R1 G1 B1 R2 G2 B2 ...] sequentially
%       for each voxel, this variable is set to FALSE by default. If you
%       would like the saved image only to be opened by AnalyzeDirect
%       Software, set old_RGB to TRUE (or 1). It will be set to 0, if it
%       is default or empty.
%
%  Tip: to change the data type, set nii.hdr.dime.datatype,
%	and nii.hdr.dime.bitpix to:
%
%     0 None                     (Unknown bit per voxel) % DT_NONE, DT_UNKNOWN
%     1 Binary                         (ubit1, bitpix=1) % DT_BINARY
%     2 Unsigned char         (uchar or uint8, bitpix=8) % DT_UINT8, NIFTI_TYPE_UINT8
%     4 Signed short                  (int16, bitpix=16) % DT_INT16, NIFTI_TYPE_INT16
%     8 Signed integer                (int32, bitpix=32) % DT_INT32, NIFTI_TYPE_INT32
%    16 Floating point    (single or float32, bitpix=32) % DT_FLOAT32, NIFTI_TYPE_FLOAT32
%    32 Complex, 2 float32      (Use float32, bitpix=64) % DT_COMPLEX64, NIFTI_TYPE_COMPLEX64
%    64 Double precision  (double or float64, bitpix=64) % DT_FLOAT64, NIFTI_TYPE_FLOAT64
%   128 uint RGB                  (Use uint8, bitpix=24) % DT_RGB24, NIFTI_TYPE_RGB24
%   256 Signed char            (schar or int8, bitpix=8) % DT_INT8, NIFTI_TYPE_INT8
%   511 Single RGB              (Use float32, bitpix=96) % DT_RGB96, NIFTI_TYPE_RGB96
%   512 Unsigned short               (uint16, bitpix=16) % DT_UNINT16, NIFTI_TYPE_UNINT16
%   768 Unsigned integer             (uint32, bitpix=32) % DT_UNINT32, NIFTI_TYPE_UNINT32
%  1024 Signed long long              (int64, bitpix=64) % DT_INT64, NIFTI_TYPE_INT64
%  1280 Unsigned long long           (uint64, bitpix=64) % DT_UINT64, NIFTI_TYPE_UINT64
%  1536 Long double, float128  (Unsupported, bitpix=128) % DT_FLOAT128, NIFTI_TYPE_FLOAT128
%  1792 Complex128, 2 float64  (Use float64, bitpix=128) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128
%  2048 Complex256, 2 float128 (Unsupported, bitpix=256) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128
%
%  Part of this file is copied and modified from:
%  http://www.mathworks.com/matlabcentral/fileexchange/1878-mri-analyze-tools
%
%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%  - "old_RGB" related codes in "save_nii.m" are added by Mike Harms (2006.06.28)
%
function save_nii(nii, fileprefix, old_RGB)

   if ~exist('nii','var') | isempty(nii) | ~isfield(nii,'hdr') | ...
	~isfield(nii,'img') | ~exist('fileprefix','var') | isempty(fileprefix)

      error('Usage: save_nii(nii, filename, [old_RGB])');
   end

   if isfield(nii,'untouch') & nii.untouch == 1
      error('Usage: please use ''save_untouch_nii.m'' for the untouched structure.');
   end

   if ~exist('old_RGB','var') | isempty(old_RGB)
      old_RGB = 0;
   end

   v = version;

   %  Check file extension. If .gz, unpack it into temp folder
   %
   if length(fileprefix) > 2 & strcmp(fileprefix(end-2:end), '.gz')

      if ~strcmp(fileprefix(end-6:end), '.img.gz') & ...
	 ~strcmp(fileprefix(end-6:end), '.hdr.gz') & ...
	 ~strcmp(fileprefix(end-6:end), '.nii.gz')

         error('Please check filename.');
      end

      if str2num(v(1:3)) < 7.1 | ~usejava('jvm')
         error('Please use MATLAB 7.1 (with java) and above, or run gunzip outside MATLAB.');
      else
         gzFile = 1;
         fileprefix = fileprefix(1:end-3);
      end
   end

   filetype = 1;

   %  Note: fileprefix is actually the filename you want to save
   %
   if findstr('.nii',fileprefix) & strcmp(fileprefix(end-3:end), '.nii')
      filetype = 2;
      fileprefix(end-3:end)='';
   end

   if findstr('.hdr',fileprefix) & strcmp(fileprefix(end-3:end), '.hdr')
      fileprefix(end-3:end)='';
   end

   if findstr('.img',fileprefix) & strcmp(fileprefix(end-3:end), '.img')
      fileprefix(end-3:end)='';
   end

   write_nii(nii, filetype, fileprefix, old_RGB);

   %  gzip output file if requested
   %
   if exist('gzFile', 'var')
      if filetype == 1
         gzip([fileprefix, '.img']);
         delete([fileprefix, '.img']);
         gzip([fileprefix, '.hdr']);
         delete([fileprefix, '.hdr']);
      elseif filetype == 2
         gzip([fileprefix, '.nii']);
         delete([fileprefix, '.nii']);
      end;
   end;

   if filetype == 1

      %  So earlier versions of SPM can also open it with correct originator
      %
      M=[[diag(nii.hdr.dime.pixdim(2:4)) -[nii.hdr.hist.originator(1:3).*nii.hdr.dime.pixdim(2:4)]'];[0 0 0 1]];
      save([fileprefix '.mat'], 'M');
   end

   return					% save_nii

%-----------------------------------------------------------------------------------
function write_nii(nii, filetype, fileprefix, old_RGB)

   hdr = nii.hdr;

   if isfield(nii,'ext') & ~isempty(nii.ext)
      ext = nii.ext;
      [ext, esize_total] = verify_nii_ext(ext);
   else
      ext = [];
   end

   switch double(hdr.dime.datatype),
   case   1,
      hdr.dime.bitpix = int16(1 ); precision = 'ubit1';
   case   2,
      hdr.dime.bitpix = int16(8 ); precision = 'uint8';
   case   4,
      hdr.dime.bitpix = int16(16); precision = 'int16';
   case   8,
      hdr.dime.bitpix = int16(32); precision = 'int32';
   case  16,
      hdr.dime.bitpix = int16(32); precision = 'float32';
   case  32,
      hdr.dime.bitpix = int16(64); precision = 'float32';
   case  64,
      hdr.dime.bitpix = int16(64); precision = 'float64';
   case 128,
      hdr.dime.bitpix = int16(24); precision = 'uint8';
   case 256
      hdr.dime.bitpix = int16(8 ); precision = 'int8';
   case 511,
      hdr.dime.bitpix = int16(96); precision = 'float32';
   case 512
      hdr.dime.bitpix = int16(16); precision = 'uint16';
   case 768
      hdr.dime.bitpix = int16(32); precision = 'uint32';
   case 1024
      hdr.dime.bitpix = int16(64); precision = 'int64';
   case 1280
      hdr.dime.bitpix = int16(64); precision = 'uint64';
   case 1792,
      hdr.dime.bitpix = int16(128); precision = 'float64';
   otherwise
      error('This datatype is not supported');
   end

   hdr.dime.glmax = round(double(max(nii.img(:))));
   hdr.dime.glmin = round(double(min(nii.img(:))));

   if filetype == 2
      fid = fopen(sprintf('%s.nii',fileprefix),'w');

      if fid < 0,
         msg = sprintf('Cannot open file %s.nii.',fileprefix);
         error(msg);
      end

      hdr.dime.vox_offset = 352;

      if ~isempty(ext)
         hdr.dime.vox_offset = hdr.dime.vox_offset + esize_total;
      end

      hdr.hist.magic = 'n+1';
      save_nii_hdr(hdr, fid);

      if ~isempty(ext)
         save_nii_ext(ext, fid);
      end
   else
      fid = fopen(sprintf('%s.hdr',fileprefix),'w');

      if fid < 0,
         msg = sprintf('Cannot open file %s.hdr.',fileprefix);
         error(msg);
      end

      hdr.dime.vox_offset = 0;
      hdr.hist.magic = 'ni1';
      save_nii_hdr(hdr, fid);

      if ~isempty(ext)
         save_nii_ext(ext, fid);
      end

      fclose(fid);
      fid = fopen(sprintf('%s.img',fileprefix),'w');
   end

   ScanDim = double(hdr.dime.dim(5));		% t
   SliceDim = double(hdr.dime.dim(4));		% z
   RowDim   = double(hdr.dime.dim(3));		% y
   PixelDim = double(hdr.dime.dim(2));		% x
   SliceSz  = double(hdr.dime.pixdim(4));
   RowSz    = double(hdr.dime.pixdim(3));
   PixelSz  = double(hdr.dime.pixdim(2));

   x = 1:PixelDim;

   if filetype == 2 & isempty(ext)
      skip_bytes = double(hdr.dime.vox_offset) - 348;
   else
      skip_bytes = 0;
   end

   if double(hdr.dime.datatype) == 128

      %  RGB planes are expected to be in the 4th dimension of nii.img
      %
      if(size(nii.img,4)~=3)
         error(['The NII structure does not appear to have 3 RGB color planes in the 4th dimension']);
      end

      if old_RGB
         nii.img = permute(nii.img, [1 2 4 3 5 6 7 8]);
      else
         nii.img = permute(nii.img, [4 1 2 3 5 6 7 8]);
      end
   end

   if double(hdr.dime.datatype) == 511

      %  RGB planes are expected to be in the 4th dimension of nii.img
      %
      if(size(nii.img,4)~=3)
         error(['The NII structure does not appear to have 3 RGB color planes in the 4th dimension']);
      end

      if old_RGB
         nii.img = permute(nii.img, [1 2 4 3 5 6 7 8]);
      else
         nii.img = permute(nii.img, [4 1 2 3 5 6 7 8]);
      end
   end

   %  For complex float32 or complex float64, voxel values
   %  include [real, imag]
   %
   if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
      real_img = real(nii.img(:))';
      nii.img = imag(nii.img(:))';
      nii.img = [real_img; nii.img];
   end

   if skip_bytes
      fwrite(fid, zeros(1,skip_bytes), 'uint8');
   end

   fwrite(fid, nii.img, precision);
%   fwrite(fid, nii.img, precision, skip_bytes);        % error using skip
   fclose(fid);

   return;					% write_nii

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [signalMovie] = createSignalBasedMovie(inputSignals,inputImages,varargin)
	% uses images and signals for sources from an original movie to create a cleaner, more binary movie
	% biafra ahanonu
	% started: 2014.07.20 [14:09:34]
	% inputs
		% inputSignals - [n t], n = number of signals, t = time
		% inputImages - [n x y], n = number of images, x/y are the dimensions of the images, use permute(inputImages,[3 1 2]) if you store z dimension last
	% outputs
		% signalMovie - [x y t] movie, reconstructed from cell traces
	% options
		% filterInputs: should the input images be automatically filtered to remove large or low SNR signals? 0 = no, 1 = yes
		% signalType 'raw' or 'peak', peaks uses a smoothed version of the detected peaks
		% inputPeaks : [n t] matrix (n = number of signals, t = time) of pre-computed peaks, should contain 1 = peak, 0 = no peak

	% changelog
		% 2017.01.14 [20:06:04] - support switched from [nSignals x y] to [x y nSignals]
        % 2019.08.26 [10:45:29] - Additional threshold, peak support
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

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

	%========================
	% should the input images be automatically filtered to remove large or low SNR signals?
	options.filterInputs = 0;
	% 'raw' or 'peak', peaks uses a smoothed version of the detected peaks
	options.signalType = 'peak';
	% [n t] matrix (n = number of signals, t = time) of pre-computed peaks, should contain 1 = peak, 0 = no peak
	options.inputPeaks = [];
	% options for peak
	options.dt = 1/5;
	options.tau = 0.5;
	options.pgam = 1-options.dt/options.tau;
	% 'yes', no'
	options.normalizeOutputMovie = 'yes';
	options.waitbarOn = 1;
    options.imageThreshold = 0.4;
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	try
		if options.filterInputs==1
			[inputImages, inputSignals, valid, imageSizes] = filterImages(inputImages, inputSignals);
		end
		inputSignalsCopy = inputSignals;
		nSignals = size(inputSignalsCopy,1);
		nPts = size(inputSignalsCopy,2);
		mSize = size(inputImages);
		movieDims = [mSize(1) mSize(2) nPts];
		signalMovie = double(NaN(movieDims));
		inputImagesBinary = thresholdImages(inputImages,'binary',1,'threshold',options.imageThreshold);
		if isempty(options.inputPeaks)
			[signalPeaks, signalPeakIdx] = computeSignalPeaks(inputSignals,'makePlots',0,'makeSummaryPlots',0,'numStdsForThresh',2.3,'detectMethod','raw');
		else
			signalPeaks = options.inputPeaks;
		end
		switch options.signalType
			case 'peak'
				inputSignalsCopy = signalPeaks;
				inputSignalsCopy = inputSignalsCopy + 0.5;
				for signalNo = 1:nSignals
					V.dt    = options.dt;  % time step size
					tau     = options.tau;    % decay time constant
			        P.gam   = 1-options.dt/options.tau; % C(t) = gam*C(t-1)
					inputSignalsCopy(signalNo,:) = filter(1,[1 -P.gam],inputSignalsCopy(signalNo,:));         % calcium concentration
				end
			case 'raw'

			otherwise
				% body
		end

		reverseStr = '';
		for signalNo = 1:nSignals
			try
				%thisSignal = zeros([1 1 nPts]);
				%thisSignal(1,1,:) = inputSignalsCopy(signalNo,:);
				thisSignal = inputSignalsCopy(signalNo,:);
				thisImage = squeeze(inputImagesBinary(:,:,signalNo));
				% get location of all cell pixels and setup linear indexing variables
				[x y] = find(thisImage==1);
				xrepmat = repmat(x,[1 nPts])';
				yrepmat = repmat(y,[1 nPts])';
				framerepmat = repmat(1:nPts,[1 length(x)]);
				% get the linear index in the movie of all points to replace, faster than looping over each frame
				linearInd = sub2ind(movieDims, xrepmat(:),yrepmat(:), framerepmat(:));
				% repeat the signal so it will fill the entire binary part of the image
				thisSignalRepmat = repmat(thisSignal,[1 length(x)]);
				%
				% [size(thisSignalRepmat(:)) size(linearInd) 10]
				%isempty(thisSignalRepmat(:))
				% [size(signalMovie)]
				signalMovie(linearInd) = thisSignalRepmat(:);
				reverseStr = cmdWaitbar(signalNo,nSignals,reverseStr,'inputStr','making trace movie','waitbarOn',options.waitbarOn,'displayEvery',1);
				% =======
				% OLD
				%for i=1:nPts
				%    signalMovie(x,y,i) = thisSignal(i);
				%end
				%playMovie(signalMovie);
				%tmpMovie = repmat(thisImage,[1 1 nPts]);
				%tmpMovie = bsxfun(@times,thisSignal,tmpMovie);
				%playMovie(tmpMovie);
				%signalMovie = signalMovie + tmpMovie;
				% =======
			catch err
				display(repmat('@',1,7))
				disp(getReport(err,'extended','hyperlinks','on'));
				display(repmat('@',1,7))
			end
		end

		switch options.normalizeOutputMovie
			case 'yes'
				% normalize movie to keep consistent
				[signalMovie] = normalizeVector(signalMovie,'normRange','zeroToOne');
			case 'no'
				%
			otherwise
				% body
		end

	catch err
		display(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		display(repmat('@',1,7))
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% R8MAT_PRINT_SOME prints out a portion of an R8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, real A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, a title.
%
  incx = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );

  for j2lo = max ( jlo, 1 ): incx : min ( jhi, n )

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col: ' );

    for j = j2lo : j2hi
      fprintf ( 1, '%7d       ', j );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row\n' );

    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, m );

    for i = i2lo : i2hi

      fprintf ( 1, '%7d  ', i );

      for j = j2lo : j2hi
        fprintf ( 1, '%12g  ', a(i,j) );
      end

      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
function bx = bernoulli_poly2 ( n, x )

%*****************************************************************************80
%
%% BERNOULLI_POLY2 evaluates the N-th Bernoulli polynomial at X.
%
%  Discussion:
%
%    Thanks to Bart Vandewoestyne for pointing out an error in the previous
%    documentation, 31 January 2008.
%
%    Special values of the Bernoulli polynomial include:
%
%      B(N,0) = B(N,1) = B(N), the N-th Bernoulli number.
%
%      B'(N,X) = N * B(N-1,X)
%
%      B(N,X+1) - B(N,X) = N * X^(N-1)
%      B(N,X) = (-1)^N * B(N,1-X)
%
%    A formula for the Bernoulli polynomial in terms of the Bernoulli
%    numbers is:
%
%      B(N,X) = sum ( 0 <= K <= N ) B(K) * C(N,K) * X^(N-K)
%
%    The first few polynomials include:
%
%      B(0,X) = 1
%      B(1,X) = X    - 1/2
%      B(2,X) = X^2 -   X      +  1/6
%      B(3,X) = X^3 - 3/2*X^2 +  1/2*X
%      B(4,X) = X^4 - 2*X^3   +      X^2 - 1/30
%      B(5,X) = X^5 - 5/2*X^4 +  5/3*X^3 - 1/6*X
%      B(6,X) = X^6 - 3*X^5   +  5/2*X^4 - 1/2*X^2 + 1/42
%      B(7,X) = X^7 - 7/2*X^6 +  7/2*X^5 - 7/6*X^3 + 1/6*X
%      B(8,X) = X^8 - 4*X^7   + 14/3*X^6 - 7/3*X^4 + 2/3*X^2 - 1/30
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 January 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the Bernoulli polynomial to
%    be evaluated.  N must be 0 or greater.
%
%    Input, real X, the value at which the polynomial is to
%    be evaluated.
%
%    Output, real BX, the value of B(N,X).
%
  fact = 1.0;

  b = bernoulli_number3 ( 0 );

  bx = b;

  for i = 1 : n
    fact = fact * ( n + 1 - i ) / i;
    b = bernoulli_number3 ( i );
    bx = bx * x + fact * b;
  end

  return
end


text:
function pass = test_matrixOutput(~)
% TEST_MATRIXOUTPUT   Chebmatrix operations resulting only in doubles should
% return a normal matrix.

dom = [-3 1];
% A chebmatrix with 5 column functions
T = chebpoly(1:5, dom);
T = chebmatrix(T);
C = T'*T;
pass(1) = isnumeric(C) && all( size(C) == [5 5]);

% Two rows, 5 columns
TT = [T; T];
CC = TT'*TT;
pass(2) = isnumeric(CC) && all( size(CC) == [5 5]);

% Operators involved, output should not be numeric
[Z, I, D, C, M] = linop.primitiveOperators(dom);
A = [I D; Z C];
f = [chebpoly(1, dom); chebpoly(5, dom)];
Af = A*f;
pass(3) = isa(Af, 'chebmatrix') && ...
    all(cellfun(@isa, Af.blocks, {'chebfun'; 'chebfun'})) && ...
    all( size(Af) == [2 1]);

% Functionals involved, output should be numeric
[Z, E, S, D] = linop.primitiveFunctionals(dom);
J = [Z-E(1), S];
Jf = J*f;
pass(4) = isnumeric(Jf) && all(size(Jf) == [1 1]);
K = [J' J'];
Kf = K*f;
pass(5) = isnumeric(Kf) && all(size(Kf) == [2 1]);
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [orbitRow] = getSingularOrbitFromKSPTOTConnect(guid)
    orbitRow = {};

    if(~isempty(guid))
        name = readStringFromKSPTOTConnect('GetVesselNameByGUID', guid, true);
        status = readStringFromKSPTOTConnect('GetVesselStatusByGUID', guid, true);
        orbitDbl = readDoublesFromKSPTOTConnect('GetVesselOrbitByGUID', guid, true);
    else
        name = readStringFromKSPTOTConnect('GetActiveVesselName', '', true);
        status = '';
        orbitDbl = readDoublesFromKSPTOTConnect('GetActiveVesselOrbit', '', true);
    end

    orbitRow{1} = name;
    orbitRow{2} = status;
    orbitRow{3} = orbitDbl(1);
    orbitRow{4} = orbitDbl(2);
    orbitRow{5} = orbitDbl(3);
    orbitRow{6} = orbitDbl(4);
    orbitRow{7} = orbitDbl(5);

    if(orbitDbl(2) > 1.0)
        orbitRow{8} = orbitDbl(8);
    else
        orbitRow{8} = AngleZero2Pi(orbitDbl(8));
    end

    orbitRow{9} = orbitDbl(9);
    orbitRow{10} = orbitDbl(13);

    if(orbitRow{6} < 0)
        orbitRow{6} = orbitRow{6} + 360;
    end

    if(orbitRow{7} < 0)
        orbitRow{7} = orbitRow{7} + 360;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function h = plus(f, g)
%+   Plus for BALLFUN.
% F + G adds F and G. F and G can be scalars or BALLFUNs.

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

fIsBallfun = isa(f, 'ballfun');
gIsBallfun = isa(g, 'ballfun');

% Empty check
if isempty( f )
    h = g;
    if ~gIsBallfun
        h = ballfun(h);
    end
    return
end

% Empty check
if isempty( g )
    h = f;
    if ~fIsBallfun
        h = ballfun(h);
    end
    return
end

if (fIsBallfun && gIsBallfun)
    [mf,nf,pf] = size(f);
    [mg,ng,pg] = size(g);
    m = max(mf,mg);
    n = max(nf,ng);
    p = max(pf,pg);
    X = zeros(m,n,p);
    X(1:mf,floor(n/2)+1-floor(nf/2):floor(n/2)+nf-floor(nf/2),floor(p/2)+1-floor(pf/2):floor(p/2)+pf-floor(pf/2)) = f.coeffs;
    X(1:mg,floor(n/2)+1-floor(ng/2):floor(n/2)+ng-floor(ng/2),floor(p/2)+1-floor(pg/2):floor(p/2)+pg-floor(pg/2)) = ...
    X(1:mg,floor(n/2)+1-floor(ng/2):floor(n/2)+ng-floor(ng/2),floor(p/2)+1-floor(pg/2):floor(p/2)+pg-floor(pg/2)) + g.coeffs;
    h = ballfun(X,'coeffs');
elseif (fIsBallfun && isnumeric(g))
    S = size(f);
    X = f.coeffs;
    % Add the constant g
    X(1,floor(S(2)/2)+1,floor(S(3)/2)+1) = X(1,floor(S(2)/2)+1,floor(S(3)/2)+1) + g;
    h = ballfun(X,'coeffs');
elseif (isnumeric(f) && gIsBallfun)
    S = size(g);
    X = g.coeffs;
    % Add the constant f
    X(1,floor(S(2)/2)+1,floor(S(3)/2)+1) = X(1,floor(S(2)/2)+1,floor(S(3)/2)+1) + f;
    h = ballfun(X,'coeffs');
else
    error('BALLFUN:plus:unknown', ...
          ['Undefined function ''plus'' for input arguments of type ' ...
           '%s and %s.'], class(f), class(g));
end
end


text:
function [B,L,U] = boundingbox(P,Options,x)
%BOUNDINGBOX Computes bounding box of a constraint
%
% If only one output is requested, only the symbolic model is returned
%  B = boundingbox(F)
%
% If three outputs are requested, numerical values are returned too
%  [B,L,U] = boundingbox(F)
%
% A second argument can be used to specificy solver settings
%  B = boundingbox(F,sdpsettings('solver','cplex'))
%
% A third argument can (should!) be used to obtain a bounding box in a
% particular set of variables, i.e. the bounding box of a projection.
% Unless you specify which variables you are computing the bounding box
% w.r.t, it will be very hard for you to understand which variables the
% values in L and U relate to
%  [B,L,U] = boundingbox(F,[],x)
% B will now be the box [L <= x <= U] (infinite bounds not included)

if nargin < 2
    Options = sdpsettings('verbose',0);
end

if nargin == 3
    [B,L,U] = boundingbox(lmi(P.Constraints),Options,x)
else
    [B,L,U] = boundingbox(lmi(P.Constraints),Options)
end
switch nargout
    case 0
        B
    case 1
        varargout{1} = B;
    case 2
        varargout{1} = B;
        varargout{2} = L;
    case 3
        varargout{1} = B;
        varargout{2} = L;
        varargout{3} = U;
    otherwise
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% Tow Thomas Bandpass Filter - Worst-Case Analysis (RSS, EVA, FMCA)
% File: c:\M_files\short_updates\TowThomaswca.m
% updated 11/16/06
clear;clc;
K=1e3;uF=1e-6; % Unit suffixes
% Component values
R1=200*K;R2=10*K;R3=20*K;R4=10*K;R5=20*K;R6=20*K;
C1=0.0159*uF;C2=C1;
%
Nom=[R1 R2 R3 R4 R5 R6 C1 C2]; % Place components in vector X
%
Tr=0.02;Tc=0.1; % 2% resistors, 10% capacitors
%
% Tolerance array T; minus in row 1, plus in row 2
% Column order follows component vector X
%
T=[-Tr -Tr -Tr -Tr -Tr -Tr -Tc -Tc;
   Tr Tr Tr Tr Tr Tr Tc Tc];
%
% Linear frequency; in Hz
%
BF=700;LF=1300;NP=301; % Number of points
F=linspace(BF,LF,NP);
%
Nc=length(Nom); % Nc = number of components
%
% Nominal output
%
Mr=1+(T(2,:)+T(1,:))/2; % For asymmetric tolerances if any
Nav=Nom.*Mr; % Average value of components; = Nom if symmetric tolerances
Tv=(T(2,:)-T(1,:))./(2*Mr); % Tv used for RSS
%
[A,B,D,E,I]=ttbpf(Nav); % call nominal arrays
%
% Arrays are constant and real, and they are calculated only once,
% not at each frequency as in admittance matrix analyses.
% This reduces run times considerably.
%
for i=1:NP;s=2*pi*F(i)*j;Vn(i)=abs(D*((s*I-A)\B)+E);end; % Nominal output
%
Nf=2^Nc;Vm1=zeros(Nf,NP); % Used in FMCA
k=1:Nf;RB=dec2bin(k-1); % Used in FMCA
%
% Sensitivities
%
dpf=0.0001; % derivative perturbation factor
Q=1+dpf;R=1-dpf;
for i=1:NP
   Qx=ones(1,Nc);Rx=Qx; % reset Q & R
   s=2*pi*F(i)*j;
   for p=1:Nc  % start sensitivity loop; p = component counter
      Qx(p)=Q;Rx(p)=R;
      if p > 1;Qx(p-1)=1;Rx(p-1)=1;end % reset previous
      %
      [A,B,D,E,I]=ttbpf(Nom.*Qx); % arrays perturbated forward
      Vr=abs(D*inv(s*I-A)*B+E); % output perturbated forward
      %
      [A,B,D,E,I]=ttbpf(Nom.*Rx); % arrays perturbated backward
      Vb=abs(D*inv(s*I-A)*B+E); % output perturbated backward
      %
      Sens(i,p)=(Vr-Vb)/(2*Vn(i)*dpf); % centered difference approximation
 % Ref:  Numerical Methods for Engineers, 3rd ed.
 % S.C. Chapra & R.P. Canale, McGraw-Hill, 1998, p.93
 %
 % EVA (Extreme Value Analysis)
 %
      if Sens(i,p) > 0
         Hi(p)=1+T(2,p);Lo(p)=1+T(1,p); % component tolerance for WC High
      else
         Hi(p)=1+T(1,p);Lo(p)=1+T(2,p); % component tolerance for WC Low
      end
   end % close p sensitivity loop
%
% Get extreme value low (EVL) and extreme value high (EVH)
%
	[A,B,D,E,I]=ttbpf(Nav.*Lo);
   VL(i)=abs(D*((s*I-A)\B)+E);
   %
	[A,B,D,E,I]=ttbpf(Nav.*Hi);
   VH(i)=abs(D*((s*I-A)\B)+E);
%
% RSS (Root-Sum-Square)
%
   STn=norm(Sens(i,:).*Tv);
   Vrss1(i)=Vn(i)*(1-STn);Vrss2(i)=Vn(i)*(1+STn);
%
% FMCA
%
% Binary array RB contains binary numbers Nc bits wide from 0 to Nf.
% Used for all possible tolerance combinations.
%
   for k=1:Nf  % start FMCA loop
      for p=1:Nc
         if RB(k,p)=='0'
            Tf(p)=1+T(1,p); % low tolerance if a logic '0'
         else
            Tf(p)=1+T(2,p); % high tolerance if a logic '1'
         end
      end % close FMCA tolerance loop
      [A,B,D,E,I]=ttbpf(Nav.*Tf);
      Vm1(k,i)=abs(D*((s*I-A)\B)+E);
   end % close FMCA loop
end % close major i loop
Vmax1=max(Vm1);Vmin1=min(Vm1);
%
% Plot sensitivities
% Number of sensitivities to be plotted = Nc
% Adjust 'g=plot(...' statements accordingly.
% All should be plotted to observe any bipolarity (non-monotonicity).
h=plot(F,Sens(:,1),'k',F,Sens(:,2),'r',F,Sens(:,3),'g',F,Sens(:,4),'b');
set(h,'LineWidth',2);grid on
axis auto
ylabel('%/%');xlabel('Freq (Hz)');
title('Sensitivities')
hold off
legend('R1','R2','R3','R4',0);
figure
%
h=plot(F,Sens(:,5),'k',F,Sens(:,6),'r',F,Sens(:,7),'g',F,Sens(:,8),'b');
set(h,'LineWidth',2);
axis auto;grid on
ylabel('%/%');xlabel('Freq (Hz)');
title('Sensitivities')
legend('R5','R6','C1','C2',0);
figure
%
% RSS
h=plot(F,Vn,'k--',F,Vrss1,'b',F,Vrss2,'r');
set(h,'LineWidth',2);grid on
axis auto
ylabel('Volts');xlabel('Freq (Hz)');
title('RSS')
legend('Nom','RSS Lo','RSS Hi',0);
figure
%
% EVA
h=plot(F,Vn,'k--',F,VH,'r',F,VL,'b');
set(h,'LineWidth',2);grid on
axis auto
xlabel('Freq(Hz)');ylabel('Volts');
title('EVA')
legend('Nom','EVA Hi','EVA Lo',0);
figure
%
% FMCA
h=plot(F,Vn,'k--',F,Vmax1,'r',F,Vmin1,'b');
set(h,'LineWidth',2);grid on
axis auto
xlabel('Freq(Hz)');ylabel('Volts');
title('FMCA');
legend('Nom','FMCA Hi','FMCA Lo',0);
%
% See comments concerning RSS, EVA, & FMCA at end of mfbpfwca.m


text:
function point = intersectPlaneLine(plane, line)
%INTERSECTPLANELINE return intersection between a plane and a line
%
%   PT = intersectPlaneLine(PLANE, LINE) return the intersection point of
%   the given line and the given plane.
%   PLANE : [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2]
%   LINE :  [x0 y0 z0 dx dy dz]
%   PT :    [xi yi zi]
%
%   Note: deprecated. Replaced by function 'intersectLinePlane'
%
%  Songbai Ji (6/23/2006). Bug fixed; also allow one plane, many lines;
% many planes one line; or N planes and N lines configuration in the input.
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 17/02/2005.
%

%   HISTORY
%   24/11/2005 add support for multiple input
%   23/06/2006 correction from Songbai Ji
%   14/12/2006 correction for parallel lines and plane normals
%   05/01/2007 fixup for parallel lines and plane normals
%   17/10/2008 add warning for deprecation

warning('IMAEL:deprecatedFunction', ...
    'This function is deprecated, use ''intersectLinePlane'' instead');

% unify sizes of data
if size(plane, 1) == 1;     % one plane possible many lines
    plane = repmat(plane, size(line,1), 1);
elseif size(line,1) == 1;   % one line and many planes
    line = repmat(line, size(plane,1), 1);
elseif (size(plane,1) ~= size(line,1)) ; % N planes and M lines, not allowed for now.
    error('input size not correct, either one/many plane and many/one line, or same # of planes and lines!');
end

% initialize empty array
point = zeros(size(plane, 1), 3);

% plane normal
n = cross(plane(:,4:6), plane(:, 7:9), 2);

% get indices of line and plane which are parallel
par = abs(dot(n, line(:,4:6), 2))<1e-14;
point(par,:) = NaN;
% old version:
% II = find(abs(dot(n, line(:,4:6), 2))<1e-14);
% if ~isempty(II)
%     point(II,:) = [NaN NaN NaN];
% end

% difference between origins of plane and line
dp = plane(:,1:3) - line(:,1:3);

% relative position of intersection on line
% Should be Array multiply, original file had a bug. (songbai ji
% 6/23/2006).
% Divide only for non parallel vectors (DL,
t = dot(n(~par,:), dp(~par,:), 2)./dot(n(~par,:), line(~par,4:6), 2);
%t = dot(n, dp, 2)./dot(n, line(:,4:6), 2);

% compute coord of intersection point
% point = line(:,1:3) + t*line(:,4:6);
% NOTE: original m file had a bug (in the above line).  It should be
% corrected as follows.  (Songbai Ji 6/23/2006).
%point = line(~par,1:3) + repmat(t,1,3).*line(~par,4:6);
%point = line(:,1:3) + repmat(t,1,3).*line(:,4:6);
point(~par, :) = line(~par,1:3) + repmat(t,1,3).*line(~par,4:6);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function output = SSIM(Reference, Clean, Target)

sim_sum=0;

[m1,n1] = size(Reference);
[m2,n2] = size(Target);
[m3,n3] = size(Clean);

r=[m1;m2;m3;];
c=[n1/3;n2/3;n3/3;];

rows=min(r);
cols=min(c);
for i=1:3

    x1=(Reference(1:rows,1:cols,i));
    y=(Target(1:rows,1:cols,i));
    x2=(Clean(1:rows,1:cols,i));

    avg = (x1 + x2)/2;

    [ssimval, ~] = ssim(y,avg);

    sim_sum = sim_sum + ssimval;

end

output = sim_sum/3;

author: lbasek

repo: image-denoising-benchmark


text:
%normu Normalize image points to be used for LS estimation.
%   A = normu(u) finds normalization matrix A so that mean(A*u)=0 and
%   mean(||A*u||)=sqrt(2). (see Hartley: In Defence of 8-point Algorithm,
%   ICCV`95).
%
%   Parameters:
%      u ... Size (2,N) or (3,N). Points to normalize.
%      A ... Size (3,3). Normalization matrix.

function A = normu(u)

if size(u,1)==3, u = p2e(u); end

m        = mean(u')'; % <=> mean (u,2)
u        = u - m*ones(1,size(u,2));
distu    = sqrt(sum(u.*u));
r        = mean(distu)/sqrt(2);
if ~r, A = []; return; end     % one of degenarate configurations
A        = diag([1/r 1/r 1]);
A(1:2,3) = -m/r;

author: strawlab

repo: MultiCamSelfCal


text:
function DEM_demo_convolution
% DEM demo for linear deconvolution:  This demo considers the deconvolution
% of the responses of a single-input-multiple output input-state-output
% model (DCM) to disclose the input or causes.  It focuses on estimating the
% causes and hidden states: The notes provide a comparative evaluation with
% extended Kalman filtering (see script after return).

% get a simple convolution model
%==========================================================================
M       = spm_DEM_M('convolution model');
M(1).V  = exp(8);                             % error precision
M(1).W  = exp(8);                             % error precision

% and generate data
%==========================================================================
N       = 32;                                 % length of data sequence
U       = exp(-([1:N] - 12).^2/(2.^2));       % this is the Gaussian cause;
DEM     = spm_DEM_generate(M,U,{},{[] 16},{16});

% invert model
%==========================================================================
DEM     = spm_DEM(DEM);

% overlay true values
%--------------------------------------------------------------------------
spm_DEM_qU(DEM.qU,DEM.pU)

return

% explore dimensions (n)
%==========================================================================
for i = 1:8
    DEM.M(1).E.n = i;
    DEM.M(1).E.d = 3;

    D{i}   = spm_DEM(DEM);
    F(i)   = D{i}.F;
    Sx(i)  = sum(sum((D{i}.qU.x{1} - DEM.pU.x{1}).^2));
    Sv(i)  = sum(sum((D{i}.qU.v{2} - DEM.pU.v{2}).^2));
end

% plot
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 1');

subplot(2,1,1)
bar(Sv)
xlabel('n - 1 (d = 2)')
ylabel('sum squared error (causal states)')
title('embedding dimension (n)','FontSize',16)
axis square

d     = [1 8];
for i = 1:length(d)
    subplot(2,length(d),i + length(d))
    hold on
    plot([1:N],D{d(i)}.qU.v{2})
    plot([1:N],DEM.pU.v{2},'linewidth',2,'color',[1 1 1]/2)
    hold off
    axis square
    xlabel('time')
    title(sprintf('n = %i',d(i)),'FontSize',16)
    if i == 1, a = axis; else, axis(a); end
end

% and d
%--------------------------------------------------------------------------
clear D F Sx Sv
for i = 1:6
    DEM.M(1).E.n = 15;
    DEM.M(1).E.d = i - 1;

    D{i}   = spm_DEM(DEM);
    F(i)   = D{i}.F;
    Sx(i)  = sum(sum((D{i}.qU.x{1} - DEM.pU.x{1}).^2));
    Sv(i)  = sum(sum((D{i}.qU.v{2} - DEM.pU.v{2}).^2));

end

% plot
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 2'); clf

subplot(2,2,1)
bar(Sx)
xlabel('d (n = 15)')
ylabel('sum squared error (hidden states)')
title('embedding dimension (d)','FontSize',16)
axis square
set(gca,'XTickLabel',[0:7])

subplot(2,2,2)
bar(F - min(F) + 32)
xlabel('d (n = 15)')
ylabel('log-evidence')
title('log-evidence','FontSize',16)
axis square
set(gca,'XTickLabel',[0:7])

d     = [1 4];
for i = 1:length(d)
    subplot(2,length(d),i + length(d))
    hold on
    plot([1:N],D{d(i)}.qU.x{1})
    plot([1:N],DEM.pU.x{1},'linewidth',2,'color',[1 1 1]/2)
    hold off
    axis square
    xlabel('time')
    title(sprintf('d = %i',d(i) - 1),'FontSize',16)
    if i == 1, a = axis; else, axis(a); end
end

% Comparison with EKF
%==========================================================================
clear SSE
for i = 1:8

    % i.i.d.
    %----------------------------------------------------------------------
    DEM      = spm_DEM_generate(M,U,{},{[] 16});
    DEM      = spm_DEM(DEM);

    % serial correlations
    %----------------------------------------------------------------------
    D0          = DEM;
    D0.M(1).E.s = 0;
    D0          = spm_DEM(D0);

    % extended kalman filter - i.i.d.
    %----------------------------------------------------------------------
    e_x      = spm_ekf(DEM.M,DEM.Y);
    d_x      = DEM.qU.x{1};
    t_x      = DEM.pU.x{1};
    i_x      = D0.qU.x{1};

    SSE(i,1) = sum(sum((e_x - t_x).^2));
    SSE(i,2) = sum(sum((i_x - t_x).^2));
    SSE(i,3) = sum(sum((d_x - t_x).^2));

end

spm_figure('GetWin','Figure 3');

subplot(2,1,1)
plot(1:N,i_x,'g',1:N,d_x,'r',1:N,e_x,'b',1:N,t_x,'k')
legend('DEM(0)',' ','DEM(0.5)',' ','EKF',' ','true',' ')
xlabel('time')
title({'hidden states','comparison with EKF'},'FontSize',16)
axis square

subplot(2,1,2)
plot(SSE','k:'), hold on
plot(SSE','k.','Markersize',16), hold off
set(gca,'Xtick',[1 2 3],'XLim',[0 4])
set(gca,'Xticklabel',{'EKF','DEM(0)','DEM(0.5)'})
title('sum of squared error (hidden states)','FontSize',16)
axis square

% Show equivalence when causes are not structured
%==========================================================================
M(1).E.s  = 0;
M(1).pE.h = eye(2);
M(1).pC   = [];

M(2).v    = [0;0];
M(2).V    = eye(2);
U         = [U; 0*U];

DEM       = spm_DEM_generate(M,U,{},{[] 16});
DEM       = spm_DEM(DEM);
e_x       = spm_ekf(DEM.M,DEM.Y);
d_x       = DEM.qU.x{1};

spm_figure('GetWin','Figure 4');

subplot(2,1,1)
plot(1:N,d_x,'r',1:N,e_x,'b')
legend('DEM(0)',' ','EKF',' ')
xlabel('time')
title({'hidden states';'no correlations'},'FontSize',16)
axis square

% Show equivalence when causes are removed
%--------------------------------------------------------------------------
DEM.M(1).pE.h = sparse(2,2);
DEM.M(1).W    = speye(2);

DEM           = spm_DEM(DEM);
e_x           = spm_ekf(DEM.M,DEM.Y);
d_x           = DEM.qU.x{1};

spm_figure('GetWin','Figure 4');

subplot(2,1,2)
plot(1:N,d_x,'r',1:N,e_x,'b')
legend('DEM(0)',' ','EKF',' ')
xlabel('time')
title({'hidden states';'no causes'},'FontSize',16)
axis square


text:
function y = vl_nnchannelshuffle(x, group, varargin)
% VL_NNCHANNELSHUFFLE Channel shuffling
%   Y = VL_NNRESHAPE(X, SHAPE) reshpaes the input data X to have
%   the dimensions specified by SHAPE. X is a SINGLE array of
%   dimension H x W x D x N where (H,W) are the height and width of
%   the map stack, D is the image depth (number of feature channels)
%   and N the number of of images in the stack. SHAPE is a 1 x 3 cell
%   array, the contents of which are passed in order to the MATLAB
%   reshape function. As a consequence, `[]` an be used to specify a
%   dimension which should be computed from the other two. The batch size
%   (the fourth dimension of the input) is left unchanged by this
%   reshaping operation.
%
%   Example:
%       Inputs: X with shape [100 100 3 5] and SHAPE = { 100 3 [] }
%       will produce an output Y with shape [100 3 100 5]
%
%   DZDX = VL_NNRESHAPE(X, SHAPE, DZDY) computes the derivatives of the
%   block projected onto DZDY. DZDX and DZDY have the same dimensions
%   as X and Y respectively.
%
%  This operation was originally described in:
%
%  Zhang, X., Zhou, X., Lin, M., & Sun, J. (2017).
%  ShuffleNet: An Extremely Efficient Convolutional Neural Network for
%  Mobile Devices. arXiv preprint arXiv:1707.01083.
%
% Copyright (C) 2017 Samuel Albanie
% Licensed under The MIT License [see LICENSE.md for details]

  [~, dzdy] = vl_argparsepos(struct(), varargin) ;
  keyboard

  batchSize = size(x, 4) ;

  if isempty(dzdy)
    y = reshape(x, shape{1}, shape{2}, shape{3}, batchSize) ;
  else
    y = reshape(dzdy{1}, size(x)) ;
  end

repo: MFT


text:
% PAD_SIZE Compute the optimal size for padding
%
% Usage
%   sz_padded = PAD_SIZE(sz, min_margin, max_ds)
%
% Input
%   sz (int): The size of the original signal.
%   min_margin (int): The minimum margin for padding.
%   max_ds (int): The maximum downsampling factor.
%
% Output
%   sz_padded (int): The minimum size of the padded signal.
%
% Description
%   Calculates the smallest multiple of 2^max_ds larger than sz by at least
%   2*min_margin. This ensures that there is enough margin on both sides of
%   the signal to avoid border effects, assuming that min_margin is equal
%   to at least half of the size of the largest filter used, while ensuring
%   that downsampling by powers of 2 up to 2^max_ds are possible through
%   periodization of the Fourier transform.
%   sz_added is also enforced to be at least 1
%
% See Also
%   PAD_SIGNAL, UNPAD_SIGNAL

function sz_padded = pad_size(sz, min_margin, max_ds)
    sz_padded = 2^max_ds * ceil( (sz + 2*min_margin)/2^max_ds );
    sz_padded = max(1, sz_padded);
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
  function [fwhm_best, costs, im_best] = ...
	fwhm_match(true_image, blurred_image, fwhms)
%|function [fwhm_best, costs, im_best] = ...
%|	fwhm_match(true_image, blurred_image, fwhms)
%|
%| given a blurred_image of a true_image, find the FHWM of a Gaussian kernel
%| that, when convolved to the true_image, yields the smoothed image
%| that best matches blurred_image.
%|
%| the set of FWHM values given in the array fwhms is tried.
%|
%| Copyright 2001-8-30, Jeff Fessler, University of Michigan

if nargin == 1 && streq(true_image, 'test'), fwhm_match_test, return, end
if nargin < 2, ir_usage, end

if nargin < 3
	fwhms = 0:0.5:4;
end

costs = zeros(size(fwhms));
cost_min = Inf;
for ii=1:length(fwhms)
	fwhm = fwhms(ii);
	kern = gaussian_kernel(fwhm);
	psf = kern * kern';
	tmp = conv2(true_image, psf, 'same');
	costs(ii) = norm(tmp(:) - blurred_image(:)) / norm(true_image(:));
	if costs(ii) < cost_min
		im_best = tmp;
	end
end

[dummy ibest] = min(costs);
if ibest == 1 || ibest == length(fwhms)
	warning 'need wider range of fwhms'
end
fwhm_best = fwhms(ibest);

% fwhm_match_test
function fwhm_match_test

% pyramidal PSF to stress the approach
psf1 = [0:5 4:-1:0]; psf1 = psf1 / sum(psf1); psf = psf1' * psf1;
true_image = zeros(128); true_image(64:96,64:96) = 1;
blurred_image = conv2(true_image, psf, 'same');
im plc 2 2
im(1, true_image, 'True Image')
im(2, blurred_image, 'Blurred Image')

fwhms = [2:0.25:8];
[fwhm_best, costs] = fwhm_match(true_image, blurred_image, fwhms);
np = length(psf);	ip = -(np-1)/2:(np-1)/2;
kern = gaussian_kernel(fwhm_best);
nk = length(kern);	ik = -(nk-1)/2:(nk-1)/2;

if im
	im subplot 3
	plot(fwhms, costs, 'c-o', fwhm_best, min(costs), 'yx')
	xlabel FWHM, ylabel Cost, title 'Cost vs FWHM'
	im subplot 4
	plot(ip, psf1, '-o', ik, kern(:), '-+')
	xlabel pixel, title 'PSF profile: actual and Gaussian fit'
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt