text:
function [out]=snip(in,arg1,varargin)

%% SNIP FUNCTION  %2nd Version, 13.Sept.2013
% Snip something out of a vector or matrix.
%
% (C) for snip.m function implementation Nicolas Ummen, under BSD2 open source copyright license.
% Comments etc. to NicolasUmmen@web.de
%
% The function snip takes the following arguments:
%
% snip(input,argument_1,argument_2)
%
% %arg_1 = what to snip, arg_2 = for matrices, to define 'r'ow or 'c'olumn
%
% It then removes either the specified kinds of elements from the input, or
% removes the element at the stated position, then compresses the input by as
% much as was deleted out of it.
%
% LIST OF POSSIBLE ARGUMENT_1s:
% 'x'   -> remove all x's from the input vector, compress it afterwards
% '1'   -> remove all 1's from the input vector, compress it afterwards
% [all other comparable items to the above are possible]
% 1     -> remove THE FIRST element, move all one up, shorten vector by 1
% i     -> remove COMPLEX elements from input, compress it afterwards
% j     -> [same as above] - note that 'i' and 'j' would remove the characters i or j.
% nan   -> remove all nan's, compress it
% inf   -> remove all infs from the input vector, compress it afterwards
%
% etc.
%
% The rest of the possible arguments follows matlab notation, such as:
% last  -> snip off the last element
% :     -> snip from:to
%
% Note that 'partial' deletions, for example getting rid of elements from
% (2,3) to (2,6) in a 7x7 matrix = A, is much better accomplished by overwriting
% them directly as A(2,3:6) = 0; or similar. A compression of the data is
% not possible, as elements remain in e.g. (2,1:2)
%
% ALLOWED DATA TYPES: Everything except self-defined structures/objects.
%
% EXAMPLES:
%
% in = [2 1 3 NaN NaN 1]; snip(in,2)   = [2 3 NaN NaN 1]
% in = [2 1 3 NaN NaN 1]; snip(in,'1') = [2 3 NaN NaN]
% in = [2 1 3 NaN NaN 1]; snip(in,nan) = [2 1 3 1]
% in = [2 1 3 inf inf 1]; snip(in,inf) = [2 1 3 1]
% compl = complex(1,2); in = [2 1 3 compl compl 1]; snip(in,i) = [2 1 3 1]
% in = [2 1 3 NaN NaN 1]; snip(in,'last') = [2 1 3 NaN NaN]
% in = ['abcdefg']; snip(in,'c')   = abdefg
% in = ['abcdefg']; snip(in,'last') = abcdef
% in = ['abcdefg']; snip(in,2)     = acdefg
%
% Whether the input is a row or column vector doesn't matter, and the output
% is returned accordingly.
%
% Advanced vector, snip out a certain, continuous length:
%
% in = [2 1 3 NaN NaN 1]; snip(in,2:4) = [2 NaN 1]
%
% MATRIX EXAMPLES: either snip out completely rows/columns, or identify rows
% and columns that have the same elements throughout and delete these.
% Note: deleting all e.g. 3's therefore leaves 3's in mixed entry columns
% or rows. In some cases using snip.m column or row-wise on the matrix
% might help.
%
% in = [1 2 3; 3 3 3; 4 3 3]; snip(in,2,'r') = [1 2 3; 4 3 3]
% in = [1 2 3; 3 3 3; 4 3 3]; snip(in,2,'c') = [1 3; 3 3; 4 3]
% in = [1 2 3; 3 3 3; 4 3 3; 5 3 3; 6 3 3]; snip(in,2:4,'r') = [1 2 3; 6 3 3]
% in = [1 2 3; 3 3 3; 4 3 3]; snip(in,'1')     = [1 2 3; 3 3 3; 4 3 3]
% Note the above doesn't detect a complete row/column full of 1's...
%
% in = [1 2 3; 3 3 3; 4 3 3]; snip(in,'3') = [1 2; 4 3]
% While this one detected a row and a column of solely 3's. To overwrite
% all e.g. 3s in the above matrix, rather use find - or logical indexing.
% E.g. in(in == 3) = 0 for the result [1 2 0; 0 0 0; 4 0 0].
%
% in = [1 2 nan; nan nan nan; 4 nan nan]; snip(in,nan) = [1 2; 4 nan]
% in = ['abc'; 'acb'; 'abc']; snip(in,'a') = ['bc','cb','bc']
%
% etc as above for vectors
%
% For more complicated choices from a matrix; like the first, third and
% seventh column; use the notation A = A(:,[1,3,7])

% Change log:
% 2nd Version: Added infinity and complex number option on request.

%% (C) for snip.m function implementation Nicolas Ummen, under BSD2 open soure copyright license.
% Comments etc. to NicolasUmmen@web.de

% main function body

%scrutinize what the user wants, operate on most likely first
flagflipped = 0;
if isstruct(in)
   out = in;
   warning('snip function: structures cannot be snipped!')
   return
end

[z s] = size(in);
if z == 0
    out = in;
    warning('snip function: 0x0 input detected!')
    return
end

if z > s
    in = in';
    flagflipped = 1;
    [z s] = size(in);
end

%operate on numerical things first

%if it is a vector
if z == 1 && strcmpi(arg1,'last') %snip end command
    out = in(1,1:end-1);
    if flagflipped == 1
        out = out';
    end
    return
end

if z == 1 && isnan(arg1(1)) %sort nans out
    out = in(~isnan(in));
    if flagflipped == 1
        out = out';
    end
    return
end

if z == 1 && isinf(arg1(1)) %sort infs out
    out = in(~isinf(in));
    if flagflipped == 1
        out = out';
    end
    return
end

if z == 1 && ~isreal(arg1(1)) %sort complex numbers out
    out = in(~imag(in));
    if flagflipped == 1
        out = out';
    end
    return
end

if z == 1 && isnumeric(arg1) %snip out a specific place/length from vector
    if arg1(end) > s
        out = in;
        if flagflipped == 1
            out = out';
        end
        warning('snip function: the vector is too short for that snip!')
        return
    end
    out = [in(1:arg1(1)-1) in(arg1(end)+1:end)];
    if flagflipped == 1
        out = out';
    end
    return
end

if z == 1 && iscell(in) %if it is a cell vector
    if strcmpi(arg1,'last') %snip end command
        out = in{1,1:end-1};
        if flagflipped == 1
            out = out';
        end
        return
    end
    arg1 = str2num(arg1);
    if arg1(end) > s
        out = in;
        if flagflipped == 1
            out = out';
        end
        warning('snip function: the vector is too short for that snip!')
        return
    end
    out = [in{1:arg1(1)-1} in{arg1(end)+1:end}];
    if flagflipped == 1
        out = out';
    end
    return
end

if z == 1 %snip out all objects, like 'a', or '1' from a vector
    if length(str2num(arg1)) == 1
        arg1 = str2num(arg1);
    end
    out = in(in~=arg1);
    if flagflipped == 1
        out = out';
    end
    return
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%if it is a matrix

if z >= 1 && isnan(arg1(1)) %snip out all nans
    in = in(:,all(isnan(in),1) ~= 1);
    out = in(all(isnan(in),2) ~= 1,:);
    if flagflipped == 1
        out = out';
    end
    return
end

if z >= 1 && isinf(arg1(1)) %snip out all infs
    in = in(:,all(isinf(in),1) ~= 1);
    out = in(all(isinf(in),2) ~= 1,:);
    if flagflipped == 1
        out = out';
    end
    return
end

if z >= 1 && ~isreal(arg1(1)) %snip out complex numbers
    in = in(:,all(imag(in),1) ~= 1);
    out = in(all(imag(in),2) ~= 1,:);
    if flagflipped == 1
        out = out';
    end
    return
end

if z >= 1 && isempty(varargin) == 1 && ischar(arg1) && ischar(in(1,1)) %snip out all equal elements
    in = in(:,all(in==arg1,1) ~= 1);
    out = in(all(in==arg1,2) ~= 1,:);
    if flagflipped == 1
        out = out';
    end
    return
end

if z >= 1 && isempty(varargin) == 1 %snip out all equal elements
    arg1 = str2num(arg1);
    in = in(:,all(in==arg1,1) ~= 1);
    out = in(all(in==arg1,2) ~= 1,:);
    if flagflipped == 1
        out = out';
    end
    return
end

if z >= 1 && ischar(varargin{1}) %snip out specific row(s) or column(s)
    arg2 = varargin{1};
    if flagflipped == 1
        if strcmp(arg2,'c')
            arg2 = 'r';
        else
            if strcmp(arg2,'r')
            arg2 = 'c';
            end
        end
    end
    switch arg2
        case 'c'
            if arg1(end) > s
                out = in;
                if flagflipped == 1
                    out = out';
                end
                warning('snip function: the matrix is not large enough for this snip!')
                return
            end
            out = [in(:,1:arg1(1)-1) in(:,arg1(end)+1:end)];
        case 'r'
            if arg1(end) > z
                out = in;
                if flagflipped == 1
                    out = out';
                end
                warning('snip function: the matrix is not large enough for this snip!')
                return
            end
            out = [in(1:arg1(1)-1,:); in(arg1(end)+1:end,:)];
    end
    if flagflipped == 1
        out = out';
        return
    end
    return
end

%catch things that didn't go in the above cases

out = in;

if flagflipped == 1
    out = out';
end
warning('snip function: function completed, but no relevant case to execute occured!')
return

end


text:
function MS = createMeshSpherical1D(varargin)
% MeshStructure = createMeshSpherical1D(Nr, Lr)
% MeshStructure = createMeshSpherical1D(facelocationR)
% creates a uniform 1D mesh (1D axial symmetry):
% Nx is the number of cells in r (radial) direction
% Width is the domain radius in r direction
%
% SYNOPSIS:
%   MeshStructure = createMeshSpherical1D(Nr, Lr)
%
% PARAMETERS:
%   Nr: number of cells in the domain
%   Lr: domain radius
%
% RETURNS:
%   MeshStructure.
%                 dimensions=1.5 (1D axial symmetry)
%                 numbering: shows the indexes of cellsn from left to right
%                 dx: cell size (=Lr/Nx)
%                 cellcenters.x: location of each cell in the x direction
%                 facecenters.x: location of interface between cells in the
%                 r direction
%                 numberofcells: Nr
%
%
% EXAMPLE:
%   R = 1.0; % length of the domain
%   Nr = 10; % number of cells in the domain
%   m = createMeshCylindrical1D(Nr, R);
%   plot(m.cellcenters.x, ones(size(m.cellcenters.x)), 'o', ...
%        m.facecenters.x, ones(size(m.facecenters.x)), '-+');
%   legend('cell centers', 'face centers');
%
% SEE ALSO:
%     createMesh1D, createMesh2D, createMesh3D, ...
%     createMeshCylindrical2D, createCellVariable, createFaceVariable

% Written by Ali A. Eftekhari
% See the license file

if nargin==2
  % uniform 1D mesh
  Nr=varargin{1};
  R=varargin{2};
  % cell size is dx
  dx = R/Nr;
  CellSize.x= dx*ones(Nr+2,1);
  CellSize.y= [0.0];
  CellSize.z= [0.0];
  CellLocation.x= [1:Nr]'*dx-dx/2;
  CellLocation.y= [0.0];
  CellLocation.z= [0.0];
  FaceLocation.x= [0:Nr]'*dx;
  FaceLocation.y= [0.0];
  FaceLocation.z= [0.0];
elseif nargin==1
  % nonuniform 1D mesh
  facelocationR=varargin{1};
  n=size(facelocationR);
  if n(1)==1
      facelocationR=facelocationR';
  end
  Nr = length(facelocationR)-1;
  CellSize.x= [facelocationR(2)-facelocationR(1); ...
    facelocationR(2:end)-facelocationR(1:end-1); ...
    facelocationR(end)-facelocationR(end-1)];
  CellSize.y= [0.0];
  CellSize.z= [0.0];
  CellLocation.x= 0.5*(facelocationR(2:end)+facelocationR(1:end-1));
  CellLocation.y= [0.0];
  CellLocation.z= [0.0];
  FaceLocation.x= facelocationR;
  FaceLocation.y= [0.0];
  FaceLocation.z= [0.0];
end

MS=MeshStructure(1.8, ...
  [Nr,1], ...
  CellSize, ...
  CellLocation, ...
  FaceLocation, ...
  [1], ...
  [1]);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function model = gpCovGradsTest(model)

% GPCOVGRADSTEST Test the gradients of the likelihood wrt the covariance.
% FORMAT
% DESC tests the gradients of the covariance to ensure they are
% correct.
% ARG model : the model to be tested.
% RETURN model : the model that was tested.
%
% SEEALSO : gpCreate, gpCovGrads
%
% COPYRIGHT : Neil D. Lawrence, 2006, 2009

% GP

% WARNING --- this isn't testing g_Lambda in gpCovGrads

changeVal = 1e-6;
switch model.approx
 case 'ftc'

 case {'dtc', 'dtcvar', 'fitc', 'pitc'}
  for i =1 :size(model.K_uu, 1)
    for j=1:i
      origK = model.K_uu(i, j);
      model.K_uu(i, j) = origK + changeVal;
      model.K_uu(j, i) = model.K_uu(i, j);
      [model.invK_uu, U] = pdinv(model.K_uu);
      model.logDetK_uu = logdet(model.K_uu, U);
      model = gpUpdateAD(model);
      objPlus = gpLogLikelihood(model);
      model.K_uu(i, j) = origK - changeVal;
      model.K_uu(j, i) = model.K_uu(i, j);
      [model.invK_uu, U] = pdinv(model.K_uu);
      model.logDetK_uu = logdet(model.K_uu, U);
      model = gpUpdateAD(model);
      objMinus = gpLogLikelihood(model);
      diffsK_uu(i, j) = (objPlus - objMinus)/(2*changeVal);
      diffsK_uu(j, i) = diffsK_uu(i, j);
      model.K_uu(i, j) = origK;
      model.K_uu(j, i) = origK;
      [model.invK_uu, U] = pdinv(model.K_uu);
      model.logDetK_uu = logdet(model.K_uu, U);
      model = gpUpdateAD(model);
    end
  end
  for i=1:size(model.K_uf, 1)
    for j=1:size(model.K_uf, 2)
      origK = model.K_uf(i, j);
      model.K_uf(i, j) = origK + changeVal;
      model = gpUpdateAD(model);
      objPlus = gpLogLikelihood(model);
      model.K_uf(i, j) = origK - changeVal;
      model = gpUpdateAD(model);
      objMinus = gpLogLikelihood(model);
      diffsK_uf(i, j) = (objPlus - objMinus)/(2*changeVal);
      model.K_uf(i, j) = origK;
      model = gpUpdateAD(model);
    end
  end

  [gK_uu, gK_uf, g_Lambda] = gpCovGrads(model, model.m);

  gK_uuMaxDiff = max(max(abs(2*(gK_uu-diag(diag(gK_uu))) ...
                             + diag(diag(gK_uu)) ...
                             - diffsK_uu)));
  gK_ufMaxDiff = max(max(abs(gK_uf - diffsK_uf)));

  fprintf('K_uu grad max diff %2.4f\n', gK_uuMaxDiff);
  if gK_uuMaxDiff > 1e-4
    disp(2*(gK_uu-diag(diag(gK_uu))) ...
         + diag(diag(gK_uu)) ...
         - diffsK_uu);
  end
  fprintf('K_uf grad max diff %2.4f\n', gK_ufMaxDiff);
  if gK_ufMaxDiff > 1e-4
    disp(gK_uf - diffsK_uf)
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef ReadingFilesTests < handle & matlab.unittest.TestCase

    properties (TestParameter)
        readingTests = {'testReadingGmsh'}
    end

    methods (Test, TestTags = {'ReadingFiles', 'Fast', 'Gmsh'})

        function testReading(testCase, readingTests)
            s.computerType    = 'GMSH';
            s.testName         = readingTests;
            s.variablesToStore = {'connec', 'coord', 'isElemInThisSet', ...
                                  'masterSlave', 'corners'};
            test = PrecomputedVariableTest(s);
            err = test.computeError();
            tol = 1e-13;
            testCase.verifyLessThanOrEqual(err, tol)
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function plotCortmagResults(cortMag, showPredicted)
% function plotCortmagResults(cortMag, [showPredicted])
%
if(~exist('showPredicted','var')) showPredicted = 1; end

nROIs = length(cortMag.bins);
pX = ceil(sqrt(nROIs));
pY = ceil(nROIs/pX);
xLower = 0;
xUpper = 100;
yLower = -pi;
yUpper = 4*pi;
for ii=1:nROIs
   mnPh = complexPh2PositiveRad(cortMag.meanPh{ii});
   %mnPh = unwrapPhases(mnPh);
   % plot unshifted data
   uDist = cortMag.corticalDist{ii}-cortMag.distanceShift(ii);
   figure(100);
   subplot(pX,pY,ii);
   plot(uDist,mnPh,'-o');
   set(gca,'xlim',[xLower xUpper]);
   set(gca,'ylim',[yLower yUpper]);
   str = sprintf('ROI %.0f',ii);
   title(str);
%    figure(99);
%    hold on;
%    plot(uDist,(mnPh-cortMag.fovealPhase)/(2*pi)*cortMag.stimulusRadius,'o');
%    hold off;
   %pause
end

figure(101);
%newGraphWin;
symbolString = 'b.';
errorbar(cortMag.allCorticalDist10deg, cortMag.allStimDeg, cortMag.allStimDegSE, symbolString);

% The predicted exponential CMF curve, when the distances have
% been adjusted so that 0 means 10 deg, is:
%
d = sort(cortMag.allCorticalDist10deg);

if showPredicted
    predictedDeg = exp(d*cortMag.fitParms.dScale + log(10));
    hold on; plot(d,predictedDeg,'k-.'); hold off
    cortMag.predictedDeg = predictedDeg;
end

% Guesses about the window parameters
%
xLower = min(cortMag.allCorticalDist10deg);
xUpper = max(cortMag.allCorticalDist10deg);
%xUpper = 15;
yLower = -4;
yUpper = 1.1*cortMag.stimulusRadius;

set(gca,'xlim',[xLower xUpper]);
set(gca,'ylim',[yLower yUpper]);
set(gca,'xtick',round([xLower:5:xUpper]),'ytick',[yLower:4:yUpper])

% Place the resulting data in the figure for later plotting
set(gca,'UserData',cortMag);
disp('Data stored in plot- retrieve it with: cortMag=get(gca,''UserData'')');

subject = eval('cortMag.subject','');
dataType = eval('cortMag.dataType','');
flatDir = eval('cortMag.subdir','');
title([subject,' (',cortMag.hemisphere,',',dataType,',',flatDir,'): ',...
      'dScale=',num2str(cortMag.fitParms.dScale),'  ',...
      'dShift=',num2str(cortMag.fitParms.dShift),'  ',...
      'foveaPh=',num2str(cortMag.fitParms.fovealPhase)]);
%plot(cortMag.allCorticalDist, complexPh2PositiveRad(cortMag.allMeanPh),'x');
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function A = sim_genVARModelFromEq(expr,morder)
%
% Generate VAR model coefficient matrix from text-based system of equations
% The model can then be realized using the tvarsim() function modified from
% ARfit [2].
% If the expr specification requests time-varying coefficients, then this
% function will return inline function placeholders for coefficients, which
% can be evaluated to actual coefficient matrices by sim_genTVARcoeffs().
%
% Inputs:
%
%   expr:       A cell vector containing each equation as a string (one
%               equation per cell element). See Examples for format.
%   morder:     The model order
%
% Outputs:
%
%   A:          VAR[p] model coefficients in format A=[A1,A2, ... Ap] where
%               p = morder. Ai is the M x M coefficient matrix for lag i.
%
% Example1: generate a static VAR[3] model from a (text-based) system of equations
%
% expr = { ...
%     'x1(t) = 0.9*x1(t-1)  + 0.3*x2(t-2)  + e1(t)' ...
%     'x2(t) = 1.3*x2(t-1)  + -0.8*x2(t-2) + e2(t)' ...
%     'x3(t) = 0.3*x1(t-2)  + 0.6*x2(t-1)  + e3(t)' ...
%     'x4(t) = -0.7*x4(t-3) + -0.7*x1(t-3) + 0.3*x5(t-3) + e4(t)' ...
%     'x5(t) = 1*x5(t-1)    + -0.4*x5(t-2) + 0.3*x4(t-2) + e5(t)' ...
%     };
% A = sim_genVARModelFromEq(expr,3)
% A =
%     0.9000         0         0         0         0         0    0.3000         0         0         0         0         0         0         0         0
%          0    1.3000         0         0         0         0   -0.8000         0         0         0         0         0         0         0         0
%          0    0.6000         0         0         0    0.3000         0         0         0         0         0         0         0         0         0
%          0         0         0         0         0         0         0         0         0         0   -0.7000         0         0   -0.7000    0.3000
%          0         0         0         0    1.0000         0         0         0    0.3000   -0.4000         0         0         0         0         0
%
% % Now generate some data from the model (this requires arsim.m or tvarsim.m)
%
% M   = size(A,1);        % number of variables
% Nl  = 1000;             % length of each trial
% Nr  = 100;              % number of trials
% C = eye(M);             % specify the covariance matrix
% data = zeros(M,Nl,Nr);
% for tr=1:Nr             % simulate data from VAR model
%    data(:,:,tr) = arsim(zeros(1,M),A,C,Nl)';
% end
% eegplot(data,'srate',1); % visualize the simulated data
%
%
% Example 2: generate a time-varying VAR[2] model
%
%
% See Also: tvarsim(), est_fitMVAR(), sim_genTVARcoeffs()
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual. Section 6.5.1
%   Available at: http://www.sccn.ucsd.edu/wiki/Sift
% [2] Schneider T, Neumaier A (2001) Algorithm 808: ARfit---a matlab package
%   for the estimation of parameters and eigenmodes of multivariate
%   autoregressive models. ACM Transactions on Mathematical Software 27:58-65
%   http://www.gps.caltech.edu/~tapio/arfit/
%
% Author: Tim Mullen 2010, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

Nvars = length(expr);

if nargin>1
    % create cell array of zeros
    A = repmat({0},Nvars,Nvars*morder);
end

fnstr ='{[\S]+}';

for i = 1:Nvars
    eq = expr{i};

    % Identify all inline functions
    % Specifically, we extract any string f(x) encapsulated by
    % curly-brackets '{f(x)}'
    [matchstr] = regexp(eq,'{[\S]+}','match');

    % replace each function with an identifier '$i' designating that the
    % ith equation will go here
    for fi=1:length(matchstr)
        eq = strtrim(strrep(eq,matchstr{fi},sprintf('$%d',fi)));
    end

    % remove all whitespace
    eq(isspace(eq))=[];
%
%     for fi=1:length(matchstr)
%         eq = [eq(1:matchstart(fi)-1) sprintf('$%d',fi) eq(matchend(fi)+1:end)];
%     end

    terms = regexp(eq,'([=+*])','split');

    row = str2double(regexp(terms{1},'(\d+)','match'));
    for k = 2:2:length(terms)-1
        vars = regexp(terms{k+1},'(\d+)','match');
        vars = str2num(char(vars));
        col = vars(1)+(vars(2)-1)*Nvars;

        if strcmpi(terms{k}(1),'$')
            % term is a function expression
            A{row,col} = matchstr{str2double(terms{k}(2))}(2:end-1);
        elseif ~isnan(str2double(terms{k}))
            % term is a number
            A{row,col} = str2double(terms{k});
        else
            % term is an expression to be directly evaluated
            A{row,col} = eval(terms{k});
        end
    end
end

% if A is all numeric, convert to standard array
if all(cellfun(@(x)isnumeric(x),A))
    A = cell2mat(A);
end

repo: REST


text:
% Distortionless data hiding based on integer wavelet transform (Watermark Embedding)
% Guorong Xuan   Jiang Zhu   Jidong Chen   Shi, Y.Q.   Zhicheng Ni   Wei Su
% Department of Computer Science, Tongji University, Shanghai
%
% This paper appears in: Electronics Letters
% Publication Date: 5 Dec 2002
% Volume: 38,  Issue: 25
% On page(s): 1646- 1648
% ISSN: 0013-5194

% August 12, 2007, 4:37pm
% updated August 19, 2007, 5:45pm
% updated 13 December, 2008
% By: Asad (asad_82@yahoo.com)

% NOTE: In case of error change the bit plane number as embedding may not be possible using that bitplane.
% Also in order to achieve more bpp, bit plane number can be changed.

clear all;
close all;

% define constants
BITPLANE_NUMBER = 4; %4 % specifies the bitplane number to embedd watermark
ERROR_NUM = 16; %16 % prevents this gray value from being used in embedding
WM_SIZE = 300;
% select the bitplane number and then the corresponding error num from the
% table below
% BITPLANE_NUMBER = [1     2   3   4  5  6  7  8]
% ERROR_NUM       = [128  64  32  16  8  4  2  1]
disp('------------------------ Embedding -------------------------------');

% read image and convert to gray scale if necessary
originalImage = imresize(imread('lena.tif'),[512 512],'bicubic');
if (isrgb(originalImage))
    originalImage = rgb2gray(originalImage);
end

% STEP 1: perform necessary preprocessing
%[valCount,X] = imhist(originalImage);

% STEP 2: perfom wavelet decomposition
LS = liftwave('cdf2.2','Int2Int');
[CA,CH,CV,CD] = lwt2(double(originalImage),LS);

% STEP 3: construct binary images from 5th bit of CH, CV and CD
for i=1:size(CH,1)
    for j=1:size(CH,2)
        % for constructing binary image using CH
        binSeq = dec2bin(abs(CH(i,j)),8);
        if binSeq(BITPLANE_NUMBER) == '1'
            bICH5(i,j) = 2;
        else
            bICH5(i,j) = 1;
        end
        % for constructing binary image using CV
        binSeq = dec2bin(abs(CV(i,j)),8);
        if binSeq(BITPLANE_NUMBER) == '1'
            bICV5(i,j) = 2;
        else
            bICV5(i,j) = 1;
        end
        % for constructing binary image using CD
        binSeq = dec2bin(abs(CD(i,j)),8);
        if binSeq(BITPLANE_NUMBER) == '1'
            bICD5(i,j) = 2;
        else
            bICD5(i,j) = 1;
        end
    end
end

% STEP 4a: Compress data in 5th Bit plane of CH
% find how many times 1 occurs in the level 5 binary bICH5 (horizontal)
[xx,yy,val] = find(bICH5 == 1);
totalVals = size(CH,1)*size(CH,2);

% code the sequence using arithmetic coding
count1 = [round((size(xx,1)/totalVals)*100) round(((totalVals - size(xx,1))/totalVals) * 100)];
seq1 = reshape(bICH5,1,size(CH,1)*size(CH,2));
arcCH5 = arithenco(seq1,count1);
str = sprintf('Original CH Bits Length = %d ------ Compressed CH Bits Length = %d',size(seq1,2),size(arcCH5,2));
disp(str);

% STEP 4b: Compress data in 5th Bit plane of CV
% find how many times 1 occurs in the level 5 binary bICV5 (vertical)
[xx,yy,val] = find(bICV5 == 1);
totalVals = size(CV,1)*size(CV,2);

% code the sequence using arithmetic coding
count2 = [round((size(xx,1)/totalVals)*100) round(((totalVals - size(xx,1))/totalVals) * 100)];
seq2 = reshape(bICV5,1,size(CV,1)*size(CV,2));
arcCV5 = arithenco(seq2,count2);
str = sprintf('Original CV Bits Length = %d ------ Compressed CV Bits Length = %d',size(seq2,2),size(arcCV5,2));
disp(str);

% STEP 4c: Compress data in 5th Bit plane of CD
% find how many times 1 occurs in the level 5 binary bICD5 (diagonal)
[xx,yy,val] = find(bICD5 == 1);
totalVals = size(CD,1)*size(CD,2);

% code the sequence using arithmetic coding
count3 = [round((size(xx,1)/totalVals)*100) round(((totalVals - size(xx,1))/totalVals) * 100)];
seq3 = reshape(bICD5,1,size(CD,1)*size(CD,2));
arcCD5 = arithenco(seq3,count3);
str = sprintf('Original CD Bits Length = %d ------ Compressed CD Bits Length = %d',size(seq3,2),size(arcCD5,2));
disp(str);

% STEP 5: read the watermark and reshape it for insertion
watermark = imresize(rgb2gray(imread('logo_im.jpg')),[WM_SIZE WM_SIZE],'bicubic');
%watermark = imresize(DMlcd5,[WM_SIZE WM_SIZE],'bicubic');
watermark = im2bw(watermark);
%figure,imshow(watermark,[]),title('Watermark');
watermark = reshape(watermark,1,WM_SIZE*WM_SIZE);

% STEP 6: Insert the watermark and compressed data into the image
% compute length of data to insert
dataLength = size(watermark,2) + size(arcCH5,2) + size(arcCV5,2) + size(arcCD5,2) + 6*8 + 4*16;%(2*8)*3;
available = size(CH,1)*size(CH,2) + size(CV,1)*size(CV,2) + size(CD,1)*size(CD,2);

if available < dataLength
    disp('Data to Embedd must be less than available limit.');
end

% allocate memory and initialize
embedData = zeros(1,dataLength);
% insert the header information
str1 = dec2bin(count1(1,1),8); str2 = dec2bin(count1(1,2),8);
header = strcat(str1,str2);
str1 = dec2bin(count2(1,1),8); str2 = dec2bin(count2(1,2),8);
header = strcat(header,str1,str2);
str1 = dec2bin(count3(1,1),8); str2 = dec2bin(count3(1,2),8);
header = strcat(header,str1,str2);
% insert length of each sequence for decoding
str1 = dec2bin(size(arcCH5,2),16); str2 = dec2bin(size(arcCV5,2),16);
str3 = dec2bin(size(arcCD5,2),16); str4 = dec2bin(size(watermark,2),32);
header = strcat(header,str1,str2,str3,str4);

for i=1:size(header,2)
    if header(1,i) == '1'
        embedData(1,i) = 1;
    else
        embedData(1,i) = 0;
    end
end
% get the header length
HL = size(header,2);
% concatenate data to get a single compressed data vector
embedData(1,HL+1:HL+size(arcCH5,2)) = arcCH5(1,:);
embedData(1,HL+size(arcCH5,2)+1:(HL+size(arcCH5,2) + size(arcCV5,2))) = arcCV5(1,:);
embedData(1,HL+size(arcCH5,2)+ size(arcCV5,2)+1:(HL + size(arcCH5,2) + size(arcCV5,2) + size(arcCD5,2))) = arcCD5(1,:);
embedData(1,HL+size(arcCH5,2)+ size(arcCV5,2)+size(arcCD5,2)+1:(HL+size(arcCH5,2) + size(arcCV5,2) + size(arcCD5,2)+ size(watermark,2))) = watermark(1,:);
originalBitsLength = HL + size(arcCH5,2)+ size(arcCV5,2)+size(arcCD5,2);

% embedd compressed data, watermark and header information
index = 1; % counter for obtaining 8 bit chunks
brk = 0;
for x=1:size(CH,1)
    for y=1:size(CH,2)
        if CH(x,y) ~= -ERROR_NUM
            neg = 0;
            if CH(x,y) < 0
                neg = 1;
            end
            binSeq = dec2bin(abs(CH(x,y)),8);
            if embedData(1,index) == 1
                binSeq(BITPLANE_NUMBER) = '1';
            else
                binSeq(BITPLANE_NUMBER) = '0';
            end
            num = bin2dec(binSeq);
            if neg == 1
                CH(x,y) = num * -1;
            else
                CH(x,y) = num;
            end

            % break the loop if all watermark bits are embedded
             if index < size(embedData,2)
                index = index + 1;
             else
                 brk = 1; break;
             end
        end
    end
    if brk == 1
        break;
    end
end

if index < size(embedData,2)
    for x=1:size(CV,1)
        for y=1:size(CV,2)
            if CV(x,y) ~= -ERROR_NUM
                neg = 0;
                if CV(x,y) < 0
                    neg = 1;
                end
                binSeq = dec2bin(abs(CV(x,y)),8);
                if embedData(1,index) == 1
                    binSeq(BITPLANE_NUMBER) = '1';
                else
                    binSeq(BITPLANE_NUMBER) = '0';
                end
                num = bin2dec(binSeq);
                if neg == 1
                    CV(x,y) = num * -1;
                else
                    CV(x,y) = num;
                end

                % break the loop if all watermark bits are embedded
                 if index < size(embedData,2)
                   index = index + 1;
                 else
                     brk = 1; break;
                 end
            end
        end
        if brk == 1
            break;
        end
    end
 end

if index < size(embedData,2)
    for x=1:size(CD,1)
        for y=1:size(CD,2)
            if CD(x,y) ~= -ERROR_NUM
                neg = 0;
                if CD(x,y) < 0
                    neg = 1;
                end
                binSeq = dec2bin(abs(CD(x,y)),8);
                if embedData(1,index) == 1
                    binSeq(BITPLANE_NUMBER) = '1';
                else
                    binSeq(BITPLANE_NUMBER) = '0';
                end
                num = bin2dec(binSeq);
                if neg == 1
                    CD(x,y) = num * -1;
                else
                    CD(x,y) = num;
                end

                % break the loop if all watermark bits are embedded
                if index < size(embedData,2)
                    index = index + 1;
                 else
                     brk = 1; break;
                end
            end
        end
        if brk == 1
            break;
        end
    end
 end

% compute inverse integer wavelet transform
watermarkedImage = ilwt2(CA,CH,CV,CD,LS);
% [tx,ty,tval] = find(watermarkedImage > 256);
% if ~isempty(tx)
%     disp('Watermarked Image values greater than 256');
% end
% [tx,ty,tval] = find(watermarkedImage <= 0);
% if ~isempty(tx)
%     disp('Watermarked Image values less than equal to 0');
% end

imwrite(uint8(watermarkedImage),'Watermarked Image.bmp','bmp');

str = sprintf('Payload(bpp) = %f -- Embedded Data(Header+Original Bits+Watermark) = %d bits -- Watermark Length = %d bits',(dataLength-originalBitsLength)/(4*(size(CH,1)*size(CH,1))),dataLength+HL,size(watermark,2));
disp(str);

[PSNR_OUT,Z] = psnr(originalImage,watermarkedImage);
str = sprintf('PSNR = %f',PSNR_OUT);
disp(str);

figure,imshow(originalImage,[]),title('Original Image');
figure,imshow(watermarkedImage,[]),title('Watermarked Image');

save WatermarkInfo watermarkedImage watermark BITPLANE_NUMBER WM_SIZE ERROR_NUM originalImage;

% WMI = imread('Watermarked Image.bmp');
% difference = double(watermarkedImage) - double(WMI);
% figure,imshow(difference,[]),title('Difference Image');

author: Sable

repo: mcbench-benchmarks


text:
function lab = mineigK(x,K)

% lab = mineigK(x,K)
%
% MINEIGK  Computes the minimum eigenvalue of a vector x with respect to a
%          self-dual homogenous cone K.
%
% See also sedumi, mat, vec, eyeK.

% New function by Michael C. Grant
% Copyright (C) 2013 Michael C. Grant.

xi = 0;
if isfield(K,'f')
    xi = xi + K.f;
end
if isfield(K,'l') && K.l > 0
    lab = min(x(xi+1:xi+K.l));
    xi = xi+K.l;
else
    lab = Inf;
end
if isfield(K,'q') && ~isempty(K.q)
    scl = sqrt(0.5);
    for k = 1:length(K.q)
        kk = K.q(k);
        lab = min(lab,scl*(x(xi+1)-norm(x(xi+2:xi+kk))));
        xi = xi + kk;
    end
end
if isfield(K,'r') && ~isempty(K.r)
    % This is a simpler formula than the one found in eigK.c. In theory
    % there could be cancellation error in the smaller eigenvalue. But
    % the rotated Lorentz vector is not used internally where this
    % cancellation error might matter.
    for k = 1:length(K.r)
        kk = K.r(k);
        x1 = xx(xi+1);
        x2 = xx(xi+2);
        lab = min(lab,0.5*(x1+x2-norm([x1-x2;2*x(xi+3:xi+kk)])));
    end
end
if isfield(K,'s') && ~isempty(K.s)
    Ks = K.s;
    Kq = K.s .* K.s;
    nc = length(Ks);
    OPTS.disp=0;
    % When used internally, Hermitian terms are broken apart into real and
    % imaginary halves, so we need to catch this.
    if isfield(K,'rsdpN')
        nr = K.rsdpN;
    else
        nr = nc;
    end
    for i = 1 : nc
        ki = Ks(i);
        qi = Kq(i);
        XX = x(xi+1:xi+qi); xi=xi+qi;
        if i > nr
            XX = XX + 1i*x(xi+1:xi+qi); xi=xi+qi;
        end
        XX = reshape(XX,ki,ki);
        XX = XX + XX';
        if ki > 500
            lab=min(lab,0.5*eigs(XX,1,'SA',OPTS));
        else
            lab=min(lab,0.5*min(eig(XX)));
        end
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function xyzl_header_print ( point_num, line_num, line_data_num )

%*****************************************************************************80
%
%% XYZL_HEADER_PRINT prints the header of an XYZL file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer LINE_NUM, the number of lines.
%
%    Input, integer LINE_DATA_NUM, the number of line items.
%
  fprintf ( 1, ' \n');
  fprintf ( 1, '  Number of points =     %d\n', point_num );
  fprintf ( 1, '  Number of lines =      %d\n', line_num );
  fprintf ( 1, '  Number of line items = %d\n', line_data_num );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function filePath = BF_WriteTempFile(dataVector,whatPrecision)
% BF_WriteTempFile   Write a temporary file in the system temporary directory
%
%---INPUTS:
% dataVector, a vector of data to write to file.

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

if nargin < 2
    % precision (number of significant digits) to write
    whatPrecision = 7;
end

% Filename is the tempname (which will be unique):
filePath = tempname;

% Write the file:
dlmwrite(filePath,dataVector,'precision',whatPrecision);

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%
% lines = my_cat(file_name)
%
% Simulate the function of "cat" in Unix systems.
% The input is a file name.
% The output is a cell array, each cell contains one line of the file
% content.
% Author: Xiong Xiao,
%         School of Computer Engineering, Nanyang Technological University
% Date Created: 23 Oct, 2008
% Last Modified: 23 Oct, 2008
%
function lines = my_cat(file_name)

FID = fopen(file_name, 'r');

lines = read_line_by_line(FID);

fclose(FID);

author: singaxiong

repo: SignalGraph


text:
function L = log_prior(CPD)
% LOG_PRIOR Return log P(theta) for a hhmm CPD
% L = log_prior(CPD)

L = log_prior(CPD.sub_CPD_trans);
if ~isempty(CPD.sub_CPD_start)
  L = L + log_prior(CPD.sub_CPD_start);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [vals,coeffs]=asinDerivVal(z,n,coeffs)
%%ASINDERIVVAL Return the value of the nth derivative of asin(x) with
%              respect to x  evaluated at x=z for real or imaginary z.
%
%INPUTS: z A matrix of real or complex values at which the nth derivative
%          of the arcsine function is desired. If an empty matrix is
%          passed, then just coeffs will be returned.
%        n The number of derivatives to take. n>=0.
%   coeffs The computation of the derivatives involves determining the
%          structure of a polynomial. If this function has been run before
%          for a given n value, then the polynomial can be passed back and
%          is not determined again. Otherwise, this can be omitted or an
%          empty matrix can be passed. This only makes a difference for
%          n>=12.
%
%OUTPUTS: vals The value of the nth derivative of the arcsine function
%              taken at all of the points in z. vals has the same
%              dimensions as z.
%       coeffs A vector of polynomial coefficients that can be passed back
%              to this function when evaluating with the same n to speed it
%              up.
%
%For n=0, vals=asin(z). Subsequent derivatives are computed essentially
%using the chain rule and combining terms to get the appropriate
%polynomial. All derivatives are of the form
%polynomial*sqrt(1-z(:).^2).^(-k). For n<12, the coefficients are tabulated
%to make the function faster.
%
%May 2018 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(n==0)
    vals=asin(z);
    coeffs=[];
    return;
end

if(nargin<3||isempty(coeffs))
    if(n<12)
        switch(n)
            case 1
                coeffs=1;
            case 2
                coeffs=[1;0];
            case 3
                coeffs=[2;0;1];
            case 4
                coeffs=[6;0;9;0];
            case 5
                coeffs=[24;0;72;0;9];
            case 6
                coeffs=[120;0;600;0;225;0];
            case 7
                coeffs=[720;0;5400;0;4050;0;225];
            case 8
                coeffs=[5040;0;52920;0;66150;0;11025;0];
            case 9
                coeffs=[40320;0;564480;0;1058400;0;352800;0;11025];
            case 10
                coeffs=[362880;0;6531840;0;17146080;0;9525600;0;893025;0];
            case 11
                coeffs=[3628800;0;81648000;0;285768000;0;238140000;0;44651250;0;893025];
            otherwise
                error('Invalid n')
        end
        k=1+2*(n-1);
    else
        %The derivative of (1-x^2)^(-(k/2)) is
        %kx(1-x^2)^(-1-k/2)
        %That is k*x*(1-x^2)^(-(k1/2)) where k1=k+2.
        xVal=[1,0];

        %This polynomial is 1-x^2.
        polyX=[-1,0,1];

        %The first derivative
        n=n-1;
        coeffs=1;
        k=1;
        while(n>0)
            %Take the derivative of the denominator.
            coeffsNum=conv(k*xVal,coeffs);
            k=k+2;

            %Now, differentiate the numerator and multiply to get a
            %common denominator.
            coeffs=polyder(coeffs);
            coeffs=conv(coeffs,polyX);
            coeffs=polySum(coeffs,coeffsNum);
            n=n-1;
        end
    end
else
    k=1+2*(n-1);
end

if(~isempty(z))
    vals=zeros(size(z));
    %Evaluate the result.
    vals(:)=polyval(coeffs,z(:)).*sqrt(1-z(:).^2).^(-k);
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [x1,lun,err] = perform_douglas_rachford(A,y,options)

% perform_douglas_rachford -
%
%   [x1,lun,err] = perform_douglas_rachford(A,y,options);
%
%   Solve noiseless Basis Pursuit
%       x1 = argmin_x \sum |x[k]|   subj.to A*x=y
%
%   A can be a matrix or an operator A(x,dir,options)
%   with dir=1 for A*x and dir=-2 for A^{+}*x (pseudo inverse).
%
%   options.x can be an initial guess
%
%   options.niter is the number of iterations
%
%   Special thanks to Jalal Fadili for useful helps on this algorithm.
%
%   Copyright (c) 2008 Gabirel Peyre

niter = getoptions(options, 'niter', 100);
mu = getoptions(options, 'mu', 3);

if isnumeric(A)
    pA = getoptions(options, 'pseudo_inv', []);
    if isempty(pA)
        pA = pinv(A);
    end
end

% initial guess
x = getoptions(options, 'x', [] );
if isempty(x)
    if isnumeric(A)
        x = pA*y;
    else
        x = feval(A, y, -2, options);
    end
end

% target result (error monitoring
x0 = getoptions(options, 'x0', [] );

drawiter = getoptions(options, 'drawiter',0);
verb = getoptions(options, 'verb', 1);
isreal = getoptions(options, 'isreal', 0);

solve_tv = getoptions(options, 'solve_tv', 0);

lun = [];
err = [];
if isnumeric(A)
    x1 = x + pA*(y-A*x);
else
    x1 = y - feval(A,x, +1, options);
    x1 = x + feval( A, x1, -2, options);
end
if drawiter
    clf;
end

nrefresh = 100;
ndisp = max(2, ceil(niter/nrefresh) );

for i=1:niter
    if verb
        progressbar(i,niter);
    end
    % compute the intermediary point
    %   x <- x-x1 + thresh( 2*x1-x )

    u = 2*x1-x;

%    sparse_weight = getoptions(options, 'sparse_weight', 0,1);
%    u = u./sparse_weight;
    if not(solve_tv)
        % L1 regularization
        u = perform_thresholding( u, mu, 'soft' );
    else
        % TV regularization
        opts.method = 'gradient';
        opts.lambda = mu;
        opts.niter = getoptions(options, 'niter_tv', 50);
        opts.niter_inner = 1;
        opts.display = 0;
        opts.verb = 0;
        [u1,err,tv,lalist,Err] = perform_tv_denoising(u,opts);
        % u1 = u - mu*div( xi, options )
%        opts.xi = grad( (u1-u)/mu, options );
        u = u1;
    end
%    u = u.*sparse_weight;

    x = x-x1 + u;

    % perform projection step on the constraints
    if isnumeric(A)
        x1 = x + pA*(y-A*x);
    else
        x1 = y - feval(A,x,+1, options);
        x1 = x + feval( A, x1, -2, options);
    end

    if isreal
 %       x1 = real(x1); % force real signal
    end
    % check for error decay
    if nargout>=3 && not(isempty(x0))
        err(end+1) = norm(x0-x1, 'fro')^2;
    end
    if nargout>=2
        if not(solve_tv)
            lun(end+1) = sum(abs(x1(:))); % mean(sum(abs(x1)));
        else
            lun(end+1) = compute_total_variation(x1,options);
        end
    end
    if drawiter && mod(i,ndisp)==1  % && ntest==1
        if nb_dims(x1)==1
            n = size(x1,1);
            plot( real(x1(:,1:min(end,3))) );
            axis([1,n,-1,1]);
            axis tight;
            drawnow;
        else
            clf; imageplot(x1); drawnow;
        end
    end
end
% x1 = real(x1);

author: gpeyre

repo: matlab-toolboxes


text:
%  Figure 3.30b      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%

clf;
num=1;
a=5;

zeta =1;
den2=[1/(zeta*a) 1];
den1=[1 2*zeta 1];
den=conv(den1,den2);
t=0:.1:8;
y1=step(num,den,t);

zeta =.7;
den2=[1/(zeta*a) 1];
den1=[1 2*zeta 1];
den=conv(den1,den2);
y2=step(num,den,t);

zeta =.5;
den2=[1/(zeta*a) 1];
den1=[1 2*zeta 1];
den=conv(den1,den2);
y3=step(num,den,t);

axis([0 5 .1 .9])
plot(t,y1,'-',t,y2,'--',t,y3,'-.'),grid
title('Fig. 3.30b Step response with extra pole, \alpha= 5')
xlabel('\omega_n t')
ylabel('y(t)')


text:
function [K Kconv] = getK(G)
	K = G.KERNEL;
	if nargout > 1,
		Kconv = G.kernel;
	end

author: thomaskuestner


text:
function g = ggKernDiagGradient(kern, x, covDiag)

% GGKERNDIAGGRADIENT Compute gradient of the diagonal of GG kernel.
% FORMAT
% DESC computes the gradient of the diagonal of the kernel matrix for the
% gaussian gaussian kernel given a design matrix of inputs.
% RETURN g  : a vector containing the gradients
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% ARG covDiag : partial derivative wrt diagonal of the covariance matrix
%
% SEEALSO : ggKernParamInit, kernDiagCompute, ggKernDiagCompute
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

grad_sigma2Latent = (kern.sensitivity^2)*sum(covDiag);
grad_sensitivity = (2*kern.sigma2Latent*kern.sensitivity)*sum(covDiag);

g = [zeros(size(kern.precisionU(:)')) zeros(size(kern.precisionG(:)')) ...
    grad_sigma2Latent grad_sensitivity];


text:
function bgui(state)

    report_this_filefun(mfilename('fullpath'));

    if nargin == 0  % Build the GUI
        % Create the Figure Window
        figure_w_normalized_uicontrolunits('Units','normal','Position',[.25 .25 .5 .5], ...
            'Name','BGUI Example','NumberTitle','off', ...
            'Color','w');
        % Set the default units for uicontrols to 'normal'.
        % This makes it easier to place the uicontrols, and
        % the increase/decrease in size if the figure is
        % resized.
        set(gcf,'DefaultUicontrolUnits','normal')
        % Add a static text uicontrol
        st = uicontrol('Style','text','String','Enter Command:', ...
            'Position',[.1 .8 .8 .1], ...
            'BackgroundColor','w','ForegroundColor','k', ...
            'HorizontalAlignment','left');
        % Add the edit uicontrol
        h = uicontrol('Style','edit','Position',[.1 .6 .8 .2], ...
            'min',0,'max',2)
        % Add the push button
        i = uicontrol('Style','push','Position',[.1 .2 .8 .3], ...
            'String','Evaluate Edit Block', ...
            'Callback','bgui(''evaluate'')');
        % Store the handle to the edit block in the figures
        % UserData propert.  This is done to avoid the use of
        % global variable.  See the UserData section in the
        % Building a Graphical User Interface Guide for more
        % information.
        set(gcf,'UserData',h)

    elseif strcmp(state,'evaluate') % Push button selected
        % Get the handle to the edit uicontrol
        h = get(gcf,'UserData');
        % Extract the command
        com = get(h,'String');
        % Evaluate the string using eval(try,catch)
        eval(com,'error(''Invalid MATLAB Command'')')

    end

author: CelsoReyes


text:
% M1QN3  Solve a UNO using M1QN3
%
% THIS IS A LOW LEVEL FUNCTION - USE opti_m1qn3() INSTEAD!
%
%   [x,fval,exitflag,iter,feval] = m1qn3(fun,grad,x0,opts)
%
%   Input arguments:
%       fun - nonlinear function handle
%       grad - nonlinear gradient function handle (required)
%       x0 - initial solution guess
%       opts - solver options (see below)
%
%   Return arguments:
%       x - solution vector
%       fval - objective value at the solution
%       exitflag - exit status (see below)
%       iter - number of iterations taken by the solver
%       feval - number of function evaluations taken by the solver
%
%   Option Fields (all optional):
%       display - solver display level [0,1,2]
%       tolafun - gradient convergence tolerance
%       maxiter - maximum solver iterations {1000}
%       maxfeval - Maximum number of function evaluations {1500}
%       maxtime - Maximum solver execution time {1000s}
%       nupdate - Number of L-BFGS updates for Hessian (5-10) {5}
%       iterfun - Iteration callback function handle
%
%   Return Status:
%       0 - user exited
%       1 - successful termination
%       2 - input argument error
%       3 - line-search is blocked
%       4 - maximum iterations reached
%       5 - maximum function evaluations reached
%       6 - stop on dxmin during the line-search
%       7 - numerical error
%       8 - maximum time reached
%
%
%   Copyright (C) 2012 Jonathan Currie (IPL)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ldlup.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [L,d,p]=ldlup(L,d,j,g)
% updates LDL^T factorization when a unit j-th row and column
% are replaced by column g
% if the new matrix is definite (signalled by p=[]);
% otherwise, the original L,d and
% a direction p of null or negative curvature are returned
%
% d contains diag(D) and is assumed positive
% Note that g must have zeros in other unit rows!!!
%
function [L,d,p]=ldlup(L,d,j,g);

p=[];

test=0;
if test,
  disp('enter ldlup')
  A=L*diag(d)*L';A(:,j)=g;A(j,:)=g';
end;

n=size(d,1);
I=1:j-1;K=j+1:n;
if j==1,
  v=zeros(0,1);
  del=g(j);
  if del<=n*eps,
    p=[1;zeros(n-1,1)];
    if test,
      A,p
      Nenner=abs(p)'*abs(A)*abs(p);
      if Nenner==0, indef1=0 ,else indef1=(p'*A*p)/Nenner, end;
      disp('leave ldlup at 1')
    end;
    return;
  end;
  w=g(K)/del;
  L(j,I)=v';
  d(j)=del;
  if test,
    A1=L*diag(d)*L',A
    quot=norm(A1-A,1)/norm(A,1),
    disp('leave ldlup at 3')
  end;
  return;
end;

% now j>1, K nonempty
LII=L(I,I);
u=LII\g(I);
v=u./d(I);
del=g(j)-u'*v;
if del<=n*eps,
  p=[LII'\v;-1;zeros(n-j,1)];
  if test,
    A,p
    indef1=(p'*A*p)/(abs(p)'*abs(A)*abs(p))
    disp('leave ldlup at 2')
  end;
  return;
end;
LKI=L(K,I);
w=(g(K)-LKI*u)/del;
[LKK,d(K),q]=ldlrk1(L(K,K),d(K),-del,w);
if isempty(q),
  % work around expensive sparse L(K,K)=LKK
  L=[L(I,:);
     v', 1,L(j,K);
     LKI,w,LKK];
  d(j)=del;
  if test,
    A1=L*diag(d)*L',A
    quot=norm(A1-A,1)/norm(A,1),
    disp('leave ldlup at 4')
  end;
else
  % work around expensive sparse L(K,K)=LKK
  L=[L(1:j,:);
     LKI,L(K,j),LKK];
  pi=w'*q;
  p=[LII'\(pi*v-LKI'*q);-pi;q];
  if test,
    indef2=(p'*A*p)/(abs(p)'*abs(A)*abs(p)),
    disp('leave ldlup at 5')
  end;
end;

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
%  INTERNAL FUNCTION: Prediction step for Kalman filter
%
%  ::
%
%    [a,P,Tt,Rt,Record]=prediction_step(T,R,att,Ptt)
%    [a,P,Tt,Rt,Record]=prediction_step(T,R,att,Ptt,MUt,OMGt,DPHI,DT,Record,ExpandedFlag)
%
%  Args:
%
%     - **T** [matrix] : m x m state matrix (autoregressive part)
%     - **R** [matrix] : m x n state matrix (shock impacts)
%     - **att** [vector] : m x 1 state update
%     - **Ptt** [matrix] : m x m covariance matrix of state update
%     - **MUt** [matrix] : k x l matrix with forward information to match
%     - **OMGt** [matrix|{[]}] : kl x kl covariance matrix of forward
%       information
%     - **DPHI** [matrix] : Matrix representing the restrictions on future
%       shocks
%     - **DT** [matrix] : Convoluted impact of initial conditions
%     - **Record** [] : Holder of invariant information
%     - **ExpandedFlag** [true|false] : if true, returns the expanded state
%       vector including the future shocks. If false, returns only the
%       endogenous variables
%
%  Returns:
%     :
%
%     - **a** [vector] : m x 1 or mm x 1 vector of predictions
%     - **P** [matrix] : m x m or mm x mm covariance of predictions
%     - **Tt** [matrix] : m x m or mm x mm time-varying matrix, modifying the
%       impact of lagged endogenous (T)
%     - **Rt** [] : m x n or mm x nn matrix of impact of shocks
%     - **Record** [] : Holder of invariant information
%
%

author: jmaih


text:
function [dat, descrip, colors, h1, s1] = bars(obj, varnames, varargin)
% Bar plot for canlab_dataset object
%
% :Usage:
% ::
%
%    [dat, descrip, colors, h1, s1] = bars(obj, varnames, [optional inputs])
%
%    Takes inputs to barplot_columns
%    e.g., 'noviolin' to suppress violin plot
%    'colors' followed by colors cell for each bar
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2013 Tor Wager
%
%     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/>.
% ..
%
%
% :Inputs:
%
%   **obj:**
%        canlab_dataset object
%
%   **varnames:**
%        Cell string of variable names to plot
%
% :Optional Inputs:
%    **colors:**
%        defined colors (default are set by scn_standard_colors.m)
%
%   **nofig:**
%        do not generate figure
%
%   Takes any optional inputs to barplot_colored.m
%
%
% :Outputs:
%
%   **dat:**
%        data matrix for each variable
%
%   **descrip:**
%        the description for this variable
%
%   **colors:**
%        selected colors (default are set by scn_standard_colors.m)
%
%   **h1:**
%        figure handle
%
%   **s1:**
%        axis handle
%
% :Examples:
% ::
%
%     create_figure('NPS values - All subjects');
%
%     varnames = {'15' '13' '11' ' 9' '16' '14' '12' '10'};
%     xvals = [1 2 4 5 8 9 11 12];
%     colors = {[1 0 0] [0 1 0] [1 0 0] [0 1 0] [1 0 0] [0 1 0] [1 0 0] [0 1 0]};
%     bars(LevoNPS, varnames, 'x', xvals, 'colors', colors, 'XTickLabels', varnames, 'within', 'nofig');
%

n = length(varnames);
colors = scn_standard_colors(n);

% -------------------------------------------------------------------------
% DEFAULTS AND INPUTS
% -------------------------------------------------------------------------

% optional inputs with default values
% -----------------------------------
% - allowable_args is a cell array of argument names
% - avoid spaces, special characters, and names of existing functions
% - variables will be assigned based on these names
%   i.e., if you use an arg named 'cl', a variable called cl will be
%   created in the workspace

allowable_args = {'colors' 'nofig'};

default_values = {colors [0]};

% define actions for each input
% -----------------------------------
% - cell array with one cell for each allowable argument
% - these have special meanings in the code below
% - allowable actions for inputs in the code below are: 'assign_next_input' or 'flag_on'

actions = {'assign_next_input' 'flag_on'};

% logical vector and indices of which inputs are text
textargs = cellfun(@ischar, varargin);
whtextargs = find(textargs);

for i = 1:length(allowable_args)

    % assign default
    % -------------------------------------------------------------------------

    eval([allowable_args{i} ' =  default_values{i};']);

    wh = strcmp(allowable_args{i}, varargin(textargs));

    if any(wh)
        % Optional argument has been entered
        % -------------------------------------------------------------------------

        wh = whtextargs(wh);
        if length(wh) > 1, warning(['input ' allowable_args{i} ' is duplicated.']); end

        switch actions{i}
            case 'assign_next_input'
                eval([allowable_args{i} ' = varargin{wh(1) + 1};']);
                %varargin{wh(1) + 1} = [];

            case 'flag_on'
                eval([allowable_args{i} ' = 1;']);

            otherwise
                error(['Coding bug: Illegal action for argument ' allowable_args{i}])
        end

    end % argument is input
end

% END DEFAULTS AND INPUTS
% -------------------------------------------------------------------------

% REST OF BAR PLOT

[dat, ~, ~, descrip] = get_var(obj, varnames, varargin{:});

%descrip = cat(1, descrip{:});
%colors = cat(1, colors{:});

xvals = 1:n;

if nofig
    % Do nothing, use existing figure
else
    create_figure('bars');
end

[h1, s1] = barplot_columns(dat, 'nofig', varargin{:});
% set(h2, 'BarWidth', .9)
%colormap(colors)

mynames = strrep(varnames, '_', ' ');
set(gca, 'XTickLabel', mynames);

% STATS OUTPUT

fprintf('T-tests against zero\n');
fprintf('Description\tt(%3.0f)\tp\n', size(dat, 1) - 1);

% Print table - t-test against zero
for i = 1:n
    [h, p, ci, stats] = ttest(dat(:, i));
    if p < .001, str = '***';
    elseif p < .01, str = '**';
    elseif p < .05, str = '*';
    elseif p < .10, str = '+';
    else str = ' ';
    end

    fprintf('%s\t%3.2f\t%3.6f\t%s\n', descrip{i}, stats.tstat, p, str);
end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function util_playMov(h, w, options, Y1, Y2,sz)
% UTIL_PLAYMOV Play a video file which is stored by rows in a 2-D matrix.
% DESC util_playMov(h, w, options, Y1, Y2,sz)
% Utility to play video files (a sequence of frames) where each row corresponds
% to one frame, serialized by columns. Each frame must then be reshaped before plotted
% and for this the function must receice its actual dimensions as an argument.
% Sanity warning: This function has been augmented in many steps and is quite messy!!!
% ARG h: the true height of the frame
% ARG w: the true width of the frame
% ARG options: further options for the function
% ARG: Y1 the matrix which keeps the video by rows
% ARG: As Y2, in case the user wants to see two videos side by side
% (as long as both have the same dimensions)
% ARG: sz: The presented frame's dimensions (cf. the true dimensions) is
% divided by this number if the argument is present.
%
% COPYRIGHT: Andreas C. Damianou, 2011
%
% SEEALSO : preprocessVideo.m, saveMovieScript.m

% VARGPLVM

if ~nargin
    fprintf(1,'# Usage: \nplayMov(h,w,[p <p_diff>],Y1,<Y2>,<figSizeFactor>)\n');
    return
end

% The following is only used if options.saveImages exists and is 1.
if ~isfield(options,'imFormat')
    imFormat = 'jpg';
else
    imFormat = options.imFormat;
end

scrsz = get(0,'ScreenSize');
fntSz=18;
 try
    if ~isstruct(options)
        p = options; % pause after each frame
    else
        if ~isfield(options,'p')
            options.p = []; % pause after each frame ([] means that next frame is displayed after pressing any key)
        end
        p = options.p;
        % indices: a vector v in {0,1} with size equal to the number of frames,
        % for which v(i)==1 means it corresponds to a training
        % frame, otherwise to a test or generated frame.
        if isfield(options, 'indices')
            indices = options.indices;
        end

        % Save movie as an avi file, also given some more options
        if isfield(options, 'saveMovie') && options.saveMovie
            saveMovie=1;
            if ~isfield(options,'movieName') || ~isfield(options,'fps')
                error('You need to provide the movieName and the fps in the options structure\n');
            end
            aviobj = avifile(options.movieName, 'fps', options.fps); % e.g: 4
            if isfield(options,'quality')
                aviobj.Quality = options.quality;
            end
            if isfield(options,'compression') && ~isempty(options.compression) % e.g. 'Cinepak'
                aviobj.Compression = options.compression;
            end
        end

        % sz will cause each frame's size to be reduced by that number.
        if isfield(options, 'sz')
             scrsz(3) = scrsz(3)/options.sz;
             scrsz(4) = scrsz(4)/options.sz;
             fntSz = fntSz/options.sz;
        end

        % Like above (sz) but you can re-size the frame with different
        % scales for the horizontal and vertical axis.
         if isfield(options, 'sz1') && isfield(options, 'sz2')
             scrsz(3) = options.sz1;
             scrsz(4) = options.sz2;
        end
    end
    if exist('sz') && sz
        scrsz(3) = scrsz(3)/sz;
        scrsz(4) = scrsz(4)/sz;
    end

    figure('Position',[scrsz(3)/4.86 scrsz(4)/6.666 1.2*scrsz(3)/1.6457 0.6*scrsz(4)/1.4682])

    % if sz is zero pause to allow user to fix manually the size
    if exist('sz') && ~sz
        pause
    end

    % Oversampling
    if nargin > 4 && (length(p)==2)
        % if p is 2-D then the first element says how much to pause for the
        % second matrix and the second elem how many frames of the second to play in
        % between from the first matrix
        k=1;
        try
            for i=1:size(Y2,1)
                subplot(1,2,1);
                fr=reshape(Y2(i,:),h,w);
                imagesc(fr);
                colormap('gray');
                title(['frame' num2str(i)])
                text(scrsz(4)/6.666,0,['frame' num2str(i)],'HorizontalAlignment','center','VerticalAlignment','top','BackgroundColor',[.7 .9 .7], 'FontSize',fntSz);
                for j=1:p(2)
                    subplot(1,2,2);
                    fr=reshape(Y1(k,:),h,w);
                    k = k+1;
                    imagesc(fr);
                    colormap('gray');
                    drawnow;
                    pause(p(1))
                end
                if isfield(options, 'saveMovie') && options.saveMovie
                    frame = getframe ( gca );
                    aviobj = addframe ( aviobj, frame );
                end
             if isfield(options, 'saveImages')
                 frame = getframe ( gca );
                   imwrite (frame.cdata,[options.saveImages num2str(i) '.' imFormat]);
                end
            end
        catch e
            e.message
        end
    elseif nargin >4 && (length(p) < 2) % Two movies
        for i=1:size(Y1,1)
            subplot(1,2,1);
            fr=reshape(Y1(i,:),h,w);
            imagesc(fr);
            colormap('gray');
            subplot(1,2,2);
            fr=reshape(Y2(i,:),h,w);
            imagesc(fr);
            colormap('gray');
            title(['frame' num2str(i)])
           % text(scrsz(4)/6.666,0,['frame' num2str(i)],'HorizontalAlignment','center','VerticalAlignment','top','BackgroundColor',[.7 .9 .7], 'FontSize',fntSz);
            if isempty(p)
                pause
            else
                pause(p);
            end
            if isfield(options, 'saveMovie') && options.saveMovie
                frame = getframe ( gca );
                aviobj = addframe ( aviobj, frame );
            end
             if isfield(options, 'saveImages')
                 frame = getframe ( gca );
                   imwrite (frame.cdata,[options.saveImages num2str(i) '.' imFormat]);
                end
        end
    elseif nargin <=4 % One movie
        for i=1:size(Y1,1)
            fr=reshape(Y1(i,:),h,w);
            if isfield(options, 'addLogo')
              %  possible arguments for the colormap: 64, 256, etc.
               % h1=imshow(grs2rgb(uint8(fr), gray(64)),'Border','tight');
                 h1=imshow(grs2rgb(uint8(fr), colormap('gray')),'Border','tight'); % Default: gray(64)
            else
                 imagesc(fr);
                 colormap('gray');
            end
            drawnow;
            if exist('indices')
                if indices(i)
                    title(['Frame: ' num2str(i) ' (Training)'], 'Color','b','FontSize',fntSz)
                    text(scrsz(4)/6.666,0,['frame' num2str(i) ' (Training)'],'VerticalAlignment','top','BackgroundColor',[.7 .9 .7], 'FontSize',fntSz);
                else
                    title(['Frame: ' num2str(i) ' (Generated)'], 'Color','r','FontSize',fntSz)
                    text(scrsz(4)/6.666,0,['frame' num2str(i) ' (Generated)'],'VerticalAlignment','top','BackgroundColor',[1 0 0], 'FontSize',fntSz);
                end
            else
                title(['frame' num2str(i)])
               % text(scrsz(4)/6.666,0,['frame' num2str(i)],'HorizontalAlignment','center','VerticalAlignment','top','BackgroundColor',[.7 .9 .7], 'FontSize',fntSz);
            end
            if isfield(options,'addLogo')
               % The following only for png and export_fig
               % set(gcf, 'color', 'none'); % Set the figure background to transparent
               % set(gca, 'color', 'none'); % Set the axes background to transparent
                hold on; h2=imshow(options.addLogo,'Border','tight'); axis off;
                if isfield(options,'logoAlpha')
                    set(h2,'AlphaData',options.logoAlpha);
                end
                hold off;
            end

            if isempty(p)
                pause
            else
                pause(p);
            end
            if isfield(options, 'saveMovie') && options.saveMovie
                frame = getframe ( gca );
                aviobj = addframe ( aviobj, frame );
            end
            if isfield(options, 'saveImages')
                frame = getframe ( gca );
                imageNo = num2str(i);
                % Pad the image number with leading zeros so that all
                % images have the same length in their filename (useful for
                % sorting frames).
                lengthAll = length(num2str(size(Y1,1)));
                padding = lengthAll - length(imageNo);
                imageNo = [repmat(['0'],1,padding) imageNo];
             %   if isfield(options, 'addLogo')
                    if ~isfield(options, 'saveImagesQuality')
                        options.saveImagesQuality = 75;
                    end
                    if ~isfield(options, 'saveImagesMagn')
                        options.saveImagesMagn =1;
                    end
                    if strcmp(imFormat,'manual')
                        % pause and wait for the user to save it manually
                        fprintf('User saves: %s \n', [options.saveImages imageNo '.' imFormat]);
                        pause
                    else
                         % Use export_fig which crops the borders from each
                        % frame.
                        % if ~strcmp(imFormat, 'eps')
                            export_fig([options.saveImages imageNo], ['-' imFormat], ['-q' num2str(options.saveImagesQuality)], ...
                                ['-m' num2str(options.saveImagesMagn)]);
                        %else
                        %    print('-depsc', [options.saveImages imageNo '.' imFormat]);
                        %    system(['epstopdf ' options.saveImages imageNo '.' imFormat]);
                        %end
                    end
              %  else
               %    imwrite (frame.cdata,[options.saveImages imageNo '.' imFormat]);
                %end
            end
        end
    else
        fprintf(1,'playMov(h,w,p,Y1,Y2)\n')
    end

 catch e
     e.getReport
     if isfield(options, 'saveMovie') && options.saveMovie
         aviobj = close (aviobj);
     end
 end

if isfield(options, 'saveMovie') && options.saveMovie
    aviobj = close ( aviobj );
end

author: SheffieldML

repo: vargplvm


text:
% Function to compute the gradient of the given optical flow
%
%    Copyright (C) 2013  Anestis Papazoglou
%
%    You can redistribute and/or modify this software for non-commercial use
%    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/>.
%
%    For commercial use, contact the author for licensing options.
%
%    Contact: a.papazoglou@sms.ed.ac.uk

function result = getFlowGradient( flow )

    if( iscell( flow ) )
        framesNumber = length( flow );

        gradients = cell( 1, framesNumber );
        for( i = 1: framesNumber )
            grad( :, :, 1 ) = gradient( single( flow{ i }( :, :, 1 ) ) );
            [ ~, grad( :, :, 2 ) ] = ...
                gradient( single( flow{ i }( :, :, 2 ) ) );

            gradients{ i } = grad;
        end
        result = gradients;
    else
        [ height, width, ~ ] = size( flow );
        grad = zeros( height, width, 2, 'single' );
        if( ~isfloat( flow ) )
            flow = single( flow );
        end

        grad( :, :, 1 ) = gradient( flow( :, :, 1 ) );
        [ ~, grad( :, :, 2 ) ] = gradient( flow( :, :, 2 ) );

        result = grad;
    end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%% Copyright (C) 2014, 2016, 2018-2019, 2022 Colin B. Macdonald
%% Copyright (C) 2022 Alex Vong
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defop  Method   @@sym rdivide {(@var{x}, @var{y})}
%% @defopx Operator @@sym {@var{x} ./ @var{y}} {}
%% Element-wise forward slash division of symbolic expressions.
%%
%% Example:
%% @example
%% @group
%% syms x
%% A = sym([1 137; 3 4])
%%   @result{} A = (sym 2×2 matrix)
%%       ⎡1  137⎤
%%       ⎢      ⎥
%%       ⎣3   4 ⎦
%% B = [x pi; 2*x 8]
%%   @result{} B = (sym 2×2 matrix)
%%       ⎡ x   π⎤
%%       ⎢      ⎥
%%       ⎣2⋅x  8⎦
%% A ./ B
%%   @result{} ans = (sym 2×2 matrix)
%%       ⎡ 1   137⎤
%%       ⎢ ─   ───⎥
%%       ⎢ x    π ⎥
%%       ⎢        ⎥
%%       ⎢ 3      ⎥
%%       ⎢───  1/2⎥
%%       ⎣2⋅x     ⎦
%% @end group
%% @end example
%%
%% Either @var{x} or @var{y} can be scalar:
%% @example
%% @group
%% A ./ 2
%%   @result{} ans = (sym 2×2 matrix)
%%       ⎡1/2  137/2⎤
%%       ⎢          ⎥
%%       ⎣3/2    2  ⎦
%% 2 ./ B
%%   @result{} ans = (sym 2×2 matrix)
%%       ⎡2   2 ⎤
%%       ⎢─   ─ ⎥
%%       ⎢x   π ⎥
%%       ⎢      ⎥
%%       ⎢1     ⎥
%%       ⎢─  1/4⎥
%%       ⎣x     ⎦
%% @end group
%% @end example
%%
%% Finally, the can both be scalar:
%% @example
%% @group
%% 2 ./ x
%%   @result{} ans = (sym)
%%       2
%%       ─
%%       x
%% @end group
%% @end example
%% @seealso{@@sym/ldivide, @@sym/mrdivide}
%% @end defop

function z = rdivide(x, y)

  %% 2022-06: TODO cannot simply call hadamard_product for sympy <1.9,
  %% see upstream: https://github.com/sympy/sympy/issues/8557

  cmd = { '(x,y) = _ins'
          'if x.is_Matrix and y.is_Matrix:'
          '    y_eltwise_recip = y.applyfunc(lambda a: 1/a)'
          '    if Version(spver) < Version("1.9"):'
          '        try:'
          '            return x.multiply_elementwise(y_eltwise_recip)'
          '        except:'
          '            pass'
          '    return hadamard_product(x, y_eltwise_recip)'
          'if not x.is_Matrix and y.is_Matrix:'
          '    return y.applyfunc(lambda a: x/a),'
          'else:'
          '    return x/y,' };

  z = pycall_sympy__ (cmd, sym(x), sym(y));

end

%!test
%! % scalar
%! syms x
%! assert (isa (x ./ 1, 'sym'))
%! assert (isa (x ./ x, 'sym'))
%! assert (isequal (x ./ 1, x))
%! assert (isequal (x ./ x, sym(1)))

%!test
%! % matrix-scalar
%! D = 2*[0 1; 2 3];
%! A = sym(D);
%! assert (isequal ( A./2 , D/2  ))
%! assert (isequal ( A./sym(2) , D/2  ))
%! assert (isequal ( D./sym(2) , D/2  ))

%!test
%! % matrix ./ matrix
%! D = [1 2; 3 4];
%! A = sym(D);
%! assert (isequal ( A./A , D./D  ))
%! assert (isequal ( A./D , D./D  ))
%! assert (isequal ( D./A , D./D  ))

%!test
%! % matrix ./ matrix with symbols
%! syms x y
%! A = [x y; x^2 2*y];
%! B = [y x; x y];
%! assert (isequal ( A./A , sym(ones(2,2)) ))
%! assert (isequal ( A./B , [x/y y/x; x 2] ))

%!test
%! % scalar ./ matrix
%! D = [1 2; 3 4];
%! A = sym(D);
%! assert (isequal ( 12./A , 12./D  ))

repo: octsympy


text:
function datwrite(filename, x, t)
%DATWRITE Write data to ascii file.
%
%	Description
%
%	DATWRITE(FILENAME, X, T) takes a matrix X of input vectors and a
%	matrix T of target vectors and writes them to an ascii file named
%	FILENAME. The file format is as follows: the first row contains the
%	string NIN followed by the number of inputs, the second row contains
%	the string NOUT followed by the number of outputs, and the third row
%	contains the string NDATA followed by the number of data vectors.
%	Subsequent lines each contain one input vector followed by one output
%	vector, with individual values separated by spaces.
%
%	See also
%	DATREAD
%

%	Copyright (c) Ian T Nabney (1996-2001)

nin = size(x, 2);
nout = size(t, 2);
ndata = size(x, 1);

fid = fopen(filename, 'wt');
if fid == -1
  error('Failed to open file.')
end

if size(t, 1) ~= ndata
  error('x and t must have same number of rows.');
end

fprintf(fid, ' nin   %d\n nout  %d\n ndata %d\n', nin , nout, ndata);
for i = 1 : ndata
  fprintf(fid, '%13e ', x(i,:), t(i,:));
  fprintf(fid, '\n');
end

flag = fclose(fid);
if flag == -1
  error('Failed to close file.')
end


text:
%% SQUARESTOKE Stokes equations on the unit square
%
%   SQUARESTOKE computes P2-P1 approximations of the Stokes equations in
%   the unit square on a sequence of meshes obtained by uniform refinement.
%   It plots the approximation error (pressue in L2 norm and velocity in H1
%   norm) vs the number of dof.
%
% See also StokesP2P1, collidingflow, squarePoisson
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

close all; clear all;
%% Set up
maxIt = 4;
N = zeros(maxIt,1); erru = zeros(maxIt,1); errp = zeros(maxIt,1);

%% Generate initial mesh
[node,elem] = squaremesh([0 1 0 1], 0.25);
bdFlag = setboundary(node,elem,'Dirichlet');

%% PDE and options
pde = Stokesdata2;

%% Finite Element Method
for k = 1:maxIt
    % solve the equation
    [u,p,edge,A] = StokesP2P1(node,elem,pde,bdFlag);
    N(k) = length(u)+length(p);
    if N(k) < 2e3 % show solution for small size
        figure(1);  showresult(node,elem,p);
    end
    % compute error
    uI = pde.exactu([node; (node(edge(:,1),:)+node(edge(:,2),:))/2]);
    erru(k) = sqrt((u-uI(:))'*A*(u-uI(:)));
    errp(k) = getL2error(node,elem,pde.exactp,p);
    % refine mesh
    [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
end

%% Plot convergence rates
figure(2); clf;
showrate2(N,erru,2,'-*','||Du_I-Du_h||',N,errp,2,'k-+','|| p - p_h||');

repo: Aero-matlab


text:
%
%  MODEL REDUCTION BY THE ALGORITHMS LRSRM AND DSPMR. THE GOAL IS TO
%  GENERATE A "NUMERICALLY MINIMAL REALIZATION" OF THE GIVEN SYSTEM
%  AS WELL AS A REDUCED SYSTEM OF RELATIVELY SMALL ORDER.
%
%  This demo program shows how the model reduction routines 'lp_lrsrm'
%  and 'lp_dspmr' work. Also, the use of 'lp_lradi', supplementary
%  routines, and user-supplied functions is demonstrated.

% -----------------------------------------------------------------------
% Generate test problem
% -----------------------------------------------------------------------
%
% As test example, we use an FEM-semidiscretized problem, which leads to
% a generalized system where M (the mass matrix) and N (the negative
% stiffness matrix) are sparse, symmetric, and definite.

load rail821   % load the matrices M N Btilde Ctilde of the generalized
               % system

%load rail3113   % Uncomment this to get an example of larger order.

disp('Problem dimensions:')

n = size(M,1)   % problem order (number of states)
m = size(Btilde,2)   % number of inputs
q = size(Ctilde,1)   % number of outputs

% -----------------------------------------------------------------------
% Initialization/generation of data structures used in user-supplied
% functions and computation of ADI shift parameters
% -----------------------------------------------------------------------
%
% See 'demo_u1', 'demo_u2', 'demo_u3', and 'demo_l1' for more detailed
% comments.

name = 'msns';
[M0,MU,N0,B0,C0,prm,iprm] = msns_pre(M,N,Btilde,Ctilde);  % preprocessing
msns_m_i(M0,MU,N0);   % initialization for multiplication with A0
msns_l_i;   % initialization for solving systems with A0

disp('Parameters for heuristic algorithm which computes ADI parameters:')
l0 = 20   % desired number of distinct shift parameters
kp = 50   % number of steps of Arnoldi process w.r.t. A0
km = 25   % number of steps of Arnoldi process w.r.t. inv(A0)

b0 = ones(n,1);   % This is just one way to choose the Arnoldi start
                  % vector.

p = lp_para(name,[],[],l0,kp,km,b0);   % computation of ADI shift
                                       % parameters

disp('Actual number of ADI shift parameters:');
l = length(p)

disp('ADI shift parameters:');
p

msns_s_i(p)   % initialization for shifted systems of linear equations
              % with A0+p(i)*I  (i = 1,...,l)

% -----------------------------------------------------------------------
% Solution of Lyapunov equations A0*X0+X0*A0' = -B0*B0' and
% A0'*X0+X0*A0 = -C0'*C0
% -----------------------------------------------------------------------

disp('Parameters for stopping criteria in LRCF-ADI iteration:')
max_it = 200   % (large value)
min_res = 0   % (avoided)
with_rs = 'S'   % ("activated")
min_in = 0   % (avoided)

zk = 'Z';
rc = 'C';
Bf = [];
Kf = [];
info = 3;

disp('... solving A0*XB0+XB0*A0'' = - B0*B0''...');
tp = 'B';

figure(1), hold off; clf;   % (lp_lradi will plot residual history.)

[ZB0,flag_B] = lp_lradi(tp,zk,rc,name,Bf,Kf,B0,p,max_it,min_res,...
               with_rs,min_in,info);
                                       % compute ZB0

title('LRCF-ADI for CALE  A_0X_{B0}+X_{B0}A_0^T = -B_0B_0^T')
disp('Termination flag:')
flag_B
disp('Size of ZB0:');
size_ZB0 = size(ZB0)

disp('... solving A0''*XC0+XC0*A0 = - C0''*C0...');
tp = 'C';

figure(2), hold off; clf;   % (lp_lradi will plot residual history.)

[ZC0,flag_C] = lp_lradi(tp,zk,rc,name,Bf,Kf,C0,p,max_it,min_res,...
               with_rs,min_in,info);
                                       % compute ZC0

title('LRCF-ADI for CALE  A_0^T X_{C0} + X_{C0} A_0 = -C_0^TC_0')
disp('Termination flag:')
flag_C
disp('Size of ZC0:');
size_ZC0 = size(ZC0)

% -----------------------------------------------------------------------
% Plot the transfer function of the system for a certain frequency range
% -----------------------------------------------------------------------

disp('... computing transfer function of original system ...');

freq = lp_lgfrq(1e-10,1e10,200);   % generate a set of 200 "frequency
                                   % sampling points" in the interval
                                   % [10^-10,10^+10].
G = lp_trfia(freq,N,Btilde,Ctilde,[],M);   % compute "transfer function
                                           % sample" for these frequency
                                           % points
nrm_G = lp_gnorm(G,m,q);   % compute norms of the "transfer function
                           % sample" for these frequency points

figure(3); hold off; clf;
loglog(freq,nrm_G,'k:');
xlabel('\omega');
ylabel('Magnitude');
t_text = 'dotted: ||G||';
title(t_text);
pause(1)

% -----------------------------------------------------------------------
% Generate reduced systems of high accuracy and possibly high order
% -----------------------------------------------------------------------

disp(' ')
disp('Generate reduced systems of high accuracy and possibly high order')
disp('-----------------------------------------------------------------')

disp('Parameters for model reduction:')
max_ord = []   % (avoided)
tol = 1e-14   % (This criterion determines the reduced order. The very
              % small value is chosen to generate a "numerically minimal
              % realization".)

disp('... computing reduced system by LRSRM ...');
[Ars,Brs,Crs] = lp_lrsrm(name,B0,C0,ZB0,ZC0,max_ord,tol);   % run LRSRM

disp('Reduced order:')
disp(length(Ars))

Grs = lp_trfia(freq,Ars,Brs,Crs,[],[]);   % compute "transfer function
                                          % sample" for reduced system
nrm_dGrs = lp_gnorm(G-Grs,m,q);   % compute norm of DIFFERENCE of
                                  % transfer function samples of original
                                  % and reduced system.
figure(3); hold on
loglog(freq,nrm_dGrs,'r-');
t_text = [t_text, ',    solid: ||G-G_{LRSRM}||'];
title(t_text); pause(1)

disp('... computing reduced system by DSPMR ...');
[Ard,Brd,Crd] = lp_dspmr(name,B0,C0,ZB0,ZC0,max_ord,tol);   % run DSPMR

disp('Reduced order:')
disp(length(Ard))

Grd = lp_trfia(freq,Ard,Brd,Crd,[],[]);   % compute "transfer function
                                          % sample" for reduced system
nrm_dGrd = lp_gnorm(G-Grd,m,q);   % compute norm of DIFFERENCE of
                                  % transfer function samples of original
                                  % and reduced system.
figure(3); hold on
loglog(freq,nrm_dGrd,'b--'); pause(1)
t_text = [t_text, ',    dashed: ||G-G_{DSPMR}||'];
title(t_text); pause(1)

% -----------------------------------------------------------------------
% Generate reduced systems of low order
% -----------------------------------------------------------------------

disp(' ')
disp('Generate reduced systems of low order')
disp('-------------------------------------')

disp('Parameters for model reduction:')
max_ord = 25   % (This criterion determines the reduced order.)
tol = 0    % (avoided)

disp('... computing reduced system by LRSRM ...');
[Ars,Brs,Crs] = lp_lrsrm(name,B0,C0,ZB0,ZC0,max_ord,tol);   % run LRSRM

disp('Reduced order:')
disp(length(Ars))

Grs = lp_trfia(freq,Ars,Brs,Crs,[],[]);   % compute "transfer function
                                          % sample" for reduced system
nrm_dGrs = lp_gnorm(G-Grs,m,q);   % compute norm of DIFFERENCE of
                                  % transfer function samples of original
                                  % and reduced system.
figure(3); hold on
loglog(freq,nrm_dGrs,'r-');

disp('... computing reduced system by DSPMR ...');
[Ard,Brd,Crd] = lp_dspmr(name,B0,C0,ZB0,ZC0,max_ord,tol);   % run DSPMR

disp('Reduced order:')
disp(length(Ard))

Grd = lp_trfia(freq,Ard,Brd,Crd,[],[]);   % compute "transfer function
                                          % sample" for reduced system
nrm_dGrd = lp_gnorm(G-Grd,m,q);   % compute norm of DIFFERENCE of
                                  % transfer function samples of original
                                  % and reduced system.
figure(3); hold on
loglog(freq,nrm_dGrd,'b--');

% -----------------------------------------------------------------------
% Destroy global data structures
% -----------------------------------------------------------------------

msns_m_d;
msns_l_d;
msns_s_d(p);


text:
function prob_test049 ( )

%*****************************************************************************80
%
%% TEST049 tests DIRICHLET_MIX_MEAN, DIRICHLET_MIX_SAMPLE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  comp_num = 2;
  elem_num = 3;
  sample_num = 1000;

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST049\n' );
  fprintf ( 1, '  For the Dirichlet Mixture PDF:\n' );
  fprintf ( 1, '  DIRICHLET_MIX_SAMPLE samples;\n' );
  fprintf ( 1, '  DIRICHLET_MIX_MEAN computes the mean;\n' );

  a(1:3,1:2) = [ 0.250, 1.500; ...
                 0.500, 0.500; ...
                 1.250, 2.000 ];

  comp_weight(1:2) = [ 1.0, 2.0 ];

  check = dirichlet_mix_check ( comp_num, elem_num, a, comp_weight );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST049 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of elements ELEM_NUM =   %6d\n', elem_num );
  fprintf ( 1, '  Number of components COMP_NUM = %6d\n', comp_num );
  fprintf ( 1, '  PDF parameters A(ELEM,COMP):\n' );
  fprintf ( 1, '\n' );
  for elem_i = 1 : elem_num
    fprintf ( 1, '  %14f  %14f\n', a(elem_i,1:comp_num) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Component weights:\n' );
  fprintf ( 1, '\n' );
  for comp_i = 1 : comp_num
    fprintf ( 1, '  %6d  %14f\n', comp_i, comp_weight(comp_i) );
  end

  mean(1:elem_num) = dirichlet_mix_mean ( comp_num, elem_num, a, comp_weight );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF mean:\n' );
  fprintf ( 1, '\n' );

  for elem_i = 1 : elem_num
    fprintf ( 1, '  %6d  %14f\n', elem_i, mean(elem_i) );
  end

  for sample_i = 1 : sample_num
    [ x(1:elem_num,sample_i), comp, seed ] = dirichlet_mix_sample ( ...
      comp_num, elem_num, a, comp_weight, seed );
  end

  xmax(1:elem_num) = r8row_max ( elem_num, sample_num, x );

  xmin(1:elem_num) = r8row_min ( elem_num, sample_num, x );

  mean(1:elem_num) = r8row_mean ( elem_num, sample_num, x );

  variance(1:elem_num) = r8row_variance ( elem_num, sample_num, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size = %6d\n', sample_num );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Observed Min, Max, Mean, Variance:\n' );
  fprintf ( 1, '\n' );

  for elem_i = 1 : elem_num
    fprintf ( 1, '  %6d  %14f  %14f  %14f  %14f\n', elem_i, xmin(elem_i), ...
      xmax(elem_i), mean(elem_i), variance(elem_i) );
  end

  return
end


text:
function out = apply_logisticRegression(C, X, varargin)
% APPLY_LOGISTICREGRESSION - Apply existing logistic regresion classifier
%
%Synopsis:
%   C = apply_logisticRegression(C, X)
%
%Arguments:
%   X: DOUBLE [NxM] - Data matrix, with N features and M samples.
%   C: STRUCT           - a logistic regression classifier structure.  Must include the
%                           field 'b'.
%   OPT: PROPLIST       - Structure or property/value list of optional
%                           properties. Options are also passed to clsutil_shrinkage.
%     'OriginalOutput'  - BOOL (default 0): If true, the orginal logistic regression output is returned (range [0 1])
%Returns:
%   out: FLOAT[]        - an array containing the classifier score for each
%                           sample, in range [-1 1] (or alternatively [0 1])
%Description:
%   APPLY_LOGISTICREGRESSION applies a logistic regresion classifier given data and a
%   trained LR classifier.
%
%
%Examples:
%   apply_logisticRegression(C, X))
%
%See also:
%   TRAIN_LOGISTICREGRESSION

props= {'OriginalOutput'      0                             'BOOL'
       };

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

misc_checkType(X, 'DOUBLE');
misc_checkType(C, 'STRUCT(b)');

pihat = mnrval(C.b, X');

pihat = pihat';
out = pihat(1,:);
if ~opt.OriginalOutput
  out = 2*(out - 0.5);
end

author: bbci


text:
%Example demonstrating peak interpolation (finds the peak of a function,
%interpolated between samples using the peak and its two neighbors).
%This example finds the peak of sin(2*pi*f*x+phi) for unevenly spaced x

N_p=10;%number of points
f=1/N_p;%frequency (relative to sampling frequency)

x=2*pi*rand(1,N_p);
x=sort(x);
y=sin(2*pi*f*x);

x_max_th=N_p/4;%theoretical x for maximum y
y_max_th=1;%theoretical maximum y

[y_max_int idx_max]=max(y);
x_max_int=x(idx_max);%integer estimate of peak position

%do a parabolic estimate
[x_max_p y_max_p A_p]=crit_interp_p(y(idx_max+(-1:1)),x(idx_max+(-1:1)));

%do a Gaussian estimate
[x_max_g y_max_g A_g]=crit_interp_g(y(idx_max+(-1:1)),x(idx_max+(-1:1)));

%plot the curves
N_plot=1000;%number of points to plot
x_plot=linspace(x(idx_max-2),x(idx_max+2),N_plot);
y_th=sin(2*pi*f*x_plot);
y_p=A_p(1)*x_plot.^2+A_p(2)*x_plot+A_p(3);
y_g=A_g(1)*exp(-A_g(2)*(x_plot-A_g(3)).^2);

fprintf('True: y_max=%f at x_max=%f\nParabolic Estimate: y_max=%f at x_max=%f\nGaussian Estimate: y_max=%f at x_max=%f\n',...
    y_max_th,x_max_th,y_max_p,x_max_p,y_max_g,x_max_g)

figure(1)
clf
subplot(2,1,1)
plot(x,y,'o',x_plot,y_th,'g',x_plot,y_p,'k',x_plot,y_g,'r')
xlabel('x')
ylabel('y')
legend('Sampled points','Theoretical Curve','Parabolic interpolation','Gaussian Interpolation','location','Best')

subplot(2,1,2)
plot(x(idx_max+(-1:1)),y(idx_max+(-1:1)),'o',x_plot,y_th,'g',x_plot,y_p,'k',x_plot,y_g,'r',x_max_p, y_max_p,'kx',x_max_g, y_max_g,'rx',x_max_th, y_max_th,'gx')
xlabel('x')
ylabel('y')
legend('Sampled points','Theoretical Curve','Parabolic interpolation','Gaussian Interpolation','location','Best')

author: Sable


text:
function [hwm,hwmm]=webmap_of_catalog(catalog,force)
    % webmap_of_catalog plot earthquakes in a browser window
    % using the geoweb toolbox
    %
    % see webmap, webmarker
    iconFolder=fullfile('resrc','img','ico');
    icons=strcat('ball_',{'cyan';'blue1';'blue2';'purp1';'puce';'pink';'red';'orange';'yellow';'lime';'green';'mint'},'_30.png');
    %icons=icons(1:2:end);
    % TODO: future idea: this can  execute scripts in links. eg: col
    % TODO: future idea: allow this to be colored by depth, magnitude, or date
    disp(['Plotting events on web map.  They will plot from smallest to largest. ', newline,...
        'Within each magnitude bin, they are plotted according to depth range']);
    MAXQUAKES=1000;
    if catalog.Count >= 1000 && (~exist('force','var') || ~force)
        errordlg('Too many events for webmap.  Reduce events to < 1000, and try again');
        return
    end
    %catalog.sort('Depth','descend')
    % dep2color= @(x) [(x - min(x)) /max(x- min(x)), repmat(0,numel(x),1), sqrt(1-(x - min(x)) /max(x- min(x)))];
    smallestEvent=min(catalog.Magnitude);
    biggestEvent=max(catalog.Magnitude);
    eventRange=biggestEvent-smallestEvent;

    scaleit = @(n) (((n-smallestEvent) ./ eventRange ).*1.4 + 0.6).^1.8 ./ 4;

    %scaleit = @(n)(( n + abs(min(n)) + 1 )/4) .^ (1/3);
    %TODO include the MagnitudeType
    %disp('desc')
    desc={};
    for i=1:catalog.Count
        desc(i)={sprintf('(%8.4f, %8.4f) depth: %4.2f km, mag %3.2f\n',...
        catalog.Latitude(i), catalog.Longitude(i), catalog.Depth(i), catalog.Magnitude(i))};
    end
    %disp('title');
    tit=cellstr(char(catalog.Date,'uuuu-MM-dd HH:mm:ss'));
    hwm=webmap('World Topographic Map');
    wmlimits(hwm,bounds2(catalog.Latitude), bounds2(catalog.Longitude))
    %disp('wmmarker');

    magrange = [10:-1:-2];
    cmp=jet(numel(magrange));

    for j=2:numel(magrange)
        name(j)={sprintf('%.1f %c M < %.1f',magrange(j),char(8804),magrange(j-1))};
        idx(j)={catalog.Magnitude >= magrange(j) & catalog.Magnitude < magrange(j-1)};
        nInRange(j)=sum(idx{j});
    end

    % get rid of empty categories
    magrange(nInRange==0)=[];
    name(nInRange==0)=[];
    idx(nInRange==0)=[];
    nInRange(nInRange==0)=[];
    tot=cumsum(nInRange);
    toomany = find(tot > MAXQUAKES,1,'first');
    nDepthBins=numel(icons);
    [bins,edgs]=discretize(catalog.Depth,nDepthBins);
    cmp=jet(nDepthBins);
    %quakeIcons={}
    for q=1:numel(edgs)-1
        depthnames(q) ={ sprintf(' %.1f %c Z < %.1f km', edgs(q), char(8804), edgs(q+1))};
    end

    %% using multi-color option for wmmarker is super slow. Either
    % add specialized icons or do a double loop.
    %for j=numel(magrange) : -1 :1
        %if j > toomany+1
        %    continue;
        %end

        %for q=nDepthBins:-1:1
        %    quakeIcons(bins==q)={fullfile( iconFolder , icons{mod(q,numel(icons))+1})};
        %end

        for q=nDepthBins:-1:1
            %thisidx = idx{j};
            thisidx = bins==q;
            %sum(thisidx)
            if ~any(thisidx), continue, end
            % thiscolor = cmp(q,:);
            %thisIcon=fullfile( iconFolder , icons{mod(q,numel(icons))+1});
            if ~exist('hwmm','var')
                hwmm=wmmarker(catalog.Latitude(thisidx),catalog.Longitude(thisidx),...
                    'OverlayName',[depthnames{q}],...
                    ... 'OverlayName',[name{j}],...
                    ...'Color',thiscolor,...
                    ...'Icon',quakeIcons(thisidx),...
                    'Icon',fullfile( iconFolder , icons{mod(q,numel(icons))+1}),...
                    'Description',desc(thisidx),...
                    'FeatureName',tit(thisidx),...
                    'Alpha',0.7,...
                    'IconScale',scaleit(catalog.Magnitude(thisidx))./1.5,...
                    'Autofit',false);
            else
                hwmm(end+1)=wmmarker(catalog.Latitude(thisidx),catalog.Longitude(thisidx),...
                    ...'OverlayName',[name{j}],...
                    'OverlayName',[depthnames{q}],...
                    ...'Color',thiscolor,...
                    ...'Icon',quakeIcons(thisidx),...
                    'Icon',fullfile( iconFolder , icons{mod(q,numel(icons))+1}),...
                    ...'Icon',quakeIcons{q},...
                    'Description',desc(thisidx),...
                    'FeatureName',tit(thisidx),...
                    'Alpha',0.7,...
                    'IconScale',scaleit(catalog.Magnitude(thisidx))./1.5,...
                    'Autofit',false);
            end
        end
    %end
    %if ~isempty(toomany)
    %    fprintf('Only the largest %d events (of %d) were shown\n',tot(toomany),sum(nInRange));
    %end
    %{
    hwmm=wmmarker(catalog.Latitude, catalog.Longitude, 'FeatureName',tit);...,...
        ...'IconScale',scaleit(catalog.Magnitude),...
        ...'Description',desc,'Color',dep2color(catalog.Depth), 'FeatureName',tit);
    %}
end