text:
function a = c8vec_unity ( n )

%*****************************************************************************80
%
%% C8VEC_UNITY returns the N roots of unity.
%
%  Discussion:
%
%    A(1:N) = exp ( 2 * PI * (0:N-1) / N )
%
%    A(1:N)^N = ( (1,0), (1,0), ..., (1,0) ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of elements of A.
%
%    Output, complex A(N,1), the N roots of unity.
%
  x = 0 : 2 : 2 * n - 2;
  theta = pi * ( x' ) / n;
  a(1:n,1) = cos ( theta ) + i * sin ( theta );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [aVal,aJacob,aHess,papt]=aCVPolar(x)
%%ACVPOLAR The drift function for a continuous-time motion model where the
%          target state is given in polar coordinates and the motion is
%          constant velocity in 2D Cartesian coordinates. This function
%          works regardless of whether the angle is measured
%          counterclockwise from the x axis or clockwise form the y axis.
%
%INPUTS: x The 4XN state vector of N targets in the order of
%          [r;theta;rDot;thetaDot]. The angle theta is given in radians.
%
%OUTPUTS: aVal The 4XN set of time-derivatives of the N state vectors under
%              the Cartesian linear motion model in polar coordinates.
%       aJacob This and higher partial derivatives can only be requested
%              if N=1. This is the 4X4  matrix of partial derivatives
%              of aVal such that aJacob(:,i) is the partial derivative of
%              aVal with respect to x(i).
%        aHess The 4X4X4  matrix of second derivatives of aVal such
%              that aHess(:,k1,k2) is the second partial derivative of
%              aVal with respect to x(k1) and x(k2).
%         papt The 4X1  partial derivative with resect to time of aVal.
%              This is all zeros, because the model is time invariant.
%
%A full derivation of the dynamic model is provided in [1]. Summarizing it
%here, let rVec be a position vector. We shall use the orthonormal basis
%vectors
%u_r=[cos(theta);sin(theta)] and u_theta=[-sin(theta);cos(theta)]
%if one is measuring theta from the x axis counter clockwise and
%u_r=[sin(theta);cos(theta)] and u_theta=[cos(theta);-sin(theta)]
%if one is measuring theta from the y axis clockwise.
%Note that u_theta is the derivative of u_r with respect to theta and the
%two vectors are orthonormal.
%A position can be expressed
%rVec=r*u_r
%Note the identities for the derivatives (dot terms):
%uDot_r=thetaDot*du_r/dtheta=dTheta*u_theta
%and uDot_theta=-thetaDor*u_r
%The velocity vector is
%rVecDot=rDot*u_r+r*uDot_r=rDot*u_r+r*thetaDot*u_theta
%Taking one more derivative to get an acceleration vectors, one has, after
%simplification:
%rDDot=(rDDot-r*thetaDot^2)*u_r+(r*thetaDDot+2*rDot*thetaDot)*u_theta
%where the coefficient of u_r is the radial acceleration term and the
%coefficient of u_theta is the angular acceleration term. Noting that for
%linear motion both accelerations must be zero, one gets the following two
%equations for the second derivative terms:
%rDDot=r*thetaDot^2
%thetaDDot=-(2/r)*rDot*thetaDot
%Thus, one has a linear dynamic model in 2D polar coordinates.
%
%EXAMPLE 1:
%Here, we verify that integrating forward with this model is equivalent to
%linear motion in Cartesian coordinates.
% xInitCart=[1000;40;-100;20];
% T=50;%The prediction time.
% F=FPolyKal(T,4,1);
% systemType=0;
% xInitSpher=stateCart2Pol(xInitCart,systemType);
% xEndCart=F*xInitCart;
% RelTol=1e-10;
% AbsTol=1e-13;
% xStepsPolar=RKAdaptiveOverRange(xInitSpher,[0;T],@(x,t)aCVPolar(x),0.1,0,[],[],RelTol,AbsTol);
% xEndPolarRK=xStepsPolar(:,end);
% xEndPolarExact=stateCart2Pol(xEndCart,systemType);
% max(abs(xEndPolarRK-xEndPolarExact))
%One will observe that the error is less than 1e-7, which is a good
%agreement.
%
%EXAMPLE 2:
%Here, we verify that the Jacobian is consistent with the numerical
%derivative of aDeriv.
% x=[1000;-0.1;10;0.01];
% [aDeriv,aJacob]=aCVPolar(x);
% AJacobNumDiff=numDiff(x,@(xState)aCVPolar(xState),4);
% err=(aJacob-AJacobNumDiff)./AJacobNumDiff;
% max(abs(err(:)))
%One will see that the maximum error is on the order of 1.6e-11, indicating
%good agreement.
%
%REFERENCES:
%[1] D. F. Crouse, "Basic Linear Dynamic Models in Local Coordinates,"
%    Naval Research Laboratory, Washington, D.C., Tech. Rep.
%    NRL/MR/5344--19-9882, 24 Aug. 2019.
%
%August 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

r=x(1,:);
rDot=x(3,:);
thetaDot=x(4,:);

rDDot=r.*thetaDot.^2;
thetaDDot=-(2./r).*rDot.*thetaDot;

aVal=[rDot;thetaDot;rDDot;thetaDDot];

if(nargout>1)
    N=size(x,2);

    if(N>1)
        error('Derivatives are only available for numPoints=1.')
    end

    aJacob=[0                    0,  1,                 0;
            0,                   0,  0,                 1;
            thetaDot^2,          0,  0,                 2*r*thetaDot;
            2*rDot*thetaDot/r^2, 0,  -((2*thetaDot)/r), -((2*rDot)/r)];

    if(nargout>2)
        aHess=zeros(4,4,4);

        aHess(:,:,1)=[zeros(2,4);
                      0,                        0,  0,                2*thetaDot;
                      -((4*rDot*thetaDot)/r^3), 0,  (2*thetaDot)/r^2, (2*rDot)/r^2];
        aHess(:,:,3)=[zeros(3,4);
                      2*thetaDot/r^2,0,0,-(2/r)];
        aHess(:,:,4)=[zeros(2,4);
                      2*thetaDot, 0,  0,        2*r;
                      2*rDot/r^2,   0,  -((2)/r), -((2*rDot)/r)];

        if(nargout>3)
            papt=zeros(4,1);
        end
    end
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:
%VGG_SEGMENT_GB  Graph-based image segmentation
%
%   S = vgg_segment_gb(A, sigma, K, min_sz[, compress])
%
% Segmentation of an image using the graph-based method described in:
%   "Efficient Graph-Based Image Segmentation.", Pedro F. Felzenszwalb and
%   Daniel P. Huttenlocher. International Journal of Computer Vision,
%   Volume 59, Number 2, September 2004.
% Code downloaded from:
%   http://people.cs.uchicago.edu/~pff/segment/
%
%IN:
%   A - HxWx3 uint8 image for segmentation.
%   sigma - scalar parameter on smoothing kernel to use prior to
%           segmentation.
%   k - scalar parameter on prefered segment size.
%   min_sz - scalar indicating the minimum number of pixels per segment.
%   compress - scalar boolean indicating whether the user wants the segment
%              indices compressed to the range [1 num_segments].
%              Default: 0.
%
%OUT:
%   S - HxW uint32 segmentation matrix, each value of which gives the index
%       of the region said pixel belongs to.

% $Id: vgg_segment_gb.m,v 1.1 2007/12/10 10:59:30 ojw Exp $

function varargout = vgg_segment_gb(varargin)
funcName = mfilename;
sd = 'seg_gb/';
sourceList = {['-I' sd], [funcName '.cxx']};
vgg_mexcompile_script; % Compilation happens in this script
return

repo: SelfExSR


text:
%TRIDEIG Eigenvalues of symmetric tridiagonal matrix.
%   E = TRIDEIG(A,B) is a vector containing the eigenvalues of a square
%   tridiagonal matrix with diagonal elements A and subdiagonal elements B.
%   A and B should be real vectors of length N and N-1, respectively, where
%   N is the size of the matrix. E is a real vector of length N.

%   P.-O. Persson <persson@mit.edu>
%   Department of Mathematics, MIT
%   October 21, 2002

error('No mex-file found for this platform.')

repo: Aero-matlab


text:
function model = VGG16_for_Fast_RCNN_VOC2007(model)
% VGG 16layers (only finetuned from conv3_1)

model.solver_def_file        = fullfile(pwd, 'models', 'fast_rcnn_prototxts', 'vgg_16layers_conv3_1', 'solver_30k40k.prototxt');
model.test_net_def_file      = fullfile(pwd, 'models', 'fast_rcnn_prototxts', 'vgg_16layers_conv3_1', 'test.prototxt');

model.net_file               = fullfile(pwd, 'models', 'pre_trained_models', 'vgg_16layers', 'vgg16.caffemodel');
model.mean_image             = fullfile(pwd, 'models', 'pre_trained_models', 'vgg_16layers', 'mean_image');

end

author: ShaoqingRen


text:
function resetTimeTicks(h, num, format)
% reset function for setTimeTicks
% SEE ALSO: setTimeTicks

%--------------------------------------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    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
%--------------------------------------------------------------------------

    if (nargin == 1)
        num = 4;
        format = 'yyyy/mm/dd HH:MM';
    end
    if isa(h, 'matlab.ui.Figure')
        h = get(h,'Children');
    end
    for i = 1 : numel(h)
        try
            h(i).TickLength = [0.005 0.01];
            ax = double(axis(h(i)));

            step = (ax(2)-ax(1))/(num);
            time_span = ax(2)-ax(1);

            % depending on the time_span round ticks
            round_val = 1 / 86400; % one second

            if time_span > 8                 % greater than 20days
                round_val =  1;               % round on 1 day
            elseif time_span > 4              % greater than 5 days
                round_val =  12 / 24;         % round on 12 hour
            elseif time_span > 1              % greater than a day
                round_val =  1 / 24;          % round on 1 hour
            elseif time_span > 12 / 24        % greater than 12 hours
                round_val =  30 / (24 * 60);  % round on 30 minutes
            elseif time_span > 5 / 24         % greater than 5 hours
                round_val =  15 / (24 * 60);  % round on 15 minutes
            elseif time_span > 1 / 24         % greater than 1 hour
                round_val =  5 / (24 * 60);   % round on 5 minutes
            elseif time_span > 10 / (24 * 60) % greater than 10 minutes
                round_val =  1 / (24 * 60);   % round on 1 minute
            elseif time_span > 5 / (24 * 60)  % greater than 5 minutes
                round_val =  30 / 86400;      % round on 30 seconds
            elseif time_span > 1 / (24 * 60)  % greater than 1 minutes
                round_val =  5 / 86400;       % round on 5 seconds
            end

            tick_pos = (ax(1) + (step/2)) : step : (ax(2) - (step/2));

            c_out = 0;
            if round_val >= 1 / (24 * 60) % rounding is greater than 1 minute
                c_out = 3;                % do not show seconds
            elseif round_val >= 1 / 24    % rounding is greater than 1 hour
                c_out = 6;                % do not show minutes
            elseif round_val >= 1         % greater than 1 day
                c_out = 9;                % do not show hours
            end

            tick_pos = round(tick_pos ./ round_val) .* round_val; % round ticks
            tick_pos((tick_pos < ax(1)) | (tick_pos > ax(2))) = []; % delete ticks outside figure;
            tick_pos = unique(tick_pos);

            if all(mod(tick_pos, 1) == 0)
                round_val = 1;
            end
            set(h(i), 'XTick', tick_pos);
            if strcmp(format, 'auto')
                last_date = [0 0 0 0 0 0];
                h.XTickLabel = cell(size(tick_pos));

                x_labels = cell(size(tick_pos));
                for l = 1 : numel(tick_pos)
                    str_time_format = 'HH:MM:SS';
                    new_date = datevec(double(tick_pos(l)));
                    if last_date(1) ~= new_date(1)        % if the year is different
                        str_time_format = [char(32 * ones(1, floor(c_out/2) + 5)) str_time_format(1:end - c_out) '-newlinemmm dd, yyyy'];
                    elseif last_date(2) ~= new_date(2)    % if the month is different
                        str_time_format = [char(32 * ones(1, floor(c_out/2) + 5)) str_time_format(1:end - c_out) '-newlinemmm dd, yyyy'];
                    elseif last_date(3) ~= new_date(3)    % if the day is different
                        str_time_format = [char(32 * ones(1, floor(c_out/2) + 5)) str_time_format(1:end - c_out) '-newlinemmm dd, yyyy'];
                    elseif last_date(4) ~= new_date(4)    % if the hour is different
                        str_time_format = str_time_format(1:end - c_out);
                    else % if last_date(5) ~= new_date(5)    % if the hour is different
                        str_time_format = str_time_format(1:end - c_out);
                    end
                    if isempty(str_time_format)
                        h.XTickLabel{l} = '';
                    else
                        if round_val == 1
                            tmp = datestr(tick_pos(l), regexp(str_time_format, '(?<=line).*','match', 'once'));
                        else
                            tmp = datestr(tick_pos(l), str_time_format);
                            tmp(tmp == '-') = '\';
                        end
                        h.XTickLabel{l} = tmp;
                    end
                    last_date = new_date;
                end
            else
                datetick(h(i),'x',format,'keepticks');
            end
            axis(h(i),ax);
        catch ex
            Core_Utils.printEx(ex);
        end
    end
end

author: goGPS-Project


text:
function t = bal_seq_unrank ( rank, n )

%*****************************************************************************80
%
%% BAL_SEQ_UNRANK unranks a balanced sequence.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer RANK, the rank of the balanced sequence.
%
%    Input, integer N, the number of 0's (and 1's) in the sequence.
%    N must be positive.
%
%    Output, integer T(2*N), a balanced sequence.
%

%
%  Check.
%
  if ( n < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BAL_SEQ_UNRANK - Fatal error!\n' );
    fprintf ( 1, '  Input N is illegal.\n' );
    error ( 'BAL_SEQ_UNRANK - Fatal error!' );
  end

  nseq = bal_seq_enum ( n );

  if ( rank < 0 || nseq < rank )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BAL_SEQ_UNRANK - Fatal error!\n' );
    fprintf ( 1, '  The input rank is illegal.\n' );
    error ( 'BAL_SEQ_UNRANK - Fatal error!' );
  end

  y = 0;
  low = 0;

  for x = 1 : 2 * n

    m = mountain ( n, x, y + 1 );

    if ( rank <= low + m - 1 )
      y = y + 1;
      t(x) = 0;
    else
      low = low + m;
      y = y - 1;
      t(x) = 1;
    end

  end

  return
end


text:
function [L,a,b] =rgb2lab(R,G,B)
% function [L, a, b] = RGB2Lab(R, G, B)
% RGB2Lab takes matrices corresponding to Red, Green, and Blue, and
% transforms them into CIELab.  This transform is based on ITU-R
% Recommendation  BT.709 using the D65 white point reference.
% The error in transforming RGB -> Lab -> RGB is approximately
% 10^-5.  RGB values can be either between 0 and 1 or between 0 and 255.
% By Mark Ruzon from C code by Yossi Rubner, 23 September 1997.
% Updated for MATLAB 5 28 January 1998.

if (nargin == 1)
  B = double(R(:,:,3));
  G = double(R(:,:,2));
  R = double(R(:,:,1));
end

if ((max(max(R)) > 1.0) | (max(max(G)) > 1.0) | (max(max(B)) > 1.0))
  R = R/255;
  G = G/255;
  B = B/255;
end

[M, N] = size(R);
s = M*N;

% Set a threshold
T = 0.008856;

RGB = [reshape(R,1,s); reshape(G,1,s); reshape(B,1,s)];

% RGB to XYZ
MAT = [0.412453 0.357580 0.180423;
       0.212671 0.715160 0.072169;
       0.019334 0.119193 0.950227];
XYZ = MAT * RGB;

X = XYZ(1,:) / 0.950456;
Y = XYZ(2,:);
Z = XYZ(3,:) / 1.088754;

XT = X > T;
YT = Y > T;
ZT = Z > T;

fX = XT .* X.^(1/3) + (~XT) .* (7.787 .* X + 16/116);

% Compute L
Y3 = Y.^(1/3);
fY = YT .* Y3 + (~YT) .* (7.787 .* Y + 16/116);
L  = YT .* (116 * Y3 - 16.0) + (~YT) .* (903.3 * Y);

fZ = ZT .* Z.^(1/3) + (~ZT) .* (7.787 .* Z + 16/116);

% Compute a and b
a = 500 * (fX - fY);
b = 200 * (fY - fZ);

L = reshape(L, M, N);
a = reshape(a, M, N);
b = reshape(b, M, N);

if ((nargout == 1) | (nargout == 0))
  L = cat(3,L,a,b);
end

author: zouchuhang

repo: LayoutNet


text:
% This example shows how to calculate the even and odd modes
% of two coupled waveguides, using the semivectorial modesolver

% Refractive indices:
n1 = 3.34;          % Lower cladding
n2 = 3.44;          % Core
n3 = 1.00;          % Upper cladding (air)

% Layer heights:
h1 = 2.0;           % Lower cladding
h2 = 1.3;           % Core thickness
h3 = 0.5;           % Upper cladding
rh = 1.1;           % Ridge height

% Horizontal dimensions:
rw = 1.0;           % Ridge half-width
d = 2.5;            % center-to-center separation
side = 2.5;         % Space on side

% Grid size:
dx = 0.0125;        % grid size (horizontal)
dy = 0.0125;        % grid size (vertical)

lambda = 1.55;      % wavelength
nmodes = 1;         % number of modes to compute

[x,y,xc,yc,nx,ny,eps,edges] = waveguidemeshfull([n1,n2,n3],[h1,h2,h3],...
                                          rh,rw,[(d/2-rw),side],dx,dy);

% First, we calculate the symmetric mode

[Ex1,neff1] = svmodes(lambda,n2,nmodes,dx,dy,eps,'000S','EX');
fprintf(1,'neff(1) = %.6f\n',neff1);

% Next, we calculate the symmetric mode

[Ex2,neff2] = svmodes(lambda,n2,nmodes,dx,dy,eps,'000A','EX');
fprintf(1,'neff(2) = %.6f\n',neff2);

subplot(211);
contourmode(x,y,Ex1);
title('Ex (Symmetric TE Mode)'); xlabel('x'); ylabel('y');
for v = edges, line(v{:}); end

subplot(212);
contourmode(x,y,Ex2);
title('Ex (Antisymmetric TE Mode)'); xlabel('x'); ylabel('y');
for v = edges, line(v{:}); end

author: Sable

repo: mcbench-benchmarks


text:
function [A, B, C, D] = Digital2Analog(Ap, Bp, Cp, Dp)

% Digital2Analog: bilinear transformation of a discrete system P to analog
% domain
%
% Usage:    [A, B, C, D] = Digital2Analog(Ap, Bp, Cp, Dp)
%
% INPUTS:   state space matrices of the input discrete system P
%
% OUTPUT:   state space matrices of the output analog system
%
% See also: designIIR, Analog2Digital

S = inv(eye(size(Ap,1)) + Ap);

A = (Ap - eye(size(Ap,1))) * S;
B = (eye(size(A))-A) * Bp;
C = Cp * S;
D = Dp - C * Bp;

sys = ss(A,B,C,D);
sys = sminreal(sys);

[A, B ,C ,D] = ssdata(sys);


text:
function  summaryWithFgVolume = dtiExtrapolateFgVolumeFromCrossectionArea(summaryWithCrossectionArea, summaryWithFiberGroupLength)
%Compute crosssection area for fiber groups/subjects
%
% summary = dti_Longitude_ComputeCrossectionalArea(summaryWithCrosSectionalArea, ...
%                                                      [summaryWithFiberGroupLength])
%
% Use case 1
% Input: a summary structure of fg properties returned by
% dtiFiberProperties
% Output:  same summary structure with field
% diberGroupVolume added or updated.
%
% Use case 2
% If two summary structures are provided as arguments, then
% the fgVolume will be extracted from the first one, and fiberLength will
% be extracted from the second one. The resulting summary will contain ONLY
% cross-sectional area estimates, and only for subjects present in both
% input summary structures. Note:it is OK if subjects order is not
% equivalent in both input summary structures. It is NOT ok if the fiber
% groups are not -- u should check that!
%
% Example:
%   load('summaryFiberPropertiesMoriGroups_Roi2RoiVolumeUniqueVoxels.mat');
%   summaryWithCrossectionArea = dtiComputeCrossectionArea(summary);
%
% See a more elaborate example after "return" statement below.
%
% See also: dtiFiberProperties, dtiComputeCrossectionalArea
%
% (c) Vistasoft

% History: ER wrote it 03/24/2010

fprintf(1, 'Computing volume of a fiber group from its crossectional area and length \n');

if nargin == 1
    summaryWithFgVolume = summaryWithCrossectionArea;
    summaryWithFiberGroupLength = summaryWithCrossectionArea;
elseif nargin == 2
    summaryWithFgVolume =[];
else
    eror ('Incorrect number of arguments');
end

for s=1:length(summaryWithCrossectionArea)

    s2 = find(arrayfun(@(x) strcmp(x.subject,summaryWithCrossectionArea(s).subject), summaryWithFiberGroupLength));

    for fg=1:numel(summaryWithCrossectionArea(s).sfg);

        summaryWithFgVolume(s).sfg(fg).fiberGroupVolume = summaryWithCrossectionArea(s).sfg(fg).crossectionArea*summaryWithFiberGroupLength(s2).sfg(fg).fiberLength(2);
        summaryWithFgVolume(s).sfg(fg).name = summaryWithCrossectionArea(s).sfg(fg).name;
        summaryWithFgVolume(s).subject = summaryWithCrossectionArea(s).subject;
    end
end

return

%% An 'elaborate' example
cd /biac3/wandell4/users/elenary/longitudinal/ANALYSES/Mori_Groups/

%load 'summary' structure with Roi2Roi cropped fiber groups' properties
load('summaryFiberPropertiesMoriGroups_Roi2RoiVolumeUniqueVoxels.mat');
%Compute crossectional areas
summaryWithCrossectionArea = dtiComputeCrossectionArea(summary);

%load 'summary' with original (full length) FG properties, incl. length
load('summaryFiberPropertiesMoriGroups_volumeUniqueVoxels.mat');
summaryWithFiberGroupLength = summary;
%Compute extrapolated volume from realLength and approximated Xarea
summaryFgVolumeExtra = dtiExtrapolateFgVolumeFromCrossectionArea(summaryWithCrossectionArea, summaryWithFiberGroupLength);

%Now your summary(1).sfg(1).fiberGroupVolume contains actual volume of fg1 for s1,
%and summaryVolumeExtra(1).sfg(1).fiberGroupVolume contains the volume
%of fg1 for s1 estimated as crossectional area of the middle (Roi to
%Roi, more compact) portion of the fg, extrapolated to the its full
%length.

%Now plot some stuff
labels = dtiGetMoriLabels;
labels_short=dtiGetMoriLabels(true);
load('../../data/subjectCodesAll4Years'); %This will load subjectCodes
subjectInitials={'ab', 'ajs', 'am', 'an', 'at', 'clr', 'crb', 'ctb', 'da', 'dh', 'dm', 'es', 'jh', 'lj', 'll', 'mb', 'md', 'mho', 'mn', 'pf', 'pt', 'rd', 'rsh', 'ss', 'tm', 'vr', 'vt', 'zs'};
[fiberMetrics_realVol, subjectCodes, year, group] = dtiGetFGProperties(summaryWithFiberGroupLength, subjectInitials, 'fiberGroupVolume', labels, '0');
[fiberMetrics_realLength, subjectCodes, year, group] = dtiGetFGProperties(summaryWithFiberGroupLength, subjectInitials, 'fiberLength', labels, '0');
[fiberMetrics_XArea, subjectCodes, year, group] = dtiGetFGProperties(summaryWithCrossectionArea, subjectInitials, 'crossectionArea', labels, '0');
[fiberMetrics_extrapVol, subjectCodes, year, group] = dtiGetFGProperties(summaryFgVolumeExtra, subjectInitials, 'fiberGroupVolume', labels, '0');
[fiberMetrics_FA, subjectCodes, year, group] = dtiGetFGProperties(summaryWithFiberGroupLength, subjectInitials, 'FA', labels, '0');

%Need to know numFibers in order to collapse data across somefiber groups.
[numFibers_fullFGs, subjectCodes, year, group] = dtiGetFGProperties(summaryWithFiberGroupLength, subjectInitials, 'numberOfFibers', labels, '0');

%Reduced 20 fiber groups to 9
%collapsingVector2 = [1 2 3 4 5 6 5 6 7 8 9 10 11 12 13 14 15 16 13 14]; %will combine cingulum cingulate and hc, and will combine slf_t and slf_fp.
collapsingVector2 = [1 1 2 2 3 3 3 3 4 5 6 6 7 7 8 8 9 9 8 8]; %will collapse symmetric fiber groups into one, will combine cingulum cingulate and hc, and will combine slf_t and slf_fp.

[fiberMetrics_realVol, labels_New] = dti_LongitudeAggregateFiberPropertiesAcrossGroups(fiberMetrics_realVol, labels_short, collapsingVector2, 'sum');
[fiberMetrics_extrapVol, labels_New] = dti_LongitudeAggregateFiberPropertiesAcrossGroups(fiberMetrics_extrapVol, labels_short, collapsingVector2, 'sum');
[fiberMetrics_realLength, labels_New] = dti_LongitudeAggregateFiberPropertiesAcrossGroups(fiberMetrics_realLength, labels_short, collapsingVector2, 'mean', numFibers_fullFGs);
[fiberMetrics_XArea, labels_New] = dti_LongitudeAggregateFiberPropertiesAcrossGroups(fiberMetrics_XArea, labels_short, collapsingVector2, 'sum');
[fiberMetrics_FA, labels_New] = dti_LongitudeAggregateFiberPropertiesAcrossGroups(fiberMetrics_FA, labels_short, collapsingVector2, 'mean', numFibers_fullFGs);

% Length -> Volume & Extrapolated volume
figure; plot(fiberMetrics_realLength(:), fiberMetrics_extrapVol(:), 'ro'); xlabel('real length'); hold on;
plot(fiberMetrics_realLength(:), fiberMetrics_realVol(:), 'bx');  legend('Extrapolated volume, mm^3', 'Real volume, mm^3');

%figure; plot(log(realVol), log(extrapolatedVol), 'ro'); xlabel('Log (real volume)'); ylabel('Log (extrapolated volume)');
f1=figure;
screen_size = get(0, 'ScreenSize'); set(f1, 'Position', [0 0 screen_size(3) screen_size(4) ] );
colordef black

markers=getUniqueMarkers(20);
for fg=1:length(labels_New)
    ax1 = subplot(2, 3, 1); hold on;
    plot(fiberMetrics_realLength(:, fg), fiberMetrics_XArea(:, fg), markers{fg}); ylabel('Crossectional area ROi2Roi, mm^2'); xlabel('real length');
    ax2 = subplot(2, 3, 2); hold on;
    plot(fiberMetrics_realVol(:, fg), fiberMetrics_extrapVol(:, fg), markers{fg}); ylabel('extrapolated volume, mm^3'); xlabel('real volume, mm^3');
    ax3 = subplot(2,3, 3); hold on;
    plot(fiberMetrics_XArea(:, fg), fiberMetrics_realVol(:, fg), markers{fg}); ylabel('Real volume'); xlabel('crossection area Roi2Roi'); h=legend(labels_New);
    ax4 = subplot(2, 3, 4); hold on;
    plot(fiberMetrics_FA(:, fg), fiberMetrics_realLength(:, fg), markers{fg}); ylabel('realLength');
    ax5 = subplot(2, 3, 5); hold on;
    plot(fiberMetrics_FA(:, fg), fiberMetrics_realVol(:, fg),  markers{fg}); ylabel('realVolume'); xlabel('FA');
    ax6 = subplot(2, 3, 6); hold on; plot(fiberMetrics_FA(:, fg), fiberMetrics_XArea(:, fg), markers{fg}); ylabel('Crossection Area'); xlabel('FA');
end
pL=get(h, 'Position');
p1 = get(ax1, 'Position');
set(h, 'Position', [p1(1) p1(4) pL(3) pL(4)]);
%colordef white

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [results]=VBA_Shapley(posterior,out,varargin)
% [results]=VBA_Shapley(posterior,out,[options])
%
% Compute the Shapley values of the model's factors (inputs or paramters).
% These scores measure the relative influence of these factor on the variaince
% explained by the model.
% -------------------------------------------------------------------------
% IN:
%  - posterior, out: model structures
%  - varargin:       list of option/values pairs (see below)
% OUT:
%  - results: structure of results with the following fields, depending on
%  the chose coalitions
%  .parameters/inputs:
%         Shapley values (first order) of paramters/inputs. This is a n x p
%         matrix where n is the number of factors of interest and p the
%         number of observations in the model.
%  - interaction: Interaction score, ie relative change in Shapley values for a
%         perturbation in respective inputs. It is an array with the same
%         size as the inputs whose elements are similar to sv.
% -------------------------------------------------------------------------
% options:
%   coalitions:         ('parameters') | 'inputs' | 'interactions'
%      > factor of interest of the analysis. For 'interaction', compute the
%      Shapley value of parameters and their relative change for selective
%      input perturbations.
%   inputPerturbation   -> ('zero') | 'average' | 'average_nonzero'
%      > type of perturbation to apply on inputs: set to zero, average
%      accross the experiment, or set non zero inputs to their average
%   paramType           -> ('phi')  | 'theta'
%      > parameters of interest.
%   paramIdx            -> (all)  | paramIdx
%      > restrict parameters of interest to the index list paramIdx
%   inputIdx            -> (all)  | inputIdx
%      > restrict inputs of interest to the index list inputIdx
%   obsIdx              -> (all)  | ObsIdx
%      > restrict observations of interest to the index list obsIdx
% -------------------------------------------------------------------------

%% Complete option structure
% -------------------------------------------------------------------------

if numel(varargin) == 1 && isstruct(varargin{1})
    options = varargin{1};
else
    options.coalitions = {'parameters','inputs','interactions'};
    options.inputPerturbation = {'zero','average','average_nonzero'};
    options.paramType = {'phi','theta'};
    options.paramIdx = 1:out.dim.n_phi ;
    options.inputIdx = 1:out.dim.u ;
    if size(out.y,2) == 1 % catch vertical data as unique observation
        options.obsIdx = 0;
    else
        options.obsIdx = 1:out.dim.p ;
    end
    parser = inputParser;
    parser.parse (varargin{:});
    options = parser.Results;
    %options = parseargs(options,varargin{:});
end

%% Prepare perturbation scheme
% -------------------------------------------------------------------------

% dimensions
nu = numel(options.inputIdx);
nw = numel(options.paramIdx);
if options.obsIdx==0 % catch vertical data as unique observation
    nResps=1;
else
    nResps = numel(options.obsIdx);
end

% factorial perturbations on coalitions of interest
switch options.coalitions
    case 'interactions'
        % factorial perturbation of parmeeters with normal inputs,
        kw = full(VBA_spm_perm_mtx(nw));
        k = [kw ones(2^nw,nu)] ;
        % factorial perturbation of parameters with each input respectively
        % pertrubed
        ku = ones(nu)-eye(nu);
        for i=1:nu
            k = vertcat(k, [kw repmat(ku(i,:),size(kw,1),1)]);
        end
        % plus factorial perturbation of inputs alone
        % (ie 2^nw + nu x 2^nw + 2^nu coalitions)
        k = vertcat(k,[ones(2^nu,nw) full(VBA_spm_perm_mtx(nu))]);
    case 'parameters'
        % factorial perturbation of paramters, normal inputs
        k = full(VBA_spm_perm_mtx(nw));
        k = [k ones(size(k,1),nu)];
    case 'inputs'
        % factorial perturbation of inputs, normal paramters
        k = full(VBA_spm_perm_mtx(nu));
        k = [ones(size(k,1),nw) k];
end
nk = size(k,1);

%% Compute explained variances
% -------------------------------------------------------------------------

% loop over coalitions
ve = nan(nk,nResps);
parfor t = 1:size(k,1)
    kt = k(t,:);
    w_perm = kt(1:nw);
    u_perm = kt(nw+(1:nu));
    ve(t,:) = explainedVar(posterior,out,options,u_perm,w_perm) ;
end

% normalize
ve1 = ve(1,:);
ve0 = ve(end,:);
for i=1:nResps
    ve(:,i) = (ve(:,i) - ve0(i) )/(ve1(i)-ve0(i)) ;
end

%% Compute Shapley values
% -------------------------------------------------------------------------

% restrict coalitions to effects of interest
switch options.coalitions
    case 'interactions'
        spl = [2^nw*ones(nu+1,1); 2^nu];
        k = mat2cell(k,spl,nw+nu);
        for i=1:nu+1
            k{i} = k{i}(:,1:nw);
        end
        k{nu+2} = k{nu+2}(:,nw+(1:nu));
        ve = mat2cell(ve,spl,nResps);
    case 'parameters'
        k = {k(:,1:nw)};
        ve = {ve};
    case 'inputs'
        k = {k(:,nw+(1:nu))};
        ve = {ve};
end

% compute first order scores
for ii = 1:numel(k)
    n = size(k{ii},2);
    nn = factorial(n);
%     v{ii} = nan(n,nResps);
    for m=1:n % loop over players
        % Shapley coeficients
        i = k{ii}(:,m);
        z = sum(k{ii},2);
        coef = (2*i-1).*factorial(z-i).*factorial(n-z-(1-i))/nn;
        % compute shapley values per se
        v{ii}(m,:) =  coef'*ve{ii};
    end
end

% compute interactions if necessary
sv = v{1};
if strcmp(options.coalitions,'interactions')
    for iu=1:nu
        svi{iu} = (v{1}-v{iu+1})./v{1};
    end
    % shapley value of inputs
    svu = v{nu+2};
else
    svi={};
end

%% Store results
% -------------------------------------------------------------------------
switch options.coalitions
    case 'parameters'
        results.parameters = sv;
    case 'inputs'
        results.inputs = sv;
    case 'interactions'
        results.parameters = sv;
        results.inputs = svu;
        results.interactions = svi;
end

end
%% Subfunctions
% =========================================================================

% -------------------------------------------------------------------------
% Compute explained variance of the model given by posterior and out, with
% the inputs and paramters pertrubed according to u_swicth and w_switch
% respectively (switch = 1 -> normal, switch=0 -> perturbation)
% -------------------------------------------------------------------------
function v=explainedVar(posterior,out,options,u_switch,w_switch)

    % prevent unecessary bells and whistles
    out.options.verbose = 0;
    out.options.DisplayWin = 0;
    out.options.inF{1}.fast = true;

    % prepare degraded model
    % .....................................................................

    % == pertub inputs
    % index of inputs to perturb
    inputIdx = options.inputIdx(u_switch==0);
    % apply perturbation
    switch options.inputPerturbation
        case 'zero'
            out.u(inputIdx,:) = 0;
        case 'average'
            out.u(inputIdx,:) = mean(out.u(inputIdx,:),2);
        case 'average_nonzero'
            for iu=inputIdx
                idxNZ = find(out.u(iu,:)~=0);
                out.u(iu,idxNZ) = mean(out.u(iu,idxNZ));
            end
    end

    % == perturb parameters
    paramIdx = options.paramIdx(w_switch==0);
    switch options.paramType
        case 'phi'
            posterior.muPhi(paramIdx) = 0*posterior.muPhi(paramIdx);
        case 'theta'
            posterior.muTheta(paramIdx) = 0*posterior.muTheta(paramIdx);
    end

    % predict data
    % .....................................................................
    [yp,~,~,~,er] = VBA_simulate (...
        out.options.dim.n_t,...
        out.options.f_fname,...
        out.options.g_fname,...
        posterior.muTheta,...
        posterior.muPhi,...
        out.u,...
        Inf,...
        Inf,...
        out.options,...
        posterior.muX0);
    g = yp-er;
    y = out.y;

    % if vertical data, transpose everything
    if options.obsIdx == 0
        g = g';
        y = y';
        options.obsIdx = 1;
        out.options.isYout = out.options.isYout';
    end

    % compute explained variance
    % .....................................................................
    v = nan(1,numel(options.obsIdx));
    for i=1:numel(options.obsIdx) % for each observation of interest
        obsIdx = options.obsIdx(i);
        in_idx = find(out.options.isYout(obsIdx,:) == 0);
        v(i) = 1-((var(y(obsIdx,in_idx)-g(obsIdx,in_idx))/var(y(obsIdx,in_idx))))  ;
    end

end

author: MBB-team

repo: VBA-toolbox


text:
function y = besseliLn(nu,x)
% Compute logarithm of besseli function (modified Bessel function of first kind).
% Written by Mo Chen (mochen80@gmail.com).
% TODO: improve precision using the method in
% Clustering on the Unit Hypersphere using von Mises-Fisher Distributions.  A. Banerjee, I. S. Dhillon, J. Ghosh, and S. Sra
[v,ierr] = besseli(nu,x);
if any(ierr ~= 0) || any(v == Inf)
    error('ERROR: logbesseli');
end
y = log(v);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
%--- help for abstvar/posterior_mode ---
%
%  Return the posterior mode from the BVAR estimation
%
%  ::
%
%     dd = posterior_mode(self)
%
%  Args:
%     self (var object): var object
%
%  Returns:
%     :
%
%     - **dd** (struct): posterior mode of the parameters
%
%

author: jmaih


text:
%% This file is the PDE simulation file of Belousov-Zhabotinsky reaction.
% Coded By: K
% Last Updated: 2019/06/24
%%
function [rhs,x_t,z_t,s_t,u_t,x,z,s,u,x_x,z_x,s_x,u_x,x_y,z_y,s_y,u_y,x_xx,z_xx,s_xx,u_xx,x_yy,z_yy,s_yy,u_yy,x_lap,z_lap,s_lap,u_lap]=...
    BZ_Reaction_PDE(t,xzsut,Kx,Kxx,Ky,Kyy,K22,n,N,Dx,Dz,Ds,Du,q,f,ksi,alpha,beta,gama,ksi2,ksi3,phi,NeedDev)

% Calculate u and v terms
xt=reshape((xzsut(1:N)),n,n);
zt=reshape((xzsut((N+1):(2*N))),n,n);
st=reshape((xzsut(2*N+1:3*N)),n,n);
ut=reshape((xzsut((3*N+1):(4*N))),n,n);

x=real(ifft2(xt));
z=real(ifft2(zt));
s=real(ifft2(st));
u=real(ifft2(ut));

% Reaction Terms
xtrhs=reshape((fft2(  (1/ksi)*( f.*z.*(q-x)./(q+x) + x - x.^2 - beta*x + s )  )),N,1);
ztrhs=reshape((fft2(  x - z - alpha*z + gama*u  )),N,1);
strhs=reshape((fft2(  (1/ksi2)*(beta*x - s + phi*u )  )),N,1);
utrhs=reshape((fft2(  (1/ksi3)*(alpha*z - gama*u )  )),N,1);

rhs=[-(Dx/Dx)*K22.*xzsut(1:N)+xtrhs
     -(Dz/Du)*K22.*xzsut(N+1:2*N)+ztrhs
     -(Ds/Du)*K22.*xzsut(2*N+1:3*N)+strhs
     -(Du/Du)*K22.*xzsut(3*N+1:4*N)+utrhs
     ];

 % If you don't need to extract the value, don't run the following code to
 % speed up
if NeedDev==1
    % Get the derivative you want
    x_x=real(ifft2(reshape((1j*Kx).*xzsut(1:N),n,n)));
    z_x=real(ifft2(reshape((1j*Kx).*xzsut(N+1:2*N),n,n)));
    s_x=real(ifft2(reshape((1j*Kx).*xzsut(2*N+1:3*N),n,n)));
    u_x=real(ifft2(reshape((1j*Kx).*xzsut(3*N+1:4*N),n,n)));
    %
    x_y=real(ifft2(reshape((1j*Ky).*xzsut(1:N),n,n)));
    z_y=real(ifft2(reshape((1j*Ky).*xzsut(N+1:2*N),n,n)));
    s_y=real(ifft2(reshape((1j*Ky).*xzsut(2*N+1:3*N),n,n)));
    u_y=real(ifft2(reshape((1j*Ky).*xzsut(3*N+1:4*N),n,n)));
    %
    x_xx=real(ifft2(reshape((1j*Kxx).^2.*xzsut(1:N),n,n)));
    z_xx=real(ifft2(reshape((1j*Kxx).^2.*xzsut(N+1:2*N),n,n)));
    s_xx=real(ifft2(reshape((1j*Kxx).^2.*xzsut(2*N+1:3*N),n,n)));
    u_xx=real(ifft2(reshape((1j*Kxx).^2.*xzsut(3*N+1:4*N),n,n)));
    %
    x_yy=real(ifft2(reshape((1j*Kyy).^2.*xzsut(1:N),n,n)));
    z_yy=real(ifft2(reshape((1j*Kyy).^2.*xzsut(N+1:2*N),n,n)));
    s_yy=real(ifft2(reshape((1j*Kyy).^2.*xzsut(2*N+1:3*N),n,n)));
    u_yy=real(ifft2(reshape((1j*Kyy).^2.*xzsut(3*N+1:4*N),n,n)));
    %
    x_t=real(ifft2(reshape(rhs(1:N),n,n)));
    z_t=real(ifft2(reshape(rhs(N+1:2*N),n,n)));
    s_t=real(ifft2(reshape(rhs(2*N+1:3*N),n,n)));
    u_t=real(ifft2(reshape(rhs(3*N+1:4*N),n,n)));
    %
    x_lap=real(ifft2(reshape(-K22.*xzsut(1:N),n,n)));
    z_lap=real(ifft2(reshape(-K22.*xzsut(N+1:2*N),n,n)));
    s_lap=real(ifft2(reshape(-K22.*xzsut(2*N+1:3*N),n,n)));
    u_lap=real(ifft2(reshape(-K22.*xzsut(3*N+1:4*N),n,n)));

end

author: dynamicslab

repo: SINDy-PI


text:
function [ algam, sgngam ] = r8_lgams ( x )

%*****************************************************************************80
%
%% R8_LGAMS evaluates the log of |gamma(x)| and sign, for an R8 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real ALGAM, the logarithm of the absolute value of
%    gamma ( X ).
%
%    Output, real SGNGAM, the sign (+1 or -1 ) of gamma ( X ).
%
  algam = r8_lngam ( x );
  sgngam = 1.0;

  if ( x <= 0.0 )

    k = floor ( mod ( - r8_aint ( x ), 2.0 ) + 0.1 );

    if ( k == 0 )
      sgngam = - 1.0;
    end

  end

  return
end


function [ux,u,m] = applypca(X)
% [ux,u,m] = applypca(X)
%
% finds principal components of
%
% input:
%  x  input data (each column is an input vector)
%
% output:
%  ux  data projected onto pca (1st row=leading principal component)
%  u   the orthogonal matrix with all principal components
%  m   mean of the data
%
% You can recover the original data with
%
% xoriginal=u*ux+repmat(m,1,size(ux,2));
%
% See also pca
%
% copyright by Kilian Q. Weinberger, 2006

[u,v]=pca(X);
m=mean(X,2);
ux=u'*(X-repmat(m,1,size(X,2)));


text:
% -----------------------------------------------------------
% file: compile_mex.m
% -----------------------------------------------------------
% authors: Sebastien Bougleux and Évariste Daller
% institution: Normandie Univ, CNRS - ENSICAEN - UNICAEN, GREYC UMR 6072
% -----------------------------------------------------------
% Execute this file in matlab to compile matlab functions for LSAP and LSAPE
% -----------------------------------------------------------
% This file is part of LSAPE.
% LSAPE is free software: you can redistribute it and/or modify
% it under the terms of the CeCILL-C License. See README for more
% details.
% -----------------------------------------------------------

disp('Compiling tools for LSAP and LSAPE ...');
files =  { 'hungarianLSAP.cpp', ...
       	   'hungarianLSAPE.cpp', ...
           'greedyLSAP.cpp'
};

str = 'mex -O -I../include ';
flags = ' ';%' -std=c++11';

for i=1:length(files)
    eval([str files{i} flags ' ']);
end

author: optimaltransport


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function tfrppagt
%TFRPPAGT Unit test for the time frequency representation TFRPPAGE.

%       O. Lemoine - March 1996.

% We test each property of the corresponding TFR :

N=128;

% Covariance by translation in time
t1=55; t2=70; f=0.3;
sig1=amgauss(N,t1).*fmconst(N,f,t1);
sig2=amgauss(N,t2).*fmconst(N,f,t2);
tfr1=tfrppage(sig1);
tfr2=tfrppage(sig2);
[tr,tc]=size(tfr1);
nu=round(f*(tc-1)*2)+1;
tfr=tfr1-tfr2(:,modulo((1:tc)-t1+t2,tc));
if any(any(abs(tfr)>sqrt(eps))),
 error('tfrppage test 1 failed');
end

% Reality of the TFR
sig=noisecg(N);
tfr=tfrppage(sig);
if sum(any(abs(imag(tfr))>sqrt(eps)))~=0,
 error('tfrppage test 2 failed');
end

% Energy conservation
sig=noisecg(N);
tfr=tfrppage(sig);
Es=norm(sig)^2;
Etfr=sum(mean(tfr));
if abs(Es-Etfr)>sqrt(eps),
 error('tfrppage test 3 failed');
end

% Time-marginal
sig=noisecg(N);
tfr=tfrppage(sig);
ip1=abs(sig).^2;
ip2=mean(tfr)';
if any(abs(ip1-ip2)>sqrt(eps)),
 error('tfrppage test 4 failed');
end

% Conservation of the time support (wide-sense)
sig=[zeros(N/4,1);noisecg(N/2);zeros(N/4,1)];
tfr=tfrppage(sig);
if sum(any(abs(tfr(:,1:N/4-1))>sqrt(eps))) | ...
   sum(any(abs(tfr(:,(3*N/4+1):N))>sqrt(eps))),
 error('tfrppage test 5 failed');
end

% time localization
t0=30; sig=((1:N)'==t0);
tfr=tfrppage(sig);
[ik,jk]=find(tfr~=0.0);
if any(jk~=t0)|any(ik'-(1:N)),
 error('tfrppage test 6 failed');
end;

% A PPAGED with a constant window is a PAGED
sig=noisecg(N);
tfr1=tfrppage(sig,1:N,N,ones(N+1,1));
tfr2=tfrpage(sig);
if any(any(abs(tfr1-tfr2)>sqrt(eps))),
 error('tfrppage test 7 failed');
end

N=129;

% Covariance by translation in time
t1=55; t2=70; f=0.3;
sig1=amgauss(N,t1).*fmconst(N,f,t1);
sig2=amgauss(N,t2).*fmconst(N,f,t2);
tfr1=tfrppage(sig1);
tfr2=tfrppage(sig2);
[tr,tc]=size(tfr1);
nu=round(f*(tc-1)*2)+1;
tfr=tfr1-tfr2(:,modulo((1:tc)-t1+t2,tc));
if any(any(abs(tfr)>sqrt(eps))),
 error('tfrppage test 8 failed');
end

% Reality of the TFR
sig=noisecg(N);
tfr=tfrppage(sig);
if sum(any(abs(imag(tfr))>sqrt(eps)))~=0,
 error('tfrppage test 9 failed');
end

% Energy conservation
sig=noisecg(N);
tfr=tfrppage(sig);
Es=norm(sig)^2;
Etfr=sum(mean(tfr));
if abs(Es-Etfr)>sqrt(eps),
 error('tfrppage test 10 failed');
end

% Time-marginal
sig=noisecg(N);
tfr=tfrppage(sig);
ip1=abs(sig).^2;
ip2=mean(tfr)';
if any(abs(ip1-ip2)>sqrt(eps)),
 error('tfrppage test 11 failed');
end

% Conservation of the time support (wide-sense)
sig=[zeros(round(N/4),1);noisecg(round(N/2));zeros(round(N/4),1)];
tfr=tfrppage(sig);
if sum(any(abs(tfr(:,1:round(N/4)-1))>sqrt(eps))) | ...
   sum(any(abs(tfr(:,(round(3*N/4)+1):N))>sqrt(eps))),
 error('tfrppage test 12 failed');
end

% time localization
t0=30; sig=((1:N)'==t0);
tfr=tfrppage(sig);
[ik,jk]=find(tfr~=0.0);
if any(jk~=t0)|any(ik'-(1:N)),
 error('tfrppage test 13 failed');
end;

% A PPAGED with a constant window is a PAGED
sig=noisecg(N);
tfr1=tfrppage(sig,1:N,N,ones(2*N+1,1));
tfr2=tfrpage(sig);
if any(any(abs(tfr1-tfr2)>sqrt(eps))),
 error('tfrppage test 14 failed');
end

author: HeLiangHIT


text:
function test_failed=test_signals;
%TEST_SIGNALS
%
%  The script ensures that all files are correctly included by loading
%  the singals and checking their sizes.

disp(' ===============  TEST_SIGNALS ===========');

test_failed=0;

[test_failed,fail]=ltfatdiditfail(numel(bat)-400,test_failed);
disp(['SIGNALS BAT       ',fail]);

[test_failed,fail]=ltfatdiditfail(numel(batmask)-1600,test_failed);
disp(['SIGNALS BATMASK   ',fail]);

[test_failed,fail]=ltfatdiditfail(numel(greasy)-5880,test_failed);
disp(['SIGNALS GREASY    ',fail]);

[test_failed,fail]=ltfatdiditfail(numel(linus)-41461,test_failed);
disp(['SIGNALS LINUS     ',fail]);

[test_failed,fail]=ltfatdiditfail(numel(gspi)-262144,test_failed);
disp(['SIGNALS GSPI      ',fail]);

[test_failed,fail]=ltfatdiditfail(numel(traindoppler)-157058,test_failed);
disp(['SIGNALS TRAINDOPPLER',fail]);

[test_failed,fail]=ltfatdiditfail(numel(otoclick)-2210,test_failed);
disp(['SIGNALS OTOCLICK  ',fail]);

[test_failed,fail]=ltfatdiditfail(numel(cameraman)-65536,test_failed);
disp(['SIGNALS CAMERAMAN ',fail]);

[test_failed,fail]=ltfatdiditfail(numel(cocktailparty)-363200,test_failed);
disp(['SIGNALS COCKTAILPARTY ',fail]);

[test_failed,fail]=ltfatdiditfail(numel(lichtenstein)-262144*3,test_failed);
disp(['SIGNALS LICHTENSTEIN ',fail]);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


%  Script file: microphone.m
%
%  Purpose:
%    This program plots the gain pattern of a cardioid
%    microphone.
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    01/15/07    S. J. Chapman        Original code
%
% Define variables:
%   g         -- Microphone gain constant
%   gain      -- Gain as a function of angle
%   theta     -- Angle from microphone axis (radians)

% Calculate gain versus angle
g = 0.5;
theta = 0:pi/20:2*pi;
gain = 2*g*(1+cos(theta));

% Plot gain
polar (theta,gain,'r-');
title ('\bfGain versus angle \theta');


text:
function varargout = rts_DynConsole(varargin)
% RTS_DYNCONSOLE MATLAB code for rts_DynConsole.fig
%      RTS_DYNCONSOLE, by itself, creates a new RTS_DYNCONSOLE or raises the existing
%      singleton*.
%
%      H = RTS_DYNCONSOLE returns the handle to a new RTS_DYNCONSOLE or the handle to
%      the existing singleton*.
%
%      RTS_DYNCONSOLE('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in RTS_DYNCONSOLE.M with the given input arguments.
%
%      RTS_DYNCONSOLE('Property','Value',...) creates a new RTS_DYNCONSOLE or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before rts_DynConsole_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to rts_DynConsole_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help rts_DynConsole

% Last Modified by GUIDE v2.5 13-Dec-2013 17:48:09

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @rts_DynConsole_OpeningFcn, ...
                   'gui_OutputFcn',  @rts_DynConsole_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT
end

% --- Executes just before rts_DynConsole is made visible.
function rts_DynConsole_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to rts_DynConsole (see VARARGIN)

% Choose default command line output for rts_DynConsole
    handles.output = hObject;

    hRtsMainGUI = varargin{1};
    vesselGUID = varargin{2};
    handles.mainGUI = hRtsMainGUI;
    handles.vesselGUID = vesselGUID;

    rts_initRTSDynPanel(hRtsMainGUI, hObject, handles.clockLabel, handles.connStatusLabel, vesselGUID,...
                        handles.attPointingTextDispLabel, handles.attRatesTextDispLabel, handles.attCmdTextDispLabel, handles.attDispAxes,...
                        handles.actStatusTextDispLabel, handles.rwaInfoUitable, handles.attDerivedTextDispLabel);

    % Update handles structure
    guidata(hObject, handles);

    % UIWAIT makes rts_DynConsole wait for user response (see UIRESUME)
    % uiwait(handles.rts_DynConsole);
end

% --- Outputs from this function are returned to the command line.
function varargout = rts_DynConsole_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;
end

% --- Executes when user attempts to close rts_DynConsole.
function rts_DynConsole_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to rts_DynConsole (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: delete(hObject) closes the figure
    rts_pageCloseStopTimers(hObject);
    delete(hObject);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function fx = p55_fun ( n, x )

%*****************************************************************************80
%
%% P55_FUN evaluates the integrand for problem 55.
%
%  Interval:
%
%    a = 0 <= x <= 1 = b
%
%  Integrand:
%
%    exp ( - c^2 * ( x - x0 )^2 )
%
%  Exact Integral:
%
%    sqrt ( pi )
%    * ( erf ( c * ( b - x0 ) ) - erf ( c * ( a - x0 ) ) )
%    / ( 2 * c )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  c = p55_param ( 'get', 'c', [] );
  x0 = p55_param ( 'get', 'x0', [] );
  x0 = x0 * ones ( size ( x ) );

  fx(1:n) = exp ( - c * c * ( x(1:n) - x0 ).^2 );

  return
end


text:
function item = initialise(item, val, dflag)

% function item = initialise(item, val, dflag)
% This is a generic initialisation function to insert values into the val
% field of a cfg_item. It is usable for all leaf cfg_item classes
% (cfg_entry, cfg_files, cfg_menu). Assignment checks are done through
% subsasgn.
% dflag is ignored for leaf entries - the passed value is assigned
% unconditionally.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: initialise.m 4898 2012-09-05 13:40:16Z volkmar $

rev = '$Rev: 4898 $'; %#ok

d1flag = false;
if ischar(val) && strcmp(val, '<DEFAULTS>')
    if isempty(item.def)
        % don't change anything if no .def field
        return;
    else
        try
            val = feval(item.def, {});
            d1flag = true;
        catch
            cfg_message('matlabbatch:initialise:defaults', ...
                    '%s: No matching defaults value found.', ...
                    subsasgn_checkstr(item,substruct('.','val')));
            val = '<UNDEFINED>';
        end;
    end;
end;

subs = substruct('.', 'val', '{}', {1});
if ischar(val) && strcmp(val, '<UNDEFINED>') % val may be <UNDEFINED>
    item.val = {};
else
    % Modify job before initialisation
    if ~d1flag && ~dflag && ~isempty(item.rewrite_job)
        val = feval(item.rewrite_job, val);
    end
    item = subsasgn(item, subs, val);
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
path(path, '../toolbox_graph_data/off/');

name = 'mushroom';
name = 'nefertiti';
clear options;
options.name = name;
[vertex,face] = read_off([name '.off']);

rep = ['results/mesh-duality/' name '/'];
if not(exist(rep))
    mkdir(rep);
end

clf;
plot_mesh(vertex,face,options);
shading faceted;
camlight; axis tight; axis equal;
saveas(gcf, [rep name '-mesh.png'], 'png');

clf;
A = triangulation2adjacency(face);
plot_graph(A,vertex);
view(2); axis tight; axis equal;
saveas(gcf, [rep name '-graph.png'], 'png');

[A1,vertex1] = compute_dual_graph(face,vertex);
clf;
plot_graph(A1,vertex1);
view(2); axis tight; axis equal;
saveas(gcf, [rep name '-dual.png'], 'png');

% 1:4 subdivision
options.sub_type = '1:4';
[vertex2,face2] = perform_mesh_subdivision(vertex',face',1, options);
[vertex3,face3] = perform_mesh_subdivision(vertex2,face2,1, options);

clf;
plot_mesh(vertex2,face2,options);
shading faceted;
camlight; axis tight; axis equal;
saveas(gcf, [rep name '-subdivide-4-once.png'], 'png');
clf;
plot_mesh(vertex3,face3,options);
shading faceted;
camlight; axis tight; axis equal;
saveas(gcf, [rep name '-subdivide-4-twice.png'], 'png');

% 1:3 subdivision
options.sub_type = '1:3';
[vertex2,face2] = perform_mesh_subdivision(vertex,face,1, options);
[vertex3,face3] = perform_mesh_subdivision(vertex2,face2,1, options);

clf;
plot_mesh(vertex2,face2,options);
shading faceted;
camlight; axis tight; axis equal;
saveas(gcf, [rep name '-subdivide-3-once.png'], 'png');
clf;
plot_mesh(vertex3,face3,options);
shading faceted;
camlight; axis tight; axis equal;
saveas(gcf, [rep name '-subdivide-3-twice.png'], 'png');

author: gpeyre

repo: matlab-toolboxes


text:
function namesPhysContrasts = tapas_physio_get_contrast_names_default()
% Returns cell of defaults names for phys contrasts used e.g. in reporting
%
%   namesPhysContrasts = tapas_physio_get_contrast_names_default()
%
% IN
%
% OUT
%
% EXAMPLE
%   tapas_physio_get_contrast_names_default
%
%   See also tapas_physio_report_contrasts tapas_physio_compute_tsnr_gains

% Author: Lars Kasper
% Created: 2016-10-03
% Copyright (C) 2016 TNU, Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public
% License (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

namesPhysContrasts = {
    'All Phys'
    'Cardiac'
    'Respiratory'
    'Card X Resp Interation'
    'HeartRateVariability'
    'RespiratoryVolumePerTime'
    'Noise Rois'
    'Movement'
    'All Phys + Move'
    };

author: translationalneuromodeling

repo: tapas


text:
% Downsampling procedure.
%
% Arguments:
%   grayscale I image
%   downsampling filter 'filter', should be a 1D separable filter.
%   'border_mode' should be 'circular', 'symmetric', or 'replicate'. See 'imfilter'.
%
% If image width W is odd, then the resulting image will have width (W-1)/2+1,
% Same for height.
%
% tom.mertens@gmail.com, August 2007
%

function R = downsample_(I, filter)

border_mode = 'symmetric';

% low pass, convolve with separable filter
R = imfilter(I,filter,border_mode);     %horizontal
R = imfilter(R,filter',border_mode);    %vertical

% decimate
r = size(I,1);
c = size(I,2);
R = R(1:2:r, 1:2:c, :);

author: mahmoudnafifi


text:
%% patchCurvature
% Below is a demonstration of the features of the |patchCurvature| function

%% Syntax
% |[Vd,Fd,Fds]=patchCurvature(V,F);|

%% Description
% Computes curvature metrics for the patch data defined by the faces F and
% the vertices V.

%% Examples

%%
clear; close all; clc;

%%
% Plot settings
cMap=warmcold(250);

%%

[F,V]=graphicsModels(9);
% [F,V]=stanford_bunny;
% [F,V]=tri2quad(F,V);

%% Compute curvature

[U_min,U_max,C_min,C_max,C_mean,C_gauss] = patchCurvature(F,V);

%% Visualize curvature on mesh

% Compute plot variables
C_min_V=faceToVertexMeasure(F,V,C_min); %Vertex data for interpolated shading
C_max_V=faceToVertexMeasure(F,V,C_max); %Vertex data for interpolated shading
VN=patchCentre(F,V); %Element centres used for vector origins
vecPlotSize=mean(patchEdgeLengths(F,V)); %Vector plotting size

% Visualize
cFigure;
subplot(1,2,1); hold on;
title('C_{min}');
hp=gpatch(F,V,C_min_V,'none',0.9);
hp.FaceColor='interp';
colormap(gca,cMap); colorbar;
quiverVec(VN,U_min,vecPlotSize,'k');
axisGeom;
c=max(abs(C_min(:)));
caxis(0.25*[-c c]);
camlight headlight;

subplot(1,2,2); hold on;
title('C_{max}');
hp=gpatch(F,V,C_max_V,'none',0.9);
hp.FaceColor='interp';
quiverVec(VN,U_max,vecPlotSize,'k');
colormap(gca,cMap); colorbar;
axisGeom;
c=max(abs(C_max(:)));
caxis(0.25*[-c c]);
camlight headlight;

drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [ X Y ] = selPts
%[ X Y ] = selPts
%	manually select several points in current figure
%   selection is stopped when ESC is pressed

X = []; Y = [];
ctrlLen = 0;
while 1
	[X1 Y1] = ginput(1);
	if isempty(X1), break; end
	X = [X;X1];
	Y = [Y;Y1];
	plot(X1,Y1,'bo');
	ctrlLen = ctrlLen+1;
% 	text(X1+3,Y1,num2str(ctrlLen));
end

end


text:
function A = multiheadAttention(Q, K, V,nvp)
% multiheadAttention   Multi-head Attention
%
%   A = multiheadAttention(Q, K, V) computes scaled dot product attention
%   for multiple attention heads as outlined in [1] (see Section 3.2.1 and
%   Figure 2). Note that this function computes the attention for multiple
%   attention heads at once for efficiency. Q is a collection of query
%   matrices, K is a collection of key matrices and V is a collection of
%   value matrices. The output A is a collection of attention matrices. See
%   below for details.
%
%   Inputs:
%       Q   - numFeatures-by-numInputSubWords-by-numHeads-by-numObs array of queries.
%       K   - numFeatures-by-numAllSubWords-by-numHeads-by-numObs array of keys.
%       V   - numFeatures-by-numAllSubWords-by-numHeads-by-numObs array of values.
%
%   Outputs:
%       A   - numFeatures-by-numInputSubWords-by-numHeads-by-numObs array of attention matrices.
%
%   A = multiheadAttention(Q, K, V, 'PARAM1', VAL1, 'PARAM2', VAL2, ...)
%   specifies the optional parameter name/value pairs:
%
%     'CausalMask'  - A scalar logical to turn causal masking on or off. Causal
%                     masking prevents tokens at time T attending to tokens
%                     at time S<T. The default is true.
%
%     'Dropout'     - The dropout probability for the attention
%                     probabilities. The default is 0.
%
%     'InputMask'   - A logical mask to mask attending to particular
%                     tokens, for example padding tokens. The default is
%                     [], interpreted as not applying any masking.
%
%   References:
%
%   [1] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion
%       Jones, Aidan N. Gomez, Lukasz Kaiser, Illia Polosukhin, "Attention
%       Is All You Need", https://arxiv.org/abs/1706.03762
arguments
    Q
    K
    V
    nvp.CausalMask (1,1) logical = true
    nvp.Dropout (1,1) double {mustBeNonnegative,mustBeLessThanOrEqual(nvp.Dropout,1)} = 0
    nvp.InputMask = []
end

% We compute attention weights by taking the product between Q and K
% matrices. W is numAllSubWords-by-numInputSubWords-by-numHeads. Each
% element of W is the dot product of a query vector from Q and a key vector
% from K.
W = dlmtimes(permute(K, [2 1 3 4]), Q);

% Divide by square root of d
W = W./sqrt(size(Q,1));

% Apply masking
W = transformer.layer.maskAttentionWeights(W,'CausalMask',nvp.CausalMask,'InputMask',nvp.InputMask);

% Apply softmax
W = softmax(W, 'DataFormat', 'CTUB');

% Apply dropout
W = transformer.layer.dropout(W,nvp.Dropout);

% We compute the attention by taking products between the attention weights
% W and V. A is numFeatures-by-numInputSubWords-by-numHeads. One
% interpretation of A is that it is the expected value of V according to
% the probability distribution given by W.
A = dlmtimes(V, W);
end

author: matlab-deep-learning

repo: transformer-models


text:
function y = str2double(x)
% Linearizes char input to a row, like Matlab does
if ischar(x)
    x = x(:).';
end
y = builtin('str2double', x);
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
function [ meta, symmetry ] = meta2sicd_tiffrs2( tiff_filename )
%META2SICD_TIFFRS2 Compile SICD metadata for RS2 data package
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

[ meta, symmetry ] = meta2sicd_tiffrs( tiff_filename, 'RS2' );

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function dfdx_sort = map_to_original(dfdx,x, x_indep);

for i = 1:length(x(:))
    x_origvar(i) = getvariables(x(i));
end
x_var = getvariables(x_indep);
dfdx_sort = [];
for i = 1:length(x_origvar)
    dfdx_sort = [dfdx_sort dfdx(:,find(x_origvar(i) == x_var))];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%STATSCHECK  Check whether the STATS toolbox is in the path

function n = statscheck

  n = ~isempty(ver('stats'));
  if (nargout == 0) & (n == 0)
		error([newline 'The Matlab STATS toolbox is missing.'])
  end

return

author: marianux

repo: ecg-kit

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


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% Images from a histologial serial sectioning,
% data courtesy Oliver Schmitt, Institute of Anatomy, University of Rostock, Germany
%==============================================================================

checkDataFile
if expfileExists, return; end;
viewPara = {'viewImage','viewImage2D','colormap','gray(256)'};
expfile = jpgs2data('','HNSP-T.jpg','HNSP-R.jpg','omega',[0,2,0,1],'viewPara',viewPara,'m',[512,256]);
checkDataFile
%==============================================================================


text:
function [W, R, P] = miso_firwiener(N, X, y)
%MISO_FIRWIENER Optimal FIR Wiener filter for multiple inputs.
%   MISO_FIRWIENER(N, X, Y) computes the optimal FIR Wiener filter of
%   order N, given any number of (stationary) random input signals as
%   the columns of matrix X, and one output signal in column vector Y.

%   Author: Keenan Pepper
%   Last modified: 2007-12-21

%   References:
%     [1] Y. Huang, J. Benesty, and J. Chen, Acoustic MIMO Signal
%     Processing, Springer−Verlag, 2006, page 48

% Number of input channels.
M = size(X, 2);

% Input covariance matrix, in abbreviated block Toeplitz form.
R = zeros(M*(N+1), M);
for m = 1:M
    for i = 1:M
         rmi = xcorr(X(:,m), X(:,i), N);
         Rmi = flipud(rmi(1:N+1));
         top = (m-1) * (N+1) + 1;
         bottom = m * (N+1);
         R(top:bottom,i) = Rmi;
    end
end

R = reshape(R, [N+1,M,M]);      % This just permutes the indices to
R = permute(R, [2,1,3]);        % change R from Toeplitz-block to
R = reshape(R, [M*(N+1),M]);    % block-Toeplitz.

% Cross−correlation vector.
P = zeros(1, M*(N+1));
for i = 1:M
    top = (i-1)*(N+1)+1;
    bottom = i * (N+1);
    p = xcorr(y, X(:,i), N);
    P(top:bottom) = p(N+1:2*N+1)';
end

P = reshape(P, [N+1,M]);        % More index rearrangement.
P = reshape(P', [M*(N+1),1]);

W = block_levinson(P, R);

W = reshape(W, [M,N+1]);        % Make output compatible with earlier
W = reshape(W', [1, M*(N+1)]);  % version.

author: Sable


text:
% ALM (Tang and Nehorai 2011)
% process_video('RPCA', 'ALM', 'dataset/demo.avi', 'output/demo_ALM.avi');
[L,S] = alm(M);

author: andrewssobral

repo: lrslibrary


text:
function row = polygonToRow(polygon, varargin)
%POLYGONTOROW Convert polygon coordinates to a row vector
%
%   ROW = polygonToRow(POLY);
%   where POLY is a N-by-2 array of points representing vertices of the
%   polygon, converts the vertex coordinates into a linear array:
%   ROW = [X1 Y1 X2 Y2 .... XN YN]
%
%   ROW = polygonToRow(POLY, TYPE);
%   Can coose another format for converting polygon. Possibilities are:
%   'interlaced' (default}, as described above
%   'packed': ROW has format [X1 X2 ... XN Y1 Y2 ... YN].
%
%   Example
%   square = [10 10 ; 20 10 ; 20 20 ; 10 20];
%   row = polygonToRow(square)
%   row =
%       10   10   20   10   20   20   10   20
%
%   % the same with different ordering
%   row = polygonToRow(square, 'packed')
%   row =
%       10   20   20   10   10   10   20   20
%
%
%   See also
%   polygons2d, rowToPolygon
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-07-23,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% determines ordering type
type = 'interlaced';
if ~isempty(varargin)
    type = varargin{1};
end

if strcmp(type, 'interlaced')
    % ordering is [X1 Y1 X2 X2... XN YN]
    Np = size(polygon, 1);
    row = reshape(polygon', [1 2*Np]);

elseif strcmp(type, 'packed')
    % ordering is [X1 X2 X3... XN Y1 Y2 Y3... YN]
    row = polygon(:)';
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function varargout = process_cohere1n( varargin )
% PROCESS_COHERE1N: Compute the coherence between all the pairs of signals, in one file.
%
% USAGE:   OutputFiles = process_cohere1n('Run', sProcess, sInputA)
%                        process_cohere1n('Test')

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Coherence NxN [Deprecated]';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Connectivity';
    sProcess.Index       = 656;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Connectivity';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data',     'results',  'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;

    % === CONNECT INPUT
    sProcess = process_corr1n('DefineConnectOptions', sProcess, 1);
    % === REMOVE EVOKED REPONSE
    sProcess.options.removeevoked.Comment = 'Remove evoked response from each trial';
    sProcess.options.removeevoked.Type    = 'checkbox';
    sProcess.options.removeevoked.Value   = 0;
    sProcess.options.removeevoked.Group   = 'input';
    % === COHERENCE METHOD
    sProcess.options.cohmeasure.Comment = {...
        ['Magnitude-squared Coherence
' ...
        '|C|^2 = |Gxy|^2/(Gxx*Gyy)'], ...
        ['Imaginary Coherence (2019)
' ...
        'IC    = |imag(C)|'], ...
        ['Lagged Coherence (2019)
' ...
        'LC    = |imag(C)|/sqrt(1-real(C)^2)'], ...
        [' Imaginary Coherence (before 2019)
' ...
        ' IC    = imag(C)^2 / (1-real(C)^2) ']; ...
        'mscohere', 'icohere2019', 'lcohere2019', 'icohere'};
    sProcess.options.cohmeasure.Type    = 'radio_label';
    sProcess.options.cohmeasure.Value   = 'mscohere';
    % === Overlap
    sProcess.options.overlap.Comment = 'Overlap for PSD estimation:';
    sProcess.options.overlap.Type    = 'value';
    sProcess.options.overlap.Value   = {50, '%', []};
    % === MAX FREQUENCY RESOLUTION
    sProcess.options.maxfreqres.Comment = 'Maximum frequency resolution:';
    sProcess.options.maxfreqres.Type    = 'value';
    sProcess.options.maxfreqres.Value   = {2,'Hz',2};
    % === HIGHEST FREQUENCY OF INTEREST
    sProcess.options.maxfreq.Comment = 'Highest frequency of interest:';
    sProcess.options.maxfreq.Type    = 'value';
    sProcess.options.maxfreq.Value   = {60,'Hz',2};
    % === OUTPUT MODE
    sProcess.options.outputmode.Comment = {'Save individual results (one file per input file)', 'Concatenate input files before processing (one file)', 'Save average connectivity matrix (one file)'};
    sProcess.options.outputmode.Type    = 'radio';
    sProcess.options.outputmode.Value   = 1;
    sProcess.options.outputmode.Group   = 'output';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputA) %#ok<DEFNU>
    % Input options
    OPTIONS = process_corr1n('GetConnectOptions', sProcess, sInputA);
    if isempty(OPTIONS)
        OutputFiles = {};
        return
    end

    % Metric options
    OPTIONS.Method = 'cohere';
    OPTIONS.RemoveEvoked  = sProcess.options.removeevoked.Value;
    OPTIONS.MaxFreqRes    = sProcess.options.maxfreqres.Value{1};
    OPTIONS.MaxFreq       = sProcess.options.maxfreq.Value{1};
    OPTIONS.CohOverlap    = 0.50;  % First pre-define the overlap
    OPTIONS.pThresh       = 0.05;
    OPTIONS.CohMeasure    = sProcess.options.cohmeasure.Value;

    % Change the overlap if it is specified
    if isfield(sProcess.options, 'overlap') && isfield(sProcess.options.overlap, 'Value') && ...
       iscell(sProcess.options.overlap.Value) && ~isempty(sProcess.options.overlap.Value) && ~isempty(sProcess.options.overlap.Value{1})
        OPTIONS.CohOverlap = sProcess.options.overlap.Value{1}/100 ;
    end

    % Compute metric
    OutputFiles = bst_connectivity({sInputA.FileName}, [], OPTIONS);
end

%% ===== TEST FUNCTION =====
function Test() %#ok<DEFNU>
    % Start a new report
    bst_report('Start');
    % Get test datasets
    sFile = process_simulate_ar('Test');
    % Loop on frequency resolutions
    for freq = [1 2 3 5 10 20]
        % Coherence process
        sTmp = bst_process('CallProcess', 'process_cohere1n', sFile, [], ...
            'timewindow',   [], ...    % All the time in input
            'cohmeasure',   1, ...     % 1=Magnitude-squared, 2=Imaginary
            'overlap',      3, ...     % 50%
            'maxfreqres',   freq, ...  % VARIES
            'maxfreq',      [], ...    % No maximum frequency
            'pThresh',      0.05, ...  % p-value thrshold
            ... 'isfreqbands',  0, ...     % No frequency bands
            ... 'freqbands',    [], ...
            'outputmode',   1);        % Save individual results (one file per input file)
        % Snapshot: spectrum
        bst_process('CallProcess', 'process_snapshot', sTmp, [], ...
            'target',       11, ...  % Connectivity matrix (image)
            'modality',     1, 'orient', 1, 'time', 0, 'contact_time', [-40, 110], 'contact_nimage', 16, ...
            'Comment',      [sTmp.Comment, ': ' sTmp.Comment]);
    end
    % Save and display report
    ReportFile = bst_report('Save', sTmp);
    bst_report('Open', ReportFile);
end

author: brainstorm-tools


text:
%% ongrid
% Below is a demonstration of the features of the |ongrid| function

%%
clear; close all; clc;

%% Syntax
% |[L1t,L2t,L3t]=ongrid(V1,V2,V3,cm,vm);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.


text:
function [fAdj,numAdj]=sortTriangFacesAroundVertices(vertInFaces,fAdj,numAdj)
%%SORTTRIANGFACESAROUNDVERTICES Consider a closed polyhedron that is made
%       up of triangles whose vertices are sorted in counterclockwise order
%       when looking from the outside the polyhedron at the triangle, given
%       a list of which vertices are in which triangle, obtain a list of
%       which faces (triangles) are adjacent to each vertex with the faces
%       specified in counterclockwise order.
%
%INPUTS: vertInFaces A 3XnumFaces set of indices of vertices in all of the
%                    triangles in the polyhedron (all faces are triangles).
%                    The vertex locations do not need to be provided; just
%                    the indices of the triangles.
%       fAdj, numAdj If available, fAdj is a maxAdjFacesXnumVert matrix,
%                    where for the kth vertex, fAdj(1:numAdj(k),k) is the
%                    set of indices of the faces that are adjacent to that
%                    vertex. It effectively selects a column of
%                    vertInFaces. numAdj is a length numVertices vector. If
%                    these two inputs are omitted, they are reconstructed
%                    from vertInFaces.
%
%OUTPUTS: fAdj,numAdj These have the same definition as fAdj and numAdj on
%                     the input but the faces are sorted in
%                     counterclockwise order around each vertex.
%
%Having the ordering of faces around vertices recorded in a
%counterclockwise order is important in algorithms such as in [1] for
%finding 3D barycentric coordinates.
%
%In a completely closed polyhedron, one can start at the vertex of interest
%and then follow the triangles around in the correct order. Assume that all
%of the trianges are given with indices in counterclockwise order. One can
%rotates the indices of each triangles adjacent to the vertex such that the
%vertex of interest is in the first index. One can then trace the triangles
%around in the correct order by following the indices. The second vertex of
%the first triangle will also be the last index traced. To trace around,
%consider the example below:
%     4-------3
%    / \     / \
%   /   \   /   \
%  /     \ /     \
% 5-------1-------2
%  \      |      /
%   \     |     /
%    \    |    /
%     \   |   /
%      \  |  /
%       \ | /
%        \|/
%         6
%
%Starting with face 1-2-3, the second triangle must share vertex 3. Thus,
%one must search for a triangle whose second index is 3. That leads to
%1-3-4. Similarly, the third index of that triangle will be 4 and the
%second index of the next triangle must be 4, so tht leads to triangle
%1-4-5. Thus, one can continue tracing triangles aorund the vertex until
%the common vertex is number 2, in which case one knows that the final
%triangle has been found.
%
%REFERENCES:
%[1] M. S. Floater, "Generalized barycentric coordinates and applications,"
%    Acta Numerica, vol. 24, pp. 161-214, 1 May 2015.
%
%September 2022 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numVert=max(vertInFaces(:));

if(nargin<2||isempty(fAdj))
    %Determine the maximum number of adjacent vertices.
    numAdj=zeros(numVert,1);
    for k=1:numVert
        idx=find(any(vertInFaces==k,1));
        numAdj(k)=length(idx);
    end

    %Store the adjacent vertices, but not in any particular order.
    maxAdj=max(numAdj);
    fAdj=zeros(maxAdj,numVert);
    for k=1:numVert
        idx=find(any(vertInFaces==k,1));
        fAdj(1:numAdj(k),k)=idx;
    end
end

for k=1:numVert
    numFaces=numAdj(k);
    %Rotate the vertices in each triangle until vertex k is first.
    for curFace=1:numFaces
        fIdx=fAdj(curFace,k);

        if(vertInFaces(2,fIdx)==k)
            vertInFaces(:,fIdx)=vertInFaces([2;3;1],fIdx);
        elseif(vertInFaces(3,fIdx)==k)
            vertInFaces(:,fIdx)=vertInFaces([3;1;2],fIdx);
        end%Otherwise, assume vertInFaces(1,fIdx)==k and no rotation is
           %needed.
    end
    newOrder=zeros(numFaces,1);
    newOrder(1)=1;%The first face remains unchanged.
    assignedFaces=false(numFaces,1);
    assignedFaces(1)=true;

    curAddedIdx=1;
    fIdx=fAdj(1,k);
    endNextIdx=vertInFaces(2,fIdx);
    curNextIdx=vertInFaces(3,fIdx);
    while(curNextIdx~=endNextIdx)
        for curFace=2:numFaces
            if(assignedFaces(curFace))
                continue;
            end
            fIdx=fAdj(curFace,k);

            if(vertInFaces(2,fIdx)==curNextIdx)
                %Found the next face.
                break;
            end
        end

        curAddedIdx=curAddedIdx+1;
        newOrder(curAddedIdx)=curFace;
        assignedFaces(curFace)=true;
        curNextIdx=vertInFaces(3,fIdx);
    end
    fAdj(1:numFaces,k)=fAdj(newOrder,k);
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 Wrot = whiteningLocal(CC, yc, xc, nRange)

Wrot = zeros(size(CC,1), size(CC,1));
for j = 1:size(CC,1)
    ds          = (xc - xc(j)).^2 + (yc - yc(j)).^2;
    [~, ilocal] = sort(ds, 'ascend');
    ilocal      = ilocal(1:nRange);

    [E, D]      = svd(CC(ilocal, ilocal));
    D           = diag(D);
    eps         = 1e-6;
    wrot0       = E * diag(1./(D + eps).^.5) * E';
    Wrot(ilocal, j)  = wrot0(:,1);
end

author: cortex-lab

repo: KiloSort

save_path: github-repos/MATLAB/cortex-lab-KiloSort


text:
%% Generate unit test data for rmathlib from MATLAB's statistics toolbox
% J.Currie October 2017
clc
clear

%% Student's t Distribution
t.x       = 0;
t.q       = 1.0589;
t.p       = 0.95;
t.df      = 1:10;
t.qt      = tinv(t.p,t.df);
t.dt      = tpdf(t.x,t.df);
t.pt      = tcdf(t.q,t.df);

%% F Distribution
f.x     = 3;
f.q     = 3;
f.p     = 0.95;
f.df1   = 5:15;
f.df2   = 10:20;
f.qf    = finv(f.p,f.df1,f.df2);
f.df    = fpdf(f.x,f.df1,f.df2);
f.pf    = fcdf(f.q,f.df1,f.df2);

%% Normal Distribution
norm.x      = [-2 -1 0 1 2];
norm.q      = [-2,-1,0,1,2];
norm.p      = [0.1,0.25,0.5,0.75,0.9];
norm.mu     = 0;
norm.sig    = 1;
norm.qnorm  = icdf('Normal',norm.p,norm.mu,norm.sig);
norm.dnorm  = pdf('Normal',norm.x,norm.mu,norm.sig);
norm.pnorm  = cdf('Normal',norm.q,norm.mu,norm.sig);

%% Uniform Distribution
unif.x      = [-2 -1 0 1 2];
unif.q      = [-2,-1,0,1,2];
unif.p      = [0.1,0.25,0.5,0.75,0.9];
unif.a      = 0.1;
unif.b      = 0.5;
unif.qunif  = icdf('Uniform',unif.p,unif.a,unif.b);
unif.dunif	= pdf('Uniform',unif.x,unif.a,unif.b);
unif.punif	= cdf('Uniform',unif.q,unif.a,unif.b);

%% Chi2 distribution
chi2.x      = 0;
chi2.q      = 1.0589;
chi2.p      = 0.95;
chi2.df     = 1:10;
chi2.qchisq	= chi2inv(chi2.p,chi2.df);
chi2.dchisq	= chi2pdf(chi2.x,chi2.df);
chi2.pchisq	= chi2cdf(chi2.q,chi2.df);

%% Exponential distribution
exp.x       = 0;
exp.q       = 1.0589;
exp.p       = 0.95;
exp.u       = 1:10;
exp.qexp	= icdf('Exponential',exp.p,exp.u);
exp.dexp	= pdf('Exponential',exp.x,exp.u);
exp.pexp	= cdf('Exponential',exp.q,exp.u);

%% Beta Distribution
beta.x      = [-2 -1 0 1 2];
beta.q      = [-2,-1,0,1,2];
beta.p      = [0.1,0.25,0.5,0.75,0.9];
beta.a      = 0.1;
beta.b      = 0.5;
beta.qbeta	= icdf('Beta',beta.p,beta.a,beta.b);
beta.dbeta	= pdf('Beta',beta.x,beta.a,beta.b);
beta.pbeta	= cdf('Beta',beta.q,beta.a,beta.b);

%% Cauchy Distribution
% cauchy.x    = [-2 -1 0 1 2];
% cauchy.q    = [-2,-1,0,1,2];
% cauchy.p    = [0.1,0.25,0.5,0.75,0.9];
% cauchy.a    = 0.1;
% cauchy.b    = 0.5;
% cauchy.qcauchy = icdf('Cauchy',cauchy.p,cauchy.a,cauchy.b);
% cauchy.dcauchy = pdf('Cauchy',cauchy.x,cauchy.a,cauchy.b);
% cauchy.pcauchy = cdf('Cauchy',cauchy.q,cauchy.a,cauchy.b);

%% Gamma Distribution
gam.x       = [-2 -1 0 1 2];
gam.q       = [-2,-1,0,1,2];
gam.p       = [0.1,0.25,0.5,0.75,0.9];
gam.a       = 0.1;
gam.b       = 0.5;
gam.qgamma	= icdf('Gamma',gam.p,gam.a,gam.b);
gam.dgamma	= pdf('Gamma',gam.x,gam.a,gam.b);
gam.pgamma	= cdf('Gamma',gam.q,gam.a,gam.b);

%% Log Normal Distribution
logn.x      = [-2 -1 0 1 2];
logn.q      = [-2,-1,0,1,2];
logn.p      = [0.1,0.25,0.5,0.75,0.9];
logn.mu     = log(20000);
logn.sig    = 1;
logn.qlnorm	= icdf('Lognormal',logn.p,logn.mu,logn.sig);
logn.dlnorm = pdf('Lognormal',logn.x,logn.mu,logn.sig);
logn.plnorm = cdf('Lognormal',logn.q,logn.mu,logn.sig);

%% Poisson distribution
pois.x      = 0;
pois.q      = 1.0589;
pois.p      = 0.95;
pois.lambda = 1:10;
pois.qpois  = icdf('Poisson',pois.p,pois.lambda);
pois.dpois	= pdf('Poisson',pois.x,pois.lambda);
pois.ppois	= cdf('Poisson',pois.q,pois.lambda);

%% Weibull Distribution
weib.x      = [-2 -1 0 1 2];
weib.q      = [-2,-1,0,1,2];
weib.p      = [0.1,0.25,0.5,0.75,0.9];
weib.shape  = 0.5;
weib.scale  = 1;
weib.qweibull = icdf('Weibull',weib.p,weib.scale,weib.shape);
weib.dweibull = pdf('Weibull',weib.x,weib.scale,weib.shape);
weib.pweibull = cdf('Weibull',weib.q,weib.scale,weib.shape);

%% Matrix Tests
mat.tx    = 0;
mat.tq    = 1.0589;
mat.tdf   = magic(3);
mat.dt    = tpdf(mat.tx,mat.tdf);
mat.pt    = tcdf(mat.tq,mat.tdf);

mat.nx     = 1;
mat.nq     = 1;
mat.mu     = 0;
mat.sig    = magic(3);
mat.dnorm  = pdf('Normal',mat.nx,mat.mu,mat.sig);
mat.pnorm  = cdf('Normal',mat.nq,mat.mu,mat.sig);

%% Upper Tests
upper.tq    = 0.1;
upper.tdf   = 0.1;
upper.pt    = tcdf(upper.tq, upper.tdf, 'upper');

upper.nq    = 0.1;
upper.mu    = 0;
upper.sig   = 1;
upper.pnorm = cdf('Normal',upper.nq, upper.mu, upper.sig, 'upper');

%% Save to mat
save Utilities/UnitTests/UnitTestData/rmathlibData

author: Arrowstar


text:
function x = sample_fs(xf, grid_sz)

% Samples the Fourier series

sz = [size(xf,1) size(xf,2)];

if nargin < 2 || all(sz == grid_sz)
    x = prod(sz) * cifft2(xf);
else
    if any(grid_sz < sz)
        error('The grid size must be larger than or equal to the signal size')
    end
    tot_pad = grid_sz - sz;
    pad_sz = ceil(tot_pad/2);
    xf_pad = padarray(xf, pad_sz);
    if any(mod(tot_pad,2) == 1)
        % Handle odd padding
        xf_pad = xf_pad(1:end-mod(tot_pad(1),2), 1:end-mod(tot_pad(2),2), :, :);
    end
    x = prod(grid_sz) * cifft2(xf_pad);
end

author: martin-danelljan

repo: ECO


text:
function n = cube_arbq_size ( degree )

%*****************************************************************************80
%
%% CUBE_ARBQ_SIZE returns the size of quadrature rule for a cube.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 July 2014
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer DEGREE, the desired degree of exactness.
%    1 <= DEGREE <= 15.
%
%    Output, integer N, the number of points in the
%    corresponding rule.
%
  n_save = [ ...
      1,   4,   6,  10,  13, ...
     22,  26,  42,  50,  73, ...
     84, 116, 130, 172, 190 ];

  if ( degree < 1 | 15 < degree )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CUBE_ARBQ_SIZE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of DEGREE.\n' );
    error ( 'CUBE_ARBQ_SIZE - Fatal error!' );
  end

  n = n_save ( degree );

  return
end

author: johannesgerer

repo: jburkardt-m