text:
function [atlasView, atlasTypeNum] = getAtlasView()
%
% [atlasView, atlasTypeNum] = getAtlasView()
%
% Author:  Brewer, Wandell
% Purpose:
%
% We check which of the currently open FLAT windows have dataTYPE set to
% Atlases-???'. If there is only one, then that FLAT view is returned. If
% there are none, an empty array is returned. If there is more than one,
% the user is prompted to select the figure showing the desired Atlas.
%
mrGlobals;
atlasViewList = {};

% Changed -- if createAtlas fails this could be the cause.
atlasTypeNum = existDataType('Atlases',[],0);

for(ii=1:length(FLAT))
    % strmatch will match if only the first part matches
    if(~isempty(FLAT{ii}) & strmatch('Atlases',getDataTypeName(FLAT{ii})))
        % These are the FLAT{} data structures with an Atlases type
        % selected.
        atlasViewList{end+1} = ii;
    end
end

if(isempty(atlasViewList))
    atlasView = [];
elseif(length(atlasViewList)==1)
    atlasView = FLAT{atlasViewList{1}};
else
    % More than one of the FLAT{} structures has an Atlases type selected.
    % We ask the user which of the FLAT{} structures they have in mind.
    % We figure this out by asking them to identify the figure number with
    % the selected atlas.  Then, we convert the figure number into the
    % index into the FLAT{} cell array.
    prompt={'Enter the figure number of the FLAT view showing the atlas'};
    def={'1'}; dlgTitle='Figure Showing Atlas'; lineNo=1;
    answer=inputdlg(prompt,dlgTitle,lineNo,def);
    atlasFig = str2num(answer{1});
    idx = fignum2FLATindex(atlasFig);
    if ~isempty(idx),  atlasView = FLAT{idx};
    else
        errordlg('Selected figure number does not have an Atlases data type.');
        atlasView = [];
    end
end

return;


text:
function K = k00(x, xp, hyp, Ubar, Ubarp, dt, i)

ubar = Ubar(:,1);
vbar = Ubar(:,2);

ubarp = Ubarp(:,1);
vbarp = Ubarp(:,2);

Ku0u0 = k00.ku0u0(x, xp, hyp, ubar, vbar, ubarp, vbarp, dt, i);
Ku0v0 = k00.ku0v0(x, xp, hyp, ubar, vbar, ubarp, vbarp, dt, i);
Kv0u0 = Ku0v0';
Kv0v0 = k00.kv0v0(x, xp, hyp, ubar, vbar, ubarp, vbarp, dt, i);

K = [Ku0u0 Ku0v0;
     Kv0u0 Kv0v0];

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
function G = mrdivide(F, a)
%/   DISKFUNV right divide.
%
% F/a divides each component of a DISKFUNV F by the scalar a.
%
% Only allowed to divide by scalars.
%
% See also MLDIVIDE.

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

% Empty check:
if ( ( isempty(F) ) || ( isempty(a) ) )
   G = diskfunv;
   return
end

% Only allow F/a, where a is a scalar:
if ( ~isa(a, 'double') )
    error('CHEBFUN:DISKFUNV:mrdivide:nonScalar', ...
        'Division must be scalar valued.');
end

% Componentwise divide.
G = F;
G.components{1} = mrdivide(F.components{1}, a);
G.components{2} = mrdivide(F.components{2}, a);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function preprocessing()
cd matconvnet;
complienn;
cd ..
end

author: jiangqy


text:
% demo script for splitting the field of view in patches and processing in
% parallel with or without memory mapping. The demo follows the file
% demo_patches.m using the CNMF class object

clear;
%% setup path to file and package
gcp;
path_to_package = '../ca_source_extraction';   % path to the folder that contains the package
addpath(genpath(path_to_package));

filename = '/Users/epnevmatikakis/Documents/Ca_datasets/Neurofinder/neurofinder.02.00/images/neurofinder0200_rig.tif';
        % path to file (assumed motion corrected)

is_memmaped = true;        % choose whether you want to load the file in memory or not

%% create object and load file

CNM = CNMF();
if is_memmaped
    CNM.readFile(filename,is_memmaped);
else
    CNM.readFile(filename,is_memmaped,1,2000); % load only a part of the file due to memory
end

%% set options and create patches

patch_size = [32,32];    % size of each patch along each dimension (optional, default: [32,32])
overlap = [6,6];         % amount of overlap in each dimension (optional, default: [4,4])
K = 10;                  % number of components to be found
gSig = 7;                % std of gaussian kernel (size of neuron)
p = 2;                   % order of autoregressive system (p = 0 no dynamics, p=1 just decay, p = 2, both rise and decay)
gnb = 3;                 % order of background
merge_thr = 0.8;         % merging threshold

options = CNMFSetParms(...
    'd1',CNM.dims(1),'d2',CNM.dims(2),...
    'search_method','dilate',...                % search locations when updating spatial components
    'deconv_method','constrained_foopsi',...    % activity deconvolution method
    'nb',1,...                                  % number of background components per patch
    'gnb',gnb,...                               % number of global background components
    'ssub',2,...
    'tsub',1,...
    'p',p,...                                   % order of AR dynamics
    'merge_thr',merge_thr,...                   % merging threshold
    'gSig',gSig,...
    'spatial_method','regularized',...
    'cnn_thr',0.2,...
    'patch_space_thresh',0.25,...
    'min_SNR',2);

CNM.optionsSet(options);
CNM.gnb = gnb;
CNM.K = K;
CNM.patch_size = patch_size;                % size of each patch along each dimension (optional, default: [32,32])
CNM.overlap = overlap;                      % amount of overlap in each dimension (optional, default: [4,4])
CNM.createPatches();                        % create patches

%% fit all patches
CNM.fitPatches();

%% component classification

CNM.evaluateComponents();   % evaluate spatial components based on their correlation with the data
CNM.CNNClassifier('')       % evaluate spatial components with the CNN classifier
CNM.eventExceptionality();  % evaluate traces
CNM.keepComponents();       % keep the components that are above certain thresholds

%% repeat processing

CNM.updateSpatial();
CNM.updateTemporal();
CNM.extractDFF();            % extract DF/F values.

%% do some plotting
figure;
CNM.correlationImage();
CNM.plotContours();
CNM.plotComponentsGUI();     % display all components

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function submitWeb(partId)
%SUBMITWEB Generates a base64 encoded string for web-based submissions
%   SUBMITWEB() will generate a base64 encoded string so that you can submit your
%   solutions via a web form

  fprintf('==\n== [ml-class] Submitting Solutions | Programming Exercise %s\n==\n', ...
          homework_id());
  if ~exist('partId', 'var') || isempty(partId)
    partId = promptPart();
  end

  % Check valid partId
  partNames = validParts();
  if ~isValidPartId(partId)
    fprintf('!! Invalid homework part selected.\n');
    fprintf('!! Expected an integer from 1 to %d.\n', numel(partNames));
    fprintf('!! Submission Cancelled\n');
    return
  end

  [login] = loginPrompt();
  if isempty(login)
    fprintf('!! Submission Cancelled\n');
    return
  end

  [result] = submitSolution(login, partId, output(partId), ...
                            source(partId));
  result = base64encode(result);

  fprintf('\nSave as submission file [submit_ex%s_part%d.txt]: ', ...
          homework_id(), partId);
  saveAsFile = input('', 's');
  if (isempty(saveAsFile))
    saveAsFile = sprintf('submit_ex%s_part%d.txt', homework_id(), partId);
  end

  fid = fopen(saveAsFile, 'w');
  if (fid)
    fwrite(fid, result);
    fclose(fid);
    fprintf('\nSaved your solutions to %s.\n\n', saveAsFile);
    fprintf(['You can now submit your solutions through the web \n' ...
             'form in the programming exercises. Select the corresponding \n' ...
             'programming exercise to access the form.\n']);

  else
    fprintf('Unable to save to %s\n\n', saveAsFile);
    fprintf(['You can create a submission file by saving the \n' ...
             'following text in a file: (press enter to continue)\n\n']);
    pause;
    fprintf(result);
  end

end

% ================== CONFIGURABLES FOR EACH HOMEWORK ==================

function id = homework_id()
  id = '2';
end

function [partNames] = validParts()
  partNames = { 'Sigmoid Function ', ...
                'Logistic Regression Cost', ...
                'Logistic Regression Gradient', ...
                'Predict', ...
                'Regularized Logistic Regression Cost' ...
                'Regularized Logistic Regression Gradient' ...
                };
end

function srcs = sources()
  % Separated by part
  srcs = { { 'sigmoid.m' }, ...
           { 'costFunction.m' }, ...
           { 'costFunction.m' }, ...
           { 'predict.m' }, ...
           { 'costFunctionReg.m' }, ...
           { 'costFunctionReg.m' } };
end

function out = output(partId)
  % Random Test Cases
  X = [ones(20,1) (exp(1) * sin(1:1:20))' (exp(0.5) * cos(1:1:20))'];
  y = sin(X(:,1) + X(:,2)) > 0;
  if partId == 1
    out = sprintf('%0.5f ', sigmoid(X));
  elseif partId == 2
    out = sprintf('%0.5f ', costFunction([0.25 0.5 -0.5]', X, y));
  elseif partId == 3
    [cost, grad] = costFunction([0.25 0.5 -0.5]', X, y);
    out = sprintf('%0.5f ', grad);
  elseif partId == 4
    out = sprintf('%0.5f ', predict([0.25 0.5 -0.5]', X));
  elseif partId == 5
    out = sprintf('%0.5f ', costFunctionReg([0.25 0.5 -0.5]', X, y, 0.1));
  elseif partId == 6
    [cost, grad] = costFunctionReg([0.25 0.5 -0.5]', X, y, 0.1);
    out = sprintf('%0.5f ', grad);
  end
end

% ========================= SUBMIT HELPERS =========================

function src = source(partId)
  src = '';
  src_files = sources();
  if partId <= numel(src_files)
      flist = src_files{partId};
      for i = 1:numel(flist)
          fid = fopen(flist{i});
          while ~feof(fid)
            line = fgets(fid);
            src = [src line];
          end
          fclose(fid);
          src = [src '||||||||'];
      end
  end
end

function ret = isValidPartId(partId)
  partNames = validParts();
  ret = (~isempty(partId)) && (partId >= 1) && (partId <= numel(partNames));
end

function partId = promptPart()
  fprintf('== Select which part(s) to submit:\n', ...
          homework_id());
  partNames = validParts();
  srcFiles = sources();
  for i = 1:numel(partNames)
    fprintf('==   %d) %s [', i, partNames{i});
    fprintf(' %s ', srcFiles{i}{:});
    fprintf(']\n');
  end
  fprintf('\nEnter your choice [1-%d]: ', ...
          numel(partNames));
  selPart = input('', 's');
  partId = str2num(selPart);
  if ~isValidPartId(partId)
    partId = -1;
  end
end

function [result, str] = submitSolution(email, part, output, source)

  result = ['a:5:{' ...
            p_s('homework') p_s64(homework_id()) ...
            p_s('part') p_s64(part) ...
            p_s('email') p_s64(email) ...
            p_s('output') p_s64(output) ...
            p_s('source') p_s64(source) ...
            '}'];

end

function s = p_s(str)
   s = ['s:' num2str(numel(str)) ':"' str '";'];
end

function s = p_s64(str)
   str = base64encode(str, '');
   s = ['s:' num2str(numel(str)) ':"' str '";'];
end

% =========================== LOGIN HELPERS ===========================

function [login] = loginPrompt()
  % Prompt for password
  [login] = basicPrompt();
end

function [login] = basicPrompt()
  login = input('Login (Email address): ', 's');
end

% =========================== Base64 Encoder ============================
% Thanks to Peter John Acklam
%

function y = base64encode(x, eol)
%BASE64ENCODE Perform base64 encoding on a string.
%
%   BASE64ENCODE(STR, EOL) encode the given string STR.  EOL is the line ending
%   sequence to use; it is optional and defaults to '\n' (ASCII decimal 10).
%   The returned encoded string is broken into lines of no more than 76
%   characters each, and each line will end with EOL unless it is empty.  Let
%   EOL be empty if you do not want the encoded string broken into lines.
%
%   STR and EOL don't have to be strings (i.e., char arrays).  The only
%   requirement is that they are vectors containing values in the range 0-255.
%
%   This function may be used to encode strings into the Base64 encoding
%   specified in RFC 2045 - MIME (Multipurpose Internet Mail Extensions).  The
%   Base64 encoding is designed to represent arbitrary sequences of octets in a
%   form that need not be humanly readable.  A 65-character subset
%   ([A-Za-z0-9+/=]) of US-ASCII is used, enabling 6 bits to be represented per
%   printable character.
%
%   Examples
%   --------
%
%   If you want to encode a large file, you should encode it in chunks that are
%   a multiple of 57 bytes.  This ensures that the base64 lines line up and
%   that you do not end up with padding in the middle.  57 bytes of data fills
%   one complete base64 line (76 == 57*4/3):
%
%   If ifid and ofid are two file identifiers opened for reading and writing,
%   respectively, then you can base64 encode the data with
%
%      while ~feof(ifid)
%         fwrite(ofid, base64encode(fread(ifid, 60*57)));
%      end
%
%   or, if you have enough memory,
%
%      fwrite(ofid, base64encode(fread(ifid)));
%
%   See also BASE64DECODE.

%   Author:      Peter John Acklam
%   Time-stamp:  2004-02-03 21:36:56 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   if isnumeric(x)
      x = num2str(x);
   end

   % make sure we have the EOL value
   if nargin < 2
      eol = sprintf('\n');
   else
      if sum(size(eol) > 1) > 1
         error('EOL must be a vector.');
      end
      if any(eol(:) > 255)
         error('EOL can not contain values larger than 255.');
      end
   end

   if sum(size(x) > 1) > 1
      error('STR must be a vector.');
   end

   x   = uint8(x);
   eol = uint8(eol);

   ndbytes = length(x);                 % number of decoded bytes
   nchunks = ceil(ndbytes / 3);         % number of chunks/groups
   nebytes = 4 * nchunks;               % number of encoded bytes

   % add padding if necessary, to make the length of x a multiple of 3
   if rem(ndbytes, 3)
      x(end+1 : 3*nchunks) = 0;
   end

   x = reshape(x, [3, nchunks]);        % reshape the data
   y = repmat(uint8(0), 4, nchunks);    % for the encoded data

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Split up every 3 bytes into 4 pieces
   %
   %    aaaaaabb bbbbcccc ccdddddd
   %
   % to form
   %
   %    00aaaaaa 00bbbbbb 00cccccc 00dddddd
   %
   y(1,:) = bitshift(x(1,:), -2);                  % 6 highest bits of x(1,:)

   y(2,:) = bitshift(bitand(x(1,:), 3), 4);        % 2 lowest bits of x(1,:)
   y(2,:) = bitor(y(2,:), bitshift(x(2,:), -4));   % 4 highest bits of x(2,:)

   y(3,:) = bitshift(bitand(x(2,:), 15), 2);       % 4 lowest bits of x(2,:)
   y(3,:) = bitor(y(3,:), bitshift(x(3,:), -6));   % 2 highest bits of x(3,:)

   y(4,:) = bitand(x(3,:), 63);                    % 6 lowest bits of x(3,:)

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Now perform the following mapping
   %
   %   0  - 25  ->  A-Z
   %   26 - 51  ->  a-z
   %   52 - 61  ->  0-9
   %   62       ->  +
   %   63       ->  /
   %
   % We could use a mapping vector like
   %
   %   ['A':'Z', 'a':'z', '0':'9', '+/']
   %
   % but that would require an index vector of class double.
   %
   z = repmat(uint8(0), size(y));
   i =           y <= 25;  z(i) = 'A'      + double(y(i));
   i = 26 <= y & y <= 51;  z(i) = 'a' - 26 + double(y(i));
   i = 52 <= y & y <= 61;  z(i) = '0' - 52 + double(y(i));
   i =           y == 62;  z(i) = '+';
   i =           y == 63;  z(i) = '/';
   y = z;

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Add padding if necessary.
   %
   npbytes = 3 * nchunks - ndbytes;     % number of padding bytes
   if npbytes
      y(end-npbytes+1 : end) = '=';     % '=' is used for padding
   end

   if isempty(eol)

      % reshape to a row vector
      y = reshape(y, [1, nebytes]);

   else

      nlines = ceil(nebytes / 76);      % number of lines
      neolbytes = length(eol);          % number of bytes in eol string

      % pad data so it becomes a multiple of 76 elements
      y = [y(:) ; zeros(76 * nlines - numel(y), 1)];
      y(nebytes + 1 : 76 * nlines) = 0;
      y = reshape(y, 76, nlines);

      % insert eol strings
      eol = eol(:);
      y(end + 1 : end + neolbytes, :) = eol(:, ones(1, nlines));

      % remove padding, but keep the last eol string
      m = nebytes + neolbytes * (nlines - 1);
      n = (76+neolbytes)*nlines - neolbytes;
      y(m+1 : n) = '';

      % extract and reshape to row vector
      y = reshape(y, 1, m+neolbytes);

   end

   % output is a character array
   y = char(y);

end

author: schneems

repo: Octave

save_path: github-repos/MATLAB/schneems-Octave


text:
function d = AngDiff(x, y)
%ANGDIFF  Compute angle difference accurately
%
%   D = ANGDIFF(X, Y) computes Y - X, reduces the result to (-180,180] and
%   rounds the result.  X and Y must be in [-180,180].  X and Y can be any
%   compatible shapes.

  [d, t] = sumx(-x, y);
  c = (d - 180) + t > 0;
  d(c) = (d(c) - 360) + t(c);
  c = (d + 180) + t <= 0;
  d(c) = (d(c) + 360) + t(c);
end

author: Sable

repo: mcbench-benchmarks


text:
function SO3F = times(SO3F1,SO3F2)
% overloads |SO3F1 .* SO3F2|
%
% Syntax
%   sF = SO3F1 .* SO3F2
%   sF = a .* SO3F2
%   sF = SO3F1 .* a
%
% Input
%  SO3F1, SO3F2 - @SO3Fun
%  a - double
%
% Output
%  SO3F - @SO3Fun
%

% uniform component .* SO3Fun
if isa(SO3F1,'SO3FunRBF') && isempty(SO3F1.center) && isa(SO3F2,'SO3Fun')
  SO3F = SO3F1.c0 .* SO3F2;
  return
end
if isa(SO3F2,'SO3FunRBF') && isempty(SO3F2.center) && isa(SO3F1,'SO3Fun')
  SO3F = SO3F1 .* SO3F2.c0;
  return
end

if isnumeric(SO3F1)
  dim = length(size(SO3F1));
  SO3F = SO3FunHandle(@(rot) permute(SO3F1,[dim+1 1:dim]) .* SO3F2.eval(rot),SO3F2.CS,SO3F2.SS);
  return
end

if isnumeric(SO3F2)
  SO3F = SO3F2 .* SO3F1;
  return
end

ensureCompatibleSymmetries(SO3F1,SO3F2);
SO3F = SO3FunHandle(@(rot) SO3F1.eval(rot) .* SO3F2.eval(rot),SO3F1.SRight,SO3F1.SLeft);

end

author: mtex-toolbox

repo: mtex

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


text:
function using_gpu()
% Manopt example on how to use GPU with manifold factories that allow it.
%
% We are still working on this feature, and so far only few factories have
% been adapted to work on GPU. But the adaptations are rather easy. If
% there is a manifold you'd like to use on GPU, let us know via the forum
% on http://www.manopt.org, we'll be happy to help!
%
% See also: spherefactory stiefelfactory grassmannfactory complexcirclefactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Aug. 3, 2018.
% Contributors:
% Change log:

    if exist('OCTAVE_VERSION', 'builtin')
        warning('manopt:usinggpu', 'Octave does not handle GPUs at this time.');
        return;
    end

    if gpuDeviceCount() <= 0
        warning('manopt:usinggpu', 'No GPU available: cannot run example.');
        return;
    end

    % Construct a large problem to illustrate the use of GPU.
    % Below, we will compute p left-most eigenvectors of A (symmetric).
    % On a particular test computer, we found that for n = 100, 1000, CPU
    % is faster, but for n = 10000, GPU tends to be 10x faster.
    p = 3;
    n = 10000;
    A = randn(n);
    A = A+A';

    inner = @(U, V) U(:)'*V(:);

    % First, setup and run the optimization problem on the CPU.
    problem.M = grassmannfactory(n, p, 1); % 1 copy of Grassmann(n, p)
    problem.cost = @(X) .5*inner(X, A*X);  % Rayleigh quotient to be minimized
    problem.egrad = @(X) A*X;              % Could use caching to save here
    problem.ehess = @(X, Xdot) A*Xdot;
    X0 = problem.M.rand();                 % Random initial guess
    tic_cpu = tic();
    X_cpu = trustregions(problem, X0);     % run any solver
    time_cpu = toc(tic_cpu);

    % Then, move the data to the GPU, redefine the problem using the moved
    % data, activate the GPU flag in the factory, and run it again.
    A = gpuArray(A);
    problem.M = grassmannfactory(n, p, 1, true); % true is the GPU flag;
    problem.cost = @(X) .5*inner(X, A*X);        % Code for cost and gradient etc.
    problem.egrad = @(X) A*X;                    % basically didn't change, but
    problem.ehess = @(X, Xdot) A*Xdot;           % operates on gpuArrays now.
    X0 = gpuArray(X0);
    tic_gpu = tic();
    X_gpu = trustregions(problem, X0);
    time_gpu = toc(tic_gpu);

    fprintf('Total time CPU: %g\nTotal time GPU: %g\nSolution difference: %g\n', ...
            time_cpu, time_gpu, norm(X_cpu - X_gpu, 'fro'));

end

author: NicolasBoumal

repo: manopt


text:
function [Rlat,Rlon,alt] = LatLongCalcSingle_ADI(msg, inputLat, inputLong)
% Calculate latitude, longitude and altitude from message bits
% Copyright 2010, The MathWorks, Inc.

persistent NL
persistent latzones
persistent Dlat0
persistent Dlat1
persistent latOffset0
persistent latOffset1

if isempty(NL)
    Dlat0 = 360/(4*15-0);
    latOffset0 = floor(inputLat/Dlat0);
    Dlat1 = 360/(4*15-1);
    latOffset1 = floor(inputLat/Dlat1);
    NL=2:59;
    latzones = [(180/pi)*acos(sqrt((1-cos(pi/2/15))./(1-cos(2*pi./NL)))) 0];
end

% Altitude calculation
q = msg(48);
if q == 0
    af = 100;
else
    af = 25;
end
altBits=[msg(41:47);msg(49:52)]';
alt = altBits*[1024;512;256;128;64;32;16;8;4;2;1]*af - 1000;

evenOdd1 = msg(54);
latBits = msg(55:71)';
longBits = msg(72:88)';
la1 = latBits*[65536;32768;16384;8192;4096;2048;1024;512;256;128;64;32;16;8;4;2;1];
lo1 = longBits*[65536;32768;16384;8192;4096;2048;1024;512;256;128;64;32;16;8;4;2;1];

% Technically you need both even and odd messages to calculate lat/long
% unambiguously. For this code, use a single message and then check to see
% if the lat/long values are reasonable. If not, change the lat/long base
% factors (LL.a1, LL.a2, etc.) and recompute.

% Latitude calculation
if evenOdd1 == 0
    Rlat = Dlat0*(latOffset0 + la1/131072);
else
    Rlat = Dlat1*(latOffset1 + la1/131072);
end

% Compare latitude to known location. If it's off by more than two degrees,
% use new base factors.
if Rlat > inputLat+2
    if strcmp(evenOdd1,'Even')
        Rlat = Dlat0*(latOffset0 - 1 + la1/131072);
    else
        Rlat = Dlat1*(latOffset1 - 1 + la1/131072);
    end
elseif Rlat < inputLat-2
    if strcmp(evenOdd1,'Even')
        Rlat = Dlat0*(latOffset0 + 1 + la1/131072);
    else
        Rlat = Dlat1*(latOffset1 + 1 + la1/131072);
    end
end

% Based on latitude, calculate longitude
NL0 = find(latzones<Rlat,1,'first');
ni0 = NL0;
ni1 = NL0 - 1;

Dlon0 = 360/ni0;
longOffset0 = floor(inputLong/Dlon0);
Dlon1 = 360/ni1;
longOffset1 = floor(inputLong/Dlon1);

if evenOdd1 == 0
    Rlon = Dlon0*(longOffset0 + lo1/131072);
else
    Rlon = Dlon1*(longOffset1 + lo1/131072);
end

% Compare longitude to known location. If it's off by more than two
% degrees, use new base factors.
if Rlon > inputLong+2
    if strcmp(evenOdd1,'Even')
        Rlon = Dlon0*(longOffset0 - 1 + lo1/131072);
    else
        Rlon = Dlon1*(longOffset1 - 1 + lo1/131072);
    end
elseif Rlon < inputLong-2
    if strcmp(evenOdd1,'Even')
        Rlon = Dlon0*(longOffset0 + 1 + lo1/131072);
    else
        Rlon = Dlon1*(longOffset1 + 1 + lo1/131072);
    end
end

% disp(sprintf('Plane is at altitude %d\nLatitude value: %d\nLongitude value: %d', alt, la1, lo1));

% GoogleMap(aircraftID, alt1, Rlat, Rlon)

author: analogdevicesinc


text:
function N = uminus(N)
%UMINUS   Unitary minus for CHEBOP2 objects.

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

% Negate all the variables coefficents:
A = N.coeffs;
if ( iscell(A) )
    for jj = 1:size(A, 1)
        for kk = 1:size(A, 2)
            A{jj,kk} = -A{jj,kk};
        end
    end
else
    A = -A;
end

% Do not negate the BCs:
if ( ~isempty(N.lbc) || ~isempty(N.rbc) || ~isempty(N.ubc) || ~isempty(N.dbc) )
    warning('CHEBFUN:CHEBOP2:unimus:BCs', ...
        'Operator has BCs. These were not negated.');
end

% Update the properties of the CHEBOP2:
N.coeffs = A;
op = N.op;
N.op = @(u) -op(u);

end


text:
%> @brief Split 'str_in' into lines and prefix each line with 'prefix'
%> @author Eric Cousineau <eacousineau@gmail.com>, AMBER Lab, under Dr.
%> Aaron Ames
function [str] = indent(str_in, prefix)
% Using sprintf for sprintf:
% http://www.mathworks.com/matlabcentral/newsreader/view_thread/249016
lines = strsplit(sprintf(str_in), sprintf('\n'));
str = sprintf([prefix, strjoin(lines, ['\n', prefix])]);
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
% STD_CENTROID - compute cluster centroid in EEGLAB dataset STUDY.
%                  Compute and store the centroid(s) (i.e., mean(s))
%                  for some combination of six measures on specified
%                  clusters in a STUDY. Possible measures include: scalp
%                  maps, ERPs, spectra, ERSPs, ITCs, dipole_locations
% Usage:
%        >> [STUDY, centroid] = std_centroid(STUDY, ALLEEG, ...
%                                              clusters, measure1, measure2, ...);
%
% Inputs:
%   STUDY        - STUDY set
%   ALLEEG       - ALLEEG dataset vector (else an EEG dataset) containing the STUDY
%                  datasets, typically created using LOAD_ALLEEG.
%   clusters     - [vector] of cluster indices. Computes measure means for the
%                  specified clusters. {deffault|[]: compute means for all
%                  STUDY clusters}
%   measure(s)   - ['erp'|'spec'|'scalp'|'dipole'|'itc'|'ersp'].
%                  The measures(s) for which to calculate the cluster centroid(s):
%                     'erp'    ->  mean ERP of each cluster.
%                     'dipole' ->  mean dipole of each cluster.
%                     'spec'   ->  mean spectrum of each cluster (baseline removed).
%                     'scalp'  ->  mean topoplot scalp map of each cluster.
%                     'ersp'   ->  mean ERSP of each cluster.
%                     'itc'    ->  mean ITC of each cluster.
%                  If [], re-compute the centroid for whichever centroids
%                  have previously been computed.
% Outputs:
%   STUDY        - input STUDY structure with computed centroids added.
%                  If the requested centroids already exist, overwrites them.
%   centroid     - cell array of centroid structures, each cell corrasponding
%                  to a different cluster requested in 'clusters' (above).
%                  fields of 'centroid' may include centroid.erp, centroid.dipole,
%                  etc. (as above). The structure is similar as the output
%                  of the STD_READDATA function (with some fields
%                  about the cluster name and index missing).
% Examples:
%
%   >> [STUDY, centroid] = std_centroid(STUDY, ALLEEG,[], 'scalp');
%   % For each of the clusters in STUDY, compute a mean scalp map.
%   % The centroids are saved in the STUDY structure as entries in array
%   % STUDY.cluster(k).centroid.scalp. The centroids are also returned in
%   % a cell array the size of the clusters (i.e., in: centroid(k).scalp).
%
%   >> [STUDY, centroid] = std_centroid(STUDY, ALLEEG,5,'spec','scalp');
%   % Same as above, but now compute only two centroids for Cluster 5.
%   % The returned 'centroid' has two fields: centroid.scalp and centroid.spec
%
% Authors: Hilit Serby & Arnaud Delorme, SCCN, INC, UCSD, Feb 03, 2005

% Copyright (C) Hilit Serby, SCCN, INC, UCSD, Feb 03, 2005, hilit@sccn.ucsd.edu
%
% 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.

% Coding notes: Useful information on functions and global variables used.

function [STUDY, centroid] = std_centroid(STUDY,ALLEEG, clsind, varargin);

 if nargin < 3
     help std_centroid;
     return
 end

 if isempty(clsind)
     for k = 2: length(STUDY.cluster) %don't include the ParentCluster
         if ~strncmpi('Notclust',STUDY.cluster(k).name,8)
             % don't include 'Notclust' clusters
             clsind = [clsind k];
         end
     end
 end
 %default values
erpC =0;
specC =0 ;
scalpC = 0;
dipoleC = 0;
itcC = 0;
erspC = 0;

commands = {};
if isempty(varargin)
    if isfield(STUDY.cluster(clsind(1)).centroid,'scalp')
        commands{end+1} = 'scalp';
    end
    if isfield(STUDY.cluster(clsind(1)).centroid,'spec')
        commands{end+1} = 'spec';
    end
    if isfield(STUDY.cluster(clsind(1)).centroid,'erp')
        commands{end+1} = 'erp';
    end
    if isfield(STUDY.cluster(clsind(1)).centroid,'ersp')
        commands{end+1} = 'ersp';
    end
    if isfield(STUDY.cluster(clsind(1)).centroid,'itc')
        commands{end+1} = 'itc';
    end
    if isfield(STUDY.cluster(clsind(1)).centroid,'dipole')
        commands{end+1} = 'dipole';
    end
else
    commands = varargin;
end

Ncond = length(STUDY.condition);
if Ncond == 0
    Ncond = 1;
end
centroid = cell(length(clsind),1);
fprintf('Computing ');
for k = 1:length(clsind)
    for l = 1:Ncond
        for ind = 1:length(commands)
            ctr = commands{ind};
            switch ctr
                case 'scalp'
                    centroid{k}.scalp = 0;
                    scalpC = 1;
                    if (l ==1) && (k ==1)
                        fprintf('scalp ');
                    end
                case 'erp'
                    centroid{k}.erp{l} = 0;
                    erpC = 1;
                    if (l ==1) && (k ==1)
                        fprintf('erp ');
                    end
                case 'spec'
                    centroid{k}.spec{l} = 0;
                    specC = 1;
                    if (l ==1) && (k ==1)
                        fprintf('spectrum ');
                    end
                case 'ersp'
                    centroid{k}.ersp{l} = 0;
                    centroid{k}.ersp_limits{l} = 0;
                    erspC =1;
                    if (l ==1) && (k ==1)
                        fprintf('ersp ');
                    end
                case 'itc'
                    centroid{k}.itc{l} = 0;
                    centroid{k}.itc_limits{l} = 0;
                    itcC = 1;
                    if (l ==1) && (k ==1)
                        fprintf('itc ');
                    end
                case 'dipole'
                    dipoleC =1;
                    if (l ==1) && (k ==1)
                        fprintf('dipole ');
                    end
            end
        end
    end
end
fprintf('centroid (only done once)\n');
if itcC || erspC || specC || erpC || scalpC
    for clust = 1:length(clsind) %go over all requested clusters
        for cond = 1:Ncond %compute for all conditions
            for k = 1:length(STUDY.cluster(clsind(clust)).comps) % go through all components
                comp  = STUDY.cluster(clsind(clust)).comps(k);
                abset = STUDY.cluster(clsind(clust)).sets(cond,k);
                if scalpC && cond == 1  %scalp centroid, does not depend on condition
                    grid = std_readtopo(ALLEEG, abset, comp);
                    if isempty(grid)
                        return;
                    end
                    centroid{clust}.scalp = centroid{clust}.scalp + grid;
                end
                if erpC %erp centroid
                    [erp, t] = std_readerp(ALLEEG, abset, comp, STUDY.preclust.erpclusttimes);
                    fprintf('.');
                    if isempty(erp)
                        return;
                    end
                    if (cond==1) && (k==1)
                        all_erp = zeros(length(erp),length(STUDY.cluster(clsind(clust)).comps));
                    end
                    all_erp(:,k) = erp';
                    if k == length(STUDY.cluster(clsind(clust)).comps)
                        [all_erp pol] = std_comppol(all_erp);
                        centroid{clust}.erp{cond} = mean(all_erp,2);
                        centroid{clust}.erp_times = t;
                    end
                end
                if specC %spec centroid
                    [spec, f] = std_readspec(ALLEEG, abset, comp, STUDY.preclust.specclustfreqs);
                    fprintf('.');
                    if isempty(spec)
                        return;
                    end
                    centroid{clust}.spec{cond} = centroid{clust}.spec{cond} + spec;
                    centroid{clust}.spec_freqs = f;
                end
                if erspC %ersp centroid
                    fprintf('.');
                    if cond == 1
                        tmpabset = STUDY.cluster(clsind(clust)).sets(:,k);
                        [ersp, logfreqs, timevals] = std_readersp(ALLEEG, tmpabset, comp, STUDY.preclust.erspclusttimes, ...
                                                                STUDY.preclust.erspclustfreqs );
                        if isempty(ersp)
                            return;
                        end
                        for m = 1:Ncond
                            centroid{clust}.ersp{m} = centroid{clust}.ersp{m} + ersp(:,:,m);
                            centroid{clust}.ersp_limits{m} = max(floor(max(max(abs(ersp(:,:,m))))), centroid{clust}.ersp_limits{m});
                        end
                        centroid{clust}.ersp_freqs  = logfreqs;
                        centroid{clust}.ersp_times = timevals;
                    end
                end
                if itcC %itc centroid
                    fprintf('.');
                    [itc, logfreqs, timevals] = std_readitc(ALLEEG, abset, comp, STUDY.preclust.erspclusttimes, ...
                                                                STUDY.preclust.erspclustfreqs );
                    if isempty(itc)
                        return;
                    end
                    centroid{clust}.itc{cond} = centroid{clust}.itc{cond} + itc;
                    centroid{clust}.itc_limits{cond} = max(floor(max(max(abs(itc)))), centroid{clust}.itc_limits{cond}); %ersp image limits
                    centroid{clust}.itc_freqs  = logfreqs;
                    centroid{clust}.itc_times = timevals;
                end
            end
        end
        if ~scalpC
            fprintf('\n');
        end
	end
end

if dipoleC %dipole centroid
    for clust = 1:length(clsind)
        max_r = 0;
        len = length(STUDY.cluster(clsind(clust)).comps);
        tmppos = 0;
        tmpmom = 0;
        tmprv = 0;
        ndip = 0;
        for k = 1:len
            fprintf('.');
            comp  = STUDY.cluster(clsind(clust)).comps(k);
            abset = STUDY.cluster(clsind(clust)).sets(1,k);
            if ~isfield(ALLEEG(abset), 'dipfit')
               warndlg2(['No dipole information available in dataset ' num2str(abset) ], 'Aborting compute centroid dipole');
               return;
            end
            if ~isempty(ALLEEG(abset).dipfit.model(comp).posxyz)
                ndip = ndip +1;
                tmppos = tmppos + ALLEEG(abset).dipfit.model(comp).posxyz;
                tmpmom = tmpmom + ALLEEG(abset).dipfit.model(comp).momxyz;
                tmprv = tmprv + ALLEEG(abset).dipfit.model(comp).rv;
                if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical')
                   if isfield(ALLEEG(abset).dipfit, 'hdmfile') %dipfit 2 spherical model
                       load('-mat', ALLEEG(abset).dipfit.hdmfile);
                       max_r = max(max_r, max(vol.r));
                   else % old version of dipfit
                       max_r = max(max_r,max(ALLEEG(abset).dipfit.vol.r));
                   end
               end
            end
        end
        centroid{clust}.dipole.posxyz =  tmppos/ndip;
        centroid{clust}.dipole.momxyz =  tmpmom/ndip;
        centroid{clust}.dipole.rv =  tmprv/ndip;
        if strcmpi(ALLEEG(abset).dipfit.coordformat, 'spherical') && (~isfield(ALLEEG(abset).dipfit, 'hdmfile')) %old dipfit
            centroid{clust}.dipole.maxr = max_r;
        end
        STUDY.cluster(clsind(clust)).centroid.dipole = centroid{clust}.dipole;
   end
end

%update STUDY
for clust =  1:length(clsind) %go over all requested clusters
    for cond  = 1:Ncond
        ncomp = length(STUDY.cluster(clsind(clust)).comps);
        if scalpC && cond == 1%scalp centroid
            centroid{clust}.scalp  = centroid{clust}.scalp/ncomp;
            STUDY.cluster(clsind(clust)).centroid.scalp = centroid{clust}.scalp ;
        end
        if erpC
            STUDY.cluster(clsind(clust)).centroid.erp{cond} = centroid{clust}.erp{cond};
		    STUDY.cluster(clsind(clust)).centroid.erp_times = centroid{clust}.erp_times;
        end
		if specC
            centroid{clust}.spec{cond} = centroid{clust}.spec{cond}/ncomp;
            STUDY.cluster(clsind(clust)).centroid.spec{cond} = centroid{clust}.spec{cond};
		    STUDY.cluster(clsind(clust)).centroid.spec_freqs = centroid{clust}.spec_freqs;
        end
        if erspC %ersp centroid
            centroid{clust}.ersp{cond} = centroid{clust}.ersp{cond}/ncomp;
            STUDY.cluster(clsind(clust)).centroid.ersp{cond} = centroid{clust}.ersp{cond};
            STUDY.cluster(clsind(clust)).centroid.ersp_limits{cond} = floor(0.75*centroid{clust}.ersp_limits{cond});
            %[round(0.9*min(cell2mat({centroid{clust}.ersp_limits{cond,:}})))  round(0.9*max(cell2mat({centroid{clust}.ersp_limits{cond,:}})))];
            STUDY.cluster(clsind(clust)).centroid.ersp_freqs = centroid{clust}.ersp_freqs;
            STUDY.cluster(clsind(clust)).centroid.ersp_times = centroid{clust}.ersp_times;
        end
        if itcC
            centroid{clust}.itc{cond} = centroid{clust}.itc{cond}/ncomp;
            STUDY.cluster(clsind(clust)).centroid.itc{cond} = centroid{clust}.itc{cond} ;
            STUDY.cluster(clsind(clust)).centroid.itc_limits{cond} = floor(0.75*centroid{clust}.itc_limits{cond});%round(0.9*max(cell2mat({centroid{clust}.itc_limits{cond,:}})));
            STUDY.cluster(clsind(clust)).centroid.itc_freqs = centroid{clust}.itc_freqs;
            STUDY.cluster(clsind(clust)).centroid.itc_times = centroid{clust}.itc_times;
        end

    end
end
fprintf('\n');

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
  function [T, reuse] = Gmri_gram(ob, W, reuse)
%|function [T, reuse] = Gmri_gram(ob, W, reuse)
%|
%| build Toeplitz-like gram-matrix object for G'WG gram matrix
%| to be called by build_gram() indirectly rather than directly by user!
%|
%| Jeff Fessler

if nargin < 1, ir_usage, end
if streq(ob, 'test'), Gmri_gram_test, return, end

if isempty(W)
	wi = 1; % default unweighted

% 1D column vector or scalar:
elseif isnumeric(W) && ndims(W) == 2 && size(W,2) == 1
	wi = W;
elseif isa(W, 'Fatrix') && streq(W.caller, 'diag_sp')
	wi = W.arg.diag;
elseif isa(W, 'fatrix2') && streq(W.caller, 'Gdiag')
	wi = W.arg.diag;
else
	fail('Gmri_gram requires W to be Gdiag or diag_sp or wi array')
end

if ~isreal(wi), fail('only real wi supported; otherwise T not Hermitian'), end

arg2.wi = wi;
arg2.Gmri = ob;
arg2.new_zmap = @Gmri_gram_new_zmap;

T = Gmri_gram_work(arg2, reuse);

% Gmri_gram_work()
function T = Gmri_gram_work(arg2, reuse)

arg1 = arg2.Gmri.arg;
wi = arg2.wi;

if isempty(arg1.zmap)
	T = build_gram(arg1.Gnufft, wi .* abs(arg1.basis.transform).^2);
else
	LL = ncol(arg1.aB);
	arg2.T = cell(LL,1);
	reuse = [];
	for ll=1:LL
		wl = arg1.aB(:,ll) .* wi .* abs(arg1.basis.transform).^2;
		[arg2.T{ll} reuse] = build_gram(arg1.Gnufft, wl, reuse);
	end

	arg2.dim = arg1.dim([2 2]); % [np np]

	switch class(arg2.Gmri)
	case 'Fatrix'
		T = Fatrix(arg2.dim, arg2, ...
			'caller', [mfilename '.Gmri_gram'], ...
			'forw', @Gmri_zmap_forw_Fatrix, ...
			'back', @Gmri_zmap_forw_Fatrix); % trick: Hermitian

	case 'fatrix2'
		T = fatrix2('arg', arg2, ...
			'imask', arg2.Gmri.mask, ...
			'omask', arg2.Gmri.mask, ...
			'idim', arg2.Gmri.Nd, ...
			'odim', arg2.Gmri.Nd, ...
			'forw', @Gmri_zmap_forw, ...
			'back', @Gmri_zmap_forw); % trick: because Hermitian

	otherwise
		fail('unknown class "%s"', class(arg2.Gmri))
	end
end

% Gmri_gram_new_zmap()
% update Toeplitz-like gram-matrix Fatrix object for new zmap
function T = Gmri_gram_new_zmap(T, varargin) % (ti, zmap, L, aL)
T.arg.Gmri = T.arg.Gmri.arg.new_zmap(T.arg.Gmri, varargin{:}); % yikes!
T = Gmri_gram_work(T.arg, []);

% Gmri_zmap_forw_Fatrix()
% y = T * x
function y = Gmri_zmap_forw_Fatrix(arg2, x)

arg1 = arg2.Gmri.arg;

if size(x,1) ~= arg2.dim(2)
	x = reshapee(x, prod(arg1.Nd), []); % [(N) (nc)] to [*N *nc]
	x = x(arg1.mask,:); % [np *nc]
end
nc = ncol(x);

LL = ncol(arg1.aB);
y = 0;
for ll=1:LL
	tmp = repmat(arg1.aCt(:,ll), [1 nc]) .* x;
	tmp = arg2.T{ll} * tmp;
	tmp = repmat(conj(arg1.aCt(:,ll)), [1 nc]) .* tmp;
	y = y + tmp;
%	y = y + conj(arg1.aCt(:,ll)) .* (arg2.T{ll} * (arg1.aCt(:,ll) .* x));
end

% Gmri_zmap_forw()
function y = Gmri_zmap_forw(arg2, x)
y = Gmri_zmap_forw_Fatrix(arg2, x);
y = embed(y, arg2.Gmri.arg.mask); % required for fatrix2

% Gmri_gram_test()
function Gmri_gram_test
ig = image_geom('nx', 6, 'ny', 8, 'dx', 1, 'offsets', 'dsp');
ig.mask = ellipse_im(ig) > 0;
x = ellipse_im(ig, 'shepplogan-emis', 'type', 'slow');
kspace = mri_trajectory('spiral1', {}, [ig.nx ig.ny], ig.fov);
ti = linspace(0, 10e-3, size(kspace,1));
zmap = 20 * ig.ones + 2i * pi * 10;
L = 8;
%A = Gmri(kspace, ig.mask, 'exact', 1, 'n_shift', ig.dim/2); % perfect gram
%A = Gmri(kspace, ig.mask, 'class', 'fatrix2'); % gram mult not so well matched
A = Gmri(kspace, ig.mask, 'L', L, 'ti', ti, 'zmap', zmap);
T = build_gram(A);
fatrix2_tests(T, 'complex', 1)
x1 = A' * (A * x);
x1 = ig.embed(x1);
x2 = T * x;
im plc 1 3
im(1, x1)
im(2, x2)
im(3, x2 - x1)
try
	equivs(x2, x1, 'thresh', 5e-3) % todo: why so big - due to nufft?
%	equivs(x2, x1)
catch
	keyboard
end
% fatrix2_tests(A, 'complex', 1)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [Ay, by]  = makeAy(baseMVA, ng, gencost, pgbas, qgbas, ybas)
%MAKEAY  Make the A matrix and RHS for the CCV formulation.
%   [AY, BY]  = MAKEAY(BASEMVA, NG, GENCOST, PGBAS, QGBAS, YBAS)
%
%   Constructs the parameters for linear "basin constraints" on Pg, Qg
%   and Y used by the CCV cost formulation, expressed as
%
%       AY * X <= BY
%
%   where X is the vector of optimization variables. The starting index
%   within the X vector for the active, reactive sources and the Y
%   variables should be provided in arguments PGBAS, QGBAS, YBAS. The
%   number of generators is NG.
%
%   Assumptions: All generators are in-service.  Filter any generators
%   that are offline from the GENCOST matrix before calling MAKEAY.
%   Efficiency depends on Qg variables being after Pg variables, and
%   the Y variables must be the last variables within the vector X for
%   the dimensions of the resulting AY to be conformable with X.
%
%   Example:
%       [Ay, by]  = makeAy(baseMVA, ng, gencost, pgbas, qgbas, ybas);

%   MATPOWER
%   Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC)
%   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;

% find all pwl cost rows in gencost, either real or reactive
iycost = find(gencost(:, MODEL) == PW_LINEAR);

% this is the number of extra "y" variables needed to model those costs
ny = size(iycost, 1);

if ny == 0
   Ay = sparse([], [], [], 0, ybas+ny-1, 0);
   by = [];
   return
end

% if p(i),p(i+1),c(i),c(i+1) define one of the cost segments, then
% the corresponding constraint on Pg (or Qg) and Y is
%                                             c(i+1) - c(i)
%  Y   >=   c(i) + m * (Pg - p(i)),      m = ---------------
%                                             p(i+1) - p(i)
%
% this becomes   m * Pg - Y   <=   m*p(i) - c(i)

% Form A matrix.  Use two different loops, one for the Pg/Qg coeffs,
% then another for the y coeffs so that everything is filled in the
% same order as the compressed column sparse format used by MATLAB;
% this should be the quickest.

m = sum(gencost(iycost, NCOST));  % total number of cost points
Ay = sparse([], [], [], m-ny, ybas+ny-1, 2*(m-ny));
by = [];
% First fill the Pg or Qg coefficients (since their columns come first)
% and the rhs
k = 1;
for i=iycost'
   ns = gencost(i, NCOST);                % # of cost points; segments = ns-1
   p = gencost(i, COST:2:COST+2*ns-1) / baseMVA;
   c = gencost(i, COST+1:2:COST+2*ns);
   m = diff(c) ./ diff(p);                % slopes for Pg (or Qg)
   if any(diff(p) == 0)
     fprintf('\nmakeAy: bad x axis data in row %i of gencost matrix\n',i);
   end
   b = m .* p(1:ns-1) - c(1:ns-1);        % and rhs
   by = [by;  b'];
   if i > ng
     sidx = qgbas + (i-ng) - 1;           % this was for a q cost
   else
     sidx = pgbas + i - 1;                % this was for a p cost
   end
   Ay(k:k+ns-2, sidx) = m';
   k = k + ns - 1;
end
% Now fill the y columns with -1's
k = 1;
j = 1;
for i=iycost'
   ns = gencost(i, NCOST);
   Ay(k:k+ns-2, ybas+j-1) = -ones(ns-1,1);
   k = k + ns - 1;
   j = j + 1;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [strika,dipa,rakea,dipdia,strikb,dipb,rakeb,dipdib,ierr] = pt2pl(trendp,plungp,trendt,plungt)

    % compute strike dip and rake (and dip direction) of two nodal planes from trend and plung of P and T axes
    %
    %     usage:
    %     call pt2pl(trendp,plungp,trendt,plungt,strika,dipa,rakea
    %    1,dipdia,strikb,dipb,rakeb,dipdib,ierr)
    %
    %     arguments:
    %     trendp         trend of P axis in degrees (INPUT)
    %     plungp         plunge of P axis in degrees (INPUT)
    %     trendt         trend of P axis in degrees (INPUT)
    %     plungt         plunge of P axis in degrees (INPUT)
    %     strika         strike angle of first nodal plane in degrees (OUTPUT)
    %     dipa           dip angle of first nodal plane in degrees (OUTPUT)
    %     rakea          rake angle of first nodal plane in degrees (OUTPUT)
    %     dipdia         dip direction angle of first nodal plane in degrees (OUTPUT)
    %     strikb         strike angle of second nodal plane in degrees (OUTPUT)
    %     dipb           dip angle of second nodal plane in degrees (OUTPUT)
    %     rakeb          rake angle of second nodal plane in degrees (OUTPUT)
    %     dipdib         dip direction angle of second nodal plane in degrees (OUTPUT)
    %     ierr           error indicator (OUTPUT)
    %
    %     errors:
    %     1              input TREND angle of P axis out of range
    %     2              input PLUNGE angle P axis out of range
    %     3              1+2
    %     4              input TREND angle of P axis out of range
    %     5              input PLUNGE angle P axis out of range
    %     6              4+5
    %     8,9,10         internal errors
    %
    % c
    %      call fpsset
    amistr=-360.;
    amastr=360.;
    amidip=0.;
    amadip=90.;
    amirak=-360.;
    amarak=360.;
    amitre=-360.;
    amatre=360.;
    amiplu=0.;
    amaplu=90.;
    orttol=2.;
    ovrtol=0.001;
    tentol=0.0001;
    dtor=0.017453292519943296;
    c360=360.;
    c90=90.;
    c0=0.;
    c1=1.;
    c2=2.;
    c3=3.;
    % c
    %      call ax2ca(trendp,plungp,px,py,pz,ierr)
    [px,py,pz,ierr] = focal_ax2ca(trendp,plungp);
    if (ierr ~= 0)
        disp(['PT2PL: ierr=' num2str(ierr)]);
        return;
    end

    %      call ax2ca(trendt,plungt,tx,ty,tz,ierr)
    [tx,ty,tz,ierr] = focal_ax2ca(trendt,plungt);
    if (ierr ~= 0)
        ierr = ierr + 3;
        disp(['PT2PL: ierr=' num2str(ierr)]);
        return;
    end

    [anx,any,anz,dx,dy,dz,ierr] = focal_pt2nd(px,py,pz,tx,ty,tz);
    if (ierr ~= 0)
        ierr = 8;
        disp(['PT2PL: ierr=' num2str(ierr)]);
        return;
    end
    [strika,dipa,rakea,dipdia,ierr] = focal_nd2pl(anx,any,anz,dx,dy,dz);
    if (ierr ~= 0)
        ierr = 9;
        disp(['PT2PL: ierr=' num2str(ierr)]);
        return;
    end
    [strikb,dipb,rakeb,dipdib,ierr] = focal_nd2pl(dx,dy,dz,anx,any,anz);
    if (ierr ~= 0)
        ierr = 10;
        disp(['PT2PL: ierr=' num2str(ierr)]);
        return;
    end
end

author: CelsoReyes


text:
function c = tapas_softmax_mu3_wld_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the softmax observation model for multinomial responses with phasic
% volatility exp(mu3) as the decision temperature and parameters accounting for win- and
% loss-distortion of state values.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2018 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Is the decision based on predictions or posteriors? Comment as appropriate.
c.predorpost = 1; % Predictions
%c.predorpost = 2; % Posteriors

% Model name
c.model = 'softmax_mu3_wld';

% Sufficient statistics of Gaussian parameter priors

% Win-distortion
c.la_wdmu = 0;
c.la_wdsa = 2^-2;

% Loss-distortion
c.la_ldmu = 0;
c.la_ldsa = 2^-2;

% Gather prior settings in vectors
c.priormus = [
    c.la_wdmu,...
    c.la_ldmu,...
         ];

c.priorsas = [
    c.la_wdsa,...
    c.la_ldsa,...
         ];

% Model filehandle
c.obs_fun = @tapas_softmax_mu3_wld;

% Handle to function that transforms observation parameters to their native space
% from the space they are estimated in
c.transp_obs_fun = @tapas_softmax_mu3_wld_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
% SATGLOBE - Draw an idealized satellite view of earth
%
% This file renders a fully manipulatable satellite view
% of earth at a resolution of two pixels per degree, with added
% international political boundaries and gridlines.
% The imagery was obtained from NASA, and the globe was
% rendered using the Matlab Mapping Toolbox.
%
% The Mapping Toolbox is not needed to use this file.
%
% In order to save storage space, this m-file loads image
% data from the file satglobe.mat, and then creates the
% graticule mesh itself. This process allows users who
% do not have the Matlab Mapping Toolbox to render the
% figure, but it does take a few moments to compute the
% mesh. Using this trick, the data storage is reduced
% considerably; however, once the figure is redered, you
% may wish to save it as a regular Matlab figure file
% to increase speed.
%
% Michael Kleder, 2004

function satglobe
load satglobe
th=repmat((0:.5:180)'*pi/180,[1 721]);
ph=repmat((-180:.5:180)*pi/180,[361 1]);
s.children(1).properties.XData = sin(th).*cos(ph);
s.children(1).properties.YData = sin(th).*sin(ph);
s.children(1).properties.ZData = cos(th);
s.children(1).properties.CData = double(c)/255;
figure;
struct2handle(s,gcf);
set(gcf,'color','k','renderer','zbuffer','inverthardcopy','off','name',...
   'Earth at 2 Pixels per Degree, by Michael Kleder');
return


text:
function [X,P] = femUnk(fe)
%+========================================================================+
%|                                                                        |
%|              OPENFEM - LIBRARY FOR FINITE ELEMENT METHOD               |
%|           openFem is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal & Francois Alouges (c) 2017-2018.          |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%|             francois.alouges@polytechnique.edu                         |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : femUnk.m                                      |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal & François Alouges            |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Unknowns and reduction matrix for constrained |
%|  `---'  |                finite elements                               |
%+========================================================================+

% Initialization with dofs
X  = fe.dof;
Nx = size(X,1);
P  = speye(Nx);

% Dirichlet
if ~isempty(fe.dir)
    % Finite element for dirichlet mesh
    fed = fem(fe.dir,fe.typ);

    % Dof indices
    I = find(~ismember(X,fed.dof,'rows'));

    % Reduction matrix
    M = speye(Nx);
    M = M(:,I);

    % Update
    P  = P * M;
    X  = X(I,:);
    Nx = length(I);
end

% Junction
if ~isempty(fe.jct)
    % Number of junction
    Njct = length(fe.jct)/2;

    % Initialization
    feJct = cell(1,Njct);
    Ijct  = cell(1,Njct);

    % Finite element for junction meshes
    for i = 1:Njct
        if isa(fe.jct{2*i-1},'msh') && isnumeric(fe.jct{2*i})
            feJct{i} = fem(fe.jct{2*i-1},fe.typ);
        else
            error('femUnknown.m : unavailable case');
        end
    end

    % Valid dof indices for junction meshes
    for i = 1:Njct
        Xi         = feJct{i}.dof;
        Ijct{i}    = zeros(size(Xi,1),1);
        [~,I,J]    = intersect(Xi,X,'rows','stable');
        Ijct{i}(I) = J;
    end
    Ijct = cell2mat(Ijct);

    % Extract dirichlet condition and multiple indices
    bool = (Ijct(:,Njct) == 0);
    for i = 1:Njct-1
        bool = bool + (Ijct(:,Njct) == Ijct(:,i));
    end
    Ijct = Ijct(~bool,:);

    % Final dof after extraction of the last indices
    I = setdiff((1:Nx)',Ijct(:,end));

    % Linear relation matrix
    M = speye(Nx);
    for  i = 1:Njct
        for j = setdiff(1:Njct,i)
            M = M + sparse(Ijct(:,i),Ijct(:,j),-fe.jct{2*j}/fe.jct{2*i},Nx,Nx);
        end
    end

    % Reduction
    M = M(:,I);

    % Update
    P  = P * M;
    X  = X(I,:);
%     Nx = length(I);
end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%% Demo: White balancing a single image
%
% Copyright (c) 2018-present, Mahmoud Afifi
% York University, Canada
% mafifi@eecs.yorku.ca | m.3afifi@gmail.com
%
% This source code is licensed under the license found in the
% LICENSE file in the root directory of this source tree.
% All rights reserved.
%
% Please cite the following work if this program is used:
% Mahmoud Afifi, Brian Price, Scott Cohen, and Michael S. Brown,
% "When color constancy goes wrong: Correcting improperly white-balanced
% images", CVPR 2019.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%

%% input and options
infileName = fullfile('..','example_images','figure3.jpg');
outfileName = fullfile('result.jpg');
device = 'cpu'; % 'cpu' or 'gpu'
gamut_mapping = 2; % use 1 for scaling, 2 for clipping (our paper's results reported using clipping).
upgraded_model = 1; % use 1 to load our new model that is upgraded with new training examples.

%%
switch lower(device)
    case 'cpu'
        if upgraded_model == 1
            load(fullfile('models','WB_model+.mat'));
        elseif upgraded_model == 0
            load(fullfile('models','WB_model.mat'));
        else
            error('Wrong upgraded_model value; please use 0 or 1');
        end
    case 'gpu'
        try
            gpuDevice();
        catch
            error('Cannot find a GPU device');
        end
        if upgraded_model == 1
            load(fullfile('models','WB_model+_gpu.mat'));
        elseif upgraded_model == 0
            load(fullfile('models','WB_model_gpu.mat'));
        else
            error('Wrong upgraded_model value; please use 0 or 1');
        end
    otherwise
        error('Wrong device; please use ''gpu'' or ''cpu''')
end
model.gamut_mapping = gamut_mapping;
fprintf('Processing image: %s\n',infileName);
I_in = imread(infileName);
tic
I_corr = model.correctImage(I_in);
disp('Done!');
toc
subplot(1,2,1); imshow(I_in); title('Input');
subplot(1,2,2); imshow(I_corr); title('Our result');
disp('Saving...');
if strcmpi(device,'cpu')
    imwrite(I_corr,outfileName);
else
    imwrite(gather(I_corr),outfileName);
end
disp('Saved!');

author: mahmoudnafifi


text:
function dataNdCombined = combine(this, dataNdArray, combineDims, tolerance)
% combines multiple n-dim. datasets into a single one along specified
% dimensions. Makes sure data is sorted into right place according to
% different dimInfos
%
% NOTE: inverse operation of MrDataNd.split
%
%   Y = MrDataNd()
%   dataNdCombined = Y.combine(dataNdArray, combineDims, tolerance)
%
% This is a method of class MrDataNd.
%
% IN
%   dataNdArray     cell(nDatasets,1) of MrDataNd to be combined
%   combineDims     [1, nCombineDims] vector of dim indices to be combined
%                       OR
%                   cell(1, nCombineDims) of dimLabels to be combined
%                   NOTE: If specified dimLabels do not exist, new
%                   dimensions are created with these names and default
%                   samplingPoints (1:nDatasets)
%                   default: all singleton dimensions (i.e. dims with one
%                   sample only within each individual dimInfo)
%                   NOTE: if a non-singleton dimension is given, images are
%                         concatenated along this dimension
%
%   tolerance                   dimInfos are only combined, if their
%                               information is equal for all but the
%                               combineDims (because only one
%                               representation is retained for those,
%                               usually from the first of the dimInfos).
%                               However, sometimes numerical precision,
%                               e.g., rounding errors, preclude the
%                               combination. Then you can increase this
%                               tolerance;
%                               default: single precision (eps('single')
%                               ~1.2e-7)
% OUT
%
% EXAMPLE
%   combine
%
%   See also MrDataNd MrDimInfo.combine MrDataNd.split

% Author:   Lars Kasper
% Created:  2018-05-16
% Copyright (C) 2018 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

if nargin < 4
    tolerance = eps('single');
end

% create new
dataNdCombined = this.copyobj();

%% 1) dimInfoCombined = Y.combine(dimInfoArray, combineDims)
doCombineSingletonDims = nargin < 3;
if doCombineSingletonDims
    indSplitDims = this.dimInfo.get_singleton_dimensions();
    combineDims = this.dimInfo.dimLabels(indSplitDims);
else
    % for 1-dim case, make cell
    if ~iscell(combineDims) && isstr(combineDims)
        combineDims = {combineDims};
    end
end

%% Create dimInfoArray from all data and combine it first
dimInfoArray = cellfun(@(x) x.dimInfo.copyobj(), dataNdArray, 'UniformOutput', false);

% if dimLabels not existing previously, add as new dimensions
for iDim = 1:numel(combineDims)
    combineDim = combineDims{iDim};
    % new dimension that did not exist in dimInfo
    if isempty(this.dimInfo.get_dim_index(combineDim))
        cellfun(@(x,y) x.add_dims(combineDim, 'units', 'nil', ...
            'samplingPoints', y), dimInfoArray, ...
            num2cell(1:size(dataNdArray, iDim))', ...
            'UniformOutput', false);
    end
end

[dimInfoCombined, indSamplingPointCombined] = dimInfoArray{1}.combine(...
    dimInfoArray, combineDims, tolerance);

%% Loop over all splits dataNd and put data into right place, as defined by combined DimInfo
if ~isempty(combineDims) % otherwise, nothing to do hear
    % dimInfo sampling points
    indSplitDims        = dimInfoArray{1}.get_dim_index(combineDims);
    nSplits             = numel(dataNdArray);
    dataMatrixCombined  = nan(dimInfoCombined.nSamples);
    for iSplit = 1:nSplits
        % write out indices to be filled in final array, e.g. tempData(:,:,sli, dyn)
        % would be {':', ':', sli, dyn}
        index = repmat({':'}, 1, dimInfoCombined.nDims);
        index(indSplitDims) = indSamplingPointCombined(iSplit,:);
        dataMatrixCombined(index{:}) = dataNdArray{iSplit}.data;
    end

    %% assemble the output object
    dataNdCombined.dimInfo = dimInfoCombined;
    dataNdCombined.data = dataMatrixCombined;
end
end


text:
function range = p04_range ( )

%*****************************************************************************80
%
%% P04_RANGE returns an interval bounding the root for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real RANGE(2), the minimum and maximum values of
%    an interval containing the root.
%
  range(1) =  0.00001;
  range(2) = 20.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Param, extra] = parsepv(Param, pvpairs, varargin)
%PARSEPV Parses parameter/value pairs
%
% NewParam = parsepv(Param, pvpairs)
% [NewParam, extra] = parsepv(Param, pvpairs, 'returnextra')
%
% This function is an extension of parse_pv_pairs.  It allows the option of
% returning unrecognized parameter/value pairs, rather than erroring.
%
% Input variables:
%
%   Param:          1 x 1 structure holding default parameters (fieldnames)
%                   and values
%
%   pvpairs:        1 x n cell array of parameter/value pairs
%
%   'returnextra':  if this string is included, the function will return a
%                   cell array holding any unrecognized parameters and the
%                   corresponding values.  Otherwise, it will error if a
%                   parameter is not recognized.
%
% Output variables:
%
%   NewParam:       1 x 1 struct identical to Param but with defaults
%                   replaced by the values from pvpairs
%
%   extra:          1 x m cell array of any unrecognized parameter/value
%                   pairs

% Copyright 2009 Kelly Kearney

if nargin == 3
    returnextra = strcmpi(varargin{1}, 'returnextra');
else
    returnextra = false;
end

npv = length(pvpairs);
n = npv/2;

if n~=floor(n)
  error 'Property/value pairs must come in PAIRS.'
end
if n<=0
  % just return the defaults
  if returnextra
      extra = cell(0);
  end
  return
end

if ~isstruct(Param)
    error 'No structure for defaults was supplied'
end

if returnextra
    extra = cell(0);
end

% there was at least one pv pair. process any supplied
propnames = fieldnames(Param);
lpropnames = lower(propnames);
for i=1:n
    p_i = lower(pvpairs{2*i-1});
    v_i = pvpairs{2*i};

    ind = strmatch(p_i,lpropnames,'exact');
    if isempty(ind)
        ind = find(strncmp(p_i,lpropnames,length(p_i)));
        if isempty(ind)
            if returnextra
                extra = [extra pvpairs(2*i-1:2*i)];
                continue;
            else
                error(['No matching property found for: ',pvpairs{2*i-1}]);
            end
        elseif length(ind)>1
            error(['Ambiguous property name: ',pvpairs{2*i-1}])
        end
    end
    p_i = propnames{ind};

    % override the corresponding default in params
    Param = setfield(Param,p_i,v_i);

end


text:
function [K,dK] = gpcov_ratquad(D, p1, p2, p3)
% f = 1 + D.^2/(2*p3^2*p2^2);
% K = p1^2*f.^(-p3^2);

f = 1 + D.^2/(2*p3^2*p2^2);
K = p1^2*f.^(-p3^2);

if nargout >= 2
  dK = zeros([size(D),3]);
  dK(:,:,1) = K .* 2 / p1;
  dK(:,:,2) = K .* (-p3^2).*f.^(-1) .* (-2).*D.^2/(2*p3^2)*p2^(-3);
  dK(:,:,3) = K .* (-2*p3*log(f) + (-p3^2)./f.*D.^2/(2*p2^2) * (-2) * p3^(-3));
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function psnr=compute_psnr(im1,im2)
if size(im1, 3) == 3,
    im1 = rgb2ycbcr(im1);
    im1 = im1(:, :, 1);
end

if size(im2, 3) == 3,
    im2 = rgb2ycbcr(im2);
    im2 = im2(:, :, 1);
end

imdff = double(im1) - double(im2);
imdff = imdff(:);

rmse = sqrt(mean(imdff.^2));
psnr = 20*log10(255/rmse);

author: tyshiwo


text:
function [DIP, teta, hr]= DipoleGeneratorAbnormal(N,fs,rr,alphai,bi,tetai,teta0)
%
% [DIP teta]= DipoleGeneratorAbnormal(N,fs,rr,alphai,bi,tetai,teta0)
% Synthetic cardiac dipole generator using the 'differential form' of the
% dipole equations. Refer to references of the toolbox for further details.
%
% inputs:
% N: signal length
% fs: sampling rate
% rr: rr interval time series
% alphai: structure contaning the amplitudes of Gaussian functions used for
%       modeling the x, y, and z coordinates of the cardiac dipole
% bi: structure contaning the widths of Gaussian functions used for
%       modeling the x, y, and z coordinates of the cardiac dipole
% tetai: structure contaning the phase of Gaussian functions used for
%       modeling the x, y, and z coordinates of the cardiac dipole
% teta0: initial phase of the synthetic dipole

%
% output:
% DIP: structure contaning the x, y, and z coordinates of the cardiac dipole
% teta: vector containing the dipole phase
%
%
% Open Source ECG Toolbox, version 2.0, April 2008
% Released under the GNU General Public License
% Copyright (C) 2008  Reza Sameni
% Sharif University of Technology, Tehran, Iran -- GIPSA-Lab, Grenoble, France
% reza.sameni@gmail.com

% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details.

% Last update Shamim Nemati Feb 13 2009

rr = rr(:); t = [cumsum(rr)]; tp = linspace(0,sum(rr),N);
rrp = spline(t,rr,tp);
f=1./rrp; % make a heart rate (per second) time-series from rr intervals resampled to the current sampling frequency
hr = f*60;
w = 2*pi*f;     dt = 1/fs;
teta = zeros(1,N); X = zeros(1,N); Y = zeros(1,N); Z = zeros(1,N);

teta(1) = teta0;

for i = 1:N-1
    teta(i+1) = teta(i) + w(i)*dt; %dtheta/dt = w
    if(teta(i+1)>pi) % beat transition ------------------------------------
        teta(i+1) = teta(i+1) - 2*pi;
    end
    %------------------------------------------
    dtetaix = mod(teta(i) - tetai.x + pi , 2*pi) - pi;
    dtetaiy = mod(teta(i) - tetai.y + pi , 2*pi) - pi;
    dtetaiz = mod(teta(i) - tetai.z + pi , 2*pi) - pi;

    if(i==1),
        X(i) = sum(alphai.x .* exp(-dtetaix .^2 ./ (2*bi.x .^ 2)));
        Y(i) = sum(alphai.y .* exp(-dtetaiy .^2 ./ (2*bi.y .^ 2)));
        Z(i) = sum(alphai.z .* exp(-dtetaiz .^2 ./ (2*bi.z .^ 2)));
    end
    %------ Eq.(1) in Clifford et al. CinC2008 ------
    X(i+1) = X(i) - dt*sum(w(i)*alphai.x ./ (bi.x) .^ 2 .* dtetaix .* exp(-dtetaix .^2 ./ (2* bi.x .^ 2)));   % x state variable
    Y(i+1) = Y(i) - dt*sum(w(i)*alphai.y ./ (bi.y) .^ 2 .* dtetaiy .* exp(-dtetaiy .^2 ./ (2* bi.y .^ 2)));   % y state variable
    Z(i+1) = Z(i) - dt*sum(w(i)*alphai.z ./ (bi.z) .^ 2 .* dtetaiz .* exp(-dtetaiz .^2 ./ (2* bi.z .^ 2)));   % z state variable
    %------------------------------------------------

end
DIP.x = X;
DIP.y = Y;
DIP.z = Z;

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function PlotRegMean(ops1,ops)

figure('position', [900 50 900 900])
ax = ceil(sqrt(numel(ops1)));
for i = 1:length(ops1)
    subplot(ax,ax,i)
    imagesc(ops1{i}.mimg)
    colormap('gray')
    title(sprintf('Plane %d', i))
end
annotation('textbox', [0 .96 1 .03], 'String', ...
    sprintf('Registration for mouse %s, date %s', ops.mouse_name, ops.date), ...
    'Interpreter', 'none', 'FontWeight', 'bold', 'LineStyle', 'none', ...
    'HorizontalAlignment', 'center');
drawnow

author: cortex-lab


text:
function trans = createScaling3d(varargin)
%CREATESCALING3D Create the 4x4 matrix of a 3D scaling
%
%   TRANS = createScaling3d(S);
%   returns the scaling transform corresponding to a scaling factor S in
%   each direction. S can be a scalar, or a 1x3 vector containing the
%   scaling factor in each direction.
%
%   TRANS = createScaling3d(SX, SY, SZ);
%   returns the scaling transform corresponding to a different scaling
%   factor in each direction.
%
%   The returned matrix has the form :
%   [SX  0  0  0]
%   [ 0 SY  0  0]
%   [ 0  0 SZ  0]
%   [ 0  0  0  0]
%
%   See also:
%   transforms3d, transformPoint3d, transformVector3d, createTranslation3d,
%   createRotationOx, createRotationOy, createRotationOz
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 20/04/2006.
%

%   HISTORY
%   25/11/2008 rename from scale3d to scaling3d
%   30/04/2009 rename to createScaling3d

% process input parameters
if isempty(varargin)
    % assert uniform scaling in each direction
    sx = 1;
    sy = 1;
    sz = 1;
elseif length(varargin)==1
    % only one argument
    var = varargin{1};
    if length(var)==1
        % same scaling factor in each direction
        sx = var;
        sy = var;
        sz = var;
    elseif length(var)==3
        % scaling is a vector, giving different scaling in each direction
        sx = var(1);
        sy = var(2);
        sz = var(3);
    else
        error('wrong size for first parameter of "createScaling3d"');
    end
elseif length(varargin)==3
    % 3 arguments, giving scaling in each direction
    sx = varargin{1};
    sy = varargin{2};
    sz = varargin{3};
else
    error('wrong number of arguments for "createScaling3d"');
end

% create the scaling matrix
trans = [sx 0 0 0;0 sy 0 0;0 0 sz 0;0 0 0 1];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function Dcross(sel)
    % Calculate the Dvalue for a volume around a grid points.
    %
    % see DCPARAIN.
    %
    %   Stefan Wiemer 1/95
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    global eq0p

    % the new data vector to be analysed is called Da, relative to the center of the x-section and already in km
    % D = [x,y,z ]
    catalog=ZG.primeCatalog; % points to same thing
    Da = [eq0p(1,:)' eq0p(2,:)' catalog.Date catalog.Date.Month catalog.Date.Day catalog.Magnitude catalog.Depth];
    Da = Da.subset(-2.99 < Da(:,7));
    if exist('sel','var')
        switch sel
            case 'ca'
                my_calculate();
            case 'lo'
                my_load();
        end
        % get the grid parameter
        % initial values
        %
        dd = 1.00;
        dx = 1.00 ;
        ni = 600;
        Nmin = 600; %on line 303 it has been replaced by ni
        stan2 = nan;
        stan = nan;
        prf = nan;
        av = nan;

        % make the interface
        %
        figure_w_normalized_uicontrolunits(...
            'Name','Grid Parameters',...
            'NumberTitle','off', ...
            'units','points',...
            'Visible','on', ...
            'Position',[ 100 200 500 200]);
        axis off
        % Francesco ...

        % creates a dialog box to input grid parameters
        %
        freq_field=uicontrol('Style','edit',...
            'Position',[.32 .57 .12 .08],...
            'Units','normalized','String',num2str(ni),...
            'callback',@callbackfun_001);

        freq_field0=uicontrol('Style','edit',...
            'Position',[.80 .57 .12 .08],...
            'Units','normalized','String',num2str(ra),...
            'callback',@callbackfun_002);

        freq_field2=uicontrol('Style','edit',...
            'Position',[.32 .44 .12 .08],...
            'Units','normalized','String',num2str(dx),...
            'callback',@callbackfun_003);

        freq_field3=uicontrol('Style','edit',...
            'Position',[.32 .31 .12 .08],...
            'Units','normalized','String',num2str(dd),...
            'callback',@callbackfun_004);

        tgl1 = uicontrol('Backgroundcolor', [0.8 0.8 0.8], 'Fontweight','bold',...
            'FontSize', 10, 'Style','checkbox',...
            'string','Number of Events:',...
            'Position',[.05 .56 .2 .10],...
            'callback',@callbackfun_005,...
            'Units','normalized');

        set(tgl1,'value',1);

        tgl2 =  uicontrol('BackGroundColor', [0.8 0.8 0.8],'Style','checkbox',...
            'string','Constant Radius:','Fontweight','bold','FontSize', 10,...
            'Position',[.55 .56 .2 .1],...
            'callback',@callbackfun_006,...
            'Units','normalized');

        set(tgl2, 'ForegroundColor', 'w');

        close_button=uicontrol('Style','Pushbutton',...
            'Position',[.45 .05 .15 .13 ],...
            'Units','normalized', 'Callback', @callbackfun_007,'String','Cancel');

        help_button=uicontrol('Style','Pushbutton',...
            'Position',[.70 .05 .15 .13 ],...
            'Units','normalized', 'Callback', @callbackfun_008,'String','Help');

        go_button1=uicontrol('Style','Pushbutton',...
            'Position',[.20 .05 .15 .13 ],...
            'Units','normalized',...
            'callback',@callbackfun_009,...
            'String','Go');

        txt3 = text(...
            'Position',[0.35 0.9 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.l ,...
            'FontWeight','bold',...
            'String',' Grid Parameters');
        txt5 = text(...
            'Position',[-0.07 0.46 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'String','Horizontal Spacing [km]:');

        txt6 = text(...
            'Position',[-0.07 0.30 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'String','Depth spacing [km]:');

        txt7 = text(...
            'Position',[0.45 0.62 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'Color', 'r',...
            'String','OR');

        set(gcf,'visible','on');
        watchoff

        % get the grid-size interactively and
        % calculate the b-value in the grid by sorting
        % the seismicity and selectiong the ni neighbors
        % to each grid point

        functin my_calculate()

        figure(xsec_fig);
        figure;
        ax = plot(Da(:,1),-Da(:,7),'o');
        xlabel('Distance in [km]')
        ylabel('Depth in [km]')

        set(gca,'NextPlot','add')
        ax=findobj(gcf,'Tag','mainmap_ax');
        [x,y, mouse_points_overlay] = select_polygon(ax);

        plos2 = plot(x,y,'b-', 'Color', 'r');        % plot outline
        sum3 = 0.;
        pause(0.3)

        %create a rectangular grid
        xvect=[min(x):dx:max(x)];
        yvect=[min(y):dd:max(y)];
        gx = xvect;gy = yvect;
        tmpgri=zeros((length(xvect)*length(yvect)),2);
        n=0;

        for i=1:length(xvect)
            for j=1:length(yvect)
                n=n+1;
                tmpgri(n,:)=[xvect(i) yvect(j)];
            end
        end

        %extract all gridpoints in chosen polygon
        XI=tmpgri(:,1);
        YI=tmpgri(:,2);

        ll = polygon_filter(x,y, XI, YI, 'inside');
        %grid points in polygon
        newgri=tmpgri(ll,:);

        % Plot all grid points
        plot(newgri(:,1),newgri(:,2),'+k')

        if length(xvect) < 2  ||  length(yvect) < 2
            errordlg('Selection too small! (not a matrix)');
            return
        end

        itotal = length(newgri(:,1));
        if length(gx) < 4  ||  length(gy) < 4
            errordlg('Selection too small! ');
            return
        end

        %  make grid, calculate start- endtime etc.  ...
        %
        [t0b, teb] = bounds(newa.Date) ;
        n = newa.Count;
        tdiff = round((teb-t0b)/ZG.bin_dur);
        loc = zeros(3, length(gx)*length(gy));

        % loop over  all points
        %
        i2 = 0.;
        i1 = 0.;
        bvg = [];
        allcount = 0.;
        wai = waitbar(0,' Please Wait ...  ');
        set(wai,'NumberTitle','off','Name','D-value grid - percent done', 'position', [250 80 270 50]);
        drawnow;
        %

        %
        for i= 1:length(newgri(:,1))
            x = newgri(i,1);y = newgri(i,2);
            allcount = allcount + 1.;
            i2 = i2+1;

            % calculate distance from center point and sort wrt distance

            l = sqrt(((Da(:,1) - x)).^2 + ((Da(:,7) + y)).^2 + (Da(:,2).^2)) ;
            [s,is] = sort(l);
            b = Da(is(:,1),:) ;       % re-orders matrix to agree row-wise

            if tgl1 == 0   % take point within r
                l3 = l <= ra;
                b = Da.subset(l3);      % new data per grid point (b) is sorted in distanc
                rd = ra;
            else
                % take first ni points
                b = b(1:ni,:);      % new data per grid point (b) is sorted in distance
                l2 = sort(l);
                rd = l2(ni);

            end

            %estimate the completeness and b-value, and take the zero depth events away.
            %ZG.newt2 = [b.Longitude b.Latitude zeros(size(b,1),1) zeros(size(b,1),1) zeros(size(b,1),1) zeros(size(b,1),1) b.Date];
            %
            ZG.newt2  = b;
            E = ZG.newt2;

            if b.Count >= ni  % enough events?

                dtokm = 0;
                [bv magco0 stan av pr] =  bvalca3(b.Magnitude,McAutoEstimate.auto);

                if range == 1 | range == 2

                    pdc3nofig(E);

                elseif range == 3

                    pdc3;
                    pause;

                end %if range = 1|2

                D = coef(1,1);
                fdallfig;

            else
                D = nan;
                bv = nan;

            end %if length >= ni

            bvg = [bvg ; D x y rd bv deltar];
            waitbar(allcount/itotal)

        end  % for  newgri

        figure(HCIfig);
        cb = colorbar('horiz');
        set(cb, 'position', [0.32 0.08 0.4 0.03], 'XTickLabel', col);
        axes('pos',[0 0 1 1]); axis off; set(gca,'NextPlot','add');
        te= text('string','D-value','pos',[0.49,0.01], 'fontsize',8,'FontWeight','bold')
        set(gcf,'visible','on');
        % save data
        %
        %  set(txt1,'String', 'Saving data...')

        drawnow
        gx = xvect;gy = yvect;

        catsave3('Dcross');
        %corrected the window positioning error
        close(wai)
        watchoff

        %
        % reshape a few matrices
        %
        normlap2=nan(length(tmpgri(:,1)),1)

        reshaper=@(v)reshape(v,length(yvect),length(xvect));
        normlap2(ll)= bvg(:,1);
        valueMap=reshaper(normlap2);

        normlap2(ll)= bvg(:,4);
        reso = reshaper(normlap2);

        normlap2(ll)= bvg(:,5);
        BM=reshaper(normlap2);

        old = valueMap;

        % View the b-value map
        view_Dv(valueMap, lab1, Da)

    end

    % Load exist D-grid
    function my_load()
        [file1,path1] = uigetfile(['*.mat'],'b-value gridfile');
        if length(path1) > 1

            load([path1 file1])
            xsecx = newa(:,end)';
            xsecy = newa(:,7);
            xvect = gx; yvect = gy;
            tmpgri=zeros((length(xvect)*length(yvect)),2);

            normlap2=nan(length(tmpgri(:,1)),1)
            normlap2(ll)= bvg(:,1);
            valueMap=reshape(normlap2,length(yvect),length(xvect));

            normlap2(ll)= bvg(:,4);
            r=reshape(normlap2,length(yvect),length(xvect));

            old = valueMap;

            nlammap
            [xsecx xsecy,  inde] =mysect(catalog.Latitude',catalog.Longitude',catalog.Depth,ZG.xsec_defaults.WidthKm,0,lat1,lon1,lat2,lon2);
            % Plot all grid points
            set(gca,'NextPlot','add')
            plot(newgri(:,1),newgri(:,2),'+k')
            view_Dv(valueMap, lab1, Da)
        else
            return
        end
    end

    function callbackfun_001(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ni=str2double(freq_field.String);
        freq_field.String=num2str(ni);
        tgl2.Value=0;
        tgl1.Value=1;
    end

    function callbackfun_002(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ra=str2double(freq_field0.String);
        freq_field0.String=num2str(ra);
        tgl2.Value=1;
        tgl1.Value=0;
    end

    function callbackfun_003(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        dx=str2double(freq_field2.String);
        freq_field2.String=num2str(dx);
    end

    function callbackfun_004(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        dd=str2double(freq_field3.String);
        freq_field3.String=num2str(dd);
    end

    function callbackfun_005(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        set(tgl2,'Value',0, 'ForegroundColor', 'w');
        set(tgl1, 'ForegroundColor', 'k');
    end

    function callbackfun_006(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        set(tgl1,'Value',0,'ForegroundColor', 'w');
        set(tgl2, 'ForegroundColor', 'k');
    end

    function callbackfun_007(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        close;

    end

    function callbackfun_008(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        close;

    end

    function callbackfun_009(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        tgl1=tgl1.Value;
        tgl2=tgl2.Value;
        close;
        gobut = 3;
        org = 1;
        startfd(1, gobut);
    end

end


text:
function [ obblist, labellist ] = visualizeRelpos3_alignment( kids, mergereps, leafreps, params, labelset )
%VISUALIZERELPOS3_ALIGNMENT
% to visualize the absolute position of objects
% from the relative position respecting to its parents
% the layout needs to be adjust based on the "same orientation" and "alignment"
%
% the default setting is 2D, one-hot label vectors
%
% INPUT: kids, mergereps, leafreps, params   - the data to be visualized
%        labelset  - the name of the categories (order matters)
%
% OUTPUT: obblist  - the 3d obblist of all the leaves
%         labellist  - labels of all the leaves

dim = 2;
BOX = 0;
SUPPORT = 1;
GROUP = 2;
SURROUND = 3;

LABELLEN = length(labelset);
labeldefs = zeros(LABELLEN, LABELLEN);
for i = 1:LABELLEN
    labeldefs(i,i) = 1;
end

% transform the layers to discrete values
layers = [0.2,0.4,0.6,0.8,1];
for i = 1:size(leafreps,2)
    dist = abs(layers - leafreps(3,i));
    [~,I] = min(dist);
    leafreps(3,i) = layers(I)*5;
end

% transform the labels to one-hot labels
labellist = {};
for i = 1:size(leafreps,2)
    label = leafreps(4:end,i);
    [~,I] = max(label);
    leafreps(4:end,i) = labeldefs(:,I);
    labellist{length(labellist)+1} = labelset{I};
end

% back propagate the positions
rcent = [0;0]; % root center
rfront = [1;0]; % root front
framelist = zeros(4,length(kids));
framelist(:,length(kids)) = [rcent;rfront];
for i = length(kids):-1:1
    k = kids{i};
    frame = framelist(:,i);
    pcent = frame(1:2);
    pfront = frame(3:4);
    paxes = [-pfront(2);pfront(1)];
    transform = [pfront(:)';paxes(:)'];
%     transform = inv(transform);

    rps = mergereps{i};
    for j = 2:length(k)
        ind = k(j);
        rp = rps(:,j-1);

        % adjust the orientation
        sameori = rp(end);
        gtori = [-1,0,1];
        dist = abs(sameori - gtori);
        [~,I] = min(dist);
        sameori = gtori(I);
        if(sameori==1)
            rp(3:4) = [1;0];
        end
        if(sameori==-1)
            rp(3:4) = [-1;0];
        end

        ccent = rp(1:2);
        cfront = rp(3:4);
        cfront = cfront/norm(cfront,2);
        ncent = transform\ccent;
        ncent = ncent + pcent;
        nfront = transform\cfront;
        nfront = nfront/norm(nfront,2);
        nframe = [ncent;nfront];
        framelist(:,ind) = nframe;
    end
end

obblist = [];
for i = 1:length(kids)
    k = kids{i};
    if(k(1)==0)
        frame = framelist(:,i);
        leaf = leafreps(:,i);
        obb = [frame(1);0;frame(2);frame(3);0;frame(4);0;1;0;leaf(1);0;leaf(2)];
        obblist = [obblist,obb];
    end
end

%% forward propagate to adjust the alignments
translationlist = cell(length(kids),1);
nobblist = zeros(size(obblist,1),length(kids));
gtparam = [0,0.25,0.5,0.75,1];
for i = 1:length(kids)
    k = kids{i};
    flag = k(1);
    mparam = params{i};
    for j = 1:length(mparam)
       dist = abs(mparam(j)- gtparam);
       [~,I] = min(dist);
       mparam(j) = gtparam(I)*4;
    end
    if(flag==BOX)
        translationlist{i} = [0;0;0];
        obb = obblist(:,i);
        nobblist(:,i) = obb;
    end
    if(flag==SUPPORT)
        translationlist{i} = [0;0;0];
        list = k(2:end);
        nobblist(:,i) = updateOBB(nobblist,list);
    end
    if(flag==GROUP)
        if(all(mparam>0))
            obb1 = nobblist(:,k(2));
            obb2 = nobblist(:,k(3));
            [trans1, trans2, pobb] = adjustGroupOBB(obb1,obb2,mparam);
            translationlist{k(2)} = translationlist{k(2)} + trans1;
            translationlist{k(3)} = translationlist{k(3)} + trans2;
            translationlist{i} = [0;0;0];
            nobblist(:,i) = pobb;
        else
            translationlist{i} = [0;0;0];
            list = k(2:end);
            nobblist(:,i) = updateOBB(nobblist,list);
        end
    end
    if(flag==SURROUND)
        if(all(mparam>0))
            obb1 = nobblist(:,k(2));
            obb2 = nobblist(:,k(3));
            obb3 = nobblist(:,k(4));
            [trans1, trans2, trans3, pobb] = adjustSurroundOBB(obb1,obb2,obb3,mparam);
            translationlist{i} = [0;0;0];
            translationlist{k(2)} = translationlist{k(2)} + trans1;
            translationlist{k(3)} = translationlist{k(3)} + trans2;
            translationlist{k(4)} = translationlist{k(4)} + trans3;
            nobblist(:,i) = pobb;
        else
            translationlist{i} = [0;0;0];
            list = k(2:end);
            nobblist(:,i) = updateOBB(nobblist,list);
        end
    end
end

for i = length(kids):-1:1
    trans = translationlist{i};
    k = kids{i};
    list = k(2:end);
    for j = 1:length(list)
        ctrans = translationlist{list(j)};
        ctrans = ctrans+trans;
        translationlist{list(j)} = ctrans;
    end
end

for i = 1:size(obblist,2)
    obb = obblist(:,i);
    trans = translationlist{i};
    obb(1) = obb(1)+trans(1);
    obb(3) = obb(3)+trans(3);
    obblist(:,i) = obb;
end

end

repo: GRAINS


text:
%ITHIN Morphological skeletonization
%
% OUT = ITHIN(IM) is the binary skeleton of the binary image IM.  Any non-zero
% region is replaced by a network of single-pixel wide lines.
%
% OUT = ITHIN(IM,DELAY) as above but graphically displays each iteration
% of the skeletonization algorithm with a pause of DELAY seconds between
% each iteration.
%
% References::
%  - Robotics, Vision & Control, Section 12.5.3,
%    P. Corke, Springer 2011.
%
% See also HITORMISS, ITRIPLEPOINT, IENDPOINT.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function out = ithin(im, delay)

    % create a binary image
    im = im > 0;

    o = im;

    Sa = [0 0 0; NaN 1 NaN; 1 1 1];
    Sb = [NaN 0 0; 1 1 0; NaN 1 NaN];

    o = im;
    while true
        for i=1:4
            r = hitormiss(im, Sa);
            im = im - r;
            r = hitormiss(im, Sb);
            im = im - r;
            Sa = rot90(Sa);
            Sb = rot90(Sb);
        end
        if nargin > 1
            idisp(im);
            pause(delay);
        end
        if all(o == im)
            break;
        end
        o = im;
    end
    o = im;
    if nargout > 0
        out = o;
    end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function G = rayxgridfRCC(F, transla, val)
%------------------------------------------------------------------------------
% Multiplies gridfunction F with element of stencil given by its value
% val and its position transla.
% Excess area is filled by Cell-Centered (CC) Reflection across boundaries.
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: June 23, 2000.
% (c) 1999-2002 Stichting CWI, Amsterdam.
%------------------------------------------------------------------------------
o=[0 0];
if ~all(size(o) == size(transla))
  error(' rayxgridfRCC - unexpected dimensions of transla ')
end
G = val * moveLRRCC(moveUDRCC(F, -transla(1)), -transla(2));
%------------------------------------------------------------------------------

author: Sable


text:
function [ xnew, ferr, q ] = roots_rc ( n, x, fx, q )

%*****************************************************************************80
%
%% ROOTS_RC solves a system of nonlinear equations using reverse communication.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 January 2013
%
%  Author:
%
%    Original FORTRAN77 version by Gaston Gonnet.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Gaston Gonnet,
%    On the Structure of Zero Finders,
%    BIT Numerical Mathematics,
%    Volume 17, Number 2, June 1977, pages 170-183.
%
%  Parameters:
%
%    Input, integer N, the number of equations.
%
%    Input, real X(N).  Before the first call, the user should
%    set X to an initial guess or estimate for the root.  Thereafter, the input
%    value of X should be the output value of XNEW from the previous call.
%
%    Input, real FX(N), the value of the function at XNEW.
%
%    Workspace, real ( kind = 8 ) Q(2*N+2,N+2).  Before the first call
%    for a given problem, the user must set Q(2*N+1,1) to 0.0.
%
%    Output, real ( kind = 8 ) XNEW(N), a new point at which a function
%    value is requested.
%
%    Output, real ( kind = 8 ) FERR, the function error, that is, the sum of
%    the absolute values of the most recently computed function vector.
%
%    Workspace, real ( kind = 8 ) Q(2*N+2,N+2).
%
  ferr = sum ( abs ( fx(1:n) ) );
%
%  Initialization if Q(2*N+1,1) = 0.0.
%
  if ( q(2*n+1,1) == 0.0 )

    for i = 1 : n
      for j = 1 : n + 1
        q(i,j) = 0.0;
        q(i+1,j) = 0.0;
      end
      q(i,i) = 100.0;
      q(i+n,i) = 1.0;
    end

    q(2*n+1,1:n) = 1.0E+30;
    q(2*n+2,1:n) = n;

    for i = 1 : n
      q(i+n,n+1) = x(i);
    end

    q(1:n,n+1) = fx(1:n);

    q(2*n+1,n+1) = ferr;
    q(2*n+2,n+1) = 0.0;
    damp = 0.99;

  else

    jsus = 1;
    for i = 2 : n + 1
      if ( 2 * n <= q(2*n+2,i) )
        q(2*n+1,i) = 1.0E+30;
      end
      if ( q(2*n+2,jsus) < floor ( ( n + 3 ) / 2 ) )
        jsus = i;
      end
      if ( ( n + 3 ) / 2 <= q(2*n+2,i) && q(2*n+1,jsus) < q(2*n+1,i) )
        jsus = i;
      end
    end

    for i = 1 : n
      q(i+n,jsus) = x(i);
      q(i,jsus) = fx(i);
    end

    q(2*n+1,jsus) = ferr;
    q(2*n+2,jsus) = 0;
    jsma = 1;
    damp = 0.0;

    for j = 1 : n + 1
      if ( 1.0E+30 / 10.0 < q(2*n+1,j) )
        damp = 0.99;
      end
      if ( q(2*n+1,j) < q(2*n+1,jsma) )
        jsma = j;
      end
    end

    if ( jsma ~= n + 1 )
      for i = 1 : 2 * n + 2
        t = q(i,jsma);
        q(i,jsma) = q(i,n+1);
        q(i,n+1) = t;
      end
    end

  end

  q(1:n,n+2) = q(1:n,n+1);
%
%  Call the linear equation solver, which should not destroy the matrix
%  in Q(1:N,1:N), and should overwrite the solution into Q(1:N,N+2).
%
  q(1:n,n+2) = q(1:n,1:n) \ q(1:n,n+2);

  sump = sum ( q(1:n,n+2) );

  if ( abs ( 1.0 - sump ) <= 1.0E-10 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'ROOT - Fatal error!\n' );
    fprintf ( 1, '  SUMP almost exactly 1.\n' );
    fprintf ( 1, '  SUMP = %g\n', sump );
    error ( 'ROOTS_RC - Fatal error!' );
  end

  for i = 1 : n
    xnew(i) = q(i+n,n+1);
    for j = 1 : n
      xnew(i) = xnew(i) - q(i+n,j) * q(j,n+2);
    end
%
%  If system not complete, damp the solution.
%
    xnew(i) = xnew(i) / ( 1.0 - sump ) * ( 1.0 - damp ) ...
      + q(i+n,n+1) * damp;

  end

  for j = 1 : n + 1
    q(2*n+2,j) = q(2*n+2,j) + 1.0;
  end

  return
end


text:
function p = quantize(p,R,minV,maxV,minS,maxS)
%
% p = quantize(p,R,type) -- "quantize" elements of KDE p to R bits
%

%p.centers   = round(1000*p.centers)/1000;
%p.means     = round(1000*p.means)/1000;
%p.ranges    = round(1000*p.ranges)/1000;
%p.bandwidth = ceil(10000*p.bandwidth)/10000;

p.centers   = roundVals(p.centers  ,minV,maxV,R);
p.means     = roundVals(p.means    ,minV,maxV,R);
p.ranges    = roundVals(p.ranges   ,minV,maxV,R);
p.bandwidth = roundVals(p.bandwidth,minS,maxS,R);

function x = roundVals(x,minV,maxV,R)
 scale = 2^R ./ repmat(maxV-minV,size(x)./size(maxV));
 minV = repmat(minV,size(x)./size(minV)); x = x - minV; x = max(x,0);
 x = x .* scale; x = min(round(x),2^R-1); x = x./scale; x = x + minV;

author: ShapeNet

repo: RenderForCNN


text:
function cMap = imfConsistency(mu, refIdx, consistencyThres)

[s1, s2, s3] = size(mu);
cMap = zeros(s1, s2, s3);
cMap(:,:,refIdx) = ones(s1, s2);

refMu = mu(:,:,refIdx);
N = 256;
for i = 1 : s3
    if i ~= refIdx
          cMu  = imhistmatch(mu(:,:,i), refMu, N);
          diff = abs(cMu - refMu);
          cMap(:,:,i) = diff <= consistencyThres;
    end
end

author: thfylsty


text:
%% Copyright (C) 2014, 2016 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
%% @deftypemethod  @@symfun {@var{d} =} size (@var{x})
%% @deftypemethodx @@symfun {[@var{n}, @var{m}] =} size (@var{x})
%% @deftypemethodx @@symfun {@var{d} =} size (@var{x}, @var{dim})
%% Return the size of a symbolic function.
%%
%% This behaves differently than for @@sym:
%% @example
%% @group
%% syms x y
%% f(x, y) = [1 x; y 2];
%%
%% size(f)
%%   @result{} 1  1
%% @end group
%%
%% @group
%% length(f)
%%   @result{} 1
%% @end group
%% @end example
%%
%% @seealso{@@sym/size}
%% @end deftypemethod

function [n, m] = size(x, dim)

  n = [1 1];
  if (nargin == 2) && (nargout == 2)
    print_usage ();
  elseif (nargout == 2)
    m = 1;
    n = 1;
  elseif (nargin == 2)
    n = 1;
  end

end

%!test
%! syms x
%! f(x) = x;
%! d = size(f);
%! assert (isequal (d, [1 1]))
%! [n, m] = size(f);
%! assert (isequal ([n m], [1 1]))
%! assert (size(f, 1) == 1)
%! assert (size(f, 2) == 1)

%!test
%! syms x
%! f(x) = [1 x];
%! d = size(f);
%! assert (isequal (d, [1 1]))
%! [n, m] = size(f);
%! assert (isequal ([n m], [1 1]))
%! assert (size(f, 1) == 1)
%! assert (size(f, 2) == 1)

repo: octsympy