text:
function value = r8_besi1 ( x )

%*****************************************************************************80
%
%% R8_BESI1 evaluates the Bessel function I of order 1 of an R8 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the Bessel function I of order 1 of X.
%
  persistent bi1cs
  persistent nti1
  persistent xmax
  persistent xmin
  persistent xsml

  if ( isempty ( nti1 ) )

    bi1cs = [ ...
      -0.19717132610998597316138503218149E-02, ...
      +0.40734887667546480608155393652014, ...
      +0.34838994299959455866245037783787E-01, ...
      +0.15453945563001236038598401058489E-02, ...
      +0.41888521098377784129458832004120E-04, ...
      +0.76490267648362114741959703966069E-06, ...
      +0.10042493924741178689179808037238E-07, ...
      +0.99322077919238106481371298054863E-10, ...
      +0.76638017918447637275200171681349E-12, ...
      +0.47414189238167394980388091948160E-14, ...
      +0.24041144040745181799863172032000E-16, ...
      +0.10171505007093713649121100799999E-18, ...
      +0.36450935657866949458491733333333E-21, ...
      +0.11205749502562039344810666666666E-23, ...
      +0.29875441934468088832000000000000E-26, ...
      +0.69732310939194709333333333333333E-29, ...
      +0.14367948220620800000000000000000E-31 ]';

    nti1 = r8_inits ( bi1cs, 17, 0.1 * r8_mach ( 3 ) );
    xmin = 2.0 * r8_mach ( 1 );
    xsml = sqrt ( 8.0 * r8_mach ( 3 ) );
    xmax = log ( r8_mach ( 2 ) );

  end

  y = abs ( x );

  if ( y <= xmin )
    value = 0.0;
  elseif ( y <= xsml )
    value = 0.5 * x;
  elseif ( y <= 3.0 )
    value = x * ( 0.875 + r8_csevl ( y * y / 4.5 - 1.0, bi1cs, nti1 ) );
  elseif ( y <= xmax )
    value = exp ( y ) * r8_besi1e ( x );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_BESI1 - Fatal error!\n' );
    fprintf ( 1, '  Result overflows.\n' );
    error ( 'R8_BESI1 - Fatal error!' )
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ centers, mincenter, mindist, q2, quality ] = kmeans_fast ( ...
  data, initcenters, method )

%*****************************************************************************80
%
%% KMEANS_FAST carries out the fast KMEANS algorithm of Charles Elkan.
%
%  Discussion:
%
%    Note that if the input quantity INITCENTERS is a scalar, it is taken
%    to be the number of clusters desired, but if it is a KxM array, it is]
%    taken to be the M-dimensional coordinates of K points, to be used as
%    initial guesses for the cluster centers.
%
%  Modified:
%
%    04 September 2013
%
%  Author:
%
%    Charles Elkan
%
%  Reference:
%
%    Charles Elkan,
%    Using the Triangle Inequality to Accelerate k-Means,
%    Proceedings of the Twentieth International Conference on Machine Learning
%    (ICML-2003),
%    Washington DC, 2003.
%
%  Parameters:
%
%    Input, real DATA(N,M), the M-dimensional coordinates of N points.
%
%    Input, integer INITCENTERS, the number K of clusters desired.
%
%    Input, real INITCENTERS(K,M), the M-dimensional coordinates of K
%    points, which the program will use as starting guesses for the
%    cluster centers.
%
%    Input, integer METHOD, selects the algorithm to be used.
%    * 0, unoptimized, using n by k matrix of distances O(nk) space;
%    * 1, vectorized, using only O(n+k) space;
%    * 2, like 1, in addition using distance inequalities (default).
%
%    Output, real CENTERS(K,M), the M-dimensional coordinates of the K
%    cluster centers.
%
%    Output, integer MINCENTER(N,1), the index, for each data point, of
%    the cluster to which it belongs.
%
%    Output, real MINDIST(N,1), an upper bound of the distance of each point to
%    the nearest center.
%
%    Output, real Q2, the mean of UDIST^2.
%
%    Output, real QUALITY, the mean of UDIST.

tic
if nargin < 3 method = 2; end
[n,dim] = size(data);

if max(size(initcenters)) == 1
    k = initcenters;
    [centers, mincenter, mindist, lower, computed] = anchors(mean(data),k,data);
    total = computed;
    skipestep = 1;
else
    centers = initcenters;
    mincenter = zeros(n,1);
    total = 0;
    skipestep = 0;
    [k,dim2] = size(centers);
    if dim ~= dim2 error('dim(data) ~= dim(centers)'); end;
end

nchanged = n;
iteration = 0;
oldmincenter = zeros(n,1);

while nchanged > 0
    % do one E step, then one M step
    computed = 0;

    if method == 0 & ~skipestep
        for i = 1:n
            for j = 1:k
                distmat(i,j) = calcdist(data(i,:),centers(j,:));
            end
        end
        [mindist,mincenter] = min(distmat,[],2);
        computed = k*n;

    elseif (method == 1 | (method == 2 & iteration == 0)) & ~skipestep
        mindist = Inf*ones(n,1);
        lower = zeros(n,k);
        for j = 1:k
           jdist = calcdist(data,centers(j,:));
           lower(:,j) = jdist;
           track = find(jdist < mindist);
           mindist(track) = jdist(track);
           mincenter(track) = j;
        end
        computed = k*n;

    elseif method == 2 & ~skipestep
        computed = 0;
%
% for each center, nndist is half the distance to the nearest center
% if d(x,center) < nndist then x cannot belong to any other center
% mindist is an upper bound on the distance of each point to its nearest center
%
        nndist = min(centdist,[],2);
% the following usually is not faster
%        ldist = min(lower,[],2);
%        mobile = find(mindist > max(nndist(mincenter),ldist));
        mobile = find(mindist > nndist(mincenter));

% recompute distances for point i and center j
%       only if j can possibly be the new nearest center
% for speed, the first check has been optimized by modifying centdist
% swapping the order of the checks is slower for data with natural clusters

        mdm = mindist(mobile);
        mcm = mincenter(mobile);

        for j = 1:k
% the following is incorrect: for j = unique(mcm)'
            track = find(mdm > centdist(mcm,j));
            if isempty(track) continue; end
            alt = find(mdm(track) > lower(mobile(track),j));
            if isempty(alt) continue; end
            track1 = mobile(track(alt));
%
% calculate exact distances to the mincenter
% recalculate separately for each jj to avoid copying too much of data
% redo may be empty, but we don't need to check this.
%
            redo = find(~recalculated(track1));
            redo = track1(redo);
            c = mincenter(redo);
            computed = computed + size(redo,1);
            for jj = unique(c)'
                rp = redo(find(c == jj));
                udist = calcdist(data(rp,:),centers(jj,:));
                lower(rp,jj) = udist;
                mindist(rp) = udist;
            end
            recalculated(redo) = 1;

            track2 = find(mindist(track1) > centdist(mincenter(track1),j));
            track1 = track1(track2);
            if isempty(track1) continue; end

            % calculate exact distances to center j
            track4 = find(lower(track1,j) < mindist(track1));
            if isempty(track4) continue; end
            track5 = track1(track4);
            jdist = calcdist(data(track5,:),centers(j,:));
            computed = computed + size(track5,1);
            lower(track5,j) = jdist;

            % find which points really are assigned to center j
            track2 = find(jdist < mindist(track5));
            track3 = track5(track2);
            mindist(track3) = jdist(track2);
            mincenter(track3) = j;
        end % for j=1:k
    end % if method

    oldcenters = centers;
%
% M step: recalculate the means for each cluster
% if a cluster is empty, its mean is left unchanged
% we minimize computations for clusters with little changed membership
%
    diff = find(mincenter ~= oldmincenter);
    diffj = unique([mincenter(diff);oldmincenter(diff)])';
    diffj = diffj(find(diffj > 0));

    if size(diff,1) < n/3 & iteration > 0
         for j = diffj
            plus = find(mincenter(diff) == j);
            minus = find(oldmincenter(diff) == j);
            oldpop = pop(j);
            pop(j) = pop(j) + size(plus,1) - size(minus,1);
            if pop(j) == 0 continue; end
            centers(j,:) = (centers(j,:)*oldpop + sum(data(diff(plus),:),1) - sum(data(diff(minus),:),1))/pop(j);
        end
    else
        for j = diffj
            track = find(mincenter == j);
            pop(j) = size(track,1);
            if pop(j) == 0 continue; end
% it's correct to have mean(data(track,:),1) but this can make answer worse!
            centers(j,:) = mean(data(track,:),1);
        end
    end

    if method == 2
        for j = diffj
            offset = calcdist(centers(j,:),oldcenters(j,:));
            computed = computed + 1;
            if offset == 0 continue; end
            track = find(mincenter == j);
            mindist(track) = mindist(track) + offset;
            lower(:,j) = max(lower(:,j) - offset,0);
        end
%
% compute distance between each pair of centers
% modify centdist to make "find" using it faster.
%
        recalculated = zeros(n,1);
        realdist = alldist(centers);
        centdist = 0.5*realdist + diag(Inf*ones(k,1));
        computed = computed + k + k*(k-1)/2;
    end

    nchanged = size(diff,1) + skipestep;
    iteration = iteration+1;
    skipestep = 0;
    oldmincenter = mincenter;

%   difference = max(max(abs(oldcenters - centers)));
%   [iteration toc nchanged computed size(diffj,2)]
    fprintf ( 1, '%4d  %g  %d  %d\n', iteration, toc, nchanged, computed );
%   [iteration toc nchanged computed]
    total = total + computed;
end % while nchanged > 0

  udist = calcdist(data,centers(mincenter,:));
  quality = mean(udist);
  q2 = mean(udist.^2);
  %[iteration toc quality q2 total]
  fprintf ( 1, '  %4d  %g  %g  %g  %d\n', iteration, toc, quality, q2, total );

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

  return
end


text:
function pot = convert_to_pot(CPD, pot_type, domain, evidence)
% CONVERT_TO_POT Convert a root CPD to one or more potentials
% pots = convert_to_pot(CPD, pot_type, domain, evidence)

assert(length(domain)==1);
assert(~isempty(evidence(domain)));
T = 1;

sz = CPD.sizes;
ns = zeros(1, max(domain));
ns(domain) = sz;

switch pot_type
 case 'u',
  pot = upot(domain, 1, T, 0);
 case 'd',
  ns(domain) = 1;
  pot = dpot(domain, ns(domain), T);
 case {'c','g'},
  ns(domain) = 0;
  pot = cpot(domain, ns(domain), 0);
 case 'cg',
  ddom = [];
  cdom = domain; % we assume the root node is cts
  %pot = cgpot(ddom, cdom, ns, {cpot([],[],0)});
  pot = cgpot(ddom, cdom, ns);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [logp,logq1,logq2] = spm_mci_switch (Pr,M,U,Y,beta)
% Return log probability of tempered model switch
% FORMAT [logp,logq1,logq2] = spm_mci_switch (Pr,M,U,Y,beta)
%
% Pr        parameters (vectorised and in M.V subspace)
% M,U,Y     as usual
% beta      inverse temperature (set to 1 to get usual posterior)
%
% logp      log prob of model switch
% logq1     log joint of model 1
% logq2     log joint of model 2
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mci_switch.m 6548 2015-09-11 12:39:47Z will $

logq1 = spm_mci_joint(Pr,M{1},U{1},Y);

% Parameters in original space
P = M{1}.V*Pr+M{1}.vpE;
[L2,tmp,st] = feval(M{2}.L,P,M{2},U{2},Y);
e = P-M{2}.vpE;
L1 = - e'*M{2}.ipC*e/2 + M{2}.log_prior_t2;
logq2 = L1+L2;

logp=(1-beta)*logq1+beta*logq2;

author: spm


text:
% Estimate spatial covariance matrix for sentences using a mask. The mask
% specifies speech presense probability at all time frequency locations,
% with a 1 means speech present and 0 means speech absent.
%
function curr_layer = F_MVDR_spatialCov(input_layer, curr_layer)
input = input_layer.a;
fs = curr_layer.fs;
freqBin = curr_layer.freqBin;
nFreqBin = length(freqBin);

[D,T,N] = size(input);
if isfield(curr_layer, 'noiseCovL2')
    noiseCovL2 = curr_layer.noiseCovL2;
else
    noiseCovL2 = 0;  % add noiseCovRegularization*\lambda*I to noise covariance, where \lambda is the maximum eigenvalue
end

D = D/2;
speechCov = input(1:D,:,:,:);
noiseCov = input(D+1:end,:,:,:);

dimTmp = size(speechCov,1) / nFreqBin;
nCh = sqrt(dimTmp);

speechCov = reshape(speechCov, nCh, nCh, nFreqBin, T, N);
noiseCov = reshape(noiseCov, nCh, nCh, nFreqBin, T, N);

speechCov_cell = num2cell(speechCov, [1 2]);       % convert to cell array and call cellfun for speed
noiseCov_cell = num2cell(noiseCov, [1 2]);

% add regularization
if noiseCovL2 > 0
    eig_val = cellfun(@GetEigVal, noiseCov_cell, 'UniformOutput', 0);
    eig_val = cell2mat(eig_val);
    noise_floor = eig_val(end,:,:,:,:) * noiseCovL2;
    noise_floor(:) = 100;
    % because our noise floor depends on noise covariance and network
    % parameters, it will not be able to pass through gradient check. This
    % is because we didn't consider this dependancy in the backpropagation.
    % But this should not have much effect on the training.
    noise_floor_cell = num2cell(noise_floor, 1);
    noiseCov_cell = cellfun(@DiagLoading, noiseCov_cell, noise_floor_cell, 'UniformOutput', 0);
end

ninv_x = cellfun(@(x,n) (inv( n )*x), speechCov_cell, noiseCov_cell, 'UniformOutput', 0);
lambda = cellfun(@(x) abs(trace(x)), ninv_x, 'UniformOutput', 0);
if 0    % only lambda
    weight = cellfun(@(x,y) ones(size(x,1),1)/y, ninv_x, lambda, 'UniformOutput', 0);
elseif 0    % no lambda
    weight = cellfun(@(x) x(:,1), ninv_x, 'UniformOutput', 0);
else
    weight = cellfun(@(x,y) x(:,1)/y, ninv_x, lambda, 'UniformOutput', 0);
end
output = cell2mat(weight);
output = permute(output, [3 1 2 4 5]);
output = reshape(output, nFreqBin*nCh, T, N);

curr_layer.a = output;
curr_layer.lambda = lambda;
curr_layer.phi_s = speechCov_cell;
curr_layer.phi_n = noiseCov_cell;
if noiseCovL2>0
    curr_layer.noise_floor = noise_floor_cell;
end
end

%%
function ev = GetEigVal(A)
[~,V] = eig(A);
ev = diag(V);
end

function B = DiagLoading(A, alpha)
D = size(A,1);
B = A + alpha * eye(D);
end

author: singaxiong

repo: SignalGraph


text:
function [rs,rp]=reflectionCoeffsVec(e1,e2,k,n,u1,u2)
%%REFLECTIONCOEFFSVEC Reflection coefficients for the reflection of an
%           electromagnetic plane wave from a LOSSLESS medium off a LOSSY
%           medium (or another lossless medium) are computed given a vector
%           in the direction of propagation of the light and a vector
%           normal to the surface. Both media must have real
%           permeabilities. For example, one might approximate the
%           troposphere as lossless while accounting for the loss in sea
%           water. Lossy media have complex permittivities. This function
%           is a different parameterization of reflectionCoeffs.
%
%INPUTS: e1 The permittivity (refraction index) of the lossless medium.
%           The incoming ray is traveling through this prior to reflecting
%           off of the surface with permittivity e2. This must be a real
%           quantity (lossless). e1 and e2 can both be either absolute
%           permittivities or relative permittivities. "Relative" means
%           that the permittivity of the medium has been divided by the
%           permittivity of free space and is a dimensionless quantity.
%        e2 The permittivity of the medium against which the ray reflects.
%           This can be complex (lossy).
%         k A 3X1 vector in the direction of propagation of the light
%           approaching the surface from which it will reflect.
%         n A 3X1 normal vector to the surface of reflection. It does not
%           matter whether n is pointing up or down from the surface.
%        u1 The permeability of the lossless medium. This must be
%           a real quantity.  u1 and u2 can both be either absolute
%           permeabilities or relative permeabilities. If omitted or an
%           empty matrix is passed, a value of 1 is used (the permeability
%           equals that of free space). That is a reasonable approximation
%           for air.
%        u2 The permeability of the lossy medium against which the incoming
%           ray reflects. medium. This must be a real quantity. If omitted
%           or an empty matrix is passed, a value of 1 is used. That is a
%           reasonable approximation for water.
%
%OUTPUTS: rs The complex reflection coefficient for s-polarized light, also
%            known as transverse-Electric (TE) polarized light.
%         rp The complex reflection coefficient for p-polarized light, also
%            known as transverse-Magnetic (TM) polarized light and as
%            tangent-plane polarized light.
%
%This function uses angBetweenVecs to get the angle between the incoming
%vector and the normal to the surface (adjusting in case one should have
%used -n instead of n) and then calls reflectionCoeffs.
%
%July 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<6||isempty(u2))
    u2=1;
end

if(nargin<5||isempty(u1))
    u1=1;
end

thetai=angBetweenVecs(k,n);

if(thetai>pi/2)
    %If -n should have been used instead of n.
    thetai=pi-thetai;
end
[rs,rp]=reflectionCoeffs(e1,e2,thetai,u1,u2);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function s = spmatrix(a)
%SPMATRIX Converts a two-way sparse tensor to sparse matrix.
%
%   SPMATRIX(X) converts a sparse tensor to a sparse matrix. The sparse
%   tensor must be two-dimensional.
%
%   See also SPTENSOR, SPTENSOR/RESHAPE, SPTENMAT
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if ndims(a) ~= 2
    error('Sparse tensor must be two dimensional.');
end

if isempty(a.subs)
    s = sparse(a.size(1), a.size(2));
else
    s = sparse(a.subs(:,1), a.subs(:,2), a.vals, a.size(1), a.size(2));
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function seedRNG(s)
%SEEDRNG   Seed the MATLAB random number generator.
%   SEEDRNG(S) seeds the MATLAB random number generator with the integer value
%   S.  This function is meant to provide an interface for seeding the RNG that
%   is independent of the currently installed version of MATLAB.  It is used
%   solely for seeding the RNG in tests which use random numbers to ensure
%   deterministic behavior between runs.
%
%   Users should NOT call this function from their own code, as it may be
%   removed without warning in a future release.

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

if ( verLessThan('matlab', '7.12') )
    % Before R2011a.
    rand('seed', s);
    randn('seed', s);
else
    % R2011a and later
    rng('default');
    rng(s);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [posterior] = tapas_linear_estimate(y, x, model, inference, pars)
%% Estimate a linear model
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

n = 3;
if nargin < n
    model = struct();
end

n = n + 1;
if nargin < n
    inference = struct();
end

n = n + 1;
if nargin < n
    pars = struct();
end

[pars] = tapas_linear_pars(pars);
[data] = tapas_linear_data(y, x, pars);
[model] = tapas_linear_model(model, pars);
[inference] = tapas_linear_inference(inference, pars);

[posterior] = tapas_linear_estimate_interface(data, model, inference);

end

author: translationalneuromodeling

repo: tapas


text:
function outputImage = times(this, otherImage)
% Multiplies data of other image with this one and writes output in new image
%
% NOTE: If voxel dimensions of 2nd image do not match - or a scalar is
% given as 2nd argument - data in the 2nd argument is automatically
% replicated to match this image geometry.
%
%
%   Y = MrImage()
%   outputImage = times(Y, otherImage, ...
%   functionHandle)
%
%   OR
%
%   outputImage = Y.*otherImage
%
% This is a method of class MrImage.
%
%
% IN
%   otherImage              image that will be multiplied with this one
%
% OUT
%   outputImage             new MrImage, product of this and otherImage
%
% EXAMPLE
%
%   % Compute difference of 2 images
%		Y = MrImage();
%		Z = MrImage();
%		X = Y.times(Z);
%
%   % OR (cool overload!):
%       X = Y.*Z
%
%   See also MrImage perform_binary_operation MrImage.plus

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

outputImage = this.perform_binary_operation(otherImage, @times);


text:
function y = vl_nnconcat(inputs, dim, dzdy, varargin)
%VL_NNCONCAT CNN concatenate multiple inputs.
%  Y = VL_NNCONCAT(INPUTS, DIM) concatenates the inputs in the cell
%  array INPUTS along dimension DIM generating an output Y.
%
%  DZDINPUTS = VL_NNCONCAT(INPUTS, DIM, DZDY) computes the derivatives
%  of the block projected onto DZDY. DZDINPUTS has one element for
%  each element of INPUTS, each of which is an array that has the same
%  dimensions of the corresponding array in INPUTS.

% Copyright (C) 2015 Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.inputSizes = [] ;
opts = vl_argparse(opts, varargin, 'nonrecursive') ;

if nargin < 2, dim = 3; end;
if nargin < 3, dzdy = []; end;

if isempty(dzdy)
  y = cat(dim, inputs{:});
else
  if isempty(opts.inputSizes)
    opts.inputSizes = cellfun(@(inp) [size(inp,1),size(inp,2),size(inp,3),size(inp,4)], inputs, 'UniformOutput', false) ;
  end
  start = 1 ;
  y = cell(1, numel(opts.inputSizes)) ;
  s.type = '()' ;
  s.subs = {':', ':', ':', ':'} ;
  for i = 1:numel(opts.inputSizes)
    stop = start + opts.inputSizes{i}(dim) ;
    s.subs{dim} = start:stop-1 ;
    y{i} = subsref(dzdy,s) ;
    start = stop ;
  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
% STD_LOADALLEEG - constructs an ALLEEG structure, given the paths and file names
%                    of all the EEG datasets that will be loaded into the ALLEEG
%                    structure. The EEG datasets may be loaded without their EEG.data
%                    (see the POP_EDITOPTIONS function), so many datasets can be
%                    loaded simultaneously. The loaded EEG datasets have dataset
%                    information and a (filename) pointer to the data.
% Usage:
%                  % Load sseveral EEG datasets into an ALLEEG structure.
%                  >> ALLEEG = std_loadalleeg(paths,datasets) ;
%                  >> ALLEEG = std_loadalleeg(STUDY) ;
% Inputs:
%   paths      - [cell array of strings] cell array with all the datasets paths.
%   datasets   - [cell array of strings] cell array with all the datasets file names.
%
% Output:
%   ALLEEG     - an EEGLAB data structure, which holds the loaded EEG sets
%                 (can also be one EEG set).
%  Example:
%          >> paths = {'/home/eeglab/data/sub1/','/home/eeglab/data/sub2/', ...
%          >>          '/home/eeglab/data/sub3/', '/home/eeglab/data/sub6/'};
%          >> datasets = { 'visattS1', 'visattS2', 'visattS3', 'visattS4'};
%          >> ALLEEG = std_loadalleeg(paths,datasets) ;
%
% See also: POP_LOADSTUDY, POP_STUDY
%
% Authors: Hilit Serby, Arnaud Delorme, SCCN, INC, UCSD, October , 2004

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

function ALLEEG = std_loadalleeg(varargin)

    if nargin < 1
        help std_loadalleeg;
        return;
    end

    genpath = '';
    oldgenpath = '';
    if isstruct(varargin{1})
        datasets = {varargin{1}.datasetinfo.filename};
        try,
            paths = {varargin{1}.datasetinfo.filepath};
        catch,
            paths = cell(1,length(datasets));
            paths(:) = { '' };
        end
        genpath = varargin{1}.filepath;
        if isfield(varargin{1}.etc, 'oldfilepath')
            oldgenpath = varargin{1}.etc.oldfilepath;
        end
    else
        paths = varargin{1};
        if nargin > 1
            datasets = varargin{2};
        else
            datasets = paths;
            paths    = cell(size(datasets));
        end
    end

    set = 1;
    ALLEEG = [];

    eeglab_options;

    % read datasets
    % -------------
    comp = computer;
    warnfold = 'off';
    if length(oldgenpath) > 1 && oldgenpath(2) == ':' && ~strcmpi(comp(1:2), 'PC')
        oldgenpath = [ filesep oldgenpath(4:end) ];
        oldgenpath(find(oldgenpath == '\')) = filesep;
    end

    for dset = 1:length(paths)
        if ~isempty(paths{dset})
            if paths{dset}(2) == ':' && ~strcmpi(comp(1:2), 'PC')
                paths{dset} = [ filesep paths{dset}(4:end) ];
                paths{dset}(find(paths{dset} == '\')) = filesep;
            end
        end

        [sub2 sub1] = fileparts(char(paths{dset}));
        [sub3 sub2] = fileparts(sub2);

        % priority is given to relative path of the STUDY if STUDY has moved
        if ~isempty(oldgenpath) && oldgenpath(end) == filesep
             oldgenpath(end) = [];
        end
        if ~isequal(genpath, oldgenpath) && ~isempty(oldgenpath) %
            oldgenpath = strrep(oldgenpath, [ filesep filesep ], filesep);
            if ~isequal(genpath, oldgenpath) && ~isempty(oldgenpath)

                disp('Warning: STUDY moved since last saved, trying to load data files using relative path');
                if  ~isempty(strfind(char(paths{dset}), oldgenpath))
                    relativePath = char(paths{dset}(length(oldgenpath)+1:end));
                    relativePath = fullfile(genpath, relativePath);
                else
                    disp('Important warning: relative path cannot calculated, make sure the correct data files are loaded');
                    relativePath = char(paths{dset});
                end

                % fix issue when datasets are in a parent folder of the STUDY
                if dset == 1
                    indCommon = 1;
                    while indCommon <= length(oldgenpath) && indCommon <= length(paths{1}) && paths{1}(indCommon) == oldgenpath(indCommon)
                        indCommon = indCommon+1;
                    end
                    indCommon = indCommon-1;
                    if indCommon > 1 && indCommon < length(genpath) % do not change path if nothing in common between the two paths
                        genpath(indCommon-length(oldgenpath)+length(genpath)+1:end) = [];
                    end
                end
            else
                relativePath = char(paths{dset});
            end
        else
            relativePath = char(paths{dset});
        end

        % load files with a variety of extension using file-io
        [~,~,ext] = fileparts(datasets{dset});
        if ~isequal(lower(ext), '.set')
            disp('Importing binary file and resaving with a .set file extension as an EEGLAB dataset');
            EEG = pop_biosig(fullfile(relativePath, datasets{dset}));
            [pathTmp,fileTmp,ext] = fileparts(fullfile(relativePath, datasets{dset}));
            EEG = pop_saveset(EEG, fullfile(pathTmp, [ fileTmp '.set' ]));
        elseif exist(fullfile(relativePath, datasets{dset})) == 2
            EEG = pop_loadset('filename', datasets{dset}, 'filepath', relativePath, 'loadmode', 'info', 'check', 'off');
        elseif exist(fullfile(char(paths{dset}), datasets{dset})) == 2
            EEG = pop_loadset('filename', datasets{dset}, 'filepath', char(paths{dset}), 'loadmode', 'info', 'check', 'off');
        elseif exist( fullfile(genpath, datasets{dset})) == 2
            [tmpp tmpf ext] = fileparts(fullfile(genpath, datasets{dset}));
            EEG = pop_loadset('filename', [tmpf ext], 'filepath',tmpp, 'loadmode', 'info', 'check', 'off');
            warnfold = 'on';
        elseif exist( fullfile(genpath, sub1, datasets{dset})) == 2
            [tmpp tmpf ext] = fileparts(fullfile(genpath, sub1, datasets{dset}));
            EEG = pop_loadset('filename', [tmpf ext], 'filepath',tmpp, 'loadmode', 'info', 'check', 'off');
            warnfold = 'on';
        elseif exist( fullfile(genpath, sub2, datasets{dset}))  == 2
            [tmpp tmpf ext] = fileparts(fullfile(genpath, sub2, datasets{dset}));
            EEG = pop_loadset('filename', [tmpf ext], 'filepath',tmpp, 'loadmode', 'info', 'check', 'off');
            warnfold = 'on';
        elseif exist( fullfile(genpath, sub2, sub1, datasets{dset}))  == 2
            [tmpp tmpf ext] = fileparts(fullfile(genpath, sub2, sub1, datasets{dset}));
            EEG = pop_loadset('filename', [tmpf ext], 'filepath',tmpp, 'loadmode', 'info', 'check', 'off');
            warnfold = 'on';
        elseif exist(lower(fullfile(char(paths{dset}), datasets{dset}))) == 2
            EEG = pop_loadset('filename', lower(datasets{dset}), 'filepath',lower(char(paths{dset})), 'loadmode', 'info', 'check', 'off');
        elseif exist(fullfile(pwd, datasets{dset})) == 2
            EEG = pop_loadset('filename', lower(datasets{dset}), 'filepath',pwd, 'loadmode', 'info', 'check', 'off');
        else
            txt = [ sprintf('The dataset %s is missing\n', datasets{dset}) 10 ...
                          'Is it possible that it might have been deleted?' 10 ...
                          'If this is the case, re-create the STUDY using the remaining datasets' ];
            error(txt);
        end

        EEG = eeg_checkset(EEG);
        if ~option_storedisk
            EEG = eeg_checkset(EEG, 'loaddata');
        elseif ~ischar(EEG.data)
            EEG.data   = 'in set file';
            EEG.icaact = [];
        end

        [ALLEEG, EEG] = eeg_store(ALLEEG, EEG, 0, 'notext');
    end
    ALLEEG = eeg_checkset(ALLEEG);

    if strcmpi(warnfold, 'on') && ~strcmpi(pwd, genpath) && ~isempty(genpath)
        disp('Changing current path to STUDY path...');
        cd(genpath);
    end
    if strcmpi(warnfold, 'on')
        disp('This STUDY has a relative path set for the datasets');
        disp('so the current path MUST remain the STUDY path');
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
clear, clc;

% This is an example for running the function sparseInverseCovariance
%
%   max log( det( Theta ) ) - < S, Theta> - lambda * ||Theta||_1
%
% For detailed description of the function, please refer to the Manual.
%
%% Related papers
%
% The implementation is based on the following paper:
%
% [1]  Jerome Friedman, Trevor Hastie, and Robert Tibshirani,
%      Sparse inverse covariance estimation with the graphical lasso, 2007
%
% This function has been used in the following paper:
%
% [2] Shuai Huang, Jing Li, Liang Sun, Jun Liu, Teresa Wu,
%     Kewei Chen, Adam Fleisher, Eric Reiman, and Jieping Ye,
%     Learning Brain Connectivity of Alzheimer's Disease
%     from Neuroimaging Data, NIPS, 2009
%
%% ------------   History --------------------
%
% First version on September 18, 2008.
%
% For any problem, please contact Jun Liu (j.liu@asu.edu)

cd ..
cd ..

root=cd;
addpath(genpath([root '/SLEP']));
                     % add the functions in the folder SLEP to the path

% change to the original folder
cd Examples/invCov;

n=100;
% problem size

% generate the data
randn('state',1);
A=randn(n,n);

mean_1=mean(A,1);
A=A-repmat(mean_1,n,1);
norm_2=sqrt( sum(A.^2,1) );
A=A./repmat(norm_2,n,1);

% S is the empirical covariance matrix
S=A'*A;

% set the maximal number of iterations
opts.maxIter=100;

lambda=0.2;

tic;
Theta=sparseInverseCovariance(S, lambda, opts);
toc;

author: zzstefan

repo: BrainNetClass


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

% This function allows you to compute optimal
% center and blur values on a small subset of
% the MIT ICCV dataset (see below) by a very simple
% grid search over these values. You can specify
% according to which criteria you want to choose
% these values. This file prints and saves optimal values
% per metric, and returns the optimal values on
% the ROC metric, but you can choose your own
% strategy to weight the metrics and select a
% set of parameters.
% Given the parameters chosen on the ICCV
% dataset, you can apply them to your saliency
% benchmark maps.

function [optCenter,optBlur] = optimizeParams(mainpath,salMapDir,resDir,targetHist,showOutput)
% mainpath is the full path to ParameterTrainingOnICCV

% PATHS
FIXATIONDB = fullfile(mainpath,'ParameterTrainingOnICCV/ICCVdata/FixationDBCleanedICCV.mat'); % fixation data for dataset
FIXATIONMAPS = fullfile(mainpath,'ParameterTrainingOnICCV/ICCVdata/ALLFIXATIONMAPSICCV/');
STIMULI = fullfile(mainpath,'ParameterTrainingOnICCV/100_ICCV_STIMULI');
HASH = fullfile(mainpath,'ParameterTrainingOnICCV/ICCVdata/hash.mat');

if ~exist(resDir,'dir'), mkdir(resDir); end

if nargin < 5,
    showOutput = 0;
end

load('center.mat');
load(FIXATIONDB);

%% PARAMETERS

whichIms = 1:60; % which images to optimize on, from the STIMULI folder
                            % to save time and computation, can only use subset

metrics = {'S','ROC','ROC_borji','sROC_borji','CC','NSS','EMD'};
metricNames = {'Similarity','AUC (Judd)',...
                          'AUC (Borji)','shuffled AUC','Cross-correlation',...
                          'Normalized Scanpath Saliency','Earth Mover Distance'};
metricOpt = [1,1,1,1,1,1,-1];
% optimal value per metric is either the maximum (1) or the minimum (-1)

sigmaVals = [0,1,10,20,30,40,60,80,100]; % blur parameters to try
centVals = [0:0.1:1]; % center weight parameters to try

nims = length(whichIms);
nsigmas = length(sigmaVals);
ncents = length(centVals);

%% Run code on multiple cores for speed-up
try
    matlabpool open
end

%% Check that files are present and of right format

if isempty(salMapDir), error('%s is empty. Terminating.\n',salMapDir); end
files = dir(fullfile(salMapDir, '*.jpg'));
if isempty(files)
    error('Could not find image files...');
end

%% Create folders of varying blur and center
% NOTE: this code creates AND STORES saliency
% maps for all levels of blur and center, in order
% to be able to visually inspect intermediate
% results created. If you don't want to save all
% these extra images, rewrite this code to
% combine this cell with the next one, and calculate
% performance numbers without without saving
% maps.

% initialize empty directories
fprintf('Creating folders of varying blur and center.\n')
for s = 1:nsigmas
    for w=1:ncents
        curdir = fullfile(resDir,sprintf('center%d_blur%d',w,s));
        mkdir(curdir);
    end
end

% compute a blurred and centered version of each image
fprintf('Populating folders with images of varying blur and center.\n')

if showOutput, figure; end
parfor i = whichIms %can optimize all files or only a subset

    fprintf('On image %d',i);

    mapOrig = im2double(imread(fullfile(salMapDir, files(i).name)));
    cent = imresize(center, size(mapOrig));

    % try different blur values
    for s = 1:nsigmas

        % try different center weights
        for w=1:ncents

            fprintf('.') % put down one dot for each computation

            % directory for image files
            curdir = fullfile(resDir,sprintf('center%d_blur%d',w,s));

            map = processMap(mapOrig,cent,targetHist,sigmaVals(s),centVals(w));

            % save result
            imwrite(map, fullfile(curdir, files(i).name));

        end

    end

    fprintf('\n')

end

%% Test and record performance of each folder (different blur, center values)

% create a file to store all calculated performances
paramsSaveFile = fullfile(resDir, 'params.mat');
if exist(paramsSaveFile,'file')
    load(paramsSaveFile);
end

folders = dir(fullfile(resDir, 'center*'));

fprintf('Evaluating performance of each folder.\n')

for ii = length(metrics); %1:length(metrics)

    metric = metrics{ii};
    params.(metric).S= nan(ncents,nsigmas);
    params.(metric).scores = nan(ncents,nsigmas,nims);

    % test original performance (no center, no blur)
    scores = scoreModel_parallel(salMapDir, metric, 'jpg', FixationDBCleaned, FIXATIONMAPS, STIMULI, HASH);
    params.(metric).orig = scores(whichIms);

    % test all variants of center and blur
    for i = 1:length(folders)

        fprintf('On metric %s, folder %s...',metric,folders(i).name); tic
        scores = scoreModel_parallel(fullfile(resDir, folders(i).name), metric, 'jpg', FixationDBCleaned, FIXATIONMAPS, STIMULI, HASH);

        temp = regexp(folders(i).name,'center(\d*)_blur(\d*)','tokens');
        temp = temp{1};
        center = str2num(temp{1}); blur = str2num(temp{2});

        params.(metric).scores(center,blur,:) = scores;
        params.(metric).S(center,blur) = mean(scores);

        save(paramsSaveFile, 'params'); % save results after each iteration, just in case
        fprintf('Done in %2.2f s\n',toc)
    end

end

%% Replace this with your own code to choose optimal blur, center weight parameters

% this code chooses the best parameters for ROC,
% and prints out the best performance PER METRIC;
% you may like to combine all this information to
% choose a single optimal set of parameters.

fprintf('Selecting optimal parameters.\n')

for ii =1:length(metrics)
    metric = metrics{ii}; metricName = metricNames{ii};

    % get the blur and center weight parameters with highest performance on
    clear S

    S = params.(metric).S;
    if metricOpt(ii)==1
        [bestVal,idx] = max(S(:));
    else
        [bestVal,idx] = min(S(:));
    end
    [center,blur] = ind2sub(size(S),idx);

    % obtain the performance with no extra blur or center
    origPerf = params.(metric).orig;

    if (metricOpt(ii)==1 && mean(origPerf) >= bestVal) || ...
       (metricOpt(ii)==-1 && mean(origPerf) <= bestVal)
        fprintf('%s performance was best for original model (no extra optimizations).\n\n',metricName);
        optCenterParam = 0;
        optBlurParam = 0;
        scores = origPerf;
    else
        optCenterParam = centVals(center);
        optBlurParam = sigmaVals(blur);
        scores = params.(metric).scores(center,blur,:);
        fprintf('Optimal center weight = %2.2f, blur radius = %2.2f\n',optCenterParam,optBlurParam)
        fprintf('%s performance improves from %2.2f to %2.2f\n\n',metricName,mean(origPerf),bestVal)
    end

    params.(metric).optCenterParam = optCenterParam;
    params.(metric).optBlurParam = optBlurParam;
    params.(metric).optScore = scores;

    % save the chosen parameter settings
    save(paramsSaveFile, 'params');

end

% replace the following with your own choices
optCenter = params.ROC.optCenterParam;
optBlur = params.ROC.optBlurParam;

%% clean-up
try
    matlabpool close
end

fprintf('Done everything.\n')

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
function imorse = ic_to_imorse ( ic )

%*****************************************************************************80
%
%% IC_TO_IMORSE converts an ASCII integer code to a Morse integer code.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer IC, the integer code for the ASCII character.
%
%    Output, integer IMORSE, the integer code for the Morse character,
%    or -1 if no corresponding Morse code is available.
%
  junk = [ ...
    -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, ...
    -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, ...
     1, -1, 45, -1, -1, -1, -1, 42, -1, -1, -1, -1, 39, 43, 38, 44, ...
    37, 28, 29, 30, 31, 32, 33, 34, 35, 36, 40, -1, -1, -1, -1, 41, ...
    -1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, ...
    17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, -1, -1, -1, -1, -1, ...
    -1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, ...
    17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, -1, -1, -1, -1, -1, ...
    -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, ...
    -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, ...
    -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, ...
    -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, ...
    -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, ...
    -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, ...
    -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, ...
    -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 ];

  if ( 0 <= ic && ic <= 255 )
    imorse = junk(ic+1);
  else
    imorse = -1;
  end

  return
end


text:
function ds = setinterpreter(ds,funchandle)
   %SETINTERPRETER changes the associated interpreter function
   % datasource = setinterpreter(datasource,function_handle)

   if ~isa(funchandle,'function_handle')
      error('Trying to set a function handle to a %s',class(funchandle));
   end
   ds.interpreter = funchandle;
   % ds = datasource(funchandle,ds.file_string,ds.file_args{:});
end

author: geoscience-community-codes

repo: GISMO


text:
function [map,num,typ,scheme] = brewermap_view(N,scheme) %#ok<*ISMAT>
% An interactive figure for ColorBrewer colorscheme selection (RGB colormaps)
%
% (c) 2014-2022 Stephen Cobeldick
%
% View Cynthia Brewer's ColorBrewer colorschemes in a figure.
%
% * Two colorbars give the colorscheme in color and grayscale.
% * A button toggles between 3D-cube and 2D-lineplot of the RGB values.
% * A button toggles an endless demo cycle through the colorschemes.
% * A button reverses the colormap.
% * 35 buttons select any ColorBrewer colorscheme.
% * Text with the colorscheme's type (Diverging/Qualitative/Sequential)
% * Text with the colorscheme's number of nodes (defining colors).
%
%%% Syntax:
% brewermap_view
% brewermap_view(N)
% brewermap_view(N,scheme)
% brewermap_view([],...)
% brewermap_view(axes/figure handles,...) % see "Adjust Colormaps"
% [map,num,typ] = brewermap_view(...)
%
% Calling the function with an output argument blocks MATLAB execution until
% the figure is deleted: the final colormap and colorscheme are then returned.
%
%% Adjust Colormaps of Figures or Axes %%
%
% Only R2014b or later. Provide axes or figure handles as the first input
% and their colormaps will be updated in real-time by BREWERMAP_VIEW.
%
%%% Example:
%
% >> S = load('spine');
% >> image(S.X)
% >> brewermap_view(gca)
%
%% Input and Output Arguments %%
%
%%% Inputs (*=default):
% N = NumericScalar, an integer to define the colormap length.
%   = *[], colormap length of two hundred and fifty-six (256).
%   = NaN, same length as the defining RGB nodes (useful for Line ColorOrder).
%   = Array of axes/figure handles. R2014b or later only.
% scheme = CharRowVector or StringScalar, a ColorBrewer colorscheme name.
%
%%% Outputs (these block code execution until the figure is closed!):
% map = NumericMatrix, the colormap defined when the figure is closed.
% num = NumericVector, the number of nodes defining the ColorBrewer colorscheme.
% typ = CharRowVector, the colorscheme type: 'Diverging'/'Qualitative'/'Sequential'.
%
% See also BREWERMAP BREWERMAP_PLOT CUBEHELIX PRESET_COLORMAP MAXDISTCOLOR
% RGBPLOT COLORMAP COLORMAPEDITOR COLORBAR UICONTROL ADDLISTENER

%% Input Wrangling %%
%
persistent ax2D ln2D ax3D pt3D txtH is2D cbAx cbIm pTxt pSld bEig bGrp bRev scm isr
%
new = isempty(ax2D)||~ishghandle(ax2D);
dfn = 256;
upd = false;
upb = false;
hgv = [];
nmr = dfn;
%
err = 'First input must be a real positive scalar numeric or [] or NaN.';
if nargin==0 || isnumeric(N)&&isequal(N,[])
	N = dfn;
elseif isnumeric(N)
	assert(isscalar(N),'SC:brewermap_view:NotScalarNumeric',err)
	assert(isnan(N)||isreal(N)&&isfinite(N)&&fix(N)==N&&N>0,...
		'SC:brewermap_view:NotRealPositiveNotNaN',err)
	N = double(N);
elseif all(ishghandle(N(:))) % R2014b or later
	assert(isgraphics(N(:),'axes')|isgraphics(N(:),'figure'),...
		'SC:brewermap_view:NotAxesNorFigureHandles',...
		'First input may be an array of figure or axes handles.')
	hgv = N(:);
	nmr = arrayfun(@(h)size(colormap(h),1),hgv);
	N = nmr(1);
else
	error('SC:brewermap_view:UnsupportedInput',err)
end
%
[mcs,nmn,pyt] = brewermap('list');
%
% Check BREWERMAP output:
tmp = find([any(diff(double(char(pyt)),1),2);1]);
assert(isequal(tmp,[9;17;35]),'SC:brewermap_view:SchemeSequence',...
	'The BREWERMAP function returned an unexpected scheme sequence.')
%
% Default pseudo-random colorscheme:
if nargin==0 || new
	isr = rand(1)>0.5;
	scm = mcs{1+mod(round(now*1e7),numel(mcs))};
end
%
% Parse input colorscheme:
if nargin==2
	scheme = bmv1s2c(scheme);
	assert(ischar(scheme) && ndims(scheme)==2 && size(scheme,1)==1,...
		'SC:brewermap_view:NotCharacterVector',...
		'Second input <scheme> must be a 1xN char vector.')
	% Check if a reversed colormap was requested:
	isr = strncmp(scheme,'-',1) || strncmp(scheme,'*',1);
	isd = strncmp(scheme,'+',1) || isr;
	scm = scheme(1+isd:end);
end
%
if isnan(N)
	N = nmn(strcmpi(scm,mcs));
end
%
%% Ensure Figure Exists %%
%
% LHS and RHS slider bounds/limits, and slider step sizes:
lbd = 1;
rbd = dfn;
stp = [1,10]; % [minor,major]
%
% Define the 3D cube axis order:
xyz = 'RGB'; % choose order
[~,xyz] = ismember(xyz,'RGB');
%
if new % Create a new figure.
	%
	% Figure parameters:
	M = 9; % buttons per column
	gap = 0.01; % gaps
	bth = 0.04; % demo button height
	btw = 0.10; % demo button width
	bgh = 0.40; % button group height
	cbw = 0.23; % colorbar width (both together)
	cbh = 1-3*gap-bth; % colorbar height
	axh = 1-bgh-2*gap; % axes height
	axw = 1-cbw-2*gap; % axes width
	bgw = axw-btw-gap; % button group width
	%
	figH = figure('HandleVisibility','callback', 'Color','white',...
		'IntegerHandle','off', 'NumberTitle','off', 'Units','normalized',...
		'Name','ColorBrewer Interactive ColorScheme Selector',...
		'MenuBar','figure', 'Toolbar','none', 'Tag',mfilename);
	%
	% Add 2D lineplot:
	ax2D = axes('Parent',figH, 'Position',[gap,bgh+gap,axw,axh], 'Box','on',...
		'ColorOrder',[1,0,0; 0,1,0; 0,0,1; 0.6,0.6,0.6], 'HitTest','off',...
		'Visible','off', 'XLim',[0,1], 'YLim',[0,1], 'XTick',[], 'YTick',[]);
	ln2D = line([0,0,0,0;1,1,1,1],[0,0,0,0;1,1,1,1], 'Parent',ax2D,...
		'Visible','off', 'Linestyle','-', 'Marker','.');
	%
	% Add 3D scatterplot:
	ax3D = axes('Parent',figH, 'OuterPosition',[0,bgh,axw+2*gap,1-bgh],...
		'Visible','on', 'XLim',[0,1], 'YLim',[0,1], 'ZLim',[0,1], 'HitTest','on');
	pt3D = patch('Parent',ax3D, 'XData',[0;1], 'YData',[0;1], 'ZData',[0;1],...
		'Visible','on', 'LineStyle','none', 'FaceColor','none', 'MarkerEdgeColor','none',...
		'Marker','o', 'MarkerFaceColor','flat', 'MarkerSize',10, 'FaceVertexCData',[1,1,0;1,0,1]);
	view(ax3D,3);
	grid(ax3D,'on')
	lbl = {'Red','Green','Blue'};
	xlabel(ax3D,lbl{xyz(1)})
	ylabel(ax3D,lbl{xyz(2)})
	zlabel(ax3D,lbl{xyz(3)})
	%
	% Add warning text:
	txtH = text('Parent',ax2D, 'Units','normalized', 'Position',[1,1],...
		'HorizontalAlignment','right', 'VerticalAlignment','top', 'Color','k');
	%
	% Add demo button:
	demo = uicontrol(figH, 'Style','togglebutton', 'Units','normalized',...
		'Position',[1-cbw/2,1-bth-gap,cbw/2-gap,bth], 'String','Demo',...
		'Max',1, 'Min',0, 'Callback',@bmvDemo); %#ok<NASGU>
	% Add 2D/3D button:
	is2D = uicontrol(figH, 'Style','togglebutton', 'Units','normalized',...
		'Position',[1-cbw/1,1-bth-gap,cbw/2-gap,bth], 'String','2D / 3D',...
		'Max',1, 'Min',0, 'Callback',@bmv2D3D);
	% Add reverse button:
	bRev = uicontrol(figH, 'Style','togglebutton', 'Units','normalized',...
		'Position',[bgw+2*gap,bgh-bth,btw,bth], 'String','Reverse',...
		'Max',1, 'Min',0, 'Callback',@bmvRevM);
	%
	% Add colorbars:
	C(1,1,:) = [1,1,1];
	cbAx(2) = axes('Parent',figH, 'Visible','off', 'Units','normalized',...
		'Position',[1-cbw/2,gap,cbw/2-gap,cbh], 'YLim',[0.5,1.5],...
		'YDir','reverse', 'HitTest','off');
	cbAx(1) = axes('Parent',figH, 'Visible','off', 'Units','normalized',...
		'Position',[1-cbw/1,gap,cbw/2-gap,cbh], 'YLim',[0.5,1.5],...
		'YDir','reverse', 'HitTest','off');
	cbIm(2) = image('Parent',cbAx(2), 'CData',C);
	cbIm(1) = image('Parent',cbAx(1), 'CData',C);
	%
	% Add parameter slider, listener, and corresponding text:
	sv = mean([lbd,rbd],2);
	pTxt = uicontrol(figH,'Style','text', 'Units','normalized',...
		'Position',[bgw+2*gap,bgh-2*bth-gap,btw,bth], 'String','X');
	pSld = uicontrol(figH,'Style','slider', 'Units','normalized',...
		'Position',[bgw+2*gap,gap,btw,bgh-2*bth-gap], 'Min',lbd(1), 'Max',rbd(1),...
		'SliderStep',stp(1,:)/(rbd(1)-lbd(1)), 'Value',sv(1));
	addlistener(pSld, 'Value', 'PostSet',@bmvSldr);
	%
	% Add colorscheme button group:
	bGrp = uibuttongroup('Parent',figH, 'BorderType','none', 'Units','normalized',...
		'BackgroundColor','white', 'Position',[gap,gap,bgw,bgh-gap]);
	% Determine button locations:
	Z = 1:numel(mcs);
	Z = Z+(Z>17);
	C = (ceil(Z/M)-1)/4;
	R = (M-1-mod(Z-1,M))/M;
	% Add colorscheme buttons to group:
	for jj = numel(mcs):-1:1
		bEig(jj) = uicontrol('Parent',bGrp, 'Style','Toggle', 'String',mcs{jj},...
			'Unit','normalized', 'Position',[C(jj),R(jj),1/4,1/M]);
	end
	set(bGrp,'SelectionChangeFcn',@bmvChgS);
	%
end
%
%% Nested Functions %%
%
	function str = bmvName()
		% Generate the colorscheme name.
		drn = '+-';
		ids = strcmpi(scm,mcs);
		if any(ids)
			str = [drn(1+isr),mcs{ids}];
		else % will throw an error in BREWERMAP:
			str = scm;
		end
	end
%
	function bmvUpDt()
		% Update all graphics objects in the figure.
		%
		% Get ColorBrewer colormap and grayscale equivalent:
		[map,num,typ] = brewermap(N,bmvName());
		mag = map*[0.298936;0.587043;0.114021];
		%
		% Update colorbar values:
		set(cbAx, 'YLim', [0,abs(N)+(N==0)]+0.5);
		set(cbIm(1), 'CData',reshape(map,[],1,3))
		set(cbIm(2), 'CData',repmat(mag,[1,1,3]))
		%
		% Update 2D line / 3D patch values:
		if  get(is2D, 'Value') % 2D
			set(ln2D, 'XData',linspace(0,1,abs(N)));
			set(ln2D,{'YData'},num2cell([map,mag],1).');
		else % 3D
			set(pt3D,...
				'XData',map(:,xyz(1)),...
				'YData',map(:,xyz(2)),...
				'ZData',map(:,xyz(3)),...
				'FaceVertexCData',map)
		end
		%
		% Update reverse button:
		set(bRev, 'Value',isr)
		%
		% Update warning text:
		str = {[bmvName(),' '];[typ,' '];sprintf('%d Nodes ',num)};
		set(txtH,'String',str);
		%
		% Update parameter value text:
		set(pTxt(1), 'String',sprintf('N = %.0f',N));
		%
		% Update external axes/figure:
		nmr(1) = N;
		for k = 1:numel(hgv)
			colormap(hgv(k),brewermap(nmr(k),bmvName()));
		end
		%
		drawnow()
	end
%
	function bmv2D3D(h,~)
		% Switch between 2D-line and 3D-cube representation.
		%
		if get(h,'Value') % 2D
			set(ax3D, 'HitTest','off', 'Visible','off')
			set(ax2D, 'HitTest','on',  'Visible','on')
			set(pt3D, 'Visible','off')
			set(ln2D, 'Visible','on')
		else % 3D
			set(ax2D, 'HitTest','off', 'Visible','off')
			set(ax3D, 'HitTest','on',  'Visible','on')
			set(ln2D, 'Visible','off')
			set(pt3D, 'Visible','on')
		end
		%
		bmvUpDt();
	end
%
	function bmvChgS(~,e)
		% Change the colorscheme.
		%
		scm = get(e.NewValue,'String');
		%
		bmvUpDt()
	end
%
	function bmvRevM(h,~)
		% Reverse the colormap.
		%
		isr = logical(get(h,'Value'));
		%
		bmvUpDt()
	end
%
	function bmvSldr(~,~)
		% Update the slider position.
		%
		if ~upd
			return
		end
		%
		N = round(get(pSld,'Value'));
		%
		bmvUpDt()
	end
%
	function bmvDemo(h,~)
		% Display all ColorBrewer colorschemes sequentially.
		%
		cnt = 0;
		while ishghandle(h)&&get(h,'Value')
			cnt = mod(cnt+1,pow2(53));
			%
			if mod(cnt,23)<1
				ids = 1+mod(find(strcmpi(scm,mcs)),numel(mcs));
				scm = mcs{ids};
				try %#ok<TRYNC>
					set(bGrp, 'SelectedObject',bEig(ids));
				end
			end
			%
			if mod(cnt,69)<1
				isr = ~isr;
			end
			%
			upb = (upb || N<=1) && N<dfn;
			N = N - 1 + 2*upb;
			%
			% Update slider position:
			upd = false;
			try %#ok<TRYNC>
				set(pSld, 'Value',N)
				bmvUpDt()
			end
			upd = true;
			%
			% Faster/slower:
			pause(0.1);
		end
		%
	end
%
%% Initialize the Figure %%
%
set(bGrp,'SelectedObject',bEig(strcmpi(scm,mcs)));
set(pSld,'Value',max(lbd,min(rbd,N)));
upd = true;
bmvUpDt()
%
if nargout
	waitfor(ax2D);
	scheme = bmvName();
else
	clear map
end
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%brewermap_view
function arr = bmv1s2c(arr)
% If scalar string then extract the character vector, otherwise data is unchanged.
if isa(arr,'string') && isscalar(arr)
	arr = arr{1};
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%bmv1s2c
%
% Copyright (c) 2014-2022 Stephen Cobeldick
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%license

author: DrosteEffect

repo: BrewerMap


text:
classdef Zmap3DGridFunction < ZmapGridFunction
    % ZMAP3DGRIDFUNCTION is a ZmapFunction that produces a grid of volumetric results
    %
    % see also ZMAPGRIDFUNCTION

    properties
        features={'borders'}; % features to show on the map, such as 'borders','lakes','coast',etc.
    end
    properties(Constant)
        Type = GridTypes.XYZ;
    end

    methods

        function obj=Zmap3DGridFunction(varargin)
            unimplemented_error()
            obj@ZmapGridFunction(varargin{:});
        end

        function plot(obj,choice, varargin)
            % plots the results on the provided axes.
            if ~exist('choice','var')
                choice=obj.active_col;
            end
            if ~isnumeric(choice)
                choice = find(strcmp(obj.Result.values.Properties.VariableNames,choice));
            end

            mydesc = obj.Result.values.Properties.VariableDescriptions{choice};
            myname = obj.Result.values.Properties.VariableNames{choice};

            f=findobj(groot,'Tag',obj.PlotTag,'-and','Type','figure');
            if isempty(f)
                f=figure('Tag',obj.PlotTag);
            end
            figure(f);
            set(f,'name',['results from bvalgrid : ', myname])
            delete(findobj(f,'Type','axes'));

            % this is to show the data
            obj.Grid.pcolor([],obj.Result.values.(myname), mydesc);
            set(gca,'NextPlot','add');

            % the imagesc exists is to enable data cursor browsing.
            obj.plot_image_for_cursor_browsing(myname, mydesc, choice);

            shading(obj.ZG.shading_style);
            set(gca,'NextPlot','add')

            obj.add_grid_centers();

            ax=gca;
            for n=1:numel(obj.features)
                ft=obj.ZG.features(obj.features{n});
                copyobj(ft,ax);
            end

            colorbar
            title(mydesc)
            xlabel('Longitude')
            ylabel('Latitude')

            dcm_obj=datacursormode(gcf);
            dcm_obj.Updatefcn=@ZmapGridFunction.mydatacursor;
            if isempty(findobj(gcf,'Tag','lookmenu'))
                add_menu_divider();
                lookmenu=uimenu(gcf,'label','graphics','Tag','lookmenu');
                shademenu=uimenu(lookmenu,'Label','shading','Tag','shading');

                % TODO: combine mapdata_viewer with this function
                exploremenu=uimenu(gcf,'label','explore');
                uimenu(exploremenu,'label','explore'    , 'MenuSelectedFcn',@(~,~)mapdata_viewer(obj.Result,obj.RawCatalog,gcf));

                uimenu(shademenu,'Label','interpolated' , 'MenuSelectedFcn',@(~,~)shading('interp'));
                uimenu(shademenu,'Label','flat'         , 'MenuSelectedFcn',@(~,~)shading('flat'));

                plottype=uimenu(lookmenu,'Label','plot type');
                uimenu(plottype,'Label','Pcolor plot','Tag','plot_pcolor',...
                    'MenuSelectedFcn',@(src,~)obj.plot(choice),'Checked','on');

                % countour-related menu items

                uimenu(plottype,'Label','Plot Contours','Tag','plot_contour',...
                    'Enable','off',...not fully unimplmented
                    'MenuSelectedFcn',@(src,~)obj.contour(choice));
                uimenu(plottype,'Label','Plot filled Contours','Tag','plot_contourf',...
                    'Enable','off',...not fully unimplmented
                    'MenuSelectedFcn',@(src,~)contourf(choice));
                uimenu(lookmenu,'Label','change contour interval',...
                    'Enable','off',...
                    'MenuSelectedFcn',@(src,~)changecontours());

                % display overlay menu items

                uimenu(lookmenu,'Label','Show grid centerpoints','Checked',char(obj.showgridcenters),...
                    'MenuSelectedFcn',@obj.togglegrid_cb);
                uimenu(lookmenu,'Label',['Show ', obj.RawCatalog.Name, ' events'],...
                    'MenuSelectedFcn',@(src,~)obj.addquakes_cb(src,obj.RawCatalog));

                uimenu(lookmenu,'Separator','on',...
                    'Label','brighten',...
                    'MenuSelectedFcn',@(~,~)colormap(ax,brighten(colormap,0.4)));
                uimenu(lookmenu,'Label','darken',...
                    'MenuSelectedFcn',@(~,~)colormap(ax,brighten(colormap,-0.4)));

            end

            update_layermenu(obj,myname);
        end % plot function

    end % Public methods

    methods(Access=protected)
        function plot_image_for_cursor_browsing(obj, myname, mydesc, choice)
            h=obj.Grid.imagesc([],obj.Result.values.(myname));
            h.AlphaData=zeros(size(h.AlphaData))+0.0;

            % add some details that can be picked up by the interactive data cursor
            h.UserData.vals= obj.Result.values;
            h.UserData.choice=choice;
            h.UserData.myname=myname;
            h.UserData.myunit=obj.Result.values.Properties.VariableUnits{choice};
            h.UserData.mydesc=obj.Result.values.Properties.VariableDescriptions{choice};
        end

        function update_layermenu(obj, myname)
            if isempty(findobj(gcf,'Tag','layermenu'))
                layermenu=uimenu(gcf,'Label','layer','Tag','layermenu');
                for i=1:width(obj.Result.values)
                    tmpdesc=obj.Result.values.Properties.VariableDescriptions{i};
                    tmpname=obj.Result.values.Properties.VariableNames{i};
                    uimenu(layermenu,'Label',tmpdesc,'Tag',tmpname,...
                        'Enable',tf2onoff(~all(isnan(obj.Result.values.(tmpname)))),...
                        'MenuSelectedFcn',@(~,~)plot_cb(tmpname));
                end
            end

            % make sure the correct option is checked
            layermenu=findobj(gcf,'Tag','layermenu');
            set(findobj(layermenu,'Tag',myname),'Checked','on');

            % plot here
            function plot_cb(name)
                set(findobj(layermenu,'type','uimenu'),'Checked','off');
                obj.plot(name);
            end
        end

    end % Protected methods
    methods(Access=protected, Static)
        function addquakes_cb(src, catalog)
            qtag=findobj(gcf,'tag','quakes');
            if isempty(qtag)
                set(gca,'NextPlot','add')
                plot(catalog.X, catalog.Y, 'o',...
                    'MarkerSize',3,...
                    'markeredgecolor',[.2 .2 .2],...
                    'tag','quakes');
                set(gca,'NextPlot','replace')
            else
                ison=qtag.Visible == "on";
                qtag.Visible=tf2onoff(~ison);
                src.Checked=tf2onoff(~ison);
                drawnow
            end
        end
    end
end

author: CelsoReyes


text:
function ye = fL(y, D, hom)
% Neumann condition at outflow
z =  2^(-2.5)*(y*y*0.25/D - 1)* exp(-y*y/(8*D));
z = z + 2^(-2.5)*((2-y)*(2-y)*0.25/D - 1)*exp(-(2-y)*(2-y)/(8*D));
ye = hom*z;

repo: Aero-matlab


text:
helpdlg('SEISMOLAP is no longer supported as a part of ZMAP. Please contact  Prof. J. Zschau , GFZ Potsdam, zschau@gzf-potsdam.de, for information on recent SeismoLap implementations')

report_this_filefun(mfilename('fullpath'));

return

ttlStr='  Introduction to Seismo Lap in ZMAP';
hlpStr= ...
    ['                                                                             '
    ' SeismoLap is a technique develloped by J. Zschau                            '
    ' (Geoforschungzentrum Potsdam, Germany) to measure the relative              '
    ' seismic activity in an area. The space-time overlap of cubical              '
    ' volumes is measured, normalized by the size of the cube (SeismoLap1)        '
    ' The negaite invers of this value Seismolap2 is calculated and the           '
    ' probabilty estimated using a Pearson type 3 distribution.                   '
    '                                                                             '
    ' The ZMAP implementation of SeismoLap offers the following options:          '
    '                                                                             '
    ' One point calculation: SeismoLap1 and the probability can be calculated     '
    '     at one point in space. It is nessesary to define the Input paramters    '
    '     lat/long of the point of interest, box siz, and interaction time.       '
    '  Grid: Both the seismolap 1 and probability values can be calculated        '
    '       using a grid of points. At each grid-point the time series            '
    '       is calcualted and the user can view different time cuts in a color    '
    '       representation. The time series at a user defined point can be viewed '
    '       by selecting the Mouse option. Since it take a while to               '
    '       calculate a grid, it is recommended to save the grid and re-load them '
    '       when desired.                                                         '
    ' Movie: A number of time cuts can be animated as a movie sequence.           '
    '        This movie file can be save and re-loaded.                           '
    '                                                                             '
    ' If   you have problems whith this software, please contact:                 '
    ' Stefan Wiemer                                                               '
    ' Geophysical Institute University of Alaska Fairbanks                        '
    ' Fairbanks, AK, 99775-7320, USA                                              '
    ' phone: 907 474 6171, Fax 907 474 7290                                       '
    ' e-mail stefan@giseis.alaska.edu                                             '
    '                                                                             '
    ' Enjoy!                                                                      '
    ' Stefan Wiemer   September . 95                                              '
    '                                                                             '];

zmaphelp(ttlStr,hlpStr)


text:
function [gb_thin_CSG, gb_thin_CS, gb_fat_CS, texton] = simpleGPb(I)

% place the current directory on the top of paths (so rgb2lab in this
% directory will get called when needed)
curr_dir = fileparts(which(mfilename));
rmpath(curr_dir);
addpath(curr_dir);

[gb_thin_CSG, gb_thin_CS, gb_fat_CS] = Gb_CSG(I);
texton = genTexton(I);
gb_thin_CSG = uint8(round(255*gb_thin_CSG));
gb_thin_CS = uint8(round(255*gb_thin_CS));
gb_fat_CS = 255*gb_fat_CS;
texton = uint8(texton);

end

author: Cloud-CV

repo: object-proposals


text:
function x = p09_start ( option, nvar )

%*****************************************************************************80
%
%% P09_START returns a starting point for problem 9.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Output, real X(NVAR), the starting point.
%
  x = zeros ( nvar, 1 );

  [ height, ival, val ] = p09_gx ( option );

  x(ival) = val;

  return
end


text:
function stopfun = stopifdeletedfile(filename)
% Create an interactive stopping criterion based on the existence of a file
%
% function stopfun = stopifdeletedfile()
% function stopfun = stopifdeletedfile(filename)
%
% Use on the options structure passed to a Manopt solver, e.g.:
%
%   problem = ... % manopt problem structure with manifold, cost function, ...
%   options.stopfun = stopifdeletedfile(); % add this option
%   trustregions(problem, x0, options); % run this or any other solver
%
% This will create a temporary file called MANOPT_DELETE_ME_TO_STOP_SOLVER
% in the present working directory. If this file is deleted at any time
% during the solver's execution, the solver will terminate gracefully and
% return its current iterate as soon as it gets to the point of evaluating
% stopping criteria. A different file name can also be specified using the
% input string filename (optional).
%
% Note: certain solvers (including trustregions) check stopping criteria
% only at outer iterations, not during inner iterations; hence, their may
% be a delay before actual termination.
%
% See also: statsfunhelper stopifclosedfigure

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Aug. 2, 2018.
% Contributors:
% Change log:

    % Default name for the temporary file.
    if ~exist('filename', 'var') || isempty(filename)
        filename = 'MANOPT_DELETE_ME_TO_STOP_SOLVER';
    end

    % Make sure the file exists, and release our handle on it.
    fid = fopen(filename, 'a');
    if fid >= 0
        fclose(fid);
        % The stopping criterion is a function handle.
        stopfun = @checkcriterion;
    else
        warning('manopt:stopifdeletedfile', ...
              'Couldn''t create the file: no stopping criterion created.');
        stopfun = @(problem, x, info, last) false;
    end

    % The function is defined as a subfunction so that it has access to
    % filename without the need for an @() construct. This makes it easier
    % for Matlab to determine the number of output arguments of the
    % function handle @checkcriterion, which ultimately helps
    % stoppingcriterion determine how to call it.
    function [stop, reason] = checkcriterion(problem, x, info, last) %#ok<INUSD>

        reason = sprintf(['Interactive stopping criterion ' ...
                     '(file %s deleted). See options.stopfun.'], filename);

        % Try to access the file.
        fid = fopen(filename, 'r');

        % If we can't, it means the file no longer exists: stop the solver.
        % Otherwise, release our handle on the file to make sure it can be
        % deleted by another program.
        if fid < 0
            stop = true;
        else
            fclose(fid);
            stop = false;
        end

    end

end

author: NicolasBoumal

repo: manopt


text:
% Implementation of the MEM-EKF* algorithm based on the article
%
% "Tracking the Orientation and Axes Lengths of an Elliptical Extended Object"
% Shishan Yang and Marcus Baum
% arXiv preprint, 2018,
% https://arxiv.org/abs/1805.03276
%
% Further information:
% http://www.fusion.informatik.uni-goettingen.de
% https://github.com/Fusion-Goettingen
%
% Source code written by Shishan Yang
% =============================

function [r,p,Cr,Cp]= time_update(r,p,Cr,Cp,Ar, Ap,Cwr, Cwp)
r = Ar*r;
Cr = Ar*Cr*Ar'+Cwr;

p = Ap*p;
Cp = Ap*Cp*Ap'+Cwp;
end

author: Fusion-Goettingen

repo: ExtendedObjectTracking


text:
function [data, clamped] = mk_mutilated_samples(bnet, ncases, max_clamp, usecell)
% GEN_MUTILATED_SAMPLES Do random interventions and then draw random samples
% [data, clamped] = gen_mutilated_samples(bnet, ncases, max_clamp, usecell)
%
% At each step, we pick a random subset of size 0 .. max_clamp, and
% clamp these nodes to random values.
%
% data(i,m) is the value of node i in case m.
% clamped(i,m) = 1 if node i in case m was set by intervention.

if nargin < 4, usecell = 1; end

ns = bnet.node_sizes;
n = length(bnet.dag);
if usecell
  data = cell(n, ncases);
else
  data = zeros(n, ncases);
end
clamped = zeros(n, ncases);

csubsets = subsets(1:n, max_clamp, 0); % includes the empty set
distrib_cset = normalise(ones(1, length(csubsets)));

for m=1:ncases
  cset = csubsets{sample_discrete(distrib_cset)};
  nvals = prod(ns(cset));
  distrib_cvals = normalise(ones(1, nvals));
  cvals = ind2subv(ns(cset), sample_discrete(distrib_cvals));
  mutilated_bnet = do_intervention(bnet, cset, cvals);
  ev = sample_bnet(mutilated_bnet);
  if usecell
    data(:,m) = ev;
  else
    data(:,m) = cell2num(ev);
  end
  clamped(cset,m) = 1;
end

author: bayesnet


text:
function [out] = evap_21(p1,p2,S,Ep,dt)
%evap_21

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Threshold-based evaporation with constant minimum rate
% Constraints:  f <= S/dt
% @(Inputs):    p1   - wilting point (1st threshold) [mm]
%               p2   - 2nd threshold as fraction of wilting point [-]
%               S    - current storage [mm]
%               Ep   - potential evapotranspiration rate [mm/d]
%               dt   - time step size [d]

out = min(max(p2,min(S./p1,1)).*Ep,S./dt);

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
executable = '"../../x64/Release/FeatureExtraction.exe"';

outputDir = 'D:\Datasets\face2face\2006_processed/';

% First collect the filenames of the data to be used
input_label_dir = 'D:\Datasets\face2face\rapport-oct-2006-all-transcriptions/';

folds = dir([input_label_dir, '*SES*']);

listener_file_labels = {};
speaker_file_labels = {};

listener_vid_files = {};
speaker_vid_files = {};

speaker_vid_dir = 'D:\Datasets\face2face\rapport-oct-2006-all-movie-speaker/';
listener_vid_dir = 'D:\Datasets\face2face\rapport-oct-2006-all-movie-listener/';
verbose = true;
for i=1:numel(folds)

    listener_file = dir([input_label_dir, folds(i).name, '/*.L.nod.eaf']);
    speaker_file = dir([input_label_dir, folds(i).name, '/*.S.nod.eaf']);

    if(~isempty(listener_file))
        % Need to find the appropriate video file if it exists
        num = listener_file.name(end-13:end-10);
        vid_file_dir = dir([listener_vid_dir, '/*', num, '*']);
        vid_file = dir([listener_vid_dir, '/', vid_file_dir.name, '/*.mp4']);
        if(~isempty(vid_file))
            listener_vid_files = cat(1, listener_vid_files, [listener_vid_dir, '/', vid_file_dir.name, '/', vid_file.name]);
            listener_file_labels = cat(1, listener_file_labels, [input_label_dir, '/' folds(i).name, '/' listener_file.name]);
        end
    end

    if(~isempty(speaker_file))
        num = speaker_file.name(end-13:end-10);
        vid_file_dir = dir([speaker_vid_dir, '/*', num, '*']);
        vid_file = dir([speaker_vid_dir, '/', vid_file_dir.name, '/*.mp4']);
        if(~isempty(vid_file))
            speaker_vid_files = cat(1, speaker_vid_files, [speaker_vid_dir, '/', vid_file_dir.name, '/', vid_file.name]);
            speaker_file_labels = cat(1, speaker_file_labels, [input_label_dir, '/' folds(i).name, '/' speaker_file.name]);
        end
    end

end

% file_labels = cat(1, listener_file_labels, speaker_file_labels);
% video_files = cat(1, listener_vid_files, speaker_vid_files);

file_labels = listener_file_labels;
video_files = listener_vid_files;

parfor i=1:numel(file_labels)

    [~,short_name,vid_ext] = fileparts(video_files{i});

    command = executable;

    inputFile = video_files{i};
    outputFile = [outputDir short_name '.txt'];

    outputEaf = [outputDir short_name '.eaf'];

    command = cat(2, command, [' -f "' inputFile '" -of "' outputFile  '"']);

    if(verbose)
        outputVideo = [outputDir short_name '.track.avi'];
        command = cat(2, command, [' -ov "' outputVideo '"']);
    end

    command = cat(2, command, [' -no2Dfp -no3Dfp -noMparams -noAUs -noGaze']);

    dos(command);

    copyfile(file_labels{i}, outputEaf);
    copyfile(video_files{i}, [outputDir '/' short_name, '.mp4']);

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function sparse_test ( )

%*****************************************************************************80
%
%% SPARSE_TEST tests the SPARSE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test MATLAB''s sparse matrix features.\n' );

  sparse_test01 ( );
  sparse_test02 ( );
  sparse_test03 ( );
  sparse_test04 ( );
  sparse_test05 ( );
  sparse_test06 ( );
  sparse_test07 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function boolVal=pointIsInTriangle(point2D,triVertices)
%%POINTISINTRIANGLE Determine whether a specified 2D points is within a
%           triangle.
%
%INPUTS: point2D A 2XnumPts set of 2D points.
%    triVertices A 2X3 set of vertices of the triangle in any order.
%
%OUTPUTS: bookVal A 1XnumPts set of boolean values that are true if the
%                 point is in the triangle or is on an edge of the
%                 triangle.
%
%The function pt2TriangCoords2D is used to convert the points to 2D
%triangular barycentric coordinates. If the point is within the triangle,
%then all of the cooridnates will be positive. Otherwise, some coordinates
%will be negative.
%
%EXAMPLE:
%In this example, the first 3 points are in or on the triangle and the
%final two points are outside the trangle. Thus, boolPoints=[1,1,1,0,0];
% v=[0,1,0.5;
%    0,0,0.5];
% points=[0.25,0,1,10,2;
%         0.25,0,0,-1,0];
% boolVals=pointIsInTriangle(points,v)
%
%September 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    numPts=size(point2D,2);
    phi=pt2TriangCoords2D(point2D,triVertices);

    boolVal=false(1,numPts);
    for k=1:numPts
        boolVal(k)=all(phi(:,k)>=0);
    end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [LL, prior, transmat, mu, Sigma, mixmat] = ...
     mhmm_em(data, prior, transmat, mu, Sigma, mixmat, varargin);
% LEARN_MHMM Compute the ML parameters of an HMM with (mixtures of) Gaussians output using EM.
% [ll_trace, prior, transmat, mu, sigma, mixmat] = learn_mhmm(data, ...
%   prior0, transmat0, mu0, sigma0, mixmat0, ...)
%
% Notation: Q(t) = hidden state, Y(t) = observation, M(t) = mixture variable
%
% INPUTS:
% data{ex}(:,t) or data(:,t,ex) if all sequences have the same length
% prior(i) = Pr(Q(1) = i),
% transmat(i,j) = Pr(Q(t+1)=j | Q(t)=i)
% mu(:,j,k) = E[Y(t) | Q(t)=j, M(t)=k ]
% Sigma(:,:,j,k) = Cov[Y(t) | Q(t)=j, M(t)=k]
% mixmat(j,k) = Pr(M(t)=k | Q(t)=j) : set to [] or ones(Q,1) if only one mixture component
%
% Optional parameters may be passed as 'param_name', param_value pairs.
% Parameter names are shown below; default values in [] - if none, argument is mandatory.
%
% 'max_iter' - max number of EM iterations [10]
% 'thresh' - convergence threshold [1e-4]
% 'verbose' - if 1, print out loglik at every iteration [1]
% 'cov_type' - 'full', 'diag' or 'spherical' ['full']
%
% To clamp some of the parameters, so learning does not change them:
% 'adj_prior' - if 0, do not change prior [1]
% 'adj_trans' - if 0, do not change transmat [1]
% 'adj_mix' - if 0, do not change mixmat [1]
% 'adj_mu' - if 0, do not change mu [1]
% 'adj_Sigma' - if 0, do not change Sigma [1]
%
% If the number of mixture components differs depending on Q, just set  the trailing
% entries of mixmat to 0, e.g., 2 components if Q=1, 3 components if Q=2,
% then set mixmat(1,3)=0. In this case, B2(1,3,:)=1.0.

if ~isstr(varargin{1}) % catch old syntax
  error('optional arguments should be passed as string/value pairs')
end

[max_iter, thresh, verbose, cov_type,  adj_prior, adj_trans, adj_mix, adj_mu, adj_Sigma] = ...
    process_options(varargin, 'max_iter', 10, 'thresh', 1e-4, 'verbose', 1, ...
		    'cov_type', 'full', 'adj_prior', 1, 'adj_trans', 1, 'adj_mix', 1, ...
		    'adj_mu', 1, 'adj_Sigma', 1);

previous_loglik = -inf;
loglik = 0;
converged = 0;
num_iter = 1;
LL = [];

if ~iscell(data)
  data = num2cell(data, [1 2]); % each elt of the 3rd dim gets its own cell
end
numex = length(data);

O = size(data{1},1);
Q = length(prior);
if isempty(mixmat)
  mixmat = ones(Q,1);
end
M = size(mixmat,2);
if M == 1
  adj_mix = 0;
end

while (num_iter <= max_iter) & ~converged
  % E step
  [loglik, exp_num_trans, exp_num_visits1, postmix, m, ip, op] = ...
      ess_mhmm(prior, transmat, mixmat, mu, Sigma, data);

  % M step
  if adj_prior
    prior = normalise(exp_num_visits1);
  end
  if adj_trans
    transmat = mk_stochastic(exp_num_trans);
  end
  if adj_mix
    mixmat = mk_stochastic(postmix);
  end
  if adj_mu | adj_Sigma
    [mu2, Sigma2] = mixgauss_Mstep(postmix, m, op, ip, 'cov_type', cov_type);
    if adj_mu
      mu = reshape(mu2, [O Q M]);
    end
    if adj_Sigma
      Sigma = reshape(Sigma2, [O O Q M]);
    end
  end

  if verbose, fprintf(1, 'iteration %d, loglik = %f\n', num_iter, loglik); end
  num_iter =  num_iter + 1;
  converged = em_converged(loglik, previous_loglik, thresh);
  previous_loglik = loglik;
  LL = [LL loglik];
end

%%%%%%%%%

function [loglik, exp_num_trans, exp_num_visits1, postmix, m, ip, op] = ...
    ess_mhmm(prior, transmat, mixmat, mu, Sigma, data)
% ESS_MHMM Compute the Expected Sufficient Statistics for a MOG Hidden Markov Model.
%
% Outputs:
% exp_num_trans(i,j)   = sum_l sum_{t=2}^T Pr(Q(t-1) = i, Q(t) = j| Obs(l))
% exp_num_visits1(i)   = sum_l Pr(Q(1)=i | Obs(l))
%
% Let w(i,k,t,l) = P(Q(t)=i, M(t)=k | Obs(l))
% where Obs(l) = Obs(:,:,l) = O_1 .. O_T for sequence l
% Then
% postmix(i,k) = sum_l sum_t w(i,k,t,l) (posterior mixing weights/ responsibilities)
% m(:,i,k)   = sum_l sum_t w(i,k,t,l) * Obs(:,t,l)
% ip(i,k) = sum_l sum_t w(i,k,t,l) * Obs(:,t,l)' * Obs(:,t,l)
% op(:,:,i,k) = sum_l sum_t w(i,k,t,l) * Obs(:,t,l) * Obs(:,t,l)'

verbose = 0;

%[O T numex] = size(data);
numex = length(data);
O = size(data{1},1);
Q = length(prior);
M = size(mixmat,2);
exp_num_trans = zeros(Q,Q);
exp_num_visits1 = zeros(Q,1);
postmix = zeros(Q,M);
m = zeros(O,Q,M);
op = zeros(O,O,Q,M);
ip = zeros(Q,M);

mix = (M>1);

loglik = 0;
if verbose, fprintf(1, 'forwards-backwards example # '); end
for ex=1:numex
  if verbose, fprintf(1, '%d ', ex); end
  %obs = data(:,:,ex);
  obs = data{ex};
  T = size(obs,2);
  if mix
    [B, B2] = mixgauss_prob(obs, mu, Sigma, mixmat);
    [alpha, beta, gamma,  current_loglik, xi, gamma2] = ...
	fwdback(prior, transmat, B, 'obslik2', B2, 'mixmat', mixmat);
  else
    B = mixgauss_prob(obs, mu, Sigma);
    [alpha, beta, gamma,  current_loglik, xi] = fwdback(prior, transmat, B);
  end
  loglik = loglik +  current_loglik;
  if verbose, fprintf(1, 'll at ex %d = %f\n', ex, loglik); end

  exp_num_trans = exp_num_trans + sum(xi,3);
  exp_num_visits1 = exp_num_visits1 + gamma(:,1);

  if mix
    postmix = postmix + sum(gamma2,3);
  else
    postmix = postmix + sum(gamma,2);
    gamma2 = reshape(gamma, [Q 1 T]); % gamma2(i,m,t) = gamma(i,t)
  end
  for i=1:Q
    for k=1:M
      w = reshape(gamma2(i,k,:), [1 T]); % w(t) = w(i,k,t,l)
      wobs = obs .* repmat(w, [O 1]); % wobs(:,t) = w(t) * obs(:,t)
      m(:,i,k) = m(:,i,k) + sum(wobs, 2); % m(:) = sum_t w(t) obs(:,t)
      op(:,:,i,k) = op(:,:,i,k) + wobs * obs'; % op(:,:) = sum_t w(t) * obs(:,t) * obs(:,t)'
      ip(i,k) = ip(i,k) + sum(sum(wobs .* obs, 2)); % ip = sum_t w(t) * obs(:,t)' * obs(:,t)
    end
  end
end
if verbose, fprintf(1, '\n'); end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function triangle_lyness_rule_test02 ( )

%*****************************************************************************80
%
%% TRIANGLE_LYNESS_RULE_TEST02 performs the weight sum test on Lyness rules.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGLE_LYNESS_RULE_TEST02\n' );
  fprintf ( 1, '  LYNESS_RULE returns the points and weights\n' );
  fprintf ( 1, '  of a Lyness rule for the triangle.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this test, we simply check that the weights\n' );
  fprintf ( 1, '  sum to 1.\n' );

  rule_num = lyness_rule_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of available rules = %d\n', rule_num );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      Rule    Sum of weights\n' );
  fprintf ( 1, '\n' );

  for rule = 0 : rule_num

    order = lyness_order ( rule );

    [ w, x ] = lyness_rule ( rule, order );

    w_sum = sum ( w(1:order) );

    fprintf ( 1, '  %8d  %25.16f\n', rule, w_sum );

  end

  return
end


text:
function imageData = DMdicomread(fileName_or_metaHeader)
    % DMdicomread Reads image data from a Dicom file
    %
    % Usage:
    %     imageData = DMdicomread(fileName)
    %         fileName: path and filename of the file to test
    %
    %     imageData = DMdicomread(metaheader)
    %         metaheader: structure returned from DMdicominfo
    %
    %  Returns:
    %     The pixel data
    %
    %
    %
    %     Licence
    %     -------
    %     Part of DicoMat. https://github.com/tomdoel/dicomat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the BSD 3-Clause license. Please see the file LICENSE for details.
    %

    if isstruct(fileName_or_metaHeader)
        fileName = fileName_or_metaHeader.Filename;
    else
        fileName = fileName_or_metaHeader;
    end

    header = DMReadDicomTags(fileName, DMDicomDictionary.EssentialDictionary);
    if isfield(header, 'PixelData')
        imageData = header.PixelData;
    else
        imageData = [];
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [ c, seed ] = c8vec_uniform_01 ( n, seed )

%*****************************************************************************80
%
%% C8VEC_UNIFORM_01 returns a unit pseudorandom C8VEC.
%
%  Discussion:
%
%    The angles should be uniformly distributed between 0 and 2 * PI,
%    the square roots of the radius uniformly distributed between 0 and 1.
%
%    This results in a uniform distribution of values in the unit circle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Second Edition,
%    Springer, 1987,
%    ISBN: 0387964673,
%    LC: QA76.9.C65.B73.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, December 1986, pages 362-376.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley, 1998,
%    ISBN: 0471134031,
%    LC: T57.62.H37.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, Number 2, 1969, pages 136-143.
%
%  Parameters:
%
%    Input, integer N, the number of values to compute.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, complex C(N), the pseudorandom complex vector.
%
%    Output, integer SEED, a seed for the random number generator.
%
  c = zeros ( n, 1 );

  i4_huge = 2147483647;

  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'C8VEC_UNIFORM_01 - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'C8VEC_UNIFORM_01 - Fatal error!' );
  end

  for j = 1 : n

    k = floor ( seed / 127773 );

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

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

    r = sqrt ( seed * 4.656612875E-10 );

    k = floor ( seed / 127773 );

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

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

    theta = 2.0 * pi * seed * 4.656612875E-10;

    c(j) = r * ( cos ( theta ) + sin ( theta ) * i );

  end

  return
end


text:
function [weight] = ons_kernel(data, weight_o, eta, beta, delta)
% This program output the final portfolio the ONS strategy on data
% Only one expert, thus, go through to the exper routine
% If multiple experts, combine them similar to BK algorithm (bk_kernel.m)
%
% function [weight] = ons_kernel(data, weight_o, eta, beta, delta)
%
% weight: final portfolio, used for next rebalance
%
% data: market sequence vectors
% weight_o: last portfolio, also can be last price relative adjusted.
% eta, beta, delta: ONS' parameters
%
% Example: [weight] = ons_kernel(data, weight_o, 0, 1, 1/8)

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

weight = ons_expert(data, weight_o, eta, beta, delta);

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text: %HORZCAT Datafile overload: horizontal concatenation

author: marianux

repo: ecg-kit

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


text:
classdef MimBrushSizeSlider < MimGuiPluginSlider
    % MimBrushSizeSlider. Gui Plugin for changing the size of the edit brush
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Brush size'
        SelectedText = 'Brush size'

        ToolTip = 'Change the size of the editing paint brush'
        Category = 'Paint'
        Visibility = 'Dataset'
        Mode = 'Edit'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
        Location = 42

        MinValue = 0
        MaxValue = 100
        SmallStep = 0.01
        LargeStep = 0.1
        DefaultValue = 50

        EditBoxPosition = 90
        EditBoxWidth = 30

        StackVertically = false
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded && gui_app.ImagePanel.OverlayImage.ImageExists && ...
                isequal(gui_app.ImagePanel.SelectedControl, 'Paint');
        end

        function is_selected = IsSelected(gui_app)
            is_selected = true;
        end

        function [value_instance_handle, value_property_name, limits_instance_handle, limits_property_name] = GetHandleAndProperty(gui_app)
            value_instance_handle = gui_app.ImagePanel;
            value_property_name = 'PaintBrushSize';
            limits_instance_handle = [];
            limits_property_name = [];
        end

    end
end


text:
%MDL_PLANAR2 Create model of a simple planar 2-link mechanism
%
% MDL_PLANAR2 is a script that creates the workspace variable p2 which
% describes the kinematic characteristics of a simple planar 2-link
% mechanism.
%
% Also defines the vector:
%   qz   corresponds to the zero joint angle configuration.
%
% Also defines the vector:
%   qz   corresponds to the zero joint angle configuration.
%
% Notes::
% - Moves in the XY plane.
% - No dynamics in this model.
%
% See also mdl_twolink, mdl_planar1, mdl_planar3, SerialLink.

% MODEL: generic, planar, 2DOF, symbolic, standard_DH

% Copyright (C) 1993-2017, 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

syms a1 a2 real;

p2 = SerialLink([
    Revolute('d', 0, 'a', a1, 'alpha', 0, 'standard')
    Revolute('d', 0, 'a', a2, 'alpha', 0, 'standard')
    ], ...
    'name', 'two link');
qz = [0 0];

author: petercorke

repo: robotics-toolbox-matlab


text:
function [output,idx,validFrameMask] = F_tmaxpool(input_layer, curr_layer)
input = input_layer.a;

% the input contains N samples, each of DxT size
[D,T,N] = size(input);
precision = class(gather(input(1)));

if isfield(curr_layer, 'context')
    context = curr_layer.context;
else
    context = 0;
end
if isfield(curr_layer, 'stride')
    stride = curr_layer.stride;
else
    stride = 0;
end

if context==0 || stride==0  % global pooling
    [output,idx] = max(input,[], 2);
    validFrameMask = zeros(1,N);
else
    nWindow = length(1:stride:(T-context+1));
    if strcmpi(class(input), 'gpuArray')
        output = gpuArray.zeros(D, length(nWindow), N, precision);
    else
        output = zeros(D, length(nWindow), N, precision);
    end
    idx = output;
    for i=1:nWindow
        offset = (i-1)*stride;
        [output(:,i,:), idx(:,i,:)] = max(input(:,offset+1:offset+context,:), [], 2);
    end
    validFrameMask = zeros(nWindow,N);
end

end

author: singaxiong

repo: SignalGraph


text:
%% fnc_SelectInput: select the input indexes from a generic list
%
% Usage:
%   index = fnc_SelectInput(pro, iset)
%
% Inputs:
%    pro                project structure
%    iset               cell array or array of inputs of the considered set, they can be selected
%                       by index (1,2,3 ...) or by name ('in1','x',..) or
%                       mixed
%
% Output:
%     index              vector of indexes corresponding to the iset inputs
%
% ------------------------------------------------------------------------
%
% Author : Flavio Cannavo'
% e-mail: flavio(dot)cannavo(at)gmail(dot)com
% Release: 1.0
% Date   : 02-02-2011
%
% History:
% 1.0  02-02-2011  First release.
%%

function index = fnc_SelectInput(pro, iset)

index = nan(size(iset));

if iscell(iset)
    for i=1:length(iset)
        in = iset{i};
        if isnumeric(in)
            index(i) = floor(in);
        else
            for j=1:length(pro.Inputs.pdfs)
                if strcmp(in,pro.Inputs.Names{j})
                    index(i) = j;
                end
            end
        end
    end
else
    index = iset;
end

if sum(isnan(index))>0
    disp('Warning: some input is not well defined');
end

index = sort(unique(index));

author: Sable

repo: mcbench-benchmarks


text:
% P(X) = GAUSSIAN(MU0,COV0)
%
% Q(X) = GAUSSIAN(MU,COV)
%
% where the parameters are actually functions
%   MU = MU(THETA),
%   COV = COV(THETA).
%
% More generally, the whole joint distribution Q(X,MU0,COV0) can be a
% function of THETA.
%
% Evaluates KL(Q||P) w.r.t. THETA up to a constant, and the gradient.
%
% vb_rotationcost_gaussian(logdet_Cov, x_invCov0_x, x_invCov0_mu, ...
%                          mu_invCov0_mu, logdet_Cov0, grad_logdet_Cov, ...
%                          grad_x_invCov0_x, grad_x_invCov0_mu, ...
%                          grad_mu_invCov0_mu, grad_logdet_Cov0)
%
% LOGDET_COV         : <LOG(DET(COV))>
% X_INVCOV0_X        : <X'*INV(COV0)*X>
% X_INVCOV0_MU       : <X'*INV(COV0)*MU>
% MU_INVCOV0_MU      : <MU'*INV(COV0)*MU>
% LOGDET_COV0        : <LOG(DET(COV0))>
% GRAD_LOGDET_COV    : D(<LOG(DET(COV))>) / D(THETA)
% GRAD_X_INVCOV0_X   : D(<X'*INV(COV0)*X>) / D(THETA)
% GRAD_X_INVCOV0_MU  : D(<X'*INV(COV0)*MU>) / D(THETA)
% GRAD_MU_INVCOV0_MU : D(<MU'*INV(COV0)*MU>) / D(THETA)
% GRAD_LOGDET_COV0   : D(<LOG(DET(COV0))>) / D(THETA)
%
% where the expectations <...> are taken over the distribution
% Q(X,MU0,COV0). Note that the terms can be evaluated up to a constant
% w.r.t. THETA, which often gives huge savings.

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

function [KL,dKL] = vb_rotationcost_gaussian(logdet_Cov,         ...
                                             x_invCov0_x,        ...
                                             x_invCov0_mu,       ...
                                             mu_invCov0_mu,      ...
                                             logdet_Cov0,        ...
                                             grad_logdet_Cov,    ...
                                             grad_x_invCov0_x,   ...
                                             grad_x_invCov0_mu,  ...
                                             grad_mu_invCov0_mu, ...
                                             grad_logdet_Cov0)

% Cost from <log q(X)>
log_qX = -0.5 * logdet_Cov;
dlog_qX = -0.5 * grad_logdet_Cov;

% Cost from <log p(x)>
log_pX = -0.5 * (logdet_Cov0 + x_invCov0_x - 2*x_invCov0_mu + mu_invCov0_mu);
dlog_pX = -0.5 * (grad_logdet_Cov0 + grad_x_invCov0_x - 2*grad_x_invCov0_mu ...
                  + grad_mu_invCov0_mu);

KL = log_qX - log_pX;
dKL = dlog_qX - dlog_pX;

% $$$ function [f,df_dR] = vb_rotationcost_gaussian(R,logdet_R,invT_R, ...
% $$$                                               XX_R_invCov,X_Mu_invCov)
% $$$
% $$$ [D,N] = size(X);
% $$$
% $$$ % TODO: Check that the sizes match
% $$$
% $$$ % Cost from <log q(X)>
% $$$ log_qX = -N*logdet_R;
% $$$ dlog_qX = -N*invT_R;
% $$$
% $$$ % Cost from <log p(x)>
% $$$ log_pX = -0.5 * (traceprod(XX_R_invCov,R) - 2*traceprod(X_Mu_invCov,R));
% $$$ dlog_pX = X_Mu_invCov - XX_R_invCov;
% $$$
% $$$ f = log_pX -

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab