text:
function [W,H]=nmfprob(X,K,maxiter,speak)
%
% Probabilistic NFM interpretating X as samples from a multinomial
%
% INPUT:
% X (N,M) : N (dimensionallity) x M (samples) non negative input matrix
% K       : Number of components
% maxiter : Maximum number of iterations to run
% speak   : prints iteration count and changes in connectivity matrix
%           elements unless speak is 0
%
% OUTPUT:
% W       : N x K matrix
% H       : K x M matrix
%
% Lars Kai Hansen, IMM-DTU (c) November 2005
%
print_iter=50;

powers=1.5+(2.5-1.5)*((1:maxiter)-1)/(maxiter-1);
% INITIALIZE
[D,N]=size(X);
X_factor = (sum(sum(X)));
X_org = X;
X=X/X_factor;

W=rand(D,K);
W=W./repmat(sum(W,1),D,1);
H=rand(K,N);
H=H./repmat(sum(H,2),1,N);
P=ones(K,1);
P=P/sum(P);
W1=W;H1=H;

% use W*H to test for convergence
Xr_old = W*H;

for n=1:maxiter,
    %E-step
    Qnorm=(W*diag(P))*H;

    for k=1:K,
        %E-step
        Q=(W(:,k)*H(k,:)*P(k))./(Qnorm+eps);
        XQ=X.*Q;
        %M-step W
        dummy=sum(XQ,2);
        W1(:,k)=dummy/(sum(dummy));
        dummy=sum(XQ,1);
        H1(k,:)=dummy/(sum(dummy));
    end

    W=W1;
    H=H1;

    %%%%%%%%%%%%%%%%%%%%%%%
    % print to screen
    %%%%%%%%%%%%%%%%%%%%%%%
    if (rem(n,print_iter)==0) & speak,
        Xr = W*H;
        diff = sum(sum(abs(Xr_old-Xr)));
        Xr_old = Xr;
        eucl_dist = nmf_euclidean_dist(X_org,W*diag(sqrt(P))*X_factor*diag(sqrt(P))*H);
        errorx = mean(mean(abs(X-W*H)))/mean(mean(X));
        disp(['Iter = ',int2str(n),...
            ', relative error = ',num2str(errorx),...
            ', diff = ', num2str(diff),...
            ', eucl dist ' num2str(eucl_dist)])
        if errorx < 10^(-5), break, end
    end
end,

W=W*diag(sqrt(P))*X_factor;
H=diag(sqrt(P))*H;

author: andrewssobral

repo: lrslibrary


text:
function [ftHeadmodel, ftLeadfield, iChannels] = out_fieldtrip_headmodel(HeadModelFile, ChannelFile, iChannels, isIncludeRef)
% OUT_FIELDTRIP_HEADMODEL: Converts a head model file into a FieldTrip structure (see ft_datatype_headmodel).
%
% USAGE:  [ftHeadmodel, ftLeadfield, iChannels] = out_fieldtrip_headmodel(HeadModelFile, ChannelFile, isIncludeRef=1);
%         [ftHeadmodel, ftLeadfield, iChannels] = out_fieldtrip_headmodel(HeadModelMat,  ChannelMat,  isIncludeRef=1);
%
% INPUTS:
%    - HeadModelFile  : Relative path to a head model file available in the database
%    - HeadModelMat   : Brainstorm head model file structure
%    - ChannelFile    : Relative path to a channel file available in the database
%    - ChannelMat     : Brainstorm channel file structure
% OUTPUTS:
%    - ftHeadmodel    : Volume conductor model, typically returned by ft_prepare_headmodel
%    - ftLeadfield    : Leadfield matrix, typically returned by ft_prepare_leadfield
%    - iChannels      : Modified list of channels (after adding channels)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Jeremy T. Moreau, Elizabeth Bock, Francois Tadel, 2015

% ===== PARSE INPUTS =====
if (nargin < 4) || isempty(isIncludeRef)
    isIncludeRef = 1;
end
ftHeadmodel = [];
ftLeadfield = [];

% ===== LOAD INPUTS =====
% Load head model file
if ischar(HeadModelFile)
    HeadModelMat = in_bst_headmodel(HeadModelFile);
elseif isstruct(HeadModelFile)
    HeadModelMat = HeadModelFile;
else
    error('Failed to load head model.');
end
% If this file was computed with FieldTrip, it should include the original FieldTrip headmodel
if isfield(HeadModelMat, 'ftHeadmodelMeg') && ~isempty(HeadModelMat.ftHeadmodelMeg)
    ftHeadmodel = HeadModelMat.ftHeadmodelMeg;
elseif isfield(HeadModelMat, 'ftHeadmodelEeg') && ~isempty(HeadModelMat.ftHeadmodelEeg)
    ftHeadmodel = HeadModelMat.ftHeadmodelEeg;
end
% Load channel file
if ischar(ChannelFile)
    ChannelMat = in_bst_channel(ChannelFile);
elseif isstruct(ChannelFile)
    ChannelMat = ChannelFile;
else
    error('Failed to load channel file.')
end
% Get sensor type
Modality = ChannelMat.Channel(iChannels(1)).Type;
% Get headmodel type
switch (Modality)
    case 'EEG',   HeadModelMethod = HeadModelMat.EEGMethod;
    case 'ECOG',  HeadModelMethod = HeadModelMat.ECOGMethod;
    case 'SEEG',  HeadModelMethod = HeadModelMat.SEEGMethod;
    case {'MEG','MEG MAG','MEG GRAD','MEG REF'}, HeadModelMethod = HeadModelMat.MEGMethod;
end

% ===== ADD MEG REF =====
if isIncludeRef && ismember(Modality, {'MEG','MEG MAG','MEG GRAD'})
    iRef = channel_find(ChannelMat.Channel, 'MEG REF');
    if ~isempty(iRef)
        iChannels = [iRef, iChannels];
    end
end

% ===== CREATE FIELDTRIP HEADMODEL =====
if isempty(ftHeadmodel)
    % Get subject
    sSubject = bst_get('SurfaceFile', HeadModelMat.SurfaceFile);
    % Headmodel type
    switch (HeadModelMethod)
        case {'meg_sphere', 'singlesphere'}
            ftHeadmodel.type = 'singlesphere';
            ftHeadmodel.r = HeadModelMat.Param(iChannels(1)).Radii(1);
            ftHeadmodel.o = HeadModelMat.Param(iChannels(1)).Center(:)';

        case {'eeg_3sphereberg', 'concentricspheres'}
            ftHeadmodel.type = 'concentricspheres';
            ftHeadmodel.r = HeadModelMat.Param(iChannels(1)).Radii(:)';
            ftHeadmodel.o = HeadModelMat.Param(iChannels(1)).Center(:)';
            % Get default conductivities
            BFSProperties = bst_get('BFSProperties');
            ftHeadmodel.c = BFSProperties(1:3);

        case {'os_meg', 'localspheres'}
            ftHeadmodel.type = 'localspheres';
            ftHeadmodel.r = [HeadModelMat.Param(iChannels).Radii]';
            ftHeadmodel.o = [HeadModelMat.Param(iChannels).Center]';

        case {'singleshell'}
            ftHeadmodel.type = HeadModelMethod;
            % Check if the surfaces are available
            if isempty(sSubject.iInnerSkull)
                error('No inner skull surface available for this subject.');
            else
                disp(['BST> ' HeadModelMethod ': Using the default inner skull surface available in the database.']);
            end
            % Load surfaces
            SurfaceFiles = {sSubject.Surface(sSubject.iInnerSkull).FileName};
            ftHeadmodel.bnd = out_fieldtrip_tess(SurfaceFiles);

        case {'openmeeg', 'dipoli', 'bemcp'}
            ftHeadmodel.type = HeadModelMethod;
            % Check if the surfaces are available
            if isempty(sSubject.iInnerSkull) || isempty(sSubject.iOuterSkull) || isempty(sSubject.iScalp)
                error('No BEM surfaces available for this subject.');
            else
                disp(['BST> ' HeadModelMethod ': Using the default surfaces available in the database (inner skull, outer skull, scalp).']);
            end
            % Load surfaces
            SurfaceFiles = {sSubject.Surface(sSubject.iScalp).FileName, ...
                            sSubject.Surface(sSubject.iOuterSkull).FileName, ...
                            sSubject.Surface(sSubject.iInnerSkull).FileName};
            ftHeadmodel.bnd = out_fieldtrip_tess(SurfaceFiles);
            % Default OpenMEEG options
            ftHeadmodel.cond         = [0.33, 0.004125, 0.33];
            ftHeadmodel.skin_surface = 1;
            ftHeadmodel.source       = 3;
            % ERROR: MISSING .mat field??

        otherwise
            error('out_fieldtrip_headmodel does not support converting this type of head model.');
    end
    % Unit and labels
    ftHeadmodel.unit  = 'm';
    ftHeadmodel.label = {ChannelMat.Channel(iChannels).Name}';
end

% ===== CREATE FIELDTRIP LEADFIELD =====
if (nargout >= 2)
    % Create FieldTrip structure
    nSources = length(HeadModelMat.GridLoc);
    ftLeadfield.pos             = HeadModelMat.GridLoc;
    ftLeadfield.unit            = 'm';
    ftLeadfield.inside          = true(nSources, 1);
    ftLeadfield.leadfielddimord = '{pos}_chan_ori';
    ftLeadfield.label           = {ChannelMat.Channel(iChannels).Name};
    ftLeadfield.leadfield       = cell(1, nSources);
    for i = 1:nSources
        ftLeadfield.leadfield{i} = HeadModelMat.Gain(iChannels, 3*(i-1)+[1 2 3]);
    end
end

author: brainstorm-tools


text:
function save_im(hf)

% rszs=get(hf,'ResizeFcn');
% set(hf,'ResizeFcn','');

%set(hf,'UserData',[status status2]);
ss2=get(hf,'UserData');
status=ss2(1);
status2=ss2(2);
status3=ss2(3);

ssz=get(0,'ScreenSize');

% set(status,'FontUnits','points','FontSize',8*ssz(4)/768);
% set(status2,'FontUnits','points','FontSize',8*ssz(4)/768);
% set(status3,'FontUnits','points','FontSize',8*ssz(4)/768);
fs=get(status,'FontSize');
set(status,'FontSize',round(0.7*fs));
fs=get(status2,'FontSize');
set(status2,'FontSize',round(0.7*fs));
fs=get(status3,'FontSize');
set(status3,'FontSize',round(0.7*fs));
drawnow;

[file,path] = uiputfile({'*.png'; '*.jpg'; '*.tif'; '*.bmp'; '*.emf'; '*.eps'; '*.fig '; '*.mat'},'Save Image As');
%set(hf,'PaperPositionMode','manual','PaperOrientation','portrait');
saveas(hf,[path file],'png');
%set(hf,'PaperPositionMode','manual','PaperOrientation','Landscape');

% set(hf,'ResizeFcn',rszs);

% set(status,'FontUnits','points','FontSize',10*ssz(4)/768);
% set(status2,'FontUnits','points','FontSize',10*ssz(4)/768);
% set(status3,'FontUnits','points','FontSize',10*ssz(4)/768);

set(status,'FontSize',fs);
set(status2,'FontSize',fs);
set(status3,'FontSize',fs);

author: Sable

repo: mcbench-benchmarks


text:
function Y = subsref(G, index)
%
%  Define field name indexing for kronMatrix object.
%
%               Y = subsref(A, index)
%
%          This is called whenever a subscripted reference to the
%          kronMatrix object is made.  Options are:
%              G.a
%              G.b
%              G.a{k}
%              G.b{k}
%          where k is an integer index.  That is, if A = G.a is a cell
%          array, then A{k} = G.a{k}.
%

%  J. Nagy  6/5/03

Y = G;

switch index(1).type

case '.'
  switch index(1).subs
  case 'a'
    Y = G.a;
  case 'b'
    Y = G.b;
  otherwise
    error('Subscript should be ''a'' or ''bb''')
  end

otherwise
  error('Array and cell array indexing are not supported.')
end

if length(index) == 2
  switch index(2).type
  case '{}'
    if isa(Y, 'cell')
      Y = Y{index(2).subs{1}};
    else
      if index(2).subs{1} == 1
        return
      else
        error('Index exceeds cell array dimension.')
      end
    end
  otherwise
    error('Only cell indexing is supported for factors.')
  end
elseif length(index) > 2
  error('Illegal reference.')
end

if length(Y) == 1
  Y = Y{1};
end

repo: IRtools


text:
function triangulation_order6_print ( node_num, triangle_num, node_xy, ...
  triangle_node, triangle_neighbor )

%*****************************************************************************80
%
%% TRIANGULATION_ORDER6_PRINT prints out information defining a Delaunay triangulation.
%
%  Discussion:
%
%    Triangulations created by R8TRIS2 include extra information encoded
%    in the negative values of TRIANGLE_NEIGHBOR.
%
%    Because some of the nodes counted in NODE_NUM may not actually be
%    used in the triangulation, I needed to compute the true number
%    of vertices.  I added this calculation on 13 October 2001.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input, real NODE_XY(2,NODE_NUM), the point coordinates.
%
%    Input, integer TRIANGLE_NODE(6,TRIANGLE_NUM), the nodes that make up the triangles.
%
%    Input, integer TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), the triangle neighbors on each side.
%    If there is no triangle neighbor on a particular side, the value of
%    TRIANGLE_NEIGHBOR should be negative.  If the triangulation data was created by
%    R8TRIS2, then there is more information encoded in the negative values.
%
  dim_num = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATION_ORDER6_PRINT\n' );
  fprintf ( 1, '  Information defining a triangulation.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of nodes is %d\n', node_num );

  r8mat_transpose_print ( dim_num, node_num, node_xy, '  Nodes:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of triangles is %d\n', triangle_num );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sets of six points are used as vertices of\n' );
  fprintf ( 1, '  the triangles.  For each triangle, the vertices\n' );
  fprintf ( 1, '  are listed in counterclockwise order, followed by\n' );
  fprintf ( 1, '  the midside nodes.\n' );

  i4mat_transpose_print ( 6, triangle_num, triangle_node, '  Triangle nodes:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  On each side of a given triangle, there is either\n' );
  fprintf ( 1, '  another triangle, or a piece of the convex hull.\n' );
  fprintf ( 1, '  For each triangle, we list the indices of the three\n' );
  fprintf ( 1, '  neighbors, or (if negative) the codes of the\n' );
  fprintf ( 1, '  segments of the convex hull.\n' );

  i4mat_transpose_print ( 3, triangle_num, triangle_neighbor, ...
    '  Triangle neighbors' );
%
%  Determine the number of vertices.  This is not the same as the
%  number of points!
%
  vertex_list(               1:  triangle_num) = ...
    triangle_node(1,1:triangle_num);
  vertex_list(  triangle_num+1:2*triangle_num) = ...
    triangle_node(2,1:triangle_num);
  vertex_list(2*triangle_num+1:3*triangle_num) = ...
    triangle_node(3,1:triangle_num);

  vertex_list = i4vec_sort_heap_a ( 3*triangle_num, vertex_list );

  [ vertex_num, vertex_list ] = i4vec_sorted_unique ( 3*triangle_num, ...
    vertex_list );
%
%  Determine the number of boundary points.
%
  boundary_num = 2 * vertex_num - triangle_num - 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of boundary points is %d\n', boundary_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The segments that make up the convex hull can be\n' );
  fprintf ( 1, '  determined from the negative entries of the triangle\n' );
  fprintf ( 1, '  neighbor list.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     #   Tri  Side    N1    N2    N3\n' );
  fprintf ( 1, '\n' );

  skip = 0;

  k = 0;

  for i = 1 : triangle_num

    for j = 1 : 3

      if ( triangle_neighbor(j,i) < 0 )

        s = -triangle_neighbor(j,i);
        t = floor ( s / 3 );

        if ( t < 1 | triangle_num < t )
          fprintf ( 1, '\n' );
          fprintf ( 1, '  Sorry, this data does not use the R8TRIS2\n' );
          fprintf ( 1, '  convention for convex hull segments.\n' );
          skip = 1;
          break
        end

        s = mod ( s, 3 ) + 1;
        k = k + 1;
        n1 = triangle_node(s,t);
        n2 = triangle_node(s+3,t);
        n3 = triangle_node(i4_wrap(s+1,1,3),t);
        fprintf ( 1, '  %4d  %4d  %4d  %4d  %4d  %4d\n', k, t, s, n1, n2, n3 );
      end

    end

    if ( skip )
      break
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function J = NLMF2Dtree(I, Options)
% This function NLMF2Dtree performs Non-Local Means noise filtering of
% 2D grey/color image. This function constructs a KD-tree of the whole
% image to find patches which have the closest intensity difference,
% instead of searching only in the local neighboorhood like in the
% NLMF function.
% Warning, despite the usage of a kd-tree, this function will take more
% then a hour for one image.
%
% Function:
%
%   J = NLMF2Dtree(I, Options);
%
% inputs,
%   I : 2D grey/color or 3D image data, of type Single or Double
%           in range [0..1]
%   Options : Struct with options
%
% outputs,
%   J : The NL-means filtered image or image volume
%
% options,
%   Options.kernelratio : Radius of local Patch (default 3)
%   Options.filterstrength : Strength of the NLMF filtering (default 0.05)
%   Options.number : The number of patches which have the closest distance
%                      to a local patch used for filtering (default 9).
%
% First Compile c-code!!!!, with :
%   mex image2vectors_double.c -v
%
% Note: This function uses the "Matlab Statistics Toolbox"
%

% Process inputs
defaultoptions=struct('kernelratio',3,'filterstrength',0.05,'number',9);
if(~exist('Options','var')), Options=defaultoptions;
else
    tags = fieldnames(defaultoptions);
    for i=1:length(tags), if(~isfield(Options,tags{i})), Options.(tags{i})=defaultoptions.(tags{i}); end, end
    if(length(tags)~=length(fieldnames(Options))),
        warning('NLMF2Dtree:unknownoption','unknown options found');
    end
end

% Detect if this is a color or grey scale image
switch(size(I,3))
    case 1, rgb=false;
    case 3, rgb=true;
    otherwise,error('NLMF2Dtree:inputs','This is not a 2D image');
end

% Pad the image, to allow simple extraction of local patches at the image
% boundary
Ipad = padarray(I,[Options.kernelratio Options.kernelratio],'symmetric');

% Get the local patches
V=image2vectors_double(double(Ipad),Options.kernelratio,2);

% Create a KDtree with all local patches
ns = createns(V','nsmethod','kdtree');

% Find the patches closest by in intensity in relation to the local patch
% itself
[VM,D] = knnsearch(ns,V','k',Options.number+1);

% Calculate the weight of the patches used for a pixel
W=exp(-(D.^2)*Options.filterstrength);
% The first column is the local patch itself, set the weight to maximum
% found in the other patches
W(:,1)=max(W(:,2:end),[],2);
% Make the sum of weights which will be used for one patch, equal to one.
W=W./repmat((sum(W,2)+eps),[1 size(W,2)]);

% Add the weighted intensities, of the center pixels of the in intensity
% closest by patches
if(rgb)
    Ir=I(:,:,1);
    Ig=I(:,:,2);
    Ib=I(:,:,3);
    Jr=sum(Ir(VM).*W,2);
    Jg=sum(Ig(VM).*W,2);
    Jb=sum(Ib(VM).*W,2);
    Jr=reshape(Jr,[size(I,1) size(I,2)]);
    Jg=reshape(Jg,[size(I,1) size(I,2)]);
    Jb=reshape(Jb,[size(I,1) size(I,2)]);
    J(:,:,1)=Jr;
    J(:,:,2)=Jg;
    J(:,:,3)=Jb;
else
    J=sum(I(VM).*W,2);
    J=reshape(J,size(I));
end

author: thomaskuestner


text:
classdef print

properties
end

methods (Static)

    function [ ] = sparse_signal( x )
    %sparse_signal Prints a sparse vector as pairs of indices and values
    N = length(x);
    K = 1;
    for i=1:N
        if x(i)
            fprintf('(%d,%0.4f) ',i, x(i));
            if mod(K, 5) == 0
                % We introduce a new line after every few values
                fprintf('\n');
            end
            K = K+1;
        end
    end
    fprintf('  N=%d, K=%d\n', N, K-1);
    end

    function [ ] = sorted_sparse_signal( x )
        %PRINTSORTEDSPARSEVECTOR Sorts non-zero values in x and prints them.
        % We identify non-zero values of x. We sort them. We print them in
        % the descending magnitude order along with their indices
        values = spx.commons.sparse.sorted_non_zero_elements(x);
        fprintf('Index:\tValue\n');
        fprintf('%4d:\t%f\n', values);
    end

    function vector(x, precision)
        % prints a vector
        if nargin < 2
            precision = 2;
        end
        n = numel(x);
        if (precision == 'e')
            format = '%.4e ';
        else
            format = sprintf('%%.%df ', precision);
        end
        for i=1:n
            fprintf(format, x(i));
            if mod(i, 20) == 0
                fprintf('\n');
            end
        end
        if mod(i, 40) ~= 0
            fprintf('\n');
        end
    end

end
end

author: indigits

repo: sparse-plex


text:
function I = mutualinfo1(A,B)
%==================================================
%This is a prog in the MutualInfo 0.9 package written by
% Hanchuan Peng.
%
% I = mutualinfo(vec1,vec2)
% calculate the mutual information of two vectors
%
% For small images this function is faster than
% function 'mutualinfo2'.
%==================================================

pA = estpa(A);
HA = estentropy(pA);

pB = estpa(B);
HB = estentropy(pB);

pAB = estpab(A,B);
HAB = estjointentropy(pAB);

I = HA + HB - HAB;
return

author: thfylsty


text:
function drawPolyFromVertices(c, color, varargin)

if isempty(c)
  return
end

k = convhull(c(1,:), c(2,:), c(3,:));
X = reshape(c(1,k'), size(k'));
Y = reshape(c(2,k'), size(k'));
Z = reshape(c(3,k'), size(k'));
plot3(X,Y,Z,color)
fill3(X,Y,Z,color, 'FaceAlpha', 0.5, varargin{:})

author: rdeits

repo: iris-distro

save_path: github-repos/MATLAB/rdeits-iris-distro


text:
function out = idwMvInterp(imgw, map, maxhw, p )
% Description:
% Fill holes using Inverse Distance Weighting interpolation
%
% Inputs:
% imgw - input image
% map - Map of the canvas with 0 indicating holes and 1 indicating pixel
% maxhw - Radius of inverse weighted interpolation
% p - power for inverse weighted interpolation
%
% Output:
% out - interpolated image
%
% Author: Fitzgerald J Archibald
% Date: 23-Apr-09

outH  = size(imgw,1);
outW  = size(imgw,2);
out = imgw;

[yi_arr, xi_arr] = find(map==0); % Find locations needing fill
if isempty(yi_arr) == false
    color = size(imgw,3);
    for ix = 1:length(yi_arr),

        xi = xi_arr(ix);
        yi = yi_arr(ix);

        % Find min window which has non-hole neighbors
        yixL=max(yi-maxhw,1);
        yixU=min(yi+maxhw,outH);
        xixL=max(xi-maxhw,1);
        xixU=min(xi+maxhw,outW);

        % use inverse distance weighting filter for filling
        mapw = map(yixL:yixU,xixL:xixU);
        if isempty(find(mapw, 1)) == false,
            wk = compWk(mapw, xi-xixL+1, yi-yixL+1, p); % compute weights for interpolation
            for colIx = 1:color,
                out(yi,xi,colIx) = idw(imgw(yixL:yixU, xixL:xixU, colIx), mapw, wk); % interpolation
            end
        end
    end
end

return;

% compute wk
function wk = compWk(map, cx, cy, p)

[h,w] = size(map);
[x,y] = meshgrid(1:h,1:w);
y=(y-cx).^2;
x=(x-cy).^2;
d2=x+y; % square of distance
wk = 1./(d2'.^(p/2));

return;

% Inverse distance weighting
function out = idw(in, map, wk)

num=sum(double(in(find(map))).*wk(find(map))); % weight the available pixel values among the neighbors
den=sum(wk(find(map))); % sum of contributing weights within radius r

out = num / den;

return


text:
function obox = orientedBox(points)
%ORIENTEDBOX Minimum-width oriented bounding box of a set of points
%
%   OBOX = orientedBox(PTS)
%   Computes the oriented bounding box of a set of points. Oriented box is
%   defined by a center, two dimensions (the length and the width), and the
%   orientation of the length axis. Orientation is counted in degrees,
%   counter-clockwise.
%
%   Example
%     % Draw oriented bounding box of an ellipse
%     elli = [30 40 40 20 30];
%     pts = ellipseToPolygon(elli, 120);
%     obox = orientedBox(pts);
%     figure; hold on;
%     drawEllipse(elli);
%     drawOrientedBox(obox, 'm');
%
%   See also
%   drawOrientedBox, orientedBoxToPolygon
%

% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2012-03-29,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2012 INRA - Cepia Software Platform.

%% initialisations

% first, compute convex hull of the polygon
inds = convhull(points(:,1), points(:,2));
hull = points(inds, :);

% if first and last points are the same, remove the last one
if inds(1) == inds(end)
    hull = hull(1:end-1, :);
end

% compute convex hull centroid, that corresponds to approximate
% location of rectangle center
center = mean(hull, 1);
hull = bsxfun(@minus, hull, center);

% number of hull vertices
nV = size(hull, 1);

% default values
rotatedAngle = 0;
minWidth = inf;
minAngle = 0;

% avoid degenerated cases
if nV < 3
    return;
end

% indices of vertices in extreme y directions
[tmp, indA] = min(hull(:, 2)); %#ok<ASGLU>
[tmp, indB] = max(hull(:, 2)); %#ok<ASGLU>

caliperA = [ 1 0];    % Caliper A points along the positive x-axis
caliperB = [-1 0];    % Caliper B points along the negative x-axis

%% Find the direction with minimum width (rotating caliper algorithm)

while rotatedAngle < pi
    % compute the direction vectors corresponding to each edge
    indA2 = mod(indA, nV) + 1;
    vectorA = hull(indA2, :) - hull(indA, :);

    indB2 = mod(indB, nV) + 1;
    vectorB = hull(indB2, :) - hull(indB, :);

    % Determine the angle between each caliper and the next adjacent edge
    % in the polygon
    angleA = vectorAngle(caliperA, vectorA);
    angleB = vectorAngle(caliperB, vectorB);

    % Determine the smallest of these angles
    angleIncrement = min(angleA, angleB);

    % Rotate the calipers by the smallest angle
    caliperA = rotateVector(caliperA, angleIncrement);
    caliperB = rotateVector(caliperB, angleIncrement);

    rotatedAngle = rotatedAngle + angleIncrement;

    % compute current width, and update opposite vertex
    if angleA < angleB
        line = createLine(hull(indA, :), hull(indA2, :));
        width = distancePointLine(hull(indB, :), line);
        indA = mod(indA, nV) + 1;

    else
        line = createLine(hull(indB, :), hull(indB2, :));
        width = distancePointLine(hull(indA, :), line);
        indB = mod(indB, nV) + 1;

    end

    % update minimum width and corresponding angle if needed
    if width < minWidth
        minWidth = width;
        minAngle = rotatedAngle;
    end
end

%% Compute box dimensions

% orientation of the main axis
theta = rad2deg(minAngle);

% pre-compute trigonometric functions
cot = cos(minAngle);
sit = sin(minAngle);

% elongation in direction of rectangle length
x = hull(:,1);
y = hull(:,2);
x2  =   x * cot + y * sit;
y2  = - x * sit + y * cot;

% compute extension along main directions
xmin = min(x2);    xmax = max(x2);
ymin = min(y2);    ymax = max(y2);

% position of the center with respect to the centroid compute before
dl = (xmax + xmin)/2;
dw = (ymax + ymin)/2;

% change  coordinate from rectangle to user-space
dx  = dl * cot - dw * sit;
dy  = dl * sit + dw * cot;

% coordinates of oriented box center
center = center + [dx dy];

% size of the rectangle
rectLength  = xmax - xmin;
rectWidth   = ymax - ymin;

% concatenate rectangle data
obox = [center rectLength rectWidth theta];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
close all;
clear all;
clc;

load('bin/omp_vs_gomp_comparison.mat');

mf = spx.graphics.Figures();

mf.new_figure('Iterations');

lengends = cell(5, 1);
lengends{1} = 'OMP';
for nl=1:num_ls
    lengends{nl+1} = sprintf('GOMP-%d', Ls(nl));
end
hold all;

plot(Ks, omp_average_iterations_with_k(Ks));
for nl=1:num_ls
    plot(Ks, gomp_average_iterations_with_k(Ks, nl));
end
xlabel('Sparsity level');
ylabel('Average iterations');
legend(lengends);
grid on;

mf.new_figure('Success rate');
hold all;
plot(Ks, omp_success_rates_with_k(Ks));
for nl=1:num_ls
    plot(Ks, gomp_success_rates_with_k(Ks, nl));
end
xlabel('Sparsity level');
ylabel('Success rates');
legend(lengends);
grid on;


text:
% GEODEMO_6 - Extracting a layer of velocity from a C-GRID model like ROMS
url = 'http://geoport.whoi.edu/thredds/dodsC/examples/bora_feb.nc';
hname = 'h';
uname = 'u';
vname = 'v';
aname = 'angle';

%% Use NCTOOLBOX ncgeodataset for data access
% Create a geodataset object that represents the dataset from the opendap endpoint
nc = ncgeodataset(url);

% Create variable objects that represent the variables of interest in the dataset
uvar = nc.geovariable(uname);
vvar = nc.geovariable(vname);
hvar = nc.geovariable(hname);
avar = nc.geovariable(aname);

%% Interpolate ROMS u and v vectors on to same grid as hvar (a rho variable)
% Use the uvar and vvar objects and the avar (angle) variable object
% to interpolate velocity vectors onto rho grid points and rotate them according
% to the angle variable
Uobj = hvar.getvectors(uvar, vvar, avar);

%% Grab the data at the indices of interest
itime = 3; % 3rd time step
klev = -1; % last (top) layer
% Grab coordinate information at our indices
g = Uobj.grid(itime, klev, :, :);

%% Plot vectors
% Get the magnitude of the new vectors at our indices and plot,
% overlay with the vector representation of the new velocities
figure;
pcolorjw(g.lon, g.lat, Uobj.magnitude(itime, klev, :, :));
colorbar;
arrows(g.lon(1:end,1:end), g.lat(1:end,1:end),...
    Uobj.vectors(itime, klev, 1:end,1:end), 0.08, 'black');
title({nc.attribute('title');nc.location;datestr(g.time)},'interpreter','none');
dasp(44);

%% Plot subset of vectors
figure;
pcolorjw(g.lon(1:58, 1:70), g.lat(1:58, 1:70), ...
    Uobj.magnitude(itime, klev, 1:58, 1:70));
colorbar;
arrows(g.lon(1:2:58,1:2:70), g.lat(1:2:58,1:2:70),...
    Uobj.vectors(itime, klev, 1:2:58,1:2:70), 0.08, 'black');
title({nc.attribute('title');nc.location;datestr(g.time)},'interpreter','none');
dasp(44);

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function a = c8mat_house ( n, x )

%*****************************************************************************80
%
%% C8MAT_HOUSE constructs a complex Householder elementary reflector matrix.
%
%  Formula:
%
%     A = I - ( 2 * X * hermitian ( X ) ) / ( conjg ( X ) * X )
%
%  Example:
%
%    N = 5, X = ( 1, 1, 1, 0, -1 )
%
%   1/2 -1/2 -1/2  0  1/2
%  -1/2  1/2 -1/2  0  1/2
%  -1/2 -1/2  1/2  0  1/2
%    0    0    0   1   0
%   1/2  1/2  1/2  0  1/2
%
%  Properties:
%
%    A is hermitian: hermitian ( A ) = A.
%
%    Because A is hermitian, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is unitary: hermitian ( A ) * A = A * hermitian ( A ) = I.
%
%    inverse ( A ) = A.
%
%    det ( A ) = -1.
%
%    LAMBDA(1) = -1.
%
%    If X is the vector used to define A, then X is an eigenvector
%    of A associated with the eigenvalue of -1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, complex X(N), the vector that defines the
%    Householder matrix.
%
%    Output, complex A(N,N), the matrix.
%
  a = c8mat_identity ( n );

  xdot = real ( x(1:n) * x(1:n)' );

  if ( 0.0 < xdot )

    for i = 1 : n
      for j = 1 : n
        a(i,j) = a(i,j) - 2.0 * x(i) * x(j)' / xdot;
      end
    end

  end

  return
end


function [sigma] = validateCovMatrix(sig)

% [sigma] = validateCovMatrix(sig)
%
% -- INPUT --
% sig:      sample covariance matrix
%
% -- OUTPUT --
% sigma:    positive-definite covariance matrix
%

EPS = 10^-6;
ZERO = 10^-10;

sigma = sig;
[r err] = cholcov(sigma, 0);

if (err ~= 0)
    % the covariance matrix is not positive definite!
    [v d] = eig(sigma);

    % set any of the eigenvalues that are <= 0 to some small positive value
    for n = 1:size(d,1)
        if (d(n, n) <= ZERO)
            d(n, n) = EPS;
        end
    end
    % recompose the covariance matrix, now it should be positive definite.
    sigma = v*d*v';

    [r err] = cholcov(sigma, 0);
    if (err ~= 0)
        disp('ERROR!');
    end
end


text:
function [gx,dgdx,dgdp] = g_classif0(x,P,u,in)
% logistic-like observation function for VBA_classification
% function [gx,dgdx,dgdp] = g_classif0(x,P,u,in)
% IN:
%   - x: [useless]
%   - P: classification weights
%   - u: [useless]
%   - in: contains design matrix (in.X) and sparsify flag (in.sparse)
% OUT:
%   - gx: E[y|P]
%   - dgdx: [useless]
%   - dgdp: gradient of E[y|P] wrt P.

if in.sparse
    [sP, dsdP] = VBA_sparsifyPrior(P);
else
    sP = P;
end
gx = sss(in.X'*sP);
dgdx = [];
dgdp = diag(gx.*(1-gx))*in.X';
dgdp = dgdp';
if in.sparse % for exploiting the analytical gradients from g_GLM
    dgdp = dsdP*dgdp;
end

function sx = sss(x)
sx = 1./(1+exp(-x));
sx(sx < 1e-8) = 1e-8;
sx(sx > 1-1e-8) = 1-1e-8;

author: MBB-team

repo: VBA-toolbox


text:
function [horiList, vertList, horiPos, vertPos] = createDeBruijnSeq(prjW, prjH)
%% Create the De Bruijn sequence we used in the paper.
% 1 = Red
% 2 = Yellow
% 3 = Lime
% 4 = Green
% 5 = Cyan
% 6 = Blue
% 7 = Purple
% 8 = Magenta
% See also: ImgProc.genStructuredLight

%% License
% ACADEMIC OR NON-PROFIT ORGANIZATION NONCOMMERCIAL RESEARCH USE ONLY
% Copyright (c) 2018 Bingyao Huang
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:

% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.

% If you publish results obtained using this software, please cite our paper.

% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
% SOFTWARE.

%%

HStripeT = 2;
VStripeT = 2;
TopBotbuffer = 10;
LeftRightbuffer = 10;
Debrujin_k = 4; % of color
Debrujin_n = 3; %size of window (seq length)
Debrujin_size = power(Debrujin_k, Debrujin_n)+(Debrujin_n-1); %length of Debruijn sequence, +n-1 for wrap

Hsr = TopBotbuffer; % start position of Horizontal strips
Vsr = LeftRightbuffer; % start position of Vertical strips
Hspace = floor((prjH - TopBotbuffer * 2 - HStripeT) / (Debrujin_size-1));  % has to match CreateDeBrujin
Vspace = floor((prjW - LeftRightbuffer * 2 - VStripeT) / (Debrujin_size-1));

% Create Debruijn sequence
sequence = zeros(Debrujin_size,1);

HL = 4;

digOne = HL - 1;
digTwo = HL;
i = 0;
for curHL = 1 : HL
    for curdigOne = 1:digOne
        for curdigTwo = 1:digTwo
            i = i+1;
            sequence(i) = (HL - (curdigOne - 1));
            i = i+1;
            sequence(i) = (HL - (curdigTwo - 1));
            i = i+1;
            sequence(i) = (curHL);
        end
    end

    if (curHL < HL)
        i = i+1;
        sequence(i) = (curHL);
    end

    digOne = digOne - 1;
    digTwo = digTwo - 1;

end

% Three additions of HL to the end
i = i+1;
sequence(i) = (HL);
i = i+1;
sequence(i) = (HL);
i = i+1;
sequence(i) = (HL);

horiList = zeros(Debrujin_size,1);
vertList = zeros(Debrujin_size,1);
horiPos = zeros(Debrujin_size,1);
vertPos = zeros(Debrujin_size,1);

% Create De Bruijn look up table(position and list)
for i = 1:Debrujin_size
    horiList(i) = 2 * sequence(i) - 1;
    horiPos(i)  = Hsr+2; % for 1280 and topbottom buffer = 20, space is 19.375 without HstripT
    Hsr = Hsr + Hspace;

    %Vertical
    vertList(i)  = 2 * sequence(i);
    vertPos(i)  = Vsr+2; % for 1280 and topbottom buffer = 20, space is 19.375 without HstripT
    Vsr = Vsr + Vspace;
end

horiList = horiList';
vertList = vertList';
horiPos = horiPos';
vertPos = vertPos';

end

%     horiList = [7,7,1,7,5,1,7,3,1,7,1,1,5,7,1,5,5,1,5,3,1,5,1,1,3,7,1,3,5,1,...
%         3,3,1,3,1,1,1,7,7,3,7,5,3,7,3,3,5,7,3,5,5,3,5,3,3,3,7,7,5,7,5,5,5,7];
%
%     vertList = [8,8,2,8,6,2,8,4,2,8,2,2,6,8,2,6,6,2,6,4,2,6,2,2,4,8,2,4,6,2,...
%         4,4,2,4,2,2,2,8,8,4,8,6,4,8,4,4,6,8,4,6,6,4,6,4,4,4,8,8,6,8,6,6,6,8];
%
%     horiPos = [10,21,32,43,54,65,76,87,98,109,120,131,142,153,164,175,186,197,...
%         208,219,230,241,252,263,274,285,296,307,318,329,340,351,362,373,384,395,...
%         406,417,428,439,450,461,472,483,494,505,516,527,538,549,560,571,582,593,604,...
%         615,626,637,648,659,670,681,692,703];
%
%     vertPos = [10,29,48,67,86,105,124,143,162,181,200,219,238,257,276,295,314,333,...
%         352,371,390,409,428,447,466,485,504,523,542,561,580,599,618,637,656,675,...
%         694,713,732,751,770,789,808,827,846,865,884,903,922,941,960,979,998,1017,...
%         1036,1055,1074,1093,1112,1131,1150,1169,1188,1207];

author: BingyaoHuang

repo: single-shot-pro-cam-calib


text:
function [ o, x, w ] = en_r2_03_2 ( n )

%*****************************************************************************80
%
%% EN_R2_03_2 implements the Stroud rule 3.2 for region EN_R2.
%
%  Discussion:
%
%    The rule has order O = 2^N.
%
%    The rule has precision P = 3.
%
%    EN_R2 is the entire N-dimensional space with weight function
%
%      w(x) = exp ( - x1^2 - x2^2 ... - xn^2 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971,
%    ISBN: 0130438936,
%    LC: QA311.S85.
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Output, integer O, the order.
%
%    Output, real X(N,O), the abscissas.
%
%    Output, real W(O), the weights.
%
  o = 2^n;
  volume = sqrt ( pi^n );

  a = volume / o;
  r = sqrt ( 1 / 2 );

  x = zeros ( n, o );
  w = zeros ( o, 1 );

  k = 0;
%
%  2^N points.
%
  k = k + 1;
  x(1:n,k) = - r;
  w(k) = a;
  more = 1;

  while ( more )
    more = 0;
    for i = n : -1 : 1
      if ( x(i,k) < 0.0 )
        k = k + 1;
        x(1:n,k) = x(1:n,k-1);
        x(i,k)     = + r;
        x(i+1:n,k) = - r;
        w(k) = a;
        more = 1;
        break;
      end
    end
  end

  return
end


text:
function [subject] = mrGetSubject(voldr)
%function [subject] = mrGetSubject(voldr)
%
% PURPOSE: Get the name of a subject on whom we have volume data.
% 07.19.97 - Created by Poirson
% 10.27.98 - Rewritten by Press to get rid of unix calls

subject = 'NotAPerson';

while ~exist([voldr '/' subject],'dir')

  disp('Current list of subject:');
  dir(voldr);
  subject = input('Subject name: ','s');
  disp('');

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [success] = example_downloadTestData(varargin)
	% Downloads CIAtah example test data from online into data folder.
	% Biafra Ahanonu
	% started: September 2018
	% inputs
		%
	% outputs
		%

	% changelog
		% 2019.09.16 [13:03:33] - Added three new imaging sessions to use for cross-day alignment and made downloading more generalized.
		% 2020.09.14 [13:17:57] - Added example two-photon dataset.
		% 2021.02.02 [11:25:34] - Function now calls data directory via standardized ciapkg.getDirPkg('data') to avoid placing data in incorrect folder.
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
		% 2022.11.07 [19:09:38] - Update links from Stanford to UCSF box (Stanford migrating away from Box, links are dead).
	% TODO
		% Switch Box downloads to Google Drive due to Stanford phase out.

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

	%========================
	options.downloadPreprocessed = 0;
	% Download extra files
	options.downloadExtraFiles = 1;
	% Binary: 1 = download meta data files.
	options.downloadMeta = 0;
	% Directory where download folder goes
	options.dataDir = ciapkg.getDirPkg('data');
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	try
		success = 0;
		downloadList = {};
			% downloadList{end+1}.folderName
			% downloadList{end}.fileName
			% downloadList{end}.fileUrl
			% downloadList{end}.metaFile

		downloadList{end+1}.folderName = '2014_04_01_p203_m19_check01';
		downloadList{end}.fileName = 'concat_recording_20140401_180333.h5';
			% downloadList{end}.fileUrl = 'https://stanford.box.com/shared/static/jmld9o9s0oemvn6oionr3lf9lwobqk9l.h5';
		downloadList{end}.fileUrl = 'https://ucsf.box.com/shared/static/dtdwpr4xjs1u67dv30qnhu7arf889ysr.h5';

		if options.downloadExtraFiles==1
			[downloadList] = localfxn_downloadExtra(downloadList);
		end

		nFiles = length(downloadList);

		for fileNo = 1:nFiles
			fileInfo = downloadList{fileNo};
			rawSavePathDownload = [options.dataDir filesep fileInfo.folderName];
			if ~exist(rawSavePathDownload,'dir');mkdir(rawSavePathDownload);fprintf('Made folder: %s\n',rawSavePathDownload);end

			rawSavePathDownload = [rawSavePathDownload filesep fileInfo.fileName];
			if exist(rawSavePathDownload,'file')~=2
				fprintf('Downloading file to %s\n',rawSavePathDownload)
				websave(rawSavePathDownload,fileInfo.fileUrl);
			else
				fprintf('Already downloaded %s\n',rawSavePathDownload)
			end

			if options.downloadMeta==1
				rawSavePathDownload = [rawSavePathDownload filesep fileInfo.metaFile];
				if exist(rawSavePathDownload,'file')~=2
					fprintf('Downloading file to %s\n',rawSavePathDownload)
					websave(rawSavePathDownload,fileInfo.metaFile);
				else
					fprintf('Already downloaded %s\n',rawSavePathDownload)
				end
			end
		end

		if options.downloadPreprocessed==1
			rawSavePathDownload = [options.dataDir filesep '2014_04_01_p203_m19_check01']
			if ~exist(rawSavePathDownload,'dir');mkdir(rawSavePathDownload);fprintf('Made folder: %s',rawSavePathDownload);end

			rawSavePathDownload = [rawSavePathDownload filesep '2014_04_01_p203_m19_check01_turboreg_crop_dfof_downsampleTime_1.h5'];
			if exist(rawSavePathDownload,'file')~=2
				fprintf('Downloading file to %s\n',rawSavePathDownload)
				websave(rawSavePathDownload,'https://stanford.box.com/shared/static/0zasceqd7b9ea6pa4rsgx1ag1mpjwmrf.h5');
			end

			rawSavePathDownload = [options.dataDir filesep '2014_04_01_p203_m19_check01' filesep '2014_04_01_p203_m19_check01_turboreg_crop_dfof_1.h5'];
			if exist(rawSavePathDownload,'file')~=2
				fprintf('Downloading file to %s\n',rawSavePathDownload)
				websave(rawSavePathDownload,'https://stanford.box.com/shared/static/azabf70oky7vriek48pb98jt2c5upj5i.h5');
			end
		end

		success = 1;
	catch err
		success = 0;
		disp('Check internet connection or download files manually, see *data* folder.')
		display(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		display(repmat('@',1,7))
	end
end

function [downloadList] = localfxn_downloadExtra(downloadList)
	% Download additional files

	downloadList{end+1}.folderName = ['batch' filesep '2014_08_05_p104_m19_PAV08'];
	downloadList{end}.fileName = 'concat_recording_20140805_162046.h5';
		% downloadList{end}.fileUrl = 'https://stanford.box.com/shared/static/iv3v3iqmji7fd5lvdcere8q6tlb008vb.h5';
	downloadList{end}.fileUrl = 'https://ucsf.box.com/shared/static/6uhqtff2qlgsrjqpt81iz5i7xbu5hfkv.h5';

	downloadList{end+1}.folderName = ['batch' filesep '2014_08_06_p104_m19_PAV09'];
	downloadList{end}.fileName = 'concat_recording_20140806_103546.h5';
		% downloadList{end}.fileUrl = 'https://stanford.box.com/shared/static/vyj69a7u2ay9uva9fbsyp34d8kvcvclw.h5';
	downloadList{end}.fileUrl = 'https://ucsf.box.com/shared/static/dh13zgf3ecsd5weolfuln4as09f6v4dg.h5';

	downloadList{end+1}.folderName = ['batch' filesep '2014_08_07_p104_m19_PAV10'];
	downloadList{end}.fileName = 'concat_recording_20140807_102507.h5';
		% downloadList{end}.fileUrl = 'https://stanford.box.com/shared/static/xmhgqx3atceq8f0zwqdctegiqo2d71su.h5';
	downloadList{end}.fileUrl = 'https://ucsf.box.com/shared/static/m2wvus76erqrng7bp13q7s3tlr727dgy.h5';

	downloadList{end+1}.folderName = ['twoPhoton' filesep '2017_04_16_p485_m487_runningWheel02'];
	downloadList{end}.fileName = 'concat_recording_20140807_102507.h5';
		% downloadList{end}.fileUrl = 'https://ucsf.box.com/shared/static/y71zhjejy7k90ov2o8wnbupi20o4q0nq.h5';
	downloadList{end}.fileUrl = 'https://ucsf.box.com/shared/static/y71zhjejy7k90ov2o8wnbupi20o4q0nq.h5';

	% downloadList{end+1}.folderName = ['batch' filesep '2014_07_31_p104_m80_PAV03'];
	% downloadList{end}.fileName = 'concat_recording_20140731_105559.h5';
	% downloadList{end}.fileUrl = 'https://stanford.box.com/shared/static/ddgy2zq3as9l3paifd1hltpojk93x9oa.h5';
	% % downloadList{end}.metaFile = 'https://stanford.box.com/shared/static/6qs3dt593uud314ins8qnxrnuj9sdm87.txt';

	% downloadList{end+1}.folderName = ['batch' filesep '2014_08_05_p104_m80_PAV04'];
	% downloadList{end}.fileName = 'concat_recording_20140805_180816.h5';
	% downloadList{end}.fileUrl = 'https://stanford.box.com/shared/static/psm56avve4gt3tquc88r7fri2xda9woe.h5';
	% % downloadList{end}.metaFile = 'https://stanford.box.com/shared/static/gf7agxqvp3ks0mr0ovy54t067agdp4tf.txt';

	% downloadList{end+1}.folderName = ['batch' filesep '2014_08_06_p104_m80_PAV05'];
	% downloadList{end}.fileName = 'concat_recording_20140806_104210.h5';
	% downloadList{end}.fileUrl = 'https://stanford.box.com/shared/static/ae8qcmkxcv8ax7g1yfs9qjp010gqyip3.h5';
	% % downloadList{end}.metaFile = 'https://stanford.box.com/shared/static/ja18hyade4jmh6czw0vxpxosiu3i9j9e.txt';
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function M = stiefelfactory(n, p, k)
% Returns a manifold structure to optimize over orthonormal matrices.
%
% function M = stiefelfactory(n, p)
% function M = stiefelfactory(n, p, k)
%
% The Stiefel manifold is the set of orthonormal nxp matrices. If k
% is larger than 1, this is the Cartesian product of the Stiefel manifold
% taken k times. The metric is such that the manifold is a Riemannian
% submanifold of R^nxp equipped with the usual trace inner product, that
% is, it is the usual metric.
%
% Points are represented as matrices X of size n x p x k (or n x p if k=1,
% which is the default) such that each n x p matrix is orthonormal,
% i.e., X'*X = eye(p) if k = 1, or X(:, :, i)' * X(:, :, i) = eye(p) for
% i = 1 : k if k > 1. Tangent vectors are represented as matrices the same
% size as points.
%
% By default, k = 1.
%
% See also: grassmannfactory rotationsfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%  July  5, 2013 (NB) : Added ehess2rhess.
%  Jan. 27, 2014 (BM) : Bug in ehess2rhess corrected.
%  June 24, 2014 (NB) : Added true exponential map and changed the randvec
%                       function so that it now returns a globally
%                       normalized vector, not a vector where each
%                       component is normalized (this only matters if k>1).

    if ~exist('k', 'var') || isempty(k)
        k = 1;
    end

    if k == 1
        M.name = @() sprintf('Stiefel manifold St(%d, %d)', n, p);
    elseif k > 1
        M.name = @() sprintf('Product Stiefel manifold St(%d, %d)^%d', n, p, k);
    else
        error('k must be an integer no less than 1.');
    end

    M.dim = @() k*(n*p - .5*p*(p+1));

    M.inner = @(x, d1, d2) d1(:).'*d2(:);

    M.norm = @(x, d) norm(d(:));

    M.dist = @(x, y) error('stiefel.dist not implemented yet.');

    M.typicaldist = @() sqrt(p*k);

    M.proj = @projection;
    function Up = projection(X, U)

        XtU = multiprod(multitransp(X), U);
        symXtU = multisym(XtU);
        Up = U - multiprod(X, symXtU);

% The code above is equivalent to, but much faster than, the code below.
%
%     Up = zeros(size(U));
%     function A = sym(A), A = .5*(A+A'); end
%     for i = 1 : k
%         Xi = X(:, :, i);
%         Ui = U(:, :, i);
%         Up(:, :, i) = Ui - Xi*sym(Xi'*Ui);
%     end

    end

    M.tangent = M.proj;

    % For Riemannian submanifolds, converting a Euclidean gradient into a
    % Riemannian gradient amounts to an orthogonal projection.
	M.egrad2rgrad = M.proj;

    M.ehess2rhess = @ehess2rhess;
    function rhess = ehess2rhess(X, egrad, ehess, H)
        XtG = multiprod(multitransp(X), egrad);
        symXtG = multisym(XtG);
        HsymXtG = multiprod(H, symXtG);
        rhess = projection(X, ehess - HsymXtG);
    end

    M.retr = @retraction;
    function Y = retraction(X, U, t)
        if nargin < 3
            t = 1.0;
        end
        Y = X + t*U;
        for i = 1 : k
            [Q, R] = qr(Y(:, :, i), 0);
            % The instruction with R assures we are not flipping signs
            % of some columns, which should never happen in modern Matlab
            % versions but may be an issue with older versions.
            Y(:, :, i) = Q * diag(sign(sign(diag(R))+.5));
        end
    end

    M.exp = @exponential;
    function Y = exponential(X, U, t)
        if nargin == 2
            t = 1;
        end
        tU = t*U;
        Y = zeros(size(X));
        for i = 1 : k
            % From a formula by Ross Lippert, Example 5.4.2 in AMS08.
            Xi = X(:, :, i);
            Ui = tU(:, :, i);
            Y(:, :, i) = [Xi Ui] * ...
                         expm([Xi'*Ui , -Ui'*Ui ; eye(p) , Xi'*Ui]) * ...
                         [ expm(-Xi'*Ui) ; zeros(p) ];
        end

    end

    M.hash = @(X) ['z' hashmd5(X(:))];

    M.rand = @random;
    function X = random()
        X = zeros(n, p, k);
        for i = 1 : k
            [Q, unused] = qr(randn(n, p), 0); %#ok<NASGU>
            X(:, :, i) = Q;
        end
    end

    M.randvec = @randomvec;
    function U = randomvec(X)
        U = projection(X, randn(n, p, k));
        U = U / norm(U(:));
    end

    M.lincomb = @lincomb;

    M.zerovec = @(x) zeros(n, p, k);

    M.transp = @(x1, x2, d) projection(x2, d);

    M.vec = @(x, u_mat) u_mat(:);
    M.mat = @(x, u_vec) reshape(u_vec, [n, p, k]);
    M.vecmatareisometries = @() true;

end

% Linear combination of tangent vectors
function d = lincomb(x, a1, d1, a2, d2) %#ok<INUSL>

    if nargin == 3
        d = a1*d1;
    elseif nargin == 5
        d = a1*d1 + a2*d2;
    else
        error('Bad use of stiefel.lincomb.');
    end

end


text:
function [theta,axis] = dquat2rot(dq)

% DQUAT2ROT  extracts the rotation axis and angle of a rotation dual
%            quaternion
%
%   [THETA,AXIS] = DQUAT2ROT(DQ) returns the rotation angle, THETA [deg],
%       and the rotation axis, AXIS, of a rotation dual quaternion DQ.
%       - DQ is a rotation dual quaternion. It is a 8-vector or an 8*N
%       array (each column is a rotation dual quaternion) where N is the
%       number of rotation dual quaternions.
%       - THETA is the rotation angle [deg]. It is comprised between 0 and
%       180deg. THETA is a scalar or an 1*N array (element i is the
%       rotation angle of dual quaternion i).
%       - AXIS is the unitary (norm equal to 1) rotation axis. It is a
%       3-vector or an 3*N array (column i is the rotation axis of
%       dual quaternion i).
%
% See also DQUAT2TRANS, DQUAT2SCREW, ROT2DQUAT

sdq = size(dq);
if sdq == [1 8], dq = dq'; sdq = size(dq); end
n = sdq(2);

% wrong format
if sdq(1) ~= 8
    error('DualQuaternion:dquat2rot:wrongsize',...
        '%d rows in the DQ array. It should be 8.',sdq(1));
end

% check that it is a rotation dual quaternion
tol = 1e-5;
[maxval,imax] = max(max(abs(dq(5:8,:))));
if maxval > tol
    warning('DualQuaternion:dquat2rot:wrongFormat',...
        'At least one dual quaternion is not a rotation dual quaternion (tol = %.1e).\n Indices of max values: %d \n Max value = %.2e',...
        tol,imax,maxval);
end
normDQ = sqrt(sum(dq(1:4,:).^2));
[maxval2,imax2] = max(abs(normDQ-1));
if maxval2 > tol
       warning('DualQuaternion:dquat2rot:wrongFormatForRotation',...
        'At least one dual quaternion is not a rotation dual quaternion (is it a rotation dual quaternion?) (tol = %.1e).\n Indices of max values: %d \n Max value = %.2e',...
        tol,imax2,maxval2);
end

indpos = find(dq(1,:) > 1);
dq(1,indpos) = ones(1,length(indpos));
indneg = find(dq(1,:) < -1);
dq(1,indneg) = ones(1,length(indneg));

% Extraction of the rotation parameters
theta = 2*acosd(dq(1,:)); % acosd: [-1,1] --> [0,180]
theta = mod(theta,360);
indthetaOK = find(theta > 0);
axis = [ones(1,n) ; zeros(2,n)];
if length(indthetaOK > 0)
    repsin = repmat(sind(theta(indthetaOK)/2),3,1);
    axis(:,indthetaOK) = dq(2:4,indthetaOK)./repsin;
end


text:
function rho_w = WeightedPearsonCorr(XI,XJ, w )
%WeightedPearsonCorr Pearson Correlation with Weights for use with
%outlier ensembles
%   Detailed explanation goes here

% Ensure that w is a column vector summing to 1
if size(w, 1) == 1
    w = w';
elseif size(w, 2) == 1
else
    error('w must be a vector');
end
w = w ./ sum(w);

m=size(XJ,1); % number of samples of p
p=size(XI,2); % dimension of samples

assert(p == size(XJ,2)); % equal dimensions
assert(size(XI,1) == 1); % pdist requires XI to be a single sample

mean_XI = dot(w, XI);
var_XI = dot(w, (XI - mean_XI).^2);

mean_XJ = sum(XJ * w, 2);
var_XJ = sum(((XJ - repmat(mean_XJ, 1, size(XJ, 2))).^2) * w, 2);

cov_XIJ = sum( (repmat((XI-mean_XI), size(XJ, 1), 1) .* (XJ - repmat(mean_XJ, 1, size(XJ, 2)))) * w, 2);

rho_w = cov_XIJ ./ sqrt(var_XJ) /  sqrt(var_XI);

% rho_w=zeros(m,1); % initialize output array
% for i=1:m
%     mean_XJ = dot(w, X
%     var_XJ
%     covariance = cov(XI, XJ(i,:)); covariance = covariance(1,2);
%     d(i,1) = 1 - covariance/stdXI/std(XJ(i,:));
% end

end

author: dsmi-lab-ntust

repo: AnomalyDetectionToolbox


text:
 function h = titlef(varargin)
%function h = titlef(varargin)
%|
%| version of title with built-in sprintf
%| and defaults to latex interpreter
%| also supports default font size from ir_fontsize()
%|
%| Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end

opt = {'fontsize', ir_fontsize('title')};
%opt = {opt{:}, 'fontname', 'times'};

if ir_is_octave
	opt = {opt{:}, 'fontname', 'Helvetica'};
end

varargin{1} = ir_strrep_tex(varargin{1}); % trick

if ir_is_octave
	tex = {};
	str = varargin{1};
	str = strrep(str, '$', ''); % remove $ needed for latex
	varargin{1} = str;
else
	tex = {'interpreter', 'latex'};
	str = varargin{1};
	str = strrep(str, '\', '\\'); % because of sprintf below
	% trick to replace "|" with "$|$" in strings
	if ~isempty(strfind(str, '|')) && isempty(strfind(str, '$'))
		str = strrep(str, '|', '$|$');
	end
	varargin{1} = str;
end

if isfreemat
	for ii=1:length(varargin)
		if streq(varargin{ii}, 'interpreter') % not supported by freemat
			varargin{ii} = {};
			varargin{ii+1} = {};
		end
	end
end

if im
	tmp = sprintf(varargin{:});
	tmp = strrep(tmp, '%', '\%'); % else latex will ignore %
	hh = title(tmp, tex{:}, opt{:});
else
	hh = [];
end

if nargout
	h = hh;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [x,info] = mgfracLapP1P2(A,b,elem,option,varargin)
%% MGFRACLAPP1P2 multigrid solvers for fractional Laplacian
%
%
% See also mg
%
% Documentation in Help browser ifem mgdoc
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

tic;

%% Size of systems
Ndof = length(b);               % number of dof
N = max(elem(:));               % number of nodes
Nydof = Ndof/N;                 % number of dof in the extend direction
Ny = (Nydof + 1)/2;             % number of nodes in the extend direction
% NT = size(elem,1);                 % number of elements
dim = size(elem,2)-1;

%% Options
% Assign default values to unspecified parameters
if ~exist('option','var'),
    option = [];
end
option = mgoptions(option,Ndof);    % parameters
x0 = option.x0;
N0 = option.N0/5;
tol = option.tol;
maxIt = option.solvermaxit;
mu = option.smoothingstep;
solver = option.solver;
% preconditioner = option.preconditioner;
preconditioner = 'none';
coarsegridsolver = option.coarsegridsolver;
printlevel = option.printlevel;
option.smoother = 'LINE';

%% Fixed dof
isFixDof = [];
if isfield(option,'freeDof') % freeDof is given
   isFreeDof = false(N,1);
   isFreeDof(option.freeDof) = true;
   NA = size(A,1);
   if NA > length(option.freeDof)
       isFixDof = true(NA,1);
       isFixDof(isFreeDof) = false;
   end
else % Find free dof and eliminate isolated dof
    deg = sum(spones(A));  % degree
    isFreeDof = false(Ndof,1);
    isFreeDof(deg>1) = true;
    isFixDof = ~isFreeDof;
    isFixDof(deg == 0) = false;
end
if any(isFixDof) % a bigger matrix is given
    xD = zeros(Ndof,1);
    xD(isFixDof) = b(isFixDof)./diag(A(isFixDof,isFixDof));
    A = A(isFreeDof,isFreeDof);
end
if length(b) > sum(isFreeDof)
    b = b(isFreeDof);
end
if length(x0) > sum(isFreeDof)
    x0 = x0(isFreeDof);
    option.x0 = x0;
end
isFreeNode = isFreeDof(1:N); % default choice for P1

%% Hierarchical Structure of Mesh
%  hierarchical structure for linear element.
if dim == 2  % 2-D
    [HB, NL, level] = HBstructure(elem,N0);
end
if dim == 3  % 3-D
    if nargin > 4
        HBmesh = varargin{end}; % need HB for bisection refinement
        [HB, NL, level] = HBstructure3(elem,N0,HBmesh);
    else % no HBmesh is given. only work for the red uniform refinement
        [HB, NL, level] = HBstructure3(elem,N0);
    end
end

%% Transfer operators between multilevel meshes for P1 element
% standard prolongation and restriction operator for P1 element
[Pro,Res] = transferoperator(HB,NL,isFreeNode);
for j = 1:level-1
    Pro{j} = kron(speye(Nydof-1,Nydof-1),Pro{j}); % for each level in y direction
    Res{j+1} = Pro{j}';
end
clear HB

%% Matrices in each level
Ai = cell(level,1);
Ai{level} = A;
for j = level:-1:2
    Ai{j-1} = Res{j}*Ai{j}*Pro{j-1};           % Ac = Res*Af*Pro
    switch option.smoother
        case 'GS'
            Bi{j} = tril(Ai{j});        % Forward Gauss-Seidel   B = D+L
            BBi{j} = triu(Ai{j});       % Backward Gauss-Seidel BB = D+U
        case 'LINE'
            NxNy = size(Ai{j},1);
            Nx = NxNy/(Nydof-1);
            di = zeros(NxNy,4);
            di(:,1) = full(diag(Ai{j}))/2; % half of the main diagoal
            temp = full(diag(Ai{j},Nx)); % next vertex in the extend direction
            di(1:length(temp),2) = temp;
            temp = full(diag(Ai{j},Nx*(Ny-2))); %
            di(1:length(temp),3) = temp;
            temp = full(diag(Ai{j},Nx*(Ny-1))); % middle point in the extend direction
            di(1:length(temp),4) = temp;
            halfBi = spdiags(di,[0 -Nx -Nx*(Ny-2) -Nx*(Ny-1)],NxNy,NxNy);
            Bi{j} = halfBi + halfBi';
            L{j} = chol(Bi{j},'lower');
            R{j} = L{j}';
%             Bi{j} = spdiags(diag(Ai{j}),0,size(Ai{j},1),size(Ai{j},1));        %
%             BBi{j} = Bi{j};       % Jacobi iteration
    end
end

%% MG cycles
% initial set up
k = 1;
x = x0;
r = b - A*x;
nb = norm(b);
err = zeros(maxIt,2);
if nb > eps  % nb is non-zero
    err(1,:) = norm(r)/nb;
else
    err(1,:) = norm(r);
end
switch solver
    case 'VCYCLE'
        if printlevel >= 1
            fprintf('Multigrid Vcycle Iteration \n')
        end
        while (max(err(k,:)) > tol) && (k <= maxIt)
            k = k + 1;
            % Step 2: Compute Br by one Vcylce MG
            Br = vcycle(r);
            % Step 3: Correct the solution
            x = x + Br;
            % Step 1: Form residual r
            r = r - A*Br;
            err(k,1) = sqrt(abs(Br'*r/(x'*b))); % approximate relative error in energy norm
            err(k,2) = norm(r)/nb; % relative error of the residual in L2-norm
            if printlevel >= 2
                fprintf('#dof: %8.0u,  #nnz: %8.0u, MG Vcycle iter: %2.0u, err = %8.4e\n',...
                         Ndof, nnz(A), k-1, max(err(k,:)));
            end
        end
        err = err(1:k,:);
        itStep = k-1;
    case 'WCYCLE'
        if printlevel >= 1
            fprintf('Multigrid Wcycle Iteration \n')
        end
        while (max(err(k,:)) > tol) && (k <= maxIt)
            k = k + 1;
            % Step 2: Compute Br by one Vcylce MG
            Br = wcycle(r);
            % Step 3: Correct the solution
            x = x + Br;
            err(k,1) = sqrt(abs(Br'*r/(x'*b))); % approximate relative error in energy norm
            % Step 1: Form residual r
            r = r - A*Br;
            err(k,2) = norm(r)/nb; % relative error of the residual in L2-norm
            if printlevel >= 2
                fprintf('#dof: %8.0u,  #nnz: %8.0u, MG Wcycle iter: %2.0u, err = %8.4e\n',...
                         Ndof, nnz(A), k-1, max(err(k,:)));
            end
        end
        err = err(1:k,:);
        itStep = k-1;
end

%% Krylov space method
% set up preconditioner
switch preconditioner
    case 'V'
        prefunc = @vcycle;
        if printlevel >= 1
            fprintf('Multigrid V-cycle Preconditioner with ')
        end
    case 'W'
        prefunc = @wcycle;
        if printlevel >= 1
            fprintf('Multigrid W-cycle Preconditioner with ')
        end
end
switch solver
    case 'CG'
        if printlevel >= 1
            fprintf('Conjugate Gradient Method\n')
        end
        [x,flag,err,itStep] = pcg(A,b,tol,maxIt,prefunc,[],x0);
    case 'MINRES'
        if printlevel >= 1
            fprintf('Minimum Residual Method \n')
        end
        [x,flag,err,itStep] = minres(A,b,tol,maxIt,prefunc,[],x0);
    case 'GMRES'
        if printlevel >= 1
            fprintf('General Minimum Residual Method\n')
        end
        if isfield(option,'restart')
            restart = option.restart;
        else
            restart = min(N,10);
        end
        [x,flag,err,itStep] = gmres(A,b,restart,tol,maxIt,prefunc,[],x0);
        itStep = (itStep(1)-1)*restart + itStep(2);
end

%% Modify x to include fix dof
if any(isFixDof)
    xD(isFreeDof) = x;
    x = xD;
end

%% Output
if k > maxIt
    flag = 1;
else
    flag = 0;
end
time = toc;
if printlevel >= 2
    fprintf('#dof: %8.0u, level: %2.0u,   coarse grid %2.0u, #nnz: %8.0u\n',...
              Ndof, level, size(Ai{1},1), nnz(Ai{1}))
end
if printlevel >= 1
    fprintf('#dof: %8.0u,  #nnz: %8.0u, smoothing: %2.0u, iter: %2.0u,   err = %8.4e,   time = %4.2g s\n',...
                 Ndof, nnz(A), mu, itStep, max(err(end,:)), time)
end
if (flag == 1) && (printlevel>0)
   fprintf('NOTE: the iterative method does not converge! \n');
end
info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag,'stopErr',max(err(end,:)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions vcycle, wcycle, fcycle, bpx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Vcycle MG
    function Br = vcycle(r,J)        % solve equations Ae = r in each level
    if nargin<=1
        J = level;
    end
    ri = cell(J,1);            % record residual in each level
    ei = cell(J,1);            % record err in each level
    ri{J} = r;
    for i = J:-1:2
        ei{i} = linesmoother(ri{i},i);   % pre-smoothing
        for s = 1:mu-1           % extra mu-1 steps smoothing
            ei{i} = ei{i} + linesmoother(ri{i}-Ai{i}*ei{i},i);
        end
        ri{i-1} = Res{i}*(ri{i} - Ai{i}*ei{i});
    end
    if strcmp(coarsegridsolver,'direct')
        ei{1} = Ai{1}\ri{1}; % direct solver in the coarest level
    else                         % iterative solver in the coarest level
        D = spdiags(diag(Ai{1}),0,size(Ai{1},1),size(Ai{1},1));
        [ei{1},flag] = pcg(Ai{1},ri{1},1/size(Ai{1},1),1000,D);
    end
    for i = 2:J
        ei{i} = ei{i} + Pro{i-1}*ei{i-1};
        ei{i} = ei{i} + linesmoother(ri{i}-Ai{i}*ei{i},i);
        for s = 1:mu-1
            ei{i} = ei{i} + linesmoother(ri{i}-Ai{i}*ei{i},i); % post-smoothing
        end
    end
    Br = ei{J};
    end

%% Wcycle MG
    function e = wcycle(r,J)        % solve equations Ae = r in each level
    if nargin<=1
        J = level;
    end
    if J == 1
        e = Ai{J}\r;   % exact solver in the coaresest grid
        return
    end
    % fine grid pre-smoothing
    e = Bi{J}\r;        % pre-smoothing
    for s = 1:mu-1        % extra mu-1 steps smoothing
        e = e + Bi{J}\(r-Ai{J}*e);
    end
    % coarse grid correction twice
    rc = Res{J}*(r - Ai{J}*e);
    ec = wcycle(rc,J-1);
    ec = ec + wcycle(rc - Ai{J-1}*ec,J-1);
    e = e + Pro{J-1}*ec;
    % fine grid post-smoothing
    e = e + BBi{J}\(r-Ai{J}*e);
    for s = 1:mu-1
        e = e + BBi{J}\(r-Ai{J}*e); % post-smoothing
    end
    end

%% Line smoother
    function e = linesmoother(r,J)
%         NxNy = length(r);
%         Nx = NxNy/(Ny-1);
%         e = zeros(NxNy,1);
%         for vi = 1:Nx
%             linedof = vi + 0:Nx:NxNy-1;
%             e(linedof) = Ai{J}(linedof,linedof)\r(linedof);
%         end
%         e = Bi{J}\r;
        e = R{J}\(L{J}\r);
        e = 0.75*e;
    end
end

repo: ifem


text:
function t = toeplitz( c, r )

%   Disciplined convex/geometric programming information for TOEPLITZ:
%      TOEPLITZ imposes no convexity restrictions on its arguments.
%      Instead of using the TOEPLITZ function, however, consider
%      creating a matrix variable using the 'toeplitz' keyword; e.g.
%          variable X(5,5) toeplitz;

%
% Check arguments
%

narginchk(1,2);
if nargin < 2,
    c    = vec( c );
    m    = length( c );
    p    = m;
    x    = [ cvx_subsref( c, p : -1 : 1 ) ; conj( cvx_subsref( c, 2 : p ) ) ];
else
    temp = cvx_subsref( r, 1 ) - cvx_subsref( c, 1 );
    if ~cvx_isconstant( temp ) || cvx_constant( temp ) ~= 0,
        warning('MATLAB:toeplitz:DiagonalConflict',['First element of ' ...
               'input column does not match first element of input row. ' ...
               '\n         Column wins diagonal conflict.'])
    end
    r = vec( r );
    c = vec( c );
    p = length( r );
    m = length( c );
    x = [ cvx_subsref( r, p : -1 : 2 ) ; c ];
end

%
% Construct matrix
%

cidx = ( 0 : m - 1 )';
ridx = p : -1 : 1;
t    = cidx( :, ones( p, 1 ) ) + ridx( ones( m, 1 ) , : );
t    = reshape( cvx_subsref( x, t ), size( t ) );

% 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:
% Log-harmonic scale
%
% Inputs
%  Hb    : Below this harmonic limit, the scale is linear
%          (similar to the mel scale which is linear below 1000Hz)
%          (e.g. 12)
%          Based on observation of the LF model, the asymptotic behavior of the
%          spectrum starts around the 12th harmonic (for the most tense voice)
%  Hmax  : The maxmimum number of harmonic considered during synthesis
%          (e.g. 256)
%  order : The reduced number of phase coefficients (e.g. 24)
%
% Outputs
%  hsl   : The log-harmonic scale
%
% Copyright (c) 2013 University of Crete - Computer Science Department(UOC-CSD)/
%                    Foundation for Research and Technology-Hellas - Institute
%                    of Computer Science (FORTH-ICS)
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%  Jonas Ballani for the bezier curve (see nested function)
%

function hsl = hlin2hlog(Hb, Hmax, order)

    %   The last compressed coef (i.e. order) corresponds to a given Hmax

    % Use linear scale below Hb (higher coefs will be overwritten)
    hsl = 1:Hmax;

    % Build a Bezier curve to start with a linear scale and finish smoothly
    % at (Hmax,order)
    p(1,:) = [Hb, Hb];
    p(2,:) = [order, order];
    p(3,:) = [Hmax, order];
    t = 0:0.01:1;
    [X,Y,p_bez] = CASTELJAU(0,1,p,t);
    hsl(Hb+1:end) = interp1(p_bez(:,1), p_bez(:,2), (Hb+1):Hmax);

    if 0
        hs = 1:Hmax;
        plot(hs, hsl, 'k');

        hold on;
        keyboard
    end

return

function [X,Y,val] = CASTELJAU(a,b,p,y)

    % function val = CASTELJAU(a,b,p,y)
    %
    % INPUT:  a   Linke Intervallgrenze
    %         b   Rechte Intervallgrenze
    %         p   Stützstellen (nx2-Matrix)
    %         y   Auswertungspunkte (Spaltenvektor)
    %
    % OUTPUT: val   Werte des Bezierpolynoms an y (mx2-Matrix)
    %
    % Date:   2007-11-05
    % Author: Jonas Ballani

    % Notes from degottex@csd.uoc.gr:
    %  From bezier.zip
    %  From http://m2matlabdb.ma.tum.de/download.jsp?MC_ID=7&SC_ID=8&MP_ID=480
    %  No license or copyright specified.

    n = size(p,1);
    m = length(y);
    T = zeros(n,n);
    val = zeros(m,2);
    X(:,1) = p(:,1);
    Y(:,1) = p(:,2);

    for j = 1:m
        for i = 2:n
            X(i:n,i) = (b-y(j))/(b-a)*X(i-1:n-1,i-1) + (y(j)-a)/(b-a)*X(i:n,i-1);
            Y(i:n,i) = (b-y(j))/(b-a)*Y(i-1:n-1,i-1) + (y(j)-a)/(b-a)*Y(i:n,i-1);
        end
        val(j,1) = X(n,n);
        val(j,2) = Y(n,n);
    end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
classdef dislocationSystem
% class representing dislocation
%
% Syntax
%   dS = dislocationSystem(b,l)
%   dS = dislocationSystem(sS)
%
% Input
%  b  - @Miller Burgers vector
%  l  - @Miller line vector
%  sS - @slipSystem
%  pr - poisson ratio
%
% Class Properties
%  b  - @Miller Burgers vector
%  l  - @Miller line vector
%  u  - energy
%  CS - @crystalSymmetry
%  isEdge  - is edge dislocation
%  isScrew - is screw dislocation
%
% See also
% DislocationSystems SlipSystems GND

  properties
    b % burgers vector
    l % line vector
    u % line energy of the dislocation system
  end

  properties (Dependent = true)
    CS
    isEdge
    isScrew
  end

  methods
    function dS = dislocationSystem(sS,l,u)

      if nargin == 0, return; end

      % adjust length burger vector edge dislocation a/2
      % for screw dislocations ??

      if isa(sS,'slipSystem')

        % define edge dislocations
        if sS.CS.lattice == 'hexagonal' %#ok<BDSCA>
          dS.b = 1/3 * sS.b;
        elseif sS.CS.lattice == 'cubic' %#ok<BDSCA>
          dS.b = 1/2 * sS.b;
        else
          dS.b = sS.b;
          warning('I could not determine the correct length of the Burgers vector. Please adjust it manually.')
        end
        dS.l = cross(sS.b,sS.n);

        % define screw dislocations
        b = 0.5*unique(sS.b,'antipodal','noSymmetry');
        dS.b = [dS.b(:);b(:)];
        dS.l = [dS.l(:);b(:)];

        % line energy
        dS.u = 1 + dS.isEdge;

      else

        b = sS;
        omega = angle(b,l,'noSymmetry');
        assert(all(omega > pi/2-1e-5 | omega<1e-5),...
          'line vector and burgers vector should be either be orthogonal! or identical')

        dS.b = sS;
        l.antipodal = true;
        dS.l = l;
        if nargin < 3, u = 1; end
        if numel(u) ~= length(dS.b), u = repmat(u,size(dS.b)); end
        dS.u = u;

      end

    end

    function CS = get.CS(sS)
      if isa(sS.b,'Miller')
        CS = sS.b.CS;
      else
        CS = specimenSymmetry;
      end
    end

    function isEdge = get.isEdge(dS)
      isEdge = angle(dS.b,dS.l,'noSymmetry') > pi/2 - 1e-5;
    end

    function isScrew = get.isScrew(dS)
      isScrew = angle(dS.b,dS.l,'noSymmetry') < 1e-5;
    end

    function display(dS,varargin)
      % standard output

      displayClass(dS,inputname(1),varargin{:});

      % display symmetry
      if isa(dS.CS,'crystalSymmetry')
        if ~isempty(dS.CS.mineral)
          disp([' mineral: ',char(dS.CS,'verbose')]);
        else
          disp([' symmetry: ',char(dS.CS,'verbose')]);
        end
      end

      toChar = @(x) char(round(x),'spaceSep');

      if any(dS.isEdge(:))

        disp([' edge dislocations : ',size2str(submatrix(ones(size(dS)),dS.isEdge))]);

        if isa(dS.CS,'crystalSymmetry')
          matrix = [arrayfun(toChar,dS.b(dS.isEdge),'UniformOutput',false),...
            arrayfun(toChar,dS.l(dS.isEdge),'UniformOutput',false),...
            vec2cell(dS.u(dS.isEdge))];

          cprintf(matrix,'-L',' ','-Lc',...
            {'Burgers vector' 'line vector' 'energy'},'-d','  ','-ic',true);
          disp(' ');

        end

      end

      if any(dS.isScrew(:))

        disp([' screw dislocations: ',size2str(submatrix(ones(size(dS)),dS.isScrew))]);

        if isa(dS.CS,'crystalSymmetry')
          matrix = [arrayfun(toChar,dS.l(dS.isScrew),'UniformOutput',false),...
            vec2cell(dS.u(dS.isScrew))];

          cprintf(matrix,'-L',' ','-Lc',...
            {'Burgers vector' 'energy'},'-d','  ','-ic',true);
          disp(' ');
        end
      end

    end

    function n = numArgumentsFromSubscript(varargin)
      n = 0;
    end

  end

  methods (Static = true)
    function dS = fcc(cs,varargin)
      dS = dislocationSystem(symmetrise(slipSystem.fcc(cs),'antipodal'));
    end

    function dS = bcc(cs,varargin)
      dS = dislocationSystem(symmetrise(slipSystem.bcc(cs),'antipodal'));
    end

    function dS = hcp(cs,varargin)
      dS = dislocationSystem(symmetrise(slipSystem.hcp(cs),'antipodal'));
    end

  end

end

author: mtex-toolbox

repo: mtex

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


text:
function value = year_length_months_julian ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_MONTHS_JULIAN returns the number of months in a Julian year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of months
%    in the year.
%
  value = 12;

  return
end


text:
%VAR Datafile overload
%
%   [V,U] = VAR(A,W)
%
% A - Datafile
% W - Vector of weights, one per object
%
% V - Vector of feature variances
% U - Mean vector
%
% Computes variance V and mean U in a single run for speed.
% Objects are assumed to have the same number of features.
% Take care that the feature size of A has been correctly set.
% The routine is useful in case the data is too large to be
% converted to a dataset first.

author: marianux

repo: ecg-kit

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


text:
function [ RGBpersp, XYZpersp ] = ht_Points2Persp( RGB, XYZ, P, H, W )
%RGB, XYZ: 3*(#Number of Points)
%P = K * [R, t]
%H, W: size of output perspective

%outputs
RGBpersp = nan(H, W, 3);
XYZpersp = nan(H, W, 3);

%XYZ->image coordinate
UV = P * [XYZ; ones(1, size(XYZ, 2))];
UV_norm = UV(3, :);
UV = bsxfun(@rdivide, UV(1:2, :), UV_norm);

% keyboard;

front_idx = UV_norm > 0;
UV = UV(:, front_idx);
RGB = RGB(:, front_idx);
XYZ = XYZ(:, front_idx);
UV_norm = UV_norm(front_idx);

visible_idx = UV(1, :) >= 0 & UV(2, :) >= 0 & UV(1, :) < W & UV(2, :) < H;
UV = UV(:, visible_idx);
RGB = RGB(:, visible_idx);
XYZ = XYZ(:, visible_idx);
UV_norm = UV_norm(visible_idx);

% UV= 1 + uint16(UV);
% UV = UV(:, UV_norm>0);
% UV_norm = UV_norm(UV_norm>0);
% RGB = RGB(:, UV_norm>0);
% XYZ = XYZ(:, UV_norm>0);
%
% UV = UV(:, UV(1, :) <=W & UV(2, :) <= H);
% UV_norm = UV_norm(UV(1, :) <=W & UV(2, :) <= H);
% RGB = RGB(:, UV(1, :) <=W & UV(2, :) <= H);
% XYZ = XYZ(:, UV(1, :) <=W & UV(2, :) <= H);

% keyboard;

%rendering
UV = 1 + uint16(floor(UV));
check_norm = ones(H, W) * max(UV_norm);
for pp = 1:1:size(UV, 2)
    if UV_norm(pp) <= check_norm(UV(2, pp), UV(1, pp))
        RGBpersp(UV(2, pp), UV(1, pp), 1) = RGB(1, pp);
        RGBpersp(UV(2, pp), UV(1, pp), 2) = RGB(2, pp);
        RGBpersp(UV(2, pp), UV(1, pp), 3) = RGB(3, pp);
        XYZpersp(UV(2, pp), UV(1, pp), 1) = XYZ(1, pp);
        XYZpersp(UV(2, pp), UV(1, pp), 2) = XYZ(2, pp);
        XYZpersp(UV(2, pp), UV(1, pp), 3) = XYZ(3, pp);
        check_norm(UV(2, pp), UV(1, pp)) = UV_norm(pp);
    end
end

RGBpersp = uint8(RGBpersp);

end

author: HajimeTaira


text:
function lagrange_interp_1d_test03 ( prob, m, nd )

%*****************************************************************************80
%
%% LAGRANGE_INTERP_1D_TEST03 tests LAGRANGE_APPROX_1D with Chebyshev spaced data.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROB, the problem index.
%
%    Input, integer M, the polynomial approximant degree.
%
%    Input, integer ND, the number of data points.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGRANGE_INTERP_1D_TEST03:\n' );
  fprintf ( 1, '  Approximate Chebyshev-spaced data from TEST_INTERP_1D problem #%d.\n', prob );
  fprintf ( 1, '  Use polynomial approximant of degree %d.\n', m );
  fprintf ( 1, '  Number of data points = %d\n', nd );

  a = 0.0;
  b = 1.0;
  xd = r8vec_chebyspace ( nd, a, b );

  yd = p00_f ( prob, nd, xd );

  if ( nd < 10 )
    r8vec2_print ( nd, xd, yd, '  Data array:' );
  end
%
%  #1:  Does interpolant match function at interpolation points?
%
  ni = nd;
  xi = xd;
  yi = lagrange_approx_1d ( m, nd, xd, yd, ni, xi );

  int_error = norm ( yi - yd ) / ni;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  L2 approximation error averaged per data node = %g\n', int_error );
%
%  #2: Plot the piecewise linear interpolant and polynomial approximants.
%
  figure ( 1 );
  clf
  ni = 501;
  xi = r8vec_linspace ( ni, a, b );
  yi = lagrange_approx_1d ( m, nd, xd, yd, ni, xi );

  hold on
  plot ( xd, yd, 'k.', 'Markersize', 30 );
  plot ( xd, yd, 'b-', 'Linewidth', 3 );
  plot ( xi, yi, 'r-', 'Linewidth', 4 );
  xlabel ( '<--- X --->' );
  ylabel ( '<--- Y --->' );
  title ( sprintf ( 'Lagrange/Chebyshev Polynomial Approximant of degree %d for %d nodes', m, nd ) )
  grid on
  hold off

  filename = sprintf ( 'p%02d_lagcheby_%02d_%02d.png', prob, m, nd );
  print ( '-dpng', filename );
  fprintf ( 1, '  Created plot file "%s".\n', filename );

  return
end

author: johannesgerer


text:
function verifyquadpts
%% VERIFYQUADPTS examples and verfication on quadrature rules.
%
%  error = verifyquadpts(n) computes the error of n-th order quadrature
%  rule in a triangle. This is an example on the usage of quadrature points
%  and verification of qudarture order for approximating integrals in a
%  triangle.
%
% See also quadpts
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

% a reference triangle
node = [0 0; 1 0; 0 1];
elem = [1 2 3];
area = 0.5;
err = zeros(9,2);
for n = 1:9
    % get quadrature points
    [lambda,weight] = quadpts(n);
    nQuad = size(lambda,1);
    t1 = 0;
    t2 = 0;
    for p = 1:nQuad
        % quadrature points in the x-y coordinate
        pxy = lambda(p,1)*node(elem(:,1),:) ...
            + lambda(p,2)*node(elem(:,2),:) ...
            + lambda(p,3)*node(elem(:,3),:);
        t1 = t1 + weight(p)*f1(pxy(1),pxy(2),n);
        t2 = t2 + weight(p)*f2(pxy(1),pxy(2));
    end
    t1 = t1*area;
    t2 = t2*area;
    err(n,1) = abs(t1 - 2/((n+1)*(n+2)));
    err(n,2) = abs(t2 - (sin(1) - cos(1)));
end
ts = zeros(9,3); ts = char(ts);
display('n     x^n + y^n     sin(x+y)');
display([num2str((1:9)') ts num2str(err(:,1),'%0.5e') ts num2str(err(:,2),'%0.5e')]);
end

function z = f1(x,y,n)
z = x.^n + y.^n;
end

function z = f2(x,y)
z = sin(x+y);
end

%% Results
% Let T be the triangle formed by (0,0), (1,0), and (0,1).
%
% Error1 is for the integral
% $\int _{T} x^n + y^n \, dxdy$.
% It should be numerically exact.
%
% Error2 is for the integral
% $\int _{T} \sin(x+y) \, dxdy$.
% It decays as n increas.
%
% See the doc for qudrature rules in <matlab:ifemdoc('quadpts') quadpts>.

repo: Aero-matlab


text:
function [out] = exchange_1(p1,p2,p3,S,fmax,dt)
%exchange_1 two-way channel exchange: linear and exponential.

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Water exchange between aquifer and channel
% Constraints:  f <= fIn
% @(Inputs):    p1   - linear scaling parameter [-]
%               p2   - linear scaling parameter [-]
%               p3   - exponential scaling parameter [-]
%               S    - current storage [mm]
%               fmax - maximum flux size [mm/d]
%               dt   - time step size [d]

out = max((p1*abs(S/dt) + p2*(1-exp(-1*p3*abs(S/dt)))).*sign(S),-1*fmax);

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT