text:
function spm_MDP_plot(MDP)
% creates a movie of hierarchical expectations and outcomes
% FORMAT spm_MDP_plot(MDP))
%
% MDP - nested MDP (and DEM) structures
%     - (requires fields to specify the labels of states and outcomes)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MDP_plot.m 6977 2016-12-24 17:48:44Z karl $

% Preliminaries
%--------------------------------------------------------------------------
cla,axis([1 20 0 8]); axis image, hold on, box on
set(gca,'Fontsize',16,'XColor','r','TickLength',[0 0],'XTick',[],'YTick',[])
title('Narrative construction','FontSize',16)

% switch off basis functions for visual stimuli
%--------------------------------------------------------------------------
try
global STIM
B = STIM.B;
STIM.B = 1;
end

% hidden factors to display (one per level)
%--------------------------------------------------------------------------
scale  = 1 - 1/16;                             % background level for text
try k1 = MDP.label.k;     catch, k1 = 1; end
try k2 = MDP.MDP.label.k; catch, k2 = 1; end

% cycle over highest level
%--------------------------------------------------------------------------
T1     = size(MDP.xn{k1},4);
F1     = fix(1024/(T1*numel(MDP.label.name{k1}{1})));
F1     = min(max(F1,8),32);
for t1 = 1:T1

    % draw stuff
    %----------------------------------------------------------------------
    p     = 1 - squeeze(MDP.xn{k1}(end,:,:,t1));
    p     = p*scale;
    str   = MDP.label.name{k1};
    ns    = size(p,1);
    nt    = size(p,2);
    for i = 1:nt
        for j = 1:ns
            try
                set(h1(i,j),'Color',[1 1 1]*p(j,i));
            catch
                h1(i,j) = text(16*i/nt,j/ns + 6,str(j),...
                    'HorizontalAlignment','center',...
                    'Color',[1 1 1]*p(j,i),'FontSize',F1);
            end
        end
    end

    % next level down
    %----------------------------------------------------------------------
    if isfield(MDP,'mdp')
        T2     = size(MDP.mdp(t1).xn{k2},4);
        F2     = fix(1024/(T2*numel(MDP.MDP.label.name{k1}{1})));
        F2     = min(max(F2,8),32);
        for t2 = 1:T2

            % draw stuff
            %--------------------------------------------------------------
            p     = 1 - squeeze(MDP.mdp(t1).xn{k2}(end,:,:,t2));
            p     = p*scale;
            str   = MDP.MDP.label.name{k2};
            ns    = size(p,1);
            nt    = size(p,2);
            for i = 1:nt
                for j = 1:ns
                    try
                        set(h2(i,j),'Color',[1 1 1]*p(j,i));
                    catch
                        h2(i,j) = text(16*i/nt,j/ns + 4,str(j),...
                            'HorizontalAlignment','center',...
                            'Color',[1 1 1]*p(j,i),'FontSize',F2);
                    end
                end
            end

            % check for extra time points
            %--------------------------------------------------------------
            try
                delete(h2((nt + 1):end,:))
            end

            % next level down
            %--------------------------------------------------------------
            if isfield(MDP.mdp(t1),'dem')

                % accumulated expectations from Bayesian filtering
                %----------------------------------------------------------
                o  = MDP.mdp(t1).dem(t2).X;
                no = numel(o);
                T  = size(o{1},2);

            else

                % observed outcomes
                %----------------------------------------------------------
                o  = MDP.mdp(t1).o;
                no = size(o,1);
                T  = 16;
            end

            % outcome modalities
            %--------------------------------------------------------------
            set(gca,'XTickLabel',MDP.MDP.label.modality)
            set(gca,'XTick',16*(1:no)/no)

            F3     = fix(1024/(no*numel(MDP.MDP.label.outcome{1}{1})));
            F3     = min(max(F3,8),32);
            for t3 = 1:T

                % draw stuff
                %----------------------------------------------------------
                for i = 1:no

                    if iscell(o)

                        % probabilistic outcomes (expectations)
                        %--------------------------------------------------
                        x     = 16*i/no;
                        p     = spm_softmax(log(o{i}(:,t3))/4);
                        str   = MDP.MDP.label.outcome{i};
                        [p,j] = sort(p,'descend');
                        col   = [1 1 1] - [0 1 1]*(p(1) - 1/numel(p));
                        try
                            set(h3(i),'String',str(j(1)),...
                                'Color',col);
                        catch
                            h3(i) = text(x,3,str(j(1)),...
                                'HorizontalAlignment','center',...
                                'Color',col,'FontSize',F3);
                        end

                        % observations in image format
                        %--------------------------------------------------
                        x     = [-1 1] + x;
                        y     = [-1 1] + 1;
                        try
                            px = MDP.mdp(t1).dem(t2).pU.x{1}(:,t3);
                            pv = MDP.mdp(t1).dem(t2).pU.v{2}(:,t3);
                            px = spm_unvec(px,MDP.MDP.DEM.G(1).x);
                            pv = spm_unvec(pv,MDP.MDP.DEM.G(2).v);
                            Y  = MDP.MDP.DEM.label{i}(px,pv);
                            try
                                delete(h4(i))
                            end
                            h4(i) = image(x,y,flipud(Y)*64);

                        end

                    else

                        % deterministic outcomes
                        %--------------------------------------------------
                        str = MDP.MDP.label.outcome{i}{o(i,t2)};
                        try
                            set(h3(i),'String',str);
                        catch
                            h3(i) = text(16*i/no,2,str,'Color','r',...
                                'HorizontalAlignment','center','FontSize',F3);
                        end
                    end
                end

                % plot timelines
                %----------------------------------------------------------
                x  = [1 1]*(t3 + (t2 - 1)*T + (t1 - 1)*T*MDP.MDP.T);
                x  = 1 + 16*x/(MDP.MDP.T*MDP.T*T);
                y  = [7.2 8];
                try
                    set(l1,'XData',x);
                catch
                    l1 = line(x,y,'Color','r','LineWidth',8);
                end

                x  = [1 1]*(t3 + (t2 - 1)*T);
                x  = 1 + 16*x/(T*T2);
                y  = [5.2 6];
                try
                    set(l2,'XData',x);
                catch
                    l2 = line(x,y,'Color','r','LineWidth',8);
                end

                x  = [1 1]*t3;
                x  = 1 + 16*x/T;
                y  = [3.2 4];
                try
                    set(l3,'XData',x);
                catch
                    l3 = line(x,y,'Color','r','LineWidth',8);
                end

                % and save graphics
                %----------------------------------------------------------
                try
                    M(end + 1) = getframe(gca);
                catch
                    M = getframe(gca);
                end

            end
        end

    else

        % observed outcomes
        %------------------------------------------------------------------
        o      = MDP.o;
        no     = size(o,1);
        T      = 16;
        F3     = fix(1024/(no*numel(MDP.label.outcome{1}{1})));
        F3     = min(max(F3,8),48);

        % outcome modalities
        %------------------------------------------------------------------
        set(gca,'XTickLabel',MDP.label.modality)
        set(gca,'XTick',16*(1:no)/no)

        for t3 = 1:T

            % draw stuff
            %--------------------------------------------------------------
            for i = 1:no
                str = MDP.label.outcome{i}{o(i,t1)};
                try
                    set(h3(i,j),'String',str);
                catch
                    h3(i,j) = text(16*i/no,2,str,'Color','r',...
                        'HorizontalAlignment','center','FontSize',F3);
                end
            end

            % plot timelines
            %--------------------------------------------------------------
            x  = [1 1]*(t3 + (t1 - 1)*T);
            x  = 1 + 16*x/(T*T1);
            y  = [7.2 8];
            try
                set(l2,'XData',x);
            catch
                l2 = line(x,y,'Color','r','LineWidth',8);
            end

            x  = [1 1]*t3;
            x  = 1 + 16*x/T;
            y  = [5.2 6];
            try
                set(l3,'XData',x);
            catch
                l3 = line(x,y,'Color','r','LineWidth',8);
            end

            % and save graphics
            %--------------------------------------------------------------
            try
                M(end + 1) = getframe(gca);
            catch
                M = getframe(gca);
            end

        end
    end
end

% save movie
%--------------------------------------------------------------------------
try, STIM.B = B; end
set(gca,'Userdata',{M,16})
set(gca,'ButtonDownFcn','spm_DEM_ButtonDownFcn')

author: spm


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Released under the MIT License.
% If you use this code, please cite the following paper:
% Mahmoud Afifi, Abdelrahman Abdelhamed, Abdullah Abuolaim, Abhijith
% Punnappurath, and Michael S Brown. CIE XYZ Net: Unprocessing Images for
% Low-Level Computer Vision Tasks. arXiv preprint, 2020.
%
% Author: Mahmoud Afifi | Email: mafifi@eecs.yorku.ca, m.3afifi@gmail.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

classdef mae_cos_LossLayer < nnet.layer.RegressionLayer
    properties
        scale
    end
    methods
        function layer = mae_cos_LossLayer(name, scaleFactor)
            layer.Name = name;
            layer.Description = ...
                'Mean abs error + (negative) cos similarity';
            layer.scale = scaleFactor;
        end

        function loss = forwardLoss(layer, Y, T)
            % Calculate MAE + (negative) cos similarity.
            Y(Y<0) = 0;

            sz = size(Y);
            if length(sz) == 3
                R = 1;
                mae_srgb = sum(reshape(abs(Y(:,:,4:6) - T(:,:,4:6)),...
                    [],3),2);
                mae_xyz = sum(reshape(abs(Y(:,:,1:3) - T(:,:,1:3)),...
                    [],3),2);

                % compute cosine simialrity
                cosSim_srgb = zeros(size(mae_srgb), 'like', Y);
                cosSim_xyz = zeros(size(mae_xyz), 'like', Y);

                y_srgb = reshape(Y(:,:,3:6),[],3);
                t_srgb = reshape(T(:,:,4:6),[],3);
                y_xyz = reshape(Y(:,:,1:3),[],3);
                t_xyz = reshape(T(:,:,1:3),[],3);

                cosSim_srgb = cosSim_srgb + sum(y_srgb.*t_srgb,2)./...
                    (sqrt(sum(y_srgb.^2,2)) .* ...
                    sqrt(sum(t_srgb.^2,2)) + eps);

                cosSim_xyz = cosSim_xyz + sum(y_xyz.*t_xyz,2)./...
                    (sqrt(sum(y_xyz.^2,2)) .* ...
                    sqrt(sum(t_xyz.^2,2)) + eps);

            else
                R = sz(4);
                mae_srgb = zeros(sz(1)*sz(2),1,'like',Y);
                cosSim_srgb = zeros(size(mae_srgb),'like',mae_srgb);

                mae_xyz = zeros(sz(1)*sz(2),1,'like',Y);
                cosSim_xyz = zeros(size(mae_xyz),'like',mae_xyz);

                for j = 1 : R

                    mae_srgb = mae_srgb + sum(reshape( ...
                        abs(Y(:,:,4:6,j) - T(:,:,4:6,j)) ,...
                        [],3),2);
                    mae_xyz = mae_xyz + sum(reshape( ...
                        abs(Y(:,:,1:3,j) - T(:,:,1:3,j)) ,...
                        [],3),2);

                    % compute cosine simialrity
                    y_srgb = reshape(Y(:,:,4:6),[],3);
                    t_srgb = reshape(T(:,:,4:6),[],3);
                    y_xyz = reshape(Y(:,:,1:3),[],3);
                    t_xyz = reshape(T(:,:,1:3),[],3);

                    cosSim_srgb = cosSim_srgb + sum(y_srgb.*t_srgb,2)./...
                        (sqrt(sum(y_srgb.^2,2)) .* ...
                        sqrt(sum(t_srgb.^2,2)) + eps);

                    cosSim_xyz = cosSim_xyz + sum(y_xyz.*t_xyz,2)./...
                        (sqrt(sum(y_xyz.^2,2)) .* ...
                        sqrt(sum(t_xyz.^2,2)) + eps);
                end
            end

            loss = (sum(mae_srgb - cosSim_srgb) + ...
                layer.scale * sum(mae_xyz - cosSim_xyz))/R;

        end
    end
end

author: mahmoudnafifi


text:
% NOLM_Attractor - Plot a chaotic attractor for the NOLM.
% Copyright Springer 2013 A.L. Steele and S. Lynch.
clear
N=50000;P=100;start=500;
lambda=1.55E-6;n2=3.2E-20;Aeff=30E-12;L=80;
E1(1)=0;phi=0*pi;
Eout(1:N)=0;
phif=0*pi; E6p=0;

kappa1=0.25;kappa2=0.8;kappa3=0.8;
rootk1 = sqrt(kappa1); irootk1 = 1i*sqrt(1-kappa1);
rootk2 = sqrt(kappa2); irootk2 = 1i*sqrt(1-kappa2);
rootk3 = sqrt(kappa3); irootk3 = 1i*sqrt(1-kappa3);

G=sqrt((1-kappa2)*(1-kappa3));
% iterate
for n=1:N

    E1 = rootk3*sqrt(P)+irootk3*E6p;
    P1 = abs(E1)^2;

    E3 = rootk1*E1;
    E4 = irootk1*E1;

    phic=2*pi*n2*L*(2-kappa1)*P1/(lambda*Aeff);
    phicc=2*pi*n2*L*(1+kappa1)*P1/(lambda*Aeff);

    E3p = E3*exp(-1i*(phi+phic));
    E4p = E4*exp(-1i*(phi+phicc));

    E5 = rootk1*E3p+irootk1*E4p;

    Eout(n) = rootk2*E5;
    x(n)=real(Eout(n));
    y(n)=imag(Eout(n));
    E6 = irootk2*E5;

    E6p = E6*exp(1i*phif);

end
axis([-10 10 -10 10])
axis equal
plot(x(start:N),y(start:N),'.','MarkerSize',1);
fsize=15;
set(gca,'xtick',-10:5:10,'FontSize',fsize)
set(gca,'ytick',-10:5:10,'FontSize',fsize)
xlabel('Real E','FontSize',fsize)
ylabel('Imag E','FontSize',fsize)

author: Sable

repo: mcbench-benchmarks


text:
function y=fitpop(thetabit,bits,loss)
% J. C. Spall, August 1999
% This function returns a vector of N fitness values.
global p N thetamin thetamax
theta=zeros(p,1);
for i=1:N
   pointer=1;
   for j=1:p
      theta(j)=bit2num(thetabit(i,pointer:pointer+bits(j)-1),...
         bits(j),thetamin(j),thetamax(j));
      pointer=pointer+bits(j);
   end
   fitness(i)=-feval(loss,theta);
end
y=fitness;


text:
function score = score_all_subid(M,sid,bid,lib)
%
% For a MotorProgram M, compute the log-score of switching
% its "sub-stroke id" to every other possibility
%
% Input
%  sid: stroke index
%  bid: sub-stroke index
%
% Output
%  score: [nclust x 1] score for each possibility, using all CPDs
%     that directly depend on it
%
    assert(numel(sid)==1);
    assert(numel(bid)==1);
    ncat = lib.N;

    seq_score = zeros(ncat,1);

    % Score the transition TO the current primitive
    if bid == 1 % first primitive
        seq_score = seq_score + lib.logStart(:);
    else % not the first
        prev_subid = M.S{sid}.ids(bid-1);
        seq_score = seq_score + vec(lib.logT(prev_subid,1:ncat));
    end

    % Score the transition AWAY from the current primitive
    if bid < M.S{sid}.nsub
        next_subid = M.S{sid}.ids(bid+1);
        seq_score = seq_score + vec(lib.logT(1:ncat,next_subid));
    end

    % Scale term
    invscale_type = M.S{sid}.invscales_type(bid);
    rep_invscale_type = repmat(invscale_type,[ncat 1]);

    % Shape term
    if ~isempty(M.S{sid}.shapes_type)
        shape_type = M.S{sid}.shapes_type(:,:,bid);
        rep_shape_type = repmat(shape_type,[1 1 ncat]);
        shape_score = CPD.score_shape_type(lib,rep_shape_type,1:ncat);
    else
        shape_token = M.S{sid}.shapes_token(:,:,bid);
        rep_shape_token = repmat(shape_token,[1 1 ncat]);
        shape_score = CPD.score_shape_marginalize(lib,rep_shape_token,1:ncat);
    end

    % compute scores
    scale_score = CPD.score_invscale_type(lib,rep_invscale_type,1:ncat);
    score = seq_score + scale_score + shape_score;
end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
 function ir_colorbar_text(str)
%function ir_colorbar_text(str)
%|
%| Add a colorbar and label it with string str.
%|
%| Needed as of octave 3.8.2 and matlab R2014b because they have incompatible
%| graphics support for colorbars.

if nargin < 1, ir_usage, end
if nargin == 1 && streq(str, 'test'), ir_colorbar_text_test, return, end

if ~im, return, end

h = colorbar;
if ir_is_octave
	subplot(h); ylabel(str)
else
	h.Label.String = str;
end

function ir_colorbar_text_test
im clf
im(eye(8))
ir_colorbar_text('test1')

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function g = linardKernGradient(kern, x, varargin)

% LINARDKERNGRADIENT Gradient of LINARD kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the
% automatic relevance determination linear
% kernel's parameters. As well as the kernel structure and the
% input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in X.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations associated with the rows of the
% kernel matrix.
% ARG x2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as X1 and the same number of columns
% as X2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO linardKernParamInit, kernGradient, linardKernDiagGradient, kernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006, 2009

% KERN

  g = zeros(1, size(x, 2)+1);
  if nargin < 4
    [k, sk] = linardKernCompute(kern, x);
  else
    [k, sk] = linardKernCompute(kern, x, varargin{1});
  end
  g(1) = sum(sum(varargin{end}.*sk));
  if nargin < 4
    for i = 1:size(x, 2)
      g(1+i) =  x(:, i)'*varargin{end}*x(:, i)*kern.variance;
    end
  else
    for i = 1:size(x, 2)
      g(1+i) =  x(:, i)'*varargin{end}*varargin{1}(:, i)*kern.variance;
    end
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
I=double(imread ('vessel.png'));
Ivessel=FrangiFilter2D(I);
figure,
subplot(1,2,1), imshow(I,[]);
subplot(1,2,2), imshow(Ivessel,[0 0.25]);

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function hf_out = lhs_operation_joint(hf, samplesf, reg_filter, init_samplef, XH, init_hf, proj_reg)

% This is the left-hand-side operation in Conjugate Gradient

hf_out = cell(size(hf));

% Extract projection matrix and filter separately
P = cellfun(@real, hf(2,1,:), 'uniformoutput',false);
hf = hf(1,1,:);

% Get sizes
num_features = length(hf);
filter_sz = zeros(num_features,2);
for k = 1:num_features
    filter_sz(k,:) = [size(hf{k},1), size(hf{k},2)];
end
[~, k1] = max(filter_sz(:,1));  % Index for the feature block with the largest spatial size
block_inds = 1:num_features;
block_inds(k1) = [];
output_sz = [size(hf{k1},1), 2*size(hf{k1},2)-1];

% Compute the operation corresponding to the data term in the optimization
% (blockwise matrix multiplications)
%implements: A' diag(sample_weights) A f

% sum over all features and feature blocks
pad_sz = cell(1,1,num_features);
sh = mtimesx(samplesf{k1}, permute(hf{k1}, [3 4 1 2]), 'speed');    % assumes the feature with the highest resolution is first
for k = block_inds
    pad_sz{k} = (output_sz - [filter_sz(k,1), 2*filter_sz(k,2)-1]) / 2;

    sh(:,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end) = ...
        sh(:,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end) + mtimesx(samplesf{k}, permute(hf{k}, [3 4 1 2]), 'speed');
end

% weight all the samples
% sh = bsxfun(@times,sample_weights,sh);

% multiply with the transpose
hf_out1 = cell(1,1,num_features);
hf_out1{k1} = permute(conj(mtimesx(sh, 'C', samplesf{k1}, 'speed')), [3 4 2 1]);
for k = block_inds
    hf_out1{k} = permute(conj(mtimesx(sh(:,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end), 'C', samplesf{k}, 'speed')), [3 4 2 1]);
end

% compute the operation corresponding to the regularization term (convolve
% each feature dimension with the DFT of w, and the tramsposed operation)
% add the regularization part
% hf_conv = cell(1,1,num_features);
for k = 1:num_features
    reg_pad = min(size(reg_filter{k},2)-1, size(hf{k},2)-1);

    % add part needed for convolution
    hf_conv = cat(2, hf{k}, conj(rot90(hf{k}(:, end-reg_pad:end-1, :), 2)));

    % do first convolution
    hf_conv = convn(hf_conv, reg_filter{k});

    % do final convolution and put toghether result
    hf_out1{k} = hf_out1{k} + convn(hf_conv(:,1:end-reg_pad,:), reg_filter{k}, 'valid');
end

% Stuff related to the projection matrix

% B * P
BP_cell = cell(1,1,num_features);
for k = 1:num_features
    BP_cell{k} = mtimesx(mtimesx(init_samplef{k}, P{k}, 'speed'), init_hf{k}, 'speed');
end

BP = BP_cell{k1};
for k = block_inds
    BP(1,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end) = ...
        BP(1,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end) + BP_cell{k};
end

% multiply with the transpose: A^H * BP
hf_out{1,1,k1} = hf_out1{k1} +  permute(bsxfun(@times, BP, conj(samplesf{k1})), [3 4 2 1]);

% B^H * BP
fBP = cell(1,1,num_features);
fBP{k1} = reshape(bsxfun(@times, conj(init_hf{k1}), BP), size(init_hf{k1},1), []).';

% Compute proj matrix part: B^H * A_m * f
shBP = cell(1,1,num_features);
shBP{k1} = reshape(bsxfun(@times, conj(init_hf{k1}), sh), size(init_hf{k1},1), []).';

for k = block_inds
    % multiply with the transpose: A^H * BP
    hf_out{1,1,k} = hf_out1{k} +  permute(bsxfun(@times, BP(1,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end), conj(samplesf{k})), [3 4 2 1]);

    % B^H * BP
    fBP{k} = reshape(bsxfun(@times, conj(init_hf{k}), BP(1,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end)), size(init_hf{k},1), []).';

    % Compute proj matrix part: B^H * A_m * f
    shBP{k} = reshape(bsxfun(@times, conj(init_hf{k}), sh(1,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end)), size(init_hf{k},1), []).';
end

% hf_out2 = cell(1,1,num_features);
for k = 1:num_features
    fi = size(hf{k},1) * (size(hf{k},2)-1) + 1;  % index where the last frequency column starts

    % B^H * BP
    hf_out2 = 2*real(XH{k} * fBP{k} - XH{k}(:,fi:end) * fBP{k}(fi:end,:)) + proj_reg * P{k};

    % Compute proj matrix part: B^H * A_m * f
    hf_out{2,1,k} = hf_out2 + (2*real(XH{k} * shBP{k} - XH{k}(:,fi:end) * shBP{k}(fi:end,:)));
end
end

author: martin-danelljan

repo: ECO


text:
function [node,elem]=readsmf(fname)
%
% [node,elem]=readsmf(fname)
%
% read simple model format (SMF)
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2007/11/21
%
% input:
%    fname: name of the	SMF data file
%
% output:
%    node: node coordinates of the mesh
%    elem: list of elements of the mesh
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%
node=[];
elem=[];
fid=fopen(fname,'rt');
while(~feof(fid))
    line=fgetl(fid);
    if(line(1)=='v')
        dd=sscanf(line,'v %f %f %f');
        if(length(dd)==3)
            node=[node;dd];
        end
    elseif(line(1)=='f')
        dd=sscanf(line,'f %d %d %d');
        if(length(dd)==3)
            elem=[elem;dd];
        end
    end
end
fclose(fid);
node=reshape(node,3,length(node)/3)';
elem=reshape(elem,3,length(elem)/3)';

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function x = At_fhp_rect(y, OMEGA, m, n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Shiqian Ma
% Date : 09/05/2007
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

K = length(y);

fx = zeros(m,n);
fx(1,1) = y(1);
% fx(OMEGA) = sqrt(2)*(y(2:(K+1)/2) + i*y((K+3)/2:K));
fx(OMEGA) = (y(2:(K+1)/2) + i*y((K+3)/2:K));
% x = reshape(real(n*ifft2(fx)), m*n, 1);
x = reshape(real(sqrt(m*n)*ifft2(fx)),m,n);

author: thomaskuestner


text:
function [meta] = getmeta(db,index)
 meta = db.enteries(index);
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function fsmplay(filename, speed)
% Play an fsm file as a sequence of single wavelength images.
%   fsmplay(filename, speed)
%       filename is the full path to a PE .fsm file.
%       speed is frames per second

data = fsmload(filename);
fsmplaydata(data, speed);


text:
function variance = bradford_variance ( a, b, c )

%*****************************************************************************80
%
%% BRADFORD_VARIANCE returns the variance of the Bradford PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, the parameters of the PDF.
%    A < B,
%    0.0 < C.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  variance = ( b - a )^2 * ...
    ( c * ( log ( c + 1.0 ) - 2.0 ) + 2.0 * log ( c + 1.0 ) ) ...
    / ( 2.0 * c * ( log ( c + 1.0 ) )^2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%Plot the original STL mesh:
figure
[stlcoords] = READ_stl('sample.stl');
xco = squeeze( stlcoords(:,1,:) )';
yco = squeeze( stlcoords(:,2,:) )';
zco = squeeze( stlcoords(:,3,:) )';
[hpat] = patch(xco,yco,zco,'b');
axis equal

%Voxelise the STL:
[OUTPUTgrid] = VOXELISE(100,100,100,'sample.stl','xyz');

%Show the voxelised result:
figure;
subplot(1,3,1);
imagesc(squeeze(sum(OUTPUTgrid,1)));
colormap(gray(256));
xlabel('Z-direction');
ylabel('Y-direction');
axis equal tight

subplot(1,3,2);
imagesc(squeeze(sum(OUTPUTgrid,2)));
colormap(gray(256));
xlabel('Z-direction');
ylabel('X-direction');
axis equal tight

subplot(1,3,3);
imagesc(squeeze(sum(OUTPUTgrid,3)));
colormap(gray(256));
xlabel('Y-direction');
ylabel('X-direction');
axis equal tight


text:
function x = dvand ( n, alpha, b )

%*****************************************************************************80
%
%% DVAND solves a Vandermonde system A' * x = b.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 November 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ake Bjorck, Victor Pereyra,
%    Solution of Vandermonde Systems of Equations,
%    Mathematics of Computation,
%    Volume 24, Number 112, October 1970, pages 893-903.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real ALPHA(N), the parameters that define the matrix.
%    The values should be distinct.
%
%    Input, real B(N), the right hand side of the linear system.
%
%    Output, real X(N), the solution of the linear system.
%
  x(1:n) = b(1:n);

  for k = 1 : n - 1
    for j = n : -1 : k + 1
      x(j) = ( x(j) - x(j-1) ) / ( alpha(j) - alpha(j-k) );
    end
  end

  for k = n - 1 : -1 : 1
    for j = k : n - 1
      x(j) = x(j) - alpha(k) * x(j+1);
    end
  end

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Setting up advection-diffusion solver
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function C = please_Update_Adv_Diff_Concentration_Unsplit(C,dt,dx,dy,uX,uY,k)

% C:     concentration
% dt:    time-step
% dx,dy: spatial steps in x and y, respectively
% uX:    x-Component of Velocity
% uY:    y-Component of Velocity
% k:     diffusion coefficient

% Compute Necessary Derivatives (Note: these calculations could be parallalized)
Cx = give_Necessary_Derivative(C,dx,uX,'x');
Cy = give_Necessary_Derivative(C,dy,uY,'y');
Cxx = DD(C,dx,'x');
Cyy = DD(C,dy,'y');

% Update Concentration
%C = C + dt * ( k*(Cxx+Cyy) - uX'.*Cx - uY'.*Cy );

C = C + dt * ( k*(Cxx+Cyy) - uX.*Cx - uY.*Cy );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Computes derivative based on sign of Velocity, u, using UPWIND
% approach
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function C_z = give_Necessary_Derivative(C,dz,uZ,string)

C_z = zeros(size(C));
len = length(uZ(:,1));
signs = sign(uZ);

if strcmp(string,'x')

    %For periodicity on ends w/ UPWIND
    for i=1:len

        %left side of grid
        if signs(i,1) <= 0
            C_z(i,1) =  ( C(i,2) - C(i,1) ) / (dz);
        else
            C_z(i,1) =  ( C(i,1) - C(i,len) ) / (dz);
        end

        %right side of grid
        if signs(i,len) <= 0
            C_z(i,len) =  ( C(i,1) - C(i,len) ) / (dz);
        else
            C_z(i,len) =  ( C(i,len) - C(i,len-1) ) / (dz);
        end

    end

    %Standard Upwind
    for i=1:len
        for j=2:len-1
            if signs(i,j) <= 0
                C_z(i,j) = ( C(i,j+1) - C(i,j) ) / (dz);
            else
                C_z(i,j) = ( C(i,j) - C(i,j-1) ) / (dz);
            end
        end
    end

    % Ends x-Direction calculation %

elseif strcmp(string,'y')

    %For periodicity on ends w/ UPWIND
    for i=1:len

        %bottom of grid
        if signs(1,i) <= 0
            C_z(1,i) =  ( C(2,i) - C(1,i) ) / (dz);
        else
            C_z(1,i) =  ( C(1,i) - C(len,i) ) / (dz);
        end

        %top of grid
        if signs(len,i) <= 0
            C_z(len,i) =  ( C(1,i) - C(len,i) ) / (dz);
        else
            C_z(len,i) =  ( C(len,i) - C(len-1,i) ) / (dz);
        end

    end

    %Standard Upwind
    for i=2:len-1
        for j=1:len
            if signs(i,j) <= 0
                C_z(i,j) = ( C(i+1,j) - C(i,j) ) / (dz);
            else
                C_z(i,j) = ( C(i,j) - C(i-1,j) ) / (dz);
            end
        end
    end

    % Ends y-Direction calculation %

else

    fprintf('\n\n\n ERROR IN FUNCTION D FOR COMPUTING 1ST DERIVATIVE\n');
    fprintf('Need to specify which desired derivative, x or y.\n\n\n');

end

clear signs; clear len;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Finds CENTERED finite difference approximation to 2ND
% DERIVATIVE in z direction, specified by input and 'string'
% Note: It automatically accounts for periodicity of the domain.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function u_zz = DD(u,dz,string)

% u:      velocity
% dz:     spatial step in "z"-direction
% string: specifies which 2ND derivative to take (to enforce periodicity)

len = length(u(:,1));

if strcmp(string,'x')

    %For periodicity on ends
    u_zz(:,1) =  ( u(:,2) - 2*u(:,1)   + u(:,len) )   / (dz^2);
    u_zz(:,len)= ( u(:,1) - 2*u(:,len) + u(:,len-1) ) / (dz^2);

    %Standard Upwind Scheme (Centered Difference)
    for j=2:len-1
        u_zz(:,j) = ( u(:,j+1) - 2*u(:,j) + u(:,j-1) ) / (dz^2);
    end

elseif strcmp(string,'y')

    %For periodicity on ends
    u_zz(1,:) =  ( u(2,:) - 2*u(1,:)   + u(len,:) )   / (dz^2);
    u_zz(len,:)= ( u(1,:) - 2*u(len,:) + u(len-1,:) ) / (dz^2);

    %Standard Upwind Scheme (Centered Difference)
    for j=2:len-1
        u_zz(j,:) = ( u(j+1,:) - 2*u(j,:) + u(j-1,:) ) / (dz^2);
    end

else

    fprintf('\n\n\n ERROR IN FUNCTION DD FOR COMPUTING 2ND DERIVATIVE\n');
    fprintf('Need to specify which desired derivative, x or y.\n\n\n');

end

author: nickabattista


text:
function [R, tau] = CholeskyMultIdentity(H)
%  Implements Cholesky with added multiple of the identity.  This attempts to
%  to find a scalar tau > 0 such that H + tau * I is sufficiently positive
%  definite, where I is the identity matrix.

global numFact  % Number of Cholesky factorizations attempted.

%  First, we intially try a Cholesky factorization.
[R, fail] = chol(H);
numFact = numFact + 1;
%  If it does not fail, we're done.
if fail == 0
    tau = 0;
    return;
end

%  If the initial Cholesky factorization fails, we attempt to find a scalar
%  tau > 0 such that H + tau * I is sufficiently positive definite.
beta = 0.001;  % Heuristic for increasing tau.
min_H_diag = min(diag(H));  % Smallest diagonal of H.

%  If smallest diagonal of H is positive, set tau to 0; otherwise, set to
%  nonnegative version of the smallest diagonal plus the beta heuristic.
if min_H_diag > 0
    tau = 0;
else
    tau = -min_H_diag + beta;
end

I = eye(size(H,1));  % Identity matrix.

%  Repeatedly add a tau-multiple of the identity to H until the Cholesky
%  factorization succeeds.  Upon each failure, double tau.
while 1
    [R, fail] = chol(H + tau * I);
    numFact = numFact + 1;
    if fail == 0
        return;
    else
        %  NOTE:  In order to decrease number of factorizations, we may want to
        %         increase tau by a factor of 10 instead of 2.
        tau = max(2*tau, beta);
    end
end
end

author: clarkzinzow

repo: Nonlinear-Optimization-Algorithms


text:
function ra = scimat_scimat2imref(scimat)
% SCIMAT_SCIMAT2IMREF  Matlab's image reference frame from SCIMAT metadata.
%
% SCIMAT_SCIMAT2IMREF produces an image frame (i.e. imref2d object) as used
% by Matlab to display images in real world coordinates. For example,
% Matlab's imshow() function uses imref2d objects as inputs.
%
% Note that imref2d objects do not support rotations, so the rotation
% matrix SCIMAT.rotmat will be ignored with a warning if it's not the
% identity.
%
% RA = SCIMAT_SCIMAT2IMREF(SCIMAT)
%
%   SCIMAT is a struct with a 2D image and its metadata (see "help
%   scimat"). The image can have several frames (2D+t), but it cannot be
%   3D. The image can also be multi-channel (e.g. RGB).
%
%   RA is the output imref2d object (see "help imref2d").
%
% Example:
%
% % create a SCIMAT image with some offset and non-isotropic pixel spacing
% offset = [2, 3];
% spacing = [0.3, 0.6];
% scimat = scimat_im2scimat(imread('cameraman.tif'), spacing, offset);
%
% % alternative 1: display it using scimat_imagesc
% subplot(2, 1, 1)
% hold off
% scimat_imagesc(scimat)
% axis equal
%
% % alternative 2: display it using a reference frame from scimat_scimat2imref
% ra = scimat_scimat2imref(scimat);
% subplot(2, 1, 2)
% hold off
% imshow(scimat.data, ra);
%
%
% See also: scimat, imref2d, scimat_imagesc, imshow.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2015 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(1, 1);
nargoutchk(0, 1);

% check input dimensions
if (size(scimat.data, 3) > 1)
    error('Implemented only for 2D images, not for 3D images')
end

% Matlab's image frame doesn't support rotation matrices
if (isfield(scimat, 'rotmat') && any(scimat.rotmat(:) ~= [1 0 0 1]'))
    warning('Matlab''s image frame doesn''t support rotation matrices. Ignoring rotation matrix')
end

% get image limits
ImageSize = [size(scimat.data, 1) size(scimat.data, 2)];
xmin = scimat_index2world([1, 1], scimat);
xmax = scimat_index2world(ImageSize, scimat);
XWorldLimits = [xmin(1) xmax(1)];
YWorldLimits = [xmin(2) xmax(2)];

ra = imref2d(ImageSize, XWorldLimits, YWorldLimits);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% OP_FDOTN_V: assemble the vector r = [r(i)], with  r(i) = (f * n, v), with n the normal exterior vector.
%  v should be a scalar-valued function.
%
%   rhs = op_fdotn_v (spv, msh, coeff);
%
% INPUT:
%
%   spv:   structure representing the function space (see sp_vector/sp_eval_boundary_side)
%   msh:   structure containing the domain partition and the quadrature rule (see msh_cartesian/msh_eval_boundary_side)
%   coeff: source function evaluated at the quadrature points
%
% OUTPUT:
%
%   rhs: assembled right-hand side
%
% Copyright (C) 2009, 2010 Carlo de Falco, Rafael Vazquez
% Copyright (C) 2011, 2014, 2017 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function rhs = op_fdotn_v (spv, msh, coeff)

  coeff = reshape (coeff, [], msh.nqn, msh.nel);
  shpv  = reshape (spv.shape_functions, spv.ncomp, msh.nqn, spv.nsh_max, msh.nel);

  rhs = zeros (spv.ndof, 1);

  for iel = 1:msh.nel
    if (all (msh.jacdet(:,iel)))
     jacdet_weights = reshape (msh.jacdet(:, iel) .* msh.quad_weights(:, iel), 1, msh.nqn);

     shpv_iel = reshape (shpv(:, :, :, iel), spv.ncomp, msh.nqn, spv.nsh_max);
     coeff_dot_n = sum (bsxfun (@times, coeff(:,:,iel), msh.normal(:,:,iel)), 1);

     coeff_times_jw = bsxfun (@times, jacdet_weights, coeff_dot_n);

     aux_val = bsxfun (@times, coeff_times_jw, shpv_iel);
     rhs_loc = sum (sum (aux_val, 1), 2);

     indices = find (spv.connectivity(:,iel));
     rhs_loc = rhs_loc(indices); conn_iel = spv.connectivity(indices,iel);
     rhs(conn_iel) = rhs(conn_iel) + rhs_loc(:);
    else
      warning ('geopdes:jacdet_zero_at_quad_node', 'op_fdotn_v: singular map in element number %d', iel)
    end
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
%TREXP Matrix exponential for so(3) and se(3)
%
% For so(3)::
%
% R = TREXP(OMEGA) is the matrix exponential (3x3) of the so(3) element OMEGA that
% yields a rotation matrix (3x3).
%
% R = TREXP(OMEGA, THETA) as above, but so(3) motion of THETA*OMEGA.
%
% R = TREXP(S, THETA) as above, but rotation of THETA about the unit vector S.
%
% R = TREXP(W) as above, but the so(3) value is expressed as a vector W
% (1x3) where W = S * THETA. Rotation by ||W|| about the vector W.
%
% For se(3)::
%
% T = TREXP(SIGMA) is the matrix exponential (4x4) of the se(3) element SIGMA that
% yields a homogeneous transformation  matrix (4x4).
%
% T = TREXP(SIGMA, THETA) as above, but se(3) motion of SIGMA*THETA, the
% rotation part of SIGMA (4x4) must be unit norm.
%
% T = TREXP(TW) as above, but the se(3) value is expressed as a twist vector TW
% (1x6).
%
% T = TREXP(TW, THETA) as above, but se(3) motion of TW*THETA, the
% rotation part of TW (1x6) must be unit norm.
%
% Notes::
% - Efficient closed-form solution of the matrix exponential for arguments that are
%   so(3) or se(3).
% - If THETA is given then the first argument must be a unit vector or a
%   skew-symmetric matrix from a unit vector.
% - Angle vector argument order is different to ANGVEC2R.
%
% References::
% - Robotics, Vision & Control: Second Edition, P. Corke, Springer 2016; p42-43.
% - Mechanics, planning and control, Park & Lynch, Cambridge, 2017.
%
% See also ANGVEC2R, TRLOG, TREXP2, SKEW, SKEWA, Twist.

%## 3d homogeneous differential

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math
function T = trexp(S, theta)

    if ishomog(S) || isvec(S,6)
        % input is se(3)

        if nargin == 1
            % twist vector 1x6 or augmented skew matrix 4x4
            if isvec(S,6)
                % it's a twist vector
                S = skewa(S);
            end
            T = expm(S);
        else
            % se(3) plus twist
            if all(size(S) == 4)
                % it's se(3) matrix
                [skw,v] = tr2rt(S);
            else
                % it's a twist vector
                v = S(1:3); v= v(:);
                skw = skew(S(4:6));
            end

            % use an efficient solution
            R = trexp(skw, theta);
            t = (eye(3,3)*theta + (1-cos(theta))*skw + (theta-sin(theta))*skw^2)*v;

            T = rt2tr(R,t);
        end
    elseif isrot(S) || isvec(S,3)
        % input is so(3)

        if isrot(S)
            % input is 3x3 skew symmetric
            w = vex(S);
        elseif isvec(S)
            % input is a 3-vector
             w = S;
        end

        % for a zero so(3) return unit matrix, theta not relevant
        if norm(w) < 10*eps
            T = eye(3,3);
            return;
        end

        if nargin == 1
            %  theta is not given, extract it
            theta = norm(w);
            w = unit(w);
        else
            % theta is given
            assert(isunit(w), 'SMTB:trexp:badarg',  'w must be a unit twist');
        end

        S = skew(w);

        T = eye(3,3) + sin(theta)*S + (1-cos(theta))*S^2;

    else
        error('SMTB:trexp:badarg', 'first argument must be so(3), 3-vector, se(3) or 6-vector');
    end
end

author: petercorke

repo: spatialmath-matlab


text:
function [fSignificanceLevel, fMu, fSigma] = kj_CalcNormSig(mValues, fTestValue)
% function [fSignificanceLevel, fMu, fSigma] = kj_CalcNormSig(mValues, fTestValue)
% --------------------------------------------------------------------------------
% Calculates the level of significance of fTestValue assuming the distribution
% to be a normal distribution
%
% Input parameters:
%   mValues               Value distribution (assumed to be a normal distribution)
%   fTestValue            Value to be tested
%
% Output parameter:
%   fSignificanceLevel    Level of significance
%   fMu                   Mu of normal distribution
%   fSigma                Sigma of normal distribution
%
% Danijel Schorlemmer
% March 13, 2002

global bDebug;
if bDebug
  report_this_filefun(mfilename('fullpath'));
end

% Select all non-NaN values of the distribution
vSelection = ~isnan(mValues);
mNoNanValues = mValues(vSelection);

% Fit the values to a normal distribution
[fMu, fSigma] = normfit(mNoNanValues);

% Return the significance level of the testvalue
fSignificanceLevel = 1 - (normcdf(fTestValue, fMu, fSigma));

author: CelsoReyes


text:
function [A, C, Z] = ldsPca(X, k, m)
% Subspace method for learning linear dynamic system.
% Input:
%   X: d x n data matrix
%   k: dimension of hidden variable
%   m: stacking order for the Hankel matrix
% Output:
%   A: k x k transition matrix
%   C: k x d emission matrix
%   Z: k x n latent variable
%   Y: d x n reconstructed data
% reference: Bayesian Reasoning and Machine Learning (BRML) chapter 24.5.3 p.507
% Written by Mo Chen (sth4nth@gmail.com).
[d,n] = size(X);
H = reshape(X(:,hankel(1:m,m:n)),d*m,[]);
[U,S,V] = svd(H,'econ');
C = U(1:d,1:k);
Z = S(1:k,1:k)*V(:,1:k)';
A = Z(:,2:end)/Z(:,1:end-1); % estimated transition
% Y = C*Z; % reconstructions

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function test_suite = testSetupExample
initTestSuite;

function fh = setup
fh = figure;

function teardown(fh)
delete(fh);

function testColormapColumns(fh)
assertEqual(size(get(fh, 'Colormap'), 2), 3);

function testPointer(fh)
assertEqual(get(fh, 'Pointer'), 'arrow');

author: ricedsp


text:
function h=ref_pconv(f,g,ctype)
%REF_PCONV  Reference PCONV
%   Usage:  h=ref_pconv(f,g)
%
%   PCONV(f,g) computes the periodic convolution of f and g.

% AUTHOR: Peter L. Søndergaard

L=length(f);
h=zeros(L,1);

switch(lower(ctype))
  case {'default'}
    for ii=0:L-1
      for jj=0:L-1
	h(ii+1)=h(ii+1)+f(jj+1)*g(mod(ii-jj,L)+1);
      end;
    end;
  case {'r'}
    for ii=0:L-1
      for jj=0:L-1
	h(ii+1)=h(ii+1)+f(jj+1)*conj(g(mod(jj-ii,L)+1));
      end;
    end;
  case {'rr'}
    for ii=0:L-1
      for jj=0:L-1
	h(ii+1)=h(ii+1)+conj(f(mod(-jj,L)+1))*conj(g(mod(jj-ii,L)+1));
      end;
    end;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [beta_gibbs sigma_gibbs]=mgibbs(It,Bu,betabar,omegabar,sigma,q)

% function [beta_gibbs sigma_gibbs]=mgibbs(It,Bu,betabar,omegabar,sigma,q)
% performs the Gibbs algorithm 2.1.2 for the Minnesota prior, and returns draws from posterior distribution
% inputs:  - integer 'It': total number of iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'Bu': number of burn-in iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - vector 'betabar': posterior mean vector (defined in 1.3.18)
%          - matrix 'omegabar': posterior covariance matrix for the VAR coefficients (defined in 1.3.17)
%          - matrix 'sigma': 'true' variance-covariance matrix of VAR residuals, for the original Minnesota prior
%          - integer 'q': total number of coefficients to estimate for the BVAR model (defined p 7 of technical guide)
% outputs: - matrix 'beta_gibbs': record of the gibbs sampler draws for the beta vector
%          - matrix'sigma_gibbs': record of the gibbs sampler draws for the sigma matrix (vectorised)

% start iterations
for ii=1:(It-Bu)
% draw a random vector beta from N(betabar,sigmabar) (step 3 of the algorithm)
beta=betabar+chol(bear.nspd(omegabar),'lower')*randn(q,1);

% record values before starting next iteration
beta_gibbs(:,ii)=beta;
% go for next iteration
end

% record the values for the variance-covariance matrix sigma
% vectorize sigma and record it (it-Bu) times
sigma_gibbs=repmat(sigma(:),1,It-Bu);

author: european-central-bank

repo: BEAR-toolbox


text:
% Implement a simple version of "Geodesics in Heat" by Crane et al. 2013

clc; clear all; close all;
[V,F] = readOBJ('../data/spot.obj');
[V,F,S] = loop(V,F,2); % upsample the mesh to see the speed difference
use_prefactorization = true; % set it to "true" or "false" to see speed difference

% precompute matrices
A = massmatrix(V,F);
Lc = cotmatrix(V,F);
t = avgedge(V,F)^2;
LHS = A - t*Lc;
G = grad(V,F); % #F*dim by #V
D = div(V,F); % #V by #F*dim

% prefactorization
preLHS = decomposition(LHS);
preLc = decomposition(Lc);

%% change heat source
heatSrcIdx = 1;

tic;

% step 1
delta = zeros(size(V,1),1);
delta(heatSrcIdx) = 1;
if use_prefactorization
    u = preLHS \ delta;
else
    u = LHS \ delta;
end

% step 2
gradu = G*u;
gradu = reshape(gradu, size(F,1), 3);
gradu_normalized = gradu ./ normrow(gradu);
X = -gradu_normalized;

% step 3
divX = D * reshape(X, size(F,1)*3, 1);
if use_prefactorization
    phi = preLc \ divX;
else
    phi = Lc \ divX;
end
phi = phi - phi(heatSrcIdx);

toc;

%% visualization
tsurf(F,V, 'CData', phi);
axis equal;
CM = cbrewer('Reds', 500);
CM = CM(size(CM,1):-1:1,:);
colormap(CM);
colorbar
shading interp
drawnow

author: odedstein

repo: sgi-introduction-course


text:
function SR = superresolve_dbrsr(slidingWindows, magFactor)

numberOfFrames = size(slidingWindows.frames,3);

% Set PSF kernels for the different binning factors.
psf{1} = fspecial('gaussian', 5, 0.8); % PSF for binning factor 2
psf{2} = fspecial('gaussian', 7, 1.2); % PSF for binning factor 3
psf{3} = fspecial('gaussian', 9, 1.6); % PSF for binning factor 4

% Specific preprocessing steps for the non-uniform interpolation based multi-frame methods.
weighting_struct = struct('mec_weight_flag',4,'dis_weight_flag',1,'qps_weight_flag',0,'mec_confMapN',[],...
    'mec_thr_scaler',2,'mec_weight_exp',2,'mec_weighting',[],'dis_rho',0.7,'dis_scaler',magFactor*10,'dis_weighting',[],...
    'qps_map',[],'qps_rho',0.7,'qps_weighting',[],'alpha',1,'beta',1,'gamma',1,'whi_mode',2);
warped_meshXN = zeros(size(slidingWindows.referenceFrame,1),size(slidingWindows.referenceFrame,2),numberOfFrames-1);
warped_meshYN = zeros(size(slidingWindows.referenceFrame,1),size(slidingWindows.referenceFrame,2),numberOfFrames-1);

lr_seq_sorted = slidingWindows.frames;
lr_seq_tmp = lr_seq_sorted(:,:,(numberOfFrames+1)/2);
lr_seq_sorted(:,:,(numberOfFrames+1)/2) = lr_seq_sorted(:,:,1);
lr_seq_sorted(:,:,1) = lr_seq_tmp;

lr_seq_sorted = permute(lr_seq_sorted,[1 2 4 3]);

permutation_vec = 1:numberOfFrames-1;
permutation_vec((numberOfFrames+1)/2:end) = permutation_vec((numberOfFrames+1)/2:end) + 1;
permutation_vec(1:(numberOfFrames-1)/2)   = circshift(permutation_vec(1:(numberOfFrames-1)/2),[0 -1]);

fw_optFlowN_xc = zeros(size(lr_seq_sorted,1),size(lr_seq_sorted,2),2,numberOfFrames-1);

for itera = 1:numberOfFrames-1
    [warped_meshXN(:,:,itera),warped_meshYN(:,:,itera)] = lmsSR_warpImgGridLocalTrans(size(slidingWindows.referenceFrame,1),size(slidingWindows.referenceFrame,2),slidingWindows.flowToReference{permutation_vec(itera)}.mvs_xc);

    % Storing the MEC-SSD-Confidence into the Weighting Struct
    weighting_struct.mec_confMapN(:,:,itera) = slidingWindows.mec_confMap{permutation_vec(itera)};

    % Forward MVF for WDBR
    fw_optFlowN_xc(:,:,:,itera) = slidingWindows.flowToReference{permutation_vec(itera)}.mvs_xc;
end

% Actual SR call
[~,wnuisr_nodeconv,~,~] = lmsSR_generateSRusingNUIv4weighted(lr_seq_sorted,warped_meshXN,warped_meshYN,[magFactor, magFactor],psf{magFactor-1},2,weighting_struct);
[SR,~] = lmsSR_generateSRusingWDBR(wnuisr_nodeconv,fw_optFlowN_xc,magFactor,numberOfFrames,psf{magFactor-1},2,weighting_struct);

author: thomas-koehler

repo: SupER


text:
function f = p02_f ( n, x )

%*****************************************************************************80
%
%% P02_F evaluates the objective function for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the argument of the objective function.
%
%    Output, real F, the value of the objective function.
%
  f = 0.0;

  for i = 1 : 13

    c = - i / 10.0;

    fi = x(3)     * exp ( c * x(1) )         - x(4) * exp ( c * x(2) ) ...
       + x(6)     * exp ( c * x(5) )         -        exp ( c ) ...
       + 5.0  * exp ( 10.0 * c ) - 3.0  * exp ( 4.0 * c );

    f = f + fi * fi;

  end

  return
end


text:
function linplus_test485 ( )

%*****************************************************************************80
%
%% TEST485 tests R8PBU_INDICATOR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  n = 9;
  mu = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST485\n' );
  fprintf ( 1, '  R8PBU_INDICATOR sets up a R8PBU indicator matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
  fprintf ( 1, '  Bandwidth MU = %d\n', mu );

  a = r8pbu_indicator ( n, mu );

  r8pbu_print ( n, mu, a, '  The R8PBU indicator matrix:' );

  return
end


text:
   function y = tang(x,y)
    % to find the exact angle in degrees
           y1 = atan2(y,x);
             if y1 >= 0
               y = y1*180.0/pi;
                 elseif  y1 < 0
                    y = y1*180.0/pi + 360.0;
                end;

author: Sable

repo: mcbench-benchmarks


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function optval = LookBackFixedStrike(S,K,C)
% S = NSim x Nt matrix of simulated prices
% K = Strike price
% C = 1 -> Call; C = 0 -> Put

    if(C == 1)
        optval = mean(max(S,2) - K, 0);
    else
        optval = mean(K - min(S,2) - K, 0);
    end
end


text:
% LTFAT - Non-stationary Gabor systems
%
%  Florent Jaillet and Peter L. Søndergaard, 2011 - 2018
%
%  Transforms
%    NSDGT                - Non-stationary DGT
%    UNSDGT               - Uniform non-stationary DGT
%    INSDGT               - Inverse NSDGT and UNSDGT
%    NSDGTREAL            - Non-stationary DGT for real-valued signals
%    UNSDGTREAL           - Uniform non-stationary DGT for real-valued signals
%    INSDGTREAL           - Inverse NSDGTREAL and UNSDGTREAL
%
%  Window construction and bounds
%    NSGABDUAL            - Non-stationary dual windows
%    NSGABTIGHT           - Non-stationary tight windows
%    NSGABFRAMEBOUNDS     - Frame bounds of an NSDGT system
%    NSGABFRAMEDIAG       - Diagonal of non-stationary Gabor frame operator
%
%  Plots
%    PLOTNSDGT            - Plot output coefficients from NSDGT
%    PLOTNSDGTREAL        - Plot output coefficients from NSDGTREAL
%
%  For help, bug reports, suggestions etc. please visit
%  http://github.com/ltfat/ltfat/issues


text:
function [plaza,v,vmax]=random_slow(plaza,v,vmax,probslow);
    [L,W]=size(plaza);
    for lanes=2:W-1;
        temp=find(plaza(:,lanes)==1);
        nn=length(temp);
        for k=1:nn;
            i=temp(k);
            if(rand<=probslow)
                v(i,lanes)=max(v(i,lanes)-1,0);
            end
        end
    end
end

author: ravenxrz

repo: Mathematical-Modeling


text:
function euler = rotMat2euler(R)
%ROTMAT2EULER Converts a rotation matrix orientation to ZYX Euler angles
%
%   euler = rotMat2euler(R)
%
%   Converts a rotation matrix orientation to ZYX Euler angles where phi is
%   a rotation around X, theta around Y and psi around Z.
%
%   For more information see:
%   http://www.x-io.co.uk/node/8#quaternions
%
%   Date          Author          Notes
%   27/09/2011    SOH Madgwick    Initial release

    phi = atan2(R(3,2,:), R(3,3,:) );
    theta = -atan(R(3,1,:) ./ sqrt(1-R(3,1,:).^2) );
    psi = atan2(R(2,1,:), R(1,1,:) );

    euler = [phi(1,:)' theta(1,:)' psi(1,:)'];
end

author: xioTechnologies

repo: Oscillatory-Motion-Tracking-With-x-IMU


text:
function theta = circle_segment_angle_from_chord_angles ( omega1, omega2 )

%*****************************************************************************80
%
%% CIRCLE_SEGMENT_ANGLE_FROM_CHORD_ANGLES computes the angle of a circle segment.
%
%  Discussion:
%
%    Begin with a circle of radius R.  Choose two points P1 and P2 on the
%    circle, and draw the chord P1:P2.  This chord divides the circle
%    into two pieces, each of which is called a circle segment.
%    Consider one of the pieces.  The "angle" of this segment is the angle
%    P1:C:P2, where C is the center of the circle.  Let Q be the point on
%    the chord P1:P2 which is closest to C.  The "height" of the segment
%    is the distance from Q to the perimeter of the circle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real OMEGA1, OMEGA2, the angles of the points P1 and P2.
%    OMEGA1 <= OMEGA2.
%
%    Output, real THETA, the angle of the circle segment.
%    Essentially, THETA = OMEGA2 - OMEGA1.
%
  while ( omega2 < omega1 )
    omega2 = omega2 + 2.0 * pi;
  end

  theta = omega2 - omega1;

  return
end


text:
function i = end(x,k,n)
%END          Overloaded functions end, specifies last index
%

% written  10/03/02     S.M. Rump
%

  if x.complex
    xx = x.mid;
  else
    xx = x.inf;
  end
  if n==1           % call as one-dimensional array
    i = length(xx(:));
  else
    i = size(xx,k);
  end

repo: OpenMAS


text:
function cdf = gp_predcdf(gp, x, y, varargin)
%GP_PREDCDF  Predictive CDF evaluated at YT
%
%  Description
%    CDF = GP_PREDCDF(GP, X, Y, XT, 'yt', YT, OPTIONS)
%    takes a GP structure together with matrix X of training
%    inputs and vector Y of training targets, and evaluates the
%    cdf of the predictive distribution at test inputs XT, YT.
%
%    CDF = GP_PREDCDF(GP, X, Y, OPTIONS) evaluates the
%    cdf of the predictive distribution at training inputs X, Y.
%
%    OPTIONS is optional parameter-value pair
%      tstind - a vector defining, which rows of X belong to which
%               training block in *IC type sparse models. Default is [].
%               See also GP_PRED.
%      z      - optional observed quantity in triplet (x_i,y_i,z_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, expected value
%               for ith case.
%      zt     - optional observed quantity in triplet (xt_i,yt_i,zt_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, the expected
%               value for the ith case.
%
%  See also
%    GP_PRED, GP_PAK, GP_UNPAK
%
% Copyright (c) 2012 Aki Vehtari

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

  ip=inputParser;
  ip.FunctionName = 'GP_PREDCDF';
  ip.addRequired('gp',@(x) isstruct(x) || iscell(x));
  ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addOptional('xt', [], @(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))))
  ip.addParamValue('yt', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('zt', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('prct', [5 50 95], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('nsamp', 5000, @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('predcf', [], @(x) isempty(x) || ...
                   isvector(x) && isreal(x) && all(isfinite(x)&x>0))
  ip.addParamValue('tstind', [], @(x) isempty(x) || iscell(x) ||...
                   (isvector(x) && isreal(x) && all(isfinite(x)&x>0)))
  if numel(varargin)==0 || isnumeric(varargin{1})
    % inputParser should handle this, but it doesn't
    ip.parse(gp, x, y, varargin{:});
  else
    ip.parse(gp, x, y, [], varargin{:});
  end
  xt=ip.Results.xt;
  yt=ip.Results.yt;
  z = ip.Results.z;
  zt = ip.Results.zt;
  prct = ip.Results.prct;
  nsamp = ip.Results.nsamp;
  predcf=ip.Results.predcf;
  tstind=ip.Results.tstind;
  if isempty(xt)
    xt=x;
    if isempty(tstind)
      if iscell(gp)
        gptype=gp{1}.type;
      else
        gptype=gp.type;
      end
      switch gptype
        case {'FULL' 'VAR' 'DTC' 'SOR'}
          tstind = [];
        case {'FIC' 'CS+FIC'}
          tstind = 1:size(x,1);
        case 'PIC'
          if iscell(gp)
            tstind = gp{1}.tr_index;
          else
            tstind = gp.tr_index;
          end
      end
    end
    if isempty(yt)
      yt=y;
    end
    if isempty(zt)
      zt=z;
    end
  end

  % pass these forward
  options=struct();
  if ~isempty(z);options.z=z;end
  if ~isempty(yt);options.yt=yt;end
  if ~isempty(zt);options.zt=zt;end
  if ~isempty(predcf);options.predcf=predcf;end
  if ~isempty(tstind);options.tstind=tstind;end

  [tn, nin] = size(x);

  if iscell(gp) || numel(gp.jitterSigma2)>1 || isfield(gp,'latent_method')
    % gp_array
    if iscell(gp)
      nGP = numel(gp);
      for i1=1:nGP
        Gp=gp{i1};
        P_TH(:,i1)=Gp.ia_weight;
        [Ef,Varf]=gp_pred(Gp,x,y,xt,options);
        cdfs(:,i1)=Gp.lik.fh.predcdf(Gp.lik, Ef, Varf, yt);
      end
      cdf=sum(bsxfun(@times,cdfs,P_TH),2);
    elseif numel(gp.jitterSigma2)>1
      % MCMC samples
      [Efs,Varfs]=gpmc_preds(gp,x,y,xt,options);
      nmc=size(gp.jitterSigma2,1);
      for i1=1:nmc
        Gp = take_nth(gp,i1);
        cdfs(:,i1)=Gp.lik.fh.predcdf(Gp.lik, Ef, Varf, yt);
      end
      cdf=mean(cdfs, 2);
    else
      [Ef,Varf]=gp_pred(gp,x,y,xt,options);
      cdf=gp.lik.fh.predcdf(gp.lik, Ef, Varf, yt);
    end

  end

end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function TP = getTransProbs (hmm)
% Returns the transition probabilies from any state to any other state,
% without considering the persistence probabilities (i.e. the probability
% to remain in the same state)
%
% Author: Diego Vidaurre, University of Oxford (2017)

TP = hmm.P;
[K,~,Q] = size(TP);
for j = 1:Q
   TPj = TP(:,:,j);
   TPj(eye(K)==1) = 0;
   for k = 1:K
      TPj(k,:) = TPj(k,:) / sum(TPj(k,:));
   end
   TP(:,:,j) = TPj;
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function [azStart,dist,azEnd,latLonWaypoints]=indirectGreatCircleProb(latLonStart,latLonEnd,r,N,algorithm,waypointType,noJumps)
%%INDIRECTGREATCIRCLEPROB Solve the indirect great circle problem. That is,
%                      given two points on a spherical Earth, find the
%                      initial bearing and distance one must travel to
%                      take the shortest (geodesic) path between the
%                      points. Additionaly, this function can compute
%                      waypoints on the path. This type of navigation is
%                      sometimes referred to as "great circle sailing."
%                      When used for navigation, it is less accurate than
%                      an ellipsoidal Earth approximation.
%
%INPUTS: latLonStart The 2XnumPts set of numPts initial points given in
%               latitude and longitude in radians in the format
%               [latitude;longitude] (on a reference sphere, latitude is
%               spherical elevation; longitude is spherical azimuth). If
%               all these points are the same but latLonEnd varies, then a
%               single 2X1 vector can be passed. Extra rows, if passed,
%               will be ignored.
%     latLonEnd The 2XnumPts final points given in latitude and longitude
%               in adians in the format [latitude;longitude]. If all these
%               points are the same but latLonStart varies, then a single
%               2X1 vector can be passed.
%             r The assumed radius of the spherical Earth model. If omitted
%               or an empty matrix is passed, the default of
%               r=osculatingSpher4LatLon(latLonStart) is used.
%             N The number of waypoints, besides the initial and final
%               points on the trajectory, to produce. If omitted or an
%               empty matrix is passed, then no waypoints are generated and
%               latLonWaypoints is an empty matrix.
%     algorithm An optional parameter selecting the algorithm to use.
%               Possible values are:
%               0 Use the algorithm of [2], which avoids issues with
%                 singularities at the poles.
%               1 (The default if omitted or an empty matrix is passed and
%                 waypointType=1) Use the COFI algorithm of [1]. Latitudes
%                 within 2^24*eps(pi/2) of +/-pi/2 (the North and South
%                 poles) will be clipped to that bound. This avoids a
%                 singularity at the poles.
%               2 (The default if omitted or an empty matrix is passed and
%                 waypointType=0) Obtain angles using Equation 5-4b of [3],
%                 which is stable at the poles. Obtaincthe distance using a
%                 cross product formula, which is also used in algorithm 1,
%                 as discussed below.
%    waypointType An optional parameter indicating how the waypoints will
%                 be spaced. Possible values are:
%                 0 (The default if omitted or an empty matrix is passed)
%                   Space the N waypoints between the start and end
%                   uniformly in distance.
%                 1 Space the N waypoints between the start and end
%                   uniformly in longitude. This option is only available
%                   for algorithm 0. If a trajectory is meridional, then
%                   only the starting and stopping points will be returned.
%
%OUTPUTS: azStart The 1XnumPts scalar forward azimuths at the starting
%                 point in radians East of true North on the reference
%                 sphere.
%            dist The 1XnumPts distances on the sphere-approximated Earth
%                 between the starting and stopping points.
%           azEnd The 1XnumPts forward azimuth at the ending point in
%                 radians East of true North on the reference sphere.
% latLonWaypoints A 2X(N+2) set of waypoints on the trajectory.
%                 latLonWaypoints(:,i) is the ith points as
%                 [latitude;longitude] in radians. The first point is
%                 latLonStart and the final points is latLonEnd, though the
%                 points are made to avoid jumps of 2*pi in longitude,
%                 which can mean that the last point in this can be off
%                 from the provided last point by a factor of 2*pi. The
%                 other points are equally spaced along the trajectory. If
%                 the trajectory is meridional and waypointType=1, then
%                 only the starting and ending points will be returned.
%
%Great circles are geodesics on a sphere. A geodesic path between two
%points on a sphere is much simpler to determine as compared to an
%ellipsoid (which is what one gets using the indirectGeodeticProb
%function). For the spherical model, the latitude is taken as the elevation
%and the longitude as the azimuth on the sphere (we assume coordinate
%system 0 in the spher2Cart function). The latitudes and longitudes
%returned as waypoints by this function can be used to piece together an
%_approximate_ geodesic trajectory on the reference ellipsoid. One can
%obtain headings and distances between the reference points using rhumb
%lines.
%
%The azimuthal values returned when one of the points is at the pole
%depends on the longitude value given for the points. The getENUAxes
%function can return East-North-Up axes anywhere on the globe (for a
%sphere, the ellipsoidal flattening is 0). When given a pole, the result is
%the equivalent one would get by taking a point with the same latitude and
%with the latitude magnitude a tiny epsilon less than pi/2.
%
%The implementation of the COFI algorithm of [1] has been modified so that
%the angular difference D is obtained using the angBetweenVecs given F and
%T rather than the technique derived from the dot product relation starting
%in Equation 13 in [1]. This is because the cross product relation used in
%angBetweenVecs is more accurate. This is also used for the distance in
%algorithm 2. Minor changes had to be made to algorithm 1 to deal with
%trajectories going West, since the derivation in [1] is for East-bound
%trajectories.
%
%EXAMPLE 1:
%This is example 2 from [1]. It is a trajectory that crosses the
%international date and and goes from the Northern hemisphere to the
%southern hemisphere. We also compute the reverse path and show that the
%start and end azimuth angles produced in each direction are consistent
%with each other (when the same Earth radius is used eahc time). We then
%plot the trajectory on an image of the spherical Earth. For better
%plotting, the radius of the Earth has been normalized to 1.
% N=500;
% latStart=degMinSec2Rad(37,47.5);
% lonStart=degMinSec2Rad(-122,-27.8);
% latEnd=degMinSec2Rad(-33,-51.7);
% lonEnd=degMinSec2Rad(151,12.7);
%
% latLonStart=[latStart;lonStart];
% latLonEnd=[latEnd;lonEnd];
%
% [azStartFwd,distFwd,azEndFwd,latLonWayPoints]=indirectGreatCircleProb(latLonStart,latLonEnd,1,N);
% [azStartRev,distRev,azEndRev,latLonWayPointsRev]=indirectGreatCircleProb(latLonEnd,latLonStart,1,N);
%
% %All four of these should be approximately zero if the algorithm works
% %forwards and backwards:
% distFwd-distRev
% wrapRange(azStartFwd-wrapRange(azEndRev+pi,-pi,pi),-pi,pi)
% wrapRange(azStartRev-wrapRange(azEndFwd+pi,-pi,pi),-pi,pi)
% max(max(abs(wrapRange(latLonWayPoints-fliplr(latLonWayPointsRev),-pi,pi))))
%
% xStartCart=ellips2Cart([latLonStart;0],1,0);
% xEndCart=ellips2Cart([latLonEnd;0],1,0);
% %Give the points a nonzero elevation so the lines are easier to see.
% pathPoints=ellips2Cart([latLonWayPoints;0.02*ones(1,N+2)],1,0);
%
% figure(1)
% clf
% hold on
% plotMapOnEllipsoid([],1,0);
% scatter3(xStartCart(1),xStartCart(2),xStartCart(3),100,'filled')
% scatter3(xEndCart(1),xEndCart(2),xEndCart(3),100,'filled')
% plot3(pathPoints(1,:),pathPoints(2,:),pathPoints(3,:),'-r','linewidth',4)
% view([-75.368365161957513,3.494544111878512])
%
%EXAMPLE 2:
%This is the same as example 1, except one point is given at the North
%pole.
% N=500;
% latStart=pi/2;%North pole.
% lonStart=degMinSec2Rad(-122,-27.8);
% latEnd=degMinSec2Rad(-33,-51.7);
% lonEnd=degMinSec2Rad(151,12.7);
%
% latLonStart=[latStart;lonStart];
% latLonEnd=[latEnd;lonEnd];
%
% [azStartFwd,distFwd,azEndFwd,latLonWayPoints]=indirectGreatCircleProb(latLonStart,latLonEnd,1,N);
% [azStartRev,distRev,azEndRev,latLonWayPointsRev]=indirectGreatCircleProb(latLonEnd,latLonStart,1,N);
%
% %All four of these should be approximately zero if the algorithm works
% %forwards and backwards:
% distFwd-distRev
% wrapRange(azStartFwd-wrapRange(azEndRev+pi,-pi,pi),-pi,pi)
% wrapRange(azStartRev-wrapRange(azEndFwd+pi,-pi,pi),-pi,pi)
% max(max(abs(wrapRange(latLonWayPoints-fliplr(latLonWayPointsRev),-pi,pi))))
%
% xStartCart=ellips2Cart([latLonStart;0],1,0);
% xEndCart=ellips2Cart([latLonEnd;0],1,0);
% pathPoints=ellips2Cart([latLonWayPoints;0.02*ones(1,N+2)],1,0);
%
% figure(1)
% clf
% hold on
% plotMapOnEllipsoid([],1,0);
% scatter3(xStartCart(1),xStartCart(2),xStartCart(3),100,'filled')
% scatter3(xEndCart(1),xEndCart(2),xEndCart(3),100,'filled')
% plot3(pathPoints(1,:),pathPoints(2,:),pathPoints(3,:),'-r','linewidth',4)
%
%REFERENCES:
%[1] C.-L. Chen, P.-F. Liu, and W.-T. Gong, "A simple approach to great
%    circle sailing: The COFI method," The Journal of Navigation, vol. 67,
%    no. 3, pp. 403-418, May 2014.
%[2] D. F. Crouse, "Singularity-free great-circle sailing," Naval Research
%    Laboratory, Washington, DC, Tech. Rep. NRL/5340/MR-2021/4, 26
%    Jul. 2021.
%[3] J. P. Snyder, "Map projections-a working manual," U.S. Geological
%    Survey, Tech. Rep. 1395, 1987.
%
%August 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<7||isempty(noJumps))
    noJumps=true;
end

if(nargin<6||isempty(waypointType))
    waypointType=0;
end

if(nargin<5||isempty(algorithm))
    if(waypointType==1)
        algorithm=1;
    else
        algorithm=2;
    end
end

if(nargin<4)
    N=[];%Generate no waypoints.
end

if(nargin<3||isempty(r))
    r=osculatingSpher4LatLon(latLonStart);
end

if(nargout<3)
    N=[];%Do not generate waypoints.
end

numPts1=size(latLonStart,2);
numPts2=size(latLonEnd,2);
numPts=max(numPts1,numPts2);

if(numPts1>numPts2)
    latLonEnd=repmat(latLonEnd,[1,numPts]);
elseif(numPts2>numPts1)
    latLonStart=repmat(latLonStart,[1,numPts]);
end

azStart=zeros(1,numPts);
dist=zeros(1,numPts);
azEnd=zeros(1,numPts);
if(~isempty(N))
    latLonWaypoints=zeros(2,N+2,numPts);
else
    latLonWaypoints=[];
end

switch(algorithm)
    case 0
        if(waypointType~=0)
            error('This algorithm only supports waypointType=0.');
        end

        for k=1:numPts
            [azStart(k),dist(k),azEnd(k),latLonWaypointsCur]=indirectGreatCircleProbCrouse(latLonStart(1:2,k),latLonEnd(1:2,k),r,N);
            if(~isempty(N))
                latLonWaypoints(:,:,k)=latLonWaypointsCur;
            end
        end
    case 1
        for k=1:numPts
            [azStart(k),dist(k),azEnd(k),latLonWaypointsCur]=indirectGreatCircleProbChen(latLonStart(1:2,k),latLonEnd(1:2,k),r,N,waypointType);
            if(~isempty(N))
                latLonWaypoints(:,:,k)=latLonWaypointsCur;
            end
        end
    case 2
        for k=1:numPts
            [azStart(k),dist(k),azEnd(k),latLonWaypointsCur]=indirectGreatCircleProbSnyder(latLonStart(1:2,k),latLonEnd(1:2,k),r,N);
            if(~isempty(N))
                latLonWaypoints(:,:,k)=latLonWaypointsCur;
            end
        end
    otherwise
        error('Unknown algorithm chosen.')
end

if(~isempty(N)&&noJumps&&N>0)
    cumOffset=0;
    for k=1:numPts
        for j=2:(N+2)
            diffVal=latLonWaypoints(2,j)+cumOffset-latLonWaypoints(2,j-1);
            if(diffVal>pi)
                cumOffset=cumOffset-2*pi;
            elseif(diffVal<-pi)
                cumOffset=cumOffset+2*pi;
            end
            latLonWaypoints(2,j)=latLonWaypoints(2,j)+cumOffset;
        end
    end
end
end

function [azStart,dist,azEnd,latLonWayPoints]=indirectGreatCircleProbChen(latLonStart,latLonEnd,r,N,waypointType)
%%INDIRECTGREATCIRCLEPROBCHEN Solve the indirect great circle problem using
%                   the algorithm of [1], except the angular difference D
%                   is obtained using the angBetweenVecs given F and T.
%
%Starting latitude values within 2^24*eps(pi/2); of the pole will be
%clipped to be 2^24*eps(pi/2); away from the pole. This reduces problems
%with singularities at the poles at the cost of a loss of accuracy compared
%to elsewhere on the globe.
%
%REFERENCES:
%[1] C.-L. Chen, P.-F. Liu, and W.-T. Gong, "A simple approach to great
%    circle sailing: The COFI method," The Journal of Navigation, vol. 67,
%    no. 3, pp. 403-418, May 2014.
%
%August 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Latitude at the start.
LF=latLonStart(1);
%Deal with locations very close to the North or South pole by offsetting
%them just enough to avoid finite precision issues.
epsVal=2^24*eps(pi/2);
if(abs(LF-pi/2)<epsVal)
    LF=pi/2-epsVal;
elseif(abs(LF+pi/2)<epsVal)
    LF=pi/2+epsVal;
end

%Latitude at the destination.
LT=latLonEnd(1);

%Difference in longitude. If DLo is positive, then the shortest path is
%going East. If it is negative, then the shortest path is going west.
DLo=wrapRange(latLonEnd(2)-latLonStart(2),-pi,pi);

goingWest=(DLo<0);

sinLF=sin(LF);
cosLF=cos(LF);
sinLT=sin(LT);
cosLT=cos(LT);

lonStart=latLonStart(2);
cosLS=cos(lonStart);
sinLS=sin(lonStart);

lonEnd=latLonEnd(2);
cosLE=cos(lonEnd);
sinLE=sin(lonEnd);

%The Cartesian starting point (unit sphere).
F=[cosLS*cosLF;
   sinLS*cosLF;
   sinLF];

%The Cartesian ending point (unit sphere).
T=[cosLE*cosLT;
   sinLE*cosLT;
   sinLT];

%This is more accurate than using the cosine formula that is implicit in
%Equations 13 and 15 in [1].
D=angBetweenVecs(F,T);
if(D==0)
   %Special case: The starting and ending points are the same.
    azStart=NaN;%Undefined direction.
    dist=0;
    azEnd=NaN;%Undefined direction.
    latLonWayPoints=repmat(latLonStart,[1,N+2]);
    return
end

cosD=cos(D);
sinD=sin(D);

%From Equation 17 in [1].
cosC=(sinLT-sinLF*cosD)/(cosLF*sinD);

%Deal with finite precision issues.
cosC=min(1,max(-1,cosC));
C=acos(cosC);

%The sign adjustment for going west is not mentioned in [1].
%Alternatively, one could get sinC=cosLT*sin(DLo)/sinD; and then use the
%atan2 function.
if(goingWest)
    azStart=-C;
else
    azStart=C;
end
dist=D*r;

%%%%%%%%%%%%%
%%%%The Azimuth at the End
%%%%%%%%%%%%%

%To get the azimuth at the end, use Equation 17 in [1] switching the
%beginning and the end; this flips the sign of D, so the sign of sinD
%flips. Then, flip the direction of the result 180 degrees if we are going
%West.
cosCEnd=-(sinLF-sinLT*cosD)/(cosLT*sinD);
%Deal with finite precision issues.
cosCEnd=max(-1,min(1,cosCEnd));
azEnd=acos(cosCEnd);
if(goingWest)
    azEnd=-azEnd;
end

%%%%%%%%%%%%%
%%%%Waypoints
%%%%%%%%%%%%%
if(~isempty(N))
    %If waypoints are desired.
    latLonWayPoints=zeros(2,N+2);
    latLonWayPoints(:,1)=latLonStart;
    latLonWayPoints(:,N+2)=latLonEnd;
    if(waypointType==0)
        %Waypoints uniformly spaced in distance.
        step=D/(N+1);
        DFX=step:step:(step*N);

        cosDFX=cos(DFX);
        sinDFX=sin(DFX);

        %Equation 19 in [1].
        sinLX=sinLF.*cosDFX+cosLF.*sinDFX.*cosC;
        cosLX=sqrt(1-sinLX.^2);%Equivalent to cosLX=cos(asin(sinLX));

        %Equation 21 in [1].
        cosDLoFX=(cosDFX-sinLF.*sinLX)./(cosLF.*cosLX);

        %Deal with finite precision issues.
        cosDLoFX=max(-1,min(1,cosDLoFX));
        if(goingWest)
            lonPoints=wrapRange(latLonStart(2)-acos(cosDLoFX),-pi,pi);
        else
            lonPoints=wrapRange(latLonStart(2)+acos(cosDLoFX),-pi,pi);
        end
        latLonWayPoints(:,2:(N+1))=[asin(sinLX);lonPoints];
    else
        %Waypoints uniformly spaced in longitude.
        DLonFT=wrapRange(lonEnd-lonStart,-pi,pi);

        if(DLonFT==0)
            %Meridonal trajectory.
            latLonWayPoints=[latLonStart,latLonEnd];
            return
        end

        deltaSign=sign(DLonFT);

        step=DLonFT/(N+1);
        DLoFX=abs(step:step:(step*N));
        cosDLoFX=cos(DLoFX);
        sinDLoFX=sin(DLoFX);

        sinC=sin(C);

        %Equation 38 in [1].
        tanLX=(cosC*sinDLoFX+sinLF*sinC*cosDLoFX)/(cosLF*sinC);

        latLonWayPoints(:,2:(N+1))=[atan(tanLX);wrapRange(lonStart+deltaSign*DLoFX,-pi,pi)];
    end
else
    latLonWayPoints=[];
end
end

function [azStart,dist,azEnd,latLonWayPoints]=indirectGreatCircleProbCrouse(latLonStart,latLonEnd,r,N)
%%INDIRECTGREATCIRCLEPROBCROUSE Solve the indirect great circle problem
%               using the algorithm of [1].
%
%REFERENCES:
%[1] D. F. Crouse, "Singularity-free great-circle sailing," Naval Research
%    Laboratory, Washington, DC, Tech. Rep. NRL/5340/MR-2021/4, 26
%    Jul. 2021.
%
%August 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Convert the latitude-longitude values to unit vectors.
u1=spher2Cart([latLonStart(2);latLonStart(1)]);
u2=spher2Cart([latLonEnd(2);latLonEnd(1)]);

%Rotate the vectors both onto the equator.
[vec1Rot,vec2Rot,R]=rot2Vecs2CommonPlane(u1,u2);

%The geodesic goes along the equator in the rotated coordinate system in
%the tangent plane. We need to find the shortest direction around the
%equator. To do this, we shall get the azimuth angles.
az1=atan2(vec1Rot(2),vec1Rot(1));
az2=atan2(vec2Rot(2),vec2Rot(1));
DLo=wrapRange(az2-az1,-pi,pi);
goingWest=(DLo<0);

%The distance traveled between the points.
dist=abs(DLo)*r;

%Unit vectors in the traveling directions are thus either the positive or
%the negative of the East unit direction vector at the starting and ending
%points.
uENU1=getENUAxes([0;az1],false,1,0);
uENU2=getENUAxes([0;az2],false,1,0);

if(goingWest)
    uDir1=-uENU1(:,1);
    uDir2=-uENU2(:,1);
else
    uDir1=uENU1(:,1);
    uDir2=uENU2(:,1);
end

%Now, we shall undo the rotations to get the unit direction vectors in the
%original Cartesian coordinate system.
uDir1=R'*uDir1;
uDir2=R'*uDir2;

%Given the unit vectors in the original Cartesian coordinate system, we
%need to get the azimuth angle. The azimuth angle is in terms of degrees
%East of North, so we need to project into the East-North plane.
%Get the new set of ENU basis vectors in the  original coordinate system
uENU1=getENUAxes(latLonStart,false,1,0);
uENU2=getENUAxes(latLonEnd,false,1,0);

EProj=dot(uDir1,uENU1(:,1));
NProj=dot(uDir1,uENU1(:,2));
azStart=atan2(EProj,NProj);

EProj=dot(uDir2,uENU2(:,1));
NProj=dot(uDir2,uENU2(:,2));
azEnd=atan2(EProj,NProj);

if(~isempty(N))%If waypoints are desired.
    distPts=linspace(0,dist,N+2);
    latLonWayPoints=directGreatCircleProb(latLonStart,azStart,distPts,r);
    %This ensures that the longitudes match even if the starting and/or
    %stopping points are at the poles.
    latLonWayPoints(:,1)=latLonStart;
    latLonWayPoints(:,N+2)=latLonEnd;
else
    latLonWayPoints=[];
end
end

function [azStart,dist,azEnd,latLonWayPoints]=indirectGreatCircleProbSnyder(latLonStart,latLonEnd,r,N)
%%INDIRECTGREATCIRCLEPRBSNYDER This function calls greatCircleDistance with
%   algorithm 0 selected, which should be the most stable solution. Then,
%   it calls greatCircleAzimuth with algorithm 2 selected, which is the
%   very simple formula in Equation 5-4b in [1].
%
%REFERENCES:
%[1] J. P. Snyder, "Map projections-a working manual," U.S. Geological
%    Survey, Tech. Rep. 1395, 1987.
%
%December 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Latitude at the start.
phi1=latLonStart(1);
%Latitude at the destination.
phi2=latLonEnd(1);

lambda1=latLonStart(2);
lambda2=latLonEnd(2);
cosPhi1=cos(phi1);
sinPhi1=sin(phi1);
cosPhi2=cos(phi2);
sinPhi2=sin(phi2);

cosLambda1=cos(lambda1);
sinLambda1=sin(lambda1);
cosLambda2=cos(lambda2);
sinLambda2=sin(lambda2);
deltaLambda=lambda2-lambda1;
cosDeltaLambda=cos(deltaLambda);
sinDeltaLambda=sin(deltaLambda);

%The Cartesian starting point (unit sphere).
F=[cosLambda1*cosPhi1;
   sinLambda1*cosPhi1;
   sinPhi1];
%The Cartesian ending point (unit sphere).
T=[cosLambda2*cosPhi2;
   sinLambda2*cosPhi2;
   sinPhi2];
dist=r*angBetweenVecs(F,T);

azStart=atan2(cosPhi2*sinDeltaLambda,cosPhi1*sinPhi2-sinPhi1*cosPhi2*cosDeltaLambda);
azEnd=atan2(-cosPhi1*sinDeltaLambda,cosPhi2*sinPhi1-sinPhi2*cosPhi1*cosDeltaLambda);
azEnd=wrapRange(azEnd+pi,-pi,pi);

if(~isempty(N))%If waypoints are desired.
    distPts=linspace(0,dist,N+2);
    latLonWayPoints=directGreatCircleProb(latLonStart,azStart,distPts,r);
    %This ensures that the longitudes match even if the starting and/or
    %stopping points are at the poles.
    latLonWayPoints(:,1)=latLonStart;
    latLonWayPoints(:,N+2)=latLonEnd;
else
    latLonWayPoints=[];
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function hammersley_step_set ( step )

%*****************************************************************************80
%
%% HAMMERSLEY_STEP_SET sets the step of the leaped Hammersley subsequence.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer STEP, the step of the leaped Hammersley subsequence.
%
  global hammersley_STEP

  step = floor ( step );

  if ( ~halham_step_check ( step ) )
    error ( 'HAMMERSLEY_STEP_SET - Fatal error!' );
  end

  hammersley_STEP = step;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y=wrap(x,sm,lg)

% y=wrap(x,sm,lg)
%
% Wrap the values of x so that they are >=sm and <lg

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

error(nargchk(3,3,nargin));
error(chk_param(x,'x','numeric','real'));
error(chk_param(sm,'sm','scalar','real'));
error(chk_param(sm,'lg','scalar','real'));
if (lg<=sm)
  error('lg must be > sm');
end

y=mod(x-sm,lg-sm)+sm;

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function vw = volume3View(vw, loc, xHairs, quickUpdate)
% vw = volume3View(vw, <loc>, <xHairs>, <quickUpdate=1>);
%
% Refresh a volume '3-view': a 'BrainVoyager'-like
% 3-window view, showing sagittal, coronal, and axial slices, at the
% specified location. 'loc' should be a 3-vector containing, respectiveley,
% the Superior/Inferior, Anterior/Posterior, and Right/Left coordinates.
% (For most of the volume anatomies used, S/I range from 1-256, A/P and
% R/L ranges from 1-124.
%
% If omitted, 'loc' is taken from the vw.loc field, or
% failing this set to the middle of the volume.
%
% xHairs: if set to nonzero, will place crosshairs on the current location,
% otherwise will omit them (0 default).
%
% 3/06/03 by ras
% 12/03 ras: using imagesc instead of imshow for now to display
% views, since there's a trickyness in accurately detecting the
% clicked-on location of images put up using imshow.
% 01/04 ras: cleaned up a lot, added 'zoom' field to ui struct allowing
% zooming in on vAnat. (see also zoom3view)
% to do: re-merge ras_recomputeImage with the regular recomputeImage.
% 11/04 ras: added flipping axes option
% 01/07 ras, bw: various ROI display options; trying to make it sensible

%TODO: Remove the input argument 'quickUpdate' - is not used

if ~exist('vw', 'var') || isempty(vw)
    vw = getSelectedVolume;
end

if ~exist('loc','var') || length(loc)<3
    loc = getLocFromUI(vw);
else
    setLocInUI(vw, loc);
end
vw.loc = loc;

% get the ui struct from the view
ui = viewGet(vw,'ui');

if ~exist('xHairs','var')
    if isfield(ui, 'crosshairs') && ui.crosshairs==1
        xHairs = 1;
    else
        xHairs = 0;
    end
end

%%%%% check if window exists
winTag = ['3VolumeWindow: ',vw.name];
winExists = findobj('Tag',winTag);

% if window doesn't exist, open the window
if isempty(winExists)
    openRaw3ViewWindow;
else
    % if one or more 3-view figures exist, make the most recent one current
    figure(winExists(end));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ***** Refresh Window ***** %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
modeStr=['ui.',ui.displayMode,'Mode'];
mode = eval(modeStr);
cmap = mode.cmap;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get images, apply phase/amp/param maps        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[vw, axi] = recompute3ViewImage(vw, loc(1), 1);
[vw, cor] = recompute3ViewImage(vw, loc(2), 2);
[vw, sag] = recompute3ViewImage(vw, loc(3), 3);

% not needed for anything, but memory is cheap and may be useful:
vw.ui.image = {axi cor sag};

% flip images if selected
if isfield(ui,'flipAP') && ui.flipAP==1
    sag = fliplr(sag);
    % update loc
end
if isfield(ui,'flipLR') && ui.flipLR==1
    axi = fliplr(axi);
    cor = fliplr(cor);
    % update loc
end

%%%%%%%%%%%%%%%%%%%%%
% put up each image %
%%%%%%%%%%%%%%%%%%%%%
if ~isempty(ui.underlayHandle) && ishandle(ui.underlayHandle(1))
    % handles exist, update image CData only
    set(ui.underlayHandle(1), 'CData', axi);
    set(ui.underlayHandle(2), 'CData', cor);
    set(ui.underlayHandle(3), 'CData', sag);
    quickUpdate = 1; % will bypass axes-related functions below
else
    % create image objects
    axes(ui.axiAxesHandle); cla; %#ok<*MAXES>
    vw.ui.underlayHandle(1) = image(axi);
    colormap(cmap); axis image, axis off

    axes(ui.corAxesHandle); cla;
    vw.ui.underlayHandle(2) = image(cor);
    colormap(cmap); axis image, axis off

    axes(ui.sagAxesHandle); cla;
    vw.ui.underlayHandle(3) = image(sag);
    colormap(cmap); axis image, axis off

    quickUpdate = 0; % will do zoom, text creation below
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If creating axes, add text annotation, zoom   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if quickUpdate==0

    % zoom in on selected coords %
    if ~isfield(ui,'zoom')
        % initialize to full range if not set
        ui.zoom = [1 size(vw.anat,1);...
            1 size(vw.anat,2);...
            1 size(vw.anat,3)];
    else
        axis(ui.axiAxesHandle, [ui.zoom(3,:) ui.zoom(2,:)]);
        axis(ui.corAxesHandle, [ui.zoom(3,:) ui.zoom(1,:)]);
        axis(ui.sagAxesHandle, [ui.zoom(2,:) ui.zoom(1,:)]);
    end

    % set a callback for each image which dynamically updates
    % the location based on where the user clicks
    % callback:
    % vw = recenter3view(vw,[orientationFlag]);
    for i=1:3
        bdf = sprintf('%s=recenter3view(%s,%i); ',vw.name, vw.name, i);
        set(vw.ui.underlayHandle(i), 'ButtonDownFcn', bdf);
    end

    % Add titles to indicate directions (esp. L-R: radiological conventions)
    % TO DO: do this at window creation, making them text children of
    % invisible axes (can't use text uicontrols, since they can't use
    % the TeX markup for the arrows)
    txtcol = 'k';
    SIsz = ui.zoom(1,2) - ui.zoom(1,1);
    APsz = ui.zoom(2,2) - ui.zoom(2,1);
    LRsz = ui.zoom(3,2) - ui.zoom(3,1);
    minS = ui.zoom(1,1);
    minA = ui.zoom(2,1);
    minL = ui.zoom(3,1);

    axes(ui.sagAxesHandle);
    if isfield(ui,'flipAP') && ui.flipAP==1
        title('Pos \leftrightarrow Ant','FontSize',10,'Color',txtcol);
    else
        title('Ant \leftrightarrow Pos','FontSize',10,'Color',txtcol);
    end

    axes(ui.corAxesHandle);
    if isfield(ui,'flipLR') && ui.flipLR==1
        title('Right \leftrightarrow Left','FontSize',10,'Color',txtcol);
    else
        title('Left \leftrightarrow Right','FontSize',10,'Color',txtcol);
    end

    % the I/S and P/A directions for the cor, axial don't
    % get flipped...
    text(minL-0.2*LRsz,minS + 0.5*SIsz, 'Inf \leftrightarrow Sup', ...
        'FontSize', 10,  'HorizontalAlignment', 'center', ...
        'Rotation', 90, 'Color', txtcol);
    axes(ui.axiAxesHandle);
    text(minL-0.2*LRsz,minA + 0.5*APsz, 'Pos \leftrightarrow Ant', ...
        'FontSize', 10,  'HorizontalAlignment', 'center', ...
        'Rotation', 90, 'Color', txtcol);

    % set tags so we can find each of these axes later
    set(ui.sagAxesHandle, 'Tag', 'sagAxes');
    set(ui.corAxesHandle, 'Tag', 'corAxes');
    set(ui.axiAxesHandle, 'Tag', 'axiAxes');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              draw ROIs if selected            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% remove any old ROI objects
hOld = findobj('Tag', sprintf('ROI_Handles_%s', viewGet(vw,'Name')));
delete(hOld);

% draw selected ROIs
if isfield(vw,'ui') & logical(ui.showROIs ~= 0) ...
        & ~isempty(vw.ROIs) & ~viewGet(vw, 'hideVolumeROIs')
	roiList = viewGet(vw, 'ROIsToDisplay');

    drawROIs3View(vw, roiList, ui.axiAxesHandle, 1, loc);
    drawROIs3View(vw, roiList, ui.corAxesHandle, 2, loc);
    drawROIs3View(vw, roiList, ui.sagAxesHandle, 3, loc);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add crosshairs to each image to specify current location %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if xHairs
    vw = renderCrosshairs(vw, 1);
end

% Draw colorbar
f = [ui.displayMode,'Mode'];  % field name of display mode
cmap = ui.(f).cmap;
numGrays = ui.(f).numGrays;
numColors = ui.(f).numColors;
setColorBar(vw, [], [numGrays+1 numGrays+numColors]);

colormap(cmap);

% update the annotation text
curScan     = viewGet(vw, 'curscan');
annotation  = viewGet(vw, 'annotation', curScan);
set(ui.annotationHandle, 'String', annotation);

% lastly, deal with the god-awful global variable issue
if strncmp(vw.name, 'VOLUME', 6)
	updateGlobal(vw);
end

return
% /----------------------------------------------------------------/ %

% /----------------------------------------------------------------/ %
function loc = getLocFromUI(vw)
% reads off the axi,cor, and sag UI edit fields to
% get the current view location.
loc(1) = str2num(get(vw.ui.sliceNumFields(1),'String')); %#ok<*ST2NM>
loc(2) = str2num(get(vw.ui.sliceNumFields(2),'String'));
loc(3) = str2num(get(vw.ui.sliceNumFields(3),'String'));
return
% /----------------------------------------------------------------/ %

% /----------------------------------------------------------------/ %
function vw = setLocInUI(vw, loc)
% sets the current view location in the axi, cor, sag UI edit fields.
set(vw.ui.sliceNumFields(1), 'String', num2str(loc(1)));
set(vw.ui.sliceNumFields(2), 'String', num2str(loc(2)));
set(vw.ui.sliceNumFields(3), 'String', num2str(loc(3)));
return
% /----------------------------------------------------------------/ %

% /----------------------------------------------------------------/ %
function drawROIs3View(vw, roiList, axs, ori, loc)
% draw the ROIs, parsing the relevant settings, on the axes
% specified by [axs]. This is implemented a little differently
% from drawROIs or drawROIsPerim, in that it doesn't select a given
% orientation each time (which caused some extra time during refreshes).
% roiList is an array into the vw.ROIs field.
ui = vw.ui;

% build prefs for the outline function.
prefs.method = ui.roiDrawMethod;
prefs.axesHandle = axs;

for r = roiList
    R = vw.ROIs(r);
    if ~isempty(R.coords)
        prefs.color = R.color;
        if r==vw.selectedROI, prefs.color=viewGet(vw,'selRoiColor'); end

        switch ori
            case 1, pts = R.coords([2 3], R.coords(1,:)==loc(1));
            case 2, pts = R.coords([1 3], R.coords(2,:)==loc(2));
            case 3, pts = R.coords([1 2], R.coords(3,:)==loc(3));
        end
		if isempty(pts), continue; end

        if isfield(ui, 'flipLR') && ui.flipLR==1 && ori<3
            % L/R flip affects columns of axi + coronal, but not sag, orientations
            dims = viewGet(vw,'Size');
            pts(2,:) = dims(3) - pts(2,:);
        end

        % restrict to voxels within zoom range
		zoom = viewGet(vw, 'zoom');
		switch ori
			case 1, zoom = zoom([2 3],:);
			case 2, zoom = zoom([1 3],:);
			case 3, zoom = zoom([1 2],:);
		end
		ok =  pts(1,:) >= zoom(1,1) & pts(1,:) <= zoom(1,2) & ...
				   pts(2,:) >= zoom(2,1) & pts(2,:) <= zoom(2,2) ;
		pts = pts(:,ok);

        h = outline(pts, prefs);
        if ishandle(h)
            set(h, 'ButtonDownFcn', sprintf('recenter3view(%s,%i);',vw.name, ori));
			set(h, 'Tag', sprintf('ROI_Handles_%s', vw.name));
        end

    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [uttName_sorted, aliID_sorted, vocab] = readKaldiAlignmentText(file_name)
lines = my_cat(file_name);

vocab = [];
for i=1:length(lines)
    PrintProgress(i,length(lines),1000);
    idx = regexp(lines{i}, ' ');
    uttName{i} = lines{i}(1:idx(1)-1);
    aliID{i} = str2num(lines{i}(idx(1)+1:end))';
    vocab = [vocab; aliID{i}];
    if mod(i,100)==0
        vocab = unique(vocab);
    end
end
vocab = unique(vocab);

[uttName_sorted, order] = sort(uttName);
aliID_sorted = aliID(order);

author: singaxiong

repo: SignalGraph


text:
%
% lellipf(phi, k, errtol)
%
% Inputs:
%
%   phi     Input angle vector size 1 or 1xN.
%   k       Input parameter vector size 1 or 1xN.
%   errtol  Error tolerance for Carlson's algorithms.
%
% Matlab function to compute Legendre's (incomplete) elliptic integral
% F(phi, k).  Uses a vectorized implementation of Carlson's Duplication Algorithms
% for symmetric elliptic integrals as found in "Computing Elliptic
% Integrals by Duplication," by B. C. Carlson, Numer. Math. 33, 1-16 (1979)
% and also found in ACM TOMS Algorithm 577.  Section 4 in the paper cited
% here describes how to convert between the symmetric elliptic integrals
% and Legendre's elliptic integrals.
%
% Returns NaN's for any argument values outside input range.
%

function f = lellipf(phi, k, errtol)

% Argument checking for vectorization:
lphi = length(phi);
lk = length(k);
errflag = logical(0);
if (lphi ~= lk)
    if (lphi==1)
        phivec = phi * ones(1,lk);
        kvec = k;
    elseif (lk==1)
        kvec = k * ones(1,lphi);
        phivec = phi;
    else
        disp('Incompatible input vector dimensions in lellipf!');
        errflag = logical(1);
    end
else
    phivec = phi;
    kvec = k;
end

if ~errflag
    snphi = sin(phivec);
    csphi = cos(phivec);
    csphi2 = csphi .* csphi;
    onesvec = ones(1,length(phivec));
    y = onesvec - kvec.*kvec .* snphi.*snphi;
    f = snphi .* rf(csphi2,  y, onesvec, errtol);
else
    f = NaN;
end