text:
function [ ltri, ledg, rtri, redg ] = vbedg ( x, y, node_num, node_xy, ...
  triangle_num, triangle_node, triangle_neighbor, ltri, ledg, rtri, redg )

%*****************************************************************************80
%
%% VBEDG determines which boundary edges are visible to a point.
%
%  Discussion:
%
%    The point (X,Y) is assumed to be outside the convex hull of the
%    region covered by the 2D triangulation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2008
%
%  Author:
%
%    Original FORTRAN77 version by Barry Joe,
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Barry Joe,
%    GEOMPACK - a software package for the generation of meshes
%    using geometric algorithms,
%    Advances in Engineering Software,
%    Volume 13, pages 325-331, 1991.
%
%  Parameters:
%
%    Input, real X, Y, the coordinates of a point outside the
%    convex hull of the current triangulation.
%
%    Input, integer NODE_NUM, the number of points.
%
%    Input, real NODE_XY(2,NODE_NUM), the coordinates of the
%    vertices.
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input, integer TRIANGLE_NODE(3,TRIANGLE_NUM), the triangle incidence list.
%
%    Input, integer TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), the triangle neighbor list; negative
%    values are used for links of a counter clockwise linked list of boundary
%    edges;
%      LINK = -(3*I + J-1) where I, J = triangle, edge index.
%
%    Input/output, integer LTRI, LEDG.  If LTRI /= 0 then these values are
%    assumed to be already computed and are not changed, else they are updated.
%    On output, LTRI is the index of boundary triangle to the left of the
%    leftmost boundary triangle visible from (X,Y), and LEDG is the boundary
%    edge of triangle LTRI to the left of the leftmost boundary edge visible
%    from (X,Y).  1 <= LEDG <= 3.
%
%    Input/output, integer RTRI.  On input, the index of the boundary triangle
%    to begin the search at.  On output, the index of the rightmost boundary
%    triangle visible from (X,Y).
%
%    Input/output, integer REDG, the edge of triangle RTRI that is visible
%    from (X,Y).  1 <= REDG <= 3.
%
  dim_num = 2;
%
%  Find the rightmost visible boundary edge using links, then possibly
%  leftmost visible boundary edge using triangle neighbor information.
%
  if ( ltri == 0 )
    ldone = 0;
    ltri = rtri;
    ledg = redg;
  else
    ldone = 1;
  end

  while ( 1 )

    l = -triangle_neighbor(redg,rtri);
    t = floor ( l / 3 );
    e = mod ( l, 3 ) + 1;
    a = triangle_node(e,t);

    if ( e <= 2 )
      b = triangle_node(e+1,t);
    else
      b = triangle_node(1,t);
    end

    lr = lrline ( x, y, node_xy(1,a), node_xy(2,a), node_xy(1,b), ...
      node_xy(2,b), 0.0 );

    if ( lr <= 0 )
      break;
    end

    rtri = t;
    redg = e;

  end

  if ( ldone )
    return
  end

  t = ltri;
  e = ledg;

  while ( 1 )

    b = triangle_node(e,t);
    e = i4_wrap ( e-1, 1, 3 );

    while ( 0 < triangle_neighbor(e,t) )

      t = triangle_neighbor(e,t);

      if ( triangle_node(1,t) == b )
        e = 3;
      elseif ( triangle_node(2,t) == b )
        e = 1;
      else
        e = 2;
      end

    end

    a = triangle_node(e,t);

    lr = lrline ( x, y, node_xy(1,a), node_xy(2,a), node_xy(1,b), ...
      node_xy(2,b), 0.0 );

    if ( lr <= 0 )
      break
    end

  end

  ltri = t;
  ledg = e;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%DRAWAXIS  Draw coordinate system axis from pose estimation
%
%     img = cv.drawAxis(img, cameraMatrix, distCoeffs, rvec, tvec, length)
%
% ## Input
% * __img__ input image. It must have 1 or 3 channels. In the output, the
%   number of channels is not altered.
% * __cameraMatrix__ input 3x3 floating-point camera matrix
%   `A = [fx 0 cx; 0 fy cy; 0 0 1]`.
% * __distCoeffs__ vector of distortion coefficients
%   `[k1,k2,p1,p2,k3,k4,k5,k6,s1,s2,s3,s4]` of 4, 5, 8 or 12 elements.
% * __rvec__ rotation vector of the coordinate system that will be drawn.
% * __tvec__ translation vector of the coordinate system that will be drawn.
% * __length__ length of the painted axis in the same unit as `tvec`
%   (usually in meters).
%
% ## Output
% * __img__ output image.
%
% Given the pose estimation of a marker or board, this function draws the
% axis of the world coordinate system, i.e. the system centered on the
% marker/board. Useful for debugging purposes.
%
% See also: cv.estimatePoseBoard, cv.estimatePoseSingleMarkers,
%  cv.estimatePoseCharucoBoard, cv.Rodrigues, cv.line, cv.projectPoints
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [A, C, Q, R, initx, initV] = ensure_AR(A, C, Q, R, initx, initV, k, obs, diagonal)
%
% Ensure that the system matrices have the right form for an autoregressive process.

ss = length(A);
if nargin<8, obs=ones(ss, 1); end
if nargin<9, diagonal=0; end

[coef, C] = SS_to_AR(A, Q, k, diagonal);
[A, C, Q, R, initx, initV] = AR_to_SS(coef, C, obs);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% runpca() -  perform principal component analysis (PCA) using singular value
%             decomposition (SVD) using Matlab svd() or svds()
%                        >> inv(eigvec)*data = pc;
% Usage:
%    >> [pc,eigvec,sv] = runpca(data);
%    >> [pc,eigvec,sv] = runpca(data,num,norm)
%
% Inputs:
%   data   - input data matrix (rows are variables, columns observations)
%   num    - number of principal comps to return  {def|0|[] -> rows in data}
%   norm   - 1/0 = do/don't normalize the eigvec's to be equivariant
%                                                {def|0 -> no normalization}
% Outputs:
%   pc     - the principal components, i.e.        >> inv(eigvec)*data = pc;
%   eigvec - the inverse weight matrix (=eigenvectors). >> data = eigvec*pc;
%   sv     - the singular values (=eigenvalues)
%
% Author: Colin Humphries, CNL / Salk Institute, 1997
%
% See also: runica()

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

% 01/31/00 renamed runpca() and improved usage message -sm
% 01-25-02 reformated help & license, added links -ad

function [pc,A,sv]= runpca2(data,K,dosym)
[chans,frames] = size(data);

if nargin < 3 || K < chans
    dosym = 1;
end
if nargin < 2
    K = chans;
end
if nargin < 1
  help runpca
  return
end

% remove the mean
for i = 1:chans
  data(i,:) = data(i,:) - mean(data(i,:));
end

% do svd
[U,S,V] = svd(data*data'/frames); % U and V should be the same since data*data' is symmetric

sv = sqrt(diag(S));
if dosym == 1
    pc = pinv(diag(sv(1:K))) * V(:,1:K)' * data;
    A = U(:,1:K) * diag(sv(1:K));
else
    pc = U * pinv(diag(sv)) * V' * data;
    A = U * diag(sv) * V';
end

author: PatternRecognition

repo: OpenBMI


text:
% SPINV_LDLCHOL - Evaluate the sparse inverse matrix given the LDL
%                 decomposition.
%
% For a sparse symmetric positive definite matrix C,
%
%   Z = SPINV(LD)
%
% where LD = LDLCHOL(C) and [Z]_ij = [inv(C)]_ij for such ij that [C]_ij
% is non-zero.
%
% See Vanhatalo and Vehtari (2008) for details.

% Copyright (c) 2008      Jarno Vanhatalo

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

function Z = spinv_ldlchol(LD)

Z = spinv(LD, 1);
return

    n = size(LD,1);

    % The mex-file was not available, so evaluate the sparse inverse here

    % TODO:
    % For now, just evaluate the matrix A and then run..
    %[L,D] = ldlsplit(LD);
    %A = L*D*L';

    if nargin < 2
      q = 1:n;
    end

    %[LD, p, q] = ldlchol(A);

    [I,J,ld] = find(LD);
    temp = [I(:) J(:) ; J(:) I(:)];
    temp = sortrows(unique(temp,'rows'),2);
    Iz = temp(:,1); Jz = temp(:,2);

    % Find the column starting points
    a1=zeros(n,1);
    a2 = cumsum(histc(J,1:n));
    a1(1) = 1; a1(2:end) = a2(1:end-1) + 1;
    az1=zeros(n,1);
    az2 = cumsum(histc(Jz,1:n));
    az1(1) = 1; az1(2:end) = az2(1:end-1) + 1;

    for j=1:n
        indaz{j} = az1(j):az2(j);
        indIz{j} = Iz(indaz{j})';
    end

    % Evaluate the sparse inverse
    z = zeros(size(Iz));
    z(end) = 1./ld(end);
    % Allocate memory
    cindit=zeros(n,1);
    for jj = n-1:-1:1
        fil = ld(a1(jj)+1:a1(jj+1)-1);
        fi = I(a1(jj)+1:a1(jj+1)-1);
        lfi = length(fi);
        Zt = zeros(lfi,lfi);
        indz = cumsum(histc(indIz{jj},[0 ; fi]));
        indz = az1(jj) + indz(1:end-1);

        i4=0;
        for i1 = 1:lfi
            cind1=indaz{fi(i1)};
            Icind1=indIz{fi(i1)};
            indfi = lfi;
            i2=length(Icind1);
            go = true;
            while go
                if Icind1(i2)==jj  % Find the indeces for the jj'th rows in fi columns
                    i4=i4+1;
                    cindit(i4)=cind1(i2);
                    go = false;
                end
                if indfi >= 1 && fi(indfi) == Icind1(i2) % Find the indeces for the fi'th rows in i2'nd columns
                    Zt(indfi,i1) = z(cind1(i2));
                    indfi = indfi-1;
                end
                i2 = i2-1;
            end
        end
        % remove extras
        cindi=cindit(1:i4);

        zij = -fil'*Zt;
        z(cindi) = zij;
        z(indz) = zij;
        zij = 1./ld(a1(jj)) - fil'*z(indz);
        z(az1(jj)-1+find(indIz{jj}==jj,1)) = zij;
    end

    Z = sparse(Iz,Jz,z);
    r(q) = 1:n;
    Z = Z(r,r);

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function model = presolve_bounds_from_domains(model);

% Sigmonial with non-integer powers must be positive
sigmonials = find((model.variabletype == 4));
for i = 1:length(sigmonials)
    j = sigmonials(i);
    involved = find(model.monomtable(j,:));
    fractional = involved(model.monomtable(j,involved)~=fix(model.monomtable(j,involved)));
    for k = 1:length(fractional)
        model.lb(fractional(k)) = max([1e-9  model.lb(fractional(k))]);
    end
end

% The evaluation based operators can communicate the domain (although
% this really should be available via domain constraints anyway)
% The operator model can however also include range constraints, which we
% just as well might extract and add to the model
for i = 1:length(model.evalVariables)
    j = model.evalVariables(i);
    model.lb(j) = max([model.lb(j) model.evalMap{i}.properties.range(1)]);
    model.ub(j) = min([model.ub(j) model.evalMap{i}.properties.range(2)]);
    j = model.evalMap{i}.variableIndex;
    model.lb(j) = max([model.lb(j) repmat(model.evalMap{i}.properties.domain(1),length(j),1)],[],2);
    model.ub(j) = min([model.ub(j) repmat(model.evalMap{i}.properties.domain(2),length(j),1)],[],2);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function varargout = pow2(varargin)
%POW2 (overloaded)

switch class(varargin{1})

    case 'double' % What is the numerical value of this argument (needed for displays etc)
        % SHOULD NEVER HAPPEN, THIS SHOULD BE CAUGHT BY BUILT-IN
        error('Overloaded SDPVAR/NORM CALLED WITH DOUBLE. Report error')

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.
        if length(varargin{1}) == 1
            varargout{1} = yalmip('addEvalVariable',mfilename,varargin{1});
        else
            y = [];
            for i = 1:length(varargin{1})
                y = [y;yalmip('addEvalVariable',mfilename,extsubsref(varargin{1},i))];
            end
            varargout{1} = y;
        end

    case 'char' % YALMIP sends 'model' when it wants the epigraph or hypograph
        switch varargin{1}
            case 'graph'
                t = varargin{2};
                X = varargin{3};

                % This is different from so called extended operators
                % Just do it!
                F = SetupEvaluationVariable(varargin{:});

                % Now add your own code, such as domain constraints.
                % Exponential does not need any domain constraint.

                % Let YALMIP know about convexity etc
                varargout{1} = F;
                varargout{2} = struct('convexity','convex','monotonicity','increasing','definiteness','positive');
                varargout{3} = X;

            case 'milp'
                    varargout{1} = [];
                    varargout{2} = [];
                    varargout{3} = [];
            otherwise
                error('SDPVAR/EXP called with CHAR argument?');
        end
    otherwise
        error('SDPVAR/EXP called with CHAR argument?');
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [C,A,e,Av] = corrStaticFC (data,T,options,use_pca)
% When running the HMM on fMRI data (with order=0, i.e. a Gaussian
% distribution with mean and covariance),
% it can happen that states are assigned to entire subjects
% with not much switching between states. This happens typically when the
% static functional connectivity (FC) is so different that states specialise
% into specific subjects, explaining these grand patterns with
% no room to capture any dynamic FC.
% This function computes the (subjects by subjects) matrix of static FC
% similarities (measured in terms of correlation) between each pair of
% subjects. If the obtained values are too low, then covtype='sharedfull'
% has a higher chance to do a good job.
%
% If options.embeddedlags is specified, then it compares the time-embedded
% covariance matrices in a similar fashion; see Vidaurre et al. (2018) Nature
% Communications.
%
% If the use_pca parameter is set to 1, then C will be
% (subjects by subjects by no. of PCA components), with one matrix of
% similarities per number of PCA components; e.g. C(:,:,3) is the
% similarity matrix when using three PC; C(:,:,4) is the
% similarity matrix when using four PCs, etc. Finally, C(:,:,end) is the
% similarity marix in the original space (no PCA).
% C(:,:,1) and C(:,:,2) are set to NaN.
% In that case, A is the PCA mixing matrix, and e the explained variance
% for each number of components; e.g. e(3) is the explained variance for 3
% PCA components.
%
% Author: Diego Vidaurre (2017)

if nargin<3, options = struct(); end
if nargin<4, use_pca = 0; end
%if nargin<5, toPlot = 1; end

if xor(iscell(data),iscell(T)), error('X and T must be cells, either both or none of them.'); end

if iscell(T)
    if size(T,1)==1, T = T'; end
    for i = 1:length(T)
        if size(T{i},1)==1, T{i} = T{i}'; end
    end
    if size(data,1)==1, data = data'; end
end
N = length(T);
Av = [];

options.K = 2; options.order = 0;
options.BIGNbatch = 2; % to avoid checkoptions complaining
options = checkspelling(options);
options = checkoptions(options,data,T,0);
do_varimax = options.varimax;

if length(options.pca) > 1 || (options.pca > 0 && options.pca ~= 1)
    warning('options.pca will be ignored; use use_pca instead')
end

if use_pca
    [A,~,e] = highdim_pca(data,T,[],...
        options.embeddedlags,options.standardise,...
        options.onpower,0,options.detrend,...
        options.filter,options.leakagecorr,options.Fs);
    npca = find(e>0.95,1);
    if do_varimax
        Av = zeros(size(A,1),npca,npca);
        Av(:,1,1) = A(:,1);
        for d = 2:npca
            try
                Av(:,1:d,d) = rotatefactors(A(:,1:d));
            catch
                disp(['Varimax cannot deal with more components; npca set to ' num2str(d)])
                npca = d-1;
                break;
            end
        end

    end
    if isfield(options,'A'), warning('options.A will be ignored; use use_pca instead'); end
else
    A = []; e = [];
end

if use_pca
    FC = cell(npca,1);
end

for j = 1:N
    if iscell(data)
        if ischar(data{j})
            dat = load(data{j});
            if isfield(dat,'X')
                X = dat.X;
            else
                X = getfield(dat,char(fieldnames(dat)));
            end
        else
            X = data{j};
        end
        Tj = T{j};
    else
        ind = (1:T(j)) + sum(T(1:j-1));
        if isstruct(data)
            X = data.X(ind,:);
        else
            X = data(ind,:);
        end
        Tj = T(j);
    end
    % Standardise data and control for ackward trials
    X = standardisedata(X,Tj,options.standardise);
    % Filtering
    if ~isempty(options.filter)
       X = filterdata(X,Tj,options.Fs,options.filter);
    end
    % Detrend data
    if options.detrend
       X = detrenddata(X,Tj);
    end
    % Leakage correction
    if options.leakagecorr ~= 0
        X = leakcorr(X,Tj,options.leakagecorr);
    end
    % Hilbert envelope
    if options.onpower
       X = rawsignal2power(X,Tj);
    end
    % Embedding
    if length(options.embeddedlags) > 1
        [X,Tj] = embeddata(X,Tj,options.embeddedlags);
    end
    if use_pca
        X_nopca = X;
        ndim = size(X,2);
        for d = 2:npca
            if j==1, FC{d} =  zeros(d*(d-1)/2,N); end
            if do_varimax
                Ad = Av(:,1:d,d);
            else
                Ad = A(:,1:d);
            end
            X = X_nopca * Ad; % * Ad';
            X = X + 1e-8 * randn(size(X)); % add some noise to avoid ill-conditioning
            % Downsampling
            if options.downsample > 0
                X = downsampledata(X,Tj,options.downsample,options.Fs);
            end
            c = corr(X);
            FC{d}(:,j) = c(triu(true(size(X,2)),1))';
        end
        X = X_nopca;
        if options.downsample > 0
            X = downsampledata(X,Tj,options.downsample,options.Fs);
        end
        c = corr(X);
    else
        % PCA transform (precomputed)
        if isfield(options,'A') && ~isempty(options.A)
            if do_varimax
                A = rotatefactors(A);
            end
            X = X * A;
        end
        % Downsampling
        if options.downsample > 0
            X = downsampledata(X,Tj,options.downsample,options.Fs);
        end
        c = corr(X);
        if j==1
            ndim = size(X,2);
            FC = zeros(ndim*(ndim-1)/2,N);
        end
        FC(:,j) =  c(triu(true(ndim),1))';
    end
end

if use_pca
    C = NaN(N,N,npca);
    for d = 3:npca
        c = corr(FC{d});
        c(eye(N)==1) = Inf;
        C(:,:,d) = c;
    end
else
    C = corr(FC);
    C(eye(N)==1) = Inf;
end

% if toPlot
%     imagesc(C,[-1 1]); colorbar
%     set(gca,'FontSize',18)
%     colormap('jet');
%     grotc=colormap;  grotc(end,:)=[.8 .8 .8];  colormap(grotc);
% end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function varargout = inthist(varargin)
% VL_INTHIST  Calculate Integral Histogram
%   INTHIST = VL_INTHIST(LABELS) computes the integral histogram for
%   the label map(s) LABELS.
%
%   LABELS is a [M,N,L] UINT32 array containing L label maps (one for
%   each MxN layer). Each label map associates to each of the MxN
%   pixels one of K labels (a value of 0 denotes no association) with
%   unitary mass.
%
%   INTHIST is the [M,N,K] integral histogram obtained by integrating
%   the label maps (an integral histogram is just an array of K
%   integral images, one for each of the K labels). All layers of
%   labels are accumulated to the same integral histogram.
%
%   The total number of labels K is obtained as the maximum value of
%   LABELS. VL_INTHIST(..., 'NUMLABELS', K) specifies the number of
%   labels explicitly.
%
%   VL_INTHIST(..., 'MASS', MASSES) specifies a mass MASSES for each
%   entry of LABELS. MASSES can be either of class UINT32 or DOUBLE,
%   and the class of INTHIST varies accordingly.
%
%   See also: VL_IMINTEGRAL(), VL_SAMPLEINTHIST(), VL_HELP().
[varargout{1:nargout}] = vl_inthist(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
% TANDEMO   This file generates the tan(x) - x example from Chapter 1.
%           See Figures 1.2 and 1.3.
%
% C. T. Kelley, October 22, 2002.
%
x0=4.5; tol=1.d-20;
%
% Solve the problem three times.
%
[x,hist]=newtsol(x0,'ftan',tol,tol,1);
[x,histc]=chordsol(x0,'ftan',tol,tol);
[x,hists]=secant(x0,'ftan',tol,tol);
%
% Iteration history for Newton.
% I use handle graphics to get the axis the way I
% want them. Try commenting out "set(p1,'XTick',[0 1 2 3 4 5]);"
% and see what the difference is.
%
maxit=6;
figure(1);
p1=subplot(1,1,1);
semilogy(hist(1:maxit,1),abs(hist(1:maxit,2)))
set(p1,'XTick',[0 1 2 3 4 5],'FontSize',14);
axis([0 5 1.d-16 1]);
xlabel('Nonlinear iterations'); ylabel('Absolute Nonlinear Residual');

%
% Plot 15 iterations for all three methods.
%
figure(2);
maxit=15;
p2=subplot(1,1,1);
semilogy(hist(1:maxit,1),abs(hist(1:maxit,2)),'-',...
histc(1:maxit,1),abs(histc(1:maxit,2)),'--',...
hists(1:maxit,1),abs(hists(1:maxit,2)),'-.');
set(p2,'XTick',[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14],'FontSize',14);
axis([0 14 1.d-16 1]);
legend('Newton','Chord','Secant');
xlabel('Nonlinear iterations'); ylabel('Absolute Nonlinear Residual');

repo: Aero-matlab


text:
function mr = mrReadVAnat(pth)
% Read a mrGray vAnatomy .dat file into an mr struct.
%
% mr = mrReadVAnat([path to file])
%
% Loads the vAnatomy.dat file specified by fileName (full path!)
% into the [rows,cols,planes] image cube 'img'.
%
% If fileName is omitted, a get file dialog appears.
%
% RETURNS:
%   * img is the [rows,cols,planes] intensity array
%   * mmPerPix is the voxel size (in mm/pixel units)
%   * fileName is the full-path to the vAnatomy.dat file. (If
%     you pass fileName in, you obviously don't need this. But
%     it may be useful when the user selects the file.)
%
% ras, 06/30/05: imported into mrVista 2.0 Test repository.

% initialize mr struct
[p f ext] = fileparts(pth);
mr = mrCreateEmpty;
mr.format = 'vanat';
mr.pth = pth;
mr.name = f;

[mr.data, mr.voxelSize, mr.dims] = readFileVAnat(pth);
mr.data = double(mr.data);
mr.spaces = mrStandardSpaces(mr);

% based on the usual format of vAnatomy files, we can
% label the directions in the pixel space with a good guess:
mr.spaces(1).dirLabels = {'S <--> I'  'A <--> P'  'L <--> R'};
mr.spaces(1).sliceLabels =  {'Axial' 'Coronal' 'Sagittal'};
mr.spaces(2).dirLabels = {'S <--> I'  'A <--> P'  'L <--> R'};
mr.spaces(2).sliceLabels =  {'Axial' 'Coronal' 'Sagittal'};

mr.spaces(end+1) = mr.spaces(2); % add I|P|R space
mr.spaces(end).name = 'I|P|R';

mr.spaces(end+1) = mr.spaces(end);
mr.spaces(end).name = 'R|A|S';
mr.spaces(end).dirLabels =  {'L <--> R' 'P <--> A' 'I <--> S'};
mr.spaces(end).sliceLabels =  {'Sagittal' 'Coronal' 'Axial'};
mr.spaces(end).xform = ipr2ras(mr.spaces(4).xform,mr.dims,mr.voxelSize);

mr.dimUnits = {'mm' 'mm' 'mm' 'sec'};
mr.dataUnits = 'Scaled T1 Intensity';

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%DEMOIMAGEGRAPHS  One-line description here, please.
%
%   output = demoImageGraphs(input)
%
%   Example
%   demoImageGraphs
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
% Created: 2022-04-01,    using Matlab 9.12.0.1884302 (R2022a)
% Copyright 2022 INRAE.

% create a simple 2D binary image

lx = -10:10;
ly = -10:10;
[x, y] = meshgrid(lx, ly);
img = hypot(x, y) < 8.5;

figure;
imshow(img, 'InitialMagnification', 1000);
print(gcf, 'sample_image.png', '-dpng');

% create adjacency graph

[v, e] = imageGraph(img);
figure; imshow(ones(size(img)), 'InitialMagnification', 1000);
hold on; drawGraph(v, e);
print(gcf, 'imageGraph_2d.png', '-dpng');

% create boundary graph

[v, e] = imageBoundaryGraph(img);
figure; imshow(ones(size(img)), 'InitialMagnification', 1000);
hold on; drawGraphEdges(v, e, 'Color', 'b', 'LineWidth', 2);
print(gcf, 'imageBoundaryGraph_2d.png', '-dpng');

% save into figure
print(gcf, 'image_graphs_2d.png', '-dpng');

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [h]=m_triplot(long,lat,tri)
%  M_CONTOURF Adds filled contours to a map
%    M_CONTOUR(LONG,LAT,DATA,...)
%
global MAP_PROJECTION

% Have to have initialized a map first

if isempty(MAP_PROJECTION)
  disp('No Map Projection initialized - call M_PROJ first!');
  return;
end
bcol=[.8,.9,1];
[X,Y]=m_ll2xy(long,lat,'clip','on');

hold on; trimesh(tri,X,Y,0*X,'facecolor',bcol,'edgecolor','k');
view(2)
%triplot(tri,X,Y,'k','facecolor',bcol);

%m_coast('patch','red') ;

end

author: CHLNDDEV


text:
%COMPUTEEXACTMARGINALSBP Runs exact inference and returns the marginals
%over all the variables (if isMax == 0) or the max-marginals (if isMax == 1).
%
%   M = COMPUTEEXACTMARGINALSBP(F, E, isMax) takes a list of factors F,
%   evidence E, and a flag isMax, runs exact inference and returns the
%   final marginals for the variables in the network. If isMax is 1, then
%   it runs exact MAP inference, otherwise exact inference (sum-prod).
%   It returns an array of size equal to the number of variables in the
%   network where M(i) represents the ith variable and M(i).val represents
%   the marginals of the ith variable.
%
% Copyright (C) Daphne Koller, Stanford University, 2012

function M = ComputeExactMarginalsBP(F, E, isMax)

% initialization
% you should set it to the correct value in your code
M = [];
vars = unique([F.var]);

N = length(vars);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR CODE HERE
%
% Implement Exact and MAP Inference.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CliqTree = CreateCliqueTree(F,E);
P = CliqueTreeCalibrate(CliqTree,isMax);
M = repmat(struct('var',[],'card',[],'val',[]),1,N);
for i = 1:N
	for j = 1:length(CliqTree.cliqueList)
		inter = intersect(vars(i),CliqTree.cliqueList(j).var);
		if length(inter)==1
			V = setdiff(CliqTree.cliqueList(j).var,vars(i));
			if isMax == 0
				M(i) = FactorMarginalization(P.cliqueList(j),V);
			else
				M(i) = FactorMaxMarginalization(P.cliqueList(j),V);
			end
			break;
		end
	end
end
if isMax==0
for i=1:N
	M(i).val = M(i).val/sum(M(i).val);
end
end
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [fx,tx,pv,fv]=fxpefac(s,fs,tinc,m,pp)
%FXPEFAC PEFAC pitch tracker [FX,TT,PV,FV]=(S,FS,TINC,M,PP)
%
% Input:   s(ns)      Speech signal
%          fs         Sample frequency (Hz)
%          tinc       Time increment between frames (s) [0.01]
%                     or [start increment end]
%          m          mode
%                     'g' plot graph showing waveform and pitch
%                     'G' plot spectrogram with superimposed pitch
%                     'x' use external files for algorithm parameter
%                         initialization: fxpefac_g and fxpefac_w
%          pp         structure containing algorithm parameters
%
% Outputs: fx(nframe)     Estimated pitch (Hz)
%          tx(nframe)     Time at the centre of each frame (seconds).
%          pv(nframe)     Probability of the frame of being voiced
%          fv             structure containing feature vectors
%                           fv.vuvfea(nframe,2) = voiced/unvoiced GMM features

% References
%  [1]  S.Gonzalez and M. Brookes,
%       A pitch estimation filter robust to high levels of noise (PEFAC), Proc EUSIPCO,Aug 2011.

% Bugs/Suggestions
% (1) do long files in chunks
% (2) option of n-best DP

%	   Copyright (C) Sira Gonzalez and Mike Brookes 2011
%      Version: $Id: fxpefac.m,v 1.2 2011/07/27 07:22:07 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

persistent w_u m_u v_u w_v m_v v_v dpwtdef
% initialize persistent variables
if ~numel(w_u)

    % voiced/unvoiced decision based on 2-element feature vector
    % (a) mean power of the frame's log-freq spectrum (normalized so its short-term average is LTASS)
    % (b) sum of the power in the first three peaks
    %===== VUV
    if nargin>3 && any(m=='x')
        fxpefac_g;    % read in GMM parameters
        fxpefac_w;     % read in Weights parameters
    else
        w_u=[0.2123723 0.207788 0.2701817 0.1293616 0.04741722 0.1328791 ]';

        m_u=[0.2220388 0.4067706 ;
            0.04567656 0.4016914 ;
            0.8415278 0.3192158 ;
            0.2194808 0.1910079 ;
            1.6347 0.5819833 ;
            1.181519 0.6996485 ];

        v_u=reshape([0.01413822 0.003357913 0.003357913 0.01786169 ;
            0.0009377269 0.0006220489 0.0006220489 0.03422057 ;
            0.1233703 0.004299293 0.004299293 0.007660504 ;
            0.01779449 0.002078821 0.002078821 0.001605052 ;
            1.110173 0.00718649 0.00718649 0.005734435 ;
            0.5477135 -0.00182316 -0.00182316 0.05659796 ]',[2 2 6]);

        w_v=[0.07758689 0.2109879 0.1856225 0.06853158 0.2701563 0.1871148 ]';

        m_v=[1.208656 0.3365564 ;
            1.216643 0.5971916 ;
            4.08585 1.240948 ;
            8.322102 1.349939 ;
            1.734108 1.168643 ;
            0.5107205 0.940308 ];

        v_v=reshape([0.06181574 0.002950501 0.002950501 0.004528442 ;
            0.2946077 0.01433284 0.01433284 0.02684239 ;
            2.508473 -0.03310555 -0.03310555 0.1098579 ;
            14.17252 -0.09009174 -0.09009174 0.07989255 ;
            0.5834894 -0.07854027 -0.07854027 0.1108958 ;
            0.05978017 0.005528601 0.005528601 0.1309329 ]',[2 2 6]);
    end
    %===== PDP
    %     dfm = -0.4238; % df mean
    %     dfv = 3.8968; % df variance (although treated as std dev here)
    %     delta = 0.15;
    %     dflpso=[dfm 0.5/(log(10)*dfv^2) -log(2*delta/(dfv*sqrt(2*pi)))/log(10)]; % scale factor & offset for df pdf
    %     dpwtdef=[1.0000, 0.8250, 1.3064, 1.9863]; % default DP weights
    dpwtdef=[1.0000, 0.8250, 0.01868, 0.006773, 98.9, -0.4238]; % default DP weights
    %===== END

end

% Algorithm parameter defaults

p.fstep=5;              % frequency resolution of initial spectrogram (Hz)
p.fmax=4000;            % maximum frequency of initial spectrogram (Hz)
p.fres = 20;            % bandwidth of initial spectrogram (Hz)
p.fbanklo = 40;         % low frequency limit of log filterbank (Hz)
p.mpsmooth = 201;       % width of smoothing filter for mean power
% p.maxtranf = 1000;      % maximum value of tranf cost term
p.shortut = 7;          % max utterance length to average power of entire utterance
p.pefact = 1.5;         % shape factor in PEFAC filter
p.numopt = 3;           % number of possible frequencies per frame
p.flim = [60 400];      % range of feasible fundamental frequencies (Hz)
p.w = dpwtdef;          % DP weights
% p.rampk = 1.1;          % constant for relative-amplitude cost term
% p.rampcz = 100;         % relative amplitude cost for missing peak
p.tmf = 2;              % median frequency smoothing interval (s)
p.tinc = 0.01;          % default frame increment (s)

% update parameters from pp argument

if nargin>=5 && isstruct(pp)
    fnq=fieldnames(pp);
    for i=1:length(fnq)
        if isfield(p,fnq{i})
            p.(fnq{i})=pp.(fnq{i});
        end
    end
end

% Sort out input arguments
if nargin>=3  && numel(tinc)>0
    p.tinc = tinc;   % 0.01 s between consecutive time frames
end
if nargin<4
    m='';
end

% Spectrogram of the mixture
fmin = 0; fstep = p.fstep; fmax = p.fmax;
fres = p.fres;  % Frequency resolution (Hz)
[tx,f,MIX]=spgrambw(s,fs,fres,[fmin fstep fmax],[],p.tinc);
nframes=length(tx);
txinc=tx(2)-tx(1);  % actual frame increment
%  ==== we could combine spgrambw and filtbankm into a single call to spgrambw or use fft directly ====
% Log-frequency scale
[trans,cf]=filtbankm(length(f),2*length(f)-1,2*f(end),p.fbanklo,f(end),'usl');
O = MIX*trans'; % Original spectrum in Log-frequency scale

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Amplitude Compression

% Calculate alpha based on LTASS ratios
ltass = stdspectrum(6,'p',cf);
auxf = [cf(1),(cf(1:end-1)+cf(2:end))./2,cf(end)];
ltass = ltass.*diff(auxf);                  % weight by bin width

% estimated ltass
O = O.*repmat(diff(auxf),nframes,1);     % weight spectrum by bin width

if tx(end)<p.shortut                        % if it is a short utterance
    eltass = mean(O,1);                     % mean power per each frequency band
    eltass = smooth(eltass,p.mpsmooth);     % smooth in log frequency
    eltass= eltass(:).';                    % force a row vector

    % Same mean power per frame as ltass
    cte = mean(ltass)/mean(eltass);
    eltass = eltass.*cte;                   % normalize to have the same mean as LTASS
    O = O.*cte;

    % Linear AC
    alpha = (ltass)./(eltass);
    alpha = alpha(:).';
    alpha = repmat(alpha,nframes,1);
    O = O.*alpha;                           % force O to have an average LTASS spectrum
    % ==== should perhaps exclude the silent portions ***
else                                        % long utterance
    tsmo = 2; % time smoothing over 1 sec
    stt = round(tsmo/txinc);
    filttime = [ones(stt,1); zeros(stt-1,1)];
    filtfreq = ones(1,p.mpsmooth);
    eltass = imfilter(O,filttime);
    eltass = imfilter(eltass,filtfreq);     % filter in time and log frequency

    % Same mean power per frame than ltass
    cte = repmat(mean(ltass),nframes,1)./mean(eltass,2);
    eltass = eltass.*repmat(cte,1,length(cf));
    O = O.*repmat(cte,1,length(cf));

    % Linear AC
    alpha = repmat(ltass,nframes,1)./(eltass);
    O = O.*alpha;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the filter to detect the harmonics
ini = find(cf>2*cf(1));
sca = cf/cf(ini(1)); % normalize bin frequencies to start at approximately 0.5
sca = sca(sca<10.5 & sca>0.5);  % restrict to 0.5 - 10.5 times fundamental
filh = -log10(p.pefact-cos(2*pi*sca));
filh = filh-mean(filh);  % force filter to be zero mean
posit = find(sca>=1);  % ==== this should just equal ini(1) ====
if ~mod(length(posit),2)
    filh = [filh 0];  % force to be an odd length after central tap
end
negat = find(sca<1);
numz = length(posit)-1-length(negat);
filh = filh./max(filh);
filh = [zeros(1,numz) filh];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Filter the log-frequency scaled spectrogram
B = imfilter(O,filh);  % ==== no good reason to use imfilter here ====

% Feasible frequency range
numopt = p.numopt; % Number of possible fundamental frequencies per frame
flim = p.flim;
pfreq = find(cf>flim(1) & cf<flim(2));
ff = zeros(nframes,numopt);
amp = zeros(nframes,numopt);
for i=1:nframes
    [pos,peak]=findpeaks(B(i,pfreq),[],5/(cf(pfreq(2))-cf(pfreq(1)))); % ==== calculate some out of loop ====
    if numel(pos)
        [peak,ind]=sort(peak,'descend');
        pos = pos(ind);                     % indices of peaks in the B array
        posff = cf(pfreq(pos));             % frequencies of peaks
        fin = min(numopt,length(posff));
        ff(i,1:fin)=posff(1:fin);           % save both frequency and amplitudes
        amp(i,1:fin)=peak(1:fin);
        %     else
        %         ff(i,:)=0;          % ==== unnecessary since they start as zeros ====
        %         amp(i,:)=0;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Probabilitly of the frame of being voiced

% voiced/unvoiced decision based on 2-element feature vector
% (a) mean power of the frame's log-freq spectrum (normalized so its short-term average is LTASS)
% (b) sum of the power in the first three peaks

pow = mean(O,2)*1e-6;
vuvfea = [pow sum(amp,2)./(pow*1e9)];

pru=gaussmixp(vuvfea,m_u,v_u,w_u);  % Probability of being unvoiced
prv=gaussmixp(vuvfea,m_v,v_v,w_v);  % Probability of being voiced

% pru = exp(pru);
% prv=exp(prv);    % Linear probability
% pv = prv./(prv+pru); % ==== better to write pv=(1+exp(pru-prv)).^(-1) ====
pv=(1+exp(pru-prv)).^(-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Dynamic programming

% w(1): relative amp, voiced local cost
% w(2): median pitch deviation cost
% w(3): df cost weight
% w(4): max df cost
% w(5): relative amp cost for missing peaks (very high)
% w(6): df mean

w = p.w;

% Relative amplitude
camp = -amp./repmat(max(amp,[],2),1,numopt);  % relative amplitude used as cost
camp(amp==0)=w(5); % If no frequency found

% Time interval for the median frequency
tmf = p.tmf; % in sec
inmf = round(tmf/txinc);

%--------------------------------------------------------------------------
% FORWARDS
% Initialize values
cost = zeros(nframes,numopt);
prev = zeros(nframes,numopt);
medfx = zeros(nframes,1);
dffact=2/txinc;

% First time frame
% cost(1,:) = w(1)*ramp(1,:);
cost(1,:) = w(1)*camp(1,:);  % only one cost term for first frame
fpos = ff(1:min(inmf,end),1);
mf=median(fpos(pv(1:min(inmf,end))>0.6));   % calculate median frequency of first 2 seconds
if isnan(mf)
    mf=median(fpos(pv(1:min(inmf,end))>0.5));
    if isnan(mf)
        mf=median(fpos(pv(1:min(inmf,end))>0.4));
        if isnan(mf)
            mf=median(fpos(pv(1:min(inmf,end))>0.3)); % ==== clumsy way of ensuring that we take the best frames ====
            if isnan(mf)
                mf=0;
            end
        end
    end
end
medfx(1)=mf;

for i=2:nframes              % main dynamic programming loop
    if i>inmf
        fpos = ff(i-inmf:i,1);  % fpos is the highest peak in each frame
        mf=median(fpos(pv(1:inmf)>0.6));  % find median frequency over past 2 seconds
        if isnan(mf)
            mf=median(fpos(pv(1:inmf)>0.5));
            if isnan(mf)
                mf=median(fpos(pv(1:inmf)>0.4));
                if isnan(mf)
                    mf=median(fpos(pv(1:inmf)>0.3));% ==== clumsy way of ensuring that we take the best frames ====
                    if isnan(mf)
                        mf=0;
                    end
                end
            end
        end
    end
    medfx(i)=mf;
    % Frequency difference between candidates and cost
    df = dffact*(repmat(ff(i,:).',1,numopt) - repmat(ff(i-1,:),numopt,1))./(repmat(ff(i,:).',1,numopt) + repmat(ff(i-1,:),numopt,1));
    costdf=w(3)*min((df-w(6)).^2,w(4));

    % Cost related to the median pitch
    if mf==0                                   % this test was inverted in the original version
        costf = zeros(1,numopt);
    else
        costf = abs(ff(i,:) - mf)./mf;
    end
    [cost(i,:),prev(i,:)]=min(costdf + repmat(cost(i-1,:),numopt,1),[],2); % ==== should we allow the possibility of skipping frames ? ====
    cost(i,:)=cost(i,:)+w(2)*costf + w(1)*camp(i,:);  % add on costs that are independent of previous path

end

% Traceback

fx=zeros(nframes,1);
best = zeros(nframes,1);

nose=find(cost(end,:)==min(cost(end,:))); % ==== bad method (dangerous) ===
best(end)=nose(1);
% ff = [ff zeros(nframes,1)];  % not clear why this was here
fx(end)=ff(end,best(end));
for i=nframes:-1:2
    best(i-1)=prev(i,best(i));
    fx(i-1)=ff(i-1,best(i-1));
end

if nargout>=4
    fv.vuvfea=vuvfea;  % voiced-unvoiced features
    fv.best=best;  % selected path
    fv.ff=ff;  % pitch candidates
    fv.amp=amp;  % pitch candidate amplitudes
    fv.medfx=medfx;  % median pitch
    fv.w=w;  % DP weights
    fv.dffact=dffact;  % df scale factor
end

if ~nargout || any(m=='g') || any(m=='G')
    nax=0;  % number of axes sets to link
    msk=pv>0.5; % find voiced frames as a mask
    fxg=fx;
    fxg(~msk)=NaN; % allow only good frames
    fxb=fx;
    fxb(msk)=NaN; % allow only bad frames
    if any(m=='G') || ~nargout && ~any(m=='g')
        clf;
        spgrambw(s,fs,'ilcwpf'); % draw spectrogram with log axes
        hold on
        plot(tx,log10(fxg),'-b',tx,log10(fxb),'-r'); % fx track
        yy=get(gca,'ylim');
        plot(tx,yy(1)+yy*[-1;1]*(0.02+0.05*pv),'-k'); % P(V) track
        hold off
        nax=nax+1;
        axh(nax)=gca;
        if any(m=='g')
            figure;   % need a new figure if plotting two graphs
        end
    end
    if any(m=='g')
        ns=length(s);
        [tsr,ix]=sort([(1:ns)/fs 0.5*(tx(1:end-1)+tx(2:end))']); % intermingle speech and frame boundaries
        jx(ix)=1:length(ix); % create inverse index
        sp2fr=jx(1:ns)-(0:ns-1);  % speech sample to frame number
        spmsk=msk(sp2fr);   % speech sample voiced mask
        sg=s;
        sg(~spmsk)=NaN;   % good speech samples only
        sb=s;
        sb(spmsk)=NaN;    % bad speech samples only
        clf;
        subplot(5,1,1);
        plot(tx,pv,'-b',(1:ns)/fs,0.5*mod(cumsum(fx(sp2fr)/fs),1)-0.6,'-b');
        nax=nax+1;
        axh(nax)=gca;
        ylabel('\phi(t), P(V)');
        set(gca,'ylim',[-0.65 1.05]);
        subplot(5,1,2:3);
        plot((1:ns)/fs,sg,'-b',(1:ns)/fs,sb,'-r');
        nax=nax+1;
        axh(nax)=gca;
        subplot(5,1,4:5);
        plot(tx,fxg,'-b',tx,fxb,'-r');
        ylabel('Pitch (Hz)');
        %         semilogy(tx,fxg,'-b',tx,fxb,'-r');
        %         ylabel(['Pitch (' yticksi 'Hz)']);
        set(gca,'ylim',[min(fxg)-30 max(fxg)+30]);
        nax=nax+1;
        axh(nax)=gca;
    end
    if nax>1
        linkaxes(axh,'x');
    end
end

function y=smooth(x,n)
nx=length(x);
c=cumsum(x);
y=[c(1:2:n)./(1:2:n) (c(n+1:end)-c(1:end-n))/n (c(end)-c(end-n+2:2:end-1))./(n-2:-2:1)];

author: decouples


text:
function saveNEVTetrodes(NEVFullFilename)

% saveNEVTetrodes
%
% Opens saves a new NEV file and splits it into smaller NEV files
% containing only the tetrodes according to the associated CCF file. The
% CCF file should have the same name as the data file.
%
% NEVFullFilename:  The full path to the NEV to be opened.
%                   DEFAULT: If the filename is not provided, the user will
%                   be prompted to select a file.
%
% Use saveNEVTetrodes(NEVFullFilename)
%
%   Example: saveNEVTetrodes('c:\datafolder\datafile.nev');
%
%   The function will open datafile.nev and datafile.ccf and based on the
%   tetrode information saved in datafile.ccf, it will split datafile.nev
%   into smaller chunks that will only contain channels associated to that
%   particular tetrode.
%
%   Kian Torab
%   ktorab@blackrockmicro.com
%   Blackrock Microsystems
%   Version 1.0.0.0
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Version History
%
% 1.0.0.0:
%   - Initial release.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Opening the file
if ~exist('NEVFullFilename', 'var')
    [dataFilename dataFolder] = getFile('*.nev');
    NEVFullFilename = [dataFolder dataFilename];
elseif exist(NEVFullFilename, 'file') ~= 2
        disp('The NEV file name does not exist.');
        return;
else
    [filepath filename fileext] = fileparts(NEVFullFilename);
    dataFolder = [filepath '/'];
    dataFilename = [filename fileext];
end

%% Openning the associated CCF file
ccfFullFilename = [dataFolder dataFilename(1:end-3) 'ccf'];
if exist(ccfFullFilename, 'file') ~= 2 % for 2.x file type
    ccfFullFilename = [dataFolder dataFilename(1:end-8) '.ccf'];
    if exist(ccfFullFilename, 'file') ~= 2 % for TOC file type
        disp('Cannot find the associated CCF file.');
        disp('This function requires the CCF file used during the recording.');
        disp('The CCF must have the same name as the original recorded file.');
        return;
    end
end
ccf = openCCF(ccfFullFilename);

% Calculating the number of NTrode groups in the file
if isfield(ccf, 'NTrodeInfo')
    if isfield(ccf.NTrodeInfo, 'NTrodeMembers')
        numberOfNTrodeGroups = size(ccf.NTrodeInfo.NTrodeMembers,2);
    else
        disp('This data file does not contain any tetrodes.');
        return;
    end
else
    disp('There is an error in the CCF file.');
    return;
end

% Splitting the data file according to the information saved in the
% associated CCF file.
for idx = 1:numberOfNTrodeGroups
    fprintf('Saving tetrode %d containing channels %s.\n', idx, num2str(ccf.NTrodeInfo.NTrodeMembers{idx}));
    saveNEVSubSpikes(ccf.NTrodeInfo.NTrodeMembers{idx}, NEVFullFilename, 'tet');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text: %% S2FunHandle

author: mtex-toolbox

repo: mtex

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


text:
function [vol] = volume_intrinsic(l)
  % VOLUME Compute volumes of tets defined intrinsically by edge lengths l
  %
  % v = volume(l)
  %
  % Inputs:
  %   l  #T by 6 list of tetrahedra side lengths of edges opposite *face* pairs
  %     [23 31 12 41 42 43]
  % Ouputs:
  %   vol  #T list of tet volumes (always positive)
  %

  % http://en.wikipedia.org/wiki/Heron%27s_formula#Heron-type_formula_for_the_volume_of_a_tetrahedron

  % U, V, W, u, v, w are lengths of edges of the tetrahedron (first three form
  % a triangle; u opposite to U and so on)
  u = l(:,1); v = l(:,2); w = l(:,3);
  U = l(:,4); V = l(:,5); W = l(:,6);
  X = (w - U + v).*(U + v + w);
  x = (U - v + w).*(v - w + U);
  Y = (u - V + w).*(V + w + u);
  y = (V - w + u).*(w - u + V);
  Z = (v - W + u).*(W + u + v);
  z = (W - u + v).*(u - v + W);
  a = sqrt(x.*Y.*Z);
  b = sqrt(y.*Z.*X);
  c = sqrt(z.*X.*Y);
  d = sqrt(x.*y.*z);
  vol = sqrt( ...
    (-a + b + c + d).* ...
    ( a - b + c + d).* ...
    ( a + b - c + d).* ...
    ( a + b + c - d))./ ...
    (192.*u.*v.*w);

end

author: alecjacobson

repo: gptoolbox


text:
function varargout = solveTSP( cities, display)
% cities = solveTSP( cities, maxItt, display)
%
% cities - An Nx2 matrix containing cartesian coordinates of the "cities"
% beeing visited. The initial trail is assumed from the first city to the
% scond and so on...
%
% display - bolean flag decide if to display the progress of the program (slows the running time).
%                       default = false;
%
% maxIteration - maximum iterations for the program
%                                   default = 10,000
%
%
% [cities ind] = solveTSP( cities, display) returns the aranged cities and
% an index vector of the visiting order
%
% [cities ind totalDist] = solveTSP( cities, display)
% totalDist is the route total distance
%
% demo1:
% cities = solveTSP( rand(100,2), true );
%
% demo2:
% t = (0:999)' /1000;
% cities = [ t.^2.*cos( t*30 ) t.^2.*sin( t*30 ) ];
% [ans ind] = sort( rand(1000,1) );
% [cities ind] = solveTSP( cities(ind,:), true );

    if nargin < 2
        display = false;
    end

    siz = size(cities);
    if siz(2) ~= 2
        error( 'The program is expecting cities to be an Nx2 matix of cartesian coordinates' );
    end
    N = siz(1);

    order = (1:N)';     % initial cities visit order

    if display
        hFig = figure;
        hAx = gca;
        updateRate = ceil( N/50 );
    end

    itt = 1;
    maxItt = min(20*N,1e5);
    noChange = 0;

    while itt < maxItt  && noChange < N

        dist = calcDistVec( cities(order,:),1 );    % travel distance between the cities

        %% ----------- Displaying current route -----------------------
        if display && ~mod(itt,updateRate) && ishandle( hFig )
            hold(hAx,'off');
            plot( hAx, cities( order,1),cities( order,2),'r.' );
            hold( hAx,'on');
            plot( hAx, cities( order,1),cities( order,2) );
            str = {[ 'iteration: ' num2str( itt ) ] ;
                           [ 'total route: ' num2str( sum( dist) ) ] };
            title(  hAx,str );
            pause(0.02)
        end

        flip = mod( itt-1, N-3 )+2 ;

        untie = dist(1:end-flip) + dist(flip+1:end);    % the distance saved by untying a loop
        shufledDist = calcDistVec( cities( order,:),flip );
        connect =  shufledDist(1:end-1) + shufledDist( 2:end);  % the distance payed by connecting the loop (after flip)
        benifit = connect - untie;  % "what's the distance benifit from this loop fliping

        %% --------------- Finding the optimal flips (most benficial) ----------------
        localMin = imerode(benifit,ones(2*flip+1,1) );
        minimasInd = find( localMin == benifit);
        reqFlips = minimasInd( benifit(minimasInd) < -eps );

        %% -------- fliping all loops found worth fliping --------------------
        prevOrd = order;
        for n=1:numel( reqFlips )
            order( reqFlips(n) : reqFlips(n)+flip-1 ) = order( reqFlips(n) +flip-1: -1 :reqFlips(n) );
        end

        %% -------  counting how many iterations there was no improvement
        if isequal( order,prevOrd )
            noChange = noChange + 1;
        else
            noChange = 0;
        end

        itt = itt+1;

    end     % while itt < maxItt  && noChange < N

    output = {cities( order,:), order, sum( dist)};
    varargout = output(1:nargout);

function dist = calcDistVec( cord,offset )
% dist = calcDistVec( cord,offset )
% offset is the number of cities to calculate the distence between
% the distance for the first city is allway 0

    dist = zeros( size(cord,1)-offset+1,1 );
    temp = cord( 1:end-offset,:) - cord( offset+1:end,:);
    dist(2:end) = sqrt( sum(temp.^2,2) );

author: Sable

repo: mcbench-benchmarks


text:
function [flag, t, lambda, S] = ray_polygon_intersect(o,d,V,E)
% RAY_POLYGON_INTERSECT 2D Ray/polygon intersection
%
% [flag, t, lambda] = ray_polygon_intersect(o,d,V,E)
%
% Input:
%    o  2D vector ray origin.
%    d  2D vector ray direction.
%    V  #V by 2 list of vertex positions
%    E  #E by 2 list of edge indices
% Output:
%    flag  #E list of bools: (false) Reject, (true) Intersect.
%    t  #E list of distances from the ray origin.
%    lambda  #E list of parameter of hit between E(:,1) and E(:,2)
%    S  #E by #o list of signs
%

  epsilon = eps;%0.00001;

  assert(size(V,2) == 2);
  assert(size(E,2) == 2);
  % number of edges
  m = size(E,1);
  %assert(numel(d) == 2);

  p1 = V(E(:,1),:);
  p2 = V(E(:,2),:);
  % edge vectors from 1 to 2
  e12 = p2-p1;
  % perpendiculars of edge vectors
  pe12 = perp(e12);

  if numel(d) == 2 && numel(o) == 2
    % make direction a row vector
    d = reshape(d,1,2);
    %d = normalizerow(d);
    d = d ./ sqrt(sum(d.^2,2));
    assert(numel(o) == 2);
    % make origin a row vector
    o = reshape(o,1,2);
    % p1 minus o
    p1mo = plusrow(p1,-o);
    % pe12 dot d
    pe12dd = dotrow(pe12,d);
    % perp d
    pd = perp(d);
    % project to edges
    % http://objectmix.com/graphics/132701-ray-line-segment-intersection-2d.html
    t = dot2(pe12,p1mo)./pe12dd;
    lambda = dotrow(plusrow(-p1,o+d),pd)./dotrow(e12,pd);
    %lambda = dotrow(perp(d),p1mo)./pe12dd;
    flag = true(m,1);
    flag(lambda > 1) = 0;
    flag(lambda < 0) = 0;
    flag(t < 0) = 0;
  else
    d = normalizerow(d);
    pd = perp(d);
    d = permute(d,[3 2 1]);
    pd = permute(pd,[3 2 1]);
    o = permute(o,[3 2 1]);
    p1mo = p1-o;
    pe12dd = sum(pe12.*d,2);
    % project to edges
    t = sum(pe12.*p1mo,2)./pe12dd;
    lambda = sum((o+d-p1).*pd,2)./sum(e12.*pd,2);
    n = max(size(o,3),size(d,3));
    flag = true(m,n);
    flag(lambda > 1) = 0;
    flag(lambda < 0) = 0;
    flag(t < 0) = 0;
    t = squeeze(t);
    S = sign(squeeze(pe12dd));
  end

    function u = perp(v)
      % [x,y] = [y,-x]
      u = [v(:,2),-v(:,1)];
    end

    function r = dot2(a,b)
      % Optimizes r = dot(a,b,2), that is it computes dot products per row
      % Faster than dot if I know that I'm calling it correctly
      r = sum(a.*b,2);
    end

    function r = dotrow(a,b)
      % Computes dot product of rows in a against b.
      % optimizes r = dot(a,repmat(b,size(a,1),1),2)
      r = a(:,1).*b(1,1) + a(:,2).*b(1,2);
    end

    function r = plusrow(a,b)
      % Computes sum rows in a with b.
      % optimizes r = a + repmat(b,size(a,1),1)
      r = [a(:,1)+b(1,1) a(:,2)+b(1,2)];
    end

end


text:
function [dlnZ_dmu, dlnZ_dvs] = heavisideNoiseGradVals(noise, mu, varsigma, y)

% HEAVISIDENOISEGRADVALS Gradient wrt mu and varsigma of log-likelihood for heaviside noise model.

% IVM

D = size(mu, 2);
c = y./sqrt(varsigma);
u = zeros(size(c));
dlnZ_dmu = zeros(size(c));
for i = 1:D
  u(:, i) = c(:, i).*(mu(:, i) + noise.bias(i));
  dlnZ_dmu(:, i) = c(:, i).*ngaussian(u(:, i))...
            ./(cumGaussian(u(:, i))+ ...
               noise.eta/(1-2*noise.eta));
end
dlnZ_dvs = -.5*c.*u.*dlnZ_dmu;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef DTLZ7 < PROBLEM
% <multi/many> <real> <large/none> <expensive/none>
% Benchmark MOP proposed by Deb, Thiele, Laumanns, and Zitzler

%------------------------------- Reference --------------------------------
% K. Deb, L. Thiele, M. Laumanns, and E. Zitzler, Scalable test problems
% for evolutionary multiobjective optimization, Evolutionary multiobjective
% Optimization. Theoretical Advances and Applications, 2005, 105-145.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            if isempty(obj.M); obj.M = 3; end
            if isempty(obj.D); obj.D = obj.M+19; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj = zeros(size(PopDec,1),obj.M);
            g      = 1+9*mean(PopDec(:,obj.M:end),2);
            PopObj(:,1:obj.M-1) = PopDec(:,1:obj.M-1);
            PopObj(:,obj.M)     = (1+g).*(obj.M-sum(PopObj(:,1:obj.M-1)./(1+repmat(g,1,obj.M-1)).*(1+sin(3*pi.*PopObj(:,1:obj.M-1))),2));
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            interval     = [0,0.251412,0.631627,0.859401];
            median       = (interval(2)-interval(1))/(interval(4)-interval(3)+interval(2)-interval(1));
            X            = UniformPoint(N,obj.M-1,'grid');
            X(X<=median) = X(X<=median)*(interval(2)-interval(1))/median+interval(1);
            X(X>median)  = (X(X>median)-median)*(interval(4)-interval(3))/(1-median)+interval(3);
            R            = [X,2*(obj.M-sum(X/2.*(1+sin(3*pi.*X)),2))];
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                x      = linspace(0,1,100)';
                y      = 2*(2-x/2.*(1+sin(3*pi*x)));
                nd     = NDSort([x,y],1)==1;
                x(~nd) = nan;
                R      = [x,y];
            elseif obj.M == 3
                [x,y]  = meshgrid(linspace(0,1,20));
                z      = 2*(3-x/2.*(1+sin(3*pi*x))-y/2.*(1+sin(3*pi*y)));
                nd     = reshape(NDSort([x(:),y(:),z(:)],1)==1,size(z));
                z(~nd) = nan;
                R      = {x,y,z};
            else
                R = [];
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%% Display Start message
disp('Running PSK SER Simulation set')
%% Plot theoretical curves
SNRs = -4:24;
[h_fig, h_lines] = PSK_SER_Curves(SNRs);
%% Run Monte Carlo simulations

PSK_SER = zeros(5,length(SNRs));
hold on
simLines = semilogy(SNRs, PSK_SER,'*');

tic
[PSK_BER, PSK_SER] = PSK_Simulate(SNRs, simLines);
toc
%% Plot simulation results
% Simulation results are plotted inside the PSK_Simulate function
%% Add Legend
% Do MATLAB graphics magic to create concise legend
% See "Controling Legends" in MATLAB doc

% group lines together
simGrp = hggroup('DisplayName','Simulation');
theoGrp = hggroup('DisplayName','Theoretical');
set(simLines,'Parent',simGrp)
set(h_lines,'Parent',theoGrp)
set(get(get(simGrp,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','on'); % Include this hggroup in the legend
set(get(get(theoGrp,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','on'); % Include this hggroup in the legend
legend show

disp('Note: Fig. 5-2-10 in Proakis 3rd Edition is labeled incorrectly.')
disp('The y-axis in that figure should be labeled from 10^-5 at the bottom')
disp('rather than 10^-6.  All other labels should be corrected accordingly.')
disp('This figure has been corrected in the 4th Edition.')


text:
function [X maxdot] = packing_on_the_sphere(d, n, epsilon, X0)
% Return a set of points spread out on the sphere.
%
% function [X maxdot] = packing_on_the_sphere(d, n, epsilon, X0)
%
% Using optimization on the oblique manifold, that is, the product of
% spheres, this function returns a set of n points with unit norm in R^d in
% the form of a matrix X of size nxd, such that the points are spread out
% on the sphere. Ideally, we would minimize the maximum inner product
% between any two points X(i, :) and X(j, :), i~=j, but that is a nonsmooth
% cost function. Instead, we replace the max function by a classical
% log-sum-exp approximation and (attempt to) solve:
%
% min_{X in OB(d, n)} log( .5*sum_{i~=j} exp( xi'*xj/epsilon ) ),
%
% with xi = X(:, i) and epsilon is some "diffusion constant". As epsilon
% goes to zero, the cost function is a sharper approximation of the max
% function (under some assumptions), but the cost function becomes stiffer
% and hence harder to optimize.
%
% The second output, maxdot, is the maximum inner product between any two
% points in the returned X. This number is the one we truly are trying to
% minimize.
%
% Notice that this cost function is invariant under rotation of X:
% f(X) = f(XQ) for all orthogonal Q in O(d).
% This calls for optimization over the set of symmetric positive
% semidefinite matrices of size n and rank d with unit diagonal, which can
% be thought of as the quotient of the oblique manifold OB(d, n) by O(d):
% See elliptopefactory.
%
% This is known as the Thomson or, more specifically, the Tammes problem:
% http://en.wikipedia.org/wiki/Tammes_problem
% An interesting page by Neil Sloane collecting best known packings is
% available here http://neilsloane.com/packings/

% This file is part of Manopt and is copyrighted. See the license file.
%
% Main author: Nicolas Boumal, July 2, 2013
% Contributors:
%
% Change log:
%   Aug. 14, 2013 (NB) : Code now compatible to experiment with both the
%                        obliquefactory and the elliptopefactory.
%
%   Jan.  7, 2014 (NB) : Added reference to Neil Sloane's page and the
%                        maxdot output.
%
%   June 24, 2014 (NB) : Now shifting exponentials to alleviate numerical
%                        trouble when epsilon is too small.
%

    if ~exist('d', 'var') || isempty(d)
        % Dimension of the embedding space: R^d
        d = 3;
    end
    if ~exist('n', 'var') || isempty(n)
        % Number n of points to place of the sphere in R^d.
        % For example, n=12 yields an icosahedron:
        % https://en.wikipedia.org/wiki/Icosahedron
        % Notice though that platonic solids are not always optimal.
        % Try for example n = 8: you don't get a cube.
        n = 24;
    end
    if ~exist('epsilon', 'var') || isempty(epsilon)
        % This value should be as close to 0 as affordable.
        % If it is too close to zero, optimization first becomes much
        % slower, than simply doesn't work anymore becomes of floating
        % point overflow errors (NaN's and Inf's start to appear).
        % If it is too large, then log-sum-exp is a poor approximation of
        % the max function, and the spread will be less uniform.
        % An okay value seems to be 0.01 or 0.001 for example. Note that a
        % better strategy than using a small epsilon straightaway is to
        % reduce epsilon bit by bit and to warm-start subsequent
        % optimization in that way. Trustregions will be more appropriate
        % for these fine tunings.
        epsilon = 0.0015;
    end

    % Pick your manifold (the elliptope factory quotients out the global
    % rotation invariance of the problem, which is more natural but
    % conceptually a bit more complicated --- for usage with the toolbox it
    % is the same though: just uncomment the appropriate line).
    manifold = obliquefactory(d, n, true);
    % manifold = elliptopefactory(n, d);

    % Generate a random initial guess if none was given.
    if ~exist('X0', 'var') || isempty(X0)
        X0 = manifold.rand();
    end

    % Define the cost function with caching system used: the store
    % structure we receive as input is tied to the input point X. Everytime
    % this cost function is called at this point X, we will receive the
    % same store structure back. We may modify the store structure inside
    % the function and return it: the changes will be remembered for next
    % time.
    function [f store] = cost(X, store)
        if ~isfield(store, 'ready')
            XXt = X*X';
            % Shift the exponentials by the maximum value to reduce
            % numerical trouble due to possible overflows.
            s = max(max(triu(XXt, 1)));
            expXXt = exp((XXt-s)/epsilon);
            % Zero out the diagonal
            expXXt(1:(n+1):end) = 0;
            u = sum(sum(triu(expXXt, 1)));
            store.XXt = XXt;
            store.s = s;
            store.expXXt = expXXt;
            store.u = u;
            store.ready = true;
        end
        u = store.u;
        s = store.s;
        f = s + epsilon*log(u);
    end

    % Define the gradient of the cost. When the gradient is called at a
    % point X for which the cost was already called, the store structure we
    % receive remember everything that the cost function stored in it, so
    % we can reuse previously computed elements.
    function [g store] = grad(X, store)
        if ~isfield(store, 'ready')
            [~, store] = cost(X, store);
        end
        % Compute the Euclidean gradient
        eg = store.expXXt*X / store.u;
        % Convert to the Riemannian gradient (by projection)
        g = manifold.egrad2rgrad(X, eg);
    end

    % Setup the problem structure with its manifold M and cost+grad
    % functions.
    problem.M = manifold;
    problem.cost = @cost;
    problem.grad = @grad;

    % For debugging, it's always nice to check the gradient a few times.
    % checkgradient(problem);
    % pause;

    % Call a solver on our problem with a few options defined. We did not
    % specify the Hessian but it is still okay to call trustregion: Manopt
    % will approximate the Hessian with finite differences of the gradient.
    opts.tolgradnorm = 1e-8;
    opts.maxtime = 1200;
    opts.maxiter = 1e5;
    % X = trustregions(problem, X0, opts);
    X = conjugategradient(problem, X0, opts);

    % Evaluate the maximum inner product between any two points of X.
    XXt = X*X';
    dots = XXt(find(triu(ones(n), 1))); %#ok<FNDSB>
    maxdot = max(dots);

    % Similarly, even though we did not specify the Hessian, we may still
    % estimate its spectrum at the solution. It should reflect the
    % invariance of the cost function under a global rotatioon of the
    % sphere, which is an invariance under the group O(d) of dimension
    % d(d-1)/2 : this translates into d(d-1)/2 zero eigenvalues in the
    % spectrum of the Hessian.
    % The approximate Hessian is not a linear operator, and is it a
    % fortiori not symmetric. The result of this computation is thus not
    % reliable. It does display the zero eigenvalues as expected though.
    if manifold.dim() < 300
        evs = real(hessianspectrum(problem, X));
        figure;
        stem(1:length(evs), sort(evs), '.');
        title(['Eigenvalues of the approximate Hessian of the cost ' ...
               'function at the solution']);
    end

    % Show how the inner products X(:, i)'*X(:, j) are distributed.
    figure;
    hist(real(acos(dots)), 20);
    title('Histogram of the geodesic distances');

    % This is the quantity we actually want to minimize.
    fprintf('Maximum inner product between two points: %g\n', maxdot);

    % Give some visualization if the dimension allows
    if d == 2
        % For the circle, the optimal solution consists in spreading the
        % points with angles uniformly sampled in (0, 2pi). This
        % corresponds to the following value for the max inner product:
        fprintf('Optimal value for the max inner product: %g\n', cos(2*pi/n));
        figure;
        t = linspace(-pi, pi, 201);
        plot(cos(t), sin(t), '-', 'LineWidth', 3, 'Color', [152,186,220]/255);
        daspect([1 1 1]);
        box off;
        axis off;
        hold on;
        plot(X(:, 1), X(:, 2), 'r.', 'MarkerSize', 25);
        hold off;
    end
    if d == 3
        figure;
        % Plot the sphere
        [sphere_x sphere_y sphere_z] = sphere(50);
        handle = surf(sphere_x, sphere_y, sphere_z);
        set(handle, 'FaceColor', [152,186,220]/255);
        set(handle, 'FaceAlpha', .5);
        set(handle, 'EdgeColor', [152,186,220]/255);
        set(handle, 'EdgeAlpha', .5);
        daspect([1 1 1]);
        box off;
        axis off;
        hold on;
        % Add the chosen points
        Y = 1.02*X';
        plot3(Y(1, :), Y(2, :), Y(3, :), 'r.', 'MarkerSize', 25);
        % And connect the points which are at minimal distance,
        % within some tolerance.
        min_distance = real(acos(maxdot));
        connected = real(acos(XXt)) <= 1.20*min_distance;
        [Ic Jc] = find(triu(connected, 1));
        for k = 1 : length(Ic)
            i = Ic(k); j = Jc(k);
            plot3(Y(1, [i j]), Y(2, [i j]), Y(3, [i j]), 'k-');
        end
        hold off;
    end

end

author: andrewssobral

repo: lrslibrary


text:
function one_pixel_result=seg_eva_one_img(predict_mask, gt_mask, class_info)

    exclude_class_idxes=class_info.void_class_idxes;
    class_num=class_info.class_num;

    assert(all(size(gt_mask)==size(predict_mask)));

    gt_mask_vector=gt_mask(:);
    predict_mask_vector=predict_mask(:);

    if ~isempty(exclude_class_idxes)
        valid_pixel_sel=~ismember(gt_mask, exclude_class_idxes);
        gt_mask_vector=gt_mask_vector(valid_pixel_sel);
        predict_mask_vector=predict_mask_vector(valid_pixel_sel);
    end

    [tmp_con_mat, label_vs]=confusionmat(gt_mask_vector, predict_mask_vector);
    confusion_mat=zeros(class_num, class_num);
    confusion_mat(label_vs,label_vs)=tmp_con_mat;

    one_pixel_result=seg_eva_gen_result_from_con_mat(confusion_mat, exclude_class_idxes);

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function [rtk,stat]=pppos(rtk,obs,nav)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% precise point positioning %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%input:  rtk - rtk control struct
%        obs - observations
%        nav - navigation message
%output: rtk - rtk control struct
%        stat - state (0:error 1:ok)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global glc
nobs=size(obs,1); MAXITER=8; iter=1; stat=0;
exc=zeros(nobs,1); dr=zeros(3,1);

% initialize rtk.sat.fix
for i=1:glc.MAXSAT
    for j=1:rtk.opt.nf
        rtk.sat(i).fix(j)=0;
    end
end

% debug tracing
% fprintf('\n \n \n');
% fprintf('time= %d',obs(1).time.time);fprintf('\n');
% fprintf('before time update');
% printx_ppp(rtk.x,rtk);
% printP(rtk.P,rtk);

% time update
rtk=udstate_ppp(rtk,obs,nav);

% debug tracing
% fprintf('after time update');
% printx_ppp(rtk.x,rtk);
% printP(rtk.P,rtk);

% cumpute satellite position,clock bias,velocity,clock drift
sv=satposs(obs,nav,rtk.opt.sateph);

% exclude measurements of eclipsing satellite (block IIA)
if rtk.opt.posopt(4)==1
    sv=testeclipse(obs,nav,sv);
end

% compute earth tides correction
if rtk.opt.tidecorr==1
    dr=tidedisp(gpst2utc(obs(1).time),rtk.sol.pos,7,nav.erp,nav.otlp);
end

while iter<=MAXITER

    % calculate residuals,measurement matrix,measurement noise matrix
    [v,H,R,~,exc,nv,rtk]=ppp_res(0,rtk.x,rtk,obs,nav,sv,dr,exc);
    if nv==0,break;end

    % debug tracing
%     printv(v);
%     printH(H);
%     printR(R);

    % measurement update
    [X,P,stat_tmp]=Kfilter_h(v,H,R,rtk.x,rtk.P);
    if stat_tmp==0
        [week,sow]=time2gpst(obs(1).time);
        fprintf('Warning:GPS week = %d sow = %.3f,filter error!\n',week,sow);
        stat=0;break;
    end

    % debug tracing
%     fprintf('after measurement update');
%     printx_ppp(X,rtk);
%     printP(P,rtk);

    % calculate posteriori residuals,validate the solution
    [~,~,~,~,exc,stat,rtk]=ppp_res(iter,X,rtk,obs,nav,sv,dr,exc);
    iter=iter+1;

    if stat==1
        rtk.x=X; rtk.P=P;
        break;
    end

end

if iter>MAXITER
    [week,sow]=time2gpst(obs(1).time);
    fprintf('Warning:GPS week=%d sow=%.3f, ppp iteration overflows',week,sow);
end

if stat==1
    % PPP ambiguity resolution (Not supported for the time being)

    % update solution
    rtk=update_stat(rtk,obs,glc.SOLQ_PPP);

    % hold fixed ambiguity (Not supported for the time being)

end

return

repo: GINav


text:
%% RRT*FN Toolbox
% *RRT*FN Toolbox* for solving path/motion planning problems.
% Toolbox provides *RRT*, *RRT** and *RRT*FN* for solving various problems
% (e.g. 2D mobile robot, planar redundant manipulator)
%
% Sampling-based motion planning algorithms got a lot of research attention
% in past decade or so. The main reason of the popularity is an ability
% to solve relatively easy motion planning problems of 2D mobile robot
% as well as hard high dimensional problems both.
%
% *Rapidly-Exporing Random Tree (RRT)* is sampling-based algorithm, solves
% the problem of motion and path planning providing feasible solutions
% however it does not take into account optimality of the solution.
%
% <<rrt_pseudocode.png>>
%
% *RRT** is probabilistically optimal variant of RRT which converges to the
% optimal solution asymtotically. Due to the fact that RRT* takes cost of
% the path into consideration, the computational diffuculty increases with
% each node added to tree, a good example of this problem is finding
% neighbors.
%
% <<rrt_star_pseudocode.png>>
%
% <<rrt_star_pseudocode_1.png>>
%
% *RRT*FN (Fixed Nodes)* is a memory efficient variant of RRT*. It
% inherents all the optimization features of RRT*, but in contrast using
% limited number of nodes i.e. memory is limited.
%
% <<rrt_star_fn_pseudocode.png>>
%
%%
%
% <<rrts_comparison.png>>
%%
%

% Please do not change anything in rrt.m,
% rrt_star.m and rrt_star_fn.m files unless it is a bug.
Everything
% you want to add (e.g. new function or modifications of a model)
% please do it creating/editing classes for models.
% e.g. FNSimple2D.m or
% FNRedundantManipulator.m
%

%

author: olzhas


text:
%NBAYESC Bayes Classifier for given normal densities
%
%   W = NBAYESC(U,G)
%
% INPUT
%   U  Dataset of means of classes
%   G  Covariance matrices (optional; default: identity matrices)
%
% OUTPUT
% 	W  Bayes classifier
%
% DESCRIPTION
% Computation of the Bayes normal classifier between a set of classes.
% The means, labels and priors are defined by the dataset U of the size
% [C x K]. The covariance matrices are stored in a matrix G of the
% size [K x K x C], where K and C correspond to the dimensionality and
% the number of classes, respectively.
%
% If C is 1, then G is treated as the common covariance matrix, yielding
% a linear solution. For G = I, the nearest mean solution is obtained.
%
% This routine gives the exact solution for the given parameters, while
% the trainable classifiers QDC and LDC give approximate solutions, based
% on the parameter estimates from a training set. For a given dataset, U
% and G can be computed by MEANCOV.
%
% EXAMPLES
% [U,G] = MEANCOV(GENDATB(25));
% W = NBAYESC(U,G);
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, QDC, LDC, NMC.

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: nbayesc.m,v 1.4 2009/01/04 21:11:07 duin Exp $

function W = nbayesc(U,G);

		[cu,ku] = size(U);		% CU is the number of classes and KU - the dimensionality
	if nargin == 1,
		prwarning(4,'Covariance matrix is not specified, the identity matrix is assumed.');
		G = eye(ku);
	end

	[k1,k2,c] = size(G);	% C = 1, if G is the common covariance matrix.

	if (c ~= 1 & c ~= cu) | (k1 ~= k2) | (k1 ~= ku)
		error('Covariance matrix or a set of means has a wrong size.')
	end

	pars.mean  = +U;
	pars.cov   = G;
	pars.prior = getprior(U);

	%W = prmapping('normal_map','trained',pars,getlablist(U),ku,cu);
	W = normal_map(pars,getlablist(U),ku,cu);
	W = setname(W,'BayesNormal');
	W = setcost(W,U);

return;

author: marianux

repo: ecg-kit

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


text:
function [dist, proj] = distancePointEllipse(point, elli)
%DISTANCEPOINTELLIPSE Distance from a point to an ellipse.
%
%   DIST = distancePointEllipse(POINT, ELLI)
%   Computes the Euclidean distance between the point POINT and the ellipse
%   ELLI.
%   POINT may also be a N-by-2 array of point coordinates. In that case the
%   result is a N-by-1 array of distances.
%
%   [DIST, PROJ] = distancePointEllipse(POINT, ELLI)
%   Also return the coordinates of the projection of the point onto the
%   ellipse. PROJ has same dimensions as the array POINT.
%
%   Example
%     % create an ellipse
%     elli = [50 50 40 30 20];
%     % generate points along a regular grid
%     [x, y] = meshgrid(1:100, 1:100);
%     pts = [x(:) y(:)];
%     % compute distance map
%     distMap = reshape(distancePointEllipse(pts, elli), size(x));
%     figure; imshow(distMap, []); colormap parula;
%
%
%   References:
%     https://blog.chatfield.io/simple-method-for-distance-to-ellipse/
%
%   See also
%     ellipses2d, projPointOnEllipse, distancePoints
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2022-07-17, using Matlab 9.12.0.1884302 (R2022a)
% Copyright 2022 INRAE - BIA Research Unit - BIBS Platform (Nantes)

proj = projPointOnEllipse(point, elli);

dist = distancePoints(point, proj, 'diag');


text:
function imdb = setup_imdb_generic(datasetDir, varargin)

opts.seed = 0 ;             % random seed generator
opts.ratio = [0.7 0.3];     % train:val ratio
opts.ext = '.jpg';          % extension of target files
opts.per_class_limit = inf; % inf indicates no limit
opts = vl_argparse(opts, varargin);

opts.ratio = opts.ratio(1:2)/sum(opts.ratio(1:2));
rng(opts.seed);
imdb.imageDir = datasetDir;

% meta
fprintf('Scanning for classes ... ');
contents = dir(imdb.imageDir);
contents_name = {contents.name};
imdb.meta.classes = setdiff(contents_name(cell2mat({contents.isdir})),{'.','..'});
imdb.meta.sets = {'train', 'val', 'test'};
fprintf('%d classes found! \n', numel(imdb.meta.classes));

% images
imdb.images.name    = {};
imdb.images.class   = [];
imdb.images.set     = [];
fprintf('Scanning for images: \n');
for c = imdb.meta.classes,
  c = c{1};
  fprintf('\t%s ...', c);
  contents = dir(fullfile(imdb.imageDir,c,['*' opts.ext]));
  curr_name = cellfun(@(s) fullfile(c,s),{contents.name},'UniformOutput',false);
  curr_name = curr_name(1:min(numel(curr_name),opts.per_class_limit));
  imdb.images.name = [imdb.images.name curr_name];
  imdb.images.class = [imdb.images.class ones(1,numel(curr_name))*find(strcmp(imdb.meta.classes,c))];
  curr_set = ones(1,floor(opts.ratio(1)*numel(curr_name)));
  curr_set = [curr_set 2*ones(1,numel(curr_name)-numel(curr_set))];
  imdb.images.set = [imdb.images.set curr_set(randperm(numel(curr_set)))];
  fprintf(' done!\n');
end

% id
imdb.images.id = 1:length(imdb.images.name);

repo: ResNet-Matconvnet


text:
function c8_tan_test ( )

%*****************************************************************************80
%
%% C8_TAN_TEST tests C8_TAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2015
%
%  Author:
%
%    John Burkardt
%
  seed = 123456678;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8_TAN_TEST\n' );
  fprintf ( 1, '  C8_TAN computes the tangent of a C8.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '       C1=C8_UNIFORM_01          C2 = C8_TAN(C1)           C3 = C8_ATAN(C1)\n' );
  fprintf ( 1, '     ---------------------     ---------------------     ---------------------\n' );
  fprintf ( 1, '\n' );

  for test = 1 : 10

    [ c1, seed ] = c8_uniform_01 ( seed );
    c2 = c8_tan ( c1 );
    c3 = c8_atan ( c2);

    fprintf ( 1, '  (%12f  %12f)  (%12f  %12f)  (%12f  %12f)\n', ...
      real ( c1 ), imag ( c1 ), real ( c2 ), imag ( c2 ), real ( c3 ), imag ( c3 ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% ChromosomePositionHistogram
% NOTE: Histograms for RNA Polymerase
% Author: Derek Macklin, macklin@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Created: 8/18/2011
classdef ChromosomePositionHistogram
    properties (Constant)
        boundIndices = struct( ...
            'complex', {{...
                82, ...                                             % SMC
                50, ...                                             % Helicase
                [180, 181, 183, 185, 187, 189, 191, 193], ...       % DnaA
                1, ...                                              % Gyrase
                78, ...                                             % Topo IV
                165, ...                                            % LuxR
                }}, ...
            'monomer', {{
                101, ...                                            % HTH regulator
                242, ...                                            % Ferric uptake repressor
                }}...
            );
        boundIndicesNames = struct( ...
            'complex', {{...
                'SMC', ...
                'Helicase', ...
                'DnaA', ...
                'Gyrase', ...
                'Topo IV', ...
                'LuxR'
                }}, ...
            'monomer', {{
                'HTH regulator', ...
                'Ferric uptake\nrepressor', ...
                }}...
            );
    end

    methods (Static = true)
        function run(simDir, fileName)
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;
            import edu.stanford.covert.cell.sim.util.DiskLogger;
            import edu.stanford.covert.cell.sim.util.PlotUtil;
            import edu.stanford.covert.cell.sim.analysis.ChromosomePositionHistogram;

            %% options
            if nargin < 1 || isempty(simDir)
                simDir = SimulationDiskUtil.getLatestSimulation();
            end

            [simDir, ~, sim] = SimulationDiskUtil.getSimulation(simDir);
            c = sim.state('Chromosome');

            %% get data
            stateNames = {
                'ProteinComplex' 'counts'
                'Chromosome' 'complexBoundSites'
                'Chromosome' 'monomerBoundSites'
                'RNAPolymerase' 'positionStrands'
                'RNAPolymerase' 'states'
                'Time' 'values'};
            states = DiskLogger.load(simDir, stateNames, [], [], [], 'extract');

            cDensityMatrix = ChromosomePositionHistogram.makeDensityMatrix(states, ChromosomePositionHistogram.boundIndices, sim);

            %% figure 1
            [~, figHandle1] = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle();
            clf(figHandle1);
            nPlots = numel(ChromosomePositionHistogram.boundIndices.complex);
            [axesHandles, xAxesHandle] = PlotUtil.multiElementPlot(...
                figHandle1, 3.5 * ones(nPlots, 1), [1 c.sequenceLen], ...
                struct(...
                'titleStr', 'Protein Density on Chromosome 1', ...
                'xlabelStr', 'Position on Chromosome 1'));
            set(xAxesHandle, 'XTick', [1 c.terCPosition c.sequenceLen], 'XTickLabel', {'oriC', 'terC', 'oriC'});

            for i = 1:size(cDensityMatrix, 2)
                PlotUtil.plotLine(axesHandles(i), 1:size(cDensityMatrix, 1), cDensityMatrix(:, i), false, true, false);
                ylabel(axesHandles(i), {ChromosomePositionHistogram.boundIndicesNames.complex{i} '(Counts)'});
            end
            if exist('fileName', 'var') && ~isempty(fileName)
                saveas(figHandle1, sprintf('%s-1.pdf', fileName));
                close(figHandle1);
            end

            %% figure 2
            [~, figHandle2] = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle();
            ChromosomePositionHistogram.plotRnaPolymeraseDensity(sim, states, figHandle2);

            if exist('fileName', 'var') && ~isempty(fileName)
                saveas(figHandle2, sprintf('%s-2.pdf', fileName));
                close(figHandle2);
            end
        end
    end

    %multiple simulations
    methods (Static = true)

        % TODO: Handle cases where we don't have one of monomers or complexes
        function [cDensityMatrix, mDensityMatrix] = makeDensityMatrix(states, boundIndices, simulation)
            complexIdxs = boundIndices.complex;
            monomerIdxs = boundIndices.monomer;

            cbs = states.Chromosome.complexBoundSites;
            mbs = states.Chromosome.monomerBoundSites;

            ch = simulation.state('Chromosome');

            cDensityMatrix = zeros(size(cbs, 1), numel(complexIdxs));
            for i = 1:numel(complexIdxs)
                ftpt = ones(ch.complexDNAFootprints(complexIdxs{i}(1)), 1);
                for j = 1:numel(complexIdxs{i})
                    complexIdx = complexIdxs{i}(j);
                    cDensityMatrix(:, i) = cDensityMatrix(:, i) + ...
                        cconv(...
                        full(sum(sum(cbs(:, 1:2, :) == complexIdx, 3), 2)), ...
                        ftpt, ...
                        size(cDensityMatrix, 1) ...
                        );
                end
            end
            cDensityMatrix = cDensityMatrix / size(cbs, 3);

            mDensityMatrix = zeros(size(mbs, 1), numel(monomerIdxs));
            for i = 1:numel(monomerIdxs)
                ftpt = ones(ch.monomerDNAFootprints(monomerIdxs{i}(1)), 1);
                for j = 1:numel(monomerIdxs{i})
                    monomerIdx = monomerIdxs{i}(j);
                    mDensityMatrix(:, i) = mDensityMatrix(:, i) + ...
                        cconv(...
                        full(sum(sum(mbs(:, 1:2, :) == monomerIdx, 3), 2)), ...
                        ftpt, ...
                        size(mDensityMatrix, 1) ...
                        );
                end
            end
            mDensityMatrix = mDensityMatrix / size(mbs, 3);

            % Hack to deal with rounding issues that lead to small
            % negative numbers
            cDensityMatrix = max(0, cDensityMatrix);
            mDensityMatrix = max(0, mDensityMatrix);
        end

        function plotRnaPolymeraseDensity(sim, states, figHandle)
            import edu.stanford.covert.cell.sim.util.DiskLogger;
            import edu.stanford.covert.cell.sim.util.PlotUtil;
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;
            import edu.stanford.covert.util.CircularSparseMat;

            c = sim.state('Chromosome');
            pol = sim.state('RNAPolymerase');
            pc = sim.state('ProteinComplex');

            nTimePoints = size(states.RNAPolymerase.positionStrands, 3);
            nRnaPol = size(states.RNAPolymerase.positionStrands, 1);
            polStates = reshape(states.RNAPolymerase.states, [nRnaPol * nTimePoints 1]);
            posStrndTimes = reshape(permute([states.RNAPolymerase.positionStrands repmat(permute(1:nTimePoints, [1 3 2]), [nRnaPol 1])], [1 3 2]), [nRnaPol * nTimePoints 3]);

            polStates = polStates(posStrndTimes(:, 1) ~= 0, :);
            posStrndTimes = posStrndTimes(posStrndTimes(:, 1) ~= 0, :);

            % Get rid of strands 3 and 4 (i.e., only look at Chromosome 1)
            polStates = polStates(posStrndTimes(:, 2) <= 2, :);
            posStrndTimes = posStrndTimes(posStrndTimes(:, 2) <= 2, :);

            boundSites = CircularSparseMat(posStrndTimes, polStates, [c.sequenceLen 4 nTimePoints], 1);

            activeDensity = cconv(...
                full(sum(sum(boundSites >= pol.activelyTranscribingValue, 3), 2)), ...
                ones(c.complexDNAFootprints(pc.getIndexs('RNA_POLYMERASE')), 1), ...
                c.sequenceLen) / nTimePoints;
            nonSpecBoundDensity = cconv(...
                full(sum(sum(boundSites == pol.nonSpecificallyBoundValue, 3), 2)), ...
                ones(c.complexDNAFootprints(pc.getIndexs('RNA_POLYMERASE')), 1), ...
                c.sequenceLen) / nTimePoints;
            specBoundDensity = cconv(...
                full(sum(sum(boundSites == pol.specificallyBoundValue, 3), 2)), ...
                ones(c.complexDNAFootprints(pc.getIndexs('RNA_POLYMERASE')), 1), ...
                c.sequenceLen) / nTimePoints;

            % Hack to deal with rounding issues that lead to small
            % negative numbers
            activeDensity = max(0, activeDensity);
            nonSpecBoundDensity = max(0, nonSpecBoundDensity);
            specBoundDensity = max(0, specBoundDensity);

            %layout axes
            clf(figHandle);
            [axesHandles, xAxesHandle] = PlotUtil.multiElementPlot(figHandle, 8 * ones(3, 1), [1 c.sequenceLen], struct( ...
                'titleStr', 'RNA Polymerase Density on Chromosome 1', ...
                'xlabelStr', 'Position on Chromosome 1'));
            set(xAxesHandle, 'XTick', [1 c.terCPosition c.sequenceLen], 'XTickLabel', {'oriC', 'terC', 'oriC'});

            %plot data
            PlotUtil.plotLine(axesHandles(1), 1:c.sequenceLen, activeDensity, false, true, false);
            ylabel(axesHandles(1), 'Active');

            PlotUtil.plotLine(axesHandles(2), 1:c.sequenceLen, specBoundDensity, false, true, false);
            ylabel(axesHandles(2), {'Specifically' 'Bound'});

            PlotUtil.plotLine(axesHandles(3), 1:c.sequenceLen, nonSpecBoundDensity, false, true, false);
            ylabel(axesHandles(3), {'Non-' 'Specifically' 'Bound'});

            % Format Y-axis
            PlotUtil.alignYAxesLabels(axesHandles);
            PlotUtil.offsetYAxes(axesHandles, 0.015);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function T = GenPencil(im, P, J)
% ==============================================
%   Compute the pencil map 'T'
%
%   Paras:
%   @im        : input image ranging value from 0 to 1.
%   @P         : the pencil texture.
%   @J         : the tone map.
%

    %% Parameters
    theta = 0.2;

    [H, W, ~] = size(im);

    %% Initialization
    P = imresize(P, [H, W]);
    P = reshape(P, H*W, 1);
    logP = log(P);
    logP = spdiags(logP, 0, H*W, H*W);

    J = imresize(J, [H, W]);
    J = reshape(J, H*W, 1);
    logJ = log(J);

    e = ones(H*W, 1);
    Dx = spdiags([-e, e], [0, H], H*W, H*W);
    Dy = spdiags([-e, e], [0, 1], H*W, H*W);

    %% Compute matrix A and b
    A = theta * (Dx * Dx' + Dy * Dy') + (logP)' * logP;
    b = (logP)' * logJ;

    %% Conjugate gradient
    beta = pcg(A, b, 1e-6, 60);

    %% Compute the result
    beta = reshape(beta, H, W);

    P = reshape(P, H, W);

    T = P .^ beta;
end

repo: PencilDrawing


text:
function v = wrap(u,N)
% WRAP Wrap a vector of indices around a torus.
% v = wrap(u,N)
%
% e.g., wrap([-1 0 1 2 3 4], 3)   =   2 3 1 2 3 1

v = mod(u-1,N)+1;


text:
filename='Arch_hexahedra_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'MMA'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.1;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function drawText(x,y, drawString, fontSize, fontStyle, color, font, fontWeight, flagCenter)
%drawText Draws text to figure
%
% drawText(x,y, drawString, fontSize, fontStyle, color, font, fontWeight, flagCenter)
%
%INPUTS
% x             x coordinate
% y             y coordinate
% drawString    String to draw
%
%OPTIONAL INPUTS
% fontSize      Font size (Default = 12)
% fontStyle     Font style (Default = 'normal')
% color         Font color (Default = [0 0 0])
% font          Font (Default = 'sans-serif')
% fontWeight    Font weight (Default = 'normal')
% flagCenter    Center text on coordinate (Default = false)
%
global CB_MAP_OUTPUT
global mapHandle

if nargin < 4
    fontSize = 12;
end
if nargin < 5
    fontStyle = 'normal';
end
if nargin < 6
    color = [0 0 0];
end
if nargin < 7
    font = 'sans-serif';
end
if nargin <8
    fontWeight = 'normal';
end
if nargin < 9
    flagCenter = false;
end

if length(fontSize)>1
    error('Fontsize should be a scalar')
end

if iscell(drawString)
    display ('whoops - this shouldnt happen.  drawText.m')
%     drawString
%     whos('drawString')
    drawString = drawString{1};
end
if isnan(x) | isnan(y) | isnan(fontSize)
%     display ('whoops - this shouldnt happen.  drawText.m')
    return;
end
if strcmp(CB_MAP_OUTPUT, 'matlab')
    if find(color>1)
        color = color/255;
    end
    if flagCenter
        text(x,-y,drawString, 'FontSize', fontSize/2, 'color', color,'HorizontalAlignment','center');
    else
        text(x,-y,drawString, 'FontSize', fontSize/2, 'color', color);
    end
elseif strcmp(CB_MAP_OUTPUT, 'java')
    % need to insert code
elseif strcmp(CB_MAP_OUTPUT, 'svg')
    %determine type of color input
    if ischar(color)
        if color=='k'
            colorStroke = 'rgb(0,0,0)';
        else
            colorStroke = color;
        end
    else if isvector(color)
            colorStroke = strcat('rgb(',num2str(color(1)),',',num2str(color(2)),',',num2str(color(3)),')');
        end
    end
    drawString = regexprep(drawString,'&','&');
    fprintf(mapHandle, '<g style="font-family: %s; font-style: %s; font-weight: %s; stroke: none;">\n',font,fontStyle, fontWeight);
    if flagCenter
        %fprintf('<text style="fill: %s; text-rendering: optimizeLegibility;" x="%.2f" y="%.2f" font-size="%dpx">%s</text>\n</g>\n',colorStroke,x-(length(drawString)*fontSize*1.25),y+(5*fontSize),5*fontSize,drawString);
        fprintf(mapHandle,'<text style="fill: %s; text-rendering: optimizeLegibility;" x="%.2f" y="%.2f" font-size="%dpx">%s</text>\n</g>\n',colorStroke,x-(length(drawString)*fontSize*1.25),y+(5*fontSize),5*fontSize,drawString);
    else
        %code for debugging
%         fprintf('<text style="fill: %s; text-rendering: optimizeLegibility;" x="%.2f" y="%.2f" font-size="%dpx">%s</text>\n</g>\n',colorStroke,x,y,5*fontSize,drawString);
%         fprintf('\n%s\n','----')
%         disp(colorStroke)
%         class(colorStroke)
%         fprintf('%s\n','----')
%         disp(x)
%         class(x)
%         fprintf('%s\n','----')
%         disp(y)
%         class(y)
%         fprintf('%s\n','----')
%         disp(5*fontSize)
%         class(5*fontSize)
%         fprintf('%s\n','----')
%         disp(drawString)
%         class(drawString)
        fprintf(mapHandle,'<text style="fill: %s; text-rendering: optimizeLegibility;" x="%.2f" y="%.2f" font-size="%dpx">%s</text>\n</g>\n',colorStroke,x,y,5*fontSize,drawString);
    end
%     fprintf(mapHandle,'<text style="stroke: none;" x="%8.2f" y="%8.2f">%s</text>\n',x,-y,drawString);
else
    display('error CB_MAP_OUTPUT in drawText');
end

author: opencobra

repo: cobratoolbox


text:
function template = PTKCreateTemplateForLobe(lobe_mask, context, reporting)
    % PTKCreateTemplateForLobe. Function for creating a template image
    %     for a lobar regon
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    template = lobe_mask.Copy;
end

author: tomdoel

repo: pulmonarytoolkit


text:
function fixed_num = p12_fixed_num ( )

%*****************************************************************************80
%
%% P12_FIXED_NUM returns the number of fixed points in problem 12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer FIXED_NUM, the number of fixed points.
%
  fixed_num = 12;

  return
end


text:
clear all; close all;
I=imread('rice.png');
I=im2double(I);
LEN=20;
THETA=10;
PSF=fspecial('motion', LEN, THETA);
J=imfilter(I, PSF, 'conv', 'circular');
figure;
subplot(121);  imshow(I);
subplot(122);  imshow(J);
noise=0.03*randn(size(I));
K=imadd(J, noise);
NP=abs(fft2(noise)).^2;
NPower=sum(NP(:))/prod(size(noise));
NCORR=fftshift(real(ifft2(NP)));
IP=abs(fft2(I)).^2;
IPower=sum(IP(:))/prod(size(I));
ICORR=fftshift(real(ifft2(IP)));
L=deconvwnr(K, PSF, NCORR, ICORR);
figure;
subplot(121);  imshow(K);
subplot(122);  imshow(L);

author: HuangCongQing