text:
% UNDER CONSTRUCTION...
%
%
% Author: Javier Lopez-Calderon

function  d = emg2pulse(ch1, fs)

d=[];

if nargin<1
        help emg2pulse
        return
end
if nargin>3
        error('ERPLAB says: emg2pulse() works until 3 inputs')
end
if size(ch1,1)>1 && size(ch1,2)>1
        error('ERPLAB says: emg2pulse() works with vector inputs')
end

d = abs(ch1);
[b,a] = butter(2,30/fs);

for t=1:size(d,3)
        d(d(1,:,t)<6) = 0;
        d(1,:,t) = filtfilt(b,a, d(1,:,t));
        d(1,:,t) = 100*d(1,:,t);
        m = mean(d(1,:,t));
        sdd = std(d(1,:,t));
        th = m+2*sdd;
        d(d(1,:,t)<th) = 0;
        d(d(1,:,t)>=th) = 100;
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function output = FeaturePipe(input, processing, input2)

output = input;
for i=1:length(processing)
    switch processing{i}.name

        % data manipulations

        case 'CPU2GPU'
            output = gpuArray(output);
        case 'GPU2CPU'
            output = gather(output);
        case 'selectDim'
            output = output(processing{i}.transform,:);
        case {'splice', 'Splice'}
            output = ExpandContext_v2(output, processing{i}.transform);
        case 'stream2'  % extra input
            output = [output; input2];
        case 'segmentation'     % this is similar to enframe. It is usually used for high dimensional input, while enframe is usually for single/multichannel waveforms
            seglen = processing{i}.transform(1);
            segshift = processing{i}.transform(2);
            % if input is 1D array, output is 2D matrix. If input is 2D matrix, output is 3D tensor.The last dimension is the number of segments
            % Any remaining frames that is shorter than seglen is discarded.
            output = DivideSent2Segments(output, seglen, segshift, 0);
        case 'transpose'
            output = output';
        case 'permute'
            output = permute(output, processing{i}.transform);
        case 'vectorize'
            output = reshape(output,numel(output),1);

        case 'removeNonspeechFrame' % assume the input is something related to energy. Do energy based VAD first, then discard those frames classified as nonspeech

        case 'upsample'

        case 'downsample'

            % signal processing

        case 'addNoise'
            output = output + randn(size(output)) * processing{i}.transform;
        case 'removeDC'
            output = DC_remove(output,0.999);
        case 'preemphasis'
            output = filter([1 -0.97],1,output);
        case 'enframe'
            output = my_enframe(output, processing{i}.transform(1), processing{i}.transform(2));
        case 'windowing'
            switch processing{i}.transform
                case 'hamming'
                    window = hamming(frame_size);
                case 'hanning'
                    window = hanning(frame_size);
                otherwise
                    window = [];
            end
            if ~isempty(window)
                output = bsxfun(@times, output, window);
            end
        case 'fft'
            output = fft(output,processing{i}.transform);
        case 'log'
            if isfield(processing{i}, 'transform')
                output = log(output+processing{i}.transform);
            else
                output = log(output+eps);
            end
        case 'power'
            output = real(output .* conj(output));
        case 'complex2realImag'
            [D,nFr] = size(output);
            output = output(1:D/2,:) + sqrt(-1)*output(D/2+1:end,:);
        case 'realImag2complex'
            output = [real(output); imag(output)];

            % language processing

        case 'seq2ngram'
            output = seq2ngram(output, processing{i}.ngram, processing{i}.vocab);
        case 'idx2vec'
            output2 = sparse(processing{i}.outputDim, length(output));
            for j=1:length(output)
                output2(output(j),j)= 1;
            end
            output = output2;
        case {'context_sum', 'Context_sum'}
            context = processing{i}.transform;
            output = SumContext(output, context);
        case 'fulltify'
            output = full(output);

            % temporal processing

        case 'delta'
            output = comp_dynamic_feature(output', processing{i}.delta_order, processing{i}.delta_order)';
        case 'dynamic'
            D = genDeltaTransform(size(output,2), 2);
            A = D*D;
            output = [output; output*D'; output*A'];

            % neural networks steps

        case {'affinetransform','AffineTransform'};
            output = processing{i}.transform * output;
            if isfield(processing{i}, 'bias')
                bias = processing{i}.bias;
                output = bsxfun(@plus, output, bias(:));
            end
        case {'sigmoid','Sigmoid'}
            output = sigmoid(output);
        case {'softmax', 'Softmax'}
            output = softmax(output);
        case {'addshift', 'AddShift'}
            output = bsxfun(@minus, output, processing{i}.transform');
        case {'rescale', 'Rescale'}
            if issparse(output)
                visible_nonzero_idx = find(sum(abs(output),2)>0);
                visible_nonzero = full(output(visible_nonzero_idx,:));
                output(visible_nonzero_idx) = bsxfun(@times, visible_nonzero, processing{i}.transform(visible_nonzero_idx)');
            else
                output = bsxfun(@times, output, processing{i}.transform');
            end
        case 'linear'
            % linear activation node, do nothings
        case 'phoneID2posterior'
            output = phoneID2posterior(output, length(processing{i}.classID), processing{i}.classID);
        case 'phoneID2classID'
            output = phoneID2classID(output, processing{i}.classID);

            % normalization methods

        case 'MVN'
            output = MVN(output')';
        case 'CMN'
            output = CMN(output')';
        case 'length_norm'
            scale = 1./sqrt(sum(output.^2,1));
            output = bsxfun(@times, output, scale);

        otherwise
            fprintf('Error: unknown processing step: %s\n', processing{i}.name);
            break;
    end
end
end

author: singaxiong

repo: SignalGraph


text:
function result = swt16rec3_GPU(dec_full, level, waveletName)
%SWT16rec3 Stationary wavelet-based 2-level frame in 3-D.
%   swt16rec3_GPU performs a multi-level 3-D stationary wavelet synthesis
%   using a specific waveletName.
%
%   result = swt16rec3_GPU(dec_full, level, waveletName) computes the
%   3-D matrix 'result' at using the coefficients in 'dec_full', level
%   'level' and using 'waveletName'. Level must be either 1 or 2. Wavelet Toolbox
%   is required to retrieve the filters using 'waveletName'. Both 'level'
%   and 'waveletName' should match those used in creating 'dec_full'.
%
%   Requires the Parallel Computing Toolbox and a compatible GPU.
%
%   See conv3DGPU and swt16dec3_GPU

%   A. Kucukelbir 12-June-2012.
%   Last Revision: 26-June-2012.

% Create filters
if ischar(waveletName)
    [~,~,LR,HR] = wfilters(waveletName);
end

% Define kernels to be used by GPU
rowsKernel    = parallel.gpu.CUDAKernel('conv3DGPU.ptx', 'conv3DGPU.cu', 'rowsKernel');
columnsKernel = parallel.gpu.CUDAKernel('conv3DGPU.ptx', 'conv3DGPU.cu', 'columnsKernel');
beamsKernel   = parallel.gpu.CUDAKernel('conv3DGPU.ptx', 'conv3DGPU.cu', 'beamsKernel');
halfSum       = parallel.gpu.CUDAKernel('conv3DGPU.ptx', 'conv3DGPU.cu', 'halfSum');

% Upsample filters and create temporary cell
x = cell(1,1);
for i = 1:level-1
   LR = dyadup(LR,1);
   HR = dyadup(HR,1);
end

% Set parameters for convolution kernels on GPU
BLOCKDIM_X = 8;
BLOCKDIM_Y = 8;
BLOCKDIM_Z = 8;

for i = 1:level

    % Grab first level of frame coeffieicents.
    dec = dec_full(1:8);
    dec_full(1:8) = [];

    % Downsample filters
    if(i==2)
        LR = dyaddown(LR,0);
        HR = dyaddown(HR,0);
    end

    % Move filters to GPU
    LR_g = gpuArray(LR);
    HR_g = gpuArray(HR);

    % Calculate kernel radius
    if( mod(length(LR),2) == 0 )
        kernelRadius = length(LR)/2;
    else
        kernelRadius = (length(LR)-1)/2;
    end

    % Ensure volume has dimensions that are a multiple of 8
    diff = 0;
    if( mod(size(dec{1,1},1),8) ~= 0 )
    diff = 8*ceil(size(dec{1,1},1)/8) - size(dec{1,1},1);
        for j=1:8
        dec{j,1} = padarray(dec{j,1},[diff diff diff],'post');
        end
    end

    % Calculate volume sizes and transfer to GPU
    volRowSize = size(dec{1,1},1);
    volColSize = size(dec{1,1},2);
    volBeaSize = size(dec{1,1},3);
    dec1_g = gpuArray(dec{1,1});
    dec2_g = gpuArray(dec{2,1});
    dec3_g = gpuArray(dec{3,1});
    dec4_g = gpuArray(dec{4,1});
    dec5_g = gpuArray(dec{5,1});
    dec6_g = gpuArray(dec{6,1});
    dec7_g = gpuArray(dec{7,1});
    dec8_g = gpuArray(dec{8,1});

    % Setup GPU grids and blocks
    rowsKernel.GridSize        = [(volRowSize / BLOCKDIM_X) * (volBeaSize / BLOCKDIM_Z) , (volColSize / BLOCKDIM_Y)];
    rowsKernel.ThreadBlockSize = [BLOCKDIM_X, BLOCKDIM_Y, BLOCKDIM_Z];

    columnsKernel.GridSize        = [(volRowSize / BLOCKDIM_X) * (volBeaSize / BLOCKDIM_Z) , (volColSize / BLOCKDIM_Y)];
    columnsKernel.ThreadBlockSize = [BLOCKDIM_X, BLOCKDIM_Y, BLOCKDIM_Z];

    beamsKernel.GridSize        = [(volRowSize / BLOCKDIM_X) * (volBeaSize / BLOCKDIM_Z) , (volColSize / BLOCKDIM_Y)];
    beamsKernel.ThreadBlockSize = [BLOCKDIM_X, BLOCKDIM_Y, BLOCKDIM_Z];

    halfSum.GridSize        = [(volRowSize / BLOCKDIM_X) * (volBeaSize / BLOCKDIM_Z) , (volColSize / BLOCKDIM_Y)];
    halfSum.ThreadBlockSize = [BLOCKDIM_X, BLOCKDIM_Y, BLOCKDIM_Z];

    % Do the SWT for this level
    tmp_g = swt3_1level(dec1_g, dec2_g, dec3_g, dec4_g, dec5_g, dec6_g, dec7_g, dec8_g,...
                      LR_g, HR_g, ...
                      volRowSize, volColSize, volBeaSize, kernelRadius,...
                      rowsKernel, columnsKernel, beamsKernel, halfSum);

    % Remove multiple of 8 padding
    tmp = gather(tmp_g);
    x{i,1} = tmp(1:end-diff, 1:end-diff, 1:end-diff);

    % Remove kernelRadius padding
    if( mod(length(LR),2) == 0 )
            x{i,1} = x{i,1}(kernelRadius:end-(kernelRadius-1)-1, kernelRadius:end-(kernelRadius-1)-1, kernelRadius:end-(kernelRadius-1)-1);
    else
            x{i,1} = x{i,1}(kernelRadius+1:end-kernelRadius, kernelRadius+1:end-kernelRadius, kernelRadius+1:end-kernelRadius);
    end

    % Take care of second level specific processing
    if (i == 2)
        x{1,1} = padarray(x{1,1},[diff diff diff],'post');
        dec1_g = gpuArray(x{1,1});
            tmp_g = swt3_1level(dec1_g, dec2_g, dec3_g, dec4_g, dec5_g, dec6_g, dec7_g, dec8_g,...
                      LR_g, HR_g, ...
                      volRowSize, volColSize, volBeaSize, kernelRadius,...
                      rowsKernel, columnsKernel, beamsKernel, halfSum);
        tmp = gather(tmp_g);
        tmp = tmp(1:end-diff, 1:end-diff, 1:end-diff);
        if( mod(length(LR),2) == 0 )
            x{1,1} = tmp(kernelRadius:end-(kernelRadius-1)-1, kernelRadius:end-(kernelRadius-1)-1, kernelRadius:end-(kernelRadius-1)-1);
        else
            x{1,1} = tmp(kernelRadius:end-kernelRadius-1, kernelRadius:end-kernelRadius-1, kernelRadius:end-kernelRadius-1);
        end
    end

end

% Combine results from all levels
result = zeros(size(x{1,1}));
for i = 1:level
    result = result + 1/level*x{i,1};
end

end

function x = swt3_1level(dec1_g, dec2_g, dec3_g, dec4_g, dec5_g, dec6_g, dec7_g, dec8_g,...
                  LR_g, HR_g, ...
                  volRowSize, volColSize, volBeaSize, kernelRadius,...
                  rowsKernel, columnsKernel, beamsKernel, halfSum)

%% FIRST PASS ALONG BEAMS
% Get the AA
aaa_Lo_Lo_Lo_g = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);
daa_Lo_Lo_Hi_g = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);
aa_Lo_Lo_g     = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);

aaa_Lo_Lo_Lo_g = feval( beamsKernel, aaa_Lo_Lo_Lo_g, dec1_g, LR_g, volRowSize, volColSize, volBeaSize, kernelRadius );
daa_Lo_Lo_Hi_g = feval( beamsKernel, daa_Lo_Lo_Hi_g, dec2_g, HR_g, volRowSize, volColSize, volBeaSize, kernelRadius );

aa_Lo_Lo_g = feval( halfSum, aa_Lo_Lo_g, aaa_Lo_Lo_Lo_g, daa_Lo_Lo_Hi_g, volRowSize, volColSize, volBeaSize );

% Get the DA
ada_Lo_Hi_Lo_g = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);
dda_Lo_Hi_Hi_g = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);
da_Lo_Hi_g     = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);

ada_Lo_Hi_Lo_g = feval( beamsKernel, ada_Lo_Hi_Lo_g, dec3_g, LR_g, volRowSize, volColSize, volBeaSize, kernelRadius );
dda_Lo_Hi_Hi_g = feval( beamsKernel, dda_Lo_Hi_Hi_g, dec4_g, HR_g, volRowSize, volColSize, volBeaSize, kernelRadius );

da_Lo_Hi_g = feval( halfSum, da_Lo_Hi_g, ada_Lo_Hi_Lo_g, dda_Lo_Hi_Hi_g, volRowSize, volColSize, volBeaSize );

% Get the AD
aad_Hi_Lo_Lo_g = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);
dad_Hi_Lo_Hi_g = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);
ad_Hi_Lo_g     = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);

aad_Hi_Lo_Lo_g = feval( beamsKernel, aad_Hi_Lo_Lo_g, dec5_g, LR_g, volRowSize, volColSize, volBeaSize, kernelRadius );
dad_Hi_Lo_Hi_g = feval( beamsKernel, dad_Hi_Lo_Hi_g, dec6_g, HR_g, volRowSize, volColSize, volBeaSize, kernelRadius );

ad_Hi_Lo_g = feval( halfSum, ad_Hi_Lo_g, aad_Hi_Lo_Lo_g, dad_Hi_Lo_Hi_g, volRowSize, volColSize, volBeaSize );

% Get the DD
add_Hi_Hi_Lo_g = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);
ddd_Hi_Hi_Hi_g = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);
dd_Hi_Hi_g     = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);

add_Hi_Hi_Lo_g = feval( beamsKernel, add_Hi_Hi_Lo_g, dec7_g, LR_g, volRowSize, volColSize, volBeaSize, kernelRadius );
ddd_Hi_Hi_Hi_g = feval( beamsKernel, ddd_Hi_Hi_Hi_g, dec8_g, HR_g, volRowSize, volColSize, volBeaSize, kernelRadius );

dd_Hi_Hi_g = feval( halfSum, dd_Hi_Hi_g, add_Hi_Hi_Lo_g, ddd_Hi_Hi_Hi_g, volRowSize, volColSize, volBeaSize );

%% SECOND PASS ALONG COLUMNS
% Get the A
aa_Lo_Lo_gt = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);
da_Lo_Hi_gt = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);
a_Lo_g      = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);

aa_Lo_Lo_gt = feval( columnsKernel, aa_Lo_Lo_gt, aa_Lo_Lo_g, LR_g, volRowSize, volColSize, volBeaSize, kernelRadius );
da_Lo_Hi_gt = feval( columnsKernel, da_Lo_Hi_gt, da_Lo_Hi_g, HR_g, volRowSize, volColSize, volBeaSize, kernelRadius );

a_Lo_g = feval( halfSum, a_Lo_g, aa_Lo_Lo_gt, da_Lo_Hi_gt, volRowSize, volColSize, volBeaSize );

% Get the D
ad_Hi_Lo_gt = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);
dd_Hi_Hi_gt = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);
d_Hi_g      = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);

ad_Hi_Lo_gt = feval( columnsKernel, ad_Hi_Lo_gt, ad_Hi_Lo_g, LR_g, volRowSize, volColSize, volBeaSize, kernelRadius );
dd_Hi_Hi_gt = feval( columnsKernel, dd_Hi_Hi_gt, dd_Hi_Hi_g, HR_g, volRowSize, volColSize, volBeaSize, kernelRadius );

d_Hi_g = feval( halfSum, d_Hi_g, ad_Hi_Lo_gt, dd_Hi_Hi_gt, volRowSize, volColSize, volBeaSize );

%% FIRST PASS ALONG ROWS
a_Lo_gt = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);
d_Hi_gt = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);
x       = parallel.gpu.GPUArray.zeros(volRowSize, volColSize, volBeaSize);

a_Lo_gt = feval( rowsKernel, a_Lo_gt, a_Lo_g, LR_g, volRowSize, volColSize, volBeaSize, kernelRadius );
d_Hi_gt = feval( rowsKernel, d_Hi_gt, d_Hi_g, HR_g, volRowSize, volColSize, volBeaSize, kernelRadius );

x = feval( halfSum, x, a_Lo_gt, d_Hi_gt, volRowSize, volColSize, volBeaSize );

end

author: Sable

repo: mcbench-benchmarks


text:
function [theta_median theta_std theta_lbound theta_ubound sigma_median]=panel6estimates(d,N,n,T,theta_gibbs,sigma_gibbs,cband)

% obtain point estimates for the structural factors
% loop over sample periods
for ii=1:T
   % loop over structural factors
   for jj=1:d
   theta_median(jj,1,ii)=[quantile(theta_gibbs(jj,:,ii),0.5)];
   theta_std(jj,1,ii)=std(theta_gibbs(jj,:,ii));
   theta_lbound(jj,1,ii)=[quantile(theta_gibbs(jj,:,ii),(1-cband)/2)];
   theta_ubound(jj,1,ii)=[quantile(theta_gibbs(jj,:,ii),1-(1-cband)/2)];
   end
end

% obtain point estimates for sigma
% loop over sample periods
for ii=1:T
   % loop over sigma entries
   for jj=1:(N*n)^2
   sigma_median(jj,1,ii)=[quantile(sigma_gibbs(jj,:,ii),0.5)];
   end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function[varargout]=cms2kmd(varargin)
%CMS2KMD  Converts centimeters per second to kilometers per day.
%
%   Y=CMS2KMD(X)   <==>  Y=X*(3600*24/100/1000)
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2012--2015 J.M. Lilly --- type 'help jlab_license' for details

for i=1:nargin
    varargout{i}=(3600*24/100/1000).*varargin{i};
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
% sgwt_randmat : Compute random  (Erdos-Renyi model) graph
%
% function A=sgwt_randmat(N,thresh)
%
% Inputs :
% N - number of vertices
% thresh - probability of connection of each edge
%
% Outputs :
% A - adjacency matrix

% This file is part of the SGWT toolbox (Spectral Graph Wavelet Transform toolbox)
% Copyright (C) 2010, David K. Hammond.
%
% The SGWT toolbox 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.
%
% The SGWT toolbox 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 the SGWT toolbox.  If not, see <http://www.gnu.org/licenses/>.

function [A]=sgwt_randmat(N,thresh)
  assert(thresh<=1 && thresh>=0);
  A=rand(N)>1-thresh;
  B=triu(A);
  A=B+B';
  for i=1:size(A,1)
    A(i,i)=0;
  end
  A=sparse(A);

repo: gspbox


text:
function varargout = squeeze(varargin)
%SQUEEZE   Squeeze a DISKFUN to one variable, if possible.
%   G = squeeze(F) returns a DISKFUN if F depends on THETA and R in
%   polar coordinates. If F depends only on the theta-variable, a row
%   CHEBFUN is returned. If it depends on just the r-variable, a
%   column CHEBFUN is returned.

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

[varargout{1:nargout}] = squeeze@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [Jul1,Jul2]=UTC2TCG(Jul1,Jul2)
%%UTC2TCG  Convert from universal coordinated time (UTC) given as a
%          two-part pseudo-Julian date to geocentric coordinate time (TCG),
%          represented as a two-part Julian date.
%
%INPUTS: Jul1, Jul2 Two parts of a pseudo-Julian date given in UTC. The
%                   units of the date are days. The full date is the sum of
%                   both terms. The date is broken into two parts to
%                   provide more bits of precision. It does not matter how
%                   the date is split.
%
%OUTPUTS: Jul1, Jul2 The time as a Julian date in TCG.
%
%The UTC date is only pseudo-Julian, because there is not a fixed number
%of seconds in a Julian day. The convention used in the IAU standard is
%that the Julian day matches the UTC day regardless of whether the UTC day
%is 86399, 86400 or 86401 SI seconds (depending on the presence of leap
%seconds).
%
%UTC began at 1960 January 1.0 (JD 2436934.5) and this function should not
%be called with an earlier date.
%
%This just calls a number of intermediate conversion functions out of the
%International Astronomical Union's (IAU) Standard's of Fundamental
%Astronomy library.
%
%Many temporal coordinate systems standards are compared in [1].
%
%REFERENCES:
%[1] D. F. Crouse, "An Overview of Major Terrestrial, Celestial, and
%    Temporal Coordinate Systems for Target Tracking," Formal Report, Naval
%    Research Laboratory, no. NRL/FR/5344--16-10,279, 10 Aug. 2016, 173
%    pages.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

[Jul1,Jul2]=UTC2TAI(Jul1,Jul2);
[Jul1,Jul2]=TAI2TT(Jul1,Jul2);
[Jul1,Jul2]=TT2TCG(Jul1,Jul2);
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:
%RandomPath Vehicle driver class
%
% Create a "driver" object capable of steering a Vehicle object through random
% waypoints within a rectangular region and at constant speed.
%
% The driver object is connected to a Vehicle object by the latter's
% add_driver() method.  The driver's demand() method is invoked on every
% call to the Vehicle's step() method.
%
% Methods::
%  init       reset the random number generator
%  demand     return speed and steer angle to next waypoint
%  display    display the state and parameters in human readable form
%  char       convert to string
%
% Properties::
%  goal          current goal/waypoint coordinate
%  veh           the Vehicle object being controlled
%  dim           dimensions of the work space (2x1) [m]
%  speed         speed of travel [m/s]
%  closeenough   proximity to waypoint at which next is chosen [m]
%
% Example::
%
%    veh = Vehicle(V);
%    veh.add_driver( RandomPath(20, 2) );
%
% Notes::
% - It is possible in some cases for the vehicle to move outside the desired
%   region, for instance if moving to a waypoint near the edge, the limited
%   turning circle may cause the vehicle to temporarily move outside.
% - The vehicle chooses a new waypoint when it is closer than property
%   closeenough to the current waypoint.
% - Uses its own random number stream so as to not influence the performance
%   of other randomized algorithms such as path planning.
%
% Reference::
%
%   Robotics, Vision & Control, Chap 6,
%   Peter Corke,
%   Springer 2011
%
% See also Vehicle.

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

% TODO
%  should be a subclass of VehicleDriver
%  Vehicle should be an abstract superclass

classdef RandomPath < handle
    properties
        goal        % current goal
        h_goal      % graphics handle for goal
        veh         % the vehicle we are driving
        dim
        speed       % speed of travel
        closeenough  % proximity to goal before
        d_prev
        randstream  % random stream just for Sensors
    end

    methods

        function driver = RandomPath(dim, varargin)
        %RandomPath.RandomPath Create a driver object
        %
        % D = RandomPath(DIM, OPTIONS) returns a "driver" object capable of driving
        % a Vehicle object through random waypoints.  The waypoints are positioned
        % inside a rectangular region bounded by +/- DIM in the x- and y-directions.
        %
        % Options::
        % 'speed',S      Speed along path (default 1m/s).
        % 'dthresh',D    Distance from goal at which next goal is chosen.
        %
        % See also Vehicle.

        % TODO options to specify region, maybe accept a Map object?

            driver.dim = dim;

            opt.speed = 1;
            opt.dthresh = 0.05 * dim;
            opt = tb_optparse(opt, varargin);

            driver.speed = opt.speed;
            driver.closeenough = opt.dthresh;
            drive.d_prev = Inf;
            driver.randstream = RandStream.create('mt19937ar');
        end

        function init(driver)
        %RandomPath.init Reset random number generator
        %
        % R.init() resets the random number generator used to create the waypoints.
        % This enables the sequence of random waypoints to be repeated.
        %
        % See also RANDSTREAM.
            driver.goal = [];
            driver.randstream.reset();
        end

        % not used
        function visualize(driver)
            clf
            d = driver.dim;
            axis([-d d -d d]);
            hold on
            xlabel('x');
            ylabel('y');
        end

        % private method, invoked from demand() to compute a new waypoint
        function setgoal(driver)
            r = driver.randstream.rand(2,1);
            driver.goal = 0.8 * driver.dim * (r - 0.5)*2;
            %fprintf('set goal: (%.1f %.1f)\n', driver.goal);
            if isempty(driver.h_goal)
                %driver.h_goal = plot(driver.goal(1), driver.goal(2), '*')
            else
                %set(driver.h_goal, 'Xdata', driver.goal(1), 'Ydata', driver.goal(2))
            end
        end

        function [speed, steer] = demand(driver)
        %RandomPath.demand Compute speed and heading to waypoint
        %
        % [SPEED,STEER] = R.demand() returns the speed and steer angle to
        % drive the vehicle toward the next waypoint.  When the vehicle is
        % within R.closeenough a new waypoint is chosen.
        %
        % See also Vehicle.
            if isempty(driver.goal)
                driver.setgoal()
            end

            speed = driver.speed;

            goal_heading = atan2(driver.goal(2)-driver.veh.x(2), ...
                driver.goal(1)-driver.veh.x(1));
            d_heading = angdiff(goal_heading, driver.veh.x(3));

            steer = d_heading;

            % if nearly at goal point, choose the next one
            d = colnorm(driver.veh.x(1:2) - driver.goal);
            if d < driver.closeenough
                driver.setgoal();
            elseif d > driver.d_prev
                driver.setgoal();
            end
            driver.d_prev = d;
        end

        function display(driver)
        %RandomPath.display Display driver parameters and state
        %
        % R.display() displays driver parameters and state in compact
        % human readable form.
        %
        % See also RandomPath.char.
            loose = strcmp( get(0, 'FormatSpacing'), 'loose');
            if loose
                disp(' ');
            end
            disp([inputname(1), ' = '])
            disp( char(driver) );
        end % display()

        function s = char(driver)
        %RandomPath.char Convert to string
        %
        % s = R.char() is a string showing driver parameters and state in in
        % a compact human readable format.
            s = 'RandomPath driver object';
            s = char(s, sprintf('  current goal=(%g,%g), dimension %.1f', ...
                driver.goal, driver.dim));
        end

    end % methods
end % classdef

author: Allopart

repo: rbpf-gmapping


text:
classdef hmx
%+========================================================================+
%|                                                                        |
%|         OPENHMX - LIBRARY FOR H-MATRIX COMPRESSION AND ALGEBRA         |
%|           openHmx is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : hmx.m                                         |
%|    #    |   VERSION    : 0.61                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 05.09.2019                                    |
%| ( === ) |   SYNOPSIS   : H-Matrix class definition and functions       |
%|  `---'  |                                                              |
%+========================================================================+

properties
    typ = [];            % LEAF TYPE (0=H-MATRIX ; 1=COMPRESSED ; 2=FULL/SPARSE)
    pos = [];            % COORDINATES POSITIONS (X,Y)
    row = cell(1,4);     % CHILDREN ROWS INDICES (M11,M12,M21,M22)
    col = cell(1,4);     % CHILDREN COLUMNS INDICES
    chd = cell(1,4);     % CHILDREN
    dat = [];            % LEAF DATA
    tol = [];            % COMPRESSORS ACCURACY
end

methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%% CONSTRUCTOR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CONSTRUCTOR
    function Mh = hmx(varargin)
        % Empty object
        if (nargin == 0)

        % Initialization other H-matrix
        elseif (nargin == 2)
            Mh = hmxCopy(varargin{1},varargin{2});

        % Initialization with dimension and accuracy
        elseif (nargin == 3)
            Mh     = hmx();
            Mh.pos = {varargin{1},varargin{2}};
            Mh.tol = varargin{3};

        % Particles builder with partial and total pivoting
        elseif (nargin == 4)
            X     = varargin{1};
            Y     = varargin{2};
            green = varargin{3};
            acc   = varargin{4};
            Mh    = hmxBuilder(X,Y,green,acc);

        % Compressed builder
        elseif (nargin == 5)
            X      = varargin{1};
            Y      = varargin{2};
            A      = varargin{3};
            B      = varargin{4};
            acc    = varargin{5};
            Mh     = hmx(X,Y,acc);
            Mh.typ = 1;
            Mh.dat = {A,B};

       % Finite element builder
        elseif (nargin == 8)
            Xunk  = varargin{1};
            Yunk  = varargin{2};
            Mx    = varargin{3};
            X     = varargin{4};
            green = varargin{5};
            Y     = varargin{6};
            My    = varargin{7};
            acc   = varargin{8};
            Mh    = hmxBuilderFem(Xunk,Yunk,Mx,X,green,Y,My,acc);

        else
            error('hmx.m : undefined constructor case')
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%% GLOBAL DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % SIZE
    function s = size(varargin)
        Mh = varargin{1};
        s  = [size(Mh.pos{1},1),size(Mh.pos{2},1)];
        if (nargin == 2)
            s = s(varargin{2});
        end
    end

    % LENGTH
    function l = length(Mh)
        l = max(size(Mh));
    end

    % ISLOWER
    function b = islower(Mh)
        if (Mh.typ == 0)
            b = (Mh.chd{2}.typ == 1) && isempty(Mh.chd{2}.dat{1});
        elseif (Mh.typ == 2)
            b = 1;
        else
            error('hmx.m : unavailable case.')
        end
    end

    % ISUPPER
    function b = isupper(Mh)
        if (Mh.typ == 0)
            b = (Mh.chd{3}.typ == 1) && isempty(Mh.chd{3}.dat{1});
        elseif (Mh.typ == 2)
            b = 1;
        else
            error('hmx.m : unavailable case.')
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%% VISUALISATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % STRUCTURE VISUALISATION
    function spy(Mh)
        hmxSpy(Mh);
    end

    % PLOT POSITIONS
    function plot3(Mh)
        plot3(Mh.pos{1}(:,1),Mh.pos{1}(:,2),Mh.pos{1}(:,3),'bo',...
            Mh.pos{2}(:,1),Mh.pos{2}(:,2),Mh.pos{2}(:,3),'*r');
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CONVERSION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % FULL
    function M = full(varargin)
        if (nargin == 1)
            M = hmxFull(varargin{1});
        else
            M = hmxFullSub(varargin{1},varargin{2},varargin{3});
        end
    end

    % SPARSE
    function M = sparse(varargin)
        if (nargin == 1)
            M = hmxSparse(varargin{1});
        elseif (nargin == 2)
            M = hmxSparsify(varargin{1},varargin{2});
        else
            M = hmxSparseSub(varargin{1},varargin{2},varargin{3});
        end
    end

    % LOW-RANK
    function [A,B] = lowrank(varargin)
        if (nargin == 1)
            [A,B] = hmxLowrank(varargin{1});
        else
            [A,B] = hmxLowrankSub(varargin{1},varargin{2},varargin{3});
        end
    end

    % DIAGONAL
    function D = diag(Mh)
        D = hmxDiag(Mh,(1:size(Mh,1))',(1:size(Mh,2))');
        D = full(diag(D));
    end

    % ZEROS
    function Mh = zeros(Mh)
        Mh.typ = 1;
        Mh.row = cell(1,4);
        Mh.col = cell(1,4);
        Mh.chd = cell(1,4);
        Mh.dat = {zeros(size(Mh,1),0),zeros(0,size(Mh,2))};
    end

    % ONES
    function Mh = ones(Mh)
        Mh.typ = 1;
        Mh.row = cell(1,4);
        Mh.col = cell(1,4);
        Mh.chd = cell(1,4);
        Mh.dat = {ones(size(Mh,1),1),ones(1,size(Mh,2))};
    end

    % SINGLE
    function Mh = single(Mh)
        Mh = hmxSingle(Mh);
    end

    % DOUBLE
    function Mh = double(Mh)
        Mh = hmxDouble(Mh);
    end

    % TRANSPOSITION
    function Mh = transpose(Mh)
        Mh = hmxTranspose(Mh);
    end

    % TRANSPOSITION CONJUGATE
    function Mh = ctranspose(Mh)
        Mh = hmxCtranspose(Mh);
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ALGEBRA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % SCALAR PRODUCT
    function Mh = times(Ml,Mr)
        Mh = hmxTimes(Ml,Mr);
    end

    % UMINUS
    function Mh = uminus(Mh)
        Mh = (-1).*Mh;
    end

    % ADDITION
    function Mh = plus(Ml,Mr)
        Mh = hmxPlus(Ml,Mr);
    end

    % SUBSTRACTION
    function Mh = minus(Ml,Mr)
       Mh = Ml + (-Mr);
    end

    % MATRIX PRODUCT
    function Mh = mtimes(Ml,Mr)
        if isnumeric(Ml) && (numel(Ml) == 1)
            Mh = Ml .* Mr;
        elseif isnumeric(Mr) && (numel(Mr) == 1)
            Mh = Ml .* Mr;
        else
            Mh = hmxMtimes(Ml,Mr);
        end
    end

    % ADDITION WITH MATRIX PRODUCT
    function Mh = plusmtimes(Mh,alpha,Ml,Mr)
        Mh = hmxPlusMtimes(Mh,alpha,Ml,Mr);
    end

    % INVERSION
    function Mh = inv(Mh)
        Mh = hmxInv(Mh);
    end

    % CHOLESKY FACTORISATION
    function Mh = chol(Mh)
        Mh = hmxChol(Mh);
    end

    % LDLt FACTORISATION
    function [Mh,Dh] = ldl(Mh)
        [Mh,Dh] = hmxLdl(Mh);
    end

    % LU FACTORISATION
    function [Lh,Uh] = lu(Mh)
        [Lh,Uh] = hmxLU(Mh);
    end

    % MLDIVIDE
    function B = mldivide(Mh,B)
       B = hmxMldivide(Mh,B);
    end

    % MRDIVIDE
    function B = mrdivide(B,Mh)
        B = (Mh.'\B.').';
    end

    % VERTICAL CONCATENATION
    function Mh = vertcat(varargin)
        Mh = varargin{1}.';
        for i = 2:nargin
            Mh = hmxHorzcat(Mh,varargin{i}.');
        end
        Mh = Mh.';
    end

    % HORIZONTAL CONCATENATION
    function Mh = horzcat(varargin)
        Mh = varargin{1};
        for i = 2:nargin
            Mh = hmxHorzcat(Mh,varargin{i});
        end
    end
end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [] = aps_spectrometer_PWV_modis(batchfile)
% aps_spectrometer_PWV_modis(batchfile)
% Scipt to load modis data, mask out clouds. The modis data is assumed to
% be structured in date folders. The batchfile contains the full path to the
% modis files in these folders. Note that the first line of the batchfile
% should read "files".
%

% INPUTS used:
% batchfile             A txt file containing the full path and file names of the
%                       meris data that needs to be processed. The first
%                       line of this file should read "files". The data
%                       should be structured in date folders.
% xlims                 Limits in the x-direction, either in degrees
% ylims                 Limits in the y-direction, either in degrees
% smpres                The output resolution, either in degrees
%                       Units needs to be consistend with xlims and ylims.
%
% By David Bekaert - University of Leeds
% August 2014

% setting the defaults and checking the input arguments
if nargin<1
    fprintf('aps_spectrometer_PWV_modis(batchfile) \n')
    error('myApp:argChk', ['Not enough input arguments...\n'])
end

smpres = getparm_aps('region_res',1);          % in degrees
xlims = getparm_aps('region_lon_range',1);
ylims = getparm_aps('region_lat_range',1);
stamps_processed = getparm_aps('stamps_processed',1);

%% the actual scripting
%bounds for all ifgms in degrees or in meters
xmin = xlims(1);
xmax = xlims(2);
ymin = ylims(1);
ymax = ylims(2);

% getting the number of files to be processed
files = char(textread(batchfile,'%s','headerlines',1));
ix_keep = [];
for counter=1:size(files,1)
    if isempty(strfind(files(counter,:),'recal'));
       ix_keep = [ix_keep; counter];
    end
end
files = files(ix_keep,:);
% exclude recalibrated files
ndates = size(files,1);

% loading the date information
if strcmp(stamps_processed,'y')
   ps = load(getparm_aps('ll_matfile',1));
   ifgs_dates = ps.day;
   fprintf('Stamps processed structure \n')
else
    ifgday_matfile = getparm_aps('ifgday_matfile',1);
    ifgs_dates = load(ifgday_matfile);
    ifgs_dates = ifgs_dates.ifgday;
    ifgs_dates = reshape(ifgs_dates,[],1);
    ifgs_dates = unique(ifgs_dates);
end

fprintf(['Lon range: ' num2str(xmin) ' -- '  num2str(xmax) ' degrees\n'])
fprintf(['Lat range: ' num2str(ymin) ' -- '  num2str(ymax) ' degrees\n'])
fprintf(['Output resolution is assumed to be ' num2str(smpres) ' degrees \n'])

% extracting the dates from the filenames
for k=1:ndates
    [path,filename_temp,ext_temp] = fileparts(files(k,:));
    clear filename_temp ext_temp
    [path_temp,date,ext_temp] = fileparts(path);
    clear path_temp ext_temp

    % save the paths as structures to allow for variable path lengths
    pathlist{k} = path;
    clear path
    % saving the date information
    datelist(k,:) =date;
    clear date
end

%start loop here to calculate atmos correction for each date
fprintf('Starting the masking and writing of the PWV for each SAR date \n')
for n = 1:ndates

    file = [files(n,:)];
    outfile_watervapor = [pathlist{n} filesep datelist(n,:) '_ZPWV_nointerp.xyz'];
    [lon_grid,lat_grid,watervap] = grdread2(file);          % g/cm^2 water vapour

    D = grdinfo2(file);
    mod_xmin = D(1);
    mod_xmax = D(2);
    mod_ymin = D(3);
    mod_ymax = D(4);
    mod_cols = length(lon_grid);
    mod_rows = length(lat_grid);
    clear D

    %% Calculate wet delay
    %mask out dodgy pixels
    grdwrite2(lon_grid,lat_grid,watervap,'tmp.grd')

    % downsample and output without interpolation
    grdsmp_cmd = ['grdsample -R',num2str(xmin),'/',num2str(xmax),'/',num2str(ymin),'/',num2str(ymax),' -I',num2str(smpres),' -F tmp.grd -Gtmp_smp.grd'];
    [a,b]=system(grdsmp_cmd);
    grd2xyz_cmd = ['grd2xyz -R',num2str(xmin),'/',num2str(xmax),'/',num2str(ymin),'/',num2str(ymax),' tmp_smp.grd -bo >' outfile_watervapor];
    [a,b]=system(grd2xyz_cmd);

    % load nointerp back in and write out
    nointfid = fopen(outfile_watervapor,'r');
    data_vector = fread(nointfid,'double');
    fclose(nointfid);
    % reshaping into the right n column matrix
    data = reshape(data_vector,3,[])';
    noint = data(:,3);
    xy = data(:,[1:2]);
    clear data data_vector

%     figure; scatter3(xy(:,1),xy(:,2),noint,15,noint,'filled'); view(0,90); axis equal; axis tight

    %output
    data_write = [xy noint]';
    fid = fopen(outfile_watervapor,'w');
    fwrite(fid,data_write,'double');
    fclose(fid);
    clear data_write

    fprintf([num2str(n) ' completed out of ' num2str(ndates) '\n'])
end

[a,b] = system('rm tmp.grd tmp.xyz tmp2.grd tmp_fil.grd tmp_smp.grd tmp_smp2.grd tmp_smp_fil.grd');

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
function At = ctranspose(A)
At.func = A.transpose;
At.transpose = A.func;
At = class(At,'A_operator');

author: thomaskuestner


text:
function varargout = size( x, dim )

%   Disciplined convex/geometric programming information for SIZE:
%       SIZE imposes no convexity restrictions on its first argument.
%       The second argument, if supplied, must be a positive integer.

s = x.size_;
if nargin > 1,
    if nargout > 1,
        error( 'Too many output arguments.' );
    elseif ~isnumeric( dim ) || length( dim ) ~= 1 || dim <= 0 || dim ~= floor( dim ),
        error( 'Dimension argument must be a positive integer scalar.' );
    elseif dim > length( s ),
        varargout{1} = 1;
    else
        varargout{1} = s(dim);
    end
elseif nargout > 1,
    ns = length( s );
    no = nargout;
    if no > ns,
        s( end+1:no ) = 1;
    elseif no < ns,
        s( no ) = prod( s( no : end ) );
    end
    for k = 1 : no,
        varargout{k} = s(k); %#ok
    end
else
    varargout{1} = s;
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
%% HELLO_BATCH_LOCAL uses the BATCH command to run the HELLO code locally.
%
%  Discussion:
%
%    The HELLO code is a function, so first we must write a script
%    called HELLO_SCRIPT that runs the function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2011
%
%  Author:
%
%    John Burkardt
%
  clear

  workers = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HELLO_BATCH_LOCAL\n' );
  fprintf ( 1, '  Run HELLO_SCRIPT locally.\n' );
%
%  The BATCH command sends the script for execution.
%
  my_job = batch ( 'hello_script', ...
    'Configuration', 'local', ...
    'CaptureDiary', true, ...
    'FileDependencies', { 'hello_fun' }, ...
    'matlabpool', workers );
%
%  WAIT pauses the MATLAB session til the job completes.
%
  wait ( my_job );
%
%  DIARY displays any messages printed during execution.
%
  diary ( my_job );
%
%  Here is where we would ordinarily use the LOAD command
%  to make the script's workspace available.
%  This job does not produce any data in the workspace,
%  so no LOAD is needed.
%
%  load ( my_job );
%
%  DESTROY cleans up data about the job we no longer need.
%
  destroy ( my_job );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HELLO_BATCH_LOCAL\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

author: johannesgerer

repo: jburkardt-m


text:
function inds = meshBoundaryVertexIndices(varargin)
%MESHBOUNDARYVERTEXINDICES Indices of boundary vertices of a mesh.
%
%   INDS = meshBoundaryVertexIndices(V, F)
%   INDS = meshBoundaryVertexIndices(V, E, F)
%
%   Example
%     % create centered icosahedron
%     [v, f] = createIcosahedron;
%     v(:,3) = v(:,3) - mean(v(:,3));
%     % convert to simili-sphere
%     [v2, f2] = subdivideMesh(v, f, 3);
%     v3 = normalizeVector3d(v2);
%     % clip with plane
%     plane = createPlane([0 0 0], [-1 -2 3]);
%     [vc, fc] = clipMeshVertices(v3, f2, plane, 'shape', 'plane');
%     figure; drawMesh(vc, fc); axis equal; view(3);
%     % draw boundary vertices
%     inds = meshBoundaryVertexIndices(vc, fc);
%     hold on; drawPoint3d(vc(inds,:), 'k*');
%
%   See also
%     meshes3d, meshBoundary, meshBoundaryEdgeIndices, meshEdgeFaces

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2019-05-01, using Matlab 8.6.0.267246 (R2015b)
% Copyright 2019-2022 INRA - Cepia Software Platform

[vertices, edges, faces] = parseMeshData(varargin{:});

% Compute edge-vertex map if not specified
if isempty(edges)
    edges = meshEdges(vertices, faces);
end

% compute edges to faces map
edgeFaces = meshEdgeFaces(vertices, edges, faces);

borderEdges = sum(edgeFaces == 0, 2) > 0;

inds = edges(borderEdges, :);
inds = unique(inds(:));

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
classdef (Abstract) ClustererX < BaseX  % Extends trackingX.BaseX
% ClustererX Abstract class
%
% Summary of ClustererX:
% This is the base class for all TrackingX clusterers.
% Any custom defined Clusterer should be derived from this ClustererX base class.
%
% ClustererX Properties:
%   ClusterList - A maintained list containing the last generated clusters
%   UnassocTrackInds - A maintained list of all tracks not associated to
%                      any measurements
%
% ClustererX Methods:
%   + ClustererX - Constructor method
%
% (+) denotes puplic properties/methods
%
% February 2018 Lyudmil Vladimirov, University of Liverpool.

    properties
        ClusterList
        UnassocTrackInds
    end

    methods (Abstract)
        cluster(this);
    end
    methods
        function this = ClustererX(varargin)
        % CLUSTERERX Constructor method
        %

        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function [V,A,flag] = symeig(S,d,m)
% Compute eigenvalues and eigenvectors of symmetric matrix
%   m == 's' smallest (default)
%   m == 'l' largest
% Written by Mo Chen (sth4nth@gmail.com).
if nargin == 2
    m = 's';
end
opt.disp = 0;
opt.issym = 1;
opt.isreal = 1;
if any(m == 'ls')
    [V,A,flag] = eigs(S,d,[m,'a'],opt);
else
    error('The third parameter must be l or s.');
end

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function kern = rbfKernParamInit(kern)

% RBFKERNPARAMINIT RBF kernel parameter initialisation.
% The radial basis function kernel (RBF) is sometimes also known as
% the squared exponential kernel. It is a very smooth non-linear
% kernel and is a popular choice for generic use.
%
% k(x_i, x_j) = sigma2 * exp(-gamma/2 *(x_i - x_j)'*(x_i - x_j))
%
% The parameters are sigma2, the process variance (kern.variance)
% and gamma, the inverse width (kern.inverseWidth). The inverse
% width controls how wide the basis functions are, the larger
% gamma, the smaller the basis functions are.
%
% There is also an automatic relevance determination version of
% this kernel provided.
%
% SEEALSO : rbfardKernParamInit
%
% FORMAT
% DESC initialises the radial basis function
%  kernel structure with some default parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006
%
% COPYRIGHT : Antti Honkela, 2012

% KERN

kern.inverseWidth = 1;
kern.variance = 1;
kern.nParams = 2;

% Constrains parameters positive for optimisation.
if (isfield(kern,'options')) && ...
   (isfield(kern.options,'boundedParam')) && kern.options.boundedParam,
  for k=1:2,
    kern.transforms(k).index = k;
    kern.transforms(k).type = optimiDefaultConstraint('bounded');
    kern.transforms(k).transformsettings = [0 1e6];
  end
else
  kern.transforms.index = [1 2];
  kern.transforms.type = optimiDefaultConstraint('positive');
end
kern.isStationary = true;
kern.isNormalised = false;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function update_num_c = update_cft(i_span,n_times,min_num_c,mid_num_c,max_num_c,num_c)

% determine the time series model
if i_span < mid_num_c*n_times
    % start with 4 coefficients model
    update_num_c = min(min_num_c,num_c);
elseif i_span < max_num_c*n_times
    % start with 6 coefficients model
    update_num_c = min(mid_num_c,num_c);
else
    % start with 8 coefficients model
    update_num_c =  min(max_num_c,num_c);
end

end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
function S=hypergeometric0F1(a,z)
%%HYPERGEOMETRIC0F1 Evaluate the confluent hypergeometric function 0F1 for
%                   real or complex parameters. This is the confluent
%                   hypergeometric limit function 0F1(;a;z).
%
%INPUTS: a A real or complex scalar value.
%        z A real or complex scalar value.
%
%OUTPUTS: S The value 0F1(;a;z).
%
%This implements the Taylor series expansion of APpendix I of [1]. It was
%noted that the series generated at least 12 digits of accuracy for
%abs(z)<=1000. A check was added for values of a that are real and negtive
%as the solution should be infinite in that case. However, neighboring
%solutions are finite.
%
%REFERENCES:
%[1] J. Pearson, "Computation of hypergeometric functions," Master's
%    thesis, University of Oxford, Worcester College, 4 Sep. 2009.
%
%July 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Check for the special case.
if(isreal(a)&&a<0)
   S=Inf;
   return;
end

A=1;
S=A;

epsVal=Inf;
k=0;
while(epsVal>eps(S))
    A=A*(1/(a+k))*(z/(k+1));
    S=S+A;

    epsVal=abs(A);

    k=k+1;
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.


text:
% TESTBED_PLOT_STATIONS - Plots the stations in the Testbed dataset.
%
% TESTBED_PLOT_STATIONS(COORDINATES, ...)
%
% COORDINATES is a matrix that has longitudes as the first row and latitudes
% as the second row. It can also be a struct which has COORDINATES as a
% field: COORDINATES.coordinates.

% Last modified 2010-06-11
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function testbed_plot_stations(coordinates, varargin)

% Default options
options = struct(...
    'Style', 'o', ...
    'MarkerSize', 4, ...
    'MarkerEdgeColor', 'k', ...
    'MarkerFillColor', [0.7 0 0.7]);
[options, errmsg, remopts] = argparse(options, varargin{:});
error(errmsg);

if nargin < 1
  % Default set of stations
  data = testbed_loaddata();
  data = testbed_preprocess(data);
  coordinates = data.coordinates;
end

if ~isnumeric(coordinates) && isstruct(coordinates)
  coordinates = coordinates.coordinates;
end

% Plot the stations
map_plot(coordinates', ...
         options.Style, ...
         'MarkerSize',      options.MarkerSize, ...
         'MarkerEdgeColor', options.MarkerEdgeColor, ...
         'MarkerFaceColor', options.MarkerFillColor, ...
         remopts{:});

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function tests = test_spm_cat_struct
% Unit Tests for spm_cat_struct
%__________________________________________________________________________
% Copyright (C) 2017 Wellcome Trust Centre for Neuroimaging

% $Id: test_spm_cat_struct.m 7074 2017-05-05 10:51:08Z guillaume $

tests = functiontests(localfunctions);

function test_spm_cat_struct_1(testCase)
s1 = struct('a',1,'b',2);
s2 = struct('a',3,'b',4);

exp = struct('a',{1,3},'b',{2,4});
act = spm_cat_struct(s1,s2);
testCase.verifyEqual(act, exp);

function test_spm_cat_struct_2(testCase)
s1 = struct('a',1);
s2 = struct('b',2);

exp = struct('a',{1,[]},'b',{[],2});
act = spm_cat_struct(s1,s2);
testCase.verifyEqual(act, exp);

function test_spm_cat_struct_3(testCase)
s1 = struct('a',{1,2});
s2 = struct('b',3);

exp = struct('a',{1,2,[]},'b',{[],[],3});
act = spm_cat_struct(s1,s2);
testCase.verifyEqual(act, exp);

exp = struct('a',{[],1,2},'b',{3,[],[]});
act = spm_cat_struct(s2,s1);
testCase.verifyEqual(act, exp);

function test_spm_cat_struct_4(testCase)
s1 = struct('a',{1;2});
s2 = struct('b',3);

exp = struct('a',{1;2;[]},'b',{[];[];3});
act = spm_cat_struct(s1,s2);
testCase.verifyEqual(act, exp);

exp = struct('a',{[],1,2},'b',{3,[],[]});
act = spm_cat_struct(s2,s1);
testCase.verifyEqual(act, exp);

function test_spm_cat_struct_5(testCase)
s1 = struct('a',{1,2});
s2 = struct('b',3);
s3 = struct('a',4,'c',5);

exp = struct('a',{1,2,[],4},'b',{[],[],3,[]},'c',{[],[],[],5});
act = spm_cat_struct(s1,s2,s3);
testCase.verifyEqual(act, exp);

function test_spm_cat_struct_6(testCase)
s1 = struct('a',{1,2;3 4},'b',{5,6;7,8});
s2 = struct('b',{0,1},'c',{3,2});

exp = struct('a',{1,3,2,4,[],[]},'b',{5,7,6,8,0,1},'c',{[],[],[],[],3,2})';
act = spm_cat_struct(s1,s2);
testCase.verifyEqual(act, exp);

author: spm


text:
% TEST_MAXWELL_SQUARE_H_DRCHLT: data function for Dirichlet boundary condition.

function h = test_maxwell_square_h_drchlt (x, y, ind)

  h = zeros (size (x));
  switch (ind)
    case 1
      h = exp(x) .* cos(y);
    case 2
      h = -exp(x) .* cos(y);
    case 3
      h = -sin(y);
    case 4
      h = sin(y);
    otherwise
      error ('h_drchlt: unknown reference number')
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function a = eq( x, y )

%Disciplined convex programming information for EQ (==):
%   Both the left- and right-hand sides of an equality constraint must
%   be affine (or constant). If either side of the constraint is complex,
%   then the real and imaginary portions are constrained separately.
%
%Disciplined geometric programming information for EQ (>):
%   Both the left- and right-hand sides of an equality constraint must
%   be log-affine, which includes positive constants and monomials.

b = newcnstr( evalin( 'caller', 'cvx_problem', '[]' ), x, y, '==' );
if nargout, a = b; end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.


text:
function variance = gumbel_variance ( )

%*****************************************************************************80
%
%% GUMBEL_VARIANCE returns the variance of the Gumbel PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VARIANCE, the variance of the PDF.
%
  variance = pi * pi / 6.0;

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Y = OneLayerAGH_Test(X, Anchor, W, s, sigma, options)
%
% One-Layer Anchor Graph Hashing Test
% Written by Wei Liu (wliu@ee.columbia.edu)
% X = nXdim input data
% Anchor = mXdim anchor points (m<<n)
% W = mXr projection matrix in spectral space
% s = number of nearest anchors
% sigma: Gaussian RBF kernel width parameter
% Y = nXr binary codes (Y_ij in {1,0})
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[n,dim] = size(X);
m = size(Anchor,1);

%% get Z
Z = zeros(n,m);
Dis = sqdist(X',Anchor');
%clear X;
%clear Anchor;

val = zeros(n,s);
pos = val;
for i = 1:s
    [val(:,i),pos(:,i)] = min(Dis,[],2);
    tep = (pos(:,i)-1)*n+[1:n]';
    Dis(tep) = 1e60;
end
clear Dis;
clear tep;

switch lower(options.CodingMethod)
    case {lower('Cosine')}
        if isfield(options,'bNormalized') && options.bNormalized
        else
            X = NormalizeFea(X);
        end
        Anchor = NormalizeFea(Anchor);
        S = full(X*Anchor');
        linidx = [1:size(pos,1)]';
        val = S(sub2ind(size(S),linidx(:,ones(1,size(pos,2))),pos));
    case {lower('Gaussian')}
        val = exp(-val/(1/1*sigma^2));
    otherwise
        error('method does not exist!');
end

val = repmat(sum(val,2).^-1,1,s).*val; %% normalize
tep = (pos-1)*n+repmat([1:n]',1,s);
Z([tep]) = [val];
Z = sparse(Z);
clear tep;
clear val;
clear pos;

%% get binary codes
Y = (Z*W>0);  %% logical format
clear Z;

author: ZJULearning

repo: MatlabFunc


text:
%
% Generate new exemplars.
%
% Rather than re-sampling from the original log-probability scores,
% take the rank score and sample according to 1/sigma(i)
% where sigma(i) is the rank score from highest to lowest score.
%
% Input
%  G: result of model fitting
%  lib: library
%  nsamp: (default=1) number of samples
%
% Output
%  samples: [nsamp x 1 cell] samples
%  types: [nsamp x 1 cell] samples before token-level is resampled
%
function [samples,types] = task_generate_exemplars_1overk(G,lib,nsamp)

    if ~exist('nsamp','var')
       nsamp = 1;
    end

    % re-score the hypotheses by their rank score instead
    % of the real score
    wts = rescore_by_rank(G.scores);
    G.scores = log(wts);

    [samples,types] = task_generate_exemplars(G,lib,nsamp);
end

%
% Input
%  G: mixture model struct
%  lib: library
%  nsamp: (default=1) number of samples
%
% Output
%  samples: [nsamp x 1 cell] samples
%  types: [nsamp x 1 cell] samples before token-level is resampled
%
function [samples,types] = task_generate_exemplars(G,lib,nsamp)

    if ~exist('nsamp','var')
       nsamp = 1;
    end

    % probability of choosing each parse
    logwts = G.scores;
    wts = exp(logwts-logsumexp(logwts(:)));
    samples = cell(nsamp,1);
    types = cell(nsamp,1);
    for iter = 1:nsamp

        % choose the parse
        [M,kindx] = rand_discrete(G.models,wts);

        % choose the type-level resampling
        Q = rand_discrete(G.samples_type{kindx});
        Q = Q.copy();
        Q.I = M.I;

        % resample at the token-level
        types{iter} = Q;
        Q = generate_exemplar(Q.copy(),lib);
        samples{iter} = Q;
    end

end

%
% Transform a score into 1/sigma(i), where
% sigma(i) is the rank (from highest to lowest)
% of the score.
%
% Input
%  scores: [K x 1] vector of scores where higher is better
%
% Output
%  wts: [K x 1] new scores (not in log space)
%
function wts = rescore_by_rank(scores)

    assert(isvector(scores));
    K = numel(scores);

    [~,rank_indx] = sort(scores(:),1,'descend');

    wts = zeros(K,1);
    for i=1:K
        wts(i) = 1 ./ rank_indx(i);
    end
    logwts = log(wts);
    wts = exp(logwts - logsumexp(logwts(:)));

    % Assure the precise normalization of the weights
    sum_wts = sum(wts);
    assert(aeq(sum_wts,1));
    wts = wts ./ sum(wts);

end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
% --------------------------------------------------------
% MDP Tracking
% Copyright (c) 2015 CVGL Stanford
% Licensed under The MIT License [see LICENSE for details]
% Written by Yu Xiang
% --------------------------------------------------------
%
% draw dres structure
function show_dres(frame_id, I, tit, dres, state, cmap)

if nargin < 5
    state = 1;
end

if nargin < 6
    cmap = colormap;
end

imshow(I);
title(tit);
hold on;

if isempty(dres) == 1
    index = [];
else
    if isfield(dres, 'state') == 1
        index = find(dres.fr == frame_id & dres.state == state);
    else
        index = find(dres.fr == frame_id);
    end
end

ids = unique(dres.id);
for i = 1:numel(index)
    ind = index(i);

    x = dres.x(ind);
    y = dres.y(ind);
    w = dres.w(ind);
    h = dres.h(ind);
    r = dres.r(ind);

    if isfield(dres, 'id') && dres.id(ind) > 0
        id = dres.id(ind);
        id_index = find(id == ids);
        str = sprintf('%d', id_index);
        index_color = min(1 + floor((id_index-1) * size(cmap,1) / numel(ids)), size(cmap,1));
        c = cmap(index_color,:);
    else
        c = 'g';
        str = sprintf('%.2f', r);
    end
    if isfield(dres, 'occluded') && dres.occluded(ind) > 0
        s = '--';
    else
        s = '-';
    end
    rectangle('Position', [x y w h], 'EdgeColor', c, 'LineWidth', 4, 'LineStyle', s);
    text(x, y-size(I,1)*0.01, str, 'BackgroundColor', [.7 .9 .7], 'FontSize', 14);

    if isfield(dres, 'id') && dres.id(ind) > 0
        % show the previous path
        ind = find(dres.id == id & dres.fr <= frame_id);
        centers = [dres.x(ind)+dres.w(ind)/2, dres.y(ind)+dres.h(ind)];
        patchline(centers(:,1), centers(:,2), 'LineWidth', 4, 'edgecolor', c, 'edgealpha', 0.3);
    end
end
hold off;

author: yuxng


text:
function [ links, params, requesturl, box] = opensearch( q )
%OPENSEARCH Queries GI-cat using the OpenSearch-geo interface
% Returns a list of available links from the matching resources
%
% Usage:  [ links, params, requesturl ] = opensearch( q )
%
%   q.endpoint = URL to opensearch server
%   q.string_text = text string to match
%   q.bbox= bounding box, vector [lon_min lon_max lat_min lat_max]
%       or string 'lon_min,lat_min,lon_max,lat_max'
%   q.time_start = start time, datenum/datestr or ISO string
%   q.time_end  = ending time, datenum/datestr or ISO string
%
% For full options and more explanation, see the GI-CAT OpenSearch doc
%    http://essi-lab.eu/cgi-bin/twiki/view/GIcat/OpenSearchGuide
%
% Example 1:
% q.endpoint='http://geoport.whoi.edu/gi-cat/services/opensearchgeo';
% q.string_text='sea_water_temperature';
% q.bbox=[-74.0 -66.0 39.0 45.];
% q.time_start='1999-01-20 04:33:45';
% q.time_end='1999-08-01 00:00:00';
% [links,params,request_url]=opensearch(q);   % make the query
% dap=links2dap(links); % find only the OPeNDAP links
%
% Example 2:
% q.endpoint='http://testbedapps.sura.org/gi-cat/services/opensearchgeo';
% q.string_text='sea_water_salinity';
% q.bbox=[-82 -73 36 40];
% q.time_start=[2000 06 20 04 33 45];
% q.time_end=[2006 08 01 00 00 00];
% [links,params,request_url]=opensearch(q);   % make the query
% dap=links2dap(links); % find only the OPeNDAP links
%
% Example 3:
% q.endpoint='http://geoport.whoi.edu/gi-cat/services/opensearchgeo';
% q.bbox=[-180 180 -30 70];  % whole globe from -30 to 70 north
% q.time_start=now_utc-28; % start 28 days ago
% q.time_end=now_utc;    % stop now
% q.string_text='relative_humidity'; %search for datasets containing 'relative_humidity'
% [links,params,request_url]=opensearch(q);   % make the query
% dap=links2dap(links); % find only the OPeNDAP links

if(nargin<1);help opensearch; return;end
if ~isfield(q,'endpoint'); disp('must define endpoint (URL)'); return; end
if ~isfield(q,'string_text'); q.string_text='';end
if ~isfield(q,'bbox')
    q.bbox='';
else
    if length(q.bbox)==4, % matlab style bbox [lon_min lon_max lat_min lat_max]
        q.bbox=sprintf('%g, %g, %g, %g',q.bbox([1 3 2 4]));
    end
end
if ~isfield(q,'time_start')
    q.time_start='';
else
    try
        dn=datenum(q.time_start); % handle datenum/datestr input
        q.time_start=datestr(dn,'yyyy-mm-ddTHH:MM:SSZ');
    catch
    end
end
if ~isfield(q,'time_end')
    q.time_end='';
else
    try
        dn=datenum(q.time_end);
        q.time_end=datestr(dn,'yyyy-mm-ddTHH:MM:SSZ');
    catch
    end
end
if ~isfield(q,'loc'); q.loc='';end
if ~isfield(q,'si'); q.si='1';end

% Default "&ct=&" returns only 10 records so specify default to be very
% large if not supplied.
if ~isfield(q,'ct'); q.ct='1000';end

% In GI-CAT 8.4, default "&rel=&" for rel:{geo:relation} returns
% records outside bounding box, so here we specify the default
% to be 'overlaps' if not supplied.  Options are
% {disjoint|overlaps|contains}
if ~isfield(q,'rel'); q.rel='overlaps';end

params = {
    'si';   q.si;
    'ct';   q.ct;
    'rel';  q.rel;
    'st';   q.string_text;
    'bbox'; q.bbox;
    'loc';  q.loc;
    'ts';   q.time_start;
    'te';   q.time_end;
    'outputFormat';'application/atom+xml';
    'enableCrawling'; 'false';
    };

% Optionally return the request URL used for the search
if nargout>2,
    requesturl=q.endpoint;
    for i=1:2:length(params)
        if (i == 1), separator = '?'; else separator = '&'; end
        param = char(java.net.URLEncoder.encode(params{i}));
        value = char(java.net.URLEncoder.encode(params{i+1}));
        requesturl = [requesturl separator param '=' value];
    end
end

% Execute the query
[results] = urlread(q.endpoint,'get', params);

% find all URL links
[~,~,~,~,links1]=regexp(results,'<gmd:URL>(.*?)</gmd:URL>','match'); % early versions of nciso/gicat
[~,~,~,~,links2]=regexp(results,'<dm:srvOp xmlns:dm="http://floraresearch.eu/sdi/services/7.0/dataModel/schema" name="(.*?)</dm:srvOp>','match'); % modern versions of nciso/gicat
for ind = 1:length(links2)
    link = links2(ind);
    st = regexp(link{1}, 'http');
    links2{ind}{1}(1:st{1}-1) = [];
end

% find all box values (bounding box)
[~,~,~,~,box]=regexp(results,'<box xmlns="http://www.georss.org/georss">(.*?)</box>','match');

if ~isempty(links2)
    links = links2;
else
    links = links1;
end

end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function [f,g,H] = taylorModel(d,f,g,H,T)

p = length(d);

fd3 = 0;
gd2 = zeros(p,1);
Hd = zeros(p);
for t1 = 1:p
    for t2 = 1:p
        for t3 = 1:p
            fd3 = fd3 + T(t1,t2,t3)*d(t1)*d(t2)*d(t3);

            if nargout > 1
                gd2(t3) = gd2(t3) + T(t1,t2,t3)*d(t1)*d(t2);
            end

            if nargout > 2
                Hd(t2,t3) = Hd(t2,t3) + T(t1,t2,t3)*d(t1);
            end
        end

    end
end

f = f + g'*d + (1/2)*d'*H*d + (1/6)*fd3;

if nargout > 1
    g = g + H*d + (1/2)*gd2;
end

if nargout > 2
    H = H + Hd;
end

if any(abs(d) > 1e5)
    % We want the optimizer to stop if the solution is unbounded
    g = zeros(p,1);
end

author: zouchuhang

repo: LayoutNet


text:
% LENGTH - length of memory mapped underlying array
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, Nov. 2008

% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function s = length(a)

    s = size(a,1);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
function sp = sp_vector_curl_transform (scalar_spaces, msh)

  warning ('geopdes:obsolete', 'The function SP_VECTOR_CURL_TRANSFORM is obsolete. Using SP_VECTOR with the transformation as last argument')
  sp = sp_vector (scalar_spaces, msh, 'curl-preserving');

end


text:
function x = spgridcc(levelseq, d)
% SPGRIDCC   Computes the Clenshaw-Curtis sparse grid points in
% each coordinate direction.
%    X = SPGRIDCC(LEVELSEQ, D)  Computes the sparse grid points for
%    the given sequence of index sets LEVELSEQ and dimension D. The
%    coordinate value of dimension i is stored in column i of the
%    matrix X. One row of matrix X represents one grid point.

% Author : Andreas Klimke
% Version: 1.1
% Date   : April 20, 2004

% Change log:
% V1.0   : July 30, 2003
%          Initial version
% V1.1   : April 20, 2004
%          Altered function header and simplified code.

% Get the number of levels
nlevels = size(levelseq,1);

% Compute number of points
totalpoints = spseqdimcc(levelseq,d);

% index contains the index of the resulting array containing all
% subdomains of the level.
index = 1;

x = zeros(totalpoints,d);

for kl = 1:nlevels
	level = levelseq(kl,:);
	for i = 1:d
		% compute the points, scaled to [0,1]
		if level(i) == 0
			% take the interval bounds if the level is zero
			c = 0.5;
		elseif level(i) == 1
			c = [0; 1];
		else
			c = ( ( ( (1:2^(level(i)-1)) *2 ) -1).*2^(-level(i)))';
		end
		% Compute the number of grid points per dimension, store it
		% in repvec.
		repvec = ones(1,d);
		for k = 1:d
			if level(k) == 0
				% repvec(k) = 1;
			elseif level(k) < 3
				repvec(k) = 2;
			else
				repvec(k) = 2^(level(k)-1);
			end
		end
		npoints = prod(repvec);
		repvec(i) = 1;
		c = repmat(shiftdim(c, 1-i), repvec);
		x(index:index+npoints-1,i) = c(:);
	end
	index = index + npoints;
end

% -----------------------------------------------------------------
function n = spseqdimcc(levelseq, d)
% SPSEQDIMCC   Compute number of grid points for given sequence of
% index sets.

% Author : Andreas Klimke
% Version: 1.0
% Date   : April 20, 2004

n = 0;
for k = 1:size(levelseq,1);
	ntemp = 1;
	for l = 1:d
		lev = levelseq(k,l);
		if lev == 0
			continue;
		elseif lev <= 2
			ntemp = ntemp * 2;
		else
			ntemp = ntemp * 2^(lev-1);
		end
	end
	n = n + ntemp;
end


text:
function stateVec=orbEls2State(orbEls,deltaT,elType,GM,epsVal)
%%ORBELS2STATE Obtain a velocity and position from a set of orbital
%              elements at a desired time in seconds since the epoch time.
%
%INPUTS: orbEls A 6XnumVec vector of orbital numVec elements to convert.
%               The elements are either Gooding's universal orbital
%               elements or equinoctial orbital elements, as specified by
%               the parameter elType. The components of the elements
%               depends on the type and are explained below.
%        deltaT The time in seconds forward that the trajectory should be
%               propagated from the epoch time of the orbital elements. If
%               this parameter is omitted, a value of 0 is used. For
%               example, if one wishes to propagate the state forward an
%               hour from the time associated with the elements, one can
%               use deltaT=60*60.
%        elType A value indicating the type of orbital elements. Possible
%               values are:
%               0 (The default if omitted) The elements are Gooding's
%                 universal orbital elements.
%               1 The elements are direct equinoctial orbital elements.
%               2 The elements are retrograde equnoctial orbital
%                 elements.
%            GM An optional value of the universal gravitational constant
%               times the mass of the Earth. If omitted, the value
%               Constants.WGS84GMWithAtmosphere is used. The units are
%               m^3/sec^2.
%        epsVal If universal orbital elements are chosen, this is a
%               precision bound used on the alpha term (GM divided by the
%               semi-major axis) for determining whether the trajectory is
%               parabolic. If omitted, a default value of eps is used.
%
%OUTPUTS: stateVec A 6XnumVec matrix of numVec state vectors consisting of
%                 position and velocity in a Cartesian inertial coordinate
%                 system where the gravitating body is at the origin. Units
%                 are meters, and then meters per second.
%
%Gooding's universal orbital elements are consist of (in order)
%alpha=GM/a where a is the semi-major axis in meters
%q=a*(1-e)  where e is the eccentricity (unitless). This is the perifocal
%           distance.
%i          inclination in radians
%Omega      longitude of the ascending node in radians
%omega      argument of periapsis in radians
%tau        the time at pericenter (seconds)
%
%The equinoctial orbital elements (for both direct and retrograde elements)
%are
%a      semi-major axis in meters
%h      first eccentricity vector component (unitless)
%k      second eccentricity vector component (unitless)
%p      first ascending node vector component (unitless)
%q      second ascending node vector component (unitless)
%lambda mean longitude in radians.
%
%The algorithm using Gooding's universal orbital elements is very robust to
%numerical errors. It is implemented using the algorithm of [1], where the
%elements are also described.
%
%Equinoctial orbital elements are discussed in Section 2 of [2] and in [3].
%The diffference between direct and retrograde elements is essentially a
%matter of handedness.
%
%Orbital elements are discussed in general in Chapter 2.2 of [4].
%
%REFERENCES:
%[1] R. H. Gooding, "On universal elements, and conversion procedures to
%    and from position and velocity," Celestial mechanics, vol. 44, no. 3,
%    pp. 283-298, 1988.
%[2] D. A. Danielson, C. P. Sagovac, B. Neta, and L. W. Early,
%    "Semianalytic satellite theory," Mathematics Department, Naval
%    Postgraduate School, Monterey, CA, Tech. Rep., 1995. [Online].
%    Available: http://oai.dtic.mil/oai/oai?verb=getRecord&metadataPrefix= html&identifier=ADA531136
%[3] R. A. Broucke and P. J. Cefola, "On the equinoctial orbit elements,"
%    Celestial Mechanics, vol. 5, no. 3, pp. 303-310, 1972.
%[4] O. Montenbruck and E. Gill, Satellite Orbits: Models, Methods
%    Applications. Berlin: Springer, 2000.
%
%January 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5)
    epsVal=eps;
end

if(nargin<4)
   GM=Constants.WGS84GMWithAtmosphere;
end

if(nargin<3)
    elType=0;
end

if(nargin<2)
    deltaT=0;
end

numVecs=size(orbEls,2);
stateVec=zeros(6,numVecs);

switch(elType)
    case 0
        for curVec=1:numVecs
            stateVec(:,curVec)=orbEls2StateUniv(orbEls(:,curVec),deltaT,GM,epsVal);
        end
    case 1
        for curVec=1:numVecs
            stateVec(:,curVec)=orbEls2StateEquinoctial(orbEls(:,curVec),deltaT,1,GM);
        end
    case 2
        for curVec=1:numVecs
            stateVec(:,curVec)=orbEls2StateEquinoctial(orbEls(:,curVec),deltaT,-1,GM);
        end
    otherwise
        error('Invalid element type provided.')
end
end

function state=orbEls2StateUniv(theEls,deltaT,GM,epsVal)
%%ORBELS2STATEUNIV  An implementation of Gooding's conversion from
%                   universal orbital elements to state vectors. This
%                   implements the ELS3PV function.
%
%The implementation is taken from Appendix B of
%R. H. Gooding, "On universal elements, and conversion procedures to
%and from position and velocity," Celestial mechanics, vol. 44, no. 3,
%pp. 283-298, 1988.
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.

%If we are given classical elements...
alpha=theEls(1);%alpha=GM/a where a is the semi-major axis in meters
q=theEls(2);%q=a*(1-e), where e=eccentricity. q is the perifocal distance
i=theEls(3);%Inclination
Omega=theEls(4);%Longitude of the ascending node
omega=theEls(5);%Argument of periapsis
tau=theEls(6)+deltaT;%The time at pericenter.

[r,u,VR,VT]=ELS2PV(alpha,q,omega,tau,GM,epsVal);

c=cos(u);
s=sin(u);
x1=r*c;
y1=r*s;
x2=VR*c-VT*s;
y2=VR*s+VT*c;

c=cos(i);
s=sin(i);
z=y1*s;
y1=y1*c;
zDot=y2*s;
y2=y2*c;
c=cos(Omega);
s=sin(Omega);
x=x1*c-y1*s;
y=x1*s+y1*c;
xDot=x2*c-y2*s;
yDot=x2*s+y2*c;

state=[x;y;z;xDot;yDot;zDot];
end

function state=orbEls2StateEquinoctial(theEls,deltaT,I,GM)
%%ORBELS2STATEEQUINOCTIAL An algorithm to convert from equinoctial orbital
%                         elements to a target state vector.
%
%The input I is the retrograde factor ond is +1 for direct equinoctial
%elements and -1 for retrograde equinoctial elements.
%
%The algorithm is taken from Section 2.1.4 of
%D. A. Danielson, C. P. Sagovac, B. Neta, and L. W. Early, "Semianalytic
%satellite theory," Mathematics Department, Naval Postgraduate School,
%Monterey, CA, Tech. Rep., 1995. [Online]. Available:
%http://oai.dtic.mil/oai/oai?verb=getRecord&metadataPrefix= html&identifier=ADA531136
%
%The equation for propagating the elements forward in time is taken from
%R. A. Broucke and P. J. Cefola, "On the equinoctial orbit elements,"
%Celestial Mechanics, vol. 5, no. 3, pp. 303-310, 1972.
%
%January 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

a=theEls(1);%Semi-major axis
h=theEls(2);%First eccentricity vector component
k=theEls(3);%Second eccentricity vector component
p=theEls(4);%First ascending node vector component
q=theEls(5);%Second ascending node vector component
lambda=theEls(6);%Mean longitude (at epoch)

n=sqrt(GM/a^3);%The mean motion
%The propagation of the mean longitude is from Equation 15 in Brouke and
%Cefola.
lambda=lambda+n*deltaT;

%Section 2.1.4 conversion from equinoctial elements to position and
%velocity.
coeff=(1/(1+p^2+q^2));
f=coeff*[1-p^2+q^2;
         2*p*q;
         -2*I*p];
g=coeff*[2*I*p*q;
         (1+p^2-q^2)*I;
         2*q];

%Solve the Equinoctial form of Kepler's Equation. to get F.
F=solveEquinoctialKeplersEq(lambda,h,k);

n=sqrt(GM/a^3);
b=1/(1+sqrt(1-h^2-k^2));
denom=1-h*sin(F)-k*cos(F);
sinL=((1-k^2*b)*sin(F)+h*k*b*cos(F)-h)/denom;
cosL=((1-h^2*b)*cos(F)+h*k*b*sin(F)-k)/denom;

r=a*(1-h^2-k^2)/(1+h*sinL+k*cosL);
X=r*cosL;
Y=r*sinL;
XDot=-n*a*(h+sinL)/sqrt(1-h^2-k^2);
YDot=n*a*(k+cosL)/sqrt(1-h^2-k^2);

r=X*f+Y*g;
rDot=XDot*f+YDot*g;

state=[r;rDot];

end

function [r,u,VR,VT]=ELS2PV(alpha,q,omega,tau,GM,epsVal)
%%ELS2PV An implementation of the ELS2PV subroutine of Gooding's algorithm
%        for universal orbital element conversion. This is a
%        two-dimensional conversion suroutine.
%
%The algorithm is taken from Appendix A of
%R. H. Gooding, "On universal elements, and conversion procedures to
%and from position and velocity," Celestial mechanics, vol. 44, no. 3,
%pp. 283-298, 1988.
%with minor changes so that alpha is compared to an epsilon value rather
%than zero for determining the parabolic case.
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.

if(abs(alpha)<=epsVal) %Parabolic case, GM cannot be zero.
    d=solveCubicEqGoodin(0.5/GM,q,1.5*GM*tau);
    r=q+(1/2)*d^2/GM;
    h=sqrt(2*GM*q);
    v=2*atan2(d,h);
else%Ellipse or hyperbola
    e1=alpha*q;
    e=GM-e1;
    ep1=GM+e;
    h=sqrt(q*ep1);

    em=tau*abs(alpha)^(3/2);
    if(alpha>0)%ellipse, GM cannot be zero
        ee2=(1/2)*solveKeplersEq(em/GM,e1/GM,2);
        s2=sin(ee2);
        c2=cos(ee2);
        r=q+2*e*s2^2/alpha;
        d=2*e*s2*c2/sqrt(abs(alpha));
        v=2*atan2(ep1*s2,h*sqrt(abs(alpha))*c2);
        emv=em/GM-v;
        v=v+4*pi*fix(abs(emv/(4*pi))+(1/2))*sign(emv);
    else%hyperbola
        s=solveKeplersEq(em/e,-e1/e,3);
        c=sqrt(1+s^2);
        s2=s^2/(c+1);
        r=q-e*s2/alpha;
        d=e*s/sqrt(abs(alpha));
        v=atan2(s*h*sqrt(abs(alpha)),-GM*s2-e1);
    end
end

%All orbits
u=omega+v;
VR=d/r;
VT=h/r;
end

function x=solveCubicEqGoodin(a,b,c)
%SOLVECUBICEQGOODING  This solves the equation a*x^3+3*b*x-2*c=0, where
%                     a>=0 and b^3+a*c^2>=0 for the real root of x.
%                     Additionally, if a and b are both zero, then zero is
%                     generated in lieu of an indeterminate solution.
%
%INPUTS: a,b,c The coefficients in the equation a*x^3+3*b*x-2*c=0 where
%              a>=0 and b^3+a*c^2>=0.
%
%OUTPUTS: x The one real solution to a*x^3+3*b*x-2*c=0.
%
%In the report
%R. H. Gooding and A. W. Odell, "The hyperbolic Kepler's equation,
%and the elliptic equations revisited," Royal Aerospace Executive,
%Procurement Executive, Ministry of Defence, Farnborough, Hants, United
%Kingdom, Tech. Rep. 369, Jul. 1989.
%the source of the algorithm is listed as being in
%R. H. Gooding, "Solution of the hyperbolic Kepler's equation," Royal
%Aerospace Executive, Procurement Executive, Ministry of Defence,
%Farnborough, Hants, United Kingdom, Tech. Rep. 87042, 1987.
%However, the 1987 report appears to be unobtainable. Nonetheless, the
%solution is given in Equation 4.47 in
%G. Beutler, Methods of Celestial Mechanics: Physical, Mathematical and
%Numerical Principles. Berlin: Springer, 2005, vol. 1.
%The modification to return zero instead of infinity is consistent with how
%Gooding used the algorithm.
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.

if(a==0&&b==0||c==0)
    x=0;
else
    d=nthroot(sqrt(a)*abs(c)+sqrt(b^3+a*c^2),3)^2;
    x=2*c/(d+b+b^2/d);
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 area = sphere_cap_area_3d ( r, h )

%*****************************************************************************80
%
%% SPHERE_CAP_AREA_3D computes the surface area of a spherical cap in 3D.
%
%  Discussion:
%
%    Draw any radius of the sphere and note the point P where the radius
%    intersects the sphere.  Consider the point on the radius line which is
%    H units from P.  Draw the circle that lies in the plane perpendicular to
%    the radius, and which intersects the sphere.  The circle divides the sphere
%    into two pieces, and the corresponding disk divides the solid sphere into
%    two pieces.  The spherical cap is the part of the solid sphere that
%    includes the point P.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the sphere.
%
%    Input, real H, the "height" of the spherical cap.
%    H must be between 0 and 2 * R.
%
%    Output, real AREA, the area of the spherical cap.
%
  if ( h <= 0.0 )
    area = 0.0;
  elseif ( 2.0 * r <= h )
    area = 4.0 * pi * r * r;
  else
    area = 2.0 * pi * r * h;
  end

  return
end


text:
function [x,y] = km_kde(data,sigma,nsteps,range_x)
% KM_KDE performs kernel density estimation (KDE) on one-dimensional data
% http://en.wikipedia.org/wiki/Kernel_density_estimation
%
% Input:	- data: input data, one-dimensional
%           - sigma: bandwidth (sometimes called "h")
%           - nsteps: optional number of abscis points. If nsteps is an
%             array, the abscis points will be taken directly from it.
%           - range_x: optional factor for abscis expansion beyond extrema
% Output:	- x: equispaced abscis points
%			- y: estimates of p(x)
% USAGE: [xx,pp] = km_kde(data,sigma,nsteps,range_x)
%
% Author: Steven Van Vaerenbergh (steven *at* gtas.dicom.unican.es), 2010.
%
% This file is part of the Kernel Methods Toolbox for MATLAB.
% https://github.com/steven2358/kmbox

N = length(data);	% number of data points

if (nargin<4), range_x = 1; end % default abscis expansion
if (nargin<3), nsteps = 100; end % default number of abscis points

% obtain full data range + extra bit
mm(1) = min(data);
mm(2) = max(data);
if length(nsteps) > 1
    x = nsteps;
else
    xm = range_x*mm;
    if (mm(1)>0), xm(1)=mm(1)/range_x; end
    if (mm(2)<0), xm(2)=mm(2)/range_x; end
    x = linspace(xm(1),xm(2),nsteps);
end
y = zeros(size(x));

% kernel density estimation
c = 1/sqrt(2*pi*sigma^2);
for i=1:N,
	y = y + 1/N*c*exp(-(data(i)-x).^2/(2*sigma^2));
end

author: epfl-lasa


text:
addpath ./propack
addpath ./yima_rpca
addpath ./lib

X = rand(1000, 1000);
k = 20;

tic;[U1 s1 V1] = svdex(X, k, struct('svd_solver', 'svd'));t1 = toc;
tic;[U2 s2 V2] = svdex(X, k, struct('svd_solver', 'svds'));t2 = toc;
tic;[U3 s3 V3] = svdex(X, k, struct('svd_solver', 'lansvd', 'init', randn(size(X,1),1)));t3 = toc;
figure; plot(1:k,s1,'r', 1:k,s2,'g^', 1:k,s3,'bv');
fprintf('Time:\n     svd=%fs\n    svds=%fs\n  lansvd=%fs\n',t1,t2,t3);

Y1 = bsxfun(@times, U1, s1')*V1';
err1 = max(abs(X(:) - Y1(:)));
Y2 = bsxfun(@times, U2, s2')*V2';
err2 = max(abs(X(:) - Y2(:)));
Y3 = bsxfun(@times, U3, s3')*V3';
err3 = max(abs(X(:) - Y3(:)));
fprintf('Error:\n     svd=%f\n    svds=%f\n  lansvd=%f\n',err1,err2,err3);

author: andrewssobral

repo: lrslibrary


text:
function SaveBinary(filename,data,varargin)

%SaveBinary - Save data in a multiplexed binary file.
%
%  USAGE
%
%    SaveBinary(filename,data,<options>)
%
%    filename       file to save to
%    data           data to save
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'precision'   sample precision (default = 'int16')
%     'mode'        either 'new' to create a new file (default) or 'append'
%                   to add to an existing file
%    =========================================================================

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% 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.

% Default values
% precision = 'int16'; % THIS SHOULD NEVER BE ASSUMED
precision = class(data);
mode = 'new';

if nargin < 2 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help SaveBinary'' for details).');
end

% Parse options
for i = 1:2:length(varargin),
  if ~ischar(varargin{i}),
    error(['Parameter ' num2str(i+3) ' is not a property (type ''help SaveBinary'' for details).']);
  end
  switch(lower(varargin{i})),
    case 'precision',
      precision = varargin{i+1};
      if ~isstring_FMAT(precision),
        error('Incorrect value for property ''precision'' (type ''help SaveBinary'' for details).');
      end
    case 'mode',
      mode = lower(varargin{i+1});
      if ~isstring_FMAT(mode,'new','append'),
        error('Incorrect value for property ''mode'' (type ''help SaveBinary'' for details).');
      end
    otherwise,
      error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help SaveBinary'' for details).']);
  end
end

sizeInBytes = 0;
switch precision,
	case {'uchar','unsigned char','schar','signed char','int8','integer*1','uint8','integer*1'},
		sizeInBytes = 1;
	case {'int16','integer*2','uint16','integer*2'},
		sizeInBytes = 2;
	case {'int32','integer*4','uint32','integer*4','single','real*4','float32','real*4'},
		sizeInBytes = 4;
	case {'int64','integer*8','uint64','integer*8','double','real*8','float64','real*8'},
		sizeInBytes = 8;
end

if strcmp(mode,'new'),
	f = fopen(filename,'w');
	status = fseek(f,0,'bof');
else
	f = fopen(filename,'a');
	status = fseek(f,0,'eof');
end
if status ~= 0,
	fclose(f);
	error('Could not start writing (possible reasons include insufficient access permissions).');
end

fwrite(f,data',precision);

fclose(f);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function mean = beta_mean ( a, b )

%*****************************************************************************80
%
%% BETA_MEAN returns the mean of the Beta PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    0.0D+00 < A,
%    0.0D+00 < B.
%
%    Output, real MEAN, the mean of the PDF.
%
  mean = a / ( a + b );

  return
end

author: johannesgerer


text:
% Codes for CVPR-15 work `Face Alignment by Coarse-to-Fine Shape Searching'
% Any question please contact Shizhan Zhu: zhshzhutah2@gmail.com
% Released on July 25, 2015

function [featPCA,pca_model,cut_id] = getPCA(featOri,cut_id)
%[featPCA,pca_model,cut_id] = getPCA(featOri,cut_id)
% pca_model: coeff(n*keep) and meanFeatOri

% fprintf('0');
pca_model.meanFeatOri = double(mean(featOri,1));
assert(size(featOri,1)>=size(featOri,2));
%-----------------------------------------------
if size(featOri,2)>5000
    featOri = single(featOri);
end
[n,p] = size(featOri);
featOri = bsxfun(@minus,featOri,mean(featOri,1));
% parfor i = 1:size(featOri,1),featOri(i,:) = featOri(i,:) - pca_model.meanFeatOri; end;
% fprintf('1');
[U,sigma,coeff] = svd(featOri,0);
clear featOri
% fprintf('2');
sigma = diag(sigma);
% fprintf('3');
score = bsxfun(@times,U,sigma');
clear U;
% fprintf('4');
sigma = sigma ./ sqrt(n-1);
latent = sigma.^2;
clear sigma;
[~,maxind] = max(abs(coeff),[],1);
d = size(coeff,2);
% fprintf('5');
colsign = sign(coeff(maxind + (0:p:(d-1)*p)));
clear maxind;
coeff = bsxfun(@times,coeff,colsign);
score = bsxfun(@times,score,colsign);
clear colsign;
% [coeff,score,latent] = princomp(featOri);

if nargin<2 || isempty(cut_id)
    cut_id = find(cumsum(latent)>sum(latent)*0.98,1);
end

pca_model.coeff = double(coeff(:,1:cut_id));
featPCA = double(score(:,1:cut_id));

end

author: zhusz


text:
function kappa = curvature(alpha)
% compute curvature tensor from a (complete) dislocation density tensor
%
% Input
%  alpha  - @dislocationDensityTensor
%
% Output
%  kappa - @curvatureTensor

kappa = curvatureTensor( alpha' - 0.5 * diag(tensor(trace(alpha),'rank',0)));

author: mtex-toolbox

repo: mtex

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


text:
% MIPS
%   Version 1.5         12-Dec-2022
%
%   The MATPOWER Interior Point Solver (MIPS) is a package of MATLAB/Octave
%   M-files for solving non-linear programming problems (NLPs) using a
%   primal dual interior point method. MIPS is based on code written in C
%   language by Hongye Wang as a graduate student at Cornell University
%   for optimal power flow applications. It was later ported to the
%   MATLAB/Octave language by Ray D. Zimmerman for use in MATPOWER.
%
%   The latest version of MIPS can be found on GitHub at:
%
%       https://github.com/MATPOWER/mips
%
%   MIPS is covered by the 3-clause BSD License (see LICENSE for details).

%   MIPS
%   Copyright (c) 2010-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MIPS.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mips for more info.

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%% POWER_WRITE writes a power table to a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POWER_WRITE\n' );
  fprintf ( 1, '  Write a table of integer powers to a file.\n' );

  file_name = 'power_write.txt';

  file_unit = fopen ( file_name, 'w' );

  fprintf ( file_unit, '\n' );
  fprintf ( file_unit, '           N   N-squared     N-cubed\n' );
  fprintf ( file_unit, '\n' );
  for i = 0 : 10
    fprintf ( file_unit, '  %10d  %10d  %10d\n', i, i^2, i^3 );
  end

  fclose ( file_unit );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  POWER_WRITE wrote the file %s\n', file_name );

  exit

repo: jburkardt-m