text:
function [net, info] = cnn_train_adagrad(net, imdb, getBatch, varargin)
% CNN_TRAIN   Demonstrates training a CNN
%    CNN_TRAIN() is an example learner implementing stochastic gradient
%    descent with momentum to train a CNN for image classification.
%    It can be used with different datasets by providing a suitable
%    getBatch function.

opts.train = [] ;
opts.val = [] ;
opts.numEpochs = 300 ;
opts.batchSize = 256 ;
opts.useGpu = false ;
opts.learningRate = 0.001 ;
opts.continue = false ;
opts.expDir = fullfile('data','exp') ;
opts.conserveMemory = false ;
opts.sync = true ;
opts.prefetch = false ;
opts.weightDecay = 0.0005 ;
opts.errorType = 'multiclass' ;
opts.plotDiagnostics = false ;
opts.delta = 1e-8;
opts.display = 1;
opts.snapshot = 1;
opts.test_interval = 1;
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir); end
if isempty(opts.train), opts.train = find(imdb.images.set==1); end
if isempty(opts.val), opts.val = find(imdb.images.set==2); end
if isnan(opts.train), opts.train = []; end

% -------------------------------------------------------------------------
%                                                    Network initialization
% -------------------------------------------------------------------------

for i=1:numel(net.layers)
  if ~strcmp(net.layers{i}.type,'conv'), continue; end
  net.layers{i}.filtersMomentum = zeros(size(net.layers{i}.filters), ...
    class(net.layers{i}.filters)) ;
  net.layers{i}.biasesMomentum = zeros(size(net.layers{i}.biases), ...
    class(net.layers{i}.biases)) ; %#ok<*ZEROLIKE>
  if ~isfield(net.layers{i}, 'filtersLearningRate')
    net.layers{i}.filtersLearningRate = 1 ;
  end
  if ~isfield(net.layers{i}, 'biasesLearningRate')
    net.layers{i}.biasesLearningRate = 1 ;
  end
  if ~isfield(net.layers{i}, 'filtersWeightDecay')
    net.layers{i}.filtersWeightDecay = 1 ;
  end
  if ~isfield(net.layers{i}, 'biasesWeightDecay')
    net.layers{i}.biasesWeightDecay = 1 ;
  end
end

if opts.useGpu
  net = vl_simplenn_move(net, 'gpu') ;
  for i=1:numel(net.layers)
    if ~strcmp(net.layers{i}.type,'conv'), continue; end
    net.layers{i}.filtersMomentum = gpuArray(net.layers{i}.filtersMomentum) ;
    net.layers{i}.biasesMomentum = gpuArray(net.layers{i}.biasesMomentum) ;
  end
end

G_f = cell(numel(net.layers), 1);
G_b = cell(numel(net.layers), 1);

for l=1:numel(net.layers)

  if ~strcmp(net.layers{l}.type, 'conv'), continue ; end

  G_f{l} = zeros(size(net.layers{l}.filters), 'single');
  G_b{l} = zeros(size(net.layers{l}.biases), 'single');

end

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

rng(0) ;

if opts.useGpu
  one = gpuArray(single(1)) ;
else
  one = single(1) ;
end

info.train.objective = [] ;
info.train.error = [] ;
info.train.topFiveError = [] ;
info.train.speed = [] ;
info.val.objective = [] ;
info.val.error = [] ;
info.val.topFiveError = [] ;
info.val.speed = [] ;

lr = opts.learningRate ;
res = [] ;
for epoch=1:opts.numEpochs

  % fast-forward to where we stopped
  modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
  modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;
  if opts.continue
    if exist(modelPath(epoch),'file')
      if epoch == opts.numEpochs
        load(modelPath(epoch), 'net', 'info') ;
      end
      continue ;
    end
    if epoch > 1
      fprintf('resuming by loading epoch %d\n', epoch-1) ;
      load(modelPath(epoch-1), 'net', 'info') ;
    end
  end

  train = opts.train(randperm(numel(opts.train))) ;
  val = opts.val ;

  info.train.objective(end+1) = 0 ;
  info.train.error(end+1) = 0 ;
  info.train.topFiveError(end+1) = 0 ;
  info.train.speed(end+1) = 0 ;
  info.val.objective(end+1) = 0 ;
  info.val.error(end+1) = 0 ;
  info.val.topFiveError(end+1) = 0 ;
  info.val.speed(end+1) = 0 ;

  for t=1:opts.batchSize:numel(train)
    % get next image batch and labels
    batch = train(t:min(t+opts.batchSize-1, numel(train))) ;
    batch_time = tic ;
    fprintf('training: epoch %02d: processing batch %3d of %3d ...', epoch, ...
            fix(t/opts.batchSize)+1, ceil(numel(train)/opts.batchSize)) ;
    [im, labels] = getBatch(imdb, batch) ;
    if opts.prefetch
      nextBatch = train(t+opts.batchSize:min(t+2*opts.batchSize-1, numel(train))) ;
      getBatch(imdb, nextBatch) ;
    end
    if opts.useGpu
      im = gpuArray(im) ;
    end

    % backprop
    net.layers{end}.class = labels ;
    res = vl_simplenn(net, im, one, res, ...
      'conserveMemory', opts.conserveMemory, ...
      'sync', opts.sync) ;

    % gradient step
    for l=1:numel(net.layers)
      if ~strcmp(net.layers{l}.type, 'conv'), continue ; end

      g_f = (net.layers{l}.filtersLearningRate) * ...
            (opts.weightDecay * net.layers{l}.filtersWeightDecay) * net.layers{l}.filters + ...
            (net.layers{l}.filtersLearningRate) / numel(batch) * res(l).dzdw{1};
      g_b = (net.layers{l}.biasesLearningRate) * ...
            (opts.weightDecay * net.layers{l}.biasesWeightDecay) * net.layers{l}.biases + ...
            (net.layers{l}.biasesLearningRate) / numel(batch) * res(l).dzdw{2};

      G_f{l} = G_f{l} + g_f .^ 2;
      G_b{l} = G_b{l} + g_b .^ 2;

      net.layers{l}.filters = net.layers{l}.filters - lr ./ (opts.delta + sqrt(G_f{l})) .* g_f;
      net.layers{l}.biases  = net.layers{l}.biases - lr ./ (opts.delta + sqrt(G_b{l})) .* g_b;
    end

    % print information
    batch_time = toc(batch_time) ;
    speed = numel(batch)/batch_time ;
    info.train = updateError(opts, info.train, net, res, batch_time) ;

    fprintf(' %.2f s (%.1f images/s)', batch_time, speed) ;
    n = t + numel(batch) - 1 ;
    switch opts.errorType
      case 'multiclass'
        fprintf(' err %.1f err5 %.1f', ...
          info.train.error(end)/n*100, info.train.topFiveError(end)/n*100) ;
        fprintf('\n') ;
      case 'binary'
        fprintf(' err %.1f', ...
          info.train.error(end)/n*100) ;
        fprintf('\n') ;
      case 'euclideanloss'
        fprintf(' err %.1f', info.train.error(end) / n);
        fprintf('\n') ;
    end

    % debug info
    if opts.plotDiagnostics
      figure(2) ; vl_simplenn_diagnose(net,res) ; drawnow ;
    end
  end % next batch

  % evaluation on validation set
  if epoch == 1 || rem(epoch, opts.test_interval) == 0 || epoch == opts.numEpochs
  for t=1:opts.batchSize:numel(val)
    batch_time = tic ;
    batch = val(t:min(t+opts.batchSize-1, numel(val))) ;
    fprintf('validation: epoch %02d: processing batch %3d of %3d ...', epoch, ...
            fix(t/opts.batchSize)+1, ceil(numel(val)/opts.batchSize)) ;
    [im, labels] = getBatch(imdb, batch) ;
    if opts.prefetch
      nextBatch = val(t+opts.batchSize:min(t+2*opts.batchSize-1, numel(val))) ;
      getBatch(imdb, nextBatch) ;
    end
    if opts.useGpu
      im = gpuArray(im) ;
    end

    net.layers{end}.class = labels ;
    res = vl_simplenn(net, im, [], res, ...
      'disableDropout', true, ...
      'conserveMemory', opts.conserveMemory, ...
      'sync', opts.sync) ;

    % print information
    batch_time = toc(batch_time) ;
    speed = numel(batch)/batch_time ;
    info.val = updateError(opts, info.val, net, res, batch_time) ;

    fprintf(' %.2f s (%.1f images/s)', batch_time, speed) ;
    n = t + numel(batch) - 1 ;
    switch opts.errorType
      case 'multiclass'
        fprintf(' err %.1f err5 %.1f', ...
          info.val.error(end)/n*100, info.val.topFiveError(end)/n*100) ;
        fprintf('\n') ;
      case 'binary'
        fprintf(' err %.1f', ...
          info.val.error(end)/n*100) ;
        fprintf('\n') ;
      case 'euclideanloss'
        fprintf(' err %.1f', info.val.error(end) / n);
        fprintf('\n') ;
    end
  end
  end

  % save
  info.train.objective(end) = info.train.objective(end) / numel(train) ;
  info.train.error(end) = info.train.error(end) / numel(train)  ;
  info.train.topFiveError(end) = info.train.topFiveError(end) / numel(train) ;
  info.train.speed(end) = numel(train) / info.train.speed(end) ;
  info.val.objective(end) = info.val.objective(end) / numel(val) ;
  info.val.error(end) = info.val.error(end) / numel(val) ;
  info.val.topFiveError(end) = info.val.topFiveError(end) / numel(val) ;
  info.val.speed(end) = numel(val) / info.val.speed(end) ;
  if epoch == 1 || rem(epoch, opts.snapshot) == 0 || epoch == opts.numEpochs
  save(modelPath(epoch), 'net', 'info') ;
  end

  if epoch == 1 || rem(epoch, opts.display) == 0 || epoch == opts.numEpochs
  figure(1) ; clf ;
  subplot(1,2,1) ;
  semilogy(1:epoch, info.train.objective, 'k') ; hold on ;
  semilogy([1 opts.test_interval : opts.test_interval : epoch epoch], info.val.objective([1 opts.test_interval : opts.test_interval : epoch epoch]), 'b') ;
  xlabel('training epoch') ; ylabel('energy') ;
  grid on ;
  h=legend('train', 'val') ;
  set(h,'color','none');
  title('objective') ;
  subplot(1,2,2) ;
  switch opts.errorType
    case 'multiclass'
      plot(1:epoch, info.train.error, 'k') ; hold on ;
      plot(1:epoch, info.train.topFiveError, 'k--') ;
      plot([1 opts.test_interval : opts.test_interval : epoch epoch], info.val.error([1 opts.test_interval : opts.test_interval : epoch epoch]), 'b') ;
      plot([1 opts.test_interval : opts.test_interval : epoch epoch], info.val.topFiveError([1 opts.test_interval : opts.test_interval : epoch epoch]), 'b--') ;
      h=legend('train','train-5','val','val-5') ;
    case 'binary'
      plot(1:epoch, info.train.error, 'k') ; hold on ;
      plot([1 opts.test_interval : opts.test_interval : epoch epoch], info.val.error([1 opts.test_interval : opts.test_interval : epoch epoch]), 'b') ;
      h=legend('train','val') ;
    case 'euclideanloss'
      plot(1 : epoch, info.train.error, 'k'); hold on;
      plot([1 opts.test_interval : opts.test_interval : epoch epoch], info.val.error([1 opts.test_interval : opts.test_interval : epoch epoch]), 'b') ;
      h = legend('train', 'val') ;
  end
  grid on ;
  xlabel('training epoch') ; ylabel('error') ;
  set(h,'color','none') ;
  title('error') ;
  drawnow ;
  print(1, modelFigPath, '-dpdf') ;
  end
end

% -------------------------------------------------------------------------
function info = updateError(opts, info, net, res, speed)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
sz = size(predictions) ;
n = prod(sz(1:2)) ;

labels = net.layers{end}.class ;
info.objective(end) = info.objective(end) + sum(double(gather(res(end).x))) ;
info.speed(end) = info.speed(end) + speed ;
switch opts.errorType
  case 'multiclass'
    [~,predictions] = sort(predictions, 3, 'descend') ;
    error = ~bsxfun(@eq, predictions, reshape(labels, 1, 1, 1, [])) ;
    info.error(end) = info.error(end) +....
      sum(sum(sum(error(:,:,1,:))))/n ;
    info.topFiveError(end) = info.topFiveError(end) + ...
      sum(sum(sum(min(error(:,:,1:5,:),[],3))))/n ;
  case 'binary'
    error = bsxfun(@times, predictions, labels) < 0 ;
    info.error(end) = info.error(end) + sum(error(:))/n ;
  case 'euclideanloss'
    error = euclideanloss(sigmoid(predictions), labels);
    info.error(end) = info.error(end) + error;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function p27_default ( dim_num )

%*****************************************************************************80
%
%% P27_DEFAULT sets default values for problem 27.
%
%  Discussion:
%
%    If a problem uses vector parameters, then the spatial dimension
%    DIM_NUM is needed as input, so that the vector parameters can be
%    properly dimensioned.
%
%    Every problem keeps a count of the number of function calls.  Calling
%    this routine causes that count to be reset to 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the problem.
%
  i4 = 0;
  p27_i4 ( 'D', '*', i4 );

  r8 = 0.0;
  p27_r8 ( 'D', '*', r8 );

  r8vec = [];
  p27_r8vec ( 'D', '*', dim_num, r8vec );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function p=Phi(x)
   p=0.5*(1.+erf(x/sqrt(2)));
end

author: Sable

repo: mcbench-benchmarks


text:
function [scans,dt] = er_getScanGroup(vw,scan)
% [scans,dt] = er_getScanGroup(view,scan);
%
% Gets the group of scans assigned as a group to the
% current scan; if no group of scans is assigned, prompts
% to assign it.
%
% See: er_groupScans.
%
% ras 04/07/04: wrote it.
global dataTYPES;

cdt = vw.curDataType;

if ieNotDefined('scan')
    scan = viewGet(vw, 'curScan');
end

% check if there's a scan group assigned to this scan
if ~isfield(dataTYPES(cdt).scanParams(scan),'scanGroup') || ...
        isempty(dataTYPES(cdt).scanParams(scan).scanGroup)
    er_groupScans(vw);
end

% parse the groupScans text into a data type series
% and a set of scans (e.g. 'Original: 1 2 3'):
txt = dataTYPES(cdt).scanParams(scan).scanGroup;
%colon = findstr(':',txt); %findstr is obsolete, using strfind instead
colon = strfind(txt,':');
dtName = txt(1:colon-1);
scans = str2num(txt(colon+2:end));
dt = cellfind({dataTYPES.name}, dtName);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function out=proj_l1_ball(x,r)
%PROJ_L1_BALL computes the orthogonal projection onto the l1 ball {x: norm(x(:),1) <= r}
%
%  Usage:
%  out = PROJ_L1_BALL(x,[r])
%  =============================================================
%  x:
%  x - point to be projected (vector/matrix)
%  r - a positive scalar [default: 1]
%  ==============================================================
%  Assumptions:
%  r > 0
%  ==============================================================
%  Output:
%  out - projection vector

% This file is part of the FOM package - a collection of first order methods for solving convex optimization problems
% Copyright (C) 2017 Amir and Nili Beck
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

if (nargin < 1)
     error ('usage: prox_l1_ball(x,[r])') ;
end

if (nargin < 2)
    %setting default value to 1
    r = 1;
end

out = sign(x) .* proj_simplex(abs(x),r,'ineq') ;

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
% State transition function for the UNGM-model.
%
% Copyright (C) 2007 Jouni Hartikainen
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.

function x_n = ungm_f(x,param)
n = param(1);
x_n = 0.5*x(1,:) + 25*x(1,:)./(1+x(1,:).*x(1,:)) + 8*cos(1.2*(n-1));
if size(x,1) > 1
   x_n = x_n + x(2,:);
end

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
function value = angle_rad_3d ( p1, p2, p3 )

%*****************************************************************************80
%
%% ANGLE_RAD_3D returns the angle in radians between two rays in 3D.
%
%  Discussion:
%
%    The routine always computes the SMALLER of the two angles between
%    two rays.  Thus, if the rays make an (exterior) angle of
%    1.5 pi radians, the (interior) angle of 0.5 pi radians will be reported.
%
%    X dot Y = Norm(X) * Norm(Y) * Cos ( Angle(X,Y) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(3), P2(3), P3(3), points defining an angle.
%    The rays are P1 - P2 and P3 - P2.
%
%    Output, real VALUE, the angle between the two rays,
%    in radians.  This value will always be between 0 and PI.  If either ray has
%    zero length, then the angle is returned as zero.
%
  dim_num = 3;

  v1norm = sqrt ( sum ( ( p2(1:dim_num) - p1(1:dim_num) ).^2 ) );

  if ( v1norm == 0.0 )
    value = 0.0;
    return
  end

  v2norm = sqrt ( sum ( ( p2(1:dim_num) - p3(1:dim_num) ).^2 ) );

  if ( v2norm == 0.0 )
    value = 0.0;
    return
  end

  dot = ( p1(1:dim_num) - p2(1:dim_num) ) * ( p3(1:dim_num) - p2(1:dim_num) )';

  value = r8_acos ( dot / ( v1norm * v2norm ) );

  return
end


text:
function [ invJ ] = vec2jaclInv( vec )
tolerance = 1e-12;

if size(vec,1) == 3

    phi = vec;

    ph = norm(phi);
    if ph < tolerance
        % If the angle is small, fall back on the series representation
        invJ = vec2jacInvSeries(phi,10);
    else
        axis = phi/norm(phi);
        ph_2 = 0.5*ph;

        invJ =   ph_2 * cot(ph_2)* eye(3)...
               + (1 - ph_2 * cot(ph_2))* axis * axis'...
               - ph_2 * hat(axis);
    end

elseif size(vec,1) == 6

    phi = vec(1:3);
    rho = vec(4:6)

    ph = norm(phi);
    if ph < tolerance
        % If the angle is small, fall back on the series representation
        invJ = vec2jaclInvSeries(phi,10);
    else
        invJsmall = vec2jaclInv( phi );
        Q = vec2Q( vec );
        invJ = [ invJsmall -invJsmall*Q*invJsmall; zeros(3) invJsmall ];
    end
end
end

author: mbrossar


text:
function [trans,rot,scale,skew] = affineDecompose(A)
%
% [trans,rot,scale,skew] = affineDecompose(A)
%
% Decomposes an affine (orthogonal linear) transformation matrix into its
% four componets: translations, rotations scales and skews.
%
% (See affineBuild for more info.)
%
% HISTORY:
%   2004.03.12 RFD (bob@white.stanford.edu) shamelessly copied the core
%   algorithm from spm99's spm_imatrix.m.
%   2005.07.08 ras (sayres at stanford edu) imported into mrVista 2.0

% Translations and zooms
%-----------------------------------------------------------------------
R         = A(1:3,1:3);
C         = chol(R'*R);
P         = [A(1:3,4)' 0 0 0  diag(C)'  0 0 0];
if det(R)<0, P(7)=-P(7);end % Fix for -ve determinants

% Shears
%-----------------------------------------------------------------------
C         = diag(diag(C))\C;
P(10:12)  = C([4 7 8]);
R0        = affineBuild([0 0 0], [0 0 0], P(7:9), P(10:12));
R0        = R0(1:3,1:3);
R1        = R/R0;

% This just leaves rotations in matrix R1
%-----------------------------------------------------------------------
%[          c5*c6,           c5*s6, s5]
%[-s4*s5*c6-c4*s6, -s4*s5*s6+c4*c6, s4*c5]
%[-c4*s5*c6+s4*s6, -c4*s5*s6-s4*c6, c4*c5]

P(5) = asin(rang(R1(1,3)));
if (abs(P(5))-pi/2).^2 < 1e-9,
	P(4) = 0;
	P(6) = atan2(-rang(R1(2,1)), rang(-R1(3,1)/R1(1,3)));
else,
	c    = cos(P(5));
	P(4) = atan2(rang(R1(2,3)/c), rang(R1(3,3)/c));
	P(6) = atan2(rang(R1(1,2)/c), rang(R1(1,1)/c));
end;
trans = P(1:3);
rot = P(4:6);
scale = P(7:9);
skew = P(10:12);
return;

% There may be slight rounding errors making b>1 or b<-1.
function a = rang(b)
a = min(max(b, -1), 1);
return;


text:
classdef prtClusterMeanShiftEuclidean < prtCluster
    % nMaxClusters   Kmeans clustering object
    %
    %    CLUSTER = prtClusterMeanShift returns a MeanShift clustering object.
    %
    %    CLUSTER = prtClusterMeanShift(PROPERTY1, VALUE1, ...) constructs a
    %    prtClusterMeanShift object CLUSTER with properties as specified by
    %    PROPERTY/VALUE pairs.
    %
    %    A prtClusterMeanShift object inherits all properties from
    %    the abstract class prtCluster. In addition is has the following
    %    properties:
    %
    %    nClusters              - Maximum Number of cluster to learn
    %
    %    A prtClusterMeanShift object inherits the TRAIN, RUN,
    %    CROSSVALIDATE and KFOLDS methods from prtAction. It also inherits
    %    the PLOT method from prtCluster.
    %
    %    Invoking the RUN method on a prtClusterMeanShift object classifies
    %    the input data by assigning each observation a label according to
    %    the cluster center it is closest to. The cluster centers are found
    %    during training.
    %
    %   Example:
    %
    %   ds = prtDataGenMary                  % Load a prtDataSet
    %   clusterAlgo = prtClusterMeanShift;   % Create a prtClusterMeanShift object
    %   clusterAlgo.nClusters = 10;       % Set the max number of desired clusters
    %
    %   % Set the internal decision rule to be MAP. Not required for
    %   % clustering, but necessary to plot the results.
    %   clusterAlgo.internalDecider = prtDecisionMap;
    %   clusterAlgo = clusterAlgo.train(ds); % Train the cluster algorithm
    %   plot(clusterAlgo);                   % Plot the results
    %
    %   See also prtCluster, prtClusterGmm

    properties (SetAccess=private)
        name = 'Mean Shift Clustering'
        nameAbbreviation = 'MeanShift'
    end

    properties
        nClusters = 5;

        nMaxIterations = 100;
        meanShiftThreshold = 1e-5;

        meanSeparationThreshold = 1; % Must be exceeded to

        membershipDistance = 1;
        distanceFunction = @(ds,mu)prtDistanceEuclidean(ds,mu);

    end
    properties (SetAccess = protected)
        clusterCenters = [];   % The cluster centers
    end

    methods

        function self = set.nClusters(self,value)
            if ~prtUtilIsPositiveScalarInteger(value)
                error('prt:prtClusterMeanShift:nClusters','value (%s) must be a positive scalar integer',mat2str(value));
            end
            self.nClusters = value;
        end

        function self = set.distanceFunction(self,value)
            if ~isa(value,'function_handle')
                error('prt:prtClusterMeanShift:distanceFunction','distanceFunction must be a function handle');
            elseif nargin(value) ~= 2
                error('prt:prtClusterMeanShift:distanceFunction','distanceFunction must be a function handle that takes two input arguments');
            end
            self.distanceMetricFn = value;
        end
        function self = prtClusterMeanShiftEuclidean(varargin)
            % Allow for string, value pairs
            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end
    end

    methods (Access=protected, Hidden = true)

        function self = trainAction(self,ds)

            X = ds.X;
            beenMappedToAMean = false(ds.nObservations,1);

            means = nan(self.nClusters,ds.nFeatures);
            for iCluster = 1:self.nClusters

                if all(beenMappedToAMean) % Visited every point
                    break
                end

                cPossibleSeedPoints = find(~beenMappedToAMean);
                cSeedPoint = cPossibleSeedPoints(prtRvUtilRandomSample(length(cPossibleSeedPoints)));

                % Pick random point as this mean;
                cMean = X(cSeedPoint,:);

                for iter = 1:self.nMaxIterations
                    cD = self.distanceFunction(X, cMean);
                    isMemberThisCluster = cD < self.membershipDistance;

                    beenMappedToAMean(isMemberThisCluster) = true;

                    newMean = mean(X(isMemberThisCluster,:),1);

                    if norm(cMean-newMean) < self.meanShiftThreshold
                        break
                    else
                        cMean = newMean;
                    end
                end

                % Throw out cluster if too close to existing
                cGoodMeans = means(~any(isnan(means),2),:);
                if ~isempty(cGoodMeans)

                    distanceFromThisMeanToOthers = self.distanceFunction(cMean, cGoodMeans);

                    if min(distanceFromThisMeanToOthers) < self.meanSeparationThreshold
                        cMean = nan(1,size(X,2)); % Set to nan so that we will throw away below.
                    end
                end

                means(iCluster,:) = cMean;
            end

            % Clean up - Remove clusters we didn't use or threw away.
            self.clusterCenters = means(~any(isnan(means),2),:);

        end

        function ds = runAction(self,ds)

            fn = self.distanceFunction;
            distance = fn(ds.getObservations,self.clusterCenters);

            if size(distance,1) ~= ds.nObservations || size(distance,2) ~= size(self.clusterCenters,1)
                error('prt:prtClusterKmeans:badDistanceMetric','Expected a matrix of size %s from the distance metric, but distance metric function output a matrix of size %s',mat2str([DataSet.nObservations, size(Obj.clusterCenters,1)]),mat2str(size(distance)));
            end

            %The smallest distance is the expected class:
            [dontNeed,clusters] = min(distance,[],2);  %#ok<ASGLU>

            binaryMatrix = zeros(size(clusters,1),self.nClusters);
            for i = 1:self.nClusters
                binaryMatrix(i == clusters,i) = 1;
            end
            ds = ds.setObservations(binaryMatrix);

        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function pass = test_rhs( pref )
% Check that non-zero righthand forcing terms are being dealt with
% correctly.
% Alex Townsend, March 2013.

if ( nargin < 1 )
    pref = chebfunpref();
end
tol = 100*pref.cheb2Prefs.chebfun2eps;

d = [0 pi 0 pi];
N = chebop2(@(u) laplacian(u), d);
N.lbc = 0;
N.rbc = @(y) pi*y.^3./6;
N.dbc = 0;
N.ubc = @(x) x*pi^3/6 + sin(x).*sinh(pi);

rhs = chebfun2(@(x,y) x.*y, d);

u = N \ rhs;

exact = chebfun2(@(x,y) x.*y.^3/6 + sin(x).*sinh(y), d);

pass(1) = ( norm(u - exact) < 10*tol );

d = [0 pi 0 pi];
N = chebop2(@(u) laplacian(u), d);
N.lbc = 0;
N.rbc = @(y) pi^4*y./12 + sinh(pi)*(cos(y) + sin(y));
N.dbc = @(x) sinh(x);
N.ubc = @(x) pi*x.^4/12 - sinh(x);

rhs = chebfun2(@(x,y) x.^2.*y, d);

u = N \ rhs;

exact = chebfun2(@(x,y) x.^4.*y./12 + sinh(x).*(cos(y) + sin(y)), d);

pass(2) = ( norm(u - exact) < 10*tol );

% right hand side when we are on a rectangular domain.
% real part of analytic plus forcing term
d = [0 pi 0 1];
exact = chebfun2(@(x,y) real(exp(x+1i*y)) + x.^3.*y.^3, d);
N = chebop2(@(u) laplacian(u), d);
N.lbc = exact(d(1),:);
N.rbc = exact(d(2),:);
N.dbc = exact(:,d(3));
N.ubc = exact(:,d(4));

u = N \ chebfun2(@(x,y) 6*x.*y.^3  + 6*y.*x.^3, d);

pass(3) = ( norm(u - exact) < 10*tol );

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [data g] = pre_detrend(varargin)
%
% Detrend each channel and trial using a least-squares linear fit.
% Alternately, one can center the data by removing the mean from each
% trial. See [1] for more details.
%
% Inputs:
%
%   data:           raw EEG data
%   SamplingRate:   data sampling rate
%
% Optional:     <'Name',Value> pairs
%
%     VerbosityLevel:   Verbosity level. 0 = no output, 1 = text, 2 = graphical
%                       Possible values: 0,1,2
%                       Default value  : 1
%                       Input Data Type: real number (double)
%
%     DetrendingMethod: Detrending options
%                       Linear: removes the least-squares fit of a straight line from each trial.
%                       Constant: removes the mean from each trial (centering)
%                       Possible values: 'linear','constant'
%                       Default value  : 'linear'
%                       Input Data Type: boolean
% Outputs:
%
%   data:       processed EEG data
%   g:          argument specification structure
%
% See Also: pop_pre_prepData(), pre_prepData(), detrend()
%
% Refences:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual. Section 6.5.1
%   Available at: http://www.sccn.ucsd.edu/wiki/Sift
%
% Author: Tim Mullen 2010, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

SegLenRange = [eps 10];
StepSizeRange = [eps 10];

g = arg_define(varargin, ...
    arg_norep('data',mandatory,[],'Raw data. This should be [num_chans x num_time x num_trials]'), ...
    arg_norep({'srate','SamplingRate'},mandatory,[],'Sampling Rate'), ...
    arg({'verb','VerbosityLevel'},1,{int32(0) int32(1) int32(2)}, ...
    'Verbosity level. 0 = no output, 1 = text, 2 = graphical'), ...
    arg({'method','DetrendingMethod'},{'linear'},{'linear','constant'}, ...
    'Detrending options. Linear: removes the least-squares fit of a straight line from each trial. Constant: removes the mean from each trial (centering)','type','logical'), ...
    arg_subtoggle({'piecewise','Piecewise'},[], ...
    { ...
    arg({'seglength','SegmentLength'},0.33,SegLenRange,'Length of each detrending segment (sec).'), ...
    arg({'stepsize','StepSize'},0.0825,StepSizeRange,'Step size between segment centers (sec). It is recommended to use at least 0.5*SegmentLength (50% overlap).'), ...
    },'Use piecewise detrending. Divide the data into (overlapping) segments and detrend each segment separately. Segment endpoints are merged and stitched together using a cubic spline function to minimize discontinuities at segment intersection points. This is useful for as an alternative to high-pass filtering for removing infraslow oscillations (e.g. SC drift) from the '), ...
    arg({'plot','Plot'},false,[],'Plot results for inspection.') ...
    );

% commit data to workspace
data = g.data;
g=rmfield(g,'data');

[nbchan pnts trials] = size(data);

% set sliding-window parameters
if ~g.piecewise.arg_selection
    % global detrending
    g.piecewise.seglength = pnts/g.srate;
    g.piecewise.stepsize = g.piecewise.seglength;
end
windowing_params = [g.piecewise.seglength g.piecewise.stepsize];

% apply each detrending method (centering,linear detrending, ...)
for i=1:length(g.method)

    % initialize verbose output
    m = g.method{i};
    if g.verb && g.piecewise.arg_selection,
        pre = 'Piecewise ';
        post = sprintf(' (using %0.4f sec segment length)',g.piecewise.seglength);
    else
        pre = '';
        post = '';
    end
    if g.verb && strcmpi(m,'mean')
        fprintf('%sCentering data%s...\n',pre,post);    end
    if g.verb && strcmpi(m,'linear')
        fprintf('%sDetrending data%s...\n',pre,post);   end

    if g.verb==2
        multiWaitbar('Detrending','Reset','Color',hlp_getNextUniqueColor);
    end
end

% variable intialization
fitlines = zeros(size(data));

% apply the detrending/centering
for ch=1:nbchan

    if g.verb==2
        multiWaitbar('Detrending',ch/nbchan);
    end

    if g.plot
        % return detrended data as well as fitted curves
        [data(ch,:,:) fitlines(ch,:,:)] = locdetrend_siftmod(squeeze(data(ch,:,:)),g.srate,windowing_params,m);
    else
        % return only detrended data (faster)
        data(ch,:,:) = locdetrend_siftmod(squeeze(data(ch,:,:)),g.srate,windowing_params,m);
    end

end

if g.verb==2
    multiWaitbar('Detrending','Close');
end

% plot results, if requested
if g.plot
    eegplot(data+fitlines,'srate',g.srate,'data2',fitlines,'title','Original Data');
    h = gcf;
    ax = findobj(gcf,'tag','eegaxis');
    plts = get(ax,'children');
    legend([plts(end) plts(1)],'original','best local-linear fit');
    eegplot(data,'srate',g.srate,'title','Detrended Data','children',h);
    ax = findobj(gcf,'tag','eegaxis');
    plts = get(ax,'children');
    legend(plts(end),'detrended data');
end

repo: REST


text:
% function to calculate lcm of two number a and b using gcd of two numbers.

function [ret] =  lcm(a, b)
    if b == 0
       re = a;
       ret=(a*b)/re;    % i.e zero

    else
       re = gcd(b, mod(a, b));
       ret=(a*b)/re;
    end

author: TheAlgorithms

repo: MATLAB-Octave


text:
function [errorCode] = pco_sdk_example_live_add(looptime,exposure_time,triggermode)
% grab and display images in a loop
%
%   [errorCode] = pco_sdk_example_live_add(looptime,triggermode)
%
% * Input parameters :
%    looptime                time the loop is running (default=10 seconds)
%    exposure_time           camera exposure time (default=10ms)
%    triggermode             camera trigger mode (default=AUTO)
%
% * Output parameters :
%    errorCode               ErrorCode returned from pco.camera SDK-functions
%
%grab images from a recording pco.camera
%using functions PCO_AddBufferEx and PCO_WaitforBuffer
%display the grabbed images
%
%
%break loop either after waittime or nr_of_images images are done
%

glvar=struct('do_libunload',1,'do_close',0,'camera_open',0,'out_ptr',[]);

if(~exist('looptime','var'))
 looptime = 10;
end

if(~exist('exposure_time','var'))
 exposure_time = 10;
end

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

%reduce_display_size=1: display only top-left corner 800x600Pixel
reduce_display_size=1;

pco_camera_load_defines();
subfunc=pco_camera_subfunction();

[errorCode,glvar]=pco_camera_open_close(glvar);
pco_errdisp('pco_camera_setup',errorCode);
disp(['camera_open should be 1 is ',int2str(glvar.camera_open)]);
if(errorCode~=PCO_NOERROR)
 commandwindow;
 return;
end

out_ptr=glvar.out_ptr;

try

act_recstate = uint16(10);
[errorCode,~,act_recstate] = calllib('PCO_CAM_SDK', 'PCO_GetRecordingState',out_ptr,act_recstate);
if(errorCode)
  pco_errdisp('PCO_GetRecordingState',errorCode);
else
 disp(['actual recording state is ',int2str(act_recstate)]);
end

%stop camera
subfunc.fh_stop_camera(out_ptr);

cam_desc=libstruct('PCO_Description');
set(cam_desc,'wSize',cam_desc.structsize);
[errorCode,~,cam_desc] = calllib('PCO_CAM_SDK', 'PCO_GetCameraDescription', out_ptr,cam_desc);
cam_desc
pco_errdisp('PCO_GetCameraDescription',errorCode);

bitpix=uint16(cam_desc.wDynResDESC);

%set bitalignment LSB
bitalign=uint16(BIT_ALIGNMENT_LSB);
errorCode = calllib('PCO_CAM_SDK', 'PCO_SetBitAlignment', out_ptr,bitalign);
pco_errdisp('PCO_SetBitAlignment',errorCode);

errorCode = calllib('PCO_CAM_SDK', 'PCO_SetRecorderSubmode',out_ptr,RECORDER_SUBMODE_RINGBUFFER);
pco_errdisp('PCO_SetRecorderSubmode',errorCode);

%enable ASCII and binary timestamp
if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_NO_TIMESTAMP)==0)
 subfunc.fh_enable_timestamp(out_ptr,TIMESTAMP_MODE_BINARYANDASCII);
end

%enable MetaData if available
if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_METADATA))
 subfunc.fh_set_metadata_mode(out_ptr,0);
end

	%{
	%Pixel Binning
	%binning funktioniert nur wenn gleichzeitig ROI gesetzt wird.
	h_binning=1; %1,2,4
	v_binning=1; %1,2,4

	v_ROI_reduction = 4
	[errorCode] = calllib('PCO_CAM_SDK', 'PCO_SetBinning', out_ptr,h_binning,v_binning); %2,4, etc.
	%ROI selection
	[errorCode] = calllib('PCO_CAM_SDK', 'PCO_SetROI', out_ptr,1,1,5120/h_binning/4,5120/v_binning/v_ROI_reduction); %does this work for the panda...?

%}

%set default Pixelrate
subfunc.fh_set_pixelrate(out_ptr,2);

subfunc.fh_set_triggermode(out_ptr,triggermode);
subfunc.fh_set_exposure_times(out_ptr,exposure_time,2,0,2);

%if PCO_ArmCamera does fail no images can be grabbed
errorCode = calllib('PCO_CAM_SDK', 'PCO_ArmCamera', out_ptr);
if(errorCode~=PCO_NOERROR)
 pco_errdisp('PCO_ArmCamera',errorCode);
 ME = MException('PCO_ERROR:ArmCamera','Cannot continue script with not armed camera');
 subfunc.fh_lasterr(errorCode);
 throw(ME);
end

%adjust transfer parameter if necessary
subfunc.fh_set_transferparameter(out_ptr);

triggermode=subfunc.fh_get_triggermode(out_ptr);

%display actual image time and maximal frequency
imatime=subfunc.fh_show_frametime(out_ptr)

%%
%% Measure time to acquire 1 image
		dwSec=uint32(0);
		dwNanoSec=uint32(0);
		[errorCode,~,dwSec,dwNanoSec] = calllib('PCO_CAM_SDK', 'PCO_GetCOCRuntime', out_ptr,dwSec,dwNanoSec);
		if(errorCode)
			pco_errdisp('PCO_GetCOCRuntime',errorCode);
		end
		disp(['Max double image capture freq: ' num2str(round(1/(double(dwNanoSec)/1000/1000/1000),3)) ' Hz.'])
%%

metadatasize=0;
if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_METADATA))
 wMetaDataMode=uint16(0);
 wMetaDataSize=uint16(0);
 wMetaDataVersion=uint16(0);

 [errorCode,~,wMetaDataMode,wMetaDataSize,wMetaDataVersion] = calllib('PCO_CAM_SDK', 'PCO_GetMetaDataMode',out_ptr,wMetaDataMode,wMetaDataSize,wMetaDataVersion);
 pco_errdisp('PCO_GetMetaDataMode',errorCode);
 if(errorCode)
  pco_errdisp('PCO_GetMetaDataMode',errorCode);
  return;
 end

 disp(['Camera supports Metadata Version: ',num2str(wMetaDataVersion)]);

 if(wMetaDataMode~=0)
  metadatasize=wMetaDataSize;
  disp(['Metadata is enabled. Number of Bytes: ',num2str(metadatasize)]);
 end
end

%calculate images to grab
nr_of_images=uint32(fix(looptime/imatime)+1);
disp(['maximal ',num2str(nr_of_images),' images will be grabbed in ',num2str(looptime),' seconds' ]);

act_xsize=uint16(0);
act_ysize=uint16(0);
ccd_xsize=uint16(0);
ccd_ysize=uint16(0);

%use PCO_GetSizes because this always returns accurate image size for next recording
[errorCode,~,act_xsize,act_ysize]  = calllib('PCO_CAM_SDK', 'PCO_GetSizes', out_ptr,act_xsize,act_ysize,ccd_xsize,ccd_ysize);
if(errorCode)
 pco_errdisp('PCO_GetSizes',errorCode);
end

disp(['sizes: horizontal ',int2str(act_xsize),' vertical ',int2str(act_ysize)]);

flags=2; %IMAGEPARAMETERS_READ_WHILE_RECORDING;
errorCode = calllib('PCO_CAM_SDK', 'PCO_SetImageParameters', out_ptr,act_xsize,act_ysize,flags,[],0);
pco_errdisp('PCO_SetImageParameters',errorCode);

lineadd=0;
if(metadatasize>0)
 xs=uint32(fix((double(bitpix)+7)/8));
 xs=xs*uint32(act_xsize);
 i=uint16(floor((metadatasize*2)/double(xs)));
 lineadd=i+1;
end

disp(['lines added: ',int2str(lineadd)]);

%allocate memory for display, 4 buffers are used
bufcount=4;
bufnum=zeros(4,1,'int16');

imas=uint32(fix((double(bitpix)+7)/8));
imas= imas*uint32(act_xsize)* uint32(act_ysize);
imasize=imas;

image_stack=zeros(act_xsize,(act_ysize+lineadd),bufcount,'uint16');

%Allocate 4 SDK buffer and set address of buffers from image_stack
ev_ptr(bufcount) = libpointer('voidPtr');
im_ptr(bufcount) = libpointer('voidPtr');

buflist=libstruct('PCO_Buflist');
names=fieldnames(buflist);
x=1:4:length(names);
bufnr=names(x);
x=3:4:length(names);
statusdll=names(x);
x=4:4:length(names);
statusdrv=names(x);

for n=1:bufcount
 sBufNri=int16(-1);
 im_ptr(n) = libpointer('uint16Ptr',image_stack(:,:,n));
 ev_ptr(n) = libpointer('voidPtr');

 [errorCode,~,sBufNri]  = calllib('PCO_CAM_SDK','PCO_AllocateBuffer', out_ptr,sBufNri,imasize,im_ptr(n),ev_ptr(n));
 if(errorCode~=PCO_NOERROR)
  pco_errdisp('PCO_AllocateBuffer',errorCode);
  ME = MException('PCO_ERROR:AllocateBuffer','Cannot continue script without allocated Buffers');
  subfunc.fh_lasterr(errorCode);
  throw(ME);
 end
 buflist.(bufnr{n})=int16(sBufNri);
 bufnum(n)=int16(sBufNri);
 e=bitor(buflist.(statusdll{n}),uint32(PCO_BUFFER_EVAUTORES),'uint32');
 buflist.(statusdll{n})=e;
end

%show figure
ima=image_stack(:,1:act_ysize,1);
if(reduce_display_size~=0)
 [xs,ys]=size(ima);
 xmax=800;
 ymax=600;
 if((xs>xmax)&&(ys>ymax))
  ima=ima(1:xmax,1:ymax);
 elseif(xs>xmax)
  ima=ima(1:xmax,:);
 elseif(ys>ymax)
  ima=ima(:,1:ymax);
 end
end
ima=ima';

imah=draw_image(ima,[0 100]);
axish=gca;
set(axish,'CLim',[0 1000]);

pause(0.0001);

%grab preimage to get actual image value range and set limits
errorCode=subfunc.fh_start_camera(out_ptr);
if(errorCode~=PCO_NOERROR)
 ME = MException('PCO_ERROR:StartCamera','Cannot continue script with stopped camera');
 subfunc.fh_lasterr(errorCode);
 throw(ME);
end

errorCode  = calllib('PCO_CAM_SDK','PCO_AddBufferEx', out_ptr,0,0,bufnum(1),act_xsize,act_ysize,bitpix);
if(errorCode~=PCO_NOERROR)
 pco_errdisp('PCO_AddBufferEx',errorCode);
 ME = MException('PCO_ERROR:AddBufferEx','Cannot continue script without added Buffers');
 subfunc.fh_lasterr(errorCode);
 throw(ME);
end

disp('get pre images');
trigdone=int16(1);
if((triggermode==1)||(triggermode==2))
 [errorCode,~,trigdone]  = calllib('PCO_CAM_SDK','PCO_ForceTrigger',out_ptr,trigdone);
 if(errorCode)
  pco_errdisp('PCO_ForceTrigger',errorCode);
 else
  disp([int2str(trigdone),' trigger done return: ',int2str(trigdone)]);
 end
elseif(triggermode>2)
 disp('send external trigger pulse within 1 second');
 pause(0.001);
end

[errorCode,~,buflist]  = calllib('PCO_CAM_SDK','PCO_WaitforBuffer', out_ptr,1,buflist,1000);
if(errorCode)
 pco_errdisp('PCO_WaitforBuffer',errorCode);
else
 disp(['PCO_WaitforBuffer done bufnr: ',int2str(bufnum(1))]);
end

if((bitand(buflist.dwStatusDll_1,hex2dec('00008000')))&&(buflist.dwStatusDrv_1==0))
 disp(['Event buf',int2str(1),' pre image done, StatusDrv 0x',num2str(buflist.dwStatusDrv_1,'%08X'),' Statusdll 0x',num2str(buflist.dwStatusDll_1,'%08X')]);
%get data and show image
 ima=get(im_ptr(1),'Value');
 if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_NO_TIMESTAMP)==0)
  [imanum,~]=subfunc.fh_get_timestamp(ima,bitalign,bitpix);
  disp(['get_timestamp imacount: ',int2str(imanum),' done ']);
  txt=subfunc.fh_print_timestamp(ima,bitalign,bitpix);
  disp(['timestamp ima:  ',txt]);
  str_time=subfunc.fh_get_struct_timestamp(ima,bitalign,bitpix);
  txt=subfunc.fh_print_struct_timestamp(str_time);
  disp(['timestamp str:  ',txt]);
 end

 if(metadatasize>0)
  metaline=ima(1:metadatasize,act_ysize+1:end);
  [errorCode_m,metastruct]=subfunc.fh_get_struct_metadata(metaline,metadatasize);
  if(errorCode_m)
   pco_errdisp('get_struct_metadata',errorCode_m);
  else
   txt=subfunc.fh_print_meta_timestamp(metastruct);
   disp(['timestamp meta: ',txt]);
  end
 end

 if(metadatasize>0)
  ima=ima(:,1:act_ysize);
 end

 if(reduce_display_size~=0)
  [xs,ys]=size(ima);
  xmax=800;
  ymax=600;
  if((xs>xmax)&&(ys>ymax))
   ima=ima(1:xmax,1:ymax);
  elseif(xs>xmax)
   ima=ima(1:xmax,:);
  elseif(ys>ymax)
   ima=ima(:,1:ymax);
  end
 end
 ima=ima';
 if(bitalign==BIT_ALIGNMENT_MSB)
  s=int16(16-bitpix);
  s=s*-1;
  ima=bitshift(ima,s);
 end

 m=max(max(ima(10:end-10,10:end-10)));
 set(axish,'CLim',[0 m+100]);
 disp(['pre image done maxvalue: ',int2str(m)]);
 set(imah,'CData',ima,'CDataMapping','scaled');
 pause(0.0001);
end

subfunc.fh_stop_camera(out_ptr);

%this will remove all pending buffers in the queue and does reset grabber
errorCode = calllib('PCO_CAM_SDK', 'PCO_CancelImages', out_ptr);
pco_errdisp('PCO_CancelImages',errorCode);

pause(0.01);
%variable to reduce amount of messages
d=10;
if(nr_of_images>100)
 if(nr_of_images<500)
  d=50;
 else
  d=100;
 end
end

errorCode=subfunc.fh_start_camera(out_ptr);
if(errorCode~=PCO_NOERROR)
 ME = MException('PCO_ERROR:StartCamera','Cannot continue script with stopped camera');
 subfunc.fh_lasterr(errorCode);
 throw(ME);
end

%setup loop
for n=1:bufcount
 errorCode = calllib('PCO_CAM_SDK','PCO_AddBufferEx', out_ptr,0,0,bufnum(n),act_xsize,act_ysize,bitpix);
 if(errorCode~=PCO_NOERROR)
  pco_errdisp('PCO_AddBufferEx',errorCode);
  ME = MException('PCO_ERROR:AddBufferEx','Cannot continue script without added buffers');
  subfunc.fh_lasterr(errorCode);
  throw(ME);
 end
end

[errorCode,~,act_recstate] = calllib('PCO_CAM_SDK', 'PCO_GetRecordingState',glvar.out_ptr,act_recstate);
if(errorCode)
  pco_errdisp('PCO_GetRecordingState',errorCode);
else
 disp(['Actual recording state is ',int2str(act_recstate)]);
end

if(act_recstate==1)
 disp('get images');

 trigdone=int16(1);
 trigcount=0;
 if(triggermode==1)
  [errorCode,~,trigdone]  = calllib('PCO_CAM_SDK','PCO_ForceTrigger',out_ptr,trigdone);
  if(errorCode)
   pco_errdisp('PCO_ForceTrigger',errorCode);
  else
   trigcount=trigcount+1;
   disp(['first trigger done return: ',int2str(trigdone)]);
  end
 elseif(triggermode>=2)
  disp('send external trigger pulses within 1 second');
 end
 pause(0.0001);
 tic;

%grab and display loop
 ima_nr=0;
 last_ok=0;
 while(ima_nr<nr_of_images)
%wait for buffers
  [errorCode,~,buflist]  = calllib('PCO_CAM_SDK','PCO_WaitforBuffer', out_ptr,bufcount,buflist,1000);
  if(errorCode)
   pco_errdisp('PCO_WaitforBuffer',errorCode);
   break;
  end
%  disp(['After wait image ',int2str(ima_nr),' last_ok ',num2str(last_ok)]);
%first image done trigger next
  if(triggermode==1)
   [errorCode,~,trigdone]  = calllib('PCO_CAM_SDK','PCO_ForceTrigger',out_ptr,trigdone);
   if(errorCode)
    pco_errdisp('PCO_ForceTrigger',errorCode);
   else
    trigcount=trigcount+1;
   end
  end

%test and display buffer
  next=last_ok+1;
  multi=0;
  for n=1:bufcount
   if(next>bufcount)
    next=1;
   end
%   disp(['Status buf',int2str(next),' StatusDll ',num2str(buflist.(statusdll{next}),'%08X'),' StatusDrv ',num2str(buflist.(statusdrv{next}),'%08X')]);
   if((bitand(buflist.(statusdll{next}),hex2dec('00008000')))&&(buflist.(statusdrv{next})==0))
%get data and show image
    last_ok=next;
    ima=get(im_ptr(next),'Value');
    ima_nr=ima_nr+1;
    multi=multi+1;

    if(rem(ima_nr,d)==0)
     disp(['Status buf',int2str(next),' ok  image ',int2str(ima_nr),' done, last_ok ',num2str(last_ok)]);
     if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_NO_TIMESTAMP)==0)
      txt=subfunc.fh_print_timestamp(ima,bitalign,bitpix);
      disp(['Timestamp of image(',num2str(ima_nr,'%04d'),'): ',txt]);
     end

     if(metadatasize>0)
      metaline=ima(1:metadatasize,act_ysize+1:end);
      [errorCode_m,metastruct]=subfunc.fh_get_struct_metadata(metaline,metadatasize);
      if(errorCode_m)
       pco_errdisp('get_struct_metadata',errorCode_m);
      else
       txt=subfunc.fh_print_meta_timestamp(metastruct);
       disp(['Timestamp of meta (',num2str(ima_nr,'%04d'),'): ',txt]);
      end
     end
    end

    if(metadatasize>0)
     ima=ima(:,1:act_ysize);
    end

    if(reduce_display_size~=0)
     [xs,ys]=size(ima);
     xmax=800;
     ymax=600;
     if((xs>xmax)&&(ys>ymax))
      ima=ima(1:xmax,1:ymax);
     elseif(xs>xmax)
      ima=ima(1:xmax,:);
     elseif(ys>ymax)
      ima=ima(:,1:ymax);
     end
    end
    ima=ima';
    if(bitalign==BIT_ALIGNMENT_MSB)
     s=int16(16-bitpix);
     s=s*-1;
     ima=bitshift(ima,s);
    end
    set(imah,'CData',ima,'CDataMapping','scaled');
	pause(0.0001);
	assignin('base','letztes_bild',ima);

    errorCode = calllib('PCO_CAM_SDK','PCO_AddBufferEx', out_ptr,0,0,bufnum(next),act_xsize,act_ysize,bitpix);
    if(errorCode)
     pco_errdisp('PCO_AddBufferEx',errorCode);
     break;
    end
   end
   next=next+1;
  end

  if(multi>1)
   disp(['Multi Buffers found after wait:',num2str(multi)]);
  end
  t=toc;
  if(t>looptime)
   break;
  end
 end

 disp(['Last image ',int2str(ima_nr),' done ']);
 disp([int2str(ima_nr),' images done in ',num2str(t),' seconds. time per image is ',num2str(t/double(ima_nr),'%.3f'),'s']);
 if(bitand(cam_desc.dwGeneralCapsDESC1,GENERALCAPS1_NO_TIMESTAMP)==0)
  subfunc.fh_print_timestamp_t(ima,1,bitpix); %image is already shifted
 end

%this will remove all pending buffers in the queue
 errorCode = calllib('PCO_CAM_SDK', 'PCO_CancelImages', out_ptr);
 pco_errdisp('PCO_CancelImages',errorCode);

 [errorCode,~,buflist] = calllib('PCO_CAM_SDK','PCO_WaitforBuffer', out_ptr,bufcount,buflist,500);
 pco_errdisp('PCO_WaitforBuffer',errorCode);

 for next=1:bufcount
  disp(['Event buf',int2str(next),' StatusDll ',num2str(buflist.(statusdll{next}),'%08X'),' StatusDrv ',num2str(buflist.(statusdrv{next}),'%08X')]);
 end

 %disp('Press "Enter" to close window and proceed')
 %pause();
 %close();
 %pause(1);
end

subfunc.fh_stop_camera(out_ptr);

for n=1:bufcount
 errorCode = calllib('PCO_CAM_SDK','PCO_FreeBuffer',out_ptr,bufnum(n));
 if(errorCode)
  pco_errdisp('PCO_FreeBuffer',errorCode);
 else
  disp(['PCO_FreeBuffer',num2str(n),' done ']);
 end
end

catch ME
 errorCode=subfunc.fh_lasterr();
 txt=blanks(101);
 txt=calllib('PCO_CAM_SDK','PCO_GetErrorTextSDK',pco_uint32err(errorCode),txt,100);

 calllib('PCO_CAM_SDK', 'PCO_CancelImages', out_ptr);
 for n=1:bufcount
  calllib('PCO_CAM_SDK','PCO_FreeBuffer',out_ptr,bufnum(n));
 end

 clearvars -except ME glvar errorCode txt;

 if(glvar.camera_open==1)
  glvar.do_close=1;
  glvar.do_libunload=1;
  pco_camera_open_close(glvar);
 end

 if strfind(ME.identifier,'PCO_ERROR:')
  msg=[ME.identifier,' ',ME.message];
  warning('off','backtrace')
  warning(msg)
  disp(txt);
  for k=1:length(ME.stack)
    disp(['from file ',ME.stack(k).file,' at line ',num2str(ME.stack(k).line)]);
  end
  close();
  clearvars -except errorCode;
  commandwindow;
  return;
 else
  close();
  clearvars -except ME;
  rethrow(ME)
 end
end

clearvars -except glvar errorCode;

if(glvar.camera_open==1)
 glvar.do_close=1;
 glvar.do_libunload=1;
 pco_camera_open_close(glvar);
end

clearvars glvar;
commandwindow;

end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function spm_cost_SHC_path(qU,A)
% plots path for cost_SHC demo's
% FORMAT spm_cost_SHC_path(qU,A)
%
% qU  - DEM condotioal esimates of states
% A.x - locations of attrcuor
% A.d - radius
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_cost_SHC_path.m 3757 2010-03-08 11:41:53Z guillaume $

% plot
%==========================================================================
[X Y Z] = sphere;

% plot path
%--------------------------------------------------------------------------
x     = qU.v{1}(1:2,:);
plot(A.x(1,:),A.x(2,:),'c.','MarkerSize',64), hold on
plot(x(1,:),x(2,:));
plot(x(1,:),x(2,:),'.','MarkerSize',8);

% plot locations
%--------------------------------------------------------------------------
% for i = 1:size(A.x,2)
%     surf(X*A.d + A.x(1,i),Y*A.d + A.x(2,i),Z*A.d - 8)
% end
% shading interp

axis([-4 4 -4 4])
axis square
title('trajectory','FontSize',16),hold off

% occupancy
%--------------------------------------------------------------------------
for i = 1:size(A.x,2)
    for j = 1:2
        d(j,:) = x(j,:) - A.x(j,i);
    end
    b(i,:) = sqrt(sum(d.^2)) < 2*A.d;
end

% plot precent occupancy
%--------------------------------------------------------------------------
p     = 100*sum(b,2)/length(b);
for i = 1:length(p)
    text(A.x(1,i)*2,A.x(2,i)*1.4,sprintf('%2.0f%%',p(i)),'Fontsize',12)
end
hold off

author: spm


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
% ##2
%==============================================================================
%
%  Example for 2D registration of hand data,
%
%   - data                 hands, omega=(0,20)x(0,25), level=3:7, m=[128,128]
%   - viewer               viewImage2D
%   - image model          splineInter
%   - distance             SSD
%   - pre-registration     affine2D
%   - regularizer          mfElastic
%   - optimization         Gauss-Newton
% ===============================================================================

close all, help(mfilename)

setup2DhandData
viewImage('reset','viewImage','viewImage2D','colormap',bone(256),'axis','off');
imgModel('reset','imgModel','splineInter','regularizer','moments','theta',1e-2);
trafo('reset','trafo','affine2D');
distance('reset','distance','SSD');
regularizer('reset','regularizer','mfElastic','alpha',1e3,'mu',1,'lambda',0);

[yc,wc,his] = MLIR(ML,'parametric',1,'plotIter',1,'plotMLiter',1);

showResults(ML,yc)

%==============================================================================


text:
function [x,w]=ffejer2(N, a, b)
%Fejer
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

N2 = N+1;
x = a+(b-a)*0.5*(cos(pi*(0:N2)/N2)+1)';
x = x(2:N+1);
w = zeros(N2,1);

for i=0:N2
    curw = 0;
    for j=1:floor(N2/2)
        curw = curw + sin((2*j-1)*i*pi/N2)/(2*j-1);
    end;
    w(i+1)=curw*4/N2*sin(i*pi/N2);
end;
w = w(2:N+1);
w = w*(b-a)/2;

end

author: oseledets

repo: TT-Toolbox


text:
function reconstructSelSearchHierarchyFromFz(varargin)
% reconstructSelSearchHierarchyFromFz(varargin)
%
% Reconstruct a Selective Search hierachy to find out which superpixels make up the image.
%
% Note that labelHistos require pixel-level ground-truth.
%
% Copyright by Holger Caesar, 2015

% Parse input
p = inputParser;
addParameter(p, 'dataset', SiftFlowDatasetMC);
addParameter(p, 'projectName', 'WeaklySupervisedLearning');
addParameter(p, 'segmentationNameSS', 'Uijlings2013-ks100-sigma0.8-colorTypesRgb'); %can also be used for Felzen regions (but then has to be equal to segmentationNameFz)
addParameter(p, 'segmentationNameFz', 'Felzenszwalb2004-k100-sigma0.8-colorTypesRgb');
addParameter(p, 'subset', 'all');
addParameter(p, 'overwrite', false); % replaces the fields added by this script, but keeps all other fields
addParameter(p, 'randomOrder', false);
parse(p, varargin{:});

dataset = p.Results.dataset;
projectName = p.Results.projectName;
segmentationNameSS = p.Results.segmentationNameSS;
segmentationNameFz = p.Results.segmentationNameFz;
subset = p.Results.subset;
overwrite = p.Results.overwrite;
randomOrder = p.Results.randomOrder;

% Make program reproducible
rng(42);

% Init
globalTimer = tic;
labelCount = dataset.labelCount; %#ok<NASGU>
imagesProcessed = 0;

% Get image list and image size
[imageList, imageCount] = dataset.getImageListSubset(subset);

% Randomize image order
if randomOrder,
    % Store old randomness
    seed = rng;

    % Create randomness and shuffle image list
    rng('shuffle');
    imageList = imageList(randperm(imageCount));

    % Reset
    rng(seed.Seed);
end;

% Create paths
global glFeaturesFolder;
segmentationFolderSS = fullfile(glFeaturesFolder, projectName, dataset.name, 'segmentations', segmentationNameSS);
segmentationFolderFz = fullfile(glFeaturesFolder, projectName, dataset.name, 'segmentations', segmentationNameFz);

for imageIdx = 1 : imageCount,
    imagesRemaining = imageCount - imageIdx;
    secsPerImage = toc(globalTimer) / imagesProcessed;
    minutesLeft = secsPerImage / 60 * imagesRemaining;
    fprintf('Reconstructing hierarchy for image %d of %d (%.1fmin left) ##\n', imageIdx, imageCount, minutesLeft);

    imageName = imageList{imageIdx};
    imageSize = dataset.getImageSize(imageName);

    % Get SS propBlobs
    segmentPathSS = fullfile(segmentationFolderSS, [imageName, '.mat']);
    segmentStructSS = load(segmentPathSS);
    propBlobsSS = segmentStructSS.propBlobs;
    if ~overwrite,
        if isfield(segmentStructSS, 'overlapList'),
            fprintf('Skipping file which has already been processed: %s\n', segmentPathSS);
            continue;
        end;
    end;

    % Get Fz propBlobs
    segmentPathFz = fullfile(segmentationFolderFz, [imageName, '.mat']);
    segmentStructFz = load(segmentPathFz, 'propBlobs');
    propBlobsFz = segmentStructFz.propBlobs;

    % Find superpixels in region proposals
    spCount = numel(propBlobsFz);
    blobSizesSS = [propBlobsSS.size]';
    for spIdx = 1 : spCount,
        candidates = find(blobSizesSS == propBlobsFz(spIdx).size);

        for candidateIdx = 1 : numel(candidates),
            regionIdx = candidates(candidateIdx);
            if isequal(propBlobsFz(spIdx), propBlobsSS(regionIdx)),
                break;
            end;
        end;
    end;

    % Find Fz segments among the SS regions
    overlapList = sparse(computeBlobOverlapAny(propBlobsFz, propBlobsSS, imageSize))';

    % Construct GT label histograms
    superPixelLabelHistos = [];

    % Store new information
    segmentStructSS.superPixelLabelHistos = superPixelLabelHistos;
    segmentStructSS.overlapList = overlapList;
    save(segmentPathSS, '-struct', 'segmentStructSS', '-v6');
    imagesProcessed = imagesProcessed + 1;
end;

author: nightrome

repo: matconvnet-calvin


text:
%% plot array map & compute eastings and northings
if make_figures
    disp('Plotting array map')
    close all
    cols = 'rwbggg';
    for c=1:length(lat)
        chan = get(w(c),'channel');;
        plot(easting(c),northing(c),'o','MarkerFaceColor',cols(c),'MarkerSize',10)
        hold on
        quiver(easting(c),northing(c),-easting(c)/100,-northing(c)/100,0); % /100 just gives arrow length
        text(easting(c)+1,northing(c),chan(1:3));
    end
    grid on
    quiver(440,1325,wind_speed*sin(deg2rad(wind_direction)), wind_speed*cos(deg2rad(wind_direction)) ,0,'k');
    text(440,1325,'wind')
    hold off
    title('Beach House array position relative to SLC40');
    xlabel('metres east');
    ylabel('metres north');
    axis equal;
    outfile = sprintf('%s/arraymap.png',figureOutDirectory);
    feval('print', '-dpng', outfile);
    close
end

author: geoscience-community-codes

repo: GISMO


text:
function H= plot_scalpEvolution(erp, mnt, ival, varargin)
%PLOT_SCALPEVOLUTION - Draws scalp topographies for specified intervals
%
%Synposis:
% H= plot_scalpEvolution(ERP, MNT, IVAL, <OPTS>)
%
%Description:
% Draws scalp topographies for specified intervals,
% separately for each each Class. For each Classes topographies are
% plotted in one row and shared the same Color map scaling. (In future
% versions there might be different options for Color scaling.)
%
%Input:
% ERP: struct of epoched EEG data.
% MNT: struct defining an electrode montage
% IVAL: [nIvals x 2]-sized array of interval, which are marked in the
%       ERP plot and for which scalp topographies are drawn.
%       When all interval are consequtive, ival can also be a
%       vector of interval borders.
% OPTS: struct or property/value list of optional fields/properties:
%  .IvalColor - [nColors x 3]-sized array of rgb-coded Colors
%                with are used to mark intervals and corresponding
%                scalps. Colors are cycled, i.e., there need not be
%                as many Colors as interval. Two are enough,
%                default [0.6 1 1; 1 0.6 1].
%  .LegendPos - specifies the position of the legend in the ERP plot,
%               default 0 (see help of legend for choices).
%
%The opts struct is passed to plot_scalpPattern
%
%Output:
% H: struct of handles to the created graphic objects.
%
%See also plot_scalpEvolutionPlusChannel, plot_scalpPatterns, plot_scalp.

% 01-2005 Benjamin Blankertz

props= {'PrintIval',    0,  'BOOL';
        'PlotChannel',  0,  'BOOL'};
props_scalpEvolutionPlusChannel= plot_scalpEvolutionPlusChannel;

if nargin==0,
  H= opt_catProps(props, props_scalpEvolutionPlusChannel);
  return
end

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

opt_scalpEvolutionPlusChannel= ...
    opt_substruct(opt, props_scalpEvolutionPlusChannel(:,1));

H= plot_scalpEvolutionPlusChannel(erp, mnt, [], ival, ...
                                  opt_scalpEvolutionPlusChannel);

if nargout<1,
  clear H
end

author: bbci


text:
function grad = B_SpatialCovMask(future_layers, prev_layers, curr_layer)
mask = prev_layers{1}.a;
data = prev_layers{2}.a;
[D,T,N] = size(mask);
[D2,T,N] = size(data);
nCh = D2/D;
data = reshape(data, D, nCh, T, N);
data = permute(data, [2 3 1 4]);
% data = abs(data);

mask_f = sum(mask,2);     % sum of mask over time

spatCov = curr_layer.a;
speechCov = reshape(spatCov(1:D*nCh^2,:,:,:), nCh^2,1,D,N);
noiseCov = reshape(spatCov(D*nCh^2+1:end,:,:,:), nCh^2,1,D,N);

future_grad = GetFutureGrad(future_layers, curr_layer);
future_grad_speech = reshape(future_grad(1:D*nCh^2,:,:), nCh^2,1,D,N);
future_grad_noise = reshape(future_grad(D*nCh^2+1:end,:,:), nCh^2,1,D,N);

% from speech covariance
if 0    % slow implementation
    data = gather(data);
    data_cell = num2cell(data, [1]);
    tic; xx = cellfun(@(x) (reshape(x*x',nCh^2,1)), data_cell, 'UniformOutput', 0);toc;
    xx2 = cell2mat(xx);
else    % really fast implementation
    xx2 = reshape(outProdND(data), nCh^2, T,D);
end

gradFromSpeech = bsxfun(@minus, xx2, speechCov);
gradFromSpeech = bsxfun(@times, gradFromSpeech, future_grad_speech);
gradFromSpeech = squeeze(sum(gradFromSpeech))';
gradFromSpeech = bsxfun(@times, gradFromSpeech, 1./mask_f);

gradFromNoise = bsxfun(@minus, xx2, noiseCov);
gradFromNoise = bsxfun(@times, gradFromNoise, future_grad_noise);
gradFromNoise = squeeze(sum(gradFromNoise))';
gradFromNoise = bsxfun(@times, gradFromNoise, 1./(T-mask_f));

grad = gradFromSpeech - gradFromNoise;
% we may get complex valued gradient for real valued mask, as the partial
% derivative formula does not know that the mask must be real valued. To
% preserve the real-valued property of the mask, we just take the real part
% of the gradient, which represents the search direction of mask in the
% real domain.
grad = real(grad);

end

author: singaxiong

repo: SignalGraph


text:
%%%% Create animation for a 2-link serial chain manipiulator

% Course: Robotic Manipulation and Mobility
% Advisor: Dr. V. Krovi
%
% Homework Number: MIDTERM
%
% Names: Sourish Chakravarty
% 	Hrishi Lalit Shah

function [aviobj]= CREATE_ANIME_1(aviobj,T,th1,x2,y2,th2,h)

global l1 lc1 l2 lc2
figure(h)
Z0=[0 0]; %%% Point of suspension
for i=1:length(T)

    c1=cos(th1(i));
    c2=cos(th2(i));
    s1=sin(th1(i));
    s2=sin(th2(i));
    %%% Plotting Link - 1
    Z1= Z0 + [l1*c1, l1*s1];% Coordinate of the hanging end link 1
    plot([Z0(1),Z1(1)],[Z0(2),Z1(2)],'b','linewidth',4);
    hold on
    Z1m= Z0 + [lc1*c1, lc1*s1];% Coordinate of the CM of link 1
    plot(Z1m(1), Z1m(2),'r*'); %Plots mid point
    plot(Z0(1),Z0(2),'k*'); % Plots point of suspension

    %%% Plotting Link - 2
    Z2m = [x2(i), y2(i)];% Coordinate of the CM of link 2
    Z2l = Z2m - [lc2*c2, lc2*s2]; % Coordinate of source-end of link - 2
    Z2r = Z2m + [(l2-lc2)*c2, (l2-lc2)*s2]; % Coordinate of effector-end of link-2
    plot([Z2l(1),Z2r(1)],[Z2l(2),Z2r(2)],'m','linewidth',4);
    plot(Z2m(1), Z2m(2),'r*'); %Plots mid point
    xlim([-4,4]);
    ylim([-4,4]);
%     title('Midterm Animation');
    xlabel('X-axis');
    ylabel('Y-axis');
    grid on
    hold off
    pause(0.01);     %Stop execution for 0.01 to make animation visible
    frame= getframe(gcf);   %Step 2: Grab the frame
    aviobj = addframe(aviobj,frame); % Step 3: Add frame to avi object
end


text:
function pass = test_integral( )
% Test INTEGRAL()

tol = 1000*chebfunpref().cheb2Prefs.chebfun2eps;

%test empty diskfun
f = diskfun();
pass(1) = (integral(f) == 0);
pass(2) = (integral(f, chebfun(@(x) x)) == 0);

%test unitcircle feature
f = diskfun(@(x,y) sin(x.^2+3*y));
pass(3) = ( abs(sum(f(:, 1))-integral(f, 'unitcircle')) < tol);
pass(4) = ( abs(integral(diskfun.harmonic(3,2, 'neumann'), 'unitcircle')) < tol);

% test integral along a contour (need more rigorous tests when complex-valued
% diskfuns are supported)
z = chebfun(@(x) .5*exp(1i*pi*x));
g = diskfun(@(x, y) exp(-2*x.^2-2*y.^2));
pass(5) = ( abs(.5*sum(g(:,.5))-integral(g, z)) < tol);
z = chebfun(@(x) x*exp(1i*pi/4));
pass(6) = ( abs(sum(g(pi/4, :))-integral(g, z) )  < tol );
end


text:
% visualizeHKS
clear all; close all; clc;

gptoolbox_path = ''; % complete
addpath(genpath(gptoolbox_path));

data_folder = '..\data\';
meshes = dir([data_folder,'*.o*']); % a list of all the meshes in the data folder
meshes = {meshes.name}; % access each mesh name by meshes{i}

cam_folder = [data_folder, 'cams\'];
cams = dir([cam_folder,'*.mat']); % a list of all the cam in the cams folder
cams = {cams.name};

addpath(data_folder); addpath(cam_folder);

%% Compute the HKS using gptoolbox

%% Display the HKS on a few meshes for a small fixed time t.

author: odedstein

repo: sgi-introduction-course


text:
% eegplugin_iirfilt() - EEGLAB plugin for importing data using IIRFILT Matlab toolbox
%
% Usage:
%   >> eegplugin_iirfilt(fig, trystrs, catchstrs);
%
% Inputs:
%   fig        - [integer]  EEGLAB figure
%   trystrs    - [struct] "try" strings for menu callbacks.
%   catchstrs  - [struct] "catch" strings for menu callbacks.
%
% Create a plugin:
%   For more information on how to create an EEGLAB plugin see the
%   help message of eegplugin_besa() or visit http://www.sccn.ucsd.edu/eeglab/contrib.html
%
% Authors: Maksym Pozdin (mpozdin.ece04@gtalumni.org, IOL/ONRC,2004),
%          with Arnaud Delorme and Scott Makeig (SCCN/INC/UCSD, La Jolla CA)

% Copyright (C) 2003 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function vers = eegplugin_iirfilt(fig, trystrs, catchstrs)

    vers = 'iirfilt1.01';
  if nargin < 3
      error('eegplugin_iirfilt requires 3 arguments');
  end;

  % add besa folder to path
  % -----------------------
  if ~exist('pop_iirfilt')
      p = which('eegplugin_iirfilt');
      p = p(1:findstr(p,'eegplugin_iirfilt.m')-1);
      if ~isdeployed
        addpath([ p vers ] ); end
  end;

  % find import data menu
  % ---------------------
  menu = findobj(fig, 'tag', 'filter');

  % menu callbacks
  % --------------
  combio = [ trystrs.no_check '[EEG LASTCOM] = pop_iirfilt(EEG);' catchstrs.new_and_hist ];

  % create menus if necessary
  % -------------------------
  uimenu( menu, 'Label', 'Short IIR filter',  'CallBack', combio, 'Separator', 'on');

author: PatternRecognition

repo: OpenBMI


text:
function writePolygonSet(polys, filename)
%WRITEPOLYGONSET Write a set of simple polygons into a file.
%
%   writePolygonSet(POLYS, FILENAME);
%   Writes the set of polygons in the file FILENAME.
%   Following format is used:
%     X11 X12 X13 ... X1N
%     Y11 Y12 Y13 ... Y1N
%     X21 X22 X23 ... X2N
%     Y21 Y22 Y23 ... Y2N
%   Each polygon may have a different number of vertices.
%
%   See also
%   polygons2d, readPolygonSet

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2013-01-14
% Copyright 2013-2022 INRA - TPV URPOI - BIA IMASTE

% open file for reading
fid = fopen(filename, 'wt');

for i = 1:length(polys)
    poly = polys{i};
    n = size(poly, 1);

    % precompute format
    format = [repmat('%g ', 1, n) '\n'];

    % write one line for x, then one line for y
    fprintf(fid, format, poly(:,1)');
    fprintf(fid, format, poly(:,2)');
end

% close file
fclose(fid);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


% sensors.m
%   Compute the output of rate gyros, accelerometers, and pressure sensors
%
%  Revised:
%   3/5/2010  - RB
%   5/14/2010 - RB

function y = sensors(uu, P)

    % relabel the inputs
%    pn      = uu(1);
%    pe      = uu(2);
    pd      = uu(3);
%    u       = uu(4);
%    v       = uu(5);
%    w       = uu(6);
    phi     = uu(7);
    theta   = uu(8);
%    psi     = uu(9);
    p       = uu(10);
    q       = uu(11);
    r       = uu(12);
    F_x     = uu(13);
    F_y     = uu(14);
    F_z     = uu(15);
%    M_l     = uu(16);
%    M_m     = uu(17);
%    M_n     = uu(18);
    Va      = uu(19);
%    alpha   = uu(20);
%    beta    = uu(21);
%    wn      = uu(22);
%    we      = uu(23);
%    wd      = uu(24);

    % simulate rate gyros (units are rad/sec)
    y_gyro_x = p + P.sigma_gyro * randn;
    y_gyro_y = q + P.sigma_gyro * randn;
    y_gyro_z = r + P.sigma_gyro * randn;

    % simulate accelerometers (units of g)
    y_accel_x = F_x / P.mass / P.gravity + sin(theta) + P.sigma_accel * randn;
    y_accel_y = F_y / P.mass / P.gravity - cos(theta) * sin(phi) + P.sigma_accel * randn;
    y_accel_z = F_z / P.mass / P.gravity - cos(theta) * cos(phi) + P.sigma_accel * randn;

    % simulate pressure sensors
    y_static_pres = P.rho * P.gravity * (-pd) + P.beta_abs_pres + P.sigma_abs_pres * randn;
    y_diff_pres = P.rho * Va^2 / 2 + P.beta_diff_pres + P.sigma_diff_pres * randn;

    % construct total output
    y = [...
        y_gyro_x;...
        y_gyro_y;...
        y_gyro_z;...
        y_accel_x;...
        y_accel_y;...
        y_accel_z;...
        y_static_pres;...
        y_diff_pres;...
    ];

end


text:
function S=generateSAbsolute( anim )
% Generate S tranformed by R and t
%
% USAGE
%  anim = anim.generateSAbsolute()
%
% INPUTS
%  anim     - Animation object (help Animation for details)
%
% OUTPUTS
%  S        - anim.S after R, t have been applied (or homographies)
%
% EXAMPLE
%
% See also GENERATETOYANIMATION
%
% Vincent's Structure From Motion Toolbox      Version 3.0
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

if isempty(anim.S)
  warning('S is not defined nor computable'); S=[]; return;
end
if (isempty(anim.R) || isempty(anim.t)) && (isempty(anim.P))
  warning('Need to define R,t or P'); S=[]; return;
end

if size(anim.P,2)==3
  % homography case
  if size(anim.S,3)==1
    S = multiTimes( anim.P, anim.S, 1 );
  else
    S = multiTimes( anim.P, anim.S, 2 );
  end
else
  % compute the projected points for normal camera and 3D points
  if isempty(anim.R)
    if size(anim.S,3)==1
      S = bsxfun( @plus, multiTimes( anim.P(:,1:3,:), anim.S, 1 ), ...
        anim.P(:,4,:) );
    else
      S = bsxfun( @plus, multiTimes( anim.P(:,1:3,:), anim.S, 2 ), ...
        anim.P(:,4,:) );
    end
  else
    if size(anim.S,3)==1
      S = bsxfun( @plus, multiTimes( anim.R, anim.S, 1 ), ...
        reshape(anim.t, 3, 1, anim.nFrame ) );
    else
      S = bsxfun( @plus, multiTimes( anim.R, anim.S, 2 ), ...
        reshape(anim.t, 3, 1, anim.nFrame ) );
    end
  end
end

author: vrabaud


text:
function [ joints ] = get_anno_joints( rect, pidxs, parts )

num_joints = length(pidxs);
joints = NaN(num_joints, 2);

points = rect.annopoints.point;
for j = 1:num_joints
    pidx = pidxs(j);
    annopoint_idxs = parts(pidx+1).pos;
    assert(annopoint_idxs(1) == annopoint_idxs(2));
    pt = util_get_annopoint_by_id(points, annopoint_idxs(1));
    if (~isempty(pt))
        joints(j, :) = [pt.x pt.y];
    end
end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [X,Y,Z]=contourRegularise(varargin)

switch nargin
    case 2
        Vcs=varargin{1};
        np=varargin{2};
        interpMethod='pchip'; %Default method
    case 3
        Vcs=varargin{1};
        np=varargin{2};
        interpMethod=varargin{3};
    otherwise
        error('Wrong number of input arguments');
end

%%
numContours=numel(Vcs);

%Allocate coordinate matrices
X=nan(numContours,np); Y=nan(numContours,np); Z=nan(numContours,np);
startDefined=0;
for q=1:1:numContours

    %Join goups if present N.B. assumes they belong to the same curve!
    Vs=[];
    for qGroup=1:1:numel(Vcs{q})
        Vss=Vcs{q}{qGroup};
        if ~isempty(Vss)
            if isPolyClockwise(Vss)
                Vss=flipud(Vss);
            end
            Vs=[Vs; Vss];
        end
    end

    %Resample curve so they have the same number of points
    if ~isempty(Vs)
        numDigitsMerge=6-numOrder(mean(sum(diff(Vs,[],1).^2,2)));
        [~,ind1,~]=unique(pround(Vs,numDigitsMerge),'rows');

        Vs=Vs(ismember(1:size(Vs,1),ind1),:); %This maintains point order

        [Vs]=evenlySampleCurve(Vs,np,interpMethod,1);

        if startDefined==0 %After 1 contour has been added
            V_start=Vs(1,:);
            startDefined=1;
        else
            %Reorder curves so start and end points match closely
            D=sqrt(sum((Vs-V_start(ones(1,size(Vs,1),1),:)).^2,2));
            [~,indMin]=min(D);
            if indMin~=1
                Vs=[Vs(indMin:size(Vs,1),:); Vs(1:indMin-1,:)];
            end
            V_start=Vs(1,:);
        end

        %Store coordinates
        [X(q,:),Y(q,:),Z(q,:)]=getColumns(Vs);
    end
end

%%
%Remove twist
for q=2:1:size(numContours,1)
    v1=[X(q-1,:)' Y(q-1,:)' Z(q-1,:)'];
    v2=[X(q,:)' Y(q,:)' Z(q,:)'];
    [v2f,~,~]=minPolyTwist(v1,v2);
    X(q,:)=v2f(:,1);
    Y(q,:)=v2f(:,2);
    Z(q,:)=v2f(:,3);
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [ meta ] = read_nitf_meta( filename )
%READ_NITF_META Read metadata from NITF into a structure
%
%   NITFHDR = READ_NITF_META(NITFFILE) calls functions to parse NITF file
%   segments.
%
% Written by: Matt Donath, NGA, matthew.b.donath@nga.ic.gov
% Modified by: Wade Schwartzkopf, NGA/IDT
% References:
%  > MIL-STD-2500A, NATIONAL IMAGERY TRANSMISSION FORMAT VERSION 2.0
%  > MIL-STD-2500C, NATIONAL IMAGERY TRANSMISSION FORMAT VERSION 2.1
%  > STDI-0001,     NATIONAL SUPPORT DATA EXTENSIONS (SDE)(VERSION 1.3/CN2)
%                   FOR THE NATIONAL IMAGERY TRANSMISSION FORMAT (NITF)
%  > STDI-0002,     THE COMPENDIUM OF CONTROLLED EXTENSIONS (CE) FOR THE
%                   NATIONAL IMAGERY TRANSMISSION FORMAT (NITF) VERSION 2.1
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% OPEN FILE
fid=fopen(filename,'r','b'); % NITF always big-endian

%% READ NITF FILE HEADER
meta.filehdr = read_nitf_filehdr(fid);

%% READ IMAGE SEGMENTS
nextImgSubhdroffset=meta.filehdr.HL;
for i = 1:meta.filehdr.NUMI
    fseek(fid,nextImgSubhdroffset,'bof'); % Jump to ith image segment
    meta.imagesubhdr{i}=read_nitf_imgsubhdr(fid);
    nextImgSubhdroffset=nextImgSubhdroffset+meta.filehdr.LISH(i)+meta.filehdr.LI(i);
end
fseek(fid,nextImgSubhdroffset,'bof'); % Jump to end of image segments

%% SKIP GRAPHIC SEGMENTS
if isfield(meta.filehdr,'NUMS') && meta.filehdr.NUMS > 0
    for seg = 1:meta.filehdr.NUMS
        fseek(fid, meta.filehdr.LSSH(seg), 'cof'); % Skip graphic segment subheaders
        fseek(fid, meta.filehdr.LS(seg), 'cof'); % Skip graphic segment data
    end
end

%% SKIP RESERVED SEGMENTS
if isfield(meta.filehdr,'NUMX') && meta.filehdr.NUMX > 0
    for seg = 1:meta.filehdr.NUMX
        fseek(fid, meta.filehdr.LXSH(seg), 'cof'); % Skip reserved segment subheaders
        fseek(fid, meta.filehdr.LX(seg), 'cof'); % Skip reserved segment data
    end
end

%% SKIP LABEL SEGMENTS
if isfield(meta.filehdr,'NUML') && meta.filehdr.NUML > 0
    for seg = 1:meta.filehdr.NUML
        fseek(fid, meta.filehdr.LLSH(seg), 'cof'); % Skip label segment subheaders
        fseek(fid, meta.filehdr.LL(seg), 'cof'); % Skip label segment data
    end
end

%% SKIP TEXT SEGMENTS
if isfield(meta.filehdr,'NUMT') && meta.filehdr.NUMT > 0
    for seg = 1:meta.filehdr.NUMT
%         % Skip text segments
%         fseek(fid, meta.filehdr.LTSH(seg), 'cof'); % Skip text segment subheaders
%         fseek(fid, meta.filehdr.LT(seg), 'cof'); % Skip text segment data
        % Read text segments
        meta.txtsubhdr{seg}=read_nitf_txtsubhdr(fid);
        meta.txtseg{seg} = fread(fid,meta.filehdr.LT(seg),'uint8=>char')';
    end
end

%% READ DATA EXTENSION SEGMENT SUBHEADERS
curr_pos = ftell(fid);
if isfield(meta.filehdr,'NUMDES') && meta.filehdr.NUMDES > 0
    for seg = 1:meta.filehdr.NUMDES
        fseek(fid, curr_pos, 'bof');
        fld = strcat('DESSubhdr',num2str(seg));
        meta = readDESegSubhdr(meta,fid,fld);
        curr_pos = curr_pos + meta.filehdr.LDSH(seg) + meta.filehdr.LD(seg);
    end
end

%% READ RESERVED EXTENSION SUBHEADERS
if isfield(meta.filehdr,'NUMRES') && meta.filehdr.NUMRES > 0
    for seg = 1:meta.filehdr.NUMRES
        fld = strcat('RESSubhdr',num2str(seg));
        meta = readRESegSubhdr(meta,fid,fld);
    end
end

%% Close file
fclose(fid);

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function generateOutput(dataTree)
global UNIT_TEST

if dataTree.IsEmpty()
    return
end

UNIT_TEST = true; %#ok<NASGU>
iG = dataTree.GetCurrElemIndexID();
banner = sprintf('Re-calculating derived data at %s with the following processing stream:\n\n', char(datetime(datetime, 'Format','HH:mm:ss, MMMM d, yyyy')));
dataTree.PrintProcStream(banner);
dataTree.groups(iG).Calc();
dataTree.groups(iG).Save();
UNIT_TEST = false;

author: BUNPC


text:
function [intersects, edgeIndices] = intersectLinePolyline(line, poly, varargin)
%INTERSECTLINEPOLYLINE Intersection points between a line and a polyline.
%
%   P = intersectLinePolyline(LINE, POLY)
%   Returns the intersection points of the lines LINE with polyline POLY.
%   LINE is a 1-by-4 row vector containing parametric representation of the
%   line (in the format [x0 y0 dx dy], see the function 'createLine' for
%   details).
%   POLY is a NV-by-2 array containing coordinates of the polyline vertices
%   P is a K-by-2 array containing the coordinates of the K intersection
%   points.
%
%   P = intersectLinePolyline(LINE, POLY, TOL)
%   Specifies the tolerance for geometric tests. Default is 1e-14.
%
%   [P INDS] = intersectLinePolyline(...)
%   Also returns the indices of edges involved in intersections. INDS is a
%   K-by-1 column vector, such that P(i,:) corresponds to intersection of
%   the line with the i-th edge of the polyline. If the intersection occurs
%   at a polyline vertex, the index of only one of the two neighbor edges
%   is returned.
%   Note that due to numerical approximations, the use of function
%   'isPointOnEdge' may give results not consistent with this function.
%
%
%   Examples
%   % compute intersections between a square and an horizontal line
%     poly = [0 0;10 0;10 10;0 10];
%     line = [5 5 1 0];
%     intersectLinePolyline(line, poly)
%     ans =
%           10     5
%     % also return indices of edges
%     [inters inds] = intersectLinePolyline(line, poly)
%     inters =
%           10     5
%     inds =
%           2
%
%   % compute intersections between a square and a diagonal line
%     poly = [0 0;10 0;10 10;0 10];
%     line = [5 5 1 1];
%     intersectLinePolyline(line, poly)
%     ans =
%            0     0
%           10    10
%
%   See also
%   lines2d, polylines2d, intersectLines, intersectLinePolygon
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2003-10-31
% Copyright 2003-2022 INRA - TPV URPOI - BIA IMASTE

% get computation tolerance
tol = 1e-14;
if ~isempty(varargin)
    tol = varargin{1};
end

% create the array of edges
N = size(poly, 1);
edges = [poly(1:N-1, :) poly(2:N, :)];

% compute intersections with supporting lines of polyline edges
supportLines = edgeToLine(edges);
intersects = intersectLines(line, supportLines, tol);

% find edges that are not parallel to the input line
inds = find(isfinite(intersects(:, 1)));

% compute position of intersection points on corresponding lines
pos = linePosition(intersects(inds, :), supportLines(inds, :), 'diag');

% and keep only intersection points located on edges
b = pos > -tol & pos < 1+tol;
inds = inds(b);
intersects = intersects(inds, :);

% remove multiple vertices (can occur for intersections located at polyline
% vertices)
[intersects, I, J] = unique(intersects, 'rows'); %#ok<ASGLU>

if nargout > 1
    % return indices of edges involved in intersection
    % (in case of intersection located at a vertex, only one of the
    % neighbor edges is returned)
    edgeIndices = inds(I);
end


text:
function [params, s_new] = denoise_fica_tanh(params, s, state)
% FastICA tanh nonlinearity as DSS denoising function
%   [params, s_new] = denoise_fica_tanh(params, s, state)
%     params  Function specific modifiable parameters
%       a     Scaling constant (default: 1)
%     state   DSS algorithm state
%     s       Source signal estimate, matrix of row vector signals
%     s_new   Denoised signal estimate

% Copyright (C) 2004, 2005 DSS MATLAB package team (dss@cis.hut.fi).
% Distributed by Laboratory of Computer and Information Science,
% Helsinki University of Technology. http://www.cis.hut.fi/projects/dss/.
% $Id$

if nargin<3 | ~isstruct(state)
    params.name = 'FastICA tanh nonlinearity';
    params.description = '';
    params.param = {'a'};
    params.param_value ={1};
    params.param_type = {'scalar'};
    params.param_desc = {'Steepness modifier'};
    params.approach = {'defl','symm'};
    params.alpha = {};
    params.beta = {};
    return;
end

if ~isfield(params, 'a')
  params.a = 1;
end
a = params.a;

hypTan = tanh(a * s);
s_new = hypTan - a * repmat(mean(1 - hypTan .^ 2, 2), 1, length(s)) .* s;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function b = r8blt_vxm ( n, ml, a, x )

%*****************************************************************************80
%
%% R8BLT_VXM multiplies a vector by a R8BLT matrix.
%
%  Discussion:
%
%    The R8BLT storage format is appropriate for a banded lower triangular matrix.
%    The matrix is assumed to be zero below the ML-th subdiagonal.
%    The matrix is stored in an ML+1 by N array, in which the diagonal
%    appears in the first row, followed by successive subdiagonals.
%    Columns are preserved.
%
%  Example:
%
%    N = 5, ML = 2
%
%    A11   0   0   0   0
%    A21 A22   0   0   0
%    A31 A32 A33   0   0
%      0 A42 A43 A44   0
%      0   0 A53 A54 A55
%                --- ---
%                    ---
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer ML, the lower bandwidth.
%
%    Input, real A(ML+1,N), the R8BLT matrix.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Output, real B(N), the product X*A.
%
  b(1:n) = 0.0;

  for i = 1 : n
    jlo = max ( 1, i - ml );
    jhi = i;
    for j = jlo : jhi
      b(j) = b(j) + x(i) * a(i-j+1,j);
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% check latlon2local

origin = [30.5298813668333; 114.350355321667; 20.8869];
lla = [30.5298266666667; 114.350367833333; 20.8];
lla = [30.5299376666667; 114.35528; 17.6];
% matlab functions
[e1, n1, u1] = geodetic2enu(lla(1), lla(2), lla(3), origin(1),...
    origin(2), origin(3), wgs84Ellipsoid, 'degrees');
[e2, n2, u2] = latlon2local(lla(1), lla(2), lla(3), origin);

assert(abs(e1 - e2) < 1e-8);
assert(abs(n1 - n2) < 1e-8);
assert(abs(u1 - u2) < 1e-8);

% our own implementation
Ce2n0=llh2dcm_v000(origin(1:2) * pi/ 180,[0;1]);
qs02e=rotro2qr(Ce2n0');
origin_xyz = ecef2geo_v000([origin(1:2) * pi / 180; origin(3)],1);
xyz = ecef2geo_v000([lla(1:2) * pi / 180; lla(3)],1);
ned = quatrot_v000(qs02e, xyz - origin_xyz ,1)';

assert(abs(e1 - ned(2)) < 1e-8);
assert(abs(n1 - ned(1)) < 1e-8);
assert(abs(u1 + ned(3)) < 1e-8);

repo: ekfmonoslam


text:
function [gps_week, gps_sow, gps_dow] = jd2gps(jd)

% SYNTAX:
%   [gps_week, gps_sow, gps_dow] = jd2gps(jd);
%
% INPUT:
%   jd = julian day
%
% OUTPUT:
%   gps_week = GPS week
%   gps_sow  = GPS seconds of week
%   gps_dow  = GPS day of week
%
% DESCRIPTION:
%   Conversion of julian day number to GPS week and
%	seconds of week.

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

deltat = jd - 2444244.5;
gps_week = floor(deltat/7);
gps_dow  = floor(deltat - gps_week*7);
gps_sow  = (deltat - gps_week*7)*86400;

author: goGPS-Project


text:
function lockfiles = check_for_lock_files(results_directory)
%Uses the unix find command to see if there are any lock files
%(directories which end in ".lock") inside the directory specified
%by results_directory
%
%Returns a cell array of files, potentially of length 0
%Tomasz Malisiewicz (tomasz@csail.mit.edu)

[a,b] = unix(sprintf('find %s -name "*.lock"',results_directory));
if a==1 || length(b) == 0
  lockfiles = {};
else
  b=textscan(b,'%s');
  lockfiles = b{1};
end

author: quantombone

repo: exemplarsvm


text:
% Test file for scribble.m.

function pass = test_scribble(pref)

if ( nargin < 1 )
    pref = chebfunpref();
end

% Save the current warning state and then clear it in case scribble() issued a
% warning previously.
[lastmsg, lastid] = lastwarn();
lastwarn('');

% To pass, this call to scribble() must not crash or generate any warnings.
c = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789.,;:?!''"`_ )([]{}-+*/^=<>\|%#~@$&';
s = scribble(c);

[wmsg, wid] = lastwarn();
if ( strcmp(wid, 'CHEBFUN:scribble:unknownchar') )
    pass(1) = false;
else
    % If this run of scribble() didn't issue a warning, pass and restore the
    % warning state.
    pass(1) = true;
    lastwarn(lastmsg, lastid);
end

end

author: chebfun

repo: chebfun


text:
function [d pred] = dag_sp(A,u,varargin)
% DAG_SP Compute the weighted single source shortest path problem.
%
% The DAG shortest path algorithm for the single source shortest path
% problem only works on directed acyclic-graphs (DAGs).
%
% If the graph is not a DAG, the results are undefined.  In the future, the
% function may throw an error if the graph is not a DAG.
%
% See the shortest_paths function for calling information.  This function
% just calls shortest_paths(...,struct('algname','dag'));
%
% This algorithm works on weighted directed acyclic graphs.
% The runtime is O(V+E)
%
% ... = clustering_coefficients(A,...) takes a set of
% key-value pairs or an options structure.  See set_matlab_bgl_options
% for the standard options.
%    There are no additional options for this function.
%
% Example:
%    load graphs/kt-3-7.mat
%    dag_sp(A,1)
%
% See also SHORTEST_PATHS

% David Gleich
% Copyright, Stanford University, 2006-2008

%% History
%  2006-04-23: Initial version
%  2008-10-07: Changed options parsing
%%

algname = 'dag';
if ~isempty(varargin),
    options = merge_options(struct(),varargin{:});
    options.algname= algname;
else options = struct('algname',algname);
end

[d pred] = shortest_paths(A,u,options);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore