text:
function partners=assignStableRoommates(prefLists)
%%ASSIGNSTABLEROOMMATES Solve the basic form of the stable roommates
%        problem. In this problem, an even number of individuals rank their
%        preferences for each other. The algorithm then finds a stable
%        pairing of individuals, if one exists. An assignment is stable if
%        there is no pair of (non-assigned) members where BOTH prefer each
%        other over their current assignment. This function finds one
%        solutions; though multiple solutions may exist. If no solutions
%        exist, an empty matrix is returned.
%
%INPUTS: prefLists An NX(N-1) matrix where prefLists(n,i) is the ith
%                  preferred roommate for person n. N must be an even
%                  number. The "roommate" is just a number from 1 to N that
%                  cannot equal n. Alternatively, an NXN matrix can be
%                  passed, where the last column is equal to the row
%                  number. This is allowed as this algorithm just augments
%                  any NX(N-1) matrix to have that format before proceeding
%                  (it plays a role in identifying infeasible assignments).
%
%OUTPUTS: partners An NX1 vector listing which each partner of the N
%                  individuals is assigned to share a room in a stable
%                  assignment (thus, partners is a permutation vector). If
%                  no stable assignment exists, then an empty matrix is
%                  returned.
%
%The O(N^2) algorithm of [1] is used to solve the problem.
%
%Example with a solution (from [1]):
% prefLists=[4,6,2,5,3;
%            6,3,5,1,4;
%            4,5,1,6,2;
%            2,6,5,1,3;
%            4,2,3,6,1;
%            5,1,4,2,3];
% partners=assignStableRoommates(prefLists)
%The solution should be partners=[6;3;2;5;4;1] for pairs of (1,6), (2,3),
%and (4,5)
%
%Example without a solution (from [1]):
% prefLists=[2,6,4,3,5;
%            3,5,1,6,4;
%            1,6,2,5,4;
%            5,2,3,6,1;
%            6,1,3,4,2;
%            4,2,5,1,3];
% partners=assignStableRoommates(prefLists)
%The solution in this instance is an empty matrix, because no stable
%solution exists.
%
%Example with multiple solutions (from [1]). Only 1 solution is found by
%this function, though.
% prefLists=[2,5,4,6,7,8,3;
%            3,6,1,7,8,5,4;
%            4,7,2,8,5,6,1;
%            1,8,3,5,6,7,2;
%            6,1,8,2,3,4,7;
%            7,2,5,3,4,1,8;
%            8,3,6,4,1,2,5;
%            5,4,7,1,2,3,6];
% partners=assignStableRoommates(prefLists)
%The solution that will be found by this algorithm is
%partners=[4;3;2;1;6;5;8;7]. However, the other two possible solutions are
%[1,2,3,4,8,7,6,5] and [5,6,7,8,1,2,3,4].
%
%REFERENCES:
%[1] R. W. Irving, "An Efficient Algorithm for the 'Stable Rommates'
%    Problem," Journal of Algorithms, vol. 6, no. 4, pp. 577-595, Dec.
%    1985.
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    numPeople=size(prefLists,1);

    %If the preference lists do not already contain "self" (sentinel) as
    %the last column, then add it.
    if(size(prefLists,2)<numPeople)
        prefLists=[prefLists,(1:numPeople)'];
    end

    firstUnmatched=1;
    firstInCycle=1;

    %Create the ranking matrix to quickly compare preferences and
    %initialize leftmost and rightmost matrices
    ranking=zeros(numPeople,numPeople);
    %Extra entries are sentinels for procedure
    %find1stwithMoteThan1Potential
    leftmost=ones(numPeople+1,1);
    rightmost=numPeople*ones(numPeople+1,1);
    for curPerson=1:numPeople
        for rank=1:numPeople
            ranking(curPerson,prefLists(curPerson,rank))=rank;
        end
    end

    [solutionPossible,leftmost,rightmost]=phase1Reduce(prefLists,ranking,leftmost,rightmost);

    %Proper initialization unnecessary
    second(1:numPeople)=leftmost(1:numPeople)+1;

    solutionFound=false;
    tail=zeros(numPeople,1);%Initially, nothing in tail (all zeros)
    cycle=zeros(numPeople,1);
    while(solutionPossible&&solutionFound==false)
        firstUnmatched=find1stwithMoteThan1Potential(firstUnmatched,leftmost,rightmost);
        if(firstUnmatched>numPeople)
            solutionFound=true;
            break;
        else
            [cycle,firstInCycle,lastInCycle,firstUnmatched,second,tail]=seekCycle(prefLists,ranking,cycle,firstInCycle,firstUnmatched,second,tail,rightmost);

            [solutionPossible,leftmost,second,rightmost]=phase2Reduce(prefLists,ranking,cycle,leftmost,second,rightmost,firstInCycle,lastInCycle);
        end
    end

    if(solutionFound)
        partners=zeros(numPeople,1);

        for curPerson=1:numPeople
            partners(curPerson)=prefLists(curPerson,leftmost(curPerson));
        end
    else
        partners=[];%No solution found.
    end
end

function [solutionPossible,leftmost,rightmost]=phase1Reduce(prefLists,ranking,leftmost,rightmost)
    numPeople=size(prefLists,1);

    %Initially, no one proposed to (all zeros).
    setProposedTo=zeros(numPeople,1);
    for curPerson=1:numPeople
        proposer=curPerson;
        while(1)
            %Best potential partner
            nextChoice=prefLists(proposer,leftmost(proposer));
            %nextChoice holds current
            current=prefLists(nextChoice,rightmost(nextChoice));

            while(ranking(nextChoice,proposer)>ranking(nextChoice,current))
                %Proposer is rejected by nextChoice
                leftmost(proposer)=leftmost(proposer)+1;
                nextChoice=prefLists(proposer,leftmost(proposer));
                current=prefLists(nextChoice,rightmost(nextChoice));
            end
            %nextChoice holds proposer
            rightmost(nextChoice)=ranking(nextChoice,proposer);
            %and rejects current
            proposer=current;

            if(setProposedTo(nextChoice)==0)
                break;
            end
        end
        setProposedTo(nextChoice)=1;
    end

    solutionPossible=(proposer==nextChoice);
end

function firstUnmatched=find1stwithMoteThan1Potential(firstUnmatched,leftmost,rightmost)
%%Find the first person with more than one potential partner.
    while(leftmost(firstUnmatched)==rightmost(firstUnmatched))
        firstUnmatched=firstUnmatched+1;
    end
end

function [cycle,firstInCycle,lastInCycle,firstUnmatched,second,tail]=seekCycle(prefLists,ranking,cycle,firstInCycle,firstUnmatched,second,tail,rightmost)
    numPeople=size(prefLists,1);

    if(firstInCycle>1)
        %Last person in previous tail.
        person=cycle(firstInCycle-1);
        %His second choice may have to be updated
        posInCycle=firstInCycle-1;
        cycleSet=tail;
    else
        cycleSet=zeros(numPeople,1);%Nothing in the cycle set
        posInCycle=1;
        person=firstUnmatched;
    end

    %Generate sequence
    while(1)
        cycleSet(person)=1;%Add person to the cycle set.

        cycle(posInCycle)=person;
        posInCycle=posInCycle+1;
        posInList=second(person);

        while(1)%Update second choice for current person
            nextChoice=prefLists(person,posInList);
            posInList=posInList+1;

            if(ranking(nextChoice,person)<=rightmost(nextChoice))
               break;
            end;
        end

        second(person)=posInList-1;
        person=prefLists(nextChoice,rightmost(nextChoice));

        %If sequence starts to cycle.
        if(cycleSet(person)==1)
            break;
        end
    end
    lastInCycle=posInCycle-1;
    tail=cycleSet;

    %Work back to the beginning of the cycle
    while(1)
       posInCycle=posInCycle-1;
       tail(cycle(posInCycle))=0;

       if(cycle(posInCycle)==person)
           break;
       end
    end
    firstInCycle=posInCycle;
end

function [solutionPossible,leftmost,second,rightmost]=phase2Reduce(prefLists,ranking,cycle,leftmost,second,rightmost,firstInCycle,lastInCycle)
    for rank=firstInCycle:lastInCycle
        %Allow next person in cycle to be rejected
        proposer=cycle(rank);
        leftmost(proposer)=second(proposer);
        %Proper update unnaecessary at this stage
        second(proposer)=leftmost(proposer)+1;

        nextChoice=prefLists(proposer,leftmost(proposer));
        rightmost(nextChoice)=ranking(nextChoice,proposer);
        %nextChoice holds proposer.
    end

    rank=firstInCycle;
    solutionPossible=true;
    while(rank<=lastInCycle&&solutionPossible)
        %Check that no one has run out of potential partners
        proposer=cycle(rank);
        solutionPossible=(leftmost(proposer)<=rightmost(proposer));
        rank=rank+1;
    end
end

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


function [lx, r0, ux] = get_rate_grid_boundaries( model, modparam, t, gamma)
% Retrive lower and upper boundaries for variance grid, based on first two moments of variance process
% model : see below
% modparam : container with each of the model parameters required
% t: time to maturity used for the vairance grid
% gamma: grid mult param, rougly centered around +/- gamma standard deviations of the variance process

if model == 1 %CIR  (eta, theta, Rho, Sigmar, r0)
    eta = modparam.eta; theta = modparam.theta; Sigmar = modparam.Sigmar; r0 = modparam.r0;

    mu_H = exp(-eta*t)*r0 + theta*(1-exp(-eta*t));
    sig2_H = Sigmar^2/eta*r0*(exp(-eta*t)-exp(-2*eta*t)) +theta*Sigmar^2/(2*eta)*(1-2*exp(-eta*t)+exp(-2*eta*t));

elseif model == 2 %VASICEK ( in SV this is STEIN STEIN / same moments as SCOTT) ... (eta, theta, Rho, Sigmar, r0)
    eta = modparam.eta; theta = modparam.theta; Sigmar = modparam.Sigmar; r0 = modparam.r0;

    mu_H = exp(-eta*t)*r0 + theta*(1-exp(-eta*t));
    sig2_H = Sigmar^2/(2*eta)*(1 -exp(-2*eta*t));

elseif model == 3 % Merton
    lambda = modparam.lambda;  Sigmar = modparam.Sigmar; r0 = modparam.r0;

    mu_H = r0 + lambda*t;
    sig2_H = Sigmar^2*t;

elseif model == 4 % Dothan
    a = modparam.a;  Sigmar = modparam.Sigmar; r0 = modparam.r0;

    mu_H = r0*exp(a*t);
    sig2_H = r0^2*exp(2*a*t)*(exp(Sigmar^2*t) - 1);

end

ux = mu_H + gamma*sqrt(sig2_H);  %variance grid upper bound
lx = mu_H - gamma*sqrt(sig2_H);

if model == 1 || model == 4
    lx = max(0.00001,lx);  %variance grid lower bound
end

end


text:
function ret = svargplvmScales(method, model, varargin)

% SVARGPLVMSCALES A small utility to get or set the ARD weights (scales) for an svargplvm model.
%
% COPYRIGHT : Andreas C. Damianou, 2011

% VARGPLVM

if ~isfield(model, 'numModels') && isfield(model, 'M')
    model.numModels = model.M;
end

if strcmp(method, 'get')
    for i=1:model.numModels
        if strcmp(model.comp{i}.kern.type, 'rbfardjit')
            scales{i} = model.comp{i}.kern.inputScales;
        else
            scales{i} = model.comp{i}.kern.comp{1}.inputScales;
        end
    end
    ret = scales;
elseif strcmp(method, 'set')
    if nargin < 2
        error('Function requires at least 3 arguments with the set method')
    end
    scales = varargin{1};
    if length(scales) == 1
        scales = repmat(median(scales{1}), 1, model.q);
    end

    for i=1:model.numModels
        if strcmp(model.comp{i}.kern.type, 'rbfardjit')
            model.comp{i}.kern.inputScales = scales;
        else
            model.comp{i}.kern.comp{1}.inputScales = scales;
        end
    end
    ret = model;
end

author: SheffieldML

repo: vargplvm


text:
function rank = bal_seq_rank ( n, t )

%*****************************************************************************80
%
%% BAL_SEQ_RANK ranks a balanced sequence.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the number of 0's (and 1's) in the sequence.
%    N must be positive.
%
%    Input, integer T(2*N), a balanced sequence.
%
%    Output, integer RANK, the rank of the balanced sequence.
%

%
%  Check.
%
  ierror = bal_seq_check ( n, t );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BAL_SEQ_RANK - Fatal error!\n' );
    fprintf ( 1, '  The input array is illegal.\n' );
    fprintf ( 1, '  IERROR = %d\n', ierror );
    error ( 'BAL_SEQ_RANK - Fatal error!' );
  end

  y = 0;
  rank = 0;

  for x = 1 : 2 * n - 1

    if ( t(x) == 0 )
      y = y + 1;
    else
      mxy = mountain ( n, x, y + 1 );
      rank = rank + mxy;
      y = y - 1;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [err, paramsOut] = showErrors(params, gt, figName, verbose)
%% Compute calibration errors given stereo parameters from synthetic calibration.
% Error types:
% 1. General errors
%   a. Reprojection error
%   b. 3D alignment error: RMS of Euclidean distance between reconstructed and
%   ground truth point cloud.
%   c. Rotation error
%   d. Translation error
% 2. Projector intrinsics errors
% 3. Camera intrinsics errors
% See also: Calibration.poseErr, Calibration.alignmentErrorSynthetic

%% License
% ACADEMIC OR NON-PROFIT ORGANIZATION NONCOMMERCIAL RESEARCH USE ONLY
% Copyright (c) 2018 Bingyao Huang
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:

% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.

% If you publish results obtained using this software, please cite our paper.

% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
% SOFTWARE.

%%

if(nargin < 4)
    verbose = false;
end

% concatenate every plane's gt points for 3D alignment error calculation
gtPtsCamImg = cell2mat(gt.synthData.gtPtsCamImg');
gtPtsPrjImg = cell2mat(gt.synthData.gtPtsPrjImg');
gtPtsWorld =  cell2mat(gt.synthData.gtPtsWorld');

%% unpack parameters
% intrinsics
camK = params.camK;
camKc = params.camKc;
prjK = params.prjK;
prjKc = params.prjKc;

% extrinsics
R = params.R;
T = params.T;

rvec = cv.Rodrigues(R)';
%% Reconstruct gt points using calibrated parameters
% VERY IMPORTANT!!!!, undistort camera image points first
gtCamPtsUndistort = ImgProc.cvUndistortPoints(gtPtsCamImg, camK, camKc);
gtCrjPtsUndistort = ImgProc.cvUndistortPoints(gtPtsPrjImg, prjK, prjKc);

% triangulate
reconPtsWorld = Reconstruct.triangulatePoints(camK, prjK, R, T, gtCamPtsUndistort, gtCrjPtsUndistort);
paramsOut.reconPtsWorld = reconPtsWorld;

% visualize
if(verbose)
    Reconstruct.visualizePts3d(reconPtsWorld, R, T, figName);
end

%% Reprojection error (sum of cam and prj)
camPts2d = cv.projectPoints(reconPtsWorld, [0 0 0], [0 0 0], camK, 'DistCoeffs', camKc);
prjPts2d = cv.projectPoints(reconPtsWorld, rvec, T', prjK, 'DistCoeffs', prjKc);

% compute residual between reprojected pts2d and captured pts2d
camReprojRes = camPts2d - gtPtsCamImg;
prjReprojRes = prjPts2d - gtPtsPrjImg;

stereoReprojRes = [camReprojRes;prjReprojRes];

% root mean square reprojection error
rmsReprojErr = sqrt( mean(sum(stereoReprojRes.^2, 2)) );

%% 3D alignment error
rmsAlignErr = Calibration.alignmentErrorSynthetic(reconPtsWorld, gtPtsWorld, verbose, figName);

%% Parameteres estimation error
% intrinsics error
camKErr = params.camK - gt.camK;
camKcErr = params.camKc - gt.camKc;
prjKErr = params.prjK - gt.prjK;
prjKcErr = params.prjKc - gt.prjKc;

% pose error (R, T)
[rotErr, transErr] = Calibration.poseErr([params.R, params.T], [gt.R, gt.T] );

%% General errors
err.name = figName;

% reconstruction, reprojection and 3d alignment error
% err.rmsReconErr = rmsReconErr;
err.rmsReprojErr = rmsReprojErr;
err.rmsAlignErr = rmsAlignErr;
% err.mEpipolarErr = mEpipolarErr;

% extrinsic err
err.rotErr = rotErr;
err.transErr = transErr;

%% Projector errors
% projection matrix err (absolute percentage value)
err.fxErrPrj = abs(prjKErr(1,1) / gt.prjK(1,1)) * 100;
err.fyErrPrj = abs(prjKErr(2,2) / gt.prjK(2,2)) * 100;

% projection matrix err (absolute value)
err.cxErrPrj = abs(prjKErr(1,3));
err.cyErrPrj = abs(prjKErr(2,3));

% distortion factors err (absolute percentage value)
err.k1ErrPrj = abs(prjKcErr(1) / gt.prjKc(1)) * 100;
err.k2ErrPrj = abs(prjKcErr(2) / gt.prjKc(2)) * 100;
err.p1ErrPrj = abs(prjKcErr(3) / gt.prjKc(3)) * 100;
err.p2ErrPrj = abs(prjKcErr(4) / gt.prjKc(4)) * 100;

%% Camera error

err.fxErrCam = abs(camKErr(1,1) / gt.camK(1,1)) * 100;
err.fyErrCam = abs(camKErr(2,2) / gt.camK(2,2)) * 100;

% projection matrix err (absolute value)
err.cxErrCam = abs(camKErr(1,3));
err.cyErrCam = abs(camKErr(2,3));

% distortion factors err (absolute percentage value)
err.k1ErrCam = abs(camKcErr(1) / gt.camKc(1)) * 100;
err.k2ErrCam = abs(camKcErr(2) / gt.camKc(2)) * 100;
err.p1ErrCam = abs(camKcErr(3) / gt.camKc(3)) * 100;
err.p2ErrCam = abs(camKcErr(4) / gt.camKc(4)) * 100;

end

author: BingyaoHuang

repo: single-shot-pro-cam-calib


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
function testall
%TESTALL Test successively all the functions of the
%	Time-Frequency Toolbox.

%	O. Lemoine  March-September 1996.

%  Choice of the Instantaneous Amplitude
   amexpo1t
   amexpo2t
   amgausst
   amrectt
   amtriant

%  Choice of the Instantaneous Frequency
   fmconstt
   fmhypt
   fmlint
   fmodanyt
   fmpart
   fmpowert
   fmsint
   gdpowert

%  Choice of Particular Signals
   altest
   anaaskt
   anabpskt
   anafskt
   anapulst
   anaqpskt
   anasingt
   anastept
   atomst
   dopnoist
   dopplert
   klaudert
   mexhatt

%  Addition of Noise
   noisecgt
   noisecut

%  Modification
   scalet

% Processing Files

%  Time-Domain Processing
   ifestart
   instfret
   loctimet

%  Frequency-Domain Processing
   fmtt
   ifmtt
   locfreqt
   sgrpdlat

%  Linear Time-Frequency Processing
   tfrgabot
   tfrstftt

% Bilinear Time-Frequency Processing in the Cohen's Class
   tfrbjt
   tfrbudt
   tfrcwt
   tfrgrdt
   tfrmhst
   tfrmht
   tfrmmcet
   tfrpaget
   tfrpmht
   tfrppagt
   tfrpwvt
   tfrridbt
   tfrridht
   tfrridnt
   tfrridtt
   tfrrit
   tfrspt
   tfrspwvt
   tfrwvt
   tfrzamt

%  Bilinear Time-Frequency Processing in the Affine Class
   tfrbertt
   tfrdflat
   tfrscalt
   tfrspawt
   tfruntet

%  Reassigned TimeFrequency Processing
   tfrrgabt
   tfrrmsct
   tfrrpmht
   tfrrppat
   tfrrpwvt
   tfrrspt
   tfrrspwt

%  Ambiguity Functions
   ambifunt
   ambifuwt

%  PostProcessing or Help to the Interpretation
   friedmat
   holdert
   htlt
   margtfrt
   midpoitt
   momftfrt
   momttfrt
   renyit
   ridgest

% Other
   dividert
   dwindowt
   integt
   integ2dt
   izakt
   kayttht
   modulot
   oddt
   sigmergt
   zakt

author: HeLiangHIT


text:
classdef SlidingWindowSegmentation < Algorithm
    properties (Access = public)
        windowSize = 512;
        iterationSize = 256;
    end

    methods (Access = public)
        function obj = SlidingWindowSegmentation(windowSize,iterationSize)
            if nargin > 0
                obj.windowSize = windowSize;
                if nargin > 1
                    obj.iterationSize = iterationSize;
                end
            end

            obj.name = 'slidingWindow';
            obj.inputPort = DataType.kSignal;
            obj.outputPort = DataType.kSegment;
        end

        function segments = compute(obj,data)
            file = Algorithm.GetSharedContextVariable(Constants.kSharedVariableCurrentDataFile);
            nSamples = length(data);
            nSegments = int32 ((nSamples - obj.windowSize + 1) / obj.iterationSize);
            segments = repmat(Segment,1,nSegments);
            segmentsCount = 1;
            for i = 1 : obj.iterationSize : nSamples - obj.windowSize
                endSample = i + obj.windowSize - 1;
                segment = Segment(file.fileName,file.data(i:endSample,:));
                segment.startSample = i;
                segment.endSample = endSample;
                segments(segmentsCount) = segment;
                segmentsCount = segmentsCount + 1;
            end
        end

        function str = toString(obj)
            str = sprintf('%s_%d_%d',obj.name,obj.windowSize,obj.iterationSize);
        end

        function metrics = computeMetrics(obj,input)
            file = Algorithm.GetSharedContextVariable(Constants.kSharedVariableCurrentDataFile);
            nSamples = length(input);
            nSegments = int32 ((nSamples - obj.windowSize + 1) / obj.iterationSize);
            flops = (nSamples - obj.windowSize ) / obj.iterationSize;
            memory = obj.windowSize * file.numColumns * Constants.kSensorDataBytes;
            outputSize = obj.windowSize * file.numColumns * nSegments * Constants.kSensorDataBytes;
            metrics = Metric(flops,memory,outputSize);
        end

        function editableProperties = getEditableProperties(obj)
            property1 = Property('windowSize',obj.windowSize,50,1024,PropertyType.kNumber);
            property2 = Property('iterationSize',obj.iterationSize,10,1024,PropertyType.kNumber);
            editableProperties = [property1,property2];
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function x = transpose( x )

%   Disciplined convex/geometric programming information for TRANSPOSE:
%      The transpose operation may be applied to CVX variables without
%      restriction.

s = x.size_;
if length( s ) > 2,
    cvx_throw( 'Transpose of an ND array is not defined.' );
elseif any( s > 1 ),
    b = x.basis_;
    if all( s > 1 ),
        b = b(:,reshape( 1 : prod( s ), s )');
    end
    x = cvx( [ s(2), s(1) ], b );
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [ceq,ceqJac] = autoGen_cst_steplength(q1m,q2m,q4m,q5m,l1,l2,l4,l5,stepLength)
%AUTOGEN_CST_STEPLENGTH
%    [CEQ,CEQJAC] = AUTOGEN_CST_STEPLENGTH(Q1M,Q2M,Q4M,Q5M,L1,L2,L4,L5,STEPLENGTH)

%    This function was generated by the Symbolic Math Toolbox version 6.3.
%    25-Oct-2015 18:36:52

t2 = cos(q1m);
t3 = l1.*t2;
t4 = cos(q2m);
t5 = l2.*t4;
t6 = cos(q4m);
t7 = cos(q5m);
t8 = sin(q1m);
t9 = sin(q2m);
t10 = sin(q4m);
t11 = l4.*t10;
t12 = sin(q5m);
t13 = l5.*t12;
ceq = [-stepLength+t11+t13-l1.*t8-l2.*t9;t3+t5-l4.*t6-l5.*t7];
if nargout > 1
    ceqJac = reshape([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-t3,-l1.*t8,-t5,-l2.*t9,0.0,0.0,l4.*t6,t11,l5.*t7,t13,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],[2,32]);
end

author: MatthewPeterKelly

repo: OptimTraj


text:
% EEGFILTFFT -  (high|low|band)-pass filter data using inverse fft
%                 (without using the Matlab signal processing toolbox)
% Usage:
%  >> [smoothdata] = eegfiltfft(data,srate,locutoff,hicutoff);
%  >> [smoothdata] = eegfiltfft(data,srate,locutoff,hicutoff,epochframes,filtorder,revfilt);
%
% Inputs:
%   data        = (channels,frames*epochs) data to filter
%   srate       = data sampling rate (Hz)
%   locutoff    = low-edge frequency in pass band (Hz)  {0 -> lowpass}
%   hicutoff    = high-edge frequency in pass band (Hz) {0 -> highpass}
%
% Optional inputs:
%   epochframes = frames per epoch (filter each epoch separately {def/0: data is 1 epoch}
%   filtorder   = argument not used (but required for symmetry with EEGFILT function).
%   revfilt     = [0|1] reverse filter (i.e. bandpass filter to notch filter). {0}
%
% Outputs:
%    smoothdata = smoothed data
%
% Known problems:
%    The signal drop off is much smaller compared to standard filtering methods
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2003
%
% See also: EEGFILT

% inspired from a google group message
% http://groups.google.com/groups?q=without+%22the+signal+processing+toolbox%22&hl=en&lr=&ie=UTF-8&oe=UTF-8&selm=f56893ae.0311141025.3069d4f8%40posting.google.com&rnum=8

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

function smoothdata = eegfiltfft(data, fs, lowcut, highcut, epochframes, filtorder, revfilt);

    if nargin < 4
        help eegfiltfft;
    end
    [chans frames] = size(data);
    if nargin < 5 || epochframes == 0
        epochframes = frames;
    end
    if nargin < 7
        revfilt = 0;
    end

    epochs = frames/epochframes;
    if epochs ~= round(epochs)
        error('Epochframes does not divide the total number of frames');
    end
    fv=reshape([0:epochframes-1]*fs/epochframes,epochframes,1); % Frequency vector for plotting

    %figure;
    %plot(fv, 20*log(abs(X))/log(10))  % Plot power spectrum in dB scale
    %xlabel('Frequency [Hz]')
    %ylabel('Signal power [dB]')

    % find closest freq in fft decomposition
    % --------------------------------------
    if lowcut ~= 0
        [tmp idxl]=min(abs(fv-lowcut));  % Find the entry in fv closest to 5 kHz
    else
        idxl = 0;
    end
    if highcut ~= 0
        [tmp idxh]=min(abs(fv-highcut));  % Find the entry in fv closest to 5 kHz
    else
        idxh = ceil(length(fv)/2);
    end

    % filter the data
    % ---------------
    smoothdata = zeros(chans,frames);
    for e = 1:epochs                % filter each epoch, channel
        for c=1:chans
            X=fft(data(c,(e-1)*epochframes+1:e*epochframes));
            if revfilt
                X(idxl+1:idxh-1)=0;
                if mod(length(X),2) == 0
                    X(end/2:end)=0;
                else
                    X((end+1)/2:end)=0;
                end
            else
                X(1:idxl)=complex(0);
                X(end-idxl:end)=complex(0);
                X(idxh:end)=complex(0);
            end;
            smoothdata(c,(e-1)*epochframes+1:e*epochframes) = 2*real(ifft(X));
            if epochs == 1
                if rem(c,20) ~= 0
                    fprintf('.');
                else
                    fprintf('%d',c);
                end
            end
        end
    end
    fprintf('\n');

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function C=Zlevelset_color(k)

if (k==1000)
    C(1)=255;
    C(2)=0;
    C(3)=0;
end

if (k==2000)
    C(1)=0;
    C(2)=255;
    C(3)=0;
end

if (k==3000)
    C(1)=0;
    C(2)=0;
    C(3)=255;
end

if (k==4000)
    C(1)=255;
    C(2)=255;
    C(3)=0;
end

if (k==5000)
    C(1)=255;
    C(2)=0;
    C(3)=255;
end

if (k==6000)
    C(1)=0;
    C(2)=255;
    C(3)=255;
end

if (k==7000)
    C(1)=128;
    C(2)=0;
    C(3)=0;
end

if (k==8000)
    C(1)=0;
    C(2)=128;
    C(3)=0;
end

if (k==9000)
    C(1)=0;
    C(2)=0;
    C(3)=128;
end

if (k==10000)
    C(1)=128;
    C(2)=128;
    C(3)=0;
end

if (k==11000)
    C(1)=128;
    C(2)=0;
    C(3)=128;
end

if (k==12000)
    C(1)=0;
    C(2)=128;
    C(3)=128;
end

if (k==13000)
    C(1)=255;
    C(2)=128;
    C(3)=128;
end

if (k==14000)
    C(1)=255;
    C(2)=128;
    C(3)=0;
end

if (k==15000)
    C(1)=255;
    C(2)=0;
    C(3)=128;
end

if (k==16000)
    C(1)=128;
    C(2)=255;
    C(3)=128;
end

if (k==17000)
    C(1)=128;
    C(2)=255;
    C(3)=0;
end

if (k==18000)
    C(1)=0;
    C(2)=255;
    C(3)=128;
end

if (k==19000)
    C(1)=128;
    C(2)=128;
    C(3)=255;
end

if (k==20000)
    C(1)=128;
    C(2)=0;
    C(3)=255;
end

if (k==21000)
    C(1)=0;
    C(2)=128;
    C(3)=255;
end

author: Sable

repo: mcbench-benchmarks


text:
function Y=plot(varargin)
%PLOT (overloaded)

% Fast version for plotting simple PWA objects
if nargin == 1
    X = varargin{1};
    if isa(varargin{1},'sdpvar')
        if length(X) == 1
            if isequal(full(getbase(X)),[0 1])
                extstruct = yalmip('extstruct',getvariables(X));
                if ~isempty(extstruct)
                    if isequal(extstruct.fcn,'pwa_yalmip') | isequal(extstruct.fcn,'pwq_yalmip')%#ok
                        switch extstruct.arg{3}
                            case ''
                            otherwise
                                Pn = polytope; Bi = {}; Ci = {};
                                index = extstruct.arg{5};
                                for i = 1:length(extstruct.arg{1})
                                    % Pick out row
                                    for j = 1:length(extstruct.arg{1}{i}.Bi)
                                        extstruct.arg{1}{i}.Bi{j} = extstruct.arg{1}{i}.Bi{j}(index,:);
                                        extstruct.arg{1}{i}.Ci{j} = extstruct.arg{1}{i}.Ci{j}(index,:);
                                    end
                                    if isempty(extstruct.arg{1}{i}.Ai{1})
                                        Pn = [Pn extstruct.arg{1}{i}.Pn];
                                        Bi = cat(2, Bi, extstruct.arg{1}{i}.Bi);
                                        Ci = cat(2, Ci, extstruct.arg{1}{i}.Ci);
                                    else
                                        if nnz([extstruct.arg{1}{i}.Ai{:}]) == 0
                                            Pn = [Pn extstruct.arg{1}{i}.Pn];
                                            Bi = cat(2, Bi, extstruct.arg{1}{i}.Bi);
                                            Ci = cat(2, Ci, extstruct.arg{1}{i}.Ci);
                                        else

                                            hold on
                                            mpt_plotPWQ(extstruct.arg{1}{i}.Pn, ...
                                                extstruct.arg{1}{i}.Ai, ...
                                                extstruct.arg{1}{i}.Bi, ...
                                                extstruct.arg{1}{i}.Ci, []);
                                            hold off
                                        end
                                    end
                                end
                                if ~isempty(Bi),
                                    mpt_plotPWA(Pn, Bi, Ci);
                                end
                                drawnow
                                return
                        end
                    end
                end
            end
        end
    end
end

% More complex expression. Get epi-graph model
% project to our variables, and extract defining facets
if nargin == 1
    [p,Bi,Ci,Pn,Pfinal] = pwa(varargin{1});%#ok
elseif isa(varargin{2},'lmi')
    [p,Bi,Ci,Pn,Pfinal] = pwa(varargin{1},varargin{2});%#ok
else
    error('Second argument should be a domain defining SET object.');
end
if nargout>0
    Y = mpt_plotPWA(Pn,Bi,Ci);
else
    mpt_plotPWA(Pn,Bi,Ci);
end

function S = extractrow(S,index)
for i = 1:length(S)
    S{i} = S{i}(index,:);
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [footPrint, nPix]= get_footprint(xs, Ly, Lx, ops)

[Mmax, imax] = max(xs, [], 2);

yp = rem(imax-1, Ly) + 1;
xp = ceil(imax/Ly);

Nk = size(xs,1);
footPrint = zeros(Nk, 1);
nPix = zeros(Nk, 1);

d = ops.diameter;
XS = repmat([-d*2:d*2], 4*d+1, 1);
YS = XS';

ds = (XS.^2 + YS.^2).^.5;

for j = 1:Nk
    yp1 = yp(j) + YS;
    xp1 = xp(j) + XS;

    badi = (yp1<1) | (yp1>Ly) | (xp1<1) | (xp1>Lx);
    yp1 = yp1(~badi);
    xp1 = xp1(~badi);
    allds = ds(~badi);

    ind = xs(j, yp1 + (xp1-1) * Ly) > Mmax(j)/2;

    footPrint(j) = mean(allds(ind));
    nPix(j) = numel(allds);
end

end

author: cortex-lab


text:
function anonAnonBeamGeometryS = defineAnonBeamgeometry()
%
% function anonAnonBeamGeometryS = defineAnonBeamgeometry()
%
% Defines the anonymous "header" data structure. The fields containing PHI have
% been left out by default. In order to leave out additional fields, remove
% them from the list below.
%
% The following three operations are allowed per field of the input data structure:
% 1. Keep the value as is: listed as 'keep' in the list below.
% 2. Allow only the specific values: permitted values listed within a cell array.
%    If no match is found, a default anonymous string will be inserted.
% 3. Insert dummy date: listed as 'date' in the list below. Date will be
%    replaced by a dummy value of 11/11/1111.
%
% APA, 1/11/2018

anonAnonBeamGeometryS = struct( ...
                   'imageType',  {{'BEAM GEOMETRY'}}, ...
                  'beamNumber',  'keep', ...
                'beamModality',  {{'PHOTON','PROTON'}}, ...
               'beamEnergyMeV',  'keep', ...
             'beamDescription',  'keep', ...
               'RxDosePerTxGy',  'keep', ...
           'MonitorUnitsPerTx',  'keep', ...
                  'numberOfTx',  'keep', ...
             'fractionGroupID',  {{'LTLUNG_D','RTLUNG_D'}}, ...
                    'beamType',  {{'DYNAMIC'}}, ...
              'planIDOfOrigin',  'keep', ...
              'collimatorType',  'keep', ...
                'apertureType',  {{'MLCX'}}, ...
         'apertureDescription',  'keep', ...
             'collimatorAngle',  'keep', ...
                 'gantryAngle',  'keep', ...
                  'couchAngle',  'keep', ...
                   'headInOut',  {{'HFS','HFP','FFS','FFP'}}, ...
    'nominalIsocenterDistance',  'keep', ...
        'numberRepresentation',  'keep', ...
                  'apertureID',  'keep', ...
                  'wedgeAngle',  'keep', ...
          'wedgeRotationAngle',  'keep', ...
                    'arcAngle',  'keep', ...
                   'machineID',  'keep', ...
                  'beamWeight',  'keep', ...
                 'weightUnits',  'keep', ...
                 'compensator',  'keep', ...
           'compensatorFormat',  'keep', ...
                     'beamUID',  'keep', ...
               'assocBeamsUID',  'keep'  ...
           );

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function y = tapas_sgm(x, a)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

y = a./(1+exp(-x));

return;

author: translationalneuromodeling

repo: tapas


text:
function [samples, E_hist, g_mean] = gpnddisimSampleHMC(model, display, iters, options);

% GPNDDISIMSAMPLEHMC Do HMC sampling for the GPNDDISIM model.
% FORMAT
% DESC performs HMC sampling for the Gaussian process single input
% motif model for a given number of iterations.
% ARG model : the model to be optimised.
% ARG display : whether or not to display while optimisation
% proceeds, set to 2 for the most verbose and 0 for the least
% verbose.
% ARG iters : number of samples to return.
% RETURN samples : the samples.
%
% SEEALSO : hmc, gpsimCreate, gpsimGradient, gpsimObjective
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% COPYRIGHT : Antti Honkela, 2007-2013
%
% COPYRIGHT : Jaakko Peltonen, 2011

% GPNDDISIM

if nargin < 3
  iters = 2000;
  if nargin < 2
    display = 1;
  end
end

[params,paramnames] = modelExtractParam(model);

if nargin < 4,
  options = hmcDefaultOptions;
end

if isfield(model, 'uniformPriors') && model.uniformPriors,
  options.iters = iters;
  bounds = gpnddisimExtractParamTransformSettings(model);
  bounds = cat(1, bounds{:})';
  [samples, E_hist, g_mean] = ...
      myhmc_bounded('modelObjective', params,  options, ...
		    'modelGradient', bounds, model);
else
  hmcopt = optOptions;
  if display
    hmcopt(1) = display;
    if length(params) <= 100
      hmcopt(9) = 1;
    end
  end

  % Momentum persistence
  hmcopt(5) = 1;

  % Leapfrog steps
  hmcopt(7) = options.tau;

  % hmcopt(18) = epsilon, step length
  hmcopt(18) = options.epsilon;

  % Number of samples to return
  hmcopt(14) = iters;

  samples = hmc('modelObjective', params,  hmcopt, ...
		'modelGradient', model);

  E_hist = [];
  g_mean = [];
end

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [Ppad, center_pad] = padIm(P, dim, center)
%PADPSF Pad an array P with zeros to make it bigger.
%
%      Ppad = padIm(PSF, dim);
%
%  Pad P with zeros to make it an m-by-n array.
%
%  Input:
%        P  Array (matrix)
%      dim  Desired dimension of padded array.
%             If dim is a scalar, then n = m.
%  Optional Input:
%    center integer array of two values indicating location of "center"
%           of P. This is used for PSFs
%
%  Output:
%        P  Padded m-by-n array.

% Reference: See Chapter 4,
%            "Deblurring Images - Matrices, Spectra, and Filtering"
%            by P. C. Hansen, J. G. Nagy, and D. P. O'Leary,
%            SIAM, Philadelphia, 2006.

%
% Set default parameters.
%
switch nargin
    case 1
        error('Need desired dimension of padded array')
    case 3
        if length(center) == 1
            error('center should be a vector with two integers')
        else
            ci = center(1); cj = center(2);
        end
end
if length(dim) == 1
    m = dim; n = dim;
else
    m = dim(1); n = dim(2);
end
%
% Pad the with zeros.
%
[mp, np] = size(P);
Ppad = zeros(m, n);
putidx = fix((size(Ppad) - size(P))/2) + 1;
Ppad(putidx(1):putidx(1)+mp-1, putidx(2):putidx(2)+np-1) = P;

if nargin == 3
    center_pad = [ci+putidx(1)-1, cj+putidx(2)-1];
else
    center_pad = [];
end

repo: IRtools


text:
% LFDispSetup - helper function used to set up a light field display
%
% Usage:
%
%     [ImageHandle, FigureHandle] = LFDispSetup( InitialFrame )
%     [ImageHandle, FigureHandle] = LFDispSetup( InitialFrame, ScaleFactor )
%
%
% This sets up a figure for LFDispMousePan and LFDispVidCirc. New in 0.5: the display is set up in
% the curently active figure if one exists. It not, this creates a new one.  The function should handle both
% mono and colour images. You can have multiple open interactive LFDisp* figures.
%
%
% Inputs:
%
%     InitialFrame : a 2D image with which to start the display
%
% Optional Inputs:
%
%     ScaleFactor : Adjusts the size of the display -- 1 means no change, 2 means twice as big, etc.
%                   Integer values are recommended to avoid scaling artifacts.
%
% Outputs:
%
%     FigureHandle, ImageHandle : handles of the created objects
%
%
% User guide: LFToolbox.pdf
% See also:  LFDispVidCirc, LFDispMousePan

% Copyright (c) 2013-2020 Donald G. Dansereau

function [ImageHandle, FigureHandle] = LFDispSetup( InitialFrame, ScaleFactor )

%---
ScaleFactor = LFDefaultVal('ScaleFactor', 1);

%---
CurSize = [size(InitialFrame), ScaleFactor];
FigureHandle = get(0, 'CurrentFigure');

%---Check if window already existed with correct size---
RedoWindow = false;
if( ~isempty(FigureHandle) )
	PrevSize = FigureHandle.UserData;
	PrevTag = get(FigureHandle, 'Tag');
	if( strcmp(PrevTag, 'LFDisplay') || ~all(size(PrevSize)==size(CurSize)) || ~all(CurSize == PrevSize) )
		% either not tagged or incorrect size; clear the figure and rebuild it below
		clf( FigureHandle );
		RedoWindow = true;
	end
end

if( isempty(FigureHandle) || RedoWindow )
	% Get screen size
	set( 0, 'units','pixels' );
	ScreenSize = get( 0, 'screensize' );
	ScreenSize = ScreenSize(3:4);

	% Create the figure if one doesn't already exist
	if( isempty(FigureHandle) )
		FigureHandle = figure(...
			'doublebuffer','on',...
			'backingstore','off',...
			...%'menubar','none',...
			...%'toolbar','none',...
			'tag','LFDisplay');
	else
		% make sure it's set up correctly
		FigureHandle.set( 'doublebuffer','on',...
			'backingstore','off',...
			...%'menubar','none',...
			...%'toolbar','none',...
			'tag','LFDisplay');
	end

	FigureHandle.UserData = [size(InitialFrame), ScaleFactor];
	FigureHandle.NextPlot = 'add';

	% Set the axis position and size within the figure
	AxesPos = [0,0,size(InitialFrame,2),size(InitialFrame,1)];
	axes( ...
		'units','pixels',...
		'Position', AxesPos,...
		'xlimmode','manual',...
		'ylimmode','manual',...
		'zlimmode','manual',...
		'climmode','manual',...
		'alimmode','manual',...
		'layer','bottom' ...
		);

	ImageHandle = imshow(InitialFrame);

	% Scaling factor sets the size of the display
	% this has the nice side-effect of rescaling and, if needed, repositining the window
	truesize(floor(ScaleFactor*size(InitialFrame(:,:,1))));

else
	ImageHandle = findobj(FigureHandle,'type','image');
	set(ImageHandle,'cdata', InitialFrame);
end

repo: LFToolbox


text:
function [gKern, gVarmeans, gVarcovars, gInd] = linard2VardistPsi2Gradient(linard2Kern, vardist, Z, covGrad, learnInducing)

% LINARD2VARDISTPSI2GRADIENT description.

% VARGPLVM

if nargin < 5
    learnInducing = 1;
end

% inverse variances
A = linard2Kern.inputScales;

Psi1 = linard2VardistPsi1Compute(linard2Kern, vardist, Z);

Amat = sparse(diag(A));
K1 = (vardist.means'*vardist.means + diag(sum(vardist.covars,1)));
K = Amat*K1*Amat;

% TYPICAL WAY
%sumS = sum(vardist.covars,1);
%for q=1:vardist.latentDimension
%   sc = sum(vardist.means(:,q)*ones(1,size(Z,1)).*Psi1,1);
%
%   ZZ = Z(:,q)*Z(:,q)';
%
%   gKern(q)  = sum(sum( (Z(:,q)*sc + (sumS(q)*A(q))*ZZ).*covGrad ));
%
%   gVarmeans(:,q) =  A(q)*sum(Psi1*((ones(size(Z,1),1)*Z(:,q)').*covGrad),2);
%
%   gVarcovars(q) = sum(sum((ZZ.*covGrad)));
%end
%gKern = 2*gKern(:)';
%gVarmeans = 2*gVarmeans(:)';
%gVarcovars = (A.^2).*gVarcovars;
%gVarcovars = ones(size(Psi1,1),1)*gVarcovars;
%gVarcovars = gVarcovars(:)';
%gInd = covGrad*Z*K;
%gInd = 2*gInd(:)';
%%% end of typical way

% FAST WAY
AZ = Z*Amat;
AZK = AZ*K1;
gKern = sum((covGrad*Z).*AZK,1);
gVarmeans = (Psi1*(covGrad*Z))*Amat;
AZ = AZ*Amat;
gVarcovars = ones(size(Psi1,1),1)*sum((covGrad*Z).*AZ,1);
if learnInducing
    gInd = covGrad*Z*K;
else
    [M Q] = size(Z);
    gInd = zeros(M, Q);
end
%
gKern = 2*gKern(:)';
gVarmeans = 2*gVarmeans(:)';
gVarcovars = gVarcovars(:)';
gInd = 2*gInd(:)';
%%%  end of FAST WAY

%sum(sum(abs(gKern1-gKern)))
%sum(sum(abs(gVarmeans1 - gVarmeans)))
%sum(sum(abs(gVarcovars1 - gVarcovars)))
%sum(sum(abs(gInd1 - gInd)))
%pause


text:
function make_eMouseData_drift(fpath, KS2path, chanMapName, useGPU, useParPool)
% this script makes a binary file of simulated eMouse recording
% written by Jennifer Colonell, based on Marius Pachitariu's original eMouse simulator for Kilosort 1
% Adds the ability to simulate simple drift of the tissue relative to the
% probe sites.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% you can play with the parameters just below here to achieve a signal more similar to your own data!!!
norm_amp  = 1.5 * 16.7; % if 0, use amplitudes of input waveforms; if > 0, set all amplitudes to norm_amp*rms_noise
mu_mean   = .75; %0.75; % mean of mean spike amplitudes. Incoming waveforms are in uV; make <1 to make sorting harder
noise_model = 'gauss'; %'gauss' or 'fromData'; 'fromData' requires a noiseModel.mat built by make_noise_model
rms_noise = 10; % rms noise in uV. Will be added to the spike signal. 15-20 uV an OK estimate from real data
t_record  = 1200; % duration in seconds of simulation. longer is better (and slower!) (1000)
fr_bounds = [1 10]; % min and max of firing rates ([1 10])
tsmooth   = 0.5; % gaussian smooth the noise with sig = this many samples (increase to make it harder) (0.5)
chsmooth  = 0.5; % smooth the noise across channels too, with this sig (increase to make it harder) (0.5)
amp_std   = .1; % standard deviation of single spike amplitude variability (increase to make it harder, technically std of gamma random variable of mean 1) (.25)
fs_rec    = 30000; % sample rate for the for the recording. Waveforms must be sampled at a rate w/in .01% of this rate
nt        = 81; % number of timepoints expected. All waveforms must have this time window
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%drift params. See the comments in calcYPos_v2 for details
drift.addDrift = 1;
drift.sType = 'rigid';  %'rigid' or 'point';
drift.tType = 'sine';   %'exp' or 'sine'
drift.y0 = 3800;        %in um, position along probe where motion is largest
                        %y = 0 is the tip of the probe
drift.halfDistance = 1000;   %in um, distance along probe over which the motion decays
drift.amplitude = 5;        %in um for a sine wave
%                             peak variation is 2Xdrift.amplitude
drift.halfLife = 2;     %in seconds
drift.period = 600;      %in seconds
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%waveform source data
% waveform files to use for the simulation
% these come in two types:
%     3A data, used for "large" units that cover > 100 um in Y
%     Data from Kampff ultradense survey, analyzed by Nick Steimetz and
%     Susu Chen, for "small" units that cover < 100 um in Y
% Interpolation of large units uses a linear "scattered interpolant"
% Interpolation of small units uses a grid interpolant.
% currently, the signal at the simulated sites is taken as the interpolated
% field value at the center of the site. This is likely appropriate for
% mapping 3A data onto simulated 3A data (or 3B data). For creating 3A data
% from the Kampff data, averaging over site area might be more realistic.

%fileCopies specifies how many repeats of each file (to make dense data
%sets). The units will be placed at random x and y on the probe, but
%using many copies can lead to too many very similar units.

useDefault = 1;     %use waveforms from eMouse folder in KS2

if useDefault
    %get waveforms from eMouse folder in KS2
     filePath{1} = fullfile(KS2path,'eMouse_drift','kampff_St_unit_waves_allNeg_2X.mat');
     fileCopies(1) = 2;
     filePath{2} = fullfile(KS2path,'eMouse_drift','121817_SU_waves_allNeg_gridEst.mat');
     fileCopies(2) = 2;
else
    %fill in paths to waveform files
    filePath = {};
    filePath{1} = 'C:\Users\labadmin\Documents\emouse_drift\eMouse\121817_single_unit_waves_allNeg.mat';
    fileCopies(1) = 1;
    filePath{2} = 'C:\Users\labadmin\Documents\emouse_drift\UltradenseKS4Jennifer\SpikeMeanWaveforms\kampff_St_unit_waves_allNeg_2X.mat';
    fileCopies(2) = 1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%other parameters
bPair     = 0; % set to 0 for randomly distributed units, 1 for units in pairs
pairDist  = 50; % distance between paired units
bPlot     = 0; %make diagnostic plots of waveforms

% Noise can either be generated from a gaussian distribution, or modeled on
% noise from real data, matching the frequency spectrum and cross channel
% correlation. The sample noise data is taken from a 3B2 recording performed
% by Susu Chen. Note that the noise data should come from a probe with
% the same geometry as the model probe.
if ( strcmp(noise_model,'fromData') )
    if useDefault
        nmPath = [KS2path,'\eMouse_drift\','SC026_3Bstag_noiseModel.mat'];
    else
        %fill in path to desired noise model.mat file
    end
    noiseFromData = load(nmPath);
end

% Add a SYNC channel to the file
% 16 bit word with a 1 Hz square wave in 7th bit
addSYNC = false;
syncOffset = 0.232; % must be between 0 and 0.5, offset to first on edge
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rng('default');

% There are 3 seeds for the randome number generator, so some parts of the
% simulation can be fixed while others vary from run to run

% For unit placment, average amplitude, and spike times

unit_seed = 101;

% For individual spike amplitudes, but still based on the same average
% Meant to simulate the same spikes showing up in different streams

amp_seed = 101;

% For noise generation

noise_seed = 101;

% set the seed of the random number generator used for unit definition
rng(unit_seed);

%bitPerUV = 0.42667; %imec 3A or 3B, gain = 500
bitPerUV = 1.3107; %for NP 2.0, fixed gain of 80

%  load channel map file built with make_eMouseChannelMap_3A_short.m

chanMapFile = fullfile(fpath, chanMapName);
load(chanMapFile);

zeroSites = find(connected == 0);

Nchan = numel(chanMap); %physical sites on the probe
%invChanMap(chanMap) = [1:Nchan]; % invert the  channel map here--create the order in which to write output

% range of y positions to place units
minY = min(ycoords(find(connected==1)));
maxY = max(ycoords(find(connected==1)));
minX = min(xcoords(find(connected==1)));
maxX = max(xcoords(find(connected==1)));

%build the complete filePath list for waveforms
nUniqueFiles = numel(filePath);
currCount = nUniqueFiles;
for i = 1:nUniqueFiles
    for j = 1:fileCopies(i)-1
        currCount = currCount + 1;
        filePath{currCount} = filePath{i};
    end
end

nFile = length(filePath);

% for i = 1: nFile
%     fprintf('%s\n',filePath{i});
% end

NN = 0; %number of units included
uF = {}; %cell array for interpolants
uR = []; %array of max/min X and Y for points that can be interpolateed
         %structure with uR.maxX, uR.maxY, uR.minX, uR.minY
origLabel = []; %array of original labels

for fileIndex = 1:nFile
    % generate waveforms for the simulation
    uData = load(filePath{fileIndex});
    fs_diff = 100*(abs(uData.fs - fs_rec)/fs_rec);
    if (fs_diff > 0.01)
        fprintf( 'Waveform file %d has wrong sample rate.\n', fileIndex );
        fprintf( 'Skipping to next file.');
        continue
    end
    if (uData.nt ~= nt)
        fprintf( 'Waveform file %d has wrong time window.\n', fileIndex );
        fprintf( 'Skipping to next file.');
        continue
    end
    if ~( strcmp(uData.dataType, '3A') || strcmp(uData.dataType, 'UD') )
        fprintf( 'Waveform file %d has unknown sample type.\n', fileIndex );
        fprintf( 'Skipping to next file.');
        continue
    end

    nUnit = length(uData.uColl);
    unitRange = NN+1:NN+nUnit;  %for now, using all the units we have

    if strcmp(uData.dataType, '3A')
        uType(unitRange) = 0;
    end
    if strcmp(uData.dataType, 'UD')
        uType(unitRange) = 1;
    end

    %vary intensity for each unit
    if uType(NN+1) == 0
        scaleIntensity = mu_mean*(1 + (rand(nUnit,1) - 0.5));
    end
    if uType(NN+1) == 1     %these are already small, so don't make them smaller.
        scaleIntensity = (1 + (rand(nUnit,1)-0.5));
    end

   if (norm_amp > 0)
        %get min and max intensity, and normalize to norm_amp X the rms noise
        targAmp = rms_noise*norm_amp;
        for i = 1:nUnit
            maxAmp = max(max(uData.uColl(i).waves)) - min(min(uData.uColl(i).waves));
            uData.uColl(i).waves = uData.uColl(i).waves*(targAmp/maxAmp);
        end
    else
        for i = 1:nUnit
            uData.uColl(i).waves = uData.uColl(i).waves*scaleIntensity(i);
        end
   end

    %deprecated
    %allowed use of scattered interpolants (non-grid sampling of the
    %waveforms -- but the calculations are 5-10X slower than gridded
    %interpolants. For non-square sampling of waveforms (e.g. Neuropixels
    %probes) -- make a scattered interpolant and then sample on a gridded
    %interpolant to feed to the simulator.
    %for these units, create an intepolant which will be used to
    %calculate the waveform at arbitrary sites

%     if (uType(unitRange(1)) == 1)
%         [uFcurr, uRcurr] = makeGridInt( uData.uColl, nt );
%     else
%         [uFcurr, uRcurr] = makeScatInt( uData.uColl, nt );
%     end

    % with both types gridded, always make a gridded interpolant
    [uFcurr, uRcurr] = makeGridInt( uData.uColl, nt );

    %append these to the array over all units
    uF = [uF, uFcurr];
    uR = [uR, uRcurr];

    for i = unitRange
        origLabel(i) = uData.uColl(i-NN).label;
    end

    NN = NN + nUnit;
end

% calculate a size for each unit
  uSize = zeros(1,NN);
  for i = 1:NN
      uSize(i) = (uR(i).maxX - uR(i).minX) * (uR(i).maxY - uR(i).minY);
  end

% distribute units along the length the probe, either in pairs separated
% by unitDist um, or fully randomly.
% for now, keep x = the original position (i.e. don't try to recenter)
uX = zeros(1,NN);
uY = zeros(1,NN);
uX = uX + ((maxX - minX)/2 + minX);

yRange = maxY - minY;

unitOrder = randperm(NN);   %really only important for the pairs

if ~bPair
    uX(unitOrder) = (maxX - minX)/2 + minX;
    uY(unitOrder) = minY + yRange * rand(NN,1);
else
    %pairs will be placed at regular spacing (so the spacing is controlled)
    %want to avoid placing units at the ends of the probe
    endBuff = 50; %min distance in um from end of probe
    if (mod(NN,2))
        %odd number of units; pair up the first NN - 1
        %place the first units between (minY + pairDist) and maxY
        nPair = (NN-1)/2;
        % (0:nPair) = nPair + 1 positions; one extra for the loner
        pairPos = minY + endBuff + (pairDist/2) + ...
            (0:nPair)*(yRange - pairDist - 2*endBuff)/(nPair);
        uY(unitOrder(1:2:NN-2)) =  pairPos(1:nPair) + (pairDist/2);
        uY(unitOrder(2:2:NN-1)) = pairPos(1:nPair) - (pairDist/2);
        uY(unitOrder(NN)) = pairPos(nPair + 1);
    else
        %even number of units; pair them all
        nPair = NN/2;
        pairPos = minY + endBuff + (pairDist/2) + ...
            (0:nPair-1)*(yRange - pairDist - 2*endBuff)/(nPair-1);
        uY(unitOrder(1:2:NN-2)) = pairPos(1:nPair) + (pairDist/2);
        uY(unitOrder(2:2:NN-1)) = pairPos(1:nPair) - (pairDist/2);
    end
end

% calculate monitor sites for these units
monSite = zeros(1,NN);
for i = 1:NN
    [currWav, uSites] = intWav( uF{i}, uX(i), uY(i), uR(i), xcoords, ycoords, connected, nt );
%         fprintf( 'site %d: ',i);
%         for uCount = 1:length(uSites)
%             fprintf( '%d,', uSites(uCount));
%         end
%         fprintf( '\n' );
    monSite(i) = pickMonSite( currWav, uSites, connected );
    if( monSite(i) < 0 )
        fprintf( 'All sites nearby unit %d are not connected. Probably an error.\n', i);
        return;
    end
end

if (bPlot)

    %write out units and positions to console
    fprintf( 'label\torig label\ty position\tmonitor site\n' );
    for i = 1:NN
        fprintf( '%d\t%d\t%.1f\t%d\n', i, origLabel(i), uY(i), monSite(i) );
    end

    %for an overview of the units included, plot waves over the whole
    %probe, at the initial position, and shifted by 1*delta, and 2*delta

    deltaUM = 10;

    %set up waves for whole probe
    testwav = zeros(Nchan,nt);

    %for each unit, find the sites with signal, calculate the waveform based on
    %the current probe position (using the interpolant) and add to wav

    for i = 1:NN
        [currWav, uSites] = intWav( uF{i}, uX(i), uY(i), uR(i), xcoords, ycoords, connected, nt );
        testwav(uSites,:) = testwav(uSites,:) + currWav;
    end

    %calculate waveforms with units shifted by deltaUM
    testwav2 = zeros(Nchan,nt);
    for i = 1:NN
        currYPos = uY(i) + deltaUM;
        [currWav, uSites] = intWav( uF{i}, uX(i), currYPos, uR(i), xcoords, ycoords, connected, nt );
        testwav2(uSites,:) = testwav2(uSites,:) + currWav;
    end

    %calculate waveforms with units shifted by 2*deltaUM
    testwav3 = zeros(Nchan,nt);

    for i = 1:NN
        currYPos = uY(i) + 2*deltaUM;
        [currWav, uSites] = intWav( uF{i}, uX(i), currYPos, uR(i), xcoords, ycoords, connected, nt );
        testwav3(uSites,:) = testwav3(uSites,:) + currWav;
    end
    figure(1);
    tRange = 1.1*nt;
    yRange = 1.1*( max(max(testwav)) - min(min(testwav)) );

    %properties of the 3A probe.
    %TODO: add derivation of these from xcoords, ycoords
    xMin = 11;
    yMin = 20;
    xStep = 16;
    yStep = 20;

    tDat = 1:nt;
    %only plotting 1-384 (rather than 1-385) because 385 is unconnected
    %digital channel

    for i = 1:Nchan-1
        currDat = testwav( i, : );
        currDat2 = testwav2(i,:);
        currDat3 = testwav3(i,:);
        xPlotPos = (xcoords(i) - xMin)/xStep;
        xOff = tRange * xPlotPos + nt/3;
        yOff = yRange * (ycoords(i) - yMin)/yStep;
        figure(1);
        %plot(tDat + xOff, currDat + yOff,'b-');
        %plot(tDat + xOff, scatDat + yOff,'b-');
        plot(tDat + xOff, currDat + yOff,'b-', tDat + xOff, currDat3 + yOff, 'r-' );
        if(xPlotPos == 0)
          msgstr=sprintf('ch%d',i);
          text(xOff,yOff+50,msgstr);
        end
        hold on
    end

    %plot signals at monitor sites at initial position.
    figure(2);

    %some plotting params
    colorStr ={};
    colorStr{1} = '-b';
    colorStr{2} = '-r';
    %fprintf('unit\tmaxAmp\n');
    for i = 1:NN
        [currWav, uSites] = intWav( uF{i}, uX(i), uY(i), uR(i), xcoords, ycoords, connected, nt );
        cM = find(uSites == monSite(i));
        currDat = currWav(cM,:);
        %fprintf( '%d\t%.2f\n', i, (max(currDat)-min(currDat)));
        currColor = colorStr{uType(i) + 1};
        plot(tDat,currDat,currColor);
        hold on;
    end

end

bContinue = 1;

% set the sample rate to that specified in the hard coded params in this
% file (independent of fs read in through channel map)
% allows simulation of multiple streams with slightly different clock rates

fs = fs_rec;
fs_std = 30000; %used to generate spike times

%same for range of firing rates (note that we haven't included any info
%about the original firign rates of the units
fr = fr_bounds(1) + (fr_bounds(2)-fr_bounds(1)) * rand(NN,1); % create variability in firing rates

% totfr = sum(fr); % total firing rate

spk_times = [];
clu = [];

if bContinue        %done with setup, now starting the time consuming stuff

for j = 1:length(fr)      %loop over neurons
    %generate a set of time differences bewteen spikes
    %random numbers from a geometric distribution with with probability
    %(sample time)*firing rate = (1/sample rate)*firing rate =
    %(1/(samplerate*firingrate))
    %geometric distribution an appropriate model for the number of trials
    %before a success (neuron firing)
    %second two params for geornd are size of the array, here 2*firing
    %rate*total time of the simulation.

    dspks = int64(geornd(1/(fs_std/fr(j)), ceil(2*fr(j)*t_record),1));
    dspks(dspks<ceil(fs_std * 2/1000)) = [];  % remove ISIs below the refractory period
    res = cumsum(dspks);
    spk_times = cat(1, spk_times, res);
    clu = cat(1, clu, j*ones(numel(res), 1));
end
% convert spike times to the requested rate
spk_times = int64( double(spk_times)*(fs/fs_std));

[spk_times, isort] = sort(spk_times);
clu = clu(isort);
clu       = clu(spk_times<t_record*fs);
spk_times = spk_times(spk_times<t_record*fs);
nspikes = numel(spk_times);

% this generates single spike amplitude with mean = 1 and std deviation
% ~amp_std, while ensuring all values are positive
% re-seed the random number generator to allow variable amplitudes
% while holding positions of units constant
rng(amp_seed);
amps = gamrnd(1/amp_std^2,amp_std^2, nspikes,1);

%
buff = 128;
NT   = 4 * fs + buff; % batch size + buffer

fidW     = fopen(fullfile(fpath, 'sim_binary.imec.ap.bin'), 'w');

t_all    = 0;

if useGPU
    %set up the random number generators for run to run reproducibility
    gpurng('default'); % start with a default set
    gpurng(noise_seed); % set the seed
    %gpurng('shuffle'); uncomment to have the seed set using the clock
end

% reset random number generator for noise
rng(noise_seed);

%record a y position for each spike time
yDriftRec = zeros( length(spk_times), 5, 'double' );
allspks = 0;

% The parpool option can speed up the calculation when using scattered
% interpolants AND running with a large number of workers (>8). With all
% gridded interpolants, the overhead is too large and
if (useParPool)
    %delete any currently running pool
    delete(gcp('nocreate'))
    %create parallel pool
    locCluster = parcluster('local');
    parpool('local',locCluster.NumWorkers);
    %get a handle to it
    p = gcp();
    %add variables to the workers
    p_xcoords = parallel.pool.Constant(xcoords);
    p_ycoords = parallel.pool.Constant(ycoords);
    p_connected = parallel.pool.Constant(connected);
    p_uF = parallel.pool.Constant(uF);
    p_uX = parallel.pool.Constant(uX);
    p_uY = parallel.pool.Constant(uY);
    p_uR = parallel.pool.Constant(uR);
end

while t_all<t_record

    if ( strcmp(noise_model,'gauss') )
        if useGPU
            enoise = gpuArray.randn(NT, Nchan, 'single');
        else
            enoise = randn(NT, Nchan, 'single');
        end
        if t_all>0
            enoise(1:buff, :) = enoise_old(NT-buff + [1:buff], :);
        end

        dat = enoise;
        dat = my_conv2(dat, [tsmooth chsmooth], [1 2]);
        %rescale the smoothed data to make the std = 1;
        dat = zscore(dat, 1, 1);
        dat = gather_try(dat);
        %multiply the final noise calculation by the expected rms in uV
        dat = dat*rms_noise;
        %fprintf( 'Noise mean = %.3f; std = %.3f\n', mean(dat(:,1)), std(dat(:,1)));
    elseif ( strcmp(noise_model,'fromData') )
        enoise = makeNoise( NT, noiseFromData, chanMap, connected, NchanTOT );
        if t_all>0
            enoise(1:buff, :) = enoise_old(NT-buff + [1:buff], :);
        end
        dat = enoise;
    end

    if addSYNC
        sync = zeros(NT,1,'int16');
        hilo = zeros(NT,1,'logical');
        % for each time point, calculate where it is in the 1Hz cycle
        currTime = zeros(NT,1,'single');
        currTime(:,1) = (0:NT-1) + t_all*fs;
        currTime(:,1) = currTime(:,1)/fs - syncOffset;  %time in seconds, relative to offset
        hilo(:,1) = currTime(:,1) - floor(currTime(:,1)) < 0.5;
        sync(hilo,1) = 64;
    end

    if t_all>0
        dat(1:buff/2, :) = dat_old(NT-buff/2 + [1:buff/2], :);
    end

    % now we add spikes all channels; ref channels zeroed out after
    ibatch = (spk_times >= t_all*fs) & (spk_times < t_all*fs+NT-buff);
    ts = spk_times(ibatch) - t_all*fs;
    ids = clu(ibatch);
    am = amps(ibatch);
    tRange = int64(1:nt);
    tic

    if (useParPool)
        %     %first run a parfor loop for the time consuming calculation of the
        %     %interpolated waves

        currWavArray = zeros(length(ts),Nchan,nt,'double'); %actual array much shorter
        currSiteArray = zeros(length(ts),Nchan,'uint16'); %record indices of the sites in currWavArray;
        currNsiteArray = zeros(length(ts),'uint16'); %record number of sites

        parfor i = 1:length(ts)
            cc = ids(i); %current cluster index
            currT = t_all + double(ts(i))/fs;
            currYPos = calcYPos_v2( currT, p_uY.Value(cc), drift );
            [currWav, uSites] = ...
                intWav( p_uF.Value{cc}, p_uX.Value(cc), currYPos, p_uR.Value(cc), p_xcoords.Value, p_ycoords.Value, p_connected.Value, nt );
            nSite = length(uSites);
            currWavArray(i,:,:) = padarray(currWav * am(i),[(Nchan-nSite),0],0,'post');
            currSiteArray(i, :) = padarray(uSites,(Nchan-nSite),0,'post')';
            currNsiteArray(i) = nSite;
        end

    end

    for i = 1:length(ts)
        %for a given time, need to calcuate the wave that correspond to
        %the current position of the probe.
        allspks = allspks + 1;
        currT = t_all + double(ts(i))/fs;
        cc = ids(i); %current cluster index

        %get current position of this unit
        currYPos = calcYPos_v2( currT, uY(cc), drift );

        %currYdrift = calcYPos( currT, ycoords );
        yDriftRec(allspks,1) = currT;
        yDriftRec(allspks,2) = currYPos;
        yDriftRec(allspks,3) = cc; %record the unit label in the drift record

        if (useParPool)
            %get the waves for this unit from the big precalculated array
            uSites = squeeze(currSiteArray(i,1:currNsiteArray(i)));
            tempWav = squeeze(currWavArray(i,1:currNsiteArray(i),:));
            dat(ts(i) + tRange, uSites) = dat(ts(i) + tRange, uSites) + tempWav';
        else
            %calculate the interpolations now
            [tempWav, uSites] = ...
                intWav( uF{cc}, uX(cc), currYPos, uR(cc), xcoords, ycoords, connected, nt );

            dat(ts(i) + tRange, uSites) = dat(ts(i) + tRange, uSites) +...
                am(i) * tempWav(:,:)';
        end

        cM = find(uSites == monSite(cc));
        %if drift has moved the monitor site outside the footprint of the
        %unit, the recorded amplitudes just stay = 0;

        if ( cM )
            yDriftRec(allspks,4) = max(tempWav(cM,:)) - min(tempWav(cM,:));
            yDriftRec(allspks,5) = max(dat(ts(i) + tRange,monSite(cc))) - min(dat(ts(i) + tRange,monSite(cc)));
        end

    end

    %zero out the unconnected channels

    dat(:, zeroSites') = 0; % these are the reference and dig channels

    dat_old    =  dat;
    %convert to 16 bit integers; waveforms are in uV
    dat = int16(bitPerUV * dat);
    if addSYNC
        %add the column of sync data
        dat = horzcat(dat, sync);
    end
    fwrite(fidW, dat(1:(NT-buff),:)', 'int16');

    t_all = t_all + (NT-buff)/fs;
    elapsedTime = toc;
    fprintf( 'created %.2f seconds of data; nSpikes %d; calcTime: %.3f\n ', t_all, length(ts), elapsedTime );
    enoise_old = enoise;
    clear currWavArray;
    clear currNSiteArray;
    clear currSiteArray;
end

fclose(fidW); % all done

gtRes = spk_times + nt/2; % add back the time of the peak for the templates (half the time span of the defined waveforms)
gtClu = clu;

save(fullfile(fpath, 'eMouseGroundTruth'), 'gtRes', 'gtClu')
save(fullfile(fpath, 'eMouseSimRecord'), 'yDriftRec' );

end

end

function [uF, uR] = makeScatInt( uColl, nt )
    %create an interpolating function for each unit in
    % the array uColl
    uF = {};
    % array of maximum radii for each unit. only sites within this radius
    % will get get contributions from this unit
    uR = [];
    for i = 1:length(uColl)
        %reshape waves to a 1D vector
        wave_pts = reshape(uColl(i).waves',[uColl(i).nChan*nt,1]);
        xpts = repelem(uColl(i).xC,nt)';
        ypts = repelem(uColl(i).yC,nt)';
        tpts = (repmat(1:nt,1,uColl(i).nChan))';
        %specify points as y, x to be congruent with row,column of grid
        %interpolant.
        uF{i} = scatteredInterpolant( ypts, xpts, tpts, wave_pts, 'linear', 'nearest' );
        uR(i).minX = min(uColl(i).xC);
        uR(i).maxX = max(uColl(i).xC);
        uR(i).minY = min(uColl(i).yC);
        uR(i).maxY = max(uColl(i).yC);
        %earlier version that used y extend
        %uR(i) = min( max(uColl(i).yC), abs(min(uColl(i).yC)) );
    end

end

function [uF, uR] = makeGridInt( uColl, nt )
    %Reshape the array of [nsite, nt] into [Y, X, nt]
    %Note that this is dependent on how the data were stored:
    %here, assumes the sites are stored in row major order.
    %Could also derive this from the X and Y coordinates for generality

    %create an interpolating function for each unit in
    % the array uColl
    uF = {};
    % array of maximum radii for each unit. only sites within this radius
    % will get get contributions from this unit
    uR = [];

    for i = 1:length(uColl)
        sCol = length(unique(uColl(i).xC));
        sRow = length(unique(uColl(i).yC));
        %reshape waveform array into [sRow x sCol x nt array]
        %the data are stored as (1,1),(1,2),(1,3)...(2,1),(2,3)
        %reshape transforms as "row fast", so need to reshape and permute
        v = permute(reshape(uColl(i).waves,[sCol,sRow,nt]),[2,1,3]);
        xVal = uColl(i).xC(1:sCol);
        %pick off first element of each column to get y values
        colOne = (0:sRow-1)*sCol + 1;
        yVal = uColl(i).yC(colOne);
        tVal = 1:nt;
        %remember: y = rows, x = columns!
        [Y,X,T] = ndgrid(yVal, xVal, tVal);
        uF{i} = griddedInterpolant(Y,X,T,v,'makima');
        uR(i).minX = min(uColl(i).xC);
        uR(i).maxX = max(uColl(i).xC);
        uR(i).minY = min(uColl(i).yC);
        uR(i).maxY = max(uColl(i).yC);
        %earlier version that just used y extent
        %uR(i) = min( max(uColl(i).yC), abs(min(uColl(i).yC)) );
    end

end

function uSites = findSites( xPos, yPos, xcoords, ycoords, connected, uR )

    %find the sites currently in range for this unit at this point in time
    %(i.e. this value of yDrift)
    % xPos, yPos are the coordinates of the com of the unit signal in the
    % coordinates of the probe at the current time
    % xcoords and ycoords are the positions of the sites, assumed constant
    % motion of the probe should be modeled as rigid motion of the units

    %calculate distance from xPos, yPos for each site
    xDist = xcoords - xPos;
    yDist = ycoords - yPos;

    uSites = find( (xDist >= uR.minX) & (xDist <= uR.maxX) & ...
                   (yDist >= uR.minY) & (yDist <= uR.maxY) & ...
           	       (connected==1) );
%     dist_sq = (xPos - xcoords).^2 + (yPos - (ycoords + yDrift)).^2;
%     % want to exclude sites that aren't connected; just add a constant so
%     % they won't pass the distance test
%     dist_sq(find(connected==0)) = dist_sq(find(connected==0)) + 10*maxRad^2;
%     uSites = find( dist_sq < maxRad^2 );

end

function [uWav, uSites] = intWav( currF, xPos, yPos, uR, xcoords, ycoords, connected, nt )

    % figure out for which sites we need to calculate the waveform
    uSites = findSites( xPos, yPos, xcoords, ycoords, connected, uR );

    % given an array of sites on the probe, calculate the waveform using
    % the interpolant determined for this unit
    % xPos and yPos are the positions of the current unit
    % nt = number of timepoints

    currX = xcoords - xPos;
    currY = ycoords - yPos;

    nSites = length(uSites);
    xq = double(repelem(currX(uSites), nt));
    yq = double(repelem(currY(uSites), nt));
    tq = (double(repmat(1:nt, 1, nSites )))';
    %remember, y = rows in the grid, and x = columns in the grid
    %interpolation, and scattered interpolation set to match.
    nVal = numel(currF.Values);
    %tic
    uWav = currF( yq, xq, tq );
    %fprintf( '%d\t%d\t%.3f\n', numel(uSites), nVal, 1000*toc);
    uWav = (reshape(uWav', [nt,nSites]))';

end
function [uWav, uSites] = dumWav( currF, xPos, yPos, uR, xcoords, ycoords, connected, nt )

    % figure out for which sites we need to calculate the waveform
    uSites = findSites( xPos, yPos, xcoords, ycoords, connected, uR );
    uWav = zeros(length(uSites),nt);
end

function monSite = pickMonSite( currWav, uSites, connected )

    % find the connected site with the largest expected signal with no
    % drift. If all uSites are unconnected, returns an error

    %calc amplitudes for each site
    currAmp = max(currWav,[],2) - min(currWav,[],2);

    %sort amplitudes
    [sortAmp, ind] = sort(currAmp,'descend');

    monSite = -1;
    i = 1;

    while ( i < length(uSites) && monSite < 0 )
        if (connected(uSites(ind(i))))
            monSite = uSites(ind(i));
        else
            i = i + 1;
        end
    end

end

function currYPos = calcYPos_v2( t, yPos0, drift  )
%   calculate current position of a unit given the current time and
%   initial position (yPos0)

%   The pattern of tissue motion in space is set by drift.sType:
%       'rigid' -- all the units move together
%       'point' -- motion is largest at a point y0 (furthest from tip) and
%         decreases exponentially for units far from y0. Need to
%         specify y0 and the halfDistance
%
%   The pattern of tissue motion in time is set by drift.tType:
%       'exp' -- initial fast transition followed by exponential
%            decay; specify stepsize (fast transition distance); halfLife (sec)
%            and period (in sec).
%       'sine' -- specify amplitude and period
%
%   Drift parameter values for 20 um p-p, uniform sine motion, with period = 300 s:
%     drift.sType = 'rigid';        %'rigid' or 'point';
%     drift.tType = 'sine';          %'exp' or 'sine'
%
%     drift.y0 = 3800;        %in um, position along probe where motion is largest
%                             %conventially, y = 0 is the bottom of the probe
%     drift.halfDistance = 1000;     %in um
%
%     drift.amplitude = 10;         %in um. For a sine wave, the peak to
%                                    peak variation is 2Xdrift.amplitude
%     drift.halfLife = 10;          %in seconds
%     drift.period = 300;           %in seconds
%
if (drift.addDrift)
    switch drift.tType
        case 'exp'
            timeIntoCycle = t - (floor(t/drift.period))*drift.period;
            delta = drift.amplitude*exp(-timeIntoCycle/drift.halfLife);
        case 'sine'
            delta = drift.amplitude*sin(t*2*pi()/drift.period);
        otherwise
            fprintf( 'unknown parameter in drift calculation \n')
            return;
    end

    switch drift.sType
        case 'rigid'
            %delta is equal for any position
        case 'point'
            %delta falls off exponentially from y0
            delta = delta * exp( -abs(drift.y0 - yPos0)/drift.halfDistance );
        otherwise
            fprintf( 'unknown parameter in drift calculation \n')
            return;
    end
    currYPos = yPos0 + delta;
else
    currYPos = yPos0;
end

end

function eNoise = makeNoise( noiseSamp,noiseModel,chanMap,connected,NchanTOT )

    %if chanMap is a short version of a 3A probe, use the first
    %nChan good channels to generate noise, then copy that array
    %into an NT X NChanTot array

    nChan = numel(chanMap);        %number of noise channels to generate
    goodChan = sum(connected);
    tempNoise = zeros( noiseSamp, goodChan, 'single' );
    nT_fft = noiseModel.nm.nt;         %number of time points in the original time series
    fftSamp = noiseModel.nm.fft;

    noiseBatch = ceil(noiseSamp/nT_fft);
    lastWind = noiseSamp - (noiseBatch-1)*nT_fft; %in samples

    for j = 1:goodChan
            for i = 1:noiseBatch-1
                tStart = (i-1)*nT_fft+1;
                tEnd = i * nT_fft;
                tempNoise(tStart:tEnd,j) = fftnoise(fftSamp(:,j),1);
            end
            %for last batch, call one more time and truncate
            lastBatch = fftnoise(fftSamp(:,j),1);
            tStart = (noiseBatch-1)*nT_fft+1;
            tEnd = noiseSamp;
            tempNoise(tStart:tEnd,j) = lastBatch(1:lastWind);
    end

    %unwhiten this array
    Wrot = noiseModel.nm.Wrot(1:goodChan,1:goodChan);
    tempNoise_unwh = tempNoise/Wrot;

    %scale to uV; will get scaled back to bits at the end
    tempNoise_unwh = tempNoise_unwh/noiseModel.nm.bitPerUV;

    %to get the final noise array, map to an array including all channels
    eNoise = zeros(noiseSamp, NchanTOT, 'single');
    %indicies of the good channels
    goodChanIndex = find(connected);
    eNoise(:,chanMap(goodChanIndex)) = tempNoise_unwh;

end

function noise=fftnoise(f,Nseries)
% Generate noise with a given power spectrum.
% Useful helper function for Monte Carlo null-hypothesis tests and confidence interval estimation.
%
% noise=fftnoise(f[,Nseries])
%
% INPUTS:
% f: the fft of a time series (must be a column vector)
% Nseries: number of noise series to generate. (default=1)
%
% OUTPUT:
% noise: surrogate series with same power spectrum as f. (each column is a surrogate).
%
%   --- Aslak Grinsted (2009)
%
if nargin<2
    Nseries=1;
end
f=f(:);     %ensures f is a column vector
N=length(f);
Np=floor((N-1)/2);
phases=rand(Np,Nseries)*2*pi;
phases=complex(cos(phases),sin(phases)); % this was the fastest alternative in my tests.
f=repmat(f,1,Nseries);
f(2:Np+1,:)=f(2:Np+1,:).*phases;
f(end:-1:end-Np+1,:)=conj(f(2:Np+1,:));
noise=real(ifft(f,[],1));

end

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
%%*******************************************************************
%% schurmat_qblk: compute schur matrix corresponding to SOCP blocks.
%%
%% HKM direction: output = schur + Ax*Ae' + Ae*Ax' - Ad*Ad'
%% NT  direction: output = schur + Ae*Ae' - Ad*Ad'
%%
%% where schur = A*D*A', and Ad is the modification to ADA'
%% so that the latter is positive definite.
%%
%% [schur,UU,EE] = schurmat_qblk(blk,At,schur,UU,EE,p,dd,ee,xx);
%%
%% UU: stores the dense columns of Ax, Ae, Ad, and possibly
%%     those of A*D^{1/2}. It has the form UU = [Ax Ae Ad].
%% EE: stores the assocaited (2,2) block matrix when the
%%     output matrix is expressed as an augmented matrix.
%%     It has the form EE = [0 -lam 0; -lam 0 0; 0 0 I].
%%
%% options = 0, HKM
%%         = 1, NT
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function [schur,UU,EE] = schurmat_qblk(blk,At,par,schur,UU,EE,p,dd,ee,xx);

   global idxdenAq nnzschur_qblk

   if (nargin == 10); options = 0; else; options = 1; end;
   iter = par.iter;

   if isempty(EE)
      count = 0;
   else
      count = max(max(EE(:,2)),max(EE(:,1)));
   end
   pblk = blk(p,:); n = sum(pblk{2}); numblk = length(pblk{2});
%%
   Ae = qprod(pblk,At{p}',ee{p});
   if (options == 0)
      Ax = qprod(pblk,At{p}',xx{p});
   end
   idxden = checkdense(Ae);
   ddsch = dd{p};
   if ~isempty(idxden);
      spcolidx = setdiff([1:numblk],idxden);
      s = 1 + [0, cumsum(pblk{2})];
      idx = s(idxden);
      tmp = zeros(n,1);
      tmp(idx) = sqrt(2*abs(ddsch(idx)));
      Ad = qprod(pblk,At{p}',tmp);
      ddsch(idx) = abs(ddsch(idx));
      if (options == 0)
         len = length(idxden);
         gamzsub = par.gamz{p}(idxden);
         lam = gamzsub.*gamzsub;
         UU = [UU, Ax(:,idxden), Ae(:,idxden)*spdiags(lam,0,len,len), Ad(:,idxden)];
         tmp = count+[1:len]';
         EE = [EE; [tmp, len+tmp, -lam; len+tmp, tmp, -lam; ...
                    2*len+tmp, 2*len+tmp, ones(len,1)] ];
         count = count+3*len;
         Ax = Ax(:,spcolidx); Ae = Ae(:,spcolidx);
         tmp = Ax*Ae';
         schur = schur + (tmp + tmp');
      else
         len = length(idxden);
         w2 = par.gamz{p}./par.gamx{p};
         lam = w2(idxden);
         UU = [UU, Ae(:,idxden)*spdiags(sqrt(lam),0,len,len), Ad(:,idxden)];
         tmp = count+[1:len]';
         EE = [EE; [tmp, tmp, -lam; len+tmp, len+tmp, ones(len,1)] ];
         count = count + 2*len;
         Ae = Ae(:,spcolidx);
         schur = schur + Ae*Ae';
      end
   else
      if (options == 0)
         tmp = Ax*Ae';
         schur = schur + (tmp+tmp');
      else
         tmp = Ae*Ae';
         schur = schur + tmp;
      end
   end
   if (iter==1)
      idxdenAq{p} = checkdense(At{p}');
   end
   if ~isempty(idxdenAq{p});
      idxden = idxdenAq{p};
      len = length(idxden);
      Ad = At{p}(idxden,:)'*spdiags(sqrt(abs(ddsch(idxden))),0,len,len);
      UU = [UU, Ad];
      tmp = count+[1:len]';
      EE = [EE; [tmp, tmp, -sign(ddsch(idxden))]];
      count = count + len;
      ddsch(idxden) = zeros(len,1);
   end
   schurtmp = At{p}' *spdiags(ddsch,0,n,n) *At{p};
   schur = schur + schurtmp;
%%*******************************************************************


text:
function fiber_curvature_vals=dtiFiberCurvature(fiber)

%For 3xn fiber computes curvature at each point as defined here
%http://en.wikipedia.org/wiki/Curvature (I know it is a shame to use Wiki
%for research -- but oh well
%
%ER 12/2007
xprime=gradient(fiber(1, :));
yprime=gradient(fiber(2, :));
zprime=gradient(fiber(3, :));

fiber_curvature_vals=sqrt((gradient(zprime).*yprime-gradient(yprime).*zprime).^2+(gradient(xprime).*zprime-gradient(zprime).*xprime).^2+(gradient(yprime).*xprime-gradient(xprime).*yprime).^2)./((xprime.^2 + yprime.^2 + zprime.^2).^(3/2));
%fiber_curvature_vals=fiber_curvature_vals./max(fiber_curvature_vals(:));

% These values are sometimes huge. what is the distribution? The original
% data are not unitless.

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function amd_install
%AMD_INSTALL compile and install amd2 for use in MATLAB
%   Your current directory must be AMD/MATLAB for this function to work.
%
% Example:
%   amd_install
%
% See also amd, amd2.

% Copyright 1994-2007, Tim Davis, University of Florida,
% Patrick R. Amestoy, and Iain S. Duff.

% This orders the same matrix as the ANSI C demo, amd_demo.c.  It includes an

amd_make
addpath (pwd)
fprintf ('\nThe following path has been added.  You may wish to add it\n') ;
fprintf ('permanently, using the MATLAB pathtool command.\n') ;
fprintf ('%s\n\n', pwd) ;
amd_demo

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

% Script chap::2::script
% Density Heston Model
%
%   uses the characteristic function of the Heston model
%
t = 10;                             % maturity
a = 600;                            % spot value
N = 512;                            % number of grid points
x = ( (0:N-1) - N/2 ) / a;          % range

figname = 'Risk Neutral Density - Heston-Hull-White';

vInst_base = 0.02;                  % instantanuous variance of base parameter set
vLong_base = 0.02;                  % long term variance of base parameter set
kappa_base = 0.1;                   % mean reversion speed of variance of base parameter set
omega_base = 0.2;                   % volatility of variance of base parameter set
rho_base = 0;                       % correlation of base parameter set

lambda_base = 0.1;
eta_base = 0.02;

icurveData = [0.999884333380315;0.996803132736937;0.993568709230647;0.990285301195274;0.986945903402709;0.983557350486521;0.980185549124449;0.976782934344041;0.973361992614499;0.969976793305220;0.966616749933289;0.962914317958160;0.959904777446077;0.920091903961326;0.882870065420196;0.847186544281939;0.812742515687365;0.779459552415061;0.747152463119429;0.715745016074346;0.685138723808460;0.655753392359115;0.627333845297308;0.599226698198774;0.572763319281569;0.547259133751455;0.523441996253080;0.499646068368557;0.477507905873099;0.456481811728753;0.436385788738282;0.417350253831050;0.399187111819286;0.381865611666566;0.365435617455498;0.349786183601181;0.334806921914717;0.320548897004994;0.306983265264429;0.294081800917050;0.282443547729164;0.269929224010243];
icurveDates = [734472;734501;734534;734562;734591;734622;734653;734683;734713;734744;734775;734807;734836;735202;735567;735931;736298;736663;737028;737393;737758;738125;738489;738854;739219;739585;739949;740316;740680;741046;741411;741776;742140;742507;742872;743237;743602;743967;744334;744698;745063;745429];
icurveInterpMethod = 'spline';
icurveType = 'Discount';
icurveSettle = 734471;
irdc = IRDataCurve(icurveType,icurveSettle,icurveDates,icurveData);
% cf (characteristic function) for the heston hull white model
f = @(x) cf_heston(x,vInst_base,vLong_base,kappa_base,omega_base,rho_base,t,0) .* cf_hullwhite(x+1i,t,0,lambda_base,eta_base,irdc);

legendname_base = 'Base parameter set';
y_base = fftdensity(f,a,N);         % density calculated from cf for base parameter set
%% Changing lambda_base
lambda_low = .005;                    % changing parameter (low value)
lambda_high = .5;                   % changing parameter (high value)

f_low = @(x) cf_heston(x,vInst_base,vLong_base,kappa_base,omega_base,rho_base,t,0).* cf_hullwhite(x+1i,t,0,lambda_low,eta_base,irdc);
f_high = @(x)  cf_heston(x,vInst_base,vLong_base,kappa_base,omega_base,rho_base,t,0).* cf_hullwhite(x+1i,t,0,lambda_high,eta_base,irdc);

y_low = fftdensity(f_low,a,N);
y_high = fftdensity(f_high,a,N);

legendname_low = 'Changing \lambda low value';
legendname_high = 'Changing \lambda high value';

% output density as figure
createfigure_density(x,y_base,y_low,y_high,...
    figname, legendname_base,legendname_low,legendname_high);
%% Changing eta_base
eta_low = .01;
eta_high = .03;

f_low =  @(x) cf_heston(x,vInst_base,vLong_base,kappa_base,omega_base,rho_base,t,0).* cf_hullwhite(x+1i,t,0,lambda_base,eta_low,irdc);
f_high =  @(x) cf_heston(x,vInst_base,vLong_base,kappa_base,omega_base,rho_base,t,0).* cf_hullwhite(x+1i,t,0,lambda_base,eta_high,irdc);

y_low = fftdensity(f_low,a,N);
y_high = fftdensity(f_high,a,N);

legendname_low = 'Changing \eta low value';
legendname_high = 'Changing \eta high value';

createfigure_density(x,y_base,y_low,y_high,...
    figname, legendname_base,legendname_low,legendname_high);

%% Changing curve_base
zero = -log(icurveData)./(icurveDates-icurveSettle)*360 - 0.01;
icurveData_low = exp(-zero .*(icurveDates-icurveSettle)/360);
zero = -log(icurveData)./(icurveDates-icurveSettle)*360 + 0.01;
icurveData_high = exp(-zero .*(icurveDates-icurveSettle)/360);

irdc_low = IRDataCurve(icurveType,icurveSettle,icurveDates,icurveData_low);
irdc_high = IRDataCurve(icurveType,icurveSettle,icurveDates,icurveData_high);

f_low =  @(x) cf_heston(x,vInst_base,vLong_base,kappa_base,omega_base,rho_base,t,0).* cf_hullwhite(x+1i,t,0,lambda_base,eta_base,irdc_low);
f_high =  @(x) cf_heston(x,vInst_base,vLong_base,kappa_base,omega_base,rho_base,t,0).* cf_hullwhite(x+1i,t,0,lambda_base,eta_base,irdc_high);

y_low = fftdensity(f_low,a,N);
y_high = fftdensity(f_high,a,N);

legendname_low = 'Changing curve parallel shift - low value';
legendname_high = 'Changing curve parallel shift - high value';

createfigure_density(x,y_base,y_low,y_high,...
    figname, legendname_base,legendname_low,legendname_high);


text:
function [ mB ] = ApplyGaussianBlur( mI, gaussianKernelStd, stdToRadiusFactor )

gaussianBlurRadius  = ceil(stdToRadiusFactor * gaussianKernelStd);

vGaussianKernel = exp(-((-gaussianBlurRadius:gaussianBlurRadius) .^ 2) / (2 * gaussianKernelStd * gaussianKernelStd));
vGaussianKernel = vGaussianKernel / sum(vGaussianKernel);

mI   = PadArrayReplicate(mI, gaussianBlurRadius);

mB = conv2(vGaussianKernel, vGaussianKernel.', mI, 'valid');

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function e = end(X,k,n)
%END Last index of indexing expression for sparse tensor.
%
%   The expression X(end,:,:) will call END(X,1,3) to determine
%   the value of the first index.
%
%   See also SPTENSOR, SPTENSOR/SUBSREF, END.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

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

if n > ndims(X)
  error('Subscript out of range.');
end
e = X.size(k);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%%
% test for linear programming and interpolation of discrete meaures.

addpath('../toolbox/');
addpath('../toolbox/mexEMD/');

test = 'weighted';
test = 'empirical';

rep = ['results/' test '/'];
[~,~] = mkdir(rep);

% helpers
SetAR = @(ar)set(gca, 'PlotBoxAspectRatio', [1 ar 1]);
myplot = @(x,y,ms,col)plot(x,y, 'o', 'MarkerSize', ms, 'MarkerEdgeColor', col, 'MarkerFaceColor', col, 'LineWidth', 2);
myplot = @(x,y,ms,col)plot(x,y, 'o', 'MarkerSize', ms, 'MarkerEdgeColor', 'k', 'MarkerFaceColor', col, 'LineWidth', 1);

% Dimensions  of the clouds.
n0 = 4000;
n1 = n0;

% Compute a first point cloud \(X_0\) that is Gaussian.
% and a second point cloud \(X_1\) that is Gaussian mixture.
randn('state', 666);
gauss = @(q,a,c)a*randn(2,q)+repmat(c(:), [1 q]);
X0 = randn(2,n0)*.3;
X1 = [gauss(n1/2,.5, [0 1.6]) gauss(n1/4,.3, [-1 -1]) gauss(n1/4,.3, [1 -1])];
% weights
normalize = @(a)a/sum(a(:));
switch test
    case 'weighted'
        p0 = normalize(rand(n0,1));
        p1 = normalize(rand(n1,1));
    case 'empirical'
        p0 = 2*ones(n0,1)/n0;
        p1 = 2*ones(n1,1)/n1;
end

%%
% Display the point clouds.
% The size of each dot is proportional to its probability density weight.

clf; hold on;
for i=1:length(p0)
    myplot(X0(1,i), X0(2,i), p0(i)*length(p0)*10, 'b');
end
for i=1:length(p1)
    myplot(X1(1,i), X1(2,i), p1(i)*length(p1)*10, 'r');
end
axis([min(X1(1,:)) max(X1(1,:)) min(X1(2,:)) max(X1(2,:))]); axis off;

%%
% Compute the cost matrix

C = repmat( sum(X0.^2)', [1 n1] ) + ...
    repmat( sum(X1.^2), [n0 1] ) - 2*X0'*X1;

%%
% Solve the linprog of OT

[cost,gamma] = mexEMD(p0,p1,C);

%%
% Compute displacement interpolation.

[I,J,gammaij] = find(gamma);
tlist = linspace(0,1,6);
clf;
for k=1:length(tlist)
    t=tlist(k);
    Xt = (1-t)*X0(:,I) + t*X1(:,J);
    % subplot(2,3,i);
    clf;
    hold on;
    for i=1:length(gammaij)
        myplot(Xt(1,i), Xt(2,i), gammaij(i)*length(gammaij)*6, [t 0 1-t]);
    end
    % title(['t=' num2str(t,2)]);
    axis([min(X1(1,:)) max(X1(1,:)) min(X1(2,:)) max(X1(2,:))]);
    % dummy points
    plot([min(X1(1,:)) max(X1(1,:))], [min(X1(2,:)) max(X1(2,:))], '.', 'MarkerFaceColor', 'w', 'MarkerEdgeColor', 'w');
    axis equal; axis square;
    axis off;
    drawnow;
    saveas(gcf, [rep 'interp-' num2str(k) '.eps'], 'epsc');
end

author: optimaltransport


text:
clear; clc; close all;

DEVICE_ID = 0; %set gpu id starting from 0

% set paths to Caffe and DeepFool
PATH_CAFFE = '/path/to/caffe/matlab';
PATH_DEEPFOOL = '/path/to/DeepFool';
PATH_IMAGENET_TRAIN = '/path/to/ILSVRC2012/train';

addpath(PATH_CAFFE);
addpath(PATH_DEEPFOOL);

caffe.set_mode_gpu();
caffe.set_device(DEVICE_ID);

model = 'caffenet';

if (strcmp(model, 'vgg_16'))
    fprintf('Loading VGG 16\n');
    net_model_path = fullfile('data', 'deploy_vgg_16.prototxt');
    net_weights_path = fullfile('data', 'vgg_16.caffemodel');
    net_url = 'http://www.robots.ox.ac.uk/~vgg/software/very_deep/caffe/VGG_ILSVRC_16_layers.caffemodel';
elseif (strcmp(model, 'vgg_19'))
    fprintf('Loading VGG 19\n');
    net_model_path = fullfile('data', 'deploy_vgg_19.prototxt');
    net_weights_path = fullfile('data', 'vgg_19.caffemodel');
    net_url = 'http://www.robots.ox.ac.uk/~vgg/software/very_deep/caffe/VGG_ILSVRC_19_layers.caffemodel';
elseif (strcmp(model, 'googlenet'))
    fprintf('Loading GoogLeNet\n');
    net_model_path = fullfile('data', 'deploy_googlenet.prototxt');
    net_weights_path = fullfile('data', 'googlenet.caffemodel');
    net_url = 'http://dl.caffe.berkeleyvision.org/bvlc_googlenet.caffemodel';
elseif (strcmp(model, 'vgg_f'))
    fprintf('Loading VGG-F\n');
    net_model_path = fullfile('data', 'deploy_vgg_f.prototxt');
    net_weights_path = fullfile('data', 'googlenet.caffemodel');
    net_url = 'http://dl.caffe.berkeleyvision.org/bvlc_vgg_f.caffemodel';
elseif (strcmp(model, 'caffenet'))
    fprintf('Loading CaffeNet\n');
    net_model_path = fullfile('data', 'deploy_caffenet.prototxt');
    net_weights_path = fullfile('data', 'caffenet.caffemodel');
    net_url = 'http://dl.caffe.berkeleyvision.org/bvlc_reference_caffenet.caffemodel';
elseif (strcmp(model, 'resnet-152'))
    fprintf('Loading ResNet-152\n');
    net_model_path = fullfile('data', 'deploy_resnet.prototxt');
    net_weights_path = fullfile('data', 'caffenet.caffemodel');
    net_url = 'https://deepdetect.com/models/resnet/ResNet-152-model.caffemodel';
else
    error('Model is not recognized!');
end

if (~exist(net_weights_path, 'file'))
    fprintf('Downloading model...\n');
    websave(net_weights_path, net_url);
end

net = caffe.Net(net_model_path, net_weights_path, 'test'); % run with phase test (so that dropout isn't applied)
fprintf('Network is loaded\n');

% Loading the data
fprintf('Loading the data...\n');
im_array = makeImagenetData(PATH_IMAGENET_TRAIN);
fprintf('Data is loaded\n');

% set options
opts.library = 'caffe';
opts.net = net;

v = universal_perturbation(im_array, opts);

% Test perturbation on a sample image

d = load(fullfile('data', 'ilsvrc_2012_mean.mat'));
mean_data = d.mean_data;

labels_file = fullfile('data', 'synset_words.txt');
fileID = fopen(labels_file);
C = textscan(fileID, '%c %d %s', 'Delimiter', '\n');
synset_names = C{3};
synset_names_short = synset_names;
for i = 1:1000
    synset_names_short(i) = strtok(synset_names(i), ',');
end

im_test = imread(fullfile('data', 'test_img.png'));
im_test = preprocess_img(im_test, mean_data);

original_label = predict_caffe(im_test, net, 1);
perturbed_label = predict_caffe(im_test+v, net, 1);

% Show original and perturbed images
im_ = im_test + mean_data(1:224, 1:224, :);
im_ = permute(im_,[2,1,3]);
im_ = im_(:,:,[3,2,1],:);

figure; imshow(uint8(im_));
xlabel(synset_names_short(original_label));
set(gca,'xtick',[],'ytick',[])
set(gcf, 'Color', 'white');
set(gca, 'FontSize', 15);
title('Original image');

v_ = permute(v,[2,1,3]);
v_ = v_(:,:,[3,2,1],:);
figure; imshow(uint8(im_ + v_));
xlabel(synset_names_short(perturbed_label));
set(gca,'xtick',[],'ytick',[])
set(gcf, 'Color', 'white');
set(gca, 'FontSize', 15);
title('Perturbed image');

repo: universal