text:
function OUT = xcorr_multisubject(data, varargin)
% Cross-correlation and partial correlation matrices for 3-D data, i.e., a cell array of subject data matrices
%
% :Usage:
% ::
%
%     OUT = xcorr_multisubject(data, [optional inputs])
%
% ..
%     Author and copyright information:
%     -------------------------------------------------------------------------
%     Copyright (C) 2009 Tor Wager
%
%     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/>.
% ..
%
% :Inputs:
%
%   **data:**
%        A cell array, one cell per subject/replicate, of n x k
%        data to be inter-correlated.
%
% :Optional Inputs:
%
%   **'partialr':**
%        Use partial correlations obtained via ridge regression (k = 1 fixed)
%
%   **'shift_by':**
%        Followed by integer value for max number of time points to shift
%
% :Output:
%
%   **OUT:**
%        A structure containing subject correlation matrices, the mean
%        matrix, and raw and FDR-thresholded group matrix
%
% :Examples:
% ::
%
%    % cl cell structure:
%    for i = 1:length(clpos_data), data{i} = cat(2, clpos_data{i}.timeseries); end
%
%    % cl structure:
%    for i = 1:size(cl(1).all_data, 2), for c = 1:length(cl), data{i}(:,c) = cl(c).all_data(:, i); end, end
%
%    % parcel_cl_avgs or clpos_data2 structure:
%    for i = 1:length(parcel_cl_avgs), for j = 1:N, data{j}(:,i) = parcel_cl_avgs(i).timeseries{j}; end, end
%
%    create_figure('Xcorr', 1, 3);
%    imagesc(OUT.stats.mean);
%    subplot(1, 3, 2);
%    imagesc(OUT.stats.sig);
%    subplot(1, 3, 3);
%    imagesc(OUT.stats.fdrsig);
%    colormap gray
%
% Example of MDS and plotting total (not decomposed) relationships:
% ::
%
%    OUT.stats.D = (1 - OUT.stats.mean) ./ 2;
%    [OUT.stats_mds.GroupSpace,OUT.stats_mds.obs,OUT.stats_mds.implied_dissim] = shepardplot(OUT.stats.D,[]);
%    OUT.stats_mds = nmdsfig_tools('cluster_solution',OUT.stats_mds, OUT.stats_mds.GroupSpace, 2:5, 1000, []);
%    nmdsfig(OUT.stats_mds.GroupSpace,'classes',OUT.stats_mds.ClusterSolution.classes,'names',OUT.stats_mds.names,'sig',OUT.stats.fdrsig);
%
% Example of MDS and plotting direct relationships:
% ::
%
%    OUT.ridge = matrix_direct_effects_ridge(data);
%    D = OUT.ridge.mean; D(find(eye(size(D)))) = 1;
%    D = (D' + D) ./ 2;
%    OUT.ridge.D = (1 - D) ./ 2;
%    [OUT.stats_mds.GroupSpace,OUT.stats_mds.obs,OUT.stats_mds.implied_dissim] = shepardplot(OUT.ridge.D,[]);
%    OUT.stats_mds = nmdsfig_tools('cluster_solution',OUT.stats_mds, OUT.stats_mds.GroupSpace, 2:10, 1000, []);
%    nmdsfig(OUT.stats_mds.GroupSpace,'classes',OUT.stats_mds.ClusterSolution.classes,'names',OUT.stats_mds.names,'sig',OUT.ridge.fdrsig);
%    hh = nmdsfig_fill(OUT.stats_mds);
%    axis image, axis equal
%
%    OUT.stats_mds = nmdsfig_tools('cluster_solution',OUT.stats_mds, OUT.stats_mds.GroupSpace, 2:5, 1000, []);
%
%    % data is cell, one cell per subject
%    [OUT.stats_mds.GroupSpace,OUT.stats_mds.obs,OUT.stats_mds.implied_dissim] = shepardplot(OUT.stats.D,[]);
%    OUT.stats_mds = nmdsfig_tools('cluster_solution',OUT.stats_mds, OUT.stats_mds.GroupSpace, 2:5, 1000, []);
%
% See also:
% ttest3d, plot_correlation_matrix -- for more compact ways of estimating
% and plotting correlation matrices. See also cellfun and canlab_mat2cell
% for conversion options.

% ..
%    Defaults
% ..
shift_by = 0;   % currently, if shift_by = 1, runs cross-correls and returns betas
robustflag = 0; % used if shift_by > 0, can do robust correlations
dopartialr = 0; % partial correlations
betaflag = 1;   % used if shift_by > 0, return betas instead of correlation
nconditions = 1; % fixed at 1 now

% optional inputs with default values
% -----------------------------------

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case 'partialr', dopartialr = 1;
            case 'shift_by', shift_by = varargin{i+1}; varargin{i+1} = [];

                %case 'basistype', basistype = varargin{i+1}; varargin{i+1} = [];

                %otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

fprintf(1,'Multi-subject cross correlation\n')

if shift_by
    fprintf(1,'Estimating time-shifted (cross-lagged) correlations\n')

elseif dopartialr
    fprintf(1,'Estimating partial correlations\n')
end

fprintf('\t   ')

warning off     % for robustfit iteration limit

numsub = length(data);

for i = 1:numsub            %size(DATA.dat,3);

    subjdat = data{i};
    m = size(subjdat, 2);

    fprintf(1,'\b\b\b%03d', i);

    warning off

    for n=1:nconditions

        if shift_by > 0

            if betaflag
                betastr = 'max_cross_lagged_beta'; % or corr
            else
                betastr = 'max_cross_lagged_correlation';
            end

            for j = 1:m - 1
                for k = j + 1 : m
                    % now can handle shift by 0 in shift_correl

                    [sval,myc, mylat, myxc] = shift_correl(subjdat(:,j), subjdat(:,k), shift_by, robustflag, betaflag);

                    if isempty(mylat), mylat = NaN; end

                    nxl(j,k) = mylat;
                    nxc(j,k) = myxc;
                end
            end
        elseif dopartialr
            betastr = 'partial_correlation';
            nxc = calc_partial_r(subjdat);
            nxl = [];

        else % %if there's no shift by; fastest
            betastr = 'correlation';
            nxc = corr(subjdat);
            nxl = [];
        end

        warning on

        if shift_by > 0
            % adjust matrices
            nxc(end+1,:) = 0;
            nxc=nxc + nxc' + eye(size(nxc,1));

            nxl(end+1,:) = 0;
            nxl=nxl+nxl';

        end

        OUT.metric_returned = betastr;
        OUT.shift_by = shift_by;
        OUT.shift_explanation = '0 for no latency est., n for cross-correlations shifting up to n time points forward/back';

        OUT.pairwise_assoc{n}(:,:,i) = nxc;
        clear nxc;

        OUT.latency{n}(:,:,i)  = nxl;
        clear nxl;

    end         % condition loop

end             % subject loop

fprintf('\n');

fprintf('stats...')
[mxc,t,sig,OUT.stats] = ttest3d(OUT.pairwise_assoc{1});

t = OUT.stats.t .* abs(OUT.stats.fdrsig);
t = (t + t') ./ 2;
OUT.stats.fdr_thresholded_tvalues = t;

fprintf('\n');

end % Main function

function [b p] = calc_partial_r(X)

X = zscore(X);

p = [];

for i = 1:size(X, 2)
    y = X(:, i);
    xx = X;
    xx(:, i) = 1;  % intercept; need it, and also placeholder

    % ols version
    b(:, i) = xx \ y;
    b(i, i) = 1;        % use 1 for self, not intercept

    % ridge version
    % note: we have double intercept...
    b(:, i) = ridge(y, xx, 1);
    b(i, i) = 1;        % use 1 for self, not intercept

    %create_figure; plot(b(:, i)); hold on; plot(b1, 'r'); drawnow; pause(.05)

end

end % function


text:
function cS = times(alpha,cS)

if ~isa(cS,'crystalShape')
  [cS,alpha] = deal(alpha,cS);
end

alpha = repmat(alpha(:).',size(cS.V,1),size(cS.V,2)/numel(alpha));
cS.V = alpha .* cS.V;

author: mtex-toolbox

repo: mtex

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


text:
%% Qudarature
%
%%
%

fun = @(v) v.x .* v.y;

sF = S2FunHarmonic.quadrature(fun)

surf(sF)


text:
function test_bug3336

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY

% there are two versions: one with error() and one with ft_error()
% these should give the same result

% compare the two ways of giving an error

try
  subfunction1err
catch
  err1 = lasterror;
end

try
  subfunction2err
catch
  err2 = lasterror;
end

assert(~isequal(err1.message, err2.message))
assert(isequal(err1.identifier, err2.identifier))

% idem for warnings
% but here I cannot get access to the stack

subfunction1warn
[msg1, id1] = lastwarn;

subfunction2warn
[msg2, id2] = lastwarn;

assert(~isequal(msg1, msg2))
assert(isequal(id1, id2))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% VERSION 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function subfunction1err
subsubfunction1err

function subsubfunction1err
error('SOME:ERROR', 'something went wrong 1');

function subfunction1warn
subsubfunction1warn

function subsubfunction1warn
warning('SOME:WARNING', 'something went wrong 1');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% VERSION 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function subfunction2err
subsubfunction2err

function subsubfunction2err
ft_error('SOME:ERROR', 'something went wrong 2');

function subfunction2warn
subsubfunction2warn

function subsubfunction2warn
warning('SOME:WARNING', 'something went wrong 2');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function msh = mrmInflate(meshFile,iter,inflatedMeshFile)
% Inflate (unfold) mesh while minimizing distance and area distortions.
%
%   msh = mrmInflate([meshFile],[iterations=400],[inflatedMeshFile]);
%
% Must be run on a Linux machine, as this mainly uses Jonas Larson's
% surfrelax (Linux executable)
%

if ieNotDefined('meshFile'),
  [filename, pathname] = uigetfile('*.mat','Select mrMesh file');
  meshFile             = fullfile(pathname,filename);
  drawnow;
end;
if ieNotDefined('iter'),
    iter = 400;
end;
if ieNotDefined('inflatedMeshFile')
    if(nargout<1)
        [p,f,e]          = fileparts(meshFile);
        inflatedMeshFile = fullfile(p, [f '_surfRelaxed_' num2str(iter) e]);
    else
        inflatedMeshFile = '';
    end
end;

surfrelaxFile = fullfile(fileparts(which('mrMesh')), 'surfrelax');

if(ischar(meshFile))
    load(meshFile);
elseif(isstruct(meshFile))
    msh = meshFile;
    clear meshFile;
else
    error('First argument must be either a filename or a mesh structure.');
end

% if data field exists this is probably an older mesh
if(isfield(msh,'data'))
    triangles = msh.data.triangles;
    origin    = msh.data.origin;
    vertices  = msh.data.vertices;
else
    triangles = msh.triangles;
    origin    = msh.origin;
    vertices  = msh.vertices;
end
if(isfield(msh,'initVertices')&&~isempty(msh.initVertices))
    vertices = msh.initVertices;
end

% Save mesh in OOGL (.off) format
inFile = [tempdir 'tempOOGL.off'];
outFile = [tempdir 'tempOOGL_relaxed.off'];

%
% WRITE .off FILE
%
% Open file
fid = fopen(inFile, 'w');
if fid == -1,
  delete(inFile);
  error('Could not open file');
end;

% Header
fprintf(fid, 'OFF\n');
fprintf(fid, '%d %d %d\n', ...
        size(vertices,2), size(triangles,2), 0);

% Vertices
for ii=1:3
  vertices(ii,:) = vertices(ii,:) + origin(ii);
end
fprintf(fid, '%f %f %f\n', vertices);

% Faces
fprintf(fid, '3 %d %d %d\n', triangles);

% Close file
fclose(fid);

%--- actual smoothing
opt = ['-iter ' num2str(iter) ' -dist 0.1'];
system([surfrelaxFile ' ' opt ' ' inFile ' ' outFile]);

% Read new vertex locations
fid = fopen(outFile, 'r', 'ieee-be');
if fid == -1,
  delete(inFile);
  delete(outFile);
  error('Could not open file');
end

%--- read .off output file
% Header
hdr = fgetl(fid); %#ok<NASGU>
n = fread(fid, 3, 'int32');

% Read the vertices
vertices = fread(fid, [3 n(1)], 'float32');
newTriangles = fread(fid, [5 n(2)], 'int32');
fclose(fid);

for ii=1:3
  vertices(ii,:) = vertices(ii,:) - origin(ii);
end

newTriangles = newTriangles(2:4,:);
if(sum(newTriangles(:)-triangles(:))>0)
    error('triangles do not match!');
end

%--- put back
if isfield(msh,'data'),
  % mrMesh will recompute the normals if we remove them
  if(isfield(msh.data,'normals'))
    msh.data =  rmfield(msh.data,'normals');
  end;
  msh.data.vertices = vertices;
else
  if(isfield(msh,'normals'))
    msh =  rmfield(msh,'normals');
  end;
  msh.vertices = vertices;
end;
msh.lights = [];
msh.id = -1;

if(~isempty(inflatedMeshFile))
    % save
    [p, msh.name] = fileparts(inflatedMeshFile);
    msh.fileName = inflatedMeshFile;
    save(inflatedMeshFile, 'msh');
end
% clean up
delete(inFile);
delete(outFile);
if(nargout==0)
    clear msh;
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [coeffs, poles, k] = residue(u, v, k)
%RESIDUE   Partial-fraction expansion (residues).
%   [R, P, K] = RESIDUE(B, A) finds the residues, poles and direct term of a
%   partial fraction expansion of the ratio of two CHEBFUN objects B(s)/A(s).
%   If there are no multiple roots,
%        B(s)       R(1)       R(2)             R(n)
%        ----  =  -------- + -------- + ... + -------- + K(s)
%        A(s)     s - P(1)   s - P(2)         s - P(n)
%   B and A are CHEBFUN objects consisting of a single fun. The residues are
%   returned in the column vector R, the pole locations in column vector P, and
%   the direct terms in the CHEBFUN K. The number of poles is n = length(A) - 1
%   = length(R) = length(P). The direct term CHEBFUN is zero if length(B) <
%   length(A), otherwise length(K) = length(B) - length(A) + 1.
%
%   If P(j) = ... = P(j+m-1) is a pole of multiplicity m, then the expansion
%   includes terms of the form
%                  R(j)        R(j+1)                R(j+m-1)
%                -------- + ------------   + ... + ------------
%                s - P(j)   (s - P(j))^2           (s - P(j))^m
%
%   [B, A] = RESIDUE(R, P, K), with 3 input arguments and 2 output arguments,
%   converts the partial fraction expansion back to the polynomials with CHEBFUN
%   representation in B and A.
%
%   Warning: Numerically, the partial fraction expansion of a ratio of
%   polynomials represents an ill-posed problem. If the denominator polynomial,
%   A(s), is near a polynomial with multiple roots, then small changes in the
%   data, including roundoff errors, can make arbitrarily large changes in the
%   resulting poles and residues. Problem formulations making use of state-space
%   or zero-pole representations are preferable.

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

if ( min(size(u)) > 1 || min(size(v)) > 1 )
    error('CHEBFUN:CHEBFUN:residue:quasi', ...
          'Residue does not support CHEBFUN objects with multiple columns.');
end

if ( nargin == 2 )

    if ( (numel(u.funs) > 1) || (numel(v.funs) > 1) )
        error('CHEBFUN:CHEBFUN:residue:multipleFuns', ...
              'Residue does not support CHEBFUNs consisting of multiple FUNs.');
    end
    if ( ~isfinite(u) || ~isfinite(v) )
        error('CHEBFUN:CHEBFUN:residue:inf', ...
              'Residue does not support functions with nonzero exponents.');
    end

    b = poly(u);
    a = poly(v);
    [coeffs, poles, k] = residue(b, a);
    k = chebfun(@(x) polyval(k, x), max(length(k), 1), 'vectorize');

elseif ( nargin == 3 )

    % Interpret an empty CHEBFUN as a zero CHEBFUN for the user's convenience.
    if ( isempty(k) )
        k = 0;
    else
        if ( numel(k.funs) > 1 )
            error('CHEBFUN:CHEBFUN:residue:multipleFuns', ...
                  'RESIDUE does not support CHEBFUNs with multiple FUNS.');
        end

        if ( ~isfinite(k) )
            error('CHEBFUN:CHEBFUN:residue:inf', ...
                  'RESIDUE does not support functions which are unbounded.');
        end

        k = poly(k);
    end

    [b, a] = residue(u, v, k);
    coeffs = chebfun(@(x) polyval(b, x), max(length(b), 1), 'vectorize');
    poles = chebfun(@(x) polyval(a, x), max(length(a), 1), 'vectorize');

else

    error('CHEBFUN:CHEBFUN:residue:args', ...
        'Residue requires either 2 or 3 arguments.');

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ax,ci,cen,ub,lb,S,e,lam,Fm,Fc,F,pval, msb] = conf_region(X,varargin)
% :Usage:
% ::
%
%    [ax,ci,cen,ub,lb,S,e,lam,Fm,Fc,F,pval, msb] = conf_region(X,[doplot])
%
% alternative conf region multivariate based on Johnson & Wichern, 4th ed., p. 236
% 2D (3D)
%
% :Outputs:
%
%   **ax:**
%        axes of confidence region, scaled to be half-length of confidence hyperellipsoid
%
%   **ci:**
%        length of axes of conf region (diagonal matrix), axes in reverse order of importance
%        as in ouput of eig
%
%   **cen:**
%        center of region (means of variables)
%
%   **ub:**
%        upper boundary coordinates for region, rows are dims (vars), cols index coordinates
%        last coordinate is on axis of greatest variation ("reverse order"), from output of eig
%
%   **lb:**
%        lower boundary coordinates
%        to plot, try: plot(ub(1,:),ub(2,:),'bo'); hold on; plot(lb(1,:),lb(2,:),'ro')
%
%   **S:**
%        covariance matrix
%
%   **e:**
%        eigenvectors
%
%   **lam:**
%        eigenvalues
%
%   **Fm:**
%        degrees of freedom multiplier, p(n-1) / n(n-p)
%
%   **Fc:**
%        critical F value based on alpha level
%
%   **F:**
%        F value for test - probably not quite right
%
%   **pval:**
%        p value for test - probably not quite right
%
% :To plot:
% Either enter 1 or color as a 2nd argument, or do it yourself:
% ::
%
%    [ax,ci,cen,ub,lb,S,e,lam] = conf_region(X);
%    b = pinv([X(:,1) ones(size(X,1),1)] ) * X(:,2);
%    theta = atan(b(1)) - pi/2;
%    [h,h2] = plot_ellipse(cen(1),cen(2),theta,ci(1,1),ci(2,2));
%
% There is also an example for rendering individual subject conf regions
%
% :Examples:
% ::
%
%    N = 250;
%    X = mvnrnd([1 2], [1 .6; .6 1], N);
%    [ax,ci,cen,ub,lb,S,e,lam] = conf_region(X);
%    b = pinv([X(:,1) ones(size(X,1),1)] ) * X(:,2);
%    theta = atan(b(1)) - pi/2;
%    %create_figure('test');
%    [h,h2] = plot_ellipse(cen(1),cen(2),theta,ci(1,1),ci(2,2));
%
%    % plot standard deviation - for boostrapping, or for individual cases
%    [h,h2] = plot_ellipse(cen(1),cen(2),theta,ci(1,1)*sqrt(N),ci(2,2)*sqrt(N));
%
% :Example: Render individual subject conf regions
% ::
%
%    X = x(wh, [3 1]);
%    dfe = size(X, 1) - size(X, 2);   % df
%    alph = .1;  % .1 for 90%, .05 for 95%
%    tcrit = tinv(1 - alph, dfe);
%    ci = ((lam) .^ .5) .* tcrit;
%    b = pinv([X(:,1) ones(size(X,1),1)] ) * X(:,2);
%    theta = atan(b(1)) - pi/2;
%    [h,h2] = plot_ellipse(cen(1),cen(2),theta,ci(1,1),ci(2,2));
%
%    set(h, 'Color', 'k', 'LineWidth', 1);
%    set(h2, 'FaceColor', [.5 1 0]);

doplot = 0;
if length(varargin) > 0, doplot = varargin{1}; end

% -----------------------------------------------------
% * determine multivariate standard deviation matrix S
% -----------------------------------------------------

% center
Xs = X - repmat(mean(X),size(X,1),1);

% covariance matrix S
S = (Xs' * Xs) ./ (size(Xs,1)-1);

% -----------------------------------------------------
% * get squared distance
% -----------------------------------------------------

% squared distance, Johnson & Wichern p. 201
% (X-mean(X))' * inv(S) * (X - mean(X)) generalized to matrices
d = Xs * inv(S) * Xs';
d = diag(d);            % what do the off-diagonals signify?

% -----------------------------------------------------
% * get critical F and coefficient of mult. for axes
% -----------------------------------------------------
p = size(X,2);
n = size(X,1);
Fm = (p * (n - 1)) ./ (n * (n - p));
a = .95;                                % alpha
Fc = finv(a,p,n - p);                   % critical F
coef = sqrt(Fm .* Fc);

[e,lam] = eig(S);

% -----------------------------------------------------
% * get  F and p-values
% ratio of squared distances (mahal?) between to within
% -----------------------------------------------------
msb = pdist([mean(X); zeros(1,size(X,2))]); % .^ 2;
msw = sum(d) ./ (n-p);
F = msb / msw;
pval = 1 - fcdf(F,p,n-p); % - this is wrong I think. do nonparam.

% nonparametric
niter = 1000;dnull = [];
signv = sign(randn(size(X,1),niter));
for i = 1:niter
    xtmp = X .* repmat(signv(:,i),1,size(X,2));
    dnull(i) = pdist([mean(xtmp); zeros(1,size(X,2))]);
end
pval = sum(msb <= dnull) ./ length(dnull);

% -----------------------------------------------------
% * confidence interval, in units - mult by axes.
% -----------------------------------------------------
ci = ((lam) .^ .5 .* coef);
ax = e * ci;

cen = mean(X)';
lb = repmat(cen,1,size(e,2)) - (e * ci);
ub = repmat(cen,1,size(e,2)) + (e * ci);

if doplot
    % 2-d plot

    % degrees or slope to angle in radians
    %ang = corrcoef(X); ang = ang(1,2);
    %theta = acos(ang-(pi/2));
    %pi * ang / 180

    b = pinv([X(:,1) ones(size(X,1),1)] ) * X(:,2);
    theta = atan(b(1)) - pi/2;
    [h,h2] = plot_ellipse(cen(1),cen(2),theta,ci(1,1),ci(2,2));
    set(h,'Color',[1 1 1],'LineWidth',.5);

    if length(doplot) == 3 , set(h2,'FaceColor',doplot); end
    if ischar(doplot), set(h2,'FaceColor',doplot(1)); end

    drawnow
end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function p00_jac_test ( problem_num )

%*****************************************************************************80
%
%% P00_JAC_TEST compares the jacobian to a finite difference estimate.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM_NUM, the number of problems.
%
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'P00_JAC_TEST\n' );
  fprintf ( 1, '  Find the maximum relative difference between the\n' );
  fprintf ( 1, '  jacobian and a finite difference estimate.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   Problem    Option      Diff               I         J\n' );

  for problem = 1 : problem_num

    option_num = p00_option_num ( problem );

    fprintf ( 1, '\n' );

    for option = 1 : option_num

      nvar = p00_nvar ( problem, option );

      x0 = p00_start ( problem, option, nvar );

      for i = 1 : nvar
        [ dx, seed ] = r8_uniform_01 ( seed );
        x0(i) = x0(i) + 0.01 * dx;
      end

      [ max_adif, max_adif_i, max_adif_j, max_rdif, max_rdif_i, max_rdif_j ] = ...
        p00_jac_check ( problem, option, nvar, x0 );

      fprintf ( 1, '  %8d  %8d  %14e  %8d  %8d\n', ...
        problem, option, max_rdif, max_rdif_i, max_rdif_j );

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef TestEigen
    %TestEigen

    methods (Static)
        function test_1
            A = randn(5); A = A.'*A;
            evals = cv.eigen(A);
            [evals,evecs,b] = cv.eigen(A);
        end

        function test_compare_against_eig
            A = gallery('lehmer',4);
            [evals,evecs] = cv.eigen(A);

            % match orientation, order, and sign
            [V,D] = eig(A);
            [D,ord] = sort(diag(D), 'descend');
            V = V(:,ord).';
            idx = (sign(evecs(:,1)) ~= sign(V(:,1)));
            V(idx,:) = -1 * V(idx,:);

            % compare
            assert(norm(D - evals) < 1e-6);
            assert(norm(V - evecs) < 1e-6);
        end

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

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function variance = i4row_variance ( m, n, a )

%*****************************************************************************80
%
%% I4ROW_VARIANCE returns the variances of an I4ROW.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of data in TABLE.
%
%    Input, integer A(M,N), the array.
%
%    Output, real VARIANCE(M), the variance of each row of TABLE.
%
  for i = 1 : m

    mean = sum ( a(i,1:n) ) / n;

    variance(i) = 0.0;
    for j = 1 : n
      variance(i) = variance(i) + ( a(i,j) - mean )^2;
    end

    if ( 1 < n )
      variance(i) = variance(i) / ( n - 1 );
    else
      variance(i) = 0.0
    end

  end

  return
end


text:
function [CC, rhs, bb, gg, Px, Py, xsplit, ysplit] = discretize(N, f, m, n, flag)
%DISCRETIZE   Given a CHEBOP2, this function converts the problem to one of the form
%
%  sum_i  kron(A_i,B_i)
%
% and computes the discretisation as a cell array:
%
%  {{
%
%       A_1 , B_1
%       A_2 , B_2
%        .  ,  .
%        .  ,  .
%        .  ,  .
%       A_k , B_k
%
%                   }}
%
%    INPUTS:
%      N = PDE (CHEBOP2).
%      f = forcing term (CHEBFUN2).
%      m = discretization size in 1st variable.
%      n = discretization size in 2nd variable.
%      flag = 0 (default) means assigned boundary conditions, flag = 1
%      means do not assign boundary conditions.
%
%    OUTPUTS:
%      CC = cell array of matrices storing the terms in the matrix
%      equation.
%      rhs = matrix discretizing forcing term.
%      bb = cell array storing the discretized linear constraints.
%      gg = cell array storing the discretized nonhomogeneous part of the
%      constraints.
%      Px, Py = store permutation matrix to ensure bcs are linear
%      dependent.
%      XSPLIT, YSPLIT = 0 if subproblems cannot be formed. XSPLIT = 1 if
%      even and odd modes decouple in 1st variable. YSPLIT = 1 if even and
%      odd modes decouple in 2nd variable.
%
% Returns RHS with degrees of freedom removed and bb which stores
% elminated boundary conditions, gg eliminated boundary rows.

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

% Parse inputs.
if nargin < 5
    flag = 0;
end

%%
% Get information of PDE and pref.
A = N.coeffs;
rect = N.domain;
pref = chebfunpref();
tol = pref.cheb2Prefs.chebfun2eps;
xorder = N.xorder;
yorder = N.yorder;

%%
% Check if the PDO was given as a variable coefficient PDO with the notation
% @(x,y,u), but is actually a constant coefficient PDO.
if ( iscell( A ) )
    doesNotDependOnXorY = all(all(cellfun(@isnumeric, N.coeffs)));
    if ( doesNotDependOnXorY )
        A = cell2mat( A ).';
    end
end

%%
% Convert matrix of coefficients to a discretization for the PDE using the
% singular value decomposition.  We find the rank of the PDE operator and
% then use the optimal low rank expansion of the operator as a way to
% discretise the PDE.
if (  isempty(N.V) || isempty(N.U) )
    if ( iscell(A) )
        counter = 1;
        U = cell(size(A, 1), 1);
        V = cell(size(A, 2), 1);
        for jj = 1:size(A, 1)
            for kk = 1:size(A, 2)
                a = A{jj,kk};
                if ( isa(a, 'chebfun2') )
                    if ( abs(vscale(a)) > tol )
                        [C, D, R] = cdr(a);
                        for col = 1:size(C, 2)
                            U{jj,counter} = C(:,col)*D(col,col);
                            V{kk,counter} = R(:,col);
                            counter = counter + 1;
                        end
                    end
                elseif ( isa(a,'double') && ( isempty(a) || abs(a) > tol ) )
                    U{jj,counter} = a;
                    V{kk,counter} = 1;
                    counter = counter + 1;
                end
            end
        end
        rk = size(U, 2);
        S = diag(ones(rk, 1));
        na = size(U, 1);
        nb = size(V, 1);
    else
        % Compute the SVD of the coefficient matrix.
        [U, S, V] = svd(A.');
        % Find the rank of A, which is also the rank of the PDE operator and
        % construct the low rank expansion for A.
        rk = find(diag(S) > tol, 1, 'last');
        U = U(:,1:rk);
        S = S(1:rk,1:rk);
        V = V(:,1:rk);
        [na, nb] = size(A.');
    end
else
    rk = size(N.S, 2);
    U = N.U;
    V = N.V;
    S = N.S;
    na = size(U, 1);
    nb = size(V, 1);
end
% LEFT = zeros(m); RIGHT = zeros(n);

% Construct the discretisation in matrix equation form.
CC = cell(rk,2);

for jj = 1 : rk

    RIGHT = unconstrainedMatrixEquation(V, jj, n, xorder, rect(1:2)); % jjth term on the right.
    LEFT = unconstrainedMatrixEquation(U, jj, m, yorder, rect(3:4)); % jjth term on the left.

    % Balance out the scaling from the singular value. This does slightly
    % improve the accuracy.
    singvalue = sqrt(S(jj,jj));
    CC{jj,2} = singvalue * RIGHT;
    CC{jj,1} = singvalue * LEFT;

end

%%
% Test to see if we can solve subproblems. This checks if the PDE operator
% contains differential terms of the same parity.
ysplit = 0; xsplit=0;
if ~iscell(U) && ~iscell(V)
    emask = 1:2:na;
    omask = 2:2:na;
    if ( min( norm(U(emask,:)), norm(U(omask,:)) ) < 10*tol )
        ysplit = 1;
    end
    emask = 1:2:nb;
    omask = 2:2:nb;
    if ( min( norm(V(emask,:)), norm(V(omask,:)) ) < 10*tol )
        xsplit = 1;
    end
end

%%
% We have a discretisation for the PDE operator, now let's find a
% discretisation for the boundary conditions.

% If no boundary conditions is prescribed then make it empty.
bcLeft = [];
leftVal = [];
bcRight = [];
rightVal = [];
bcUp = [];
upVal = [];
bcDown = [];
downVal = [];

if ( ~isempty(N.lbc) ) % Left boundary conditions.
    [bcLeft, leftVal] = chebop2.constructBC(N.lbc, -1, m, n, rect(3:4), rect(1:2), xorder);
end
if ( ~isempty(N.rbc) ) % Right boundary conditions.
    [bcRight, rightVal] = chebop2.constructBC(N.rbc, 1, m, n, rect(3:4), rect(1:2), xorder);
end
if ( ~isempty(N.ubc) ) % Top boundary conditions.
    [bcUp, upVal] = chebop2.constructBC(N.ubc, 1, n, m, rect(1:2), rect(3:4), yorder);
end
if ( ~isempty(N.dbc) ) % Bottom boundary conditions.
    [bcDown, downVal] = chebop2.constructBC(N.dbc, -1, n, m, rect(1:2), rect(3:4), yorder);
end

%%

% For the down and up BCs we have B^TX = g^T.
By = [ bcUp.'; bcDown.' ];
Gy = [ upVal.'; downVal.' ];
[By, Gy, Py] = canonicalBC(By, Gy);

% For the left and right BCs we have X*B = g. We do the LU to B^T.
Bx = [ bcLeft.'; bcRight.' ];
Gx = [ leftVal.'; rightVal.' ];
[Bx, Gx, Px] = canonicalBC(Bx, Gx);
Bx = Bx.';
Gx = Gx.'; % Now transpose so that X*B = g;

%%
% Construct the RHS of the Sylvester matrix equation.

% Complete the RHS (part of the RHS could have been in the operator):
if ( ~isempty(N.rhs) && ( isa(N.rhs, 'chebfun2') || isa(N.rhs, 'double') ) )
    f = f + N.rhs;
end
E = zeros(m, n);
[n2, n1] = length(f);
F = chebcoeffs2(f);

% Map the RHS to the right ultraspherical space.
lmap = ultraS.convertmat(n1, 0, yorder-1);
rmap = ultraS.convertmat(n2, 0, xorder-1);
F = lmap * F * rmap.';

% Place those coefficients of the forcing function onto the RHS.
n1 = min(n1, m);
n2 = min(n2, n);
E(1:n1,1:n2) = F(1:n1,1:n2);

if ( ~flag ) % Impose boundary conditions.

    % Use the eliminated boundary condition to place zeros in the columns of
    % the matrix equation discretization. There are rk columns to zero out.

    for jj = 1:rk  % For term in the matrix equation.
        [C, E] = zeroDOF(CC{jj,1}, CC{jj,2}, E, By, Gy);
        CC{jj,1} = C;
        [C, E] = zeroDOF(CC{jj,2}, CC{jj,1}, E.', Bx.', Gx.');
        CC{jj,2} = C;
        E = E.';
    end

    % Remove degrees of freedom.
    nn = n - max(xorder, yorder);
    mm = m - max(xorder, yorder);
    df1 = max(0, xorder - yorder);
    df2 = max(0, yorder - xorder);
    for jj = 1:rk
        CC{jj,1} = CC{jj,1}(1:mm, yorder+1:m-df1);
        CC{jj,2} = CC{jj,2}(1:nn, xorder+1:n-df2);
    end
    % Truncation of righthand side.
    rhs = E(1:mm, 1:nn);

else
    rhs = E;
end

% Pass back the eliminated boundary conditions.
bb = {bcLeft bcRight bcUp bcDown};
gg = {leftVal rightVal upVal downVal};

%%
% Check boundary continunity conditions.

% Check BCs at corners:
allbc = 0;
if ( ~isempty(bcUp) && ~isempty(upVal) && ~isempty(bcRight) && ~isempty(rightVal) )
    if ( norm(rightVal(end-4:end),inf) < sqrt(tol) && norm(upVal(end-4:end),inf) < sqrt(tol) )
        allbc = allbc + norm(upVal.'*bcRight - bcUp.'*rightVal);
    end
end
if ( ~isempty(bcUp) && ~isempty(upVal) && ~isempty(bcLeft) && ~isempty(leftVal) )
    if ( norm(leftVal(end-4:end),inf) < sqrt(tol) && norm(upVal(end-4:end),inf) < sqrt(tol) )
        allbc = allbc + norm(upVal.'*bcLeft - bcUp.'*leftVal);
    end
end
if ( ~isempty(bcDown) && ~isempty(downVal) && ~isempty(bcRight) && ~isempty(rightVal) )
    if ( norm(rightVal(end-4:end),inf) < sqrt(tol) && norm(downVal(end-4:end),inf) < sqrt(tol) )
        allbc = allbc + norm(downVal.'*bcRight - bcDown.'*rightVal);
    end
end
if ( ~isempty(bcDown) && ~isempty(downVal) && ~isempty(bcLeft) && ~isempty(leftVal) )
    if ( norm(leftVal(end-4:end),inf)<sqrt(tol) && norm(downVal(end-4:end),inf)<sqrt(tol) )
        allbc = allbc + norm(downVal.'*bcLeft - bcDown.'*leftVal);
    end
end

% [TODO]: Should there be an error if the compatibility conditions do not
% match?
%
% if allbc >= 100*sqrt(tol)
%     s = sprintf('Boundary conditions differ by %1.4f', allbc');
%     warning('CHEBFUN:CHEBOP2:discretize:BCs', s)
% end

end

function B = unconstrainedMatrixEquation(ODE, jj, n, order, dom)
% Construct the unconstrained Matix Equation. Adding in the constraints later.

B = spalloc(n, n, 3*n);
for kk = 1:size(ODE, 1)

    % Get conversion and differentiation matrices:
    S = ultraS.convertmat(n, kk-1, order-1);
    D = ((2./diff(dom))^(kk-1)) * ultraS.diffmat(n, kk-1);

    if ( iscell(ODE(kk,jj)) && isa(ODE{kk,jj}, 'chebfun') )
        % Variable coefficient term:
        c = ODE{kk,jj}.coeffs;
        M = ultraS.multmat(n, c, kk-1);
        A = S * M * D;

    elseif ( iscell(ODE(kk,jj)) && ~isempty(ODE{kk,jj}) )
        % Constant coefficient term in a variable coefficient ODE:
        A = ODE{kk,jj}.* S * D;

    elseif ( isa(ODE(kk,jj),'double') )
        % Constant coefficient term in a constant coefficient ODE:
        A = ODE(kk,jj).* S * D;

    else
        % Empty cell in array, so no term:
        A = zeros(n);

    end

    % Form the ODE operator:
    B = B + A;

end

end

function [B, G, P] = canonicalBC(B, G)
%CANONICALBC   Form a linear combintation of the boundary conditions
%so that they can be used for imposing on the PDE.

P = nonsingularPermute(B);
B = B*P;
[L, B] = lu(B);
G = L \ G;

% Scale so that B is unit upper triangular.
if ( min(size(B)) > 1 )
    D = diag(1./diag(B));
elseif ( ~isempty(B) )
    D = 1./B(1,1);
else
    D = []; % No boundary conditions.
end
B = D*B;
G = D*G;

end

function P = nonsingularPermute(B)
%NONSINGULARPERMUTE   Permute the columns of B to ensure that the principal
%m*m submatrix of B is nonsingular, where m = size(B, 1).
%
% Note: This is needed for solving the matrix equations with linear
% constraints, see DPhil thesis of Alex Townsend (section 6.5).

m = size(B, 1);
k = 1;

% [TODO]: improve this check.
% Try each mxm block in a linear fashion:
while ( rank(B(:,k:m+k-1)) < m )
    k = k+1;
    if ( m+k > size(B, 2) )
        error('CHEBFUN:CHEBOP2:discretize:nonsingularPermute:BCs', ...
            'Boundary conditions are linearly dependent.');
    end
end

P = speye(size(B, 2));
P = P(:,[k:m+k-1, 1:k-1, m+k:end]);

end

function [C1, E] = zeroDOF(C1, C2, E, B, G)
%ZERODOF   Eliminate so degrees of freedom in the matrix equation can be
%removed.

for ii = 1:size(B, 1) % For each boundary condition, zero a column.
    for kk = 1:size(C1, 1)
        if ( abs(C1(kk,ii)) > 10*eps )
            c = C1(kk, ii); % Constant required to zero entry out.
            C1(kk,:) = C1(kk,:) - c*B(ii,:);
            E(kk,:) = E(kk,:) - c*G(ii,:)*C2.';
        end
    end
end

end


text:
function [f] = spm_fx_NMDA(x_V,x_G,P,M)

VN = 60;
%% Pyamidal Cells & interneuron NMDA receptos

mag_block = 1/(1 + 0.2*exp(-0.062*(exp(P.scale_NMDA)*x_V)));
f =  (x_G*(VN - x_V)*mag_block);

author: spm


text:
function p=get_required_precision(d,epsilon)
% Finds the number of decimal positions required to print all data
% elements.
%
% P=GET_REQUIRED_PRECISION(D) returns the least number of decimals
% 'after the period' to represent all numbers D without loss of precision.
% D should be an array with real numbers and can be of any size.
%
% P=GET_REQUIRED_PRECISION(D,EPSILON) allows for an imprecision of
% at most EPSILON.
%
% In the current implementation, this function always returns a value
% between 0 and 7 (inclusive).
%
% Examples: GET_REQUIRED_PRECISION([2 5])       returns 0
%           GET_REQUIRED_PRECISION([2 5.532])   returns 3
%
% NNO Jan 2010

if nargin<2
    epsilon=0; % we want high precision of course
end

% max precision is 7, because that seems more or less the limit of 32 bit
% floats

for p=0:7
    if max(abs(rem(d(:),10^(-p))))<=epsilon
        return
    end
end


text:
function asa103_test01 ( )

%*****************************************************************************80
%
%% TEST01 demonstrates the use of DIGAMMA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01:\n' );
  fprintf ( 1, '  DIGAMMA computes the Digama or Psi function.\n' );
  fprintf ( 1, '  Compare the result to tabulated values.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '          X         ' );
  fprintf ( 1, 'FX                        FX2\n' );
  fprintf ( 1, '                    ' );
  fprintf ( 1, '(Tabulated)               (DIGAMMA)               DIFF\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = psi_values ( n_data );

    if ( n_data == 0 )
      break
    end

    [ fx2, ifault ] = digamma ( x );

    fprintf ( 1, '  %12.4f  %24.16e  %24.16e  %10.4e\n', ...
    x, fx, fx2, abs ( fx - fx2 ) );

  end

  return
end


text:
function tec_io_test01 ( )

%*****************************************************************************80
%
%% TEC_IO_TEST01 tests TEC_READ.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2006
%
%  Author:
%
%    John Burkardt
%
  tec_file_name = 'ell.dat';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEC_IO_TEST01\n' );
  fprintf ( 1, '  TEC_READ can read a TECPLOT ASCII file that contains\n' );
  fprintf ( 1, '  finite element data.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this example, we will read "%s".\n', tec_file_name );

  [ dim_num, node_num, element_num, element_order, node_data_num, ...
    node_coord, element_node, node_data ] = tec_read ( tec_file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension         = %d\n', dim_num );
  fprintf ( 1, '  Number of nodes           = %d\n', node_num );
  fprintf ( 1, '  Number of elements        = %d\n', element_num );
  fprintf ( 1, '  Element order             = %d\n', element_order );
  fprintf ( 1, '  Number of node data items = %d\n', node_data_num );

  r8mat_transpose_print_some ( dim_num, node_num, node_coord, 1, 1, dim_num, 10, ...
    '  Coordinates of first 10 nodes:' );

  i4mat_transpose_print_some ( element_order, element_num, element_node, 1, 1, ...
    element_order, 10, '  Nodes of first 10 elements:' );

  r8mat_transpose_print_some ( node_data_num, node_num, node_data, 1, 1, ...
    node_data_num, 10, '  Node data for first 10 nodes:' );

  return
end


text:
%% Stereo Image Matching
%
% Example of stereo image matching to produce a disparity map and point cloud
% generation.
%
% Resulting .ply file can also be viewed using
% <http://meshlab.sourceforge.net/ MeshLab>.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/python/stereo_match.py>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/stereo_match.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/gpu/stereo_match.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/calib3d/stereoBM/SBM_Sample.cpp>
% * <https://docs.opencv.org/3.2.0/dd/d53/tutorial_py_depthmap.html>
%

%% Theory
%
% Previously, we saw basic concepts like epipolar constraints and other
% related terms. We also saw that if we have two images of same scene, we can
% get depth information from that in an intuitive way. Below is an image and
% some simple mathematical formulas which prove that intuition:
%
% <<https://docs.opencv.org/3.2.0/stereo_depth.jpg>>
%
% The above diagram contains equivalent triangles. Writing their equivalent
% equations will yield us following result:
%
% $$disparity = x - x' = \frac{Bf}{Z}$$
%
% $x$ and $x'$ are the distance between points in image plane corresponding to
% the scene point 3D and their camera center. $B$ is the distance between two
% cameras (which we know) and $f$ is the focal length of camera (already
% known). So in short, the above equation says that the depth of a point in a
% scene is inversely proportional to the difference in distance of
% corresponding image points and their camera centers. So with this
% information, we can derive the depth of all pixels in an image.
%
% So it finds corresponding matches between two images. We have already seen
% how epiline constraint make this operation faster and accurate. Once it
% finds matches, it finds the disparity.
%

%% Code

function stereo_match_demo()
    %% Images
    % load pair of images
    % (SGBM works with either grayscale or color images, BM only grayscale)
    imgL = cv.imread(fullfile(mexopencv.root(),'test','aloeL.jpg'), 'Color',true);
    imgR = cv.imread(fullfile(mexopencv.root(),'test','aloeR.jpg'), 'Color',true);
    subplot(121), imshow(imgL), title('Left')
    subplot(122), imshow(imgR), title('Right')

    %%
    % downscale for faster processing
    scale = 0.5;
    [imgL, imgR] = scale_images(imgL, imgR, scale);
    whos imgL imgR
    [h,w,cn] = size(imgL);

    %%
    % disparity-to-depth mapping 4x4 matrix, used to compute point cloud
    if true
        % manually enter Q matrix
        % (turns points 180 deg around x-axis, so that y-axis looks up)
        f = 0.8*w;  % guess for local length
        Q = [1 0 0 -0.5*w; 0 -1 0 0.5*h; 0 0 0 -f; 0 0 1 0];
    elseif false
        % images must be rectified, if not we rectify and get Q matrix
        % (requires calibrated stereo camera, see stereo_calibration_demo.m)
        intrinsicFile = fullfile(tempdir(), 'stereo_intrinsic.yml');
        extrinsicFile = fullfile(tempdir(), 'stereo_extrinsic.yml');
        [imgL, imgR, Q] = rectify_images(imgL, imgR, scale, ...
            intrinsicFile, extrinsicFile);
    else
        % when empty, point cloud generation is skipped
        Q = [];
    end
    display(Q)

    %% Stereo Matching
    % create stereo matcher, with params tuned for 'aloe' image pair
    win_size = 3;
    min_disp = 16;             % 0
    num_disp = 112 - min_disp; % fix(w/8) + 15
    num_disp = double(bitand(int32(num_disp), int32(-16))); % divisible by 16
    stereo = cv.StereoSGBM('MinDisparity',min_disp, 'NumDisparities',num_disp, ...
        'BlockSize',16, 'P1',8*cn*win_size^2, 'P2',32*cn*win_size^2, ...
        'Disp12MaxDiff',1, 'UniquenessRatio',10, ...
        'SpeckleWindowSize',100, 'SpeckleRange',32, 'Mode','SGBM');
    display(stereo)

    %%
    % compute disparity map (from a pair of rectified stereo images)
    tic, D = stereo.compute(imgL, imgR); toc
    fprintf('16-bit disparity map: min=%d, max=%d\n', min(D(:)), max(D(:)));
    D = single(D) / 16; % fixed-point numbers with 4 fractional bits -> float
    DD = (D - stereo.MinDisparity) / stereo.NumDisparities; % normalized [0,1]
    figure, imshow(DD), colorbar, title('Disparity')

    %% Point Cloud
    if ~isempty(Q)
        % generate 3d point cloud
        ptc = create_point_cloud(D, Q, imgL, false);
        fprintf('%d point cloud\n', ptc.Count);

        % write point cloud to PLY file
        plyfile = fullfile(tempdir(), 'aloe.ply');
        write_point_cloud(ptc, plyfile);
        fprintf('Point cloud saved to: %s\n', plyfile);

        % visualize point cloud
        if ~mexopencv.isOctave()
            %HACK: Octave hangs if we plot too many scatter points
            figure, vis_point_cloud(ptc);
            title('Point Cloud'); xlabel('X'); ylabel('Y'); zlabel('Z');
            %axis([-10 10 -10 10 -20 0])
        end
    end
end

%% Helper functions

function [imgL, imgR] = scale_images(imgL, imgR, scale)
    if scale ~= 1
        if scale < 1
            interpo = 'Area';
        else
            interpo = 'Cubic';
        end
        imgL = cv.resize(imgL, scale, scale, 'Interpolation',interpo);
        imgR = cv.resize(imgR, scale, scale, 'Interpolation',interpo);
    end
end

function [imgL, imgR, Q] = rectify_images(imgL, imgR, scale, intrinsicFile, extrinsicFile)
    % load params from previously calibrated stereo camera
    assert(exist(intrinsicFile, 'file') == 2, 'missing intrinsic params');
    assert(exist(extrinsicFile, 'file') == 2, 'missing extrinsic params');
    I = cv.FileStorage(intrinsicFile);  % intrinsic: M1, D1, M2, D2
    E = cv.FileStorage(extrinsicFile);  % extrinsic: R, T

    % account for new image size by scaling camera matrix: fx, fy, cx, cy
    I.M1 = I.M1 * scale;
    I.M2 = I.M2 * scale;

    % Note: we assume that cameras were calibrated using images of same size
    % as the original image size here, i.e: [I.width, I.height] == sz/scale)
    sz = [size(imgL,2) size(imgL,1)];

    % re-rectify using scaled image size
    RCT = cv.stereoRectify(I.M1, I.D1, I.M2, I.D2, sz, E.R, E.T);
    Q = RCT.Q;

    % apply rectification
    [map11, map12] = cv.initUndistortRectifyMap(I.M1, I.D1, sz, ...
        'P',RCT.P1, 'R',RCT.R1);
    [map21, map22] = cv.initUndistortRectifyMap(I.M2, I.D2, sz, ...
        'P',RCT.P2, 'R',RCT.R2);
    imgL = cv.remap(imgL, map11, map12, 'Interpolation','Linear');
    imgR = cv.remap(imgR, map21, map22, 'Interpolation','Linear');
end

function ptc = create_point_cloud(D, Q, imgL, cvst)
    xyz = cv.reprojectImageTo3D(D, Q);
    mask = repmat(D > min(D(:)), [1 1 3]);  % where disparity was not computed
    xyz = reshape(xyz(mask), [], 3);
    if size(imgL,3) == 3
        clr = reshape(imgL(mask), [], 3);
    else
        clr = repmat(imgL(mask(:,:,1)), 1, 3);
    end

    if nargin < 4, cvst = true; end
    if cvst && ~mexopencv.isOctave() && mexopencv.require('vision')
        ptc = pointCloud(xyz, 'Color',clr);
    else
        %HACK: pointCloud and related functions are not implemented in Octave
        ptc = struct('Location',xyz, 'Color',clr, 'Count',size(xyz,1));
    end
end

function vis_point_cloud(ptc)
    if isobject(ptc)
        pcshow(ptc);
    else
        scatter3(ptc.Location(:,1), ptc.Location(:,2), ptc.Location(:,3), ...
            6, single(ptc.Color)/255, '.')
        axis tight vis3d
        rotate3d on
    end
end

function write_point_cloud(ptc, fname)
    if isobject(ptc)
        pcwrite(ptc, fname, 'Encoding','ascii');
    else
        fid = fopen(fname, 'wt');
        fprintf(fid, 'ply\n');
        fprintf(fid, 'format ascii 1.0\n');
        fprintf(fid, 'element vertex %d\n', ptc.Count);
        fprintf(fid, 'property float x\n');
        fprintf(fid, 'property float y\n');
        fprintf(fid, 'property float z\n');
        fprintf(fid, 'property uchar red\n');
        fprintf(fid, 'property uchar green\n');
        fprintf(fid, 'property uchar blue\n');
        fprintf(fid, 'end_header\n');
        fprintf(fid, '%f %f %f %d %d %d\n', [ptc.Location single(ptc.Color)].');
        fclose(fid);
    end
end


text:
% ABHISHEK MEENA
% Department of Electrical Engineering
% IIT KANPUR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc
clear all

% Gaussian Kernel is taken

%%
% VARIABLE PARAMETERS CHANGE TO GET DIFFERENT RESULTS
hs = 8; %spatial bandwidth
hr= 7;   % range bandwidth
threshold_convergence_mean = 0.25;
bandwidth=[hs,hr];
%%
i=imread('42049.jpg');
[height,width,frame] = size(i);
x=zeros(5,height*width);
%% Going from RGB space to Luv using the function RGB2LUV
for j=1:height
    for l=1:width
        x(1,l+width*(j-1)) = j;
        x(2,l+width*(j-1)) = l;
        [x(3,l+width*(j-1)),x(4,l+width*(j-1)),x(5,l+width*(j-1))] = RGB2LUV(i(j,l,1),i(j,l,2),i(j,l,3));
    end
end
%%
%%
% finding the clusters  and plotting the clusters with their data points

% centres_clusters = centres of all clusters obtained

[centres_clusters,data2cluster,datapoints_cluster_no] = mean_shift_algorithm(x,bandwidth,threshold_convergence_mean);

no_clusters = length(datapoints_cluster_no);
figure(1);

hold on
color_vector = 'bgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrmykbgrcmykbgrcmykbgrcmykbgrcmy';
for k = 1:no_clusters
    clusters_dataset_members = datapoints_cluster_no{k};
    cluster_centres = centres_clusters(:,k);
    plot(x(1,clusters_dataset_members),x(2,clusters_dataset_members),[color_vector(k) '.'])
    plot(cluster_centres(1),cluster_centres(2),'o','MarkerEdgeColor','k','MarkerFaceColor',color_vector(k), 'MarkerSize',13')
end
title(['data points with their cluster centres Gaussian Kernel (hs,hr)=',num2str(hs),',',num2str(hr)]);
hold off
%%
% creating a 3-D RGB matrix and then plotting the 3-D matrix to get the segmented image
[h2,w2] = size(centres_clusters);
zfilter=zeros(5,height*width);
for i12=1:w2
    mem=datapoints_cluster_no{i12,1};
    p1=size(mem);

    for s1=1:p1(1,2)
        zfilter(:,mem(s1))=centres_clusters(:,i12);
    end
end
zluv(:,:,1)=(reshape(zfilter(3,:),width,height))';
zluv(:,:,2)=(reshape(zfilter(4,:),width,height))';
zluv(:,:,3)=(reshape(zfilter(5,:),width,height))';
zrgb = colorspace('Luv->RGB',zluv);
zrgb1=round(zrgb*255);
figure(2)
imshow(zrgb);
title(['Segmented Image (hs,hr)=',num2str(hs),',',num2str(hr)]);
%%

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%IBLOBS	Blob features
%
% F = IBLOBS(IM, OPTIONS) is a vector of RegionFeature objects that
% describe each connected region in the image IM.
%
% Options::
%  'aspect',A        set pixel aspect ratio, default 1.0
%  'connect',C	     set connectivity, 4 (default) or 8
%  'greyscale'	     compute greyscale moments 0 (default) or 1
%  'boundary'        compute boundary (default off)
%  'area',[A1,A2]    accept only blobs with area in the interval A1 to A2
%  'shape',[S1,S2]   accept only blobs with shape in the interval S1 to S2
%  'touch',T         accept only blobs that touch (1) or do not touch (0)
%                    the edge (default accept all)
%  'class',C         accept only blobs of pixel value C (default all)
%
% The RegionFeature object has many properties including:
%
%  uc            centroid, horizontal coordinate
%  vc            centroid, vertical coordinate
%  p             centroid (uc, vc)
%  umin          bounding box, minimum horizontal coordinate
%  umax          bounding box, maximum horizontal coordinate
%  vmin          bounding box, minimum vertical coordinate
%  vmax          bounding box, maximum vertical coordinate
%  area          the number of pixels
%  class         the value of the pixels forming this region
%  label         the label assigned to this region
%  children      a list of indices of features that are children of this feature
%  edgepoint     coordinate of a point on the perimeter
%  edge          a list of edge points 2xN matrix
%  perimeter     edge length (pixels)
%  touch         true if region touches edge of the image
%  a             major axis length of equivalent ellipse
%  b             minor axis length of equivalent ellipse
%  theta         angle of major ellipse axis to horizontal axis
%  shape         aspect ratio b/a (always <= 1.0)
%  circularity   1 for a circle, less for other shapes
%  moments       a structure containing moments of order 0 to 2
%
% References::
%  - Robotics, Vision & Control, Section 13.1,
%    P. Corke, Springer 2011.
% - METHODS TO ESTIMATE AREAS AND PERIMETERS OF BLOB-LIKE OBJECTS: A COMPARISON
%   Luren Yang, Fritz Albregtsen, Tor Lgnnestad and Per Grgttum
%   IAPR Workshop on Machine Vision Applications Dec. 13-15, 1994, Kawasaki
% - Area and perimeter measurement of blobs in discrete binary pictures.
%   Z.Kulpa.
%   Comput. Graph. Image Process., 6:434-451, 1977.
%
% Notes::
% - The RegionFeature objects are ordered by the raster order of the top most
%   point (smallest v coordinate) in each blob.
% - Circularity is computed using the raw perimeter length scaled down by Kulpa's
%   correction factor.
%
% See also RegionFeature, ILABEL, IDISPLABEL, IMOMENTS.

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

function [features,labimg] = iblobs(im, varargin)

	[nr,nc] = size(im);

    opt.area = [0 Inf];
    opt.shape = [0 Inf];
    opt.class = NaN;
    opt.touch = NaN;
    opt.aspect = 1;
    opt.connect = 4;
    opt.greyscale = false;
    opt.moments = false;
    opt.boundary = false;

    opt = tb_optparse(opt, varargin);

    % HACK ilabel should take int image
	%[li,nl,parent,color,edge] = ilabel(im, opt.connect);

    [b,li,nb,a] = bwboundaries(im);
    nl = max(li(:));

	blob = 0;
	for i=0:nl
		binimage = (li == i);

		% determine the blob extent
		[y,x] = find(binimage);
		umin = min(x); umax = max(x);
		vmin = min(y); vmax = max(y);

        % it touches the edge if its parent is 0
		%t = (parent(i) == 0);
        try
        t = (umin == 1) || (vmin == 1) || (umax == nc) || (vmax == nr);
        catch
            t = false;
            if i==0
                continue
            end
        end

        % compute the moments
		if opt.greyscale
			F = imoments(binimage .* im, 'aspect', opt.aspect);
		else
			F = imoments(binimage, 'aspect', opt.aspect);
		end

        % compute shape property, accounting for degenerate case
		if F.a == 0,
			shape = NaN;
		else
			shape = F.b / F.a;
        end

        color = (i <= nb) && (i > 0);

		% apply various filters
		if 	((t == opt.touch) || isnan(opt.touch)) && ...
            ((color == opt.class) || isnan(opt.class)) && ...
			(F.area_ >= opt.area(1)) && ...
			(F.area_ <= opt.area(2)) && ...
			(					...
				isnan(shape) ||			...
				(               ...
					(shape >= opt.shape(1)) &&	...
					(shape <= opt.shape(2))	...
				)				...
			)

            % this blob matches the filter

            % record a perimeter point
            %[y,x] = ind2sub(size(im), edge(i));
            if i==0
                % special work for the background blob, it's edge is not traced by
                % bwboundaries

                % first find a point on the edge that belongs to background
                k = find( (x==1) |(y==1) | (x==nc) | (y==nr) );
                if length(k) >0
                    % blob 0 does touch the edge
                    k = k(1);
                    ep = edgelist(binimage, [x(k) y(k)]); % trace the perim
                else
                    ep = [];
                end
            else
                ep = b{i};
            end

            if ~isempty(ep)
                F.edgepoint = [ep(1,2) ep(1,1)];    % a point on the perimeter
                e = [ep(:,2) ep(:,1)];
                F.edge = e;

                e = diff([e; e(1,:)])';

                % compute length:
                %   - 1 for horizontal/vertical segment
                %   - sqrt(2) for diagnonal
                %
                % Apply Kulpa's correction factor when computing
                % circularity
                kulpa = pi/8*(1+sqrt(2));
                F.perimeter_ = sum( colnorm(e) );

                F.circularity_ = 4*pi*F.area_/ (F.perimeter_*kulpa)^2;
            end
            % set object properties
			F.umin = umin;
			F.umax = umax;
			F.vmin = vmin;
			F.vmax = vmax;
			F.touch_ = t;
            if i > 0
                F.parent = find(a(i,:))+1;
                if isempty(F.parent)
                    F.parent = 1;
                end
                F.children = find(a(:,i))'+1;
            else
                F.parent = 1;
                F.children = [];
            end
			F.shape_ = shape;
            F.label_ = i+1;
            F.class_ = color;

            % save it in the feature vector
			blob = blob+1;
			features(blob) = F;
		end
	end

    if blob == 0
        features = [];
    end

%     % add children property
%     % the numbers in the children property refer to elements in the feature vector
%     for i=1:length(features)
%         parent = features(i).parent;
%         for F=features
%             if F.label == parent
%                 F.children = [F.children i];
%             end
%         end
%     end

	%fprintf('%d blobs in image, %d after filtering\n', nl, blob);
    if nargout > 1
        labimg = li+1;
    end

author: petercorke

repo: machinevision-toolbox-matlab


text:
% version from 102901
function [w,h]=nmf(v,r,verbose,winit,hinit,w_fixvec,h_fixvec)
%
% Jean-Philippe Brunet
% Cancer Genomics
% The Broad Institute
% brunet@broad.mit.edu
%
% This software and its documentation are copyright 2004 by the
% Broad Institute/Massachusetts Institute of Technology. All rights are reserved.
% This software is supplied without any warranty or guaranteed support whatsoever.
% Neither the Broad Institute nor MIT can not be responsible for its use, misuse,
% or functionality.
%
% NMF divergence update equations :
% Lee, D..D., and Seung, H.S., (2001), 'Algorithms for Non-negative Matrix
% Factorization', Adv. Neural Info. Proc. Syst. 13, 556-562.
%
% v (n,m) : N (genes) x M (samples) original matrix
%           Numerical data only.
%           Must be non negative.
%           Not all entries in a row can be 0. If so, add a small constant to the
%           matrix, eg.v+0.01*min(min(v)),and restart.
%
% r       : number of desired factors (rank of the factorization)
%
% verbose : prints iteration count and changes in connectivity matrix elements
%           unless verbose is 0
%
% Note : NMF iterations stop when connectivity matrix has not changed
%        for 10*stopconv interations. This is experimental and can be
%        adjusted.
%
% w    : N x r NMF factor
% h    : r x M NMF factor
%
% winit - initial value for w
% hinit - initial value for h
% h_fixvec & w_fixvec:
% fixvec - binary vector which component should be fixed: eg [2 3] will
% fix second and third component while varying the first...

if ~exist ('w_fixvec', 'var')
    w_fixvec = [];
end

if ~exist ('h_fixvec', 'var')
    h_fixvec = [];
end

% test for negative values in v
if min(min(v)) < 0
    error('matrix entries can not be negative');
    return
end
if min(sum(v,2)) == 0
    error('not all entries in a row can be zero');
    return
end

[n,m]=size(v);
stopconv=40;      % stopping criterion (can be adjusted)
niter = 10000;     % maximum number of iterations (can be adjusted)

cons=zeros(m,m);
consold=cons;
inc=0;
j=0;

%
% initialize random w and h
%
w=rand(n,r);
if exist('winit', 'var')
    if ~isempty(winit)
        w = winit;
        fprintf ('initial values of ''w'' used...\n');
    end
else
    winit = zeros(n,r);
end

h=rand(r,m);
if exist('hinit', 'var')
    if ~isempty(hinit)
        h = hinit;
        fprintf ('initial values of ''h'' used...\n');
    end
else
    hinit = zeros(r,m);
end

if ~isempty(w_fixvec)
    fprintf('Fixing [ ');
    fprintf('%g ', w_fixvec);
    fprintf('] component of ''w''\n');
end

if ~isempty(h_fixvec)
    fprintf('Fixing [ ');
    fprintf('%g ', h_fixvec);
    fprintf('] component of ''h''\n');
end

% oldpar = load('/afs/inf.ed.ac.uk/user/s08/s0880377/project/MATLAB/qdots/h5.mat');
% w = oldpar.w;
% h=oldpar.h;

for i=1:niter

    % divergence-reducing NMF iterations

    x1=repmat(sum(w,1)',1,m);
    h=h.*(w'*(v./(w*h)))./x1;
    h(h_fixvec,:) = hinit(h_fixvec,:);

    x2=repmat(sum(h,2)',n,1);
    w=w.*((v./(w*h))*h')./x2;
    w(:,w_fixvec) = winit(:,w_fixvec);

    % test convergence every 10 iterations

    if(mod(i,10)==0)
        j=j+1;

        % adjust small values to avoid undeflow
        h=max(h,eps);w=max(w,eps);

        % construct connectivity matrix
        [y,index]=max(h,[],1);   %find largest factor
        mat1=repmat(index,m,1);  % spread index down
        mat2=repmat(index',1,m); % spread index right
        cons=mat1==mat2;

        if(sum(sum(cons~=consold))==0) % connectivity matrix has not changed
            inc=inc+1;                     %accumulate count
        else
            inc=0;                         % else restart count
        end
        if verbose                     % prints number of changing elements
            fprintf('\t%d\t%d\t%d\n',i,inc,sum(sum(cons~=consold))),
        end

        if(inc>stopconv)
            break,                % assume convergence is connectivity stops changing
        end

        consold=cons;

    end
end

if i == niter
    fprintf('\nMAximum number of iteration (%g) reached! \n', niter)
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% FGPLVM toolbox
% Version 0.132		Friday 24 Feb 2006 at 19:26
% Copyright (c) 2006 Neil D. Lawrence
%
% CMDSROADDATA This script uses classical MDS to visualise some road distance data.
% DEMOIL1 Oil data with fully independent training conditional.
% DEMOIL2 Oil data with fully independent training conditional, and MLP back constraints.
% DEMOIL3 Oil data with deterministic training conditional.
% DEMOIL4 Oil data with deterministic training conditional, and MLP back constraints.
% DEMOIL5 Oil data with partially independent training conditional.
% DEMOIL6 Oil data with partially independent training conditional, and MLP back constraints.
% DEMROBOTWIRELESSNAVIGATE Take some test data for the robot and navigate with it.
% DEMROBOTTRACES1 Wireless Robot data from University of Washington, with tailored dynamics.
% DEMROBOTWIRELESS1 Wireless Robot data from University of Washington, without dynamics and without back constraints.
% DEMROBOTWIRELESS2 Wireless Robot data from University of Washington, without dynamics and without back constraints.
% DEMROBOTWIRELESS4 Wireless Robot data from University of Washington with dynamics and no back constraints.
% DEMROBOTWIRELESS4 Wireless Robot data from University of Washington with dynamics and back constraints.
% DEMROBOTWIRELESSNAVIGATE Take some test data for the robot and navigate with it.
% DEMSPGP1D1 Do a simple 1-D regression after Snelson & Ghahramani's example.
% DEMSPGP1D2 Do a simple 1-D regression after Snelson & Ghahramani's example.
% DEMSPGP1D3 Do a simple 1-D regression after Snelson & Ghahramani's example.
% DEMSPGP1D4 Do a simple 1-D regression after Snelson & Ghahramani's example.
% DEMSPGP1DPLOT Plot results from 1-D sparse GP.
% DEMSTICK1 Model the stick man using an RBF kernel.
% DEMSTICK2 Model the stick man using an RBF kernel and dynamics.
% DEMSTICK3 Model the stick man using an RBF kernel and mlp back constraints.
% DEMSTICK4 Model the stick man using an RBF kernel and 3-D latent space.
% DEMVOWELS1 Model the vowels data with a 2-D FGPLVM using RBF kernel.
% DEMVOWELS2 Model the vowels data with a 2-D FGPLVM using RBF kernel and back constraints.
% DEMVOWELS3 Model the vowels data with a 2-D FGPLVM using RBF kernel and back constraints, but without PCA initialisation.
% DEMVOWELSISOMAP Model the vowels data with a 2-D FGPLVM using RBF kernel.
% DEMVOWELSLLE Model the vowels data with a 2-D FGPLVM using RBF kernel.
% DEMWALKSITJOG1 Model the stick man using an RBF kernel.
% DEMWALKSITJOG2 Model the stick man using an RBF kernel.
% DEMWALKSITJOG3 Model the stick man using an RBF kernel.
% DEMWALKSITJOG4 Model the stick man using an RBF kernel.
% DEMWALKSITJOGDYNAMICSLEARN Learn the stick man dynamics.
% DYNAMICSTEST Run some tests on the specified dynamics model.
% FGPLVMADDDYNAMICS Add a dynamics kernel to the model.
% FGPLVMCLASSVISUALISE Callback function for visualising data in 2-D.
% FGPLVMCOVGRADSTEST Test the gradients of the covariance.
% FGPLVMCREATE Create a GPLVM model with inducing varibles.
% FGPLVMDISPLAY Display an FGPLVM model.
% FGPLVMDYNAMICSFIELDPLOT 2-D field plot of the dynamics.
% FGPLVMDYNAMICSPLOT 2-D scatter plot of the latent points.
% FGPLVMDYNAMICSPOSTERIORMEANVAR Mean and variances of the posterior at points given by X.
% FGPLVMDYNAMICSRUN Visualise the manifold.
% FGPLVMDYNAMICSSAMPLE Sample a field from the GP.
% FGPLVMEXPANDPARAM Expand a parameter vector into a GP-LVM model.
% FGPLVMEXTRACTPARAM Extract a parameter vector from a GP-LVM model.
% FGPLVMFIELDPLOT 2-D field plot of the dynamics.
% FGPLVMGRADIENT GP-LVM gradient wrapper.
% FGPLVMKERNDYNAMICSSAMPLE Sample a field from a given kernel.
% FGPLVMLOADRESULT Load a previously saved result.
% FGPLVMLOGLIKEGRADIENTS Compute the gradients of the EZFT sparse covariance.
% FGPLVMLOGLIKELIHOOD Log-likelihood for a GP-LVM.
% FGPLVMNEARESTNEIGHBOUR Give the number of errors in latent space for 1 nearest neighbour.
% FGPLVMOBJECTIVE Wrapper function for GPLVM objective.
% FGPLVMOPTIMISE Optimise the inducing variable based kernel.
% FGPLVMOPTIMISEPOINT Optimise the postion of a point.
% FGPLVMOPTIONS Return default options for FGPLVM model.
% FGPLVMPOINTGRADIENT Wrapper function for gradient of a single point.
% FGPLVMPOINTLOGLIKEGRADIENT Log-likelihood gradient for of a point of the GP-LVM.
% FGPLVMPOINTLOGLIKELIHOOD Log-likelihood of a point for the GP-LVM.
% FGPLVMPOINTOBJECTIVE Wrapper function for objective of a single point.
% FGPLVMPOSTERIORMEANVAR Mean and variances of the posterior at points given by X.
% FGPLVMPRINTPLOT Print latent space for learnt model.
% FGPLVMREADFROMFID Load from a FID produced by the C++ implementation.
% FGPLVMREADFROMFILE Load a file produced by the c++ implementation.
% FGPLVMRESULTSDYNAMIC Load a results file and visualise them.
% FGPLVMTEST Test the gradients of the gpCovGrads function and the fgplvm models.
% FGPLVMVISUALISE Visualise the manifold.
% GPBLOCKINDICES Return indices of given block.
% GPCOMPUTEALPHA Update the vector `alpha' for computing posterior mean quickly.
% GPCOMPUTEM Compute the matrix m given the model.
% GPCOVGRADS Sparse objective function gradients wrt Covariance functions for inducing variables.
% GPCOVGRADSTEST Test the gradients of the covariance.
% GPCREATE Create a GP model with inducing varibles/pseudo-inputs.
% GPDATAINDICES Return indices of present data.
% GPDISPLAY Display a Gaussian process model.
% GPDYNAMICSCREATE Create the dynamics model.
% GPDYNAMICSDISPLAY Display a GP dynamics model.
% GPDYNAMICSEXPANDPARAM Place the parameters vector into the model for GP dynamics.
% GPDYNAMICSEXTRACTPARAM Extract parameters from the GP dynamics model.
% GPDYNAMICSLATENTGRADIENTS Gradients of the X vector given the dynamics model.
% GPDYNAMICSLOGLIKEGRADIENTS Gradients of the GP dynamics wrt parameters.
% GPDYNAMICSLOGLIKELIHOOD Give the log likelihood of the dynamics part.
% GPDYNAMICSSAMP Sample from the dynamics for a given input.
% GPDYNAMICSSETLATENTVALUES Set the latent values inside the model.
% GPEXPANDPARAM Expand a parameter vector into a GP model.
% GPEXTRACTPARAM Extract a parameter vector from a GP model.
% GPGRADIENT Gradient wrapper for a GP model.
% GPLOGLIKEGRADIENTS Compute the gradients for the parameters and X.
% GPLOGLIKELIHOOD Compute the log likelihood of a GP.
% GPOBJECTIVE Wrapper function for GP objective.
% GPOPTIMISE Optimise the inducing variable based kernel.
% GPOPTIONS Return default options for GP model.
% GPOUT Evaluate the output of an Gaussian process model.
% GPPOSTERIORGRADMEANVAR Gadient of the mean and variances of the posterior at points given by X.
% GPPOSTERIORMEANCOVAR Mean and covariances of the posterior at points given by X.
% GPPOSTERIORMEANVAR Mean and variances of the posterior at points given by X.
% GPREVERSIBLEDYNAMICSCREATE Create the dynamics model.
% GPREVERSIBLEDYNAMICSDISPLAY Display a GP dynamics model.
% GPREVERSIBLEDYNAMICSEXPANDPARAM Place the parameters vector into the model for GP dynamics.
% GPREVERSIBLEDYNAMICSEXTRACTPARAM Extract parameters from the GP reversible dynamics model.
% GPREVERSIBLEDYNAMICSLATENTGRADIENTS Gradients of the X vector given the dynamics model.
% GPREVERSIBLEDYNAMICSLOGLIKEGRADIENTS Gradients of the GP reversible dynamics wrt parameters.
% GPREVERSIBLEDYNAMICSLOGLIKELIHOOD Give the log likelihood of the dynamics part.
% GPREVERSIBLEOPTIONS Return default options for GP reversible dynamics model.
% GPREVERSIBLEDYNAMICSSAMP Sample from the dynamics for a given input.
% GPREVERSIBLEDYNAMICSSETLATENTVALUES Set the latent values inside the model.
% GPSCALEBIASGRADIENT Compute the gradient of the scale and bias.
% GPUPDATEKERNELS Update the kernels that are needed.
% MODELLATENTGRADIENTS Gradients of the latent variables for dynamics models in the GPLVM.
% MODELSETLATENTVALUES Set the latent variables for dynamics models in the GPLVM.
% ROBONEDYNAMICSCREATE Create the dynamics model.
% ROBONEDYNAMICSDISPLAY Display the robot dynamics model.
% ROBONEDYNAMICSEXPANDPARAM Place the parameters vector into the model for first robot dynamics.
% ROBONEDYNAMICSEXTRACTPARAM Extract parameters from the robot one dynamics model.
% ROBONEDYNAMICSLATENTGRADIENTS Gradients of the X vector given the dynamics model.
% ROBONEDYNAMICSLOGLIKEGRADIENTS Gradients of the robot one dynamics wrt parameters.
% ROBONEDYNAMICSLOGLIKELIHOOD Give the log likelihood of the robot one dynamics part.
% ROBONEDYNAMICSSETLATENTVALUES Set the latent values inside the model.
% ROBTHREEDYNAMICSCREATE Create the dynamics model.
% ROBTHREEDYNAMICSDISPLAY Display the robot dynamics model.
% ROBTHREEDYNAMICSEXPANDPARAM Place the parameters vector into the model for first robot dynamics.
% ROBTHREEDYNAMICSEXTRACTPARAM Extract parameters from the robot three dynamics model.
% ROBTHREEDYNAMICSLATENTGRADIENTS Gradients of the X vector given the dynamics model.
% ROBTHREEDYNAMICSLOGLIKEGRADIENTS Gradients of the robot three dynamics wrt parameters.
% ROBTHREEDYNAMICSLOGLIKELIHOOD Give the log likelihood of the robot three dynamics part.
% ROBTHREEDYNAMICSSETLATENTVALUES Set the latent values inside the model.
% ROBTHREEDYNAMICSSETLATENTVALUES Set the latent values inside the model.
% ROBTWODYNAMICSCREATE Create the dynamics model.
% ROBTWODYNAMICSDISPLAY Display the robot dynamics model.
% ROBTWODYNAMICSEXPANDPARAM Place the parameters vector into the model for first robot dynamics.
% ROBTWODYNAMICSEXTRACTPARAM Extract parameters from the robot two dynamics model.
% ROBTWODYNAMICSLATENTGRADIENTS Gradients of the X vector given the dynamics model.
% ROBTWODYNAMICSLOGLIKEGRADIENTS Gradients of the robot two dynamics wrt parameters.
% ROBTWODYNAMICSLOGLIKELIHOOD Give the log likelihood of the robot one dynamics part.
% ROBTWODYNAMICSSETLATENTVALUES Set the latent values inside the model.

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function s = sharpness(im, metric)
%SHARPNESS   Calculate image sharpness metric.
%
% DESCRIPTION:
%       sharpness returns a scalar metric related the sharpness of the 2D
%       or 3D image matrix defined by im. By default, the metric is based
%       on the Brenner gradient which returns the sum of the centered
%       finite-difference at each matrix element in each Cartesian
%       direction. Metrics calculated using the Sobel operator or the
%       normalised variance can also be returned by setting the input
%       paramater metric.
%
%       For further details, see B. E. Treeby, T. K. Varslot, E. Z. Zhang,
%       J. G. Laufer, and P. C. Beard, "Automatic sound speed selection in
%       photoacoustic image reconstruction using an autofocus approach," J.
%       Biomed. Opt., vol. 16, no. 9, p. 090501, 2011.
%
% USAGE:
%       s = sharpness(im)
%       s = sharpness(im, metric)
%
% INPUTS:
%       im          - 2D or 3D image data to evaluate
%
% OPTIONAL INPUTS:
%       metric      - sharpness metric. Supported values are:
%                       'Brenner'           (default)
%                       'Tenenbaum'
%                       'NormVariance'
%
% OUTPUTS:
%       s           - computed sharpness metric
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 15th January 2012
%       last update - 15th January 2012
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

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

% check for metric input
if nargin == 1
   metric = 'Brenner';
end

switch metric

    case 'Brenner'

        % compute sharpness metric based on the Brenner gradient
        switch numDim(im)
            case 2

                % compute metric
                bren_x = (im(1:end-2, :) - im(3:end, :)).^2;
                bren_y = (im(:, 1:end-2) - im(:, 3:end)).^2;
                s = sum(bren_x(:)) + sum(bren_y(:));

            case 3

                % compute metric
                bren_x = (im(1:end-2, :, :) - im(3:end, :, :)).^2;
                bren_y = (im(:, 1:end-2, :) - im(:, 3:end, :)).^2;
                bren_z = (im(:, :, 1:end-2) - im(:, :, 3:end)).^2;
                s = sum(bren_x(:)) + sum(bren_y(:)) + sum(bren_z(:));

        end

    case 'Tenenbaum'

        % compute sharpness metric based on the Tenenbaum gradient
        switch numDim(im)
            case 2

                % define the 2D sobel gradient operator
                sobel = [-1 0 1; -2 0 2; -1 0 1];

                % compute metric
                s = conv2(sobel, im).^2 + conv2(sobel.', im).^2;
                s = sum(s(:));

            case 3

                % define the 3D sobel gradient operator
                sobel3D(:, :, 1) = [1 2 1; 2 4 2; 1 2 1];
                sobel3D(:, :, 2) = zeros(3);
                sobel3D(:, :, 3) = -sobel3D(:, :, 1);

                % compute metric
                s = convn(im, sobel3D).^2 + convn(im, permute(sobel3D, [3 1 2])).^2 +  convn(im, permute(sobel3D, [2 3 1])).^2;
                s = sum(s(:));

        end

    case 'NormVariance'

        % compute sharpness metric based on the normalised variance
        mu = mean(im(:));
        s = sum((im(:) - mu).^2)./mu;

end

repo: Aero-matlab


text:
%%
% Copyright (c) 2019-present, Mahmoud Afifi
%
% This source code is licensed under the license found in the
% LICENSE file in the root directory of this source tree.
%
% Please, cite the following paper if you use this code:
%
% Mahmoud Afifi and Michael S. Brown. What else can fool deep learning?
% Addressing color constancy errors on deep neural network performance.
% ICCV, 2019
%
% Email: mafifi@eecs.yorku.ca | m.3afifi@gmail.com
%%

%%
% This demo generates ten images from the given image. Each generated file
% will be in the following format: originalName_WB_CS.ext, where
% originalName is the original filename, WB refers to white balance (WB)
% settings, and CS is the camera style.
% WB settings are:
%   - T: Tungsten WB [2850 Kelvin (K)]
%   - F: Fluorescent WB [3800K]
%   - D: Daylight WB [5500K]
%   - C: Cloudy WB [6500K]
%   - S: Shade WB [7500K]
% Camera styes are:
%   - AS: Adobe Standard
%   - CS: Camera Standard
% The generated images and a copy of the input image (optional) will be
% saved in the output directory (output_dir)

%%
clear;
clc

imagename = fullfile('..','images','image1.jpg'); % image filename
output_dir = fullfile('..','results'); % output directory to save the
% generated images and a copy of input image
useGPU = false; %to use GPU
NumOfImgs = 10; % should be less than or equal 10
saveOrig = true; %to save a copy of the original image in output_dir

if NumOfImgs > 10
    error('Cannot generate more than 10 images for each input image');
end

if exist(output_dir,'dir') == 0
    mkdir(output_dir);
end

if useGPU
    load('synthWBmodel_GPU.mat'); % load WB_emulator GPU model
else
    load('synthWBmodel.mat'); % load WB_emulator CPU model
end

I_in = imread(imagename); % read the image
[~,name,ext] = fileparts(imagename);
if saveOrig == true % save a copy of the original image
    imwrite(I_in,fullfile(output_dir,sprintf('%s%s%s',name,'_original',ext)));
end

%%
disp('processing...');
tic
% generate images with synthetic WB effects
out = WB_emulator.generate_wb_srgb(I_in, NumOfImgs);
toc
if useGPU
    out = gather(out); % if GPU is used, convert images to a double tensor
end
disp('done!');
disp('saving...');
for i =1 : size(out,4)
    % save generated images
    imwrite(out(:,:,:,i),fullfile(output_dir,sprintf('%s%s%s',...
        name,WB_emulator.wb_photo_finishing{i},ext)));
end

author: mahmoudnafifi


text:
function [distances,v,w] = convolutionalDistance(p0,p1,areaWeights,kernel,kernelTranspose,options)

% convolutionalDistance - compute entropy-OT ditance and rescaling factors
%   [distances,v,w] = convolutionalDistance(p0,p1,kernel,kernelTranspose, options)
%
% The normalized coupling reads
%       pi = diag(v)*K*diag(w)
% and it should satisfy
%       pi*a  = p1
%   and pi'*a = p0

n = size(p0,1);
options.null = 0;
niter = getoptions(options, 'niter', 100);
tol = getoptions(options, 'tol', 1e-6);
verb = getoptions(options, 'verb', 1);
displayFunction = getoptions(options, 'disp', @(x,y) disp(''));
disp_rate = getoptions(options, 'disp_rate', 10);
disp_time = getoptions(options, 'disp_time', 0);

if nargin<5 || isempty(kernelTranspose);
    kernelTranspose = kernel; % assume symmetry
end

if isempty(areaWeights)
    areaWeights = ones(n,1);
end

p0 = p0+eps;
p1 = p1+eps;

v = ones(size(p0));
w = ones(size(p1));

distances = zeros(size(v,2),1);

A = sum(areaWeights);
aw = bsxfun(@times,areaWeights,w);

for i=1:niter
    if disp_time
        tic
    end

    v = p1 ./ kernel(aw);
    av = bsxfun(@times,areaWeights,v);
    w = p0 ./ kernelTranspose(av);
    aw = bsxfun(@times,areaWeights,w);

    oldDistances = distances;

    ll = @(x) real(log(x));
    lv = av.*ll(v);
    lw = aw.*ll(w);
    distances = sum((log(A)*av+lv).*kernel(aw),2) + sum(av.*kernel(lw),2);

    change = norm(oldDistances-distances,'fro');

    if verb==1, fprintf('Iteration %d:  %g\n', i, change);
    elseif verb==2, progressbar(i,niter); end

    if isa(displayFunction,'function_handle') && mod(i,disp_rate)==1
        displayFunction(v,w); drawnow;
    end

    if change<tol && i > 2
        if verb==2, progressbar(niter,niter); end
        break;
    end

    if disp_time, toc, end
end

distances = sqrt(max(distances,0));
fprintf('\n');

author: gpeyre


text:
function x_new = sls1mm(f, x_0, opt, gradf, varargin)
%SLS1MM  Markov Chain Monte Carlo sampling using Slice Sampling
%
%  Description
%    SLS1MM is faster streamlined version of SLS for generating samples
%    from one dimensional distribution with known min-max limits.
%
%    SAMPLE = SLS1MM(F, X, OPTIONS) uses slice sampling to generate
%      single value from the *one dimensional* distribution P ~ EXP(-F),
%      where F is the first argument to SLS1MM. Markov chain starts from
%      point X and the sampling is made using min-max slice sampling.
%      See SLS1MM_OPT for details.
%
%    SAMPLES = SLS1MM(F, X, OPTIONS, [], P1, P2, ...) allows additional
%      arguments to be passed to F(). The fourth argument is ignored,
%      but included for compatibility with HMC2 and the optimisers.
%
%  See SLS1MM_OPT for the optional parameters in the OPTIONS structure.
%  Note, that unlike in SLS, missing fields give an error.
%
%  See also
%    SLS1MM_OPT, SLS

%  Based on "Slice Sampling" by Radford M. Neal in "The Annals of Statistics"
%  2003, Vol. 31, No. 3, 705-767, (c) Institute of Mathematical Statistics, 2003

%       Copyright (c) 2003-2004 Toni Auranen
%       Copyright (c) 2004 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% Set up some variables
maxiter = opt.maxiter;
l = opt.mmlimits(1);
r = opt.mmlimits(2);

% Generate sample
y = -f(x_0,varargin{:}) + log(rand);
x_new = x_0;
for iter=1:maxiter
  x_new = l + (r-l).*rand;
  y_new = -f(x_new,varargin{:});
  if y < y_new
    return;
  end
  if x_new < x_0
    l = x_new;
  else
    r = x_new;
  end
end
if iter+1 > maxiter
  fprintf('Maximum number (%d) of iterations reached during shrinkage.\n',maxiter);
  error('Check function F, decrease the interval ''mmlimits'' or increase the value of ''maxiter''.');
end


text:
function sct_merge_nii(flist,output)
% sct_merge_nii(flist,output)
% merge niftis in the 4th dimenstion.
% sct_merge_nii(sct_tools_ls('*dwi.nii'),'DWI.nii.gz') --> always .nii.gz

sct_unix(['fslmerge -t ' output ' ' strjoin(flist)])

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [NEIGH,posInfo] = comp_filterbankneighbors(a,M,N,do_real)

%this function is called by filterbankconstphase

chanStart = [0;cumsum(N)];

NEIGH = zeros(6,chanStart(end));

%Horizontal neighbors
for kk = 1:M
  NEIGH(1,chanStart(kk)+1) = chanStart(kk)+2;
  NEIGH(1,chanStart(kk+1)) = chanStart(kk+1)-1;
  NEIGH(1:2,chanStart(kk)+(2:N(kk)-1)) = chanStart(kk)+[(1:N(kk)-2);(3:N(kk))];
end

%Vertical neighbors
%Set time distance limit
LIM = .8;

%One channel higher
for kk = 1:M-1
  aTemp = a(kk)/a(kk+1);
  POSlow = chanStart(kk+1)+min(max(0,ceil(((0:N(kk)-1)-LIM)*aTemp)),N(kk+1)-1);
  POShigh = chanStart(kk+1)+max(0,min(floor(((0:N(kk)-1)+LIM)*aTemp),N(kk+1)-1));

%   for ll = 1:N(kk)
%     tmpIdx = (POSlow(ll):POShigh(ll))+1;
%     NEIGH((5:4+numel(tmpIdx)),chanStart(kk)+ll) = tmpIdx.';
%   end

NEIGH(5,chanStart(kk)+(1:N(kk))) = POSlow + 1;
NEIGH(6,chanStart(kk)+(1:N(kk))) = POShigh + 1;
end
if ~do_real
    aTemp = a(M)/a(1);
    POSlow = chanStart(1)+min(max(0,ceil(((0:N(M)-1)-LIM)*aTemp)),N(1)-1);
    POShigh = chanStart(1)+max(0,min(floor(((0:N(M)-1)+LIM)*aTemp),N(1)-1));

    NEIGH(5,chanStart(M)+(1:N(M))) = POSlow + 1;
    NEIGH(6,chanStart(M)+(1:N(M))) = POShigh + 1;
end
NEIGH(6,NEIGH(6,:)==NEIGH(5,:)) = 0;

%One channel lower
for kk = 2:M
  aTemp = a(kk)/a(kk-1);
  POSlow = chanStart(kk-1)+min(max(0,ceil(((0:N(kk)-1)-LIM)*aTemp))',N(kk-1)-1);
  POShigh = chanStart(kk-1)+max(0,min(floor(((0:N(kk)-1)+LIM)*aTemp),N(kk-1)-1)');

%   for ll = 1:N(kk)
%     tmpIdx = (POSlow(ll):POShigh(ll))+1;
%     NEIGH((3:2+numel(tmpIdx)),chanStart(kk)+ll) = tmpIdx.';
%   end
NEIGH(3,chanStart(kk)+(1:N(kk))) = POSlow + 1;
NEIGH(4,chanStart(kk)+(1:N(kk))) = POShigh + 1;
end
if ~do_real
    aTemp = a(1)/a(M);
    POSlow = chanStart(M)+min(max(0,ceil(((0:N(1)-1)-LIM)*aTemp)),N(M)-1);
    POShigh = chanStart(M)+max(0,min(floor(((0:N(1)-1)+LIM)*aTemp),N(M)-1));

    NEIGH(3,chanStart(1)+(1:N(1))) = POSlow + 1;
    NEIGH(4,chanStart(1)+(1:N(1))) = POShigh + 1;
end
NEIGH(4,NEIGH(4,:)==NEIGH(3,:)) = 0;

posInfo = zeros(chanStart(end),2);
for kk = 1:M
    posInfo(chanStart(kk)+(1:N(kk)),:) = [(kk-1)*ones(N(kk),1),(0:N(kk)-1)'.*a(kk)];
end
posInfo = posInfo.';

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function p_data = p01_data ( dim_num, data_num )

%*****************************************************************************80
%
%% P01_DATA returns the data for problem p01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension of the dependent
%    variables.
%
%    Input, integer DATA_NUM, the number of data points.
%
%    Output, real P_DATA(DIM_NUM,DATA_NUM), the data.
%
  p_data = [ ...
    0.0, 4.0; ...
    1.0, 5.0; ...
    2.0, 6.0; ...
    4.0, 6.0; ...
    5.0, 5.0; ...
    6.0, 3.0; ...
    7.0, 1.0; ...
    8.0, 1.0; ...
    9.0, 1.0; ...
   10.0, 3.0; ...
   11.0, 4.0; ...
   12.0, 4.0; ...
   13.0, 3.0; ...
   14.0, 3.0; ...
   15.0, 4.0; ...
   16.0, 4.0; ...
   17.0, 3.0; ...
   18.0, 0.0 ]';

 return
end


text:
function v = rotate_outer(m,rot,varargin)
% rotate crystal directions
%
% Input
%  m - @Miller
%  ori - @orientation
%
% Output
%  v - vector3d
%

% ensure that the rotations have the right reference frame
if isa(rot,'orientation') && nargin == 2
  rot = m.CS.ensureCS(rot);
end

v = rotate_outer@vector3d(m,rot);


text:
%% Remove duplicate segments
function smaps2 = msPruneSegments(smaps)

smaps2 = smaps;
segment_list = {};

for m = 1:size(smaps, 2)
    for k = 1:max(smaps(:, m))
        ind = smaps(:, m)==k;
        [segment_list, isnew] = prune_add2list(segment_list, find(ind));
        if ~isnew
            smaps2(ind, m) = 0;
        end
    end
    uid = unique(smaps2(:, m));
    uid = setdiff(uid, 0);
    for k = 1:numel(uid)
        ind = smaps2(:, m)==uid(k);
        smaps2(ind, m) = k;
    end
end

function [seglist, isnew] = prune_add2list(seglist, ids)

isnew = true;
for k = 1:numel(seglist)
    seg = seglist{k};
    if numel(seg)==numel(ids) && all(seg==ids)
        isnew = false;
        return;
    end
end
seglist{end+1} = ids;

author: zouchuhang

repo: LayoutNet


text:
function FitInfo = GLG_EM_wrapper( theta, tree, varargin )

% Fit GLG model to 1D wavelet tree with composite EM algorithms
% Run the two EM algorithms for the GLG model until effective convergence.
%
% Syntax:
%   FitInfo = GLG_EM_wrapper( theta, tree, ... )
%
% Input:
%   theta      : Parameter array
%
%   tree       : Wavelet tree from DWT2_TO_CELL
%
%
% The optional arguments are specified as NAME, VALUE.
%
%   conv_thres : Convergence threshold for the log-likelihood function
%
%   conv       : How should CONV_THRES be applied: 'absolute' og 'relative'
%                Default: 'absolute'
%
%   max_itr    : The maximum number of iterations
%                Default: 50
%
%   nodes      : The number of nodes used in quadrature rule
%                Default: 20
%
%   noise_dev  : Standard dev of noise; 0 if not relevant.
%                Default: 0
%
%
% Output:
%   FitInfo    : Structure with model, initial model and likelihood
%
%
% See also: EM_ROOT, EM_TREE

% --------------------------------------------------------------------
% 													Set default values
% --------------------------------------------------------------------

p = inputParser;

p.addRequired( 'theta' );
p.addRequired( 'tree' );

p.addOptional( 'conv_thres', 0.01, @isnumeric );
p.addOptional( 'conv', 'absolute', @ischar );
p.addOptional( 'max_itr', 50, @isnumeric );
p.addOptional( 'nodes', 20, @isnumeric );
p.addOptional( 'noise_dev', 0, @isnumeric );

parse( p, tree, theta, varargin{:} );

conv_thres = p.Results.conv_thres;
conv       = p.Results.conv;
max_itr    = p.Results.max_itr;
nodes      = p.Results.nodes;
noise_dev  = p.Results.noise_dev;

% --------------------------------------------------------------------
% 													 Initialize output
% --------------------------------------------------------------------

no_levels = length(tree)-1;
no_dir = size( theta, 3 );

ll = cell(no_levels, no_dir);

% Initialize output
FitInfo = struct(...
    'init_model', theta, ...
    'model', zeros(no_levels, 5, no_dir), ...
    'll', [] ...
    );

% --------------------------------------------------------------------
%                                    Run EM algorithm on the top level
% --------------------------------------------------------------------

% Inform user
fprintf('Running EM algorithm...\n');

if no_dir > 1
    fprintf(repmat( ' ', 1, 10 ));

    num_itr_space = numel(num2str(max_itr))-numel(num2str(no_dir))+3;

    itr_space = repmat( ' ', 1, num_itr_space );
    for d = 1:no_dir
        fprintf( ['Direction %u' itr_space], d );
    end

    fprintf('\n');
end

fprintf('Level 1   ');

for d = 1:no_dir
    w = tree{2}{d}(:);

    % Inform user
    fprintf('Iteration:  ');

    % Run EM algorithm
    [theta, ll{1,d}] = EM_root( theta, w, d, nodes, max_itr, noise_dev, conv_thres, conv );

    % Prepare for iteration count for next direction
    if no_dir > 1
        count = length(ll{1,d});
        fprintf( repmat(' ', 1, numel(num2str(max_itr))-numel(num2str(count))+2) );
    end
end

% --------------------------------------------------------------------
%                             Run EM algorithm on the remaining levels
% --------------------------------------------------------------------

for l = 2:no_levels
    fprintf('\nLevel %u   ', l);

    for d = 1:no_dir
        % Arrange coefficients
        w_parent = tree{l}{d}(:);

        % For images!
        w_child = zeros( numel(w_parent), 4 );
        w_child(:,1) = reshape( tree{l+1}{d}(1:2:end, 1:2:end), 1, [] );
        w_child(:,2) = reshape( tree{l+1}{d}(1:2:end, 2:2:end), 1, [] );
        w_child(:,3) = reshape( tree{l+1}{d}(2:2:end, 1:2:end), 1, [] );
        w_child(:,4) = reshape( tree{l+1}{d}(2:2:end, 2:2:end), 1, [] );

        % Inform user
        fprintf('Iteration:  ');

        % Run EM algorithm
        [theta, ll{l,d}] = EM_tree( theta, w_parent, w_child, l, d, nodes, max_itr, noise_dev, conv_thres, conv );

        % Prepare for iteration count for next direction
        if no_dir > 1
            count = length(ll{l,d});
            fprintf( repmat(' ', 1, numel(num2str(max_itr))-numel(num2str(count))+2) );
        end
    end
end

fprintf('\n\n');

% Final model
FitInfo.model = theta;
FitInfo.ll = ll;

% Test if any likelihood function has decreased
[row, col] = find(cellfun( @(ll) any(ll(2:end) < ll(1:end-1)), ll ));

if ~isempty( row )
    fprintf('Likelihood functions decreased on level and direction:\n');
    disp( [row col] );
end

end

author: Sable

repo: mcbench-benchmarks


text:
%------------------------------ PolyMesher -------------------------------%
% Ref: C Talischi, GH Paulino, A Pereira, IFM Menezes, "PolyMesher: A     %
%      general-purpose mesh generator for polygonal elements written in   %
%      Matlab," Struct Multidisc Optim, DOI 10.1007/s00158-011-0706-z     %
%-------------------------------------------------------------------------%
function [x] = SuspensionDomain(Demand,Arg)
  BdBox = [-2 24 -2 24];
  switch(Demand)
    case('Dist');  x = DistFnc(Arg,BdBox);
    case('BC');    x = BndryCnds(Arg{:},BdBox);
    case('BdBox'); x = BdBox;
    case('PFix');  x = FixedPoints(BdBox);
  end
%----------------------------------------------- COMPUTE DISTANCE FUNCTIONS
function Dist = DistFnc(P,BdBox)
  d1  = dRectangle(P,0,18.885,0,14.56);
  d2  = dLine(P,18.885,1.3030,4,0);
  d3  = dLine(P,3.92,14.56,6.1699,6.88);
  d4  = dLine(P,9.8651,4.0023,18.885,3.70);
  d5  = dLine(P,4,0,0,4);
  d13 = dLine(P,0,14,3.92,14.56);
  d14 = dCircle(P,10,8,4);
  d15 = dLine(P,9.8651,4.0023,6.1699,6.88);
  d   = dDiff(dDiff(dDiff(dDiff(d1,d2),d5),d13),...
        dUnion(dDiff(dIntersect(d3,d4),d15),d14));
  d6  = dCircle(P,2,2,2);
  d7  = dCircle(P,4,2,2);
  d8  = dCircle(P,2,4,2);
  d   = dUnion(d,dUnion(d6,dUnion(d7,d8)));
  d9  = dCircle(P,2,14,2);
  d10 = dCircle(P,2,16,2);
  d11 = dCircle(P,18.885,2.5,1.2);
  d12 = dCircle(P,20,2.5,1.2);
  Dist = dUnion(d,dUnion(d9,dUnion(d10,dUnion(d11,d12))));
%---------------------------------------------- SPECIFY BOUNDARY CONDITIONS
function [x] = BndryCnds(Node,Element,BdBox)
  CornerCircle = sqrt((Node(:,1)-2.0).^2+(Node(:,2)-2.0).^2);
  [foo,CornerCircle] = sort(CornerCircle);
  UpperCircle = sqrt((Node(:,1)- 2.0).^2+(Node(:,2)-16.0).^2);
  [foo,UpperCircle] = sort(UpperCircle);
  Supp = ones(2,3);
  Supp(1,:) = [CornerCircle(1) 1 1];
  Supp(2,:) = [UpperCircle(1)  1 0];
  RightCircle = sqrt((Node(:,1)-20.0).^2+(Node(:,2)-2.5).^2);
  [foo,RightCircle] = sort(RightCircle);
  Load = ones(1,3);
  Load(1,:) = [RightCircle(1) -8 -1];
  x = {Supp,Load};
%----------------------------------------------------- SPECIFY FIXED POINTS
function [PFix] = FixedPoints(BdBox)
  PFix = [ 2   2;
           2  16;
          20 2.5];
%-------------------------------------------------------------------------%

repo: ifem


text:
function test_bug1529

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_plot_ortho

[ftver, ftpath] = ft_version;

%% Load the SPMtemplate
mritem = ft_read_mri([ftpath '/external/spm8/templates/T1.nii']);
% mrichecked = ft_determine_coordsys(mritem);
mritem.coordsys = 'spm'; % define the coordinate system (why is this not part of ft_read_mri?)

%% Display the Anatomy via plot_ortho (defaults)
figure('name','MRI TEM');
ft_plot_ortho(mritem.anatomy);

%% Display the Anatomy via plot_ortho (at a certain location )
location_voxel = [46 64 37];
location_str   = sprintf('x: %.0f y: %.0f z:%.0f',location_voxel(1:3));
figure('name',['MRI TEM VOX' location_str ]);
ft_plot_ortho(mritem.anatomy,'location',location_voxel(1:3));

%% Display the Anatomy via plot_ortho (at the same certain location, using external transform )
location_mri = [0 0 0];
location_voxel = mritem.transform\[location_mri 1]';
location_str   = sprintf('x: %.0f y: %.0f z:%.0f',location_voxel(1:3));
figure('name',['MRI TEM MRI EXTERN ' location_str]);
ft_plot_ortho(mritem.anatomy,'location',location_voxel(1:3));

%% Display the Anatomy via plot_ortho (at the same certain location, using internal transform )
location_mri = [0 0 0];
location_str   = sprintf('x: %.0f y: %.0f z:%.0f',location_mri(1:3));
figure('name',['MRI TEM MRI INTERN ' location_str]);
ft_plot_ortho(mritem.anatomy,'location',location_mri(1:3),'transform',mritem.transform);

%% Display the Anatomy via plot_ortho (at the same certain location, using inverse internal transform )
location_mri = [0 0 0];
location_str   = sprintf('x: %.0f y: %.0f z:%.0f',location_mri(1:3));
figure('name',['MRI TEM MRI INTERN ' location_str]);
ft_plot_ortho(mritem.anatomy,'location',location_mri(1:3),'transform',inv(mritem.transform));

author: fieldtrip


text:
%SF_TRI_P5 Fifth order Lagrange shape functions for triangles (P5).
%
%   [ VBASE, NLDOF, XLDOF, SFUN ] = SF_TRI_P5( I_EVAL, N_SDIM, N_VERT, I_DOF, XI, AINVJAC, VBASE )
%   Evaluates conforming fifth order P5 Lagrange shape functions on 2D triangular elements
%   with values defined in the nodes, edges, and center. XI are Barycentric coordinates.
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       i_eval      scalar:  1             Evaluate function values
%                           >1             Evaluate values of derivatives
%       n_sdim      scalar: 2              Number of space dimensions
%       n_vert      scalar: 3              Number of vertices per cell
%       i_dof       scalar: 1-21           Local basis function to evaluate
%       xi          array [3,1]            Local coordinates of evaluation point
%       aInvJac     [n,6]                  Inverse of transformation Jacobian
%       vBase       [n]                    Preallocated output vector
%                                                                                .
%       Output      Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       vBase       [n]                    Evaluated function values
%       nLDof       [4]                    Number of local degrees of freedom on
%                                          vertices, edges, faces, and cell interiors
%       xLDof       [3,n_ldof]             Local coordinates of local dofs
%       sfun        string                 Function name of called shape function
%
%   See also SF_TRI_P1

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function r8mat_normal_ab_test ( )

%*****************************************************************************80
%
%% R8MAT_NORMAL_AB_TEST tests R8MAT_NORMAL_AB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 March 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT_NORMAL_AB_TEST\n' );
  fprintf ( 1, '  R8MAT_NORMAL_AB computes a scaled pseudonormal matrix.\n' );

  m = 5;
  n = 4;
  mu = 100.0;
  sigma = 5.0;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  MU = %g\n', mu );
  fprintf ( 1, '  SIGMA = %g\n', sigma );
  fprintf ( 1, '  SEED = %d\n', seed );

  [ r, seed ] = r8mat_normal_ab ( m, n, mu, sigma, seed );

  r8mat_print ( m, n, r, '  The matrix:' );

  return
end