text:
function [ywarp,dwarp_dt,dwarp_dtheta,d2warp_dthetadt] = outwarp_negpowc1(hyp,y,invflag)
%GPLITE_NOISEFUN Noise function for lite Gaussian Process regression.
%   SN2 = GPLITE_NOISEFUN(HYP,X,NOISEFUN) computes the GP noise function
%   NOISEFUN, that is the variance of observation noise evaluated at test
%   points X. HYP is a single column vector of noise function
%   hyperparameters. NOISEFUN is a numeric array whose elements specify
%   features of the noise function, as follows:
%
%   See also GPLITE_COVFUN, GPLITE_MEANFUN.

if nargin < 2; y = []; end
if nargin < 3 || isempty(invflag); invflag = false; else; invflag = true; end

if invflag && nargout > 1
    error('outwarp_fun:InverseOnly', ...
        ['When calling for the inverse output warping function, only one function output is expected.']);
end

%--------------------------------------------------------------------------
% CUSTOM: Number of hyperparameters
Noutwarp = 2;       % # hyperparameters of the output warping function
%--------------------------------------------------------------------------

N = size(y,1);      % Number of training points

% Return number of output warping function hyperparameters and additional info
if ischar(hyp)
    ywarp = Noutwarp;
    if nargout > 1

        if isempty(y); y = [0;1]; end

        % Initialize bounds for all hyperparameters
        outwarp_info.LB = -Inf(1,Noutwarp);
        outwarp_info.UB = Inf(1,Noutwarp);
        outwarp_info.PLB = -Inf(1,Noutwarp);
        outwarp_info.PUB = Inf(1,Noutwarp);
        outwarp_info.x0 = NaN(1,Noutwarp);

        %------------------------------------------------------------------
        % CUSTOM: Initialize hyperparameter bounds and other details

        % Threshold parameter
        outwarp_info.LB(1) = min(y);
        outwarp_info.UB(1) = max(y);
        outwarp_info.PLB(1) = min(y);
        outwarp_info.PUB(1) = max(y);
        outwarp_info.x0(1) = NaN;

        % Power exponent k (log space)
        outwarp_info.LB(2) = -Inf;
        outwarp_info.UB(2) = Inf;
        outwarp_info.PLB(2) = -3;
        outwarp_info.PUB(2) = 3;
        outwarp_info.x0(2) = 0;

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

        % Assign handle of current output warping function
        outwarp_info.outwarpfun = str2func(mfilename);

        % Plausible starting point
        idx_nan = isnan(outwarp_info.x0);
        outwarp_info.x0(idx_nan) = 0.5*(outwarp_info.PLB(idx_nan) + outwarp_info.PUB(idx_nan));

        dwarp_dt = outwarp_info;

    end

    return;
end

[Nhyp,Ns] = size(hyp);      % Hyperparameters and samples

if Nhyp ~= Noutwarp
    error('outwarp_fun:WrongLikHyp', ...
        ['Expected ' num2str(Noutwarp) ' output warping function hyperparameters, ' num2str(Nhyp) ' passed instead.']);
end
if Ns > 1
    error('outwarp_fun:nosampling', ...
        'Output warping function output is available only for one-sample hyperparameter inputs.');
end

%--------------------------------------------------------------------------
% CUSTOM: Compute output warping function and gradients

% Read hyperparameters
y0 = hyp(1);
k = exp(hyp(2));

% Compute output warping or inverse warping
ywarp = y;
idx = y < y0;
if invflag      % Inverse output warping
    ywarp(idx) = y0 + 1 - (1 + k*y0 - k*y(idx)).^(1/k);
else            % Direct output warping
    delta = (1 + y0 - y(idx));
    deltak = delta.^k;
    ywarp(idx) = y0 - deltak/k + 1/k;
end

if nargout > 1
    % First-order derivative of output warping function in output space
    dwarp_dt = ones(size(y));
    deltakm1 = delta.^(k-1);

    dwarp_dt(idx) = deltakm1;

    if nargout > 2
        % Gradient of output warping function wrt hyperparameters
        dwarp_dtheta = zeros(N,Noutwarp);

        dwarp_dtheta(idx,1) = 1 - deltakm1;                         % y0
        dwarp_dtheta(idx,2) = -deltak.*log(delta) + deltak/k - 1/k; % log(k)

        if nargout > 3
            % Gradient of derivative of output warping function
            d2warp_dthetadt = zeros(N,Noutwarp);

            d2warp_dthetadt(idx,1) = (k-1)*delta.^(k-2);        % y0
            d2warp_dthetadt(idx,2) = k*deltakm1.*log(delta);    % log(k)

        end

    end
end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function  display(grains,varargin)
% standard output

displayClass(grains,inputname(1));

disp(' ')
%disp(char(dynOption(grains)));

% generate phase table
matrix = cell(numel(grains.phaseMap),6);

for ip = 1:numel(grains.phaseMap)

  ind = grains.phaseId == ip;

  % phase
  matrix{ip,1} = num2str(grains.phaseMap(ip)); %#ok<*AGROW>

  % grains
  matrix{ip,2} = int2str(nnz(ind));

  % grains
  matrix{ip,3} = int2str(sum(grains.grainSize(ind)));

  % abort in special cases
  if isempty(grains.CSList{ip})
    continue
  elseif ischar(grains.CSList{ip})
    matrix{ip,4} = grains.CSList{ip};
    continue
  else
    % mineral
    matrix{ip,4} = char(grains.CSList{ip}.mineral);
  end

  % symmetry
  matrix{ip,5} = grains.CSList{ip}.pointGroup;

  % reference frame
  matrix{ip,6} = option2str(grains.CSList{ip}.alignment);

end

% remove empty rows
matrix(accumarray(full(grains.phaseId),1,[size(matrix,1) 1])==0,:) = [];

if ~isempty(grains)
  cprintf(matrix,'-L',' ','-Lc',...
    {'Phase' 'Grains' 'Pixels' 'Mineral'  'Symmetry' 'Crystal reference frame'},...
    '-d','  ','-ic',true);
else
  disp('  no grains here!')
end

disp(' ')

% show boundary and triple points
disp([' ' varlink([inputname(1),'.boundary'],'boundary segments') ': ',int2str(length(grains.boundary))])
disp([' ' varlink([inputname(1),'.innerBoundary'],'inner boundary segments') ': ',int2str(length(grains.innerBoundary))])
disp([' ' varlink([inputname(1),'.triplePoints'],'triple points') ': ',int2str(length(grains.triplePoints))])
disp(' ');

if isempty(grains), return; end

% show properties
disp(char(dynProp(grains.prop),...
  'Id',grains.id,'Phase',grains.phase,'Pixels',grains.grainSize))
disp(' ')

author: mtex-toolbox

repo: mtex

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


text:
function [xdot] = nldetect_anex1_v1_eom(t,x); % ,junk
% Note:  Include ',junk' before eom if this is to be called as an m-file,
%   i.e. [T,Y] = ode45('nldetec_anex1_v1_eom',...)
% Exlude ',junk' if this is to be called as:
%   sol = ode45(@nldetect_anex1_v1_eom, ...)

% State Space representation:
% x = [q; qdot], xdot = [qdot; qdoubledot]
global eom % more efficient to treat as global than to pass in every time.

% Define nonlinear force:
if strcmp(eom.NLType,'bang');
    % Bang (Contact Nonlinearity)
    % Assumes that stiffness of contacted parts is k4mult times higher.
    % Damping is also assumed higher by a factor of c4mult
    del4 = (x(eom.at_ns(1))-x(eom.at_ns(2)));
    if del4 < eom.delcont % No Contact Occurs
        Fnl = eom.kat*del4 ...
            + eom.cfactk*eom.kat*(x(eom.at_ns(1)+eom.Ntot)-x(eom.at_ns(2)+eom.Ntot)); % Add damping terms
    else                % Contact Occurs
        Fnl = eom.kat*eom.k4mult*del4 ...
            + eom.c4mult*eom.cfactk*eom.kat*(x(eom.at_ns(1)+eom.Ntot)-x(eom.at_ns(2)+eom.Ntot)); % Add damping terms
    end
elseif strcmp(eom.NLType,'cubic');
    % Cubic Spring
    del4 = (x(eom.at_ns(1))-x(eom.at_ns(2)));
    Fnl = eom.kat*del4*(1 + eom.katnl*del4^2) ...
        + eom.cfactk*eom.kat*(x(eom.at_ns(1)+eom.Ntot)-x(eom.at_ns(2)+eom.Ntot)); % Add damping terms
else
    error('NLType not recognized');
end

% Define input force - normalized to have unit area (energy) * Afnl
% This is applied to mass 5 below
if t <= eom.tfp
    fext = eom.Afnl*(2*eom.tfp/pi)*sin((pi/eom.tfp)*t);
else
    fext = 0;
end

% Equations of Motion
xdot(eom.dnodes,1) = x(eom.vnodes);
xdot(eom.vnodes,1) = (-eom.Ctot*x(eom.vnodes) - eom.Ktot*x(eom.dnodes) + fext*eom.fext_vec + ...
    Fnl*eom.fnl_vec)./diag(eom.Mtot);

author: Sable

repo: mcbench-benchmarks


text:
function datapt = ma_GASpacecraftMassTasks(stateLogEntry, subTask, celBodyData)
%ma_GASpacecraftMassTasks Summary of this function goes here
%   Detailed explanation goes here

    switch subTask
        case 'fuelOx'
            datapt = stateLogEntry(10);
        case 'monoprop'
            datapt = stateLogEntry(11);
        case 'xenon'
            datapt = stateLogEntry(12);
        case 'dry'
            datapt = stateLogEntry(9);
        case 'total'
            datapt = sum(stateLogEntry(9:12));
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function TestXSum(doSpeed)
% Automatic test: XSum
% This is a routine for automatic testing. It is not needed for processing and
% can be deleted or moved to a folder, where it does not bother.
%
% TestXSum(doSpeed)
% INPUT:
%   doSpeed: Optional logical flag to trigger time consuming speed tests.
%            Default: TRUE. If no speed test is defined, this is ignored.
% OUTPUT:
%   On failure the test stops with an error.
%
% Tested: Matlab 6.5, 7.7, 7.8, WinXP
% Author: Jan Simon, Heidelberg, (C) 2010 matlab.THISYEAR(a)nMINUSsimon.de

% $JRev: R0p V:015 Sum:yqlcwUee8BQZ Date:24-Feb-2010 23:50:29 $
% $License: BSD (see Docs\BSD_License.txt) $
% $File: Published\XSum\TestXSum.m $

% Initialize: ==================================================================
if nargin == 0
   doSpeed = true;
end
if doSpeed
   TestTime = 1;  % [sec]
else
   TestTime = 0.2;
end

% Hello:
ErrID = ['JSim:', mfilename];
whichXSum = which('XSum');

disp(['==== Test XSum:  ', datestr(now, 0), char(10), ...
      'Version: ', whichXSum, char(10)]);

disp('Ask XSum for compilation settings: ------');
Id = evalc('XSum');
fprintf(Id);
disp('-----------------------------------------');

% Methods in most likely (depends on compiler) increasing accuracy:
Method = {'Double', 'Long', 'Kahan', 'Knuth', 'KnuthLong', 'Knuth2', 'QFloat'};

% I've created a version with 384 bit QFloats also, but this needs the LCC v3.8
% compiler, which is not completely compatible with Matlab 2009a:
hasQFloat = any(findstr(lower(Id), 'qfloat: yes'));
if ~hasQFloat
   Method(strcmpi(Method, 'QFloat'))= [];
end

hasLongDouble = any(findstr(lower(Id), 'long double: yes'));
if ~hasLongDouble
   fprintf(['\n!!! The compiler does not support long doubles !!!\n', ...
         '  1. KnuthLong is forwarded to Knuth2.\n', ...
         '  2. XSum(A, [], ''LONG'') can be more precise than SUM(A), if ', ...
         'the compiler \n     accumulates the sum in a 80 bit register, ', ...
         'e.g. Open Watcom 1.8.\n']);
   Method(strcmpi(Method, 'KnuthLong'))= [];
end

nMethod = length(Method);

% Start tests: -----------------------------------------------------------------
% At first proove that some trivial summing work well:
disp([char(10), '== Test small input without rounding errors:']);
disp('  This would stop with an error on problems.');

for iMethod = 1:nMethod
   aMethod = Method{iMethod};
   disp(['-- ', aMethod, ':']);

   S = XSum([], [], aMethod);
   if isempty(S)
      disp('  ok: XSum([])');
   else
      error(ErrID, 'XSum([])');
   end

   S = XSum([], 1, aMethod);
   if isempty(S)
      disp('  ok: XSum([], 1)');
   else
      error(ErrID, 'XSum([], 1)');
   end

   S = XSum([], 2, aMethod);
   if isempty(S)
      disp('  ok: XSum([], 2)');
   else
      error(ErrID, 'XSum([], 2)');
   end

   S = XSum(1, [], aMethod);
   if isequal(S, 1)
      disp('  ok: XSum(1)');
   else
      error(ErrID, 'XSum(1)');
   end

   S = XSum(1, 1, aMethod);
   if isequal(S, 1)
      disp('  ok: XSum(1, 1)');
   else
      error(ErrID, 'XSum(1, 1)')
   end

   S = XSum(1, 2, aMethod);
   if isequal(S, 1)
      disp('  ok: XSum(1, 2)');
   else
      error(ErrID, 'XSum(1, 2)');
   end

   S = XSum([1, 1], [], aMethod);
   if isequal(S, 2)
      disp('  ok: XSum([1, 1], [])');
   else
      error(ErrID, 'XSum([1, 1], [])');
   end

   S = XSum([1, 1], 1, aMethod);
   if isequal(S, [1, 1])
      disp('  ok: XSum([1, 1], 1)');
   else
      error(ErrID, 'XSum([1, 1], 1)');
   end

   S = XSum([1, 1], 2, aMethod);
   if isequal(S, 2)
      disp('  ok: XSum([1, 1], 2)');
   else
      error(ErrID, 'XSum([1, 1], 2)');
   end

   S = XSum([1; 1], [], aMethod);
   if isequal(S, 2)
      disp('  ok: XSum([1; 1], [])');
   else
      error(ErrID, 'XSum([1; 1], [])');
   end

   S = XSum([1; 1], 1, aMethod);
   if isequal(S, 2)
      disp('  ok: XSum([1; 1], 1)');
   else
      error(ErrID, 'XSum([1; 1], 1)');
   end

   S = XSum([1; 1], 2, aMethod);
   if isequal(S, [1; 1])
      disp('  ok: XSum([1; 1], 2)');
   else
      error(ErrID, 'XSum([1; 1], 2)');
   end

   S = XSum([1, 1, 1], [], aMethod);
   if isequal(S, 3)
      disp('  ok: XSum([1, 1, 1], [])');
   else
      error(ErrID, 'XSum([1, 1, 1], [])');
   end

   S = XSum([1, 1, 1], 1, aMethod);
   if isequal(S, [1, 1, 1])
      disp('  ok: XSum([1, 1, 1], 1)');
   else
      error(ErrID, 'XSum([1, 1, 1], 1)');
   end

   S = XSum([1, 1, 1], 2, aMethod);
   if isequal(S, 3)
      disp('  ok: XSum([1, 1, 1], 2)');
   else
      error(ErrID, 'XSum([1, 1, 1], 2)');
   end

   S = XSum([1; 1; 1], [], aMethod);
   if isequal(S, 3)
      disp('  ok: XSum([1; 1; 1], [])');
   else
      error(ErrID, 'XSum([1; 1; 1], [])');
   end

   S = XSum([1; 1; 1], 1, aMethod);
   if isequal(S, 3)
      disp('  ok: XSum([1; 1; 1], 1)');
   else
      error(ErrID, 'XSum([1; 1; 1], 1)');
   end

   S = XSum([1; 1; 1], 2, aMethod);
   if isequal(S, [1; 1; 1])
      disp('  ok: XSum([1; 1; 1], 2)');
   else
      error(ErrID, 'XSum([1; 1; 1], 2)');
   end

   S = XSum([1, 2; 3, 4], [], aMethod);
   if isequal(S, [4, 6])
      disp('  ok: XSum([1,2;3,4], [])');
   else
      error(ErrID, 'XSum([1,2;3,4], [])');
   end

   S = XSum([1, 2; 3, 4], 1, aMethod);
   if isequal(S, [4, 6])
      disp('  ok: XSum([1,2;3,4], 1)');
   else
      error(ErrID, 'XSum([1,2;3,4], 1)');
   end

   S = XSum([1, 2; 3, 4], 2, aMethod);
   if isequal(S, [3; 7])
      disp('  ok: XSum([1,2;3,4], 2)');
   else
      error(ErrID, 'XSum([1,2;3,4], 2)');
   end

   X = reshape(1:2*3*5*7, [2, 3, 5, 7]);
   for iDim = 1:4
      Ssum = sum(X, iDim);
      S    = XSum(X, iDim, aMethod);
      if isequal(S, Ssum)
         fprintf('  ok: XSum([2 x 3 x 5 x 7], %d)\n', iDim);
      else
         error(ErrID, '  ok: XSum([2 x 3 x 5 x 7], %d)\n', iDim);
      end
   end
end

% Measure error compensation - this cannot be "tested", because it is not known,
% what is "correct" or "wrong".
disp([char(10), '== Test error compensation:']);
disp('  Standard truncation error ignores small numbers');
disp('-- [10^k, 1, 1, -10^k] = 2 ?');

Data = [15:21, 104:106];  % 15:120
fprintf('       k =');
fprintf(' %3d', Data);
fprintf('\n');
for iMethod = 0:nMethod
   switch iMethod
      case 0,    aMethod = 'Matlab';
      otherwise, aMethod = Method{iMethod};
   end
   fprintf('%-9s:  ', aMethod);
   for iN = Data;
      X = [10^iN, 1, 1, -10^iN];
      switch iMethod
         case 0,    S = sum(X, 2);
         otherwise, S = XSum(X, 2, aMethod);
      end
      fprintf('%2g  ', S);
   end
   fprintf('\n');
   drawnow;
end

fprintf('\n');
disp('-- [10^2k, 1, 10^k, 1, -10^k, 1, -10^2k] = 3 ?');
Data = [7:11, 15:21, 51:54];   % 4:120;
fprintf('       k =');
fprintf(' %3d', Data);
fprintf('\n');
for iMethod = 0:nMethod
   switch iMethod
      case 0,    aMethod = 'Matlab';
      otherwise, aMethod = Method{iMethod};
   end
   fprintf('%-9s:  ', aMethod);
   for iN = Data;
      X = [10^(2*iN), 1, 10^iN, 1, -10^iN, 1, -10^(2*iN)];
      switch iMethod
         case 0,    S = sum(X, 2);
         otherwise, S = XSum(X, 2, aMethod);
      end
      fprintf('%2g  ', S);
   end
   fprintf('\n');
   drawnow;
end

% ------------------------------------------------------------------------------
nLoop = 100;  % A small averaging
Len   = 1000;
fprintf('\n-- Test with random data:\n')
fprintf('  X = +-rand(1E3, 1) * 10 .^ rand(1E3, 1) * Exp,  exact sum = 0.0\n');
fprintf('  Average over absolute results of %d runs:\n', nLoop);
fprintf('  Smaller results are better!\n');

cVector = 6:32;
if hasQFloat
   cVector = [cVector, 99];
end
fprintf('     Exp = ');
fprintf('       1E%.2d', cVector);
fprintf('\n');

for iMethod = 0:nMethod
   switch iMethod
      case 0,    aMethod = 'Matlab';
      otherwise, aMethod = Method{iMethod};
   end
   fprintf('%-9s: ', aMethod);

   for c = cVector
      S = 0;
      for j = 1:nLoop
         X = RandTestData(Len, c);
         switch iMethod
            case 0,    S = S + abs(sum(X, 1)) / nLoop;
            otherwise, S = S + abs(XSum(X, 1, aMethod)) / nLoop;
         end
      end
      fprintf('%11.3G', S);
   end
   fprintf('\n');
   drawnow;
end

% ------------------------------------------------------------------------------
fprintf('\n-- Compare accuracy of results:\n')
nLoop = 100;  % A small averaging
if hasQFloat && 0
   c = 120;
else
   c = 48;   % > 16, e.g. with 48 all outputs have a similar format
end

fprintf('  X = +-rand(N, 1) * 10 .^ rand(N, 1) * Exp,  exact sum = 0.0\n');
fprintf('  Average over absolute results of %d runs:\n', nLoop);
fprintf('  Smaller results are better!\n');

for LenExp = 3:5
   fprintf('  N = 1E%d,  Exp = 1E%.2d\n', LenExp, c);
   Len = 10 ^ LenExp;

   TestM = zeros(Len, nLoop);
   for iLoop = 1:nLoop
      TestM(:, iLoop) = RandTestData(Len, c);
   end

   for iMethod = -1:nMethod
      switch iMethod
         case -1,    aMethod = 'SUM';
         case  0,    aMethod = 'SUM(SORT)';
         otherwise,  aMethod = Method{iMethod};
      end

      S = 0;
      for iLoop = 1:nLoop
         X = TestM(:, iLoop);
         switch iMethod
            case -1,   S = S + abs(sum(X, 1)) / nLoop;
            case  0,   S = S + abs(sum(sort(X), 1)) / nLoop;
            otherwise, S = S + abs(XSum(X, 1, aMethod)) / nLoop;
         end
      end
      fprintf('    %-9s: %.6G\n', aMethod, S);
      drawnow;
   end
end

% RANDN:
nLoop = 20;
fprintf(['\n  RANDN: Absolute results, sum estimated by Knuth2, ', ...
      'average over %d runs:\n'], nLoop);
fprintf('  Smaller results are better!\n');
fprintf('  0 must be replied for Knuth2, but even Knuth should be exact.\n');

for LenExp = 4:6
   fprintf('  X = RANDN(1, 1E%d):\n', LenExp);
   Len = 10 ^ LenExp;

   MethodSum = zeros(1, nMethod + 2);
   for iLoop = 1:20
      TestM    = randn(Len, 1);
      EstimSum = XSum(TestM, 1, 'Knuth2');

      for iMethod = -1:nMethod
         switch iMethod
            case -1
               S = abs(sum(TestM, 1) - EstimSum);
            case 0
               S = abs(sum(sort(TestM), 1) - EstimSum);
            otherwise
               aMethod = Method{iMethod};
               S = abs(XSum(TestM, 1, aMethod) - EstimSum);
         end
         MethodSum(iMethod + 2) = MethodSum(iMethod + 2) + S / nLoop;
      end
   end

   for iMethod = -1:nMethod
      switch iMethod
         case -1,   aMethod = 'SUM';
         case 0,    aMethod = 'SUM(SORT)';
         otherwise, aMethod = Method{iMethod};
      end
      fprintf('    %-9s: %.8G\n', aMethod, MethodSum(iMethod + 2));
      drawnow;
   end
end

% Speed: -----------------------------------------------------------------------
fprintf('\n-- Compare times:\n');
LenList = [1E3, 1E4, 1E5, 1E6, 1E7];

fprintf('  SUM(RAND(N, 1)):\n');
fprintf('       N =');
fprintf('     1E%d', round(log10(LenList)));
fprintf('\n');
Time_SUM = zeros(1, length(LenList));
for iMethod = -1:nMethod
   switch iMethod
      case -1,   aMethod = 'SUM';
      case 0,    aMethod = 'SUM(SORT)';
      otherwise, aMethod = Method{iMethod};
   end
   fprintf('%-9s:', aMethod);

   Time_Method = zeros(1, length(LenList));

   for iLen = 1:length(LenList)
      Len = LenList(iLen);
      X   = rand(Len, 1);

      % Get number of loops:
      N     = 0;
      iTime = cputime;
      while cputime - iTime < 0.5
         S = sum(X);  clear('S');
         N = N + 1;
      end
      nLoop = N / (cputime - iTime);  % Loops per second
      nLoop = max(1, ceil(nLoop * TestTime));

      switch iMethod
         case -1
            tic;
            for iLoop = 1:nLoop
               S = sum(X);  clear('S');
            end
            nPerSec = nLoop / toc;
            Time_SUM(iLen) = nPerSec;

         case 0
            tic;
            mLoop = nLoop / 10;  % Reduce test time!
            for iLoop = 1:mLoop
               S = sum(sort(X));  clear('S');
            end
            nPerSec = mLoop / toc;

         otherwise
            if strcmpi(aMethod, 'qfloat')  % QFloat is very slow...
               mLoop = nLoop / 10;
            else
               mLoop = nLoop;
            end

            tic;
            for iLoop = 1:mLoop
               S = XSum(X, 1, aMethod);  clear('S');   %#ok<*NASGU>
            end
            nPerSec = mLoop / toc;
      end
      Time_Method(iLen) = nPerSec;

      PrintLoop(nPerSec);
   end
   fprintf('   loops/sec\n');

   if iMethod ~= -1
      fprintf('          ');
      fprintf('%8.1f', Time_SUM ./ Time_Method);
      fprintf('   time for SUM\n');
   end
end

% Goodbye:
fprintf('\nXSum seems to work fine.\n');

return;

% ******************************************************************************
function PrintLoop(N)
if N > 10
   fprintf('  %6.0f', N);
else
   fprintf('  %6.1f', N);
end

return;

% ******************************************************************************
function X = RandTestData(n, c)
% Test vector with exact sum is 0.
% INPUT:  n: Length of output vector, must be even.
%         c: The output has a maximum absolute value of 10^c.
% OUTPUT: X: Random [1 x n] double vector with exact sum 0.
%
% Simple approach: Create random numbers with random exponent and append a copy
% with opposite sign. Finally the vector is mixed.
% Data created with this method are definitely not valid to test the SUM over
% sorted absolute values

m = n / 2;
X = rand(m, 1) .* 10 .^ ceil(rand(m, 1) * c);
X = [X; -X];
X = randmix(X);

return;

% ******************************************************************************
function X = randmix(X)
% Faster RANDPERM, Knuth's shuffle
% Author: Jan Simon, Heidelberg, (C) 2010 matlab.THISYEAR(a)nMINUSsimon.de

for i = 1:numel(X)
   w    = ceil(rand * i);
   t    = X(w);
   X(w) = X(i);
   X(i) = t;
end

return;


text:
% AUTHOR:         Aaron Nicolson
% AFFILIATION:    Signal Processing Laboratory, Griffith University
%
% This Source Code Form is subject to the terms of the Mozilla Public
% License, v. 2.0. If a copy of the MPL was not distributed with this
% file, You can obtain one at http://mozilla.org/MPL/2.0/.

clear all; close all; clc;

%% GET MATLAB_FEAT REPOSITORY
addpath('~/Dropbox/GitHub/matlab_feat/feat')
addpath('./deepxi')

%% PARAMETERS
T_d = 32; % window duration (ms).
T_s = 16; % window shift (ms).
f_s = 16000; % sampling frequency (Hz).
s.N_d = round(f_s*T_d*0.001); % window duration (samples).
s.N_s = round(f_s*T_s*0.001); % window shift (samples).
s.f_s = f_s; % sampling frequency (Hz).
s.NFFT = 2^nextpow2(s.N_d); % frequency bins (samples).
d = s; x = s;
SNR_avg = -5:5:15; % SNR levels used to compute average SD level.

%% DIRECTORIES
s.dir = '/home/aaron/set/deep_xi_test_set/test_clean_speech';
d.dir = '/home/aaron/set/deep_xi_test_set/test_noise';
gamma_hat_dir = input('gamma_hat path:', 's');
gamma_hat_dir_split = strsplit(gamma_hat_dir, '/');
ver = [gamma_hat_dir_split{end-2}, '_', gamma_hat_dir_split{end-1}];

%% FILE LISTS
gamma_hat_paths = dir([gamma_hat_dir, '/*.mat']); % noise file paths.

results = MapNested();
noise_src_set = {};
SNR_set = {};
for i = 1:length(gamma_hat_paths)

    load([gamma_hat_paths(i).folder, '/', gamma_hat_paths(i).name])

    if any(isnan(gamma_hat(:))) || any(isinf(gamma_hat(:)))
        error('NaN or Inf value in gamma_hat: %s.', gamma_hat_paths(i).name)
    end

    split_basename = strsplit(gamma_hat_paths(i).name,'_');
    noise_src = split_basename{end-1};
    SNR = split_basename{end};
    clean_speech = extractBefore(gamma_hat_paths(i).name, ['_', noise_src, '_', SNR]);
    SNR = SNR(1:end-6);

    s.wav = audioread([s.dir, '/', clean_speech, '_', noise_src, '.wav']); % clean speech.
    d.src = audioread([d.dir, '/', clean_speech, '_', noise_src, '.wav']); % noise.
    [x.wav, d.wav] = add_noise(s.wav, d.src(1:length(s.wav)), str2double(SNR)); % noisy speech.

    x = analysis_stft(x, 'polar'); % noisy-speech STMS.
    d = analysis_stft(d, 'polar'); % noise STMS.

    gamma = (x.STMS.^2)./(d.STMS.^2); % instantaneous a posteriori SNR.

    gamma_hat = gamma_hat(1:size(gamma, 1), :);
    D = spectral_distortion(gamma, gamma_hat);

    if any(isnan(D(:))) || any(isinf(D(:)))
        error('NaN or Inf value in D.')
    end

    if ~any(strcmp(SNR_set, SNR))
        SNR_set{end+1} = SNR;
    end

    if ~any(strcmp(noise_src_set, noise_src))
        noise_src_set{end+1} = noise_src;
    end

    if results.isKey(noise_src, SNR)
        results(noise_src, SNR) = [D; results(noise_src, SNR)];
    else
        results(noise_src, SNR) = D;
    end
    clc;
    fprintf('%.2f%%\n', 100*i/length(gamma_hat_paths));
end

% there has to be a better way to do this.
for i=1:length(SNR_set)
    SNR_set{i} = str2num(SNR_set{i});
end
tmp = sort(cell2mat(SNR_set));
for i=1:length(SNR_set)
    SNR_set{i} = num2str(tmp(i));
end

res_dir = 'log/results/spectral_distortion_gamma';
if ~exist(res_dir, 'dir')
    mkdir(res_dir)
end

fileID = fopen([res_dir, '/', ver, '.csv'],'w');
fprintf(fileID, 'noise, snr_db, SD\n');
avg = [];
for i=1:length(noise_src_set)
    for j=1:length(SNR_set)
        D = mean(results(noise_src_set{i}, SNR_set{j}));
        fprintf('%s, %s: %.2f.\n', noise_src_set{i}, SNR_set{j}, D);
        fprintf(fileID, '%s, %s, %.2f\n', noise_src_set{i}, SNR_set{j}, D);
        if ismember(j, SNR_avg)
            avg = [avg; results(noise_src_set{i}, SNR_set{j})];
        end
    end
end
fclose(fileID);

avg_path = [res_dir, '/average.csv'];
if ~exist(avg_path, 'file')
    fileID = fopen(avg_path, 'w');
    fprintf(fileID, 'ver, SD\n');
    fclose(fileID);
end
fileID = fopen(avg_path, 'a');
fprintf(fileID, '%s, %.2f\n', ver, mean(avg));
fclose(fileID);

author: anicolson

repo: DeepXi

save_path: github-repos/MATLAB/anicolson-DeepXi


text:
function Y = invariantMultiply_new( A, X)
%
%           Y = invariantMultiply( A, X );
%
%  This function computes the multiplication of a spatially invariant
%  point spread function (PSF) times an image
%                y = A*x
%
%  Here we use a Kronecker product approach to better handle reflexive
%  boundary conditions
%
%  Input:
%            A  -  psfMatrix object
%            X  -  array containing the image to which the psfMatrix
%                  is to be multiplied.
%
%  Output:
%            Y  -  contains the result after PSF multiplication.
%

%  J. Nagy  7/5/2016

imsize = size( X );
psfMatData = A.matdata;

switch A.boundary
    case 'periodic'
        if A.transpose
            Y = real(ifft2(conj(psfMatData).*fft2(X)));
        else
            Y = real(ifft2(psfMatData.*fft2(X)));
        end
    case 'zero'
        p = A.p;
        padSize = [size(psfMatData,1),size(psfMatData,2)] - size(X);
        Xpad = padarray(X, padSize, 'post');
        Xtilde = fft2(Xpad);
        Xtilde = Xtilde(:,p(:,1));
        if A.transpose
            Y = conj(psfMatData(:,:,1)).*Xtilde;
        else
            Y = psfMatData(:,:,1).*Xtilde;
        end
        Y = real(ifft2(Y));
        Y = Y(:, p(:,1));
        Y = Y(1:imsize(1), 1:imsize(2));
    case 'reflexive'
        p = A.p;
        padSize = [size(psfMatData,1),size(psfMatData,2)] - size(X);
        Xpad = padarray(X, padSize, 'post');
        Xtilde = fft2(Xpad);
        Xtilde = Xtilde(:,p(:,1));
        if A.transpose
            term1 = conj(psfMatData(:,:,1)).*Xtilde;
            term2 = conj(psfMatData(:,:,2)).*Xtilde(:,p(:,1));
            term3 = conj(psfMatData(:,:,3)).*Xtilde(p(:,2),:);
            term4 = conj(psfMatData(:,:,4)).*Xtilde(p(:,2),p(:,1));
        else
            term1 = psfMatData(:,:,1).*Xtilde;
            term2 = psfMatData(:,:,2).*Xtilde;
            term2 = term2(:,p(:,1));
            term3 = psfMatData(:,:,3).*Xtilde;
            term3 = term3(p(:,2),:);
            term4 = psfMatData(:,:,4).*Xtilde;
            term4 = term4(p(:,2),p(:,1));
        end
        Y = real(ifft2(term1 + term2 + term3 + term4));
        Y = Y(:,p(:,1));
        Y = Y(1:imsize(1), 1:imsize(2));
    otherwise
        error('still working on this')
end

repo: IRtools


text:
function Dataset2XLS(DATA, filename)
%DATASET2XLS - Save Dataset format to EXCEL File
% This function wil help you to save simulation results to Excel files.
% INPUTS:
% DATA - The data to be saved, in DataSet format (Simulink).
% Filename - optional file name to save in. if not supplied, will use
% example_Dataset.xlsx.
% usage:
% 1. Run a Simulink model, and save the result in a Dataset format. The
% defaults filename is 'logsout'.
% 2. Run this function: Dataset2XLS(logsout, 'test.xls');
%

%% Check inputs:
if (nargin < 1 || nargin >3)
    error('Wrong number of arguments. Enter Dataset parameter Name');
end
if nargin < 2
    filename = 'example_Dataset.xlsx';
    warning('### No filename supplied, using example_Dataset.xlsx instead ###');
end
if ~ischar(filename)
    error('Illegal file name');
end
if isempty(DATA)
    error('Illegal Data to save');
end
if ~isa(DATA,'Simulink.SimulationData.Dataset')
    error('Data is not in Simulink.SimulationData.Dataset Format');
end
if ~exist(filename,'file')
    existingFile = 0;
else
    existingFile = 1;
end

%% Write data to file:
disp('### Writing to file ###')
warning('off'); %#ok<*WNOFF>
for i = 1:DATA.numElements
    fprintf('### Writing element #%d ###\n',i);
    locSig = DATA.getElement(i);
    stat = xlswrite(filename, [locSig.Values.Time locSig.Values.Data], locSig.Name);
    if ~stat
        error('Problem writing to file, check if it''s not open');
    end
end
warning('on'); %#ok<*WNON>
fprintf('### Finished writing to file %s ###\n',filename)

%% Delete unwanted excel sheets (If file did not exist beforehand):
if ~existingFile
    disp('### Removing extra sheets ###')
    % Retrieve sheet names
    [~, sheetNames] = xlsfinfo(filename);
    % Open Excel as a COM Automation server
    hExcel = actxserver('Excel.Application');
    % Open Excel workbook
    Workbook = hExcel.Workbooks.Open(fullfile(pwd,filename));
    % Clear the content of the sheets (from the second onwards)
    cellfun(@(x) hExcel.ActiveWorkBook.WorkSheets.Item(x).Delete, sheetNames(1:3));
    % Now save/close/quit/delete
    Workbook.Save;
    hExcel.Workbook.Close;
    hExcel.Quit;
    hExcel.delete;
end


text:
% ismember_bc - ismember backward compatible with Matlab versions prior to 2013a

% Copyright (C) 2013 Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [C,IA] = ismember_bc(A,B,varargin)

errorFlag = error_bc;

v = version;
indp = find(v == '.');
v = str2num(v(1:indp(2)-1));
if v > 7.19, v = floor(v) + rem(v,1)/10; end

if nargin > 2
    ind = strmatch('legacy', varargin);
    if ~isempty(ind)
        varargin(ind) = [];
    end
end

if v >= 7.14
    [C,IA] = ismember(A,B,varargin{:},'legacy');
    if errorFlag
        [C2,IA2] = ismember(A,B,varargin{:});
        if (~isequal(C, C2) || ~isequal(IA, IA2))
            warning('backward compatibility issue with call to ismember function');
        end
    end
else
    [C,IA] = ismember(A,B,varargin{:});
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function value = p18_f ( dim_num, point_num, x )

%*****************************************************************************80
%
%% P18_F evaluates the integrand for problem 18.
%
%  Discussion:
%
%    This is the characteristic function of the interior of the
%    N sphere of radius R and center Z, to be integrated within the
%    unit hypercube [0,1]^N.  If the user picks a combination of R
%    and Z that causes the volume of the sphere to lie at least
%    partially outside the unit hypercube, the formula for the
%    exact integral will no longer be correct.
%
%  Dimension:
%
%    DIM_NUM is arbitrary.
%
%  Region:
%
%    0 <= X(1:DIM_NUM) <= 1
%
%  Integral Parameters:
%
%    R defaults to 0.50.
%    You can change R by calling P18_R8.
%
%    Z(1:DIM_NUM) defaults to (0.5,0.5,...0.5).
%    You can change Z by calling P18_R8VEC.
%
%  Integrand:
%
%    f(x) = 1 if X(1:DIM_NUM) is less than R from Z(1:DIM_NUM),
%           0 otherwise.
%
%  Exact Integral:
%
%    The volume of the DIM_NUM sphere of radius R.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the evaluation points.
%
%    Output, real VALUE(POINT_NUM), the integrand values.
%
  z = [];
  z = p18_r8vec ( 'G', 'Z', dim_num, z );

  r = 0.0;
  r = p18_r8 ( 'G', 'R', r );

  value(1:point_num) = 0.0;

  for point = 1 : point_num

    d = sqrt ( sum ( ( x(1:dim_num,point) - z(1:dim_num)' ).^2 ) );

    if ( d <= r )
      value(point) = 1.0;
    else
      value(point) = 0.0;
    end

  end

  p18_i4 ( 'I', '#', point_num );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function boundary = p02_boundary_nearest ( m, n, point )

%*****************************************************************************80
%
%% P02_BOUNDARY_NEAREST returns a nearest boundary point in problem 02.
%
%  Discussion:
%
%    The given input point need not be inside the region.
%
%    In some cases, more than one boundary point may be "nearest",
%    but only one will be returned.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real POINT(M,N), the coordinates
%    of the points.
%
%    Output, real BOUNDARY(M,N), points on the boundary
%    that are nearest to each point.
%
  center = [ 0.0, 0.0 ];
  r1 = 1.0;
  r2 = 0.4;

  for j = 1 : n

    if ( point(1:2,j) == center(1:2)' )

      boundary(1,j) = r2;
      boundary(2,j) = center(2);

    else

      r = sqrt ( ( point(1,j) - center(1) ).^2 ...
               + ( point(2,j) - center(2) ).^2 );

      if ( r1 - r < r - r2 )
        boundary(1,j) = center(1) + r1 * ( point(1,j) - center(1) ) / r;
        boundary(2,j) = center(2) + r1 * ( point(2,j) - center(2) ) / r;
      else
        boundary(1,j) = center(1) + r2 * ( point(1,j) - center(1) ) / r;
        boundary(2,j) = center(2) + r2 * ( point(2,j) - center(2) ) / r;
      end

    end

  end

  return
end


text:
function [Torr] = ftH2O2Torr(ftH2O)
% Convert pressure from feet of water column at 4 degrees to torr
% Chad Greene 2012
Torr = ftH2O*22.4198;


text:
function M = diag(f)
%DIAG   Multiplication operator.
%   M = DIAG(F) returns a multiplication operator M so that M*u is equivalent to
%   F.*u.
%
%   This method is mainly provided for backwards compatibility.

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

M = operatorBlock.mult(f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ a, b, c ] = line_exp2imp_2d ( p1, p2 )

%*****************************************************************************80
%
%% LINE_EXP2IMP_2D converts an explicit line to implicit form in 2D.
%
%  Discussion:
%
%    The explicit form of a line in 2D is:
%
%      ( P1, P2 ) = ( (X1,Y1), (X2,Y2) ).
%
%    The implicit form of a line in 2D is:
%
%      A * X + B * Y + C = 0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2,1), P2(2,1), two points on the line.
%
%    Output, real A, B, C, the implicit form of the line.
%

%
%  Take care of degenerate cases.
%
  if ( p1(1:2,1) == p2(1:2,1) )
    fprintf ( 1, '\n' );
    fprintf ( 1,  'LINE_EXP2IMP_2D - Fatal error!\n' );
    fprintf ( 1,  '  P1 = P2\n' );
    fprintf ( 1,  '  P1 = %f  %f\n', p1(1:2,1) );
    fprintf ( 1,  '  P2 = %f  %f\n', p2(1:2,1) );
    error ( 'LINE_EXP2IMP_2D - Fatal error!' );
  end

  a = p2(2,1) - p1(2,1);
  b = p1(1,1) - p2(1,1);
  c = p2(1,1) * p1(2,1) - p1(1,1) * p2(2,1);

  norm = a * a + b * b + c * c;

  if ( 0.0 < norm )
    a = a / norm;
    b = b / norm;
    c = c / norm;
  end

  if ( a < 0.0 )
    a = -a;
    b = -b;
    c = -c;
  end

  return
end


text:
function s=openRaw3ViewWindow
%
% Opens a volume 3-view (sag,axi,cor) window and initializes the corresponding data structure.
%
% VOLUME is a cell array of volume structures.
% s is the index of the new one.
%
% This function should ONLY be called by either:
% openVolumeWindow or openGrayWindow.
%
% Modifications:
%
% djh and baw, 7/98
% djh, 2/99 modified to enable the window to be opened under
% either volume or gray mode
%
% djh, 4/99
% - Eliminate overlayClip sliders
% - Added mapWin sliders to show overlay only for pixels with parameter
%   map values that are in the appropriate range.
% bw, 12/29/00
% - scan slider instead of buttons
% djh 2/13/2001
% - open multiple volume windows simultaneously
% ras 10/28/02
% - adds 'Rory' menu for Rory's analysis; local version
% ras 03/06/03
% revamp: now makes a '3-view' showing sagittal, coronal, and axial views
% all at once, instead of paging through the slices

% Make sure the global variables exist
mrGlobals

disp('Initializing Volume view')

% s is the index of the new volume structure.
s = getNewViewIndex(VOLUME); %#ok<NODEF>

% Set name, viewType, & subdir
VOLUME{s} = [];
VOLUME{s} = viewSet(VOLUME{s}, 'name', ['VOLUME{',num2str(s),'}']);
VOLUME{s} = viewSet(VOLUME{s}, 'viewType',  'Volume');
VOLUME{s} = viewSet(VOLUME{s}, 'subdir', 'Volume');

% Refresh function, gets called by refreshScreen
VOLUME{s} = viewSet(VOLUME{s}, 'refreshFn', 'volume3View');

%%%%%%%%%%%%%%%%%%
% Load Anatomies %
%%%%%%%%%%%%%%%%%%
VOLUME{s} = loadAnat(VOLUME{s});
%TODO: Change the functionality of loadAnat to bring it in line with the
%inplane load anat functionality

%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize data slots %
%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize slots for co, amp, and ph
VOLUME{s}.co = [];
VOLUME{s}.amp = [];
VOLUME{s}.ph = [];
VOLUME{s}.map = [];
VOLUME{s}.mapName = '';
VOLUME{s}.mapUnits = '';

% Initialize slots for tSeries
VOLUME{s}.tSeries = [];
VOLUME{s}.tSeriesScan = NaN;
VOLUME{s}.tSeriesSlice = NaN;

% Initialize ROIs
VOLUME{s}.ROIs = [];
VOLUME{s}.selectedROI = 0;

% Initialize curDataType / curScan
VOLUME{s}.curDataType = 1;
VOLUME{s}.curScan = 1;

% Initialize location in volume space
dims = viewGet(VOLUME{s},'Size');
VOLUME{s}.loc = round(dims./2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make displayModes and color maps %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize displayModes
VOLUME{s} = resetDisplayModes(VOLUME{s});
VOLUME{s}.ui.displayMode = 'anat';

%%%%%%%%%%%%%%%%%%%
% Open the window %
%%%%%%%%%%%%%%%%%%%
% Figure number for volume view window
figName = sprintf('Volume %i', s);
if isfield(mrSESSION,'sessionCode') && isfield(mrSESSION,'description'),
  figName = sprintf('%s: %s  [%s]', figName, mrSESSION.sessionCode, ...
                    mrSESSION.description);
end

VOLUME{s}.ui.figNum = figure('Name', figName, ...
                             'MenuBar', 'none', ...
                             'Units', 'Normalized',...
                             'Color', [.9 .9 .9], ...
                             'Position', [0.25 0.25 0.5 0.5],...
                             'NumberTitle', 'off',  ...
                             'Tag', ['3VolumeWindow: ',VOLUME{s}.name]);

% Handle for volume view window
VOLUME{s}.ui.windowHandle = gcf;

% set up 3 axes, for each view of the volume:
VOLUME{s}.ui.sagAxesHandle = axes('position',[.15 .5 .4 .4]);
VOLUME{s}.ui.corAxesHandle = axes('position',[.5 .5 .4 .4]);
VOLUME{s}.ui.axiAxesHandle = axes('position',[.5 .15 .4 .4]);

% The 'main' axis will arbitarily refer to the axial window
VOLUME{s}.ui.mainAxisHandle = VOLUME{s}.ui.axiAxesHandle;

% ras 05/06: for faster updating, keep track of image/obj handles:
VOLUME{s}.ui.underlayHandle = []; % for ML7+, will make overlays
VOLUME{s}.ui.overlayHandle = [];  % a separate img, and use alpha maps
VOLUME{s}.ui.roiHandles = {};

% Set minColormap property so there's potentially room for 128
% colors
%set(VOLUME{s}.ui.windowHandle, 'minColormap', 128);

% Set closeRequestFcn so we can clean up when the window is closed
set(gcf, 'CloseRequestFcn', 'closeVolumeWindow');

% Set selectedVOLUME when click in this window
set(gcf, 'WindowButtonDownFcn', ['selectedVOLUME =',num2str(s),';']);

%%%%%%%%%%%%%
% Add Menus %
%%%%%%%%%%%%%
disp('Attaching menus')

VOLUME{s} = filesMenu(VOLUME{s});
VOLUME{s} = editMenu(VOLUME{s});
VOLUME{s} = windowMenu(VOLUME{s});
VOLUME{s} = analysisMenu(VOLUME{s});
VOLUME{s} = viewMenu(VOLUME{s});
VOLUME{s} = roiMenu(VOLUME{s});
VOLUME{s} = plotMenu(VOLUME{s});
VOLUME{s} = colorMenu(VOLUME{s});
VOLUME{s} = xformVolumeMenu(VOLUME{s});
VOLUME{s} = grayMenu(VOLUME{s});
VOLUME{s} = eventMenu(VOLUME{s});
VOLUME{s} = helpMenu(VOLUME{s}, 'Volume');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Add Annotation String; move to nicer location %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
VOLUME{s} = makeAnnotationString(VOLUME{s});

%%%%%%%%%%%%%%%%%
% Add Color Bar %
%%%%%%%%%%%%%%%%%
% % Make color bar and initialize it to 'off'
VOLUME{s}.ui.colorbarHandle = makeColorBar(VOLUME{s});
setColorBar(VOLUME{s},'off');
VOLUME{s}.ui.cbarRange = [];

%%%%%%%%%%%%%%%%%%%
% Add buttons     %
%%%%%%%%%%%%%%%%%%%
disp('Attaching buttons')

% Buttons and editable text fields for choosing slice number and
% slice orientation
VOLUME{s} = make3ViewSliceUI(VOLUME{s});
setCurSliceOri(VOLUME{s},1); % default to axials

% Buttons for choosing gray mode vs volume mode
VOLUME{s} = makeGrayVolButtons(VOLUME{s});

% add some zoom buttons
VOLUME{s} = makeZoomButtons(VOLUME{s});

% now add a checkbox to toggle the crosshairs
xHairColor = [1 .5 .5];
if ispref('VISTA', 'xHairColor')
    xHairColor = getpref('VISTA', 'xHairColor', xHairColor);
end
callback = [VOLUME{s}.name ' = toggle3ViewCrossHairs(' VOLUME{s}.name ', gcbo);'];
htmp = uicontrol('Style', 'checkbox', 'String', 'Crosshairs', ...
                 'Value', 0, 'Callback', callback, ...
                 'BackgroundColor', get(gcf,'Color'), ...
                 'ForegroundColor', xHairColor ,...
                 'Units', 'Normalized', 'Position', [0 0.1 0.1 0.05]);
VOLUME{s}.ui.xHairToggleHandle = htmp;
VOLUME{s}.ui.crosshairs = 0; % crosshairs initialized off

%%%%%%%%%%%%%%%%%%%
% Add popup menus %
%%%%%%%%%%%%%%%%%%%
disp('Attaching popup menus')

VOLUME{s} = makeROIPopup(VOLUME{s});
VOLUME{s} = makeDataTypePopup(VOLUME{s});

%%%%%%%%%%%%%%%
% Add sliders %
%%%%%%%%%%%%%%%
disp('Attaching sliders')

% Scan number slider
w = 0.12; h = 0.03; l = 0; b = 0.95;
VOLUME{s} = makeSlider(VOLUME{s},'scan',[],[l b w h]);
VOLUME{s} = initScanSlider(VOLUME{s},1);
VOLUME{s} = selectDataType(VOLUME{s},VOLUME{s}.curDataType);

% correlation threshold:
VOLUME{s} = makeSlider(VOLUME{s},'cothresh',[0,1],[.85,.85,.15,.03]);
setCothresh(VOLUME{s},0);

% phase window:
VOLUME{s} = makeSlider(VOLUME{s},'phWinMin',[0,2*pi],[.85,.75,.15,.03]);
VOLUME{s} = makeSlider(VOLUME{s},'phWinMax',[0,2*pi],[.85,.65,.15,.03]);
setPhWindow(VOLUME{s},[0 2*pi]);

% parameter map window:
VOLUME{s} = makeSlider(VOLUME{s},'mapWinMin',[0,1],[.85,.55,.15,.03]);
VOLUME{s} = makeSlider(VOLUME{s},'mapWinMax',[0,1],[.85,.45,.15,.03]);
setMapWindow(VOLUME{s},[0 1]);

% Brightness / Contrast: (Replaces older Anat Clip).
VOLUME{s} = makeBrightnessSlider(VOLUME{s},[.85 .2 .15 .03]);
VOLUME{s} = makeSlider(VOLUME{s},'contrast',[0 1],[.85 .1 .15 .03]);
VOLUME{s} = viewSet(VOLUME{s},'brightness',0.5);
VOLUME{s} = viewSet(VOLUME{s},'contrast',0.05);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize display parameters %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize image field
VOLUME{s}.ui.image = [];

% New ROI display fields (ras 01/07)
VOLUME{s}.ui.showROIs = -2;        % list of ROIs to show (0 = hide, -1 = selected)
VOLUME{s}.ui.roiDrawMethod = 'perimeter'; % can be 'boxes', 'perimeter', 'patches'
VOLUME{s}.ui.filledPerimeter = 0; % filled perimeter toggle

disp('Done initializing Volume view')

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ac] = printResult(X, label, K, kmeansFlag)

if kmeansFlag == 1
    indic = litekmeans(X, K, 'Replicates',20);
else
    [~, indic] = max(X, [] ,2);
end
result = bestMap(label, indic);
[ac, nmi_value, cnt] = CalcMetrics(label, indic);
disp(sprintf('ac: %0.4f\t%d/%d\tnmi:%0.4f\t', ac, cnt, length(label), nmi_value));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function fun = makelinefun(x1,y1,x2,y2,o)
%ref http://stackoverflow.com/questions/13209373/matlab-straight-line-between-2-points-with-n-points-between
if nargin < 5
    o = 1;
end
if o == 2,
    fun  = @(N) deal(linspace(x1,x2,N), linspace(y1,y2,N));
else
	fun  = @(N) [linspace(x1,x2,N) ; linspace(y1,y2,N)];
end

%
% USAGE:
%
% f = makelinefun(0,0,6,9);
% xy = f(4)
%
%   xy =
%       0     2     4     6
%       0     3     6     9
%
%       f = makelinefun(0,0,6,9);
% [x,y] = f(4)
%
%   x =
%       0     2     4     6
%   y =
%       0     3     6     9
%


text:
%% Copyright (C) 2014, 2016, 2019 Colin B. Macdonald
%%
%% 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
%% @defmethod  @@sym repmat (@var{A}, @var{n}, @var{m})
%% @defmethodx @@sym repmat (@var{A}, [@var{n} @var{m}])
%% Build symbolic block matrices.
%%
%% Examples:
%% @example
%% @group
%% repmat([1 2 sym(pi)], 2, 3)
%%   @result{} (sym 2×9 matrix)
%%       ⎡1  2  π  1  2  π  1  2  π⎤
%%       ⎢                         ⎥
%%       ⎣1  2  π  1  2  π  1  2  π⎦
%%
%% repmat(sym(pi), [1 3])
%%   @result{} (sym) [π  π  π]  (1×3 matrix)
%% @end group
%% @end example
%%
%% @seealso{@@sym/vertcat, @@sym/horzcat}
%% @end defmethod

function B = repmat(A, n, m)

  if (nargin == 2)
    m = n(2);
    n = n(1);
  elseif (nargin == 3)
    % no-op
  else
    print_usage ();
  end

  cmd = { '(A, n, m) = _ins'
          'if n == 0 or m == 0:'
          '    return sp.Matrix(n, m, [])'
          'if A is None or not A.is_Matrix:'
          '    A = sp.Matrix([A])'
          'L = [A]*m'
          'B = sp.Matrix.hstack(*L)'
          'L = [B]*n'
          'B = sp.Matrix.vstack(*L)'
          'return B' };

  B = pycall_sympy__ (cmd, sym(A), int32(n), int32(m));

end

%!test
%! % simple
%! syms x
%! A = [x x x; x x x];
%! assert (isequal (repmat(x, 2, 3), A))

%!test
%! % block cf double
%! A = [1 2 3; 4 5 6];
%! B = sym(A);
%! C = repmat(A, 2, 3);
%! D = repmat(B, 2, 3);
%! assert (isequal (C, D))

%!test
%! % empty
%! A = repmat(sym([]), 2, 3);
%! assert (isempty(A));
%! assert (isequal (size(A), [0 0]))

%!test
%! % more empties
%! A = repmat(sym(pi), [0 0]);
%! assert (isequal (size(A), [0 0]))
%! A = repmat(sym(pi), [0 3]);
%! assert (isequal (size(A), [0 3]))
%! A = repmat(sym(pi), [2 0]);
%! assert (isequal (size(A), [2 0]))

repo: octsympy


text:
function kam = KAM(ebsd,varargin)
% intragranular average misorientation angle per orientation
%
% Syntax
%
%   plot(ebsd,ebsd.KAM ./ degree)
%
%   % ignore misorientation angles > threshold
%   kam = KAM(ebsd,'threshold',10*degree);
%   plot(ebsd,kam./degree)
%
%   % ignore grain boundary misorientations
%   [grains, ebsd.grainId] = calcGrains(ebsd)
%   plot(ebsd, ebsd.KAM./degree)
%
%   % consider also second order neigbors
%   kam = KAM(ebsd,'order',2);
%   plot(ebsd,kam./degree)
%
% Input
%  ebsd - @EBSD
%
% Options
%  threshold - ignore misorientation angles larger then threshold
%  order     - consider neighbors of order n
%  max       - take not the mean but the maximum misorientation angle
%
% See also
% grain2d.GOS

%weights = get_option(varargin,'weights',ones(3));

if check_option(varargin,'max')
  fun = @(a,b) nanmax(a,[],b);
else
  fun = @(a,b) nanplus(a,b);
end

weights = get_option(varargin,'weights',[]);

if isempty(weights)

  % get order
  order = get_option(varargin,'order',1);
  [i,j] = meshgrid(-order:order,-order:order);
  weights = (abs(i) + abs(j)) <= order;

  psi = get_option(varargin,'kernel');
  if ~isempty(psi)
    weights = psi(sqrt(ebsd.dy^2 * i.^2 + ebsd.dx^2 * j.^2));
  end
else

  order = (min(size(weights)) -1)/2;

end

% set center to zero
weights(order+1,order+1) = 0;

% get threshold
threshold = get_option(varargin,'threshold',inf);

% prepare the result
kam = zeros(size(ebsd));

% for taking the mean count the non nan values
count = zeros(size(kam));

% extract grainIds and make them a bit larger
if isfield(ebsd.prop,'grainId')
  grainId = nan(size(ebsd)+2*order);
  grainId(order+1:end-order,order+1:end-order) = ebsd.grainId;
end

% for all phases
for id = ebsd.indexedPhasesId

  % extract orientations and make them a bit larger
  rot = ebsd.rotations;
  rot(ebsd.phaseId ~= id) = NaN;
  ori = orientation.nan(size(ebsd)+2*order,ebsd.CSList{id});
  ori(order+1:end-order,order+1:end-order) = rot;

  % take the mean
  for i = -order:order
    for j = -order:order

      if  weights(order+1+i,order+1+j) == 0, continue; end

      % compute misorientation angles
      omega = angle(ori(order+1:end-order,order+1:end-order),...
        ori((order+1:end-order)+i,(order+1:end-order)+j));

      % apply angle threshold
      omega(omega > threshold) = NaN;

      % avoid grain boundaries
      if isfield(ebsd.prop,'grainId')
        omega( grainId(order+1:end-order,order+1:end-order) ~= ...
          grainId((order+1:end-order)+i,(order+1:end-order)+j) ) = NaN;
      end

      kam = fun(kam, omega .* weights(order+1+i,order+1+j));
      count = fun(count,~isnan(omega) .* weights(order+1+i,order+1+j));

    end
  end

end

kam = kam ./ count;


text:
%%  kmeans algorithm for an image
%---input---------------------------------------------------------
%   Y: 2D image
%   k: number of clusters
%   g: number of GMM components
%---output--------------------------------------------------------
%   X: 2D labels
%   GMM: Gaussian mixture model parameters

function [X GMM]=image_kmeans(Y,k,g)
[m n temp]=size(Y);
y=reshape(Y,[m*n 3]);
x=kmeans(y,k);
X=reshape(x,[m n]);

GMM=get_GMM(X,Y,g);


text:
function z = nmi(x, y)
% Compute normalized mutual information I(x,y)/sqrt(H(x)*H(y)) of two discrete variables x and y.
% Input:
%   x, y: two integer vector of the same length
% Ouput:
%   z: normalized mutual information z=I(x,y)/sqrt(H(x)*H(y))
% Written by Mo Chen (sth4nth@gmail.com).
assert(numel(x) == numel(y));
n = numel(x);
x = reshape(x,1,n);
y = reshape(y,1,n);

l = min(min(x),min(y));
x = x-l+1;
y = y-l+1;
k = max(max(x),max(y));

idx = 1:n;
Mx = sparse(idx,x,1,n,k,n);
My = sparse(idx,y,1,n,k,n);
Pxy = nonzeros(Mx'*My/n); %joint distribution of x and y
Hxy = -dot(Pxy,log2(Pxy));

% hacking, to elimative the 0log0 issue
Px = nonzeros(mean(Mx,1));
Py = nonzeros(mean(My,1));

% entropy of Py and Px
Hx = -dot(Px,log2(Px));
Hy = -dot(Py,log2(Py));

% mutual information
MI = Hx + Hy - Hxy;

% normalized mutual information
z = sqrt((MI/Hx)*(MI/Hy));
z = max(0,z);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function orientation  = FlyOrient(subset_frame, threshold)

%FLYORIENT
%   Usage:
%       orientation = FlyOrient(subset_frame, threshold)
%
% This function takes in a subset frame around the fly (calculated by
% FindFly) and discards the 3D data by placing points where the pixel intensity
% is larger than the user chosen threshold.  Then, Principal Components
% Analysis (by way fot the pca1 function) is performed on the resulting
% scatter plot to find the direction of maximum variance --- this direction
% is taken to be the fly's (ambiguous) orientation.
% orientation is a vector consisting of two angles (complements) that comprise
% the body axis. The first element is an angle in the upper half plane; the
% second element is an angle in the lower half plane.

% Written by Dan Valente
% 11 October 2006

%Normalize frame data by pixel of maximum intensity
subset_frame = subset_frame/max(max(subset_frame));

% Put dots where fly is and do PCA on reduced data set
[rows, cols] = find(subset_frame >= threshold);
rows = length(subset_frame(:,1))-rows+1;
x = [cols';rows'];
[xnew, PC, V, data] = pca1(x);

% Find orientation vectors (two, mirrored across diagonal), and group into
% upper half and lower half planes.
a1 = PC(1,1);
b1 = PC(2,1);
a2 = -PC(1,1);
b2 = -PC(2,1);
if (b1 >= 0 );
    orientUHP = atan2(b1,a1);
    orientLHP = atan2(b2,a2);
elseif (b2 >=0);
    orientUHP = atan2(b2,a2);
    orientLHP = atan2(b1,a1);
else
end

% The vector we will return
orientation = [orientUHP orientLHP];

return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
clc, clear
a=[-1 -2 0;-1 0 0];b=[-1;0];
[x,y]=ga(@ycfun1,3,a,b,[],[],[],[],@ycfun2);
x, y=-y

author: Eurus-Holmes


text:
function exact = p42_exact ( )

%*****************************************************************************80
%
%% P42_EXACT returns the exact integral for problem 42.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  alpha = p42_param_get ( );

  exact = 2.0^( alpha - 2.0 ) * ( gamma ( alpha / 2.0 ) )^2 / gamma ( alpha );

  return
end


text:
%
% A variational approach to SPCP (Aravkin et al. 2014)
%
% RPCA | SPCP-sum-SPG | Stable PCP-sum solved by Spectral Projected Gradient (Aravkin et al. 2014)
% process_video('RPCA', 'flip-SPCP-sum-SPG', 'dataset/demo.avi', 'output/demo_flip-SPCP-sum-SPG.avi');

alg_path_aux = fullfile(lrs_conf.rpca_path,'SPGL1');
addpath(genpath(alg_path_aux));

nFrames     = size(M,2);
lambda      = 1/sqrt(max(size(M,1),size(M,2)));
L0          = repmat(median(M,2), 1, nFrames);
S0          = M - L0;
epsilon     = 5e-3*norm(M,'fro'); % tolerance for fidelity to data

% Flip-Flop version of SPCP-sum solved by Spectral Projected Gradient
opts = struct('sum',true,'L0',L0,'S0',S0,'max',false,...
  'tau0',3e5,'SPGL1_tol',1e-1,'tol',1e-3);
[L,S] = solver_RPCA_SPGL1(M,lambda,epsilon,[],opts);

rmpath(genpath(alg_path_aux));

author: andrewssobral

repo: lrslibrary


text:
%-----------------------------------------------------------------------
% Job configuration created by cfg_util (rev $Rev: 4899 $)
%-----------------------------------------------------------------------
matlabbatch{1}.menu_cfg.menu_struct.conf_choice.type = 'cfg_choice';
matlabbatch{1}.menu_cfg.menu_struct.conf_choice.name = 'File/Dir Operations';
matlabbatch{1}.menu_cfg.menu_struct.conf_choice.tag = 'file_dir';
matlabbatch{1}.menu_cfg.menu_struct.conf_choice.values = {
                                                          '<UNDEFINED>'
                                                          '<UNDEFINED>'
                                                          '<UNDEFINED>'
                                                          }';
matlabbatch{1}.menu_cfg.menu_struct.conf_choice.check = [];
matlabbatch{1}.menu_cfg.menu_struct.conf_choice.help = {''};

author: spm


text:
function geometry_test0243 ( )

%*****************************************************************************80
%
%% TEST0243 tests R8VEC_ANY_NORMAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 10;
  test_num = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0243\n' );
  fprintf ( 1, '  R8VEC_ANY_NORMAL computes a vector V2 that is normal\n' );
  fprintf ( 1, '  to a given vector V1.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Test    ||V1||      ||V2||        V1.V2\n' );
  fprintf ( 1, '\n' );

  seed = 123456789;

  for test = 1 : test_num

    [ v1, seed ] = r8vec_uniform_01 ( dim_num, seed );
    v1_length = r8vec_norm ( dim_num, v1 );
    v2 = r8vec_any_normal ( dim_num, v1 );
    v2_length = r8vec_norm ( dim_num, v2 );
    v1v2_dot = v1(1:dim_num) * v2(1:dim_num)';
    fprintf ( 1, '  %6d  %10f  %10f  %10f\n', ...
      test, v1_length, v2_length, v1v2_dot );
  end

  return
end


text:
function [fdMc, fProbability, fAICc] = calc_llhd_dMc2(mCat1, mCat2)
% function [fdMc, fProbability, fAICc] = calc_llhd_dMc2(mCat1, mCat2);
% ----------------------------------------------------------------------------------------------
% Calculate log-likelihood estimation of a magnitude of completness change between the two periods
%
% Incoming variable
% mCat1 : EQ catalog period 1 (Catalog to be modified)
% mCat2 : EQ catalog period 2 (Observed catalog)
%
% Outgoing variable
% fdMc         : Change in the magnitude of completeness
% fProbability : log-likelihood probabilty
% fAICc         : Corrected Akaike Information Criterion
%
% Author: J. Woessner, woessner@seismo.ifg,.ethz.ch
% last update: 21.10.03

fBinning = 0.1;
nSample = 100;
nMethod = 6;

% Determine exact time period
fPeriod1 = max(mCat1(:,3)) - min(mCat1(:,3));
fPeriod2 = max(mCat2(:,3)) - min(mCat2(:,3));

% Initialize values
fMinMag = min([min(mCat1(:,6)) min(mCat2(:,6))]);
fMaxMag = max([max(mCat1(:,6)) max(mCat2(:,6))]);

try
%% Calculate model for best fitting Mc, both periods
[mResult1, fMls1, fMc1, fMu1, fSigma1, mDataPredBest1, vPredBest1, fBValue1, fAvalue1] = calc_McCdfnormal(mCat1, fBinning);
[mResult, fMls, fMc2, fMu, fSigma, mDataPredBest, vPredBest, fBValue2, fAvalue2] = calc_McCdfnormal(mCat2, fBinning);

% [fMc1, fStd_Mc1, fBvalue1, fStd_B1, fAvalue1, fStd_A1, vMc1, mBvalue1] = calc_McBboot(mCat1, fBinning, nSample, nMethod)
% [fMc2, fStd_Mc2, fBvalue2, fStd_B2, fAvalue2, fStd_A2, vMc2, mBvalue2] = calc_McBboot(mCat2, fBinning, nSample, nMethod)

% Mc difference
fdMc = fMc2-fMc1;

% Select part of catalog
vSel = mCat1(:,6) >= fMc2-fBinning/2;
mCat1tmp=mCat1(vSel,:);

% Create Model distribution
vMagnitudes = [fMc2:0.1:floor(max(mCat2(:,6)))+0.1];
% Productuvity for fMc2
nNumberEvents = length(mCat1tmp(:,1));
vNumbers = 10.^(log10(nNumberEvents) - fBValue1*(vMagnitudes-fMc2));
%vNumbers = 10.^(fAvalue1- fBValue1*(vMagnitudes-fMc2));
vNumbers = round(vNumbers);
vNCumFMD = round(-diff(vNumbers));
% Calculate synthetic data below Mc
fMinMag = min(mCat2(:,6));
vMagstep = fMinMag:0.1:fMc2-0.1;
vProb = normcdf(vMagstep,fMu1, fSigma1);
vProb = vProb';
vMagstep = vMagstep';

% Calculate number of EQs in bins
fN_Mc = vNCumFMD(1,1);
vN = round(vProb(:,1)*fN_Mc);
%mNonCumModel = [vN vMagstep; vNCumFMD' vBin'];
mNonCumModel = [vN vMagstep; vNCumFMD' vMagnitudes(:,1:end-1)'];

% vMags = mDataPredBest(:,2)';
% FMD to be modeled (second period)
[vObsFMD,vBin2] = hist(mCat2(:,6),roundn(min(mCat2(:,6)),-1):0.1:floor(max(mCat2(:,6))));

% Select bins to calculate loglikelihood
vSel = (mNonCumModel(:,2) < min(vBin2) | mNonCumModel(:,2) > max(vBin2));
vPredFMD = mNonCumModel(~vSel,:);
% Normalize
vObsFMD = ceil(vObsFMD./fPeriod2);
vPredFMD(:,1) = vPredFMD(:,1)./fPeriod1;

if length(vObsFMD') ~= length(vPredFMD(:,1))
    disp('warning')
end
% Calculate the likelihood
vProb_ = calc_log10poisspdf2(vObsFMD', vPredFMD(:,1));
% Sum the probabilities
fProbability = (-1) * sum(vProb_);

nDegFree = 1; % degree of freedom
n_samples = length(mCat2(:,6));
%% Corrected Akaike Information Criterion (AICc)
fAICc = -2*(-fProbability)+2*nDegFree+2*nDegFree*(nDegFree+1)/(n_samples-nDegFree-1);

catch
    fdMc = nan;
    fProbability = nan;
    fAICc = nan;
end
% figure_w_normalized_uicontrolunits(30)
% %plot(mNonCumModel(:,2),mNonCumModel(:,1),'r*',vBin2,vObsFMD,'o')
% histogram(mCat1(:,6),fMinMag:0.1:fMaxMag);
% hold on
% histogram(mCat2(:,6),fMinMag:0.1:fMaxMag)
% h = findobj(gca,'Type','patch');
% set(h,'FaceColor','r','EdgeColor','w')
% plot(vPredFMD(:,2),vPredFMD(:,1),'g*',vBin2,vObsFMD,'o')

author: CelsoReyes


text:
function [prob,p] = coinRead(filename,type,print)
%coinRead  Read a Mathematical File and convert to Matlab Values
%
% prob = coinRead(filename,type) reads the file specified by filename
% and problem type specified by type and converts it to an optiprob. If you
% do not specify a file extension it will default to the type specified.
% The returned structure is solver independent so the user can manually
% extract matrices as required or supply it directly to opti().
%
% prob = coinRead(filename,type,print) specifies to display information
% and error messages in MATLAB when print > 0. The default is 0.
%
%   Available File Types:
%       - MPS  [.mps]
%       - QPS  [.qps]
%       - LP   [.lp]  (Appears to be CPLEX LP version - Simplified)
%       - GMPL [.mod] (Subset of AMPL, implemented via GLPK)
%       - GAMS [.gms] (Only first COIN version... Will probably not work)
%
% You may specify a full path to the file, or if you specify a filename
% only, it must exist on the MATLAB path.
%
% The routines underneath use COIN-OR & GLPK utilities for File IO. See
% attached EPL License for COIN-OR and GPL for GLPK.

%   Copyright (C) 2011 Jonathan Currie (I2C2)

%Quick Checks
if(~exist('type','var') || isempty(type))
    %See if contained within file name
    dots = strfind(filename,'.');
    if(isempty(dots))
        error('You must supply the file type to this function');
    else %hopefully got a good ext, will check below
        type = filename(dots(end)+1:end);
    end
end
if(~exist('print','var') || isempty(print))
    print = 0;
end
if(~ischar(filename))
    error('Filename must be a char array!');
end

%Determine file type
switch(lower(type))
    case 'mps'
        filetype = 'mps';
    case 'qps'
        filetype = 'qps';
    case 'lp'
        filetype = 'lp';
    case {'gmpl','mod'}
        filetype = 'mod';
    case {'gams','gms'}
        filetype = 'gms';
    otherwise
        error('Unknown file type %s',type);
end

%Check if extension is specified
if(isempty(strfind(filename,'.')))
    filename = [filename '.' filetype];
end

% If filename includes absolute path, we can skip which
if(~isempty(strfind(filename,':')))
    p = filename;
    %Check it exists
    if(~exist(filename,'file'))
        error('Cannot locate file %s!',filename);
    end
else %Locate the full path to the file (must be on matlab path)
    p = which(filename);
    %Check it exists
    if(isempty(p))
        error('Cannot locate file %s!',filename);
    end
end

%Call coin-or mex
[f,A,rl,ru,lb,ub,ivars,H,name,sostype,sosind,soswt,objc] = coinR(p,filetype,print);

%No guarantee matrices have been created with correct order (T.Kelman, 6/9/12), double transpose to ensure order)
H = H'';
A = A'';

%Remove H if all zero
if(~nnz(H))
    H = [];
%Assume we need to make symmetric if e.g. [8 0; 2 10]
elseif(nnz(triu(H,1)) == 0 && ~(nnz(tril(H,-1)) == 0))
    optiwarn('OPTI:COINR_QUADOBJ','Assumed QOBJ has been specified as lower triangular only, forcing to symmetric');
    H = H + tril(H,-1)';
%Check for upper tri only
elseif(nnz(tril(H,-1)) == 0 && ~(nnz(triu(H,1)) == 0))
    optiwarn('OPTI:COINR_QUADOBJ','Assumed QOBJ has been specified as upper triangular only, forcing to symmetric');
    H = H + triu(H,1)';
end

%Build Integer String
int = repmat('C',size(f')); ivars = logical(ivars);
int(ivars) = 'I'; %assume just integer variables (no binary - even with MPS)
%Build SOS
if(~isempty(sostype))
    sostype = num2str(sostype);
    if(length(sostype) == 1)
        sosind = sosind{1}; %remove cells
        soswt = soswt{1};
    end
else
    sostype = []; sosind = []; soswt = [];
end

%Build return problem
prob = optiprob('H',H,'f',f,'objbias',objc,'lin',A,rl,ru,'bounds',lb,ub,'int',int,'name',name,'sos',sostype,sosind,soswt);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [ w, x ] = gm_rule_set ( rule, dim_num, point_num )

%*****************************************************************************80
%
%% GM_RULE_SET sets a Grundmann-Moeller rule.
%
%  Discussion:
%
%    This is a revised version of the calculation which seeks to compute
%    the value of the weight in a cautious way that avoids intermediate
%    overflow.  Thanks to John Peterson for pointing out the problem on
%    26 June 2008.
%
%    This rule returns weights and abscissas of a Grundmann-Moeller
%    quadrature rule for the DIM_NUM-dimensional unit simplex.
%
%    The dimension POINT_NUM can be determined by calling GM_RULE_SIZE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Axel Grundmann, Michael Moeller,
%    Invariant Integration Formulas for the N-Simplex
%    by Combinatorial Methods,
%    SIAM Journal on Numerical Analysis,
%    Volume 15, Number 2, April 1978, pages 282-290.
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule.
%    0 <= RULE.
%
%    Input, integer DIM_NUM, the spatial dimension.
%    1 <= DIM_NUM.
%
%    Input, integer POINT_NUM, the number of points in the rule.
%
%    Output, real W(POINT_NUM), the weights.
%
%    Output, real X(DIM_NUM,POINT_NUM), the abscissas.
%
  s = rule;
  d = 2 * s + 1;
  k = 0;
  n = dim_num;
  one_pm = 1;

  for i = 0 : s

    weight =  one_pm;

    for j = 1 : max ( n, d, d + n - i )

      if ( j <= n )
        weight = weight *  ( j );
      end
      if ( j <= d )
        weight = weight * ( d + n - 2 * i );
      end
      if ( j <= 2 * s )
        weight = weight / 2.0;
      end
      if ( j <= i )
        weight = weight / j;
      end
      if ( j <= d + n - i )
        weight = weight / j;
      end

    end

    one_pm = - one_pm;

    beta_sum = s - i;
    more = 0;
    beta = [];
    h = 0;
    t = 0;

    while ( 1 )

      [ beta, more, h, t ] = comp_next ( beta_sum, dim_num + 1, ...
        beta, more, h, t );

      k = k + 1;

      w(k) = weight;

      x(1:dim_num,k) = ( 2 * beta(2:dim_num+1)' + 1 ) / ( d + n - 2 * i );

      if ( ~more )
        break
      end

    end

  end

  return
end


text:
function y = cvx_isconcave( x, full )
error( nargchk( 1, 2, nargin ) );
if nargin == 2,
    y = ~imag( x );
elseif isreal( x ),
    y = true;
else
    y = nnz(imag(x)) == 0;
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function varargout=rem(varargin)
%REM (overloaded)

switch class(varargin{1})

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.

        if ~isa(varargin{2},'double')
            error('REM is currently only supported for DOUBLE second argument');
        end

        x = varargin{1};
        y = varargin{2};
        % Some boring code for scalarization
        if prod(size(x))==1 & prod(size(y))>1
            x = x*ones(size(y));
        elseif prod(size(y))==1 & prod(size(x))>1
            y = y*ones(size(x));
        end
        if ~all(size(x) == size(y))
            error('Matrix dimensions must agree.');
        end
        dim = size(x);
        x = reshape(x,prod(dim),1);
        y = reshape(y,prod(dim),1);
        z = [];
        % Create one variable for each element
        for i = 1:length(x)
            xi = extsubsref(x,i);
            yi = extsubsref(y,i);
            inarg = {xi,yi};
            z = [z;yalmip('define',mfilename,inarg{:})];
        end
        z = reshape(z,dim);
        varargout{1} = z;

    case 'char' % YALMIP send 'graph' when it wants the epigraph or hypograph

        % Description using epigraphs
        t = varargin{2};
        x = varargin{3};
        y = varargin{4};

        % t = rem(x,y), i.e. t = x - n*y, n = fix(x/y)
        d1 = binvar(1,1);
        d2 = binvar(1,1);
        [M,m] = derivebounds(x);
        n = intvar(1,1);
        F = [t == x - y*n, d1+d2 == 1];
        F = [F,x>=m*(1-d1), x<=M*(1-d2),(x/y)-d1 <= n <= (x/y)+d2];
        if ~isinf(m) && isa(y,'double')
            F = [F, fix(m/y) <= n];
        end
        if ~isinf(M) && isa(y,'double')
            F = [F, n <= fix(M/y)];
        end

        varargout{1} = F;
        varargout{2} = struct('convexity','none','monotonicity','none','definiteness','none','model','integer');
        varargout{3} = [x(:);y(:)];

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%% (Internal) Replace NaN from PRdatasets
%
%   dsOut = deNaN_dataset(dsIn, type)
%
% Arguments:
%
%      + dsIn: paths to add
%
%      + type: Choose the type of replacement:
%
%             - 'remove', remove the whole vector which includes NaN
%
%             - 'change', change NaNs for repeated values of the same
%             feature, added an small variance, estimated ignoring this nan
%             values.
%
%             - 'change_same_class' change NaNs for repeated values of the
%             same feature, same class, added an small variance, estimated
%             ignoring this nan values.
%
%      + func_ptr: A pointer to an isX function. Default @isnan
%
% Output:
%
%      + dsOut: the clean PRdataset
%
% Example:
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
%
function dsOut = deNaN_dataset(dsIn, type, func_ptr)

    if( nargin < 2 || isempty(type) )
        type = 'remove';
    end

    if( nargin < 3 || isempty(func_ptr) )
        func_ptr = @isnan;
    end

    features_with_nans = find(any(func_ptr(+dsIn)));
    rows_without_nans = ~any(func_ptr(+dsIn),2);

    if( any(~rows_without_nans) )

        rows_without_nans_idx = find(rows_without_nans);

        switch(type)

            case 'remove'
                % just remove nans

                dsOut = dsIn(rows_without_nans_idx,:);

            case 'change'

                % change NaNs for repeated values of the same feature, added an small variance, estimated ignoring this nan
                % values.
                dsOut = dsIn;

                robust_std = nanmeda(+dsIn(rows_without_nans_idx,:));

                for kk = rowvec(features_with_nans)

                    nan_idx = find(func_ptr(+(dsIn(:,kk))));
                    not_nan_idx = find(~func_ptr(+(dsIn(:,kk))));
                    laux_idx = length(nan_idx);

                    dsIn(nan_idx,kk) = randsample( +dsIn(not_nan_idx,kk), laux_idx, true ) + 1/100*robust_std(kk)*rand(laux_idx,1);

                    dsOut(:,kk) = dsIn(:,kk);

                end

            case 'change_same_class'

                % change NaNs for repeated values of the same feature, same
                % class, added an small variance, estimated ignoring this nan
                % values.
                dsOut = dsIn;

                dsIn = seldat(dsIn);

                ds_aux = dsIn(rows_without_nans_idx,:);
                ds_aux = setprior(ds_aux, 0);

                w_aux = qdc(ds_aux);
                w_aux = +w_aux;
                Class_indexes = getnlab(dsIn);

                for kk = rowvec(features_with_nans)
                    nan_idx = find(func_ptr(+(dsIn(:,kk))));

                    classes_involved = unique(Class_indexes(nan_idx));

                    for ll = rowvec(classes_involved)

                        aux_idx = find(Class_indexes(nan_idx) == ll);
                        this_class_not_nan_idx = find(Class_indexes == ll);
                        this_class_not_nan_idx = setdiff(this_class_not_nan_idx, nan_idx(aux_idx) );

                        laux_idx = length(aux_idx);

                        dsIn(nan_idx(aux_idx),kk) = randsample( +dsIn(this_class_not_nan_idx,kk), laux_idx, true ) + 1/100*w_aux.cov(kk,kk,ll)*rand(laux_idx,1);

                    end

                    dsOut(:,kk) = dsIn(:,kk);

                end

        end
    else

        dsOut = dsIn;

    end

author: marianux

repo: ecg-kit

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


text:
function [ point_num, edge_num, face_num, face_order_max ] = ...
  soccer_size_3d ( )

%*****************************************************************************80
%
%% SOCCER_SIZE_3D gives "sizes" for a truncated icosahedron in 3D.
%
%  Discussion:
%
%    The shape is a truncated icosahedron, which is the design used
%    on a soccer ball.  There are 12 pentagons and 20 hexagons.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    http://polyhedra.wolfram.com/uniform/u25.html
%
%  Parameters:
%
%    Output, integer POINT_NUM, the number of points.
%
%    Output, integer EDGE_NUM, the number of edges.
%
%    Output, integer FACE_NUM, the number of faces.
%
%    Output, integer FACE_ORDER_MAX, the maximum order of any face.
%
  point_num = 60;
  edge_num = 90;
  face_num = 32;
  face_order_max = 6;

  return
end


text:
function h = drawPlane3d(plane, varargin)
%DRAWPLANE3D Draw a plane clipped by the current axes.
%
%   drawPlane3d(PLANE) draws a plane of the format:
%       [x0 y0 z0  dx1 dy1 dz1  dx2 dy2 dz2]
%
%   drawPlane3d(...,'PropertyName',PropertyValue,...) sets the value of the
%   specified patch property. Multiple property values can be set with
%   a single statement. See function patch for details.
%
%   drawPlane3d(AX,...) plots into AX instead of GCA.
%
%   H = drawPlane3d(...) returns a handle H to the patch object.
%
%   Example
%
%     p0 = [1 2 3];
%     v1 = [1 0 1];
%     v2 = [0 -1 1];
%     plane = [p0 v1 v2];
%     axis([-10 10 -10 10 -10 10]);
%     drawPlane3d(plane)
%     drawLine3d([p0 v1])
%     drawLine3d([p0 v2])
%     set(gcf, 'renderer', 'zbuffer');
%
%   See also
%   planes3d, createPlane, patch

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% INRA - TPV URPOI - BIA IMASTE
% created the 17/02/2005.
%
%   HISTORY
%   2008-10-30 replace intersectPlaneLine by intersectLinePlane, add doc
%   2010-10-04 fix a bug for planes touching box by one corner
%   2011-07-19 fix a bug for param by Xin KANG (Ben)
%

% Parse and check inputs
valFun = @(x) size(x,1)==1 && isPlane(x);
defOpts.FaceColor = 'm';
[hAx, plane, varargin] = ...
    parseDrawInput(plane, valFun, 'patch', defOpts, varargin{:});

% extract axis bounds to crop plane
lim = get(hAx, 'xlim');
xmin = lim(1);
xmax = lim(2);
lim = get(hAx, 'ylim');
ymin = lim(1);
ymax = lim(2);
lim = get(hAx, 'zlim');
zmin = lim(1);
zmax = lim(2);

% create lines corresponding to cube edges
lineX00 = [xmin ymin zmin 1 0 0];
lineX01 = [xmin ymin zmax 1 0 0];
lineX10 = [xmin ymax zmin 1 0 0];
lineX11 = [xmin ymax zmax 1 0 0];

lineY00 = [xmin ymin zmin 0 1 0];
lineY01 = [xmin ymin zmax 0 1 0];
lineY10 = [xmax ymin zmin 0 1 0];
lineY11 = [xmax ymin zmax 0 1 0];

lineZ00 = [xmin ymin zmin 0 0 1];
lineZ01 = [xmin ymax zmin 0 0 1];
lineZ10 = [xmax ymin zmin 0 0 1];
lineZ11 = [xmax ymax zmin 0 0 1];

% compute intersection points with each plane
piX00 = intersectLinePlane(lineX00, plane);
piX01 = intersectLinePlane(lineX01, plane);
piX10 = intersectLinePlane(lineX10, plane);
piX11 = intersectLinePlane(lineX11, plane);
piY00 = intersectLinePlane(lineY00, plane);
piY01 = intersectLinePlane(lineY01, plane);
piY10 = intersectLinePlane(lineY10, plane);
piY11 = intersectLinePlane(lineY11, plane);
piZ00 = intersectLinePlane(lineZ00, plane);
piZ01 = intersectLinePlane(lineZ01, plane);
piZ10 = intersectLinePlane(lineZ10, plane);
piZ11 = intersectLinePlane(lineZ11, plane);

% concatenate points into one array
points = [...
    piX00;piX01;piX10;piX11; ...
    piY00;piY01;piY10;piY11; ...
    piZ00;piZ01;piZ10;piZ11;];

% check validity: keep only points inside window (with tolerance)
ac = sqrt (eps);
ivx = points(:,1) >= xmin-ac & points(:,1) <= xmax+ac;
ivy = points(:,2) >= ymin-ac & points(:,2) <= ymax+ac;
ivz = points(:,3) >= zmin-ac & points(:,3) <= zmax+ac;
valid = ivx & ivy & ivz;
pts = unique(points(valid, :), 'rows');

% If there is no intersection point, escape.
if size(pts, 1) < 3
    disp('plane is outside the drawing window');
    if nargout > 0
        h = [];
    end
    return;
end

% the two spanning lines of the plane
d1 = plane(:, [1:3 4:6]);
d2 = plane(:, [1:3 7:9]);

% position of intersection points in plane coordinates
u1 = linePosition3d(pts, d1);
u2 = linePosition3d(pts, d2);

% reorder vertices in the correct order
ind = convhull(u1, u2);
ind = ind(1:end-1);

% draw the patch
htmp = patch( ...
    'XData', pts(ind,1), ...
    'YData', pts(ind,2), ...
    'ZData', pts(ind,3), ...
    'Parent', hAx, varargin{:});

% Do not return axis if not requested
% avoids output when called without semicolon
if nargout > 0
    h = htmp;
end


text:
function sol = SolveBP(A, y, N, maxIters, lambda, OptTol)
% SolveBP: Solves a Basis Pursuit problem
% Usage
%	sol = SolveBP(A, y, N, maxIters, lambda, OptTol)
% Input
%	A           Either an explicit nxN matrix, with rank(A) = min(N,n)
%               by assumption, or a string containing the name of a
%               function implementing an implicit matrix (see below for
%               details on the format of the function).
%	y           vector of length n.
%   N           length of solution vector.
%	maxIters    maximum number of PDCO iterations to perform, default 20.
%   lambda      If 0 or omitted, Basis Pursuit is applied to the data,
%               otherwise, Basis Pursuit Denoising is applied with
%               parameter lambda (default 0).
%	OptTol      Error tolerance, default 1e-3
% Outputs
%	 sol        solution of BP
% Description
%   SolveBP solves the basis pursuit problem
%      min ||x||_1 s.t. A*x = y
%   by reducing it to a linear program, and calling PDCO, a primal-dual
%   log-barrier algorithm. Alternatively, if lambda ~= 0, it solves the
%   Basis Pursuit Denoising (BPDN) problem
%      min lambda*||x||_1 + 1/2||y - A*x||_2^2
%   by transforming it to an SOCP, and calling PDCO.
%   The matrix A can be either an explicit matrix, or an implicit operator
%   implemented as a function. If using the implicit form, the user should
%   provide the name of a function of the following format:
%     y = OperatorName(mode, m, n, x, I, dim)
%   This function gets as input a vector x and an index set I, and returns
%   y = A(:,I)*x if mode = 1, or y = A(:,I)'*x if mode = 2.
%   A is the m by dim implicit matrix implemented by the function. I is a
%   subset of the columns of A, i.e. a subset of 1:dim of length n. x is a
%   vector of length n is mode = 1, or a vector of length m is mode = 2.
% See Also
%   SolveLasso, SolveOMP, SolveITSP
%

if nargin < 6,
	OptTol = 1e-3;
end
if nargin < 5,
	lambda = 0;
end
if nargin < 4,
    maxIters = 20;
end

n = length(y);

n_pdco = 2*N;    % Input size
m_pdco = n;      % Output size

% upper and lower bounds
bl = zeros(n_pdco,1);
bu = Inf .* ones(n_pdco,1);

% generate the vector c
if (lambda ~= 0)
    c = lambda .* ones(n_pdco,1);
else
    c = ones(n_pdco,1);
end

% Generate an initial guess
x0 = ones(n_pdco,1)/n_pdco;       % Initial x
y0 = zeros(m_pdco,1);             % Initial y
z0 = ones(n_pdco,1)/n_pdco;       % Initial z

d1 = 1e-4;                 % Regularization parameters
if (lambda ~= 0) % BPDN
    d2 = 1;
else
    d2 = 1e-4;
end

xsize = 1;                 % Estimate of norm(x,inf) at solution
zsize = 1;                 % Estimate of norm(z,inf) at solution

options = pdcoSet;         % Option set for the function pdco
options = pdcoSet( options, ...
                     'MaxIter    ', maxIters  , ...
                     'FeaTol     ', OptTol    , ...
                     'OptTol     ', OptTol    , ...
                     'StepTol    ', 0.99      , ...
                     'StepSame   ', 0         , ...
                     'x0min      ', 0.1       , ...
                     'z0min      ', 1.0       , ...
                     'mu0        ', 0.01      , ...
                     'method     ', 1         , ...
                     'LSQRMaxIter', 20        , ...
                     'LSQRatol1  ', 1e-3      , ...
                     'LSQRatol2  ', 1e-15     , ...
                     'wait       ', 0    );

if (ischar(A) || isa(A, 'function_handle'))
    [xx,yy,zz,inform,PDitns,CGitns,time] = ...
        pdco(c, @pdcoMat, y, bl, bu, d1, d2, options, x0, y0, z0, xsize, zsize);
else
    Phi = [A -A];
    [xx,yy,zz,inform,PDitns,CGitns,time] = ...
        pdco(c, Phi, y, bl, bu, d1, d2, options, x0, y0, z0, xsize, zsize);
end

% Extract the solution from the output vector x
sol = xx(1:N) - xx((N+1):(2*N));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function y = pdcoMat(mode,m,n,x)
        if (mode == 1) % Direct operator
            % Decompose input
            n2 = n/2;
            u = x(1:n2);
            v = x(n2+1:n);

            % Apply matrix A
            Au = feval(A,1,m,n2,u,1:n2,n2);
            Av = feval(A,1,m,n2,v,1:n2,n2);

            y = Au-Av;
        else % Adjoint operator
            n2 = n/2;
            Atx = feval(A,2,m,n2,x,1:n2,n2);
            y = [Atx; -Atx];
        end
    end

end

%
% Copyright (c) 2006. Yaakov Tsaig
%

%
% Part of SparseLab Version:100
% Created Tuesday March 28, 2006
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail sparselab@stanford.edu
%

author: gpeyre

repo: matlab-toolboxes


text:
function obj = initializeObj(obj)
	% DESCRIPTION
	% Biafra Ahanonu
	% started: 2014.07.31 - branch from calciumImagingAnalysis 2020.05.07 [15:47:29]
	% inputs
		%
	% outputs
		%

	% changelog
		% 2020.06.05 [11:29:24] - Improved toolbox check.
		% 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
		% 2021.11.16 [11:13:07] - Wrap calls to diverging_map to avoid potential ciatah initialization errors on macOS.
		% 2022.06.27 [12:48:22] - Ensure all save paths (pics, data, log, settings, video, etc.) use absolute CIAtah private directory to match CIAtah default properties initialization and to avoid instances where users have non-CIAtah root path working directory.
	% TODO
		%

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

	% load dependencies.
	% loadBatchFxns();
	cnmfVersionDirLoad('none','displayOutput',0);
	% [success] = cnmfVersionDirLoad('cnmfe');

	% try
	% 	% Attempt to load the VERSION from file
	% 	verN = importdata('ciapkg\VERSION');
	% 	obj.classVersion = verN{1};
	% catch

	% end

	% Load colormaps
	obj.colormap = customColormap({[0 0 1],[1 1 1],[0.5 0 0],[1 0 0]});
	obj.colormapAlt = customColormap({[0 0 0.7],[1 1 1],[0.7 0 0]});
	% obj.colormapAlt2 = diverging_map(linspace(0,1,100),[0 0 0.7],[0.7 0 0]);

	% try-catch to handle cases in which diverging_map is not in the path.
	try
		obj.colormapAlt2 = diverging_map(linspace(0,1,100),[0 0 0.7],[0.7 0 0]);
	catch
		obj.colormapAlt2 = customColormap({[0 0.1 0.7],[1 1 1],[0.7 0.1 0]});
	end

	% Check required toolboxes are available, warn if not
	disp(repmat('*',1,42))
	ciapkg.io.dependencyCheck();
	disp(repmat('*',1,42))

	% Ensure date paths are up to date
	obj.picsSavePath = [ciapkg.getDir() filesep 'private' filesep 'pics' filesep datestr(now,'yyyymmdd','local') filesep];
	obj.dataSavePath = [ciapkg.getDir() filesep 'private' filesep 'data' filesep datestr(now,'yyyymmdd','local') filesep];
	obj.logSavePath = [ciapkg.getDir() filesep 'private' filesep 'logs' filesep datestr(now,'yyyymmdd','local') filesep];
	obj.settingsSavePath = [ciapkg.getDir() filesep 'private' filesep 'settings'];
	obj.videoSaveDir = [ciapkg.getDir() filesep 'private' filesep 'vids' filesep datestr(now,'yyyymmdd','local') filesep];

	% ensure private folders are set
	if ~exist(obj.picsSavePath,'dir');mkdir(obj.picsSavePath);fprintf('Creating directory: %s\n',obj.picsSavePath);end
	if ~exist(obj.dataSavePath,'dir');mkdir(obj.dataSavePath);fprintf('Creating directory: %s\n',obj.dataSavePath);end
	if ~exist(obj.logSavePath,'dir');mkdir(obj.logSavePath);fprintf('Creating directory: %s\n',obj.logSavePath);end
	if ~exist(obj.settingsSavePath,'dir');mkdir(obj.settingsSavePath);fprintf('Creating directory: %s\n',obj.settingsSavePath);end
	if ~exist(obj.videoSaveDir,'dir');mkdir(obj.videoSaveDir);fprintf('Creating directory: %s\n',obj.videoSaveDir);end

	% load user specific settings
	loadUserSettings = [obj.settingsSavePath filesep 'calciumImagingAnalysisInitialize.m'];
	if exist(loadUserSettings,'file')~=0
		run(loadUserSettings);
	else
		% create privateLoadBatchFxns.m
	end

	% if use puts in a single folder or a path to a txt file with folders
	if ~isempty(obj.rawSignals)&ischar(obj.rawSignals)
		if isempty(regexp(obj.rawSignals,'.txt'))&exist(obj.rawSignals,'dir')==7
			% user just inputs a single directory
			obj.rawSignals = {obj.rawSignals};
		else
			% user input a file linking to directories
			fid = fopen(obj.rawSignals, 'r');
			tmpData = textscan(fid,'%s','Delimiter','\n');
			obj.rawSignals = tmpData{1,1};
			fclose(fid);
		end
		obj.inputFolders = obj.rawSignals;
		obj.dataPath = obj.rawSignals;
	end
	% add subject information to object given datapath
	if ~isempty(obj.dataPath)
		obj.modelGetFileInfo();
	else
		disp('No folder paths input, run modelAddNewFolders method to add new folders.');
		% warning('Input data paths for all files!!! option: dataPath')
	end
	if ~isempty(obj.discreteStimulusTable)&~strcmp(class(obj.discreteStimulusTable),'table')
		obj.modelReadTable('table','discreteStimulusTable');
		obj.modelTableToStimArray('table','discreteStimulusTable','tableArray','discreteStimulusArray','nameArray','stimulusNameArray','idArray','stimulusIdArray','valueName',obj.stimulusTableValueName,'frameName',obj.stimulusTableFrameName);
	end
	if ~isempty(obj.continuousStimulusTable)&~strcmp(class(obj.continuousStimulusTable),'table')
		obj.delimiter = ',';
		obj.modelReadTable('table','continuousStimulusTable','addFileInfoToTable',1);
		obj.delimiter = ',';
		obj.modelTableToStimArray('table','continuousStimulusTable','tableArray','continuousStimulusArray','nameArray','continuousStimulusNameArray','idArray','continuousStimulusIdArray','valueName',obj.stimulusTableValueName,'frameName',obj.stimulusTableFrameName,'grabStimulusColumnFromTable',1);
	end
	% load behavior tables
	if ~isempty(obj.behaviorMetricTable)&~strcmp(class(obj.behaviorMetricTable),'table')
		obj.modelReadTable('table','behaviorMetricTable');
		obj.modelTableToStimArray('table','behaviorMetricTable','tableArray','behaviorMetricArray','nameArray','behaviorMetricNameArray','idArray','behaviorMetricIdArray','valueName','value');
	end
	% modify stimulus naming scheme
	if ~isempty(obj.stimulusNameArray)
		obj.stimulusSaveNameArray = obj.stimulusNameArray;
		obj.stimulusNameArray = strrep(obj.stimulusNameArray,'_',' ');
	end
	% load all the data
	if ~isempty(obj.rawSignals)&ischar(obj.rawSignals{1})
		disp('paths input, going to load files')
		obj.guiEnabled = 0;
		obj = modelVarsFromFiles(obj);
		obj.guiEnabled = 1;
	end
	% check if signal peaks have already been calculated
	if isempty(obj.signalPeaks)&~isempty(obj.rawSignals)
		% obj.computeSignalPeaksFxn();
	else
		disp('No folder data specified, load data with modelVarsFromFiles method.');
		% warning('no signal data input!!!')
	end
	% load stimulus tables

	disp(repmat('*',1,42))
	% Check java heap size
	try
		javaHeapSpaceSizeGb = java.lang.Runtime.getRuntime.maxMemory*1e-9;
		javaHeapPrefGb = 2;
		javaHelpUrl = "https://www.mathworks.com/help/matlab/matlab_external/java-heap-memory-preferences.html";
		javaTmpStr = sprintf('\n\nGo to Preferences->General->"Java Heap Memory" and increase to >%d MB. See more on how to change Java heap memory online at:Java Heap Memory Preferences (%s).',javaHeapPrefGb*1024,javaHelpUrl,javaHelpUrl,javaHelpUrl);
		if javaHeapSpaceSizeGb<javaHeapPrefGb
			javaErrorStr = @(x) sprintf('Java max heap memory is %0.3f GB. This might cause Miji errors when loading videos due to insufficient memory.%s\n\nOtherwise put "java.opts" (located in "ciapkg" sub-folder) file in the MATLAB start-up path or change MATALB start-up folder to the calciumImagingAnalysis root folder then restart MATLB before continuing.\n',javaHeapSpaceSizeGb,x);
			warning(javaErrorStr(javaTmpStr));
			ciapkg.overloaded.msgbox(javaErrorStr(sprintf('\n\nSee command line for URL to change Java heap memory.')),'Note to user','modal')
		else
			fprintf('Java max heap memory is %0.3f Gb, this should be sufficient to run Miji without problems. Please change "java.opts" file to increase heap space if run into Miji memory errors.\n',javaHeapSpaceSizeGb);
		end
	catch err
		disp(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		disp(repmat('@',1,7))
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function out = cfg_run_named_dir(job)

% Return selected dirs as separate output arguments.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_run_named_dir.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

out.dirs = job.dirs;


text:
function [img,seg] = pre_process_isotropic(img_nii,seg_nii,use_isotropic,id)
% Pre-process the data
% Normalize the intensity of the images and also do some re-sample.
%
% If use_isotropic is True, we will resize all samples to same resolution,
% otherwise, we will only resize the sample 4 and 5 (which have much low
% resolution than other samples)
%------------------------------------------------------------------------

    if nargin <4
        id = 1;
    end

    if use_isotropic
        scale = img_nii.hdr.dime.pixdim(2:4)./[0.80 0.7588 0.7588];
        new_size = round(img_nii.hdr.dime.dim(2:4).*scale);
        img = imresize3d(img_nii.img,[],new_size,'cubic','bound');
        if ~isempty(seg_nii)
            seg = imresize3d(seg_nii.img,[],new_size,'nearest','bound');
        end
    else
        if id==4 || id==5
            % resize img in up and down plane
            img = permute(img_nii.img,[3,2,1]);
            img = imresize(img,1.5,'bilinear');
            img = permute(img,[3,2,1]);
            if ~isempty(seg_nii)
                seg = permute(seg_nii.img,[3,2,1]);
                seg = imresize(seg,1.5,'nearest');
                seg = permute(seg,[3,2,1]);
            end
        else
            img = img_nii.img;
            if ~isempty(seg_nii)
                seg = seg_nii.img;
            end
        end
    end

    %Normalize the intensity of images
    img = single(img);
    mask = img>0;
    mean_value = mean(img(mask));
    std_value  = std(img(mask),1);
    img = bsxfun(@rdivide, bsxfun(@minus, img, mean_value), std_value);
end

author: yulequan

repo: HeartSeg

save_path: github-repos/MATLAB/yulequan-HeartSeg


text:
function [y] = crossprodw(x, w)

% expand, compute cross-product and whiten
y = w*crossprod(x);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function transformParams = buildAntsTransform(antsWarpProducts,inverseFlag)

if ~exist('inverseFlag','var')
    inverseFlag = 0;
end

% generate transform command
transformParams ='';

if ~inverseFlag
    if ~isempty(antsWarpProducts.Warp)
        if exist(antsWarpProducts.Warp,'file')
            transformParams = [' -t ' antsWarpProducts.Warp ' '];
        else
            error(['Unable to complete transform, warp product missing: ' antsWarpProducts.Warp]);
        end
    end
    if ~isempty(antsWarpProducts.Affine)
        if exist(antsWarpProducts.Affine, 'file')
            transformParams = [transformParams ' -t ' antsWarpProducts.Affine ' '];
        else
            error(['Unable to complete transform, affine product missing: ' antsWarpProducts.Affine]);
        end
    end
    if isempty(transformParams)
        error('ANTs: No transformation products specified');
    end
else
    if ~isempty(antsWarpProducts.InverseWarp)
        if exist(antsWarpProducts.InverseWarp,'file')
            transformParams = [' -t ' antsWarpProducts.InverseWarp ' '];
        else
            error(['Unable to complete transform, warp product missing: ' antsWarpProducts.Warp]);
        end
    end
    if ~isempty(deformS.algorithmParamsS.antsWarpProducts.Affine)
        if exist(deformS.algorithmParamsS.antsWarpProducts.Affine, 'file')
            transformParams = [transformParams ' -t  [ ' antsWarpProducts.Affine ',1] '];
        else
            error(['Unable to complete transform, affine product missing: ' antsWarpProducts.Affine]);
        end
    end
    if isempty(transformParams)
        error('ANTs: No transformation products specified');
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR