text:
function [zPred, PzPred,otherInfo]=cubKalMeasPred(xPred,PPred,zDim,h,xi,w,innovTrans,measAvgFun,stateDiffTrans,stateTrans)
%%CUBKALMEASPRED Perform the measurement prediction part of the measurement
%           update step of the cubature Kalman filter with additive
%           measurement noise. The function cubKalUpdateWithPred can be
%           used to complete the measurement update. Separating the
%           measurement prediction step from the rest of the update step
%           can make the creation of multiple measurement association
%           hypotheses from a single target prediction more efficient.
%           The full measurement update function is cubKalUpdate.
%
%INPUTS: xPred The xDimXnumComp predicted target states.
%        PPred The xDimXxDimXnumComp predicted state covariance matrices.
%         zDim The dimensionality of the output of h.
%            h A function handle for the measurement function that takes
%              the state as its argument.
%           xi An xDimXnumCubPoints matrix of cubature points. If this
%              and the next parameter are omitted or empty matrices are
%              passed, then fifthOrderCubPoints(xDim) is used. It is
%              suggested that xi and w be provided to avoid needless
%              recomputation of the cubature points.
%            w A numCubPointsX1 vector of the weights associated with the
%              cubature points.
%   innovTrans An optional function handle that computes and optionally
%              transforms the value of the difference between the
%              observation and any predicted points. This is called as
%              innovTrans(a,b) and the default if omitted or an empty
%              matrix is passed is @(a,b)bsxfun(@minus,a,b). This must be
%              able to handle sets of values. For a zDimX1 measurement,
%              either of the inputs could be zDimXN in size while one of
%              the inputs could be zDimX1 in size.  This only needs to be
%              supplied when a measurement difference must be restricted
%              to a certain range. For example, the innovation between two
%              angles will be 2*pi if one angle is zero and the other
%              2*pi, even though they are the same direction. In such an
%              instance, a function handle to the
%              wrapRange(bsxfun(@minus,a,b),-pi,pi) function with the
%              appropriate parameters should be passed for innovTrans.
%   measAvgFun An optional function handle that, when given N measurement
%              values with weights, produces the weighted average. This
%              function only has to be provided if the domain of the
%              measurement is not linear. For example, when averaging
%              angular values, then the function meanAng should be used.
% stateDiffTrans An optional function handle that takes an xDimXN matrix of
%              N differences between states and transforms them however
%              might be necessary. If not transformation is necessary, this
%              parameter can be omitted or an empty matrix passed.
%   stateTrans An optional function that takes a state estimate and
%              transforms it. This is useful if one wishes the elements of
%              the state to be bound to a certain domain. For example, if
%              an element of the state is an angle, one might generally
%              want to bind it to the region +/-pi.
%
%OUTPUTS: zPred The zDimXnumComp measurement predictions from the filter.
%        PzPred The zDimXzDimXnumComp covariance matrices associated with
%               zPred.
%     otherInfo A structure containing members of intermediate results of
%               this function that can be passed to cubKalUpdateWithPred
%               when updating with a measurement.
%
%The mathematics behind the cubature Kalman filter are described in more
%detail in Section IX of [1] and in [2]. See the comments to cubKalUpdate
%for more information.
%
%EXAMPLE:
%With this example, we demonstrate that one gets the same result using
%cubKalUpdate in one step as with using cubKalMeasPred followed by
%cubKalUpdateWithPred.
% xPred=[1e3;-2e3;100;200];
% h=@(x)([sum(x.^2);x(1)-x(4)^(3/2)]);
% PPred=[28,   3.5,    6,  8.5;
%      3.5,    23,  8.5,   11;
%        6,   8.5,   18, 13.5;
%      8.5,    11, 13.5,   13];
% z=1e6*[5.050000548964568;
%       -0.001829553054023];
% zDim=size(z,1);
% R=eye(zDim,zDim);%Measurement covariance matrix.
% %The update in one step.
% [xUpdate,PUpdate,innov,Pzz,W]=cubKalUpdate(xPred,PPred,z,R,h);
% %The update in two steps.
% [zPred, PzPred,otherInfo]=cubKalMeasPred(xPred,PPred,zDim,h);
% [xUpdate1,PUpdate1,innov1,Pzz1,W1]=cubKalUpdateWithPred(z,R,zPred,PzPred,otherInfo);
% %One will see that the one and two step updates agree.
% max(abs([xUpdate1-xUpdate;PUpdate1(:)-PUpdate(:);innov1(:)-innov;Pzz1(:)-Pzz(:);W1(:)-W(:)]))
%
%REFERENCES:
%[1] D. F. Crouse , "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems Magazine,
%    vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%[2] I. Arasaratnam and S. Haykin, "Cubature Kalman filters," IEEE
%    Transactions on Automatic Control, vol. 54, no. 6, pp. 1254-1269,
%    Jun. 2009.
%
%June 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    xDim=size(xPred,1);
    numComp=size(xPred,2);

    if(nargin<5||isempty(xi))
        [xi,w]=fifthOrderCubPoints(xDim);
    end

    if(nargin<7||isempty(innovTrans))
        %The function just returns the input.
        innovTrans=@(a,b)bsxfun(@minus,a,b);
    end

    if(nargin<8||isempty(measAvgFun))
        measAvgFun=@(zPoints,w)calcMixtureMoments(zPoints,w);
    end

    if(nargin<9||isempty(stateDiffTrans))
        stateDiffTrans=@(x)x;
    end

    if(nargin<10||isempty(stateTrans))
        stateTrans=@(x)x;
    end

    numCubPoints=size(xi,2);

    zPred=zeros(zDim,numComp);
    PzPred=zeros(zDim,zDim,numComp);
    xPredCenPoints=zeros(xDim,numCubPoints,numComp);
    zPredCenPoints=zeros(zDim,numCubPoints,numComp);
    Pxz=zeros(xDim,zDim,numComp);

    for k=1:numComp
        %cholSemiDef is used instead of chol in case a positive semi-definite
        %covariance matrix is passed.
        SPred=cholSemiDef(PPred,'lower');
        %Predicted cubature state points
        xPredPoints=stateTrans(transformCubPoints(xi,xPred,SPred));

        %Predicted cubature measurement points
        zPredPoints=zeros(zDim,numCubPoints);
        for curP=1:numCubPoints
            zPredPoints(:,curP)=h(xPredPoints(:,curP));
        end

        %Measurement prediction.
        zPred=measAvgFun(zPredPoints,w);

        %Centered, predicted cubature measurement points, transformed as
        %necessary to keep the values within a desired range.
        zPredCenPoints(:,:,k)=innovTrans(zPredPoints,zPred);
        xPredCenPoints(:,:,k)=stateDiffTrans(bsxfun(@minus,xPredPoints,xPred));
        for curP=1:numCubPoints
            diff=zPredCenPoints(:,curP);
            PzPred(:,:,k)=PzPred(:,:,k)+w(curP)*(diff*diff');
            Pxz(:,:,k)=Pxz(:,:,k)+w(curP)*xPredCenPoints(:,curP)*diff';
        end
        %Pxz is not needed for the measurement prediction, but we compute it
        %here, so that it need not be recomputed again and again if
        %cubKalUpdateWithPred is called for multiple measurements.
    end

    otherInfo.innovTrans=innovTrans;
    otherInfo.stateDiffTrans=stateDiffTrans;
    otherInfo.stateTrans=stateTrans;
    otherInfo.xPredCenPoints=xPredCenPoints;
    otherInfo.zPredCenPoints=zPredCenPoints;
    otherInfo.xPred=xPred;
    otherInfo.w=w;
    otherInfo.Pxz=Pxz;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function y=sqr(x)
y = x*x;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function b = isdir_private(dirname)
%
% isdir_private() is a backward compatible version of matlab's isdir()
% function.
%
% isdir() is a new matlab function that is an improvment over exist() to
% tell if a pathname is a directory or not. But it didn't exist prior to R2017.
% Therefore we use try/catch to still be able to use isdir when it exists

try
    b = isdir(dirname);
catch
    b = (exist(dirname,'dir') == 7);
end

author: BUNPC


text:
function disp(obj)

% function disp(obj)
% Disp a configuration dependency
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: disp.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

sz = size(obj);
fprintf('%s object: ', class(obj));
if length(sz)>4,
    fprintf('%d-D\n',length(sz));
else
    for i=1:(length(sz)-1),
        fprintf('%d-by-',sz(i));
    end;
    fprintf('%d\n',sz(end));
end;
if prod(sz)==1,
    so = struct(obj);
    disp(so);
end;
return;

author: spm


text:
function [imo,rois,labels,targets] = fast_rcnn_train_get_batch(images, imdb, batch, opts)
% FAST_RCNN_GET_BATCH_TRAIN  Generates mini-batches for Fast-RCNN train

% opts.numFgRoisPerImg = 128;
% opts.numRoisPerImg = 64;
% opts.maxScale = 1000;
% opts.bgLabel = 21;
% opts.visualize = 0;
% opts.scale = 600;
% opts.interpolation = 'bicubic';
% opts.averageImage = [];
% opts.numThreads = 2;
% opts.prefetch = true;
%
% Copyright (C) 2016 Hakan Bilen.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if isempty(images)
  imo = [] ;
  rois = [] ;
  labels = [] ;
  targets = [] ;
  return ;
end

% fetch is true if images is a list of filenames (instead of
% a cell array of images)
fetch = ischar(images{1}) ;

% prefetch is used to load images in a separate thread
prefetch = fetch & opts.prefetch ;

if prefetch
  vl_imreadjpeg(images, 'numThreads',opts.numThreads,'prefetch') ;
  imo = [] ;
  rois = [] ;
  labels = [] ;
  targets = [] ;
  return ;
end

if fetch
  ims = vl_imreadjpeg(images,'numThreads',opts.numThreads) ;
else
  ims = images ;
end

maxW = 0;
maxH = 0;

% labels = imdb.images.label(:,batch);

pboxes   = cell(1,numel(batch));
plabels  = cell(1,numel(batch));
ptargets = cell(1,numel(batch));

% get fg and bg rois
for b=1:numel(batch)
  pbox   = imdb.boxes.pbox{batch(b)};
  plabel = imdb.boxes.plabel{batch(b)};
  ptarget = imdb.boxes.ptarget{batch(b)};

  if size(pbox,2)~=4
    error('wrong box size');
  end

  % get pos boxes
  pos = find((plabel~=opts.bgLabel) & (plabel > 0)) ;
  npos = numel(pos);
  % get neg boxes
  neg = find((plabel==opts.bgLabel)) ;
  nneg = numel(neg);

    bbox = [];
    label = [];
    target = [];

    opts.numFgRoisPerImg = min(npos,opts.numFgRoisPerImg);
    nBneg = min(nneg,opts.numRoisPerImg-opts.numFgRoisPerImg);

    if npos>0
      r = randperm(npos);
      p = pos(r(1:opts.numFgRoisPerImg));
      bbox = pbox(p,:);
      label = plabel(p);
      target = ptarget(p,:);
    end
    if nneg>0
      r = randperm(nneg);

      n = neg(r(1:nBneg));
      bbox = [bbox ; pbox(n,:)];
      label = [label ; plabel(n)];
      target = [target ; ones(size(ptarget(n,:)))];
    end
  pboxes{b} = bbox;
  plabels{b} = label;
  ptargets{b} = target;
end

if isempty(pboxes)
  warning('No gt box\n');
end

labels = vertcat(plabels{:});
targets = vertcat(ptargets{:});

% rescale images and rois
rois = [];
imre = cell(1,numel(batch));
for b=1:numel(batch)
  imSize = size(ims{b});

  h = imSize(1);
  w = imSize(2);

  factor = max(opts.scale(1)/h,opts.scale(1)/w);

  if any([h*factor,w*factor]>opts.maxScale)
    factor = min(opts.maxScale/h,opts.maxScale/w);
  end

  if abs(factor-1)>1e-3
    imre{b} = imresize(ims{b},factor,'Method',opts.interpolation);
  else
    imre{b} = ims{b};
  end

  if imdb.boxes.flip(batch(b))
    im = imre{b};
    imre{b} = im(:,end:-1:1,:);
  end

  imreSize = size(imre{b});

  maxH = max(imreSize(1),maxH);
  maxW = max(imreSize(2),maxW);

  % adapt bounding boxes into new coord
  bbox = pboxes{b};
  if any(bbox(:)<=0)
    error('bbox error');
  end

  nB = size(bbox,1);
  tbbox = bbox_scale(bbox,factor,[imreSize(2) imreSize(1)]);
  if any(tbbox(:)<=0)
    error('tbbox error');
  end

  rois = [rois [b*ones(1,nB) ; tbbox' ] ];
end

imo = zeros(maxH,maxW,size(imre{1},3),numel(batch),'single');
for b=1:numel(batch)
  % subtract mean
  if ~isempty(opts.averageImage)
    imre{b} = single(bsxfun(@minus,imre{b},opts.averageImage));
  end
  sz = size(imre{b});
  imo(1:sz(1),1:sz(2),:,b) = single(imre{b});
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function S=PSNR(sss,aaa)

[m n p]=size(sss);
A=double(sss);
B=double(aaa);
sumaDif=0;
maxI=m*n*max(max(A.^2));
for u=1:m
    for v=1:n
        sumaDif=sumaDif+(A(u,v)-B(u,v))^2;
    end
end
if  (sumaDif==0)
    sumaDif=1;
end
S=maxI/sumaDif;
S=10*log10(S);

author: decouples


text:
function M_rice=imrician(M,s)

% function M_rice=imrician(M,s)
% ------------------------------------------------------------------------
% IMCRICIAN Random samples from the Rice/Rician probability distribution.
%
% R ~ Rice(v, s) if R = sqrt(X^2 + Y^2), where X ~ N(v*cos(a), s^2) and
% Y ~ N(v*sin(a), s^2) are independent normal distributions (any real a).
%
% Reference: http://en.wikipedia.org/wiki/Rice_distribution
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 15/04/2009
% ------------------------------------------------------------------------

X_GAUSS=s.*randn(size(M))+M;
Y_GAUSS=s.*randn(size(M));
M_rice = hypot(X_GAUSS,Y_GAUSS);

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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 3 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, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function getdepthRefined_structureIO(id)
%cd /n/fs/sun3d/code/SiftFuv2StructureIO
%/n/fs/vision/ionicNew/starter.sh getdepthRefined_structureIO 20000mb 120:00:00 1 300 1 /n/fs/modelnet/log/

basicSetup;
interval =20;
OuputPath = '/n/fs/sun3d/data/SUNRGBDv2/';
webtemplate = 'https://sun3d.cs.princeton.edu/player/?name=SUNRGBDv2/%s/&box3D=true&write=true&annotation=annotation3Dfinal&width=640&height=480&R=true&&highlight=false&rect=only';
%directory ='/net/pvd00/p/sunrgbd/mingrub/data_capture/11082015/single_image/2015-11-08T15.25.26.491/';
directory ='/Users/shurans/Downloads/2015-11-08T15.25.26.491/';
load(fullfile(directory,'metadata.mat'))
frameIDtarget = metadata.main_Id;
filename = metadata.main_frame(1:end-4);
folder = [fullfile(OuputPath,filename) '/'];
%% get web page link
webpagelink = sprintf(webtemplate,filename);
display(webpagelink)

if ~metadata.bad
    %% refine the depth
    [depthRefined, image] = SiftFuv2warpFast(directory,frameIDtarget, interval);

    %% put it into folder
    mkdir([folder 'annotation/']);
    mkdir([folder 'annotation3D/']);
    system(sprintf('chmod -R 777 %sannotation/',folder));
    system(sprintf('chmod -R 777 %sannotation3D/',folder));
    mkdir([folder 'depth/']);
    mkdir([folder 'extrinsics/']);
    mkdir([folder 'image/']);
    depth = uint16(depthRefined*1000);
    depth = bitor(bitshift(depth,3), bitshift(depth,-13));
    imwrite(depth,sprintf('%s/depth/%s.png',folder,filename));
    imwrite(image,sprintf('%s/image/%s.jpg',folder,filename));

    data = loadStructureIOdata(directory,[]);
    fid = fopen([folder 'intrinsics.txt'],'w');
    K = data.K';
    fprintf(fid,'%f %f %f\n%f %f %f\n%f %f %f\n',K(1),K(2),K(3),K(4),K(5),K(6),K(7),K(8),K(9));
    fclose(fid);
    %% get extrinsics

    depth = depthRefined;
    [x,y] = meshgrid(1:size(depth,2), 1:size(depth,1));
    XYZcamera(:,:,1) = (x-data.K(1,3)).*depth/data.K(1,1);
    XYZcamera(:,:,2) = (y-data.K(2,3)).*depth/data.K(2,2);
    XYZcamera(:,:,3) = depth;
    XYZcamera(repmat(depth==0,[1,1,3])) = NaN;
    X = XYZcamera(:,:,1);Y = XYZcamera(:,:,3);Z = -XYZcamera(:,:,2);

    [Rtilt,R] = rectify(cat(3,X,Y,Z));
    cameraRt =[eye(3) zeros(3,1)];
    cameraRt(1:3,1:3) =[1 0 0; 0 0 -1 ;0 1 0]*Rtilt*[1 0 0; 0 0 -1 ;0 1 0]';
    timeStamp = clock;
    timeStamp = sprintf('%.4d%.2d%.2d%.2d%.2d%.2d',timeStamp(1),timeStamp(2),timeStamp(3),timeStamp(4),timeStamp(5),round(timeStamp(6)));
    delete([folder '/extrinsics/*.txt'])
    fp = fopen([folder '/extrinsics/' timeStamp '.txt'],'w');
    for rowID=1:size(cameraRt,1)
        fprintf(fp, '%f %f %f %f\n',cameraRt(rowID,:));
    end

    system(sprintf('chmod -R 775 %s',folder));

end

%%
%{
if 0
    rgb = [reshape(image(:,:,1),[],1),reshape(image(:,:,2),[],1),reshape(image(:,:,3),[],1)];
    figure,
    XYZnew = Rtilt*[X(:),Y(:),Z(:)]';
    vis_point_cloud(XYZnew',double(rgb),40,4000);
end
points2ply('refined.ply', [X(:),Y(:),Z(:)], rgb)
%}
end

author: jianxiongxiao

repo: ProfXkit


text:
function tests = LocnTest
  tests = functiontests(localfunctions);
  clc
end

function setupOnce(testCase)
    testCase.TestData.Duration = 50;
end

function Vehicle_test(tc)
    %%
    randinit
    V = diag([0.005, 0.5*pi/180].^2);

    v = Bicycle('covar', V);
    v.add_driver( RandomPath(10) );

    v.run(tc.TestData.Duration);
    v.plot_xy();
    s = v.char();

    J = v.Fx(v.x, [.1 .2]);
    J = v.Fv(v.x, [.1 .2]);
end

function DeadReckoning_test(tc)
    %%
    randinit
    V = diag([0.005, 0.5*pi/180].^2);
    P0 = diag([0.005, 0.005, 0.001].^2);

    v = Bicycle('covar', V);
    v.add_driver( RandomPath(10) );
    s = char(v);

    ekf = EKF(v, V, P0);
    ekf.run(tc.TestData.Duration);

    clf
    ekf.plot_xy
    hold on
    v.plot_xy('r')
    grid on
    xyzlabel

    ekf.plot_ellipse('g')
    ekf.plot_P()
end

function MapLocalization_test(tc)
    %%
    randinit
    W = diag([0.1, 1*pi/180].^2);
    P0 = diag([0.005, 0.005, 0.001].^2);
    V = diag([0.005, 0.5*pi/180].^2);

    map = LandmarkMap(20);
    map = LandmarkMap(20, 'verbose');
    map = LandmarkMap(20, 10, 'verbose');
    map = LandmarkMap(20, 10);
    s = char(map);

    veh = Bicycle('covar', V);
    veh.add_driver( RandomPath(10) );
    sensor = RangeBearingSensor(veh, map, 'covar', W);
    sensor.interval = 5;
    ekf = EKF(veh, W, P0, sensor, W, map);

    ekf.run(tc.TestData.Duration);

    clf
    map.plot()
    veh.plot_xy('b');
    ekf.plot_xy('r');
    ekf.plot_ellipse('k')
    grid on
    xyzlabel

    clf
    ekf.plot_P()
end

function Mapping_test(tc)
    %%
    randinit
    W = diag([0.1, 1*pi/180].^2);
    V = diag([0.005, 0.5*pi/180].^2);

    map = LandmarkMap(20, 10);

    veh = Bicycle('covar', V);
    veh.add_driver( RandomPath(10) );

    sensor = RangeBearingSensor(veh, map, 'covar', W);
    sensor.interval = 5;

    ekf = EKF(veh, [], [], sensor, W, []);
    ekf.run(tc.TestData.Duration);

    clf
    map.plot()
    veh.plot_xy('b');
    ekf.plot_map('g');
    grid on
    xyzlabel

    %%
    verifyEqual(tc, numcols(ekf.landmarks), 20);

end

function SLAM_test(tc)
    %%
    randinit
    W = diag([0.1, 1*pi/180].^2);
    P0 = diag([0.005, 0.005, 0.001].^2);
    V = diag([0.005, 0.5*pi/180].^2);

    map = LandmarkMap(20, 10);

    veh = Bicycle(V);
    veh.add_driver( RandomPath(10) );

    sensor = RangeBearingSensor(veh, map, 'covar', W);
    sensor.interval = 1;

    ekf = EKF(veh, V, P0, sensor, W, []);
    ekf
    ekf.verbose = false;
    ekf.run(tc.TestData.Duration);

    clf
    map.plot()
    veh.plot_xy('b');
    ekf.plot_xy('r');
    ekf.plot_ellipse('k')
    grid on
    xyzlabel

    clf
    ekf.plot_P()

    clf
    map.plot();
    ekf.plot_map('g');

    %%
    verifyEqual(tc, numcols(ekf.landmarks), 20);

end

function ParticleFilter_test(tc)
    %%
    randinit
    map = LandmarkMap(20);

    W = diag([0.1, 1*pi/180].^2);
    v = Bicycle('covar', W);
    v.add_driver( RandomPath(10) );
    V = diag([0.005, 0.5*pi/180].^2);
    sensor = RangeBearingSensor(v, map, 'covar', V);

    Q = diag([0.1, 0.1, 1*pi/180]).^2;
    L = diag([0.1 0.1]);
    pf = ParticleFilter(v, sensor, Q, L, 1000);
    pf
    pf.run(tc.TestData.Duration);

    plot(pf.std)
    xlabel('time step')
    ylabel('standard deviation')
    legend('x', 'y', '\theta')
    grid

    clf
    pf.plot_pdf();
    clf
    pf.plot_xy();
end

function posegraph_test(tc)
    pg = PoseGraph('pg1.g2o')
    tc.verifyClass(pg, 'PoseGraph');
    tc.verifyEqual(pg.graph.n, 4);

    clf
    pg.plot()
    pg.optimize('animate')
    close all

    pg = PoseGraph('killian-small.toro')
    tc.verifyClass(pg, 'PoseGraph');
    tc.verifyEqual(pg.graph.n, 1941);

    pg = PoseGraph('killian.g2o', 'laser')
    tc.verifyClass(pg, 'PoseGraph');
    tc.verifyEqual(pg.graph.n, 3873);

    [r,theta] = pg.scan(1);
    tc.verifyClass(r, 'double');
    tc.verifyLength(r, 180);
    tc.verifyClass(theta, 'double');
    tc.verifyLength(theta, 180);

    [x,y] = pg.scan(1);
    tc.verifyClass(x, 'double');
    tc.verifyLength(x, 180);
    tc.verifyClass(y, 'double');
    tc.verifyLength(y, 180);

    pose = pg.pose(1);
    tc.verifyClass(pose, 'double');
    tc.verifySize(pose, [3 1]);

    t = pg.time(1);
    tc.verifyClass(t, 'double');
    tc.verifySize(t, [1 1]);

    w = pg.scanmap('ngrid', 3000);
    tc.verifyClass(w, 'int32');
    tc.verifySize(w, [3000 3000]);

    tc.assumeTrue(exist('idisp', 'file'));  %REMINDER
    clf
    pg.plot_occgrid(w);
    close all

end

function makemap_test(tc)
        tc.assumeTrue(false);  %REMINDER

end

function chi2inv_test(tc)
    tc.verifyEqual( chi2inv_rtb(0,2), 0);
    tc.verifyEqual( chi2inv_rtb(1,2), Inf);
    tc.verifyEqual( chi2inv_rtb(3,2), NaN);

    tc.verifyError( @() chi2inv_rtb(1,1), 'RTB:chi2inv_rtb:badarg');

end

author: petercorke

repo: robotics-toolbox-matlab


text:
function x=pentsolve(A,b)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% pentsolve.m
%
% Solve a pentadiagonal system Ax=b where A is a strongly nonsingular matrix
%
% If A is not a pentadiagonal matrix, results will be wrong
%
% Reference: G. Engeln-Muellges, F. Uhlig, "Numerical Algorithms with C"
%               Chapter 4. Springer-Verlag Berlin (1996)
%
% Written by Greg von Winckel 3/15/04
% Contact: gregvw@chtm.unm.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[M,N]=size(A);

% Check dimensions
if M~=N
    error('Matrix must be square');
    return;
end

if length(b)~=M
    error('Matrix and vector must have the same number of rows');
    return;
end

x=zeros(N,1);

% Check for symmetry
if A==A'    % Symmetric Matrix Scheme

    % Extract bands
    d=diag(A);
    f=diag(A,1);
    e=diag(A,2);

    alpha=zeros(N,1);
    gamma=zeros(N-1,1);
    delta=zeros(N-2,1);
    c=zeros(N,1);
    z=zeros(N,1);

    % Factor A=LDL'
    alpha(1)=d(1);
    gamma(1)=f(1)/alpha(1);
    delta(1)=e(1)/alpha(1);

    alpha(2)=d(2)-f(1)*gamma(1);
    gamma(2)=(f(2)-e(1)*gamma(1))/alpha(2);
    delta(2)=e(2)/alpha(2);

    for k=3:N-2
        alpha(k)=d(k)-e(k-2)*delta(k-2)-alpha(k-1)*gamma(k-1)^2;
        gamma(k)=(f(k)-e(k-1)*gamma(k-1))/alpha(k);
        delta(k)=e(k)/alpha(k);
    end

    alpha(N-1)=d(N-1)-e(N-3)*delta(N-3)-alpha(N-2)*gamma(N-2)^2;
    gamma(N-1)=(f(N-1)-e(N-2)*gamma(N-2))/alpha(N-1);
    alpha(N)=d(N)-e(N-2)*delta(N-2)-alpha(N-1)*gamma(N-1)^2;

    % Update Lx=b, Dc=z

    z(1)=b(1);
    z(2)=b(2)-gamma(1)*z(1);

    for k=3:N
        z(k)=b(k)-gamma(k-1)*z(k-1)-delta(k-2)*z(k-2);
    end

    c=z./alpha;

    % Backsubstitution L'x=c
    x(N)=c(N);
    x(N-1)=c(N-1)-gamma(N-1)*x(N);

    for k=N-2:-1:1
        x(k)=c(k)-gamma(k)*x(k+1)-delta(k)*x(k+2);
    end

else        % Non-symmetric Matrix Scheme

    % Extract bands
    d=diag(A);
    e=diag(A,1);
    f=diag(A,2);
    h=[0;diag(A,-1)];
    g=[0;0;diag(A,-2)];

    alpha=zeros(N,1);
    gam=zeros(N-1,1);
    delta=zeros(N-2,1);
    bet=zeros(N,1);

    c=zeros(N,1);
    z=zeros(N,1);

    % Factor A=LR
    alpha(1)=d(1);
    gam(1)=e(1)/alpha(1);
    delta(1)=f(1)/alpha(1);
    bet(2)=h(2);
    alpha(2)=d(2)-bet(2)*gam(1);
    gam(2)=( e(2)-bet(2)*delta(1) )/alpha(2);
    delta(2)=f(2)/alpha(2);

    for k=3:N-2
        bet(k)=h(k)-g(k)*gam(k-2);
        alpha(k)=d(k)-g(k)*delta(k-2)-bet(k)*gam(k-1);
        gam(k)=( e(k)-bet(k)*delta(k-1) )/alpha(k);
        delta(k)=f(k)/alpha(k);
    end

    bet(N-1)=h(N-1)-g(N-1)*gam(N-3);
    alpha(N-1)=d(N-1)-g(N-1)*delta(N-3)-bet(N-1)*gam(N-2);
    gam(N-1)=( e(N-1)-bet(N-1)*delta(N-2) )/alpha(N-1);
    bet(N)=h(N)-g(N)*gam(N-2);
    alpha(N)=d(N)-g(N)*delta(N-2)-bet(N)*gam(N-1);

    % Update b=Lc
    c(1)=b(1)/alpha(1);
    c(2)=(b(2)-bet(2)*c(1))/alpha(2);

    for k=3:N
        c(k)=( b(k)-g(k)*c(k-2)-bet(k)*c(k-1) )/alpha(k);
    end

    % Back substitution Rx=c
    x(N)=c(N);
    x(N-1)=c(N-1)-gam(N-1)*x(N);

    for k=N-2:-1:1
        x(k)=c(k)-gam(k)*x(k+1)-delta(k)*x(k+2);
    end

end

author: Sable

repo: mcbench-benchmarks


text:
function [okay] = inspect(mesh,varargin)
%INSPECT helper routine to safely query a MESH structure.

%-----------------------------------------------------------
%   Darren Engwirda
%   github.com/dengwirda/jigsaw-matlab
%   19-Dec-2018
%   darren.engwirda@columbia.edu
%-----------------------------------------------------------
%

    okay = false; base = ''; item = '';

    if (nargin >= +2), base = varargin{1}; end
    if (nargin >= +3), item = varargin{2}; end

    if (nargin <= +1 || nargin >= +4)
        error('inspect:incorrectNumbInputs' , ...
            'Incorrect number of arguments!') ;
    end

    if (~isstruct(mesh))
        error('inspect:incorrectInputClass' , ...
            'MESH must be a valid struct.') ;
    end
    if (~isempty(base) && ~ischar(base))
        error('inspect:incorrectInputClass' , ...
            'BASE must be a valid string!') ;
    end
    if (~isempty(item) && ~ischar(item))
        error('inspect:incorrectInputClass' , ...
            'ITEM must be a valid string!') ;
    end

    if (isempty(item))
%-- default ITEM kinds given BASE types
    switch (lower(base))
        case 'point', item = 'coord' ;
        case 'edge2', item = 'index' ;
        case 'tria3', item = 'index' ;
        case 'quad4', item = 'index' ;
        case 'tria4', item = 'index' ;
        case 'hexa8', item = 'index' ;
        case 'wedg6', item = 'index' ;
        case 'pyra5', item = 'index' ;
        case 'bound', item = 'index' ;
    end
    end

    if (isempty(item))
%-- check whether MESH.BASE exists
    okay = isfield(mesh,base) && ...
            ~isempty(mesh.(base)) ;
    else
%-- check whether MESH.BASE.ITEM exists
    okay = isfield(mesh,base) && ...
        isfield(mesh.(base),item) && ...
            ~isempty(mesh.(base).(item)) ;
    end

end

author: dengwirda


text:
function om_write_cond(condfile,c,names)
%   OM_WRITE_COND
%       [] = OM_WRITE_COND(CONDFILE,C)
%
%   Write conductivity file for OpenMEEG
%
%   Authors: Alexandre Gramfort alexandre.gramfort@inria.fr

% Copyright (C) 2010-2017, OpenMEEG developers

ndomains = length(c);

if nargin < 3
    names = {};
    for k=1:ndomains
        names{k} = ['domain', num2str(k)];
    end
end

if ndomains ~= length(names)
    error('Number of conductivities is not equal to the number of domain names.');
end

cfid = fopen(condfile, 'w');
if cfid == -1
    error(['Failed to open file ''',condfile,'''.']);
end

fprintf(cfid,'# Properties Description 1.0 (Conductivities)\n\n');
fprintf(cfid,'air         0.0\n');

for k=1:ndomains
    fprintf(cfid, '%s       %f\n', names{k}, c(k));
end

fclose(cfid);

end %  function

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ output_args ] = test_gmm(dataset)
%TEST_GMM Summary of this function goes here
%   Detailed explanation goes here
addpath('../NCM');
addpath('../common');

close all

% select dataset - '21' for Pavia University, '62' for the gulfport dataset
% '00' for generating unsupervised synthetic dataset
% '001' or '002' for running on the unsupervised synthetic image 1 or 2
if nargin < 1
    dataset = '001';
end

%% load dataset
SNR = 60;
% asphalt, basalt, concrete, conifer, grass, limestone, quartzite
chosen = [6,2,3,1];
rows = 60;
cols = 60;

options = struct('reduced_dim',10);
options.convergence_thresh = 0.0001;
options.show_fig = 0;

switch dataset
    case '00'
        Y_noise = 1e-3;
        Y_noise = rand(200,1) * Y_noise;

        gen_im_opts = [];
        gen_im_opts.chosen = chosen;
        gen_im_opts.noise_sigma = Y_noise;

        if 0 % synthetic unsupervised A
%             [I,E_gt,A_gt,names,wl,ws_gt,mus_gt,sigmas_gt,Y_noise] = ...
%                 generate_toy_image_end_var(rows, cols, Y_noise, chosen, 1);
            [I,E_gt,A_gt,names,wl,Y_noise,extra] = generate_toy_image_unified(...
                rows, cols, 'gmm', 'quadrant', 0, gen_im_opts);
            ws_gt = extra.ws;
            mus_gt = extra.mus;
            sigmas_gt = extra.sigmas;
            save('toy_image_end_var_3.mat','I','E_gt','A_gt','names','wl',...
                'ws_gt','mus_gt','sigmas_gt','Y_noise');
        elseif 0 % synthetic unsupervised B
            gen_im_opts.material_gaussian_centers_normalized = [0.5,0.25;...
                0.25,0.75; 0.75,0.75]; % (x,y)
            gen_im_opts.material_gaussian_sigma_normalized = 0.25;
            gen_im_opts.material_num_samples = 150;
            gen_im_opts.material_shape_width = 1.5;
            gen_im_opts.material_shape_width_sd = 0.5;
            gen_im_opts.material_gaussian_coefficient = 10;

            [I,E_gt,A_gt,names,wl,Y_noise,extra] = generate_toy_image_unified(...
                rows, cols, 'gmm', 'gaussian_dotted', 0, gen_im_opts);
            ws_gt = extra.ws;
            mus_gt = extra.mus;
            sigmas_gt = extra.sigmas;
            save('toy_image_end_var_3_B.mat','I','E_gt','A_gt','names','wl',...
                'ws_gt','mus_gt','sigmas_gt','Y_noise');
        else
%             load('toy_image_end_var_3.mat');
            load('toy_image_end_var_3_B.mat');
        end

%         noise_sigmas = noise_est_roger(I);
%         disp('The error of noise estimation is ');
%         mdiff(Y_noise,noise_sigmas);
%
%         noise_sigmas = noise_est_mlr(I);
%         mdiff(Y_noise,noise_sigmas);

        M = 4;
        options.beta1 = 5;
        options.beta2 = 5; % 0.005
        options.shrink_size = 5;
        options.beta2_decay = 0.05;
    case '001'
        load('toy_image_end_var_3.mat');

        M = 4;
        options.beta1 = 5;
        options.beta2 = 5; % 0.005
        options.shrink_size = 5;
        options.beta2_decay = 0.05;
    case '002'
        load('toy_image_end_var_3_B.mat');

        M = 4;
        options.beta1 = 0.1;
        options.beta2 = 0.1; % 0.005
        options.shrink_size = 5;
        options.beta2_decay = 0.05;
    case '01'
        load('toy_image_end_var_snr60_1111.mat');
        M = 4;
        options.beta1 = 2000;
        options.beta2 = 50; % 2e4
        options.eta = 0.05;
        options.rho1 = 5;
    case '2'
%         [I,R_gt,A_gt,names,wl] = load_pavia_university;
        load('../../data/PaviaUniversity_corrected.mat');
        % remove gravel and bitumen
        R_gt([3,7],:) = [];
        names([3,7]) = [];
        if 1
            A_gt(:,:,[3,7]) = [];
        else
            to_be_merged = [8,1];
            merged = [3,7];
            for j = 1:2
                A1 = A_gt(:,:,to_be_merged(j));
                A2 = A_gt(:,:,merged(j));
                A3 = A2 + A1;
                A_gt(:,:,to_be_merged(j)) = A3;
            end
            A_gt(:,:,[3,7]) = [];
        end
        M = 7;

        options.beta1 = 0.5;
        options.beta2 = 0.5; % 2e4
        options.rho1 = 0.01;
        options.rho2 = 0;
        options.shrink_size = 2;
        options.beta2_decay = 0.5;

    case '21'
        load('../../data/PaviaUniversity_A.mat');
        A_gt = double(A_gt);
        M = 5;

        options.beta1 = 5;
        options.beta2 = 5; % 2e4
        options.rho1 = 0;
        options.rho2 = 0;
        options.shrink_size = 2;
        options.beta2_decay = 0.05;
    case '3'
        [I,R_gt,A_gt,names,wl] = load_neon;
    case '6'
        [I,R_gt,A_gt,names,wl] = load_gulfport;
    case '62'
        load('../../data/muufl_gulfport_B.mat');
        M = 5;

%         remove_bands = [1:5,72-9:72];
%         [I,R_gt,wl] = remove_noisy_bands(I,R_gt,wl,remove_bands);

        options.beta1 = 5;
        options.beta2 = 5; % 2e4
        options.rho1 = 0;
        options.rho2 = 0;
        options.shrink_size = 1;
        options.beta2_decay = 0.05;
    otherwise
end
if exist('rgb','var') && ~isempty(rgb)
    I1 = rgb;
else
    I1 = retrieve_rgb(I,wl);
end
figure('name','RGB image of the original image');
imshow(I1);

[Y,A_gt,rows,cols] = reshape_hsi(I,A_gt);
[N,B] = size(Y);

if exist('ws_gt','var') && ~isempty(ws_gt)
    endmember_scatter_plot_end_var(Y,ws_gt,mus_gt,sigmas_gt,names);
    R_gt = zeros(M,B);
    for j = 1:M
        R_gt(j,:) = mean(Y(A_gt(:,j)==1,:), 1);
    end
else
    endmember_scatter_plot(Y,R_gt,{'Ground Truth'});
end
set(gcf,'name','Scatter plot of the original image with ground truth');

if exist('ws_gt','var') && ~isempty(ws_gt)
    opts = struct('show_approx',1,'w_jk',{ws_gt},'mu_jk',{mus_gt},'sigma_jk',...
        {sigmas_gt},'legend_names',{{'Ground truth GMM'}});
else
    opts = struct('legend_names',{{}});
end
hist_end_var(Y,A_gt,names,1,opts);
set(gcf,'name','Histogram of the ground truth pure pixels vs GMM');

I1 = reshape(Y, [rows, cols, B]);

D = 0.001^2 * eye(B);
options.D = D;

% [N,B] = size(Y);
% options.A_gt = A_gt;
% options.fix_A = 1;
% options.fix_sigma_jk = 1;

[A,R,w_jk,mu_jk,sigma_jk,extra] = gmm_hu(I1,M,options);
E = gmm_hu_endmember(I1,A,D,w_jk,mu_jk,sigma_jk);

%% Permute the results to accord with the GT for comparison
[error_M,error_A,best_p] = compare_2_endmembers(R_gt, R, A_gt, A, ...
    rows,cols,names,wl,1);

A = A*best_p';
R = best_p*R;
for i = 1:size(E,3)
    E(:,:,i) = best_p * E(:,:,i);
end
w_jk = w_jk(best_p*(1:M)');
mu_jk = mu_jk(best_p*(1:M)');
sigma_jk = sigma_jk(best_p*(1:M)');

if options.show_fig % Play the movie
    replay_scatter_abund(extra.frames_scatter, extra.frames_abund);
end

save(['result_gmm_',dataset,'.mat'],'A','R','E','w_jk','mu_jk','sigma_jk');

opts = struct('show_approx',1,'w_jk',{w_jk},'mu_jk',{mu_jk},'sigma_jk',...
    {sigma_jk},'legend_names',{{'Estimated distribution'}});
hist_end_var(Y,A_gt,names,1,opts);

function [ws,mus,sigmas] = create_endmember_params_1_1_1_1(E)
[M,B] = size(E);
ws = cell(1,M);
mus = cell(1,M);
sigmas = cell(1,M);

for j = 1:M
    ws{j} = 1;
    mus{j} = E(j,:);
    sigmas{j} = 0.01^2*eye(B);
end

% endmember variability of basalt
idx = 2;
major_dir = ones(B,1);
major_dir = major_dir/sqrt(sum(major_dir.^2));
sigmas{idx} = 0.01^2*eye(B) + 0.05^2*major_dir*major_dir';

% endmember variability of concrete
idx = 3;
major_dir = ones(B,1);
major_dir = major_dir/sqrt(sum(major_dir.^2));
sigmas{idx} = 0.01^2*eye(B) + 0.03^2*major_dir*major_dir';

author: zhouyuanzxcv

repo: Hyperspectral


text:
function [rect center corners]=calcRectCenter(varargin)
% function drawbox(width,height, param, properties)
%                 ([width,height], param, properties)
%
%   param, properties are optional
%

%% Copyright (C) Jongwoo Lim and David Ross.
%% All rights reserved.

%----------------------------------------------------------
% Process the input.
%----------------------------------------------------------
if (length(varargin{1}) == 2)
  w = varargin{1}(1);
  h = varargin{1}(2);
  varargin(1) = [];
else
  [w,h] = deal(varargin{1:2});
  varargin(1:2) = [];
end

if (length(varargin) < 1 || any(length(varargin{1}) ~= 6))
  M = [0,1,0; 0,0,1];
else
  p = varargin{1};
  if (length(varargin) > 1 && strcmp(varargin{2},'geom'))
    p = affparam2mat(p);
    varargin(1:2) = [];
  else
    varargin(1) = [];
  end
  M = [p(1) p(3) p(4); p(2) p(5) p(6)];
end

%----------------------------------------------------------
% Draw the box.
%----------------------------------------------------------

corners = [ 1,-w/2,-h/2; 1,w/2,-h/2; 1,w/2,h/2; 1,-w/2,h/2; 1,-w/2,-h/2 ]';
corners = M * corners;

result_corners = floor(corners(:,1:4));
x=result_corners(1,1);
y=result_corners(2,1);
w=result_corners(1,3)-result_corners(1,1);
h=result_corners(2,3)-result_corners(2,1);
rect = [x y w h];

center = mean(corners(:,1:4),2);

author: flyers

repo: drone-tracking


text:
classdef ParadigmSSRCSP < ParadigmBase
    % Selected-Subjects Regularized CSP (SSRCSP)
    %
    % This paradigm implements SSRCSP [1], which is a generalization of the Common Spatial Patterns
    % algorithms to calibration data comprising multiple subjects (or recordings). To train a model
    % for a particular "goal" (or target) subject using auxiliary data from other subjects, this
    % algorithm attempts to find a subset of other subjects such that, when data from those subjects
    % is combined with data from the goal subject, the performance on the goal subject's data is
    % optimal. The combination of data from multiple subjects and from the goal subject is done by
    % shrinking the goal-subject covariance matrix towards the average covariance matrix of the
    % other subjects (within each class), using a regularization parameter. The subset selection
    % algorithm employed is Sequential Floating Forward Selection (SFFS) [2].
    %
    % References:
    % [1] Lotte, F., & Guan, C.
    %     "Regularizing common spatial patterns to improve BCI designs: unified theory and new algorithms."
    %     Biomedical Engineering, IEEE Transactions on, 58(2), 355-362, 2011.
    %
    % [2] Pudil, P., Ferri, F. J., Novovicova, J., & Kittler, J.
    %     "Floating search methods for feature selection with nonmonotonic criterion functions."
    %     In Pattern Recognition, Proceedings of the 12th IAPR International. Conference on (Vol. 2, pp. 279-283). (1994)
    %
    % Name:
    %   Selected-Subjects Regularized Common Spatial Patterns
    %
    %                            Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
    %                            2014-02-06

    methods

        function defaults = preprocessing_defaults(self)
            defaults = {'FIRFilter',{'Frequencies',[6 8 28 32],'Type','minimum-phase'}, 'EpochExtraction',[0.5 3.5], 'Resampling',100};
        end

        function defaults = machine_learning_defaults(self)
            % set up the default parameters for machine learning
            defaults = {'lda',0.01,'regularization','shrinkage'};
        end

        function model = calibrate(self,varargin)
            % calibrate an SSRCSP model from a corpus of training sets
            args = arg_define(varargin, ...
                arg_norep({'collection','Collection'}), ...
                arg_norep({'goal_identifier','GoalIdentifier'}), ...
                arg({'patterns','PatternPairs'},3,uint32([1 1 64 10000]),'Number of CSP patterns (times two).'),...
                arg({'beta','CovarianceShrinkage'},0.5,[0 1],'Covariance shrinkage. This is the degree to which data from the goal subject is shrunken towards that of other subjects.'),...
                arg_sub({'flt','SignalProcessing'}, self.preprocessing_defaults(), @flt_pipeline, 'Signal processing stages. These parameters control filter stages that run on the signal level; they can be enabled, disabled and configured for the given paradigm. The prediction operates on the outputs of this stage.'), ...
                arg_sub({'ml','MachineLearning'},{'Learner',self.machine_learning_defaults()},@ml_train,'Machine learning stage of the paradigm. Operates on the feature vectors that are produced by the feature-extraction stage.'),...
                arg({'arg_dialogsel','ConfigLayout'},self.dialog_layout_defaults(),[],'Parameters displayed in the config dialog. Cell array of parameter names to display (dot-notation allowed); blanks are translated into empty rows in the dialog. Referring to a structure argument lists all parameters of that struture, except if it is a switchable structure - in this case, a pulldown menu with switch options is displayed.','type','cellstr','shape','row'));

            % get the data of the reference subject
            [reference,remaining] = utl_collection_closest(args.collection,args.goal_identifier);
            % preprocess each recording in the reference collection and concatenate them across epochs into a single set
            for r=1:length(reference)
                refsets{r} = exp_eval_optimized(flt_pipeline('signal',reference{r}.streams{1}, args.flt)); end
            refdata = exp_eval(set_joinepos(refsets{:}));

            % pre-process data of all other subjects
            otherdata = {};
            for s=1:length(remaining)
                if length(remaining{s}.streams) > 1
                    disp_once('Note: ParadigmMKLCSP will use only the first data stream of a recording (no support for multi-modal data).'); end
                % preprocess
                otherdata{s} = exp_eval_optimized(flt_pipeline('signal',remaining{s}.streams{1}, args.flt)); %#ok<*NODEF>
                if otherdata{s}.nbchan < args.patterns
                    error('CSP requires at least as many channels as you request output patterns. Please reduce the number of pattern pairs.'); end
            end

            % get the best subjects
            model.best_subjects = hlp_diskcache('featuremodels',@self.find_best_subjects,otherdata,refdata,args.patterns,args.ml);

            % calculate composite CSP
            covar = self.class_covariances(refdata);
            other_covar = self.class_covariances(exp_eval(set_joinepos(otherdata{model.best_subjects})));
            for k=1:2
                covar{k} = (1-args.beta)*covar{k} + args.beta*other_covar{k}; end
            [V,D] = eig(covar{1},covar{1}+covar{2}); P = inv(V); %#ok<NASGU>
            model.featuremodel.filters = V(:,[1:args.patterns end-args.patterns+1:end]);
            model.featuremodel.patterns = P([1:args.patterns end-args.patterns+1:end],:);

            % train predictive model
            model.predictivemodel = ml_train('data',{self.feature_extract(refdata,model.featuremodel),set_gettarget(refdata)}, args.ml);
            % set the filter graph based on the reference data
            model.tracking.filter_graph = refsets{end};
            % also store channel locations for model visualization
            model.chanlocs = refdata.chanlocs;
        end

        function predictions = predict(self,bundle,model)
            % extract features
            features = self.feature_extract(bundle.streams{1},model.featuremodel);
            % apply classifier
            predictions = ml_predict(features,model.predictivemodel);
        end

        function best_subjects = find_best_subjects(self,otherdata,refdata,patterns,ml)
            % find set of best subjects to include
            selected = {[]};
            remaining = {1:length(otherdata)};
            accuracy = {-Inf};
            n = 0;
            while n < length(otherdata)
                % find best subject to add
                best_accuracy = -Inf;
                best_index = NaN;
                for k = remaining{1+n}
                    acc = self.evaluate_subset(otherdata([selected{1+n} k]),refdata,patterns,ml);
                    if acc > best_accuracy
                        best_accuracy = acc;
                        best_index = k;
                    end
                end
                selected{1+n+1} = [selected{1+n} best_index];
                remaining{1+n+1} = setdiff(remaining{1+n},best_index);
                accuracy{1+n+1} = best_accuracy;
                n = n+1;
                % remove subjects
                while n > 2
                    % find best subject to remove
                    best_accuracy = -Inf;
                    best_index = NaN;
                    for k=selected{1+n}
                        acc = self.evaluate_subset(otherdata(setdiff(selected{1+n},k)),refdata,patterns,ml);
                        if acc > best_accuracy
                            best_accuracy = acc;
                            best_index = k;
                        end
                    end
                    if best_accuracy > accuracy{1+n-1}
                        selected{1+n-1} = setdiff(selected{1+n},best_index);
                        remaining{1+n-1} = [remaining{1+n} best_index];
                        accuracy{1+n-1} = best_accuracy;
                        n = n-1;
                    else
                        break;
                    end
                end
            end
            best_n = argmax([accuracy{:}]);
            best_subjects = selected{best_n};
        end

        function accuracy = evaluate_subset(self,trainset,testset,patterns,ml)
            % note: we are here implicitly weighting by the amount of training data per subject
            if iscell(trainset)
                trainset = exp_eval(set_joinepos(trainset{:})); end
            % train CSP on training set
            covar = self.class_covariances(trainset);
            [V,D] = eig(covar{1},covar{1}+covar{2}); %#ok<NASGU>
            featuremodel.filters = V(:,[1:patterns end-patterns+1:end]);
            % train classifier
            classifier = ml_train('data',{self.feature_extract(trainset,featuremodel),set_gettarget(trainset)},ml);
            % test on test set
            accuracy = -ml_calcloss('auto',set_gettarget(testset),ml_predict(self.feature_extract(testset,featuremodel),classifier));
        end

        function covar = class_covariances(self,dataset)
            % calculate per-class covariance matrices
            for k=1:2
                classdata = exp_eval(set_picktrials(dataset,'rank',k));
                covar{k} = (classdata.data(:,:) * classdata.data(:,:)') / (size(classdata.data,2)*size(classdata.data,3));
                covar{k}(~isfinite(covar{k})) = 0;
            end
        end

        function features = feature_extract(self,signal,featuremodel)
            % extract log-variance features from an epoched and preprocessed recording
            features = zeros(size(signal.data,3),size(featuremodel.filters,2));
            for t=1:size(signal.data,3)
                features(t,:) = sum((signal.data(:,:,t)'*featuremodel.filters).^2,1); end
            features = log(features/size(signal.data,2));
        end

        function visualize(self,varargin) %#ok<*INUSD>
            % visualize an mklCSP model
            args = arg_define(varargin, ...
                arg_norep({'model','Model'},[],[],'BCI Model to visualize.'), ...
                arg({'patterns','PlotPatterns'},true,[],'Plot patterns instead of filters. Whether to plot spatial patterns (forward projections) rather than spatial filters.'), ...
                arg({'paper','PaperFigure'},false,[],'Use paper-style font sizes. Whether to generate a plot with font sizes etc. adjusted for paper.'));

            f = figure;
            % get number of pairs, and index of pattern per subplot
            np = size(args.model.featuremodel.patterns,1)/2;
            idx = [1:np 2*np:-1:np+1];
            % for each CSP pattern...
            for p=1:np*2
                subplot(2,np,p,'Parent',f);
                if args.patterns
                    topoplot(args.model.featuremodel.patterns(idx(p),:),args.model.featuremodel.chanlocs);
                else
                    topoplot(args.model.featuremodel.filters(:,idx(p)),args.model.featuremodel.chanlocs);
                end
                t = title(['CSP Pattern ' num2str(idx(p))]);
                if args.paper
                    set(t,'FontUnits','normalized');
                    set(t,'FontSize',0.1);
                end
            end
        end

        function layout = dialog_layout_defaults(self)
            % define the default configuration dialog layout
            layout = {'SignalProcessing.Resampling.SamplingRate', 'SignalProcessing.FIRFilter.Frequencies', ...
                'SignalProcessing.FIRFilter.Type', 'SignalProcessing.EpochExtraction', '', ...
                'PatternPairs', 'CovarianceShrinkage', '', 'MachineLearning.Learner'};
        end

    end
end

% (turn off a few editor warnings because some actual implementations are missing in this file)
%#ok<*INUSD,*STOUT,*MANU>

repo: REST


text:
clear all

dt = 0.1;
t  = 0:dt:10;

Nsamples = length(t);

Xsaved   = zeros(Nsamples, 2);
DeXsaved = zeros(Nsamples, 2);
Zsaved   = zeros(Nsamples, 1);

for k=1:Nsamples
  z = GetPos();
  [pos vel]   = DvKalman(z);
  [dpos dvel] = DeDvKalman(z);

  Xsaved(k, :)   = [pos vel];
  DeXsaved(k, :) = [dpos dvel];
  Zsaved(k)      = z;
end

figure
hold on
plot(t, Xsaved(:, 1), 'ro')
plot(t, DeXsaved(:,1))

figure
hold on
plot(t, Xsaved(:, 2), 'ro')
plot(t, DeXsaved(:,2))

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
classdef lsl_streaminfo < handle
	% The stream_info object stores the declaration of a data stream.
	% Represents the following information:
	%  a) stream data format (#channels, channel format)
	%  b) core information (stream name, content type, sampling rate)
	%  c) optional meta-data about the stream content (channel labels, measurement units, etc.)
	%
	% Available to recipients of a stream and also written to disk when being recorded
	% (corresponds to the file header).

    properties
        LibHandle = 0;      % this is a handle to the liblsl library (to call its functions)
        InfoHandle = 0;     % this is a handle to an lsl_streaminfo object within the library.
    end

    methods

        function self = lsl_streaminfo(lib,name,type,channelcount,samplingrate,channelformat,sourceid)
            % Construct a new stream_info object.
            % StreamInfo = lsl_streaminfo(LibHandle,Name,Type,ChannelCount,SamplingRate,ChannelFormat,SourceId)
            %
            % Core stream information is specified here. Any remaining meta-data can be added later.
            %
            % In:
            % 	Name : Name of the stream. Describes the device (or product series) that this stream makes available
            %		   (for use by programs, experimenters or data analysts).
            %
            %   Type : Content type of the stream. Please see Table of Content Types in the documentation
            %          for naming recommendations.
            %
            %   ChannelCount : Number of channels per sample. This stays constant for the lifetime of the stream.
            %
            %   NominalSrate : The sampling rate (in Hz) as advertised by the data source,
            %                  if regular (otherwise set to IRREGULAR_RATE).
            %
            %   ChannelFormat : Format/type of each channel. If your channels have different formats,
            %                   consider supplying multiple streams or use the largest type that can
            %                   hold them all (such as cf_double64). The possible formats are:
            %                   'cf_float32', 'cf_double64', 'cf_string', 'cf_int32', 'cf_int16',
            %                   'cf_int8'
            %
            %   SourceId: Unique identifier of the source or device, if available (such as the serial number).
            %             Allows recipients to recover from failure even after the serving app or device crashes.
            %			  May in some cases also be constructed from device settings.

            self.LibHandle = lib;

            if ~exist('type','var') && ~ischar(name)
                % assume that the second argument is a handle
                self.InfoHandle = name;
            else
                % parse inputs
                if ~exist('channelcount','var') || isempty(channelcount) channelcount = 1; end
                if ~exist('samplingrate','var') || isempty(samplingrate) samplingrate = 0; end
                if ~exist('channelformat','var') || isempty(channelformat) channelformat = 'cf_float32'; end
                if ~exist('sourceid','var') || isempty(sourceid) sourceid = ''; end

                % translate channel formats
                if ischar(channelformat)
                    switch(channelformat)
                        case 'cf_float32'       % For up to 24-bit precision measurements in the appropriate physical unit
                            channelformat = 1;  % (e.g., microvolts). Integers from -16777216 to 16777216 are represented accurately.
                        case 'cf_double64'      % For universal numeric data as long as permitted by network & disk budget.
                            channelformat = 2;  % The largest representable integer is 53-bit.
                        case 'cf_string'        % For variable-length ASCII strings or data blobs, such as video frames,
                            channelformat = 3;  % complex event descriptions, etc.
                        case 'cf_int32'         % For high-rate digitized formats that require 32-bit precision. Depends critically on
                            channelformat = 4;  % meta-data to represent meaningful units. Useful for application event codes or other coded data.
                        case 'cf_int16'         % For very high rate signals (40Khz+) or consumer-grade audio
                            channelformat = 5;  % (for professional audio float is recommended).
                        case 'cf_int8'          % For binary signals or other coded data.
                            channelformat = 6;  % Not recommended for encoding string data.
                        case 'cf_int64'         % For future compatibility. Support for this type is not yet widespread in programming languages
                            channelformat = 7;  % and therefore not exposed in the API for all languages. Some builds of liblsl will not be able to
                                                % send or receive data of this type.
                        otherwise
                            error('Unsupported channel format.');
                    end
                end

                % create info object
                self.InfoHandle = lsl_create_streaminfo(lib,name,type,channelcount,samplingrate,channelformat,sourceid);
            end
        end

        function delete(self)
            % Destroy the streaminfo when it is deleted.
            lsl_destroy_streaminfo(self.LibHandle,self.InfoHandle);
        end

		% === Core Information ===
		% (these fields are assigned at construction)

        function result = name(self)
            % Get the name of the stream.
            % This is a human-readable name. For streams offered by device modules, it refers to the type of device or product series
            % that is generating the data of the stream. If the source is an application, the name may be a more generic or specific identifier.
            % Multiple streams with the same name can coexist, though potentially at the cost of
            % ambiguity (for the recording app or experimenter).

            result = lsl_get_name(self.LibHandle,self.InfoHandle);
        end

        function result = type(self)
            % Get the content type of the stream.
            % The content type is a short string such as "EEG", "Gaze" which describes the content carried by the channel (if known).
            % If a stream contains mixed content this value need not be assigned but may instead be stored in the description of channel types.
            % To be useful to applications and automated processing systems using the recommended content types is preferred.
            % See Table of Content Types in the documentation.

            result = lsl_get_type(self.LibHandle,self.InfoHandle);
        end

        function result = channel_count(self)
            % Get the number of channels of the stream.
            % A stream has at least one channel; the channel count stays constant for all samples.

            result = lsl_get_channel_count(self.LibHandle,self.InfoHandle);
        end

        function result = nominal_srate(self)
            % Get the sampling rate of the stream, according to the source (in Hz).
            % If a stream is irregularly sampled, this should be set to IRREGULAR_RATE.
            %
            % Note that no data will be lost even if this sampling rate is incorrect or if a device has temporary
            % hiccups, since all samples will be recorded anyway (except for those dropped by the device itself). However,
            % when the recording is imported into an application, a good importer may correct such errors more accurately
            % if the advertised sampling rate was close to the specs of the device.

            result = lsl_get_nominal_srate(self.LibHandle,self.InfoHandle);
        end

        function result = channel_format(self)
            % Get the channel format of the stream.
            % All channels in a stream have the same format. However, a device might offer multiple time-synched streams
            % each with its own format.

            channel_format_strings = {'cf_float32','cf_double64','cf_string','cf_int32','cf_int16','cf_int8','cf_int64'};
            result = channel_format_strings{lsl_get_channel_format(self.LibHandle,self.InfoHandle)};
        end

        function result = source_id(self)
            % Get the unique identifier of the stream's source, if any.
            % The unique source (or device) identifier is an optional piece of information that, if available, allows that
            % endpoints (such as the recording program) can re-acquire a stream automatically once it is back online.

            result = lsl_get_source_id(self.LibHandle,self.InfoHandle);
        end

		% === Additional Network Info ===
		% (these fields are auto-assigned once the streaminfo is bound to an outlet)

        function result = version(self)
            % Get the protocol version used to deliver the stream.

            result = lsl_get_version(self.LibHandle,self.InfoHandle);
        end

        function result = created_at(self)
            % Get the creation time stamp of the stream.
            % This is the time stamp when the stream was first created
            % (as determined via lsl_local_clock() on the providing machine).

            result = lsl_get_created_at(self.LibHandle,self.InfoHandle);
        end

        function result = uid(self)
            % Get the unique ID of the stream outlet (once assigned).
            % This is a unique identifier of the stream outlet, and is guaranteed to be different
            % across multiple instantiations of the same outlet (e.g., after a re-start).

            result = lsl_get_uid(self.LibHandle,self.InfoHandle);
        end

        function result = session_id(self)
            % Get the session ID for the given stream.
            % The session id is an optional human-assigned identifier of the recording session.
            % While it is rarely used, it can be used to prevent concurrent recording activitites
            % on the same sub-network (e.g., in multiple experiment areas) from seeing each other's streams
            % (assigned via a configuration file by the experimenter, see Configuration File in the docs).

            result = lsl_get_session_id(self.LibHandle,self.InfoHandle);
        end

        function result = desc(self)
            % Get access to the extended description of the stream.
            % It is highly recommended that at least the channel labels are described here.
            % See code examples in the documentation. Other information, such as amplifier settings,
            % measurement units if deviating from defaults, setup information, subject information, etc.,
            % can be specified here, as well. See Meta-Data Recommendations in the docs.

            result = lsl_xml_ptr(self.LibHandle,lsl_get_desc(self.LibHandle,self.InfoHandle));
        end

        function result = as_xml(self)
            % Retrieve the entire stream_info in XML format.
            % This yields an XML document (in string form) whose top-level element is <info>. The info element contains
            % one element for each field of the stream_info class, including:
            %  a) the core elements <name>, <type>, <channel_count>, <nominal_srate>, <channel_format>, <source_id>
            %  b) the misc elements <version>, <created_at>, <uid>, <session_id>, <v4address>, <v4data_port>, <v4service_port>, <v6address>, <v6data_port>, <v6service_port>
            %  c) the extended description element <desc> with user-defined sub-elements.

            result = lsl_get_xml(self.LibHandle,self.InfoHandle);
        end

    end
end

author: PatternRecognition

repo: OpenBMI


text:
%% ParticleFilterX Demo
% ----------------------
% * This script demonstrates the process of configuring an running a
%   ParticleFilterX object to perform single-target state estimation.
%
% * A toy single-target scenario is considered, for a target that generates
%   regular reports of it's position in the absence of clutter and/or
%   missed detection.
%
%% Extract the GroundTruth data from the example workspace
load('single-target-tracking.mat');
NumIter = size(TrueTrack.Trajectory,2);

%% Models
% Instantiate a Transitionamic model
transition_model = ConstantHeadingX('VelocityErrVariance',0.01^2, 'HeadingErrVariance', (pi/20)^2);

% Instantiate an Observation model
measurement_model = LinearGaussianX('NumMeasDims',2,'NumStateDims',4,'MeasurementErrVariance',0.02,'Mapping',[1 3]);
%measurement_model = RangeBearing2CartesianX('NumStateDims',4,'MeasurementErrVariance',[0.001,0.02],'Mapping',[1 3]);

% Compile the State-Space model
model = StateSpaceModelX(transition_model,measurement_model);

%% Simulation
% Data Simulator
dataSim = SingleTargetMeasurementSimulatorX(model);

% Simulate some measurements from ground-truth data
MeasurementScans = dataSim.simulate(TrueTrack);
measurements = [MeasurementScans.Vectors];
measurement_model = LinearGaussianX('NumMeasDims',2,'NumStateDims',4,'MeasurementErrVariance',0.02,'Mapping',[1 2]);
model.Measurement = measurement_model;

%% Initiation
% Use the first measurement scan to perform single-point initiation
measurement = MeasurementScans(1).Measurements;
timestamp = measurement.Timestamp;

% Setup prior
xPrior = measurement_model.finv(measurement.Vector);
PPrior = 10*transition_model.covar();
dist = GaussianDistributionX(xPrior,PPrior);
StatePrior = ParticleStateX(dist,5000,timestamp);

% Initiate a track using the generated prior
track = TrackX(StatePrior, TagX(1));

%% Estimation
% Instantiate a filter objects
filter = ParticleFilterX('Model',model);
figure;
for t = 2:NumIter

    % Provide filter with the new measurement
    MeasurementList = MeasurementScans(t);

    % Perform filtering
    prior = track.State;
    prediction = filter.predict(prior, MeasurementList{1}.Timestamp);
    posterior = filter.update(prediction, MeasurementList);

    % Log the data
    track.Trajectory(end+1) = posterior;

    clf;
    hold on;
    meas = measurement_model.finv(measurements(:,1:t));
    true_means = [TrueTrack.Trajectory(1:t).Vector];
    track_means = [track.Trajectory(1:t).Mean];
    plot(true_means(1,1:t), true_means(2,1:t),'.-k', track_means(1,1:t), track_means(2,1:t), 'b-', meas(1,:), meas(2,:), 'rx');
    plot_gaussian_ellipsoid(track.Trajectory(t).Mean([1,2],1), track.Trajectory(t).Covar([1,2],[1,2]));
    legend('GroundTrouth','Estimated Mean','Measurements', 'Estimated Covariance');
    xlabel("x coordinate (m)");
    ylabel("y coordinate (m)");
    axis([2 9 1 9]);
    drawnow();
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function prob_test133 ( )

%*****************************************************************************80
%
%% TEST133 tests RECIPROCAL_CDF, RECIPROCAL_CDF_INV, RECIPROCAL_CDF;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST133\n' );
  fprintf ( 1, '  For the Reciprocal PDF:\n' );
  fprintf ( 1, '  RECIPROCAL_CDF evaluates the CDF.\n' );
  fprintf ( 1, '  RECIPROCAL_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  RECIPROCAL_PDF evaluates the PDF.\n' );

  a = 1.0;
  b = 3.0;

  check = reciprocal_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST133 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =         %14f\n', a );
  fprintf ( 1, '  PDF parameter B =         %14f\n', b );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '       X            PDF           CDF            CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ x, seed ] = reciprocal_sample ( a, b, seed );

    pdf = reciprocal_pdf ( x, a, b );

    cdf = reciprocal_cdf ( x, a, b );

    x2 = reciprocal_cdf_inv ( cdf, a, b );

    fprintf ( 1, ' %14f  %14f  %14f  %14f\n', x, pdf, cdf, x2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [AddedRxns] = submitGapFill(swift, modelFile, dbFile, dictionaryFile, prepareFGFResults, weightsPerRxnFile, forceRerun, epsilon, blackList, listCompartments)
% A test function for `prepareFastGapFill`, `swiftGapFill`, and
% `fastGapFill`, allowing all files to be optionally specified and not
% requiring `prepareFastGapFill` to be rerun (if relevant variables are saved
% in a specified workspace file) if it has already been run.
%
% This function performs three runs of `fastGapFill` or `swiftGapFill on the
% same prepared model, with weights as specified in the 'Define Runs' section
% of the code.  The 'runs' variable can be appended to or changed to enable
% the running of `fastGapFill` or `swiftGapFill` for different combinations
% of weights (and different `weightsPerRxnFile` files).
%
% This function draws inspiration heavily from the example already
% available for `fastGapFill`, i.e. `runGapFill_example.m`
%
% N.B. The defaults are set for testing all of the functionality of the
% algorithm, including setting a (toy) weight file, which is not necessary
% in running `fastGapFill` or `swiftGapFill` so should be specified as empty
% if running without a weight file.
%
% All inputs are optional, defaults to E. coli model iAF1260 and example files)
%
% USAGE:
%
%    [AddedRxns] = submitGapFill(swift, modelFile, dbFile, dictionaryFile, prepareFGFResults, weightsPerRxnFile, forceRerun, epsilon, blackList, listCompartments)
%
% OPTIONAL INPUTS:
%    swift:                Boolean specifying whether to run swiftGapFill or
%                          fastGapFill (default: false)
%    modelFileIn:          File containing model, either .mat or .xml
%                          (default: 'examples/iAF1260.mat')
%    dbFileIn:             File containing universal database
%                          (default: 'AuxillaryFiles/reaction.lst')
%    dictionaryFileIn:     File containing db metabolite IDs and model
%                          counterparts, either .xls or .tsv
%                          (default: 'AuxillaryFiles/KEGG_dictionary.xls')
%    workspaceFileIn:      File for storing `prepareFastGapFill` results
%                          (default: 'examples/defaultWorkspace.mat')
%    weightsPerRxnFile:    File containing individual weights for reactions
%                          (default: 'examples/sampleWeights.tsv')
%    forceRerun:           Boolean specifying whether to rerun
%                          `prepareFastGapFill` even if it has already been
%                          run. N.B. If this is set to 'true' it will overwrite
%                          the precalculated default results files, unless
%                          `prepareFGFResults` is specified outside the examples
%                          directory (default: false)
%    epsilon:              Float, a fastCore parameter (default: getCobraSolverParams('LP', 'feasTol')*100)
%    blackList:            List of excluded universal DB reactions
%                          (default: none)
%    listCompartments:     List of compartments in the model to be gapFilled N.B.
%                          the default is set within `prepareGapFill`
%                          (default: all compartments specified in the model)
%
% OUTPUT:
%    AddedRxns:            Reactions that have been added from UX matrix to `S`
%
% .. Author:
%   - Will Bryant, Jan 2016
%   - Mojtaba Tefagh, May 2019

global CBTDIR

warning('off','MATLAB:load:variableNotFound') % Preparation - load data. Suppress load warnings and deal with in the function

% Get folder for finding default files relative to submitGapFill
runDir = [CBTDIR filesep 'src' filesep 'reconstruction' filesep 'fastGapFill' filesep];

% fastGapFill results files will be created in the same directory as prepareFGFResults
resultsDir = [CBTDIR filesep 'tutorials' filesep 'reconstruction' filesep 'fastGapFill' filesep];

% Get input files where specified
if ~exist('swift','var') || isempty(swift)
    swift=false;
end
if ~exist('modelFile','var') || isempty(modelFile)
    modelFile = strcat(runDir, ['examples' filesep 'iAF1260.mat']);
end
if ~exist('dbFile','var') || isempty(dbFile)
    dbFile = strcat(runDir, ['AuxillaryFiles' filesep 'reaction.lst']);
end
if ~exist('dictionaryFile','var') || isempty(dictionaryFile)
    dictionaryFile = strcat(runDir, ['AuxillaryFiles' filesep 'KEGG_dictionary.xls']);
end
if ~exist('prepareFGFResults','var') || isempty(prepareFGFResults)
    prepareFGFResults =  strcat(resultsDir, ['example' filesep 'prepareFGFResultsDefault.mat']);
end
if ~exist('weightsPerRxnFile','var') || isempty(weightsPerRxnFile)
    weightsPerRxnFile = strcat(resultsDir, ['example' filesep 'sampleWeights.tsv']);
end

if ~exist('forceRerun','var') || isempty(forceRerun)
    forceRerun=false;
end
if ~exist('epsilon','var') || isempty(epsilon)
    epsilon=getCobraSolverParams('LP', 'feasTol')*100;
end
if ~exist('blackList','var') || isempty(blackList)
    blackList=[];
end
if ~exist('listCompartments','var') || isempty(listCompartments)
    listCompartments=[];
end

%% prepareFastGapFill
% If workspaceFile is present, check for all variables; if consistModel,
% consistMatricesSUX and BlockedRxns are present, do not rerun
% prepareFastGapFill, otherwise run and save workspace; check for but do
% not require prepStats
try
    a = load(prepareFGFResults,'consistModel','consistMatricesSUX','BlockedRxns','model');
catch
    fprintf('Workspace file not found, proceeding with prepareGapFill ...\n')
    a = struct();
end
length_a = length(fieldnames(a));
clear a;
if (length_a == 4) && ~forceRerun
    fprintf('prepareGapFill already run, and forceRerun set to "false", prepareGapFill will not be rerun\n');
else

    % load model using relevant load function
    if regexp(modelFile,'.mat$')
        model = readMlModel(modelFile);
    elseif regexp(modelFile,'.xml$')
        model = readCbModel(modelFile);
        % If subSystems is empty, create a dummy subSystems for mergeTwoModels
        if ~exist('model.subSystems') || length(model.subSystems) ~= length(model.rxnNames)
            model.subSystems = repmat({{''}},length(model.rxns),1);
        end
        if ~exist('model.genes')
            model.genes = repmat({'no_gene'},1);
        end
        if ~exist('model.rxnGeneMat')
            model.rxnGeneMat = zeros(length(model.rxns),1);
        end
        if ~exist('model.grRules')
            model.grRules = repmat({''},length(model.rxns),1);
        end
    end

    % remove constraints from exchange reactions
    EX = strncmp('EX_',model.rxns,3);
    model.lb(EX)=-100;
    model.ub(EX)=100;
    clear EX

%     % Switch all model IDs to lower-case
%     model.mets = cellfun(@lower,model.mets,'UniformOutput',false);

    tic;
    [consistModel,consistMatricesSUX,BlockedRxns] = prepareFastGapFill(model, listCompartments, epsilon, dbFile, dictionaryFile, blackList);
    tpre=toc;

    % Prepare the output table with statistics
    cnt = 1;
    prepStats{cnt,1} = 'Model name';cnt = cnt+1;
    prepStats{cnt,1} = 'Size S (original model)';cnt = cnt+1;
    prepStats{cnt,1} = 'Number of compartments';cnt = cnt+1;
    prepStats{cnt,1} = 'List of compartments';cnt = cnt+1;
    prepStats{cnt,1} = 'Number of blocked reactions';cnt = cnt+1;
    prepStats{cnt,1} = 'Number of solvable blocked reactions';cnt = cnt+1;
    prepStats{cnt,1} = 'Size S (flux consistent)';cnt = cnt+1;
    prepStats{cnt,1} = 'Size SUX (including solvable blocked reactions)';cnt = cnt+1;
    prepStats{cnt,1} = 'Time preprocessing [min]';

    % get stats
    cnt = 1;
    prepStats{cnt,2} = modelFile;cnt = cnt+1;
    [a,b] = size(model.S);
    prepStats{cnt,2} = strcat(num2str(a),'x',num2str(b));cnt = cnt+1;
    [~,rem] = strtok(model.mets,'\[');
    rem = unique(rem);
    prepStats{cnt,2} = num2str(length(rem));cnt = cnt+1;
    Rem = rem{1};
    for j = 2:length(rem)
        Rem = strcat(Rem,',',rem{j});
    end
    prepStats{cnt,2} = Rem;cnt = cnt+1;
    clear Rem rem;
    prepStats{cnt,2} = num2str(length(BlockedRxns.allRxns));cnt = cnt+1;
    prepStats{cnt,2} = num2str(length(BlockedRxns.solvableRxns));cnt = cnt+1;
    [a,b] = size(consistModel.S);
    prepStats{cnt,2} = strcat(num2str(a),'x',num2str(b));cnt = cnt+1;
    [a,b] = size(consistMatricesSUX.S);
    prepStats{cnt,2} = strcat(num2str(a),'x',num2str(b));cnt = cnt+1;
    prepStats{cnt,2} = num2str(tpre/60);cnt = cnt+1;

    % Save data
    save(prepareFGFResults,'consistModel','consistMatricesSUX','BlockedRxns','prepStats','model');
    fprintf('prepareGapFill finished\n')
end

if swift
    fprintf('Running swiftGapFill ...\n')
else
    fprintf('Running fastGapFill ...\n')
end

%% Define Runs
% define weights for reactions to be added - the lower the weight the
% higher the priority
% WARNING: weights should not be 1, as this impacts on the algorithm
% performance

runs = {};

% RUN 1
run.weightsPerReactionFile = weightsPerRxnFile;
run.weights.MetabolicRxns = 0.1; % Kegg metabolic reactions
run.weights.ExchangeRxns = 0.5; % Exchange reactions
run.weights.TransportRxns = 10; % Transport reactions
run.name = 'initial';
runs = [runs; run];

% RUN 2
run.weightsPerReactionFile = weightsPerRxnFile;
run.weights.MetabolicRxns = 0.1; % Kegg metabolic reactions
run.weights.ExchangeRxns = 0.9; % Exchange reactions
run.weights.TransportRxns = 40; % Transport reactions
run.name = 'hw';
runs = [runs; run];

% RUN 3
run.weightsPerReactionFile = '';
run.weights.MetabolicRxns = 0.01; % Kegg metabolic reactions
run.weights.ExchangeRxns = 2.1; % Exchange reactions
run.weights.TransportRxns = 20; % Transport reactions
run.name = 'vhw';
runs = [runs; run];

clear run

%% fastGapFill/swiftGapFill

cnt = 1;
clear Stats;
Stats{cnt,1} = 'Run name';cnt = cnt+1;
Stats{cnt,1} = 'Number of added reactions (all)';cnt = cnt+1;
Stats{cnt,1} = 'Number of added metabolic reactions ';cnt = cnt+1;
Stats{cnt,1} = 'Number of added transport reactions ';cnt = cnt+1;
Stats{cnt,1} = 'Number of added exchange reactions ';cnt = cnt+1;
if swift
    Stats{cnt,1} = 'Time swiftGapFill';cnt = cnt+1;
else
    Stats{cnt,1} = 'Time fastGapFill';cnt = cnt+1;
end
Stats{cnt,1} = 'SFilename';cnt = cnt+1;

try
    a = load(prepareFGFResults,'consistModel','consistMatricesSUX','BlockedRxns','model');
    assert(length(fieldnames(a)) == 4);
    clear a;
    load(prepareFGFResults,'consistModel','consistMatricesSUX','BlockedRxns','model');
    fprintf('Variables loaded from .mat file\n');
    try
        b = load(prepareFGFResults,'prepStats');
        assert(~isempty(b))
        load(prepareFGFResults,'prepStats');
        clear b;
    catch
        fprintf('Warning: no stats for prepareFastGapfill\n');
    end
catch
    fprintf('Using the preloaded variables\n');
end

for i = 1:length(runs)
    run_idx = i+1;
    fprintf('\nRun %i (%s)\n',i,runs{i}.name);
    weights = runs{i}.weights;
    weightsPerReactionFile = runs{i}.weightsPerReactionFile;
    resultsFile = strcat(resultsDir,'results_',runs{i}.name,'.mat');

    % If specified, import weights from weights file (tsv)
    if ~isempty(weightsPerReactionFile)
        file_handle = fopen(weightsPerReactionFile);
        try
            u = textscan(file_handle,'%s\t%s');
        catch
            fprintf('File "%s" could not be found, exiting\n',weightsPerReactionFile)
            return
        end
        weightsPerReaction.rxns = {};
        weightsPerReaction.weights = {};
        for k = 1:length(u{1})
            weightsPerReaction.rxns{k} = u{1}{k};
            weightsPerReaction.weights{k} = str2double(u{2}{k});
        end
        fclose(file_handle);
    else
        weightsPerReaction = [];
    end

    cnt = 1;
    Stats{cnt,run_idx} = runs{i}.name;cnt = cnt+1;

    % fastGapFill/swiftGapFill
    tic;
    if swift
        [~, AddedRxns] = swiftGapFill(consistMatricesSUX,1e-10,weights,weightsPerReaction);
    else
        [AddedRxns] = fastGapFill(consistMatricesSUX,epsilon,weights,weightsPerReaction);
    end
    tgap=toc;
    Stats{cnt,run_idx} = num2str(length(AddedRxns.rxns));cnt = cnt+1;
    save(resultsFile);

    % Postprocessing
    [AddedRxnsExtended] = postProcessGapFillSolutions(AddedRxns,model,BlockedRxns,0);

    try
        Stats{cnt,run_idx} = num2str(AddedRxnsExtended.Stats.metabolicSol);cnt = cnt+1;
        Stats{cnt,run_idx} = num2str(AddedRxnsExtended.Stats.transportSol);cnt = cnt+1;
        Stats{cnt,run_idx} = num2str(AddedRxnsExtended.Stats.exchangeSol);cnt = cnt+1;
    catch
        Stats{cnt,run_idx} = num2str(0);cnt = cnt+1;
        Stats{cnt,run_idx} = num2str(0);cnt = cnt+1;
        Stats{cnt,run_idx} = num2str(0);cnt = cnt+1;
    end
    Stats{cnt,run_idx} = num2str(tgap);cnt = cnt+1;
    Stats{cnt,run_idx} = resultsFile;cnt = cnt+1;
    clear a b

    % Reaction List
    col = 1;
    RxnList={};
    try
        RxnList{1,col}=resultsFile;RxnList(2:length(AddedRxnsExtended.rxns)+1,col) = AddedRxnsExtended.rxns; col = col + 1;
        RxnList{1,col}=resultsFile;RxnList(2:length(AddedRxnsExtended.rxns)+1,col) = AddedRxnsExtended.rxnFormula; col = col + 1;
        RxnList{1,col}=resultsFile;RxnList(2:length(AddedRxnsExtended.rxns)+1,col) = AddedRxnsExtended.subSystem; col = col + 1;
    catch
        RxnList{1,col}=resultsFile; col = col + 1;
        RxnList{1,col}=resultsFile; col = col + 1;
        RxnList{1,col}=resultsFile; col = col + 1;
    end

    save(resultsFile);
end

author: opencobra

repo: cobratoolbox


text:
function [options,data,ndim] = checkoptions_spectra (options,data,T,parametric)

if nargin>=3
    if iscell(T)
        if size(T,1)==1, T = T'; end
        for i = 1:length(T)
            if size(T{i},1)==1, T{i} = T{i}'; end
        end
        T = cell2mat(T); T = T(:);
    end
end

% if isfield(options,'pca') && options.pca~=0 && options.pca~=1
%     options.pca = 0;
% end

if ~isempty(data)
    if iscell(data) && ischar(data{1})
        fsub = data{1};
        loadfile_sub;
        data = X;
        ndim = size(data,2);
    elseif iscell(data)
        data = data{1};
        ndim = size(data,2);
    elseif isstruct(data)
        ndim = size(data.X,2);
    else
        ndim = size(data,2);
    end
else
    if isfield(options,'S'), ndim = size(options.S,1);
    else, ndim = 1;
    end
    data = rand(10,ndim); T = 10; % useless
end

% MT and common
if ~isfield(options,'p'), options.p = 0; end
if ~isfield(options,'removezeros'), options.removezeros = 0; end
if ~isfield(options,'completelags'), options.completelags = 0; end
if ~isfield(options,'rlowess'), options.rlowess = 0; end
if ~isfield(options,'numIterations'), options.numIterations = 100; end
if ~isfield(options,'tol'), options.tol = 1e-18; end
if ~isfield(options,'pad'), options.pad = 0; end
if ~isfield(options,'Fs'), options.Fs=1; end
if ~isfield(options,'fpass'),  options.fpass=[options.Fs/200 options.Fs/2]; end
mfs = max(options.Fs/200, options.fpass(1));
if ~isfield(options,'win'),  options.win = min(4*options.Fs/mfs, min(T));  end
if ~isfield(options,'tapers'), options.tapers = [4 7]; end
if ~isfield(options,'verbose'), options.verbose = 0; end
if ~isfield(options,'to_do')
    if ndim>1, options.to_do = ones(2,1); options.to_do(2) = parametric;
    else, options.to_do = zeros(2,1);
    end
end
if nargin==3 && any(T<options.win)
    error('The specified window is larger than some of the segments')
end

options.updateGamma = 1;
options.DirichletDiag = 100; % to avoid a warning later
[options,data] = checkoptions(options,data,[]);
options.leida = 0;
options.onpower = 0;
options.pca = 0;
options.pca_spatial = 0;
options.embeddedlags = 0;
if isfield(options,'As'), options = rmfield(options,'As'); end
if isfield(options,'A'), options = rmfield(options,'A'); end
if isstruct(data), data = data.X; end

% MAR
% if ~isfield(options,'loadings'), options.loadings=eye(ndim); end;
if ~isfield(options,'Nf'), options.Nf = 256; end
if ~isfield(options,'completelags'), options.completelags = 1; end

end

author: OHBA-analysis

repo: HMM-MAR


text:
% Sample of Linear Blend Skinning in Matlab
% This code is non optimized nor necessarily correct.
% Also, it has only been tested on the data in this folder so I dont know
% how it will behave on other data.

% Aaron Wetzler, 2013

%% -------------------------------------------------------------------------
clf
close all
clear all
clc

% Load the faces, vertices and weights
faces = csvread('handfaces.txt');
V0 = csvread('handverts.txt');
V0=V0';
N = size(V0,1);
weights = csvread('handweights.txt');
sk = loadbvh('HandBase.bvh');

% Mapping to correct the ordering of the weights to the bones for the
% HandBase.bvh file
wmap = [1 2 3 4 0 15 16 17 0 5 6 7 0 8 9 10 0 11 12 13 0 14 0];
wzero = wmap==0;

% This will likely crash a machine if there are two many bones and too many
% animation frames
W = sparse(weights(:,1)+1,weights(:,2)+1,weights(:,3));
W = full(W);%W = W./repmat(sum(W,2),[1 size(W,2)]);

W = W(:,wmap+wzero);
W(:,wzero)=0;

% Currently there are some rows in W which have up to 6 nonzero weights.
% The industry accepted number is 4 so should probably zero the smaller
% ones.

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

%Display the mesh

cla
for displayWeightIndex=find(~wzero)
    h = trisurf(faces+1,V0(1,:),V0(2,:),V0(3,:),W(:,displayWeightIndex));
    axis equal; axis vis3d;
    shading interp;
    set(h,'FaceAlpha',0.5)
    axis([-200 200 -100 100 -100 500]);
    pause(0.2);
end
hold on

%% -------------------------------------------------------------------------
% Display the skinned armature animation

for j=2:258

    V=V0*0;

    cla
    for i=1:length(sk)
        % -------------------------------------------------------------------------
        % Display the bones
        if sk(i).parent
            s=sk(i).t_xyz(:,j); sp = sk(sk(i).parent).t_xyz(:,j);
            plot3([s(1) sp(1)],[s(2) sp(2)],[s(3) sp(3)],'r','LineWidth',5);
        end
        plot3(sk(i).t_xyz(1,j),sk(i).t_xyz(2,j),sk(i).t_xyz(3,j),'g.','MarkerSize',10);
        if sk(i).Nchannels
            plotax(sk(i).t_xyz(:,j)',sk(i).T(1:3,1:3,j)*sk(i).R0,15)
        end
        % -------------------------------------------------------------------------

        if sk(i).Nchannels
            % Takes a local bone coordinate and puts it into the new pose global position
            poseMatrix(:,:,i,j) = [(sk(i).T(1:3,1:3,j)*sk(i).R0) sk(i).t_xyz(:,j);0 0 0 1];
            % Takes a local bone coordinate and puts it into the base pose global position
            restMatrix(:,:,i) = [((sk(i).R0)) sk(i).head0 ;0 0 0 1];

            % The matrix T for a bone takes a point from global bind space to its new
            % position also in global space as dictated by that bone

            % -------------------------------------------------------------------------------
            % To go from base pose to new pose position we
            % undo the restmatrix and then apply the posematrix i.e.
            % if we have a vertex v from the bind pose then to see where
            % it would be sent by bone i for pose j we would apply
            % vj = M(:,:,j)*[v;1]

            M(:,:,j) = poseMatrix(:,:,i,j)*inv(restMatrix(:,:,i));
            R = M(1:3,1:3,j);
            t = M(1:3,4,j);

            % We add the influence from every bone
            V = V + (R*V0 + repmat(t,1,length(V0))) * spdiags(W(:,i),0,length(V),length(V));

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

        end
    end

    % -------------------------------------------------------------------------
    % Display nicely
    h = trisurf(faces+1,V(1,:),V(2,:),V(3,:),W(:,2)*0);
    axis off;
    axis equal
    shading interp;
    light
    lighting phong
    set(h,'FaceAlpha',0.8)
    drawnow

end
%%

% Code to extract specific joint angles
% Bone i, animation frame j
i = 20; j=2;
T = poseMatrix(1:3,1:3,i,j); Tp = poseMatrix(1:3,1:3,i-1,j); T_0 = restMatrix(1:3,1:3,i); Tp_0 = restMatrix(1:3,1:3,i-1);
R = (T_0')*(Tp_0)*(Tp')*(T);
[Rx Ry Rz] = dcm2angle(R,'XYZ');
Rx = rad2deg(-Rx); Ry = rad2deg(-Ry); Rz = rad2deg(-Rz);

[Rx Ry Rz]


text:
% Test file for SPINPREFSPHERE:

function pass = test_spinprefsphere()

% Construction from inputs:
pref = spinprefsphere('Clim', [0 10], 'dataplot', 'abs', 'dealias', 'off');
pass(1) = isequal(pref.Clim, [0 10]);
pass(2) = strcmpi(pref.dataplot, 'abs');
pass(3) = strcmpi(pref.dealias, 'off');

pref = spinprefsphere('iterplot', 10);
pass(4) = isequal(pref.iterplot, 10);

pref = spinprefsphere('Nplot', 2, 'plot', 'movie');
pass(5) = isequal(pref.Nplot, 2);
pass(6) = strcmpi(pref.plot, 'movie');

pref = spinprefsphere('view', [10 20]);
pass(7) = isequal(pref.view, [10 20]);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% load motor seeds from VAR, save into file and save into GUIdata for each
% fish % careful!!! overwrite!!!
hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);

%%
range_fish = GetFishRange;%[1:3,5:18];
for i_fish = 9:18 %range_fish,
    %% load
    save_masterdir = GetCurrentDataDir();
    save_dir = fullfile(save_masterdir,['subject_' num2str(i_fish)]);
    % load 'data'
    load(fullfile(save_dir,'data_full.mat'),'data'); % struct with many fields

    %% get full-length motorseed regressor
%     ClusterIDs = [11,1];
%     [~,stimrange] = GetStimRange('F',i_fish); % full range
%     setappdata(hfig,'isMotorseed',0);
%     [cIX,gIX,M] = LoadSingleFishDefault(i_fish,hfig,ClusterIDs,stimrange);
%     Behavior_full_motorseed = FindClustermeans(gIX,M);
    % need to do this in raw time instead! (note that full range is not the
    % same sequence as CellResp for fishset=2)
    LoadFullFish(hfig,i_fish);
    ClusterIDs = [11,1];
    [cIX,gIX] = LoadCluster_Direct(i_fish,ClusterIDs(1),ClusterIDs(2));

    CellRespZ = getappdata(hfig,'CellRespZ');
    CellRespAvrZ = getappdata(hfig,'CellRespAvrZ');

    Behavior_full_motorseed = FindClustermeans(gIX,CellRespZ(cIX,:));
    BehaviorAvr_motorseed = FindClustermeans(gIX,CellRespAvrZ(cIX,:));

    %% compute stim-rep averaged motor regressor

%     if length(data.timelists_names)==1, % fish with single stimulus set
%         period = data.periods;
%         BehaviorAvr_motorseed = mean(reshape(Behavior_full_motorseed,size(Behavior_full_motorseed,1),period,[]),3);
%
%     else % fish with multiple stimulus sets
%         nTypes = length(data.timelists_names);
%         periods = data.periods;
%         BehaviorAvr_motorseed = [];
%         for i = 1:nTypes,
%             M_behav = Behavior_full_motorseed(:,data.timelists{i});
%             if mod(numel(M_behav),size(M_behav,1)*periods(i))==0,
%                 avr = mean(reshape(M_behav,size(M_behav,1),periods(i),[]),3);
%             else % for patterns with dummy periods where the stimulus is constant, like 'spont'
%                 nrep = floor(numel(M_behav)/(size(M_behav,1)*periods(i)));
%                 M2_behav = M_behav(:,1:nrep*periods(i));
%                 avr = mean(reshape(M2_behav,size(M_behav,1),periods(i),[]),3);
%             end
%             BehaviorAvr_motorseed = horzcat(BehaviorAvr_motorseed,avr); %#ok<AGROW>
%         end
%     end

    %% save into GUIdata % careful!!! overwrite!!!
    data.Behavior_full_motorseed = Behavior_full_motorseed;
    data.BehaviorAvr_motorseed = BehaviorAvr_motorseed;
    save(fullfile(save_dir,'data_full.mat'),'data');
end

author: xiuyechen

repo: FishExplorer


text:
% [covmatrix] = covfunc(theta)
% [covmatrix, theta] = rand_theta(theta, covmatrix, Y, covfunc,
%                                 logposterior)
% [rand_y] = get_randfunc_y(covmatrix)
% Y = rand_y(covmatrix, Y, I)
% [y, dy] = logprior_theta(theta)
% [y, dy] = loglikelihood_theta(Y, covmatrix)

% COVMATRIX:
% L1/LD1
% L2/LD2
% linsolve1
% linsolve2
% K1
% K2
% logdet1
% logdet2

function res = test_gp_kron(N1,N2,N_samples,seed)

%
% DATA
%

%randn('state', 10)
%rand('state', 10)

if nargin == 0
  N1 = 200;
  N2 = 300;
end
if nargin < 3
  N_samples = 100;
end
if nargin >= 4
  randn('state', seed);
  rand('state', seed);
else
  randn('state', 10);
  rand('state', 10);
end

% Inputs
x1 = 1:N1;
x2 = 1:N2;

% Covariance models

d1 = sqrt(sq_dist(x1(1), x1));
covfunc1 = gp_cov_pp(d1,1);
covfunc1 = gp_cov_toeplitz(covfunc1);
theta1 = 10;
K1 = covfunc1(theta1);

d2 = sqrt(sq_dist(x2(1), x2));
covfunc2 = gp_cov_pp(d2,1);
covfunc2 = gp_cov_toeplitz(covfunc2);
theta2 = 10;
K2 = covfunc2(theta2);

% $$$ nz1 = nnz(K1)
% $$$ nz2 = nnz(K2)
% $$$ K = kron(K1,K2);
% $$$ whos
% $$$ t = cputime();
% $$$   LD = ldlchol(K);
% $$$   ldlsolve(LD,randn(N1*N2,1));
% $$$ time = cputime() - t
% $$$ return

% $$$ % Plot gradients
% $$$ [K1, dK1] = covfunc1(theta1);
% $$$ figure
% $$$ subplot(3,1,1)
% $$$ imagesc(K1)
% $$$ subplot(3,1,2)
% $$$ imagesc(dK1{1})
% $$$ subplot(3,1,3)
% $$$ imagesc(dK1{2})
% $$$ return

disp('Generate data..')

% Generate noisy data
a = 1;
s = 1;
Y = a * kronprod(lchol(K1), lchol(K2), randn(N2,N1));
% $$$ Y = a * (lchol(K2) * randn(N2,N1) * lchol(K1)');
save('Y_noiseless', 'Y');
Y = Y + s*randn(N2,N1);
save('Y_noisy', 'Y');

% Missing values
pmv = 0.5;
%Imv = randperm(N2*N1);
%Imv = Imv(1:floor(pmv*(N2*N1)));
Imv = rand(N2,N1) < pmv;
%Imv = (rand(N2,N1) < 0.2);
ind1 = ceil(N1/2) + (1:(2*theta1));
ind2 = ceil(N2/2) + (1:(2*theta2));
Imv(ind2,ind1) = true;
Y(Imv) = nan;

%
% INFERENCE
%

a = 1e-3;
b = 1e-3;
logprior_theta = @(theta) sum(gamma_logpdf(theta, a, b));
dlogprior_theta = @(theta) gamma_dlogpdf(theta, a, b);

% Gibbs sampling for Y(missing) and covariance parameters
burnin = floor(N_samples/2);

% Initial guess for covariance parameters
theta_init = [2.0        ... % total magnitude
              2.0*theta1 ... % 1) length scale
              0.5*theta2 ... % 2) length scale
              2.0]';         % noise magnitude

res = gp_kron(Y, covfunc1, covfunc2, N_samples, theta_init, ...
              logprior_theta, dlogprior_theta, burnin)

save(sprintf('/home/jluttine/matlab/gp/results_test_gp_kron_%s', ...
             datestr(now,'yyyymmdd')), ...
     '-struct', 'res');

cl = max(abs(Y(:))) + 0.1;
clim = [-cl cl];

Y(Imv) = cl;

figure()
clf

subplot(2,3,1)
imagesc(Y,clim)

load('Y_noiseless')

subplot(2,3,4)
imagesc(Y,clim)
title('true')

subplot(2,3,5)
imagesc(res.F,clim)
title('mean')

subplot(2,3,6)
F_var = sqrt(res.FF - res.F.*res.F);
clim_var = [0, max(F_var(:))];
imagesc(F_var,clim_var)
title('std')

map_colormap();
cm = colormap();
cm(end,:) = [0.5 0.5 0.5];
colormap(cm);

rmse_F = rmse(Y,res.F)

subplot(2,3,2)
semilogy(res.theta');

subplot(2,3,3)
%lag = 200;
plot(acorr(res.theta(:,burnin:end)'));

figure()
plot_scatterhist(res.theta(:,burnin:end)');

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function vw = getCurView
%Returns the selected mrVista view.
%
%    vw = getCurView;
%
% If it can't find a selected one, it initializes a hidden inplane. A sort
% of super-function to getSelectedInplane, getSelectedGray, etc.
%
% How it looks: first, it checks if any views are open in the global
% INPLANE, VOLUME, or FLAT variables. If not, opens a hidden inplane.
% If there are open views, it checks if the current figure is the handle to
% one of the views. If this fails, then tries to find the selected inplane;
% then the selected volume; then the selected flat.
%
% ras, 05/01/06.
mrGlobals;

vw = [];

if isempty(cellfind(INPLANE)) && isempty(cellfind(VOLUME)) && ...
        isempty(cellfind(FLAT))
    disp('No open mrVista UIs. Initializing a hiddden inplane...')
    try
        vw = initHiddenInplane;
    catch
        disp('Couldn''t start a hidden inplane -- not in a mrVista session!')
    end
    return
end

% check if gcf is a view handle
for i = cellfind(INPLANE)
    if checkfields(INPLANE{i}, 'ui', 'windowHandle') && ...
            isequal(INPLANE{i}.ui.windowHandle, gcf)
        vw = INPLANE{i};
        return
    end
end
for i = cellfind(VOLUME)
    if checkfields(VOLUME{i}, 'ui', 'windowHandle') && ...
            isequal(VOLUME{i}.ui.windowHandle, gcf)
        vw = VOLUME{i};
        return
    end
end
for i = cellfind(FLAT)
    if checkfields(FLAT{i}, 'ui', 'windowHandle') && ...
            isequal(FLAT{i}.ui.windowHandle, gcf)
        vw = FLAT{i};
        return
    end
end

% None of the above worked; try to get selected inplane/volume/flat
if isempty(vw), vw = getSelectedInplane; end
if isempty(vw), vw = getSelectedVolume;  end
if isempty(vw), vw = getSelectedFlat;    end

% Report on failure
if isempty(vw)
    disp('Failed to find any views ... returning empty')
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function structNum = getStructNum(str,planC,indexS)
%Get the structure number for a structure with a given name.
%Returns 0 if not found.
%14 Oct 02, JOD
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

structNum = 0;
for i = 1 : length(planC{indexS.structures})
  if strcmp(str,planC{indexS.structures}(i).structureName)
    structNum = i;
    return
  end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR