text:
function directories = XomicsToMultipleModels(modelGenerationConditions, param)
% Variations of the xomicstomodel function can be generated using this
% function.
%
% USAGE:
%
%    directories = XomicsToMultipleModels(modelGenerationConditions, param)
%
% INPUT:
%    modelGenerationConditions: Options to vary or to save the data
%
%       * .activeGenesApproach -The different approached to identify the active
%          genes (Possible options: 'deleteModelGenes' and 'oneRxnPerActiveGene';
%          default: 'deleteModelGenes');
%       * .boundsToRelaxExoMet - The type of bounds that can be relaxed, upper bounds,
%          lower bounds or both ('b'; possible options: 'u', 'l' and 'b';
%          default: 'b');
%       * .closeIons - Indicate whether the ions are open or closed (Possible options:
%          true and false; default: false);
%       * .cobraSolver - Optimisation solvers supported by the function. Possible
%          options: 'glpk', 'gurobi', 'ibm_cplex', 'matlab'; default: 'gurobi';
%       * .curationOverOmics - indicates whether curated data should take priority
%          over omics data ; default: false;
%       *. genericModel: Generic COBRA model(s)
%       * .inactiveGenesTranscriptomics - Use inactive transcriptomic genes or not
%          (Possible options: true and false; default: false);
%       * .specificData - Specific data variations (Default: empty)
%       * .limitBounds - Boundary on the model (Default: 1000).
%       * .metabolomicsBeforeExtraction - Indicate whether the metabolomic
%          data is included before or after the extraction (Possible options:
%          true and false; default: true);
%       * .tissueSpecificSolver - Extraction solver (Possible options: 'fastCore' and
%          'thermoKernel'; default: 'thermoKernel')
%       * .outputDir - Directory where the models will be generated (Default: current
%          directory)
%       * .transcriptomicThreshold - Transcriptomic thresholds that are defined by the
%          user (Default: log2(2));
%
%   param:           Variable with fixed parameters (Default: empty struct array)
%
% OUTPUTS:
%
%	directories - Array with the name of the new directories

if nargin < 2 || isempty(param)
    param = struct;
end

% Directory to save the models
if isfield(modelGenerationConditions, 'outputDir')
    outputDir = modelGenerationConditions.outputDir;
    outputDir = [regexprep(char(outputDir),'(/|\\)$',''), filesep];
else
    outputDir = [pwd filesep];
end

directories = [];

%% Set conditions

% Cobra solver
if isfield(modelGenerationConditions, 'cobraSolver')
    cobraSolver = modelGenerationConditions.cobraSolver;
else
    cobraSolver = {'gurobi'};
end
% genericModel
if isfield(modelGenerationConditions, 'genericModel')
    models = modelGenerationConditions.genericModel;
    modelsLabels = fieldnames(models);
else
    error('A generic model is needed in modelGenerationConditions.genericModel')
end
% Context-specific input data
if isfield(modelGenerationConditions, 'specificData')
    specificDataforXomics = modelGenerationConditions.specificData;
    specificDataLabels = fieldnames(specificDataforXomics);
else
    specificDataforXomics = struct;
end
% Tissue specific solver
if isfield(modelGenerationConditions, 'tissueSpecificSolver')
    tissueSpecificSolver = modelGenerationConditions.tissueSpecificSolver;
elseif ~isfield(modelGenerationConditions, 'tissueSpecificSolver') && isfield(param, 'tissueSpecificSolver')
    tissueSpecificSolver = {param.tissueSpecificSolver};
else
    tissueSpecificSolver = {'thermoKernel'};
end
% Active genes approach
if isfield(modelGenerationConditions, 'activeGenesApproach')
    activeGenesApproach = modelGenerationConditions.activeGenesApproach;
elseif ~isfield(modelGenerationConditions, 'activeGenesApproach') && isfield(param, 'activeGenesApproach')
    activeGenesApproach = {param.activeGenesApproach};
else
    activeGenesApproach = {'oneRxnPerActiveGene'};
end
% Transcriptomic threshold
if isfield(modelGenerationConditions, 'transcriptomicThreshold')
    transcriptomicThreshold = modelGenerationConditions.transcriptomicThreshold;
elseif ~isfield(modelGenerationConditions, 'transcriptomicThreshold') && isfield(param, 'transcriptomicThreshold')
    transcriptomicThreshold = param.transcriptomicThreshold;
else
    transcriptomicThreshold = 2;
end
% Limit bounds
if isfield(modelGenerationConditions, 'limitBounds')
    limitBounds = modelGenerationConditions.limitBounds;
elseif ~isfield(modelGenerationConditions, 'limitBounds') && isfield(param, 'limitBounds')
    limitBounds = param.limitBounds;
else
    limitBounds = 10e4;
end
% Use inactive genes from transcriptomics
if isfield(modelGenerationConditions, 'inactiveGenesTranscriptomics')
    inactiveGenesTranscriptomics = modelGenerationConditions.inactiveGenesTranscriptomics;
elseif ~isfield(modelGenerationConditions, 'inactiveGenesTranscriptomics') && isfield(param, 'inactiveGenesTranscriptomics')
    inactiveGenesTranscriptomics = param.inactiveGenesTranscriptomics;
else
    inactiveGenesTranscriptomics = false;
end
% Ions exchange
if isfield(modelGenerationConditions, 'closeIons')
    closeIons = modelGenerationConditions.closeIons;
elseif ~isfield(modelGenerationConditions, 'closeIons') && isfield(param, 'closeIons')
    closeIons = param.closeIons;
else
    closeIons = false;
end
% boundsToRelaxExoMet
if isfield(modelGenerationConditions, 'boundsToRelaxExoMet')
    boundsToRelaxExoMet = modelGenerationConditions.boundsToRelaxExoMet;
else
    boundsToRelaxExoMet = {'b'};
end
% curationOverOmics
if isfield(modelGenerationConditions, 'curationOverOmics')
    curationOverOmics = modelGenerationConditions.curationOverOmics;
elseif ~isfield(modelGenerationConditions, 'curationOverOmics') && isfield(param, 'curationOverOmics')
    curationOverOmics = param.curationOverOmics;
else
    curationOverOmics = false;
end
% metabolomicsBeforeExtraction
if isfield(modelGenerationConditions, 'metabolomicsBeforeExtraction')
    metabolomicsBeforeExtraction = modelGenerationConditions.metabolomicsBeforeExtraction;
elseif ~isfield(modelGenerationConditions, 'metabolomicsBeforeExtraction') && ...
        isfield(param, 'metabolomicsBeforeExtraction')
    metabolomicsBeforeExtraction = param.metabolomicsBeforeExtraction;
else
    metabolomicsBeforeExtraction = true;
end

disp('Generating models ...');

% Data with more than one variation will be included in the folder name
conditionsBool = false(11, 1);
if length(cobraSolver) > 1
    conditionsBool(1) = true;
end
if length(modelsLabels) > 1
    conditionsBool(2) = true;
end
if length(specificDataLabels) > 1
    conditionsBool(3) = true;
end
if length(tissueSpecificSolver) > 1
    conditionsBool(4) = true;
end
if length(activeGenesApproach) > 1
    conditionsBool(5) = true;
end
if length(transcriptomicThreshold) > 1
    conditionsBool(6) = true;
end
if length(limitBounds)> 1
    conditionsBool(7) = true;
end
if length(inactiveGenesTranscriptomics) > 1
    conditionsBool(8) = true;
end
if length(closeIons) > 1
    conditionsBool(9) = true;
end
if length(boundsToRelaxExoMet) > 1
    conditionsBool(10) = true;
end
if length(curationOverOmics) > 1
    conditionsBool(11) = true;
end
if length(metabolomicsBeforeExtraction) > 1
    conditionsBool(12) = true;
end

%% Check if there are any fields in modelGenerationConditions that are not specified above
% Because these need to be added to param and cannot yet be used here.
unknownmodelGenerationConditions = setdiff(fieldnames(modelGenerationConditions), ...
    {'cobraSolver'
    'genericModel'
    'specificData'
    'tissueSpecificSolver'
    'activeGenesApproach'
    'transcriptomicThreshold'
    'limitBounds'
    'inactiveGenesTranscriptomics'
    'closeIons'
    'boundsToRelaxExoMet'
    'curationOverOmics'
    'metabolomicsBeforeExtraction'
    'outputDir'});

if ~isempty(unknownmodelGenerationConditions)
    msg = ['Unknown parameter(s) specified in modelGenerationConditions: ' strjoin(unknownmodelGenerationConditions) '. Specify these in param'];
    error(msg)
end

%% Generate models

% Select cobraSolver
dirIdx = 0;
for firstGruoup = 1:length(cobraSolver)
    switch cobraSolver{firstGruoup}
        case 'gurobi'
            solverLabel = 'gurobi';
            [solverOK, solverInstalled] = changeCobraSolver('gurobi','all');
        case 'ibm_cplex'
            solverLabel = 'ibmCplex';
            [solverOK, solverInstalled] = changeCobraSolver('ibm_cplex','all');
        case 'mosek'
            solverLabel = 'mosek';
            [solverOK, solverInstalled] = changeCobraSolver('mosek','all');
    end

    % Input model
    for secondGroup = 1:length(modelsLabels)
        modelLabel = modelsLabels{secondGroup};
        genericModel = models.(modelsLabels{secondGroup});

        % specificDataInfo
        for thirdGruoup = 1:length(specificDataLabels)
            specificDataLabel = specificDataLabels{thirdGruoup};
            specificData = specificDataforXomics.(specificDataLabels{thirdGruoup});

            % Extraction algorithm
            for fourthGroup = 1:length(tissueSpecificSolver)
                param.tissueSpecificSolver = tissueSpecificSolver{fourthGroup};

                % Active genes approach
                for fifthGroup = 1:length(activeGenesApproach)
                    param.activeGenesApproach = activeGenesApproach{fifthGroup};
                    param.activeGenesApproach = regexprep(activeGenesApproach{fifthGroup}, 'oneRxnsPerActiveGene', 'oneRxnPerActiveGene');

                    % Transcriptomic Threshold
                    for sixthGroup = 1:length(transcriptomicThreshold)
                        param.transcriptomicThreshold = transcriptomicThreshold(sixthGroup);

                        % Load input data
                        for seventhGroup = 1:length(limitBounds)
                            param.TolMinBoundary = -limitBounds(seventhGroup);
                            param.TolMaxBoundary = limitBounds(seventhGroup);

                            % Option to identify genes as inactive if
                            % expressed below transcriptomic threshold
                            for eigthGroup = 1:length(inactiveGenesTranscriptomics)
                                param.inactiveGenesTranscriptomics = inactiveGenesTranscriptomics(eigthGroup);
                                if inactiveGenesTranscriptomics(eigthGroup)
                                    inactiveGenesTLabel = 'inactiveGenesT';
                                else
                                    inactiveGenesTLabel = 'NoInactiveGenesT';
                                end

                                % closeIons
                                for ninethGroup = 1:length(closeIons)
                                    param.closeIons = closeIons(ninethGroup);
                                    if closeIons(ninethGroup)
                                        ionsLabel = 'closedIons';
                                    else
                                        ionsLabel = 'openIons';
                                    end

                                    % Bounds to relax metabolomicsToModel (TODO)
                                    for tenthGroup = 1:length(boundsToRelaxExoMet)
                                        param.relaxOptions.bounds = boundsToRelaxExoMet{tenthGroup};
                                        switch boundsToRelaxExoMet{tenthGroup}
                                            case 'upper'
                                                boundsLabel = 'mediaUBRelaxed';
                                            case 'lower'
                                                boundsLabel = 'mediaLBRelaxed';
                                            case 'both'
                                                boundsLabel = 'mediaBBRelaxed';
                                        end

                                        for eleventhGroup = 1:length(curationOverOmics)
                                            param.curationOverOmics = curationOverOmics(eleventhGroup);
                                            if curationOverOmics(eleventhGroup)
                                                priorityLabel = 'curationOverOmics';
                                            else
                                                priorityLabel = 'omicsOverCuration';
                                            end

                                            for twelfthGroup = 1:length(metabolomicsBeforeExtraction)
                                                param.metabolomicsBeforeExtraction = metabolomicsBeforeExtraction(twelfthGroup);
                                                if metabolomicsBeforeExtraction(twelfthGroup)
                                                    exoMetLabel = 'exoMetBeforeExtraction';
                                                else
                                                    exoMetLabel = 'exoMetAfterExtraction';
                                                end

                                                % Establish the working directory
                                                conditions = [solverLabel; ...
                                                    modelLabel; ...
                                                    specificDataLabel; ...
                                                    tissueSpecificSolver{fourthGroup}; ...
                                                    activeGenesApproach{fifthGroup}; ...
                                                    {['transcriptomicsT' num2str(transcriptomicThreshold(sixthGroup))]}; ...
                                                    {['limitBoundary.' num2str(limitBounds(seventhGroup))]}; ...
                                                    inactiveGenesTLabel; ...
                                                    ionsLabel; ...
                                                    boundsLabel;
                                                    priorityLabel;
                                                    exoMetLabel];

                                                workingDirectory = [outputDir ...
                                                    strjoin(conditions(conditionsBool), '_')];
                                                param.workingDirectory = workingDirectory;

                                                dirContent = dir(workingDirectory);
                                                if ~isfolder(workingDirectory)
                                                    mkdir(workingDirectory)
%                                                 elseif any(~cellfun(@isempty, regexp({dirContent.name}, '.txt')))
%                                                     continue
                                                end

                                                % Start the diary
                                                if isunix()
                                                    name = getenv('USER');
                                                else
                                                    name = getenv('username');
                                                end
                                                param.diaryFilename = [workingDirectory filesep datestr(now,30) '_' name '_diary.txt'];

                                                % Model generation only if no model already present in folder
                                                fileNames = what(workingDirectory);

                                                if any(ismember(fileNames.mat,'Model.mat'))
                                                    load([workingDirectory filesep 'Model.mat'])
                                                    %test properties of previously computed model to decide if recomputation is necessary
                                                    if size(Model.S,2)==length(Model.expressionRxns)
                                                        recompute = 0;
                                                    else
                                                        fprintf('%s\n', ['Recomputing because prexisting inconsistent model in ' workingDirectory ])
                                                        recompute = 1;
                                                    end
                                                else
                                                    recompute = 1;
                                                end

                                                if ~recompute && any(ismember(fileNames.mat,'Model.mat'))
                                                    fprintf('%s\n', ['Prexisting model in ' workingDirectory ])
                                                else
                                                    fprintf('%s\n', ['Computing new model in ' workingDirectory ])
                                                    try
                                                        [omicsModel, modelGenerationReport] = XomicsToModel(genericModel, specificData, param);
                                                        % Save the model with the correct name
                                                        Model = omicsModel;
                                                        save([workingDirectory filesep 'Model.mat'], 'Model', 'modelGenerationReport')
                                                    catch ME
                                                        warning('XomicsToModel failed to run')
                                                        disp(param)
                                                        disp(ME)
                                                        msgText = getReport(ME)
                                                        % Close the diary if the run crashed
                                                        fprintf('%s\n', ['Diary written to: ' param.diaryFilename])
                                                        diary off
                                                    end
                                                end

                                                if any(conditionsBool)
                                                    dirIdx = dirIdx + 1;
                                                    directories{dirIdx} = strjoin(conditions(conditionsBool), '_');
                                                else
                                                    directories = pwd;
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
classdef MimMarkerPoint
    %MimMarkerPoint Used in exporting of marker points to disk
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %
    properties
        X
        Y
        Z
        Label
    end

    methods
        function obj = MimMarkerPoint(x, y, z, label)
            if nargin > 0
                obj.X = x;
                obj.Y = y;
                obj.Z = z;
                obj.Label = label;
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% FROM: LabelMe Toolbox
function h=subplottight(Ny, Nx, j, margin)
% General utility function
%
% This function is like subplot but it removes the spacing between axes.
%
% subplottight(Ny, Nx, j)

if nargin <4
    margin = 0;
end

j = j-1;
x = mod(j,Nx)/Nx;
y = (Ny-fix(j/Nx)-1)/Ny;
h=axes('position', [x+margin/Nx y+margin/Ny 1/Nx-2*margin/Nx 1/Ny-2*margin/Ny]);

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
function progress_line=cosmo_show_progress(clock_start, ratio_done, msg, prev_progress_line)
% Shows a progress bar, and time elapsed and expected to complete.
%
% progress_line=cosmo_show_progress(clock_start, progress[, msg[, prev_progress_line]])
%
% Inputs:
%   clock_start         The time the task started (from clock()).
%   ratio_done          0 <= ratio_done <= 1, where 0 means nothing
%                       completed and 1 means fully completed.
%   msg                 String with a message to be shown next to the
%                       progress bar (optional).
%   prev_progress_line  The output from the previous call to this
%                       function, if applicable (optional). If provided
%                       then invoking this function prefixes the output
%                       with numel(prev_progress_msg) backspace characters,
%                       which deletes the output from the previous call
%                       from the console. In other words, this allows for
%                       showing a progress message at a fixed location in
%                       the console window.
%
% Output:
%   progress_line       String indicating current progress using a bar,
%                       with time elapsed and expected time to complete
%                       (using linear extrapolation).
%
% Notes:
%   - As a side effect of this function, progress_msg is written to standard
%     out (the console).
%   - The use of prev_progress_line may not work properly if output is
%     written to standard out without using this function.
%
% Example:
%   % this code takes just over 3 seconds to run, and fills a progress bar.
%   prev_msg='';
%   clock_start=clock();
%   for k=1:100
%       pause(.03);
%       status=sprintf('done %.1f%%', k);
%       ratio_done=k/100;
%       prev_msg=cosmo_show_progress(clock_start,ratio_done,status,prev_msg);
%   end
%   % output:
%   > +00:00:03 [####################] -00:00:00  done 100.0%
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    if nargin<4 || isempty(prev_progress_line)
        delete_count=0; % nothing to delete
    elseif ischar(prev_progress_line)
        delete_count=numel(prev_progress_line); % count the characters
    end % if not a string, die ungracefully

    if nargin<3 || isempty(msg)
        msg='';
    end
    if ratio_done<0 || ratio_done>1
        error('illegal progress %d: should be between 0 and 1', ratio_done);
    end

    if ratio_done==0
        ratio_to_do=Inf;
    else
        ratio_to_do=(1-ratio_done)/ratio_done;
    end

    clock_now=clock();
    took=etime(clock_now, clock_start);
    eta=ratio_to_do*took; % 'estimated time of arrival'

    % set number of backspace characters
    delete_str=repmat(sprintf('\b'),1,delete_count);

    % define the bar
    bar_width=20;
    bar_done=round(ratio_done*bar_width);
    bar_eta=bar_width-bar_done;
    bar_str=[repmat('#',1,bar_done) repmat('-',1,bar_eta)];

    % because msg may contain the '%' character (which is not to be
    % interpolated) care is needed to ensure that neither building the
    % progress line nor printing it to standard out applies interpolation.
    progress_line=[sprintf('+%s [%s] -%s  ', secs2str(took), bar_str, ...
                                        secs2str(-eta)),...
                   msg];

    if ratio_done==1
        postfix=sprintf('\n');
    else
        postfix='';
    end

    fprintf('%s%s%s',delete_str,progress_line,postfix);

function [m,d]=moddiv(x,y)
    % helper function that does mod and div together so that m+d*y==x
    m=mod(x,y);
    d=(x-m)/y;

function str=secs2str(secs)
    % helper function that formats the number of seconds as
    % human-readable string

    % make secs positive (calling function should add '+' or '-')
    secs=abs(secs);

    if ~isfinite(secs)
        str='oo'; % attempt to look like 'infinity' symbol
        return
    end

    secs=round(secs); % do not provide sub-second precision

    % compute number of seconds, minutes, hours, and days
    [s,secs]=moddiv(secs,60);
    [m,secs]=moddiv(secs,60);
    [h,d]=moddiv(secs,24);

    % add prefix for day, if secs represents at least one day
    if d>0
        daypf=sprintf('%dd+',d);
    else
        daypf='';
    end

    str=sprintf('%s%02d:%02d:%02d', daypf, h, m, s);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
% ASIN   Inverse sine, result in radians.
%    ASIN(X) is the arcsine of the elements of X. Complex
%    results are obtained if ABS(x) > 1.0 for some element.
%
%    See also SIN, ASIND.
%
%    Reference page in Doc Center
%       doc asin
%
%    Other functions named asin
%
%       codistributed/asin    gpuArray/asin    sym/asin    ts/asin
%

author: jmaih


text:
% ir_example_kurtosis1
% example of kurtosis for edge-preserving denoising
% Copyright 2012-07-28, Jeff Fessler, University of Michigan

if ~exist('xtrue')
	N = 2^12;
	xtrue = ones(N,1);
	xtrue(end/2+1:end) = 2;
	plot(xtrue)

	sig = 0.2;
	rng(3)
	yi = xtrue + sig * randn(N,1);

	ax = [1 N 0 3];
	plot(yi, '.-'), axis(ax)

	ah = [0 3];
	hh = linspace(0, 3, 121);
	hp = @(y, c) bar(hh, hist(y, hh), c);

	sf = @(xh) text(.1, get(gca, 'ylim') * [0; 0.7], ...
		sprintf('\\sigma = %3.2f', std(xh(1:round(0.9*N/2)))));

	k1f = @(xh) mean((xh-mean(xh)).^4) / std(xh)^4 - 3;
	kf = @(xh) text(.1, get(gca, 'ylim') * [0; 0.4], ...
		sprintf('\\gamma_2 = %3.2f', k1f(xh(1:round(0.9*N/2)))));
end

% use local psf to help select beta
if ~isvar('R1'), printm 'R1'
	f.l2b = 2;
	mask = true(N,1);
	R1 = Reg1(mask, 'beta', 2^f.l2b);
%	qpwls_psf(1, R1, 1, mask, 1, 'loop', 0); % choose beta
end

if ~isvar('init'), printm 'init'
	psf = qpwls_psf(1, R1, 1, mask, 1);
	init = conv2(psf, yi, 'same');
end

	nn = N/2 + (-100:102);
%	nn = 1:N;
	ax = [minmax(nn)' 0 3];
	clf, pl = @(i,j) subplot(420+(i-1)*2+j);

	pl(1,1)
	plot(nn, xtrue(nn), 'g.-'), axis(ax), xtick(N/2)
	title 'true signal'

	pl(2,1)
	plot(nn, yi(nn), 'r.-'), axis(ax), xtick(N/2)
	title 'data with additive white gaussian noise'
	pl(2,2)
	hp(yi, 'r'), axisx(ah), title 'Histogram'
	sf(yi), kf(yi)

	pl(4,1)
	plot(nn, init(nn), 'm.-'), axis(ax), xtick(N/2)
	pl(4,2)
	hp(init, 'm'), axisx(ah)
	sf(init), kf(init)

if 0
	pl(1,2)
	ii = N/2+1+[-13:13];
	plot(ii, xtrue(ii), 'g.-', ii, init(ii), 'm.-')
	axis([minmax(ii)' 0.5 2.5])
end

if ~isvar('xpwls1'), printm 'pwls1'
	f.niter1 = 20;
	xpwls1 = pwls_pcg1(init(mask), 1, 1, yi, R1, 'niter', f.niter1);
	plot(xpwls1), axis(ax)
end

	pl(3,1)
	plot(nn, xpwls1(nn), 'y.-'), axis(ax), xtick(N/2)
	title 'PWLS with quadratic regularization'
	pl(3,2)
	hp(xpwls1, 'y'), axisx(ah)
	sf(xpwls1), kf(xpwls1)

if ~isvar('R2'), printm 'R2' % use small delta to show kurtosis
	R2 = Reg1(mask, 'beta', 6*2^f.l2b, 'pot_arg', {'huber', 0.01});
%	qpwls_psf(1, R2, 1, mask, 1);
end

if 0 || ~isvar('xpwls2'), printm 'pwls2'
	f.niter2 = 400;
	xpwls2 = pwls_pcg1(init(mask), 1, 1, yi(:), R2, 'niter', f.niter2);
end

	pl(4,1)
	plot(nn, xpwls2(nn), 'c.-'), axis(ax), xtick(N/2)
	title 'PWLS with TV (strongly edge-preserving) regularization'
	pl(4,2)
	hp(xpwls2, 'c'), axisx(ah)
	axis([0 3 0 300])
	sf(xpwls2), kf(xpwls2)

if 0
	pl(1,2)
	ii = N/2+1+[-13:13];
	plot(ii, xtrue(ii), 'g.-', ii, xpwls1(ii), 'y.-', ii, xpwls2(ii), 'c.-')
	axis([minmax(ii)' 0.5 2.5])
end

%	ir_savefig c ir_example_kurtosis1a

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function tsat_values_test ( )

%*****************************************************************************80
%
%% TSAT_VALUES_TEST demonstrates the use of TSAT_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TSAT_VALUES_TEST:\n' );
  fprintf ( 1, '  TSAT_VALUES stores values of\n' );
  fprintf ( 1, '  the saturation temperature\n' );
  fprintf ( 1, '  as a function of pressure.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      P           Tsat(P)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, p, tc ] = tsat_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', p, tc );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Y = sample_cond_multinomial(X, M)
% SAMPLE_MULTINOMIAL Sample Y(i) ~ M(X(i), :)
% function Y = sample_multinomial(X, M)
%
% X(i) = i'th sample
% M(i,j) = P(Y=j | X=i) = noisy channel model
%
% e.g., if X is a binary image,
% Y = sample_multinomial(softeye(2, 0.9), X)
% will create a noisy version of X, where bits are flipped with probability 0.1

if any(X(:)==0)
  error('data must only contain positive integers')
end

Y = zeros(size(X));
for i=min(X(:)):max(X(:))
  ndx = find(X==i);
  Y(ndx) = sample_discrete(M(i,:), length(ndx), 1);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ mD ] = CreateGradientOperator( numRows, numCols )
% ----------------------------------------------------------------------------------------------- %
% [ mD ] = CreateGradientOperator( numRows, numCols )
% Generates a Convolution Matrix for the 2D Gradient of the form [-1, 1].
% Input:
%   - numRows           -   Number of Rows.
%                           Number of rows of the image to be convolved.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
%   - numCols           -   Number of Columns.
%                           Number of columns of the image to be convolved.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
% Output:
%   - mD                -   Convolution Matrix.
%                           The output convolution matrix. The product of
%                           the matrix 'mD' and and image 'mI' in its
%                           column stack form ('mD * mI(:)') is equivalent
%                           to the convolution of 'mI' with the kernel
%                           using the valid convolution shape
%                           ('conv2(mI, mH, 'valid')').
%                           Structure: Matrix (Sparse).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  See my notes - 'Matrix Form of Image Gradient.md'.
% Remarks:
%   1.  The function basically calculates the Convolution Matrix for the
%       kernel [-1, 1] with operation mode of Valid Convolution. See
%       'CreateGradientOperatorUnitTest()' for the exact operation.
%   2.  The matrix 'mT' is the template for creating Vertical / Horizontal
%       derivative operator.
% TODO:
%   1.
%   Release Notes:
%   -   1.0.000     25/03/2020  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

% Vertical Operator - T(numRows)
mT  = spdiags([ones(numRows - 1, 1), -ones(numRows - 1, 1)], [0, 1], numRows - 1, numRows);
mDv = kron(eye(numCols), mT);

% Vertical Operator - T(numCols)
mT  = spdiags([ones(numCols, 1), -ones(numCols, 1)], [0, 1], numCols - 1, numCols);
mDh = kron(mT, eye(numRows));

mD = [mDv; mDh];

end

author: RoyiAvital

repo: StackExchangeCodes


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function FAIRplots(task,varargin);
%
% FAIRplots: generates FAIR plots, see E6_FAIRplots for an example
%
% Create the following 2-by-3 plot:
% _________________________________________________________________
% | R0                 | T0                 |  Tk                 |
% | titleR0            | titleT0            |  title(Tk)          |
% | showImage(R0)      | showImage(T0)      |  showImage(Tk)      |
% |                    |                    |                     |
% |_______________________________________________________________|
% |                    | D0                 |  Dk                 |
% | titleGrid          | titleT0            |  title(Tk)          |
% | showImage(T0) G0   | showDifference(D0) |  showDifference(Tk) |
% | and showGrid  Gk   |                    |                     |
% |_______________________________________________________________|
%
% The function uses a persitent struct plotOptn controling its behaviour
% and perform the following tasks:
%
%------------------------------------------------------------------------------------
% - FAIRplots('clear')
%   clears all persistent variables and returns
%------------------------------------------------------------------------------------
% - FAIRplots('reset',varargin), FAIRplots('set',varargin)
%   initializes various variables and plot-handles (see details below)
%   overwrites defaults by the varargin list
%   note: no plots are shown using set
%   Example: FAIRplots('reset','mode','NPIR','fig',level);
%
%------------------------------------------------------------------------------------
% - FAIRplots('init',struct('Tc',dataT,'Rc',dataR,'omega',omega,'m',m));
%   initializes the figure 'fig' at position 'position,
%     R(xc)            -> (2,3,1) title: Rname
%     T(xc)            -> (2,3,2) title: 'T(xc)'
%     T(xc)            -> (2,3,4) title: 'T(xc)'
%     T(xc)-R(xc)      -> (2,3,5) title: '|T(xc)-R|'
%------------------------------------------------------------------------------------
%   para  = struct('Tc',Tstop,'Rc',Rc,'omega',omega,'m',m,'yc',yStop,'Jc',Dc(Tstop));
% - FAIRplots('stop',para); T = Tstop
%     T(yStop)         -> (2,3,3) title: 'Tstop'
%     T(yStop)-R(xc)   -> (2,3,6) title: '|Tstop-R|=100%'
%     yStop            -> (2,3,4) (update grid)
%------------------------------------------------------------------------------------
%  para.Tc = T0; para.yc = Y0; para.Jstop = para.Jc; para.Jc = Dc(T0);
% - FAIRplots('start',para);
%     T(y0)            -> (2,3,2) title: Tname
%     T(y0)-R(xc)      -> (2,3,5) title: Dname
%     y0               -> (2,3,4) (update grid)
%------------------------------------------------------------------------------------
%   para.Tc = Tk; para.yc = yk; para.Jc = Dc(Tk);
%   para.normdY = norm(yk-Y0)/norm(Ystop);
% - FAIRplots(k,para); (note: k is an integer)
%   T(yc)            -> (2,3,3) title: Tanme
%   T(yc)-R(xc)      -> (2,3,6) title: Dname
%   yc               -> (2,3,4) (replace)
%
%------------------------------------------------------------------------------------
%
%  Details on initialization (default values in []):
%
%   the figure is controlled by FIG[=[]] (handle), position[=screen dependent] (positioning),
%   subplots are controllable via handles [TRG]?HANDLE
%   PLOTS[=1] is used to disable the functionality (set plots=0)
%   MODE[=name of calling function] is used to control the default figurename and reference title:
%     mode==PIR*:
%       figname: PIR*: inter/distance/trafo, dimension/size
%       Rname:   R, m=value of m, length(w)=number of parameters
%     mode==NPIR*:
%       figname: NPIR*: inter/distance/regularizer, dimension/size
%       Rname:   R, m=value of m, alpha==value of alpha
%     otherwise
%       figname: mode, dimension/size
%       Rname:   R, R, m=value of m
%
%    omega and m are used at various places (computation and output)
%
% The final plots and their entiteling  are performed using reconfigurable function
%
%    Function(arguments)  default
%    Tshow(T,omega,m)     @viewImage
%    Tname(k)             @(iter) sprintf('T(%d)',iter)
%
%    Rshow(R,omega,m)     @viewImage
%    Rname(m)             context dependent,
%                         e.g. @(m)sprintf('R, %s, length(w)=%d',dimstr(m),wLen);
%
%    Dshow(T,R,omega,m)   @(T,R,omega,m) viewImage(255-abs(T-R),omega,m)
%    Dname(k,Jc,Jstop)    @(k,Jc,Jstop) ...
%                         sprintf('|J(%d)/Jstop|=%s%%',k,num2str(abs(100*Jc/Jstop)))
%
%    Gshow(yc,omega,m)    context dependent, e.g.
%                         @(yc,omega,m)  plotGrid(yc,omega,m,'spacing',ceil(m/32))
%    Gname(normdY)        @(normdY) sprintf('T(xc), |dY|= %s',num2str(normdY))
%
% run('FAIRplots') or  - see also E6_FAIRplots for an example.
%==============================================================================

function FAIRplots(task,varargin)

if nargin==0
    help(mfilename)
    runMinimalExample;
    return;
end

persistent plotOptn
task  = lower(task);

Error  = @(str) error(sprintf('[%s-task=%s] %s\n',mfilename,task,str));
dimstr = @(m) sprintf('[%s]',sprintf(' %d',m));

% ----- CLEAR plotOptn ----------------------------------------------
if strcmp(task,'clear')
  if nargin == 1,
    clear plotOptions;
    disp('cleared plotOptn')
  else
    Error('nargin>0')
  end;
  return;
end;

% ----- INITIALIZE --------------------------------------------------
if strcmp(task,'reset') | strcmp(task,'set')
  % RESET plotOptn
  if strcmp(task,'reset'),
    plotOptn = [];
  end;

  % INITIALIZE plotOptn
  fields = {
    'fig'
    'plots'
    'figname'
    'position'
    'mode'
    'omega'
    'm'
    'T0handle'
    'Tkhandle'
    'Tshow'
    'Tname'
    'R0handle'
    'Rshow'
    'Rname'
    'D0handle'
    'Dkhandle'
    'Dshow'
    'Dname'
    'G0handle'
    'Gkhandle'
    'Gshow'
    'Gname'
    };

  for j=1:length(fields)
    if ~isfield(plotOptn,fields{j}),
      plotOptn = setfield(plotOptn,fields{j},[]);
    end;
  end;

  % prepare the figname
  J = find(strcmp(varargin,'mode'));
  if isempty(J),
    mode = dbstack;
    mode = mode(min(length(mode),2)).name;
  else
    mode = varargin{min(J)+1};
  end;

  % set defaults
  defaults = {
    'mode',     mode,...
    'plots',    1,...
    'position', [],...
    'Rshow',    @viewImage,...
    'Tshow',    @viewImage,...
    'Tname',    @(iter) sprintf('T(%d)',iter),...
    'Dshow',    @(T,R,omega,m) viewImage(255-abs(T-R),omega,m),...
    'Dname',    @(j,Jc,Jstop) ...
        sprintf('|J(%d)/Jstop|=%s%%',j,num2str(abs(100*Jc/Jstop))),...
    'Gshow',    @(yc,omega,m)  plotGrid(yc,omega,m,'spacing',ceil(m/32)),...
    'Gname',    @(normdY) sprintf('T(xc), |dY|= %s',num2str(normdY))
    };

% percentage of ratio Jc to Jstop
%     'Dname',    @(j,Jc,Jstop) ... % percentage of ratio difference Jstop and Jc to Jstop
%         sprintf('|J(%d)/Jstop|=%s%%',j,num2str(abs(100*(Jstop-Jc)/abs(Jstop)))),...

  for j=1:2:length(defaults),
    if ~isfield(plotOptn,defaults{j}),
      plotOptn.(defaults{j}) = defaults{j+1};
    elseif isempty(plotOptn.(defaults{j})),
      plotOptn.(defaults{j}) = defaults{j+1};
    end;
  end;

  % collect information about current status of the toolbox
  try, wLen = length(trafo('w0')); catch, wLen        = nan;                end;
  try, imgModelStr = imgModel;     catch, imgModelStr = 'imgModel:void';    end;
  try, trafoStr = trafo;           catch, trafoStr    = 'trafo:void';       end;
  try, distStr  = distance;        catch, distStr     = 'distance:void';    end;
  try, regStr   = regularizer;     catch, regStr      = 'regularizer:void'; end;

  % if necessary TUNE figname
  if ~any(strcmp(varargin,'figname')),

    if length(plotOptn.mode) > 2 & strcmp(lower(plotOptn.mode(1:3)),'pir'),
      str = sprintf('%s: %s/%s/%s, ',plotOptn.mode,imgModelStr,distStr,trafoStr);
    elseif length(plotOptn.mode) > 3 & strcmp(lower(plotOptn.mode(1:4)),'npir'),
      str = sprintf('%s: %s/%s/%s, ',plotOptn.mode,imgModelStr,distStr,regStr);
    else
      str = sprintf('%s: ',mode);
    end;
    plotOptn.figname = ...
      @(omega,m) sprintf('%s%dD, m=%s',str,length(omega)/2,dimstr(m));;
  end;

  % if necessary TUNE Rname
  if ~any(strcmp(varargin,'Rname')),
    if length(plotOptn.mode) > 2 & strcmp(lower(plotOptn.mode(1:3)),'pir'),
      plotOptn.Rname = @(m)sprintf('R, %s, length(w)=%d',dimstr(m),wLen);
    elseif length(plotOptn.mode) > 3 & strcmp(lower(plotOptn.mode(1:4)),'npir'),
      plotOptn.Rname = @(m) sprintf('R, %s, \\alpha=%s',...
        dimstr(m),num2str(regularizer('get','alpha')));
    else
      plotOptn.Rname = @(m) sprintf('R, %s',dimstr(m));
    end;
  end;

  % ----- OVERWRITE DEFAULTS ------------------------------------------
  for k=1:2:length(varargin), % overwrites defaults
    if ~isfield(plotOptn,varargin{k}),
      Error(sprintf('field %s not in use',varargin{k}));
    end;
    plotOptn.(varargin{k}) = varargin{k+1};
  end;
  % -------------------------------------------------------------------
  return;
end;
% -------------------------------------------------------------------

if ~plotOptn.plots, return; end;
pause = @(l) builtin('pause',l);
fig = plotOptn.fig;

% extract the parameters for plots
if nargin>1,
  T      = getField(varargin{1},'Tc');
  R      = getField(varargin{1},'Rc');
  omega  = getField(varargin{1},'omega');
  m      = getField(varargin{1},'m');
  yc     = getField(varargin{1},'yc');
  normdY = getField(varargin{1},'normdY');
  Jc     = getField(varargin{1},'Jc');
end;

switch task,

  case 'clear', return;
  case 'reset', return;
  case 'set',   return;

  case 'init',
    % activate figure, set colordef and figname
    fig = FAIRfigure(fig);
%     if isempty(fig), fig = figureh; else fig=figureh(fig); end;
    if ~isnumeric(fig),
        fig = fig.Number;
    end;

    plotOptn.fig = fig;
    % extract variables
    xc = getCellCenteredGrid(omega,m);
    if any(size(T)==1),
      T = reshape(T,m);
      R = reshape(R,m);
    end;
    T  = imgModel(T,omega,xc);
    R  = imgModel(R,omega,xc);

    if ~isempty(plotOptn.position),
      pos = plotOptn.position;
    else
      pos = FAIRposition('fig',fig);
    end;
    figure(fig)
    FAIRfigure(fig,...
      'figname',plotOptn.figname(omega,m),...
      'position',pos);

    if size(omega,2) == 6, % disable Gshow
      plotOptn.Gshow = @(yc,omega,m) [];
    end;

    % plot
    % _____________________________
    % | R       | T     |         |
    % |_________|_______|_________|
    % | T +grid | T -R  |         |
    % |_________|_______|_________|
    figureh(fig); clf
      subplot(2,3,1);
        plotOptn.R0handle = plotOptn.Rshow(R,omega,m);
        title(plotOptn.Rname(m));
      subplot(2,3,2);
        plotOptn.T0handle = plotOptn.Tshow(T,omega,m);
        title('T(xc)');
      subplot(2,3,4);
        plotOptn.G0handle = plotOptn.Tshow(T,omega,m);
        plotOptn.Gkhandle = [];
        title('T(xc)&grid'); hold on;
        axis manual
      subplot(2,3,5);
        plotOptn.D0handle = plotOptn.Dshow(T,R,omega,m);
        plotOptn.Dkhandle = [];
        title('|T(xc)-R|');
        FAIRpause(1/100);
        drawnow;

  case 'stop',
    % plot
    % _____________________________
    % |         |       | Tstop   |
    % |_________|_______|_________|
    % |   +grid |       | Tstop-R |
    % |_________|_______|_________|
    figureh(fig);
      subplot(2,3,3);
        plotOptn.Tkhandle = plotOptn.Tshow(T,omega,m);
        title('T^{stop}');
      subplot(2,3,6);
        plotOptn.Dkhandle = plotOptn.Dshow(T,R,omega,m);
        title('|T^{stop}-R|, J^{stop}=100%');
      subplot(2,3,4);
        set(plotOptn.Gkhandle,'visible','off');
        plotOptn.Gkhandle = plotOptn.Gshow(yc,omega,m);
    pause(1/100); drawnow;
    plotOptn.Jstop = Jc;

  case 'start',
    % plot
    % _____________________________
    % |         | T0    |         |
    % |_________|_______|_________|
    % |   +grid | T0-R  |         |
    % |_________|_______|_________|
    figureh(fig);
      subplot(2,3,2);
        plotOptn.T0handle = plotOptn.Tshow(T,omega,m);
        title(plotOptn.Tname(0));
      subplot(2,3,5);
        plotOptn.D0handle = plotOptn.Dshow(T,R,omega,m);
        title(plotOptn.Dname(0,Jc,plotOptn.Jstop));
     subplot(2,3,4);
       set(plotOptn.Gkhandle,'visible','off');
       plotOptn.Gkhandle = plotOptn.Gshow(yc,omega,m);
    pause(1/100); drawnow;
    plotOptn.J0 = Jc;
  otherwise,
    if ~isnumeric(task),
      warning(['don''t no how to deal task <',task,'>!']);
      return;
    end;

    % plot
    % _____________________________
    % |         |       | Tc      |
    % |_________|_______|_________|
    % |   +grid |       | Tc-R    |
    % |_________|_______|_________|
    figureh(fig)
      subplot(2,3,4);
        set(plotOptn.Gkhandle,'visible','off');
        plotOptn.Gkhandle = plotOptn.Gshow(yc,omega,m);
        title(plotOptn.Gname(normdY));
      subplot(2,3,3);
        plotOptn.Tkhandle = plotOptn.Tshow(T,omega,m);
        title(plotOptn.Tname(task));
      subplot(2,3,6);
        plotOptn.Dkhandle = plotOptn.Dshow(T,R,omega,m);
        title(plotOptn.Dname(task,Jc,plotOptn.Jstop));
    drawnow;
end;

drawnow
%------------------------------------------------------------------------------

function v = getField(s,field);
if isfield(s,field), v = s.(field); else v = []; end;

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

function s = setDefault(s,varargin);
for j=1:2:length(varargin),
  if ~isfield(s,varargin{j}),
    s.(varargin{j}) = varargin{j+1};
  elseif isempty(s.(varargin{j})),
    s.(varargin{j}) = varargin{j+1};
  end;
end;

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

function runMinimalExample
setup2DhandData
imgModel('reset','imgModel','splineInter','regularizer','moments','theta',1e0);
level = 4; omega = ML{level}.omega; m = ML{level}.m;
[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega(1,:),'out',0);
x0  = getCellCenteredGrid(omega(1,:),m);
x1  = rotation2D( 25*pi/180,x0,'c',(omega(2:2:end)-omega(1:2:end))'/2);
x2  = rotation2D(-25*pi/180,x0,'c',(omega(2:2:end)-omega(1:2:end))'/2);
Rc  = imgModel(R,omega(end,:),x0);
T0  = imgModel(T,omega(end,:),x0);
T1  = imgModel(T,omega(end,:),x1);
T2  = imgModel(T,omega(end,:),x2);

FAIRfigure(2);
subplot(1,2,1); viewImage2Dsc(T0,omega,m,'title','template','colormap','gray')
subplot(1,2,2); viewImage2Dsc(Rc,omega,m,'title','reference','colormap','gray')

FAIRplots('reset','mode','testing-PIR','fig',10,'plots',1);
FAIRplots('init',struct('Tc',T,'Rc',R,'omega',omega,'m',m));

para = struct('Tc',T0,'Rc',Rc,'omega',omega,'m',m,'yc',x0,'Jc',100);
FAIRplots('stop',para);
pause(3);

para = struct('Tc',T1,'Rc',Rc,'omega',omega,'m',m,'yc',x1,'Jc',14.3);
FAIRplots('start',para);
pause(3);

para = struct('Tc',T2,'Rc',Rc,'omega',omega,'m',m,'yc',x2,'Jc',0.3);
FAIRplots(14,para);
%==============================================================================


text:
% dbn - training a DBN with up-down algorithm
% Copyright (C) 2011 KyungHyun Cho, Tapani Raiko, Alexander Ilin
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
%
function [D] = dbn(D, patches);

actual_lrate = D.learning.lrate;

n_samples = size(patches, 1);

layers = D.structure.layers;
n_layers = length(layers);

if layers(1) ~= size(patches, 2)
    error('Data is not properly aligned');
end

minibatch_sz = D.learning.minibatch_sz;
n_minibatches = ceil(n_samples / minibatch_sz);

n_epochs = D.iteration.n_epochs;

momentum = D.learning.momentum;
weight_decay = D.learning.weight_decay;

rec.biases_grad = cell(n_layers-1, 1);
rec.W_grad = cell(n_layers-1, 1);
rec.biases_grad_old = cell(n_layers-1, 1);
rec.W_grad_old = cell(n_layers-1, 1);
for l = 1:n_layers-1
    rec.biases_grad{l} = zeros(size(D.rec.biases{l}))';
    if l < n_layers
        rec.W_grad{l} = zeros(size(D.rec.W{l}));
    end
    rec.biases_grad_old{l} = zeros(size(D.rec.biases{l}))';
    if l < n_layers
        rec.W_grad_old{l} = zeros(size(D.rec.W{l}));
    end
end

gen.biases_grad = cell(n_layers-1, 1);
gen.W_grad = cell(n_layers-1, 1);
gen.biases_grad_old = cell(n_layers-1, 1);
gen.W_grad_old = cell(n_layers-1, 1);
for l = 1:n_layers-1
    gen.biases_grad{l} = zeros(size(D.gen.biases{l}))';
    if l < n_layers
        gen.W_grad{l} = zeros(size(D.gen.W{l}));
    end
    gen.biases_grad_old{l} = zeros(size(D.gen.biases{l}))';
    if l < n_layers
        gen.W_grad_old{l} = zeros(size(D.gen.W{l}));
    end
end

top.W_grad = zeros(size(D.top.W));
top.vbias_grad = zeros(size(D.top.vbias));
top.hbias_grad = zeros(size(D.top.hbias));

top.W_grad_old = zeros(size(D.top.W));
top.vbias_grad_old = zeros(size(D.top.vbias));
top.hbias_grad_old = zeros(size(D.top.hbias));

min_recon_error = Inf;
min_recon_error_update_idx = 0;
stopping = 0;

anneal_counter = 0;
actual_lrate0 = actual_lrate;

if D.debug.do_display == 1
    figure(D.debug.display_fid);
end

try
    use_gpu = gpuDeviceCount;
catch errgpu
    use_gpu = false;
    disp(['Could not use CUDA. Error: ' errgpu.identifier])
end

for step=1:n_epochs
    if D.verbose
        fprintf(2, 'Epoch %d/%d: ', step, n_epochs)
    end
    if use_gpu
        % push
        for l = 1:n_layers-1
            if l < n_layers-1
                D.rec.W{l} = gpuArray(single(D.rec.W{l}));
                D.gen.W{l} = gpuArray(single(D.gen.W{l}));
            end
            D.rec.biases{l} = gpuArray(single(D.rec.biases{l}));
            D.gen.biases{l} = gpuArray(single(D.gen.biases{l}));
        end

        D.top.W = gpuArray(single(D.top.W));
        D.top.vbias = gpuArray(single(D.top.vbias));
        D.top.hbias = gpuArray(single(D.top.hbias));
    end

    for mb=1:n_minibatches
        D.iteration.n_updates = D.iteration.n_updates + 1;

        v0 = patches((mb-1) * minibatch_sz + 1:min(mb * minibatch_sz, n_samples), :);
        mb_sz = size(v0,1);

        if use_gpu > 0
            v0 = gpuArray(single(v0));
        end

        % learning rate
        if D.learning.lrate_anneal > 0 && (step >= D.learning.lrate_anneal * n_epochs)
            anneal_counter = anneal_counter + 1;
            actual_lrate = actual_lrate0 / anneal_counter;
        else
            if D.learning.lrate0 > 0
                actual_lrate = D.learning.lrate / (1 + D.iteration.n_updates / D.learning.lrate0);
            else
                actual_lrate = D.learning.lrate;
            end
            actual_lrate0 = actual_lrate;
        end

        D.signals.lrates = [D.signals.lrates actual_lrate];

        % recognition (wake)
        h0r = cell(n_layers-1, 1);
        h0r{1} = v0;
        if D.learning.ffactored == 0
            h0r{1} = binornd(1, h0r{1});
        end

        for l = 2:(n_layers-1)
            h0r{l} = sigmoid(bsxfun(@plus, h0r{l-1} * D.rec.W{l-1}, D.rec.biases{l}'));
            if D.learning.ffactored == 0
                h0r{l} = binornd(1, h0r{l});
            end
        end

        r0 = cell(n_layers-1, 1);
        for l = 1:(n_layers-2)
            r0{l} = sigmoid(bsxfun(@plus, h0r{l+1} * D.gen.W{l}', D.gen.biases{l}'));
        end

        % update generative weight
        for l = 1:(n_layers-2)
            gen.W_grad{l} = ((h0r{l} - r0{l})' * h0r{l+1}) / size(v0, 1);
            gen.W_grad_old{l} = (1 - momentum) * gen.W_grad{l} + momentum * gen.W_grad_old{l};
            D.gen.W{l} = D.gen.W{l} + actual_lrate * (gen.W_grad_old{l} - D.learning.weight_decay * D.gen.W{l});

            gen.biases_grad{l} = mean(h0r{l} - r0{l}, 1);
            gen.biases_grad_old{l} = (1 - momentum) * gen.biases_grad{l} + momentum * gen.biases_grad_old{l};
            D.gen.biases{l} = D.gen.biases{l} + actual_lrate * (gen.biases_grad_old{l}' - D.learning.weight_decay * D.gen.biases{l});
        end

        % contrastive steps
        if D.learning.persistent_cd == 0 || D.iteration.n_updates == 1
            vtop = h0r{end};
        end
        for cs = 1:D.learning.contrastive_step
            htop = sigmoid(bsxfun(@plus, vtop * D.top.W, D.top.hbias'));
            htop = binornd(1, htop);
            vtop = sigmoid(bsxfun(@plus, htop * D.top.W', D.top.vbias'));
            vtop = binornd(1, vtop);
        end

        tv0 = h0r{end};
        th0 = sigmoid(bsxfun(@plus, tv0 * D.top.W, D.top.hbias'));
        tv1 = vtop;
        th1 = sigmoid(bsxfun(@plus, tv1 * D.top.W, D.top.hbias'));

        top.W_grad = (tv0' * th0)/size(tv0,1) - (tv1' * th1)/size(tv1,1);
        top.vbias_grad = mean(tv0, 1) - mean(tv1, 1);
        top.hbias_grad = mean(th0, 1) - mean(th1, 1);

        top.W_grad_old = (1 - momentum) * top.W_grad + momentum * top.W_grad_old;
        top.vbias_grad_old = (1 - momentum) * top.vbias_grad' + momentum * top.vbias_grad_old;
        top.hbias_grad_old = (1 - momentum) * top.hbias_grad' + momentum * top.hbias_grad_old;

        D.top.W = D.top.W + actual_lrate * (top.W_grad_old - D.learning.weight_decay * D.top.W);
        D.top.vbias = D.top.vbias + actual_lrate * (top.vbias_grad_old - D.learning.weight_decay * D.top.vbias);
        D.top.hbias = D.top.hbias + actual_lrate * (top.hbias_grad_old - D.learning.weight_decay * D.top.hbias);

        % generation (sleep)
        h0r{n_layers-1} = tv1;
        for l = (n_layers-2):-1:1
            h0r{l} = sigmoid(bsxfun(@plus, h0r{l+1} * D.gen.W{l}', D.gen.biases{l}'));
            if D.learning.ffactored == 0
                h0r{l} = binornd(1, h0r{l});
            end
        end

        for l = 2:(n_layers - 1)
            r0{l} = sigmoid(bsxfun(@plus, h0r{l-1} * D.rec.W{l-1}, D.rec.biases{l}'));
        end

        % update recognition weight
        for l = 1:(n_layers-1)
            if l < n_layers - 1
                rec.W_grad{l} = (h0r{l}' * (h0r{l+1} - r0{l+1})) / size(h0r{1}, 1);
                rec.W_grad_old{l} = (1 - momentum) * rec.W_grad{l} + momentum * rec.W_grad_old{l};
                D.rec.W{l} = D.rec.W{l} + actual_lrate * (rec.W_grad_old{l} - D.learning.weight_decay * D.rec.W{l});
            end

            if l > 1
                rec.biases_grad{l} = mean(h0r{l} - r0{l}, 1);
                rec.biases_grad_old{l} = (1 - momentum) * rec.biases_grad{l} + momentum * rec.biases_grad_old{l};
                D.rec.biases{l} = D.rec.biases{l} + actual_lrate * (rec.biases_grad_old{l}' - D.learning.weight_decay * D.rec.biases{l});
            end
        end

        clear h0r r0 htop tv0 th0 tv1 th1;
        if D.learning.persistent_cd  == 0
            clear vtop;
        end

        % compute reconstruction error
        hr = v0;
        for l = 2:n_layers-1
            hr = sigmoid(bsxfun(@plus, hr * D.rec.W{l-1}, D.rec.biases{l}'));
        end
        hr = sigmoid(bsxfun(@plus, hr * D.top.W, D.top.hbias'));
        hr = sigmoid(bsxfun(@plus, hr * D.top.W', D.top.vbias'));
        for l = n_layers-1:-1:2
            hr = sigmoid(bsxfun(@plus, hr * D.gen.W{l-1}', D.gen.biases{l-1}'));
        end

        rerr = mean(sum((v0 - hr).^2,2));
        if use_gpu > 0
            rerr = gather(rerr);
        end
        D.signals.recon_errors = [D.signals.recon_errors rerr];

        if D.verbose == 1
            fprintf(2, '.');
        end

        if use_gpu > 0
            clear v0 h0d h0e v0_clean vr hr deltae deltad
        end

        if D.stop.criterion > 0
            if D.stop.criterion == 1
                if min_recon_error > D.signals.recon_errors(end)
                    min_recon_error = D.signals.recon_errors(end);
                    min_recon_error_update_idx = D.iteration.n_updates;
                else
                    if D.iteration.n_updates > min_recon_error_update_idx + D.stop.recon_error.tolerate_count
                        fprintf(2, '\nStopping criterion reached (recon error) %f > %f\n', ...
                            D.signals.recon_errors(end), min_recon_error);
                        stopping = 1;
                        break;
                    end
                end
            else
                error ('Unknown stopping criterion %d', D.stop.criterion);
            end
        end

        if length(D.hook.per_update) > 1
            err = D.hook.per_update{1}(D, D.hook.per_update{2});

            if err == -1
                stopping = 1;
                break;
            end
        end

        if D.debug.do_display == 1 && mod(D.iteration.n_updates, D.debug.display_interval) == 0
            D.debug.display_function (D.debug.display_fid, D, v0, v1, h0, h1, W_grad, vbias_grad, hbias_grad);
            drawnow;
        end
    end

    if use_gpu > 0
        % pull
        for l = 1:n_layers-1
            if l < n_layers-1
                D.rec.W{l} = gather(D.rec.W{l});
                D.gen.W{l} = gather(D.gen.W{l});
            end
            D.rec.biases{l} = gather(D.rec.biases{l});
            D.gen.biases{l} = gather(D.gen.biases{l});
        end

        D.top.W = gather(D.top.W);
        D.top.vbias = gather(D.top.vbias);
        D.top.hbias = gather(D.top.hbias);
    end

    if length(D.hook.per_epoch) > 1
        err = D.hook.per_epoch{1}(D, D.hook.per_epoch{2});

        if err == -1
            stopping = 1;
        end
    end

    if stopping == 1
        break;
    end

    if D.verbose == 1
        fprintf(2, '\n');
    end

    fprintf(2, 'Epoch %d/%d - recon_error: %f\n', step, n_epochs, ...
        D.signals.recon_errors(end));
end

if use_gpu > 0
    % pull
    for l = 1:n_layers-1
        if l < n_layers-1
            D.rec.W{l} = gather(D.rec.W{l});
            D.gen.W{l} = gather(D.gen.W{l});
        end
        D.rec.biases{l} = gather(D.rec.biases{l});
        D.gen.biases{l} = gather(D.gen.biases{l});
    end

    D.top.W = gather(D.top.W);
    D.top.vbias = gather(D.top.vbias);
    D.top.hbias = gather(D.top.hbias);
end

author: kyunghyuncho

repo: deepmat


text:
function [thisFrame,movieFileID,inputMovieDims] = readFrame(inputMoviePath,frameNo,varargin)
	% Fast reading of frame from files on disk. This is an alternative to loadMovieList that is much faster when only a single frame needs to be read.
	% Biafra Ahanonu
	% started: 2020.10.19 [11:45:24]
	% inputs
		% inputMoviePath | Str: path to a movie. Supports TIFF, AVI, HDF5, NWB, or Inscopix ISXD.
		% frameNo | Int: frame number.
	% outputs
		%
	% Usage
		% For non-HDF5 file types, need to open a link to the
			% [thisFrame,movieFileID,inputMovieDims] = ciapkg.io.readFrame(inputMoviePath,frameNo);
			% Then for the second call, feed in movieFileID and inputMovieDims to improve read speed.
			% [thisFrame] = ciapkg.io.readFrame(inputMoviePath,frameNo,'movieFileID',movieFileID,'inputMovieDims',inputMovieDims);

	% changelog
		% 2021.06.30 [01:26:29] - Updated handling of no file path character input.
		% 2021.07.03 [09:02:14] - Updated to have backup read method for different tiff styles.
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
		% 2022.02.24 [10:24:28] - AVI now read(...,'native') is faster.
		% 2022.07.27 [13:48:43] - Ensure subfxn_loadFrame nested function outputs frame data (even if file does not contain the requested frame, e.g. empty matrix).
	% TODO
		%

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

	% ========================
	% Str: hierarchy name in hdf5 where movie data is located
	options.inputDatasetName = '/1';
	% Int vector: list of specific frames to load.
	options.frameList = [];
	% Pointer created by Tiff (tif), VideoReader (avi), or isx.Movie.read (Inscopix ISXD files).
	options.movieFileID = [];
	% Vector: [1st 2nd 3rd] dimension of movie.
	options.inputMovieDims = [];
	% Str: default NWB hierarchy names in HDF5 file where movie data is located, will look in the order indicates
	options.defaultNwbDatasetName = {'/acquisition/TwoPhotonSeries/data'};
	% Binary: 1 = whether to display info on command line.
	options.displayInfo = 1;
	% Binary: Whether to display diagnostic information
	options.displayDiagnosticInfo = 0;
	% Binary: allow RGB display, e.g. for AVI
	options.rgbDisplay = 0;
	% Int: [] = do nothing, 1-3 indicates R,G,B channels to take from multicolor RGB AVI
	options.rgbChannel = [];
	% 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
		movieFileID = [];
		thisFrame = [];
		inputMovieDims = [];

		if ~isempty(options.frameList)
			frameNo = options.frameList(frameNo);
        end

        if iscell(inputMoviePath)==1&&length(inputMoviePath)==1
            inputMoviePath = inputMoviePath{1};
        end

		% Check if path to movie or uder
		if ischar(inputMoviePath)
			[movieType, supported, movieTypeSpecific] = ciapkg.io.getMovieFileType(inputMoviePath);
		else
			disp('Please input a file path.');
			return;
		end

		% Get the movie file identifier to faster access/reading during future calls.
		if isempty(options.movieFileID)
			options.movieFileID = subfxn_getMovieObj(inputMoviePath,options);
			movieFileID = options.movieFileID;
		end

		[thisFrame, inputMovieDims] = subfxn_loadFrame(inputMoviePath,options);

	catch err
		disp(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		disp(repmat('@',1,7))
	end
	function [thisFrame, inputMovieDims] = subfxn_loadFrame(inputMoviePathHere,options)
		inputMovieDims = options.inputMovieDims;
		thisFrame = [];
		switch movieType
			case 'hdf5'
				% Much faster to input the existing movie dimensions.
				if isempty(options.inputMovieDims)
					[inputMovieDims] = ciapkg.io.getMovieInfo(inputMoviePath,'frameList',options.frameList,'inputDatasetName',options.inputDatasetName,'displayInfo',0,'displayDiagnosticInfo',0);
					inputMovieDims = [inputMovieDims.one inputMovieDims.two inputMovieDims.three];
				else
					inputMovieDims = options.inputMovieDims;
				end
				% inputMovieDims
				thisFrame = h5read(inputMoviePathHere,options.inputDatasetName,[1 1 frameNo],[inputMovieDims(1) inputMovieDims(2) 1]);
			case 'tiff'
				warning off;
				try
					if isempty(options.movieFileID)
						% Much slower, avoid.
						tiffID = subfxn_getMovieObj(inputMoviePathHere,options)
					else
						tiffID = options.movieFileID;
					end
					tiffID.setDirectory(frameNo);
					thisFrame = read(tiffID);
				catch
					try
						% thisFrame = imread(inputMoviePathHere, 'tif', 'Index', frameNo, 'Info', imgInfo);
						thisFrame = imread(inputMoviePathHere, 'tif', 'Index', frameNo);
					catch err
						disp(repmat('@',1,7))
						disp(getReport(err,'extended','hyperlinks','on'));
						disp(repmat('@',1,7))
						% fileInfo = imfinfo(filename);
					end
					% err
					% display(repmat('@',1,7))
					% disp(getReport(err,'extended','hyperlinks','on'));
					% display(repmat('@',1,7))
				end
				warning on;
			case 'avi'
				if isempty(options.movieFileID)
					% Much slower, avoid.
					xyloObj = subfxn_getMovieObj(inputMoviePathHere,options)
				else
					xyloObj = options.movieFileID;
				end
				thisFrame = read(xyloObj, frameNo,'native');
				thisFrame = thisFrame.cdata;
				if options.rgbDisplay==0
					if size(thisFrame,3)==3&isempty(options.rgbChannel)
						thisFrame = squeeze(thisFrame(:,:,1));
					elseif ~isempty(options.rgbChannel)
						thisFrame = squeeze(thisFrame(:,:,options.rgbChannel));
					end
				end
			case 'isxd'
				if isempty(options.movieFileID)
					% Much slower, avoid.
					inputMovieIsx = subfxn_getMovieObj(inputMoviePathHere,options)
				else
					inputMovieIsx = options.movieFileID;
				end
				thisFrame = inputMovieIsx.get_frame_data(frameNo-1);
			otherwise
				thisFrame = loadMovieList(inputMoviePathHere,'inputDatasetName',options.inputDatasetName,'displayInfo',0,'displayDiagnosticInfo',0,'displayWarnings',0,'frameList',frameNo);
		end
	end
	function movieFileID = subfxn_getMovieObj(inputMoviePathHere,options)
		% Setup connection to file to reduce I/O for file types that need it.
		movieFileID = [];
		switch movieTypeSpecific
			case 'hdf5'
				%
			case 'nwb'
				try
					h5info(inputMoviePathHere,options.inputDatasetName);
				catch
					for zCheck = 1:length(options.defaultNwbDatasetName)
						try
							options.inputDatasetName = options.defaultNwbDatasetName{zCheck};
							h5info(inputMoviePathHere,options.inputDatasetName);
							fprintf('Correct NWB dataset name: %s\n',options.inputDatasetName)
						catch
							fprintf('Incorrect NWB dataset name: %s\n',options.inputDatasetName)
						end
					end
				end
			case 'tiff'
				warning off;
				try
					tiffID = Tiff(inputMoviePathHere,'r');
					movieFileID = tiffID;
				catch
				end
				warning on;
			case 'avi'
				xyloObj = VideoReader(inputMoviePathHere);
				movieFileID = xyloObj;
			case 'isxd'
				try
					inputMovieIsx = isx.Movie.read(inputMoviePathHere);
				catch
					ciapkg.inscopix.loadIDPS();
					inputMovieIsx = isx.Movie.read(inputMoviePathHere);
				end
				movieFileID = inputMovieIsx;
			otherwise
				%
		end
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function stroud_test2075 ( )

%*****************************************************************************80
%
%% STROUD_TEST2075 tests the rules for EPN with GLG weight on monomials.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2010
%
%  Author:
%
%    John Burkardt
%
  test_num = 5;

  alpha_test = [ - 0.5, 0.0, 0.5, 1.0, 2.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'STROUD_TEST2075\n' );
  fprintf ( 1, '  Demonstrate the use of quadrature rules for the region\n' );
  fprintf ( 1, '  EPN_GLG, that is, the positive half space [0,+oo)^N, with the\n' );
  fprintf ( 1, '  weight W(ALPHA;X) = product ( 1 <= I <= N ) X(I)^ALPHA exp ( -X(I) )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We use the formulas to integrate various monomials of\n' );
  fprintf ( 1, '  the form X(1)^E(1) * X(2)^E(2) * ... X(N)^E(N)\n' );
  fprintf ( 1, '  and compare to the exact integral.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The precision of each formula is known, and we only use\n' );
  fprintf ( 1, '  a formula if its precision indicates it should be able to\n' );
  fprintf ( 1, '  produce an exact result.\n' );

  for n = 1 : 6

    expon = zeros ( n, 1 );
    for test = 1 : test_num
      alpha = alpha_test(test);
      epn_glg_test ( n, expon, alpha );
    end

    expon = zeros ( n, 1 );
    expon(n) = 1;
    for test = 1 : test_num
      alpha = alpha_test(test);
      epn_glg_test ( n, expon, alpha );
    end

    if ( 2 <= n )
      expon = zeros ( n, 1 );
      expon(1) = 1;
      expon(2) = 1;
      for test = 1 : test_num
        alpha = alpha_test(test);
        epn_glg_test ( n, expon, alpha );
      end
    end

    expon = zeros ( n, 1 );
    expon(1) = 2;
    for test = 1 : test_num
      alpha = alpha_test(test);
      epn_glg_test ( n, expon, alpha );
    end

  end

  return
end


text:
function plot_histogram (samples, pd, x_label, x_title)
% plot_histogram: plots histogram from samples (empirical PDF)  and
% compares to inferred probability density function (reference PDF). It
% also plots mean and median.
%
% INPUT
%   samples: Nx1 samples.
%   pd: probality distribution object from ProbabilityDistribution class.
%   x_label: label for X axis (string).
%   x_title: title for the figure (string).
%
% OUTPUT
%   figure with histogram, reference PDF, mean and median.
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% Reference:
%
%
% Version: 004
% Date:    2021/03/02
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego`

N = length(samples);

%% REFERENCE PDF

x = linspace(min(samples), max(samples), N );
ref_pdf = pdf(pd, x);

%% STATISTIC ANALYSIS

mu = mean(samples);
med = median(samples);

i = -5;
idx1 = [];
idx2 = [];

while (isempty(idx1) || isempty(idx2))

    EPS = 10^i;
    idx1 = find( x >= mu - EPS   & x < mu + EPS );
    idx2 = find( x >= med - EPS & x < med + EPS );
    i = i + 1;
end

if ( isempty(idx1) || isempty(idx2) )
    error('plot_histogram: no match for mean or median')
end

% Middlepoints
idx1 = idx1( ceil(end/2) );
idx2 = idx2( ceil(end/2) );

%% PLOT

bins = 100;
blue_new = [0 0.4470 0.7410];
orange_new = [0.8500 0.3250 0.0980];

figure

% Plot histogram from dataquiq
histogram(samples, bins, 'Normalization', 'pdf', 'FaceColor', [.9 .9 .9]);
hold on

% Plot the reference pdf
p0 = plot(x, ref_pdf, '-',  'LineWidth', 2);

% Plot lines
y = ref_pdf (idx1);
l1 = line( [mu, mu] , [0, y], 'Color', blue_new, 'LineWidth', 2, 'LineStyle','-.');

y = ref_pdf (idx2);
l2 = line( [med, med] , [0, y], 'Color', orange_new, 'LineWidth', 2, 'LineStyle','-.' );

legend([p0, l1, l2], 'Reference PDF', 'Mean', 'Median')

xlabel(x_label);
% ylabel('Probability density function (PDF)');
title(x_title)

grid on
hold off

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function [phi,theta,r] = sofa_get_head_orientations(sofa,idx)
%SOFA_GET_HEAD_ORIENTATIONS returns phi, theta and r from the given SOFA data set
%
%   Usage: [phi,theta,r] = sofa_get_head_orientations(sofa,[idx])
%
%   Input parameters:
%       sofa   - impulse response data set (SOFA struct/file)
%       idx    - index of measurement for which head orientation should be
%                returned (default: return all head orientations)
%
%   Output parameters:
%       phi    - head orientations in the horizontal plane / rad
%       theta  - head orientations in the median plane / rad
%       r      - head orientations radii / m
%
%   SOFA_GET_HEAD_ORIENTATIONS(sofa,idx) returns head orientation [phi,theta,r]
%   as defined in the given SOFA file or struct, specified by idx. If no idx is
%   specified, all head orientations are returned.
%
%   See also: get_ir, sofa_get_header, sofa_get_secondary_sources

%*****************************************************************************
% 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 = 1;
nargmax = 2;
narginchk(nargmin,nargmax)
if nargin==nargmax-1
    idx = ':';
else
    isargvector(idx);
end

%% ===== Computation ====================================================
header = sofa_get_header(sofa);
listener_view = SOFAconvertCoordinates(header.ListenerView, ...
                                       header.ListenerView_Type,'spherical');
phi = correct_azimuth(rad(listener_view(idx,1)));
theta = correct_elevation(rad(listener_view(idx,2)));
r = listener_view(idx,3);

author: sfstoolbox

repo: sfs-matlab


text:
function vis_grammar(model)

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-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.
% -------------------------------------------------------

% visualize random derivations...forever
while true
  f = vis_grammar_rand(model);
  visualizeHOG(max(0, f));
  pause;
end

function f = vis_grammar_rand(model, s, p, f)

conf = voc_config();

if nargin < 2
  s = model.start;
  p = [0 0 0];
  f = zeros([0 0 conf.features.dim]);
end

if model.symbols(s).type == 'T'
  w = model_get_block(model, model.filters(model.symbols(s).filter));
  wsz = size(w);
  fsz = size(f);
  req_fsz = [p(2) + wsz(1), p(1) + wsz(2), wsz(3)];
  pad = max(0, req_fsz - fsz);
  f = padarray(f, pad, 0, 'post');
  f(1+p(2):1+p(2)+wsz(1)-1, 1+p(1):1+p(1)+wsz(2)-1, :) = ...
    f(1+p(2):1+p(2)+wsz(1)-1, 1+p(1):1+p(1)+wsz(2)-1, :) + w;
else
  % sample a rule weighted by production score
  len = length(model.rules{s});
  z = zeros(len,1);
  for i = 1:len
    z(i) = model_get_block(model, model.rules{s}(i).offset);
  end
  z = exp(z);
  Z = sum(z);
  if Z ~= 0
    r = find(mnrnd(1, z./Z) == 1);
  else
    r = ceil(rand*length(model.rules{s}));
  end

  if model.rules{s}(r).type == 'D'
    cs = model.rules{s}(r).rhs(1);
    f = vis_grammar_rand(model, cs, p, f);
  else
    for i = 1:length(model.rules{s}(r).rhs)
      cs = model.rules{s}(r).rhs(i);
      anchor = model.rules{s}(r).anchor{i};
      f = vis_grammar_rand(model, cs, p + anchor, f);
    end
  end
end

author: rbgirshick

repo: voc-dpm

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


text:
classdef DenoisingProblem < handle

    properties (GetAccess = public, SetAccess = private)
        optimizedImage
        plottingData
    end

    properties (Access = private)
        originalImage
        imageSize
        discreteGradient
        noisyImage
        noisyImageNorm

        optimizer
        designVariable
        iterator

        energyFunction
        l1Proximal
        dualGap

        lambda
    end

    methods (Access = public)

        function obj = DenoisingProblem(cParams)
            obj.init(cParams);
            obj.createDiscreteGradient();
            obj.createNoisyImage(cParams);
            obj.computeNoisyImageNorm();
            obj.createDesignVariable();
            obj.createIterator(cParams);
            obj.createEnergyFunction(cParams);
            obj.createL1Proximal(cParams);
            obj.createOptimizer(cParams);
        end

        function solve(obj)
            while obj.iterator.hasNotFinished()
                obj.optimizer.update();
                obj.designVariable.update();
                obj.energyFunction.computeCost();
                obj.computeU();
                obj.computeDualGap();
                obj.updatePlottingData();
                obj.iterator.update;
            end
        end

    end

    methods (Access = private)

        function updatePlottingData(obj)
            p = obj.plottingData;
            i = obj.iterator.value;
            p.cost(i) = obj.energyFunction.value;
            p.dualGap(i) = obj.dualGap;
            obj.plottingData = p;
        end

        function init(obj,cParams)
            im = obj.readImage(cParams);
            obj.computeImageSize(im);
            obj.transformImageInVectorForm(im)
            obj.lambda = cParams.totalVariationWeigth;
        end

        function computeImageSize(obj,image)
            [m,n] = size(image);
            obj.imageSize.rows = m;
            obj.imageSize.columns = n;
            obj.imageSize.rowsTimesColumns = m*n;
        end

        function transformImageInVectorForm(obj,image)
            obj.originalImage = image(:);
        end

        function createDiscreteGradient(obj)
            m  = obj.imageSize.rows;
            n  = obj.imageSize.columns;
            mn = obj.imageSize.rowsTimesColumns;
            I  = reshape(1:m*n,m,n);
            east  = [I(:,2:end), I(:,end)];
            north = [I(2:end,:); I(end,:)];
            D1 = sparse(I,east,1,mn,mn)  -speye(mn,mn);
            D2 = sparse(I,north,1,mn,mn) -speye(mn,mn);
            obj.discreteGradient = [D1 ; D2];
        end

        function createNoisyImage(obj,cParams)
            u0    = obj.originalImage;
            L     = cParams.lipschitzConstant;
            a     = cParams.noiseAmplitud;
            sigma = a/L;
            uN    = u0 + sigma*randn(size(u0));
            obj.noisyImage = uN;
        end

        function computeNoisyImageNorm(obj)
            g = obj.noisyImage;
            obj.noisyImageNorm = 0.5*(g'*g);
        end

        function createDesignVariable(obj)
            mn = 2*obj.imageSize.rowsTimesColumns;
            s.xLength = mn;
            obj.designVariable = DesignImagVariable(s);
        end

        function computeU(obj)
            g = obj.noisyImage;
            D = obj.discreteGradient;
            p = obj.designVariable.value;
            obj.optimizedImage = g(:) - D'*p;
        end

        function createEnergyFunction(obj,cParams)
            s.lipschitzConstant = cParams.lipschitzConstant;
            s.A = obj.discreteGradient';
            s.b = obj.noisyImage;
            s.designVariable = obj.designVariable;
            c = QuadraticFunction(s);
            obj.energyFunction = c;
        end

        function createL1Proximal(obj,cParams)
            s.lambda = cParams.totalVariationWeigth;
            s.imageSize = obj.imageSize;
            s.designVariable = obj.designVariable;
            obj.l1Proximal = L1VectorNormProximal(s);
        end

        function createIterator(obj,cParams)
            s.maxIter = cParams.maxIter;
            obj.iterator = Iterator(s);
        end

        function computeDualGap(obj)
            lam = obj.lambda;
            ut  = obj.optimizedImage;
            p   = obj.designVariable.value;
            D   = obj.discreteGradient;
            gN  = obj.noisyImageNorm;
            q   = D'*p;
            gap = lam*sum(abs(D*ut)) + q'*ut;
            obj.dualGap = gap/gN;
        end

        function createOptimizer(obj,cParams)
            sg.designVariable = obj.designVariable;
            sg.differentiableFunction = obj.energyFunction;
            sg.designVariable = obj.designVariable;
            s.gradientMethodParams = sg;

            sm.designVariable = obj.designVariable;
            sm.iterator       = obj.iterator;
            s.momentumParams  = sm;

            s.proximal = obj.l1Proximal;
            s.type     = cParams.optimizer;
            s.iterator = obj.iterator;
            obj.optimizer = SplittingAlgorithm.create(s);
        end

    end

    methods (Access = private, Static)

        function im = readImage(cParams)
            image = cParams.imageFile;
            im = double(imread(image));
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function out=prox_quadratic(x,A,b,alpha)
%PROX_QUADRATIC computes the proximal operator of the function alpha*(0.5*x'*A*x+b'*x)
%
%  Usage:
%  out = PROX_QUADRATIC(X,A,b,alpha)
%  ===========================================
%  INPUT:
%  x - vector to be projected
%  alpha - positive scalar
%  A - positive semidefinite matrix
%  b - vector
%  ===========================================
%  Assumptions:
%  A is psd
%  ===========================================
%  Output:
%  out - proximal operator at x

% This file is part of the FOM package - a collection of first order methods for solving convex optimization problems
% Copyright (C) 2017 Amir and Nili Beck
%
% 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 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 < 4)
    error ('usage: prox_quadratic(x,A,b,alpha)') ;
end

if (alpha < 0)
    error('usage: prox_quadratic(x,A,b,alpha) - alpha should be positive')
end

eps = 1e-10 ;

if (norm( A - A') > eps)
    error('usage: prox_quadratic(x,A,b,alpha) - A should be a symmetric matrix') ;
end
A = A + eps * eye(size(A,1)) ;
A = (A + A') / 2;

[~,p] = chol(A) ;

if (p > 0)
    error('usage: prox_quadratic(x,A,b,alpha) - A should be positive semidefinite"') ;
end

n = length(x) ;

out = (eye(n) + alpha * A) \ (x - alpha * b) ;

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function rad = GB_deg2rad(deg)
%converts from degrees to rad
rad = deg*pi/180;

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtStkRadiation.m                             |
%|    #    |   VERSION    : 0.50                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 31.10.2018                                    |
%|  / 0 \  |   LAST MODIF : 25.11.2018                                    |
%| ( === ) |   SYNOPSIS   : Stoke radiation of a unit sphere using BEM    |
%|  `---'  |                with stokeslet G and stresslet T              |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Parameters
N  = 300;
n  = 100;
x0 = [0.1 0.2 0.3];

% Mesh unit sphere
mesh = mshSphere(N,1);

% Quadrature
gamma = dom(mesh,3);

% Finite element and unknowns
phi = fem(mesh,'P1');
unk = phi.unk;

% Radiation particles
mesh2 = mshSphere(n,2);
X     = mesh2.vtx;

% Graphical rep
figure
plot(mesh)
hold on
plotNrm(mesh)
plot(gamma)
plot(phi,'*r')
plot(msh(X))
alpha(0.5)
axis equal

% Initalization of block matrix for each coordinate (x,y,z)
G      = cell(3,3);
T      = cell(3,3);
mu     = cell(3,1);
lambda = cell(3,1);

% Loop for each coordinate
for i = 1:3
    for j = 1:3
        % Using r = x - y
        % Single layer : G = \int_gamma 1/(8pi) (\delta_ij/r + r_i*r_j/|r|^3)
        name   = ['[ij/r+rirj/r^3]',num2str(i),num2str(j)];
        green  = @(X,Y) 1/(8*pi) .* femGreenKernel(X,Y,name,[]);
        G{i,j} = integral(X,gamma,green,phi);

        % Double layer : T = \int_gamma -6/(8pi) (r_i*r_j*(r.n)/|r|^5)
        T{i,j} = 0;
        for k = 1:3
            name   = ['[rirjrk/r^5]',num2str(i),num2str(j),num2str(k)];
            green  = @(X,Y) -6/(8*pi) .* femGreenKernel(X,Y,name,[]);
            T{i,j} = T{i,j} + integral(X,gamma,green,ntimes(phi,k));
        end
    end

    % mu = [u] = u_int - u_ext = - Gi1 = - Gi1(x0,y)
    name  = ['[ij/r+rirj/r^3]',num2str(i),num2str(1)];
    mu{i} = -1/(8*pi) * femGreenKernel(x0,phi.unk,name,[]);

    % lambda = [sigma] = - Ti1, with nk = xk (because boundary is unit sphere)
    lambda{i} = 0;
    for k = 1:3
        name      = ['[rirjrk/r^5]',num2str(i),num2str(1),num2str(k)];
        lambda{i} = lambda{i} - (-6/(8*pi)) * femGreenKernel(x0,phi.unk,name,[]) .* unk(:,k);
    end
end

% Convert cells to full matrix
G      = cell2mat(G);
T      = cell2mat(T);
mu     = cell2mat(mu);
lambda = cell2mat(lambda);

% Stokes radiation in domain : ui(x) = - sum_j \int_gamma Gij(x,y) lambda_j dy ...
%              + sum_j \int_gamma Tij(x,y).n(y) mu_j dy
sol = -G*lambda + T*mu ;

% Analytic solution : Gi1
ref = cell(3,1);
for i = 1:3
    name   = ['[ij/r+rirj/r^3]',num2str(i),num2str(1)];
    ref{i} = 1/(8*pi) * femGreenKernel(X,x0,name,[]);
end
ref = cell2mat(ref);

% Relative error L2 and inf
norm(ref-sol)/norm(ref)
norm(ref-sol,'inf')/norm(ref,'inf')

% Graphical representation of solution
figure(2)
for i = 1:3
    subplot(1,3,i)
    ind = (i-1)*n + (1:n);
    plot(mesh2,sol(ind))
    axis equal
    colorbar
    title(['Component ',num2str(i)])
end

disp('~~> Michto gypsilab !')


text:
function test_suite = test_chunkize
% tests for cosmo_chunkize
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_chunkize_basis
    ds=cosmo_synthetic_dataset('type','timelock','nreps',8);
    ds.sa.chunks=reshape(repmat((1:8),6,1),[],1);
    ds=cosmo_slice(ds,randperm(48));

    chunks=cosmo_chunkize(ds,8);
    assertEqual(chunks,ds.sa.chunks);

    for j=1:2:7
        chunks=cosmo_chunkize(ds,j);
        eq_chunks=bsxfun(@eq,chunks,chunks');
        eq_ds=bsxfun(@eq,ds.sa.chunks,ds.sa.chunks');

        m=eq_ds & ~eq_chunks;
        assert(~any(m(:)));
    end

    assertExceptionThrown(@()cosmo_chunkize(ds,9),'');
    ds=rmfield(ds.sa,'chunks');
    assertExceptionThrown(@()cosmo_chunkize(ds,2),'');

function test_chunkize_imbalance()
    ds=struct();
    ds.samples=(1:5)';
    assertExceptionThrown(@()cosmo_chunkize(ds,2),'');
    ds.sa.chunks=2+[1 1 2 2 2]';
    assertExceptionThrown(@()cosmo_chunkize(ds,2),'');
    ds.sa.targets=10+[1 2 1 2 2]';
    assertExceptionThrown(@()cosmo_chunkize(ds,3),'');

    count=2;
    res=cosmo_chunkize(ds,count);
    assert_chunkize_ok(ds,res,count);

    ds2=cosmo_stack({ds,ds});
    res2=cosmo_chunkize(ds2,count);
    assert_chunkize_ok(ds2,res2,count);

function test_all_unique_chunks_tiny()
    ds=struct();
    ds.samples=(1:5)';
    ds.sa.targets=2+[1 1 2 2 2]';
    ds.sa.chunks=10+[1 2 3 4 5]';

    for count=2:5
        res=cosmo_chunkize(ds,count);
        assert_chunkize_ok(ds,res,count);
    end

    assertExceptionThrown(@()cosmo_chunkize(ds,6),'');

function test_chunkize_very_unbalanced_chunks_big()
    % all chunks are unique, want a similar number of targets in each
    % output chunk
    ds=cosmo_synthetic_dataset('nreps',6,'ntargets',5);

    nsamples=size(ds.samples,1);
    ds.sa.chunks(:)=repmat((1:nsamples/10),1,10);

    n_combis=max(ds.sa.chunks)*max(ds.sa.targets);

    targets=ds.sa.targets;
    n_swap=5;

    while true
        rp=randperm(nsamples);
        ds.sa.targets=targets;
        ds.sa.targets(rp(1:n_swap))=ds.sa.targets(rp(n_swap:-1:1));
        idxs=cosmo_index_unique({ds.sa.targets,ds.sa.chunks});
        n=cellfun(@numel,idxs);
        if min(n)>=1 && max(n)<=3 && std(n)<.1 && numel(n)==n_combis
            % not too unbalanced
            break;
        end
    end

    nchunks=ceil(3+rand()*4);
    res=cosmo_chunkize(ds,nchunks);
    assert_chunkize_ok(ds,res,nchunks);

function test_chunkize_slight_unbalanced_chunks_big()
    % all chunks are unique, want a similar number of targets in each
    % output chunk
    ds=cosmo_synthetic_dataset('nreps',6,'ntargets',5);

    nsamples=size(ds.samples,1);
    ds.sa.chunks(:)=repmat((1:nsamples/2),1,2);
    ds.sa.targets(1:5)=ds.sa.targets(2:6); % slight imbalance

    nchunks=ceil(rand()*5);
    res=cosmo_chunkize(ds,nchunks);
    assert_chunkize_ok(ds,res,nchunks);

function test_chunkize_all_unique_independent_chunks()
% each sample has its own unique chunk value
    ds=cosmo_synthetic_dataset('ntargets',2,'nchunks',6*6);
    nsamples=size(ds.samples,1);
    ds.sa.chunks(:)=ceil(rand()*10)+(1:nsamples);
    ds.sa.targets=ds.sa.targets(randperm(nsamples));

    nchunks_candidates=[1 2 3 4 6 12 18];
    for nchunks=nchunks_candidates
        chunks=cosmo_chunkize(ds,nchunks);
        assert_chunkize_ok(ds,chunks,nchunks);

        idxs=cosmo_index_unique([ds.sa.targets chunks]);
        n=cellfun(@numel,idxs);

        % require full balance
        assert(all(n(1)==n(2:end)));
    end

function test_chunkize_dependent_balanced_chunks()
% each combination of chunks and targets occurs equally often
    ntargets=ceil(2+rand()*4);
    nreps=ceil(2+rand()*4);
    nchunks=36;
    ds=cosmo_synthetic_dataset('ntargets',ntargets,...
                                'nchunks',nchunks,'nreps',nreps);
    nsamples=size(ds.samples,1);
    ds=cosmo_slice(ds,randperm(nsamples));

    rep_idxs=cosmo_index_unique({ds.sa.chunks,ds.sa.targets});
    assert(all(cellfun(@numel,rep_idxs)==nreps));

    nchunks_candidates=[1 2 3 4 6 12 18];
    for nchunks=nchunks_candidates
        chunks=cosmo_chunkize(ds,nchunks);
        assert_chunkize_ok(ds,chunks,nchunks);

        idxs=cosmo_index_unique([ds.sa.targets chunks]);
        n=cellfun(@numel,idxs);

        % require full balance
        assert(all(n(1)==n(2:end)));
    end

function assert_chunkize_ok(src_ds,chunks,count)
    % number of items must match input dataset
    assertEqual(numel(src_ds.sa.chunks),numel(chunks));

    % must be balanced
    assert_chunks_targets_balanced(src_ds,chunks);

    % cannot have double dipping
    assert_no_double_dipping(src_ds,chunks);

    % must have the proper number of chunks
    assertEqual(numel(unique(chunks)),count);

    assert_chunks_targets_nonzero(src_ds,chunks);

function assert_chunks_targets_balanced(src_ds,chunks)
    idxs=cosmo_index_unique([src_ds.sa.targets chunks]);

    n=cellfun(@numel,idxs);

    % cannot test for 'optimal' balance due to combinatorial explosion;
    % this is a decent approach to make sure that chunks are not too
    % imbalanced
    assert(std(n)<=1.5);
    assert(min(n)+2>=max(n));

function assert_chunks_targets_nonzero(src_ds,chunks)
    [unused,unused,t_idxs]=unique(src_ds.sa.targets);
    [unused,unused,c_idxs]=unique(chunks);

    nt=max(t_idxs);
    nc=max(c_idxs);
    h=zeros(nt,nc);

    ns=numel(chunks);
    for k=1:ns
        t=t_idxs(k);
        c=c_idxs(k);
        h(t,c)=h(t,c)+1;
    end

    assert(all(max(h,[],1)>0));
    assert(all(max(h,[],2)>0));

function assert_no_double_dipping(src_ds,chunks)
    % samples that were in different chunks in src_ds must not be in the
    % same chunk in trg_ds
    [unq_src,unused,src_ids]=unique(src_ds.sa.chunks);
    [unq_trg,unused,trg_ids]=unique(chunks);

    n_src=numel(unq_src);
    n_trg=numel(unq_trg);
    n_samples=numel(src_ds.sa.chunks);
    chunk_count=zeros(n_src,n_trg);

    for k=1:n_samples
        i=src_ids(k);
        j=trg_ids(k);
        chunk_count(i,j)=chunk_count(i,j)+1;
    end

    assert(all(sum(chunk_count>0,2)==1));


text:
function out = isPeriodicTech(f)
%ISPERIODICTECH   Test if a CHEBFUN object is built upon a TECH of
%periodic functions.
%   out = ISPERIODICTECH(F) returns logical true if F has at least one FUN
%   which is made of a TECH of periodic functions and false otherwise.

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

% out = all(cellfun(@(f) isPeriodicTech(f), f.funs));

% Since all funs have the same technology, it is sufficient and faster to test
% only the first.
out = isPeriodicTech(f.funs{1});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun