text:
% Mathematics Q2706108
% https://math.stackexchange.com/questions/2706108
% Lasso ADMM with Positive Constraint
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     25/03/2018
%   *   First release.

%% General Parameters

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

paramLambda = 0.5;

numRows = 12;
numCols = 8;

numIterations   = 250;
stepSize        = 0.0075;

%% Generate Data

mA = randn([numRows, numCols]);
vB = 10 * randn([numRows, 1]);

hObjFun = @(vX) (0.5 * sum(((mA * vX) - vB) .^ 2)) + (paramLambda * norm(vX, 1));
hL2SubGrad = @(vX) vX ./ max(norm(vX, 2), 1e-9);

hSoftThresholdL1 = @(vX, paramLambda) sign(vX) .* max(abs(vX) - paramLambda, 0);
hProjectRPlus = @(vX, paramLambda) max(vX, 0);

%% Solution by CVX

cvx_begin('quiet')
    cvx_precision('best');
    variable vX(numCols)
    minimize( (0.5 * square_pos(norm((mA * vX) - vB, 2))) + (paramLambda * norm(vX, 1)) );
    subject to
        vX >= 0;
cvx_end

disp([' ']);
disp(['CVX Solution Summary']);
disp(['The CVX Solver Status - ', cvx_status]);
disp(['The Optimal Value Is Given By - ', num2str(cvx_optval)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Solution by Sub Gradient Method

vObjValSgm = zeros([numIterations, 1]);

vX = zeros([numCols, 1]);
vObjValSgm(1) = hObjFun(vX);

mAA = mA.' * mA;
vAb = mA.' * vB;

for ii = 2:numIterations
    vG = (mAA * vX) - vAb + (paramLambda * sign(vX));
    vX = vX - (stepSize * vG);
    vX = max(vX, 0); %<! Project onto R+

    vObjValSgm(ii) = hObjFun(vX);
end

disp([' ']);
disp(['Sub Gradient Method Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(vObjValSgm(numIterations))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Solution Proximal Projected Descent

vObjValPgm = zeros([numIterations, 1]);

[vX, mX] = SolveLsL1Prox(mA, vB, paramLambda, numIterations);

for ii = 1:numIterations
    vObjValPgm(ii) = hObjFun(mX(:, ii));
end

disp([' ']);
disp(['Proximal Gradient Method Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(vObjValPgm(numIterations))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Solution by ADMM

vObjValAdmm = zeros([numIterations, 1]);

[vX, mX] = SolveLsL1Admm(mA, vB, paramLambda, numIterations);

for ii = 1:numIterations
    vObjValAdmm(ii) = hObjFun(mX(:, ii));
end

disp([' ']);
disp(['ADMM Method Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(vObjValAdmm(numIterations))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Display Reesults

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);
hAxes       = axes();
hLineSeries = plot([1:numIterations], [cvx_optval * ones([numIterations, 1]), vObjValSgm, vObjValPgm, vObjValAdmm]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(hLineSeries(2:end), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', ['Least Squares with {L}_{1} Regularization (LASSO) with Positivity Constraint'], ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', 'Iteration Index', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', 'Objective Function Value', ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['CVX'], ['Sub Gradient Method'], ['Proximal Gradient Method'], ['ADMM']});
set(hAxes, 'LooseInset', [0.07, 0.07, 0.07, 0.07]);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
A = regressor;

An=bsxfun(@minus,A,mean(A,1));
Bn=bsxfun(@minus,M_0,mean(M_0,1));
tic
An=bsxfun(@times,An,1./sqrt(sum(An.^2,1)));
Bn=bsxfun(@times,Bn,1./sqrt(sum(Bn.^2,1)));
C=sum(An.*Bn,1);
toc

A = rand(3,3);
B = rand(3,3);

author: xiuyechen

repo: FishExplorer


text:
function f=middlepad(f,L,varargin)
%MIDDLEPAD  Symmetrically zero-extends or cuts a function
%   Usage:  h=middlepad(f,L);
%           h=middlepad(f,L,dim);
%           h=middlepad(f,L,...);
%
%   `middlepad(f,L)` cuts or zero-extends *f* to length *L* by inserting
%   zeros in the middle of the vector, or by cutting in the middle
%   of the vector.
%
%   If *f* is whole-point even, `middlepad(f,L)` will also be whole-point
%   even.
%
%   `middlepad(f,L,dim)` does the same along dimension *dim*.
%
%   If *f* has even length, then *f* will not be purely zero-extended, but
%   the last element will be repeated once and multiplied by 1/2.
%   That is, the support of *f* will increase by one!
%
%   Adding the flag `'wp'` as the last argument will cut or extend whole point
%   even functions.  Adding `'hp'` will do the same for half point even
%   functions.
%
%   See also:  isevenfunction, fir2long, fftresample

%   AUTHOR : Peter L. Søndergaard
%   TESTING: OK
%   REFERENCE: OK

if nargin<2
  error('Too few input parameters.');
end;

if  (numel(L)~=1 || ~isnumeric(L))
  error('L must be a scalar');
end;

if rem(L,1)~=0
    error('L must be an integer.');
end;

if L<1
  error('L must be larger than 0.');
end;

% Define initial value for flags and key/value pairs.
definput.flags.centering = {'wp','hp'};
definput.keyvals.dim     = [];

[flags,keyvals,dim]=ltfatarghelper({'dim'},definput,varargin);

[f,L,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(f,L,dim,'MIDDLEPAD');

Lorig=Ls;

% Skip the main section if there is nothing to do. This is necessary
% because some of the code below cannot handle the case of 'nothing to do'
if L~=Ls
  if flags.do_wp

    % ---------------   WPE case --------------------------------------

    if Lorig==1
      % Rather trivial case
      f=[f(1,:);zeros(L-1,W,assert_classname(f))];

    else
      if Lorig>L
        % Cut

        if mod(L,2)==0

          % L even. Use average of endpoints.
          f=[f(1:L/2,:);(f(L/2+1,:)+f(Lorig-L/2+1,:))/2;f(Lorig-L/2+2:Lorig,:)];

        else

          % No problem, just cut.
          f=[f(1:(L+1)/2,:);f(Lorig-(L-1)/2+1:Lorig,:)];

        end;

      else

        d=L-Lorig;

        % Extend
        if mod(Lorig,2)==0

          % Lorig even. We must split a value.

          f=[f(1:Lorig/2,:);...
             f(Lorig/2+1,:)/2;...
             zeros(d-1,W,assert_classname(f));...
             f(Lorig/2+1,:)/2;...
             f(Lorig/2+2:Lorig,:)];

        else
          % Lorig is odd, we can just insert zeros.
          f=[f(1:(Lorig+1)/2,:);zeros(d,W,assert_classname(f));f((Lorig+3)/2:Lorig,:)];

        end;

      end;
    end;

  else

    % ------------------ HPE case ------------------------------------

    if Lorig==1
        f=[f(1,:);zeros(L-1,W,assert_classname(f))];
    else
      if Lorig>L

        d=Lorig-L;
        % Cut

        if mod(L,2)==0
          % L even

          % No problem, just cut.
          f=[f(1:L/2,:);...
             f(Lorig-L/2+1:Lorig,:);];

        else

          % Average of endpoints.
          f=[f(1:(L-1)/2,:);(f((L+1)/2,:)+f(Lorig-(L-1)/2,:))/2;...
             f(Lorig-(L-1)/2+1:Lorig,:);];

        end;

      else

        d=L-Lorig;

        % Extend
        if mod(Lorig,2)==0

          % Lorig even. We can just insert zeros in the middle.

          f=[f(1:Lorig/2,:);...
             zeros(d,W,assert_classname(f));...
             f(Lorig/2+1:Lorig,:)];

        else
          % Lorig odd. We need to split a value in two
          f=[f(1:(Lorig-1)/2,:);...
             f((Lorig+1)/2,:)/2;...
             zeros(d-1,W,assert_classname(f));...
             f((Lorig+1)/2,:)/2;...
             f((Lorig-1)/2+2:Lorig,:)];

        end;

      end;

    end;
  end;

end;

f=assert_sigreshape_post(f,dim,permutedsize,order);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%

% analytical computation of the convolution between
% x(t)=1,0<t<2 and h(t)=1-t,0<t<1

 tx1=-2:.1:0;
 tx2=0:.1:2;
 tx3=2:.1:4;
 tx=[tx1 tx2 tx3];
 x1=zeros(size(tx1));
 x2=ones(size(tx2));
 x3=zeros(size(tx3));
 x=[x1 x2 x3];

 th1=-2:.1:0;
 th2=0:.1:1;
 th3=1:.1:4;
 th=[th1 th2 th3];
 h1=zeros(size(th1));
 h2=1-th2;
 h3=zeros(size(th3));
 h=[h1 h2 h3];

 plot(tx,x,th,h,':*')
 ylim([-.1 1.1])
 legend('x(\tau)','h(\tau)')
 grid

 %reflection
 figure
 plot(tx,x,-th,h,':*')
 legend('x(\tau)','h(-\tau)')
 ylim([-.1 1.1])

 %shifting
 figure
 t=-2;
 plot(tx,x,-th+t,h,':*')
 ylim([-.1 1.1])
 legend('x(\tau)','h(t-\tau)')

 %zero-overlap
 figure
 t=-2;
 plot(tx,x,-th+t,h,':*')
 ylim([-.1 1.1])
 legend('x(\tau)','h(t-\tau)')

 %partial-overlap
 figure
 t= 0.5;
 plot(tx,x,-th+t,h,':*')
 ylim([-.1 1.1])
 legend('x(\tau)','h(t-\tau)')

 T=1;
 r=0:.1:t;
 a=1/T*r+1-t/T;
 hold on; area(r,a);
 hold off;

%complete-overlap
 figure
 t=1.6;
 plot(tx,x,-th+t,h,':*')
 ylim([-.1 1.1])
 legend('x(\tau)','h(t-\tau)')

 r=t-T:.1:t;
 a=1/T*r+1-t/T;
 hold on; area(r,a);
 hold off;

%partial-overlap
 figure
 t=2.4;
 plot(tx,x,-th+t,h,':*')
 ylim([-.1 1.1])
 legend('x(\tau)','h(t-\tau)')

 r=t-T:.1:2;
 a=1/T*r+1-t/T;
 hold on; area(r,a);
 hold off;

 %zero-overlap
 figure
 t=3.6;
 plot(tx,x,-th+t,h,':*')
 ylim([-.1 1.1])
 legend('x(\tau)','h(t-\tau)')

%convolution calculation
 syms t r
 f=1-t+r;
 y=int(f,r,0,t)

 y=int(f,r,t-1,t)
 simplify(y)

 y=int(f,r,t-1,2)

 t1=0:.1:1;
 t2=1:.1:2;
 t3=2:.1:3;
 y1=t1-(t1.^2)/2;
 y2=0.5*ones(size(t2));
 y3=0.5*((3-t3).^2);

 figure
 plot(t1,y1,t2,y2,'.',t3,y3,':')
 ylim([0 0.6])
 title('Output signal y(t)');

author: Sable

repo: mcbench-benchmarks


text:
%SimulationDiskUtil
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Author: Derek Macklin, macklin@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 7/6/2011
classdef SimulationDiskUtil
    methods (Static = true)
        function [simDir, revision, sim] = getSimulation(simDir)
            import edu.stanford.covert.cell.sim.util.CachedSimulationObjectUtil;
            import edu.stanford.covert.cell.sim.util.DiskLogger;

            if ~((isunix && simDir(1) == '/') || (ispc && simDir(2) == ':'))
                simDir = [edu.stanford.covert.cell.sim.util.SimulationDiskUtil.getBaseDir() filesep simDir];
            end

            if nargout >= 2
                metadata = DiskLogger.loadMetadata(simDir);
                revision = metadata.revision;
            end

            if nargout >= 3
                sim = CachedSimulationObjectUtil.load(revision);
            end
        end

        function [simGroupDir, revision, sim] = getLatestSimulationGroup()
            import edu.stanford.covert.cell.sim.util.CachedSimulationObjectUtil;
            import edu.stanford.covert.cell.sim.util.DiskLogger;
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;

            baseDir = edu.stanford.covert.cell.sim.util.SimulationDiskUtil.getBaseDir(); %#ok<*PROP>

            files = dir([baseDir filesep '201*']);
            files = files([files.isdir]);
            fileIdx = [];
            for i = numel(files):-1:1
                if ~SimulationDiskUtil.getNumSimulations(files(i).name)
                    continue;
                end

                tmp = true;
                for j = 1:SimulationDiskUtil.getNumSimulations(files(i).name)
                    if ~exist([baseDir filesep files(i).name filesep num2str(j) filesep 'metadata.mat'], 'file')
                        tmp = false;
                        break;
                    end
                end

                if ~tmp
                    continue;
                end

                for j = 1:SimulationDiskUtil.getNumSimulations(files(i).name)
                    try
                        load([baseDir filesep files(i).name filesep num2str(j) filesep 'metadata.mat'])
                    catch %#ok<CTCH>
                        tmp = false;
                        break;
                    end
                end

                if ~tmp
                    continue;
                end

                fileIdx = i;
                break;
            end
            if isempty(fileIdx)
                simGroupDir = [];
                revision = [];
                sim = [];
                return;
            end
            simTimeStamp = files(fileIdx).name; %simulation group
            simGroupDir = sprintf('%s%s%s', baseDir, filesep, simTimeStamp);

            if nargout >= 2
                simIdx = 1; %simulation within group
                simDir = sprintf('%s%s%d', simGroupDir, filesep, simIdx);
                metadata = DiskLogger.loadMetadata(simDir);
                revision = metadata.revision;
            end

            if nargout >= 3
                sim = CachedSimulationObjectUtil.load(revision);
            end
        end

        function [simDir, revision, sim] = getLatestSimulation()
            import edu.stanford.covert.cell.sim.util.CachedSimulationObjectUtil;
            import edu.stanford.covert.cell.sim.util.DiskLogger;
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;

            simGroupDir = SimulationDiskUtil.getLatestSimulationGroup();

            simIdx = 1; %simulation within group
            simDir = sprintf('%s%s%d', simGroupDir, filesep, simIdx);

            if nargout >= 2
                metadata = DiskLogger.loadMetadata(simDir);
                revision = metadata.revision;
            end

            if nargout >= 3
                sim = CachedSimulationObjectUtil.load(revision);
            end
        end

        function simGroup = getLatestWildTypeSimulationGroup()
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;

            baseDir = SimulationDiskUtil.getBaseDir();
            files = dir([baseDir filesep '201*']);
            files = files([files.isdir]);
            [~, order] = sort({files.name});
            files = files(order(end:-1:1));
            for i = 1:numel(files)
                for j = 1:SimulationDiskUtil.getNumSimulations(files(i).name)
                    simDir = [baseDir filesep files(i).name filesep num2str(j)];
                    simGroup = files(i).name;

                    if ...
                            ~exist([simDir filesep 'metadata.mat'], 'file') || ...
                            ~exist([simDir filesep 'options.mat'], 'file') || ...
                            ~exist([simDir filesep 'parameters.mat'], 'file') || ...
                            exist([simDir filesep 'skipAnalysis'], 'file')
                        continue;
                    end

                    options = load([simDir filesep 'options.mat']);
                    if isempty(options.geneticKnockouts)
                        return;
                    end
                end
            end
        end

        function simGroup = getLatestSingleGeneDeletionSimulationGroup()
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;

            baseDir = SimulationDiskUtil.getBaseDir();
            files = dir([baseDir filesep '201*']);
            files = files([files.isdir]);
            [~, order] = sort({files.name});
            files = files(order(end:-1:1));
            for i = 1:numel(files)
                for j = 1:SimulationDiskUtil.getNumSimulations(files(i).name)
                    simDir = [baseDir filesep files(i).name filesep num2str(j)];
                    simGroup = files(i).name;

                    if ...
                            ~exist([simDir filesep 'metadata.mat'], 'file') || ...
                            ~exist([simDir filesep 'options.mat'], 'file') || ...
                            ~exist([simDir filesep 'parameters.mat'], 'file') || ...
                            exist([simDir filesep 'skipAnalysis'], 'file')
                        continue;
                    end

                    options = load([simDir filesep 'options.mat']);
                    if isscalar(options.geneticKnockouts)
                        return;
                    end
                end
            end
        end

        function [metaData, options, parameters] = getSimulations(simBatchDir, selectedSimulations)
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;
            baseDir = SimulationDiskUtil.getBaseDir(); %#ok<*PROP>

            if nargin >= 1
                simDir = SimulationDiskUtil.getSimulation([simBatchDir filesep '1']);
                simBatchDir = simDir(1:find(simDir == filesep, 1, 'last')-1);
                files = struct('name', simBatchDir, 'isdir', ~~exist(simBatchDir, 'dir'));
            else
                files = dir([baseDir filesep '201*']);
            end
            files = files([files.isdir]);
            metaData = repmat(struct(...
                'firstName', [], ...
                'lastName', [], ...
                'email', [], ...
                'affiliation', [], ...
                'userName', [], ...
                'hostName', [], ...
                'ipAddress', [], ...
                'revision', [], ...
                'differencesFromRevision', [], ...
                'shortDescription', [], ...
                'longDescription', [], ...
                'knowledgeBaseWID', [], ...
                'startTime', [], ...
                'endTime', [], ...
                'lengthSec', [], ...
                'outputDirectory', [], ...
                'segmentSizeStep', [], ...
                'downsampleStepSec', [], ...
                'baseDirectory', [], ...
                'directory', [], ...
                'simGroup', [], ...
                'simIdx', [], ...
                'name', [], ...
                'stateNames', [], ...
                'dependentStateNames' ,[]), 0, 1);
            options = repmat(struct, 0, 1);
            parameters = repmat(struct, 0, 1);
            for i = 1:numel(files)
                if nargin < 2
                    selectedSimulations = 1:SimulationDiskUtil.getNumSimulations(files(i).name);
                end
                for j = 1:numel(selectedSimulations)
                    if (isunix && files(i).name(1) == '/') || (ispc && length(files(i).name) >= 2 && files(i).name(2) == ':')
                        simDir = [files(i).name filesep num2str(selectedSimulations(j))];
                        simDir = strrep(simDir, '/', filesep);
                        simDir = strrep(simDir, '\', filesep);
                        idxs = find(files(i).name == '/' | files(i).name == '\');
                        simGroup = files(i).name(idxs(end)+1:end);
                    else
                        simDir = [baseDir filesep files(i).name filesep num2str(selectedSimulations(j))];
                        simGroup = files(i).name;
                    end

                    if exist([simDir filesep 'skipAnalysis'], 'file')
                        continue;
                    end

                    if nargout >= 1 && exist([simDir filesep 'metadata.mat'], 'file')
                        try
                            tmpMetadata = load([simDir filesep 'metadata.mat']);
                        catch %#ok<CTCH>
                            continue;
                        end
                        tmpMetadata.baseDirectory = baseDir;
                        tmpMetadata.directory = simDir;
                        tmpMetadata.simGroup = simGroup;
                        tmpMetadata.simIdx = selectedSimulations(j);
                        tmpMetadata.name = [files(i).name ' - ' num2str(selectedSimulations(j)) ' - ' tmpMetadata.shortDescription];

                        tmp = setdiff(fieldnames(metaData), fieldnames(tmpMetadata));
                        for k = 1:numel(tmp)
                            tmpMetadata.(tmp{k}) = [];
                        end

                        metaData = [metaData; tmpMetadata]; %#ok<AGROW>
                    end

                    if nargout >= 2 && exist([simDir filesep 'options.mat'], 'file')
                        tmpOptions = load([simDir filesep 'options.mat']);

                        if isempty(options)
                            options = tmpOptions;
                        else
                            tmp = setdiff(fieldnames(options), fieldnames(tmpOptions));
                            for k = 1:numel(tmp)
                                tmpOptions.(tmp{k}) = [];
                            end

                            tmp = setdiff(fieldnames(tmpOptions), fieldnames(options));
                            for k = 1:numel(tmp)
                                for l = 1:numel(options)
                                    options(l).(tmp{k}) = [];
                                end
                            end

                            options = [options; tmpOptions]; %#ok<AGROW>
                        end
                    end

                    if nargout >= 3 && exist([simDir filesep 'parameters.mat'], 'file')
                        tmpParameters = load([simDir filesep 'parameters.mat']);

                        if isempty(parameters)
                            parameters = tmpParameters;
                        else
                            tmp = setdiff(fieldnames(parameters), fieldnames(tmpParameters));
                            for k = 1:numel(tmp)
                                tmpParameters.(tmp{k}) = [];
                            end

                            tmp = setdiff(fieldnames(tmpParameters), fieldnames(parameters));
                            for k = 1:numel(tmp)
                                for l = 1:numel(parameters)
                                    parameters(l).(tmp{k}) = [];
                                end
                            end

                            parameters = [parameters; tmpParameters]; %#ok<AGROW>
                        end
                    end
                end
            end
        end

        function value = getBaseDir()
            config = getConfig();
            value = config.outputPath;
        end

        function value = getSimulationIndex(simDirAbs)
            % Get rid of any trailing slashes
            while simDirAbs(end) == '/'
                simDirAbs(end)=[];
            end

            p = find(simDirAbs == filesep, 1, 'last');
            value = str2double(simDirAbs(p+1:end));
        end

        function value = getSimulationTimeStamp(simDirAbs)
            % Get rid of any trailing slashes
            while simDirAbs(end) == filesep
                simDirAbs(end) = [];
            end

            p = find(simDirAbs == filesep, 2, 'last');
            value = simDirAbs(p(1) + 1:p(2) - 1);

			value = strrep(value, '_', ':');
			value(11) = ' ';
        end

        function value = getSimulationBatchDir(simDirAbs)
            % Get rid of any trailing slashes
            while simDirAbs(end) == '/'
                simDirAbs(end) = [];
            end

            p = find(simDirAbs == filesep, 2, 'last');
            value = simDirAbs(1:p(2) - 1);
            while value(end) == filesep
                value(end) = [];
            end
        end

        function value = getNumSimulations(simBatchDir)
			if (isunix && simBatchDir(1) == '/') || (ispc && simBatchDir(2) == ':')
				baseDir = simBatchDir;
			else
				baseDir = [edu.stanford.covert.cell.sim.util.SimulationDiskUtil.getBaseDir filesep simBatchDir];
			end
            files = dir(baseDir);
            files = files([files.isdir] & ~isnan(str2double({files.name})));
            value = numel(files);
        end

        function value = getCompleteSimulations(simBatchDir)
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;

            if (isunix && simBatchDir(1) == '/') || (ispc && simBatchDir(2) == ':')
                baseDir = simBatchDir;
            else
                baseDir = [edu.stanford.covert.cell.sim.util.SimulationDiskUtil.getBaseDir filesep simBatchDir];
            end

            nSims = SimulationDiskUtil.getNumSimulations(simBatchDir);
            tfs = false(nSims, 1);
            for i = 1:nSims
                tfs(i) = exist([baseDir filesep num2str(i) filesep 'summary.mat'], 'file');
            end
            value = find(tfs);
        end

        function warnings = getSimulationWarnings(simBatchDir, simIdx)
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;

            fid = fopen([SimulationDiskUtil.getSimulation([simBatchDir filesep num2str(simIdx)]) filesep 'out.log'], 'r');
            if fid == -1
                throw(MException('SimulationDiskUtil:error', 'Unable to open file'))
            end

            warningTemplate = struct();
            warningTemplate.message = [];
            warningTemplate.files = {};
            warningTemplate.lineNumbers = [];
            warningTemplate.times = [];
            warnings = repmat(warningTemplate, 0, 1);
            isWarn = false;
            isWarnTrace = false;
            time = 0;
            while ~feof(fid)
                line = fgetl(fid);

                if isWarnTrace && (length(line) < 5 || ~isequal(line(1:5), '  In '))
                    isWarn = false;
                    isWarnTrace = false;

                    tmpWarning.message = tmpWarning.message(12:end-4);

                    repeatedWarning = 0;
                    for j = 1:numel(warnings)
                        if isequal(warnings(j).message, tmpWarning.message) && ...
                                isequal(warnings(j).files, tmpWarning.files) && ...
                                isequal(warnings(j).lineNumbers, tmpWarning.lineNumbers)
                            repeatedWarning = j;
                            break;
                        end
                    end

                    if repeatedWarning
                        warnings(repeatedWarning).times = ...
                            [warnings(repeatedWarning).times; tmpWarning.times];
                    else
                        warnings = [warnings; tmpWarning]; %#ok<AGROW>
                    end
                end
                if numel(line) >= 11 && isequal(line(3:11), 'Warning: ')
                    isWarn = true;
                    isWarnTrace = false;
                    tmpWarning = warningTemplate;
                end
                if isWarn && line(1) == '>'
                    isWarnTrace = true;
                end

                if ~isWarn && length(line) >= 1 && ~all(line == ' ') && any(line == ' ')
                    idx1 = find(line ~= ' ', 1, 'first');
                    idx2 = find(line(idx1:end) == ' ', 1, 'first') - idx1;
                    if ~isempty(idx1) && ~isempty(idx2) && idx2 >= idx1
                        time = max(time, str2double(line(idx1:idx2)));
                    end
                elseif isWarn && ~isWarnTrace
                    tmpWarning.times = time;
                    tmpWarning.message = [tmpWarning.message line ' '];
                elseif isWarnTrace
                    tmp = line(6:end);
                    tmpWarning.files{end+1, 1} = tmp(1:strfind(tmp, ' at ')-1);
                    tmpWarning.lineNumbers(end+1, 1) = str2double(tmp(strfind(tmp, ' at ')+4:end));
                end
            end

            assert(~isWarn);
            assert(~isWarnTrace);

            fclose(fid);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [gx,dG_dX,dG_dPhi] = g_fullDCM4fmri(Xt,Phi,ut,inG)
% full DCM observation function (embedding dynamic Ballon model)

% function [gx,dG_dX,dG_dPhi] = g_fullDCM4fmri(Xt,Phi,ut,inG)
% This function evaluates the observation function of the neuronal level in
% DCM for fMRI. Note that this "generalized" observation function includes
% the HRF convolution model (balloon model).

persistent xh dxhdx dxhdp ph ti ko
% These variables are described bellow:
%   - xh{i}: current hemodynamic states for region i.
%   - dxhdx{i}: derivatives of the hemodynamic states of region i w.r.t.
%   neuronal state of region i.
%   - dxhdp{i}: derivatives of the hemodynamic states of region i w.r.t.
%   hemodynamic parameters of region i.
%   - ph: the last remembered hemodynamic parameters. This is only used
%   when this function is called from f_fullDCM4fmri.m. It is updated each
%   time this function is called outside f_fullDCM4fmri.m.
%   - ti: the index within the time series. This is used to reinitialize
%   the persistent variables.
% These are reinitialized at the beginning of the time series.

try % just to get lastly updated phi

    inG.getPhi;
    ph = Phi;

catch

    n = size(Xt,1);

    % Check whether the system is at initial state
    if isempty(xh)
        ti = 0;
        xh = kron(ones(n,1),[0;0;0;0]);
        dxhdx = zeros(n,n*4);
        dxhdp = zeros(size(inG.Phi,1),n*4);
        ko = kron(eye(n),[1 0 0 0]);
    end

    if isempty(Phi)    %- Function is called from f_fullDCM4fmri

        % Get current observation parameters mode
        if ~isempty(ph)
            Phi = ph;
        else % only the first iteration
            Phi = inG.Phi;
        end

        %- hemodynamic response convolution operation
        [fx,dfdxh,dfdp] = f_HRF2(xh,Phi,Xt,inG);

        xh = fx;
        % update gradients
        dxhdx = dxhdx*dfdxh + ko;
        dxhdp = dfdp + dxhdp*dfdxh;

    else    %- Function is called outside of f_fullDCM4fmri

        % Update persistent observation parameters
        ph = Phi;
        ti = ti+1;

        %- Get predicted observations and gradients
        [gx,dgdxh,dgdp] = g_HRF3(xh,Phi,ut,inG);
        %
        %     dgdxh = numericDiff('g_HRF',1,xh,Phi,ut,inG);
        %     dgdp = numericDiff('g_HRF',2,xh,Phi,ut,inG);

        dG_dX = dxhdx*dgdxh;
        dG_dPhi = dgdp + dxhdp*dgdxh;
        %         % empty gradients w.r.t. states
        %         dxhdx = zeros(n,n*4);

    end

    if ti==inG.n_t
        xh = [];
    end
end

author: MBB-team

repo: VBA-toolbox


text:
function afni_fig_interface(src, evnt)
   cf = gcf;
   figure(src);
   % evnt
   ud = get(src,'UserData');
   if (isempty(ud) | ~isfield(ud,'axes_lock')),
      ud = default_ud();
      set(src,'UserData',ud);
   end

   if (  lower(evnt.Character) == 'l'  ),
      ud.axes_lock = ~ud.axes_lock;
      set (src,'UserData',ud);
      if ud.axes_lock, set (src,'Name','L');
      else set (src,'Name','');
      end
   elseif (  lower(evnt.Character) == 'a'  ),
      if (ud.axes_lock),
         reset_axis(evnt.Character, get_all_axes(src));
      else
         reset_axis(evnt.Character,gca); %can send arrays for characters and axes
      end
   elseif (  lower(evnt.Character) == 'x'  |...
         lower(evnt.Character) == 'y'  |...
         lower(evnt.Character) == 'z' ),
      if (ud.axes_lock),
         zoom_axis(evnt.Character, get_all_axes(src));
      else
         zoom_axis(evnt.Character,gca); %can send arrays for characters and axes
      end
   elseif ( strcmp(evnt.Key , 'leftarrow') |...
            strcmp(evnt.Key , 'rightarrow') |...
            strcmp(evnt.Key , 'uparrow') |...
            strcmp(evnt.Key , 'downarrow')),
      if (ud.axes_lock),
         shift_axis(evnt, get_all_axes(src));
      else
         shift_axis(evnt, gca);
      end
   end
   figure(cf);
return;

function reset_axis(a,h)
   ca = gca;
   for (ih=1:1:length(h)),
      axes(h(ih));
      axis auto;
   end
   axes(ca);
return

function zoom_axis(a,h)
   ca = gca;
   for (ia=1:1:length(a)),
      for (ih=1:1:length(h)),
         axes(h(ih));
         dv = 4;
         if (  ~(a(ia) - lower(a(ia))) ),
            dv = 1;
         end
         if (~(lower(a(ia)) - 'x')), se = [1 2];
         elseif (~(lower(a(ia)) - 'y')), se = [3 4];
         elseif (~(lower(a(ia)) - 'z')), se = [5 6];
         else
            return;
         end
         v = axis;
         if (length(v) < max(se)), return; end
         vc = (v(se(1))+v(se(2)))/2;
         vd = v(se(2))-v(se(1));
         v(se(1)) = vc-vd/dv;
         v(se(2)) = vc+vd/dv;
         axis(v);
      end
   end
   axes(gca);
return;

function shift_axis(e,h)
   ca = gca;
   for (ie=1:1:length(e)),
      for (ih=1:1:length(h)),
         axes(h(ih));
         if (strcmp(e(ie).Key , 'leftarrow')),
            se = [1 2]; dv = +0.1;
         elseif (strcmp(e(ie).Key , 'rightarrow')),
            se = [1 2]; dv = -0.1;
         elseif (strcmp(e(ie).Key , 'downarrow')),
            se = [3 4]; dv = -0.1;
         elseif (strcmp(e(ie).Key , 'uparrow')),
            se = [3 4]; dv = 0.1;
         else
            return;
         end
         v = axis
         vd = v(se(2))-v(se(1));
         v(se(1)) = v(se(1))+vd*dv;
         v(se(2)) = v(se(2))+vd*dv;
         axis(v);
      end
   end
   axes(gca);
return;

function v = get_all_axes(f)
   l = get(f, 'Children');
   v = [];
   for (i=1:1:length(l)),
      if (strcmp(get(l(i),'Type'),'axes')), v = [v l(i)]; end
   end
return

function ud = default_ud()
   ud = struct('axes_lock', 0);
return

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function shepard_interp_nd_test02 ( prob, p, m, n1d )

%*****************************************************************************80
%
%% SHEPARD_INTERP_ND_TEST02 tests SHEPARD_INTERP_ND on a regular N1D^M grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROB, the problem number.
%
%    Input, real P, the power used in the distance weighting.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N1D, the number of points in 1D.
%

%
%  Set problem parameters:
%
  seed = 123456789;
  [ c, seed ] = r8vec_uniform_01 ( m, seed );
  [ w, seed ] = r8vec_uniform_01 ( m, seed );

  nd = n1d^m;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SHEPARD_INTERP_ND_TEST02:\n' );
  fprintf ( 1, '  Interpolate data from TEST_INTERP_ND problem #%d\n', prob );
  fprintf ( 1, '  using Shepard interpolation with P = %g\n', p );
  fprintf ( 1, '  spatial dimension M = %d\n', m );
  fprintf ( 1, '  and a regular grid of N1D^M = %d data points.\n', nd );

  a = 0.0;
  b = 1.0;
  x1d = r8vec_linspace ( n1d, a, b );

  xd = zeros ( m, nd );
  for i = 1 : m
    xd = r8vec_direct_product ( i, n1d, x1d, m, nd, xd );
  end

  zd = p00_f ( prob, m, c, w, nd, xd );
%
%  #1:  Does interpolant match function at interpolation points?
%
  ni = nd;
  xi = xd;
  zi = shepard_interp_nd ( m, nd, xd, zd, p, ni, xi );

  int_error = norm ( zi - zd ) / ni;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  L2 interpolation error averaged per interpolant node = %g\n', int_error );
%
%  #2: Approximation test.  Estimate the integral (f-interp(f))^2.
%
  ni = 1000;
  [ xi, seed ] = r8mat_uniform_01 ( m, ni, seed );
  zi = shepard_interp ( m, nd, xd, zd, p, ni, xi );
  ze = p00_f ( prob, m, c, w, ni, xi );

  app_error = norm ( zi - ze ) / ni;

  fprintf ( 1, '  L2 approximation error averaged per 1000 samples = %g\n', app_error );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [tc,relres,iter,xrec] = gablasso(x,g,a,M,lambda,varargin)
%GABLASSO  LASSO regression in Gabor domain
%   Usage: [tc,xrec] = gablasso(x,a,M,lambda,C,tol,maxit)
%
%   `gablasso` has been deprecated. Please use |franalasso| instead.
%
%   A call to `gablasso(x,g,a,M,lambda)` can be replaced by ::
%
%     F=frame('dgt',[],g,a,M);
%     tc=franalasso(F,lambda);
%
%   Any additional parameters passed to `gablasso` can be passed to
%   |franalasso| in the same manner.
%
%   See also: frame, franalasso

warning(['LTFAT: GABLASSO has been deprecated, please use FRANALASSO ' ...
         'instead. See the help on GABLASSO for more details.']);

F=newframe('dgt',[],g,a,M);
[tc,relres,iter,xrec] = franalasso(F,lambda,varargin{:});


text:
%XML_LOAD Loads XML file and return contents as Matlab structure or variable
%
% SYNTAX     V = xml_load( filename )
%            V = xml_load( filename, attswitch )
%
% INPUT
%  filename   filename of xml file to load (if extension .xml is omitted,
%             xml_load tries to append it if the file cannot be found).
%             This file should have been saved with xml_save().
%             If it's a non-Matlab XML file, please use e.g.
%             xmlstr = fileread(filename); V = xml_parseany(xmlstr);
%  attswitch  optional, default='on':
%             'on' takes attributes idx,size,type into account for
%                  creating corresponding Matlab data types
%             'off' ignores attributes in XML element headers
%
% OUTPUT
%   V        Matlab structure or variable
%
% SEE ALSO
%   xml_save, xml_format, xml_formatany, xml_parse, xml_parseany, xml_help, (xmlread, xmlwrite)

% Copyright (C) 2002-2005, University of Southampton
% Author: Dr Marc Molinari <m.molinari@soton.ac.uk>
% $Revision: 1.1 $ $Date: 2005/04/15 17:12:14 $ $Tag$

author: Cloud-CV

repo: object-proposals


text:
function c8mat_uniform_01_test ( )

%*****************************************************************************80
%
%% C8MAT_UNIFORM_01_TEST tests C8MAT_UNIFORM_01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2013
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8MAT_UNIFORM_01_TEST\n' );
  fprintf ( 1, '  C8MAT_UNIFORM_01 computes a "random" complex matrix.\n' );

  seed = 123456789;

  [ a, seed ] = c8mat_uniform_01 ( m, n, seed );

  c8mat_print ( m, n, a, '  The matrix:' );

  return
end


text:
function [fboxes,cs] = transfer_friends(models, bboxes)
%transfer friends onto detections
error('deprecated function.. see VOCinit');

VOCinit;
% if isfield(models{i},'detect_add_flip') && models{i}.detect_add_flip == 1
%   models{i}.gt_box = flip_box(models{i}.gt_box,[recs.size.height ...
%                     recs.size.width]);
% end

hasfriends = cellfun(@(x)size(x.friendbb,1),models);
hasfriends = find(hasfriends);

fboxes = cell(1,length(bboxes));
cs = cell(1,length(bboxes));

for i = 1:length(bboxes)
  curb = bboxes{i};
  goods = ismember(curb(:,6),hasfriends);
  curb = curb(goods,:);

  fboxes{i} = zeros(0,7);
  cs{i} = cell(0,1);
  for j = 1:size(curb,1)
    exid = curb(j,6);
    curg = models{exid}.gt_box;
    others = models{exid}.friendbb;

    if isfield(models{exid},'detect_add_flip') && models{exid}.detect_add_flip == 1
      sizeI = models{exid}.sizeI;
      %curg = flip_box(curg,sizeI);
      others = flip_box(others,sizeI);
    end

    xform_cd = find_xform(curg, curb(j,:));

    %news = others;

    if isfield(models{exid},'detect_add_flip') && models{exid}.detect_add_flip == 1
      others = flip_box(others,sizeI);
    end

    news = apply_xform(others,xform_cd);

    %if isfield(models{exid},'detect_add_flip') && models{exid}.detect_add_flip == 1
    %  news = flip_box(news,sizeI);
    %end

    news(:,5) = j;
    news(:,6) = curb(j,6);
    news(:,7) = curb(j,end);
    fboxes{i} = [fboxes{i}; news];

    cs{i} = [cs{i} models{exid}.friendclass];

  end

end

author: quantombone

repo: exemplarsvm


text:
% Estimate sinusoid parameters using peak picking or LS solutions
%
% Octave compatible
%
% Description
%  Model a given waveform using sinusoids whose parameters are estimated
%  using various methods. Depending on the combinations of the options:
%  fharmonic, fquasiharm, fadapted and use_ls, the methods below are available:
%                    | fharmonic | fquasiharm | fadapted  | use_ls
%  Peak Peaking [1]     false       false        false      false
%  Peak Peaking [1]     true        false        false      false
%  Peak Peaking [1]     true        true         false      false
%  HM [2] **            true        false        false      true    (default)
%  QHM (partly [3]*)    true        true         false      true
%  aHM [5]              true        false        true       true
%  (iQHM [3] to come)
%  (aQHM [4] to come)
%
%  ** This is a full-band HM, this is NOT the Harmonic+Noise Model (HNM)
%  *  Note that [3] is NOT fully implemented. The QHM estimation is available,
%     but NOT the iterative algorithm presented in [3].
%
%  The analysis instants need to be always provided through the f0s argument.
%
%  The window length is always odd so as the sample at the middle of the
%  window correspond to the sample of the analysis instant.
%
% Inputs
%  wav    : The waveform
%  fs     : [Hz] The sampling frequency
%  f0s    : [s, Hz] [Nx2] A temporal vector with time instants and fundamental
%           frequency f0 estimated at the given time instants.
%           The sinusoid parameters are always estimated at these time instants.
%           For the Peak Peaking method, the f0 (second column) can be omitted.
%  [opt]  : Additional options (see code below)
%
% Outputs
%  frames : N structures containing the estimated sinusoid parameters and extra
%           information (e.g. window length, the f0 used).
%           For each frame, the sinusoid parameters are in a matrix with format:
%              [5xK] for each column: the frequency [Hz], the linear amplitude,
%              the instantaneous phase [rad], the harmonic number of each
%              sinusoidal component and a boolean specifying if the sinusoidal
%              parameters are from a spectral peak or through sampling.
%              The DC is ALWAYS included at the beginning of the matrix.
%  syn    : if asked, the resynthesized waveform using an Overlap-Add method.
%  opt    : The options structure which might have been altered for consistency
%           purpose.
%
% Example
%  Please se the HOWTO_sinusoidal example
%
% References
%  [1] McAulay, R., Quatieri, T.: Speech analysis/Synthesis based on a sinusoidal
%      representation, IEEE Transactions on Acoustics, Speech and Signal
%      Processing 34(4):744-754, 1986.
%  [2] Stylianou, Y.: Harmonic plus Noise Models for Speech combined with
%      Statistical Methods, for Speech and Speaker Modification, TelecomParis,
%      PhD Thesis, 1996.
%  [3] Pantazis, Y., Rosec, O., Stylianou, Y.: Iterative Estimation of Sinusoidal
%      Signal Parameters, Signal Processing Letters, IEEE 17(5):461-464, 2010.
%  [4] Pantazis, Y., Rosec, O., Stylianou, Y.: Adaptive AM-FM Signal Decomposition
%      With Application to Speech Analysis, IEEE Transactions on Audio, Speech,
%      and Language Processing 19(2):290-300, 2010.
%  [5] Degottex, G., Stylianou, Y.: Analysis and Synthesis of Speech using an
%      Adaptive Full-band Harmonic Model, IEEE Transactions on Acoustics, Speech
%      and Language Processing, 21(10):2085-2095, 2013.
%
% Copyright (c) 2012 University of Crete - Computer Science Department
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file 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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function [frames, syn, opt] = sin_analysis(wav, fs, f0s, opt)

    if nargin<4
        % Options

        % Window
        opt.win_durf0sync = true;
        opt.win_durnbper  = 3;       % Number of period per window (def. 3)
                                     % (used only if win_durf0sync==true)
        opt.win_dur       = 30/1000; % [s] Duration of the analysis window
                                     % (used only if win_durf0sync==false)
        opt.win_fn        = @blackman; % Window type
        opt.win_dropoutside=true;    % Drop windows which are partly outside of
                                     % the signal (and thus drop the
                                     % corresponding analysis instants).

        % Partials estimation
        opt.fharmonic = true; % Use harmonic or quasi-harmonic frequencies
        opt.fquasiharm= false;% Use quasi-harmonic frequencies
                              % It is possible to estimate the parameters using
                              % the QHM model. HOWEVER, the iterative algorithm
                              % in [3] is not implemented yet !
        opt.fadapted = false; % Adapt the frequency basis to the f0 curve
                              % (currently, works only with LS solution)

        % For LS solution [2-5]
        opt.use_ls   = true;  % Use the Least Square solution (LS) [2,3]
                              % Otherwise, use Peak Picking (PP) [1]
        opt.win_ls_marg = 0.1;% For the LS solution, add 10% to the
                              % theoretical minimum window length to ensure the
                              % stability of the LS solution

        % For Peak Picking [1]
        opt.normstr       = 'sum(win)'; % Normalize the window using the given
                                % expression.
                                % Among the available terms, the following exist
                                % win, winlen and dftlen
        opt.dftlen        = []; % Force the DFT to a given length
                                % if empty, the DFT length is adapted to the
                                % window length + the oversampling factor below.
                                % (ignored when using the LS solution)
        opt.osf           = 2;  % Frequency OverSampling Factor, according to:
                                %    dftlen=2^(nextpow2(winlen)+opt.osf)
        opt.frames_keepspec = false; % Keep the halfspec
        opt.funique       = true;% Keep only unique frequencies
                                 % Considered only if harmonic=false

        % Resynthesis
        opt.resyn    = false; % Do OverLap-Add (OLA) resynthesis

        opt.debug    = 1;

        opt.fmapin{1}  = {{'wav', 'fs'}, 'snd'};
        opt.fmapin{2}  = {{'f0s'}, 'yin'};
        opt.fmapout{1} = {{'frames'}, 'mat'};
    end
    if nargin==0; frames=opt; return; end

    if opt.debug>0; disp('Sinusoidal Analysis'); end
    if opt.fquasiharm; opt.fharmonic=true; end
    if opt.use_ls; opt.fharmonic=true; end % Currently, the LS solution is only possible with a harmonic model
%      if ~opt.use_ls && opt.fharmonic==true; opt.fquasiharm=true; end % commented: otherwise we can't do harmonic PP
    if ~opt.use_ls && opt.fquasiharm==true; opt.fharmonic=true; end
    if size(f0s,2)>1 && any(f0s(:,2)<=0); error('If a fundamental frequency is specified, it cannot be zero.'); end
    if opt.fharmonic
        if size(f0s,2)<2; error('An input f0 curve is necessary for harmonic model'); end
        opt.win_durf0sync=true;
    end
    if opt.fadapted
        if size(f0s,2)<2; error('An input f0 curve is necessary for adaptivity'); end
        if ~opt.use_ls; error('Adaptivity with peak picking is not possible. Should use LS solution.'); end
    end
    if nargout<2; opt.resyn=false; end
    if opt.debug>0; disp(opt); end

    wav = wav(:);
    if opt.resyn
        syn  = zeros(size(wav));
        wins = zeros(size(wav));
    else
        syn = [];
    end

    T = f0s(:,1); % Get analysis time instants

    if opt.fadapted
        times = (0:length(wav)-1)'/fs;
        % Sample the first harmonic all along the signal
        f1 = interp1(f0s(:,1), f0s(:,2), times, 'spline');
        f1 = interp1_extrapbounds(f1); % Check if bounds are defined, and replace nan values
        % Compute the fundamental phase
        p1 = filter(1, [1 -1], 2*pi*f1/fs);
        % p1 = 2*pi*cumtrapz(f1)/fs; % Create Bad conditionned matrices !
    end

    % Use a constant winlen if winlen is not pitch sync
    if ~opt.win_durf0sync
        if size(f0s,2)>1
            f0m = exp(median(log(f0s(:,2))));
            winlen = round(opt.win_durnbper*fs/f0m/2)*2+1;
        else
            winlen = round(opt.win_dur*fs/2)*2+1;
        end
        win = opt.win_fn(winlen);
        if ~isempty(opt.normstr)
            eval(['d = ' opt.normstr ';']);
            win = win./d;
        end
        if ~isempty(opt.dftlen); dftlen=opt.dftlen;
        else                     dftlen=2^(nextpow2(winlen)+opt.osf); end
        W = delay2spec((winlen-1)/2, dftlen);
    end

    if opt.debug>0; pb = progressbar(length(T)); end
    for ind=1:length(T)

        % Be sure the analysis instant is on the sample of the window center
        T(ind) = round(T(ind)*fs)/fs;
        fr.t = T(ind);

        if opt.win_durf0sync
            f0 = f0s(ind,2);
            fr.f0 = f0;

            if opt.fadapted; winlen = get_optimal_winlen(f0s, fs, ind, opt);
            else             winlen = round(opt.win_durnbper*fs/f0/2)*2+1; end
            win = opt.win_fn(winlen);

            if ~isempty(opt.dftlen); dftlen=opt.dftlen;
            else                     dftlen=2^(nextpow2(winlen)+opt.osf); end
            W = delay2spec((winlen-1)/2, dftlen);
        end

        fr.winlen = winlen;
        fr.dftlen = dftlen;

        winids = -(winlen-1)/2:(winlen-1)/2; % Indices relative to the center
        idscenter = round(T(ind)*fs)+1;
        ids = idscenter + winids; % Indices of the window in the signal

        if opt.win_dropoutside
            if ids(1)<1 || ids(end)>length(wav);
                T(ind) = NaN;
                continue;
            end
            iddx = (1:winlen);
            idsb = ids;
            wavsel = wav(ids);
        else
            iddx = find(ids>=1 & ids<=length(wav)); % Valid indices of the window
            idsb = ids(iddx); % Indices of the win in the sig bounded by the sig limits
            wavsel = zeros(winlen,1);
            wavsel(iddx) = wav(idsb);
        end

        if opt.use_ls
            % Use Least Squares (LS) solution [2]
            if opt.fadapted
                % Adapt the frequency basis to the f0 curve
                Ho = floor(((fs/2)-max(f1(idsb))/2)/max(f1(idsb)));

                % Use interpolated phase (through interpolated frequencies)
                if length(idsb)<winlen;
                    % Extrapolate p1 if the window is partly outside of the sig
                    p1sel = interp1(idsb, p1(idsb), ids, 'nearest', 'extrap')';
                else
                    p1sel = p1(ids);
                end
                pm = p1sel - p1(idscenter);
                pm = pm*(-Ho:Ho);
                Nbk = size(pm,2);

                if opt.fquasiharm % Adaptive Quasi-Harmonic Model (aQHM) [4]
                    error('Adaptivity + Quasi-harmonicity (aQHM) not implemented ! ... yet');

                else % Adaptive Harmonic Model (aHM) [5]
                    % Build matrices to compute the LS solution of ak
                    E = cos(pm)+1j*sin(pm); % dimension of E: (2N+1)x(Ho*2+1)
                    Ew = repmat(win,1,Nbk).*E;
                    R = Ew'*Ew;
                    fr.RCN = rcond(R); % Estimate the matrix condition number

                    x = R\(Ew'*(wavsel.*win)); % The LS solution
                    ak = x(Ho+1:end);     % amplitudes (skip the negative freqs)
                end

                % Get the f0 and the center of the window
                cf0 = f1(idscenter);
                fr.sins = [cf0*(0:Ho); abs(ak)'; angle(ak)'; (0:Ho); zeros(1,Ho+1)];

            else
                % Use stationary components
                Ho = floor(((fs/2)-f0/2)/f0);
                fk = f0*(-Ho:Ho);
                Nbk = length(fk);

                if opt.fquasiharm % Quasi-harmonics
                    pm = winids'*2*pi*fk/fs;
                    E = cos(pm)+1j*sin(pm);
                    E = [E repmat(winids',1,Nbk).*E];
                    Ew = repmat(win,1,2*Nbk).*E;
                    R = Ew'*Ew;

                    fr.RCN = rcond(R);   % Estimate the matrix condition number

                    x = R\(Ew'*(wavsel.*win));   % The LS solution
                    ak = x(Ho+1:Nbk);       % amplitudes

                else % simple-harmonics [2]
                    pm = winids'*2*pi*fk/fs;
                    E = cos(pm)+1j*sin(pm);
                    Ew = repmat(win,1,Nbk).*E;
                    R = Ew'*Ew;

                    fr.RCN = rcond(R); % Estimate the matrix condition number

                    ak = R\(Ew'*(wavsel.*win)); % The LS solution

                    ak = ak(Ho+1:end);     % amplitudes (skip the negative freqs)

                    %  y = real(E*ak); % reconstructed signal
                end

                fr.sins = [fk(Ho+1:end); abs(ak)'; angle(ak)'; (0:Ho); zeros(1,Ho+1)];
            end
        else
            % Use Peak Picking (PP) from a spectrum [1]
            if opt.win_durf0sync && ~isempty(opt.normstr)
                eval(['d = ' opt.normstr ';']);
                win = win./d;
            end

            % Window the signal segment
            s = wavsel.*win;

            % Compute the spectrum and compansate the window delay
            S = fft(s, dftlen).';
            S = S.*W;

            % TODO add FChT

            if opt.fharmonic
                % Select only peaks around harmonic frequencies
                fr.sins = spec_getsins_f0(S, fs, f0);

                if ~opt.fquasiharm
                    % Force harmonic frequencies after partial estimation
                    mf0 = median(diff(fr.sins(1,:)));
                    fr.sins(1,:) = mf0*(0:size(fr.sins,2)-1);
                end
            else
                % If not quasi-harmonic => unconstrained sinusoids
                fr.sins = spec_getsins(S, fs);

                if opt.funique
                    [~, idx] = unique(fr.sins(1,:));
                    fr.sins = fr.sins(:,idx);
                end
            end

            if opt.frames_keepspec; fr.S=spec2hspec(S); end
        end

        if opt.resyn
            y = sin2sig(fr.sins, fs, winlen);
            y = 2*y;
            fr.SNR = mag2db(std(wavsel)/std(wavsel-y));
            syn(idsb) = syn(idsb) + y(iddx).*win(iddx);
            wins(idsb) = wins(idsb) + win(iddx);
        end

        if ind==1
            frames(length(T)) = fr; % pre-allocate with correct fieldnames
        end
        frames(ind)=fr;

        if 0 && T(ind)>0.3
            hold off;
            plot(wavsel, 'k');
            hold on;
            plot(y, 'b');
            keyboard
        end

        if 0
            % Compute the spectrum and compansate the window delay
            % TODO Check the behavior of the harmonic structure with the aDFT.
            win = win./sum(win);
            S = fft(wavsel.*win, dftlen);
            F = fs*(0:length(S)-1)/length(S);
            hold off;
            plot(F, ld(S), 'k');
            hold on;
            stem(fr.sins(1,:), ld(fr.sins(2,:)), 'xr');
            xlim([0 fs/2]);
%              ylim([-140 40]);
%              keyboard
            pause
        end

        if 0 && opt.debug>1 && T(ind)>0.3
            V3clear();
            V3spec(S, fs, 'k');
            V3part(frames(end).sins, fs);
            keyboard
        end

        if opt.debug>0; pb = progressbar(pb, ind); end
    end
    if opt.debug>0; pb = progressbar(pb, length(T)); end

    % Drop the necessary frames
    idx = find(~isnan(T));
    if length(idx)<length(T) && opt.debug>0
        disp(['    Some windows were outside of the signal. ' num2str(length(T)-length(idx)) ' frames dropped. (use opt.win_dropoutside=false if you want to keep all windows and zero-pad the necessary ones at signal boundaries).']);
    end
    T = T(idx);
    frames = frames(idx);

    if opt.resyn
        idx = find(wins>0);
        syn(idx) = syn(idx)./wins(idx);
        if opt.debug>0; disp(['Mean SNR=' num2str(mean([frames.SNR]))]); end
    end

    if opt.debug>1
        Hmax = 0;
        for ind=1:numel(frames)
            Hmax = max(Hmax,size(frames(ind).sins,2));
        end
        harmstruct = NaN*ones(numel(frames), Hmax);
        for ind=1:numel(frames)
            M = min(length(frames(ind).sins(1,:)),Hmax);
            harmstruct(ind,1:M) = frames(ind).sins(1,1:M);
        end
        hold off;
        [mf0, winlen, dftlen] = spec_info(fs, f0s);
        [X, Fs, Ts] = spectrogram(wav, blackman(winlen), round(0.9*winlen), dftlen, fs);
        imagesc(Ts, Fs, lin2db(X));
        axis xy;
        hold on;
        plot([0 (length(wav)-1)/fs], 0.5*fs*[1 1], '--k');
        plot(T, harmstruct, 'k');
        keyboard
    end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
% Run a simulation with modified enzyme kinetics
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 1/24/2013
classdef ModifiedMetabolicEnzymeKinetics < edu.stanford.covert.cell.sim.runners.SimulationRunner
    properties
        %Structure containing modified metabolic enzyme kinetics.
        %kCatFor, kCatRev should have units (rxn/enzyme/s).
        %Set kCatFor, kCatRev to NaN to use default value.
        modifiedKinetics = repmat(struct(...
            'reactionWholeCellModelID', '', ...
            'kCatFor', NaN, ...
            'kCatRev', NaN ...
            ), 0, 1);
    end

    methods
        function this = ModifiedMetabolicEnzymeKinetics(varargin)
            this = this@edu.stanford.covert.cell.sim.runners.SimulationRunner(varargin{:});
        end
    end

    methods (Access = protected)
        function modifyNetworkParameters(this, sim)
            %get handles
            met = sim.process('Metabolism');

            %calculate indices of modified reactions
            ids = {this.modifiedKinetics.reactionWholeCellModelID}';
            kCatFor = [this.modifiedKinetics.kCatFor]';
            kCatRev = [this.modifiedKinetics.kCatRev]';

            if numel(unique(ids)) < numel(ids)
                throw(MException('ModifiedMetabolicEnzymeKinetics:invalidReactionList', ...
                    'Modified reaction list cannot contain multiple entries for a given reaction.'));
            end

            idxs = met.reactionIndexs(ids);

            %modify kcats
            met.enzymeBounds(idxs(~isnan(kCatRev)), 1) = kCatRev(~isnan(kCatRev));
            met.enzymeBounds(idxs(~isnan(kCatFor)), 2) = kCatFor(~isnan(kCatFor));
            met.fbaEnzymeBounds(met.fbaReactionIndexs_metabolicConversion, :) = met.enzymeBounds(met.reactionIndexs_fba, :);
        end
    end
end


text:
function [qt]=matrix(qt, sz1, sz2)
%Converts a vector to a matrix in the QTT-Tucker format
%   [QT]=MATRIX(QT, SZ1, SZ2) Makes factors of QTT-Tucker to be TT-matrices
%   SZ1 and SZ2 are cell arrays with N, M dimensions for each factor
%   If they are absent, the square matrices are assumed
%
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

if (nargin<3)
    sz1 = [];
    sz2 = [];
end;

d = qt.dphys;
for i=1:d
    if (~isempty(sz1))||(~isempty(sz2))
        qt.tuck{i} = tt_matrix(qt.tuck{i}, sz1{i}, sz2{i});
    else
        cursz = qt.tuck{i}.n;
        cursz = sqrt(cursz);
        qt.tuck{i} = tt_matrix(qt.tuck{i}, cursz, cursz);
    end;
end;

end

author: oseledets

repo: TT-Toolbox


text:
function [CVA] = spm_cva_compare (Y,X,c)
% Model comparison for probabilistic CVA
% FORMAT [CVA] = spm_cva_compare (Y,X,c)
%
% Y  [N x d1] data matrix
% X  [N x d2] design matrix
% c  Contrast vector (if specified)
%
% CVA has fields:
%
% .order        number of canonical vectors (latent space dimension)
% .bic          BIC for each order
% .aic          AIC for each order
%
% and
%
% .U1,.U2      Canonical vectors
% .W1,.W2      Factor matrices
%
% for the highest order model.
%
% See spm_cva_prob.m for more details
%___________________________________________________________________________
% Copyright (C) 2011 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_cva_compare.m 4687 2012-03-14 18:15:49Z will $

if nargin > 2
    % get null-space of contrast
    X0  = X - X*c*pinv(c);
    X   = full(X*c);
    X0  = spm_svd(X0);

    Y     = Y - X0*(X0'*Y);
    X     = X - X0*(X0'*X);
end

[N1,p1]=size(Y);
[N2,p2]=size(X);
if ~(N1==N2)
    disp('X and Y are of incompatible size');
    return
end

m=min([p1,p2]);
for i=1:m+1,
    order(i)=i-1;
    CVA = spm_cva_prob (Y',X',order(i));
    bic(i)=CVA.bic;
    aic(i)=CVA.aic;
    L(i)=CVA.L;
end

CVA.order=order;
CVA.bic=bic;
CVA.aic=aic;
CVA.L=L;

author: spm


text:
function [tt]=tt_kron7(tt1,tt2,s)
%Kronecker product of two TT-tensors in the reverse order
%   [TT]=TT_KRON7(TT1,TT2,S)
% May 26, 2011
% Vladimir Kazeev
% vladimir.kazeev@gmail.com
% INM RAS
% Moscow, Russia

d1 = size(tt1,1);
d2 = size(tt2,1);

tt1{1}=permute(tt1{1},[(1:s),s+2,s+1]);
%tt2{1}=permute(tt2{1},[(1:s),s+2,s+1]);

tt = cell(d1+d2,1);
for i=1:d2
	tt{i}=tt2{i};
end
for i=1:d1
	tt{d2+i}=tt1{i};
end

%tt{1}=permute(tt{1},[(1:s),s+2,s+1]);

return
end


text:
function out = ConvertToCursor(i)

% Converts an RGB image to a matlab compatiable cursor CData.
%
% White Pixels and black pixels will be kept. All other values
% will be converted to NaN (transparent)
%
% Author: Richard Medlock. (2003)

nRows = size(i,1);
nCols = size(i,2);

for r = 1:nRows

    for c = 1:nCols

        if i(r,c,1) == 255 & i(r,c,2) == 255 & i(r,c,3) == 255 % White

            out(r,c) = 2;

        elseif i(r,c,1) == 0 & i(r,c,2) == 0 & i(r,c,3) == 0 % Black

            out(r,c) = 1;

        else % Any other colour

            out(r,c) = NaN;

        end
    end
end


text:
function catalogobj = associate(obj, maxTimeDiff, sites, source)
%ASSOCIATE Associate detections into events
% catalogobj = associate(detectionObj, maxTimeDiff) will scan through an
% Detection object and look
% for times where there are at least 2 detections on
% within maxTimeDiff seconds of each other, and declare an event. No
% checking is done to see if arrivals are on different channels
%
% catalogobj = associate(arrivals, maxTimeDiff, sites) will
% reduce the detection times by traveltimes in the sites structure. This
% allows a smaller maxTimeDiff to be used (since all airwaves should now
% line up, and all P waves etc.).
%
% Example: Imagine you have run dbdetect. You also notice that the time to cross
% the network is less than 30s. Then the following code will load the detections, subset
% for "D" states, and associate into events wherever 2 channels have an "D"
% detection within 30s. The final line returns a Catalog object - the GISMO
% container for multiple events.
%
%   dbpath = '/home/t/thompsong/NEWTON/pavlof2007';
%   detobj = Detection.retrieve(dbpath, 'status=~/D/');
%   catalogobj = detobj.associate(30);
%
% A significant problem is that there may be more than 1 event in a
% 30-second window. So a different approach is to apply a differential
% travel time correction to each detection, to reduce it to a known source
% location. Then a wave might propagate across the reduced network in a
% time of 1-second, for example, allowing multiple events to be detected in
% a 30-seconds time window.
%
%   dbpath = '/home/t/thompsong/NEWTON/pavlof2007/db';
%   detobj = Detection.retrieve(dbpath, 'status=~/D/');
%   sites = antelope.dbget_site_locations(dbpath, unique(detobj.channelinfo))
%   source.lat = 55.4203;
%   source.lon = -161.8931;
%   wavespeed = 330; % m/s
%   difftt = calc_traveltime(source, sites, wavespeed);
%   source.lat = 55.4173; source.lon = -161.8937; elev = 2.518;
%   chantags = ChannelTag(unique(detobj.channelinfo))
%   sites = antelope.dbget_site_locations(dbpath, chantags, startTime, endTime);
%   seismicspeed = 330; infrasoundspeed = 330; % m/s
%   sites = compute_travel_times(source, sites, seismicspeed, infrasoundspeed);
%   association_time_window = 2; % max seconds from seismic to infrasound arrival
%   catalogobj = detobj.associate(maxtimediff, sites);

    % break up detection types
    detections_on = obj.subset('state', 'ON');
    detections_d = obj.subset('state', 'D');
    detections_off = obj.subset('state', 'OFF');
    if detections_d.numel > 0
        obj = detections_d;
    elseif detections_on.numel > 0
        obj = detections_on;
    end

    %% REDUCE BY SUBTRACTING TRAVEL TIME
    % If sites exist, let's correct the arrival times first
    if exist('sites', 'var')
        disp('Reducing travel times')
        obj.traveltime = NaN(size(obj.time));
        for c=1:numel(sites)
            thissite = sites(c);
            thischanstr = thissite.channeltag.string();
            i = strcmp(obj.channelinfo, thischanstr)==1;
            obj.traveltime(i) = sites(c).traveltime;
            obj.time(i) = obj.time(i) - sites(c).traveltime/86400;
        end
    end
    % sort obj in ascending time order
    disp('Sorting')
    [detection_time,indices]=sort(obj.time);
    obj = obj.subset(indices);

    % First we begin with each detection and find how many detections occur in the
    % maxTimeDiff seconds that follow. This way each arrival is assigned a
    % weight equal to that number of detections. The point of this is to help us
    % identify the beginning of each event.
    fprintf('\nFinding how many detections within %e sec of each detection\n',maxTimeDiff)
    for c=1:numel(detection_time)
        associated_indices{c} = find(detection_time>=detection_time(c) & detection_time<=detection_time(c)+maxTimeDiff/86400);
        numdetections(c) = numel(associated_indices{c});
    end

    % remove decrementing series, e.g. replace sequence
    % like 6 5 4 3 2 1 with 6 0 0 0 0 0
    disp('Removing decrementing series')
    numdetections = [numdetections(1) diff(numdetections)+1];
    %numdetections(numdetections<=1)=0;

    % Now loop over numdetections and create events
    fprintf('Making events')
    eventnum = 0;
    duplicatecount = 0;
    for c=1:numel(numdetections)
        fprintf('.')
        if numdetections(c)>0
            fprintf('*')
            eventnum = eventnum + 1;
            otime(eventnum) = detection_time(c);
    %         lastDetectionTime(eventnum) = detection_time(c + numdetections(c) - 1);

            try
                eventdetobj = obj.subset([c : c + numdetections(c) - 1]);
            catch
                numdetections(c)
                c+numdetections(c)-1
                numel(numdetections)
                rethrow();
            end

            % remove the reduced time
            if exist('sites','var')
                for cc=1:numel(eventdetobj.time)
                    eventdetobj.time(cc) = eventdetobj.time(cc) + eventdetobj.traveltime(cc)/86400;
                end
            end

            % remove duplicate channels
            [uc,ia]  = unique(eventdetobj.channelinfo, 'stable');
            duplicatecount = duplicatecount + numel(eventdetobj.channelinfo) - numel(uc);
            eventdetobj = eventdetobj.subset(ia);

            % set the arrivalobj for this event
            arrivalobj{eventnum} = detection2arrival(eventdetobj);
            firstDetectionTime(eventnum) = min(eventdetobj.time);
            lastDetectionTime(eventnum) = max(eventdetobj.time);

        end

        if mod(c,30) == 0
            fprintf('\nProcessed %d out of %d\n',c, numel(numdetections));
        end
    end

    %% CREATE CATALOG
    if numel(otime)==0
        % no events
        catalogobj = Catalog();
        return
    end

    fprintf('\nCreating Catalog\n')
    olon=[];
    olat=[];
    if exist('source','var')
        olon = source.lon*ones(size(otime));
        olat = source.lat*ones(size(otime));
    end
    catalogobj = Catalog(otime, olon, olat, [], [], {}, {}, 'ontime', firstDetectionTime, 'offtime', lastDetectionTime);
    catalogobj.arrivals = arrivalobj;
    fprintf('%d detections were determined to be duplicates using a time window of %.1f seconds\n',duplicatecount, maxTimeDiff);

end

function arrivalobj = detection2arrival(detectionobj)
    if ~isa(detectionobj.channelinfo,'ChannelTag')
        ctag = ChannelTag(detectionobj.channelinfo);
    else
        ctag = detectionobj.channelinfo;
    end
    sta = get(ctag, 'station');
    chan = get(ctag, 'channel');
    arrivalobj = Arrival(cellstr(sta), cellstr(chan), detectionobj.time, ...
        cellstr(detectionobj.state), 'signal2noise', detectionobj.signal2noise);

end

% function result = alreadyHaveArrivalFromThisChannel(ctaglist, thisctag)
%   result = sum(cellfun(@(s) ~isempty(strfind(this
% end

author: geoscience-community-codes

repo: GISMO


text:
function test_bug2093

% WALLTIME 00:40:00
% MEM 3gb
% DEPENDENCY ft_read_header ft_read_data ft_read_event read_nex_header read_nex_data read_plexon_nex

% Compares output of same test NEX file from ft_read_header,
% ft_read_data, and ft_read_event using old (master branch) and new
% (nexhandling branch) versions of the FieldTrip code

testnexfile = dccnpath('/home/common/matlab/fieldtrip/data/test/original/lfp/plexon/p213parall.nex');
testmatfile = dccnpath('/home/common/matlab/fieldtrip/data/test/bug2093.mat');

%% first get output of old code
% if ~exist(testmatfile, 'file')
%   % run only once, this is to be done using FT version older than 9/9/2016,
%   % last commit before changes:  c8de94a3e13df7d10737381f1070a0d467c9bfbb
%   old.hdr = ft_read_header(testnexfile);
%   old.dat = ft_read_data(testnexfile);
%   old.evt = ft_read_event(testnexfile);
%
%   save(testmatfile, 'old')
% else
  % all future test executions with newer FT code
  load(testmatfile)
% end

%% then get output of new code
% new.hdr = ft_read_header(testnexfile);
% new.dat = ft_read_data(testnexfile);
new.evt = ft_read_event(testnexfile);

%% then compare the 2
% make sure the old header and data are identical, since I haven't changed them yet
% assert(isequal(old.hdr,new.hdr))
% assert(isequaln(old.dat,new.dat)) % data padded with NaNs should still be equal

%% make sure all old events are still present, despite adding more
assert(isequal(fieldnames(old.evt), fieldnames(new.evt)))
fnames = fieldnames(old.evt);
oldcell = cell(1,numel(fnames));
newcell = cell(1,numel(fnames));
mtnew = false(numel(new.evt),numel(fnames));
for fieldn = 1:numel(fnames)
  if ischar(old.evt(1).(fnames{fieldn})) && ...
      ischar(new.evt(1).(fnames{fieldn}))
    oldcell{fieldn} = {old.evt.(fnames{fieldn})}';
    newcell{fieldn} = {new.evt.(fnames{fieldn})}';
  else
    oldcell{fieldn} = vertcat(old.evt.(fnames{fieldn}));
    newcell{fieldn} = vertcat(new.evt.(fnames{fieldn}));
  end
  assert(isequal(class(oldcell{fieldn}), class(newcell{fieldn})))

  if ~isequal(size(oldcell{fieldn},1), numel(old.evt))
    error('This strategy won''t work.')
  end
  if ~isequal(size(newcell{fieldn},1), numel(new.evt))
    for evn = 1:numel(new.evt)
      mtnew(evn,fieldn) = isempty(new.evt(evn).(fnames{fieldn}));
    end
    if ~isequal(size(newcell{fieldn},1), sum(~mtnew(:,fieldn)))
      error('Empty fields don''t account for difference in event #s.')
    end
  end
end

oldevttable = struct2table(old.evt);
% Empty values in numeric variables cause the column to be represented as a
% cell-array when the structure is converted to a table, which is not
% supported by ismember.  The empty values should only be in the new event
% types I added support for, so they can be safely excluded from the table
% when verifying that all old event rows are contained in the new event
% data.
newrows = ~any(mtnew,2);
newevttable = table; %#ok<*AGROW>
for fieldn = 1:numel(fnames)
  if isequal(size(newcell{fieldn},1), numel(newrows))
    newevttable = [newevttable table(newcell{fieldn}(newrows,:), ...
      'VariableNames',fnames(fieldn))];
  elseif isequal(size(newcell{fieldn},1), sum(newrows))
    newevttable = [newevttable table(newcell{fieldn}, ...
      'VariableNames',fnames(fieldn))];
  else
    error('There may be different empty cells in different event fields.')
  end
end

assert(all(ismember(oldevttable,newevttable)))
disp('All old events found in new event output.')


text:
function [node,face,elem]=meshcylinders(c0, v, len, varargin)
%
% [node,face]=meshcylinders(c0, v, len, r,tsize,maxvol,ndiv)
%    or
% [node,face,elem]=meshacylinder(c0, v, len, r, tsize,maxvol,ndiv)
% [nplc,fplc]=meshacylinder(c0, v, len,r,0,0,ndiv);
%
% create the surface and (optionally) tetrahedral mesh of a 3D cylinder
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%   c0, cylinder list axis's starting point
%   v: directional vector of the cylinder
%   len: a scalar or a vector denoting the length of each
%        cylinder segment along the direction of v
%   tsize, maxvol, ndiv: please see the help for meshacylinder for details
%
% output:
%   node, face, elem: please see the help for meshacylinder for details
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

len=cumsum(len);
[ncyl,fcyl]=meshacylinder(c0,c0+v*len(1),varargin{:});

for i=2:length(len)
   [ncyl1,fcyl1]=meshacylinder(c0+v*len(i-1),c0+v*len(i),varargin{:});
   fcyl1=cellfun(@(x) {x{1}+size(ncyl,1),x{2}}, fcyl1, 'UniformOutput', false);
   ncyl=[ncyl; ncyl1];
   if(i==1)
       fcyl1=fcyl1(1:end-1);
   else
       fcyl1={fcyl1{1:end-2},fcyl1{end}};
   end
   fcyl={fcyl{:}, fcyl1{:}};
end

[ncyl,I,J]=unique(round(ncyl*1e10),'rows');
ncyl=ncyl*1e-10;
fcyl=cellfun(@(x) {J(x{1})',x{2}}, fcyl, 'UniformOutput', false);

tsize=varargin{2};
maxvol=varargin{3};

if(nargout==2 && tsize==0.0 && maxvol==0.0)
    node=ncyl;
    face=fcyl;
    return;
end
if(nargin==3)
    tsize=len/10;
end
if(nargin<5)
    maxvol=tsize*tsize*tsize;
end

centroid=cumsum([0 len(1:end-1)])+len/2;   % define the centroids of each cylinder segment
seeds=repmat(c0(:)',length(len),1)+repmat(v(:)',length(len),1).*repmat(centroid(:),1,3);
[node,elem,face]=surf2mesh(no,fc,min(no),max(no),1,maxvol,seeds,[],0);


text:
function pass = test_repeatedArithmetic(pref)

if ( nargin < 1 )
    pref = chebfunpref();
end
tol = 1e4*pref.cheb3Prefs.chebfun3eps;

% Add 50 functions together:
f = chebfun3(@(x,y,z) cos(x.*y.*z));
g = chebfun3(0);
for jj = 1:50
    g = g + f;
end
pass(1) = norm(g - 50*f ) < 10 * tol;

% Multiply 10 functions together:
f = chebfun3(@(x,y,z) cos(x.*y.*z));
g = f;
for jj = 1:10
    g = g .* f;
end
pass(2) = norm(g - f.^11) < tol;

% Multiply 10 functions together by calling the constructor in different
% ways:
f = chebfun3(@(x,y,z) cos(x.*y.*z));
g = f;
for jj = 1:10
    g = chebfun3(@(x,y,z) g(x,y,z) .* f(x,y,z));
end
pass(3) = norm(g - f.^11) < tol;

g = f;
for jj = 1:10
    g = chebfun3(@(x,y,z) g(x,y,z) .* f(x,y,z), 'fiberDim', 1);
end
pass(4) = norm(g - f.^11) < tol;

g = f;
for jj = 1:10
    g = chebfun3(@(x,y,z) g(x,y,z) .* f(x,y,z), 'fiberDim', 2);
end
pass(5) = norm(g - f.^11) < tol;

g = f;
for jj = 1:10
    g = chebfun3(@(x,y,z) g(x,y,z) .* f(x,y,z), 'fiberDim', 3);
end
pass(6) = norm(g - f.^11) < tol;

% Multiply 20 functions together:
f = chebfun3(@(x,y,z) sin(x.*y.*z));
g = f;
for jj = 1:20
    g = g .* f;
end
pass(7) = norm(g - f.^21) < tol;

% Multiply 20 functions together:
f = chebfun3(@(x,y,z) sin(x.*y.*z));
g = f;
for jj = 1:20
    g = chebfun3(@(x,y,z) g(x,y,z) .* f(x,y,z));
end
pass(8) = norm(g - f.^21) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function OUT = times(varargin)

funname = 'times';
narginchk(2,2)
OUT = builtincaller(funname,varargin{:});

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
% Philipp Berens
% CircStat: A Matlab Toolbox for Circular Statistics
% Submitted to Journal of Statistical Software
%
% Example 2
% An Application to Neuroscience
%
%
% In this example, we assess the orientation tuning properties of three
% neurons recorded from the primary visual cortex of awake macaques. the
% number of action potentials such neurons fire is modulated by the
% orientation of a visual stimulus such as an oriented grating.
%
% We thus consider two variables: The stimulus orientations ori spaced 22.5
% deg apart and the number of spikes w fired in response to each
% orientation of the stimulus.

%% part 1: load and plot data
clear
load neurodata

% orientation of bins -> convert two directions
ori = circ_axial(circ_ang2rad(ori),2);

% spacing of bins
dori = diff(ori(1:2));

% summed spikes per orientation
w;

% plot the activity of the three neurons
figure
for j = 1:3
  subplot(1,3,j)

  % compute and plot mean resultant vector length and direction
  mw = max(w(j,:));
  r = circ_r(ori,w(j,:),dori) * mw;
  phi = circ_mean(ori,w(j,:));
  hold on;
  zm = r*exp(i*phi');
  plot([0 real(zm)], [0, imag(zm)],'r','linewidth',1.5)

  % plot the tuning function of the three neurons
  polar([ori ori(1)], [w(j,:) w(j,1)],'k')

  % draw a unit circle
  zz = exp(i*linspace(0, 2*pi, 101)) * mw;
  plot(real(zz),imag(zz),'k:')
  plot([-mw mw], [0 0], 'k:', [0 0], [-mw mw], 'k:')

  formatSubplot(gca,'ax','square','box','off','lim',[-mw mw -mw mw])
  set(gca,'xtick',[])
  set(gca,'ytick',[])

end

%% part 2: descriptive statistics

stats = zeros(3,10);
for i=1:3

  spk = w(i,:);

  % circular mean angle
  stats(i,1) = circ_mean(ori,spk,2);

  % circular variance
  stats(i,2) = circ_var(ori,spk,dori,2);

  % circular standard deviation
  [stats(i,3) stats(i,4)] = circ_std(ori,spk,dori,2);

  % circular skewness
  [stats(i,5) stats(i,6)] = circ_skewness(ori,spk,2);

  % circular skewness
  [stats(i,7) stats(i,8)] = circ_kurtosis(ori,spk,2);

  % confidence limits on mean angle
  t = circ_confmean(ori,[],spk,dori,2);
  stats(i,9) = stats(i,1) + t;
  stats(i,10) = stats(i,1) - t;
end

% stats contains all data reported in table 1

%% part 3: inferential statistics

% A: tests for uniformity of distribution around the circle
% rejecting the null hypothesis allows us to assert that the neurons are
% indeed tuned to the orientation of the stimulus and fire preferentially
% at a particular orientation

uniform = zeros(3,2);
for i=1:3

  spk = w(i,:);

  % rayleigh test
  uniform(i,1) = circ_rtest(ori,spk,dori);

  % omnibus test
  uniform(i,2) = circ_otest(ori,[],spk);

  % rao's spacing test is not possible with binned data
end

% B: test for differences in preferred orientation between neurons

% differences between all groups
alpha = [ori ori ori];
idx = [ones(1,length(ori)) 2* ones(1,length(ori)) 3* ones(1,length(ori))];
spk = reshape(w',1,numel(w));

fprintf('TESTING FOR DIFFERENCES BETWEEN ANY CELLS\n')
circ_wwtest(alpha,idx,spk);

% all pairwise differences
for i=1:3
  for j=(i+1):3
    % differences between cells i and cell j
    alpha = [ori ori];
    idx = [i* ones(1,length(ori)) j* ones(1,length(ori))];
    spk = reshape(w([i j],:)',1,numel(w([i j],:)));

    fprintf('TESTING FOR DIFFERENCES BETWEEN CELLS %d AND %d\n',i,j)
    watson(i,j) = circ_wwtest(alpha,idx,spk); %#ok<AGROW>
  end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function pot = convert_to_pot(CPD, pot_type, domain, evidence)
% CONVERT_TO_POT Convert a tabular utility node to one or more potentials
% pot = convert_to_pot(CPD, pot_type, domain, evidence)

switch pot_type
 case 'u',
  sz = [CPD.sizes 1]; % the utility node itself has size 1
  pot = upot(domain, sz, 1*myones(sz), myreshape(CPD.T, sz));
 otherwise,
  error(['can''t convert a utility node to a ' pot_type ' potential']);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%  Footnote
%
%  ::
%
%     document.footnote('option_name',option_value);
%
%  Args:
%
%     varargin: arguments need to come in pairs
%
%        - **text** : footnote text
%
%

author: jmaih


text:
function [beta_median beta_std beta_lbound beta_ubound sigma_median sigma_t_median sigma_t_lbound sigma_t_ubound gamma_median]=stvol2estimates(beta_gibbs,sigma_gibbs,sigma_t_gibbs,gamma_gibbs,n,T,cband)

% compute the median, variance, and credibility intervals for the posterior distribution of beta
beta_median=quantile(beta_gibbs,0.5,2);
beta_std=std(beta_gibbs,0,2);
beta_lbound=quantile(beta_gibbs,(1-cband)/2,2);
beta_ubound=quantile(beta_gibbs,1-(1-cband)/2,2);

% compute the results for sigma (long-run value)
sigma_median=reshape(quantile(sigma_gibbs,0.5,2),n,n);

% compute the rsults for sigma (sample values)
sigma_t_median=cell(n,n);
sigma_t_lbound=cell(n,n);
sigma_t_ubound=cell(n,n);
% loop over periods and entries
for ii=1:T
   for jj=1:n
      for kk=1:jj
      sigma_t_median{jj,kk}(ii,1)=quantile(sigma_t_gibbs{ii,1}(jj,kk,:),0.5,3);
      sigma_t_lbound{jj,kk}(ii,1)=quantile(sigma_t_gibbs{ii,1}(jj,kk,:),(1-cband)/2,3);
      sigma_t_ubound{jj,kk}(ii,1)=quantile(sigma_t_gibbs{ii,1}(jj,kk,:),1-(1-cband)/2,3);
      end
   end
end

% compute the estimates for gamma
gamma_median=quantile(gamma_gibbs,0.5,1);

author: european-central-bank

repo: BEAR-toolbox


text:
function [M,MU,N,B,C,prm,iprm] = msns_pre(M,N,B,C)
%
%  Preprocessing of the system
%      .
%    M x  =  N x + B u
%                                                             (1)
%      y  =  C x,
%
%  where both M and N are REAL, SYMMETRIC and SPARSE. Moreover, M must be
%  positive definite and N must be negative definite.
%
%  The preprocessing consists of a double transformation of the state:
%
%    x <-- MU * P * x .
%
%  The first transformation with the permutation matrix P for bandwidth
%  reduction results in "overwriting" the system matrices as
%
%    M <-- P * M * P',  N <-- P * N * P',  B <-- P * B,  C <-- C * P'.
%
%  The bandwidth of the reordered matrices M and N is often much smaller
%  than that of the original matrices.
%
%  By the second transformation, the generalized system (1) is transformed
%  into a standard system
%      .
%      x  =  A x + B u
%                                                             (2)
%      y  =  C x.
%
%  To this end, the Cholesky factorization of M is computed: M = MU'*MU,
%  where MU is upper triangular. This results in:
%
%    A := inv(MU')*N*inv(MU),  B <-- inv(MU')*B,  C <-- C*inv(MU).
%
%  The matrix A, which is dense in general, is not formed explicitely.
%  It is implicitely given by N and MU.
%
%  Note that the systems (1) and (2) have an identical input-output
%  mapping.
%
%  Calling sequence:
%
%    [M,MU,N,B,C,prm,iprm] = msns_pre(M,N,B,C)
%
%  Input:
%
%    M, N      n-x-n system matrices;
%    B         n-x-m system matrix;
%    C         q-x-n system matrix.
%
%  Output:
%
%    M         permuted matrix M;
%    MU        Cholesky factor of (permuted) matrix M;
%    N         permuted matrix N;
%    B, C      transformed system matrices;
%    prm       the permutation that has been used in the first
%              transformation step;
%    iprm      the inverse permutation (needed to re-reorder certain data
%              in postprocessing).
%
%
%  LYAPACK 1.0 (Thilo Penzl, May 1999)

% Input data not completely checked!

if any(any(imag(M))) | any(any(imag(N))) | any(any(imag(B))) | ...
    any(any(imag(C)))
  disp('WARNING in ''msns_pre'': M, N, B, and C must be real matrices.');
  pause(10);
end

if norm(M-M','fro')~=0
  error('M is not symmetric!');
end

if norm(N-N','fro')~=0
  error('N is not symmetric!');
end

[prm,iprm] = lp_prm(M,N);

M = M(prm,prm);
N = N(prm,prm);

[MU,t] = chol(M);

if t~=0
  error('M is not (numerically) positive definite!');
end

if length(B)
  B = MU'\B(prm,:);
end

if length(C)
  C = C(:,prm)/MU;
end

author: Sable


text:
function optimizeFFT( dimension, dim, kernelSize, fft_planner_method, fftdimPerformed, lambdaCalib, flagZeropadding, sPrecision )
%OPTIMIZEFFT optimize the fft command
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

fprintf('Optimizing fft command\n');
if(strcmp(dimension,'2D'))
    if(dim(4) == 1)
        dimImg = [dim(1), dim(2), dim(5)]; % y-x-cha
        dimKernel = [dimImg(1)+kernelSize(1)-1, dimImg(2)+kernelSize(2)-1, dimImg(3)];
    else
        dimImg = [dim(4), dim(1), dim(2), dim(5)]; % t-y-x-cha
        dimKernel = [dimImg(1), dimImg(2)+kernelSize(1)-1, dimImg(3)+kernelSize(2)-1, dimImg(4)];
    end
elseif(strcmp(dimension,'3D'))
    dimImg = [dim(1), dim(3), dim(2), dim(5)]; % y-z-x-cha
    dimKernel = [dimImg(1)+kernelSize(1)-1, dimImg(2)+kernelSize(2)-1, dimImg(3)+kernelSize(3)-1, dimImg(4)];
elseif(strcmp(dimension,'4D'))
    dimImg = [dim(4), dim(1), dim(3), dim(2), dim(5)]; % t - y - z - x - cha
    dimKernel = [dimImg(1), dimImg(2)+kernelSize(1)-1, dimImg(3)+kernelSize(2)-1, dimImg(4)+kernelSize(3)-1, dimImg(5)];
elseif(strcmp(dimension,'5D'))
    dimImg = [dim(4), dim(1), dim(3), dim(2), dim(6), dim(5)]; % t - y - z - x - g - cha
    dimKernel = [dimImg(1), dimImg(2)+kernelSize(1)-1, dimImg(3)+kernelSize(2)-1, dimImg(4)+kernelSize(3)-1, dimImg(5), dimImg(6)];
end

% check fft along dimensions
% maximal fftdims = {1, 2, 3, [1 2], [1 3], [2 3], 1:3};
fftdims = {fftdimPerformed};
if(lambdaCalib > 0)
    fftdims{end+1} = 1:3; % due to corrKernelKspace -> convKernelImgspace
end

% check if wisdom file already exists
if(evalin('base','exist(''fftw_wisdom'',''var'')'))
    sizesCalculated = evalin('base', 'fftw_wisdom.sizesCalculated');
    doOptimization = false;
    if(size(sizesCalculated,3) ~= length(fftdims) || isemtpy(ismember(reshape(ipermute(sizesCalculated,[1 3 2]),2*size(sizesCalculated,3),size(sizesCalculated,2)),dimImg,'rows')) || isemtpy(ismember(reshape(ipermute(sizesCalculated,[1 3 2]),2*size(sizesCalculated,3),size(sizesCalculated,2)),dimKernel,'rows')))
        evalin('base', 'clear ''fftw_wisdom''');
        doOptimization = true;
    end
else
    doOptimization = true;
end
% if(exist([currpath,filesep,'utils',filesep,'general',filesep,'fftw_wisdom.mat'],'file'))
%     load([currpath,filesep,'utils',filesep,'general',filesep,'fftw_wisdom.mat']);
%     doOptimization = false;
%     if(size(sizesCalculated,3) ~= length(fftdims) || isemtpy(ismember(reshape(ipermute(sizesCalculated,[1 3 2]),2*size(sizesCalculated,3),size(sizesCalculated,2)),dimImg,'rows')) || isemtpy(ismember(reshape(ipermute(sizesCalculated,[1 3 2]),2*size(sizesCalculated,3),size(sizesCalculated,2)),dimKernel,'rows')))
%         delete([currpath,filesep,'utils',filesep,'general',filesep,'fftw_wisdom.mat']);
%         doOptimization = true;
%     end
% else
%     doOptimization = true;
% end

if(doOptimization)
    testImg = complex(randn(dimImg,sPrecision),randn(dimImg,sPrecision));
    if(flagZeropadding && lambdaCalib > 0)
        testKernel = complex(randn(dimKernel,sPrecision),randn(dimKernel,sPrecision));
    end

    sizesCalculated = zeros(2,length(dimImg),7);
    wisdom_str = cell(2,7);
    for i=1:length(fftdims)
        depth = sum(fftdims{i});
        if(length(fftdims{i}) > 1), depth = depth + 1; end;
        % optimize fft for image dimension
        fftw('planner',fft_planner_method);
        tmp = fftnshift(testImg,fftdims{i});
        wisdom_str{1,depth} = fftw('wisdom');
        sizesCalculated(1,:,depth) = dimImg;

        if(flagZeropadding && lambdaCalib > 0)
            % optimize fft for image kernel dimension
            fftw('planner',fft_planner_method);
            tmp = fftnshift(testKernel,fftdims{i});
            wisdom_str{2,depth} = fftw('wisdom');
            sizesCalculated(2,:,depth) = dimKernel;
        end
    end

    fftw_wisdom.sizesCalculated = sizesCalculated;
    fftw_wisdom.wisdom_str = wisdom_str;
    assignin('base', 'fftw_wisdom', fftw_wisdom);
%     save([currpath,filesep,'utils',filesep,'general',filesep,'fftw_wisdom.mat'], 'wisdom_str', 'sizesCalculated');
    clear 'testimg' 'testKernel';
end

author: thomaskuestner


text:
% Polynomial discrimination
% Section 8.6.2, Boyd & Vandenberghe "Convex Optimization"
% Original by Lieven Vandenberghe
% Adapted for CVX by Joelle Skaf - 10/23/05
% (a figure is generated)
%
% The goal is to find the polynomial of degree 4 on R^n that separates
% two sets of points {x_1,...,x_N} and {y_1,...,y_N}. We are trying to find
% the coefficients of an order-4-polynomial P(x) that would satisfy:
%           minimize    t
%               s.t.    P(x_i) <= t  for i = 1,...,N
%                       P(y_i) >= t   for i = 1,...,M

% Data generation
rand('state',0);
N = 100;
M = 120;

% The points X lie within a circle of radius 0.9, with a wedge of points
% near [1.1,0] removed. The points Y lie outside a circle of radius 1.1,
% with a wedge of points near [1.1,0] added. The wedges are precisely what
% makes the separation difficult and interesting.
X = 2 * rand(2,N) - 1;
X = X * diag(0.9*rand(1,N)./sqrt(sum(X.^2)));
Y = 2 * rand(2,M) - 1;
Y = Y * diag((1.1+rand(1,M))./sqrt(sum(Y.^2)));
d = sqrt(sum((X-[1.1;0]*ones(1,N)).^2));
Y = [ Y, X(:,d<0.9) ];
X = X(:,d>1);
N = size(X,2);
M = size(Y,2);

% Construct Vandermonde-style monomial matrices
p1   = [0,0,1,0,1,2,0,1,2,3,0,1,2,3,4]';
p2   = [0,1,1,2,2,2,3,3,3,3,4,4,4,4,4]'-p1;
np   = length(p1);
op   = ones(np,1);
monX = X(op,:) .^ p1(:,ones(1,N)) .* X(2*op,:) .^ p2(:,ones(1,N));
monY = Y(op,:) .^ p1(:,ones(1,M)) .* Y(2*op,:) .^ p2(:,ones(1,M));

% Solution via CVX
fprintf(1,'Finding the optimal polynomial of order 4 that separates the 2 classes...');

cvx_begin
    variables a(np) t(1)
    minimize ( t )
    a'*monX <= t;
    a'*monY >= -t;
    % For normalization purposes only
    norm(a) <= 1;
cvx_end

fprintf(1,'Done! \n');

% Displaying results
nopts = 2000;
angles = linspace(0,2*pi,nopts);
cont = zeros(2,nopts);
for i=1:nopts
   v = [cos(angles(i)); sin(angles(i))];
   l = 0;  u = 1;
   while ( u - l > 1e-3 )
      s = (u+l)/2;
      x = s * v;
      if a' * ( x(op,:) .^ p1 .* x(2*op) .^ p2 ) > 0,
          u = s;
      else
          l = s;
      end
   end;
   s = (u+l)/2;
   cont(:,i) = s*v;
end;

graph = plot(X(1,:),X(2,:),'o', Y(1,:), Y(2,:),'o', cont(1,:), cont(2,:), '-');
set(graph(2),'MarkerFaceColor',[0 0.5 0]);
title('Optimal order-4 polynomial that separates the 2 classes')
% print -deps min-deg-discr.eps

%%%% Dual infeasible ?????

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function g = NLP_gradFunction(x, Prob)
% g = NLP_gradFunction(x, Prob)
% From tomlab quickguide

if isempty(Prob.uP)
    alpha = 100;
else
    alpha = Prob.uP(1);
end

g = [-4 * alpha * x(1) * (x(2) - x(1) ^ 2) - 2 * (1 - x(1)); 2 * alpha * (x(2) - x(1) ^ 2)];

author: opencobra

repo: cobratoolbox


text:
%  Make Quotation
%
%  ::
%
%     document.quotation('option_name',option_value);
%
%  Args:
%
%     varargin: arguments need to come in pairs
%
%        - **log** : [char] the quote to insert
%        - **reference** : [char] reference of the quotation
%
%