text:
function seg=bbxflatsegment(node,loop)
%
% seg=bbxflatsegment(node,loop)
%
% decompose edge loops into flat segments along the x/y/z
% planes of the bounding box
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2008/04/08
%
% input:
%    node:  x,y,z coordinates of each node of the mesh
%    loop:  input, a single vector separated by NaN, each segment
%             is a close-polygon consisted by node IDs
% output:
%    seg:   output, a single vector separated by NaN, each segment
%             is a close-polygon on x/y/z plane
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

pos=node(loop,:);

% get the bounding box
mi=min(pos);
ma=max(pos);

% extract nodes on the bounding box
idx0=find(abs(pos(:,1)-mi(1))<1e-6)';
idx1=find(abs(pos(:,1)-ma(1))<1e-6)';

idy0=find(abs(pos(:,2)-mi(2))<1e-6)';
idy1=find(abs(pos(:,2)-ma(2))<1e-6)';

idz0=find(abs(pos(:,3)-mi(3))<1e-6)';
idz1=find(abs(pos(:,3)-ma(3))<1e-6)';

% need to be more than 3 points to make a flat polygon

if(length(idx0)<=3) idx0=[]; end
if(length(idx1)<=3) idx1=[]; end
if(length(idy0)<=3) idy0=[]; end
if(length(idy1)<=3) idy1=[]; end
if(length(idz0)<=3) idz0=[]; end
if(length(idz1)<=3) idz1=[]; end

nn=length(loop);

% if the original is a flat polygon, return

if(unique(length(idx0))==nn || unique(length(idx1))==nn ...
  |unique(length(idy0))==nn || unique(length(idy1))==nn ...
  |unique(length(idz0))==nn || unique(length(idz1))==nn)
    seg=loop(:)';
    return;
end

% otherwise, find the combination that split the loop

if(length(unique([idx0 idy0 idz0]))==nn)
   seg= [loop(idx0),nan,loop(idy0),nan,loop(idz0)];
elseif(length(unique([idx0 idy1 idz0]))==nn)
   seg= [loop(idx0),nan,loop(idy1),nan,loop(idz0)];
elseif(length(unique([idx0 idy0 idz1]))==nn)
   seg= [loop(idx0),nan,loop(idy0),nan,loop(idz1)];
elseif(length(unique([idx0 idy1 idz1]))==nn)
   seg= [loop(idx0),nan,loop(idy1),nan,loop(idz1)];
elseif(length(unique([idx1 idy0 idz0]))==nn)
   seg= [loop(idx1),nan,loop(idy0),nan,loop(idz0)];
elseif(length(unique([idx1 idy1 idz0]))==nn)
   seg= [loop(idx1),nan,loop(idy1),nan,loop(idz0)];
elseif(length(unique([idx1 idy0 idz1]))==nn)
   seg= [loop(idx1),nan,loop(idy0),nan,loop(idz1)];
elseif(length(unique([idx1 idy1 idz1]))==nn)
   seg= [loop(idx1),nan,loop(idy1),nan,loop(idz1)];
else
    seg=[];
end

% remove pattern [ ... nan nan ...] in the result

if(length(seg) && any(isnan(seg)))
    id=regexp(sprintf('%d',isnan(seg)),'11');
    if(length(id))
        seg(id+1)=[];
    end
end

author: fangq


text:
%% This script takes a directory of files that have been processed by PREP
% and produces reports.

%% Read in the file and set the necessary parameters
dataDir = 'F:\TempData';
summaryFolder = 'F:\TempDataReports';
publishOn = true;

%% Get the directory list
inList = dir(dataDir);
inNames = {inList(:).name};
inTypes = [inList(:).isdir];
inNames = inNames(~inTypes);

%% Setup up the names
basename = 'vep';
summaryReportName = [basename '_summary.html'];
sessionFolder = '.';
summaryFileName = [summaryFolder filesep summaryReportName];
if exist(summaryFileName, 'file')
   delete(summaryFileName);
end

%% Publish the reports
for k = 1:length(inNames)
    [~, theName, theExt] = fileparts(inNames{k});
    if ~strcmpi(theExt, '.set') && ~strcmpi(theExt, '.mat')
        continue;
    end
    sessionReportName = [theName '.pdf'];
    fname = [dataDir filesep inNames{k}];
    load(fname, '-mat');
    sessionFileName = [summaryFolder filesep sessionReportName];
    consoleFID = 1;
    publishPrepReport(EEG, summaryFileName, sessionFileName, consoleFID, publishOn);
end

author: VisLab

repo: EEG-Clean-Tools


text:
% DFTFILT3 - discrete complex wavelet filters
%
% Usage:
%   >> [wavelet,cycles,freqresol,timeresol] = dftfilt3( freqs, cycles, srate, varargin)
%
% Inputs:
%   freqs    - vector of frequencies of interest.
%   cycles   - cycles array. If cycles=0, then the Hanning tapered Short-term FFT is used.
%              If one value is given and cycles>0, all wavelets have
%              the same number of cycles. If two values are given, the
%              two values are used for the number of cycles at the lowest
%              frequency and at the highest frequency, with linear or
%              log-linear interpolation between these values for intermediate
%              frequencies
%   srate    - sampling rate (in Hz)
%
% Optional Inputs: Input these as 'key/value pairs.
%   'cycleinc' - ['linear'|'log'] increase mode if [min max] cycles is
%              provided in 'cycle' parameter. {default: 'linear'}
%   'winsize'  Use this option for Hanning tapered FFT or if you prefer to set the length of the
%              wavelets to be equal for all of them (e.g., to set the
%              length to 256 samples input: 'winsize',256). {default: [])
%              Note: the output 'wavelet' will be a matrix and it may be
%              incompatible with current versions of timefreq and newtimef.
%   'timesupport' The number of temporal standard deviation used for wavelet lengths {default: 7)
%
% Output:
%   wavelet - cell array or matrix of wavelet filters
%   timeresol - temporal resolution of Morlet wavelets.
%   freqresol - frequency resolution of Morlet wavelets.
%
% Note: The length of the window is always made odd.
%
% Authors: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 3/28/2003
%          Rey Ramirez, SCCN/INC/UCSD, La Jolla, 9/26/2006

% Copyright (C) 3/28/2003 Arnaud Delorme 8, SCCN/INC/UCSD, arno@salk.edu
%
% 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.

%
% Revision 1.12 2006/09/25  rey r
% Almost complete rewriting of dftfilt2.m, changing both Morlet and Hanning
% DFT to be more in line with conventional implementations.
%
% Revision 1.11  2006/09/07 19:05:34  scott
% further clarified the Morlet/Hanning distinction -sm
%
% Revision 1.10  2006/09/07 18:55:15  scott
% clarified window types in help msg -sm
%
% Revision 1.9  2006/05/05 16:17:36  arno
% implementing cycle array
%
% Revision 1.8  2004/03/04 19:31:03  arno
% email
%
% Revision 1.7  2004/02/25 01:45:55  arno
% sinus test
%
% Revision 1.6  2004/02/15 22:23:08  arno
% implementing morlet wavelet
%
% Revision 1.5  2003/05/09 20:55:10  arno
% adding hanning function
%
% Revision 1.4  2003/04/29 16:02:54  arno
% header typos
%
% Revision 1.3  2003/04/29 01:09:16  arno
% debug imaginary part
%
% Revision 1.2  2003/04/28 23:01:13  arno
% *** empty log message ***
%
% Revision 1.1  2003/04/28 22:46:49  arno
% Initial revision
%

function [wavelet,cycles,freqresol,timeresol] = dftfilt3( freqs, cycles, srate, varargin);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rey fixed all input parameter sorting.
if nargin < 3
    error(' A minimum of 3 arguments is required');
end
numargin=length(varargin);
if rem(numargin,2)
    error('There is an uneven number key/value inputs. You are probably missing a keyword or its value.')
end
varargin(1:2:end)=lower(varargin(1:2:end));

% Setting default parameter values.
cycleinc='linear';
winsize=[];
timesupport=7;  % Setting default of 7 temporal standard deviations for wavelet's length.

for n=1:2:numargin
    keyword=varargin{n};
    if strcmpi('cycleinc',keyword)
        cycleinc=varargin{n+1};
    elseif strcmpi('winsize',keyword)
        winsize=varargin{n+1};
        if ~mod(winsize,2)
            winsize=winsize+1; % Always set to odd length wavelets and hanning windows;
        end
    elseif strcmpi('timesupport',keyword)
        timesupport=varargin{n+1};
    else
        error(['What is ' keyword '? The only legal keywords are: type, cycleinc, winsize, or timesupport.'])
    end
end
if isempty(winsize) && cycles(1)==0
    error('If you are using a Hanning tapered FFT, please supply the winsize input-pair.')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute number of cycles at each frequency
% ------------------------------------------
type='morlet';
if length(cycles) == 1 && cycles(1)~=0
    cycles = cycles*ones(size(freqs));
elseif length(cycles) == 2
    if strcmpi(cycleinc, 'log') % cycleinc
         cycles = linspace(log(cycles(1)), log(cycles(2)), length(freqs));
         cycles = exp(cycles);
         %cycles=logspace(log10(cycles(1)),log10(cycles(2)),length(freqs)); %rey
    else
        cycles = linspace(cycles(1), cycles(2), length(freqs));
    end
end
if cycles==0
    type='sinus';
end

sp=1/srate; % Rey added this line (i.e., sampling period).
% compute wavelet
for index = 1:length(freqs)
    fk=freqs(index);
    if strcmpi(type, 'morlet') % Morlet.
        fk=fk/srate; % Normalize frequency for textbook equations as in TB97
        sigf=fk/cycles(index); % Computing time and frequency standard deviations, resolutions, and normalization constant.
        sigt=1./(2*pi*sigf);
        A=1./sqrt(sigt*sqrt(pi));
        timeresol(index)=2*sigt/srate; % sec
        freqresol(index)=2*sigf*srate; % Hz
        if isempty(winsize) % bases will be a cell array.
%             tneg=[-sp:-sp:-sigt*timesupport/2];
%             tpos=[0:sp:sigt*timesupport/2];
%             t=[fliplr(tneg) tpos];
            t = (0:floor(sigt*timesupport/2)*2)-floor(sigt*timesupport/2); % Always odd; backward compatible
            psi=A.*(exp(-(t.^2)./(2*(sigt^2))).*exp(2*i*pi*fk*t));
            wavelet{index}=psi;  % These are the wavelets with variable number of samples based on temporal standard deviations (sigt).
        else % bases will be a matrix.
%             tneg=[-sp:-sp:-sp*winsize/2];
%             tpos=[0:sp:sp*winsize/2];
%             t=[fliplr(tneg) tpos];
            t = (0:floor(winsize/2)*2)-floor(winsize/2); % Always odd; backward compatible
            psi=A.*(exp(-(t.^2)./(2*(sigt^2))).*exp(2*i*pi*fk*t));
            wavelet(index,:)=psi; % These are the wavelets with the same length.
            % This is useful for doing time-frequency analysis as a matrix vector or matrix matrix multiplication.
        end
    elseif strcmpi(type, 'sinus') % Hanning
        tneg=[-sp:-sp:-sp*winsize/2];
        tpos=[0:sp:sp*winsize/2];
        t=[fliplr(tneg) tpos];
        win = exp(2*i*pi*fk*t);
        wavelet(index,:) = win .* hanning(winsize)';
        %wavelet{index} = win .* hanning(winsize)';
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end
end

% symmetric hanning function
function w = hanning(n)
if ~rem(n,2)
    w = .5*(1 - cos(2*pi*(1:n/2)'/(n+1)));
    w = [w; w(end:-1:1)];
else
    w = .5*(1 - cos(2*pi*(1:(n+1)/2)'/(n+1)));
    w = [w; w(end-1:-1:1)];
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function ParforProgressStressTest2(N)
% Stress test for 'ParforProgressStarter2'. In case of timeouts, you can
% use this function to determine how many simultaneous connections your
% computer can handle, and adjust the ppm.increment() call accordingly.
%
% Copyright (c) 2010-2012, Andreas Kotowicz
%
%%
    if nargin < 1
        N = 10000;
    end

    do_debug = 1;

    %% initialize ParforProgress monitor

    try % Initialization
        ppm = ParforProgressStarter2('test task - ParforProgressStarter2', N, 0.1, do_debug);
    catch me % make sure "ParforProgressStarter2" didn't get moved to a different directory
        if strcmp(me.message, 'Undefined function or method ''ParforProgressStarter2'' for input arguments of type ''char''.')
            error('ParforProgressStarter2 not in path.');
        else
            % this should NEVER EVER happen.
            msg{1} = 'Unknown error while initializing "ParforProgressStarter2":';
            msg{2} = me.message;
            print_error_red(msg);
            % backup solution so that we can still continue.
            ppm.increment = nan(1, N);
        end
    end

    %% execute dummy loop - replace 'rand(1, 10000)' with your function.

    t0 = tic();

    parfor i = 1 : N
        rand(1, 10000);
        ppm.increment(i); %#ok<PFBNS>
    end

    total_time = toc(t0);

    %% clean up
    try % use try / catch here, since delete(struct) will raise an error.
        delete(ppm);
    catch me %#ok<NASGU>
    end

    %% show runtime

    disp(['ParforProgressStressTest2 running time: ' num2str(total_time) 's.']);

end
%% EOF

author: thomaskuestner


text:
% This file is part of libDAI - http://www.libdai.org/
%
% Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
%
% Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.

% This example program illustrates how to construct a factorgraph
% by means of the sprinkler network example discussed at
% http://www.cs.ubc.ca/~murphyk/Bayes/bnintro.html
% using the SWIG octave wrapper of libDAI

dai

C = dai.Var(0, 2);  % Define binary variable Cloudy (with label 0)
S = dai.Var(1, 2);  % Define binary variable Sprinkler (with label 1)
R = dai.Var(2, 2);  % Define binary variable Rain (with label 2)
W = dai.Var(3, 2);  % Define binary variable Wetgrass (with label 3)

% Define probability distribution for C
P_C = dai.Factor(C);
P_C(0) = 0.5;            % C = 0
P_C(1) = 0.5;            % C = 1

% Define conditional probability of S given C
P_S_given_C = dai.Factor(dai.VarSet(S,C));
P_S_given_C(0) = 0.5;    % C = 0, S = 0
P_S_given_C(1) = 0.9;    % C = 1, S = 0
P_S_given_C(2) = 0.5;    % C = 0, S = 1
P_S_given_C(3) = 0.1;    % C = 1, S = 1

% Define conditional probability of R given C
P_R_given_C = dai.Factor(dai.VarSet(R,C));
P_R_given_C(0) = 0.8;    % C = 0, R = 0
P_R_given_C(1) = 0.2;    % C = 1, R = 0
P_R_given_C(2) = 0.2;    % C = 0, R = 1
P_R_given_C(3) = 0.8;    % C = 1, R = 1

% Define conditional probability of W given S and R
SRW = dai.VarSet(S,R);
SRW.append(W);
P_W_given_S_R = dai.Factor(SRW);
P_W_given_S_R(0) = 1.0;  % S = 0, R = 0, W = 0
P_W_given_S_R(1) = 0.1;  % S = 1, R = 0, W = 0
P_W_given_S_R(2) = 0.1;  % S = 0, R = 1, W = 0
P_W_given_S_R(3) = 0.01; % S = 1, R = 1, W = 0
P_W_given_S_R(4) = 0.0;  % S = 0, R = 0, W = 1
P_W_given_S_R(5) = 0.9;  % S = 1, R = 0, W = 1
P_W_given_S_R(6) = 0.9;  % S = 0, R = 1, W = 1
P_W_given_S_R(7) = 0.99; % S = 1, R = 1, W = 1

% Build factor graph consisting of those four factors
SprinklerFactors = dai.VecFactor();
SprinklerFactors.append(P_C);
SprinklerFactors.append(P_R_given_C);
SprinklerFactors.append(P_S_given_C);
SprinklerFactors.append(P_W_given_S_R);
SprinklerNetwork = dai.FactorGraph(SprinklerFactors);

% Write factorgraph to a file
SprinklerNetwork.WriteToFile('sprinkler.fg');
fprintf('Sprinkler network written to sprinkler.fg\n');

% Output some information about the factorgraph
fprintf('%d variables\n', SprinklerNetwork.nrVars());
fprintf('%d factors\n', SprinklerNetwork.nrFactors());

% Calculate joint probability of all four variables
P = dai.Factor();
for I = 0:(SprinklerNetwork.nrFactors()-1)
    P *= SprinklerNetwork.factor(I);
end
P.normalize(); % Not necessary: a Bayesian network is already normalized by definition

% Calculate some probabilities
denom = P.marginal(dai.VarSet(W))(1);
fprintf('P(W=1) = %f\n', denom);
fprintf('P(S=1 | W=1) = %f\n', P.marginal(dai.VarSet(S,W))(3) / denom);
fprintf('P(R=1 | W=1) = %f\n', P.marginal(dai.VarSet(R,W))(3) / denom);

author: anhncs

repo: Probabilistic-Graphical-Models


text:
% This function is designed to determine whether a sphere collides with a
% give triangle

function [haveCollided] = CheckSphereTriangleIntersection(p,r,pa,pb,pc)

% INpUT pcHEpcK
assert(numel(p) == 3,'pcentroid must be a column vector [3x1]');
assert(numel(r) == 1,'Radius must be a scalar');
assert(numel(pa) == 3 && numel(pb) == 3 && numel(pc) == 3,...
       'Vertex point must be a column vector [3x1]');

% pcHEpcK pLpaNpaR SEppaRpaTION
pa = pa - p;
pb = pb - p;
pc = pc - p;
rr = r^2;
V = cross(pb-pa,pc-pa);
d = dot(pa,V);
e = dot(V,V);
sep1 = d^2 > rr*e;
% pcHEpcK VERTIpcES
aa = dot(pa,pa);
ab = dot(pa,pb);
ac = dot(pa,pc);
bb = dot(pb,pb);
bc = dot(pb,pc);
cc = dot(pc,pc);
sep2 = (aa > rr) & (ab > aa) & (ac > aa);
sep3 = (bb > rr) & (ab > bb) & (bc > bb);
sep4 = (cc > rr) & (ac > cc) & (bc > cc);
% pcHEpcK TRIpaNGLE EDGES
papb = pb - pa;
pbpc = pc - pb;
pcpa = pa - pc;
d1 = ab - aa;
d2 = bc - bb;
d3 = ac - cc;
e1 = dot(papb,papb);
e2 = dot(pbpc,pbpc);
e3 = dot(pcpa,pcpa);
Q1 = pa*e1 - d1*papb;
Q2 = pb*e2 - d2*pbpc;
Q3 = pc*e3 - d3*pcpa;
Qpc = pc*e1 - Q1;
Qpa = pa*e2 - Q2;
Qpb = pb*e3 - Q3;
sep5 = dot(Q1,Q1) > (rr*e1*e1) & dot(Q1,Qpc) > 0;
sep6 = dot(Q2,Q2) > (rr*e2*e2) & dot(Q2,Qpa) > 0;
sep7 = dot(Q3,Q3) > (rr*e3*e3) & dot(Q3,Qpb) > 0;
isSeparated = sep1 || sep2 || sep3 || sep4 || sep5 || sep6 || sep7;
haveCollided = ~isSeparated;
end

repo: OpenMAS


text:
% RunMCMCSimForBPHMM
%  Generic harness for running many iterations of MCMC,
%   allows sensible reporting/saving of samples and diagnostics
%USAGE
%  Usually called from more "user-friendly" function "runBPHMM"
%  but if specific data and initial configuration Psi are available,
%  >> RunMCMCForBPHMM( data, Psi, algP, outP )
%INPUT
%  data : SeqData object, defining a collection of sequences to fit model
%  Psi  : initial model state
%           usually generated by a function in the "init" folder
%  algParams : specifies MCMC behavior (# iterations, proposal distribs)
%           see defaults/defaultMCMCParams_BPHMM.m for details
%  outParams : specifies MCMC output behavior
%           how often to save samples, write to disk, etc.
%           see defaults/defaultOutputParams_BPHMM.m for details
%OUTPUT
%  Markov Chain state variables saved at preset frequency to hard drive
%     at filepath location specified in outParams.saveDir

function [ChainHist] = RunTimedMCMCSimForBPHMM( data, Psi, algParams, outParams, model )
tic;

if isfield( Psi, 'F' )
    % Stating chain from scratch
    n = 0;
    logPr = calcJointLogPr_BPHMMState( Psi, data );
    ChainHist = recordMCMCHistory_BPHMM( 0, outParams, [], Psi, logPr  );

    fprintf( 'Initial Config: \n' );
    printTimedMCMCSummary_BPHMM( 0, Psi, logPr, algParams);
else
    ChainHist = Psi;
    Psi = unpackBPHMMState(  ChainHist.Psi(end), data, model );
    logPr = calcJointLogPr_BPHMMState( Psi, data );
    n = ChainHist.iters.Psi(end );
    fprintf( 'Resumed Config: \n' );
    printTimedMCMCSummary_BPHMM( 0, Psi, logPr, algParams);
end

fprintf( 'Running MCMC Sampler %d : %d ... \n', outParams.jobID, outParams.taskID );

while toc < algParams.TimeLimit
    n = n + 1;
    Psi.iter = n;

    % Perform 1 iteration of MCMC, moving to next Markov state!
    [Psi, Stats] = BPHMMsample( Psi, data, algParams );

    % Diagnose convergence by calculating joint log pr. of all sampled vars
    if n == 1 || rem(n, outParams.logPrEvery)==0
        % NB: not passing "data" as arg means Psi stores all X suff stats
        logPr = calcJointLogPr_BPHMMState( Psi );
    end

    %Record current sampler state
    %  NB: internally only records at preset frequency
    ChainHist = recordMCMCHistory_BPHMM( n, outParams, ChainHist, Psi, logPr, Stats );

    doSaveToDisk = n==1 || rem(n, outParams.saveEvery)==0 || toc > algParams.TimeLimit;
    if doSaveToDisk
        filename = fullfile( outParams.saveDir,  'SamplerOutput.mat' );
        save(filename, '-struct', 'ChainHist');
    end

    if n == 1 || rem(n, outParams.printEvery)==0
       printTimedMCMCSummary_BPHMM( n, Psi, logPr, algParams);
    end

end % loop over sampler iterations

fprintf( '<<<<< --------------------------------------------------- \n');

end % main function

author: michaelchughes

repo: NPBayesHMM


text:
% This is selgp3dB_x

figure;

plot(Da(:,1),-Da(:,7),'.k','MarkerSize',1);
xlabel('Distance in [km]')
ylabel('Depth in [km]')
axis image
set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.s,...
    'LineWidth',1.0,'TickDir','out','Ticklength',[0.02 0.02],...
    'Box','on')

ax = findobj('Tag','main_map_ax');
[x,y, mouse_points_overlay] = select_polygon(ax);

zmap_message_center.set_info('Message',' Thank you .... ')

plos2 = plot(x,y,'b-','era','xor', 'Color', 'r');        % plot outline
sum3 = 0.;
pause(0.3)

%create a rectangular grid
xvect=[min(x):dx:max(x)];
yvect=[min(y):dy:max(y)];
zvect=[z2:dz:z1];

gx = xvect;gy = yvect;
tmpgri=zeros((length(xvect)*length(yvect)),2);
tmpgri2=zeros((length(xvect)*length(zvect)),2);

n=0;

for i=1:length(xvect)
    for j=1:length(yvect)
        n=n+1;
        tmpgri(n,:)=[xvect(i) yvect(j)];
    end
end

%extract all gridpoints in chosen polygon
XI=tmpgri(:,1);
YI=tmpgri(:,2);

    ll = polygon_filter(x,y, XI, YI, 'inside');
%grid points in polygon
newgri=tmpgri(ll,:);

newgri=tmpgri(ll,:);

n2 = repmat(newgri,length(zvect),1);
k = repmat(zvect',length(newgri),1);
t3 = [n2(:,1)  k n2(:,2) ];

%n3 = repmat(tmpgri2,length(yvect),1);
%k = repmat((1:1:length(yvect)),length(tmpgri),1);
%k = reshape(k,length(k)*length(yvect),1);
%t4 = [t3 n3 k];

%l = t4(:,4) == 1;
%t5 = t4(l,:);

% Plot all grid points
plot(newgri(:,1),newgri(:,2),'+k','era','back')

figure
plot3(Da(:,1),Da(:,2),-Da(:,7),'.k','MarkerSize',1)
hold on
box on ; %axis image

plot3(t3(:,1),t3(:,2),t3(:,3),'+')

if length(xvect) < 2  ||  length(yvect) < 2
    errordlg('Selection too small! (not a matrix)');
    return
end

itotal = length(newgri(:,1));
if length(gx) < 4  ||  length(gy) < 4
    errordlg('Selection too small! ');
    return
end

author: CelsoReyes


text:
function [] = chirpDetector(audioFile)

[b a]=butter(3,[25000./150000 75000./150000],'bandpass'); % create filter for chirp range
c=1;

for i = 1:150000:length(y)
   spectrogram(double(y(i:i+150000)),1000,[],[],Fs,'yaxis');
   pause(.1)

   meanPower(c) = mean(filtfilt(b,a,double(y(i:i+150000))));
   stdPower(c) = std(filtfilt(b,a,double(y(i:i+150000))));
   c=c+1;
end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function timer_id = stic
%STIC TIC equivalent using precise system time.
% Usage:
%   stic
%   timer_id = tic
%
% See also: stoc, systime

global stic_timers;

timer_id = numel(stic_timers) + 1;
stic_timers(timer_id) = systime;

if nargout < 1
    clear timer_id
end

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function mtrCreateConTracParallelScript(machineList, binDirName, exeConTracFile, scriptFileName, scriptOptsFile, pathsRoot)

% createConTrackParallelScript(machineList, binDirName, exeConTracFile,
% scriptFileName, scriptOptsFile, pathsRoot)

% Create script for parallel processing of pathway tracing on DLs
NN = length(machineList); % compute nodes
fid = fopen(scriptFileName,'wt');
fprintf(fid,'#!/bin/bash\n');
offset = 0;
for nn = 1:NN
    mach_num = floor((nn-1)/2) + 1;
    sub_n = 2 - mod(nn,2);
    paramsFile = fullfile(binDirName,'conTrack',scriptOptsFile);
    pathsFile = fullfile(binDirName,'conTrack',sprintf('%s_%d.dat',pathsRoot,nn+offset));
    fprintf(fid,'echo ''Compute %s''\n',machineList{nn});
    fprintf(fid,'ssh %s %s -i %s -p %s &\n', machineList{nn}, exeConTracFile, paramsFile, pathsFile);
    fprintf(fid,'sleep 2\n');
end
fclose(fid);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% MINIMAL_REGRESSION
%   Demonstrates automatic differentiation of a least-squares objective.

run('../../setup_autonn.m') ;  % add AutoNN to the path
rng(0) ;  % set random seed

% load simple data (4 features, 3 classes)
s = load('fisheriris.mat') ;
data_x = single(s.meas.') ;  % features-by-samples matrix
[~, ~, data_y] = unique(s.species) ;  % convert strings to class labels

% define inputs and parameters
x = Input() ;
y = Input() ;
w = Param('value', 0.01 * randn(3, 4, 'single')) ;
b = Param('value', 0.01 * randn(3, 1, 'single')) ;

% combine them using math operators, which define the prediction
prediction = w * x + b ;

% compute least-squares loss
loss = sum(sum((prediction - y).^2)) ;

% use workspace variables' names as the layers' names, and compile net
Layer.workspaceNames() ;
net = Net(loss) ;

% simple SGD
learningRate = 1e-5 ;
outputs = zeros(1, 100) ;

for iter = 1:100
  % draw minibatch
  idx = randperm(numel(data_y), 50) ;

  % evaluate network to compute gradients
  net.eval({x, data_x(:,idx), y, data_y(idx)'}) ;

  % update weights
  net.setValue(w, net.getValue(w) - learningRate * net.getDer(w)) ;
  net.setValue(b, net.getValue(b) - learningRate * net.getDer(b)) ;

  % plot loss
  outputs(iter) = net.getValue(loss) ;
end

figure(3) ;
plot(outputs) ;
xlabel('Iteration') ; ylabel('Loss') ;

loss
net

repo: MFT


text:
function [ r c ] = corner_ST( img,ptNum )
%[ r c ] = corner_ST( img )
%	Find corner points in image using Shi-Tomasi method
%   IMG should be gray.R and C are the coordinates of corners(col vec)

img = im2double(img);
h1 = fspecial('gauss',[6 6],1.5);
hd = imfilter(h1,[-1 1],'rep');
hd = hd(1:5,1:5); % use the HOG gradient
% hd = [-1 0 1];
imgdx = imfilter(img,hd,'rep');
imgdy = imfilter(img,hd','rep');

ker = fspecial('gauss',[5 5],1.5); % weight of the patch
imgdx2 = imfilter(imgdx.^2,ker,'rep');
imgdxdy = imfilter(imgdx.*imgdy,ker,'rep');
imgdy2 = imfilter(imgdy.^2,ker,'rep');

A = zeros(2,2,numel(img));
A(1,1,:) = imgdx2(:);
A(1,2,:) = imgdxdy(:);
A(2,1,:) = imgdxdy(:);
A(2,2,:) = imgdy2(:);
ev = real(eig2(A));

evTh = .1;
minEv = min(ev,[],2);
minEv = reshape(minEv,size(img));
R = (minEv > evTh*max(max(minEv)));
% fis(R)

Rmax = ordfilt2(minEv,9,ones(3));
Rcorner = (minEv == Rmax) & R; % non-maximal suppression
% fis(Rcorner)
I = find(Rcorner);

if exist('ptNum','var') && nnz(Rcorner) > ptNum
	selEvs = minEv(I); % find the ptNum corners with the largest ev
	[~,I1] = sort(selEvs,'descend');
	I = I(I1(1:ptNum));
end
[r c] = ind2sub(size(img),I); % points' positions

end

function ev = eig2(a) % eigenvalue of a 2*2 matrix
b = -(a(1,1,:)+a(2,2,:));
detA = a(1,1,:).*a(2,2,:)-a(1,2,:).*a(2,1,:);
delta = b.^2-4*detA;
ev = [-b+sqrt(delta),-b-sqrt(delta)]/2;
ev = reshape(ev,2,[])'; % N-by-2
end

author: Sable

repo: mcbench-benchmarks


text:
function k = gaussianwhiteKernDiagCompute(kern, x)

% GAUSSIANWHITEKERNDIAGCOMPUTE Compute diagonal of gaussian white kernel.
% FORMAT
% DESC computes the diagonal of the kernel matrix for the gaussian white kernel
% given a design matrix of inputs.
% RETURN K : a vector containing the diagonal of the kernel matrix computed
% at the given points.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG X : input data matrix in the form of a design matrix.
%
% SEEALSO : gaussianwhiteKernParamInit, kernDiagCompute, kernCreate,
% gaussianwhiteKernCompute
%
% COPYRIGHT : Mauricio A. Alvarez and Neil D. Lawrence, 2008
%
% MODIFICATIONS : Mauricio A. Alvarez, 2009.

% KERN

k = kern.sigma2Noise*ones(size(x,1),1);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, 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
function qp = plus(q1, q2)
%PLUS Add two quaternion objects
%
% Q1+Q2 is the element-wise sum of quaternion elements.

    if isa(q1, 'Quaternion') & isa(q2, 'Quaternion')

        qp = Quaternion(double(q1) + double(q2));
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [p,f] = spm_powell(p,xi,tolsc,func,varargin)
% Powell optimisation method
% FORMAT [p,f] = spm_powell(p,xi,tolsc,func,varargin)
%   p        - Starting parameter values
%   xi       - columns containing directions in which to begin searching
%   tolsc    - stopping criteria, optimisation stops when
%                sqrt(sum(((p-p_prev)./tolsc).^2))<1
%   func     - name of evaluated function
%   varargin - remaining arguments to func (after p)
%
%   p        - final parameter estimates
%   f        - function value at minimum
%__________________________________________________________________________
%
% Method is based on Powell's optimisation method described in
% Numerical Recipes (Press, Flannery, Teukolsky & Vetterling).
%__________________________________________________________________________
% Copyright (C) 2001-2017 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_powell.m 7252 2018-01-31 15:56:56Z john $

p = p(:);
f = feval(func,p,varargin{:});
for iter=1:512
    %if numel(p)>1, fprintf('iteration %d...\n', iter); end;            %-#
    ibig = numel(p);
    pp   = p;
    fp   = f;
    del  = 0;
    for i=1:length(p)
        ft = f;
        [p,junk,f] = min1d(p,xi(:,i),func,f,tolsc,varargin{:});
        if abs(ft-f) > del
            del  = abs(ft-f);
            ibig = i;
        end
    end
    if numel(p)==1 || sqrt(sum(((p(:)-pp(:))./tolsc(:)).^2))<1 || abs((f-fp)/(f+fp))<1e-6, return; end
    ft = feval(func,2.0*p-pp,varargin{:});
    if ft < f
        [p,xi(:,ibig),f] = min1d(p,p-pp,func,f,tolsc,varargin{:});
    end
end
warning('Too many optimisation iterations');

%==========================================================================
% function [p,pi,f] = min1d(p,pi,func,f,tolsc,varargin)
%==========================================================================
function [p,pi,f] = min1d(p,pi,func,f,tolsc,varargin)
% Line search for minimum.

global lnm % used in funeval
lnm      = struct('p',p,'pi',pi,'func',func,'args',[]);
lnm.args = varargin;

min1d_plot('Init', 'Line Minimisation','Function','Parameter Value');
min1d_plot('Set', 0, f);

tol      = 1/sqrt(sum((pi(:)./tolsc(:)).^2));
t        = bracket(f);
[f,pmin] = search(t,tol);
pi       = pi*pmin;
p        = p + pi;

%if length(p)<12,
%   for i=1:length(p), fprintf('%-8.4g ', p(i)); end;                  %-#
%   fprintf('| %.5g\n', f);                                            %-#
%else
%    fprintf('%.5g\n', f);                                              %-#
%end
min1d_plot('Clear');

%==========================================================================
% function f = funeval(p)
%==========================================================================
function f = funeval(p)
% Reconstruct parameters and evaluate.

global lnm % defined in min1d
pt = lnm.p+p.*lnm.pi;
f  = feval(lnm.func,pt,lnm.args{:});
min1d_plot('Set',p,f);

%==========================================================================
% function t = bracket(f)
%==========================================================================
function t = bracket(f)
% Bracket the minimum (t(2)) between t(1) and t(3)

gold   = (1+sqrt(5))/2; % Golden ratio

t(1)   = struct('p',0,'f',f);
t(2).p = 1;
t(2).f = funeval(t(2).p);

% if t(2) not better than t(1) then swap
if t(2).f > t(1).f
    t(3) = t(1);
    t(1) = t(2);
    t(2) = t(3);
end

t(3).p = t(2).p + gold*(t(2).p-t(1).p);
t(3).f = funeval(t(3).p);

while t(2).f > t(3).f

    % fit a polynomial to t
    tmp = cat(1,t.p)-t(2).p;
    pol = pinv([ones(3,1) tmp tmp.^2])*cat(1,t.f);

    % minimum is when gradient of polynomial is zero
    % sign of pol(3) (the 2nd deriv) should be +ve
    if pol(3)>0
        % minimum is when gradient of polynomial is zero
        d    = -pol(2)/(2*pol(3)+eps);

        % A very conservative constraint on the displacement
        if d > (1+gold)*(t(3).p-t(2).p)
            d = (1+gold)*(t(3).p-t(2).p);
        end
        u.p  = t(2).p+d;
    else
        % sign of pol(3) (the 2nd deriv) is not +ve
        % so extend out by golden ratio instead
        u.p  = t(3).p+gold*(t(3).p-t(2).p);
    end

    % FUNCTION EVALUATION
    u.f  = funeval(u.p);

    if (t(2).p < u.p) == (u.p < t(3).p)

        % u is between t(2) and t(3)
        if u.f < t(3).f
            % minimum between t(2) and t(3) - done
            t(1) = t(2);
            t(2) = u;
            return
        elseif u.f > t(2).f
            % minimum between t(1) and u - done
            t(3) = u;
            return;
        end
    end

    % Move all 3 points along
    t(1) = t(2);
    t(2) = t(3);
    t(3) = u;
end

%==========================================================================
% function [f,p] = search(t, tol)
%==========================================================================
function [f,p] = search(t, tol)
% Brent's method for line searching - given that minimum is bracketed

gold1 = 1-(sqrt(5)-1)/2;

% Current and previous displacements
d     = Inf;
pd    = Inf;

% sort t into best first order
[junk,ind] = sort(cat(1,t.f));
t   = t(ind);
brk = [min(cat(1,t.p)) max(cat(1,t.p))];

for iter=1:128
    % check stopping criterion
    if abs(t(1).p - 0.5*(brk(1)+brk(2)))+0.5*(brk(2)-brk(1)) <= 2*tol
        p = t(1).p;
        f = t(1).f;
        return;
    end

    % keep last two displacents
    ppd = pd;
    pd  = d;

    % fit a polynomial to t
    tmp = cat(1,t.p)-t(1).p;
    pol = pinv([ones(3,1) tmp tmp.^2])*cat(1,t.f);

    % minimum is when gradient of polynomial is zero
    d   = -pol(2)/(2*pol(3)+eps);
    u.p = t(1).p+d;

    % check so that displacement is less than the last but two,
    % that the displaced point is between the brackets
    % and that the solution is a minimum rather than a maximum
    eps2 = 2*eps*abs(t(1).p)+eps;
    if abs(d) > abs(ppd)/2 || u.p < brk(1)+eps2 || u.p > brk(2)-eps2 || pol(3)<=0
        % if criteria are not met, then golden search into the larger part
        if t(1).p >= 0.5*(brk(1)+brk(2))
            d = gold1*(brk(1)-t(1).p);
        else
            d = gold1*(brk(2)-t(1).p);
        end
        u.p = t(1).p+d;
    end

    % FUNCTION EVALUATION
    u.f = funeval(u.p);

    % Insert the new point into the appropriate position and update
    % the brackets if necessary
    if u.f <= t(1).f
        if u.p >= t(1).p, brk(1)=t(1).p; else brk(2)=t(1).p; end
        t(3) = t(2);
        t(2) = t(1);
        t(1) = u;
    else
        if u.p < t(1).p, brk(1)=u.p; else brk(2)=u.p; end
        if u.f <= t(2).f
            t(3) = t(2);
            t(2) = u;
        elseif u.f <= t(3).f
            t(3) = u;
        end
    end
end

%==========================================================================
% function min1d_plot(action,arg1,arg2,arg3)
%==========================================================================
function min1d_plot(action,arg1,arg2,arg3)
% Visual output for line minimisation
persistent min1dplot

if ~nargin, action = 'Init'; end

% Find the Interactive window and exit if not
%--------------------------------------------------------------------------
fg = spm_figure('FindWin','Interactive');
if isempty(fg), return; end

%-Initialize
%--------------------------------------------------------------------------
if strcmpi(action,'init')
    if nargin<4, arg3 = 'Function';          end
    if nargin<3, arg2 = 'Value';             end
    if nargin<2, arg1 = 'Line minimisation'; end

    min1dplot = struct('pointer',get(fg,'Pointer'),...
                       'name',   get(fg,'Name'),...
                       'ax',     []);
    min1d_plot('Clear');
    set(fg,'Pointer','Watch');
    min1dplot.ax = axes('Position', [0.15 0.1 0.8 0.75],...
                        'Box',      'on',...
                        'Parent',   fg);
    lab = get(min1dplot.ax,'Xlabel');
    set(lab,'string',arg3,'FontSize',10);
    lab = get(min1dplot.ax,'Ylabel');
    set(lab,'string',arg2,'FontSize',10);
    lab = get(min1dplot.ax,'Title');
    set(lab,'string',arg1);
    line('Xdata',[], 'Ydata',[],...
        'LineWidth',2,'Tag','LinMinPlot',...
        'LineStyle','-','Marker','o',...
        'Parent',min1dplot.ax);
    drawnow;

%-Reset
%--------------------------------------------------------------------------
elseif strcmpi(action,'set')
    br = findobj(fg,'Tag','LinMinPlot');
    if ~isempty(br)
        [xd,indx] = sort([get(br,'Xdata') arg1]);
        yd = [get(br,'Ydata') arg2];
        yd = yd(indx);
        set(br,'Ydata',yd,'Xdata',xd);
        drawnow;
    end

%-Clear
%--------------------------------------------------------------------------
elseif strcmpi(action,'clear')
    fg = spm_figure('FindWin','Interactive');
    if isstruct(min1dplot)
        if ishandle(min1dplot.ax), delete(min1dplot.ax); end
        set(fg,'Pointer',min1dplot.pointer);
        set(fg,'Name',min1dplot.name);
    end
    spm_figure('Clear',fg);
    drawnow;
end

author: spm


text:
function saveSameSize(h, varargin)
%SAVESAMESIZE print/export figure h at same size as on screen.
%   SAVESAMESIZE will adjust printing/exporting related properties so the
%   figure (h) will print/export the same size as the onscreen display. It
%   will print/export the figure and then restore the changed properties to
%   the original values
%
%   saveSameSize accepts some optional parameter-value pairs. If used you must
%   specify BOTH the parameter name (e.g. 'format') AND the parameter value
%   (e.g. 'png'):
%       parameter      value/meaning
%       ----------     --------------------------------------
%        'format'      output format to generate.
%                      Must be one of the formats supported by the print command
%
%                      If not specified, jpeg is used as the default
%
%        'file'        path/name of output file to create
%
%                      If not specified, 'figure<N>.<ext>' will be used
%                      where  <N> is the figure handle and ext
%                             <ext> is the extension corresponding to the
%                                   output format chosen
%
%        'renderer'    Renderer to use when creating the output.
%                      Must be one renderers supported by MATLAB
%
%                      If not specified, the figure's current renderer will be used
%
%   Example usage:
%       save current figure to jpeg file named figure<n>.jpg using the
%       figure's current renderer
%           saveSameSize(gcf);
%
%       save current figure to png file named test.png using the
%       OpenGL renderer
%           saveSameSize(gcf, 'format', 'png', 'renderer', 'opengl',...
%                          'file', 'test.png');
%
% See also PRINT

%   Copyright 2007 The MathWorks, Inc.

if ishandle(h)
    setprops = {'PaperType', 'PaperSize', 'PaperUnits', 'PaperPosition', 'PaperPositionMode'};
    try
        settings = get(h, setprops);
    catch
        error('saveSameSize unable to get figure properties');
    end

    oldSettings = cell2struct(settings, setprops, 2); % to use when resetting values
    set(h, 'PaperPositionMode', 'auto');

    paramSettings = processArgs(h, varargin{:});
    if isstruct(paramSettings)
        try
            print(h, '-r0', ...
                 paramSettings.renderer, ...
                 paramSettings.device, ...
                 paramSettings.outfile);
        catch ex
            % print had a problem - restore settings and rethrow
            set(h, oldSettings);
            rethrow(ex);
        end
        set(h, oldSettings);
    end
else
    error('saveSameSize requires a handle to a figure');
end
end

% helper function to deal with parsing command line arguments
%   expecting   ..., 'format', ['jpeg', 'eps', 'png', ...]
%                    'file', 'pathToOutputFileToCreate'
%                    'renderer', 'opengl|painters|zbuffer'
function param_settings = processArgs(h, varargin)

    if mod(length(varargin),2) ~= 0
       error('saveSameSize:InvalidOptionalParameters', ...
             'saveSameSize optional parameters must be specified as parameter-value pairs');
    end
    %setup defaults
    param_settings.outfile = '';
    ext = 'jpg';
    param_settings.device = '-djpeg';
    param_settings.renderer = lower(['-' get(h, 'renderer')]);

    for i = 1 : 2 : length(varargin)-1
        param_arg = varargin{i};
        % caller specified output device to use
        if ischar(param_arg) && strcmpi(param_arg, 'format')
           val_arg = varargin{i+1};
           % remove '-d' if caller specified it
           if length(val_arg) >= 3 && strcmpi(val_arg(1:2),'-d')
               val_arg = val_arg(3:end);
           end
           % check device against list of supported devices
           [device, ext] = localIsGoodDevice(val_arg);
           param_settings.device = ['-d' device];
        end

        % caller specified output file to use
        if ischar(param_arg) && strcmpi(param_arg, 'file')
           val_arg = varargin{i+1};
           if ischar(val_arg)
    		  param_settings.outfile = val_arg;
           else
              error('saveSameSize:InvalidOutputFile', ...
                     'saveSameSize invalid output file specified');
           end
        end

        % caller specified renderer to use
        if ischar(param_arg) && strcmpi(param_arg, 'renderer')
           val_arg = varargin{i+1};
           if ischar(val_arg)
              if strcmpi(val_arg, 'OpenGL')
                 renderer = 'opengl';
              elseif strcmpi(val_arg, 'painters')
                 renderer = 'painters';
              elseif strcmpi(val_arg, 'zbuffer')
                 renderer = 'zbuffer';
              else
                 error('saveSameSize:InvalidRenderer', ...
                     'saveSameSize invalid renderer <%s> specified', val_arg);
              end
              param_settings.renderer = ['-' renderer];
           elseif isnumeric(val_arg)
               error('saveSameSize:InvalidRenderer', ...
                     'saveSameSize invalid renderer <%s> specified', num2str(val_arg));
           else
               error('saveSameSize:InvalidRenderer', ...
                     'saveSameSize invalid renderer specified');
           end
        end
    end
    if isempty(param_settings.outfile)
       param_settings.outfile = ['figure' num2str(h) '.' ext];
    end
end

% helper function to check specified device (dev) against the list of
% supported devices.
% Returns device/format string for the supported device and the extension
% that corresponds to output of that format.
function [device, ext] = localIsGoodDevice(dev)
  pj = printjob; %#ok
  device = '';
  [options, devices, extensions ] = printtables( pj ); %#ok
  devIndex = strmatch( dev, devices, 'exact' );
  if length(devIndex) == 1 % one exact match
     device = dev;
  elseif length(devIndex) > 1
     devIndex = strmatch( dev, devices);
     if length(devIndex) == 1 % one partial match
        device = devices{devIndex};
     end
  end
  if isempty(device)
    error('saveSameSize:InvalidFormat', ...
          'saveSameSize invalid format specified ');
  end
  ext = extensions{devIndex};
  if isempty(ext)
    error('saveSameSize:NotFileFormat', ...
          'saveSameSize format specified does not support exporting to a file ');
  end
end


text:
function [S, f] = calculateSegmentSpectrum(data, lineNoise)
% Multi-taper segmented spectrum for a univariate continuous process
%
% Usage:
%     [S, f] = calculateSegmentSpectrum(data, lineNoise)
%
% Parameters:
%      data      (single channel) -- required
%      lineNoise   structure with various parameters set
%
% The lineNoise structure has the following fields set:
%       fPassBand       Frequency band used
%       Fs 	            Sampling frequency
%       pad             FFT padding factor
%       tapers          Precomputed tapers from dpss
%       taperWindowSize Taper sliding window length
%
% Output:
%       S       Spectrum
%       f       Frequencies
%

%% Check input arguments for consistency
if nargin < 2
    error('calculateSegmentSpectrum:NotEnoughArguments', ...
        'Need to provide data and segment information arguments');
end

data = change_row_to_column(data);
if size(data, 2) ~= 1;
    error('calculateSegmentSpectrum:DataNot1Dim', ...
        'Data must beunivariate time series');
end
%% Extract argument values
win = getStructureParameters(lineNoise, 'taperWindowSize', 4);
Fs = getStructureParameters(lineNoise, 'Fs', 1);
pad = getStructureParameters(lineNoise, 'pad', 0);
fpass = getStructureParameters(lineNoise, 'fPassBand', [0 lineNoise.Fs/2]);

%% Create the segmented data for the calculation of the spectrum
N = size(data, 1); % length of segmented data
dt = 1/Fs; % sampling interval
T = N*dt; % length of data in seconds
E = 0:win:(T - win); % fictitious event triggers
win = [0, win]; % use window length to define left and right limits of windows around triggers
data = createdatamatc(data, E, Fs, win); % segmented data
N = size(data,1); % length of segmented data
nfft = max(2^(nextpow2(N) + pad), N);
[f, findx] = getfgrid(Fs, nfft, fpass);
tapers = lineNoise.tapers;
J = mtfftc(data, tapers, nfft, Fs); % compute tapered fourier transforms
J = J(findx, :, :); % restrict to specified frequencies
S = squeeze(mean(conj(J).*J, 2)); % spectra of non-overlapping segments (average over tapers)
S = squeeze(mean(S, 2)); % Mean of the spectrum averaged across segments


text:
function edge2 = clipEdge(edge, box)
%CLIPEDGE Clip an edge with a rectangular box
%
%   EDGE2 = clipEdge(EDGE, BOX);
%   EDGE: [x1 y1 x2 y2],
%   BOX : [xmin xmax ymin ymax], or [xmin xmax ; ymin ymax].
%   return :
%   EDGE2 = [xc1 yc1 xc2 yc2];
%
%   If clipping is null, return [0 0 0 0];
%
%   if EDGE is a N-by-4 array, return an N-by-4 array, corresponding to each
%   clipped edge.
%
%   See also
%   edges2d, boxes2d, clipLine
%
% ---------
% author : David Legland
% created the 14/05/2005.
% Copyright 2010 INRA - Cepia Software Platform.

%   HISTORY
%   2007-01-08 sort points according to position on edge, not to x coord
%       -> this allows to return edges with same orientation a source, and
%       to keep first or end points at the same position if their are not
%       clipped.
%   01/10/2010 fix bug due to precision, thanks to Reto Zingg.

% process data input
if size(box, 1) == 2
    box = box';
end

% get limits of window
xmin = box(1);
xmax = box(2);
ymin = box(3);
ymax = box(4);

% convert window limits into lines
lineX0 = [xmin ymin xmax-xmin 0];
lineX1 = [xmin ymax xmax-xmin 0];
lineY0 = [xmin ymin 0 ymax-ymin];
lineY1 = [xmax ymin 0 ymax-ymin];

% compute outcodes of each vertex
p11 = edge(:,1) < xmin; p21 = edge(:,3) < xmin;
p12 = edge(:,1) > xmax; p22 = edge(:,3) > xmax;
p13 = edge(:,2) < ymin; p23 = edge(:,4) < ymin;
p14 = edge(:,2) > ymax; p24 = edge(:,4) > ymax;
out1 = [p11 p12 p13 p14];
out2 = [p21 p22 p23 p24];

% detect edges totally inside window -> no clip.
inside = sum(out1 | out2, 2) == 0;

% detect edges totally outside window
outside = sum(out1 & out2, 2) > 0;

% select edges not totally outside, and process separately edges totally
% inside window
ind = find(~(inside | outside));

% allocate memroty for all clipped edges
edge2 = zeros(size(edge));

% copy result of edges totally inside clipping box
edge2(inside, :) = edge(inside, :);

for i = 1:length(ind)
    % current edge
    iedge = edge(ind(i), :);

    % compute intersection points with each line of bounding window
    px0 = intersectLineEdge(lineX0, iedge);
    px1 = intersectLineEdge(lineX1, iedge);
    py0 = intersectLineEdge(lineY0, iedge);
    py1 = intersectLineEdge(lineY1, iedge);

    % create array of points
    points  = [px0; px1; py0; py1; iedge(1:2); iedge(3:4)];

    % remove infinite points (edges parallel to box edges)
	points  = points(all(isfinite(points), 2), :);

    % sort points by x then y
    points = sortrows(points);

    % get center positions between consecutive points
    centers = (points(2:end,:) + points(1:end-1,:))/2;

    % find the centers (if any) inside window
    inside = find(  centers(:,1) >= xmin & centers(:,2) >= ymin & ...
                    centers(:,1) <= xmax & centers(:,2) <= ymax);

    % if multiple segments are inside box, which can happen due to finite
    % resolution, only take the longest segment
    if length(inside) > 1
        % compute delta vectors of the segments
        dv = points(inside+1,:) - points(inside,:);
        % compute lengths of segments
        len = hypot(dv(:,1), dv(:,2));
        % find index of longest segment
        [a, I] = max(len); %#ok<ASGLU>
        inside = inside(I);
    end

    % if one of the center points is inside box, then the according edge
    % segment is indide box
    if length(inside) == 1
         % restore same direction of edge
        if iedge(1) > iedge(3) || (iedge(1) == iedge(3) && iedge(2) > iedge(4))
            edge2(ind(i), :) = [points(inside+1,:) points(inside,:)];
        else
            edge2(ind(i), :) = [points(inside,:) points(inside+1,:)];
        end
    end

end % end of loop over edges

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open 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.
%
% dslam_open 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 dslam_open. If not, see <http://www.gnu.org/licenses/>.

function arrow = poseToQuiver(T)
arrow = zeros(1, 6);
arrow(1) = T(1,4);
arrow(2) = T(2,4);
arrow(3) = T(3,4);
arrow(4) = T(1,1);
arrow(5) = T(2,1);
arrow(6) = T(3,1);
end

author: uzh-rpg


text:
function results = run_for_parameters(f, params, varargin)
%RUN_FOR_PARAMETERS Evaluate the given function for different paramters.
%Syntax:
%   results = RUN_FOR_PARAMETERS(f, params[, ...]);
%Inputs:
%   f - Function handle whose must return a single value of the same type
%       each time it is called.
%   params - An n x k matrix or cell array where each row represents a
%            parameter vector. The given function handle should accept
%            k parameters.
%   ... - Options:
%           'Parallel' - Whether to evaluate the function in parallel.
%                        Default value is false.
%           'ParallelLevel' - Set to 0 to run each task for each parameter
%                             vector. Set to 1 to run each task for each
%                             repetition. Default value is 0.
%           'Repeat' - Specify how many times the function will be
%                      evaluated for each parameter vector. This is useful
%                      in running Monte Carlo simulations. Default value is
%                      one.
%           'ShowProgress' - Shows a progress bar. Default value is false.
%Output:
%   results - Depending on the return value of the function. Let r be the
%             number of repeats.
%             1. If f returns a scalar, the returned value will be an n x r
%                matrix;
%             2. If f returns a struct, the returned value will be an n x r
%                struct;
%             3. Otherwise, the returned value will be an n x r cell array.
use_parallel = false;
parallel_level = 0;
show_progress = false;
repeat = 1;
for ii = 1:2:nargin-2
    option_name = varargin{ii};
    option_value = varargin{ii+1};
    switch lower(option_name)
        case 'parallel'
            use_parallel = option_value;
        case 'parallellevel'
            parallel_level = option_value;
        case 'repeat'
            repeat = option_value;
        case 'showprogress'
            show_progress = option_value;
        otherwise
            error('Unknown option ''%s''.', option_name);
    end
end

results = [];
n = size(params, 1);
is_params_cell = iscell(params);
first_run = true;
if use_parallel
    pool = gcp();
    if show_progress
        if parallel_level == 0
            progressbar('reset', n);
        else
            progressbar('reset', n*repeat);
        end
    end
    first_run = true;
    if parallel_level == 0
        % one task per parameter
        % TODO: split into batches when n is large
        for nn = 1:n
            if is_params_cell
                cur_args = params(nn, :);
            else
                cur_args = num2cell(params(nn, :));
            end
            ff(nn) = parfeval(pool, @run_repetitions , 1, f, cur_args, repeat);
        end
        % fetch the results
        for nn = 1:n
            [idx, cur_results] = fetchNext(ff);
            % we need to store a batch of results here
            if first_run
                results = allocate_results(cur_results(1), n, repeat);
                first_run = false;
            end
            results(idx,:) = cur_results;
            if show_progress
                progressbar('advance');
            end
        end
    else
        % one task per repetition
        for nn = 1:n
            if is_params_cell
                cur_args = params(nn, :);
            else
                cur_args = num2cell(params(nn, :));
            end
            for rr = 1:repeat
                ff(rr) = parfeval(pool, f, 1, cur_args{:});
            end
            % we need to fetch the results for each repetition here
            for rr = 1:repeat
                [idx, cur_result] = fetchNext(ff);
                if first_run
                    results = allocate_results(cur_result, n, repeat);
                    first_run = false;
                end
                if isnumeric(results) || isstruct(results)
                    results(nn,idx) = cur_result;
                else
                    results{nn,idx} = cur_result;
                end
                if show_progress
                    progressbar('advance');
                end
            end
        end
    end
    if show_progress
        progressbar('end');
    end
else
    if show_progress
        progressbar('reset', n * repeat);
    end
    for nn = 1:n
        if is_params_cell
            cur_args = params(nn, :);
        else
            cur_args = num2cell(params(nn, :));
        end
        for rr = 1:repeat
            cur_result = f(cur_args{:});
            if first_run
                results = allocate_results(cur_result, n, repeat);
                first_run = false;
            end
            if isnumeric(results) || isstruct(results)
                results(nn,rr) = cur_result;
            else
                results{nn,rr} = cur_result;
            end
            if show_progress
                progressbar('advance');
            end
        end
    end
    if show_progress
        progressbar('end');
    end
end

end

function results = run_repetitions(f, args, repeat)
% A wrapper function to evaluate f multiple times.
first_run = true;
for rr = 1:repeat
    cur_result = f(args{:});
    if first_run
        results = allocate_results(cur_result, 1, repeat);
        first_run = false;
    end
    if isnumeric(results) || isstruct(results)
        results(rr) = cur_result;
    else
        results{rr} = cur_result;
    end
end
end

function results = allocate_results(sample_result, n, r)
% determine return type
if isnumeric(sample_result) && isscalar(sample_result)
    results = zeros(n, r);
elseif isstruct(sample_result)
    field_names = fieldnames(sample_result);
    args = reshape([field_names'; cell(1, length(field_names))], 1, []);
    results = repmat(struct(args{:}), n, r);
else
    results = cell(n, r);
end
end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
function varargout = plotcoeffs(varargin)
%PLOTCOEFFS   Display the PLOTCOEFFS of the column and row slices.
%   PLOTCOEFFS(F) plots the coefficients on a semilogy scale of the
%   underlying basis used in constructing the one-dimensional slices that
%   form F. It returns two figures one for the row slices and one for the
%   column slices.
%
%   PLOTCOEFFS(F, S) allows further plotting options, such as linestyle,
%   linecolor, etc. If S contains a string 'LOGLOG', the coefficients will be
%   displayed on a log-log scale.
%
% See also DISKFUN/PLOTCOEFFS2 and DISKFUN/COEFFS2.

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

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

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% DEMO1 - Basic ncdataset data access
echo('on')
% STARTING DEMO1 ----------------------------------------------------------
% Basic usage of ncdataset

%% ---- Open your data
% Can be NetCDF, NetCDF on webserver, or OpenDAP
ds = ncdataset('http://geoport.whoi.edu/thredds/dodsC/examples/OS_M1_20081008_TS.nc');
%% ---- You can access a list of the variables available to you
ds.variables

%% Lets fetch time in Matlab's native format
t = ds.time('TIME');

%% ---- Now lets get the data
temp = double(ds.data('TEMP')); % Convert to Double!!
depth = ds.data('DEPTH');

%% ---- Plot the data
figure;
surf(t, depth, temp.');...
view(2);shading interp;...
datetick('x', 2);set(gca, 'YDir', 'reverse');...
grid('on');ch = colorbar;...
set(get(ch, 'YLabel'), 'String', '^oC');...
title({'Temperature at M1 Mooring in Monterey Bay',ds.attribute('references'),ds.location},'interpreter','none')

shg

echo('off') % ENDING DEMO1 ------------------------------------------------

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function y = logMvGamma(x, d)
% Compute logarithm multivariate Gamma function
% which is used in the probability density function of the Wishart and inverse Wishart distributions.
% Gamma_d(x) = pi^(d(d-1)/4) \prod_(j=1)^d Gamma(x+(1-j)/2)
% log(Gamma_d(x)) = d(d-1)/4 log(pi) + \sum_(j=1)^d log(Gamma(x+(1-j)/2))
% Input:
%   x: m x n data matrix
%   d: dimension
% Output:
%   y: m x n logarithm multivariate Gamma
% Written by Michael Chen (sth4nth@gmail.com).
y = d*(d-1)/4*log(pi)+sum(gammaln(x(:)+(1-(1:d))/2),2);
y = reshape(y,size(x));

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
% eegplugin_eepimport() - EEGLAB plugin for importing ANT EEProbe data files.
%                         With this menu it is possible to import and export a continuous CNT file (*.cnt) or
%                         an averaged file (*.avr), linked with an event file (*.trg).
%
% Usage:
%   >> eegplugin_eepimport(menu);
%   >> eegplugin_eepimport(menu, trystrs, catchstrs);
%
% Inputs:
%   menu       - [float]  EEGLAB menu handle
%   trystrs    - [struct] "try" strings for menu callbacks. See notes on EEGLab plugins.
%                (http://www.sccn.ucsd.edu/eeglab/contrib.html)
%   catchstrs  - [struct] "catch" strings for menu callbacks. See notes on EEGLab plugins.
%                (http://www.sccn.ucsd.edu/eeglab/contrib.html)
%
%
% Notes:
%   This plugins consist of the following Matlab files:
%   pop_loadeep.m           pop_loadeep_avg.m
%   loadeep.m               loadeep_avg.m
%   read_eep_cnt.m          read_eep_avr.m
%   read_eep_cnt.mexglx     read_eep_avr.mexglx
%   read_eep_cnt.dll        read_eep_avr.dll
%
% Authors:
% Maarten-Jan Hoeve, ANT, The Netherlands / www.ant-software.nl, 3 October 2003
% Maarten van de Velde, ANT, The Netherlands / www.ant-neuro.com, 21 July 2005
%
% See also: eeglab(), pop_loadeep(), loadeep(), read_eep_cnt(), pop_loadeep_avg(), loadeep_avg(), read_eep_avr()

%123456789012345678901234567890123456789012345678901234567890123456789012

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

% $Log: eegplugin_eepimport.m,v $
% Revision 1.6  2006-09-25 14:25:25  mvelde
% updated version number
%
% Revision 1.5  2006-09-25 14:25:25  mvelde
% updated version number
%
% Revision 1.4  2006-09-25 14:04:02  mvelde
% updated for EEGLAB 5.03
%
% fix reading data for the new EEGLAB version - Arno August 2006
%
% Revision 1.3  2005/07/21 10:11:40  mvelde
% updated menu structure
%
% Revision 1.2  2005/06/08 08:16:37  mvelde
% converted files to unix format
%
% Revision 1.1  2004/11/26 13:22:22  jwiskerke
% Added eeglab to cvs.
%
% Revision 1.3  2003/10/24 13:24:54  Maarten-Jan Hoeve
% Added GNU Licence and updated revision history
%
%
% Advanced Neuro Technology (ANT) BV, The Netherlands, www.ant-neuro.com / info@ant-neuro.com
%

function vers = eegplugin_eepimport(fig, trystrs, catchstrs)

    vers = 'eepimport1.06';

    if nargin < 3
        error('eegplugin_eepimport requires 3 arguments');
    end;

    % add folder to path
    % ------------------
    if ~exist('pop_loadeep')
        p = which('eegplugin_eepimport.m');
        p = p(1:findstr(p,'eegplugin_eepimport.m')-1);
        addpath( p );
    end;

    % find import data menu
    % ---------------------
    menu = findobj(fig, 'tag', 'import data');

    % menu callbacks
    % --------------
    comcnt = [ trystrs.no_check '[EEG LASTCOM] = pop_loadeep;' catchstrs.new_non_empty ];
    comavr = [ trystrs.no_check '[EEG LASTCOM] = pop_loadeep_avg;' catchstrs.new_non_empty ];

    uimenu( menu, 'label', 'From ANT EEProbe .CNT file', 'callback', comcnt, 'separator', 'on' );
    uimenu( menu, 'label', 'From ANT EEProbe .AVR file' , 'callback', comavr );

%   future export
%   uimenu( submenu, 'label', 'To EEProbe CNT file format' , 'callback', expcnt );
%   uimenu( submenu, 'label', 'To EEProbe AVR file format' , 'callback', expavr );

author: PatternRecognition

repo: OpenBMI


text:
classdef (Abstract) MeasurementReaderX < ReaderX
% MeasurementReaderX Abstract class
%
% Summary of MeasurementReaderX
%  This is the base class for all TrackingX measurement readers. Must be
%  used as the superclass of any custom defined TrackingX measurement readers.
%
% MeasurementReaderX Properties:
%   + Measurements - The measurements at a given timestep.
%
% MeasurementReaderX Methods:
%   ~ read - Read and return measurements at current timestep
%
% (+) denotes puplic properties/methods
% (~) denotes abstract properties/methods
%
% The above parameters and methods are accessed by the majority of
% existing TrackingX library components and are COMPULSORY to guarantee
% compatibility with the TrackingX framework.

    properties
        Measurements
    end

    methods (Abstract)
        read(this);
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function dist = sphere_distance2 ( lat1, lon1, lat2, lon2, r )

%*****************************************************************************80
%
%% SPHERE_DISTANCE2 computes great circle distances on a sphere.
%
%  Discussion:
%
%    This computation is written in terms of haversines, and can be more
%    accurate when measuring small angular distances.  It can be somewhat
%    inaccurate when the two points are antipodal.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    "Great-circle distance",
%    Wikipedia.
%
%  Parameters:
%
%    Input, real LAT1, LON1, the latitude and longitude of
%    the first point.
%
%    Input, real LAT2, LON2, the latitude and longitude of
%    the second point.
%
%    Input, real R, the radius of the sphere.
%
%    Output, real DIST, the great circle distance between
%    the points, measured in the same units as R.
%
  s = ( sin ( ( lat1 - lat2 ) / 2.0 ) ).^2 ...
    + cos ( lat1 ) * cos ( lat2 ) * ( sin ( ( lon1 - lon2 ) / 2.0 ) ).^2;

  s = sqrt ( s );

  dist = 2.0 * r * asin ( s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [D] = compute_diffusion_coefficient(I)

alpha = 5.0;
beta = 0.5;

wSize = 5;
sigma = wSize/6;

gI = imfilter(I, fspecial('gaussian', [wSize wSize], sigma),'replicate');
mask = [-1 0 1];
Ix = imfilter(gI, mask, 'replicate');
Iy = imfilter(gI, mask', 'replicate');

norm = sqrt(Ix.^2 + Iy.^2);

D = max(1e-06, exp(-alpha*norm.^beta));

end

author: thomas-koehler

repo: SupER


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

%% Syntax
% |[VN]=subCurve(Vt,np,closeLoopOpt);|

%% Description
% The |subCurve| function can be used to increase the point density of the
% input curve by adding evenly spaced points between each of the curve
% segments.

%% Examples

clear; close all; clc;

%%
% PLOT SETTINGS
fontSize=15;
markerSize1=45;
lineWidth1=2;
lineWidth2=5;
lineWidth3=2;
faceAlpha=0.5;

%% Example: Linearly upsample an open ended curve with intermediate points

%%
% Simulating a curve
Vt=[0 0 0; 10 0 0; 5 10 0; 10 0 10; 0 10 10; ];

%%
% Upsampling the curve
np=3; %Number of desired intermediate points to be added
[VN]=subCurve(Vt,np); %Using subcurve to upsample curve

%%
% Plotting results
hf1=cFigure;
title('A linearly upsampled curve','FontSize',fontSize);
xlabel('X','FontSize',fontSize);ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

plotV(Vt,'k.-.','lineWidth',lineWidth1,'MarkerSize',markerSize1);
plotV(VN,'r.-','lineWidth',lineWidth1/2,'MarkerSize',markerSize1/2);

axis equal; view(3); axis tight;  grid on;  set(gca,'FontSize',fontSize);
drawnow;

%% Example: Linearly upsample a closed curve with intermediate points

%%
% Upsampling the curve
np=3; %Number of desired intermediate points to be added
closeLoopOpt=1; %Enable closed loop option
[VN]=subCurve(Vt,np,closeLoopOpt); %Using subcurve to upsample curve

%%
% Plotting results
hf1=cFigure;
title('A linearly upsampled curve with closed end condition','FontSize',fontSize);
xlabel('X','FontSize',fontSize);ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

plotV(Vt,'k.-.','lineWidth',lineWidth1,'MarkerSize',markerSize1);
plotV(VN,'r.-','lineWidth',lineWidth1/2,'MarkerSize',markerSize1/2);

axis equal; view(3); axis tight;  grid on;  set(gca,'FontSize',fontSize);
drawnow;

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

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [ a, seed ] = r8vec_uniform ( n, b, c, seed )

%*****************************************************************************80
%
%% R8VEC_UNIFORM returns a scaled pseudorandom R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real B, C, the range of the pseudorandom values.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real A(N), the vector of pseudorandom values.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_UNIFORM - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R8VEC_UNIFORM - Fatal error!' );
  end

  for i = 1 : n

    k = floor ( seed / 127773 );

    seed = 16807 * ( seed - k * 127773 ) - k * 2836;

    if ( seed < 0 )
      seed = seed + 2147483647;
    end

    a(i) = b + ( c - b ) * seed * 4.656612875E-10;

  end

  return
end


text:
function N = nees(x,m,P)

% NEES  Normalized Estimation Error Squared.
%   N = NEES(X,M,P) computes the Normalized Estimation Error Squared given
%   true state X and a Gaussian estimate of mean M and covariance P:
%
%       NEES = (X-M)'*inv(P)*(X-M);
%
%   See also MAHALANOBIS.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

N = (x-m)'/P*(x-m);

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
%%*******************************************************************
%% detect_ublk: search for implied free variables in linear
%%              block.
%% [blk2,At2,C2,ublkinfo] = detect_ublk(blk,At,C);
%%
%% i1,i2: indices corresponding to splitting of unrestricted varaibles
%% i3   : remaining indices in the linear block
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*******************************************************************

function [blk2,At2,C2,ublkinfo,parbarrier2,X2,Z2] = ...
    detect_ublk(blk,At,C,parbarrier,X,Z,printlevel)

if (nargin < 7); printlevel = 1; end

blk2 = blk; At2 = At; C2 = C;
if (nargin >= 6)
    parbarrier2 = parbarrier;
    X2 = X; Z2 = Z;
else
    X2 = []; Z2 = [];
end
numblk = size(blk,1);
ublkinfo = cell(size(blk,1),3);
tol = 1e-14;
%%
numblknew = numblk;
%%
for p = 1:numblk
    pblk = blk(p,:);
    m = size(At{p},2);
    if strcmp(pblk{1},'l')
        r = randmat(1,m,0,'n');
        % stime = cputime;
        Ap = At{p}'; Cp = C{p};
        ApTr = (r*Ap)';
        [sApTr,perm] = sort(abs(ApTr));
        idx0 = find(abs(diff(sApTr)) < tol);
        if ~isempty(idx0)
            n = pblk{2};
            i1 = perm(idx0); i2 = perm(idx0+1);
            Api1 = Ap(:,i1);
            Api2 = Ap(:,i2);
            Cpi1 = Cp(i1)';
            Cpi2 = Cp(i2)';
            idxzr = abs(Cpi1+Cpi2) < tol & sum(abs(Api1+Api2),1) < tol;
            if any(idxzr)
                i1 = i1(idxzr');
                i2 = i2(idxzr');
                blk2{p,1} = 'u';
                blk2{p,2} = length(i1);
                At2{p} = Ap(:,i1)';
                C2{p}  = Cp(i1);
                if (printlevel)
                    fprintf('\n %1.0d linear variables from unrestricted variable.\n',...
                        2*length(i1));
                end
                if (nargin >= 6)
                    parbarrier2{p} = parbarrier{p}(i1);
                    X2{p} = X{p}(i1)-X{p}(i2);
                    Z2{p} = zeros(length(i1),1);
                end
                i3 = setdiff(1:n,union(i1,i2));
                if ~isempty(i3)
                    numblknew = numblknew + 1;
                    blk2{numblknew,1} = 'l';
                    blk2{numblknew,2} = length(i3);
                    At2{numblknew,1}  = Ap(:,i3)';
                    C2{numblknew,1}   = Cp(i3);
                    if (nargin >= 6)
                        parbarrier2{numblknew,1} = parbarrier{p}(i3);
                        X2{numblknew,1} = X{p}(i3); Z2{numblknew,1} = Z{p}(i3);
                    end
                end
                ublkinfo{p,1} = i1; ublkinfo{p,2} = i2; ublkinfo{p,3} = i3;
            end
        end
    end
end
%%*******************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function varargout = surface( f, varargin )
%SURFACE  Plot surface of a SEPARABLEAPPROX.
%   SURFACE(X, Y, Z, C) adds the surface in X,Y,Z,C to the current axes.
%
%   SURFACE(X, Y, Z) uses C = Z, so color is proportional to surface height.
%
%   See SURF for a complete description of the various forms that X,Y,Z,C can
%   take.
%
% See also SURF.

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

if ( isempty( f ) )
     surface( [ ] );
     return
end

holdState = ishold;
numpts = 200;
dom = f.domain;

% Data points:
x = linspace(dom(1), dom(2), numpts);
y = linspace(dom(3), dom(4), numpts);
[xx, yy] = meshgrid( x, y );
val = feval(f, xx, yy);

% Options:
defaultopts = {'facecolor', 'interp', 'edgealpha', .5, 'edgecolor', 'none'};

if ( isempty(varargin) )
    h1 = surface( xx, yy, val, defaultopts{:} );
    hold on
    h2 = surface( xx.', yy.', val.', defaultopts{:} );
else
    h1 = surface( xx, yy, val, defaultopts{:}, varargin{:} );
    hold on
    h2 = surface( xx.', yy.', val.', defaultopts{:}, varargin{:} );
end

if ( ~holdState )
    hold off
end

if ( nargout > 0 )
    varargout = {h1 h2};
    return
end

end


text:
% --------------------------------------------------------------------------- %
% MarrRevisited - Surface Normal Estimation
% Copyright (c) 2016 Adobe Systems Incorporated and Carnegie Mellon University.
% All rights reserved.[see LICENSE for details]
% -------------------------------------------------------------------------- %

% Written by Aayush Bansal. Please contact ab.nsit@gmail.com

% use the network to predict the values of nx, ny, and nz.
% function comp_nxyz(PREDN, NET_NAME, MODEL_ITER)
clc; clear all;

% setup the task --
DUMP_CACHE = './cachedir/best_model/';

% load the nyu-dataset --
img_set = 'test';
imgLabs = load(['./dataset/NYU/img_set/',img_set,'.mat'], 'img_set');
imgLabs = imgLabs.img_set;

% evaluate the predicted sn --
[nums_e] = eval_pred_sn(DUMP_CACHE, imgLabs);

author: aayushbansal

repo: MarrRevisited


text:
function noise = gaussianNoiseExpandParam(noise, params)

% GAUSSIANNOISEEXPANDPARAM Create noise structure from GAUSSIAN noise's parameters.
% FORMAT
% DESC returns a Gaussian noise structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG noise : the noise structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% noise structure.
% RETURN noise : noise structure with the given parameters in the
% relevant locations.
%
% SEEALSO : gaussianNoiseParamInit, gaussianNoiseExtractParam, noiseExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

noise.bias = params(1:end-1);
noise.sigma2 = params(end);

author: SheffieldML


text:
function r8vm_print ( m, n, a, title )

%*****************************************************************************80
%
%% R8VM_PRINT prints a R8VM matrix.
%
%  Discussion:
%
%    The R8VM storage format is used for an M by N Vandermonde matrix.
%    An M by N Vandermonde matrix is defined by the values in its second
%    row, which will be written here as X(1:N).  The matrix has a first
%    row of 1's, a second row equal to X(1:N), a third row whose entries
%    are the squares of the X values, up to the M-th row whose entries
%    are the (M-1)th powers of the X values.  The matrix can be stored
%    compactly by listing just the values X(1:N).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, real A(N), the R8VM matrix.
%
%    Input, string TITLE, a title to be printed.
%
  r8vm_print_some ( m, n, a, 1, 1, n, n, title );

  return
end


text:
% Test script for solving the 2D advection
Globals2D;

% Generate simple mesh
filename = 'Maxwell025.neu';

[Nv, VX, VY, K, EToV] = MeshReaderGambit2D(filename);
Nfaces = 3;
[EToE,EToF] = tiConnect2D(EToV);
BCType = Wall*(EToE==((1:K)'*ones(1,Nfaces)));

% Build mesh
Norder = ceil(10*rand(K,1));

% Set up arbitrary order elements mesh
[pinfo] = BuildPNonCon2D(Norder, K, VX, VY, EToV, BCType);

% Set initial conditions
mmode = 1; nmode = 1;

x = []; y = [];
for N1=1:max(Norder)
  pinf = pinfo(N1);
  x(pinf.ids) = pinf.x;
  y(pinf.ids) = pinf.y;
end
x = x'; y = y';
Ez = sin(mmode*pi*x).*sin(nmode*pi*y);
Hx = zeros(size(x)); Hy = zeros(size(x));

% Solve Problem for exactly one period
FinalTime = 1;
[Hx,Hy,Ez,time] = MaxwellPNonCon2D(pinfo, Hx,Hy,Ez,FinalTime);

repo: Aero-matlab


text:
function [ cost,gnd_objects,hyp_objects ] = roomLossFunction3D( gnd_objects, hyp_objects, H2G_R )
%ROOMLOSSFUNCTION Summary of this function goes here
%   Detailed explanation goes here

%% define cost
gnd_loss_cost = 2;
hyp_false_cost = 1;

%% convert hypothesis point
for i = 1:length(hyp_objects)
    hyp_objects(i).align = rotatePoint(hyp_objects(i).align, H2G_R);
end

%% greedy matching
test_type = [1:12 29];
test_type_num = length(test_type);
cost = 0;
hyp_objects(1).matchID = [];
hyp_objects(1).cost = [];

gnd_objects(1).matchID = [];
gnd_objects(1).cost = [];

for i = 1:test_type_num
    tid = test_type(i);

    hyp_rep_valid = [hyp_objects.objtype]==tid;
    hyp_rep_num = sum(hyp_rep_valid);
    hyp_rep = hyp_objects(hyp_rep_valid);
    hyp_rep_id = reshape(find(hyp_rep_valid),1,[]);

    gnd_rep_valid = [gnd_objects.objtype]==tid;
    gnd_rep_num = sum(gnd_rep_valid);
    gnd_rep = gnd_objects(gnd_rep_valid);
    gnd_rep_id = reshape(find(gnd_rep_valid),1,[]);

    if hyp_rep_num==0
        cost = cost + gnd_rep_num * gnd_loss_cost;
        for k = gnd_rep_id
            gnd_objects(k).cost = gnd_loss_cost;
            gnd_objects(k).matchID = 0;
        end

    elseif gnd_rep_num==0
        cost = cost + hyp_rep_num * hyp_false_cost;
        for k = hyp_rep_id
            hyp_objects(k).cost = hyp_false_cost;
            hyp_objects(k).matchID = 0;
        end

    else
        hyp_gnd_dist = zeros(hyp_rep_num, gnd_rep_num);
        hyp_gnd_maxs = zeros(hyp_rep_num, gnd_rep_num);
        for m = 1:hyp_rep_num
            for n = 1:gnd_rep_num

                hyp_gnd_dist(m,n) = verticesDistance(hyp_rep(m).align, gnd_rep(n).align);
                hyp_gnd_maxs(m,n) = 1;
%                 d = volumeIntersection(hyp_rep(m).point3D, gnd_rep(n).point3D);
%                 if d>=0
%                     hyp_gnd_dist(m,n) = 2*d;
%                     hyp_gnd_maxs(m,n) = 2;
%                 else
%                     hyp_gnd_dist(m,n) = verticesDistance(hyp_rep(m).point3D, gnd_rep(n).point3D);
%                     hyp_gnd_maxs(m,n) = 1;
%                 end

%                 [v1, v2, vi] = obbInteVolume(hyp_rep(m).point3D, gnd_rep(n).point3D);
%                 hyp_gnd_dist(m,n) = vi/(v1+v2-vi);
%                 p1 = hyp_rep(m).point3D;
%                 p2 = gnd_rep(n).point3D;
%                 meansize1 = (norm(p1(1,:)-p1(2,:)) + norm(p1(1,:)-p1(3,:)) + norm(p1(1,:)-p1(5,:)))/3;
%                 meansize2 = (norm(p2(1,:)-p2(2,:)) + norm(p2(1,:)-p2(3,:)) + norm(p2(1,:)-p2(5,:)))/3;
%                 meandist = mean(sqrt(sum((hyp_rep(m).point3D-gnd_rep(n).point3D).^2,2))) ...
%                     /sqrt(meansize1*meansize2);
%                 hyp_gnd_dist(m,n) = sigmoidFunc(meandist, 0.45, 4); %exp(-2*meandist);%

%                 if hyp_gnd_dist(m,n)>0.5
%                     fprintf('type: %d, hypid: %d, gndid: %d, score: %f\n', tid, hyp_rep_id(m), gnd_rep_id(n), hyp_gnd_dist(m,n));
%                 end

            end
        end

        sel_hyp = false(hyp_rep_num,1);
        sel_gnd = false(gnd_rep_num,1);
        for k = 1:min( hyp_rep_num, gnd_rep_num)
            [r,c] = find(hyp_gnd_dist==max(hyp_gnd_dist(:)));
            r = r(1); c = c(1);
            cost = cost + (hyp_gnd_maxs(r,c)-hyp_gnd_dist(r,c));
            hyp_objects(hyp_rep_id(r)).cost = (hyp_gnd_maxs(r,c)-hyp_gnd_dist(r,c))/2;
            hyp_objects(hyp_rep_id(r)).matchID = gnd_rep_id(c);
            gnd_objects(gnd_rep_id(c)).cost = (hyp_gnd_maxs(r,c)-hyp_gnd_dist(r,c))/2;
            gnd_objects(gnd_rep_id(c)).matchID = hyp_rep_id(r);

            sel_hyp(r) = true;
            sel_gnd(c) = true;
            hyp_gnd_dist(:,c) = -1;
            hyp_gnd_dist(r,:) = -1;
        end

        if hyp_rep_num>gnd_rep_num
            cost = cost + (hyp_rep_num-gnd_rep_num)*hyp_false_cost;
        else
            cost = cost + (gnd_rep_num-hyp_rep_num)*gnd_loss_cost;
        end
        rr = find(~sel_hyp);
        rc = find(~sel_gnd);
        for k = 1:length(rr)
            hyp_objects(hyp_rep_id(rr(k))).cost = hyp_false_cost;
            hyp_objects(hyp_rep_id(rr(k))).matchID = 0;
        end
        for k = 1:length(rc)
            gnd_objects(gnd_rep_id(rc(k))).cost = gnd_loss_cost;
            gnd_objects(gnd_rep_id(rc(k))).matchID = 0;
        end
    end

end

cost = cost/length(gnd_objects);
end

function d = verticesDistance(p1, p2)
% p1 = hyp_rep(m).point3D;
% p2 = gnd_rep(n).point3D;
meansize1 = (norm(p1(1,:)-p1(2,:)) + norm(p1(1,:)-p1(4,:)) + norm(p1(1,:)-p1(5,:)))/3;
meansize2 = (norm(p2(1,:)-p2(2,:)) + norm(p2(1,:)-p2(4,:)) + norm(p2(1,:)-p2(5,:)))/3;
meandist = mean(sqrt(sum((p1-p2).^2,2))) ...
    /sqrt(meansize1*meansize2);
d = sigmoidFunc(meandist, 0.45, 4); %exp(-2*meandist);%
end

% function d = volumeIntersection(p1, p2)
% [v1, v2, vi] = obbInteVolume(p1, p2);
% if v1==-1
%     d = -1;
% else
%     d = min(v1,v2)/vi;
% end
% end

function p = sigmoidFunc(v, m, a)
p = 1./(1+exp(a/m*(v-m)));
end

% gnd_objects_id = zeros(length(gnd_objects),1);
% gnd_valid = true(length(gnd_objects),1);
% gnd_xyz = zeros(length(gnd_objects),6);
% gnd_point = cell(length(gnd_objects),1);
% for i = 1:length(gnd_objects)
%     gnd_objects_id(i) = get_object_type({gnd_objects(i).name});
%     x = gnd_objects(i).x_w;
%     ad = abs(round(x(7)/(pi/2))*(pi/2) - x(7));
%     if ad>0.1
%         gnd_valid(i) = false;
%     end
%     out_points_w = gnd_objects(i).out_points_w;
%     gnd_xyz(i,:) = [min(out_points_w, [], 1) max(out_points_w, [], 1)];
%     gnd_point{i} = anno2point(out_points_w, annorule(rule_type==gnd_objects(i).type));
% end
%
%
% hyp_objects_id = zeros(length(hyp_objects),1);
% hyp_xyz = zeros(length(hyp_objects),6);
% hyp_point = zeros(8,3,length(hyp_objects));
% for i = 1:length(hyp_objects)
%     hyp_objects_id(i) = get_object_type({hyp_objects(i).name});
%     out_points_w = hyp_objects(i).out_points_w;
%     hyp_xyz(i,:) = [min(out_points_w, [], 1) max(out_points_w, [], 1)];
% end
% hyp_valid = true(length(hyp_objects),1);
%
% subset_category = [1:12 29];
% XYZ2POINTLIST = [1 2 3; 4 2 3; 1 5 3; 4 5 3; ...
%                  1 2 6; 4 2 6; 1 5 6; 4 5 6];
%
%
% %% check precision
% categorycost = inf*ones(length(subset_category),1);
% categorynumb = zeros(length(subset_category),2);
% hyp_hit_id = zeros(length(hyp_objects),2);
% for k = 1:length(subset_category)
%     oid = subset_category(k);
%     gnd_ids = find(gnd_valid & gnd_objects_id==oid);
%     hyp_ids = find(hyp_objects_id==oid);
%     num_gnd = length(gnd_ids);
%     num_hyp = length(hyp_ids);
%     categorynumb(k,:) = [num_hyp num_gnd];
%     if num_hyp == 0 && num_gnd~=0;
%         categorycost(k) = 860530;
%         continue;
%     elseif num_gnd == 0 && num_hyp~=0
%         categorycost(k) = 860530;
%         continue;
%     elseif num_gnd == 0 && num_hyp == 0
%         categorycost(k) = 0;
%         continue;
%     end
%
%     distance = zeros(length(hyp_ids), length(gnd_ids));
%     gnd_pl = zeros(length(gnd_ids)*8,3);
%     for i = 1:length(gnd_ids)
%         xyz = gnd_xyz(gnd_ids(i),:);
%         gnd_pl(i*8-7:i*8,:) = xyz(XYZ2POINTLIST);
%     end
%     for i = 1:length(hyp_ids)
%         xyz = hyp_xyz(hyp_ids(i),:);
%         hyp_pl = xyz(XYZ2POINTLIST);
%         dist = sqrt(sum((repmat(hyp_pl, length(gnd_ids), 1) - gnd_pl).^2,2));
%         dist = sum(reshape(dist, 8, []),1);
%         distance(i,:) = dist;
%     end
%
%     assigned_id = zeros(num_hyp, 1);
%     assigned_gnd = false(num_gnd,1);
%     assigned = false(num_hyp, 1);
%     assigned_dist = inf*ones(num_hyp, 1);
%     while ~all(assigned) && ~all(assigned_gnd)
%         [B,I] = min(distance, [], 2);
%         [C,J] = min(B);
%         assigned(J) = true;
%         assigned_gnd(I(J)) = true;
%         assigned_id(J) = I(J);
%         assigned_dist(J) = C;
%         hyp_hit_id(hyp_ids(J),:) = [gnd_ids(I(J)) C];
%
%         distance(J,:) = inf;
%         distance(:,I(J)) = inf;
%     end
%
%     max_penalty = 2*max(assigned_dist(~isinf(assigned_dist)));
%     assigned_dist(isinf(assigned_dist)) = max_penalty;
%     miss_gnd_num = max(0, num_gnd-num_hyp);
% %     try
%     categorycost(k) = mean([assigned_dist; max_penalty*ones(miss_gnd_num, 1)]);
% %     catch
% %         fprintf('');
% %     end
% end
%
% cost = mean(categorycost(~isinf(categorycost)));
%
% end
%
% function points = anno2point(out_points, annorule)
% inv_vertex = annorule.inv_vertex;
% for i = 1:size(inv_vertex,1)
%     out_points(end+1,:) = out_points(inv_vertex(i,1),:) ...
%                         + out_points(inv_vertex(i,2),:) ...
%                         - out_points(inv_vertex(i,3),:);
% end
% points = out_points;
% end

author: zouchuhang

repo: LayoutNet


text:
function [rtk,bias,xa,nb]=resamb(rtk)

global glc
bias=zeros(rtk.nx,1); xa=zeros(rtk.nx,1); nb=0;
x=rtk.x; P=rtk.P; opt=rtk.opt; na=rtk.na;
rtk.sol.ratio=0;

if opt.mode<=glc.PMODE_DGNSS||opt.modear==glc.ARMODE_OFF...
        ||opt.thresar(1)<1||opt.ionoopt==glc.IONOOPT_IFLC
    return;
end

% single-difference to double-difference transformation matrix
D=zeros(rtk.nx,rtk.nx);
[rtk,D,nb]=ddmat(rtk,D);
if nb<=0,return;end

% trnasform single-difference ambiguity to double-difference ambiguity
ny=na+nb; D=D(1:ny,:);
y=D*x; Qy=D*P*D';
Qxx=Qy(1:rtk.ib,1:rtk.ib);
Qxb=Qy(1:rtk.ib,rtk.ib+1:end);
Qbb=Qy(rtk.ib+1:end,rtk.ib+1:end);
x_float=y(1:rtk.ib);
bias_float=y(rtk.ib+1:end);

%Is the Q-matrix symmetric?
if ~isequal(Qbb-Qbb'<1E-8,ones(size(Qbb)));
  fprintf('Warning:Variance-covariance matrix of ambiguity is not symmetric!\n');
  nb=0; return;
end

%Is the Q-matrix positive-definite?
if sum(eig(Qbb)>0) ~= size(Qbb,1);
  fprintf ('Warning:Variance-covariance matrix is not positive definite!');
  nb=0; return;
end;

% LAMBDA algorithm for ambiguty resolution
if opt.LAMBDAtype==glc.LAMBDA_ALL

    % resolve all ambiguity
    [bias_fix,sqnorm,~,~,~,~,~]=LAMBDA(bias_float,Qbb,6,'ncands',2,'P0',...
                                            0.001,'MU',1/opt.thresar(1));

    % ratio test
    if sqnorm(1)>0
        rtk.sol.ratio=sqnorm(2)/sqnorm(1);
    else
        rtk.sol.ratio=0;
    end
    if rtk.sol.ratio>999.9,rtk.sol.ratio=999.9;end

    if sqnorm(1)<=0||rtk.sol.ratio>=opt.thresar(1)
        rtk.xa(1:na)=rtk.x(1:na);
        rtk.Pa(1:na,1:na)=rtk.P(1:na,1:na);
        bias(1:nb)=bias_fix(:,1);
        if det(Qbb)
            rtk.xa=x_float-Qxb*cholinv(Qbb)*(bias_float-bias_fix(:,1));
            rtk.Pa=Qxx-Qxb*Qbb^-1*Qxb';
            xa=restamb(rtk,bias,nb,xa);
%             [week,sow]=time2gpst(rtk.sol.time);
%             fprintf('Info:GPS week = %d sow = %.3f,AR success\n',week,sow);
        else
            nb=0;
%             [week,sow]=time2gpst(rtk.sol.time);
%             fprintf('Info:GPS week = %d sow = %.3f,AR failed!!!\n',week,sow);
        end
    else
        nb=0;
%         [week,sow]=time2gpst(rtk.sol.time);
%         fprintf('Info:GPS week = %d sow = %.3f,AR failed!!!\n',week,sow);
    end

elseif opt.LAMBDAtype==glc.LAMBDA_PART
    % resolve partial ambiguity
    [bias_fix,~,Ps,~,~,~,~]=LAMBDA(bias_float,Qbb,5,'ncands',2,'P0',...
                                               opt.thresar(2),'MU',0.6);

    % success rate test
    rtk.sol.ratio=0;

    if Ps>=opt.thresar(2)
        rtk.xa(1:na)=rtk.x(1:na);
        rtk.Pa(1:na,1:na)=rtk.P(1:na,1:na);
        bias(1:nb)=bias_fix(:,1);
        if det(Qbb)
            rtk.xa=x_float-Qxb*cholinv(Qbb)*(bias_float-bias_fix(:,1));
            rtk.Pa=Qxx-Qxb*Qbb^-1*Qxb';
            xa=restamb(rtk,bias,nb,xa);
%             [week,sow]=time2gpst(rtk.sol.time);
%             fprintf('Info:GPS week = %d sow = %.3f,PAR success\n',week,sow);
        else
            nb=0;
%             [week,sow]=time2gpst(rtk.sol.time);
%             fprintf('Info:GPS week = %d sow = %.3f,PAR failed!!!\n',week,sow);
        end
    else
        nb=0;
%         [week,sow]=time2gpst(rtk.sol.time);
%         fprintf('Info:GPS week = %d sow = %.3f,PAR failed!!!\n',week,sow);
    end
end

return

repo: GINav


text:
function [x, resnorm, residual, exitflag, output, lambda] = cplexlsqlin (C, d, Aineq, bineq, Aeq, beq, lb, ub, x0, options)
%%
% Purpose
% Solve constrained least squares problems.
%
% Syntax
%    x = cplexlsqlin(C, d, Aineq, bineq)
%    x = cplexlsqlin(C, d, Aineq, bineq, Aeq, beq)
%    x = cplexlsqlin(C, d, Aineq, bineq, Aeq, beq, lb, ub)
%    x = cplexlsqlin(C, d, Aineq, bineq, Aeq, beq, lb, ub, x0)
%    x = cplexlsqlin(C, d, Aineq, bineq, Aeq, beq, lb, ub, x0, options)
%    x = cplexlsqlin(problem)
%    [x, resnorm] = cplexlsqlin(...)
%    [x, resnorm, residual] = cplexlsqlin(...)
%    [x, resnorm, residual, exitflag] = cplexlsqlin(...)
%    [x, resnorm, residual, exitflag, output] = cplexlsqlin(...)
%    [x, resnorm, residual, exitflag, output, lambda] = cplexlsqlin(...)
%
% Description
% Finds the minimum of a problem specified by
%    min      norm(C*x-d)^2
%    st.      Aineq*x <= bineq
%             Aeq*x    = beq
%             lb <= x <= ub
%
% d, bineq, beq, lb, and ub are column vectors.
% C, Aineq, and Aeq are matrices.
%
% x = cplexlsqlin(C, d, Aineq, bineq) solves the constrained least squares
% problem min norm(C*x-d)^2 such that Aineq*x <= bineq.
%
% x = cplexlsqlin(C, d, Aineq, bineq, Aeq, beq) solves the preceding
% problem with the additional equality constraints Aeq*x = beq.
% If no inequalities exist, set Aineq=[] and bineq=[].
%
% x = cplexlsqlin(C, d, Aineq, bineq, Aeq, beq, lb, ub) defines a set of
% lower and upper bounds on the design variables, x, so that the solution
% is always in the range lb <= x <= ub. If no equalities exist, set Aeq=[]
% and beq=[].
%
% x = cplexlsqlin(C, d, Aineq, bineq, Aeq, beq, lb, ub, x0) sets the
% starting point for the algorithm to x0. If no bounds exist, set lb=[]
% and ub=[].
%
% x = cplexlsqlin(C, d, Aineq, bineq, Aeq, beq, lb, ub, x0, options)
% minimizes with the default optimization options replaced by values in the
% structure options, which can be created using the function cplexoptimset.
% If you do not wish to give an initial point, set x0=[].
%
% x = cplexlsqlin(problem) where problem is a structure.
%
% [x, resnorm] = cplexlsqlin(...) returns the value of the objective
% function at the solution x: resnorm = norm(C*x-d)^2.
%
% [x, resnorm, residual] = cplexlsqlin(...) returns the residual at the
% solution: C*x-d.
%
% [x, resnorm, residual, exitflag] = cplexlsqlin(...) returns a value
% exitflag that describes the exit condition of cplexlsqlin.
%
% [x, resnorm, residual, exitflag, output] = cplexlsqlin(...) returns a
% structure output that contains information about the optimization.
%
% [x, resnorm, residual, exitflag, output, lambda] = cplexlsqlin(...)
% returns a structure lambda whose fields contain the Lagrange multipliers
% at the solution x.
%
% Input Arguments
% problem   Structure containing the following fields:
%           C         Double matrix for objective function
%           d         Double column vector for objective function
%           Aineq     Double matrix for linear inequality constraints
%           bineq     Double column vector for linear inequality
%                     constraints
%           Aeq       Double matrix for linear equality constraints
%           beq       Double column vector for linear equality constraints
%           lb        Double column vector of lower bounds
%           ub        Double column vector of upper bounds
%           x0        Double column vector of initial point of x
%           options   Options structure created with cplexoptimset
%
% Output Arguments
% x         Solution found by the optimization function. If exitflag > 0,
%           then x is a solution; otherwise, x is the value of the
%           optimization routine when it terminated prematurely.
% resnorm   Value of the objective function at the solution x
% residual  Residual at the solution
% exitflag  Integer identifying the reason the optimization algorithm
%           terminated
% output    Structure containing information about the optimization. The
%           fields of the structure are:
%           iterations         Number of iterations
%           algorithm          Optimization algorithm used
%           message            Exit message
%           time               Execution time of the algorithm
%           cplexstatus        Status code of the solution
%           cplexstatusstring  Status string of the solution
% lambda    Structure containing the Lagrange multipliers at the solution x
%           (separated by constraint type). The fields of the structure
%           are:
%           lower              Lower bounds lb
%           upper              Upper bounds ub
%           ineqlin            Linear inequalities
%           eqlin              Linear equalities
%
%
%  See also cplexoptimset
%

% ---------------------------------------------------------------------------
% File: cplexlsqlin.m
% ---------------------------------------------------------------------------
% Licensed Materials - Property of IBM
% 5725-A06 5725-A29 5724-Y48 5724-Y49 5724-Y54 5724-Y55
% Copyright IBM Corporation 2008, 2010. All Rights Reserved.
%
% US Government Users Restricted Rights - Use, duplication or
% disclosure restricted by GSA ADP Schedule Contract with
% IBM Corp.
% ---------------------------------------------------------------------------

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell