text:
function newTree = expandTree(tree, maxOrder)
%EXPANDTREE   Convert expressions like 5*(diff(u) + u) to 5*diff(u) + 5*u.
%   The role of the EXPANDTREE method is to split up syntax trees, so that
%   expressions involving the highest derivative appearing in the equation stand
%   alone, that is, not inside parenthesis. For example, EXPANDTREE will convert
%   expressions like 5*(diff(u) + u) to 5*diff(u) + 5*u. This is necessary in
%   order to be able to put the correct expression on the right-hand side when
%   calling MATLAB's built-in solvers.
%
%   Calling sequence:
%      NEWTREE = EXPANDTREE(TREE, MAXORDER)
%   where the inputs are:
%      TREE:       The syntax tree we want to split.
%      MAXORDER:   A vector containing the maximum differential order of each
%                  variable that appear in the problem under consideration.
%   and the output is:
%      NEWTREE:    A syntax tree where the highest order deriative has been
%                  taken out of any parenthesis where other variables with lower
%                  differential order appear.

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

if ( ~isstruct(tree) || tree.height <= 1 || all(tree.diffOrder < maxOrder) ||...
    ( ~tree.hasTerms && sum(tree.ID) <= 1) )
    % We don't need to worry about expanding the tree any of the following
    % conditions apply:
    %   * The input is not a tree
    %   * It is a very short tree
    %   * Its diffOrder is less than the maxOrder we consider
    %   * It doesn't any operators that split it up in terms (that is, + or -),
    %     and it doesn't include more than one variable (which would imply
    %     nonlinearity, such as u.*w).
    newTree = tree;

elseif ( tree.numArgs == 1 )
    % If we have a unary operator, we expand the center node recursively.
    newTree = expandTree(tree.center, maxOrder);

elseif ( any(strcmp(tree.method, {'plus', 'minus'})) )
    % We're at a + or a -, so expand the syntax tree recursively:
    tree.left = treeVar.expandTree(tree.left, maxOrder);
    tree.right = treeVar.expandTree(tree.right, maxOrder);
    newTree = tree;

else
    % Must be at .*, ./ or .^.

    % If both left and right arguments to the operator are syntax trees, and the
    % highest order derivative appear in any, we know that we must have
    % nonlinearities in the highest order derivative, e.g. in the case of a
    % first order IVP, u.*diff(u). This is not supported, so throw an error:
    if ( isstruct(tree.left) && isstruct(tree.right) && ...
            ( any(tree.diffOrder == maxOrder) ) )
        error('CHEBFUN:TREEVAR:expandTree:nonlinearity', ...
            ['Nonlinearity in highest order derivative detected.\n' ...
            'Unable to convert to first order format.']);
    end

    % Find out what kind of argument we're dealing with for the left argument
    % of the operator
    if ( ~isstruct(tree.left) )
        % Left argument is not a syntax tree, must be either a CHEBFUN or a
        % scalar. In either case, the left tree has zero children.
        leftNumArgs = 0;
    else
        % Left tree is indeed a syntax tree, store its number of children.
        leftNumArgs = tree.left.numArgs;
    end

    % Find out what kind of argument we're dealing with for the right argument
    % of the operator:
    if ( ~isstruct(tree.right) )
        % Right argument is not a syntax tree, must be either a CHEBFUN or a
        % scalar. In either case, the right tree has zero children.
        rightDiffOrder = 0;
        rightNumArgs = 0;
    else
        % Left tree is indeed a syntax tree, store its number of children.
        rightNumArgs = tree.right.numArgs;
    end

    % Split up the left tree depending on how many children it has. In case we
    % don't need to split at the current level, we store the differential
    % order and the height of the left tree in the variables LEFTDIFFORDER
    % and LEFTHEIGHT. Finally, we introduce the Boolean variables SPLITLEFT
    % which indicates whether the left tree had to be split:
    if ( leftNumArgs == 0 || ...
            (strcmp(tree.left.method,'diff') && tree.left.height <= 1 )  )
        % We're either dealing with a CHEBFUN/scalar left tree, or it's
        % differentiating one of the basis variables directly (as that's the
        % only case that tree.left.height <= 1). In this case, no splitting is
        % needed.
        leftTree = tree.left;
        leftDiffOrder = 0*tree.diffOrder;
        leftHeight = 0;
        splitLeft = false;
    elseif ( leftNumArgs == 1 )
        % Here we're dealing with the unary operator case. In this case, we
        % don't need to split at the current level, however, we recursively
        % split up the child tree.
        leftTree = treeVar.expandTree(tree.left, maxOrder);
        leftDiffOrder = leftTree.diffOrder;
        leftHeight = leftTree.height;
        splitLeft = false;
    else
        % Here we're dealing with the binary operator case. In this case, we do
        % need to split at the current level, as well as recursively splitting
        % up both the child trees. LEFTLEFT and LEFTRIGHT will be the syntax
        % trees stored in the left child of the current level, once we've split
        % up the child trees of the current level.
        leftLeft = treeVar.expandTree(tree.left.left, maxOrder);
        leftRight = treeVar.expandTree(tree.left.right, maxOrder);
        splitLeft = true;
    end

    % Split up the right tree depending on how many children it has. In case we
    % don't need to split at the current level, we store the differential order
    % and the height of the right tree in the variables RIGHTDIFFORDER and
    % RIGHTHEIGHT. Finally, we introduce the Boolean variables SPLITRIGHT which
    % indicates whether the right tree had to be split:
    if ( rightNumArgs == 0 || ...
            (strcmp(tree.right.method,'diff') && tree.right.height <= 1 )  )
        % We're either dealing with a CHEBFUN/scalar right tree, or it's
        % differentiating one of the basis variables directly (as that's the
        % only case that tree.right.height <= 1). In this case, no splitting is
        % needed.
        rightTree = tree.right;
        rightDiffOrder = 0*tree.diffOrder;
        rightHeight = 0;
        splitRight = false;
    elseif ( rightNumArgs == 1 )
        % He're we're dealing with the unary operator case. In this case, we
        % don't need to split at the current level, however, we recursively
        % split up the child tree.
        rightTree = treeVar.expandTree(tree.right, maxOrder);
        rightDiffOrder = rightTree.diffOrder;
        rightHeight = 0;
        splitRight = false;
    else
        % He're we're dealing with the binary operator case. In this case, we do
        % need to split at the current level, as well as recursively splitting
        % up both the child trees. RIGHTLEFT and RIGHTRIGHT will be the syntax
        % trees stored in the right child of the current level, once we've split
        % up the child trees of the current level.
        rightLeft  = treeVar.expandTree(tree.right.left, maxOrder);
        rightRight = treeVar.expandTree(tree.right.right, maxOrder);
        % Only need to worry about splitting if both left and right trees are
        % syntax trees, not if we have CHEBFUN/scalars
        splitRight = true;
    end

    % If we had to split either the left tree or the right tree, construct new
    % syntax trees which we'll then glue together below. Note that we should
    % never expect to split both the left and right trees, as that indicates
    % nonlinearities in the highest order derivative.
    if ( ~splitLeft && ~splitRight )
        % The simple case, no splitting required, so the new left and right
        % syntax trees will be the current ones:
        newLeft = leftTree;
        newRight = rightTree;

    elseif ( splitLeft && ~splitRight )
        % Had to split on the left, not the right.

        % Multiply together the new left factor of the left tree, and the
        % current right tree:
        newLeft = struct('method','times', 'numArgs', 2, ...
            'left', leftLeft, ...
            'right', rightTree,...
            'diffOrder', max(getDifforder(leftLeft), rightDiffOrder), ...
            'height', max(getHeight(leftLeft), rightHeight) + 1, ...
            'hasTerms', getHasTerms(leftLeft) || getHasTerms(rightTree) );

        % Multiply together the new right factor of the left tree, and the
        % current right tree:
        newRight = struct('method','times', 'numArgs', 2, ...
            'left', leftRight, ...
            'right', rightTree, ...
            'diffOrder', max(getDifforder(leftRight), rightDiffOrder) + 1, ...
            'height', max(getHeight(leftRight), rightHeight) + 1, ...
            'hasTerms', getHasTerms(leftRight) || getHasTerms(rightTree) );

    elseif ( ~splitLeft && splitRight )
        % Had to split on the right, not the left.

        % Multiply together the current left tree, and the new left factor of
        % the right tree:
        newLeft = struct('method','times', 'numArgs', 2, ...
            'left', leftTree, ...
            'right', rightLeft,...
            'diffOrder', max(leftDiffOrder, getDifforder(rightLeft)), ...
            'height', max(leftHeight, getHeight(rightLeft)) + 1, ...
            'hasTerms', getHasTerms(leftTree) || getHasTerms(rightLeft) );

        % Multiply together the current left tree, and the new right factor of
        % the right tree:
        newRight = struct('method','times', 'numArgs', 2, ...
            'left', leftTree, ...
            'right', rightRight,...
            'diffOrder', max(leftDiffOrder, getDifforder(rightRight)), ...
            'height', max(leftHeight, getHeight(rightRight) + 1), ...
            'hasTerms', getHasTerms(leftTree) || getHasTerms(rightRight) );

    end

    % If the left or right trees still have multiple terms left in them, we need
    % to split them up recursively:
    if ( newLeft.hasTerms )
        newLeft = treeVar.expandTree(newLeft, maxOrder);
    end
    if ( newRight.hasTerms )
        newRight = treeVar.expandTree(newRight, maxOrder);
    end

    % Add together the new left and right trees, splitting complete!
    newTree = struct('method', 'plus', 'numArgs', 2, ...
        'left', newLeft, 'right', newRight, ...
        'diffOrder', max(newLeft.diffOrder, newRight.diffOrder), ...
        'height',max(newLeft.height, newRight.height) + 1, ...
        'hasTerms', newLeft.hasTerms || newRight.hasTerms);
end
end

function out = getDifforder(treeIn)
%GETDIFFORDER   The diffOrder of a TREEVAR.
%   We often encounter child nodes that may or may not be a syntax tree, as
%   opposed to chebfuns or scalars. We are usually interested in the heights
%   and diffOrders of the syntax trees, however, if the tree is actually a
%   CHEBFUN/scalar, they won't have the necessary fields. This method takes
%   care of checking the class of the argument, and return the correct
%   information.
if ( isstruct(treeIn) )
    out = treeIn.diffOrder;
else
    out = 0;
end

end

function out = getHeight(treeIn)
%GETHEIGHT   The height of a TREEVAR.
%   Same as GETDIFFORDER() method but for the height of a syntax tree.
if ( isstruct(treeIn) )
    out = treeIn.height;
else
    out = 0;
end

end

function out = getHasTerms(treeIn)
%GETHASTERMS   Test if a TREEVAR contains a 'plus' or 'minus' node.
%   Same as GETDIFFORDER() method but for whether the syntax tree consists of
%   multiple terms.
if ( isstruct(treeIn) )
    out = treeIn.hasTerms;
else
    out = 0;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function r = gedrnd(v,varargin)
% Generate random variables from a Generalized Error Distribution (GED) with V degrees of freedom
%
% USAGE:
%   R=gedrnd(V)
%   R=gedrnd(V,S1)
%   R=gedrnd(V,S1,S2,S3,...,SN)
%   R=gedrnd(V,[S1 S2 S3 ... SN])
%
% INPUTS:
%   V     - Degree of freedom parameter.  Either scalar or matrix.
%   Sx    - [OPTIONAL] Size of output.
%           R will be:
%               1 by 1 if V is scalar and no S are entered
%               size(V) if V is non-scalar and no S are entered
%               S1 by S1 if V is scalar and only S1 is entered
%               [S1 S2 ... SN] otherwise
%               If V is non-scalar and S are provided, size(V) must equal [S1 S2 ... SN]
%
% OUTPUTS:
%   R     - GED distributed random variables
%
% COMMENTS:
%   V>=1
%   A scalar GED r.v. with variance normalized to 1 has probability
%   density given by:
%       f(x,v) = [v/(lda*2^(1+1/v)*gamma(1/v))]*exp(-0.5*|x/lda|^v)
%       lda = [2^(-2/v)*gamma(1/v)/gamma(3/v)]^0.5
%   If X~GED(v) then abs(X)^v = Y~Gamma(1/v) (cf. Tadikamalla 1980).
%   GAMRND does the computational work.
%
% REFERENCES:
%   [1] Tadikamalla (1980), J.Am.Stat.Assoc. (75)
%   [2] Nelson (1991), Econometrica
%
% See also GEDPDF, GEDCDF, GEDINV, GEDLL, GAMRND

% Copyright: Ivana Komunjer
% komunjer@hss.caltech.edu
% Modifications Copyright:
% Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 9/1/2005

if nargin < 1,
    error('Requires at least one input argument.');
end

v(v<1)=NaN;

[err, errtext, sizeOut] = iscompatible(1,v,varargin{:});

if err
    error(errtext)
end

%   Return NaN if the argument V is outside its limit.
r = gamrnd(1./v,1,sizeOut);
r = r.^(1./v);
rndsgn = 2*((rand(sizeOut)>0.5)-0.5);  % generates a random sign -1 or +1
r = r.*rndsgn;

scalex = (gamma(3./v)./gamma(1./v)).^0.5; % standardizes the obtained values
r = r./scalex;

author: bashtage

repo: mfe-toolbox


text:
function r = esvm_apply_M(x, boxes, M)
% Apply boosting "co-occurrence" matrix M to the boxes
% function r = esvm_apply_M(x, boxes, M)
%
% Apply the multiplexer matrix M which boosts the scores of a
% window based on its friends and their scores embedded in the
% context feature vector x
%
% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

if prod(size(x))==0
  r = zeros(1,0);
  return;
end
exids = boxes(:,6);

%Because the co-occurrence matrix treats exemplar flips as separate
%exemplars, we need to check if an exemplar has been flipped at
%update its exemplar index
exids(boxes(:,7)==1) = exids(boxes(:,7)==1) + size(x,1)/2;
r = zeros(1,size(boxes,1));

for i = 1:size(boxes,1)
  r(i) = (M.w{exids(i)}'*x(:,i) + sum(x(:,i)))-M.b{exids(i)};
end

author: quantombone

repo: exemplarsvm


text:
function [ferns,hsPr] = fernsClfTrain( data, hs, varargin )
% Train random fern classifier.
%
% See "Fast Keypoint Recognition in Ten Lines of Code" by Mustafa Ozuysal,
% Pascal Fua and Vincent Lepetit, CVPR07.
%
% Dimensions:
%  M - number ferns
%  S - fern depth
%  F - number features
%  N - number input vectors
%  H - number classes
%
% USAGE
%  [ferns,hsPr] = fernsClfTrain( data, hs, [varargin] )
%
% INPUTS
%  data     - [NxF] N length F feature vectors
%  hs       - [Nx1] target output labels in [1,H]
%  varargin - additional params (struct or name/value pairs)
%   .S        - [10] fern depth (ferns are exponential in S)
%   .M        - [50] number of ferns to train
%   .thrr     - [0 1] range for randomly generated thresholds
%   .bayes    - [1] if true combine probs using bayes assumption
%   .ferns    - [] if given reuse previous ferns (recompute pFern)
%
% OUTPUTS
%  ferns    - learned fern model w the following fields
%   .fids     - [MxS] feature ids for each fern for each depth
%   .thrs     - [MxS] threshold corresponding to each fid
%   .pFern    - [2^SxHxM] learned log probs at fern leaves
%   .bayes    - if true combine probs using bayes assumption
%   .inds     - [NxM] cached indices for original training data
%   .H        - number classes
%  hsPr     - [Nx1] predicted output labels
%
% EXAMPLE
%  N=5000; H=5; d=2; [xs0,hs0,xs1,hs1]=demoGenData(N,N,H,d,1,1);
%  fernPrm=struct('S',4,'M',50,'thrr',[-1 1],'bayes',1);
%  tic, [ferns,hsPr0]=fernsClfTrain(xs0,hs0,fernPrm); toc
%  tic, hsPr1 = fernsClfApply( xs1, ferns ); toc
%  e0=mean(hsPr0~=hs0); e1=mean(hsPr1~=hs1);
%  fprintf('errors trn=%f tst=%f\n',e0,e1); figure(1);
%  subplot(2,2,1); visualizeData(xs0,2,hs0);
%  subplot(2,2,2); visualizeData(xs0,2,hsPr0);
%  subplot(2,2,3); visualizeData(xs1,2,hs1);
%  subplot(2,2,4); visualizeData(xs1,2,hsPr1);
%
% See also fernsClfApply, fernsInds
%
% Piotr's Image&Video Toolbox      Version 2.50
% Copyright 2010 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Lesser GPL [see external/lgpl.txt]

% get additional parameters and check dimensions
dfs={'S',10,'M',50,'thrr',[0 1],'bayes',1,'ferns',[]};
[S,M,thrr,bayes,ferns]=getPrmDflt(varargin,dfs,1);
[N,F]=size(data); assert(length(hs)==N);
H=max(hs); assert(all(hs>0)); assert(S<=20);

if( isempty(ferns) )
  % create ferns model and compute inds (w/o field pFern, counts)
  thrs=rand(M,S)*(thrr(2)-thrr(1))+thrr(1);
  fids=uint32(floor(rand(M,S)*F+1)); inds=fernsInds(data,fids,thrs);
  ferns=struct('fids',fids,'thrs',thrs,'bayes',bayes,'H',H,'inds',inds);
else
  % re-use cached model (will need to recompute pFern)
  ferns.H=H; ferns.pFern=[]; ferns.counts = []; inds=ferns.inds; assert(size(inds,1)==N);
end

% get counts for each leaf for each class for each fern
% KB: use hist to avoid inner loop
pFern = zeros(2^S,H,M); edges = 1:2^S;
for h=1:H, inds1=inds(hs==h,:);
  for m=1:M, pFern(:,h,m)=histc(inds1(:,m),edges); end
end
pFern = pFern + bayes;

% old code
% pFern = bayes(ones(2^S,H,M));
% for n=1:N, h=hs(n);
%   for m=1:M, ind=inds(n,m);
%     pFern(ind,h,m)=pFern(ind,h,m)+1;
%   end
% end

% KB: store the unnormalized counts
ferns.counts = pFern;

% convert fern leaf class counts into probabilities
if( bayes<=0 )
  norm = 1./sum(pFern,2);
  % KB: use bsxfun instead of loop
  pFern = bsxfun(@times,pFern,norm);
  %for h=1:H, pFern(:,h,:)=pFern(:,h,:).*norm; end
else
  norm = 1./sum(pFern,1);
  % KB: use bsxfun instead of loop
  pFern = bsxfun(@times,pFern,norm);
  %for s=1:2^S, pFern(s,:,:)=pFern(s,:,:).*norm; end
  pFern=log(pFern);
end

% store pFern and compute output values
ferns.pFern=pFern; clear pFern;
if(nargout==2), hsPr=fernsClfApply([],ferns,inds); end

end

author: kristinbranson

repo: JAABA


text:
function ellipseXSection(lmsCov, planeDims, nPts, nSD, ThreeDflag)
% Plot the predictions within a plane for an ellipsoid
%
%    ellipseXSection(lmsCov, planeDims, nPts, nSD)
%
% For the L,M, set planeDims = [1,2], and for L,S set planeDims [1,3]
%
% Example:
%
%
% See also:
%
% Copyright HH, Vistalab 2010

if ieNotDefined('lmsCov'), lmsCov = eye(3,3); end
if ieNotDefined('planeDims'), planeDims = [1,2]; end
if ieNotDefined('nSD'), nSD = 2; end
if ieNotDefined('nPts'), nPts = 45; end
if ieNotDefined('ThreeDflag'), ThreeDflag = false; end

% Make points on a circle
[x,y] = circlePoints(2*pi/nPts);
spherePoints = zeros(length(x),3);
spherePoints(:,planeDims(1)) = x(:);
spherePoints(:,planeDims(2)) = y(:);

eVec = zeros(size(spherePoints,1),2);
for ii=1:size(spherePoints,1)
    l = spherePoints(ii,:)*(lmsCov\spherePoints(ii,:)');
    eVec(ii,:) = nSD*(spherePoints(ii,planeDims)/sqrt(l)); % + mn;
end

if ThreeDflag == false,
    plot(eVec(:,1),eVec(:,2),'r.'); axis equal; grid on

elseif ThreeDflag == true,
    z = zeros(size(spherePoints,1));

    switch num2str(planeDims)
        case '1  2'
            plot3(eVec(:,1),eVec(:,2),z,'r.');
        case '1  3'
            plot3(eVec(:,1),z,eVec(:,2),'r.');
        case '2  3'
            plot3(z,eVec(:,1),eVec(:,2),'r.');
        otherwise
            error('planeDims should be [1 2],[1 3],[2 3].');
    end
    axis equal; grid on
end

% Not sure if we need a square root to compute the ratio of the longest to
% shortest axes.
% axisLengths = sqrt(eig(cov(eVec)));
% cNum = axisLengths(end)/axisLengths(1);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [scores, candidateIds] = cnn_computeScores_pairwiseModel( net, imdb, getBatch, varargin )
%cnn_computeScores_pairwiseModel computes the scores using the structure network

if ~exist('varargin', 'var')
    varargin = {};
end

%% compute results w.r.t. patches
opts = struct;
opts.conserveMemory = true;
opts.sync = true;
opts.imageSet = 1 : size( imdb.imageFiles, 4 );
opts.scoreMode = 'maxMarginals';
% parse input
opts = vl_argparse(opts, varargin);

%% do the job
numImages = length( opts.imageSet );
scores = cell( max(opts.imageSet), 1);
candidateIds = cell( max(opts.imageSet), 1);

for iImageId = 1 : numImages
    iImage = opts.imageSet( iImageId );
    fprintf('Image %d/%d: ', iImageId, numImages);
    tImageStart = tic;

    [patchData, patchLabels, patchInfo] = getBatch( imdb, iImage );
    curNumCandidates = size( patchData, 4 );
    fprintf('patches=%d, ', curNumCandidates );

    [~, ~, predictions] = vl_structuredNetwork_pairwiseModel(net, patchData, [], patchLabels, [], ...
        'computeMaxMarginals', true, ...
        'disableDropout', true, ...
        'conserveMemory', opts.conserveMemory, ...
        'sync', opts.sync);
    maxMarginals = predictions{1}.maxMarginals;

    candidateIds{ iImage } = int32( patchInfo.candidateIds{1}(:) );

    switch opts.scoreMode
        case 'maxMarginals'
            curScores = maxMarginals(:, 1) - maxMarginals(:, 2);
        otherwise
            error('cnn_computeScores_pairwiseModel:unknownScoreMode', 'options scoreMode is of incorrect value');
    end

    curScores = curScores(:);
    if numel(curScores) ~= numel( candidateIds{iImage} )
        error('cnn_computeScores_pairwiseModel:inconsistentCandidates', 'something went wrong')
    end
    scores{ iImage } = single( curScores );

    fprintf('time: %fs\n', toc(tImageStart) );
end

end

author: aosokin


text:
function y = h2y(h)

% Y = h2y(H)
%
% Hybrid-H to Admittance transformation
%
% H and Y are matrices of size [2,2,F]
% where F is the number of frequencies

% 30.09.2011    - better 0-protection, faster division
%                 add multiple frequencies (old: martie 27)

if any(h(1,1,:) == 0)
    disp('h2y > correspondent admittance matrix non-existent');
else
    y = zeros(size(h));
    den = 1./h(1,1,:);
    y(1,1,:) = den;
    y(1,2,:) = -h(1,2,:).*den;
    y(2,1,:) = h(2,1,:).*den;
    y(2,2,:) = h(2,2,:) - h(1,2,:).*h(2,1,:).*den;
end

author: Sable

repo: mcbench-benchmarks


text:
function [phi, Dphi, J] = hexbasis(node,elem,GaussPt)
%% HEXBASIS the basis and their derivatives and the det of Jacobi.
%
%   [phi, Dphi, J] = hexbasis(node,elem,GaussPts) compute the values of
%   bilinear basis and their derivatives and Jacobi det at a gauss point GaussPt in a
%   hex.
%
% Author: Huayi Wei < huayiwei1984@gmail.com>
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

NT = size(elem,1);

xi = GaussPt(1);
eta = GaussPt(2);
theta = GaussPt(3);

%% basis function and their derivatives
% for a point (x, y) in real element,  there exist a point in reference
% element (xi, eta), satisfy
%    x = x_1*phi_1 + x_2*phi_2 + x_3*phi_3 + x_4 *phi_4 +
%        x_5*phi_5 +x_6*phi_6 + x_7*phi_7 + x_8 *phi_8
%    y = y_1*phi_1 + y_2*phi_2 + y_3*phi_3 + y_4 *phi_4+
%        y_5*phi_5 +y_6*phi_6 + y_7*phi_7 + y_8 *phi_8
%    z = z_1*phi_1 + z_2*phi_2 + z_3*phi_3 + z_4 *phi_4+
%        z_5*phi_5 +z_6*phi_6 + z_7*phi_7 + z_8 *phi_8

g1 = [xi, 1 - xi];
g2 = [eta, 1 - eta];
g3 = [theta, 1 - theta];

phi = zeros(8,1);
phi(1) = g1(2)*g2(2)*g3(2);
phi(2) = g1(1)*g2(2)*g3(2);
phi(3) = g1(1)*g2(1)*g3(2);
phi(4) = g1(2)*g2(1)*g3(2);
phi(5) = g1(2)*g2(2)*g3(1);
phi(6) = g1(1)*g2(2)*g3(1);
phi(7) = g1(1)*g2(1)*g3(1);
phi(8) = g1(2)*g2(1)*g3(1);

rDphi = zeros(8,3);
rDphi(1,:) = [ -g2(2)*g3(2), -g1(2)*g3(2), -g1(2)*g2(2)];
rDphi(2,:) = [  g2(2)*g3(2), -g1(1)*g3(2), -g1(1)*g2(2)];
rDphi(3,:) = [  g2(1)*g3(2),  g1(1)*g3(2), -g1(1)*g2(1)];
rDphi(4,:) = [ -g2(1)*g3(2),  g1(2)*g3(2), -g1(2)*g2(1)];
rDphi(5,:) = [ -g2(2)*g3(1), -g1(2)*g3(1),  g1(2)*g2(2)];
rDphi(6,:) = [  g2(2)*g3(1), -g1(1)*g3(1),  g1(1)*g2(2)];
rDphi(7,:) = [  g2(1)*g3(1),  g1(1)*g3(1),  g1(1)*g2(1)];
rDphi(8,:) = [ -g2(1)*g3(1),  g1(2)*g3(1),  g1(2)*g2(1)];

X = node(:,1);
Y = node(:,2);
Z = node(:,3);

JM = zeros(NT,3,3);
JM(:,:,1) = X(elem)*rDphi;
JM(:,:,2) = Y(elem)*rDphi;
JM(:,:,3) = Z(elem)*rDphi;

J = JM(:,1,1).*JM(:,2,2).*JM(:,3,3) + JM(:,2,1).*JM(:,3,2).*JM(:,1,3)+ ...
    JM(:,3,1).*JM(:,1,2).*JM(:,2,3) - JM(:,3,1).*JM(:,2,2).*JM(:,1,3)- ...
    JM(:,1,1).*JM(:,3,2).*JM(:,2,3) - JM(:,2,1).*JM(:,1,2).*JM(:,3,3);

InvJM = zeros(NT,3,3);
InvJM(:,:,1) = [JM(:,2,2).*JM(:,3,3) - JM(:,3,2).*JM(:,2,3), JM(:,3,1).*JM(:,2,3)-JM(:,2,1).*JM(:,3,3), JM(:,2,1).*JM(:,3,2) - JM(:,3,1).*JM(:,2,2)];
InvJM(:,:,2) = [JM(:,3,2).*JM(:,1,3) - JM(:,1,2).*JM(:,3,3), JM(:,1,1).*JM(:,3,3)-JM(:,3,1).*JM(:,1,3), JM(:,3,1).*JM(:,1,2) - JM(:,1,1).*JM(:,3,2)];
InvJM(:,:,3) = [JM(:,1,2).*JM(:,2,3) - JM(:,2,2).*JM(:,1,3), JM(:,2,1).*JM(:,1,3)-JM(:,1,1).*JM(:,2,3), JM(:,1,1).*JM(:,2,2) - JM(:,2,1).*JM(:,1,2)];

InvJM(:,:,1) = InvJM(:,:,1)./[J,J,J];
InvJM(:,:,2) = InvJM(:,:,2)./[J,J,J];
InvJM(:,:,3) = InvJM(:,:,3)./[J,J,J];

Dphi = zeros(NT, 3, 8);

for i = 1:8
  Dphi(1:NT,:,i) = [InvJM(:,:,1)*rDphi(i,:)', InvJM(:,:,2)*rDphi(i,:)', InvJM(:,:,3)*rDphi(i,:)'];
end

repo: ifem


text:
function [params, names] = ardKernExtractParam(kern)

% ARDKERNEXTRACTPARAM Extract parameters from the ARD kernel structure.
% FORMAT
% DESC Extract parameters from the pre-built RBF and linear ARD
% kernel structure into a vector of parameters for optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel matrix, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
%
% FORMAT
% DESC Extract parameters and names of parameters from the
% pre-built RBF and linear ARD kernel structure.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel matrix, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
% RETURN names : celly array of strings containing parameter names.
%
% SEEALSO ardKernParamInit, ardKernExpandParam, kernExtractParam, scg, conjgrad
%
% COPYRIGHT : Neil D. Lawrence, 2004
%
% KERN

params = [kern.inverseWidth kern.rbfVariance ...
          kern.biasVariance kern.whiteVariance ...
          kern.linearVariance kern.inputScales];
if nargout > 1
  names{1} = 'ard inverse width';
  names{2} = 'ard rbf variance';
  names{3} = 'ard bias variance';
  names{4} = 'ard white variance';
  names{5} = 'ard linear variance';
  for i = 1:length(kern.inputScales)
    names{5+i} = ['ard input scale ' num2str(i)];
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function y = time_varying_estimate(meth,data,varargin)
% Performs correlation operation (default) or any function you pass in (as a function handle) on symmetric moving average of data
%
% :Usage:
% ::
%
%     y = time_varying_estimate(meth,data,[window width],[function handle])
%
% - Works on all columns of the data input matrix together, so function handles can
%   operate on multivariate data (e.g. cond())
% - For moving averages, use the Matlab built-in smoothdata.m
%
% :Inputs:
%
%   **meth:**
%   Window-type options, either:
%   - 'gaussian'
%   - 'tukey'
%
% :Optional Inputs:
%
%   **stepby:**
%
%   Sometimes input data is very high resolution, and it would take too much
%   time to work element by element across the inputs.  You can enter an
%   option here to compute estimates at every n-th lag.
%
% 	**Any function handle**
%   Any function handle to a function to evaluate on local rows of [data]
%
%   e.g., @corr
%
% :Examples:
% ::
%
%    y =  time_varying_estimate('gaussian',data,20);
%
%    % Generate sample data:
%    x = mvnrnd([0 0], [1 .6; .6 1], 100);
%
%    % Correlation between columns of x:
%    r = time_varying_estimate('tukey', x, 20);
%
%    % St. deviation of first column
%    mystd = time_varying_estimate('tukey', x(:, 1), 20, @(y) std(y));
%
% ..
%    By Tor Wager
%    Last updated: Dec 2008
% ..

nshift = 0;
stepby = 1;                 % step size for shift
center_local_data = true;   % mean-center local data: good for corr, bad for moving average

% set up the kernel
% -------------------------------------------
switch meth

    case 'gaussian'

         ntrials = varargin{1};
         kern = normpdf(-3:6/ntrials:3);
         kern = kern./max(kern);            % norm to max of 1

         mymax = find(kern == max(kern));
         nshift = mymax - 1;  % kernel shifts by n points; adjust

         kern = kern';  % column

    case 'tukey'
        % Window length is the zero-influence to zero-influence time
        kern = tukeywin(varargin{1});
        nshift = round(varargin{1} ./ 2);

        % kludgey adjust in case we need extra element
        %if length(i - shift : i + nshift) > length(kern)
            kern = [kern; 0];
        %end
    otherwise error('Unknown method.')

end

% set up function handle
% -------------------------------------------
if length(varargin) > 1
    fhandle = varargin{2};
else
    fhandle = @(y) my_corrcoef(y);
end

if length(varargin) > 2
    stepby = varargin{3};
end

% set up data
% -------------------------------------------
[nobs,ncols] = size(data);

% replicate kernel for each column
kern = repmat(kern,1,ncols);

if nobs < nshift, error('Not enough observations to support kernel.'); end
y = zeros(nobs,1);

% pad data at ends to avoid edge artifacts
% -------------------------------------------
paddat = data(end:-1:end-nshift,:);

data = [data; paddat];

paddat = data(nshift:-1:1,:);

data = [paddat; data];

% execute
% -------------------------------------------
for i = [(nshift+1):stepby:(nobs + nshift) (nobs + nshift)]
    % start at nshift to avoid ends; data is padded

    % set up windowed data

    dati = data(i - nshift : i + nshift,:);

    if center_local_data
        dati = scale(dati,1) .* kern;   % center and multiply by kern so data taper towards mean
    else
        dati = dati .* kern;
    end

    % execute: IN DEVELOPMENT: THIS is hard-coded for weighted correlation
    %r = weighted_corrcoef(dati,kern(:,1));
    %r = corrcoef(dati);

    y(i, :) = fhandle(dati);  %r(1,2);

    %y(:,i) = tmpy(nshift+1:nshift+nobs);

end

if stepby == 1
    y = y( (nshift+1):(nobs + nshift) );

else
    indx = [(nshift+1):stepby:(nobs + nshift) (nobs + nshift)];
    y = interp1(indx, y(indx), (nshift+1):(nobs + nshift));

end

end

function y = my_corrcoef(dati)
r = corrcoef(dati);
y = r(1,2);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function errors = test_gsp_remove_mean()

X = rand(100,10);

X1 = gsp_remove_mean(X);

errors = gsp_assert_test(0,sum(X1),eps(1000), 'remove mean dim 1');

X2 = gsp_remove_mean(X,2);

errors = errors + gsp_assert_test(0,sum(X2,2),eps(1000), 'remove mean dim 2');

X = rand(100,10,3);

X3 = gsp_remove_mean(X,3);

errors = errors + gsp_assert_test(0,sum(X3,3),eps(1000), 'remove mean dim 3');

end

repo: gspbox


text:
%CLOSEMESS Close progress message
%
%   CLOSEMESS(FID,N)
%
% Closes a progress message of length N on file-id FID
%
% This routine is obsolete now and just preserved to get
% old code running.

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function closemess(fid,n)

	prprogress(fid,'\n');

return

author: marianux

repo: ecg-kit

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


text:
% SLAMTB_GRAPH  A graph-SLAM algorithm with simulator and graphics.
%
%   This script performs multi-robot, multi-sensor, multi-landmark 6DOF
%   graph-SLAM with simulation and graphics capabilities.
%
%   Please read slamToolbox.pdf and courseSLAM.pdf in the root directory
%   thoroughly before using this toolbox.
%
%   - Beginners should not modify this file, just edit USERDATA_GRAPH.M and
%   enter and/or modify the data you wish to simulate.
%
%   See also USERDATAGRAPH, SLAMTB.
%
%   Also consult slamToolbox.pdf and courseSLAM.pdf in the root directory.

%   Created and maintained by
%   Copyright 2008, 2009, 2010 Joan Sola @ LAAS-CNRS.
%   Copyright 2011, 2012, 2013 Joan Sola.
%   Copyright 2015-     Joan Sola @ IRI-UPC-CSIC.
%   Programmers (for parts of the toolbox):
%   Copyright David Marquez and Jean-Marie Codol @ LAAS-CNRS
%   Copyright Teresa Vidal-Calleja @ ACFR.
%   See COPYING.TXT for full copyright license.

%% OK we start here

% clear workspace and declare globals
clear
global Map

%% I. Specify user-defined options - EDIT USER DATA FILE userDataGraph.m

userDataGraph;           % user-defined data. SCRIPT.

%% II. Initialize all data structures from user-defined data
% SLAM data
[Rob,Sen,Raw,Lmk,Obs,Trj,Frm,Fac,Tim] = createGraphStructures(...
    Robot,...
    Sensor,...
    Time,...
    Opt);

% Pre-allocate all relative-motion robots:
factorRob = Rob;

% Simulation data
[SimRob,SimSen,SimLmk,SimOpt] = createSimStructures(...
    Robot,...
    Sensor,...      % all user data
    World,...
    SimOpt);

% Graphics handles
[MapFig,SenFig]               = createGraphicsStructures(...
    Rob, Sen, Lmk, Obs,...      % SLAM data
    Trj, Frm, Fac, ...
    SimRob, SimSen, SimLmk,...  % Simulator data
    FigOpt);                    % User-defined graphic options

% Clear user data - not needed anymore
clear Robot Sensor World Time   % clear all user data

%% III. Initialize data logging
% TODO: Create source and/or destination files and paths for data input and
% logs.
% TODO: do something here to collect data for post-processing or
% plotting. Think about collecting data in files using fopen, fwrite,
% etc., instead of creating large Matlab variables for data logging.

%% IV. Startup
% TODO: Possibly put in initRobots and createFrames, createFactors, createTrj...
for rob = [Rob.rob]

    % Reset relative motion robot
    factorRob(rob) = resetMotion(Rob(rob));

    % Add first keyframe with absolute factor
    Rob(rob).state.P = 1e-6 * eye(7); % Give 1cm error
    [Rob(rob),Lmk,Trj(rob),Frm(rob,:),Fac] = addKeyFrame(...
        Rob(rob),       ...
        Lmk,            ...
        Trj(rob),       ...
        Frm(rob,:),     ...
        Fac,            ...
        factorRob(rob), ...
        'absolute');

    for sen = Rob(rob).sensors

        % Initialize new landmarks
        ninits = Opt.init.nbrInits(1);
        for i = 1:ninits

            % Observe simulated landmarks
            Raw(sen) = simObservation(SimRob(rob), SimSen(sen), SimLmk, SimOpt) ;

            % Init new lmk
            fac = find([Fac.used] == false, 1, 'first');

            % Compute and allocate lmk
            [Lmk,Obs(sen,:),Frm(rob,Trj(rob).head),Fac(fac),lmk] = initNewLmk(...
                Rob(rob),   ...
                Sen(sen),   ...
                Raw(sen),   ...
                Lmk,        ...
                Obs(sen,:), ...
                Frm(rob,Trj(rob).head), ...
                Fac(fac),        ...
                Opt) ;

            if isempty(lmk)
                break
            end

        end

    end % end process sensors

end

%% V. Main loop
for currentFrame = Tim.firstFrame : Tim.lastFrame

    % 1. SIMULATION
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Simulate robots
    for rob = [SimRob.rob]

        % Robot motion
        SimRob(rob) = simMotion(SimRob(rob),Tim);

        % Simulate sensor observations
        for sen = SimRob(rob).sensors

            % Observe simulated landmarks
            Raw(sen) = simObservation(SimRob(rob), SimSen(sen), SimLmk, SimOpt) ;

        end % end process sensors

    end % end process robots

    % 2. ESTIMATION
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % 2.a. Robot motion prediction

    % Process robots
    for rob = [Rob.rob]

        % Robot motion

        % NOTE: in a regular, non-simulated SLAM, this line is not here and
        % noise just comes from the real world. Here, the estimated robot
        % is noised so that the simulated trajectory can be made perfect
        % and act as a clear reference. The noise is additive to the
        % control input 'u'.
        Rob(rob).con.u = ...
            SimRob(rob).con.u + Rob(rob).con.uStd.*randn(size(Rob(rob).con.uStd));

        Rob(rob) = simMotion(Rob(rob),Tim);

        % Integrate odometry for relative motion factors
        factorRob(rob).con.u = Rob(rob).con.u;
        factorRob(rob) = integrateMotion(factorRob(rob),Tim);

    end

    % Advance time
    Map.t = Map.t + Tim.dt;

    % 2.b. Graph construction and solving

    if mod(currentFrame - Tim.firstFrame + 1, Opt.map.kfrmPeriod) == 0

        % Process robots
        for rob = [Rob.rob]

            % Add key frame
            [Rob(rob),Lmk,Trj(rob),Frm(rob,:),Fac] = addKeyFrame(...
                Rob(rob),       ...
                Lmk,            ...
                Trj(rob),       ...
                Frm(rob,:),     ...
                Fac,            ...
                factorRob(rob), ...
                'motion');

            % Process sensor observations
            for sen = Rob(rob).sensors

                % Observe knowm landmarks
                [Rob(rob),Sen(sen),Lmk,Obs(sen,:),Frm(rob,Trj(rob).head),Fac] ...
                    = addKnownLmkFactors( ...
                    Rob(rob),   ...
                    Sen(sen),   ...
                    Raw(sen),   ...
                    Lmk,        ...
                    Obs(sen,:), ...
                    Frm(rob,Trj(rob).head), ...
                    Fac,        ...
                    Opt) ;

                % Initialize new landmarks
                ninits = Opt.init.nbrInits(1 + (currentFrame ~= Tim.firstFrame));
                for i = 1:ninits

                    % Init new lmk
                    [Lmk,Obs(sen,:),Frm(rob,Trj(rob).head),Fac,lmk] = initNewLmk(...
                        Rob(rob),   ...
                        Sen(sen),   ...
                        Raw(sen),   ...
                        Lmk,        ...
                        Obs(sen,:), ...
                        Frm(rob,Trj(rob).head), ...
                        Fac,        ...
                        Opt) ;

                    if isempty(lmk) % Did not find new lmks
                        break
                    end

                end % for i = 1:ninits

            end % end process sensors

        end % end process robots

        % Solve graph
        [Rob,Sen,Lmk,Obs,Frm,Fac] = solveGraph(Rob,Sen,Lmk,Obs,Frm,Fac,Opt);

        % Reset odometer and sync robot with graph
        for rob = [Rob.rob]

            % Update robots with Frm info
            Rob(rob) = frm2rob(Rob(rob),Frm(rob,Trj(rob).head));

            % Reset motion robot
            factorRob(rob) = resetMotion(Rob(rob));
        end

    end

    % 3. VISUALIZATION
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if currentFrame == Tim.firstFrame ...
            || currentFrame == Tim.lastFrame ...
            || mod(currentFrame,FigOpt.rendPeriod) == 0

        % Figure of the Map:
        MapFig = drawMapFig(MapFig,  ...
            Rob, Sen, Lmk,  ...
            Trj, Frm, Fac, ...
            SimRob, SimSen, ...
            FigOpt);

        if FigOpt.createVideo
            makeVideoFrame(MapFig, ...
                sprintf('map-%04d.png',currentFrame), ...
                FigOpt, ExpOpt);
        end

        % Figures for all sensors
        for sen = [Sen.sen]
            SenFig(sen) = drawSenFig(SenFig(sen), ...
                Sen(sen), Raw(sen), Obs(sen,:), ...
                FigOpt);

            if FigOpt.createVideo
                makeVideoFrame(SenFig(sen), ...
                    sprintf('sen%02d-%04d.png', sen, currentFrame),...
                    FigOpt, ExpOpt);
            end

        end

        % Do draw all objects
        drawnow;
    end

    % 4. DATA LOGGING
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % TODO: do something here to collect data for post-processing or
    % plotting. Think about collecting data in files using fopen, fwrite,
    % etc., instead of creating large Matlab variables for data logging.

end

%% VI. Post-processing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Enter post-processing code here

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
clear all; close all;

addpath Function_SPC

%% load data (you can change the file addres in this sentence.)
%addres = 'Dataset/Original/lena.bmp';
addres = 'Dataset/lena_missing_80.png';
X0 = imread(addres);
II = size(X0);

%sig = 10.0;
%noise = sig*randn(II);
%X0 = double(X0) + noise;

%% set missing indexes
Q = (X0 ~= 255);
%Q = true(II);
T = zeros(II);
T(Q) = double(X0(Q));

%T = X0;

%% hyperparameters and run SPC

%TVQV    = 'tv';        % 'tv' or 'qv' ;
%rho     = [1 1 0]; % smoothness (0.1 - 1.0) for 'qv' and (0.01 - 0.5) for 'tv' is recommended.

TVQV    = 'qv';        % 'tv' or 'qv' ;
rho     = [0.5 0.5 0]; % smoothness (0.1 - 1.0) for 'qv' and (0.01 - 0.5) for 'tv' is recommended.
K       = 10;          % Number of components which are updated in one iteration. (typically 10)
SNR     = 25;          % error bound
nu      = 0.01;        % threshold for R <-- R + 1.
maxiter = 1000;        % maximum number of iteration
tol     = 1e-5;        % tolerance
out_im  = 1;           % you can monitor the process of 'image' completion if out == 1. 'saved' directory is necessary to save the individual rank images.

[X Z G U histo histo_R] = SPC(T,Q,TVQV,rho,K,SNR,nu,maxiter,tol,out_im);
%[X Z G U histo histo_R] = SPC_beta(T,Q,'sv',rho,K,SNR,nu,maxiter,tol,out_im);
%[X Z G U histo histo_R] = SPC_old(T,Q,TVQV,rho,K,SNR,nu,maxiter,tol,out_im);
imwrite(uint8(X),'result.png');

%% visualize distribution of G and 40 parts of rank-one tensors
[val idd] = sort(G,'descend');
UU = [];
for ii = 1:length(idd)

  Z1 = outerprod(U,idd(ii));
  Z1 = Z1/max(abs(Z1(:)))*256;
  imwrite(uint8(Z1),['saved_parts/' num2str(ii) '.png']);
  UU = [UU Z1(:)];

end
figure;
subplot(2,1,1)
plot(val,'linewidth',2);
subplot(2,1,2)
P = reshape(uint8(UU(:,1:40)),[256 256 3 10 4]);
P = permute(P,[1 5 2 4 3]);
P = reshape(P,[256*4 256*10 3]);
imagesc(P)
print -depsc parts_SPC.eps

%%%%%%%%%%%%%% test for STADARD PARAFAC
addpath('/home/yokota/toolboxes/tensorlab/');
%Tl = double(T);
%Tl(~Q) = NaN;
addres = 'Dataset/Original/lena.bmp';
X0 = imread(addres);
Tl = double(X0);
Tl = fmt(Tl);

Ufr{1} = randn(256,139);
Ufr{2} = randn(256,139);
Ufr{3} = randn(3,139);
options.Display = true;
sol = cpd_als(Tl,Ufr,options);
Xfr = cpdgen(sol);
imagesc(uint8(Xfr));
%imwrite(uint8(Xfr),'result_CPD.png');

for ii = 1:139
  Gfr(ii) = norm(sol{1}(:,ii))*norm(sol{2}(:,ii))*norm(sol{3}(:,ii));
end

author: andrewssobral


text:
%% Set up paths

path(path,'mesh_functions/');
path(path,'../data');

%% Read shape

[X,T] = readOff('../data/meshes/moomoo_s0.off');
M = getMeshData(X,T,10); % compute 10 LB eigenfunctions for fun

%% Set up Gaussian blur function

blurTime = .001; % if this gets too small, distances get noisy
blurSteps = 3;

% h = blurTime/blurSteps;

blur = @(x) blurOnMesh(x,M,blurTime,blurSteps); % faster than pre-factored?
blurTranspose = @(x) blurOnMesh(x,M,blurTime,blurSteps,1);

%% Design a few functions to average

centerVerts = [300 100 600];
nFunctions = length(centerVerts);

distributions = zeros(M.numVertices,nFunctions);

for i=1:nFunctions
    distributions(centerVerts(i),i) = 1./M.areaWeights(centerVerts(i));
end

distributions = blur(distributions);

close all
for i=1:nFunctions
    f = subplot(1,nFunctions,i);
    showDescriptor(M,distributions(:,i),[],[],[],f);
    colorbar off;
    title(sprintf('Distribution %d',i));
end

%% Take the barycenter

euclideanBarycenter = sum(distributions,2)/nFunctions;

f = subplot(1,2,1);
showDescriptor(M,euclideanBarycenter,[],[],[],f);
title('Euclidean barycenter');
colorbar off;

alpha = [1 1 1];
barycenter = convolutionalBarycenter(distributions,alpha,M.areaWeights,blur,blurTranspose);

f = subplot(1,2,2);
showDescriptor(M,barycenter,[],[],[],f);
title('Wasserstein barycenter');
colorbar off;

%% Test different entropy limits

averageEntropy = -mean(sum(bsxfun(@times,distributions.*log(distributions),M.areaWeights),1));

entropyChanges = [-1.5 -1 -.5 0 .5 1 1.5];

euclideanBarycenter = sum(distributions,2)/nFunctions;

f = subplot(1,length(entropyChanges)+1,1);
showDescriptor(M,euclideanBarycenter,[],[],[],f);
title('Euclidean barycenter');
colorbar off;

for i=1:length(entropyChanges)
    targetEntropy = averageEntropy + entropyChanges(i);

    alpha = [1 1 1];
    barycenter = convolutionalBarycenter(distributions,alpha,M.areaWeights,blur,blurTranspose,targetEntropy);

    f = subplot(1,length(entropyChanges)+1,i+1);
    showDescriptor(M,barycenter,[],[],[],f);
    title(sprintf('entropy < average+(%g)',entropyChanges(i)));
    colorbar off;

    drawnow;
end

%% Try displacement interpolation

nTimeSteps = 100;

p1 = distributions(:,1);
p2 = distributions(:,2);

interp = zeros(M.numVertices,nTimeSteps);
for i=1:nTimeSteps
    fprintf('i = %d/%d\n',i,nTimeSteps);
    t = (i-1)/(nTimeSteps-1);
    alpha = [t 1-t]; % is this right?
    interp(:,i) = convolutionalBarycenter([p1 p2],alpha,M.areaWeights,blur,blurTranspose,averageEntropy);
end

%% Animate the result

f = figure;
for i=1:nTimeSteps
    clf;
    showDescriptor(M,interp(:,i),[],[],[],f);
    colorbar off;
    drawnow;
end

author: gpeyre


text:
function [x, y] = circle(senter, radius, n_points)

	%Funksjon brukes til å lage sirkler ved å bruke polarkoordinater

    s = senter;
    r = radius;
    n = n_points;

    theta = linspace(0,2*pi, n);
    rho = ones(1, n)*r;
    [x, y] = pol2cart(theta, rho);
    x = x + s(1);
    y = y + s(2);

end

author: kennydl

repo: Reinforcment-Learning-With-Q-Learning


text:
function vol = read_asa_vol(fn)

% READ_ASA_VOL reads an ASA volume conductor file
%
% all data is converted to the following units
%   vertices        mm
%   conductivities  S/m

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

Nbnd  = read_asa(fn, 'NumberBoundaries=', '%d');
UnitC = read_asa(fn, 'UnitConduct', '%s');
UnitP = read_asa(fn, 'UnitPosition', '%s');
cond  = read_asa(fn, 'Conductivities', '%f');
radii = read_asa(fn, 'Radii', '%f');
pos   = read_asa(fn, 'Positions', '%f');
bnd1  = read_asa(fn, 'Boundary1', '%s');
bnd2  = read_asa(fn, 'Boundary2', '%s');
bnd3  = read_asa(fn, 'Boundary3', '%s');
bnd4  = read_asa(fn, 'Boundary4', '%s');

if ~isempty(radii) || ~isempty(pos)
  % this appears to be a spherical volume conductor
  if strcmpi(UnitP,'mm')
    radii = 1*radii;
    pos   = 1*pos;
  elseif strcmpi(UnitP,'cm')
    radii = 100*radii;
    pos   = 100*pos;
  elseif strcmpi(UnitP,'m')
    radii = 1000*radii;
    pos   = 1000*pos;
  else
    ft_error(sprintf('Unknown unit of distance for volume (%s)', UnitP));
  end
end

if strcmpi(UnitC,'s/m')
  cond = cond/1;
elseif strcmpi(UnitC,'s/cm')
  cond = cond/100;
elseif strcmpi(UnitC,'s/mm')
  cond = cond/1000;
else
  ft_error(sprintf('Unknown unit of conductivity for volume (%s)', UnitC));
end

if ~isempty(radii)
  % this appears to be a spherical volume conductor
  vol.radius = radii;
  vol.cond   = cond;
  vol.center = pos;
else
  % this appears to be a realistical volume conductor
  [path, name, ext] = fileparts(fn);
  if Nbnd>=1
    vol.bnd(1) = read_asa_bnd(fullfile(path, bnd1));
  end
  if Nbnd>=2
    vol.bnd(2) = read_asa_bnd(fullfile(path, bnd2));
  end
  if Nbnd>=3
    vol.bnd(3) = read_asa_bnd(fullfile(path, bnd3));
  end
  if Nbnd>=4
    vol.bnd(4) = read_asa_bnd(fullfile(path, bnd4));
  end
  if Nbnd>=5
    ft_error('cannot read more than 4 boundaries');
  end

  % if there is a precomputed matrix, read it from an external file
  mat_file = read_asa(fn, 'Matrix', '%s');
  if ~isempty(mat_file)
    nr = read_asa(fullfile(path, mat_file), 'NumberRows=', '%d');
    nc = read_asa(fullfile(path, mat_file), 'NumberColumns=', '%d');
    mab_file = read_asa(fullfile(path, mat_file), 'Matrix', '%s');
    fid = fopen_or_error(fullfile(path, mab_file), 'rb', 'ieee-le');
    vol.mat = fread(fid, [nr nc], 'float32');
    fclose(fid);
    % remove the factor 2 that ASA assumes in the system matrix
    vol.mat = vol.mat/2;
    % scale the system matrix corresponding to vertex coordinates in mm
    vol.mat = vol.mat*100;
  end

  vol.cond   = cond;
end

% remove all empty fields
field=fieldnames(vol);
for i=1:length(field)
  if isempty(getfield(vol, field{i}))
    vol = rmfield(vol, field{i});
  end
end

author: spm


text:
function fx = p42_fun ( n, x )

%*****************************************************************************80
%
%% P42_FUN evaluates the integrand for problem 42.
%
%  Discussion:
%
%    The problem has a parameter ALPHA that can be set by calling
%    P42_PARAM_SET.
%
%    The integrand has a singularity at X = 0 if ALPHA < 1.
%
%    The suggested range for ALPHA is 0.1 through 2.
%
%  Interval:
%
%    0 <= x <= pi/2
%
%  Integrand:
%
%    ( sin(x) )^( alpha - 1 )
%
%  Exact Integral:
%
%    2^( alpha - 2 ) * ( Gamma(alpha/2) )^2 / Gamma(alpha)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Piessens, Elise de Doncker-Kapenga,
%    Christian Ueberhuber, David Kahaner,
%    QUADPACK: A Subroutine Package for Automatic Integration,
%    Springer, 1983, page 84.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  alpha = p42_param_get ( );

  fx = sin ( x ).^( alpha - 1.0 );

  if ( alpha < 1.0 )
    i = find ( x == 0.0 );
    fx(i) = 0.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f = p30_fun ( n, x )

%*****************************************************************************80
%
%% P30_FUN evaluates the integrand for problem 30.
%
%  Interval:
%
%    2 <= x <= 7
%
%  Integrand:
%
%          cos (       x )
%    + 5 * cos ( 1.6 * x )
%    - 2 * cos ( 2.0 * x )
%    + 5 * cos ( 4.5 * x )
%    + 7 * cos ( 9.0 * x )
%
%  Antiderivative:
%
%          sin (       x )
%    + 5 * sin ( 1.6 * x ) / 1.6
%    - 2 * sin ( 2.0 * x ) / 2.0
%    + 5 * sin ( 4.5 * x ) / 4.5
%    + 7 * sin ( 9.0 * x ) / 9.0
%
%  Exact Integral:
%
%    -4.5275696251606720278
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dianne OLeary,
%    Scientific Computing with Case Studies,
%    SIAM, 2008,
%    ISBN13: 978-0-898716-66-5,
%     LC: QA401.O44.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  f =     cos (       x ) ...
    + 5 * cos ( 1.6 * x ) ...
    - 2 * cos ( 2.0 * x ) ...
    + 5 * cos ( 4.5 * x ) ...
    + 7 * cos ( 9.0 * x );

  return
end


text:
function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
    = bcrp(fid, data, varargins, opts)
% This file is an entry for the bcrp strategy.
%
% function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
%           = bcrp_start(fid, data, varargins, opts)
%
% cum_ret: cumulative wealth achived at the end of a period.
% cumprod_ret: cumulative wealth achieved till the end each period.
% daily_ret: daily return achieved by a strategy.
%
% data: market sequence vectors
% fid: handle for write log file
% varargins: variable parameters
% opts: option parameter for behvaioral control
%
% Example: [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
%          = bcrp_start(fid, data, {0}, opts);

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

% Extract the parameters
tc = varargins{1};      % transaction cost fee rate

% Run the bcrp simulation
[cum_ret, cumprod_ret, daily_ret, daily_portfolio]...
    = bcrp_run(fid, data, tc, opts);

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function jed = yjf_to_jed_common ( y, j, f )

%*****************************************************************************80
%
%% YJF_TO_JED_COMMON converts a Common YJF date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, J, real F, the YJF date.
%
%    Output, real JED, the Julian Ephemeris Date.
%

%
%  Copy the input.
%
  y1 = y;
  j1 = j;
  f1 = f;
%
%  Check the input.
%
  [ y1, j1, f1, ierror ] = yjf_check_common ( y1, j1, f1 );

  if ( ierror ~= 0 )
    jed = -1.0;
    return
  end
%
%  Convert the input.
%
  [ y2, m2, d2, f2 ] = yjf_to_ymdf_common ( y1, j1, f1 );

  jed = ymdf_to_jed_common ( y2, m2, d2, f2 );

  return
end


text:
function [doy,yr]=jd2doy(jd);
% JD2DOY  Converts Julian date to year and day of year.
% . Non-vectorized version. See also CAL2JD, DOY2JD,
%   GPS2JD, JD2CAL, JD2DOW, JD2GPS, JD2YR, YR2JD.
% Version: 24 Apr 99
% Usage:   [doy,yr]=jd2doy(jd)
% Input:   jd  - Julian date
% Output:  doy - day of year
%          yr  - year

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

if nargin ~= 1
  warning('Incorrect number of arguments');
  return;
end
if jd < 0
  warning('Julian date must be greater than or equal to zero');
  return;
end

[yr,mn,dy] = jd2cal(jd);
doy = jd - cal2jd(yr,1,0);


text:
function [cfg, M] = ft_sourcemovie(cfg, source, source2)

% FT_SOURCEMOVIE displays the source reconstruction on a cortical mesh
% and allows the user to scroll through time with a movie.
%
% Use as
%   ft_sourcemovie(cfg, source)
% where the input source data is obtained from FT_SOURCEANALYSIS, or a
% a parcellated source structure (i.e. contains a brainordinate field) and
% cfg is a configuration structure that should contain
%
%  cfg.funparameter    = string, functional parameter that is color coded (default = 'avg.pow')
%  cfg.maskparameter   = string, functional parameter that is used for opacity (default = [])
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure.
%
% See also FT_SOURCEPLOT, FT_SOURCEINTERPOLATE, FT_SOURCEPARCELLATE

% Copyright (C) 2011-2015, Robert Oostenveld
% Copyright (C) 2012-2014, Jorn Horschig
% Copyright (C) 2018, Jan-Mathijs Schoffelen
%
% $Id$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the initial part deals with parsing the input options and data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
if ft_abort
  return
end
ft_preamble debug
ft_preamble loadvar source
ft_preamble provenance source
ft_preamble trackconfig

% ensure that the input data is valid for this function, this will also do
% backward-compatibility conversions of old data that for example was
% read from an old *.mat file
source = ft_checkdata(source, 'datatype', 'source', 'feedback', 'yes');

% get the options
xlim          = ft_getopt(cfg, 'time');
ylim          = ft_getopt(cfg, 'frequency');
zlim          = ft_getopt(cfg, 'funcolorlim');
olim          = ft_getopt(cfg, 'opacitylim');                     % don't use alim as variable name
cfg.xparam    = ft_getopt(cfg, 'xparam', 'time');                 % use time as default
cfg.yparam    = ft_getopt(cfg, 'yparam');                         % default is dealt with below
cfg.funparameter  = ft_getopt(cfg, 'funparameter', 'avg.pow');        % use power as default
cfg.funcolormap   = ft_getopt(cfg, 'funcolormap',  'jet');
cfg.maskparameter = ft_getopt(cfg, 'maskparameter');

if isempty(cfg.yparam) && isfield(source, 'freq')
  % the default is freq (if present)
  cfg.yparam = 'freq';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the actual computation is done in the middle part
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin==2
  fun = getsubfield(source, cfg.funparameter); % might be avg.pow
elseif nargin>2 && isfield(source2, 'pos')
  % in this case two conditions of data seem to have been inputted
  fun  = getsubfield(source, cfg.funparameter); % might be avg.pow
  fun2 = getsubfield(source2, cfg.funparameter);
elseif nargin>2
  % assume the first data argument to be a parcellation, and the second a
  % parcellated structure

  % THIS DOES NOT EXIST ANYMORE: IF PARCELLATED DATA IS IN THE INPUT, IT
  % WILL BE UNPARCELLATED AUTOMATICALLY BY FT_CHECKDATA
  ft_error('bogus');
end
if size(source.pos)~=size(fun,1)
  ft_error('inconsistent number of vertices in the cortical mesh');
end

if ~isfield(source, 'tri')
  ft_error('source.tri missing, this function requires a triangulated cortical sheet as source model');
end

if ~isempty(cfg.maskparameter) && ischar(cfg.maskparameter)
  mask = double(getsubfield(source, cfg.maskparameter));
else
  mask = 0.5*ones(size(fun));
end

xparam = source.(cfg.xparam);
if length(xparam)~=size(fun,2)
  error('inconsistent size of "%s" compared to "%s"', cfg.funparameter, cfg.xparam);
end

if ~isempty(cfg.yparam)
  yparam = source.(cfg.yparam);
  if length(cfg.yparam)~=size(fun,3)
    error('inconsistent size of "%s" compared to "%s"', cfg.funparameter, cfg.yparam);
  end
end

if isempty(xlim)
  xlim(1) = min(xparam);
  xlim(2) = max(xparam);
end

xbeg = nearest(xparam, xlim(1));
xend = nearest(xparam, xlim(2));
% update the configuration
cfg.xlim = xparam([xbeg xend]);

if ~isempty(cfg.yparam)
  if isempty(ylim)
    ylim(1) = min(yparam);
    ylim(2) = max(yparam);
  end
  ybeg = nearest(yparam, ylim(1));
  yend = nearest(yparam, ylim(2));
  % update the configuration
  cfg.ylim = yparam([ybeg yend]);
  hasyparam = true;
else
  % this allows us not to worry about the yparam any more
  yparam = nan;
  ybeg = 1;
  yend = 1;
  cfg.ylim = [];
  hasyparam = false;
end

% make a subselection of the data
xparam  = xparam(xbeg:xend);
yparam  = yparam(ybeg:yend);
fun     = fun(:,xbeg:xend,ybeg:yend);
if exist('fun2', 'var')
  fun2 = fun2(:,xbeg:xend,ybeg:yend);
end
mask    = mask(:,xbeg:xend,ybeg:yend);
clear xbeg xend ybeg yend

if isempty(zlim)
  zlim(1) = min(fun(:));
  zlim(2) = max(fun(:));
  % update the configuration
  cfg.funcolorlim = zlim;
end

if isempty(olim)
  olim(1) = min(mask(:));
  olim(2) = max(mask(:));
  if olim(1)==olim(2)
    olim(1) = 0;
    olim(2) = 1;
  end
  % update the configuration
  cfg.opacitylim = olim;
end

% collect the data and the options to be used in the figure
opt.cfg     = cfg;
opt.xparam  = xparam;
opt.yparam  = yparam;
opt.xval    = 0;
opt.yval    = 0;
opt.dat     = fun;
opt.mask    = mask;
opt.pos     = source.pos;
opt.tri     = source.tri;
opt.vindx   = source.inside(:);
opt.speed   = 1;
opt.record  = 0;
opt.threshold = 0;
opt.frame   = 0;
opt.cleanup = false;
if isfield(source, 'parcellation')
  opt.parcellation = source.parcellation;
  opt.parcellationlabel = source.parcellationlabel;
end

% add functional data of optional third input to the opt structure
% FIXME here we should first check whether the meshes correspond!
if nargin>2 && isfield(source2, 'pos')
  opt.dat2 = fun2;
end

% get a handle to a figure
h  = figure;
set(h, 'color', [1 1 1]);
set(h, 'toolbar', 'figure');
set(h, 'visible', 'on');
set(h, 'CloseRequestFcn', @cb_quitbutton);
set(h, 'position', [100 50 1000 600]);
set(h, 'windowbuttondownfcn', @cb_getposition);

% get timer object
t = timer;
set(t, 'timerfcn', {@cb_timer, h}, 'period', 0.1, 'executionmode', 'fixedSpacing');

% make the user interface elements
cambutton    = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'light', 'userdata', 'C');
playbutton   = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'play',   'userdata', 'p');
recordbutton = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'record', 'userdata', 'r');
quitbutton   = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'quit',   'userdata', 'q');

thrminmin    = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '-', 'userdata', 'downarrow');
thrminplus   = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '+', 'userdata', 'shift+downarrow');
thr          = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'threshold', 'userdata', 't');
thrplusmin   = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '-', 'userdata', 'shift+uparrow');
thrplusplus  = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '+', 'userdata', 'uparrow');

spdmin       = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '<', 'userdata', 'shift+downarrow');
spd          = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'speed','userdata', 's');
spdplus      = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '>', 'userdata', 'shift+uparrow');
climminmin   = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '-', 'userdata', 'leftarrow');
climminplus  = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '+', 'userdata', 'shift+leftarrow');
clim         = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', 'colorlim', 'userdata', 'z');
climplusmin  = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '-', 'userdata', 'shift+rightarrow');
climplusplus = uicontrol('parent', h, 'units', 'normalized', 'style', 'pushbutton', 'string', '+', 'userdata', 'rightarrow');
sliderx      = uicontrol('parent', h, 'units', 'normalized', 'style', 'slider',     'string', sprintf('%s = ', cfg.xparam));
stringx      = uicontrol('parent', h, 'units', 'normalized', 'style', 'text');
slidery      = uicontrol('parent', h, 'units', 'normalized', 'style', 'slider',     'string', sprintf('%s = ', cfg.yparam));
stringy      = uicontrol('parent', h, 'units', 'normalized', 'style', 'text');
stringz      = uicontrol('parent', h, 'units', 'normalized', 'style', 'text');

set(cambutton,    'position', [0.095 0.28  0.09 0.05 ], 'callback', @cb_keyboard);
set(quitbutton,   'position', [0.005 0.28  0.09 0.05 ], 'callback', @cb_keyboard);
set(playbutton,   'position', [0.005 0.22  0.09 0.05 ], 'callback', @cb_keyboard);
set(recordbutton, 'position', [0.095 0.22  0.09 0.05 ], 'callback', @cb_keyboard);
set(thrminmin,    'position', [0.005 0.16  0.03 0.025], 'callback', @cb_keyboard);
set(thrminplus,   'position', [0.005 0.185 0.03 0.025], 'callback', @cb_keyboard);
set(thr,          'position', [0.035 0.16  0.12 0.05 ], 'callback', @cb_keyboard);
set(thrplusmin,   'position', [0.155 0.16  0.03 0.025], 'callback', @cb_keyboard);
set(thrplusplus,  'position', [0.155 0.185 0.03 0.025], 'callback', @cb_keyboard);
set(climminmin,   'position', [0.005 0.10  0.03 0.025], 'callback', @cb_keyboard);
set(climminplus,  'position', [0.005 0.125 0.03 0.025], 'callback', @cb_keyboard);
set(clim,         'position', [0.035 0.10  0.12 0.05 ], 'callback', @cb_keyboard);
set(climplusmin,  'position', [0.155 0.10  0.03 0.025], 'callback', @cb_keyboard);
set(climplusplus, 'position', [0.155 0.125 0.03 0.025], 'callback', @cb_keyboard);
set(spdmin,       'position', [0.005 0.04  0.03 0.05 ], 'callback', @cb_keyboard);
set(spd,          'position', [0.035 0.04  0.12 0.05 ], 'callback', @cb_keyboard);
set(spdplus,      'position', [0.155 0.04  0.03 0.05 ], 'callback', @cb_keyboard);
set(sliderx,      'position', [0.02  0.35  0.38 0.03 ], 'callback', @cb_slider);
set(slidery,      'position', [0.200 0.005 0.78 0.03 ], 'callback', @cb_slider);
set(stringx,      'position', [0.750 0.90  0.18 0.03 ]);
set(stringy,      'position', [0.750 0.85  0.18 0.03 ]);
set(stringz,      'position', [0.60  0.95  0.33 0.03 ]);

set(stringx, 'string', sprintf('%s = ', cfg.xparam));
set(stringy, 'string', sprintf('%s = ', cfg.yparam));
set(stringz, 'string', sprintf('location = '));
set(stringx, 'horizontalalignment', 'right', 'backgroundcolor', [1 1 1]);
set(stringy, 'horizontalalignment', 'right', 'backgroundcolor', [1 1 1]);
set(stringz, 'horizontalalignment', 'right', 'backgroundcolor', [1 1 1]);

% create axes object to contain the mesh
hx = axes;
set(hx, 'position', [0.4 0.08 0.6 0.8]);
set(hx, 'tag', 'mesh');
hs = ft_plot_mesh(source, 'edgecolor', 'none', 'facecolor', [0.5 0.5 0.5], 'vertexcolor', 0.*opt.dat(:,1,1), 'facealpha', 0.*opt.mask(:,1,1), 'clim', [0 1], 'alphalim', [0 1], 'alphamap', 'rampup', 'colormap', cfg.funcolormap, 'maskstyle', 'colormix');

lighting gouraud
material dull
cam1 = camlight('left');
cam2 = camlight('right');
caxis(cfg.funcolorlim);

% create axis object to contain a time course
hy = axes;
set(hy, 'position', [0.02 0.45 0.38 0.5]);
set(hy, 'yaxislocation', 'right');

if ~hasyparam
  tline = plot(opt.xparam, mean(opt.dat(opt.vindx,:))); hold on;
  abc = axis;
  axis([opt.xparam(1) opt.xparam(end) abc(3:4)]);
  vline = plot(opt.xparam(1)*[1 1], abc(3:4), 'r');
  hline1 = plot(opt.xparam([1 end]), [0 0], 'k');
  hline2 = plot(opt.xparam([1 end]), [0 0], 'k');

  if nargin>2 && isfield(source2, 'pos')
    tline2 = plot(opt.xparam, mean(opt.dat2(opt.vindx,:)), 'r'); hold on;
  end

else
  error('not yet implemented');
end
set(hy, 'tag', 'timecourse');

% remember the various handles
opt.h   = h;  % handle to the figure
opt.hs  = hs; % handle to the mesh
opt.hx  = hx; % handle to the axes containing the mesh
opt.hy  = hy; % handle to the axes containing the timecourse
opt.cam = [cam1 cam2]; % handles to the light objects
opt.vline = vline; % handle to the vertical line in the ERF plot
opt.tline = tline; % handle to the ERF
if nargin>2 && isfield(source2, 'pos')
  opt.tline2 = tline2;
end
opt.hline1 = hline1; % handle for the horizontal line for upper threshold
opt.hline2 = hline2; % handle for hte horizontal line for lower threshold
opt.playbutton   = playbutton; % handle to the playbutton
opt.recordbutton = recordbutton; % handle to the recordbutton
opt.quitbutton   = quitbutton; % handle to the quitbutton
opt.threshold    = [0 0];

opt.t   = t;
opt.sliderx  = sliderx;
opt.slidery  = slidery;
opt.stringx  = stringx;
opt.stringy  = stringy;
opt.stringz  = stringz;

if ~hasyparam
  set(opt.slidery, 'visible', 'off');
  set(opt.stringy, 'visible', 'off');
end

setappdata(h, 'opt', opt);

while opt.cleanup==0
  uiwait(h);
  opt = getappdata(h, 'opt');
end
stop(opt.t);

if nargout
  M = opt.movie;
end

delete(h);

% do the general cleanup and bookkeeping at the end of the function
% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous source
ft_postamble provenance

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_slider(h, eventdata)

persistent previous_valx previous_valy previous_vindx previous_clim previous_thr

if isempty(previous_valx)
  previous_valx = 0;
end
if isempty(previous_valy)
  previous_valy = 0;
end
if isempty(previous_clim)
  previous_clim = [0 1];
end
if isempty(previous_thr)
  previous_thr = [0 0];
end

h    = getparent(h);
opt  = getappdata(h, 'opt');
valx = get(opt.sliderx, 'value');
valx = round(valx*(size(opt.dat,2)-1))+1;
valx = min(valx, size(opt.dat,2));
valx = max(valx, 1);

valy = get(opt.slidery, 'value');
valy = round(valy*(size(opt.dat,3)-1))+1;
valy = min(valy, size(opt.dat,3));
valy = max(valy, 1);

mask = squeeze(opt.mask(:,valx,valy));

% update stuff
if previous_valx~=valx || previous_valy~=valy || ~isequal(previous_clim, opt.cfg.funcolorlim) || ~isequal(previous_thr, opt.threshold)
  % update strings
  set(opt.stringx, 'string', sprintf('%s = %3.3f\n', opt.cfg.xparam, opt.xparam(valx)));
  set(opt.stringy, 'string', sprintf('%s = %3.3f\n', opt.cfg.yparam, opt.yparam(valy)));

  dat = opt.dat(:,valx,valy);
  mask(dat>opt.threshold(1)&dat<opt.threshold(2)) = 0;

  % convert the color-data + opacity into rgb for robust rendering
  bgcolor = repmat([0.5 0.5 0.5], [numel(mask) 1]);
  rgb     = bg_rgba2rgb(bgcolor, dat, opt.cfg.funcolormap, opt.cfg.funcolorlim, mask, 'rampup', opt.cfg.opacitylim);

  % update data in mesh
  set(opt.hs, 'FaceVertexCData', rgb, 'facecolor', 'interp');

  set(opt.vline, 'xdata', [1 1]*opt.xparam(valx));
end

% update ERF-plot
set(opt.hy,    'ylim',   opt.cfg.funcolorlim);
set(opt.vline, 'ydata',  opt.cfg.funcolorlim);
if ~(numel(previous_vindx)==numel(opt.vindx) && all(previous_vindx==opt.vindx))
  tmp = mean(opt.dat(opt.vindx,:,valy),1);
  set(opt.tline, 'ydata', tmp);
  %set(opt.hy,    'ylim',  [min(tmp(:)) max(tmp(:))]);
  %set(opt.vline, 'ydata', [min(tmp(:)) max(tmp(:))]);

  if isfield(opt, 'dat2')
    tmp = mean(opt.dat2(opt.vindx,:,valy),1);
    set(opt.tline2, 'ydata', tmp);
  end

  set(opt.hy,    'yaxislocation', 'right');
  if isfield(opt, 'parcellation')
    set(opt.stringz, 'string', sprintf('location = %s', opt.parcellationlabel{opt.parcellation(opt.vindx)}));
  else
    set(opt.stringz, 'string', sprintf('location = [%2.1f, %2.1f, %2.1f]', opt.pos(opt.vindx,:)));
  end
end

if opt.record
  tmp = get(opt.h, 'position');
  opt.frame = opt.frame + 1;
  opt.movie(opt.frame) = getframe(opt.h,[1 1 tmp(3:4)-1]);
end
setappdata(h, 'opt', opt);

previous_valx = valx;
previous_valy = valy;
previous_vindx = opt.vindx;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_playbutton(h, eventdata)

opt = getappdata(h, 'opt');
if strcmp(get(opt.playbutton, 'string'), 'pause')
  stop(opt.t);
  set(opt.playbutton, 'string', 'play');
else
  start(opt.t);
  set(opt.playbutton, 'string', 'pause');
end
setappdata(h, 'opt', opt);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_quitbutton(h, eventdata)

opt = getappdata(h, 'opt');
opt.cleanup = 1;
setappdata(h, 'opt', opt);
uiresume(h);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_recordbutton(h, eventdata)

opt = getappdata(h, 'opt');
if strcmp(get(opt.recordbutton, 'string'), 'stop')
  opt.record = 0;
  set(opt.recordbutton, 'string', 'record');
else
  opt.record = 1;
  set(opt.recordbutton, 'string', 'stop');
end
setappdata(h, 'opt', opt);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_timer(obj, info, h)

opt   = getappdata(h, 'opt');
delta = opt.speed/size(opt.dat,2);
val = get(opt.sliderx, 'value');
val = val + delta;
if val>1
  val = val-1;
end
set(opt.sliderx, 'value', val);
setappdata(h, 'opt', opt);
cb_slider(h);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_alim(h, eventdata)
if ~ishandle(h)
  return
end
opt = guidata(h);
switch get(h, 'String')
  case '+'
    alim(alim*sqrt(2));
  case '-'
    alim(alim/sqrt(2));
end % switch
guidata(h, opt);
uiresume;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_getposition(h, eventdata)

h   = getparent(h);
opt = getappdata(h, 'opt');
if strcmp(get(get(h, 'currentaxes'), 'tag'), 'timecourse')
  % get the current point
  %pos = get(opt.hy, 'currentpoint');
  %set(opt.sliderx, 'value', nearest(opt.xparam, pos(1)));
elseif strcmp(get(get(h, 'currentaxes'), 'tag'), 'mesh')
  % get the current point, which is defined as the intersection through the
  % axis-box (in 3D)
  pos       = get(opt.hx, 'currentpoint');

  % get the intersection with the mesh
  [ipos, d] = intersect_line(opt.pos, opt.tri, pos(1,:), pos(2,:));
  [md, ix]  = min(abs(d));

  if ~isempty(ix)
    dpos      = opt.pos - ipos(ix*ones(size(opt.pos,1),1),:);
    opt.vindx = nearest(sum(dpos.^2,2),0);
  end
end
setappdata(h, 'opt', opt);
cb_slider(h);
uiresume;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cb_keyboard(h, eventdata)

if (isempty(eventdata) && ft_platform_supports('matlabversion', -Inf, '2014a')) || isa(eventdata, 'matlab.ui.eventdata.ActionData')
  % determine the key that corresponds to the uicontrol element that was activated
  key = get(h, 'userdata');
else
  % determine the key that was pressed on the keyboard
  key = parsekeyboardevent(eventdata);
end
% get focus back to figure
if ~strcmp(get(h, 'type'), 'figure')
  set(h, 'enable', 'off');
  drawnow;
  set(h, 'enable', 'on');
end

h = getparent(h);
opt = getappdata(h, 'opt');

switch key
  case 'leftarrow' % change colorlim
    cval = opt.cfg.funcolorlim(1);
    if cval<0
      cval = cval.*1.1;
    else
      cval = cval./1.1;
    end
    opt.cfg.funcolorlim(1) = cval;
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hx, 'Clim', opt.cfg.funcolorlim);
  case 'rightarrow' % change colorlim
    cval = opt.cfg.funcolorlim(2);
    if cval<0
      cval = cval.*1.1;
    else
      cval = cval./1.1;
    end
    opt.cfg.funcolorlim(2) = cval;
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hx, 'Clim', opt.cfg.funcolorlim);
  case 'shift+leftarrow'
    cval = opt.cfg.funcolorlim(1);
    if cval<0
      cval = cval./1.1;
    else
      cval = cval.*1.1;
    end
    opt.cfg.funcolorlim(1) = cval;
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hx, 'Clim', opt.cfg.funcolorlim);
  case 'shift+rightarrow'
    cval = opt.cfg.funcolorlim(2);
    if cval<0
      cval = cval./1.1;
    else
      cval = cval.*1.1;
    end
    opt.cfg.funcolorlim(2) = cval;
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hx, 'Clim', opt.cfg.funcolorlim);
  case 'uparrow' % enhance threshold
    thrval = opt.threshold(2);
    if thrval==0 && opt.cfg.funcolorlim(2)>0
      opt.threshold(2) = opt.cfg.funcolorlim(2).*0.1;
    elseif thrval==0
      opt.threshold(2) = max(opt.dat(:)).*0.1;
    else
      opt.threshold(2) = thrval.*1.1;
    end
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hline1, 'YData', [1 1].*opt.threshold(2));
  case 'shift+uparrow' % enhance threshold
    thrval = opt.threshold(2);
    if thrval>0
      opt.threshold(2) = thrval.*0.9;
    end
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hline1, 'YData', [1 1].*opt.threshold(2));
  case 'downarrow' % lower threshold
    thrval = opt.threshold(1);
    if thrval==0 && opt.cfg.funcolorlim(1)<0
      opt.threshold(1) = opt.cfg.funcolorlim(1).*0.1;
    elseif thrval==0
      opt.threshold(1) = min(opt.dat(:)).*0.1;
    else
      opt.threshold(1) = thrval.*1.1;
    end
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hline2, 'YData', [1 1].*opt.threshold(1));
  case 'shift+downarrow' % lower threshold
    thrval = opt.threshold(1);
    if thrval<0
      opt.threshold(1) = thrval.*0.9;
    end
    setappdata(h, 'opt', opt);
    cb_slider(h);
    set(opt.hline2, 'YData', [1 1].*opt.threshold(1));
  case 'ctrl+uparrow' % change speed
    opt.speed = opt.speed*sqrt(2);
    setappdata(h, 'opt', opt);
  case 'ctrl+downarrow'
    opt.speed = opt.speed/sqrt(2);
    opt.speed = max(opt.speed, 1); % should not be smaller than 1
    setappdata(h, 'opt', opt);
  case 'ctrl+rightarrow' % change channel
  case 'C' % update camera position1.373e-14
    camlight(opt.cam(1), 'left');
    camlight(opt.cam(2), 'right');
  case 'p'
    cb_playbutton(h);
  case 'q'
    cb_quitbutton(h);
  case 'r'
    cb_recordbutton(h);
  case 's'
    % select the speed
    response = inputdlg('speed', 'specify', 1, {num2str(opt.speed)});
    if ~isempty(response)
      opt.speed = str2double(response);
      setappdata(h, 'opt', opt);
    end
  case 't'
    % select the threshold
    response = inputdlg('threshold', 'specify', 1, {num2str(opt.threshold)});
    if ~isempty(response)
      tok = tokenize(response{1}, ' ');
      tok(cellfun(@isempty,tok)) = [];
      for k = 1:numel(tok)
        opt.threshold(1,k) = str2double(tok{k});
      end
      setappdata(h, 'opt', opt);
    end
    set(opt.hline1, 'YData', [1 1].*opt.threshold(2));
    set(opt.hline2, 'YData', [1 1].*opt.threshold(1));
  case 'z'
    % select the colorlim
    response = inputdlg('colorlim', 'specify', 1, {num2str(opt.cfg.funcolorlim)});
    if ~isempty(response)
      tok = tokenize(response{1}, ' ');
      tok(cellfun(@isempty,tok)) = [];
      for k = 1:2
        opt.cfg.funcolorlim(1,k) = str2double(tok{k});
      end
      setappdata(h, 'opt', opt);
      cb_slider(h);
    end
  case 'control+control'
    % do nothing
  case 'shift+shift'
    % do nothing
  case 'alt+alt'
    % do nothing
  otherwise
    setappdata(h, 'opt', opt);
    %cb_help(h);
end
cb_slider(h);
uiresume(h);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = getparent(h)
p = h;
while p~=0
  h = p;
  p = get(h, 'parent');
end


text:
function start_num = p02_start_num ( )

%*****************************************************************************80
%
%% P02_START_NUM returns the number of starting points for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer START_NUM, the number of starting points.
%
  start_num = 4;

  return
end


text:
% Multi-Frame Analysis based on Discrete Cepstral Envelope (DCE-MFA)
%
% Input
%  af            : Array of structures with matrices containing sinusoidal
%                  parameters (as the output of sin_analysis.m).
%                  Each matrix is made of:
%                   The 1st line for the frequency of each sinusoid [Hz]
%                   The 2nd line for their amplitude (linear scale)
%                  The DC has to be included (in the first column)
%  fs            : [Hz] Signal's sampling frequency
%  order         : Cepstral order
%  [extrap_dcny] : If true (default), alleviate stability problems of the DCE
%                     by replacing the DCE value and extrapolating sinusoidal
%                     components up to Nyquist.
%                  If false, use the sinusoidal components as they are.
%  [scale]       : empty  : Frequency linear scale (default)
%                 'mel'  : see frq2mel (for MFCC computation)
%                 'bark' : see frq2bark
%                 'erb'  : see frq2erb
%  [Bw]          : [Hz] Standard-deviation of the Gaussian used as weighting
%                       function (for emphasizing the importance of the low
%                       frequencies in the solution).
%                       Bw As to be big enough to have the weights still
%                       significant up to Nyquist.
%  [lr]          : Regularization parameter (as in [2]) (def. 0)
%  [dftlen]      : DFT's length, if the 4th output argument is requested.
%
% Output
%  cc             : Cepstral coefficients
%  Dk             : [log] Log energy corrections
%  af             : As in input, plus the aligned amplitudes values in the
%                   'a' field.
%  E              : The amplitude cepstral envelope
%
% References
%  [1] Y. Shiga and S. King, "Estimation of voice source and vocal tract
%      characteristics based on multi-frame analysis," EUROSPEECH, 2003.
%  [2] M. Campedel-Oudot, O. Cappe and E. Moulines, "Estimation of the Spectral
%      Envelope of Voiced Sounds Using a Penalized Likelihood Approach"
%
% Copyright (c) Yannis Stylianou, 2011, Bilbao
%
% License
%  This file is part of libphoni. libphoni 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. libphoni 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.
%
% This function is part of the COVAREP project: http://covarep.github.io/covarep
%
% Authors
%  Yannis Stylianou <yannis@csd.uoc.gr>
%  Gilles Degottex <degottex@csd.uoc.gr> (regularization term)
%

function [cc Dk af E] = env_dce_mfa(af, fs, order, extrap_dcny, scale, Bw, lr, dftlen)

    debug = 0; % 0:Do nothing; 1:Plot iterations info; 2:Plot results

    % Input parameters
    if nargin<4; extrap_dcny=true; end
    if nargin<5; scale = []; end
    if nargin<6 || isempty(Bw); Bw = 2000*(fs/16000); end% 2kHz with 16kHz[1]
    if nargin<7 || isempty(lr); lr = 0; end
    if nargin<8; dftlen=4096; end

    if ~isempty(scale)
        Bw = 100*fs; % To make it as similar as possible to the DCE-SFA
        lr = 3.5e-2; % To make it as similar as possible to the DCE-SFA [2]
        eval(['fnscale=@frq2' scale ';']);
        for n=1:numel(af)
            af(n).f = 0.5*fs*fnscale(af(n).f)/fnscale(fs/2);
        end
    end

    % If asked, extrapolate components at DC and up to Nyquist
    if extrap_dcny; af = env_extrap_sins_dcny(af, fs); end

    for n=1:numel(af)
        af(n).f = af(n).sins(1,:);
        af(n).a = log(af(n).sins(2,:));
    end

    if debug>1; subplot(211); hold off; end

    M = length(af);

    BWB = 0;  % [1](9)
    erI = 0;  % Error of previous step (initialization error)
    iter = 0; % Iteration number
    cc = zeros(order+1,1);  % initialization for ceps
    rt = 0;

    for k=1:M
        fk = af(k).f;
        ak = af(k).a;
        Nk = length(fk);
        fk = fk(:);
        ak = ak(:);
        Bk = [ones(Nk,1) 2*cos(2*pi*fk/fs*(1:order))];

        wk = exp(-fk.^2 / (2 * Bw * Bw))'; % Bottom-right paragraph p.3
        Wk = diag(wk)/Nk; % (8)

        % keep the main matrices
        BWB = BWB+(Bk'*Wk*Bk); % order by order
        af(k).Bk = Bk;        % keep Bk
        af(k).Wk = Wk;        % keep Wk
        uk = ones(Nk,1);

        dk = uk'*Wk*(ak)/(uk'*Wk*uk); % (10) with c=0
        erI = erI + ((ak-dk*uk)'*Wk*(ak-dk*uk));% (7) with c=0
        rt =  rt + (Bk'*Wk*(ak-dk*uk)); % Right-hand term of (9)
    end

    if debug>0; disp(['iter:' num2str(iter) ' error=' num2str(erI) 'log']); end

    % first estimation of ceps
    if lr==0
        cc = BWB\rt; % Solution of (9)
    else
        cc = (BWB+lr*diag(ones(size(BWB,1),1)))\rt; % Solution of (9) + Regul term
    end
    cc(1) = 0;

    iter = 1;
    while(1)
        er = 0;
        rt = 0; % Right-hand term of (9)
        Dk = zeros(M,1); % log energy corrections

        for k=1:M
            ak = af(k).a;ak= ak(:);
            Bk = af(k).Bk;
            Wk = af(k).Wk;

            Nk = length(ak);
            uk = ones(Nk,1);

            dk = uk'*Wk*(ak-Bk*cc)/(uk'*Wk*uk);% (10)
            h = (ak-dk*uk-Bk*cc);
            er = er + (h'*Wk*h);% (7)
            rt =  rt + (Bk'*Wk*(ak-dk*uk)); % Right-hand term of (9)

            Dk(k) = dk;
        end
        if lr==0
            cc = BWB\rt; % Solution of (9)
        else
            cc = (BWB+lr*diag(ones(size(BWB,1),1)))\rt; % Solution of (9) + Regul term
        end
        cc(1) = 0;

        if debug>0; disp(['iter:' num2str(iter) ' error=' num2str(er) 'log']); end
    %      disp(['reldiff=' num2str(abs((erI-er)/er))]);

        if debug>1
            % plot
            subplot(211);
                plot(iter-1,log(erI), 'o');
                hold on
                plot(iter,log(er), 'x');
                title(num2str(iter));
            subplot(212);
                hold off;
                for k=1:M
                    plot(af(k).f,af(k).a-Dk(k),'+');
                    hold on;
                end
                dftlen = 2048;
                fv = fs*(0:dftlen/2)/dftlen;
                lF = 2*cos(2*pi*fv'/fs*(1:order))*cc(2:order+1);
                plot(fv,lF,'r', 'LineWidth', 2);
            keyboard
        end

        if( abs((erI-er)/er)<0.001 ) % Stop if error doesn't improve more than 0.1%
            break;
        else
            erI = er;
            iter = iter+1;
        end
    end

    % Align the gain corrections with respect to the central frame
    ci = floor((numel(af)-1)/2)+1;
    cc(1) = cc(1)+Dk(ci);
    Dk = Dk - Dk(ci);
    cc(2:end) = 2*cc(2:end);

    % Include the log energy corrections into the output
    for fi=1:numel(af)
        af(fi).a = af(fi).a - Dk(fi);
    end

    % If asked, compute the envelope
    if nargout>2
        if isempty(scale)
            E = exp(fft(cc, dftlen));
            E = E(1:end/2+1);
        else
            if strcmp(scale,'bark')
                E = barkcc2spec(cc, fs, dftlen);
                E = E(1:end/2+1);
            elseif strcmp(scale,'mel')
                E = exp(fft(cc, dftlen));
                E = E(1:dftlen/2+1);
                E = cc2hspec(cc, fs, dftlen);
                E = fwcep2hspec(cc, fs, dftlen);
            end
        end
    end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function varargout=berk2prob(varargin)
%
% For calling details please see v_berk2prob.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_berk2prob.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: berk2prob.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_berk2prob(varargin{:});
else
    v_berk2prob(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function[fs]=give_me_source_model(CL,Nc)

% Model for simple heat

fs=2+0*CL;

end

author: nickabattista


text:
function [out] = interflow_10(S,p1,p2,p3)
%interflow_10

% 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:  Scaled linear interflow if storage exceeds a threshold
% Constraints:  f = 0, for S < p2
% @(Inputs):    p1   - time coefficient [d-1]
%               p2   - threshold for flow generation [mm]
%               p3   - linear scaling parameter [-]
%               S    - current storage [mm]

out = p1*max(0,S-p2)/(p3);

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function x_it=m_st_block(x_it,mparams)
%m_st_block   block preconditioner for Stokes equations
%   x_it = m_st_block(x_it,mparams);
%   input
%          x_it         operand for preconditioning operator
%          mparams      structure defining block preconditioning operator
%   output
%          x_it         result of ilu preconditioning operation
%
%   IFISS function: DJS, HCE; 9 April 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

nv=size(mparams.Ast,1); np=size(mparams.Q,1);
rv=x_it(1:nv); rp=x_it(nv+1:nv+np);
zv=mparams.Ast\rv; zp=mparams.Q\rp;
x_it = [zv;zp];


text:
%-----------------------------------------
% UNIVERSITY OF TENAGA NASIONAL
% PROGRAMMER: MOHD ANUAR ABD AZIZ
% STUDENT ID: ME071642
% SUBJECT   : COMPUTATIONAL FLUID DYNAMICS
% ASSIGNMENT 1
% MATLAB 6.5
%16 Feb 2005
% ------------------------------------------
clc
clear all
N=200;
L=2;
dx=L/N;
n2=25;
tb=400+273.15;
ta=34+273.15;
omega=1.78;
eps=1e-6; %minimum error being set

global N omega eps L dx n2 tb ta
warning off
%-----------------------------------------
% GET NUMBER OF ITERATION FOR EVERY METHOD
% ------------------------------------------
[iter,T]=JACOBBI;
fprintf('Jacobbi method\t\t: %d iteration \n',iter);
[iter,T]=GS;
fprintf('Gauss-Seidel method\t: %d iteration \n',iter);
[iter,T]=SOR;
fprintf('SOR method\t\t\t: %d iteration \n',iter);
%----------------------------------------------------
% GENERATE ACTUAL/THEORETICAL TEMPERATURE DISTRIBUTION
% ----------------------------------------------------
x=mesh;
for i=1:N+2
    A=cosh(sqrt(n2)*(L-x(i)))/cosh(sqrt(n2)*L);
    T(i)=(tb-ta)*A+ta;
end
Ta=T(101);
plot(x,T,'b')
hold
%-----------------------------------------
% GET RESULT AND ERROR FOR EVERY NUMBER OF
% MESH AND PLOT THE GRAPH
% ------------------------------------------
N=25;
x=mesh;
[iter,T]=SOR;
plot(x,T,'g')
[y,i]=min(abs(ones(size(x))-x));
err(1)=abs(T(i)-ta)*100/ta;
errm(1)=N;

N=50;
x=mesh;
[iter,T]=SOR;
plot(x,T,'r')
[y,i]=min(abs(ones(size(x))-x));
err(2)=abs(T(i)-ta)*100/ta;
errm(2)=N;

N=75;
x=mesh;
[iter,T]=SOR;
plot(x,T,'c')
[y,i]=min(abs(ones(size(x))-x));
err(3)=abs(T(i)-ta)*100/ta;
errm(3)=N;

N=100;
x=mesh;
[iter,T]=SOR;
plot(x,T,'m')
[y,i]=min(abs(ones(size(x))-x));
err(4)=abs(T(i)-ta)*100/ta;
errm(4)=N;

N=200;
x=mesh;
[iter,T]=SOR;
plot(x,T,'y')
[y,i]=min(abs(ones(size(x))-x));
err(5)=abs(T(i)-ta)*100/ta;
errm(5)=N;

legend('Actual', '25 mesh','50 mesh','75 mesh','100 mesh','200 mesh');
hold
xlabel('Distance, X, in meter');
ylabel('Temperature, teta, in Kelvin');
title('Temperature across fin from the base x=0');
grid
figure
plot(errm,err)
xlabel('Number of mesh, N');
ylabel('% difference/error');
title('error vs number of mesh at x= 1.0 m');
grid
fprintf('\n\nNote:-\n\tOptimum omega value is %g\n\tObtained from mfile optimum.m for N=200 mesh',omega);


text:
function [ a, seed ] = wathen_gb ( nx, ny, n, seed )

%*****************************************************************************80
%
%% WATHEN_GB returns the Wathen matrix, using general banded (GB) storage.
%
%  Discussion:
%
%    The Wathen matrix is a finite element matrix which is sparse.
%
%    The entries of the matrix depend in part on a physical quantity
%    related to density.  That density is here assigned random values between
%    0 and 100.
%
%    The matrix order N is determined by the input quantities NX and NY,
%    which would usually be the number of elements in the X and Y directions.
%    The value of N is
%
%      N = 3*NX*NY + 2*NX + 2*NY + 1,
%
%    The matrix is the consistent mass matrix for a regular NX by NY grid
%    of 8 node serendipity elements.
%
%    The local element numbering is
%
%      3--2--1
%      |     |
%      4     8
%      |     |
%      5--6--7
%
%    Here is an illustration for NX = 3, NY = 2:
%
%     23-24-25-26-27-28-29
%      |     |     |     |
%     19    20    21    22
%      |     |     |     |
%     12-13-14-15-16-17-18
%      |     |     |     |
%      8     9    10    11
%      |     |     |     |
%      1--2--3--4--5--6--7
%
%    For this example, the total number of nodes is, as expected,
%
%      N = 3 * 3 * 2 + 2 * 2 + 2 * 3 + 1 = 29
%
%    The matrix is symmetric positive definite for any positive values of the
%    density RHO(X,Y).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Nicholas Higham,
%    Algorithm 694: A Collection of Test Matrices in MATLAB,
%    ACM Transactions on Mathematical Software,
%    Volume 17, Number 3, September 1991, pages 289-305.
%
%    Andrew Wathen,
%    Realistic eigenvalue bounds for the Galerkin mass matrix,
%    IMA Journal of Numerical Analysis,
%    Volume 7, Number 4, October 1987, pages 449-457.
%
%  Parameters:
%
%    Input, integer NX, NY, values which determine the size of the matrix.
%
%    Input, integer N, the number of variables.
%
%    Input/output, integer SEED, the random number seed.
%
%    Output, real A(9*NX+13,N), the matrix.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'WATHEN_GB - Fatal error!\n' );
    fprintf ( 1, '  Not enough input.\n' );
    error ( 'WATHEN_GB - Fatal error!' );
  end

  if ( nargin < 2 )
    ny = nx;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  NY was not supplied.  Setting NY = NX = %d.\n', ny );
  end

  if ( nargin < 3 )
    n = wathen_order ( nx, ny );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  N was not supplied.  N = %d\n', n );
  end

  if ( nargin < 4 )
    seed = 123456789;
  end

  ml = 3 * nx + 4;
  mu = 3 * nx + 4;

  a = zeros ( 2 * ml + mu + 1, n );

  em = [
     6.0, -6.0,  2.0, -8.0,  3.0, -8.0,  2.0, -6.0;
    -6.0, 32.0, -6.0, 20.0, -8.0, 16.0, -8.0, 20.0;
     2.0, -6.0,  6.0, -6.0,  2.0, -8.0,  3.0, -8.0;
    -8.0, 20.0, -6.0, 32.0, -6.0, 20.0, -8.0, 16.0;
     3.0, -8.0,  2.0, -6.0,  6.0, -6.0,  2.0, -8.0;
    -8.0, 16.0, -8.0, 20.0, -6.0, 32.0, -6.0, 20.0;
     2.0, -8.0,  3.0, -8.0,  2.0, -6.0,  6.0, -6.0;
    -6.0, 20.0, -8.0, 16.0, -8.0, 20.0, -6.0, 32.0 ]';

  node = zeros(8,1);

  for j = 1 : ny

    for i = 1 : nx
%
%  For the element (I,J), determine the indices of the 8 nodes.
%
      node(1) = 3 * j * nx + 2 * j + 2 * i + 1;
      node(2) = node(1) - 1;
      node(3) = node(1) - 2;

      node(4) = ( 3 * j - 1 ) * nx + 2 * j + i - 1;
      node(8) = node(4) + 1;

      node(5) = ( 3 * j - 3 ) * nx + 2 * j + 2 * i - 3;
      node(6) = node(5) + 1;
      node(7) = node(5) + 2;

      [ rho, seed ] = r8_uniform_01 ( seed );
      rho = 100.0 * rho;

      for krow = 1 : 8
        for kcol = 1 : 8
          ii = node(krow);
          jj = node(kcol);
          a(ii-jj+ml+mu+1,jj) = a(ii-jj+ml+mu+1,jj) + rho * em(krow,kcol);
        end
      end

    end
  end

  return
end


text:
function success = createMhaMask(scanNum, maskFileName, planC, mask3M, threshold_bone)

if ~exist('planC','var')
    global planC
end

indexS = planC{end};

success = 1;

try
    if exist(maskFileName,'file')
        delete(maskFileName);
    end

    if ~isempty(threshold_bone)

        % Write .mha for scanNum
        [uniformCT, uniformScanInfoS] = getUniformizedCTScan(0,scanNum,planC);
        uniformCT = permute(uniformCT, [2 1 3]);
        uniformCT = flipdim(uniformCT,3);

        % Change data type to int16 to allow (-)ve values
        uniformCT = int16(uniformCT) - int16(planC{indexS.scan}(scanNum).scanInfo(1).CTOffset);

        uniformCTMask3M = uniformCT >= threshold_bone;

    else

        [~, uniformScanInfoS] = getUniformizedCTScan(0,scanNum,planC);

    end

    if ~isempty(mask3M)
        mask3M = permute(mask3M, [2 1 3]);
        mask3M = flipdim(mask3M,3);
    end

    if ~isempty(mask3M) && ~isempty(threshold_bone)
        mask3M = mask3M | uniformCTMask3M;
    elseif isempty(mask3M) && ~isempty(threshold_bone)
        mask3M = uniformCTMask3M;
    end

    % [dx, dy, dz]
    resolution = [uniformScanInfoS.grid2Units, uniformScanInfoS.grid1Units, uniformScanInfoS.sliceThickness] * 10;

    [xVals, yVals, zVals] = getUniformScanXYZVals(planC{indexS.scan}(scanNum));

    offset = [xVals(1) -yVals(1) -zVals(end)] * 10;

    % Write .mha file for scanNum1
    writemetaimagefile(maskFileName, mask3M, resolution, offset)

catch

    success = 0;

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR