text:
function [P,varargout] = laplace(xs,A,b,AEq,bEq,varargin)
% LAPLACE Smooth compact 1,2 or 3 dimensional atomic probability distribution subject to linear and moment constraints
%
% LAPLACE('demo') or LAPLACE() runs univariate, bivariate and trivariate examples
%
% P = LAPLACE({xs1,xs2},[],[],[],[],[mu1 mu2],[mu11,mu12,mu22]) returns an
% matrix of size length(xs1) x length(xs2). The value P(i,j) is interpreted
% as a probability assigned to lattice point (xs1(i), xs2(j)). This 2d atomic
% probability minimizes the square of the discrete laplacian assuming zero
% boundary condition, subject to moment constraints E[x1]=mu1, E[x2]=mu2,
% E[x1.*x2] = mu11, E[x1.*x2]=mu12, E[x2.*x2]=mu22 where E denotes the expectation
% operator E(f) = sum(sum(f.*P)) and x1,x2,x3 are arrays of n-d coordinates
% function generated by [x1,x2,x3] = ndgrid(xs{:}).
%
% P = LAPLACE({xs1,xs2},A,b,AEq,bEq) returns a matrix of size length(xs1) x length(xs2).
% The value P(i,j) is interpreted as a probability assigned to lattice point (xs1(i), xs2(j)).
% The 2-d atomic probability minimizes the square of the discrete laplacian subject to inequality
% constraint A*P(:) < b and/or equality constraint AEq*P(:) = bEq. Supplied vectors
% b and bEq should be of length length(xs1) x length(xs2). A,b, AEq and bEq
% may be empty. In particular ...
%
% P = LAPLACE(xs1,[],[],[f1;f2;...],[Ef1;Ef2;...]) where f1,f2,... are row vectors
% of size 1 x length(xs1) representing functions on grid points xs1, and Ef1, Ef2
% are scalars representing their means returns a smooth compact probability
% distribution with E[f1] = Ef1, E[f2] = Ef2 and so forth.
%
% P = LAPLACE(xs1,[],[],[],[],mu1,mu11,mu111,mu1111) returns a
% vector of size length(xs1) x 1. The element P(i) is interpreted
% as a probability assigned to lattice point xs1(i). The 1-d atomic
% probability minimizes the square of the discrete laplacian assuming zero
% boundary condition, subject to moment constraints E[x1]=mu1,
% E[x1.*x1]=mu11, E[x1.*x1.*x1] = mu111 and E[x1.*x1.*x1.*x1]=mu1111.
% Addtional moment constraints may be supplied.
%
% P = LAPLACE({xs1,xs2,xs3},[],[],[],[],varargin) returns an
% matrix of size length(xs1) x length(xs2) x length(xs3). The value P(i,j,k) is interpreted
% as a probability assigned to lattice point (xs1(i), xs2(j), xs3(j)). The 3d atomic
% probability minimizes the sums squares of projections of the discrete laplacian assuming zero
% boundary condition, subject to moment constraints contained in varargin.
% Each constraint should be a vector of length NCHOOSEK(K+2,2) and contain
% moment values (as defined above) for all polynomials in x1, x2 and x3 of
% degree k sorted by reverse lexicographic ordering of their coefficients.
% For example, P = LAPLACE({xs1,xs2,xs3},[],[],[],[],[mu11,mu12,mu13,mu22,mu23,mu33])
% minimizes E[x1.^2] = mu11, E[x1.x2]= mu12 and so forth. The degree of the
% implied polynomials is inferred from the length of the moment vector.
%
% P = LAPLACE(...,A,b,AEq,bEq,...) applies additional user supplied
% inequality and equality constraints on P(:), as per the convention in
% most optimization functions including QUADPROG
%
%
% ARGUMENTS:
%   xs      n x 1 cell array  of vectors  - Lattice points at which the probability distribution is supported (compatible with ndgrid(xs{:}))
%   A,b,AEq,bEq                           - Inequality and equality constraints (can be empty), as per QUADPROG
%   varargin                              - Listing of desired values for moments (see explanation below)
%   P       m1 x m2 x ... x mn            - Lattice point probabilities (n-dimensional)
%   varargout                             - Evaluation of moments inferred by P in the same order as they are supplied
%
% MOMENT CONSTRAINTS:
%   If varargin{i} has length nchoosek(k+n-1,n-1) is it intepreted as a set of moment
%   constraints on polynomials of degreee k in the n coordinate functions, listed by
%   reverse lexicographic ordering. [A shortcoming of this convention is
%   that *all* degree k moments must be supplied].
%
%           if n=2 then [0.1 0.3] means E[x1]=0.1 and E[x2]=0.3
%           if n=2 then [0.5 0.2 0.1] means E[x1.^2] = 0.5, E[x1.*x2] = 0.2 and E[x2.*x2] = 0.1
%           if n=3 then [0.5 0.2 0.1] means E[x1] = 0.5, E[x2] = 0.2 and E[x3] = 0.1
%
% LIMITATIONS:
%   The discrete laplacian and implicit boundary condition may not be the most sensible
%   objective function, especially in the 3-d case as currently implemented (I'm hoping someone
%   improves this). One might also replace the atomic measure with mixture of gaussians
%   centred at the lattice points.
%
% RELATED:
%   ksdensity.m - stats toolbox
%
% Peter Cotton

momentConstraintMultiplier = 1;

if nargin==0 || ischar(xs) && strcmpi(xs,'demo'),
    laplaceDemo;
else

    %% Limitations and checks
    [m,n,xs] = laplaceArgChecks(xs,A,b,AEq,bEq,varargin{:});

    %% Establish grid points
    X = cell(n,1);
    switch n
        case 1,
            if ~iscell(xs),
                error('laplace: Violation of programmer intent - xs should have been set to cell by laplaceArgChecks');
            elseif iscell(xs),
                X = xs;
            else
                error('Not sure how to interpret xs');
            end
        otherwise
            [X{:}] = ndgrid(xs{:});
    end

    %% Ensure probs add to one.
    if isempty(AEq),
        AEq = ones(1,prod(m));
        bEq = 1;
    else
        AEq = [AEq;ones(1,prod(m))];
        bEq = [bEq;1];
    end

    %% Interpret additional arguments as equality constraints
    nsumk_ = nan(n,1);
    for k=1:10,
        nsumk_(k) = nsumk(n,k);
    end
    nArgs = length(varargin);
    moments = cell(nArgs,1);
    for argNo=1:length(varargin),
        arg_ = varargin{argNo};
        if isvector(arg_),
            switch n
                case 1,
                    degree_ = argNo;
                    expon = degree_;
                    nExpon= 1;
                otherwise
                    degree_ = find(nsumk_==length(arg_),1);
                    if isempty(degree_),
                        error([num2str(k),'th argument appears to have the wrong number of moments']);
                    else
                        [~,expon] = nsumk(n,degree_); % Listing of all non-negs adding to degree_
                        nExpon = size(expon,1);
                        expon = expon(nExpon:-1:1,:); % Reverse lexicographic ordering is more natural for moments
                    end
            end
            moments_ = nan(prod(m),nExpon);
            for i=1:nExpon,
                % Tabulate pointwise value of polynomial in x1, x2 and x3
                xi = X{1}(:);
                moments_(:,i) = xi.^expon(i,1);
                for j=2:n,
                    xj = X{j}(:);
                    moments_(:,i) = moments_(:,i).*(xj.^expon(i,j));
                end
            end
            AEq = [AEq;momentConstraintMultiplier*moments_'];
            bEq = [bEq;momentConstraintMultiplier*arg_(:)];
            moments{argNo} = moments_;
        end
    end

    %% Optimize discrete laplacian
    switch n
        case 1,
            f = @(p) del(p,xs);
        otherwise
            unrol = @(P) P(:);
            f = @(p) unrol(del(reshape(p,m'),xs));
    end
    [C,shouldBeZero] = inferAffineMap(f,prod(m));
    if any(abs(shouldBeZero)),
        error('laplace: Violation of programmer intent (bug). Discrete laplacian should not have constant term');
    end
    H = full(C'*C);
    options = optimset('LargeScale','off','MaxIter',500,'TolX',1e-2,'MaxFunEvals',1000);
    lb = zeros(prod(m),1);
    ub = ones(prod(m),1);
    p = quadprog(H,zeros(prod(m),1),A,b,AEq,bEq,lb,ub,zeros(prod(m),1),options);
    if n>1,
        P = reshape(p,m');
    else
        P = p;
    end

    %% Return moments implied by the atomic measure
    vargout = cell(length(varargin),1);
    for argNo=1:length(varargin),
        if ~isempty(moments{argNo}),
            vargout{argNo} = (p'*moments{argNo})';
        end
    end
    varargout = cell(length(varargin),1);
    [varargout{:}] = vargout{:};
end

end

function L = del(P,xs)
% function L = del(P,xs)
%
% Augmented discrete laplacian for 1 or 2 dimensions
% Projected, augmented discrete laplacian for 3 dimensions

% Tabulate grid steps
hs = cell(1,length(xs));
for dimNo=1:length(xs),
    hs{dimNo} = xs{dimNo}(2)-xs{dimNo}(1);
end

if isvector(P),
    L = del1augmented(P,hs);
elseif length(size(P))==2,
    L = del2augmented(P,hs);
elseif length(size(P))==3,
    L = del3augmented(P,hs);
else
    error('laplace is not implemented for dimension four'); % TODO:
end

    function d = del1augmented(P,hs)
        d = del2([0;0;P(:);0;0],hs{:});
    end

    function d = del2augmented(P,hs)
        % Enlarges the area and then computes del2
        P_ = zeros(size(P)+4);
        P_(3:end-2,3:end-2) = P;
        d = del2(P_,hs{:});
    end

    function d = del3augmented(P,hs)
        % Enlarges the area and then computes del2
        P_ = zeros(size(P)+4);
        P_(3:end-2,3:end-2,3:end-2) = P;
        d = del3(P_,hs{:});
    end
    function L = del3(X,varargin)
        % FIXME: For now this returns the laplacians of 2-d marginals
        L = zeros(0,1);
        for k=1:3,
            idx_ = setdiff(1:3,k);
            d_ = del2(squeeze(sum(X,k)),varargin{idx_});
            if k==1,
                L = d_(:);
            else
                L = [L;d_(:)];
            end
        end
    end

end

function [m,n,xs] = laplaceArgChecks(xs,A,b,AEq,bEq,varargin)
% function [n,m,xs] = laplaceArgChecks(xs,A,b,AEq,bEq,varargin)
%
% Ensure arguments make sense, and infer dimension where necessary

%% Allow passing of single vector for lattice rather than cell array
if ~iscell(xs) && isvector(xs),
    % Try to infer dimension
    if ~isempty(A),
        n_ = log(size(A,2))/log(length(xs));
        n = round(n_);
        if abs(n-n_)>1e-6,
            error('error infering dimension of problem, A looks to be inconsistent');
        end
    elseif ~isempty(AEq),
        n_ = log(size(AEq,2))/log(length(xs));
        n = round(n_);
        if abs(n-n_)>1e-6,
            error('error infering dimension of problem, AEq looks to be inconsistent');
        end
    else
        n = length(varargin{1});
        if n>3,
            error('inferred dimension exceeds 3 - better to be explicit');
        end
    end
    oldLattice = xs;
    xs = cell(n,1);
    for k1=1:n,
        xs{k1} = oldLattice;
    end
end

%% Check dimensions
n = length(xs);
m = nan(n,1);
for k2=1:n,
    latticeValues = xs{k2};
    m(k2) = length(latticeValues);
    if ~isvector(latticeValues) || any(any(diff(latticeValues)<eps)) || any(any(abs(diff(diff(latticeValues)))>1e-6)),
        error(['laplace anticipates an evenly spaced vector of support points in increasing order: e.g. (-4:0.5:3). Problem with lattice{',num2str(k2),'}']);
    end
end
if ~isempty(A) && size(A,2)~=prod(m),
    error('laplace anticipates size(A,2)=prod(m), where m(k) is the number of lattice points in dimension k');
end
if ~isempty(AEq) && size(AEq,2)~=prod(m),
    error('laplace anticipates size(AEq,2)=prod(m), where m(k) is the number of lattice points in dimension k');
end
if ~isempty(b) && size(b,1)~=size(bEq,1),
    error('laplace anticipates size(b,1)=prod(m), where m(k) is the number of lattice points in dimension k');
end
if ~isempty(bEq) && size(bEq,1)~=size(AEq,1),
    error('laplace anticipates size(bEq,1)=prod(m), where m(k) is the number of lattice points in dimension k');
end

end

function [A,b] = inferAffineMap(f,n)
% f  function_handle
% A  m x n
%
% Deduce affine mapping given its function handle

b = feval(f,zeros(n,1));
m = length(b);
A = sparse(m,n);
for k=1:n,
    ek = zeros(n,1);
    ek(k) = 1;
    A(:,k) = feval(f,ek)-b;
end
end

function [m,x] = nsumk(n,k)
% NSUMK Number and listing of non-negative integer n-tuples summing to k
%    M = NSUMK(N,K) where N and K are positive integers returns M=nchoosek(K+N-1,N-1)
%    This is the number of ordered N-tuples of non-negative integers summing to K
%
%    [M,X] = NSUMK(N,K) produces a matrix X with
%    nchoosek(K+N-1,N-1) rows and n columns. Each row comprises
%    non-negative integers summing to k. The ordering of rows follows the
%    same convention as NCHOOSEK, which is undocumented but with some
%    reliability appears to be lexicographic. The reverse of this presumed ordering
%    is a natural way to list coefficients of polynomials in N variables of degree K.
%    As per nchoosek, this syntax is only practical for situations where N is
%    less than about 15.
%
%  EXAMPLES:   m = nsumk(5,2)
%              [~,x] = nsumk(5,2) returns a 15 x 5 matrix x in which rows sum to 2

if isscalar(n) && isscalar(k) && nargout<=1
    m = nchoosek(k+n-1,n-1);
elseif isscalar(n) && isscalar(k) && nargout==2,
    m = nchoosek(k+n-1,n-1);
    dividers = [zeros(m,1),nchoosek((1:(k+n-1))',n-1),ones(m,1)*(k+n)];
    x = diff(dividers,1,2)-1;
else
    error('nsumk anticipates scalar k and n');
end
end

function [y,expon] = kmoments(x,k)
%  KMOMENTS - Evaluate sample mean of all k'th degree polynomials in N-variables
%    Y = KMOMENTS(X,K) where X is an NSAMPLES by N matrix produces a NSUMK(N,K) x 1
%    vector Y comprising sample means for all K-degree polynomials in N variables
%    Samples for the j'th variable comprise the j'th column of X, and the
%    implied listing of k-degree polynomials follows reverse lexicographic ordering
%
%    [Y,expon] = KMOMENTS(X,K) also returns the reverse-lex ordered polynomial exponents
%    in a matrix of size NCHOOSEK(K+N-1,N-1) x N

%% Generate reverse lex ordering
n = size(x,2);
switch n
    case 1,
        expon = k; m=1;
    otherwise
        [~,expon] = nsumk(n,k);
        m = size(expon,1);
        expon = expon(m:-1:1,:);
end

y = ones(m,1);

%% Compute k-th degree polys for a sample x
for i=1:m,
    ys = ones(size(x,1),1);
    for j=1:n,
        ys = ys.*x(:,j).^expon(i,j);
    end
    y(i) = mean(ys);
end

end

function laplaceDemo
close all;

%% UNIVARIATE EXAMPLE:
% An example motivated by finance
% We are inferring a risk-neutral density from market prices
% Prices of call options are linear constraints on the atomic density
% We throw in a few moment constraints as well
% Don't take this too seriously in the tails

matchExactMoments = true;  % <-- Set to false to verify the dangers in fitting the wrong moments
compareShortAndLongSample = false; % <-- Set to true to use moments and 'prices' from 'true' distribution
xs = (-5:0.05:5)';

% Generate an imagined risk-neutral density sample
shortSample_ = [-0.1+0.2*randn(10,1);0.3*randn(50,1);0.5+0.5*randn(100,1)];
if compareShortAndLongSample,
   longSample_ = [-0.1+0.2*randn(10*10000,1);0.3*randn(50*10000,1);0.5+0.5*randn(100*10000,1)];
else
   longSample_ = shortSample_;
end
pShort = hist(shortSample_,xs);pShort = pShort/sum(pShort);
pLong = hist(longSample_,xs);pLong = pLong/sum(pLong);

% Generate imagined market prices for securities based on the risk-neutral density sample
strikes = [-1:0.5:1];
security = @(x) max(0,x*ones(1,length(strikes))-ones(size(x,1),1)*strikes); % (Call option) Takes  column vector x of length nSamples ->  nSamples x nStrikes matrix
securityMarketPrice = mean(security(longSample_))';
securityPayoff = security(xs)';

% Infer the density subject to price and moment constraints
allMoments = cell(5,1);
for mNo=1:length(allMoments),
    if matchExactMoments,
        allMoments{mNo} = kmoments(longSample_,mNo);
    else
        allMoments{mNo} = kmoments(shortSample_,mNo);
    end
end
for nMomentsUsed = 0:5,
    P = laplace(xs,[],[],securityPayoff,securityMarketPrice,allMoments{1:nMomentsUsed});
    subplot(2,3,nMomentsUsed+1);plot(xs,P,xs,pLong,xs,pShort);legend('laplace','long sample','short sample');ylabel(['Matched ',num2str(nMomentsUsed), ' moments']);
    grid; plotxs = axis;axis([-2,2,plotxs(3),plotxs(4)]);
    drawnow; pause(1);
end

%% BIVARIATE EXAMPLE
% Inferring a bivariate distribution from moments
figure;
xs = {(-4:0.5:4)',(-4:0.5:4)};
sample_ = mvnrnd([-0.1 0.3],[1.2,0.5;0.5,0.98],500);
sample_(:,1) = sample_(:,1);
sample_(:,2) = -0.5+sample_(:,2)+0.5./(1+abs(sample_(:,1)));
pEmpirical = hist3(sample_,xs);
pEmpirical = pEmpirical/sum(pEmpirical(:));
moments1 = kmoments(sample_,1);
moments2 = kmoments(sample_,2);
moments3 = kmoments(sample_,3);
moments4 = kmoments(sample_,4);
[P,m1,m2,m3,m4] = laplace(xs,[],[],[],[],moments1,moments2,moments3,moments4);
subplot(1,2,1);surf(pEmpirical);zlabel('Empirical');
subplot(1,2,2);surf(P);zlabel('Laplace minimizing');
bivariateMomentComparison = [[moments1;moments2;moments3;moments4],[m1;m2;m3;m4]],
drawnow; pause(1);

%% TRIVARIATE EXAMPLE
% Inferring a tri-variate distribution from moments
figure;
disp('Trivariate example takes a little while...');
nTrivariateMomentsToUse = 3;
xs = {(-3:1:3)',(-3:1:3),(-3:1:3)};
sample_ = mvnrnd([-0.1 0.3 0.1],[1.2,0.5,0.2;0.5,0.98,0.35;0.2,0.35,1.2],50);
triMoments = cell(nTrivariateMomentsToUse,1);
for mNo=1:length(triMoments),
   triMoments{mNo}= kmoments(sample_,mNo);
end
P = laplace(xs,[],[],[],[],triMoments{:});
Pxy = squeeze(sum(P,3));
pxyEmpirical = hist3(sample_(:,[1:2]),xs(1:2)); pxyEmpirical = pxyEmpirical/sum(pxyEmpirical(:));
subplot(1,2,1);surf(pxyEmpirical);zlabel('Empirical (xy projection)');
subplot(1,2,2);surf(Pxy);zlabel('Laplace minimizing (xy projection)');

end


text:
function varargout = slogfrac(varargin)
%logfrac
%
% y = SLOGFRAC(x)
%
% Computes/declares log(1 + x(1)/x(2))

switch class(varargin{1})
    case 'double'
        x = varargin{1};
        % Safe version with defined negative values (helps fmincon when
        % outside feasible region)
        if all(x==0)
            varargout{1} = log(2);% ?definition...
        elseif x(1)==0
            varargout{1} = log(1);
        else
            aux = 1 + (x(1)./(x(2)+sqrt(eps)));
            if aux <= 0
                varargout{1} = -15;
            else
                varargout{1} = log(1 + (x(1)./(x(2)+sqrt(eps))));
            end

        end

    case 'sdpvar'
        if min(size(varargin{1}))>1
            error('SLOGFRAC only defined for vector arguments');
        else
            varargout{1} = yalmip('define',mfilename,varargin{1});
        end

    case 'char'
        X = varargin{3};
        Y = X(2);
        X = X(1);
        F = [];
        operator = struct('convexity','none','monotonicity','none','definiteness','none','model','callback');
        operator.range = [-inf inf];
        operator.domain = [-inf inf];
        operator.bounds = @bounds;
        operator.derivative = @(x) ([1./(x(1)+x(2)+sqrt(eps));1./(x(1)+x(2)+sqrt(eps))-(x(2)+sqrt(eps)).^-1]);
        varargout{1} = F;
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error('SDPVAR/SLOGFRAC called with CHAR argument?');
end

function [L,U] = bounds(xL,xU)

% Derive bounds on x1/x2
z = [xL(1)./xL(2) xU(1)./xL(2) xL(1)./xU(2) xU(1)./xU(2)];
fL = min(z);
fU = max(z);

if fL <= -1
    L = -inf;
else
    L = log(1 + fL);
end
if fU <= -1
    U = -inf;
else
    U = log(1 + fU);
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function deg=rad2deg(rad)
% RAD2DEG  Converts radians to decimal degrees. Vectorized.
% Version: 8 Mar 00
% Useage:  deg=rad2deg(rad)
% Input:   rad - vector of angles in radians
% Output:  deg - vector of angles in decimal degrees

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

deg=rad.*180./pi;
%ind=(deg<0);
%deg(ind)=deg(ind)+360;

author: Sable

repo: mcbench-benchmarks


text:
function binomvtot = genbinomtab(nmax,~)

%this function generates all the possible binomial coefficients of the form
% nchoosek(n,k) with n <= nmax note that only k <= n are considered
% n and k both start from zero but matlab arrays start from one hence
% nchoosek(n,k) = binomvtot(n+1,k+1)
binomvtot = zeros(nmax+1);

% if a second argument is given it will generate the log of this function

if nargin==1

 % This is no longer used as it isn't as good!

%  for n = 1:nmax  %note n=a but index is offset by 1
%     binomvtot(n+1,1:n+1) = [binomvtot(n,1:n)./(1:-1/n:1/n),1];
%  end

 %use instead n choose k = [(n-1) choose k-1 ] + [n choose k -1 ]
binomvtot(:,1)=1;
 for n = 2:nmax+1  %note n=a but index is offset by 1
        binomvtot(n,2:n) = binomvtot(n-1,2:n)+binomvtot(n-1,1:n-1);
 end

else %take logs, note that biomial coeffs are given by exp of this
     %this uses the recurrsion relation valid for k \le n-1
 % nchoosek(n,k) = nchoosek(n-1,k) *n/(n-k)    n/(n-k) = 1/(1-k/n)
 % and nchoosek(a,a) == 1
 for n = 1:nmax  %note n=a but index is offset by 1
    binomvtot(n+1,1:n) = binomvtot(n,1:n)-log(1:-1/n:1/n);
 end
end

 %uncomment to test, I have found the relative differences, i.e. elements
 %of (binomvtot - binomvtot2)./binomvtot2 to be of order eps
% tic
% binomvtot2 = zeros(nmax+1);
% for n = 0:nmax
%     for nn=0:n
%     binomvtot2(n+1,nn+1) = nchoosek(n,nn);
%     end
%  end
% toc
% diff = (binomvtot - binomvtot2)./binomvtot2;
% diff(isnan(diff))=0
% max(abs(diff))


text:
% take the covariance matrix of input trajectories

function output = F_logdet(input)

precision = class(gather(input(1)));
if ~strcmpi(precision, 'double')    % we need to use double precision
    input = double(input);
end

[D,M,N] = size(input);
if N==1
    output = log(det(input));
else
    % to be implemented
end

if strcmpi(precision, 'single')
    output = single(output);
end

end

author: singaxiong

repo: SignalGraph


text:
%% t_meshShowContrast
%
% Illustrates how to project a contrast map onto a 3D mesh
%
% Before proceeding with tutorial:
%	* Run t_initGrayAndVolume.
%	* Run t_glmRun.
%	* Run t_glmComputeContrastMap.
%	* Run t_convertMapInplaneToGray
%
% See also T_INITGRAYANDVOLUME, T_GLMRUN, T_GLMCOMPUTECONTRASTMAP, T_CONVERTMAPINPLANETOGRAY.
%
% JW (c) Stanford VISTA
%

%% Initialize the key variables and data path:
dataDir     = fullfile(mrvDataRootPath,'functional','vwfaLoc');
pathToMesh  = fullfile(mrvDataRootPath, 'anatomy','anatomyNIFTI', 'leftMesh.mat');
pathToMap   = fullfile(dataDir, 'Gray', 'GLMs', 'FixVWordScrambleWord.mat');
displayMode = 'co'; % Co-thresholded display mode
threshold   = .05; % Co-threshold

%% Retain original directory, change to data directory
curDir = pwd;
cd(dataDir);

%% Load data structure, parameter map, and 3D mesh
vw_gray = initHiddenGray();
vw_gray = loadParameterMap(vw_gray, pathToMap);
vw_gray = meshLoad(vw_gray, pathToMesh, 1);

%% Select display mode and threshold
vw_gray = viewSet(vw_gray, 'displaymode', 'co');
vw_gray = viewSet(vw_gray, 'cothresh', threshold);

%% Set bicolor colormap (neg and pos values)
vw_gray = bicolorCmap(vw_gray);

%% Overlay contrast map onto mesh
vw_gray = meshColorOverlay(vw_gray);

%% Restore original directory
cd(curDir);

%% END

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function gui_showmenu(varargin)
% Links the BCILAB menu into another menu, or creates a new root menu if necessary.
% gui_showmenu(Options...)
%
% In:
%   Options... : optional name-value pairs; names are:
%                 'parent': parent menu to link into
%
%                 'shortcuts': whether to enable keyboard shortcuts
%
%                 'forcenew': whether to force creation of a new menu
%
% Example:
%   % bring up the BCILAB main menu if it had been closed
%   gui_showmenu;
%
% See also:
%   env_startup
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-10-29

% parse options...
hlp_varargin2struct(varargin,'parent',[], 'shortcuts',true,'forcenew',false);

% check if we're an EEGLAB plugin
folders = hlp_split(fileparts(mfilename('fullpath')),filesep);
within_eeglab = length(folders) >= 5 && strcmp(folders{end-3},'plugins') && ~isempty(strfind(folders{end-4},'eeglab'));

% highlight the menu if already there
f = findobj('Tag','bcilab_menu');
if ~isempty(f) && ~forcenew
    for k=1:length(f)
        try
            close(get(f(k),'Parent'));
        catch e
            disp('Error: previous figure handle invalid. You are in debug mode. Press dbquit to move on.');
            keyboard;
        end
    end
end

if isempty(parent) %#ok<NODEF>
    if within_eeglab && ~forcenew
        % try to link into the EEGLAB main menu
        try
            toolsmenu = findobj(0,'tag','tools');
            if ~isempty(toolsmenu)
                parent = uimenu(toolsmenu, 'Label','BCILAB');
                set(toolsmenu,'Enable','on');
            end
        catch %#ok<CTCH>
            disp('Unable to link BCILAB menu into EEGLAB menu.');
        end
    end
    if isempty(parent)
        % create new root menu, if no parent
        from_left = 70;
        from_bottom = 600;
        width = 500;
        height = 1;
        % determine position on primary monitor
        import java.awt.GraphicsEnvironment
        ge = GraphicsEnvironment.getLocalGraphicsEnvironment();
        gd = ge.getDefaultScreenDevice();
        scrheight = gd.getDisplayMode().getHeight();
        pos = [from_left, from_bottom, width, height];
        % create figure
        release = version; try release = release(find(release=='(')+1:find(release==')')-1); catch end
        figtitle = sprintf('BCILAB %s (%s on %s)',env_version,release,hlp_hostname);
        parent = figure('DockControls','off','NumberTitle','off','Name',figtitle,'Resize','off','MenuBar','none','Position',pos,'Tag','bcilab_toolwnd');
    end
end

% Data Source menu
source = uimenu(parent, 'Label','Data Source','Tag','bcilab_menu');
uimenu(source,'Label','Load recording(s)...','Accelerator',char(shortcuts*'l'),'Callback','gui_loadset');
wspace = uimenu(source,'Label','Workspace','Separator','on');
uimenu(wspace,'Label','Load...','Callback','io_loadworkspace');
uimenu(wspace,'Label','Save...','Callback','io_saveworkspace');
uimenu(wspace,'Label','Clear...','Callback','clear');
uimenu(source,'Label','Run script...','Separator','on','Callback',@invoke_script);
if isdeployed
    uimenu(source,'Label','Quit','Separator','on','Callback','exit'); end

% Offline Analysis menu
offline = uimenu(parent, 'Label','Offline Analysis');
uimenu(offline,'Label','New approach...','Accelerator',char(shortcuts*'n'),'Callback','gui_newapproach');
uimenu(offline,'Label','Modify approach...','Accelerator',char(shortcuts*'m'),'Callback','gui_configapproach([],true);');
uimenu(offline,'Label','Review/edit approach...','Accelerator',char(shortcuts*'r'),'Callback','gui_reviewapproach([],true);');
uimenu(offline,'Label','Save approach...','Accelerator',char(shortcuts*'s'),'Callback','gui_saveapproach');
uimenu(offline,'Label','Train new model...','Accelerator',char(shortcuts*'t'),'Callback','gui_calibratemodel','Separator','on');
uimenu(offline,'Label','Apply model to data...','Accelerator',char(shortcuts*'a'),'Callback','gui_applymodel');
uimenu(offline,'Label','Visualize model...','Accelerator',char(shortcuts*'v'),'Callback','gui_visualizemodel');
uimenu(offline,'Label','Run batch analysis...','Accelerator',char(shortcuts*'b'),'Callback','gui_batchanalysis','Separator','on');
uimenu(offline,'Label','Review results...','Accelerator',char(shortcuts*'i'),'Callback','gui_selectresults','Separator','on');

% Online Analysis menu
online = uimenu(parent,'Label','Online Analysis');
pipe = uimenu(online,'Label','Process data within...');
read = uimenu(online,'Label','Read input from...');
write = uimenu(online,'Label','Write output to...');
cm_read = uicontextmenu('Tag','bcilab_cm_read');
cm_write = uicontextmenu('Tag','bcilab_cm_write');

% for each plugin sub-directory...
dirs = dir(env_translatepath('functions:/online_plugins'));
for d={dirs(3:end).name}
    % find all files, their names, identifiers, and function handles
    files = dir(env_translatepath(['functions:/online_plugins/' d{1} '/run_*.m']));
    names = {files.name};
    idents = cellfun(@(n)n(1:end-2),names,'UniformOutput',false);
    % for each entry...
    for f=1:length(idents)
        try
            if ~exist(idents{f},'file') && ~isdeployed
                addpath(env_translatepath(['functions:/online_plugins/' d{1}])); end
            % get properties...
            props = arg_report('properties',str2func(idents{f}));
            % get category
            if strncmp(idents{f},'run_read',8);
                cats = [read,cm_read];
            elseif strncmp(idents{f},'run_write',9);
                cats = [write,cm_write];
            elseif strncmp(idents{f},'run_pipe',8);
                cats = pipe;
            end
            if isfield(props,'name')
                % add menu entry
                for cat=cats
                    uimenu(cat,'Label',[props.name '...'],'Callback',['arg_guidialog(@' idents{f} ');'],'Enable','on'); end
            else
                warning('env_showmenu:missing_guiname','The online plugin %s does not declare a GUI name; ignoring...',idents{f});
            end
        catch %#ok<CTCH>
            disp(['Could not integrate the online plugin ' idents{f} '.']);
        end
    end
end
uimenu(online,'Label','Clear all online processing','Callback','onl_clear','Separator','on');

% Settings menu
settings = uimenu(parent, 'Label','Settings');
uimenu(settings,'Label','Directory settings...','Callback','gui_configpaths');
uimenu(settings,'Label','Cache settings...','Callback','gui_configcache');
uimenu(settings,'Label','Cluster settings...','Callback','gui_configcluster');
uimenu(settings,'Label','Clear memory cache','Callback','env_clear_memcaches','Separator','on');

% Help menu
helping = uimenu(parent,'Label','Help');
uimenu(helping,'Label','BCI Paradigms...','Callback','env_doc code/paradigms');
uimenu(helping,'Label','Filters...','Callback','env_doc code/filters');
uimenu(helping,'Label','Machine Learning...','Callback','env_doc code/machine_learning');
scripting = uimenu(helping,'Label','Scripting');
uimenu(scripting,'Label','File input/output...','Callback','env_doc code/io');
uimenu(scripting,'Label','Dataset editing...','Callback','env_doc code/dataset_editing');
uimenu(scripting,'Label','Offline scripting...','Callback','env_doc code/offline_analysis');
uimenu(scripting,'Label','Online scripting...','Callback','env_doc code/online_analysis');
uimenu(scripting,'Label','BCILAB environment...','Callback','env_doc code/environment');
uimenu(scripting,'Label','Cluster handling...','Callback','env_doc code/parallel');
uimenu(scripting,'Label','Keywords...','Separator','on','Callback','env_doc code/keywords');
uimenu(scripting,'Label','Helpers...','Callback','env_doc code/helpers');
uimenu(scripting,'Label','Internals...','Callback','env_doc code/utils');
authoring = uimenu(helping,'Label','Plugin authoring');
uimenu(authoring,'Label','Argument declaration...','Callback','env_doc code/arguments');
uimenu(authoring,'Label','Expression functions...','Callback','env_doc code/expressions');
uimenu(authoring,'Label','Online processing...','Callback','env_doc code/online_analysis');
uimenu(helping,'Label','About...','Separator','on','Callback',@about);
uimenu(helping,'Label','Save bug report...','Separator','on','Callback','io_saveworkspace([],true)');
uimenu(helping,'Label','File bug report...','Callback','arg_guidialog(@env_bugreport);');

% toolbar (if not linked into the EEGLAB menu)
if ~(within_eeglab && ~forcenew)
    global tracking; %#ok<TLEV>
    cluster_requested = isfield(tracking,'cluster_requested') && ~isempty(tracking.cluster_requested);
    cluster_requested = hlp_rewrite(cluster_requested,false,'off',true,'on');

    ht = uitoolbar(parent,'HandleVisibility','callback');
    uipushtool(ht,'TooltipString','Load recording(s)',...
        'CData',load_icon('bcilab:/resources/icons/file_open.png'),...
        'HandleVisibility','callback','ClickedCallback','gui_loadset');
    uipushtool(ht,'TooltipString','New approach',...
        'CData',load_icon('bcilab:/resources/icons/approach_new.png'),...
        'HandleVisibility','callback','ClickedCallback','gui_newapproach','Separator','on');
    uipushtool(ht,'TooltipString','Load Approach',...
        'CData',load_icon('bcilab:/resources/icons/approach_load.png'),...
        'HandleVisibility','callback','ClickedCallback',@load_approach);
    uipushtool(ht,'TooltipString','Save approach',...
        'CData',load_icon('bcilab:/resources/icons/approach_save.png'),...
        'HandleVisibility','callback','ClickedCallback','gui_saveapproach');
    uipushtool(ht,'TooltipString','Modify approach',...
        'CData',load_icon('bcilab:/resources/icons/approach_edit.png'),...
        'HandleVisibility','callback','ClickedCallback','gui_configapproach([],true);');
    uipushtool(ht,'TooltipString','Review/edit approach',...
        'CData',load_icon('bcilab:/resources/icons/approach_review.png'),...
        'HandleVisibility','callback','ClickedCallback','gui_reviewapproach([],true);');
    uipushtool(ht,'TooltipString','Train new model',...
        'CData',load_icon('bcilab:/resources/icons/model_new.png'),...
        'HandleVisibility','callback','ClickedCallback','gui_calibratemodel','Separator','on');
    uipushtool(ht,'TooltipString','Load Model',...
        'CData',load_icon('bcilab:/resources/icons/model_load.png'),...
        'HandleVisibility','callback','ClickedCallback',@load_model);
    uipushtool(ht,'TooltipString','Save Model',...
        'CData',load_icon('bcilab:/resources/icons/model_save.png'),...
        'HandleVisibility','callback','ClickedCallback','gui_savemodel');
    uipushtool(ht,'TooltipString','Apply model to data',...
        'CData',load_icon('bcilab:/resources/icons/model_apply.png'),...
        'HandleVisibility','callback','ClickedCallback','gui_applymodel');
    uipushtool(ht,'TooltipString','Visualize model',...
        'CData',load_icon('bcilab:/resources/icons/model_visualize.png'),...
        'HandleVisibility','callback','ClickedCallback','gui_visualizemodel');
    uipushtool(ht,'TooltipString','Run batch analysis',...
        'CData',load_icon('bcilab:/resources/icons/batch_analysis.png'),...
        'HandleVisibility','callback','ClickedCallback','gui_batchanalysis');
    uipushtool(ht,'TooltipString','Read input from (online)',...
        'CData',load_icon('bcilab:/resources/icons/online_in.png'),...
        'HandleVisibility','callback','Separator','on','ClickedCallback',@click_read);
    uipushtool(ht,'TooltipString','Write output to (online)',...
        'CData',load_icon('bcilab:/resources/icons/online_out.png'),...
        'HandleVisibility','callback','ClickedCallback',@click_write);
    uipushtool(ht,'TooltipString','Clear online processing',...
        'CData',load_icon('bcilab:/resources/icons/online_clear.png'),...
        'HandleVisibility','callback','ClickedCallback','onl_clear');
    uitoggletool(ht,'TooltipString','Request cluster availability',...
        'CData',load_icon('bcilab:/resources/icons/acquire_cluster.png'),'HandleVisibility','callback','Separator','on','State',cluster_requested,'OnCallback','env_acquire_cluster;','OffCallback','env_release_cluster;');
    uipushtool(ht,'TooltipString','About BCILAB',...
        'CData',load_icon('bcilab:/resources/icons/help.png'),'HandleVisibility','callback','Separator','on','ClickedCallback',@about);
end

if within_eeglab && forcenew
    mainmenu = findobj('Tag','EEGLAB');
    % make the EEGLAB menu current again
    if ~isempty(mainmenu)
        figure(mainmenu); end
end

function about(varargin)
infotext = strvcat(...
    'BCILAB is an open-source toolbox for Brain-Computer Interfacing research.', ...
    'It is being developed by Christian Kothe at the Swartz Center for Computational Neuroscience,',...
    'Institute for Neural Computation (University of California San Diego).', ...
    ' ',...
    'Development of this software was supported by the Army Research Laboratories under', ...
    'Cooperative Agreement Number W911NF-10-2-0022, as well as by a gift from the Swartz Foundation.', ...
    ' ',...
    'The design was inspired by the preceding PhyPA toolbox, written by C. Kothe and T. Zander', ...
    'at the Berlin Institute of Technology, Chair Human-Machine Systems.', ...
    ' ',...
    'BCILAB connects to the following toolboxes/libraries:', ...
    '* AWS SDK (Amazon)', ...
    '* Amica (SCCN/UCSD)', ...
    '* Chronux (Mitra Lab, Cold Spring Harbor)', ...
    '* CVX (Stanford)', ...
    '* DAL (U. Tokyo)', ...
    '* DataSuite (SCCN/UCSD)', ...
    '* EEGLAB (SCCN/UCSD)', ...
    '* BCI2000import (www.bci2000.org)', ...
    '* Logreg (Jan Drugowitsch)', ...
    '* FastICA (Helsinki UT)', ...
    '* glm-ie (Max Planck Institute for Biological Cybernetics, Tuebingen)', ...
    '* glmnet (Stanford)', ...
    '* GMMBayes (Helsinki UT)', ...
    '* HKL (Francis Bach, INRIA)', ...
    '* KernelICA (Francis Bach, Berkeley)', ...
    '* LIBLINEAR (National Taiwan University)', ...
    '* matlabcontrol (Joshua Kaplan)', ...
    '* mlUnit (Thomas Dohmke)', ...
    '* NESTA (Caltech)', ...
    '* OSC (Andy Schmeder) and LibLO (Steve Harris)', ...
    '* PROPACK (Stanford)', ...
    '* PropertyGrid (Levente Hunyadi)', ...
    '* SparseBayes (Vector Anomaly)', ...
    '* SVMlight (Thorsten Joachims)', ...
    '* SVMperf (Thorsten Joachims)', ...
    '* Talairach (UTHSCSA)', ...
    '* Time-frequency toolbox (CRNS / Rice University)', ...
    '* t-SNE (TU Delft)', ...
    '* UnLocBox (Nathanael Perraudin et al.)', ...
    '* VDPGM (Kenichi Kurihara)'); %#ok<REMFF1>

warndlg2(infotext,'About');

function click_read(varargin)
% pop up a menu when clicking the "read input from" toolbar button
tw = findobj('tag','bcilab_toolwnd');
cm = findobj('tag','bcilab_cm_read');
tpos = get(tw,'Position');
ppos = get(0,'PointerLocation');
set(cm,'Position',ppos - tpos(1:2), 'Visible', 'on');
set(cm,'Position',ppos - tpos(1:2), 'Visible', 'on');

function click_write(varargin)
% pop up a menu when clicking the "write output to" toolbar button
tw = findobj('tag','bcilab_toolwnd');
cm = findobj('tag','bcilab_cm_write');
tpos = get(tw,'Position');
ppos = get(0,'PointerLocation');
set(cm,'Position',ppos - tpos(1:2), 'Visible', 'on');
set(cm,'Position',ppos - tpos(1:2), 'Visible', 'on');

% run a script
function invoke_script(varargin)
[filename,filepath] = uigetfile({'*.m', 'MATLAB file'},'Select script to run',env_translatepath('bcilab:/'));
if ~isnumeric(filename)
    run_script([filepath filename],true); end

% load a toolbar icon
function cols = load_icon(filename)
[cols,palette,alpha] = imread(env_translatepath(filename));
if ~isempty(palette)
    error('This function does not handle palettized icons.'); end
cls = class(cols);
cols = double(cols);
cols = cols/double(intmax(cls));
cols([alpha,alpha,alpha]==0) = NaN;

% load a BCI model from disk
function load_model(varargin)
[filename,filepath] = uigetfile({'*.mdl', 'BCI Model'},'Select BCI model to load',env_translatepath('home:/.bcilab/models'));
if ~isnumeric(filename)
    contents = io_load([filepath filename],'-mat');
    for fld=fieldnames(contents)'
        tmp = contents.(fld{1});
        if isstruct(tmp) && isfield(tmp,'timestamp')
            tmp.timestamp = now;
            assignin('base',fld{1},tmp);
        end
    end
end

% load a BCI approach from disk
function load_approach(varargin)
[filename,filepath] = uigetfile({'*.apr', 'BCI Approach'},'Select BCI approach to load',env_translatepath('home:/.bcilab/approaches'));
if ~isnumeric(filename)
    contents = io_load([filepath filename],'-mat');
    for fld=fieldnames(contents)'
        tmp = contents.(fld{1});
        if isstruct(tmp) && isfield(tmp,'paradigm')
            tmp.timestamp = now;
            assignin('base',fld{1},tmp);
        end
    end
end

repo: REST


text:
function indx = multigrid_index_one ( dim_num, order_1d, order_nd )

%*****************************************************************************80
%
%% MULTIGRID_INDEX_ONE returns an indexed multidimensional grid.
%
%  Discussion:
%
%    For dimension DIM, the number of points is ORDER_1D(DIM).
%
%    We index the points as
%      1,   2,   3,   ...,  ORDER_1D(DIM).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension of the points.
%
%    Input, integer ORDER_1D(DIM_NUM), the order of the
%    rule in each dimension.
%
%    Input, integer ORDER_ND, the product of the entries of ORDER_1D.
%
%    Output, integer INDX(DIM_NUM,ORDER_ND), the indices of the points in
%    the grid.  The second dimension of this array is equal to the
%    product of the entries of ORDER_1D.
%
  a = [];
  more = 0;
  p = 0;

  while ( 1 )

    [ a, more ] = vec_colex_next2 ( dim_num, order_1d, a, more );

    if ( ~more )
      break
    end

    p = p + 1;
%
%  The values of A(DIM) are between 0 and ORDER_1D(DIM)-1 = N - 1 = 2 * M.
%  Subtracting M sets the range to -M to +M, as we wish.
%
    indx(1:dim_num,p) = a(1:dim_num) + 1;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [fix_up,fix_down] = presolve_fixvariables(A,b,c,lb,ub,monotinicity)

% These are optimally (or w.l.o.g) set to upper bound
not_in_obj = find(c<=0);
% Setting to 1 makes Ax<b increase feasible set for these variables
constrained_blow = all(-A(:,not_in_obj)>=0,1);
% and they enter via a psd matrix in all sdp constraints
sdp_positive = monotinicity(not_in_obj) == -1;
% these variables satisffy all constraints
can_fix = not_in_obj(find(constrained_blow & sdp_positive));
% these variables are still not fixed
still_on = find(lb==0 & ub==1);
% so we can fix these
fix_up = intersect(can_fix,still_on);

% These are optimally (or w.l.o.g) set to lower bound
not_in_obj = find(c>=0);
% Setting to 1 makes Ax<b increase feasible set for these variables
constrained_blow = all(A(:,not_in_obj)>=0,1);
% and they enter via a psd matrix in all sdp constraints
sdp_positive = monotinicity(not_in_obj) == 1;
% these variables satisffy all constraints
can_fix = not_in_obj(find(constrained_blow & sdp_positive));
% these variables are still not fixed
still_on = find(lb==0 & ub==1);
% so we can fix these
fix_down = intersect(can_fix,still_on);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% Test file for @classicfun/max.m

function pass = test_max(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebfunpref();
end

% Set a domain for BNDFUN.
data.domain = [-2 7];

%%
% Spot-check the extrema for a few BNDFUN.
pass(1) = test_spotcheck_max(@(x) sin(10*x), data, 1, pref);
pass(2) = test_spotcheck_max(@airy, data, 0.535656656015700, pref);
pass(3) = test_spotcheck_max(@(x) -1./(1 + x.^2), data, -.02, pref);
pass(4) = test_spotcheck_max(@(x) (x/10).^3.*cosh(x/10), ...
    data, 0.7^3*cosh(0.7), pref);

%%
% Check operation for array-valued BNDFUN inputs.
fun_op = @(x) [sin(10*x) airy(x) (x/10).^3.*cosh(x/10)];
f = bndfun(fun_op, data, pref);
[y, x] = max(f);
exact_max = [1 0.535656656015700 0.7^3*cosh(0.7)];
fx = [sin(10*x(1)) airy(x(2)) (x(3)/10).^3.*cosh(x(3)/10)];
tol = 10*get(f, 'vscale')*eps;
pass(5) = (all(abs(y - exact_max) < 10*tol) && ...
    all(abs(fx - exact_max) < tol));

%%
% Test for complex-valued BNDFUN.
pass(6) = test_spotcheck_max(@(x) (x/2).*(exp(1i*(x/2))+1i*sin(x/2)), ...
    data, -3.277598405517787 - 2.455482593827339i, pref);

fun_op = @(x) [((x-2).^2/4+1).*exp(1i*(x/2)) ...
    -((x+1).^2/4+1).*exp(1i*(x/2))];
f = bndfun(fun_op, data, pref);
[y, x] = max(f);
exact_max = [-6.789310982858273-2.543178400749744i 15.919763683943538+5.963314870723537i];
fx = [((x(1)-2).^2/4+1).*exp(1i*(x(1)/2)) ...
    -((x(2)+1).^2/4+1).*exp(1i*(x(2)/2))];
tol = get(f, 'vscale')*eps;
pass(7) = (all(abs(y - exact_max) < 10*tol) && ...
    all(abs(fx - exact_max) < tol));

%% Test for UNBNDFUN:

% Functions on [a inf]:

% Set the domain:
data.domain = [1 Inf];

op = @(x) x.*exp(-x);
f = unbndfun(op, data);
[y, x] = max(f);
yExact = exp(-1);
xExact = 1;
errY = y - yExact;
errX = x - xExact;
pass(8) = norm([errY errX], inf) < 1e3*eps*get(f,'vscale');

end

%%
% Spot-check the results for a given BNDFUN.
function result = test_spotcheck_max(fun_op, data, exact_max, pref)

f = bndfun(fun_op, data, pref);
[y, x] = max(f);
fx = fun_op(x);
tol = 100*get(f, 'vscale')*eps;

result = (all(abs(y - exact_max) < tol) && ...
          all(abs(fx - exact_max) < tol));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [imgw, imgwr, map] = tpswarp(img, outDim, Zp, Zs, interp)
%
% Description: Thin-Plane spline warping of the input image (img) to
% output image (imgw). The warping is defined by a set of reference points
% (Zp=[Xp Yp]) on the [img] image and a set of corresponding points (Zs)
% on the (imgw) image. The warping will translate Zp in img to Zs in imgw.
%
% Input:
% img - input image
% outDim - Output canvas ([W H])
% Zp - landmark in img
% Zs - landmark in canvas
% interp.method - interpolation mode('nearest', 'invdist', 'none')
% interp.radius - Max radius for nearest neighbor interpolation or
%                 Radius of inverse weighted interpolation
% interp.power - power for inverse weighted interpolation
%
% Output:
% imgw - warped image with no holes
% imgwr - warped image with holes
% map - Map of the canvas with 0 indicating holes and 1 indicating pixel
%
% Reference:
% F.L. Bookstein, "Principal Warps: Thin-Plate splines and the
% decomposition of deformations", IEEE Transaction on Pattern Analysis and
% Machine Intelligence, Vol. 11, No. 6, June 1989
%
% Author: Fitzgerald J Archibald
% Date: 07-Apr-09

%% Initialization
NPs = size(Zp,1); % number of landmark points

imgH = size(img,1); % height
imgW = size(img,2); % width

outH = outDim(2);
outW = outDim(1);

% landmark in input
Xp = Zp(:,1)';
Yp = Zp(:,2)';

% landmark in output (homologous)
Xs = Zs(:,1)';
Ys = Zs(:,2)';

%% Algebra of Thin-plate splines

% Compute thin-plate spline mapping [W|a1 ax ay] using landmarks
[wL]=computeWl(Xp, Yp, NPs);
wY = [Xs(:) Ys(:); zeros(3,2)]; % Y = ( V| 0 0 0)'   where V = [G] where G is landmark homologous (nx2) ; Y is col vector of length (n+3)
wW = inv(wL)*wY; % (W|a1 ax ay)' = inv(L)*Y

% Thin-plate spline mapping (Map all points in the plane)
% f(x,y) = a1 + ax * x + ay * y + SUM(wi * U(|Pi-(x,y)|)) for i = 1 to n
[Xw, Yw]=tpsMap(wW, imgH, imgW, Xp, Yp, NPs);

%% Warping

% input grid for warping
[X Y] = meshgrid(1:imgH,1:imgW); % HxW

% Nearest neighbor or inverse distance weighted based interpolation
[imgw,imgwr,map] = interp2d(X(:), Y(:), img, Xw, Yw, outH, outW, interp);

return

%% [L] = [[K P];[P' 0]]
% np - number of landmark points
% (xp, yp) - coordinate of landmark points
function [wL]=computeWl(xp, yp, np)

rXp = repmat(xp(:),1,np); % 1xNp to NpxNp
rYp = repmat(yp(:),1,np); % 1xNp to NpxNp

wR = sqrt((rXp-rXp').^2 + (rYp-rYp').^2); % compute r(i,j)

wK = radialBasis(wR); % compute [K] with elements U(r)=r^2 * log (r^2)
wP = [ones(np,1) xp(:) yp(:)]; % [P] = [1 xp' yp'] where (xp',yp') are n landmark points (nx2)
wL = [wK wP;wP' zeros(3,3)]; % [L] = [[K P];[P' 0]]

return

%% Mapping: f(x,y) = a1 + ax * x + ay * y + SUM(wi * U(|Pi-(x,y)|)) for i = 1 to n
% np - number of landmark points
% (xp, yp) - coordinate of landmark points
function [Xw, Yw]=tpsMap(wW, imgH, imgW, xp, yp, np)

[X Y] = meshgrid(1:imgH,1:imgW); % HxW
X=X(:)'; % convert to 1D array by reading columnwise (NWs=H*W)
Y=Y(:)'; % convert to 1D array (NWs)
NWs = length(X); % total number of points in the plane

% all points in plane
rX = repmat(X,np,1); % Np x NWs
rY = repmat(Y,np,1); % Np x NWs

% landmark points
rxp = repmat(xp(:),1,NWs); % 1xNp to Np x NWs
ryp = repmat(yp(:),1,NWs); % 1xNp to Np x NWs

% Mapping Algebra
wR = sqrt((rxp-rX).^2 + (ryp-rY).^2); % distance measure r(i,j)=|Pi-(x,y)|

wK = radialBasis(wR); % compute [K] with elements U(r)=r^2 * log (r^2)
wP = [ones(NWs,1) X(:) Y(:)]'; % [P] = [1 x' y'] where (x',y') are n landmark points (nx2)
wL = [wK;wP]'; % [L] = [[K P];[P' 0]]

Xw  = wL*wW(:,1); % [Pw] = [L]*[W]
Yw  = wL*wW(:,2); % [Pw] = [L]*[W]

return

%% k=(r^2) * log(r^2)
function [ko]=radialBasis(ri)

r1i = ri;
r1i(find(ri==0))=realmin; % Avoid log(0)=inf
ko = 2*(ri.^2).*log(r1i);

return


text:
function value = file_line_width ( file_name )

%*****************************************************************************80
%
%% FILE_LINE_WIDTH reports the length of the longest line in a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_NAME, the name of the input file.
%
%    Output, integer VALUE, the length of the longest line.
%
  value = -1;

  input = fopen ( file_name );

  if ( input < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_LINE_WIDTH - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', file_name );
    error ( 'FILE_LINE_WIDTH - Error!' );
    return;
  end

  value = 0;
  line_width = 0;

  while ( 1 )

    [ c, count ] = fread ( input, 1 );

    if ( count < 1 )
      break;
    end

    if ( c == 10 || c == 13 )
      line_width = 0;
    else
      line_width = line_width + 1;
      value = max ( value, line_width );
    end

  end

  fclose ( input );

  return
end


text:
function g = sigmoidGradient(z)
%SIGMOIDGRADIENT returns the gradient of the sigmoid function
%evaluated at z
%   g = SIGMOIDGRADIENT(z) computes the gradient of the sigmoid function
%   evaluated at z. This should work regardless if z is a matrix or a
%   vector. In particular, if z is a vector or matrix, you should return
%   the gradient for each element.

g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the gradient of the sigmoid function evaluated at
%               each value of z (z can be a matrix, vector or scalar).

g = sigmoid(z) .* (1 - sigmoid(z));

% =============================================================

end

author: atinesh-s

repo: Coursera-Machine-Learning-Stanford


text:
%................................................................

% MATLAB codes for Finite Element Analysis
% problem16vibrationsSchultz.m
% Timoshenko beam in free vibrations
% Lee/Schultz problem
% antonio ferreira 2008

% clear memory
clear all

% E; modulus of elasticity
% G; shear modulus
% I: second moments of area
% L: length of beam
% thickness: thickness of beam
E=10e7; poisson = 0.30;L  = 1;thickness=0.1;
I=thickness^3/12;
EI=E*I;
kapa=5/6;
rho=1;
A=1*thickness;
%

P = -1; % uniform pressure
% constitutive matrix
G=E/2/(1+poisson);
C=[   EI   0; 0    kapa*thickness*G];

% mesh
numberElements     = 40;
nodeCoordinates=linspace(0,L,numberElements+1);
xx=nodeCoordinates';x=xx';
for i=1:size(nodeCoordinates,2)-1
    elementNodes(i,1)=i;
    elementNodes(i,2)=i+1
end
% generation of coordinates and connectivities
numberNodes=size(xx,1);

% GDof: global number of degrees of freedom
GDof=2*numberNodes;

% computation of the system stiffness, force, mass
[stiffness,force,mass]=...
    formStiffnessMassTimoshenkoBeam(GDof,numberElements,...
    elementNodes,numberNodes,xx,C,P,rho,I,thickness);

% boundary conditions (CC)
fixedNodeW =find(xx==min(nodeCoordinates(:))...
    | xx==max(nodeCoordinates(:)));
fixedNodeTX=fixedNodeW;
prescribedDof=[fixedNodeW; fixedNodeTX+numberNodes];

% solution
displacements=solution(GDof,prescribedDof,stiffness,force);

% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness,...
    GDof,prescribedDof)

% free vibration problem
activeDof=setdiff([1:GDof]',[prescribedDof]);
[V,D]=eig(stiffness(activeDof,activeDof),...
    mass(activeDof,activeDof));
D = diag(sqrt(D)*L*L*sqrt(rho*thickness/E/I));D=sqrt(D) ;
[D,ii] = sort(D);

% lee,schultz paper
modeNumber=4;
V1=zeros(GDof,1);
V1(activeDof,1:modeNumber)=V(:,1:modeNumber);

% drawing eigenmodes
drawEigenmodes1D(modeNumber,numberNodes,V1,xx,x)

repo: Aero-matlab


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

%% -*- texinfo -*-
%% @defun valid_sym_assumptions ()
%% Return list of valid assumptions.
%%
%% @end defun

%% Reference: http://docs.sympy.org/dev/modules/core.html

%% FIXME: Maybe in the future remove this file and replace with something else.

function L = valid_sym_assumptions()

  persistent List

  if (isempty(List))

    cmd = {'from sympy.core.assumptions import _assume_defined'
           'return list(_assume_defined),'};
    List = pycall_sympy__ (cmd);

  end

  L = List;

end

repo: octsympy


text:
function [Btuph] = hpe2Btuph(hpe)
% Convert power from electric horsepower to British
% thermal units per hour.
% Chad A. Greene 2012
Btuph = hpe*2545.457658313;


text:
function output = ctranspose(input)

output = cellfun(@ctranspose, input, 'uniformoutput', false);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Gets the 3D coordinates of the corners of a 3D bounding box.
%
% Args:
%   bb3d - 3D bounding box struct.
%
% Returns:
%   corners - 8x3 matrix of 3D coordinates.
%
% See:
%   create_bounding_box_3d.m
%
% Author: Nathan Silberman (silberman@cs.nyu.edu)
function corners = get_corners_of_bb3d(bb3d)
  corners = zeros(8, 3);

  % Order the bases.
  [~, inds] = sort(abs(bb3d.basis(:,1)), 'descend');
  basis = bb3d.basis(inds, :);
  coeffs = bb3d.coeffs(inds);

  [~, inds] = sort(abs(basis(2:3,2)), 'descend');
  if inds(1) == 2
    basis(2:3,:) = flipdim(basis(2:3,:), 1);
    coeffs(2:3) = flipdim(coeffs(2:3), 2);
  end

  % Now, we know the basis vectors are orders X, Y, Z. Next, flip the basis
  % vectors towards the viewer.
  basis = flip_towards_viewer(basis, repmat(bb3d.centroid, [3 1]));

  coeffs = abs(coeffs);

  corners(1,:) = -basis(1,:) * coeffs(1) + basis(2,:) * coeffs(2) + basis(3,:) * coeffs(3);
  corners(2,:) = basis(1,:) * coeffs(1) + basis(2,:) * coeffs(2) + basis(3,:) * coeffs(3);
  corners(3,:) = basis(1,:) * coeffs(1) + -basis(2,:) * coeffs(2) + basis(3,:) * coeffs(3);
  corners(4,:) = -basis(1,:) * coeffs(1) + -basis(2,:) * coeffs(2) + basis(3,:) * coeffs(3);

  corners(5,:) = -basis(1,:) * coeffs(1) + basis(2,:) * coeffs(2) + -basis(3,:) * coeffs(3);
  corners(6,:) = basis(1,:) * coeffs(1) + basis(2,:) * coeffs(2) + -basis(3,:) * coeffs(3);
  corners(7,:) = basis(1,:) * coeffs(1) + -basis(2,:) * coeffs(2) + -basis(3,:) * coeffs(3);
  corners(8,:) = -basis(1,:) * coeffs(1) + -basis(2,:) * coeffs(2) + -basis(3,:) * coeffs(3);

  corners = corners + repmat(bb3d.centroid, [8 1]);
end

function normals = flip_towards_viewer(normals, points)
  points = points ./ repmat(sqrt(sum(points.^2, 2)), [1, 3]);

  proj = sum(points .* normals, 2);

  flip = proj > 0;
  normals(flip, :) = -normals(flip, :);
end

author: thusiyuan


text:
function [INI,S,E] = detectmuscle(S, iter, Mode)
% Muscle detection with inverse filtering
% Artifacts are indicated with NaN's.
%
% [INI,S,E] = detectmuscle(S [, iter [,1]])
% [INI,S,E] = detectmuscle(S, Fs, Mode)  with Mode>1
% [INI,S,E] = detectmuscle(S, arg2, Mode)
%
% iter		number of iterations [default:1]
% INI.MU	mean of S
% INI.InvFilter	coefficients of inverse filter
% S		outlier replaced by NaN
% E		isnan(E) indicates muscle artifact
% Mode	1: [default] inverse filtering
%	2: based on gradient, range and amplitude [BrainVision method]
% 	3: slope > 11 uV/sample [1]
% 	4: beta2 > 0.9 uV^2/Hz [1]
%
%
% References:
% [1]	Van de Velde, M., Van Erp, G., Cluitmans, P., 1998.
%	Detection of muscle artefact in the normal human awake EEG.
% 	Electroencephalography and Clinical Neurophysiology 107 (2), 149-158.
%

%	$Id: detectmuscle.m 2202 2009-10-27 12:06:45Z schloegl $
%	Copyright (C) 2003,2008,2009 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/
%
%    BioSig 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.
%
%    BioSig 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 BioSig.  If not, see <http://www.gnu.org/licenses/>.

if nargin<2,
        iter=[];
end;
INI = [];

% Muscle Detection
if Mode==1,
	%% TODO: Validation
	if isempty(iter) iter = 1; end;
	%% inverse filter
	TH = 5;
	[se,INI.MU] = sem(S);
	if TH*se<abs(INI.MU),
		[S] = center(S);
	end;

	while iter>0,
        	[AR,RC,PE]= lattice(S',10);
		INI.InvFilter = ar2poly(AR);
        	E = zeros(size(S));
        	for k = 1:size(S,2),
                	E(:,k) = filter(INI.InvFilter(k,:),1,S(:,k));
	        end;
        	INI.V  = std(E);
	        INI.TH = INI.V * TH;
        	iter   = iter-1;

	        for k = 1:size(S,2),
        	        S(E(:,k)>(INI.TH(k)) | E(:,k)<(-INI.TH(k)),k) = NaN;
	        end;
	end;
	% the following part demonstrates a possible correction
	for k = 1:size(S,2),
                E(:,k) = filter(INI.InvFilter(k,:),1,S(:,k));
	end;
	% isnan(E), returns Artifactselse

elseif Mode==2,
	Fs = iter;
	% Criterion for bad gradient:
	% Maximal allowed voltage step / sampling point: 100.00 µV
	% Mark as bad before event: 100.00 ms
	% Mark as bad after event: 100.00 ms
	ix1 = [abs(diff(S,[],1))>100;zeros(1,size(S,2))];

	% Criterion for bad max-min:
	% Maximal allowed absolute difference: 150.00 µV
	% Interval Length: 100.00 ms
	% Mark as bad before event: 100.00 ms
	% Mark as bad after event: 100.00 ms

	ix2 = abs(filter([-1;zeros(Fs/10-1,1);1],1,S))>150;
	ix2 = [ix2(Fs/20+1:end,:);zeros(Fs/20+1,size(S,2))];

	% Criterion for bad amplitude:
	% Minimal allowed amplitude: -75.00 µV
	% Maximal allowed amplitude: 75.00 µV
	% Mark as bad before event: 100.00 ms
	% Mark as bad after event: 100.00 ms
	ix3 = abs(S)>75;

	ix = filtfilt(ones(Fs/10,1),1,double(ix1|ix2|ix3))>0;
	S(ix) = NaN;
	INI=[];

elseif Mode==3,
	Fs = iter;
	E  = filter(ones(1,Fs),Fs,abs(diff(S))/Fs*250);
	Threshold = 11; % uV/sample [1]
	S(E>Threshold) = NaN;

elseif Mode==4,
	%% TODO: free butter
	Fs = iter;
	if Fs>= 250,
		[A,B]=butter(5,25*2/Fs,'high');
	else
		[A,B]=butter(5,[25,125]*2/Fs);
	end;
	E = filtfilt(ones(1,Fs),Fs,filtfilt(A,B,S).^2);
	Threshold = 0.9; % uV^2/Hz ??? [1]
 	S(E>Threshold) = NaN;

end;

if nargout<3, return; end;

author: PatternRecognition

repo: OpenBMI


text:
function h = times(f,g)
% .*    Pointwise multiplication for SPHEREFUN objects.
%
%   F.*G multiplies SPHEREFUN objects F and G. Alternatively F or G could be a
%   double.

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

if ( isa(f, 'spherefun') && isa(g, 'spherefun') )
    % Grady's faster times for rank 1 functions:
    if ( length( f ) == 1 )
        [C, D, R] = cdr( f );
        h = g;
        onesForC = sqrt(abs(D))*ones(1,length(g));
        onesForR = sign(D)*onesForC;
        h.cols = (C*onesForC).*g.cols;
        h.rows = (R*onesForR).*g.rows;
        % Switch the parity terms if needed:
        % plus*plus = plus -> do nothing
        % plus*minus = minus -> do nothing
        % minus*minus = plus -> switch minus to plus
        % minus*plus = minus -> switch plus to minus
        if ( isempty( f.idxPlus ) )  % f is a minus
            h.idxMinus = g.idxPlus;
            h.idxPlus = g.idxMinus;
        end
        % Both terms have to be non-zero at the poles for the product to be
        % non-zero at the poles.
        h.nonZeroPoles = f.nonZeroPoles & g.nonZeroPoles;
    elseif ( length( g ) == 1 )
         h = times(g, f);
    else
        % Let separableApprox deal with the multiplication
        h = times@separableApprox(f,g);
    end
else
    % Let separableApprox deal with the multiplication
    h = times@separableApprox(f,g);
end

end


text:
function [sta_index, orid_index, eqtimes, eqindices, nstations, norigins] = arrivals_rearrange(arrivals, stations)
%% Arrange info from earthquake catalogue into a more useful format

% We want a column for each station and a row for each earthquake (origin).
% Store arrival times in this format for now, so we can later replace them
% with loaded waveforms.

% Number of stations.  Use user-specified list and order (rather than
% unique values from arrivals table) for greater flexibility.
nstations = length(stations);

% Number of earthquakes
norigins = length(arrivals.unique_orids);

% Allocate an array to store arrival times
eqtimes = zeros(norigins,nstations);

% Allocate an additional array to store the indices, so we can retrieve
% other data (e.g. event type) later
eqindices = zeros(norigins,nstations);

% Now go through each arrival and place it in the array
for n = 1:length(arrivals.orid)
	sta_index = find(strcmp(arrivals.sta{n},stations),1);
	orid_index = find(unique_orids == arrivals.orid(n),1); % Could do something more efficient here; orid is sorted in ascending order so a find is not really necessary
	eqtimes(orid_index,sta_index) = arrivals.atime(n);
	eqindices(orid_index,sta_index) = n;
end

author: geoscience-community-codes

repo: GISMO


text:
function pass = test_minus( pref )
% This tests the basic arithmetic operations on chebfun2 objects.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e5 * pref.cheb2Prefs.chebfun2eps;
j = 1;

D = [-1 1 -1 1; -2 2 -2 2; -1 pi 0 2*pi];

for r = 1 : size(D,1)
    f = chebfun2(@(x,y) cos(x.*y), D(r,:));
    g = chebfun2(@(x,y) x + y + x.*y, D(r,:));

    FminusG = chebfun2(@(x,y) cos(x.*y) - (x + y + x.*y), D(r,:) );

    tolr = norm(D(r,:),inf)*tol;

    pass(j) = ( norm( (f-g) - FminusG ) < 10*tolr ); j = j + 1;

end

end


text:
function c = solver(a)
%
% Copyright 2003 The MathWorks, Inc.

c = [1 1];