text:
function [source] = ft_datatype_source(source, varargin)

% FT_DATATYPE_SOURCE describes the FieldTrip MATLAB structure for data that is
% represented at the source level. This is typically obtained with a beamformer of
% minimum-norm source reconstruction using FT_SOURCEANALYSIS.
%
% An example of a source structure obtained after performing DICS (a frequency domain
% beamformer scan) is shown here
%
%           pos: [6732x3 double]       positions at which the source activity could have been estimated
%        inside: [6732x1 logical]      boolean vector that indicates at which positions the source activity was estimated
%           dim: [xdim ydim zdim]      if the positions can be described as a 3D regular grid, this contains the
%                                       dimensionality of the 3D volume
%     cumtapcnt: [120x1 double]        information about the number of tapers per original trial
%          time: 0.100                 the latency at which the activity is estimated (in seconds)
%          freq: 30                    the frequency at which the activity is estimated (in Hz)
%           pow: [6732x120 double]     the estimated power at each source position
%     powdimord: 'pos_rpt'             defines how the numeric data has to be interpreted,
%                                       in this case 6732 dipole positions x 120 repetitions (i.e. trials)
%           cfg: [1x1 struct]          the configuration used by the function that generated this data structure
%
% Required fields:
%   - pos
%
% Optional fields:
%   - inside, pow, coh, eta, mom, ori, leadfield, filter, or any other field with dimensions that are consistent with pos or dim
%   - dim, transform, unit, coordsys, time, freq, cumtapcnt, dimord
%
% Deprecated fields:
%   - method, outside
%
% Obsoleted fields:
%   - xgrid, ygrid, zgrid, transform, latency, frequency
%
% Revision history:
%
% (2014) The subfields in the avg and trial fields are now present in the
% main structure, e.g. source.avg.pow is now source.pow. Furthermore, the
% inside is always represented as logical vector.
%
% (2011) The source representation should always be irregular, i.e. not
% a 3-D volume, contain a "pos" field and not contain a "transform".
%
% (2010) The source structure should contain a general "dimord" or specific
% dimords for each of the fields. The source reconstruction in the avg and
% trial substructures has been moved to the toplevel.
%
% (2007) The xgrid/ygrid/zgrid fields have been removed, because they are
% redundant.
%
% (2003) The initial version was defined
%
% See also FT_DATATYPE, FT_DATATYPE_COMP, FT_DATATYPE_DIP, FT_DATATYPE_FREQ,
% FT_DATATYPE_MVAR, FT_DATATYPE_RAW, FT_DATATYPE_SOURCE, FT_DATATYPE_SPIKE,
% FT_DATATYPE_TIMELOCK, FT_DATATYPE_VOLUME

% Copyright (C) 2013-2014, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% FIXME: I am not sure whether the removal of the xgrid/ygrid/zgrid fields
% was really in 2007

% get the optional input arguments, which should be specified as key-value pairs
version = ft_getopt(varargin, 'version', 'latest');

if strcmp(version, 'latest') || strcmp(version, 'upcoming')
  version = '2014';
end

if isempty(source)
  return;
end

% old data structures may use latency/frequency instead of time/freq. It is
% unclear when these were introduced and removed again, but they were never
% used by any FieldTrip function itself.
if isfield(source, 'frequency')
  source.freq = source.frequency;
  source      = rmfield(source, 'frequency');
end
if isfield(source, 'latency')
  source.time = source.latency;
  source      = rmfield(source, 'latency');
end

switch version
  case '2014'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % ensure that it has individual source positions
    source = fixpos(source);

    if isfield(source, 'inside')
      % ensure that it is always logical
      source = fixinside(source, 'logical');
    end

    if isfield(source, 'coordsys')
      % ensure that it is in lower case
      source.coordsys = lower(source.coordsys);
    end

    if isfield(source, 'unit')
      % ensure that it is in lower case
      source.unit = lower(source.unit);
    end

    % remove obsolete fields
    if isfield(source, 'method')
      source = rmfield(source, 'method');
    end
    if isfield(source, 'transform')
      source = rmfield(source, 'transform');
    end
    if isfield(source, 'xgrid')
      source = rmfield(source, 'xgrid');
    end
    if isfield(source, 'ygrid')
      source = rmfield(source, 'ygrid');
    end
    if isfield(source, 'zgrid')
      source = rmfield(source, 'zgrid');
    end

    if isfield(source, 'avg') && isstruct(source.avg) && isfield(source, 'trial') && isstruct(source.trial) && ~isempty(intersect(fieldnames(source.avg), fieldnames(source.trial)))
      % it is not possible to convert both since they have the same field names
      ft_warning('removing ''avg'', keeping ''trial''');
      source = rmfield(source, 'avg');
    end

    if isfield(source, 'avg') && isstruct(source.avg)
      % move the average fields to the main structure
      fn = fieldnames(source.avg);
      for i=1:length(fn)
        dat = source.avg.(fn{i});
        if isequal(size(dat), [1 size(source.pos,1)])
          source.(fn{i}) = dat';
        else
          source.(fn{i}) = dat;
        end
        clear dat
      end % j
      source = rmfield(source, 'avg');
    end

    if isfield(source, 'trial') && isstruct(source.trial)
      npos = size(source.pos,1);

      % concatenate the fields for each trial and move them to the main structure
      fn = fieldnames(source.trial);

      for i=1:length(fn)
        % some fields are descriptive and hence identical over trials
        if any(strcmp(fn{i}, {'csdlabel' 'label' 'filterdimord' 'leadfielddimord'}))
          source.(fn{i}) = source.trial(1).(fn{i});
          continue
        end

        % start with the first trial
        dat    = source.trial(1).(fn{i});
        datsiz = getdimsiz(source, fn{i});
        nrpt   = datsiz(1);
        datsiz = datsiz(2:end);

        if iscell(dat)
          datsiz(1) = nrpt; % swap the size of pos with the size of rpt
          val  = cell(npos,1);
          if isfield(source, 'inside')
            indx = find(source.inside);
          else
            indx = 1:npos;
          end
          for k=1:length(indx)
            val{indx(k)}          = nan(datsiz);
            val{indx(k)}(1,:,:,:) = dat{indx(k)};
          end
          % concatenate all data as {pos}_rpt_etc
          for j=2:nrpt
            dat = source.trial(j).(fn{i});
            for k=1:length(indx)
              val{indx(k)}(j,:,:,:) = dat{indx(k)};
            end

          end % for all trials
          source.(fn{i}) = val;

        else
          % concatenate all data as pos_rpt_etc
          val = nan([datsiz(1) nrpt datsiz(2:end)]);
          val(:,1,:,:,:) = dat(:,:,:,:);
          for j=2:length(source.trial)
            dat = source.trial(j).(fn{i});
            val(:,j,:,:,:) = dat(:,:,:,:);
          end % for all trials
          source.(fn{i}) = val;

          %         else
          %           siz = size(dat);
          %           if prod(siz)==npos
          %             siz = [npos nrpt];
          %           elseif siz(1)==npos
          %             siz = [npos nrpt siz(2:end)];
          %           end
          %           val = nan(siz);
          %           % concatenate all data as pos_rpt_etc
          %           val(:,1,:,:,:) = dat(:);
          %           for j=2:length(source.trial)
          %             dat = source.trial(j).(fn{i});
          %             val(:,j,:,:,:) = dat(:);
          %           end % for all trials
          %           source.(fn{i}) = val;

        end
      end % for each field

      source = rmfield(source, 'trial');

    end % if trial

    % ensure that it has a dimord (or multiple for the different fields)
    source = fixdimord(source);

    if isfield(source, 'inside')
      % ensure that for positions outside the brain it is [], not nan
      if isfield(source, 'leadfield')
        source.leadfield(~source.inside) = {[]};
      end
      if isfield(source, 'filter')
        source.filter(~source.inside) = {[]};
      end
    end

    if isfield(source, 'leadfield') && ~isfield(source, 'label') && isfield(source, 'cfg')
      % try to determine the channel labels from the cfg
      label = ft_findcfg(source.cfg, 'channel');
      if ~isempty(label)
        source.label = label;
      end
    end

    if isfield(source, 'filter') && ~isfield(source, 'label') && isfield(source, 'cfg')
      % try to determine the channel labels from the cfg
      label = ft_findcfg(source.cfg, 'channel');
      if ~isempty(label)
        source.label = label;
      end
    end

    % ensure that all data fields have the correct dimensions
    fn = getdatfield(source);
    for i=1:numel(fn)
      dimord = getdimord(source, fn{i});
      dimtok = tokenize(dimord, '_');
      dimsiz = getdimsiz(source, fn{i}, numel(dimtok));
      if numel(dimsiz)>=3 && strcmp(dimtok{1}, 'dim1') && strcmp(dimtok{2}, 'dim2') && strcmp(dimtok{3}, 'dim3')
        % convert it from voxel-based representation to position-based representation
        try
          source.(fn{i}) = reshape(source.(fn{i}), [prod(dimsiz(1:3)) dimsiz(4:end) 1]);
        catch
          ft_warning('could not reshape %s to the expected dimensions', fn{i});
        end
      end
    end

    % ensure that the structure has all required fields
    for required={'pos'}
      assert(isfield(source, required), 'required field "%s" is missing', required{:});
    end

  case '2011'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % ensure that it has individual source positions
    source = fixpos(source);

    % remove obsolete fields
    if isfield(source, 'xgrid')
      source = rmfield(source, 'xgrid');
    end
    if isfield(source, 'ygrid')
      source = rmfield(source, 'ygrid');
    end
    if isfield(source, 'zgrid')
      source = rmfield(source, 'zgrid');
    end
    if isfield(source, 'transform')
      source = rmfield(source, 'transform');
    end

    % ensure that it has a dimord (or multiple for the different fields)
    source = fixdimord(source);

  case '2010'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % ensure that it has individual source positions
    source = fixpos(source);

    % remove obsolete fields
    if isfield(source, 'xgrid')
      source = rmfield(source, 'xgrid');
    end
    if isfield(source, 'ygrid')
      source = rmfield(source, 'ygrid');
    end
    if isfield(source, 'zgrid')
      source = rmfield(source, 'zgrid');
    end

    % ensure that it has a dimord (or multiple for the different fields)
    source = fixdimord(source);

  case '2007'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % ensure that it has individual source positions
    source = fixpos(source);

    % remove obsolete fields
    if isfield(source, 'dimord')
      source = rmfield(source, 'dimord');
    end
    if isfield(source, 'xgrid')
      source = rmfield(source, 'xgrid');
    end
    if isfield(source, 'ygrid')
      source = rmfield(source, 'ygrid');
    end
    if isfield(source, 'zgrid')
      source = rmfield(source, 'zgrid');
    end

  case '2003'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if isfield(source, 'dimord')
      source = rmfield(source, 'dimord');
    end

    if ~isfield(source, 'xgrid') || ~isfield(source, 'ygrid') || ~isfield(source, 'zgrid')
      if isfield(source, 'dim')
        minx = min(source.pos(:,1));
        maxx = max(source.pos(:,1));
        miny = min(source.pos(:,2));
        maxy = max(source.pos(:,2));
        minz = min(source.pos(:,3));
        maxz = max(source.pos(:,3));
        source.xgrid = linspace(minx, maxx, source.dim(1));
        source.ygrid = linspace(miny, maxy, source.dim(2));
        source.zgrid = linspace(minz, maxz, source.dim(3));
      end
    end

  otherwise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ft_error('unsupported version "%s" for source datatype', version);
end

function pos = grid2pos(xgrid, ygrid, zgrid)
[X, Y, Z] = ndgrid(xgrid, ygrid, zgrid);
pos = [X(:) Y(:) Z(:)];

function pos = dim2pos(dim, transform)
[X, Y, Z] = ndgrid(1:dim(1), 1:dim(2), 1:dim(3));
pos = [X(:) Y(:) Z(:)];
pos = ft_warp_apply(transform, pos, 'homogenous');

author: fieldtrip

repo: fieldtrip


text:
% sample = samplep(p,n)
% returns n samples, where the probability of i is p(i) for 1 <=i
% <= length(p), 0 otherwise.

function s = samplep(p,n)

% normalize, just in case
p = p(:);
p = normalizep(p);
p(end) = 1;
m = length(p);

% compute the cumulative sum of probabilities
c = cumsum(p);

% choose n random numbers uniformly between 0 and 1
r = rand(1,n);

% find the first element of c >= r
s = findfirstdim(repmat(c,[1,n])>=repmat(r,[m,1]),1);

author: kristinbranson

repo: JAABA


text:
function [ p, l, u ] = plu_plu ( n, pivot )

%*****************************************************************************80
%
%% PLU_PLU returns the PLU factors of the PLU matrix.
%
%  Example:
%
%    Input:
%
%      N = 5
%      PIVOT = ( 1, 3, 3, 5, 5 )
%
%    Output:
%
%      P:
%
%          1             0            0             0            0
%          0             0            1             0            0
%          0             1            0             0            0
%          0             0            0             0            1
%          0             0            0             1            0
%
%      L:
%
%         1              0            0             0            0
%         0.25           1            0             0            0
%         0.125          0.375        1             0            0
%         0.0625         0.1875       0.3125        1            0
%         0.03125        0.09375      0.15625       0.21875      1
%
%      U:
%
%        11             12           13            14           15
%         0             22           23            24           25
%         0              0           33            34           35
%         0              0            0            44           45
%         0              0            0             0           55
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer PIVOT(N), the list of pivot rows.  PIVOT(I)
%    must be a value between I and N, reflecting the choice of
%    pivot row on the I-th step.  For no pivoting, set PIVOT(I) = I.
%
%    Output, real P(N,N), L(N,N), U(N,N), the P, L and U factors
%    of A, as defined by Gaussian elimination with partial pivoting.
%    P is a permutation matrix, L is unit lower triangular, and U
%    is upper trianguler.
%

%
%  Check that the pivot vector is legal.
%
  for i = 1 : n

    if ( pivot(i) < i )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PLU_PLU - Fatal error!\n' );
      fprintf ( 1, '  PIVOT(%d) = %d\n', i, pivot(i) );
      fprintf ( 1, '  but PIVOT(I) must be no less than I!\n' );
      error ( 'PLU_PLU - Fatal error!' )
    elseif ( n < pivot(i) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PLU_PLU - Fatal error!\n' );
      fprintf ( 1, '  PIVOT(%d) = %d\n', i, pivot(i) );
      fprintf ( 1, '  but PIVOT(I) must be no greater than N\n' );
      fprintf ( 1, '  and N = %d\n', n );
      error ( 'PLU_PLU - Fatal error!' )
    end

  end
%
%  Compute U.
%
  u = zeros ( n, n );
  for i = 1 : n
    for j = 1 : n

      if ( i <= j )
        u(i,j) = 10 * i + j;
      else
        u(i,j) = 0.0;
      end

    end
  end
%
%  Compute L.
%
  l = zeros ( n, n );
  for i = 1 : n
    for j = 1 : n

      if ( i < j )
        l(i,j) = 0.0;
      elseif ( j == i )
        l(i,j) = 1.0;
      else
        l(i,j) = ( 2 * j - 1 ) / 2^i;
      end

    end
  end
%
%  Compute P.
%
  p = zeros ( n, n );
  for i = 1 : n
    for j = 1 : n
      if ( i == j )
        p(i,j) = 1.0;
      else
        p(i,j) = 0.0;
      end
    end
  end
%
%  Apply the pivot permutations, in reverse order.
%
  for i = n : -1 : 1

    if ( pivot(i) ~= i )
      for j = 1 : n
        temp          = p(i,j);
        p(i,j)        = p(pivot(i),j);
        p(pivot(i),j) = temp;
      end
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


function v = getNonUniformGrid(m_0, lx, ux, gridMethod, center, manualPoint, gridMultParam)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

nx = m_0; dx = (ux-lx)/nx;

if gridMethod == 1  % uniform grid , doesn't include lower bound
    v = lx + (1:nx)'*dx;
elseif gridMethod == 2  %%Mijatovic and Pistorious
    v = zeros(m_0,1);
    v(1) = lx; v(m_0) = ux;
    mid = floor(m_0/2);
    for k = 2: mid
        v(k) = center + sinh( (1-(k-1)/(mid -1))*asinh(v(1)-center) );
    end
    for k = mid+1:m_0-1
        v(k) = center + sinh( ((k - mid)/(mid))*asinh(v(m_0) - center) );
    end
elseif gridMethod == 3 %%Carr and Mayo
    x=0:1/(m_0-1):1; %map [0,1] to new grid
    alpha = 0.8*(ux-lx);
    c1 = asinh((lx-center)/alpha);
    c2 = asinh((ux-center)/alpha);
    v  = center + alpha*(c2*sinh(c2.*x+c1.*(1-x)));
elseif gridMethod == 4 %% Tavella and Randall
    v = zeros(m_0,1);
    v(1)   = lx;
    v(m_0) = ux;
    alpha = gridMultParam*(v(m_0) - v(1));
    c1 = asinh((v(1) - center)/alpha);
    c2 = asinh((v(m_0) - center)/alpha);
    v(2:m_0-1) = center + alpha*sinh(c2/m_0*(2:m_0-1) + c1*(1 - (2:m_0-1)/m_0));
elseif gridMethod == 5 %% Tavella and Randall PLUS we put manualPoint (e.g. v0) on grid
    tol = 1e-7;
    v = zeros(m_0,1);
    v(1)   = lx;
    v(m_0) = ux;
    alpha = gridMultParam*(v(m_0) - v(1));
    c1 = asinh((v(1) - center)/alpha);
    c2 = asinh((v(m_0) - center)/alpha);
    vtil = zeros(m_0-1,1);
    vtil(2:m_0-2) = center + alpha*sinh(c2/(m_0 - 1)*(2:m_0-2) + c1*(1 - (2:m_0-2)/(m_0-1)));
    nnot_til = 1;
    while vtil(nnot_til) < manualPoint %Find the index left of the manualPoint point
        nnot_til = nnot_til + 1;
    end
    nnot_til = nnot_til - 1;
    v(2:nnot_til) = vtil(2:nnot_til);
    v(nnot_til + 2: m_0 - 1) = vtil(nnot_til + 1: m_0 - 2);
    if manualPoint - vtil(nnot_til)<tol %if you already have a point very close to manualPoint
        v(nnot_til) = manualPoint;
        v(nnot_til+1) = (manualPoint + vtil(nnot_til+1))/2;
    elseif vtil(nnot_til+1) - manualPoint <tol
        v(nnot_til+2) = manualPoint;
        v(nnot_til + 1) = (v(nnot_til+2) + v(nnot_til))/2;
    else
        v(nnot_til + 1) = manualPoint;
    end
end

end


text:
function entry = findindict(c,dcode)
% Look up an entry in the dictionary
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

entry = [];
d = getdict;
d = d.(dcode);
if ischar(c)
    for i=1:length(d),
        if strcmpi(d(i).label,c),
            entry = d(i);
            break;
        end;
    end;
elseif isnumeric(c) && numel(c)==1
    for i=1:length(d),
        if d(i).code==c,
            entry = d(i);
            break;
        end;
    end;
else
    error(['Inappropriate code for ' dcode '.']);
end;
if isempty(entry)
    fprintf('\nThis is not an option.  Try one of these:\n');
    for i=1:length(d)
        fprintf('%5d) %s\n', d(i).code, d(i).label);
    end;
    %fprintf('\nNO CHANGES MADE\n');
end;

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function y = pow_pos( x, p )

%POW_POS   Power of positive part.
%   POW_POS(X,P) = POS(X).^P = MAX(X,0).^P.
%   Both P and X must be real, and P must be greater than or equal to 1.
%
%   Disciplined convex programming information:
%       POW_POS(X,P) is convex and nondecreasing in X; so when used in CVX
%       expressions, X must be convex. P must be constant, real, and
%       greater than or equal to 1.

narginchk(2,2);
if ~isnumeric( x ) || ~isreal( x ) || ~isnumeric( p ) || ~isreal( p ),
    error( 'Arguments must be real.' );
elseif any( p(:) <= 1 ),
    error( 'Second argument must be greater than or equal to 1.\nFor other exponents, use POW_P instead.', 1 ); %#ok
end
y = max(x,0).^p;

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function idxSelection = uniformSampling(X, voxelSize)

% Preparations -----------------------------------------------------------------

% No. of points
noPoi = size(X,1);

% Logical indices for selection of points (true = selected)
idxSelection = false(noPoi,1);

% Find voxel centers -----------------------------------------------------------

% Point with smallest coordinates
minPoi = min(X, [], 1);

% Rounded local origin for voxel structure (voxels of different pcs have coincident voxel centers if mod(100, voxelSize) == 0)
localOrigin = (floor(minPoi/100))*100;

% Find 3-dimensional indices of voxels in which points are lying
idxVoxel = [floor( (X(:,1)-localOrigin(1)) / voxelSize ) ...
            floor( (X(:,2)-localOrigin(2)) / voxelSize ) ...
            floor( (X(:,3)-localOrigin(3)) / voxelSize )];

% Remove multiple voxels
[idxVoxelUnique, ~, ic] = unique(idxVoxel, 'rows'); % ic contains "voxel index" for each point

% Coordinates of voxel centers
XVoxelCenter = [localOrigin(1) + voxelSize/2 + idxVoxelUnique(:,1)*voxelSize ...
                localOrigin(2) + voxelSize/2 + idxVoxelUnique(:,2)*voxelSize ...
                localOrigin(3) + voxelSize/2 + idxVoxelUnique(:,3)*voxelSize];

% No. of voxel (equal to no. of selected points)
noVoxel = size(XVoxelCenter,1);

% Select points nearest to voxel centers ---------------------------------------

% Sort indices and points (in order to find points inside of voxels very fast in the next loop)
[ic, idxSort] = sort(ic);
X = X(idxSort,:);

idxJump = find(diff(ic));

% Example (3 voxel)
% ic         = [1 1 1 2 2 2 3]';
% diff(ic)   = [ 0 0 1 0 0 1 ]';
% idxJump    = [     3     6 ]';
%
% idxInVoxel = [1 2 3]; for voxel 1
% idxInVoxel = [4 5 6]; for voxel 2
% idxInVoxel = [7    ]; for voxel 3

for i = 1:noVoxel

    % Find indices of points inside of voxel (very, very fast this way)
    if i == 1
        idxInVoxel = [1:idxJump(i)];
    elseif i == noVoxel
        idxInVoxel = [idxJump(i-1)+1:noPoi];
    else
        idxInVoxel = [idxJump(i-1)+1:idxJump(i)];
    end

    % Distance of points to voxel center
    dist2voxelCenter = sqrt((X(idxInVoxel,1)-XVoxelCenter(i,1)).^2 + (X(idxInVoxel,2)-XVoxelCenter(i,2)).^2 + (X(idxInVoxel,3)-XVoxelCenter(i,3)).^2);

    % Find index of point with smallest distance to voxel center
    [~, idxSelectedPoi] = min(dist2voxelCenter);

    % Select point
    idxSelection(idxSort(idxInVoxel(idxSelectedPoi))) = true;

    % 4Debug
    % idxInVoxel
    % idxSelectedPoi
    % XVoxelCenter(i,:)
    % X(idxInVoxel,:)
    % dist2voxelCenter
    % maxDist = sqrt(3*(voxelSize/2)^2);
    % correct(i) = all(dist2voxelCenter <= maxDist);

end

end

author: pglira


text:
function v = m2v(h,k,l,cs)
% Miller-indece --> cartesian coordinates
% Input
%  h,k,l -
%  cs - crystal symmetry (optional)
%
% Output
%  v - @vector3d

if any(h == 0 & k == 0 & l ==0)
  error('(0,0,0) is not a valid Miller index');
end

a = cs.axes;
V  = dot(a(1),cross(a(2),a(3)));
a_star = cross(a(2),a(3)) ./ V;
b_star = cross(a(3),a(1)) ./ V;
c_star = cross(a(1),a(2)) ./ V;

% reciprocal space
v = h * a_star + k * b_star + l * c_star;
% v = v ./ norm(v);

author: mtex-toolbox

repo: mtex

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


text:
classdef FreeFemInputPrinter < FreeFemFilePrinter

    properties (Access = private)
        inputLines
        beforeInputLines
        afterInputLines

        startingLines
        endingLines
        inputData
    end

    methods (Access = public)

        function setInputData(obj,d)
            obj.inputData = d;
        end
    end

    methods (Access = protected)

        function printingLines(obj)
            obj.obtainBeforeAfterAndInputLines();
            obj.printBeforeInputLines();
            obj.printInputLines();
            obj.printAfterInputLines();
        end

    end

    methods (Access = private)

        function obtainBeforeAfterAndInputLines(obj)
            nlines = numel(obj.linesToPrint);
            [sL,eL] = obj.obtainStartingEndingInputLines();
            obj.beforeInputLines = (1:sL-1);
            obj.inputLines = (sL:eL);
            obj.afterInputLines = (eL+1:nlines);
        end

        function [sL,eL] = obtainStartingEndingInputLines(obj)
            nlines = numel(obj.linesToPrint);
            obj.startingLines = false(nlines,1);
            obj.endingLines   = false(nlines,1);
            for iline = 1:nlines
                line = obj.linesToPrint{iline};
                obj.startingLines(iline) = obj.isInputStarting(line);
                obj.endingLines(iline)   = obj.isInputEnding(line);
            end
            sL = obj.obtainStartingLine();
            eL = obj.obtainEndingLine();
        end

        function sL = obtainStartingLine(obj)
            sLines = obj.startingLines;
            sLines = obj.booleanIndex2integerIndex(sLines);
            firstSL = sLines(1);
            sL = firstSL + 1;
        end

        function eL = obtainEndingLine(obj)
            eLines = obj.endingLines;
            eLines = obj.booleanIndex2integerIndex(eLines);
            firstEL = eLines(1);
            eL = firstEL(1) - 1;
        end

        function printBeforeInputLines(obj)
            lines = obj.beforeInputLines;
            obj.printNonInputLines(lines)
        end

        function printAfterInputLines(obj)
            lines = obj.afterInputLines;
            obj.printNonInputLines(lines)
        end

        function printNonInputLines(obj,nonInputLines)
            nLines = length(nonInputLines);
            for iline = 1:nLines
                lineNumber = nonInputLines(iline);
                line = obj.linesToPrint{lineNumber};
                obj.printLine(line);
            end
        end

        function printInputLines(obj)
            obj.printRealInputData();
            obj.printStringInputData();
        end

        function printRealInputData(obj)
            inVar = obj.inputData.reals;
            for iData = 1:numel(inVar)
                varName  = inVar{iData}{1,1};
                varValue = num2str(inVar{iData}{1,2});
                line = ['real ',varName,' = ',varValue,';'];
                obj.printLine(line);
            end
        end

        function printStringInputData(obj)
            inVar = obj.inputData.strings;
            for iData = 1:numel(inVar)
                varName  = inVar{iData}{1,1};
                varValue = ['"',inVar{iData}{1,2},'"'];
                line = ['string ',varName,' = ',varValue,';'];
                obj.printLine(line);
            end
        end

        function printInputData(obj)
            inVar = obj.inputData.reals;
            for iData = 1:numel(inVar)
                varName  = inVar{iData}{1,1};
                varValue = num2str(inVar{iData}{1,2});
                line = ['real ',varName,' = ',varValue,';'];
                obj.printLine(line);
            end
        end

        function printInputLine(obj)
            fprintf(obj.fileID,'%s\r\n',obj.line);
        end

        function printLine(obj,line)
            fprintf(obj.fileID,'%s\r\n',line);
        end

    end

    methods (Access = private, Static)

        function itIs = isInputStarting(line)
            itIs = strcmp(line,'// --- Input ---');
        end

        function itIs = isInputEnding(line)
            itIs = strcmp(line,'// --- End Input ---');
        end

        function index = booleanIndex2integerIndex(index)
            allValues = 1:length(index);
            index = allValues(index);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function grd= util_gridForLaplacian
% UTIL_GRIDFORLAPLACIAN - Channel grid to define neighbors for Laplacians
%
%Synopsis:
%  grd= util_gridForLaplacian

% 09-2012 Benjamin Blankertz

str= [
      '  _     _     _     _   Fp1   Fpz   Fp2     _     _     _     _\n' ...
      '  <     _     _     _     _  AFp1  AFp2     _     _     _     _\n' ...
      'AF9   AF7   AF5   AF3   AF1   AFz   AF2   AF4   AF6   AF8  AF10\n' ...
      '  <  FAF9  FAF7  FAF5  FAF3  FAF1  FAF2  FAF4  FAF6  FAF8 FAF10\n' ...
      ' F9    F7    F5    F3    F1    Fz    F2    F4    F6    F8   F10\n' ...
      '  <  FFC9  FFC7  FFC5  FFC3  FFC1  FFC2  FFC4  FFC6  FFC8 FFC10\n' ...
      'FT9   FT7   FC5   FC3   FC1   FCz   FC2   FC4   FC6   FT8  FT10\n' ...
      '  <  CFC9  CFC7  CFC5  CFC3  CFC1  CFC2  CFC4  CFC6  CFC8 CFC10\n' ...
      ' T9    T7    C5    C3    C1    Cz    C2    C4    C6    T8   T10\n' ...
      '  <  CCP9  CCP7  CCP5  CCP3  CCP1  CCP2  CCP4  CCP6  CCP8 CCP10\n' ...
      'TP9   TP7   CP5   CP3   CP1   CPz   CP2   CP4   CP6   TP8  TP10\n' ...
      '  <  PCP9  PCP7  PCP5  PCP3  PCP1  PCP2  PCP4  PCP6  PCP8 PCP10\n' ...
      ' P9    P7    P5    P3    P1    Pz    P2    P4    P6    P8   P10\n' ...
      '  <  PPO9  PPO7  PPO5  PPO3  PPO1  PPO2  PPO4  PPO6  PPO8 PPO10\n' ...
      'PO9   PO7   PO5   PO3   PO1   POz   PO2   PO4   PO6   PO8  PO10\n' ...
      '  <     _     _     _     _  OPO1  OPO2     _     _     _     _\n' ...
      '  _     _     _    O9    O1    Oz    O2   O10     _     _     _\n' ...
      '  <     _     _     _     _   OI1   OI2     _     _     _     _\n' ...
      '  _     _     _     _     I1   Iz    I2     _     _     _     _'];

nCols= 11;
G= cell(1, nCols);
[G{:}]= strread(sprintf(str), repmat('%s', [1 nCols]));
grd= cat(2, G{:});

author: bbci


text:
function vscl = vscale(f)
%VSCALE   Vertical scale of a SEPARABLEAPPROX.
%
% VSCL = VSCALE(F) returns the vertial scale of a SEPARABLEAPPROX as determined
% by evaluating on a coarse tensor-product grid.

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

% TODO: Should this also be taking the maximum along the edges when we are
% evaluating at 1st kind grids.

% If f is an empty SEPARABLEAPPROX, VSCL = 0:
if ( isempty( f ) )
    vscl = 0;
    return
end

% Get the degree of the SEPARABLEAPPROX:
[m, n] = length(f);

% If F is of low degree, then oversample:
m = min(max(m, 9),2000);
n = min(max(n, 9),2000); % cannot afford to go over 2000x2000.

% Calculate values on a tensor grid:
vals = sample(f, m, n);

% Take the absolute maximum:
vscl = max(abs(vals(:)));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [y,ev] = dmrg_eigb(a,k,eps,varargin)
%Find several minimal eigenvalues of a TT-matrix using DMRG method
%   [Y,EV]=DMRG_EIGB(A,K,EPS,OPTIONS) Attempts to find K minimal
%   eigenvalues of a TT-matrix A with accuracy EPS. We use minimization of
%   Block-Rayleigh quotient to do this. The solution is returned a block
%   TT-tensor (i.e, r(d+1) is equal to K).
%   Options are provided in form
%   'PropertyName1',PropertyValue1,'PropertyName2',PropertyValue2 and so
%   on. The parameters are set to default (in brackets in the following)
%   The list of option names and default values are:
%       o y0 - initial approximation [random rank-5 tensor]
%       o rmax - maximal  TT-rank of the (block) solution [2500]
%       o nswp - maximal number of sweeps [4]
%       o kick_rank - stabilization parameter, the larger, the better
%       accuracy but the higher complexity [5]
%       o verb - print debug info [ {true} | false ]
%       o msize - the size of local matrix when the iterative solver is
%       used
%   Example:
%       d=8; f=3;
%       mat=tt_qlaplace_dd(d*ones(1,f)); %Laplace in the QTT-format
%       [v,ev]=dmrg_eigb(mat,5,1e-6); %5 lowest eigenvalues
%
%
% 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
%---------------------------

warning('This realisation is old and not optimal. Consider using dmrg_eig');

%Default parameters
y0=[];
nswp=4;
msize=1;
max_l_steps=200;
kickrank=4;
verb=1;
rmax = Inf;
for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'nswp'
            nswp=varargin{i+1};
        case 'kickrank'
            kickrank=lower(varargin{i+1});
        case 'x0'
            y0=varargin{i+1};
        case 'msize'
            msize=varargin{i+1};
        case 'verb'
            verb=varargin{i+1};
        case 'rmax'
            rmax=varargin{i+1};
        otherwise
            error('Unrecognized option: %s\n',varargin{i});

    end
end

n=a.n;
d=a.d;
if ( isempty(y0) )
    kk=max(5,k);
    r=kk*ones(1,d+1);
    r(d+1)=k;
    r(1)=1;
    y0=tt_random(n,d,r);
end
y=round(y0,0);

local_eps = eps/sqrt(d);

%We start from the orthogonalization of the y vector from left-to-right
%(it does not influence the TT-ranks)

psy=y.ps;
ry=y.r;
tta=a.tt;
psa=tta.ps;
ra=tta.r;
cra=tta.core;
cry=y.core;
rm=1;
%Also we will need to compute phi-matrices
phi=cell(d+1,1);
phi{1}=1; %This one is right now unchanged
phi{d+1}=1; %Seems to be unchanged also
pos1=1;
for i=1:d
    cr=cry(psy(i):psy(i+1)-1);
    cr=reshape(cr,[ry(i),n(i)*ry(i+1)]);
    cr=rm*cr;
    ry(i)=size(cr,1);
    cr=reshape(cr,[ry(i)*n(i),ry(i+1)]);
    [u,rm]=qr(cr,0);
    rnew=size(u,2);
    %ry(i+1)=rnew;
    cry(pos1:pos1+ry(i)*n(i)*rnew-1)=u(:);
    pos1=pos1+ry(i)*n(i)*rnew;
    %With this new core compute (Ax,x) to the phi matrix
    %u=reshape(u,[ry(i),n(i),ry(i+1)]);
    crm=cra(psa(i):psa(i+1)-1);
    crm=reshape(crm,[ra(i),n(i),n(i),ra(i+1)]);
    phx=phi{i};
    phx=reshape(phx,[ra(i),ry(i),ry(i)]); phx=permute(phx,[1,3,2]);
    x0=u;
    x0=reshape(x0,[ry(i),n(i)*rnew]);
    phx=reshape(phx,[ra(i)*ry(i),ry(i)]);
    phx=phx*x0; %phx is ra(i)*ry(i)*n(i)*ry(i+1)
    phx=reshape(phx,[ra(i),ry(i),n(i),rnew]);
    crm=permute(crm,[1,3,2,4]); crm=reshape(crm,[ra(i)*n(i),n(i)*ra(i+1)]);
    %Convolve over (ak-1) j with the matrix
    phx=permute(phx,[2,4,1,3]);
    phx=reshape(phx,[ry(i)*rnew,ra(i)*n(i)]);
    phx=phx*crm; %ry(i)*ry(i+1)*n(i)*ra(i+1)
    phx=reshape(phx,[ry(i),rnew,n(i),ra(i+1)]);
    phx=permute(phx,[2,4,1,3]);
    phx=reshape(phx,[rnew*ra(i+1),ry(i)*n(i)]);
    x0=u;
    x0=reshape(x0,[ry(i)*n(i),rnew]);
    phx=phx*x0;
    phx=reshape(phx,[rnew,ra(i+1),rnew]);
    phx=permute(phx,[2,1,3]);
    phi{i+1}=phx;
    %phi{i+1}=permute(phi{i+1},[1,3,2]);
end

ry(d+1)=rnew; %This value should not be modified :(((
%Truncate the core
cry=cry(1:pos1-1);
y.r=ry;
y.core=cry;
y.ps=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
 %Test back transform
%  y0=y;
%  y0.n=[2*ones(d,1);k];
%  y0.r=[y.r;1];
%  y0.d=d;
%  e0=eye(ry(d));
%  cry=[cry,e0(:)'];
%  y0.core=cry;
%  ww1=full(y0); ww1=reshape(ww1,[numel(ww1)/k,k]);
%  bm=ww0'*ww1;
%  norm(ww1-ww0*bm)
%keyboard;
phi{d+1}=1; %Bydlocode, but seems necessary
y.r=ry;
y.core=cry;
y.ps=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
psy=y.ps;
%Now start the main DMRG sweep
swp=1;
not_converged=true;
err_max = 0;
%ry(1)=1;
ry(d+1)=1;
dir='rl';
i=d-1;
cry_left=cry(1:psy(d)-1);
cry_right=cry(psy(d):psy(d+1)-1);
while ( swp <= nswp && not_converged )
       %As usual, compute (local) B-matrix in the TT-format
       cra1=cra(psa(i):psa(i+1)-1);
       cra2=cra(psa(i+1):psa(i+2)-1);
       px1=phi{i}; px1=reshape(px1,[ra(i),ry(i),ry(i)]);
       %px1=permute(px1,[1,2,3]);
       px1=permute(px1,[1,3,2]);
       px1=reshape(px1,[ra(i),ry(i)*ry(i)]);
       px2=phi{i+2};
       px2=reshape(px2,[ra(i+2),ry(i+2),ry(i+2)]);
       %px2=permute(px2,[1,3,2]);
       %Compute the local matrix just by putting px1 into cra1, px2 into
       %cra2
       cra2=reshape(cra2,[ra(i+1)*n(i+1)*n(i+1),ra(i+2)]);
       px2=reshape(px2,[ra(i+2),ry(i+2)*ry(i+2)]);
       b2=cra2*px2; %
       b1=px1.'*reshape(cra1,[ra(i),numel(cra1)/ra(i)]);
       b1=reshape(b1,[ry(i),ry(i),n(i),n(i),ra(i+1)]); b1_save=b1; %Save for phi computations
       b1=permute(b1,[1,3,2,4,5]);
       b1=reshape(b1,[ry(i)*n(i),ry(i)*n(i),ra(i+1)]);
       b2=reshape(b2,[ra(i+1),n(i+1),n(i+1),ry(i+2),ry(i+2)]);
       b2_save=b2; %Save for phi computations
       b2=permute(b2,[2,4,3,5,1]);
       b2=reshape(b2,[n(i+1)*ry(i+2),n(i+1)*ry(i+2),ra(i+1)]);
       mm=cell(2,1); mm{1}=b1; mm{2}=b2;
       cry=[cry_left(:);cry_right(:)];
       y.r=ry;
       y.core=cry;
       y.ps=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
       %mm1=tt_matrix(mm);
       %Now setup the initial guess: i the core
       if ( strcmp(dir,'rl') ) %The block index is in the second core
         pos=numel(cry_left);
         cry1=cry_left(pos-ry(i)*n(i)*ry(i+1)+1:pos);
         cry2=cry_right(1:ry(i+1)*n(i+1)*k*ry(i+2));
         cry1=reshape(cry1,[ry(i)*n(i),ry(i+1)]);
         cry2=reshape(cry2,[ry(i+1),n(i+1)*k*ry(i+2)]);
         w=cry1*cry2; w=reshape(w,[ry(i),n(i),n(i+1),k,ry(i+2)]);
         w=permute(w,[1,2,3,5,4]); w=reshape(w,[numel(w)/k,k]);
       else %The block index is in the first core
           pos=numel(cry_left);
           cry1=cry_left(pos-ry(i)*n(i)*k*ry(i+1)+1:pos);
           cry2=cry_right(1:ry(i+1)*n(i+1)*ry(i+2));
           cry1=reshape(cry1,[ry(i)*n(i)*k,ry(i+1)]);
           cry2=reshape(cry2,[ry(i+1),n(i+1)*ry(i+2)]);
           w=cry1*cry2; w=reshape(w,[ry(i),n(i),k,n(i+1),ry(i+2)]);
           w=permute(w,[1,2,4,5,3]); w=reshape(w,[numel(w)/k,k]);
       end
       %Now run the eigenvalue solver
        bw=bfun(mm,w); ev=bw'*w;
           %er0=norm(bw-w*ev,'fro')/norm(w,'fro');
           er_old=bw-w*ev;
        erc=zeros(1,size(w,2));
        for j=1:size(w,2)
            erc(j)=norm(er_old(:,j));
        end
           er0=norm(er_old,'fro')/norm(w,'fro');
       if ( size(w,1) >= max(5*k,msize) )
           matvec='bfun';
           [wnew,ev,fail_flag]=lobpcg(w,@(x) bfun(mm,x),local_eps,max_l_steps);
       else
          matvec='full';
          fm=full(tt_matrix(mm));
          [v,dg]=eig(fm);
          %[v,dg]=eigs(sparse(fm),k+1);
          %keyboard;
          ev=diag(dg);
          [ev,ind]=sort(ev,'ascend');
          v=v(:,ind);
          wnew=v(:,1:k);
          ev=ev(1:k);
       end
       er1=norm(bfun(mm,wnew)-wnew*diag(ev),'fro')/norm(wnew,'fro');

       % Check the distance between iterations
       dx = norm(wnew*(wnew'*w)-w);

       fv=sum(ev); %The functional we minimize;
       if ( verb>1 )
         fprintf('sweep=%d block=%d fv=%10.15f loc solve=%3.2e old_solve=%3.2e \n',swp,i,fv,er1,er0);
         disp(erc);
       end
       if ( strcmp(dir,'rl') ) %Implant the auxiliary core into the i-th core
           %(a1,i1,a2,a2,i2*g,a3)-> (a1,i1*g,a2,a2,i2,a3)
           %Delete old block from the core_left, add new block to the core
           %right

           %Prepare the truncation block
           rhs=wnew*diag(ev);
            if (strcmp(matvec,'full'))
                res_true = norm(fm*wnew-rhs)/norm(rhs);
            else
                res_true = norm(bfun(mm,wnew)-rhs)/norm(rhs);
            end;

           wnew=reshape(wnew,[ry(i),n(i),n(i+1),ry(i+2),k]);
           wnew=permute(wnew,[1,2,5,3,4]); wnew=reshape(wnew,[ry(i)*n(i)*k,n(i+1)*ry(i+2)]);
           [u,s,v]=svd(wnew,'econ'); s=diag(s);
           %Truncation block
           rnew=my_chop2(s,local_eps*norm(s));

           %%%
           rnew = min(rnew, numel(s));
           rnew = min(rnew, rmax);
           %%%%

           u=u(:,1:rnew); s=s(1:rnew); v=v(:,1:rnew);% v=v';
           u=u*diag(s); %u has to be reshaped

           % Residual truncation
%            r0=1; r1=min(numel(s),rmax);
%            r=1;
%            while ( (r ~= r0 || r ~= r1) && r0 <= r1)
%             r=min(floor((r0+r1)/2),rmax);
%             %er0=norm(s(r+1:numel(s)));
%             u1=u(:,1:r)*diag(s(1:r));
%             %Sonate u1
%             u1=reshape(u1,[ry(i)*n(i),k,r]);
%             u1=permute(u1,[1,3,2]);
%             u1=reshape(u1,[numel(u1)/k,k]);
%             [u2,~,v2]=svd(u1,'econ');
%             u1=u2*v2';
%             u1=reshape(u1,[ry(i)*n(i),r,k]);
%             u1=permute(u1,[1,3,2]);
%             u1=reshape(u1,[numel(u1)/r,r]);
%             sol = u1*(v(:,1:r))';
%
%             sol = reshape(sol,[ry(i),n(i),k,n(i+1),ry(i+2)]);
%             sol=permute(sol,[1,2,4,5,3]); sol=reshape(sol,[numel(sol)/k,k]);
%             %if ( norm(sol'*sol-eye(k))>1e-3 )
%             %  keyboard
%             %end
%             if (strcmp(matvec,'full'))
%                 resid = norm(fm*sol-rhs)/norm(rhs);
%             else
%                 resid = norm(bfun(mm,sol)-rhs)/norm(rhs);
%             end;
%             if ( verb )
%             fprintf('sweep %d, block %d, r0=%d, r1=%d, r=%d, resid=%g, MatVec=%s\n', swp, i, r0, r1, r, resid,matvec);
%             end
%             if ((resid<max(res_true*1.2, eps)) ) %Value of the rank is OK
%               r1=r;
%             else %Is not OK.
%               r0=min(r+1,rmax);
%             end;
%
%            end
%            rnew=r;
%            u=u1; v=v(:,1:rnew);

%            if ( norm(sol'*sol-eye(k)) > 1e-7 )
%              keyboard;
%            end
           %u=u(:,1:rnew); s=s(1:rnew); v=v(:,1:rnew);% v=v';
           %u=u*diag(s); %u has to be reshaped

           %Random restart block
           radd=min(kickrank,size(v,1)-rnew);
           rnew=rnew+radd;
           if ( radd >  0 )
             vr=randn(size(v,1),radd);
             ur=zeros(size(u,1),radd);
             %Orthogonalize vr to v by Golub-Kahan reorth
             mvr=v'*vr; vnew=vr-v*mvr;
             reort_flag=false;
             for j=1:radd
                if ( norm(vnew(:,j)) <= 0.5*norm(vr(:,j)))
                   reort_flag=true;
                end
             end
             if (reort_flag)
                 sv=v'*vnew;
                 %mvr=mvr+v'*vnew;
                 vnew=vnew-v*sv;
             end
             [vnew,~]=qr(vnew,0);

             v=[v,vnew];
             u=[u,ur];
             %norm(u*v'-u1*v1');
             %keyboard;
             %keyboard;
           end
           v=v';

           %Memory stuff
           pos=numel(cry_left);
           cry_left(pos-ry(i)*n(i)*ry(i+1)+1:pos)=[];
           cry_right(1:ry(i+1)*n(i+1)*k*ry(i+2))=[]; %Delete the top core from cry_right
           cry_right=[u(:)', v(:)',cry_right]; %Simply add new block to cry_right
           ry(i+1)=rnew;

           %Recalculate phi block; we need to recalculate phi{i+1} here
           %using phi{i+2}
           %px2(ra(i+2),ry(i+2),ry(i+2))*cra2(ra(i+1),n(i),m(i),ra(i+2)*v(ry(i+1),n(i+1),ry(i+2))*v(ry(i+1),n(i+1),ry(i+2))
           %we already have b2 ---  (b2_save)
           %(ra(i+1),n(i+1),n(i+1),ry(i+2),ry(i+2)), convolve over the
           %n(i+1),ry(i+2)
           phx=reshape(b2_save,[ra(i+1),n(i+1),n(i+1),ry(i+2),ry(i+2)]);
           phx=permute(phx,[2,4,1,3,5]);
           %phx=permute(phx,[3,4,1,2,5]);
           %phx=permute(phx,[3,5,1,2,4]);
           % phx=permute(phx,[2,5,1,3,4]);
           phx=reshape(phx,[n(i+1)*ry(i+2),ra(i+1)*n(i+1)*ry(i+2)]);
           v0=v;
           v0=reshape(v0,[ry(i+1),n(i+1)*ry(i+2)]);
           phx=v0*phx;
           phx=reshape(phx,[ry(i+1),ra(i+1),n(i+1),ry(i+2)]);
           phx=permute(phx,[3,4,1,2]);
           phx=reshape(phx,[n(i+1)*ry(i+2),ry(i+1)*ra(i+1)]);
           phx=v0*phx;
           phx=reshape(phx,[ry(i+1),ry(i+1),ra(i+1)]);
           %phx=permute(phx,[3,2,1]);
           phx=permute(phx,[3,2,1]);
           phi{i+1}=phx;
       else %Implant the auxiliary core from the left block to the right block

           %Prepare the truncation block
           rhs=wnew*diag(ev);
            if (strcmp(matvec,'full'))
                res_true = norm(fm*wnew-rhs)/norm(rhs);
            else
                res_true = norm(bfun(mm,wnew)-rhs)/norm(rhs);
            end;
           wnew=reshape(wnew,[ry(i),n(i),n(i+1),ry(i+2),k]);
           wnew=permute(wnew,[1,2,3,5,4]); wnew=reshape(wnew,[ry(i)*n(i),n(i+1)*k*ry(i+2)]);
           [u,s,v]=svd(wnew,'econ'); s=diag(s);
           rnew=my_chop2(s,local_eps*norm(s));

           %%%%
           rnew = min(rnew, rmax);
           rnew = min(rnew, numel(s));
           %%%%

           u=u(:,1:rnew); s=s(1:rnew); v=v(:,1:rnew);% v=v';
           v=v*diag(s); %u has to be reshaped

           % Residual truncation
%            r0=1; r1=min(size(s,1),rmax);
%            r=1;
%            while ( (r ~= r0 || r ~= r1) && r0 <= r1)
%             r=min(floor((r0+r1)/2),rmax);
%             %er0=norm(s(r+1:numel(s)));
%             v1=v(:,1:r)*diag(s(1:r));
%             %Sonate v1
%             v1=reshape(v1,[n(i+1),k,ry(i+2),r]);
%             v1=permute(v1,[1,3,4,2]);
%             v1=reshape(v1,[numel(v1)/k,k]);
%             [u2,~,v2]=svd(v1,'econ');
%             v1=u2*v2';
%             v1=reshape(v1,[n(i+1),ry(i+2),r,k]);
%             v1=permute(v1,[1,4,2,3]);
%             v1=reshape(v1,[numel(v1)/r,r]);
%             sol=u(:,1:r)*v1';
%             sol = reshape(sol,[ry(i),n(i),n(i+1),k,ry(i+2)]);
%             sol=permute(sol,[1,2,3,5,4]); sol=reshape(sol,[numel(sol)/k,k]);
%             if (strcmp(matvec,'full'))
%                 resid = norm(fm*sol-rhs)/norm(rhs);
%             else
%                 resid = norm(bfun(mm,sol)-rhs)/norm(rhs);
%             end;
%             if ( verb )
%             fprintf('sweep %d, block %d, r0=%d r1=%d r=%d, resid=%g, MatVec=%s\n', swp, i, r0, r1, r, resid,matvec);
%             end
%             if ((resid<max(res_true*1.2, eps)) ) %Value of the rank is OK
%               r1=r;
%             else %Is not OK.
%               r0=min(r+1,rmax);
%             end;
%            end
%            rnew=r;
%            %Truncation block
%            %rnew=my_chop2(s,eps*norm(s));
%            u=u(:,1:rnew); v=v1;

           %Random restart block
           radd=min(kickrank,size(u,1)-rnew);
           rnew=rnew+radd;
           if ( radd >  0 )
             vr=zeros(size(v,1),radd);
             ur=randn(size(u,1),radd);
             %Orthogonalize vr to v by Golub-Kahan reorth
             mur=u'*ur; unew=ur-u*mur;
             reort_flag=false;
             for j=1:radd
                if ( norm(unew(:,j)) <= 0.5*norm(ur(:,j)))
                   reort_flag=true;
                end
             end
             if (reort_flag)
                 sv=u'*unew;
                 %mvr=mvr+v'*vnew;
                 unew=unew-u*sv;
             end
             [unew,~]=qr(unew,0);

             u=[u,unew];
             v=[v,vr];
           end
           v=v';

           cry_right(1:ry(i+1)*n(i+1)*ry(i+2))=[];
           pos=numel(cry_left);
           cry_left(pos-ry(i)*n(i)*k*ry(i+1)+1:pos)=[];
           cry_left=[cry_left,u(:)',v(:)'];

           ry(i+1)=rnew;

           %Recalculate phi block; we need to recalculate phi{i+1} using
           %phi{i}
           phx=reshape(b1_save,[ry(i),ry(i),n(i),n(i),ra(i+1)]);
           u0=u; u0=reshape(u0,[ry(i)*n(i),ry(i+1)]);
           phx=permute(phx,[1,3,5,2,4]);
           %phx=permute(phx,[2,4,5,1,3]);
           %phx=permute(phx,[1,4,5,2,3]);
           %phx=permute(phx,[2,3,5,1,4]);
           phx=reshape(phx,[ry(i)*n(i)*ra(i+1),ry(i)*n(i)]);
           phx=phx*u0;
           phx=reshape(phx,[ry(i),n(i),ra(i+1),ry(i+1)]);
           phx=permute(phx,[3,4,1,2]);
           phx=reshape(phx,[ra(i+1)*ry(i+1),ry(i)*n(i)]);
           phx=phx*u0; phx=reshape(phx,[ra(i+1),ry(i+1),ry(i+1)]);
           phx=permute(phx,[1,2,3]);
           phi{i+1}=phx;
       end

   err_max = max(err_max, dx);
   %Choose the next direction block; now implement the simple case
   if ( strcmp(dir,'rl') )
     if ( i > 1 )
       i=i-1;
     else %Change direction %The last optimization was for (1,2) core
       dir='lr';
       %One block should go from cry_right to cry_left
       cry_left=cry_right(1:ry(1)*n(1)*ry(2)*k); %This seems correct
       cry_right(1:ry(1)*n(1)*ry(2)*k)=[];
       if (verb>0)
        fprintf('--- sweep=%d, err_max=%3.3e \n',swp, err_max);
       end;
       swp=swp+1;
       not_converged = (err_max>eps);
       err_max = 0;
     end
   else
    if ( i < d-1 )
        i=i+1;
     else
       dir='rl';
       pos=numel(cry_left);
       cry_right=cry_left(pos-ry(d)*n(d)*ry(d+1)*k+1:pos);
       cry_left(pos-ry(d)*n(d)*ry(d+1)*k+1:pos)=[];
       if (verb>0)
        fprintf('--- sweep=%d, err_max=%3.3e \n',swp, err_max);
       end;
       swp=swp+1;
       not_converged = (err_max>eps);
       err_max = 0;
       %One block should go from cry_left to cry_right (?) --- seems no :)
     end
   end

end
  %Gather the solution

   ry(d+1)=k; %This is the final
   y.r=ry;
   cry=[cry_left,cry_right];
   y.core=cry;
   y.r=ry;
   y.d=d;
   y.ps=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
end

    function [x]=bfun(a,x)
    %[Y]=BFUN(A,X,K)
    %a is given as U(i1,j1,s)*V(i2,j2,s)
    %\sum_{j1,j2,s} U(i1,j1,s)*V(i2,j2,s)*X(j1,j2,q)
    re=size(x,2);
    ul=a{1}; vl=a{2};
    n1=size(ul,1);m1=size(ul,2); ral=size(ul,3);
    n2=size(vl,1); m2=size(vl,2);
    ul=permute(ul,[1,3,2]); ul=reshape(ul,[numel(ul)/m1,m1]);
    x=reshape(x,[m1,numel(x)/m1]);
    x=ul*x; %x is i1xsxj2xq s,j2
    x=reshape(x,[n1,ral,m2,re]);
    x=permute(x,[3,2,1,4]);
    x=reshape(x,[m2*ral,n1*re]);
    vl=reshape(vl,[n2,m2*ral]);
    x=vl*x; %is n2*n1*k
    x=reshape(x,[n2,n1,re]);
    x=permute(x,[2,1,3]);
    x=reshape(x,[numel(x)/re,re]);
    return
    end

author: oseledets

repo: TT-Toolbox


text:
% This script stores in a file the standard configuration of RP for using
% just one segmentation in LAB colorspace. The parameters were trained with
% the VOC 2007 dataset. Further details can be found in the paper:
%
%"S. Manen, M. Guillaumin, and L. Van Gool. Prime Object Proposals with
%Randomized Prim's Algorithm. In ICCV, 2013."
function generateRPConfig(rppath)
	%% Parameter specification:
	params.approxFinalNBoxes = 10000;                                    %Approximate number of proposals
	params.rSeedForRun = -1;                                            %Random seed to be used (-1 to generate it with a hashing function)
	params.q = 10;                                                      %Parameter to eliminate near duplicates (raise it to eliminate more duplicates)

	% LAB segmentation
	params.segmentations{1}.colorspace = 'LAB';                         %Colorspace: 'RGB', 'LAB', 'opponent', 'rg', 'HSV'

	% --> Segmentation parameters:
	params.segmentations{1}.superpixels.sigma = 0.8;
	params.segmentations{1}.superpixels.c = 100;
	params.segmentations{1}.superpixels.min_size = 100;
	% --> Parameters trained from VOC07:

	%   --> Feature weights:
	params.segmentations{1}.simWeights.wBias = 3.0017;
	params.segmentations{1}.simWeights.wCommonBorder = -1.0029;
	params.segmentations{1}.simWeights.wLABColorHist = -2.6864;
	params.segmentations{1}.simWeights.wSizePer = -2.3655;

	%   --> Size term alpha, as explained in the paper sec. 4.2.
	%       It is quantized to contain exactly 2^16 elements for speed
	%       purposes.
	params.segmentations{1}.alpha = dlmread(fullfile(rppath, 'config', 'alpha/alpha_voc07.dat'));
	params.segmentations{1}.verbose = false;                            %Set to true to display more information during execution

	%% Save parameters:
	save(fullfile(rppath, 'config', 'rp.mat'),'params');
end

author: Cloud-CV

repo: object-proposals


text:
function [past_track_plot, future_track_plot, other_eddy_plots] = ...
    plot_tracks( hdls, past_tracks, future_tracks, past_track_color, future_track_color )
%PLOT_TRACKS Plot past tracks, future tracks and eddies in the past and future (not the current ones)

if isempty(past_tracks)
    past_tracks = [NaN NaN NaN]; % Plotm requires the input to be not empty
end
past_track_plot =  plotm(past_tracks(:, 1:2), 'Color', past_track_color, 'Parent', hdls.mapAx);

if isempty(future_tracks)
    future_tracks = [NaN NaN NaN]; % Plotm requires the input to be not empty
end
future_track_plot =  plotm(future_tracks(:, 1:2), 'Color', future_track_color, 'Parent', hdls.mapAx);

% Plot other eddies by types
other_eddy_plots = plot_eddy_by_type(hdls, [past_tracks;future_tracks], hdls.constants.eddy_types, ...
    hdls.constants.eddy_markers, hdls.constants.eddy_colors, hdls.constants.OTHER_EDDY_SIZE, false);

end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
function startmen(parent_fig)
    % startmen adds menus for the basic zmap functions (Data, help)

    report_this_filefun();
    %
    %  This file display the original menu
    %
    %  Stefan Wiemer 12/94

    %  Create new figure
    % Find out if figure already exists
    %

    % Set up the Seismicity Map window Enviroment
    %
    if nargin==1
        disp(parent_fig)
        genmen = uimenu(parent_fig,'Label','Data');

        uimenu(genmen, ...
            'Label','Load Catalog (*.mat file)',...
            'MenuSelectedFcn', @(~,~)cb_load_file);%ZmapImportManager(@load_zmapfile));
        uimenu(genmen, ...
            'Label','Import Catalog from other formatted file',... %was Data ImportFilters
            'MenuSelectedFcn', @(~,~)ZmapImportManager(@zdataimport));
        uimenu(genmen, ...
            'Label','FDSN web fetch',... %TODO
            'MenuSelectedFcn', @(~,~)cb_load_web);%ZmapImportManager(@get_fdsn_data_from_web_callback));
        uimenu(genmen, ...
            'Label', 'Create or Modify *.mat datafile',...
            'MenuSelectedFcn', @(~,~)think_and_do('setup'));
        uimenu(genmen, ...
            'Label', 'Current Dataset Info',...
            'Enable','off',... % may no longer be relevent
            'MenuSelectedFcn', @(~,~)datinf());

        genmen = uimenu(parent_fig,'Label','Help');

        uimenu(genmen, ...
            'Label','Introduction and Help',...
            'MenuSelectedFcn', @(s,e) showweb('new'));

        uimenu(genmen, ...
            'Label','Sample Slide Show',...
            'MenuSelectedFcn',@(s,e) slshow());

    end
    function cb_load_file
        ok=ZmapImportManager(@load_zmapfile);
        if ok
            ZG= ZmapGlobal.Data;
            ZmapMainWindow([],ZG.primeCatalog);
        end
    end
    function cb_load_web
        ok=ZmapImportManager(@get_fdsn_data_from_web_callback);
        if ok
            ZG= ZmapGlobal.Data;
            ZmapMainWindow([],ZG.primeCatalog);
        end
    end
end

function think_and_do(~, ~, f_handle, varargin)

    if ischar(f_handle) && nargin==3
        evalin('base',f_handle); %use evalin??
    else
        f_handle(varargin{:});
    end
end

author: CelsoReyes


text:
function testLoading(rootdir)

% rootdir = '/misc/public/Adamtokristin'

matfiles = dir(fullfile(rootdir,'*.mat'));
for i = 1:numel(matfiles),
  matfile = fullfile(rootdir,matfiles(i).name);
  tmp = load(matfile);
  fprintf('Loaded from %s:\n',matfile);
  disp(tmp);
end

subdirs = mydir(rootdir,'isdir',true);
for i = 1:numel(subdirs),
  testLoading(subdirs{i});
end


text:
function [gKern, gVarmeans, gVarcovars, gInd] = kernVardistPsi2Gradient(kern, vardist, Z, covGrad, learnInducing)

% KERNVARDISTPSI2GRADIENT description.

% VARGPLVM

if nargin < 5
    % If learnInducing == 0, then gInd will be returned as a a vector of
    % zeros. It would be better if nargout was checked and gInd was not
    % returned at all if not needed, but the code is very hierarchical and
    % recursive and that would require more complicated code than just
    % including the "learnInducing" flag.
    learnInducing = 1;
end

    %% compute first the "square" terms of Psi2
    %
    if ~strcmp(kern.type,'cmpnd')
       %
       fhandle = str2func([kern.type 'VardistPsi2Gradient']);
       [gKern, gVarmeans, gVarcovars, gInd] = fhandle(kern, vardist, Z, covGrad, learnInducing);

       % Transformations
       gKern = paramTransformPsi2(kern, gKern);
       %
    else % the kernel is cmpnd
       %
       fhandle = str2func([kern.comp{1}.type 'VardistPsi2Gradient']);
       [gKern, gVarmeans, gVarcovars, gInd] = fhandle(kern.comp{1}, vardist, Z, covGrad, learnInducing);
       % Transformations
       gKern = paramTransformPsi2(kern.comp{1}, gKern);
       %
       for i = 2:length(kern.comp)
           %
           fhandle = str2func([kern.comp{i}.type 'VardistPsi2Gradient']);
           [gKerni, gVarmeansi, gVarcovarsi, gIndi] = fhandle(kern.comp{i}, vardist, Z, covGrad, learnInducing);

           % Transformations
           gKerni = paramTransformPsi2(kern.comp{i}, gKerni);

           gVarmeans = gVarmeans + gVarmeansi;
           gVarcovars = gVarcovars + gVarcovarsi;
           gInd = gInd + gIndi;
           gKern = [gKern gKerni];
           %
       end
       %
    end

    %% compute the cross-kernel terms of Psi2
    %
    if strcmp(kern.type,'cmpnd') & (length(kern.comp)>1)
      %
      index{1}=1:kern.comp{1}.nParams;
      for i=2:length(kern.comp)
           index{i} = (index{i-1}(end)+1):index{i-1}(end)+kern.comp{i}.nParams;
      end
      %
      for i=1:length(kern.comp)
        for j=i+1:length(kern.comp)

           [gKerni, gKernj, gVarm, gVarc, gI] = kernkernVardistPsi2Gradient(kern.comp{i}, kern.comp{j}, vardist, Z, covGrad, learnInducing);

           % Transformations
           gKerni = paramTransformPsi2(kern.comp{i}, gKerni);
           gKernj = paramTransformPsi2(kern.comp{j}, gKernj);

           %index{i}
           %index{j}
           %gKerni
           %gKernj
           %gKern
           %pause

           gKern(index{i}) = gKern(index{i}) + gKerni;
           gKern(index{j}) = gKern(index{j}) + gKernj;

           gVarmeans = gVarmeans + gVarm;
           gVarcovars = gVarcovars + gVarc;
           gInd = gInd + gI;
           %
       end
      end
    end
end

% variational variances are positive (This should rather go to
% vargplvmLogLikeGradients)
% gVarcovars = (gVarcovars(:).*vardist.covars(:))';

%-----------------------------------------------------
% This applies transformations
% This must be done similar to kernGradient at some point
% function gKern = paramTransformPsi2(kern, gKern)
% %
% %
%
% if strcmp(kern.type,'rbfard2') | strcmp(kern.type,'rbfardjit')
%    gKern(1) = gKern(1)*kern.variance;
%    gKern(2:end) = gKern(2:end).*kern.inputScales;
% elseif strcmp(kern.type,'linard2')
%    gKern(1:end) = gKern(1:end).*kern.inputScales;
% elseif strcmp(kern.type,'bias')
%    gKern = gKern*kern.variance;
% end

function gKern = paramTransformPsi2(kern, gKern)

    fhandle = str2func([kern.type 'KernExtractParam']);
    params  = fhandle(kern);

    if length(kern.transforms) > 1
        for i=1:length(kern.transforms)
            fhandle = str2func([kern.transforms(i).type 'Transform']);
            if ~isfield(kern.transforms(i), 'transformsettings')
                gKern(kern.transforms(i).index) = gKern(kern.transforms(i).index).*fhandle(params(kern.transforms(i).index), 'gradfact');
            else
                gKern(kern.transforms(i).index) = gKern(kern.transforms(i).index).*fhandle(params(kern.transforms(i).index), 'gradfact', kern.transforms(i).transformsettings);
            end
        end
    else
        fhandle = str2func([kern.transforms.type 'Transform']);
        if ~isfield(kern.transforms, 'transformsettings')
            gKern(kern.transforms.index) = gKern(kern.transforms.index).*fhandle(params(kern.transforms.index), 'gradfact');
        else
            gKern(kern.transforms.index) = gKern(kern.transforms.index).*fhandle(params(kern.transforms.index), 'gradfact', kern.transforms.transformsettings);
        end
    end

end

author: SheffieldML

repo: vargplvm


text:
function phi = basis_mn_tet10 ( t, n, p )

%*****************************************************************************80
%
%% BASIS_MN_TET10: all bases at N points for a TET10 element.
%
%  Discussion:
%
%    The routine is given the coordinates of the vertices of a tetrahedron.
%
%    It works directly with these coordinates, and does not refer to a
%    reference element.
%
%    P1 through P4 are vertices.
%
%    P1 <= P5  <= P2
%    P2 <= P6  <= P3
%    P1 <= P7  <= P3
%    P1 <= P8  <= P4
%    P2 <= P9  <= P4
%    P3 <= P10 <= P4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Olgierd Zienkiewicz,
%    The Finite Element Method,
%    Sixth Edition,
%    Butterworth-Heinemann, 2005,
%    ISBN: 0750663200,
%    LC: TA640.2.Z54.
%
%  Parameters:
%
%    Input, real T(3,4), the coordinates of the vertices.
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real P(3,N), the points where the basis functions
%    are to be evaluated.
%
%    Output, real PHI(10,N), the value of the basis functions
%    at the evaluation points.
%
  phi_linear(1:4,1:n) = basis_mn_tet4 ( t, n, p );

  phi( 1,1:n) = ( 2.0 * phi_linear(1,1:n)  - 1.0 ) .* phi_linear(1,1:n);
  phi( 2,1:n) = ( 2.0 * phi_linear(2,1:n)  - 1.0 ) .* phi_linear(2,1:n);
  phi( 3,1:n) = ( 2.0 * phi_linear(3,1:n)  - 1.0 ) .* phi_linear(3,1:n);
  phi( 4,1:n) = ( 2.0 * phi_linear(4,1:n)  - 1.0 ) .* phi_linear(4,1:n);
  phi( 5,1:n) =   4.0 * phi_linear(1,1:n)          .* phi_linear(2,1:n);
  phi( 6,1:n) =   4.0 * phi_linear(2,1:n)          .* phi_linear(3,1:n);
  phi( 7,1:n) =   4.0 * phi_linear(1,1:n)          .* phi_linear(3,1:n);
  phi( 8,1:n) =   4.0 * phi_linear(1,1:n)          .* phi_linear(4,1:n);
  phi( 9,1:n) =   4.0 * phi_linear(2,1:n)          .* phi_linear(4,1:n);
  phi(10,1:n) =   4.0 * phi_linear(3,1:n)          .* phi_linear(4,1:n);

  return
end


text:
function genTlbx(examples)

if nargin==0
    examples = 0;
end

version = '19.1';
ml = ver('MATLAB');
ml = ml.Release(2:end-1);

%% Unpack verilog source
if ~examples
cd('../../');
cd('hdl_wa_bsp/vendor/AnalogDevices/vivado');
!find -name '*.zip' -exec sh -c 'unzip -d "${1%.*}" "$1"' _ {} \;
!find . -name "*.zip" -type f -delete
cd('../../../..');
cd('CI/scripts');
end

%%
cd(fileparts((mfilename('fullpath'))));
cd('../..');
p = pwd;
cd(fileparts((mfilename('fullpath'))));

if examples
    fid  = fopen('bsp.tmpl','r');
else
    fid  = fopen('bsp_noexamples.tmpl','r');
end
f=fread(fid,'*char')';
fclose(fid);

f = strrep(f,'__REPO-ROOT__',p);
f = strrep(f,'__VERSION__',version);
f = strrep(f,'__ML-RELEASE__',ml);

fid  = fopen('../../bsp.prj','w');
fprintf(fid,'%s',f);
fclose(fid);

cd('../..');
addpath(genpath(matlabshared.supportpkg.getSupportPackageRoot));
addpath(genpath('.'));
rmpath(genpath('.'));
if examples
    ps = {'doc','hdl_wa_bsp','hil_models','targeting_models','deps'};
else
    ps = {'doc','hdl_wa_bsp'};
end
paths = '';
for p = ps
    pp = genpath(p{:});
    ppF = pp;
    pp = pp(1:end-1);
    pp = strrep(pp,':','</matlabPath><matlabPath>');
    paths = [paths,['<matlabPath>',pp,'</matlabPath>']]; %#ok<AGROW>
    addpath(ppF);
end
rehash
projectFile = 'bsp.prj';
currentVersion = matlab.addons.toolbox.toolboxVersion(projectFile);
if examples
    outputFile = ['AnalogDevicesBSP_v',currentVersion];
else
    outputFile = ['AnalogDevicesBSP_noexamples_v',currentVersion];
end
matlab.addons.toolbox.packageToolbox(projectFile,outputFile)

if ~usejava('desktop')
%% Update toolbox paths
mkdir other
movefile([outputFile,'.mltbx'], ['other/',outputFile,'.zip']);
cd other
unzip([outputFile,'.zip'],'out');
cd('out')
cd('metadata');
fid  = fopen('configuration.xml','r');
f=fread(fid,'*char')';
fclose(fid);

s = '</matlabPaths>';
sections = strsplit(f,s);
s1 = sections{1};
s2 = sections{2};
newfile = [s1,paths,s,s2];

fid  = fopen('configuration.xml','w');
fprintf(fid,'%s',newfile);
fclose(fid);

%% Repack
cd('..');
zip([outputFile,'.zip'], '*');
movefile([outputFile,'.zip'],['../../',outputFile,'.mltbx']);
cd('../..');
rmdir('other','s');
end

delete bsp.prj

author: analogdevicesinc


text:
function p04_title ( )

%*****************************************************************************80
%
%% P04_TITLE prints a title for problem 04.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Problem 04\n' );
  fprintf ( 1, '  Name:       HexSum\n' );
  fprintf ( 1, '              Davis, Rabinowitz, page 370, #2.\n' );
  fprintf ( 1, '  Region:     0 <= X(i) <= 1\n' );
  fprintf ( 1, '  Integrand:  F(X) = ( sum ( 2 * X(i) - 1 ) )^6\n' );

  return
end


text:
function F  = obj_actuator_freq_resp_param(x,amplitude_40, frequency_40, amplitude_100, frequency_100)
% Computes objective function for determination of basic parameters of
% the proportional valve actuator that match the required frequency
% response. The frequency response of a nonlinear system is obtained by
% using frestimate.  The required and actual frequency responses are
% compared.
% Copyright 2010 MathWorks, Inc.

% Manufacturer's characterstic
%RefFR_100_Frq =     [7 10 20 30 40 43 50 54];
%RefFR_100_Phs = -1*[25 35 59 75 87 90 96 100];
%RefFR_40_Frq =     [7 10 20 30 40 50 57 70];
%RefFR_40_Phs = -1*[21 30 50 63 75 85 90 100];

model = 'actuator_freq_resp';
load_system(model);

assignin('base','act_gain', x(1));
assignin('base','time_const', x(2));
assignin('base','act_saturation', x(3));

% Computing phase angle at current values of variable parameters

% Determine inputs and outputs
ios = getlinio(model);

% Generate input at test frequencies
in100 = frest.Sinestream(...
   'Frequency',frequency_100, ...
   'Amplitude', amplitude_100*ones(size(frequency_100)), ...
   'SimulationOrder', 'Sequential', ...
   'FreqUnits', 'Hz');

in40 = frest.Sinestream(...
   'Frequency',frequency_40, ...
   'Amplitude', amplitude_40*ones(size(frequency_40)), ...
   'SimulationOrder', 'Sequential', ...
   'FreqUnits', 'Hz');

% Test model, determine phase
sys100 = frestimate(model,ios,in100);
R100 = sys100.ResponseData(:);
Lag = angle(R100)/pi*180;
phase_100_20 = Lag(1); phase_100_43 = Lag(2);

% Test model, determine phase
sys40 = frestimate(model,ios,in40);
R40 = sys40.ResponseData(:);
Lag = angle(R40)/pi*180;
phase_40_20 = Lag(1); phase_40_57 = Lag(2);

% Calculate difference from manufacturer's characteristic
F = (phase_100_20 + 59)^2 + (phase_100_43 + 90)^2 + ...
    (phase_40_20 + 50)^2 + (phase_40_57 + 90)^2;

end

author: Sable

repo: mcbench-benchmarks


text:
function [engine, loglik] = enter_evidence(engine, evidence, varargin)
% ENTER_EVIDENCE Add the specified evidence to the network (hmm)
% [engine, loglik] = enter_evidence(engine, evidence, ...)
%
% evidence{i,t} = [] if if X(i,t) is hidden, and otherwise contains its observed value (scalar or column vector)
%
% The following optional arguments can be specified in the form of name/value pairs:
% [default value in brackets]
%
% maximize - if 1, does max-product (not yet supported), else sum-product [0]
% filter   - if 1, does filtering, else smoothing [0]
% oneslice - 1 means only compute marginals on nodes within a single slice [0]
%
% e.g., engine = enter_evidence(engine, ev, 'maximize', 1)

maximize = 0;
filter = 0;
oneslice = 0;

% parse optional params
args = varargin;
nargs = length(args);
if nargs > 0
  for i=1:2:nargs
    switch args{i},
     case 'maximize', maximize = args{i+1};
     case 'filter',  filter = args{i+1};
     case 'oneslice', oneslice = args{i+1};
     otherwise,
      error(['invalid argument name ' args{i}]);
    end
  end
end

[ss T] = size(evidence);
engine.maximize = maximize;
engine.evidence = evidence;
bnet = bnet_from_engine(engine);
engine.node_sizes = repmat(bnet.node_sizes_slice(:), [1 T]);

obs_bitv = ~isemptycell(evidence(:));
bitv = reshape(obs_bitv, ss, T);
for t=1:T
  onodes = find(bitv(:,t));
  if ~isequal(onodes, bnet.observed(:))
    error(['dbn was created assuming observed nodes per slice were '...
	   num2str(bnet.observed(:)')  ' but the evidence in slice ' num2str(t) ...
	   ' has observed nodes ' num2str(onodes(:)')]);
  end
end

obslik = mk_hmm_obs_lik_matrix(engine, evidence);

%[alpha, beta, gamma, loglik, xi] = fwdback(engine.startprob, engine.transprob, obslik, ...
[alpha, beta, gamma, loglik, xi] = fwdback_twoslice(engine, engine.startprob,...
                                                    engine.transprob, obslik, ...
                                                    'maximize', maximize, 'fwd_only', filter, ...
                                                    'compute_xi', ~oneslice);

engine.one_slice_marginal = gamma; % gamma(:,t) for t=1:T
if ~oneslice
  Q = size(gamma,1);
  engine.two_slice_marginal = reshape(xi, [Q*Q T-1]); % xi(:,t) for t=1:T-1
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function r = isfinite(a)
%ISFINITE     Array of 1's for finite components
%
%   r = isfinite(a)
%

% written  08/07/02     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if a.complex
    r = isfinite(a.mid) & isfinite(a.rad) ;
  else
    r = isfinite(a.inf) & isfinite(a.sup) ;
  end

repo: OpenMAS


text:
function unprocessed = find_unprocessed(dag)
% unprocessed = find_unprocessed(dag)
%
% francois.olivier.c.h@gmail.com

undirected_edges = dag.*dag';
[unprocessed, tmp] = find(undirected_edges)
unprocessed = unique(unprocessed);
%  N = size(dag,1);
%  unprocessed = [];
%  for i=1:(N-1)
%      for j=(i+1):N
%          if dag(i,j)==1 & dag(j,i)==1
%             unprocessed = [unprocessed,i,j];
%          end
%      end
%  end
%  unprocessed = unique(unprocessed);


text:
classdef TestBitwiseAnd
    %TestBitwiseAnd

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','img001.jpg');
    end

    methods (Static)
        function test_rgb_images
            img1 = cv.imread(TestBitwiseAnd.im, 'ReduceScale',2);
            img2 = cv.cvtColor(img1, 'RGB2HSV');

            % rectangular mask
            [h,w,~] = size(img1);
            mask = false([h w]);
            mask(100:h-100,100:w-100) = true;

            out = cv.bitwise_and(img1, img2);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_and(img1, img2);
            assert(isequal(out, expected));

            out = cv.bitwise_and(img1, img2, 'Mask',mask);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_and(img1, img2, mask);
            assert(isequal(out, expected));

            out = cv.bitwise_and(img1, img2, 'Mask',mask, 'Dest',img1);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_and(img1, img2, mask, img1);
            assert(isequal(out, expected));
        end

        function test_float_images
            img1 = cv.imread(TestBitwiseAnd.im, 'ReduceScale',2);
            img1 = single(img1) ./ 255;
            img2 = cv.cvtColor(img1, 'RGB2HSV');
            img2(:,:,1) = img2(:,:,1) ./ 360;

            % circular mask
            [h,w,~] = size(img1);
            [X,Y] = meshgrid(1:w,1:h);
            c = fix([w h]/2); r = 50;
            mask = ((X-c(1)).^2 + (Y-c(2)).^2) < r^2;

            out = cv.bitwise_and(img1, img2);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_and(img1, img2);
            assert(isequaln(out, expected));

            out = cv.bitwise_and(img1, img2, 'Mask',mask);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_and(img1, img2, mask);
            assert(isequaln(out, expected));

            out = cv.bitwise_and(img1, img2, 'Mask',mask, 'Dest',img1);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_and(img1, img2, mask, img1);
            assert(isequaln(out, expected));
        end

        function test_vectors
            A = uint8([250 253 255]);
            B = uint8([15 15 240]);
            mask = [true false true];

            % uint8([10 13 240])
            out = cv.bitwise_and(A, B);
            validateattributes(out, {class(A)}, {'size',size(A)});
            expected = bitand(A, B);
            assert(isequal(out, expected));

            % uint8([10 0 240])
            out = cv.bitwise_and(A, B, 'Mask',mask);
            validateattributes(out, {class(A)}, {'size',size(A)});
            expected = bitand(A, B);
            expected(~mask) = 0;
            assert(isequal(out, expected));

            % uint8([10 253 240])
            out = cv.bitwise_and(A, B, 'Mask',mask, 'Dest',A);
            validateattributes(out, {class(A)}, {'size',size(A)});
            expected = bitand(A, B);
            expected(~mask) = A(~mask);
            assert(isequal(out, expected));
        end

        function test_error_argnum
            try
                cv.bitwise_and();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

function out = my_bitwise_and(src1, src2, mask, dst)
    %MY_BITWISE_AND  Similar to cv.bitwise_and using core MATLAB functions

    if nargin < 3, mask = true(size(src1,1), size(src1,2)); end
    if nargin < 4, dst = zeros(size(src1), class(src1)); end

    % calculate bitwise AND
    if isinteger(src1)
        out = bitand(src1, src2);
    elseif isfloat(src1)
        if isa(src1, 'double')
            klass = 'uint64';
        elseif isa(src1, 'single')
            klass = 'uint32';
        end
        out = bitand(typecast(src1(:), klass), typecast(src2(:), klass));
        out = reshape(typecast(out, class(src1)), size(src1));
    end

    % apply masking
    for k=1:size(out,3)
        out_slice = out(:,:,k);
        dst_slice = dst(:,:,k);
        out_slice(~mask) = dst_slice(~mask);
        out(:,:,k) = out_slice;
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function mse = eval_MTL_mse (Y, X, W)
%% FUNCTION eval_MTL_mse
%   computation of mean squared error given a specific model.
%   the value is the lower the better.
%
%% FORMULATION
%
%  multi-task mse = sum_t (mse(t) * N_t) / sum_t N_t
%
%  where
%     mse(t) = sum((Yt_pred - Y{t})^2))/ N_t
%     Yt_pred = X{t} * W(:, t)
%     N_t     = length(Y{t})
%
%  which can be simplified as:
%
%  multi-task mse = sum_t sum((Yt_pred - Y{t})^2)) / sum_t N_t
%
%
%% INPUT
%   X: {n * d} * t - input matrix
%   Y: {n * 1} * t - output matrix
%   percent: percentage of the splitting range (0, 1)
%
%% OUTPUT
%   X_sel: the split of X that has the specifid percent of samples
%   Y_sel: the split of Y that has the specifid percent of samples
%   X_res: the split of X that has the 1-percent of samples
%   Y_res: the split of Y that has the 1-percent of samples
%   selIdx: the selection index of for X_sel and Y_sel for each task
%
%% LICENSE
%   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/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on Jan 6, 2016.
%
    task_num = length(X);
    mse = 0;

    total_sample = 0;
    for t = 1: task_num
        y_pred = X{t} * W(:, t);
        %mse = mse + (sum((y_pred - Y{t}).^2)/length(y_pred)) * length(y_pred);
        mse = mse + sum((y_pred - Y{t}).^2); % length of y cancelled.
        total_sample = total_sample + length(y_pred);
    end
    mse = mse/total_sample;
end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function layout=cosmo_meeg_read_layout(fn)
% Read FieldTrip layout
%
% layout=cosmo_meeg_read_layout(fn)
%
% Inputs:
%   fn                  Filename of layout file, or a string containing the
%                       layout. In the latter case fn must contain at least
%                       one newline ('\n') character.
%                       A layout file is a text file with one line per
%                       sensor, with each line containing the following
%                       data separated by white-space:
%                       1) sensor number (integer)
%                       2) horizontal position (float)
%                       3) vertical position (float)
%                       4) width (float)
%                       5) height (float)
%                       6) label (string)
%
% Output:
%   layout              struct with fields containing data for N sensors:
%     .pos              Nx2 matrix with x and y position
%     .width            Nx1 vector
%     .height           Nx1 vector
%     .label            Nx1 cell string with channel labels
%
% Notes:
%   - whitespace is trimmed from the labels
%   - the sensor number is not used; the order of the sensors in layout is
%     the same as in the layout file
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    check_layout_input(fn)

    if string_contains_newline(fn)
        lay_string=fn;
        fn_descr=@()sprintf('input:\n''%s''',lay_string);
    else
        lay_string=read_lay_string_from_file(fn);
        fn_descr=@()sprintf('file %s',fn);
    end

    layout=parse_layout(lay_string, fn_descr);

function lay_string=read_lay_string_from_file(fn)
    if ~exist(fn,'file')
        error('layout file %s does ont exist', fn);
    end

    % read FT layout (.lay) file
    fid=fopen(fn);
    file_closer=onCleanup(@()fclose(fid));
    lay_string=fread(fid,inf,'char=>char')';

function check_layout_input(fn)
    if ~ischar(fn)
        error('first argument must be string, found %s', class(fn));
    end

function tf=string_contains_newline(fn)
    tf=any(fn==sprintf('\n'));

function layout=parse_layout(lay_string, fn_descr)
    % pattern to match is integer, then 4 numeric values followed by a
    % string that can contain whitespaces and plus characters, followed by
    % newline
    integer='(\d+)';
    float='([\d\.-]+)';
    space='\s+';
    channel_label='([\w \t\r\f\v\+\-]+)';
    single_newline='\n';

    pat=[integer, space, ...
         float, space, ...
         float, space, ...
         float, space, ...
         float, space, ...
         channel_label, single_newline];

    matches=regexp(sprintf('%s\n',lay_string),pat,'tokens');
    if isempty(matches)
        error('No valid layout definition found in %s', fn_descr());
    end

    % convert to (nchannel x 6) matrix
    layout_matrix=cat(1,matches{:});

    % convert values in first five columns to numeric
    num_values_cell=layout_matrix(:,1:5)';

    str_values=sprintf('%s %s %s %s %s; ', num_values_cell{:});
    num_values=str2num(str_values);

    % store layout information (omit channel number in first column)
    layout.pos    = num_values(:,2:3);
    layout.width  = num_values(:,4);
    layout.height = num_values(:,5);

    % trim whitespace around channel names
    label=layout_matrix(:,6);
    label=regexprep(label,'^\s*','');
    label=regexprep(label,'\s*$','');
    layout.label  = label;

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [xi,w]=thirdOrder2DCubPoints()
%%THIRDORDER2DCUBPOINTS Generate third-order cubature points
%               for integration over a 2-dimensional cube with bounds of
%               (-1,-1), (-1,1), (1,1), and (1,-1).
%
%INPUTS: None
%
%OUTPUTS: xi This is a 2XnumCubPoints set of points for the standard
%            square.
%          w A 1XnumCubPoints set of cubature weights. This sums to the
%            volume of the standard square (4).
%
%This function implements the points given in [1] (4 points).
%
%EXAMPLE:
%We compare a 2nd-order moment computed using these cubature points
%to one computed using monomialIntCube (a 3rd order moment would have just
%been 0). The results are the same within typical finite precision limits.
% [xi,w]=thirdOrder2DCubPoints();
% alpha=[2;0];
% theMoment=findMomentFromSamp(alpha,xi,w);
% intVal=monomialIntCube(alpha);
% RelErr=(theMoment-intVal)/intVal
%
%REFERENCES:
%[1] F. D. Witherden and P. E. Vincent, "On the identification of symmetric
%    quadrature rules for finite element methods," Computer and Mathematics
%    with Applications, vol. 69, no. 10, pp. 1232-1241, May 2015.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

M=[0.57735026918962576450914878050195745565,   0.57735026918962576450914878050195745565, 1;
   0.57735026918962576450914878050195745565,  -0.57735026918962576450914878050195745565, 1;
  -0.57735026918962576450914878050195745565,   0.57735026918962576450914878050195745565, 1;
  -0.57735026918962576450914878050195745565,  -0.57735026918962576450914878050195745565, 1];
w=M(:,3);
xi=M(:,1:2)';

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
classdef TestAdaptiveThreshold
    %TestAdaptiveThreshold

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','sudoku.jpg');
    end

    methods (Static)
        function test_1
            img = cv.imread(TestAdaptiveThreshold.im, ...
                'Grayscale',true, 'ReduceScale',2);
            out = cv.adaptiveThreshold(img);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_2
            img = cv.imread(TestAdaptiveThreshold.im, ...
                'Grayscale',true, 'ReduceScale',2);
            out = cv.adaptiveThreshold(img, 'Method','Gaussian', ...
                'Type','BinaryInv', 'BlockSize',7, 'C',1, 'MaxValue',255);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_error_argnum
            try
                cv.adaptiveThreshold();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end


text:
function varargout = nifti_stats(varargin)
% Conversion among various statistics
% FORMAT P = nifti_stats(VAL,CODE,OPT,PARAM)
%   CODE can be one of
%     'CORREL'      'TTEST'       'FTEST'       'ZSCORE'
%     'CHISQ'       'BETA'        'BINOM'       'GAMMA'
%     'POISSON'     'NORMAL'      'FTEST_NONC'  'CHISQ_NONC'
%     'LOGISTIC'    'LAPLACE'     'UNIFORM'     'TTEST_NONC'
%     'WEIBULL'     'CHI'         'INVGAUSS'    'EXTVAL'
%     'PVAL'
%   With only one input argument, CODE defaults to 'ZSCORE'
%
%   OPT can be one of
%     '-p' ==> output P = Prob(statistic < VAL).
%     '-q' ==> output is 1-p.
%     '-d' ==> output is probability density.
%     '-1' ==> output is X such that Prob(statistic < x) = VAL.
%     '-z' ==> output is Z such that Normal cdf(Z) = p(VAL).
%     '-h' ==> output is Z such that 1/2-Normal cdf(Z) = p(VAL).
%   With less than three input arguments, OPT defaults to '-p'.
%
%   PARAM are up to three distribution parameters.
%   These default to zero if unspecified.
%
%   P is an array with the same dimensions as VAL.
%
%__________________________________________________________________________
% 99.99% of the work by RW Cox - SSCC/NIMH/NIH/DHHS/USA/EARTH - March 2004
%  0.01% of the work (the mex wrapper) by John Ashburner - FIL/ION/UCL
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

%
% $Id: nifti_stats.m 7147 2017-08-03 14:07:01Z spm $

fprintf('******************************************\n');
fprintf('Compile the nifti_stats function with\n');
fprintf('    mex nifti_stats.c nifti_stats_mex.c -O\n');
fprintf('******************************************\n');

error('nifti_stats is not compiled.');

author: spm


text:
function v = thickVehicle(size)

%THICKVEHICLE Create a vehicle graphics 3D object.
%
%   THICKVEHICLE creates a thick triangle representing a vehicle
%   This triangle is heading towards vehicle principal axis
%   and is located at vehicle position. Initial configuration
%   is: position at origin and principal axis aligned with
%   world x axis.
%
%   THICKVEHICLE(SIZE) allows for choosing a vehicle size.
%   Default is 0.5

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if nargin == 0
    size = 1;
end

v.vert0 = [ 1.0  0.0 0.0
           -0.5  0.5 0.0
           -0.5 -0.5 0.0
            1.0  0.0 0.5
           -0.5  0.5 0.5
           -0.5 -0.5 0.5] * size/1.5;

v.vert = v.vert0;
v.faces = [1 2 3 1
           4 5 6 4
           1 2 5 4
           1 3 6 4
           2 3 6 5];

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
classdef DebugLevel
    properties
        possibleValues
        value
    end

    methods

        % ----------------------------------------------------------------------------------
        function obj = DebugLevel(value) %#ok<*INUSD>
            obj.possibleValues = struct('none',0,'simulateBadData',1);
            switch(lower(value))
                case 'none'
                    obj.value = obj.possibleValues.none;
                case 'simulatebaddata'
                    obj.value = obj.possibleValues.simulateBadData;
                otherwise
                    obj.value = obj.possibleValues.none;
            end
        end

        % ----------------------------------------------------------------------------------
        function value = Get(obj)
            value = obj.value;
        end

        % ----------------------------------------------------------------------------------
        function value = SimulateBadData(obj)
            value = obj.possibleValues.simulateBadData;
        end

        % ----------------------------------------------------------------------------------
        function value = None(obj)
            value = obj.possibleValues.None;
        end

    end
end

author: BUNPC


text:
2     % problem
4     % grid parameter
1     % uniform/stretched/shishkin grid
1     % outflow boundary
0.005 % viscosity parameter
inf   % SUPG parameter (inf==>optimal)

%% Data file for test problem CD2

author: johannesgerer

repo: jburkardt-m


text:
function S = cylinderSurfaceArea(cyl)
%CYLINDERSURFACEAREA  Surface area of a cylinder.
%
%   S = cylinderSurfaceArea(CYL)
%   Computes the surface area of the cylinder defined by:
%   CYL = [X1 Y1 Z1  X2 Y2 Z2  R],
%   where [X1 Y1 Z1] and [X2 Y2 Z2] are the coordinates of the cylinder
%   extremities, and R is the cylinder radius.
%   The surface area of the cylinder comprises the surface area of the two
%   disk-shape end caps.
%
%   Example
%     cyl = [0 0 0  1 0 0  1];
%     cylinderSurfaceArea(cyl)
%     ans =
%        12.5664
%     % equals to 4*pi
%
%   See also
%     geom3d, ellipsoidSurfaceArea, intersectLineCylinder

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2017-11-02, using Matlab 9.3.0.713579 (R2017b)
% Copyright 2017-2022 INRA - Cepia Software Platform

H = distancePoints3d(cyl(:, 1:3), cyl(:, 4:6));
R = cyl(:,7);

S1 = 2*pi*R .* H;
S2 = 2 * (pi * R.^2);

S = S1 + S2;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function F = atand(F, varargin)
%ATAND   Inverse tangent of a CHEBFUN, result in degrees.
%   ATAND(F) computes the inverse tangent (in degrees) of the CHEBFUN F.
%
%   ATAND(F, PREF) does the same but uses the CHEBFUNPREF object PREF when
%   computing the composition.
%
% See also TAND, ATAN2D, ATAN.

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

% Call the compose method:
F = compose(F, @atand, varargin{:});

end

author: chebfun


text:
function jewelersLoupe(hFig)
%JEWELERSLOUPE modifies a figure to enable a Jeweler's Loupe tool.
%   JEWELERSLOUPE(FIG) will modify the figure with handle FIG to allow a
%   jeweler's loupe tool, which magnifies the axes by a factor of 4, to
%   appear whenever a user clicks on an axes. Dragging will cause the loupe
%   to follow the cursor. The tool assumes that any axes clicked will be a
%   linear 2-D axes. Performance will degrade as the amount of data in the
%   axes increates. If no figure is given, the loupe will appear in the
%   current figure.
%
%   Example:
%   plot(magic(3));
%   jewelersLoupe(gcf);
%
%   See also PAN, ZOOM.

%   Copyright 2008 The MathWorks, Inc.

% If a figure has not been passed in,
if nargin < 1
    hFig = gcf;
end

% Make sure we have a valid figure:
if ~ishandle(hFig) || ~strcmpi(get(hFig,'Type'),'figure')
    error('JEWELERSLOUPE:INVALIDFIGURE',...
        'The first input argument must be a valid figure.');
end

% Set up the callback for the loupe:
set(hFig,'WindowButtonDownFcn',@localCreateLoupe);

%------------------------------------------------------------------------%
function localCreateLoupe(hFig,evd) %#ok<INUSD>
% The loupe will appear for the current axes.
hAx = get(hFig,'CurrentAxes');

% If the axes is empty, we won't do anything.
if isempty(hAx)
    return;
end

% The loupe is implemented in terms of a separate axes:
hLoupe = axes('Parent',hFig,'Tag','loupe','HandleVisibility','off',...
    'XTickLabel','','YTickLabel','','XTick',[],'YTick',[],'Box','on');
% Set up the initial loupe position:
localUpdateLoupePosition(hFig,hLoupe);

% Copy the contents of the axes into the loupe:
copyobj(get(hAx,'Children'),hLoupe);

% We will zoom by a factor of 4 into the loupe based on the current point
% in the axes.
localDoZoom(hLoupe,hAx);

% Set callbacks for motion and button up:
set(hFig,'WindowButtonMotionFcn',{@localMoveLoupe,hAx,hLoupe});
set(hFig,'WindowButtonUpFcn',{@localCleanUp,hLoupe});

%------------------------------------------------------------------------%
function localMoveLoupe(hFig,evd,hAx,hLoupe) %#ok<INUSL>
% Move the loupe based on the current point:
localUpdateLoupePosition(hFig,hLoupe);

% Perform the zoom operation on the axes:
localDoZoom(hLoupe,hAx);

%------------------------------------------------------------------------%
function localCleanUp(hFig,evd,hLoupe) %#ok<INUSL>
% Clean up after the loupe has finished execution:

% Remove the callbacks:
set(hFig,'WindowButtonMotionFcn','');
set(hFig,'WindowButtonUpFcn','');

% Delete the loupe:
delete(hLoupe);

%------------------------------------------------------------------------%
function localUpdateLoupePosition(hFig,hLoupe)
% Moves the loupe based on the current point in the figure:

% The loupe should be centered on the current location and take up 25% of
% the figure space.
oldUnits = get(hFig,'Units');
set(hFig,'Units','Normalized');
mousePoint = get(hFig,'CurrentPoint');
set(hFig,'Units',oldUnits);
loupePosition = [(mousePoint-.125) .25 .25];
set(hLoupe,'Position',loupePosition);

%------------------------------------------------------------------------%
function localDoZoom(hLoupe,hAx)
% First, calculate the center based on the current point in the axes:
newCenter = get(hAx,'CurrentPoint');
% We are assuming a 2-D axes, so we will take the first two coordinates.
newCenter = newCenter(1,1:2);

% Zoom in by a factor of 4 around the given center point.
center_x = newCenter(1);
center_y = newCenter(2);

% Start by zooming the X-limits:
origXLim = get(hAx,'XLim');
xmin = origXLim(1);
xmax = origXLim(2);
dx = diff(origXLim);
newdx = dx * (1/4);
newdx = min(newdx,xmax-xmin);
newXLim = [center_x-newdx/2 center_x+newdx/2];

% Next, zoom the Y-limits:
origYLim = get(hAx,'YLim');
ymin = origYLim(1);
ymax = origYLim(2);
dy = diff(origYLim);
newdy = dy * (1/4);
newdy = min(newdy,ymax-ymin);
newYLim = [center_y-newdy/2 center_y+newdy/2];

% Update the limits:
set(hLoupe,'XLim',newXLim,'YLim',newYLim);


text:
function varargout = cvtIterate(germs, funcPtr, funcArgs, N)
%CVTITERATE Update germs of a CVT using random points with given density
%
%   G2 = cvtIterate(G, FPTR, FARGS, N)
%   G: inital germs
%   FPTR: pointer to a function which accept a scalar M and return M random
%       points with a given distribution
%   FARGS: arguments to be given to the FPTR function (can be empty)
%   N: number of random points to generate
%
%   Example
%   P = randPointDiscUnif(50);
%   P2 = cvtIterate(P, @randPointDiscUnif, [], 1000);
%   P3 = cvtIterate(P2, @randPointDiscUnif, [], 1000);
%
%   See also
%
%
%   Rewritten from programs found in
%   http://people.scs.fsu.edu/~burkardt/m_src/cvt/cvt.html
%
%  Reference:
%    Qiang Du, Vance Faber, and Max Gunzburger,
%    Centroidal Voronoi Tessellations: Applications and Algorithms,
%    SIAM Review, Volume 41, 1999, pages 637-676.
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2007-10-10,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

%% Init

% format input
if isempty(funcArgs)
    funcArgs = {};
end

% number of germs
Ng = size(germs, 1);

% initialize centroids with values of germs
centroids = germs;

% number of updates of each centroid
count = ones(Ng, 1);

%% random points

% generate N random points
pts = feval(funcPtr, N, funcArgs{:});

% for each point, determines which germ is the closest ones
[dist, ind] = minDistancePoints(pts, germs); %#ok<ASGLU>

h = zeros(Ng, 1);
for i = 1:Ng
    h(i) = sum(ind==i);
end

%% Centroids update

% add coordinate of each point to closest centroid
energy = 0;
for j = 1:N
    centroids(ind(j), :) = centroids(ind(j), :) + pts(j, :);
    energy = energy + sum ( ( centroids(ind(j), :) - pts(j, :) ).^2);
    count(ind(j)) = count(ind(j)) + 1;
end

% estimate coordinate by dividing by number of counts
centroids = centroids ./ repmat(count, 1, size(germs, 2));

% normalizes energy by number of sample points
energy = energy / N;

%% format output

varargout{1} = centroids;
if nargout > 1
    varargout{2} = energy;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function imOut = visPatches(W, invXForm, transIms, cLims, borderPix)
% imOut = visZCAWhitenedPatches(W, invXForm, transIms, cLims, borderPix)
%-------------------------------------------------------------------

if notDefined('invXForm'); invXForm = eye(size(W,1)); end
if notDefined('transIms'); transIms = 0; end
if notDefined('borderPix'); borderPix = 1; end

W = invXForm * W;

minW=min(W(:));
maxW=max(W(:));

[nDim,nUnits]=size(W);

nDPix = floor(sqrt(nDim)+0.999);
nUPix = floor(sqrt(nUnits)+0.999);

imOut = zeros(((nDPix+borderPix)*nUPix+borderPix));

if (nUnits/nUPix<=nUPix-1),
    imOut = imOut(:,1:(nDPix+borderPix)*(nUPix-1)+borderPix);
end

scale = 127/max(abs(minW), abs(maxW));

try
	for iW=1:nUnits;
	    if (transIms)
		    imOut(mod(iW-1,nUPix)*(nDPix+borderPix) + ...
		    borderPix+1:mod(iW-1,nUPix)*(nDPix+borderPix) + ...
		    borderPix+nDPix,floor((iW-1)/nUPix)*(nDPix+borderPix)+borderPix + ...
		    1:floor((iW-1)/nUPix)*(nDPix+borderPix)+borderPix+nDPix)...
		    = (reshape(W(:,iW),nDPix,nDPix)'*scale + 128);
	    else
		    imOut(mod(iW-1,nUPix)*(nDPix+borderPix) + ...
		    borderPix+1:mod(iW-1,nUPix)*(nDPix+borderPix) + ...
		    borderPix+nDPix,floor((iW-1)/nUPix)*(nDPix+borderPix)+borderPix + ...
		    1:floor((iW-1)/nUPix)*(nDPix+borderPix)+borderPix+nDPix)...
		    = reshape(W(:,iW),nDPix,nDPix)*scale + 128;
	    end
	end;
catch % IF ALL ELSE FAILS, JUST SHOW INVERSE-TRANSFORMED WEIGHTS
	fprintf('\nVisualization failed, displaying weights...\n')
	imOut = W;
end
if notDefined('cLims')cLims = [minW,maxW]; end

if ~nargout
	try
		imagesc(imOut/255); colormap(gray);
%  		set(gca,'clim',cLims)
		axis image;
		axis off;
	catch
	end

end

author: dustinstansbury

repo: medal