text:
function [beta,Mu_c,Sigma_y_tmp] = GMC(Priors, Mu, Sigma, x, in, out)
%GMC Gaussian Mixture Conditional. Returns the conditional P(out|in=x) of
% the Gaussian Mixture Model P(out,in) at the point in=x.
%
%
%
% Inputs -----------------------------------------------------------------
%   o Priors:  1 x K array representing the prior probabilities of the K GMM
%              components.
%   o Mu:      D x K array representing the centers of the K GMM components.
%   o Sigma:   D x D x K array representing the covariance matrices of the
%              K GMM components.
%   o x:       P x N array representing N datapoints of P dimensions.
%   o in:      1 x P array representing the dimensions to consider as
%              inputs.
%   o out:     1 x Q array representing the dimensions to consider as
%              outputs (D=P+Q)
%
% Output------------------------------------------------------------------------
%
%   o beta      :   N x K
%
%   o Mu_c      :   Q x N x K
%
%   o Sigma_c   :   Q x Q x K
%
%
%

nbData = size(x,2);
nbVar = size(Mu,1);
nbStates = size(Sigma,3);

%% Compute the influence of each GMM component, given input x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:nbStates
%     if i == 1
%         Mu(in,i)
%         Sigma(in,in,i)
%         inv(Sigma(in,in,i))
%         det(Sigma(in,in,i))
%         gaussPDF(x, Mu(in,i), Sigma(in,in,i))
%     end
    Pxi(:,i) = Priors(i).*gaussPDF(x, Mu(in,i), Sigma(in,in,i));
end
beta = Pxi./repmat(sum(Pxi,2)+realmin,1,nbStates);

%% Compute expected means y, given input x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:nbStates
    %   (N x D x K)
    %                   (D x N)                        (D x D)                               (D x N)
    Mu_c(:,:,j) = (repmat(Mu(out,j),1,nbData) + Sigma(out,in,j)*inv(Sigma(in,in,j)) * (x-repmat(Mu(in,j),1,nbData)))';
   %  Mu_c(:,:,j) = (repmat(Mu(out,j),1,nbData))';

end

%% Compute Marginal covariance matrices Sigma_y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:nbStates
    Sigma_y_tmp(:,:,j) = Sigma(out,out,j) - (Sigma(out,in,j)*inv(Sigma(in,in,j))*Sigma(in,out,j));
    Sigma_y_tmp(:,:,j) = 0.5 * (Sigma_y_tmp(:,:,j) + Sigma_y_tmp(:,:,j)');
end

end

author: epfl-lasa


text:
function [acc,estimated]=sim_NN(W1,W2,b1,b2,X,t)
[L,Col]=size(X);
error=0;
estimated=zeros(1,Col);
for k=1:Col
    estimated(k)=W2*logsig(W1*X(:,k)+b1)+b2;
    estimated_aprox=sign(estimated(k));
    if not(estimated_aprox==t(k))
        error=error+1;
    end
end
acc=(Col-error)/Col;

author: Sable

repo: mcbench-benchmarks


text:
function y = hadamard(x)
% y = hadamard(x)
%   x is a column vector or matrix
%   where the number of rows, m, must
%   be a power of 2
%
%   The transformation is symmetric and orthogonal.
%   To invert, apply again and then divide by m
%
%   This is the Help file for the mex file.
%   Mex file written by Peter Stobbe, stobbe@acm.caltech.edu, Aug 2008
%
%   This is NOT the builtin Matlab code "hadamard"
%
%   Note: in R2008b, Matlab added "fwht" and "ifwht" (the Fast Walsh-
%       Hadamart Transform and the inverse) to its Signal Processing
%       Toolbox.  With the default ordering and scaling, it's not
%       equivalent to this, but you can change this with the following:
%       y = length(x) * fwht( x, [], 'hadamard' );
%       Then y should be the same as hadamard(x) up to roundoff.
%       However, it appears that this code is faster than fwht.

error('You shouldn''t see this code -- mex file is not installed!');

repo: REST


text:
function [h, icons, plots, legend_text] = addLegendItem(item, pos)
% Add an item to the legend
% SYNTAX:
%   addLegendItem(<item>, <pos>);
%
% EXAMPLE:
%   addLegendItem("item", 2);;
%
% INPUT:
%   item = new entry
%   pos  = position in the list of inputs
%
% DEFAULT VALUES:
%   fig_handle = gcf;
%

%--------------------------------------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Andrea Gatti
%  Contributors:     Andrea Gatti, ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%    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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------
    l = legend;
    outm = get(l,'String');
    if nargin == 1
        pos = numel(outm) + 1;
    end
    if isempty(outm)
        [h, icons, plots, legend_text] = legend(item, 'Location', 'eastoutside');
        n_entry = 1;
        icons = icons(n_entry + 2 : 2 : end);
    else
        outm{pos} = item;
        [h, icons] = legend(outm, 'Location', 'eastoutside');
        n_entry = numel(outm);
        icons = icons(n_entry + 2 : 2 : end);
    end
end

author: goGPS-Project


text:
% this function compute the derivatives of the static coefficients. Its
% implementation follows that of the HMM Toolkit 3.2

function delta_coef = comp_delta(static_coef, DELTAWINDOW)

[N_vec, N_cep] = size(static_coef);

first_vec = static_coef(1,:)';
last_vec = static_coef(N_vec,:)';

static_coef = static_coef';
for i = 1:DELTAWINDOW
    static_coef = [first_vec static_coef];    % append the first feature vector DELTAWINDOW times in the front
    static_coef = [static_coef last_vec];    % append the last feature vector DELTAWINDOW times in the back
end
static_coef = static_coef';

delta_coef = zeros(size(static_coef));
% compute the delta coefficients
for i= DELTAWINDOW+1 : N_vec+DELTAWINDOW
    for j = 1:DELTAWINDOW
        delta_coef(i,:) = delta_coef(i,:) + j*(static_coef(i+j,:) - static_coef(i-j,:));
    end
end
delta_coef = delta_coef(DELTAWINDOW+1 : N_vec+DELTAWINDOW,:);

delta_coef = delta_coef / sum((1:DELTAWINDOW).^2) / 2;

author: singaxiong

repo: SignalGraph


text:
function findedge3(node,edge,range,varargin)
%% FINDEDGE3 highlights edges in certain range in 3-D
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

hold on
% set up range
if (nargin<=2) || (isempty(range))
    range = (1:size(edge,1))';
end
if islogical(range)
    range = find(range);
end
if size(range,2)>size(range,1)
    range = range';
end
midEdge = (node(edge(range,1),:)+node(edge(range,2),:))/2;
if length(range)<size(edge,1)
    h = plot3([node(edge(range,1),1)'; node(edge(range,2),1)'],...
              [node(edge(range,1),2)'; node(edge(range,2),2)'],...
              [node(edge(range,1),3)'; node(edge(range,2),3)'],...
              'LineWidth',1.5,'Color','r');
else
    h = plot3(midEdge(:,1),midEdge(:,2),midEdge(:,3),'r.','MarkerSize', 18);
end
if nargin > 3
    if strcmp(varargin{1},'noindex') || strcmp(varargin{1},'index')
        startidx = 2;
    else
        startidx = 1;
    end
    set(h,varargin{startidx:end});
end
if (nargin <=3) || ~(strcmp(varargin{1},'noindex'))
    text(midEdge(:,1)+0.025,midEdge(:,2)-0.025,midEdge(:,3),int2str(range), ...
         'FontSize',12,'FontWeight','bold','Color','k');
end
%% TODO: help and example

repo: ifem


text:
%% computing the bounds from the mean and standard deviations from dynare
prob=0.9;
Params={
	'alp'  ,  'beta'     ,    0.356000, 0.02
	'bet'  ,  'beta'     ,    0.993000, 0.002
	'gam'  ,  'normal'   ,    0.008500, 0.003
	'mst'  ,  'normal'   ,    1.000200, 0.007
	'rho'  ,  'beta'     ,    0.129000, 0.223
	'psi'  ,  'beta'     ,    0.650000, 0.05
	'del'  ,  'beta'     ,    0.010000, 0.005
	'sig_a',  'inv_gamma',    0.035449, inf
	'sig_m',  'inv_gamma',    0.008862, inf
	};
npars=size(Params,1);
bounds=nan(npars,2);
for ii=1:npars
    bounds(ii,:)=distributions.find_bounds(Params{ii,2},Params{ii,3},Params{ii,4},prob);
end
disp(bounds)
%% recover means and standard deviations
MOMS=nan(npars,2);
for ii=1:npars
    [a,b,moments,fval]=distributions.(Params{ii,2})(bounds(ii,1),bounds(ii,2),prob);
    MOMS(ii,1)=moments.mean;
    MOMS(ii,2)=moments.sd;
end
%% display and compare
disp('================ Means and Standard Deviations ================ ')
disp(num2cell(MOMS))
disp('================ Discrepancies ================')
disp(abs(cell2mat(Params(:,3:end))-MOMS))

author: jmaih


text:
function D = driving_function_mono_wfs_ls(x0,nx0,xs,f,conf)
%DRIVING_FUNCTION_MONO_WFS_LS driving signal for a line source in WFS
%
%   Usage: D = driving_function_mono_wfs_ls(x0,nx0,xs,nxs,f,conf)
%
%   Input parameters:
%       x0          - position of the secondary sources / m [nx3]
%       xs          - position and orientation of virtual line source / m [nx3]
%                     or [nx6]
%       f           - frequency of the monochromatic source / Hz
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       D           - driving function signal [nx1]
%
%   DRIVING_FUNCTION_MONO_WFS_LS(x0,nx0,xs,f,src,conf) returns WFS driving
%   signals for the given secondary sources, the virtual line source position,
%   its orientation xs(:,4:6), which is parallel to the line source, and the
%   frequency f. If no explicit orientation is given, [0 0 1] is assumed.
%
%   See also: driving_function_mono_wfs, driving_function_imp_wfs_ps

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% 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.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 5;
nargmax = 5;
narginchk(nargmin,nargmax);
isargmatrix(x0,nx0,xs);
isargpositivescalar(f);
isargstruct(conf);

%% ===== Configuration ==================================================
xref = conf.xref;
c = conf.c;
dimension = conf.dimension;
driving_functions = conf.driving_functions;

%% ===== Computation ====================================================
% Calculate the driving function in time-frequency domain

omega = 2*pi*f;
[xs,nxs] = get_position_and_orientation_ls(xs,conf);

if strcmp('2D',dimension)

    % === 2-Dimensional ==================================================

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        % D using a line source
        %
        %              iw (x0-xs) nx0   (2)/ w         \
        % D(x0,w) =  - -- -----------  H1  | - |x0-xs| |
        %              2c   |x0-xs|        \ c         /
        %
        % https://sfs.rtfd.io/en/3.2/d_wfs/#equation-fd-wfs-line
        %
        % r = |x0-xs|
        r = vector_norm(x0-xs,2);
        % Driving signal
        D = -1i.*omega./(2.*c) ...
            .* vector_product(x0-xs,nx0,2) ./ r ...
            .* besselh(1,2,omega./c.*r);
        %
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a line source.'],upper(mfilename),driving_functions);
    end

elseif strcmp('2.5D',dimension)

    % === 2.5-Dimensional ================================================

    % Reference point
    xref = repmat(xref,[size(x0,1) 1]);

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        % 2.5D correction factor
        %        ______________
        % g0 = \| 2pi |xref-x0|
        %
        g0 = sqrt(2*pi*vector_norm(xref-x0,2));
        %
        % D_2.5D using a line source
        %                         ___
        %                   g0   |i w  (x0-xs) nx0   (2)/ w         \
        % D_2.5D(x0,w) =  - -- _ |---  -----------  H1  | - |x0-xs| |
        %                   2   \| c    |x0-xs|         \ c         /
        %
        % https://sfs.rtfd.io/en/3.2/d_wfs/#equation-fd-wfs-line-25d
        %
        % r = |x0-xs|
        r = vector_norm(x0-xs,2);
        % Driving signal
        D = -g0./2 .* sqrt(i.*omega./c) ...
            .* vector_product(x0-xs,nx0,2) ./ r ...
            .* besselh(1,2,omega./c.*r);
        %
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a 2.5D line source.'],upper(mfilename),driving_functions);
    end

elseif strcmp('3D',dimension)

    % === 3-Dimensional ==================================================

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        % D using a line source
        %
        %              iw   v nx0    (2)/ w     \
        % D(x0,w) =  - -- --------  H1  | - |v| | ,
        %              2c   |v|         \ c     /
        %
        % where v = x0-xs - <x0-xs,nxs > nxs,
        % and |nxs| = 1.
        %
        % https://sfs.rtfd.io/en/3.2/d_wfs/#equation-fd-wfs-line
        %
        % v = (I - nxs'nxs)(x0-xs)
        % r = |v|
        nxs = nxs(1,:);
        v = (x0 - xs)*(eye(3) - nxs'*nxs);
        r = vector_norm(v,2);
        % Driving signal
        D = -1i*omega/(2*c) .* vector_product(v,nx0,2) ./ r .* ...
            besselh(1,2,omega/c.*r);
        %
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a line source.'],upper(mfilename),driving_functions);
    end

else
    error('%s: the dimension %s is unknown.',upper(mfilename),dimension);
end

author: sfstoolbox

repo: sfs-matlab


text:
function [Hq,tq,hq,Dq,Fq] = MFDFA1(signal,scale,q,m,Fig)
% Multifractal detrended fluctuation analysis (MFDFA)
%
% [Hq,tq,hq,Dq,Fq]=MFDFA(signal,scale,q,m,Fig);
%
% INPUT PARAMETERS---------------------------------------------------------
%
% signal:       input signal
% scale:        vector of scales
% q:            q-order that weights the local variations
% m:            polynomial order for the detrending
% Fig:          1/0 flag for output plot of Fq, Hq, tq and multifractal
%               spectrum (i.e. Dq versus hq).
%
% OUTPUT VARIABLES---------------------------------------------------------
%
% Hq:           q-order Hurst exponent
% tq:           q-order mass exponent
% hq:           q-order singularity exponent
% Dq:           q-order dimension
% Fq:           q-order scaling function
%
% EXAMPLE------------------------------------------------------------------
%
% load fractaldata
% scmin=16;
% scmax=1024;
% scres=19;
% exponents=linspace(log2(scmin),log2(scmax),scres);
% scale=round(2.^exponents);
% q=linspace(-5,5,101);
% m=1;
% signal1=multifractal;
% signal2=monofractal;
% signal3=whitenoise;
% [Hq1,tq1,hq1,Dq1,Fq1]=MFDFA1(signal1,scale,q,m,1);
% [Hq2,tq2,hq2,Dq2,Fq2]=MFDFA1(signal2,scale,q,m,1);
% [Hq3,tq3,hq3,Dq3,Fq3]=MFDFA1(signal3,scale,q,m,1);
%--------------------------------------------------------------------------
warning off
X=cumsum(signal-mean(signal));
if min(size(X))~=1||min(size(scale))~=1||min(size(q))~=1;
    error('Input arguments signal, scale and q must be a vector');
end
if size(X,2)==1;
   X=transpose(X);
end
if min(scale)<m+1
   error('The minimum scale must be larger than trend order m+1')
end
for ns=1:length(scale),
    segments(ns)=floor(length(X)/scale(ns));
    for v=1:segments(ns),
        Index=((((v-1)*scale(ns))+1):(v*scale(ns)));
        C=polyfit(Index,X(Index),m);
        fit=polyval(C,Index);
        RMS_scale{ns}(v)=sqrt(mean((X(Index)-fit).^2));
    end
    for nq=1:length(q),
        qRMS{nq,ns}=RMS_scale{ns}.^q(nq);
        Fq(nq,ns)=mean(qRMS{nq,ns}).^(1/q(nq));
    end
    Fq(q==0,ns)=exp(0.5*mean(log(RMS_scale{ns}.^2)));
end
for nq=1:length(q),
    C = polyfit(log2(scale),log2(Fq(nq,:)),1);
    Hq(nq) = C(1);
    qRegLine{nq} = polyval(C,log2(scale));
end
tq = Hq.*q-1;
hq = diff(tq)./(q(2)-q(1));
Dq = (q(1:end-1).*hq)-tq(1:end-1);

%OUTPUT FIGURE-------------------------------------------------------------
if Fig==1,
   qindex=[1,round(length(q)/2),length(q)];
   qindex2=[1,round((length(q)-1)/2),length(q)-1];

   %Variable settings---------
   qstart1=q(qindex(1));
   qmid1=q(qindex(2));
   qstop1=q(qindex(3));
   Hqstart1=Hq(qindex(1));
   Hqmid1=Hq(qindex(2));
   Hqstop1=Hq(qindex(3));
   tqstart1=tq(qindex(1));
   tqmid1=tq(qindex(2));
   tqstop1=tq(qindex(3));
   hqstart2=hq(qindex2(1));
   hqmid2=hq(qindex2(2));
   hqstop2=hq(qindex2(3));
   Dqstart1=Dq(qindex2(1));
   Dqmid1=Dq(qindex2(2));
   Dqstop1=Dq(qindex2(3));
   for nq=1:length(qindex),
       qRegFit(nq,:)=qRegLine{qindex(nq)};
   end
   X1=log2(scale);
   YMatrix1=[log2(Fq(qindex,:));qRegFit];
   X2=q;
   Y1=Hq;
   Y3=tq;
   X4=hq;
   Y5=Dq;
   %---------------------------
    q_end=num2str(max(q));
    if length(find(q==0))==1
       q_middle=0;
    else
       qm0=min(q)+((max(q)-min(q))/2);
       qm1=[q(find(q<qm0, 1, 'last' )),q(find(q>qm0, 1 ))];
       qm2=qm1-qm0;
       midind= qm2==min(qm2);
       q_middle=qm2(midind(1));
    end
    q_mid=num2str(q_middle);
    q_start=num2str(min(q));

    % Create figure
    figure1 = figure('PaperSize',[20.98 29.68],'Color',[1 1 1]);

    scaleInd=floor(min(log2(scale))):ceil(max(log2(scale)));
    for ns=1:length(scaleInd),
        scaletick{ns}=num2str(2^scaleInd(ns));
    end
    Fqind=floor(log2(min(min(Fq)))):ceil(log2(max(max(Fq))));
    for nf=1:length(Fqind),
        Fqtick{nf}=num2str(Fqind(nf));
    end

    % Create subplot
    subplot1 = subplot(2,2,1,'Parent',figure1,...
        'YTickLabel',Fqtick,...
        'XTickLabel',scaletick,...
        'XTick',scaleInd,...
        'LineWidth',2,...
        'FontSize',14);
    hold(subplot1,'all');

    % Create multiple lines using matrix input to plot
    plot1 = plot(X1,YMatrix1,'Parent',subplot1);
    set(plot1(1),'MarkerFaceColor',[1 0 0],'MarkerEdgeColor',[0 0 0],...
        'Marker','o',...
        'LineStyle','none',...
        'Color',[1 0 0],...
        'DisplayName',strcat('q = ',num2str(min(q))));
    set(plot1(2),'MarkerFaceColor',[0 0 1],'MarkerEdgeColor',[0 0 0],...
        'Marker','o',...
        'LineStyle','none',...
        'Color',[0 0 1],...
        'DisplayName',strcat('q = ',num2str(q_middle)));
    set(plot1(3),'MarkerFaceColor',[0 0.498 0],'MarkerEdgeColor',[0 0 0],...
        'Marker','o',...
        'LineStyle','none',...
        'Color',[0 0.498 0],...
        'DisplayName',strcat('q = ',num2str(max(q))));
    set(plot1(4),'LineWidth',2,'Color',[1 0 0],'DisplayName',strcat('q = ',num2str(min(q))));
    set(plot1(5),'MarkerFaceColor',[1 0 0],'MarkerEdgeColor',[1 1 1],...
        'LineWidth',2,...
        'Color',[0 0 1],...
        'DisplayName',strcat('q = ',num2str(q_middle)));
    set(plot1(6),'LineWidth',2,'Color',[0 0.498 0],'DisplayName',strcat('q = ',num2str(max(q))));

    % Create xlabel
    xlabel('Scale (segment sample size)','FontSize',14);

    % Create ylabel
    ylabel('log_2(Fq)','FontSize',14);

    % Create title
    title('Scaling function Fq (q-order RMS)','FontSize',14);

    qInd=floor(min(q)):ceil(max(q));
    for nq=1:length(qInd),
        qtick{nq}=num2str(qInd(nq));
    end

    % Create subplot
    subplot2 = subplot(2,2,2,'Parent',figure1,...
        'XTickLabel',qtick,...
        'XTick',qInd,...
        'LineWidth',2,...
        'FontSize',14);
    hold(subplot2,'all');

    % Create plot
    plot(X2,Y1,'Parent',subplot2,'Marker','o','LineStyle','none',...
        'Color',[0 0 1]);

    % Create xlabel
    xlabel('q','FontSize',16);

    % Create ylabel
    ylabel('Hq','FontSize',16);

    % Create title
    title('q-order Hurst exponent','FontSize',14);

    % Create plot
    plot(qstart1,Hqstart1,'Parent',subplot2,'MarkerFaceColor',[1 0 0],...
        'MarkerEdgeColor',[0 0 0],...
        'MarkerSize',8,...
        'Marker','o',...
        'LineWidth',2,...
        'LineStyle','none',...
        'Color',[1 0 0],...
        'DisplayName',strcat('Hq(',num2str(min(q)),') = ',num2str(Hq(q==min(q)))));

    % Create plot
    plot(qmid1,Hqmid1,'Parent',subplot2,'MarkerFaceColor',[0 0 1],...
        'MarkerEdgeColor',[0 0 0],...
        'MarkerSize',8,...
        'Marker','o',...
        'LineWidth',2,...
        'LineStyle','none',...
        'Color',[0 0 1],...
        'DisplayName',strcat('Hq(',num2str(q_middle),') = ',num2str(Hq(q==q_middle))));

    % Create plot
    plot(qstop1,Hqstop1,'Parent',subplot2,'MarkerFaceColor',[0 0.498 0],...
        'MarkerEdgeColor',[0 0 0],...
        'MarkerSize',8,...
        'Marker','o',...
        'LineWidth',2,...
        'LineStyle','none',...
        'Color',[0 0.498 0],...
        'DisplayName',strcat('Hq(',num2str(max(q)),') = ',num2str(Hq(q==max(q)))));

    % Create subplot
    subplot3 = subplot(2,2,3,'Parent',figure1,...
        'XTickLabel',qtick,...
        'XTick',qInd,...
        'LineWidth',2,...
        'FontSize',14);
    hold(subplot3,'all');

    % Create plot
    plot(X2,Y3,'Parent',subplot3,'Marker','o','LineStyle','none',...
        'Color',[0 0 1]);

    % Create xlabel
    xlabel('q','FontSize',16);

    % Create ylabel
    ylabel('tq','FontSize',16);

    % Create title
    title('q-order Mass exponent','FontSize',14);

    % Create plot
    plot(qstop1,tqstop1,'Parent',subplot3,'MarkerFaceColor',[0 0.498 0],...
        'MarkerEdgeColor',[0 0 0],...
        'MarkerSize',8,...
        'Marker','o',...
        'LineWidth',2,...
        'LineStyle','none',...
        'Color',[0 0.498 0],...
        'DisplayName',strcat('tq(',num2str(max(q)),') = ',num2str(tq(q==max(q)))));

    % Create plot
    plot(qstart1,tqstart1,'Parent',subplot3,'MarkerFaceColor',[1 0 0],...
        'MarkerEdgeColor',[0 0 0],...
        'MarkerSize',8,...
        'Marker','o',...
        'LineWidth',2,...
        'LineStyle','none',...
        'Color',[1 0 0],...
        'DisplayName',strcat('Hq(',num2str(min(q)),') = ',num2str(Hq(q==min(q)))));

    % Create plot
    plot(qmid1,tqmid1,'Parent',subplot3,'MarkerFaceColor',[0 0 1],...
        'MarkerEdgeColor',[0 0 0],...
        'MarkerSize',8,...
        'Marker','o',...
        'LineWidth',2,...
        'LineStyle','none',...
        'Color',[0 0 1],...
        'DisplayName',strcat('Hq(',num2str(q_middle),') = ',num2str(Hq(q==q_middle))));

    % Create subplot
    subplot4 = subplot(2,2,4,'Parent',figure1,'LineWidth',2,'FontSize',12);
    % Uncomment the following line to preserve the X-limits of the axes
    % xlim(subplot4,[0.2 1.8]);
    hold(subplot4,'all');

    % Create plot
    plot(X4,Y5,'Parent',subplot4,'Marker','o','LineStyle','none',...
        'Color',[0 0 1]);

    % Create xlabel
    xlabel('hq','FontSize',16);

    % Create ylabel
    ylabel('Dq','FontSize',16);

    % Create title
    title('Multifractal spectrum','FontSize',14);

    % Create plot
    plot(hqstart2,Dqstart1,'Parent',subplot4,'MarkerFaceColor',[1 0 0],...
        'MarkerEdgeColor',[0 0 0],...
        'MarkerSize',8,...
        'Marker','o',...
        'LineWidth',2,...
        'LineStyle','none',...
        'Color',[1 0 0],...
        'DisplayName',[strcat('Dq(',num2str(min(q)),') = ',num2str(Dq(q==min(q)))),sprintf('\n'),strcat('hq(',num2str(min(q)),') = ',num2str(hq(q==min(q))))]);

    % Create plot
    plot(hqmid2,Dqmid1,'Parent',subplot4,'MarkerFaceColor',[0 0 1],...
        'MarkerEdgeColor',[0 0 0],...
        'MarkerSize',8,...
        'Marker','o',...
        'LineWidth',2,...
        'LineStyle','none',...
        'Color',[0 0 1],...
        'DisplayName',[strcat('Dq(',num2str(q_middle),') = ',num2str(Dq(q==q_middle))),sprintf('\n'),strcat('hq(',num2str(q_middle),') = ',num2str(hq(q==q_middle)))]);

    % Create plot
    plot(hqstop2,Dqstop1,'Parent',subplot4,'MarkerFaceColor',[0 0.498 0],...
        'MarkerEdgeColor',[0 0 0],...
        'MarkerSize',8,...
        'Marker','o',...
        'LineWidth',2,...
        'LineStyle','none',...
        'Color',[0 0.498 0],...
        'DisplayName',[strcat('Dq(',num2str(max(q)),') = ',num2str(Dq(find(q==max(q))-1))),sprintf('\n'),strcat('hq(',num2str(max(q)),') = ',num2str(hq(find(q==max(q))-1)))]);

    % Create legend
    %legend1 = legend(subplot1,'show');
    %set(legend1,'Position',[0.4174 0.6072 0.09288 0.2055]);
    % Create legend
    legend1 = legend(subplot1,'show');
    set(legend1,'Position',[0.4126 0.6072 0.1024 0.2055]);

    % Create legend
    legend2 = legend(subplot3,'show');
    set(legend2,'Position',[0.3414 0.1576 0.1458 0.1714]);

    % Create legend
    legend3 = legend(subplot2,'show');
    set(legend3,'Position',[0.8171 0.7456 0.1528 0.1714]);

    % Create legend
    legend4 = legend(subplot4,'show');
    set(legend4,'Position',[0.8338 0.2436 0.1354 0.2764]);

    % Create textbox
    annotation(figure1,'textbox',[0.2726 0.638 0.1111 0.05249],...
        'String',{'Slope = Hq'},...
        'HorizontalAlignment','center',...
        'FontSize',14,...
        'FitBoxToText','off');

    % Create textbox
    annotation(figure1,'textbox',[0.6355 0.1669 0.2013 0.05518],...
        'String',{strcat('hq_m_a_x - hq_m_i_n = ',num2str(max(hq)-min(hq)))},...
        'HorizontalAlignment','center',...
        'FontSize',14,...
        'FitBoxToText','off',...
        'LineStyle','none');
end


text:
classdef MimOrientationSagittal < MimGuiPlugin
    % MimOrientationSagittal. Gui Plugin for using a preset bone window/level
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimOrientationSagittal is a Gui Plugin for the MIM Toolkit.
    %     The gui will create a button for the user to run this plugin.
    %     Running this plugin will change the window and level of the viewing
    %     panel to standard soft tissue values.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Sagittal'
        SelectedText = 'Sagittal'
        ToolTip = 'Changes the image orientation to a sagittal view'
        Category = 'View'
        Visibility = 'Dataset'
        Mode = 'Toolbar'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
        Icon = 'sagittal_thumb.png'
        Location = 4
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.Orientation = GemImageOrientation.YZ;
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded;
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.Orientation == GemImageOrientation.YZ && ~strcmp(gui_app.ImagePanel.Mode, char(MimModes.View3DMode));
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%% === compare timings for MULT

% number of repeats (to increase measurement accuracy)
R  = 10;

% what sizes of (square) matrices do we want measure? (log spacing)
nn = unique(round(exp(linspace(log(4),log(100),10))));

% number of different sizes to try
N  = length(nn);

% size of the output C, in Megabytes
Mbytes = 10;

time   = zeros(N,3,R);
ops    = zeros(N,1);

candidates = {'A*B',  'mmx_naive',  @(a,b)mmx_naive('m',a,b);
              'A''*B',  'mmx_naive',  @(a,b)mmx_naive('m',a,b,'tn');
              'A*B''',  'mmx_naive',  @(a,b)mmx_naive('m',a,b,'nt');
              'A''*B''',  'mmx_naive',  @(a,b)mmx_naive('m',a,b,'tt');
              'A*Bs',  'mmx_naive',  @(a,b)mmx_naive('m',a,b,'ns');
              'As*B',  'mmx_naive',  @(a,b)mmx_naive('m',a,b,'sn');
};

funcs = cell(0,0);
for i=1:size(candidates,1)
   found_it = ~isempty(which(candidates{i,2}));
   if found_it
      funcs(end+1,:) = candidates(i,[1 3]); %#ok<SAGROW>
   else
      fprintf('Could not find %s in your path, ignoring.\n',candidates{i,2});
   end
end
nf    = length(funcs);

for i = 1:N

   n = nn(i);

   r1 = n;
   c1 = n;
   r2 = c1;
   c2 = n;

   pages = round(1e6*Mbytes / (8*r1*c2));

   A     = rand(r1,c1,pages);
   B     = rand(r2,c2,pages);

   B(B>0.5) = 0;
   A(A>0.5) = 0;

   ops(i)   = r1*c2*(2*c1-1)*pages;
   fprintf('multiplying %d*%d matrix pairs of dimension [%dx%d]\n',R*nf,pages,n,n)

   for r = 1:R
      for f = 1:nf %nf:-1:1
         tstart = tic;
         C = funcs{length(funcs)+1-f,2}(A,B);
         time(i,length(funcs)+1-f,r)  = toc(tstart);
         pause(1/10000);
      end
   end
end

%% graphics

gflops   = 1e-9*bsxfun(@times, ops, 1./time);
gFLP     = mean(gflops,3);

clf
hold on
cols = get(gca,'ColorOrder');
cols = [cols; cols];
for i = 1:size(gflops,2)
   plot(nn,gFLP(:,i),'color',cols(i,:),'linewidth',3);
end
for i = 1:size(gflops,2)
   plot(nn,squeeze(gflops(:,i,:)),'.','color',cols(i,:))
end

grid on
set(gca,'xlim',[1 max(nn)])
ylabel('\bf Gigaflops')
xlabel('\bf dimension')

legend(funcs{:,1},'location','northwest')

title('\bf Comparison between mmx, ndfun, and mtimesx. (in Gflops, bigger is better)')

%% make PDF
% set(gcf,'color','w')
% export_fig 'comparison' -png

author: NicolasBoumal

repo: manopt


text:
function triangle_num_test ( )

%*****************************************************************************80
%
%% TRIANGLE_NUM_TEST tests TRIANGLE_NUM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGLE_NUM_TEST\n' );
  fprintf ( 1, '  TRIANGLE_NUM computes the triangular numbers.\n' );
  fprintf ( 1, '\n' );

  for n = 1 : 10
    fprintf ( 1, '  %2d  %6d\n', n, triangle_num ( n ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [hdr,be] = read_hdr_raw(fname)
% Read a NIFTI-1 hdr file
% FORMAT [hdr,be] = read_hdr_raw(fname)
% fname - filename of image
% hdr   - a structure containing hdr info
% be    - whether big-endian or not
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: read_hdr_raw.m 1143 2008-02-07 19:33:33Z spm $

hdr = [];
be  = [];
ok  = true;

% Name of header file
[pth,nam,ext] = fileparts(fname);
switch ext
case {'.img','.hdr'}
    hname = fullfile(pth,[nam '.hdr']);
case {'.nii'}
    hname = fullfile(pth,[nam '.nii']);
otherwise
    hname = fullfile(pth,[nam '.hdr']);
end;

% Open it if possible
fp  = fopen(hname,'r','native');
if fp==-1
    hdr = [];
    return;
end;

% Sort out endian issues of header
[unused,unused,mach] = fopen(fp);
if strcmp(mach,'ieee-be') || strcmp(mach,'ieee-be.l64')
    be = true;
else
    be = false;
end;
fseek(fp,0,'bof');
fseek(fp,40,'bof');
nd = fread(fp,1,'int16')';
if isempty(nd),
    fclose(fp);
    hdr = [];
    return;
elseif nd<1 || nd>7
    be = ~be;
    fclose(fp);
    if be, mach = 'ieee-be';
    else   mach = 'ieee-le';
    end;
    fp = fopen(hname,'r',mach);
    if fp==-1
        hdr = [];
        return;
    end;
end;

% Is it NIFTI or not
fseek(fp,0,'bof');
fseek(fp,344,'bof');
mgc = deblank(char(fread(fp,4,'uint8')'));
switch mgc
case {'ni1','n+1'}
    org = niftistruc;
otherwise
    org = mayostruc;
end;
fseek(fp,0,'bof');
% Read the fields
for i=1:length(org)
    tmp = fread(fp,org(i).len,['*' org(i).dtype.prec])';
    if length(tmp) ~= org(i).len
disp([length(tmp) org(i).len]);
        tmp = org(i).def;
        ok  = false;
    end;
    tmp = feval(org(i).dtype.conv,tmp);
    hdr.(org(i).label) = tmp;
end;

fclose(fp);
if ~ok,
     fprintf('There was a problem reading the header of\n');
     fprintf('"%s".\n', fname);
     fprintf('It may be corrupted in some way.');
end;
return;

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function h = assignColumns(f, colIdx, g)
%ASSIGNCOLUMNS   Extract columns (or rows) of an array-valued TRIGTECH.
%   G = ASSIGNCOLUMNS(F, COLIDX, G) assigns the columns specified by the row
%   vector COLIDX from the FUN F so that F(:, COLIDX) = G. COLIDX need not be
%   increasing in order or unique, but must contain only integers in the range
%   [1, M] (where F has M columns) and satisfy LENGTH(COLIDX) = SIZE(G, 2) or
%   ISEMPTY(G).
%
% See also EXTRACTCOLUMNS, MAT2CELL.

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

% g is empty - remove columns:
if ( isempty(g) )
    h = f;
    h.values(:,colIdx) = [];
    h.coeffs(:,colIdx) = [];
    h.isReal(:,colIdx) = [];
    return
end

% Prolong so that f and g have the same length:
if ( length(f) > length(g) )
    g = prolong(g, length(f));
elseif ( length(g) > length(f) )
    f = prolong(f, length(g));
end

% Assign the columns of h.values and h.coeffs:
h = f;
h.values(:, colIdx) = g.values;
h.coeffs(:, colIdx) = g.coeffs;

% Update ishappy:
h.ishappy = f.ishappy && g.ishappy;
h.isReal(colIdx) = g.isReal;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function p = exactness_1d ( family, growth, level )

%*****************************************************************************80
%
%% EXACTNESS_1D returns the exactness of a 1D quadrature rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 March 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FAMILY, indicates the family of 1D quadrature rules used
%    to provide factors.
%    'CC': Clenshaw-Curtis;
%    'L':  Gauss-Legendre
%
%    Input, string GROWTH, indicates the growth rate.
%    'E': exponential, for CC, L;
%    'HC': Hyperbolic Cross;
%    'L': linear growth, for CC, HC, L;
%    'LO', linear odd growth, for CC, L;
%    'SE': slow exponential, for CC, L.
%
%    Input, int LEVEL, the level of the sparse grid to be constructed.
%    0 <= LEVEL <= 5.
%
%    Output, integer P, the exactness of the rule.
%
  if ( upper ( family ) == 'CC' )

    o = order_1d ( family, growth, level );
    p = o;

  elseif ( upper ( family ) == 'L' )

    o = order_1d ( family, growth, level );
    p = 2 * o - 1;

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'EXACTNESS_1D - Fatal error!\n' );
    fprintf ( 1, '  Unknown family.\n' );
    error ( 'EXACTNESS_1D - Fatal error!' );

  end

  return
end


text:
function s = DotXBLAS(x,y)
%DOTXBLAS     Dot product 'as if' computed in 2-fold precision
%
%   res = DotXBLAS(x,y)
%
%On return, res approximates x'*y with accuracy as if computed
%  in 2-fold precision.
%
%Implements algorithm BLAS_ddot_x from
%  X. Li, J. Demmel, D. Bailey, G. Henry, Y. Hida, J. Iskandar,
%    W. Kahan, S. Kang, {S.}, A. Kapur, M. Martin, B. Thompson, {B.},
%    T. Tung, {T.}, D. Yoo: Design, Implementation and Testing of
%    Extended and Mixed Precision BLAS, ACM Trans. Math. Software,
%    2(28), p. 152-205, 2002.
%Requires 37n flops.
%
%Reference implementation! Slow due to interpretation!
%

% written  03/03/07     S.M. Rump
% modified 05/09/09     S.M. Rump  rounding to nearest, complex input
%

  if ( ~isreal(x) ) | ( ~isreal(y) )
    error('DotXBLAS for real input only')
  end

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  s = 0;
  t = 0;
  for i=1:length(x)
    [h,r] = TwoProduct(x(i),y(i));
    [s1,s2] = TwoSum(s,h);
    [t1,t2] = TwoSum(t,r);
    s2 = s2 + t1;
    [t1,s2] = FastTwoSum(s1,s2);
    t2 = t2 + s2;
    [s,t] = FastTwoSum(t1,t2);
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function dUpsilonS = lfmvvGradientSigmaUpsilonMatrix(gamma, sigma2, ...
    t1, t2, mode)

% LFMVVGRADIENTSIGMAUPSILONMATRIX Gradient of upsilon matrix vv wrt sigma
% FORMAT
% DESC computes the gradient wrt sigma of a portion of the LFMVV kernel.
% ARG gamma : Gamma value for system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG mode : operation mode, according to the derivative (mode 0,
% derivative wrt t1, mode 1 derivative wrt t2)
% RETURN upsilon : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio Alvarez, 2010
%
% SEEALSO : lfmvvComputeUpsilonMatrix.m

% KERN

gridt1 = repmat(t1, 1, length(t2));
gridt2 = repmat(t2', length(t1), 1);
timeGrid = gridt1 - gridt2;

dUpsilon = lfmvpGradientSigmaUpsilonMatrix(gamma, sigma2, t1, t2, mode);

if mode == 0
    dUpsilonS = gamma*dUpsilon - 4*timeGrid/(sqrt(pi)*sigma2^2).* ...
        exp(-(timeGrid.^2)./sigma2).*(3 - 2*(timeGrid.^2)/sigma2) ...
        + 2*gamma/(sqrt(pi)*sigma2)*exp(-gamma*t1)*((1-2*t2.^2/sigma2).* ...
        exp(-(t2.^2)/sigma2)).';
else
    dUpsilonS = -gamma*dUpsilon - 4*timeGrid/(sqrt(pi)*sigma2^2).* ...
        exp(-(timeGrid.^2)./sigma2).*(3 - 2*(timeGrid.^2)/sigma2);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% Test file for @deltafun/restrict.m

function pass = test_restrict(pref)

if (nargin < 1)
    pref = chebfunpref();
end
%%
d = deltafun();
pass(1) = isempty(restrict(d, [-.5, .5]));

f = bndfun(@sin);
d = deltafun(f, struct('deltaMag', 1, 'deltaLoc', 0));
pass(2) = ~isa(restrict(d, [-1, -.5]), 'deltafun');
pass(3) = ~isa(restrict(d, [.5, 1]), 'deltafun');
pass(4) = anyDelta(restrict(d, [-.5, .5]));

f = fun.constructor(@(x) exp(x), struct('domain', [-1, 1]));
d = deltafun(f, struct('deltaMag', [1 1 1 1], 'deltaLoc', [-.5, -.25, 0, 1]));
A = restrict(d, [-1 0 .5 1]);
d1 = A{1};
d2 = A{2};
d3 = A{3};

pass(5) = all(d1.deltaLoc == [-.5 -.25 0] );
pass(6) = all(d2.funPart.domain == [0, .5] );
pass(7) = all(d3.deltaMag == 1 );

%% make sure interior delta functions exactly at break points
% get divided equally in each adjacent deltafun.
f = fun.constructor(@(x) exp(x), struct('domain', [-1, 1]));
data.deltaMag = [1 1 1 1 1 1 1];
data.deltaLoc = [-1, -.5, -.25, 0, .25, .5, 1];
d = deltafun(f, data);
A = restrict(d, [-1 0 .5 1]);
d1 = A{1};
d2 = A{2};
d3 = A{3};

pass(8) = all(d1.deltaLoc == [-1 -.5 -.25 0] ) && all(d1.deltaMag == [1 1 1 .5] );
pass(9) = all(d2.deltaLoc == [0 .25 .5] ) && all(d2.deltaMag == [.5 1 .5] );
pass(10) = all(d3.deltaLoc == [.5 1] ) && all(d3.deltaMag == [.5 1] );
end


text:
function [Estimation_X] = RocEKF_propagate(Estimation_X, OdometryFromThis2Next, Sigma_ODO )

% retrieve v and w
v = OdometryFromThis2Next(1:3);
w = OdometryFromThis2Next(4:6);

NumberOfLandmarks = size(Estimation_X.landmarks, 2);
Jrw = jaco_r(-w);
ExpMinusM = so3_exp(-w);

temp = repmat({Estimation_X.orientation}, 2,1 );
adRobot = blkdiag(temp{:});
 adRobot(4:6,1:3)=skew(Estimation_X.position)*Estimation_X.orientation;

B2=zeros(3*NumberOfLandmarks ,  6 );
if NumberOfLandmarks>0
    for i=1:NumberOfLandmarks
    %B2( 3*i-2:3*i  , :  )= ExpMinusM*[  skew( Estimation_X.landmarks(1:3,i) -v  )*Jrw , -eye(3) ];
    B2( 3*i-2:3*i  , :  )= ExpMinusM*[  skew( Estimation_X.landmarks(1:3,i)   ) , -eye(3) ];
    %B2=B2*[eye(3) zeros(3,3); -skew(v) eye(3) ];
    end
end

%B1=-[Jrw  zeros(3,3); skew(v)*Jrw eye(3)];
B1=[Estimation_X.orientation  zeros(3,3); -skew( Estimation_X.orientation*v ) Estimation_X.orientation];
adA= [  B1;  B2    ];
%adA= [  -adRobot* B1; -B2    ];

odoCov=diag([w.^2;v.^2])*Sigma_ODO^2;

temp = repmat({ ExpMinusM  }, 2+NumberOfLandmarks,1 );
AA = blkdiag(temp{:});
AA(1:6,1:6)=eye(6);

% final update the covariance
Estimation_X.cov = AA*Estimation_X.cov*AA'+ adA*odoCov*adA';

% update position and orientation
Estimation_X.position = Estimation_X.position+Estimation_X.orientation*v;
Estimation_X.orientation = Estimation_X.orientation*so3_exp(w);

% update the local coordinates of landmarks
if NumberOfLandmarks>0

landmarks=Estimation_X.landmarks(1:3,:);
vv=repmat(v, 1, NumberOfLandmarks);
landmarks=ExpMinusM*(landmarks - vv  );
Estimation_X.landmarks(1:3,:)=landmarks;

end

end

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
function showgraph(node,edge)
%% SHOWGRAPH displays a planar graph
%
%    showgraph(node,edge) displays a planar undirected graph.
%
%   See also showmesh, findedge
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

dim = size(node,2);
if dim == 2 % 2-D
    line([node(edge(:,1),1)'; node(edge(:,2),1)'],...
         [node(edge(:,1),2)'; node(edge(:,2),2)'],...
         'LineWidth',1,'Color',[0.125 0.5 0.125]);
    hold on
    plot(node(:,1),node(:,2),'k.', 'MarkerSize', 12);
elseif dim == 3 % 3-D
    line([node(edge(:,1),1)'; node(edge(:,2),1)'],...
         [node(edge(:,1),2)'; node(edge(:,2),2)'],...
         [node(edge(:,1),3)'; node(edge(:,2),3)'],...
         'LineWidth',1,'Color',[0.125 0.5 0.125]);
    hold on
    plot3(node(:,1),node(:,2),node(:,3),'k.', 'MarkerSize', 22);
end
axis equal; axis off


text:
function [f,g,H,T] = autoTensor(x,useComplex,funObj,varargin)
% [f,g,H,T] = autoTensor(x,useComplex,funObj,varargin)
% Numerically compute Tensor of 3rd-derivatives of objective function from Hessian values

p = length(x);

if useComplex % Use Complex Differentials
    mu = 1e-150;

    diff = zeros(p,p,p);
    for j = 1:p
        e_j = zeros(p,1);
        e_j(j) = 1;
        [f(j) g(:,j) diff(:,:,j)] = funObj(x + mu*i*e_j,varargin{:});
    end
    f = mean(real(f));
    g = mean(real(g),2);
    H = mean(real(diff),3);
    T = imag(diff)/mu;
else % Use finite differencing
    mu = 2*sqrt(1e-12)*(1+norm(x))/norm(p);

    [f,g,H] = funObj(x,varargin{:});
    diff = zeros(p,p,p);
    for j = 1:p
        e_j = zeros(p,1);
        e_j(j) = 1;
        [junk1 junk2 diff(:,:,j)] = funObj(x + mu*e_j,varargin{:});
    end
    T = (diff-repmat(H,[1 1 p]))/mu;
end

author: huashiyiqike

repo: LSTM-MATLAB


text:
filename='Cantilever_triangle_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'PROJECTED GRADIENT';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.4;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

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

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [x, L]= time_step_RBF_new(W,x,lambda,timestep,time_step_value)

% build the degree function ( d_i = sum_j w_ij)
W=sparse(W);
D=sum(W,2);
laplacian = 1;
[dim num] = size(x);
if nargin<5

            time_step_value=0.25;
            time_step_value=0.1;
            time_step_value=0.25;
             time_step_value=0.1;
              time_step_value=0.2;%% for depth
               time_step_value=0.1;
end

% checks if elements of the degree function are zero and correct this
if(sum(D==0)>0)
    disp('Warning, Elements of d are zero');
    for i=1:num
        if(D(i)==0), D(i)=1; end
    end
end
% type of Laplacian (normalized, unnormalized)

% build the final weight matrix - reweighted by some power of the degree
% function \tilde{k}_ij = k_ij / pow(d_i d_j, lambda)
if(lambda~=0)
    f=spdiags(1./(D.^lambda), 0, num, num);
    W=1/(num)*f*W*f;
    clear f;
end
clear d

% final degree function of weights \tilde{k}
e=sum(W,2);
if(sum(e==0)>0)
    disp('Warning, Elements of e are zero');
    for i=1:num
        if(e(i)==0), e(i)=1/(num); end
    end
end

% for the time step it is easier to have the transpose of x
x=x';

tic
if(timestep==0)
    % explicit timestep
    D=spdiags(e,0,num,num);
    L=D-W;
    step=-0.1*L*x;
    x = x + step;
else
    % implicit timestep
    D=spdiags(e,0,num,num);

    L=sparse(D-W);
    %L=sparse(L);
    switch laplacian
        case 0,
            %normalized, solve L*x_new =x_old, where L=Id - E^{-1}W

            z=(D+0.25*L)\(D*x);
            % z=(D+0.1*L)\(D*x);
            diff = z - x;
            for i=1:num
                x(i,:) = x(i,:) + diff(i,:);
            end
        case 1,
            %unnormalized
            % z=cholmod(speye(num) + 0.5*L,x);

            z = (speye(num) + time_step_value*L)\x;
            diff = z - x;

            parfor i=1:num
                x(i,:) = x(i,:) + diff(i,:);
            end

    end
end

x=x'; % transform the data back in column format (one data point=one column in x)
t=toc;
% disp(['Time for time step: ', num2str(t),' seconds']);
% time_step_value

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% Demo for partially observable oscillatory system.
% This demo inverts a model of a linear oscillatory system, which is
% observed through a nonlinear sigmoid function.

clear variables
close all

% Choose basic settings for simulations
n_t = 2e2;
delta_t = 2e-1;         % integration time step (Euler method)
f_fname = @f_lin2D;
g_fname = @g_sigmoid;

u       = zeros(1,n_t);

% Build options structure for temporal integration of SDE
inF.deltat      = delta_t;
inF.a           = 0.1;
inF.b           = 0.9e-2;
inG.scale = 2;
inG.y0          = -1;
inG.ind         = 1; % only x(1) is partially observable
options.inF     = inF;
options.inG     = inG;

% Parameters of the simulation
alpha   = 1e1;
sigma   = 1e1;
theta   = 1;
phi     = [];

% Build priors for model inversion
priors.muX0 = zeros(2,1);
priors.SigmaX0 = 1e0*eye(2);
priors.muTheta = 0*ones(1,1);
priors.SigmaTheta = 1e0*eye(1);
priors.a_alpha = 1e0;
priors.b_alpha = 1e0;
priors.a_sigma = 1e0;
priors.b_sigma = 1e0;

for t=1:n_t
    priors.iQx{t} = eye(2);
    priors.iQx{t}(1,1) = 1e2;
end

% Build options and dim structures for model inversion
options.priors      = priors;
options.backwardLag  = 16;
dim.n_theta         = 1;
dim.n_phi           = 0;
dim.n               = 2;

% options.checkGrads = 1;

% Build time series of hidden states and observations
[y,x,x0,eta,e] = VBA_simulate (n_t,f_fname,g_fname,theta,phi,u,alpha,sigma,options);

% display time series of hidden states and observations
displaySimulations(y,x,eta,e);
% disp('--paused--')
% pause

% Call inversion routine
% [posterior,out] = VBA_onlineWrapper(y,u,f_fname,g_fname,dim,options);
[posterior,out] = VBA_NLStateSpaceModel(y,u,f_fname,g_fname,dim,options);

% Display results
displayResults(posterior,out,y,x,x0,theta,phi,alpha,sigma)

% Make predictions
try
    options = out.options;
    [xs,ys,xhat,vx,yhat,vy] = VBA_comparePredictions(...
        n_t,theta,phi,u,alpha,sigma,options,posterior,dim);
catch
    disp('------!!Unable to form predictions!!------')
end

author: MBB-team

repo: VBA-toolbox


text:
%

% @
% Copyright (C) 2017 Jonas Ruesch
%
% 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/>.
% @

load('../output/coefficients');

plotCoefficients(alpha, beta, CX, CY, CZ, Cl, Cm, Cn, CX_d, CY_d, CZ_d, Cl_d, Cm_d, Cn_d, ...
    CX_P, CY_P, CZ_P, Cl_P, Cm_P, Cn_P, CX_Q, CY_Q, CZ_Q, Cl_Q, Cm_Q, Cn_Q, CX_R, CY_R, CZ_R, Cl_R, Cm_R, Cn_R);

author: jrgenerative

repo: fixed-wing-sim


text:
function [U,output] = cpd_rnd(size_tens,R,options)
%CPD_RND Pseudorandom initialization for CPD.
%   U = cpd_rnd(size_tens,R) generates pseudorandom factor matrices U{1},
%   ..., U{N} of dimensions size_tens(n)-by-R that can be used to
%   initialize algorithms that compute the CPD of an N-th order tensor.
%
%   cpd_rnd(T,R) is shorthand for cpd_rnd(size(T),R) if T is real. If T is
%   complex, then U{n} will be generated as complex matrices by default
%   (cf. options).
%
%   cpd_rnd(size_tens,R,options) and cpd_rnd(T,R,options) may be used to
%   set the following options:
%
%      options.Real =        - The type of random number generator used to
%      [{@randn}|@rand|0]      generate the real part of each factor
%                              matrix. If 0, there is no real part.
%      options.Imag =        - The type of random number generator used to
%      [@randn|@rand|0|...     generate the imaginary part of each factor
%       {'auto'}]              matrix. If 0, there is no imaginary part.
%                              On 'auto', options.Imag is 0 unless the
%                              first argument is a complex tensor T, in
%                              which case it is equal to options.Real.
%      options.Orth =        - If true, the generated factor matrices are
%      [true|false|{'auto'}]   orthogonalized using a QR factorization.
%                              On 'auto', options.Orth is false if the
%                              first argument is a vector size_tens, and
%                              true if the first argument is a tensor T.
%
%   See also btd_rnd, lmlra_rnd, cpdgen.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)

% Process input.
isSizeVector = isnumeric(size_tens) && isvector(size_tens);
if ~isSizeVector
    T = size_tens;
    if isstruct(T), size_tens = T.size;
    else size_tens = size(T); end
end
N = length(size_tens);

% Check the options structure.
isfunc = @(f)isa(f,'function_handle');
if nargin < 3, options = struct; end
if ~isfield(options,'Real'), options.Real = @randn; end
if ~isfunc(options.Real), options.Real = @zeros; end
if ~isfield(options,'Imag'), options.Imag = 'auto'; end
if ischar(options.Imag) && strcmpi(options.Imag,'auto')
    if ~isSizeVector && ((isstruct(T) && ~isreal(T.val)) || ...
            (isnumeric(T) && ~isreal(T)))
        options.Imag = options.Real;
    else
        options.Imag = 0;
    end
end
if ~isfield(options,'Orth'), options.Orth = 'auto'; end
if ischar(options.Orth) && strcmpi(options.Orth,'auto')
	options.Orth = ~isSizeVector;
end

% Generate factor matrices.
U = arrayfun(@(n)options.Real(size_tens(n),R),1:N,'UniformOutput',0);
if isfunc(options.Imag)
    Ui = arrayfun(@(n)options.Imag(size_tens(n),R),1:N,'UniformOutput',0);
    U = cellfun(@(ur,ui)ur+ui*1i,U,Ui,'UniformOutput',0);
end
for n = 1:N*options.Orth
    if size(U{n},1) >= size(U{n},2), [U{n},~] = qr(U{n},0);
    else [Q,~] = qr(U{n}.',0); U{n} = Q.'; end
end
output = struct;

author: andrewssobral

repo: lrslibrary


text:
function D = spm_mrdivide(A, B)
% Regularised variant of mrdivide(A, B) or A / B, similar to B * spm_inv(A)
% FORMAT D = spm_mrdivide(A, B)
%
% D = B * inv(A), or if A is near singular D = B * inv(A + TOL*eye(size(A))
%
% where TOL is adaptively increased if necessary.
%
% This function should be preferable to B * spm_inv(A) if A is large and
% sparse or if B has few rows, since the inverse need not be explicitly
% computed (the linear system can be solved with the backslash operator).
%
% See also: spm_mldivide
%__________________________________________________________________________
% Copyright (C) 2011 Wellcome Trust Centre for Neuroimaging

% Ged Ridgway
% $Id: spm_mrdivide.m 4360 2011-06-14 16:46:37Z ged $

D = spm_mldivide(B', A')';

author: spm


text:
function krn = spm_smoothkern(fwhm,x,t)
% Generate a Gaussian smoothing kernel
% FORMAT krn = spm_smoothkern(fwhm,x,t)
% fwhm - full width at half maximum
% x    - position
% t    - either 0 (nearest neighbour) or 1 (linear).
%        [Default: 1]
%
% krn  - value of kernel at position x
%__________________________________________________________________________
%
% For smoothing images, one should really convolve a Gaussian with a sinc
% function. For smoothing histograms, the kernel should be a Gaussian
% convolved with the histogram basis function used. This function returns
% a Gaussian convolved with a triangular (1st degree B-spline) basis
% function (by default). A Gaussian convolved with a hat function (0th
% degree B-spline) can also be returned.
%__________________________________________________________________________
% Copyright (C) 2005-2011 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_smoothkern.m 4419 2011-08-03 18:42:35Z guillaume $

if nargin<3, t = 1; end

% Variance from FWHM
s = (fwhm/sqrt(8*log(2)))^2+eps;

% The simple way to do it. Not good for small FWHM
% krn = (1/sqrt(2*pi*s))*exp(-(x.^2)/(2*s));

if t==0
    % Gaussian convolved with 0th degree B-spline
    % int(exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t= -0.5..0.5)
    w1  = 1/sqrt(2*s);
    krn = 0.5*(erf(w1*(x+0.5))-erf(w1*(x-0.5)));
    krn(krn<0) = 0;

elseif t==1
    % Gaussian convolved with 1st degree B-spline
    %  int((1-t)*exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t= 0..1)
    % +int((t+1)*exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t=-1..0)
    w1  =  0.5*sqrt(2/s);
    w2  = -0.5/s;
    w3  = sqrt(s/2/pi);
    krn = 0.5*(erf(w1*(x+1)).*(x+1) + erf(w1*(x-1)).*(x-1) - 2*erf(w1*x   ).* x)...
          +w3*(exp(w2*(x+1).^2)     + exp(w2*(x-1).^2)     - 2*exp(w2*x.^2));
    krn(krn<0) = 0;

else
    error('Only defined for nearest neighbour and linear interpolation.');
    % If anyone knows a nice formula for a sinc function convolved with a
    % a Gaussian, then that could be quite useful.
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ algorithmJob ] = jobInit()
% Construct an instance of an algorithm analyser job

    % Read configuration
    load ../GUI/config/config.mat;
    algorithmJob.algorithmId    = 1;
    algorithmJob.datasetId      = 1;
    AL                          = char(algorithmList);
    algorithmJob.algorithm      = AL(algorithmJob.algorithmId,:);
    DL                          = char(dataList);
    algorithmJob.dataset        = DL(algorithmJob.datasetId,:);
    algorithmJob.parameters     = cell2mat(defaultParameters(algorithmJob.algorithmId,:));
end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function out = SY_RangeEvolve(y)
% SY_RangeEvolve    How the time-series range changes across time.
%
% Measures of the range of the time series as a function of time,
% i.e., range(x_{1:i}) for i = 1, 2, ..., N, where N is the length of the time
% series.
%
%---INPUT:
% y, the time series
%
%---OUTPUTS: based on the dynamics of how new extreme events occur with time.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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/>.
% ------------------------------------------------------------------------------

doPlot = false; % whether to plot outputs
N = length(y); % length of the time series
cums = zeros(N,1);

for i = 1:N
    cums(i) = range(y(1:i));
end
% cums = cums/range(y);

if doPlot
    figure('color','w');
    plot(cums);
end

fullr = range(y);

% Define an anonymous function for the number of unique entries in a vector, x:
lunique = @(x) length(unique(x));

out.totnuq = lunique(cums);

% How many of the unique extrema are in first <proportion> of time series?
cumtox = @(x) lunique(cums(1:floor(N*x)))/out.totnuq;
out.nuqp1 = cumtox(0.01);
out.nuqp10 = cumtox(0.1);
out.nuqp20 = cumtox(0.2);
out.nuqp50 = cumtox(0.5);

% (**1**) how many unique extrema are in first <length> of time series

Ns = [10, 50, 100, 1000];
for i = 1:length(Ns)
    if N >= Ns(i)
        out.(sprintf('nuql%u',Ns(i))) = lunique(cums(1:Ns(i)))/out.totnuq;
    else
        out.(sprintf('nuql%u',Ns(i))) = NaN;
    end
end

% (**2**) Actual proportion of full range captured at different points

out.p1 = cums(ceil(N*0.01))/fullr;
out.p10 = cums(ceil(N*0.1))/fullr;
out.p20 = cums(ceil(N*0.2))/fullr;
out.p50 = cums(ceil(N*0.5))/fullr;

Ns = [10, 50, 100, 1000];
for i = 1:length(Ns)
    if N >= Ns(i)
        out.(sprintf('l%u',Ns(i))) = cums(Ns(i))/fullr;
    else
        out.(sprintf('l%u',Ns(i))) = NaN;
    end
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function C = mtimes(A,B)
%MTIMES Multiplies two tenmat objects.
%
%  C = MTIMES(A,B) computes the product of A and B. The result is a
%  TENMAT object and can be transformed into a tensor.
%
%  C = MTIMES(A,B) is called for the syntax 'A * B' when A or B is a
%  TENMAT object.
%
%   See also TENMAT.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

% Handle scalar input
if ~isa(B,'tenmat') && numel(B) == 1
    C = A;
    C.data = C.data * B;
    return;
end
if ~isa(A,'tenmat') && numel(A) == 1
    C = B;
    C.data = C.data * A;
    return;
end

% Handle matrix input
if ~isa(A,'tenmat')
    A = tenmat(A,1);
end

if ~isa(B,'tenmat')
    B = tenmat(B,1);
end

% Error check
if size(A,2) ~= size(B,1)
    error(['Size mismatch: Number of columns in A is not equal to' ...
	   ' the number of rows in B']);
end

tsiz = [A.tsize(A.rindices) B.tsize(B.cindices)];

if ~isempty(tsiz)
    C = tenmat;
    C.tsize = tsiz;
    C.rindices = 1:length(A.rindices);
    C.cindices = (1:length(B.cindices)) + length(A.rindices);
    C.data = A.data * B.data;
else
    C = A.data * B.data;
end

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [X_out] = etrf2000itrf2008(X_in, date)

% SYNTAX:
%   [X_out] = etrf2000itrf2008(X_in, date);
%
% INPUT:
%   X_in = input coordinates (XYZ)
%   date = reference date
%
% OUTPUT:
%   X_out = output coordinates (XYZ)
%
% DESCRIPTION:
%   ETRF2000 to ITRF2008 coordinate converter.

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

%decimal year (YYYY.DDD)
[~, frac] = date2doy(datenum(date));
t = date(1) + frac;

%ITRF2008 - ETRF2000 parameters
% (http://etrs89.ensg.ign.fr/memo-V8.pdf)

%translation [mm]
T1 =  52.1; T1dot =  0.1;
T2 =  49.3; T2dot =  0.1;
T3 = -58.5; T3dot = -1.8;

%scale factor
D = 1.34e-9; Ddot = 0.08e-9;

%rotation [mas]
R1 =  0.891; R1dot =  0.081;
R2 =  5.390; R2dot =  0.490;
R3 = -8.712; R3dot = -0.792;

%propagate parameters to epoch t
T1 = T1 + T1dot*(t - 2000.0);
T2 = T2 + T2dot*(t - 2000.0);
T3 = T3 + T3dot*(t - 2000.0);
D  = D  + Ddot *(t - 2000.0);
R1 = R1 + R1dot*(t - 2000.0);
R2 = R2 + R2dot*(t - 2000.0);
R3 = R3 + R3dot*(t - 2000.0);

%convert translation parameters to [m]
T1 = T1 * 1e-3;
T2 = T2 * 1e-3;
T3 = T3 * 1e-3;

%convert rotation parameters to [rad]
R1 = R1 * 4.848136e-9;
R2 = R2 * 4.848136e-9;
R3 = R3 * 4.848136e-9;

%translation vector
T = [T1; T2; T3];

%rotation/scale matrix
R = [1 -R3 R2; R3 1 -R1; -R2 R1 1];

%7-parameters Helmert inverse transformation
X_out = inv(R)/(1+D)*(X_in - T);

author: goGPS-Project


text:
% Studies the recovery probability of different algorithms
% at a fixed sparsity level and varying number of signals.

close all;
clear all;
clc;
rng('default');
% Create the directory for storing images
[status_code,message,message_id] = mkdir('bin');
% Signal space
N = 256;
% Number of measurements
M = 48;
% Number of signals
Ss = 1:32;
% Sparsity level
K = 16;
num_trials = 100;
num_ss = length(Ss);
success_with_s.somp = zeros(num_ss, 1);
success_with_s.ra_omp = zeros(num_ss, 1);
success_with_s.ra_ormp = zeros(num_ss, 1);
success_with_s.cosamp_mmv = zeros(num_ss, 1);
success_with_s.ra_cosamp = zeros(num_ss, 1);

snr_threshold = 100;

for ns=1:num_ss
    % Current number of signals
    S = Ss(ns);
    num_successes.somp = 0;
    num_successes.ra_omp = 0;
    num_successes.ra_ormp = 0;
    num_successes.cosamp_mmv = 0;
    num_successes.ra_cosamp = 0;
    for nt=1:num_trials
        % Sensing matrix
        Phi = spx.dict.simple.gaussian_dict(M, N);
        % Sparse signal generator
        gen  = spx.data.synthetic.SparseSignalGenerator(N, K, S);
        % Gaussian distributed non-zero samples
        X = gen.gaussian;
        % Measurement vectors
        Y = Phi * X;

        % Create the solver for simultaneous orthogonal matching pursuit
        solver = spx.pursuit.joint.OrthogonalMatchingPursuit(Phi, K, 2);
        result = solver.solve(Y);
        % Solution vectors
        X_Rec = result.Z;
        % Comparison
        cs = spx.commons.SparseSignalsComparison(X, X_Rec, K);
        % Reconstruction SNR
        snr = cs.cum_signal_to_noise_ratio;
        success = snr > snr_threshold;
        num_successes.somp = num_successes.somp + success;

        % Create the solver for rank aware orthogonal matching pursuit
        solver = spx.pursuit.joint.RankAwareOMP(Phi, K);
        result = solver.solve(Y);
        % Solution vectors
        X_Rec = result.Z;
        % Comparison
        cs = spx.commons.SparseSignalsComparison(X, X_Rec, K);
        % Reconstruction SNR
        snr = cs.cum_signal_to_noise_ratio;
        success = snr > snr_threshold;
        num_successes.ra_omp = num_successes.ra_omp + success;

        % Create the solver for rank aware order recursive matching pursuit
        solver = spx.pursuit.joint.RankAwareORMP(Phi, K);
        result = solver.solve(Y);
        % Solution vectors
        X_Rec = result.Z;
        % Comparison
        cs = spx.commons.SparseSignalsComparison(X, X_Rec, K);
        % Reconstruction SNR
        snr = cs.cum_signal_to_noise_ratio;
        success = snr > snr_threshold;
        num_successes.ra_ormp = num_successes.ra_ormp + success;

        % Create the solver for CoSaMP MMV
        solver = spx.pursuit.joint.CoSaMP(Phi, K);
        result = solver.solve(Y);
        % Solution vectors
        X_Rec = result.Z;
        % Comparison
        cs = spx.commons.SparseSignalsComparison(X, X_Rec, K);
        % Reconstruction SNR
        snr = cs.cum_signal_to_noise_ratio;
        success = snr > snr_threshold;
        num_successes.cosamp_mmv = num_successes.cosamp_mmv + success;

        % Create the solver for Rank Aware CoSaMP MMV
        solver = spx.pursuit.joint.CoSaMP(Phi, K);
        solver.RankAwareResidual = true;
        result = solver.solve(Y);
        % Solution vectors
        X_Rec = result.Z;
        % Comparison
        cs = spx.commons.SparseSignalsComparison(X, X_Rec, K);
        % Reconstruction SNR
        snr = cs.cum_signal_to_noise_ratio;
        success = snr > snr_threshold;
        num_successes.ra_cosamp = num_successes.ra_cosamp + success;

        fprintf('S: %d, K=%d, trial=%d\n', S, K, nt);
    end
    success_with_s.somp(ns) = num_successes.somp / num_trials;
    success_with_s.ra_omp(ns) = num_successes.ra_omp / num_trials;
    success_with_s.ra_ormp(ns) = num_successes.ra_ormp / num_trials;
    success_with_s.cosamp_mmv(ns) = num_successes.cosamp_mmv / num_trials;
    success_with_s.ra_cosamp(ns) = num_successes.ra_cosamp / num_trials;
end

save ('bin/success_with_s_comparison.mat');

author: indigits

repo: sparse-plex


text:
function K = lfmaXlfmKernCompute(lfmKern1, lfmKern2, t1, t2)

% LFMAXLFMKERNCOMPUTE Acceleration and position LFM kernel
% FORMAT
% DESC computes cross kernel terms between acceleration and position LFM
% kernels.
% ARG lfmKern1 : the kernel structure associated with the first LFM
% kernel.
% ARG lfmKern2 : the kernel structure associated with the second LFM
% kernel.
% ARG t : inputs for which kernel is to be computed.
% RETURN K : block of values from kernel matrix.
%
% FORMAT
% DESC computes cross kernel terms between two LFM kernels for
% the multiple output kernel.
% ARG lfmKern1 : the kernel structure associated with the first LFM
% kernel.
% ARG lfmKern2 : the kernel structure associated with the second LFM
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% RETURN K : block of values from kernel matrix.
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin < 4
    t2 = t1;
end
if size(t1, 2) > 1 || size(t2, 2) > 1
    error('Input can only have one column');
end

if lfmKern1.inverseWidth ~= lfmKern2.inverseWidth
    error('Kernels cannot be cross combined if they have different inverse widths.')
end

% Get length scale out.
sigma2 = 2/lfmKern1.inverseWidth;
sigma = sqrt(sigma2);

% Parameters of the kernel
alpha(1) = lfmKern1.damper./(2*lfmKern1.mass);
alpha(2) = lfmKern2.damper./(2*lfmKern2.mass);
omega(1) = sqrt(lfmKern1.spring./lfmKern1.mass - alpha(1)*alpha(1));
omega(2) = sqrt(lfmKern2.spring./lfmKern2.mass - alpha(2)*alpha(2));

% Precomputations to increase the speed
preExp1 = zeros(length(t1),2);
preExp2 = zeros(length(t2),2);
gamma1_p = alpha(1) + j*omega(1);
gamma1_m = alpha(1) - j*omega(1);
gamma2_p = alpha(2) + j*omega(2);
gamma2_m = alpha(2) - j*omega(2);
preGamma(1) = gamma1_p + gamma2_p;
preGamma(2) = gamma1_p + gamma2_m;
preGamma(3) = gamma1_m + gamma2_p;
preGamma(4) = gamma1_m + gamma2_m;
preConst = 1./preGamma;
preFactors(1) = preConst(2) - preConst(1);
preFactors(2) = preConst(3) - preConst(4);
preFactors(3) = preConst(3) - preConst(1);
preFactors(4) = preConst(2) - preConst(4);
preExp1(:,1) = (gamma1_p^2).*exp(-gamma1_p*t1);
preExp1(:,2) = (gamma1_m^2).*exp(-gamma1_m*t1);
preExp2(:,1) = exp(-gamma2_p*t2);
preExp2(:,2) = exp(-gamma2_m*t2);
% Actual computation of the kernel
sK = lfmComputeH3AP(gamma1_p, gamma1_m, sigma2, t1,t2,preFactors([1 2]), 0) + ...
    lfmComputeH3AP(gamma2_p, gamma2_m, sigma2, t2,t1,preFactors([3 4]), 1).' + ...
    lfmComputeH4AP(gamma1_p, gamma1_m, sigma2, t1, preGamma([1 2 4 3]), preExp2, 0 ) + ...
    lfmComputeH4AP(gamma2_p, gamma2_m, sigma2, t2, preGamma([1 3 4 2]), preExp1, 1 ).';

if lfmKern1.isNormalised
    K0 =  lfmKern1.sensitivity*lfmKern2.sensitivity/(8*sqrt(2)*lfmKern1.mass*lfmKern2.mass*prod(omega));
else
    K0 =  sigma*sqrt(pi)*lfmKern1.sensitivity*lfmKern2.sensitivity/(8*lfmKern1.mass*lfmKern2.mass*prod(omega));
end

K = K0*sK;


text:
classdef PatternSearchUseParallelEnum < matlab.mixin.SetGet
    %PatternSearchUseParallelEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        UseParallel('Compute Functions in Parallel',true)
        DoNotUseParallel('Compute Functions in Serial',false);
    end

    properties
        name char = '';
        optionVal(1,1) logical = true
    end

    methods
        function obj = PatternSearchUseParallelEnum(name, optionVal)
            obj.name = name;
            obj.optionVal = optionVal;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('PatternSearchUseParallelEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('PatternSearchUseParallelEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('PatternSearchUseParallelEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function Ath = threshold_components(A,options)

% post processing of spatial components
% for each component perform the following:
%   (i)     perform median filtering
%   (ii)    keep only pixels that contibute up to a level of total energy
%   (iii)   perform morphological closing
%   (iv)    extract largest connected component

% Written by:
% Eftychios A. Pnevmatikakis, Simons Foundation, 2015

    defoptions.nrgthr = 0.9999;              % energy threshold
    defoptions.clos_op = strel('square',3);  % morphological operator for closing
    defoptions.medw = [3,3];                 % size of median filter

    if ~isfield(options,'nrgthr') || isempty(options.nrgthr); options.nrgthr = defoptions.nrgthr; end
    if ~isfield(options,'clos_op') || isempty(options.clos_op); options.clos_op = defoptions.clos_op; end
    if ~isfield(options,'medw') || isempty(options.medw); options.medw = defoptions.medw; end
    if ~isfield(options,'d3') || isempty(options.d3); options.d3 = 1; end

    [d,nr] = size(A);
    Ath = spalloc(d,nr,nnz(A));
    Ath(:,nr-options.nb+1:nr) = A(:,nr-options.nb+1:nr);
    indf = cell(nr,1);
    valf = cell(nr,1);
    parfor i = 1:nr-options.nb
        A_temp = reshape(full(A(:,i)),options.d1,options.d2,options.d3);
        for z = 1:options.d3
            A_temp(:,:,z) = medfilt2(A_temp(:,:,z),options.medw);
        end
        A_temp = A_temp(:);
        [temp,ind] = sort(A_temp(:).^2,'ascend');
        temp =  cumsum(temp);
        ff = find(temp > (1-options.nrgthr)*temp(end),1,'first');
        BW = zeros(options.d1,options.d2,options.d3);
        BW(ind(ff:d)) = 1;
        for z = 1:options.d3
            BW(:,:,z) = imclose(BW(:,:,z),options.clos_op);
        end
        [L,NUM] = bwlabeln(BW,8*(options.d3==1) + 6*(options.d3~=1));
        if NUM > 0
            nrg = zeros(NUM,1);
            for l = 1:NUM
                ff = (L==l);
                nrg(l) = sum(A_temp(ff).^2);
            end
            [~,indm] = max(nrg);
            ff = find(L==indm);
            %Ath(ff,i) = A(ff,i);
            indf{i} = ff;
            valf{i} = A_temp(ff);
        else
            valf{i} = 0;
        end
    end
    for i = 1:nr-options.nb
        Ath(indf{i},i) = valf{i};
    end
end

author: zhoupc