text:
seed = 107;
rng(seed);
L = 20;
N = 6;
T = 3000;
W = zeros(N,1,L);
for ni = 0:N-1
    W(ni+1,1,(mod(1:L,N+1)-1)==ni)=1;
end
H = (rand(1,T)<.005);
H(3*L) = 1;
imagesc(squeeze(W));shg
SimpleWHPlot(W,H);
%%
set(gca, 'xscale', 'log', 'ytick', [], 'color', 'none')
set(gca,'color','none','tickdir','out','ticklength', [0.025, 0.025])
savedir = 'C:\Users\emackev\Dropbox (MIT)\SeqNMF\Figures';
saveTitle = 'L1Example';
papersize = [4 2];

[W,H] = seqNMF(helper.reconstruct(W,H)+.0*rand(N,T), 'K', 3, 'L', 2*L,...
    'lambdaL1H', 10, 'lambda', .001, 'maxiter', 150);
SimpleWHPlot(W,H,[],0)
set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
saveas(gcf, fullfile(savedir, [saveTitle '_L1H1.pdf']));

[W,H] = seqNMF(helper.reconstruct(W,H)+.0*rand(N,T), 'K', 3, 'L', 2*L,...
    'lambdaL1W', 100, 'lambda', .001, 'maxiter', 150);
SimpleWHPlot(W,H,[],0)
set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
saveas(gcf, fullfile(savedir, [saveTitle '_L1W1.pdf']));

author: FeeLab

repo: seqNMF

save_path: github-repos/MATLAB/FeeLab-seqNMF


text:
function [s,it]=display(d,tab,it)

if nargin>1
    [s,it]=display(d.algorithm,tab,it);
    return;
end

disp(d.algorithm.name);
s=(['X = ' num2str(size(d.X,1)) 'x' num2str(size(d.X,2))  ]);
% if size(d.X,1)<=10 & size(d.X,2)<=10  & size(d.Y,1)<=10 & size(d.Y,2)<=10
%     if size(d.X,2)==1 & size(d.Y,2)==1
%         if ~isstruct(d.X(1,1))
%             disp('     X     Y'); disp([d.X d.Y]);
%         end
%     else
%         if ~isempty(d.X) & prod(size(d.X))>1  disp('X:'); disp(d.X); end
%         if ~isempty(d.Y) & prod(size(d.X))>1  disp('Y:'); disp(d.Y); end
%     end
%  else
     disp(['data dimensions:'])
    disp([s '  Y = ' num2str(size(d.Y,1)) 'x' num2str(size(d.Y,2)) ]);
% end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


img=imread('peppers.png');

imgS=single(img);
[imgH, imgW, nCh] = size(img);
if nCh==1
    img=repmat(img,[1, 1, 3]);
end

net = load('imagenet-vgg-verydeep-19.mat');
net = vl_simplenn_tidy(net) ;

imgC = imresize(imgS, net.meta.normalization.imageSize(1:2));
imgC = imgC - net.meta.normalization.averageImage;

res = vl_simplenn(net, imgC);

load projMatrix.mat

j=3; % visualizing conv3

featColor=showColorLayer(res(projMatrix(j).numLayer).x, projMatrix(j).meanFeat', projMatrix(j).V);
featColor=imresize(featColor, [imgH, imgW]);

figure, imagesc(featColor)
figure, imagesc(255-featColor);


text:
function [p] = ft_connectivity_plm(inputdata, varargin)

% FT_CONNECTIVITY_PLM computes the phase linearity measurement from a cell array of
% time-domain data, where each cell is an epoch. This implements the metric described
% in Baselice et al. "Phase Linearity Measurement: a novel index for brain functional
% connectivity", IEEE Transactions on Medical Imaging, 2018.
%
% Use as
%   [p] = ft_connectivity_plm(inputdata, ...)
%
% The input data input should be organized as a cell-array, one element for each
% epoch/repetition. Each cell should be a matrix of of nchan x nsamples values.
%
% Additional optional input arguments come as key-value pairs:
%   'bandwidth'	=	scalar, half-bandwidth parameter: the frequency range across which to integrate
%   'fsample'   = sampling frequency, needed to convert bandwidth to number of bins
%
% The output p contains the phase linearity measurement in the [0, 1] interval. It is
% organized as a 3D matrix of Nrpt x Nchan x Nchan dimensions.
%
% See also CONNECTIVITY, FT_CONNECTIVITYANALYSIS

% Copyright (C) 2018, Fabio Baselice, Pierpaolo Sorrentino, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% the sequence of steps is as follows:
%  - Hilbert transformation
%  - multiply with complex conjugate
%  - fft
%  - remove volume conduction component
%  - integrate over bandwidth

% NOTE BY JM: if the user inputs data with different length trials, the fft per trial is going
% to have different frequency resolutions, which is not good. Better to throw an error in that
% case.
fs = ft_getopt(varargin, 'fsample');
B = ft_getopt(varargin, 'bandwidth');
if isempty(fs)
  error('sampling rate is not defined');
end
if isempty(B)
  warning('bandwidth parameter is not defined, assumed 1Hz');
  B=1;
end

nsmp = cellfun('size', inputdata, 2);
assert(all(nsmp==nsmp(1)), 'currently there is no support for input, where the trials are of different length');

nrpt=numel(inputdata);
for k = 1:numel(inputdata)
  inputdata{k} = hilbert(inputdata{k}')';
end
% NOTE by JM: Is it expected that the data has been bandpassfiltered at
% this point? How would this be checked?

nchan=size(inputdata{1},1);
trial_length=size(inputdata{1},2);
ph_min=0.1;        % Eps of Eq.(17) of the manuscript
f=(fs/trial_length)*(0:(trial_length-1));
f_integr=(abs(f)<B) | (abs(f-fs)<B);
p=zeros(nchan, nchan, nrpt);

for ktime=1:nrpt
  for kchan1=1:(nchan-1)
    for kchan2=(kchan1+1):nchan
      temp=fft(inputdata{ktime}(kchan1,:).*conj(inputdata{ktime}(kchan2,:)));    % NOTE BY FB: The inner cycle can be vectorized
      temp(1)=temp(1).*(abs(angle(temp(1)))>ph_min);  % Volume conduction suppression
      temp=(abs(temp)).^2;
      p_temp=sum(temp(f_integr))./sum(temp);
      p(kchan1, kchan2, ktime)=p_temp;
      p(kchan2, kchan1, ktime)=p_temp;
    end
  end
end

p = permute(p, [3 1 2]); % permute to adhere to the conventional matrix shape of the ft_connectivity_* codebase

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [err,time,solver,eqn,node,elem] = afemPoisson(node,elem,pde,bdFlag,option,varargin)

%% Check input arguments
if nargin >=1 && ischar(node)
    option.elemType = node;
    clear node
end
if ~exist('node','var') || ~exist('elem','var')
    % default mesh: Lshape
    [node,elem] = squaremesh([-1,1,-1,1],1);
    [node,elem] = delmesh(node,elem,'x>0 & y<0');
end
if ~exist('option','var'), option = []; end
if ~exist('pde','var')
    pde = Lshapedata;                          % default data
end
if ~exist('bdFlag','var')
    bdFlag = setboundary(node,elem,'Dirichlet');
end

%% Parameters
option = afemoption(option,2);
maxIt = option.maxIt;
refType = option.refType;
elemType = option.elemType;
theta = option.theta;

%% Initialize err
errL2 = zeros(maxIt,1);   errH1 = zeros(maxIt,1); erreta = zeros(maxIt,1);
erruIuh = zeros(maxIt,1); errMax = zeros(maxIt,1);
errTime = zeros(maxIt,1); solverTime = zeros(maxIt,1);
assembleTime = zeros(maxIt,1); meshTime = zeros(maxIt,1);
itStep = zeros(maxIt,1);  stopErr = zeros(maxIt,1); flag = zeros(maxIt,1);
N = zeros(maxIt,1);

%% Generate an initial mesh
for k = 1:option.L0
    if strcmp(refType,'red')
        [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
    elseif strcmp(refType,'bisect')
        [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
    end
end

%%  Adaptive Finite Element Method
% *SOLVE* -> *ESTIMATE* -> *MARK* -> *REFINE*
for k = 1:maxIt
    % Step 1: SOLVE
    switch elemType
        case 'P1'     % piecewise linear function P1 element
            [u,Du,eqn,info] = Poisson(node,elem,pde,bdFlag,option);
        case 'CR'     % piecewise linear function CR element
            [u,Du,eqn,info] = PoissonCR(node,elem,pde,bdFlag,option);
        case 'P2'     % piecewise quadratic function
            [u,Du,eqn,info] = PoissonP2(node,elem,pde,bdFlag,option);
        case 'WG'     % weak Galerkin element
            [u,Du,eqn,info] = PoissonWG(node,elem,pde,bdFlag,option);
    end
    % compute error
    tic;
    if isfield(pde,'Du')
        if ~isfield(pde,'d')
            pde.d = [];
        end
        if ~isempty(Du)
            errH1(k) = getH1error(node,elem,pde.Du,Du,pde.d);
        else
            errH1(k) = getH1error(node,elem,pde.Du,u,pde.d);
        end
    end
    if isfield(pde,'exactu')
        errL2(k) = getL2error(node,elem,pde.exactu,u);
        % interpolation
        if strcmp(elemType,'P1')
            uI = Lagrangeinterpolate(pde.exactu,node,elem);
        else
            uI = Lagrangeinterpolate(pde.exactu,node,elem,elemType,eqn.edge);
        end
        erruIuh(k) = sqrt((u-uI)'*eqn.A*(u-uI));
        errMax(k) = max(abs(u-uI));
    end
    errTime(k) = toc;
    % record time
    solverTime(k) = info.solverTime;
    assembleTime(k) = info.assembleTime;
    if option.printlevel>1
        fprintf('Time to compute the error %4.2g s \n H1 err %4.2g    L2err %4.2g \n',...
            errTime(k),errH1(k), errL2(k));
    end
    % record solver information
    itStep(k) = info.itStep;
    stopErr(k) = info.stopErr;
    flag(k) = info.flag;
    % plot
    N(k) = size(node,1);
    if option.plotflag && N(k) < 2e3 % show mesh and solution for small size
       figure(1);  showresult(node,elem,u);
    end
    % Step 2: ESTIMATE
    switch option.estType
        case 'recovery' % recovery type
            eta = estimaterecovery(node,elem,u);
        case 'residual' % residual type
            switch elemType
                case 'P1'
                    eta = estimateresidual(node,elem,u,pde,bdFlag);
                case 'WG'
                    eta = estimateresidualWG(node,elem,u,Du,pde);
            end
    end
    erreta(k) = sqrt(sum(eta.^2));
    % Step 3: MARK
    switch option.markType
        case 'L2'
            markedElem = mark(elem,eta,theta);
        case 'MAX'
            markedElem = mark(elem,eta,theta,'MAX');
    end
    % Step 4: REFINE
    if N(k) > option.maxN
        break;
    end
    [node,elem,bdFlag,HB,tree] = bisect(node,elem,markedElem,bdFlag); %#ok<ASGLU>
    % Step 4.2: COARSEN
    if option.coarsenflag
        eta = eleminterpolate(eta,tree);
        markedElem = mark(elem,eta,0.25*theta,'COARSEN');
        [node,elem,bdFlag] = coarsen(node,elem,markedElem,bdFlag);
    end
end

%% Plot convergence rates
if option.rateflag
    figure;
    set(gcf,'Units','normal');
    set(gcf,'Position',[0.25,0.25,0.55,0.4]);
    subplot(1,2,1)
    showrate2(N(1:k),errH1(1:k),10,'-*','||Du-Du_h||',...
              N(1:k),errL2(1:k),10,'k-+','||u-u_h||');
    subplot(1,2,2)
    showrate2(N(1:k),erruIuh(1:k),10,'m-+','||Du_I-Du_h||',...
              N(1:k),errMax(1:k),10,'r-*','||u_I-u_h||_{\infty}');
end

%% Output
err = struct('N',N(1:k),'H1',errH1(1:k),'L2',errL2(1:k),...
             'uIuhH1',erruIuh(1:k),'uIuhMax',errMax(1:k),'eta',erreta(1:k));
time = struct('N',N(1:k),'err',errTime(1:k),'solver',solverTime(1:k), ...
              'assmble',assembleTime(1:k),'mesh',meshTime(1:k));
solver = struct('N',N(1:k),'itStep',itStep(1:k),'time',solverTime(1:k),...
                'stopErr',stopErr(1:k),'flag',flag(1:k));

%% Display error
ts = zeros(k,3); ts = char(ts);
display('#Dof   ||u-u_h||     ||Du-Du_h||   ||DuI-Du_h||  ||uI-u_h||_{max}    eta');
display([num2str(err.N) ts num2str(err.L2,'%0.5e') ts num2str(err.H1,'%0.5e')...
         ts num2str(err.uIuhH1,'%0.5e') ts num2str(err.uIuhMax,'%0.5e') ...
         ts num2str(err.eta,'%0.5e')]);

repo: Aero-matlab


text:
function [varargout] = slick(coords)
    % turned into MATLAB from C by Celso G Reyes
    % accepts input of [DipDir, Dip, Rake]
    % returns different items, based on varargout.

    % answer is CLOSE to original, but does not match.  Sig Figs? Technique?

    TORADS = 57.29577951;

    % COORDINATES ARE EAST,NORTH,UP
    % this version does no statistics
    % and therefore makes no plot

    assert(size(coords,2) == 3); % [ ddir, dip, rake ] as ENU

    % ddir = zeros(MAXDATA,1);  % dip direction for data
    % dip = zeros(MAXDATA,1);   % dip of data
    % rake = zeros(MAXDATA,1);  % rake of data
    % amat = zeros(MAX3,5);  % coefficient matrix for normal equation
    % stress = zeros(6,1);  % stress tensor in vector form, element order is:
    % xx,xy,xz,yy,yz,zz
    % slick_vec_el_vec  slickenside vector elements vector
    % norm  % storage of n1,n2,n3
    char name(20);      % output file name
    %FILE *fpin;   % input file pointer
    %FILE *fpout;  % output file pointer
    %FILE *fplot;  % plot file pointer
    % sigma = 0;  % for use with leasq subr
    % a2i = zeros(5,5);  % to get covariance mtrix
    line='';  % character line
    t = zeros(3,1);  % shear stress vector
    % iso = 0;  % isotropic stress mag
    %angavg = 0; angstd = 0;  % average and standard deviation of fit angle
    %isoavg = 0; isostd = 0;  % same for isotropic stress size
    magavg = 0; magstd = 0;  % same for tangential stress size
    % tf = zeros(3,1), tnorm;  % full traction vector
    % and normal traction

    % get file pointers
    %{
    -- argc;
    ++argv;
    if argc == 0
        printf("usage: slick data_file\n");
        return;
    end
    fpin = fopen(argv,"r");
    if fpin==NULL
        printf("unable to open %s.\n",argv);
        return;
    end
    fprintf(name,"%s.oput",argv);
    %}
    % read and write comment line from data file to output file
    % fgets(line,80,fpin);
    line = 'Inversion data';
    fpout = string(line);
    ddir = coords(:,1);
    dip = coords(:,2);
    rake = coords(:,3);
    % loop to get data and make up equation
    for i=1:size(coords,1)
        %ddir = coords(i,1); dip = coords(i,2); rake = coords(i,3);
        j = 3*i;%?

        z = ddir(i)/TORADS;
        z2 = dip(i)/TORADS;
        z3 = rake(i)/TORADS;

        %n1 to n3 are normal vector elements
        n1 = sin(z)*sin(z2);  % normal vector to fault plane
        n2 = cos(z)*sin(z2);
        n3 = cos(z2);

        norm(i,1:3) = [n1 n2 n3];

        % slickenside vector calculation
        slick_vec_el_vec(j,1)= -cos(z3)*cos(z) - sin(z3)*sin(z)*cos(z2);
        slick_vec_el_vec(j+1,1)= cos(z3)*sin(z) - sin(z3)*cos(z)*cos(z2);
        slick_vec_el_vec(j+2,1)= sin(z3)*sin(z2);

        % find the matrix elements
        amat(j:j+2, 1:5) = [...
            n1-n1*n1*n1+n1*n3*n3,  n2-2.*n1*n1*n2, n3-2.*n1*n1*n3, -n1*n2*n2+n1*n3*n3,    -2.*n1*n2*n3 ;...
            -n2*n1*n1+n2*n3*n3,    n1-2.*n1*n2*n2, -2.*n1*n2*n3,   n2-n2*n2*n2+n2*n3*n3,  n3-2.*n2*n2*n3;...
            -n3*n1*n1-n3+n3*n3*n3, -2.*n1*n2*n3,   n1-2.*n1*n3*n3, -n3*n2*n2-n3+n3*n3*n3, n2-2.*n2*n3*n3];

        % check to see if all possible data has been read
    end  % end of data read loop

	% solve equations via linear least squares
	[stress, sigma]=leasq(amat,slick_vec_el_vec);
	% correct zz element by using trace = 0
	stress(6)= -(stress(1)+stress(4));

	% put stress tensor into tensor form
    strten = [  stress(1),    stress(2),   stress(3) ;
                stress(2),    stress(4),   stress(5) ;
                stress(3),    stress(5),   stress(6)];

    %fpout(end+1)=sprintf("\nCOORDINATES ARE EAST,NORTH,UP.");
    %fpout(end+1)=sprintf("stress tensor is:");
    for i=1:3
    %    fpout(end+1)=sprintf("%g  %g  %g  ",strten(i,:));
    end

	% find  eigenvalues and eigenvectors
    [vecs,lam] = eig(strten,"vector"); %LAM is eigenvalues, VECS is eigenvectors
    % eigen(strten,lam,vecs);
    %fpout(end+1)=sprintf("eigenvalue   vector: E,N,UP,direction,plunge");
    for i = 3:-1:1
        [v_direction(i), v_plunge(i)] = dirplg(vecs(1,i),vecs(2,i),vecs(3,i));
        %[z, z2] = dirplg(vecs(1,i),vecs(2,i),vecs(3,i));
        %fpout(end+1)=sprintf("%g  ",lam(i)) +...
        %    sprintf("%g  %g  %g  ",vecs(:,i)) +...
        %    sprintf("%f  %f",z,z2);
    end
    %fpout(end+1)=sprintf("variance= %g",sigma);

	% order eigenvalues and compute phi
    %lam=sort(lam,'descend');

    if lam(1) ~= lam(3)
        phi = (lam(2)-lam(3)) / (lam(1)-lam(3));
        %fpout(end+1)=sprintf("phi value= %g",phi);
    else
        phi=nan;
    end
	% output data and fit angle

	angavg = 0.;
	angstd = 0.;
	isoavg = 0.;
	isostd = 0.;
	iso = 0.;
	%fpout(end+1)=sprintf("\ndip direction, dip, rake, fit angle, mag tau");
    nobs=size(coords,1);
    for i=1:nobs %from 0

        for j= 1 : 3 %from 0  % compute shear traction
            t(j)=0;
            tf(j)=0;
            myt(j) = sum(amat(3*(i-1)+j) * stress(1:5));
            for k=1:5
                t(j) = t(j)+ amat(3*(i-1)+j,k) * stress(k);
            end
            for k=1:3
                tf(j) = tf(j) + strten(j,k) * norm(i,k);
            end
        end
        tnorm = 0;
        for k=1:3
            tnorm = tnorm + tf(k) * norm(i,k);
        end
        % find angle between t and slickenside
        z = 0.;
        for j=1:3
            z = z + t(j)*slick_vec_el_vec(3*(i-1)+j);
        end
        z2 = 0.;
        for j=1:3
            z2 = z2 + t(j)*t(j);
        end
        z2 = sqrt(z2);
        z3 = 0.;
        for j=1:3
            z3 = z3 + slick_vec_el_vec(3*(i-1)+j)*slick_vec_el_vec(3*(i-1)+j);
        end
        z3 = sqrt(z3);
        z = z/(z2*z3);
        z = acos(z)*TORADS;
        angavg = angavg + z;
        angstd = angstd + z*z;
        z3= (z2/(-0.8)) - tnorm;
        iso = iso + abs(tnorm);
        isoavg = isoavg +z3;
        isostd = isostd + z3*z3;
        magavg = magavg +z2;
        magstd = magstd + z2*z2;
        %fpout(end+1)=sprintf("%7.1f  %7.1f  %7.1f  %7.1f %7.2f", ddir(i),dip(i),rake(i),z,z2);
    end
    z3 = nobs-1;
    angstd = angstd-(angavg*angavg/nobs);
    angstd = angstd/z3;
    angstd = sqrt(angstd);
    angavg = angavg/nobs;

    isostd = isostd-(isoavg*isoavg/nobs);
    isostd = isostd/z3;
    isostd = sqrt(isostd);
    isoavg = isoavg/nobs;
    iso = iso / nobs;
    isoavg = isoavg / iso;
    isostd = isostd / iso;

    magstd = magstd-(magavg*magavg/nobs);
    magstd = magstd/z3;
    magstd = sqrt(magstd);
    magavg = magavg/nobs;

    %fpout(end+1)=sprintf("fit angle mean= %f standard deviation= %f",angavg,angstd);
    %fpout(end+1)=sprintf("for f=0.8 I= %f , std. dev.= %f D norm= %f", isoavg,isostd,iso);
    %fpout(end+1)=sprintf("avg tau= %f , std. dev.= %f",magavg,magstd);

    if nargout==1
        fpout(end+1)=sprintf("fit angle mean= %f standard deviation= %f",angavg,angstd);
        fpout(end+1)=sprintf("for f=0.8 I= %f , std. dev.= %f D norm= %f", isoavg,isostd,iso);
        fpout(end+1)=sprintf("avg tau= %f , std. dev.= %f",magavg,magstd);

    varargout = {strjoin(fpout,newline)};
    elseif nargout == 5
        varargout={angavg, angstd, magstd/magavg, magavg, magstd}; %[fBeta2, fStdBeta2, fTauFit2, fAvgTau2, fStdTau2]
    elseif nargout == 9
        varargout=fastoutput();
    end

    function output = fastoutput()
        %%line 1
        % variance
        output = {... line 1 of output file
            sigma,... variance
            stress,...      stress tensor upper triangle
            ... line 2 of output file
            phi,...
            round(v_direction(1),1),...
            round(v_plunge(1),1),...
            round(v_direction(2),1),...
            round(v_plunge(2),1),...
            round(v_direction(3),1),...
            round(v_plunge(3),1)...
            };
    end

end

function [pdir, pplg]  = dirplg(e,n,u)
    % dirplb to find direction and plunge of a vector
    % double e,n,u; /* the vector in east,north,up coordinates
    % double *pdir,*pplg are pointers to the direction in east of north
    % and the plunge down the direction

    TORADS = 57.29577951;

    z=e*e+n*n;
    z=sqrt(z);
    pplg=atan2(-u,z) * TORADS;
    if pplg<0
        pplg= -pplg;
        e= -e;
        n= -n;
    end
    pdir=atan2(e,n) * TORADS;
end

function [x, psis] = leasq(a,b)
% /* finds the least squares solution of ax=b */
%{
double a[]; /* the coefficients matrix with n rows and m columns */
double x[]; /* the solution vector of length m */
double b[]; /* the constant vector of length n */
double a2[]; /* a square matrix of size m for internal use */
double c[]; /* vector of length m for internal use */
double *psis; /*pointer to the variance */

/* steps 1 a2= a transpose a */
/*       2 c= a transpose b */
/*       3 solve a2x=c by gaussian elimination */
%}

	% a2=atransa(a,a2); % computes b=a transpose*a
    a2 = a' * a;
    c = a' * b;
    x = a2 \ c;
	% x = gaus(a2, c); % solves ax=b for x by gaussian elimination
	psis = sigsq(a,x,b);
end
%{
function x = gaus(a,b)
    % /* solves ax=b for x by gaussian elimination */

%double a[];  /* a square matrix of size m */
%double b[];  /* a vector of length m */
%double x[];  /* a vector of length m */

    m=length(a);
    %  take care of special cases */
    if m<2
        x=0;
        if m==1
            x=b/a;
        end
        return
    end

	for i=1:m % (i=0;i<m;++i)   %  /* loop for each pivot */
		for i2=i+1:m %  /* loop for each row below a pivot */

			% /* see if element below pivot is 0 */
			if(a(i2,i)==0.)
                continue
            end

			% /* if element below pivot > pivot flop rows  */
			if(abs(a(i2,i))>abs(a(i,i)))
				hold=b[i];
				b[i]=b[i2];
				b[i2]=hold;
				for(i3=i;i3<m;++i3)
					hold=a(i,i3);
					a(i,i3)=a(i2,i3);
					a(i2,i3)=hold;
                end
            end

			% /* do the elimination */
			fact=a(i2,i)/a(i,i);
			a(i2,i)=0.;
			for(i3=i+1;i3<m;++i3)
                a(i2,i3)=a(i2,i3)-fact*a(i,i3);
            end
			b[i2]=b[i2]-fact*b[i];

        end
    end

	/* solve the equations */
	x[m-1]=b[m-1]/a[m*m-1];
	for(i=m-2;i> -1;--i){
		d=b[i];
		for(i2=i+1;i2<m;++i2)
            d=d-x[i2]*a(i,i2);
        end
		x[i]=d/a(i,i);
	}
	return;
}

end
%}

function psis = sigsq(a,x,b)
% SIGSQ  computes the variance of a single observation */
% double a[];    /* matrix of n rows and m columns */%
% double b[];    /* data vector length n*/
% double x[];    /* solution vector length m */
% double *psis;  /* where to put answer */
% short m,n;     /* see above */

	% double y,z,z2;     /* sum variables */

    [n,m] = size(a);
    allY = a * x; % nx1
    z=sum(b-allY).^2;

	if n ~= m
		z2 = n - m;
		z = z/z2;
    end
	psis=z;
end

author: CelsoReyes


text:
function write_filter(weight, file_name);
[filter_len, N_ch] = size(weight);
FID = fopen(file_name,'w');

fprintf(FID, '%f\n',filter_len);
fprintf(FID, '%f\n',N_ch);

for i=1:N_ch,
    for j=1:filter_len
        fprintf(FID, '%f ',weight(j,i));
    end
    fprintf(FID,'\n');
end
fclose(FID);

author: singaxiong

repo: SignalGraph


text:
function a = clement1_inverse ( n )

%*****************************************************************************80
%
%% CLEMENT1_INVERSE returns the inverse of the CLEMENT1 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.  N must not be odd%
%
%    Output, real A(N,N), the matrix.
%
  if ( mod ( n, 2 ) == 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CLEMENT1_INVERSE - Fatal error!\n' );
    fprintf ( 1, '  The matrix is singular for odd N.\n' );
    error ( 'CLEMENT1_INVERSE - Fatal error!' );
  end

  a = zeros ( n, n );

  for i = 1 : n

    if ( mod ( i, 2 ) == 1 )

      for j = i : 2 : n-1

        if ( j == i )
          prod = 1.0 / sqrt ( ( j * ( n - j ) ) );
        else
          prod = - prod ...
            * sqrt ( ( j - 1 ) * ( n + 1 - j ) ) ...
            / sqrt ( j * ( n - j ) );
        end

        a(i,j+1) = prod;
        a(j+1,i) = prod;

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Disparity Filtering Demo
% In this tutorial you will learn how to use the disparity map post-filtering
% to improve the results of |cv.StereoBM| and |cv.StereoSGBM| algorithms.
%
% Sources:
%
% * <https://github.com/opencv/opencv_contrib/blob/3.1.0/modules/ximgproc/samples/disparity_filtering.cpp>
% * <https://docs.opencv.org/3.1.0/d3/d14/tutorial_ximgproc_disparity_filtering.html>
%

%% Introduction
% Stereo matching algorithms, especially highly-optimized ones that are
% intended for real-time processing on CPU, tend to make quite a few errors on
% challenging sequences. These errors are usually concentrated in uniform
% texture-less areas, half-occlusions and regions near depth discontinuities.
% One way of dealing with stereo-matching errors is to use various techniques
% of detecting potentially inaccurate disparity values and invalidate them,
% therefore making the disparity map semi-sparse. Several such techniques are
% already implemented in the StereoBM and StereoSGBM algorithms. Another way
% would be to use some kind of filtering procedure to align the disparity map
% edges with those of the source image and to propagate the disparity values
% from high- to low-confidence regions like half-occlusions. Recent advances
% in edge-aware filtering have enabled performing such post-filtering under
% the constraints of real-time processing on CPU.

%%
% The provided example has several options that yield different trade-offs
% between the speed and the quality of the resulting disparity map. Both the
% speed and the quality are measured if the user has provided the ground-truth
% disparity map. In this tutorial we will take a detailed look at the default
% pipeline, that was designed to provide the best possible quality under the
% constraints of real-time processing on CPU.

%% Options

% left/right views of the stereopair
%left_im = 'ambush_5_left.jpg';
%right_im = 'ambush_5_right.jpg';
left_im = fullfile(mexopencv.root(),'test','aloeL.jpg');
right_im = fullfile(mexopencv.root(),'test','aloeR.jpg');

% optional ground-truth disparity (MPI-Sintel or Middlebury format),
% set it to empty string if not available
%GT_path = '';
GT_path = fullfile(mexopencv.root(),'test','aloeGT.png');

% stereo matching method: 'bm' or 'sgbm'
algo = 'bm';

% used post-filtering: 'wls_conf' or 'wls_no_conf'
filt = 'wls_conf';

% force stereo matching on full-sized views to improve quality
no_downscale = false;

% parameter of stereo matching: max disparity and window size
max_disp = 160;
wsize = -1;  % -1 to get appropriate default value

% parameter of post-filtering: wls_lambda and wls_sigma
lambda = 8000.0;
sigma = 1.5;

% coefficient used to scale disparity map visualizations
vis_mult = 1.0;

%%
% check user-provided values

algo = validatestring(algo, {'bm', 'sgbm'});
filt = validatestring(filt, {'wls_conf', 'wls_no_conf'});

if wsize < 0
    if strcmp(algo, 'sgbm')
        % default window size for SGBM
        wsize = 3;
    elseif ~no_downscale && strcmp(algo, 'bm') && strcmp(filt, 'wls_conf')
        % default window size for BM on downscaled views
        % (downscaling is performed only for wls_conf)
        wsize = 7;
    else
        % default window size for BM on full-sized views
        wsize = 15;
    end
end
assert(wsize>0 && mod(wsize,2)==1, ...
    'Incorrect window size value: must be positive and odd');

assert(max_disp>0 && mod(max_disp,16)==0, ...
    'Incorrect max disparity value: must be positive and divisible by 16');

%% Source Stereoscopic Image
% We start by loading the source stereopair. For this tutorial we will take a
% somewhat challenging example from the MPI-Sintel dataset with a lot of
% texture-less regions.

left = cv.imread(left_im, 'Color',true);
right = cv.imread(right_im, 'Color',true);
assert(~isempty(left) && ~isempty(right), 'Cannot read image files');

%%
% load ground-truth disparity if supplied

if ~isempty(GT_path)
    GT_disp = cv.DisparityWLSFilter.readGT(GT_path);
    assert(~isempty(GT_disp), 'Cannot read ground truth image file');
else
    GT_disp = [];
end

%% Prepare the views for matching
% We perform downscaling of the views to speed-up the matching stage at the
% cost of minor quality degradation. To get the best possible quality
% downscaling should be avoided.

if strcmp(filt, 'wls_conf') && ~no_downscale
    % downscale the views to speed-up the matching stage, as we will need to
    % compute both left and right disparity maps for confidence map computation
    max_disp = max_disp / 2;
    if mod(max_disp,16)~=0
        max_disp = max_disp + 16-mod(max_disp,16);
    end
    left_for_matcher = cv.resize(left, 0.5, 0.5);
    right_for_matcher = cv.resize(right, 0.5, 0.5);
else
    left_for_matcher = left;
    right_for_matcher = right;
end

if strcmp(algo, 'bm')
    left_for_matcher = cv.cvtColor(left_for_matcher,  'RGB2GRAY');
    right_for_matcher = cv.cvtColor(right_for_matcher, 'RGB2GRAY');
end

%% Process
% We are using StereoBM for faster processing. If speed is not critical,
% though, StereoSGBM would provide better quality. The filter instance is
% created by providing the StereoMatcher instance that we intend to use.
% Another matcher instance is returned by the createRightMatcher function.
% These two matcher instances are then used to compute disparity maps both for
% the left and right views, that are required by the filter.
%
% Next, disparity maps computed by the respective matcher instances, as well
% as the source left view are passed to the filter. Note that we are using the
% original non-downscaled view to guide the filtering process. The disparity
% map is automatically upscaled in an edge-aware fashion to match the original
% view resolution. The result is stored in filtered_disp.

if strcmp(filt, 'wls_conf')
    % filtering with confidence (significantly better quality than wls_no_conf)

    % Create the matching instances
    if strcmp(algo, 'bm')
        left_matcher = cv.StereoBM('NumDisparities',max_disp, 'BlockSize',wsize);
    elseif strcmp(algo, 'sgbm')
        left_matcher  = cv.StereoSGBM('NumDisparities',max_disp, 'BlockSize',wsize, ...
            'MinDisparity',0);
        left_matcher.P1 = 24*wsize*wsize;
        left_matcher.P2 = 96*wsize*wsize;
        left_matcher.PreFilterCap = 63;
        left_matcher.Mode = 'SGBM3Way';
    end
    right_matcher = cv.DisparityWLSFilter.createRightMatcher(left_matcher);

    % Perform matching
    fprintf('Matching time: '); tic
    left_disp = left_matcher.compute(left_for_matcher, right_for_matcher);
    right_disp = right_matcher.compute(right_for_matcher, left_for_matcher);
    toc

    % Create the filter instance
    wls_filter = cv.DisparityWLSFilter(left_matcher);
    wls_filter.Lambda = lambda;
    wls_filter.SigmaColor = sigma;

    % Perform filtering
    fprintf('Filtering time: '); tic
    filtered_disp = wls_filter.filter(left_disp, right_disp, left);
    toc

    % Get the confidence map that was used in the last filter call
    conf_map = wls_filter.getConfidenceMap();

    % Get the ROI that was used in the last filter call
    ROI = wls_filter.getROI();
    if ~no_downscale
        % upscale raw disparity and ROI back for a proper comparison:
        left_disp = 2.0 * cv.resize(left_disp, 2.0, 2.0);
        ROI = 2 * ROI;
    end

elseif strcmp(filt, 'wls_no_conf')
    % There is no convenience function for the case of filtering with no
    % confidence, so we will need to set the ROI and matcher parameters manually

    % Create the matching instance
    if strcmp(algo, 'bm')
        matcher = cv.StereoBM('NumDisparities',max_disp, 'BlockSize',wsize);
        matcher.TextureThreshold = 0;
        matcher.UniquenessRatio = 0;
        ddr = 0.33;
    elseif strcmp(algo, 'sgbm')
        matcher = cv.StereoSGBM('NumDisparities',max_disp, 'BlockSize',wsize, ...
            'MinDisparity',0);
        matcher.UniquenessRatio = 0;
        matcher.Disp12MaxDiff = 1000000;
        matcher.SpeckleWindowSize = 0;
        matcher.P1 = 24*wsize*wsize;
        matcher.P2 = 96*wsize*wsize;
        matcher.Mode = 'SGBM3Way';
        ddr = 0.5;
    end

    % Perform matching
    fprintf('Matching time: '); tic
    left_disp = matcher.compute(left_for_matcher, right_for_matcher);
    toc

    % Create the filter instance
    wls_filter = cv.DisparityWLSFilter(false);
    wls_filter.Lambda = lambda;
    wls_filter.SigmaColor = sigma;
    wls_filter.DepthDiscontinuityRadius = ceil(ddr*wsize);

    % manually compute ROI
    xmin = matcher.MinDisparity + matcher.NumDisparities - 1 + matcher.BlockSize/2;
    xmax = size(left_for_matcher,2) + matcher.MinDisparity - matcher.BlockSize/2;
    ymin = matcher.BlockSize/2;
    ymax = size(left_for_matcher,1) - matcher.BlockSize/2;
    ROI = [xmin, ymin, xmax - xmin, ymax - ymin];

    % Perform filtering
    fprintf('Filtering time: '); tic
    filtered_disp = wls_filter.filter(left_disp, [], left, 'ROI',ROI);
    toc

    % no confidence map
    conf_map = [];

end

%% Stats
% We compare against the ground-truth disparity

if ~isempty(GT_disp)
    MSE_before = cv.DisparityWLSFilter.computeMSE(GT_disp, left_disp, 'ROI',ROI);
    MSE_after = cv.DisparityWLSFilter.computeMSE(GT_disp, filtered_disp, 'ROI',ROI);
    percent_bad_before = cv.DisparityWLSFilter.computeBadPixelPercent(GT_disp, left_disp, 'ROI',ROI);
    percent_bad_after = cv.DisparityWLSFilter.computeBadPixelPercent(GT_disp, filtered_disp, 'ROI',ROI);

    fprintf('MSE before filtering: %.5f\n', MSE_before);
    fprintf('MSE after filtering: %.5f\n', MSE_after);
    fprintf('Percent of bad pixels before filtering: %.3f\n', percent_bad_before);
    fprintf('Percent of bad pixels after filtering: %.3f\n', percent_bad_after);
end

%% Visualize the disparity maps
% We use a convenience function getDisparityVis to visualize the disparity
% maps. The second parameter defines the contrast (all disparity values are
% scaled by this value in the visualization).
%
% Compare the raw result of StereoBM against the result of StereoBM on
% downscaled views with post-filtering

if ~isempty(GT_disp)
    GT_disp_vis = cv.DisparityWLSFilter.getDisparityVis(GT_disp, 'Scale',vis_mult);
else
    GT_disp_vis = [];
end
raw_disp_vis = cv.DisparityWLSFilter.getDisparityVis(left_disp, 'Scale',vis_mult);
filtered_disp_vis = cv.DisparityWLSFilter.getDisparityVis(filtered_disp, 'Scale',vis_mult);

% left view of the stereopair
subplot(231), imshow(left), title('left')
% right view of the stereopair
subplot(232), imshow(right), title('right')
% ground-truth disparity
subplot(233), imshow(GT_disp_vis), title('ground-truth disparity')
% disparity map before filtering
subplot(234), imshow(raw_disp_vis), title('raw disparity')
% resulting filtered disparity map
subplot(235), imshow(filtered_disp_vis), title('filtered disparity')
% confidence map used in filtering
subplot(236), imshow(conf_map), title('confidence map')

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function writedb(c,dbOut,varargin)

   %WRITEDB write antelope database tables
   %  WRITEDB(C,DBOUT) writes an arrival table into the database DBOUT. The
   %  arrival contains the station and channel names. The phase arrival time
   %  is set to the trigger field in the correlation object. If the cluster
   %  field is filled in C, then the phase names are assigned as cXX where XX
   %  is the cluster number. If the cluster field is not filled, a default
   %  phase name of P is assigned. If the table DBOUT.arrival already exists,
   %  new entries are appended to it. Duplicate entries to the database will
   %  NOT be skipped.
   %
   %  WRITEDB(C,DBOUT,PHASENAME) Assign an explicit PHASENAME. PHASENAME is
   %  a cell array containing character phase names. The cell array must be
   %  of length 1 or of length N, where N is the number of traces in the
   %  corelation object.
   %
   %  WRITEDB(C,DBOUT,'detection') Writes a detection table instead of an
   %  arrival table. The filter field is left null by default. Rows for
   %  detection ON and OFF states are not included. On state='D' rows are
   %  written.
   %
   %  WRITEDB(C,DBOUT,'detection','filterString') Fills the filter field on
   %  the detection table with the string given by filterString.
   %
   %  WRITEDB checks for the Antelope toolbox for Matlab. If it does not
   %  exist, WRITEDB writes a simple text file DBOUT.txt which contains
   %  station, channel, arrival times and phase names.

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   % READ & CHECK ARGUMENTS
   if numel(varargin)>2
      error('Wrong number of inputs');
   end;

   if ~ischar(dbOut)
      disp('Second input parameter must be a text string');
   end

   if ~isempty(c.lags)
      disp('Note: Time corrections from LAG field have not been applied to traces. Consider using ADJUSTTRIG prior to writing out');
      disp(' ');
   end;

   % SELECT TABLE TYPE
   ISDETECTION = 0;
   if numel(varargin)>=1 || ~admin.antelope_exists
      if strcmpi(varargin{1},'detection')
         ISDETECTION = 1;
      end
   end

   % ASSIGN ARGUMENTS FOR ARRIVAL TABLE
   if ~ISDETECTION
      if numel(varargin)==1
         phaseList = varargin{1};
         if numel(phaseList)== c.ntraces
            phaseList = phaseList;
         elseif numel(phaseList)~= c.ntraces && numel(phaseList)~=1
            error('phase list argument must be cell array of length 1 or N, where N is the number of traces in the correlation object');
         elseif numel(phaseList)==1
            if ischar(phaseList)
               phaseList = {phaseList};
            end
            phaseList = repmat(phaseList, c.ntraces,1);
         end
      elseif ~isempty(c.clust)
         disp('assigning phase names based on cluster field');
         clust = c.clust;
         phaseList = cell(size(clust));
         for n = 1:numel(clust)
            phaseList(n) = {['c' num2str(clust(n),'%02.0f')]};
         end
      else
         disp('assigning default phase name "P"');
         phaseList = repmat({'P'},c.ntraces,1);
      end;

   end

   % ASSIGN ARGUMENTS FOR DETECTION TABLE
   if ISDETECTION
      if numel(varargin)>=2
         if ischar(varargin{2})
            filterString = varargin(2);
         elseif iscell(varargin{2}) && (numel(varargin{2})==1)
            filterString = varargin{2};
         end
      else
         filterString = 'not_specified';
      end
   end

   % WRITE ARRIVAL TABLE OR TEXT FILE
   if ~ISDETECTION
      trig = c.trig;
      sta = c.stations; %sta = sta{1};
      chan = c.channels; %chan = chan{1};
      if admin.antelope_exists
         % TODO: check for database table
         disp(['Writing database table ' dbOut '.arrival']);
         db = dbopen(dbOut,'r+');
         db = dblookup(db,'','arrival','','');
         nMax = length(trig);
         disp('     ');
         for n = 1:nMax
            arid = dbnextid( db,'arid');
            db.record = dbaddv(db,'sta',sta{n},'chan',chan{n},'time',datenum2epoch(trig(n)),'iphase',phaseList{n},'auth','corr/writedb');
            fprintf('\b\b\b\b\b\b%5.0f%%',n/nMax*100);
         end
         fprintf('\n');
         dbclose(db);
      else
         % write out a simple text file
         disp(['Antelope toolbox not found. Writing text file ' dbOut '.txt instead']);
         fid = fopen([dbOut '.txt'],'w');
         for n = 1:numel(trig)
            fprintf(fid,'%-7s %-7s %-s %5s\n',sta,chan,datestr(trig(n),'yyyy/mm/dd HH:MM:SS.FFF'),phaseList{n});
         end
         fclose(fid);
      end
   end

   % WRITE DETECTION TABLE
   if ISDETECTION
      trig = c.trig;
      sta = c.stations;
      chan = c.channels;
      % TODO: check for database table
      disp(['Writing database table ' dbOut '.detection']);
      db = dbopen(dbOut,'r+');
      db = dblookup(db,'','detection','','');
      nMax = length(trig);
      disp('     ');
      for n = 1:nMax
         arid = dbnextid( db,'arid');
         db.record = dbaddv(db,'sta',sta{n},'chan',chan{n},'time',datenum2epoch(trig(n)),'state','D','filter',filterString);
         fprintf('\b\b\b\b\b\b%5.0f%%',n/nMax*100);
      end
      fprintf('\n');
      dbclose(db);

   end
end

author: geoscience-community-codes

repo: GISMO


text:
% created by: Zoya Bylinskii, Sept 2014
% based on code by Tilke Judd

% This function applies the chosen center bias
% and blur parameters to the saliency maps
% in salDir to produce the saliency maps in resDir.

function adjustSalMaps(salDir,resDir,centWeight,blurSigma,targetHist,showOutput)
% salDir contains all the saliency map files
% resDir is the output directory for the adjusted saliency maps
% centWeight: the chosen center weight value [0,1]
% blurSigma: the chosen blur value

assert(centWeight>=0 && centWeight<= 1, 'centWeight must be between 0 and 1');

if nargin < 6,
    showOutput = 0;
end

load('center.mat');

if exist(resDir,'dir')
    fprintf('Directory %s already exists. Will overwrite.\n',resDir);
else
    mkdir(resDir)
end

allfiles = dir(fullfile(salDir,'*.jpg'));

if showOutput, figure; end

for i = 1:length(allfiles)

    fprintf('On %s\n',allfiles(i).name)
    mapOrig = im2double(imread(fullfile(salDir,allfiles(i).name)));
    cent = imresize(center, size(mapOrig));
    map = processMap(mapOrig,cent,targetHist,blurSigma,centWeight,showOutput);
    imwrite(map,fullfile(resDir,allfiles(i).name));

end

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
classdef TSLDirectionKey < directionColorKey
  % converts directions to rgb values

  methods

    function dM = TSLDirectionKey(varargin)
      dM@directionColorKey(varargin{:});
      dM.sym = dM.sym.Laue;
      dM.sR = dM.sym.fundamentalSector;

      if ismember(dM.sym.id,[2,5,8,11,18,21,24,27,35,42])
        warning('Not a topological correct colormap! Green to blue colorjumps possible');
      end
    end

    function rgb = direction2color(dM,h,varargin)
      % in TSL all fundamental sectors are colorized with white in the center

      % project to fundamental region
      h = h.project2FundamentalRegion(dM.sym);

      % this should become white if not stated differently
      center = dM.sR.center;

      % compute angle of the points "sh" relative to the center point "center"
      % this should be between 0 and 1
      v = dM.sR.vertices;
      if ~isempty(v)
        v = v(1);
      else
        v  = dM.sym.aAxisRec;
      end
      [radius,rho] = polarCoordinates(dM.sR,h,center,v);

      % white center
      radius = 0.5+radius./2;

      % compute the color vector on the sphere
      v = vector3d('rho',rho,'theta',radius.*pi);

      % post processing of the color vector
      % by default we have white at the z, black at the -z, red
      % at x and green and blue at 120 and 240 degree accordingly
      % post rotate the color
      v = dM.colorPostRotation * v;

      % compute rgb values
      rgb = ar2rgb(mod(v.rho./ 2 ./ pi,1),v.theta./pi,get_option(varargin,'grayValue',1),'noHueCorrection');

    end
  end

end

author: mtex-toolbox

repo: mtex

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


text:
function [r,u,t,p] = macroproperties1d(n,j_x,E,nx,nv,theta)
%% Recover Macroscopic Properties
% compute back, fugacity, macroscopic velocities, temperature and pressure.
    % Computing first velocite,s from the momentum:
    u = j_x./n;

% to compute fugacity, temperature and pressure, we need to rely on the
% distribution fucntion that we where are using: MB, FD, BE.

switch theta
    case{-1} % BE
    % If BE: we apply bisection method to the approx BE distribution Eq.
        r_a = 0.001; r_b = 0.99; tol = 1e-7;
        for i = 1:nx
        psi = @(r_x) 2*E(i)- BE(r_x,1.5)*(n(i)/BE(r_x,0.5))^3/(2*pi) ...
        - n(i)*(u(i)^2);
        r_p = bisection(psi,r_a,r_b,tol);
        r(i) = r_p;
        t(i) = n(i)^2/(pi*(BE(r_p,0.5))^2);
        p(i) = E(i) - 1/2*n(i)*(u(i)^2);
        end

    case{1} % FD
    % if FD: we apply bisection method to the approx FD distribution Eq.
        r_a = 0.001; r_b = 0.99; tol = 1e-7;
        for i = 1:nx
        psi = @(r_x) 2*E(i)- FD(r_x,1.5)*(n(i)/FD(r_x,0.5))^3/(2*pi) ...
        - n(i)*(u(i)^2);
        r_p = bisection(psi,r_a,r_b,tol);
        r(i) = r_p;
        t(i) = n(i)^2/(pi*(FD(r_p,0.5))^2);
        p(i) = E(i) - 1/2*n(i)*(u(i)^2);
        end

    case{0} % MB
    % IF MB: the task is much simple.
        t = 4*E./n - 2*u.^2;
        r = n./sqrt(pi.*t);
        p = E - 1/2.*n.*u.^2;
    otherwise
        error('theta can only be: -1, 0, +1 ');
end

% Using Discrete Ordinate Method:
%     r = repmat(r,nv,1);     u = repmat(u,nv,1);
%     t = repmat(t,nv,1);     p = repmat(p,nv,1);
%     n = repmat(n,nv,1);


text:
function [lat2,lon2,h2]=direct(lat1,lon1,h1,az,va,d,a,e2)
% DIRECT  Computes direct (forward) geodetic problem.
%   Determines coordinates of 2nd station given ellipsoidal
%   coordinates of 1st station and azimuth, vertical angle
%   and distance from 1st to 2nd station. If az,va are local
%   astronomic, lat,lon must also be astronomic. If az,va
%   are local geodetic, lat,lon must be local geodetic.
%   Non-vectorized. See also INVERSE.
% Version: 2011-02-19
% Useage:  [lat2,lon2,h2]=direct(lat1,lon1,h1,az,va,d,a,e2)
%          [lat2,lon2,h2]=direct(lat1,lon1,h1,az,va,d)
% Input:   lat1 - ellipsoidal latitude of 1st station (rads)
%          lon1 - ellipsoidal longitude of 1st station (rads)
%          h1   - ellipsoidal ht. of 1st station (m)
%          az   - azimuth from station 1 to 2 (rads)
%          va   - vertical angle from 1 to 2 (rads)
%          d    - distance from 1 to 2 (m)
%          a    - ref. ellipsoid major semi-axis (m); default GRS80
%          e2   - ref. ellipsoid eccentricity squared; default GRS80
% Output:  lat2 - ellipsoidal latitude of 2nd station (rads)
%          lon2 - ellipsoidal longitude of 2nd station (rads)
%          h2   - ellipsoidal ht. of 2nd station (m)

% Copyright (c) 2011, Michael R. Craymer
% All rights reserved.
% Email: mike@craymer.com

if nargin ~= 6 & nargin ~= 8
  warning('Incorrect number of input arguments');
  return
end
if nargin == 6
  [a,b,e2]=refell('grs80');
end

[X1,Y1,Z1]=ell2xyz(lat1,lon1,h1,a,e2);
[dx,dy,dz]=sph2xyz(az,va,d);
[dX,dY,dZ]=lg2ct(dx,dy,dz,lat1,lon1);
X2=X1+dX;
Y2=Y1+dY;
Z2=Z1+dZ;
[lat2,lon2,h2]=xyz2ell(X2,Y2,Z2,a,e2);

author: Sable

repo: mcbench-benchmarks


text:
function [Y, G, A] = mlpOut(model, X);

% MLPOUT Output of an MLP model.
% FORMAT
% DESC gives the output of a multi-layer perceptron model, for
% single hidden layer models the function is a wrapper for mlpfwd.
% ARG model : the model for which the output is required.
% ARG X : the input data for which the output is required.
% RETURN Y : the output.
%
% FORMAT
% DESC gives the output of a multi-layer perceptron model.
% ARG model : the model for which the output is required.
% ARG X : the input data for which the output is required.
% RETURN Y : the output.
% RETURN G : the hidden layer activations.
%
% SEEALSO : mlpfwd, mlp, modelOut
%
% COPYRIGHT : Neil D. Lawrence, 2006, 2007

% MLTOOLS

if length(model.hiddenDim)==1
  if nargout > 1
    if nargout > 2
      [Y, G, A] = mlpfwd(model, X);
    else
      [Y, G] = mlpfwd(model, X);
    end
  else
    Y = mlpfwd(model, X);
  end
else
  ndata = size(x, 1);
  G{1} = tanh(X*model.w{1} + repmat(model.b{1}, numData, 1));
  A{1} = G{1}*model.w{2} + repmat(model.b{2}, numData, 1);
  for i = 2:length(model.numHidden)
    G{i} = tanh(A{i-1});
    A{i} = G{i}*model.w{i+1} + repmat(model.b{i+1}, numData, 1);
  end
  switch model.outfn

   case 'linear'
    y = A{end};

   otherwise
    error('Output function not implemented in multiple hidden layer model.')

  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function f=loglikGaPfromVWH(varargin)
% f=loglikGaPfromVWH(varargin)
% complete log likellihood of the GaP model
% Vxt =     varargin{1};        %data (XxT)
% Wxk =     varargin{2};        %basis functions (XxK)
% Hkt =     varargin{3};        %intensities (KxT)
% alpha =   varargin{4};        %parameters of the Gamma blinking (Kx1)
% beta =    varargin{5};        %parameters of the Gamma blinking (Kx1)
% peval =   varargin{6};        %parameters
% Gamma distribution with mean alpha/beta!
% X: number of pixels
% T: number of images (time - slices)
% K: number of components

Vxt =     varargin{1};        %data (XxT)
Wxk =     varargin{2};        %basis functions (XxK)
Hkt =     varargin{3};        %intensities (KxT)
alpha =   varargin{4};        %parameters of the Gamma blinking (Kx1)
beta =    varargin{5};        %parameters of the Gamma blinking (Kx1)
peval =   varargin{6};        %parameters

k=size(Wxk,2);
if length(alpha) ==1; alpha = repmat(alpha,k, 1); end
if length(beta) ==1; beta = repmat(beta,k, 1); end
%[Wxkbg,Hktbg]=addbg(Wxk, Hkt, peval.bg);
%P=Wxkbg*Hktbg; %current approximation (XxT)
P=Wxk*Hkt; %current approximation (XxT)

%Poisson contribution
t1=Vxt.*log(P) - P - factorialapprox(Vxt); %(XxT)
%Gamma contribution
t2a=bsxfun(@times,(alpha-1),log(Hkt)) - bsxfun(@times,beta,Hkt); % (KxT)
t2b=-alpha.*log(beta)-log(gamma(alpha)); %(1xK)

f=sum(t1(:))+sum(t2a(:))+sum(t2b(:));
% Conjugate gradient is mimimizing! for this it needs to be f=-f;

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function pass = test_nonlinSys1Breaks_C2(pref)
% Test 2x2 system (sin/cos). This is pecewiseificaion of the test
%       test_nonlinearSystem1
%
% Asgeir Birkisson, April 2014.

if ( nargin == 0 )
    pref = cheboppref;
end

tol = 1e-10;

d = [-pi 0 pi];
x = chebfun('x',d);
f = [ 0*x ; 0*x ];

%% Piecewise (chebcolloc2):
pref.discretization = @chebcolloc2;

A = chebop(@(x,u,v) [u - diff(v,2) + u.^2; diff(u) + sin(v)],d);
A.lbc = @(u,v) u-1;
A.rbc = @(u,v) [v-1/2; diff(v)];

u = mldivide(A, f, pref);
u1 = u{1}; u2 = u{2};

% Want to check BCs as well.
bcFunLeft = A.lbc(u1,u2);
bcFunRight = chebfun(A.rbc(u1,u2));

% And check that we're continuous over breakpoint
u1jump = jump(u1, 0);
u2jump = jump(u2, 0);

pass(1) = norm( chebfun(A(x, u1, u2))) < tol;
pass(2) = norm(bcFunLeft(d(1))) < tol && norm(bcFunRight(d(end))) < tol;
pass(3) = norm(u1jump) < tol && norm(u2jump) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% ------------------------------------------------------
% SwarmOps - Heuristic optimization for Matlab
% Copyright (C) 2003-2010 Magnus Erik Hvass Pedersen.
% Please see the file license.txt for license details.
% SwarmOps on the internet: http://www.Hvass-Labs.org/
% ------------------------------------------------------

% Behavioural parameters for Many Optimizing Liaisons (MOL)
% tuned by Pedersen (1). The parameter-array consists of
% the following parameters:
% - Swarm-size (denoted s)
% - Inertia weight (denoted omega)
% - Swarm's best weight (denoted phiG)
%
% Select the parameters that most closely match the
% characteristics of your optimization problem.
% For example, if you want to optimize a problem where
% the search-space has 7 dimensions and you can perform
% 5000 evaluations, then you could first try using the
% parameters MOL_5DIM_10000EVALS. If that does not yield
% satisfactory results then you could try
% MOL_10DIM_2000EVALS or perhaps MOL_10DIM_20000EVALS.
% If that does not work then you will either need to tune
% the parameters for the problem at hand, or you should
% try using another optimizer.
%
% Literature references:
% (1) M.E.H. Pedersen.
%     Good parameters for Particle Swarm Optimization.
%     Technical Report HL1001, Hvass Laboratories, 2010.

% Parameters for non-parallel version:

MOL_2DIM_400EVALS_A       = [ 23, -0.3328, 2.8446];
MOL_2DIM_400EVALS_B       = [ 50,  0.2840, 1.9466];
MOL_2DIM_4000EVALS_A      = [183, -0.2797, 3.0539];
MOL_2DIM_4000EVALS_B      = [139,  0.6372, 1.0949];
MOL_5DIM_1000EVALS        = [ 50, -0.3085, 2.0273];
MOL_5DIM_10000EVALS       = [ 96, -0.3675, 4.1710];
MOL_10DIM_2000EVALS       = [ 60, -0.2700, 2.9708];
MOL_10DIM_20000EVALS      = [116, -0.3518, 3.8304];
MOL_20DIM_40000EVALS      = [228, -0.3747, 4.2373];
MOL_20DIM_400000EVALS_A   = [125, -0.2575, 4.6713];
MOL_20DIM_400000EVALS_B   = [125, -0.2575, 4.6713];
MOL_30DIM_6000EVALS       = [ 67, -0.4882, 2.7923];
MOL_30DIM_60000EVALS      = [198, -0.2723, 3.8283];
MOL_30DIM_600000EVALS     = [134, -0.4300, 3.0469];
MOL_50DIM_100000EVALS     = [290, -0.3067, 3.6223];
MOL_100DIM_200000EVALS    = [219, -0.1685, 3.9162];
MOL_DEFAULT               = MOL_30DIM_60000EVALS;

% Parameters for parallel version (above may also work):

MOL_PAR_5DIM_10000EVALS   = [ 32, -0.3319, 5.9100];
MOL_PAR_30DIM_60000EVALS  = [164, -0.4241, 2.7729];
MOL_PAR_DEFAULT           = MOL_PAR_5DIM_10000EVALS;

% ------------------------------------------------------


text:
function [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ....
    = bk(fid, data, varargins, opts)
% This program starts the BK algorithm
%
% function [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ....
%           = bk_start(fid, data, varargins, opts)
%
% cum_ret: a number representing the final cumulative wealth.
% cumprod_ret: cumulative return until each trading period
% daily_ret: individual returns for each trading period
% daily_portfolio: individual portfolio for each trading period
%
% data: market sequence vectors
% fid: handle for write log file
% varargins: variable parameters
% opts: option parameter for behvaioral control
%
% Example: [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ...
%          = bk_start(fid, data, {0}, opts);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract the parameters
K  = varargins{1};
L  = varargins{2};
c  = varargins{3};
tc = varargins{4};      % transaction cost fee rate

% Call BK's run algorithm
[cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ...
    = bk_run(fid, data, K, L, c, tc, opts);

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function [a3,a2,a1,a0] = createTraj3(theta0,thetaf,thetad0,thetadf,tstart,tfinal)
	% inputs : initial position, velocity + final position, velocity + initial and final times
	% output : a vector specifying the polynomial and can be used with poly functions such as : polyder, polyval, etc.
	% create a 3rd order trajectory
	% example:
	% createTraj3(10,30,0,0,0,1)
	%
	%
	% By: Reza Ahmadzadeh - Matlab/Octave - 2013
	T = tfinal - tstart;
	a0 = theta0;
	a1 = thetad0;
	a2 = (-3 * (theta0 - thetaf) - (2 * thetad0+thetadf )*T)/ T ^ 2;
	a3 = (2 * (theta0 - thetaf) + (thetad0+thetadf )*T)/ T ^ 3;

end


text:
function ID = spm_data_id(varargin)
% generates a specific real number in a deterministic way
% from any data structure
% FORMAT ID = spm_data_id(X);
% X  - numeric, character, cell or stucture array[s]
% ID - specific ID
%__________________________________________________________________________
% Copyright (C) 2009 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak (based on Karl's spm_vec)
% $Id: spm_data_id.m 6712 2016-02-04 15:12:25Z peter $

X     = varargin;

if length(X) == 1
    X = X{1};
end

ID = 0;

if ischar(X) % For now strings are not taken into account
    ID = 0;
elseif isnumeric(X)
    Y = double(X(:));
    ID = sum(abs(Y(~isnan(Y) & ~isinf(Y))));
elseif isstruct(X) || isobject(X)
    X = struct(X);
    f = fieldnames(X);
    X = X(:);
    ID = 0;
    for i = 1:length(f)
        ID = ID + spm_data_id({X.(f{i})});
    end
elseif iscell(X)
    X     = X(:);
    ID = 0;
    for i = 1:length(X)
        ID = ID + spm_data_id(X{i});
    end
end

if ID > 0
    ID = 10^-(floor(log10(ID))-2)*ID;
end

author: spm


text:
% DIRECT JACOBIAN DEMO
function jacobian_planar_3gdl

%a)
q1=0
q2 = pi/2
q3 = pi/2
J = [-sin(q1)-sin(q1+q2)-sin(q1+q2+q3)  -sin(q1+q2)-sin(q1+q2+q3) -sin(q1+q2+q3);
     cos(q1)+cos(q1+q2)+cos(q1+q2+q3)    cos(q1+q2)+cos(q1+q2+q3)  cos(q1+q2+q3);
      0 0 0;
      0 0 0;
      0 0 0;
      1 1 1];
 f = [ 1 1 1 1 1 1]';

 tau = J'*f

 %b)
q1=pi/4
q2 = pi/4
q3 = pi/4
J = [-sin(q1)-sin(q1+q2)-sin(q1+q2+q3)  -sin(q1+q2)-sin(q1+q2+q3) -sin(q1+q2+q3);
     cos(q1)+cos(q1+q2)+cos(q1+q2+q3)    cos(q1+q2)+cos(q1+q2+q3)  cos(q1+q2+q3);
      0 0 0;
      0 0 0;
      0 0 0;
      1 1 1];

 tau = [1 1 1]';

 f = J*inv(J'*J)*tau

%c)
q1=pi/2
q2 = pi/2
q3 = pi/2
syms q1 q2 q3 real
J = [-sin(q1)-sin(q1+q2)-sin(q1+q2+q3)  -sin(q1+q2)-sin(q1+q2+q3) -sin(q1+q2+q3);
     cos(q1)+cos(q1+q2)+cos(q1+q2+q3)    cos(q1+q2)+cos(q1+q2+q3)  cos(q1+q2+q3)];

 %Jp = pinv(J);
 Jp = J'*inv(J*J');
 P = (eye(3) - Jp*J)
 P = simplify(P)
 qdn = P*[0 0 1]'
 qdn = simplify(qdn)

repo: ARTE


text:
function p = eliminateBinary(p,binaries)

vars = p.lmi_variables;
[mt,vt] = yalmip('monomtable');

mt = mt(vars,:);
mt(:,binaries) = min(mt(:,binaries),1);

used_variables = find(any(mt,1));
x = recover(used_variables)';
new_monoms = [];
mt = mt(:,used_variables);

y = recovermonoms(mt,x);
y = p.basis*[1;y];
if isa(y,'sdpvar')
    % copy data to p
    p.basis = y.basis;
    p.lmi_variables = y.lmi_variables;
else
    p = y;
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


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 q = scale(Q, r)
%Quaternion.scale Interpolate rotations expressed by quaternion objects
%
% QI = Q.scale(R) is a unit-quaternion that interpolates between identity for R=0
% to Q for R=1.  This is a spherical linear interpolation (slerp) that can
% be interpretted as interpolation along a great circle arc on a sphere.
%
% If R is a vector QI is a cell array of quaternions, each element
% corresponding to sequential elements of R.
%
% See also ctraj, Quaternion.interp.

    q2 = double(Q);

    if any(r<0) || (r>1)
        error('r out of range');
    end
    q1 = [1 0 0 0];         % identity quaternion
    theta = acos(q1*q2');

    if length(r) == 1
        if theta == 0
            q = Q;
        else
            q = unit(Quaternion( (sin((1-r)*theta) * q1 + sin(r*theta) * q2) / sin(theta) ));
        end
    else
        count = 1;
        for R=r(:)'
            if theta == 0
                qq = Q;
            else
                qq = Quaternion( (sin((1-r)*theta) * q1 + sin(r*theta) * q2) / sin(theta) ).unit;
            end
            q(count) = qq;
            count = count + 1;
        end
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function out = zoom(varargin)
    %ZOOM   Zoom in and out on a 2-D plot.
    %   ZOOM with no arguments toggles the zoom state.
    %   ZOOM(FACTOR) zooms the current axis by FACTOR.
    %       Note that this does not affect the zoom state.
    %   ZOOM ON turns zoom on for the current figure.
    %   ZOOM OFF turns zoom off in the current figure.
    %   ZOOM OUT returns the plot to its initial (full) zoom.
    %   ZOOM XON or ZOOM YON turns zoom on for the x or y axis only.
    %   ZOOM RESET clears the zoom out point.
    %
    %   When zoom is on, click the left mouse button to zoom in on the
    %   point under the mouse.  Click the right mouse button to zoom out
    %   (shift-click on the Macintosh).  Each time you click, the axes
    %   limits will be changed by a factor of 2 (in or out).  You can also
    %   click and drag to zoom into an area.  Double clicking zooms out to
    %   the point at which zoom was first turned on for this figure.  Note
    %   that turning zoom on, then off does not reset the zoom point.
    %   This may be done explicitly with ZOOM RESET.
    %
    %   ZOOM(FIG,OPTION) applies the zoom command to the figure specified
    %   by FIG. OPTION can be any of the above arguments.

    %   ZOOM FILL scales a plot such that it is as big as possible
    %   within the axis position rectangle for any azimuth and elevation.

    %   Clay M. Thompson 1-25-93
    %   Revised 11 Jan 94 by Steven L. Eddins
    %   Copyright (c) 1984-97 by The MathWorks, Inc.
    %   $Revision: 1399 $  $Date: 2006-08-11 11:19:27 +0200 (Fr, 11 Aug 2006) $

    %   Note: zoom uses the userdata of the zlabel of the axis and
    %   the figure buttondown and buttonmotion functions
    %
    %   ZOOM XON zooms x-axis only
    %   ZOOM YON zooms y-axis only

    switch nargin

        %%%%%%%%%%%%%%%%%%%%%%%%%%
        %%% No Input Arguments %%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%
        case 0
            fig=get(groot,'currentfigure');
            if isempty(fig), return, end
            zoomCommand='toggle';

            %%%%%%%%%%%%%%%%%%%%%%%%%%
            %%% One Input Argument %%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%
        case 1

            % If the argument is a string, the argument is a zoom command
            % (i.e. (on, off, down, xdown, etc.).  Otherwise, the argument is
            % assumed to be a figure handle, in which case all we do is
            % toggle the zoom status.

            if ischar(varargin{1})
                fig=get(groot,'currentfigure');
                if isempty(fig), return, end

                zoomCommand=varargin{1};
            else
                scale_factor=varargin{1};
                zoomCommand='scale';
                fig = gcf;
            end % if

            %%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%% Two Input Arguments %%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%
        case 2
            fig=varargin{1};
            zoomCommand=varargin{2};

        otherwise
            narginchk(0, 2);

    end % switch nargin

    %
    % handle 'off' commands first
    %
    if strcmp(zoomCommand,'off')
        %
        % turn off zoom, and take a hike
        %
        fcns = getappdata(fig,'ZOOMFigureFcns');
        if ~isempty(fcns)
            set(fig,'windowbuttondownfcn',fcns.wbdf,'windowbuttonupfcn',fcns.wbuf,...
                'windowbuttonmotionfcn',fcns.wbmf,'buttondownfcn',fcns.bdf);
            rmappdata(fig,'ZOOMFigureFcns');
        end
        return
    end % if

    ax=get(fig,'currentaxes');

    rbbox_mode = 0;
    zoomx = 1; zoomy = 1; % Assume no constraints

    zoomCommand=lower(zoomCommand);

    if ~isempty(isempty(ax)) &&  any(get(ax,'view')~=[0 90])  &&  ...
            ~(strcmp(zoomCommand,'scale') | strcmp(zoomCommand,'fill'))
        return % Do nothing
    end

    if strcmp(zoomCommand,'toggle')
        fcns = getappdata(fig,'ZOOMFigureFcns');
        if isempty(fcns)
            zoom(fig,'on');
        else
            zoom(fig,'off');
        end
        return
    end % if

    % Catch constrained zoom
    if strcmp(zoomCommand,'xdown')
        zoomy = 0; zoomCommand = 'down'; % Constrain y
    elseif strcmp(zoomCommand,'ydown')
        zoomx = 0; zoomCommand = 'down'; % Constrain x
    end

    if strcmp(zoomCommand,'down')
        % Activate axis that is clicked in
        allAxes = findobj(get(fig,'Children'),'flat','type','axes');
        ZOOM_found = 0;
        for i=1:length(allAxes)
            ax=allAxes(i);
            ZOOM_Pt1 = get(ax,'CurrentPoint');
            xlim = get(ax,'xlim');
            ylim = get(ax,'ylim');
            if (xlim(1) <= ZOOM_Pt1(1,1) & ZOOM_Pt1(1,1) <= xlim(2) & ...
                    ylim(1) <= ZOOM_Pt1(1,2) & ZOOM_Pt1(1,2) <= ylim(2))
                ZOOM_found = 1;
                set(fig,'currentaxes',ax);
                break
            end % if
        end % for

        if ZOOM_found==0, return, end

        % Check for selection type
        selection_type = get(fig,'SelectionType');
        if (strcmp(selection_type, 'normal'))
            % Zoom in
            m = 1;
            scale_factor = 2; % the default zooming factor
        elseif (strcmp(selection_type, 'open'))
            % Zoom all the way out
            zoom(fig,'out');
            return;
        else
            % Zoom partially out
            m = -1;
            scale_factor = 2;
        end

        ZOOM_Pt1 = get_currentpoint(ax);
        ZOOM_Pt2 = ZOOM_Pt1;
        center = ZOOM_Pt1;

        if (m == 1)
            % Zoom in
            units = get(fig,'units'); set(fig,'units','pixels')
            rbbox([get(fig,'currentpoint') 0 0],get(fig,'currentpoint'));
            ZOOM_Pt2 = get_currentpoint(ax);
            set(fig,'units',units)

            % Note the currentpoint is set by having a non-trivial up function.
            if min(abs(ZOOM_Pt1-ZOOM_Pt2)) >= ...
                    min(.01*[diff(get_xlim(ax)) diff(get_ylim(ax))])
                % determine axis from rbbox
                a = [ZOOM_Pt1;ZOOM_Pt2]; a = [min(a);max(a)];

                % Undo the effect of get_currentpoint for log axes
                if strcmp(get(ax,'XScale'),'log')
                    a(1:2) = 10.^a(1:2);
                end
                if strcmp(get(ax,'YScale'),'log')
                    a(3:4) = 10.^a(3:4);
                end
                rbbox_mode = 1;
            end
        end
        limits = zoom(fig,'getlimits');

    elseif strcmp(zoomCommand,'scale')
        if all(get(ax,'view')==[0 90]), % 2D zooming with scale_factor

            % Activate axis that is clicked in
            ZOOM_found = 0;
            ax = gca;
            xlim = get(ax,'xlim');
            ylim = get(ax,'ylim');
            ZOOM_Pt1 = [sum(xlim)/2 sum(ylim)/2];
            ZOOM_Pt2 = ZOOM_Pt1;
            center = ZOOM_Pt1;

            if (xlim(1) <= ZOOM_Pt1(1,1) & ZOOM_Pt1(1,1) <= xlim(2) & ...
                    ylim(1) <= ZOOM_Pt1(1,2) & ZOOM_Pt1(1,2) <= ylim(2))
                ZOOM_found = 1;
            end % if

            if ZOOM_found==0, return, end

            if (scale_factor >= 1)
                m = 1;
            else
                m = -1;
            end

        else % 3D
            old_CameraViewAngle = get(ax,'CameraViewAngle')*pi/360;
            ncva = atan(tan(old_CameraViewAngle)*(1/scale_factor))*360/pi;
            set(ax,'CameraViewAngle',ncva);
            return;
        end

        limits = zoom(fig,'getlimits');

    elseif strcmp(zoomCommand,'on')
        fcns = getappdata(fig,'ZOOMFigureFcns');
        if isempty(fcns)
            fcns.wbdf = get(fig,'windowbuttondownfcn');
            fcns.wbuf = get(fig,'windowbuttonupfcn');
            fcns.wbmf = get(fig,'windowbuttonmotionfcn');
            fcns.bdf  = get(fig,'buttondownfcn');
            setappdata(fig,'ZOOMFigureFcns',fcns);
        end
        set(fig,'windowbuttondownfcn','zoom down', ...
            'windowbuttonupfcn','ones;', ...
            'windowbuttonmotionfcn','','buttondownfcn','', ...
            'interruptible','on');
        set(ax,'interruptible','on')
        return

    elseif strcmp(zoomCommand, 'reset')
        hZlabel = get(ax, 'Zlabel');
        ZlabelUserData = get(hZlabel, 'UserData');

        if IsZoomData(ZlabelUserData)
            set(hZlabel, 'UserData', []);
        end
        return

    elseif strcmp(zoomCommand,'xon')
        zoom(fig,'on') % Set up userprop
        set(fig,'windowbuttondownfcn','zoom xdown', ...
            'windowbuttonupfcn','ones;', ...
            'windowbuttonmotionfcn','','buttondownfcn','',...
            'interruptible','on');
        set(ax,'interruptible','on')
        return

    elseif strcmp(zoomCommand,'yon')
        zoom(fig,'on') % Set up userprop
        set(fig,'windowbuttondownfcn','zoom ydown', ...
            'windowbuttonupfcn','ones;', ...
            'windowbuttonmotionfcn','','buttondownfcn','',...
            'interruptible','on');
        set(ax,'interruptible','on')
        return

    elseif strcmp(zoomCommand,'out')
        limits = zoom(fig,'getlimits');
        center = [sum(get_xlim(ax))/2 sum(get_ylim(ax))/2];
        m = -inf; % Zoom totally out

    elseif strcmp(zoomCommand,'getlimits'), % Get axis limits
        limits = get(get(ax,'ZLabel'),'UserData');
        % Do simple checking of userdata
        if size(limits,2)==4  &&  size(limits,1)<=2
            if all(limits(1,[1 3])<limits(1,[2 4]))
                getlimits = 0; out = limits(1,:); return   % Quick return
            else
                getlimits = -1; % Don't munge data
            end
        else
            if isempty(limits), getlimits = 1; else getlimits = -1; end
        end

        % If I've made it to here, we need to compute appropriate axis
        % limits.

        if isempty(get(get(ax,'ZLabel'),'userdata'))
            % Use quick method if possible
            xlim = get_xlim(ax); xmin = xlim(1); xmax = xlim(2);
            ylim = get_ylim(ax); ymin = ylim(1); ymax = ylim(2);

        elseif strcmp(get(ax,'xLimMode'),'auto')  &&  ...
                strcmp(get(ax,'yLimMode'),'auto')
            % Use automatic limits if possible
            xlim = get_xlim(ax); xmin = xlim(1); xmax = xlim(2);
            ylim = get_ylim(ax); ymin = ylim(1); ymax = ylim(2);

        else
            % Use slow method only if someone else is using the userdata
            h = get(ax,'Children');
            xmin = inf; xmax = -inf; ymin = inf; ymax = -inf;
            for i=1:length(h)
                t = get(h(i),'Type');
                if ~strcmp(t,'text')
                    if strcmp(t,'image'), % Determine axis limits for image
                        x = get(h(i),'Xdata'); y = get(h(i),'Ydata');
                        x = [min(min(x)) max(max(x))];
                        y = [min(min(y)) max(max(y))];
                        [ma,na] = size(get(h(i),'Cdata'));
                        if na>1, dx = diff(x)/(na-1); else dx = 1; end
                        if ma>1, dy = diff(y)/(ma-1); else dy = 1; end
                        x = x + [-dx dx]/2; y = y + [-dy dy]/2;
                    end
                    xmin = min(xmin,min(min(x)));
                    xmax = max(xmax,max(max(x)));
                    ymin = min(ymin,min(min(y)));
                    ymax = max(ymax,max(max(y)));
                end
            end

            % Use automatic limits if in use (override previous calculation)
            if strcmp(get(ax,'xLimMode'),'auto')
                xlim = get_xlim(ax); xmin = xlim(1); xmax = xlim(2);
            end
            if strcmp(get(ax,'yLimMode'),'auto')
                ylim = get_ylim(ax); ymin = ylim(1); ymax = ylim(2);
            end
        end

        limits = [xmin xmax ymin ymax];
        if getlimits~=-1, % Don't munge existing userdata.
            % Store limits in ZLabel userdata
            set(get(ax,'ZLabel'),'UserData',limits);
        end

        out = limits;
        return

    elseif strcmp(zoomCommand,'getconnect'), % Get connected axes
        limits = get(get(ax,'ZLabel'),'UserData');
        if all(size(limits)==[2 4]), % Do simple checking
            out = limits(2,[1 2]);
        else
            out = [ax ax];
        end
        return

    elseif strcmp(zoomCommand,'fill')
        old_view = get(ax,'view');
        view(45,45);
        set(ax,'CameraViewAngleMode','auto');
        set(ax,'CameraViewAngle',get(ax,'CameraViewAngle'));
        view(old_view);
        return

    else
        error(['Unknown option: ',zoomCommand,'.']);
    end

    %
    % Actual zoom operation
    %

    if ~rbbox_mode
        xmin = limits(1); xmax = limits(2);
        ymin = limits(3); ymax = limits(4);

        if m==(-inf)
            dx = xmax-xmin;
            dy = ymax-ymin;
        else
            dx = diff(get_xlim(ax))*(scale_factor.^(-m-1)); dx = min(dx,xmax-xmin);
            dy = diff(get_ylim(ax))*(scale_factor.^(-m-1)); dy = min(dy,ymax-ymin);
        end

        % Limit zoom.
        center = max(center,[xmin ymin] + [dx dy]);
        center = min(center,[xmax ymax] - [dx dy]);
        a = [max(xmin,center(1)-dx) min(xmax,center(1)+dx) ...
            max(ymin,center(2)-dy) min(ymax,center(2)+dy)];

        % Check for log axes and return to linear values.
        if strcmp(get(ax,'XScale'),'log')
            a(1:2) = 10.^a(1:2);
        end
        if strcmp(get(ax,'YScale'),'log')
            a(3:4) = 10.^a(3:4);
        end

    end

    % Check for v4-type equal
    fillequal = strcmp(get(ax,'plotboxaspectratiomode'),'manual') & ...
        strcmp(get(ax,'dataaspectratiomode'),'manual');
    pbar = get(ax,'plotboxaspectratio');
    % Update circular list of connected axes
    list = zoom(fig,'getconnect'); % Circular list of connected axes.
    if zoomx
        if a(1)==a(2), return, end % Short circuit if zoom is moot.
        if fillequal & (pbar(1) < pbar(2))
            set(ax,'xlimmode','auto')
        else
            set(ax,'xlim',a(1:2))
        end
        h = list(1);
        while h ~= ax
            if fillequal & zoomx & zoomy & (pbar(1) < pbar(2))
                set(h,'xlimmode','auto')
            else
                set(h,'xlim',a(1:2))
            end
            % Get next axes in the list
            next = get(get(h,'ZLabel'),'UserData');
            if all(size(next)==[2 4]), h = next(2,1); else h = ax; end
        end
    end
    if zoomy
        if a(3)==a(4), return, end % Short circuit if zoom is moot.
        if fillequal & (pbar(1) >= pbar(2))
            set(ax,'ylimmode','auto')
        else
            set(ax,'ylim',a(3:4))
        end
        h = list(2);
        while h ~= ax
            if fillequal & zoomx & zoomy & (pbar(1) >= pbar(2))
                set(h,'ylimmode','auto')
            else
                set(h,'ylim',a(3:4))
            end
            % Get next axes in the list
            next = get(get(h,'ZLabel'),'UserData');
            if all(size(next)==[2 4]), h = next(2,2); else h = ax; end
        end
    end

function bZoomData = IsZoomData(data)
    % Return 1 if the data represents zoom data
    % Return 0 if someone else is using user data

    if size(data,2)==4  &&  size(data,1)<=2
        if all(data(1,[1 3])<data(1,[2 4]))
            bZoomData = 1;
        else
            bZoomData = 0;
        end
    else
        bZoomData = 0;
    end

function p = get_currentpoint(ax)
    %GET_CURRENTPOINT Return equivalent linear scale current point
    p = get(ax,'currentpoint'); p = p(1,1:2);
    if strcmp(get(ax,'XScale'),'log')
        p(1) = log10(p(1));
    end
    if strcmp(get(ax,'YScale'),'log')
        p(2) = log10(p(2));
    end

function xlim = get_xlim(ax)
    %GET_XLIM Return equivalent linear scale xlim
    xlim = get(ax,'xlim');
    if strcmp(get(ax,'XScale'),'log')
        xlim = log10(xlim);
    end

function ylim = get_ylim(ax)
    %GET_YLIM Return equivalent linear scale ylim
    ylim = get(ax,'ylim');
    if strcmp(get(ax,'YScale'),'log')
        ylim = log10(ylim);
    end


text:
function [HR, HRS, HlR, HlRS, HiR, HiRS, ChiR, HshR, HshRS] = gaussian_method_v7_1_0(R, pars)
%GAUSSIAN_METHOD_V7_1_0 Estimate entropy values using the Gaussian Method.
%
% ---------
% ARGUMENTS
% ---------
% R    - response matrix
% pars - parameters structure

%   Copyright (C) 2009 Cesare Magri
%   Version 7.1.0

% -------
% LICENSE
% -------
% This software is distributed free under the condition that:
%
% 1. it shall not be incorporated in software that is subsequently sold;
%
% 2. the authorship of the software shall be acknowledged and the following
%    article shall be properly cited in any publication that uses results
%    generated by the software:
%
%      Magri C, Whittingstall K, Singh V, Logothetis NK, Panzeri S: A
%      toolbox for the fast information analysis of multiple-site LFP, EEG
%      and spike train recordings. BMC Neuroscience 2009 10(1):81;
%
% 3.  this notice shall remain in place in each source file.

persistent previous_Nt previous_mask;

% HiR, HiRS and ChiR are not computed for the gaussian case, null values
% are returned instead:
HiR  = NaN;
HiRS = NaN;
ChiR = NaN;

% We have to recompute Nc, Nt, maxNt (in the case this routine is called by
% quadratic extrapolation)
Nt    = pars.Nt;
Nc    = size(R,1);
maxNt = size(R,2);
totNt = sum(Nt);
Ns    = pars.Ns;

logOf2times2 = 2*log(2);
logOf2piTimesExp1 = log(2*pi*exp(1));

if ~isequal(previous_Nt, Nt)
  mask = build_logical_mask(Nt, maxNt, Ns);
else
	mask = previous_mask;
end
R(:,~mask) = nan;

% H(R) and H_lin(R) -------------------------------------------------------
if pars.doHR || pars.doHlR

    % Removing mean:
    %meanRvec = mean(R(:,mask),2);
		meanRvec = nansum(nansum(R,3),2)./totNt; %JM: this is faster, some juggling is needed to account for the -1's

	  %Rvec = R(:,mask) - meanRvec(:,ones(totNt,1));
    Rvec = R(:,mask) - meanRvec*ones(1,totNt);

    covPr = (Rvec * Rvec') ./ (totNt-1);
end

% H(R)
if pars.doHR
    HR = (Nc * logOf2piTimesExp1 + log(det(covPr))) / logOf2times2;

    % Bias correction:
    if pars.biasCorrNum==3
        HRbias = gaussian_bias(totNt, Nc);
        HR = HR - HRbias;
    end
else
    HR = 0;
end

% H_lin(R)
if pars.doHlR
    HlR = (Nc * logOf2piTimesExp1 + nansum(log(diag(covPr)), 1)) / logOf2times2;

    % Bias correction:
    if pars.biasCorrNum==3
        HlRbias = Nc * gaussian_bias(totNt, 1);
        HlR = HlR - HlRbias;
    end
else
    HlR = 0;
end

% H(R|S) and H_lin(R|S) ---------------------------------------------------
if pars.doHRS || pars.doHlRS

    % Removing mean
    meanR = nansum(R.*permute(mask(:,:,ones(Nc,1)), [3 1 2]), 2) ./ permute(Nt(:, ones(Nc,1)), [2 3 1]);
    Rzeromean = R - meanR(:, ones(maxNt,1), :);

    if Nc==1
        % HlRS is never computed for Nc==1, if we got here it's cause we
        % need to compute H(R|S):

				%covPrsMat = covPrsForSingletonNc(Rzeromean, Nt, mask, maxNt, Ns);
        covPrsMat = reshape(nansum(Rzeromean.*Rzeromean,2),1,[])./(Nt'-1); %JM:

        HRS  = (Nc * logOf2piTimesExp1 + log(covPrsMat(:))) / logOf2times2;

        HlRS = 0;

        % Bias correction:
        if pars.biasCorrNum==3
            HRSbias = gaussian_bias(Nt, 1);
            HRS = HRS - HRSbias;
        end
    else
        HRS  = zeros(Ns,1);
        HlRS = zeros(Ns,1);

        for s=1:Ns
            covPrs = Rzeromean(:, 1:Nt(s), s) * Rzeromean(:, 1:Nt(s), s).' ./ (Nt(s)-1);

            if pars.doHRS
                detCovPrs = det(covPrs);

                HRS(s) = (Nc * logOf2piTimesExp1 + log(detCovPrs)) / logOf2times2;

                % Bias correction:
                if pars.biasCorrNum==3
                    HRSbias = gaussian_bias(Nt(s), Nc);
                    HRS(s) = HRS(s) - HRSbias;
                end
            end;

            if pars.doHlRS
                diagCovPrs = diag(covPrs);

                HlRS(s) = (length(diagCovPrs) * logOf2piTimesExp1 + sum(log(diagCovPrs))) / logOf2times2;

                % Bias correction:
                if pars.biasCorrNum==3
                    HlRSbias = Nc * gaussian_bias(Nt(s), 1);
                    HlRS(s) = HlRS(s) - HlRSbias;
                end
            end
        end
    end

    if ~pars.doHRS
        HRS = 0;
    end

    if ~pars.doHlRS
        HlRS = 0;
    end

else
    HRS  = 0;
    HlRS = 0;
end

% H_sh(R) and H_sh(R|S) ---------------------------------------------------
if pars.doHshR || pars.doHshRS

    Rsh = shuffle_R_across_cells(R, Nt);

    if pars.doHshR
        % Removing mean:
        meanRshvec = mean(Rsh(:,mask),2);
        Rshvec = Rsh(:,mask) - meanRshvec(:,ones(totNt,1));

        covPshr = (Rshvec * Rshvec') ./ (totNt-1);

        HshR = (Nc * logOf2piTimesExp1 + log(det(covPshr))) / logOf2times2;

        % Bias correction:
        if pars.biasCorrNum==3
            HshRbias = gaussian_bias(totNt, Nc);
            HshR = HshR - HshRbias;
        end
    else
        HshR = 0;
    end

    % Remark: HshRS not computed for Nc==1
    if pars.doHshRS
        % Removing mean
        meanRsh = sum(Rsh.*permute(mask(:,:,ones(Nc,1)), [3 1 2]), 2) ./ permute(Nt(:, ones(Nc,1)), [2 3 1]);
        Rsh = Rsh - meanRsh(:,ones(maxNt,1),:);

        HshRS = zeros(Ns,1);
        for s=1:Ns
            covPshrs = Rsh(:, 1:Nt(s), s) * Rsh(:, 1:Nt(s), s).' ./ (Nt(s)-1);
            detCovPshrs = det(covPshrs);

            HshRS(s) = (Nc * logOf2piTimesExp1 + log(detCovPshrs)) / logOf2times2;
        end

        % Bias correction:
        if pars.biasCorrNum==3
            HshRSbias = gaussian_bias(Nt(:), Nc);
            HshRS = HshRS - HshRSbias;
        end
    else
        HshRS = 0;
    end

else
    HshR  = 0;
    HshRS = 0;
end

previous_Nt   = Nt;
previous_mask = mask;

function covPrs = covPrsForSingletonNc(R, Nt, mask, maxNt, Ns)

% Removing first singleton dimension:
R2D = zeros(maxNt, Ns);
R2D(mask) = R(1,mask);

% Computing the variance (R.*R faster than R.^2):
covPrs = (sum(R2D.*R2D,1)) ./ (Nt.' - 1);


text:
% DEMROBOTWIRELESSFGPLVM2 Wireless Robot data from University of Washington, without dynamics and without back constraints.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'robotWireless';
experimentNo = 2;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('ftc');
options.back = 'mlp';
options.backOptions = mlpOptions;

latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
modelWriteResult(model, dataSetName, experimentNo);

if exist('printDiagram') & printDiagram
  lvmPrintPlot(model, lbls, dataSetName, experimentNo);
end

% Load the results and display dynamically.
lvmResultsDynamic(model.type, dataSetName, experimentNo, 'vector')


text:
function datapt = lvd_AeroTasks(stateLogEntry, subTask, dummy)
%lvd_AeroTasks Summary of this function goes here
%   Detailed explanation goes here
    arguments
        stateLogEntry(1,1) LaunchVehicleStateLogEntry
        subTask(1,:) char
        dummy
    end

    ut = stateLogEntry.time;
    rVect = stateLogEntry.position;
    vVect = stateLogEntry.velocity;
    bodyInfo = stateLogEntry.centralBody;
    mass = stateLogEntry.getTotalVehicleMass();
    aero = stateLogEntry.aero;
    attState = stateLogEntry.attitude;
    celBodyData = bodyInfo.celBodyData;

    maStateLogEntry = [ut, rVect(:)', vVect(:)', bodyInfo.id, mass, 0, 0, 0, -1];
    altitude = ma_GALongLatAltTasks(maStateLogEntry, 'alt', celBodyData);
    vVectEcefMag = ma_GALongLatAltTasks(maStateLogEntry, 'bodyFixedVNorm', celBodyData);
    [lat, long, ~, ~, ~, ~, ~, ~] = getLatLongAltFromInertialVect(ut, rVect, bodyInfo, vVect);
    [density, pressureKPA, ~] = getAtmoDensityAtAltitude(bodyInfo, altitude, lat, ut, long);
    [~,angOfAttack,angOfSideslip,totalAoA] = attState.getAeroAngles(ut, rVect, vVect, bodyInfo);

    switch subTask
        case 'dragCoeff'
            datapt = stateLogEntry.aero.getDragCoeff(ut, rVect, vVect, bodyInfo, mass, altitude, pressureKPA, density, vVectEcefMag, totalAoA, angOfAttack, angOfSideslip);
        case 'dragForce'
            dragForceModel = DragForceModel();
            dragForceVect = dragForceModel.getForce(ut, rVect, vVect, mass, bodyInfo, aero, [], [], [], [], [], [], [], [], [], [], attState);
            dragForceMag = norm(dragForceVect); %mT*km/s^2
            datapt = dragForceMag*1000; %kN
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%plot_vehicle Draw ground vehicle pose
%
% PLOT_VEHICLE(X,OPTIONS) draws a representation of ground robot as an
% oriented triangle with pose X (1x3) = [x,y,theta] or X (3x3) as an SE(2)
% homogeneous transform.
%
% Options::
% 'scale',S    Draw vehicle with length S x maximum axis dimension (default
%              1/60)
% 'size',S     Draw vehicle with length S
%
% See also Vehicle.plot.

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

% TODO needs to work for 3D point

function plot_vehicle(x, varargin)

    opt.scale = 1/60;
    opt.size = [];

    [opt,args] = tb_optparse(opt, varargin);

    % get the current axes
    a = axis;

    % compute the dimensions of the robot
    if ~isempty(opt.size)
        d = opt.size;
    else
        d = (a(2)+a(4) - a(1)-a(3)) * opt.scale;
    end

    if numel(x) == 3
        % convert vector form of pose to SE(2)
        T = se2(x(1), x(2), x(3));
    else
        T = x;
    end

    % draw it
%     points = [
%         d 0 1
%         -d -0.6*d 1
%         -d 0.6*d 1
%         d 0 1]';

        points = [
        d 0
        -d -0.6*d
        -d 0.6*d]';

    points = homtrans(T, points);

    plot_poly(points, args{:});

    end

author: Allopart

repo: rbpf-gmapping


text:
function R = project2SO3(M)

if size(M,1) ~= 3 || size(M,2) ~= 3
    error('project2SO3 requires a 3x3 matrix as input')
end

[U,S,V] = svd(M);
R = U*V';
if(det(R)<0)
  R = U * diag([1 1 -1]) * V';
end

end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function p = mvnormpdfln(x, m, S, V)
% MVNORMPDFLN    log of multivariate normal density.
%   See MVNORMPDF for argument description.

log2pi = 1.83787706640935;
[d, n] = size(x);
if nargin == 1
  dx = x;
elseif isempty(m)
  dx = x;
else
  % m specified
  sz = size(m);
  if sz(1) ~= d
    error('rows(m) ~= rows(x)')
  end
  nm = sz(2);
  if nm == 1
    dx = x - repmat(m,1,n);
  elseif n == 1
    dx = repmat(x,1,nm) - m;
  elseif nm == n
    dx = x - m;
  else
    error('incompatible number of columns in x and m')
  end
end
if nargin < 3
  % unit variance
  p = -0.5*(d*log2pi + col_sum(dx.*dx));
  return
end
have_inv = 0;
if nargin == 3
  % standard deviation given
  if d == 1
    dx = dx./S;
    p = (-log(S) -0.5*log2pi) - 0.5*(dx.*dx);
    return;
  end
  if S(2,1) ~= 0
    error('S is not upper triangular')
  end
  if any(size(S) ~= [d d])
    error('S is not the right size')
  end
else
  if ischar(V)
    if strcmp(V,'inv')
      % inverse stddev given
      iS = S;
      have_inv = 1;
    else
      error('unknown directive')
    end
  elseif ischar(S)
    if strcmp(S,'inv')
      % inverse variance given
      if d == 1
	iS = sqrt(V);
      else
	iS = chol(V);
      end
      have_inv = 1;
    else
      error('unknown directive')
    end
  else
    % variance given
    if d == 1
      S = sqrt(V);
    else
      S = chol(V);
    end
  end
end
if have_inv
  if d == 1
    dx = iS .* dx;
    logdetiS = log(iS);
  else
    dx = iS*dx;
    logdetiS = sum(log(diag(iS)));
  end
else
  if d == 1
    dx = dx./S;
    logdetiS = -log(S);
  else
    dx = solve_tril(S',dx);
    %dx = S'\dx;
    logdetiS = -sum(log(diag(S)));
  end
end
p = (logdetiS -0.5*d*log2pi) -0.5*col_sum(dx.*dx);

author: andrewssobral

repo: lrslibrary


text:
fid = fopen('G:/100_25/zzzpairs.txt');
CC = fscanf(fid,'%d %d');
n_set = CC(1);n_num=CC(2);

same_pair = cell(n_set*n_num,2);
diff_pair = cell(n_set*n_num,2);
lfw_label = zeros(n_set*n_num * 2,2);
label_ht = java.util.Hashtable;

for i=1:n_set
    for j = 1 : n_num
        CC = textscan(fid, '%s %d %d\n');
        p = CC{1};id1=CC{2};id2=CC{3};
        same_pair((i-1)*n_num + j,1) = {sprintf('G:/100_25/%s_%04d-l.bmp',p{1},id1)};
        same_pair((i-1)*n_num + j,2) = {sprintf('G:/100_25/%s_%04d-l.bmp',p{1},id2)};
        if label_ht.containsKey(p{1})
            lfw_label((i-1)*n_num + j,1) = label_ht.get(p{1});
            lfw_label((i-1)*n_num + j,2) = label_ht.get(p{1});
        else
            idx = label_ht.size()+1;
            lfw_label((i-1)*n_num + j,1) = idx;
            lfw_label((i-1)*n_num + j,2) = idx;
            label_ht.put(p{1},idx);
        end;
    end;
    for j = 1 : n_num
         CC = textscan(fid, '%s %d %s %d\n');
         p1 = CC{1};id1=CC{2};p2=CC{3};id2=CC{4};
        diff_pair((i-1)*n_num + j,1) = {sprintf('G:/100_25/%s_%04d-l.bmp',p1{1},id1)};
        diff_pair((i-1)*n_num + j,2) = {sprintf('G:/100_25/%s_%04d-l.bmp',p2{1},id2)};
        if label_ht.containsKey(p1{1})
            lfw_label((i-1)*n_num + j + n_set*n_num,1) = label_ht.get(p1{1});
        else
            idx = label_ht.size()+1;
            label_ht.put(p1{1},idx);
            lfw_label((i-1)*n_num + j + n_set*n_num,1) = idx;
        end;
        if label_ht.containsKey(p2{1})
            lfw_label((i-1)*n_num + j + n_set*n_num,2) = label_ht.get(p2{1});
        else
            idx = label_ht.size()+1;
            label_ht.put(p2{1},idx);
            lfw_label((i-1)*n_num + j + n_set*n_num,2) = idx ;
        end;
    end;
end;
fclose(fid);

author: happynear

repo: FaceVerification


text:
function dat =  gen(algo)
%
%           data =  gen(algorithm)
%
%  Data will be generated from a generative model. This should be
%  equivalent to
%                    test(algo)
%  (no data is supplied to the test function)

  dat = generate(algo);

author: canlab

repo: CanlabCore


text:
function [y a]=lpredict(x, np, npred, pos)
% LPREDICT estimates the values of a data set before/after the observed set.
%
% LPREDICT uses linear prediction to extrapolate data, typically a
% timeseries. Note that this is not the same as linear extrapolation.
% A window of autocorrelation coefficients is moved beyond the data
% limits to extrapolate the data. For a discussion, see Press et. al. [1].
%
% The required coefficients are derived from a call to LPC in MATLAB's
% Signal Processing Toolbox
%
% Example:
% y=LPREDICT(x, np, npred, pos)
% [y, a]=LPREDICT(x, np, npred, pos)
%      x:       the input data series as a column vector or a matrix
%                   with series organized in columns
%      np:      the number of predictor coefficients to use (>=2)
%      npred:   the number of data values to return in the output
%      pos:     a string 'pre' or 'post' (default: post)
%                   This determines whether extrapolation occurs before or
%                   after the observed series x.
%
%      y:       the output, appropriately sequenced for concatenation with
%                   input x
%      a:       the coefficients returned by LPC (organized in rows).
%                   These can be used to check the quality/stability of the
%                   fit to the observed data as described in the
%                   documenation to the LPC function.
%
% The output y is given by:
%       y(k) = -a(2)*y(k-1) - a(3)*y(k-2) - ... - a(np)*y(k-np)
%                where y(n) => x(end-n) for n<=0
%
% Note that sum(-a(2:end))is always less than unity. The output will
% therefore approach zero as npred increases. This may be a problem if x
% has a large DC offset. Subtract the the column mean(s) of x from x on
% input and add them to the output column(s) to restore DC. For a more
% accurate DC correction, see [1].
%
% To pre-pend data, the input sequence is reversed and the output is
% similarly reversed before being returned. The output may always be
% vertically concatenated with the input to extend the timeseries e.g:
%       k=(1:100)';
%       x=exp(-k/100).*sin(k/5);
%       x=[lpredict(x, 5, 100, 'pre'); x; lpredict(x, 5, 100, 'post')];
%
%
% See also LPC
%
% References:
% [1] Press et al. (1992) Numerical Recipes in C. (Ed. 2, Section 13.6).
%
% Toolboxes Required: Signal Processing
%
% Revisions:    10.07 renamed to avoid filename clash with System ID
%                     Toolbox
%                     DC correction help text corrected.
%
% -------------------------------------------------------------------------
% Author: Malcolm Lidierth 10/07
% Copyright Š The Author & King's College London 2007
% -------------------------------------------------------------------------

% ---------------Argument checks---------------
if nargin<3
    error('Not enough input arguments');
end

if np<2
    error('np must be >=2');
end

if nargin<4
    pos='post';
end
%---------------------------------------------

%------------------MATRIX--------------------
% Deal with matrix input.
% Apply function to each column via recursive calls
cols=size(x,2);
if cols>1
    y=zeros(npred, cols);
    a=zeros(cols, np+1);
    for k=1:size(x,2)
        [y(:,k) a(k,:)]=lpredict(x(:,k), np, npred, pos);
    end
    return
end
%---------------------------------------------

% ---------------MAIN FUNCTION---------------
% Order input sequence
if nargin==4 && strcmpi(pos,'pre')
    x=x(end:-1:1);
end

% Get the forward linear predictor coefficients via the LPC
% function
try
    a=lpc(x,np);
catch
    % LPC missing?
    m=lasterror();
    if strcmp(m.identifier, 'MATLAB:UndefinedFunction')
        error('Requires the LPC function from the Signal Processing Toolbox');
    else
        rethrow(lasterror);
    end
end

% Negate coefficients, and get rid of a(1)
cc=-a(2:end);

% Pre-allocate output
y=zeros(npred,1);
% Seed y with the first value
y(1)=cc*x(end:-1:end-np+1);
% Next np-1 values
for k=2:min(np,npred)
    y(k)=cc*[y(k-1:-1:1); x(end:-1:end-np+k)];
end
% Now do the rest
for k=np+1:npred
    y(k)=cc*y(k-1:-1:k-np);
end

% Order the output sequence if required
if nargin==4 && strcmpi(pos,'pre')
    y=y(end:-1:1);
end

return
end