text:
function test_suite = test_clipLine3d
%testClipLine3d  One-line description here, please.
%   output = testClipLine3d(input)
%
%   Example
%   testClipLine3d
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-06-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testOx(testCase) %#ok<*DEFNU>
% line parallel to Ox axis
box     = [0 100 0 100 0 100];
line    = [10 20 30 10 0 0];

clipped = clipLine3d(line, box);
edge    = [0 20 30 100 20 30];
testCase.assertEqual(edge, clipped, 'AbsTol', .001);

function testOx_outside(testCase)

% line parallel to Ox axis
box     = [0 100 0 100 0 100];

line    = [10 -20 30 10 0 0];
edge    = NaN(1, 6);
clipped = clipLine3d(line, box);
testCase.assertEqual(edge, clipped, 'AbsTol', .001);

line    = [-10 120 30 10 0 0];
edge    = NaN(1, 6);
clipped = clipLine3d(line, box);
testCase.assertEqual(edge, clipped, 'AbsTol', .001);

line    = [10 20 -30 10 0 0];
edge    = NaN(1, 6);
clipped = clipLine3d(line, box);
testCase.assertEqual(edge, clipped, 'AbsTol', .001);

line    = [10 20 130 10 0 0];
edge    = NaN(1, 6);
clipped = clipLine3d(line, box);
testCase.assertEqual(edge, clipped, 'AbsTol', .001);

function testOy(testCase)
% line parallel to Ox axis
box     = [0 100 0 100 0 100];
line    = [10 20 30 0 10 0];

clipped = clipLine3d(line, box);
edge    = [10 0 30 10 100 30];
testCase.assertEqual(edge, clipped, 'AbsTol', .001);

function testOy_outside(testCase)

% line parallel to Ox axis
box     = [0 100 0 100 0 100];

line    = [-10 20 30 0 10 0];
edge    = NaN(1, 6);
clipped = clipLine3d(line, box);
testCase.assertEqual(edge, clipped, 'AbsTol', .001);

line    = [110 20 30 0 10 0];
edge    = NaN(1, 6);
clipped = clipLine3d(line, box);
testCase.assertEqual(edge, clipped, 'AbsTol', .001);

line    = [10 20 -30 0 10 0];
edge    = NaN(1, 6);
clipped = clipLine3d(line, box);
testCase.assertEqual(edge, clipped, 'AbsTol', .001);

line    = [10 20 130 0 10 0];
edge    = NaN(1, 6);
clipped = clipLine3d(line, box);
testCase.assertEqual(edge, clipped, 'AbsTol', .001);

function testOz(testCase)
% line parallel to Ox axis
box     = [0 100 0 100 0 100];
line    = [10 20 30 0 0 10];

clipped = clipLine3d(line, box);
edge    = [10 20 0 10 20 100];
testCase.assertEqual(edge, clipped, 'AbsTol', .001);

function testOz_outside(testCase)
% line parallel to Ox axis
box     = [0 100 0 100 0 100];

line    = [-10 20 30 0 0 10];
edge    = NaN(1, 6);
clipped = clipLine3d(line, box);
testCase.assertEqual(edge, clipped, 'AbsTol', .001);

line    = [110 20 30 0 0 10];
edge    = NaN(1, 6);
clipped = clipLine3d(line, box);
testCase.assertEqual(edge, clipped, 'AbsTol', .001);

line    = [10 -20 30 0 0 10];
edge    = NaN(1, 6);
clipped = clipLine3d(line, box);
testCase.assertEqual(edge, clipped, 'AbsTol', .001);

line    = [10 120 30 0 0 10];
edge    = NaN(1, 6);
clipped = clipLine3d(line, box);
testCase.assertEqual(edge, clipped, 'AbsTol', .001);

function testArray(testCase)

% test for several lines with multiple directions
box     = [0 100 0 100 0 100];
lineOx  = [10 20 30 10 0 0];
lineOy  = [10 20 30 0 10 0];
lineOz  = [10 20 30 0 0 10];
lines   = [lineOx;lineOy;lineOz];
clipped = clipLine3d(lines, box);

testCase.assertEqual(3, size(clipped, 1));

% the same, but with some lines outside the box
box     = [0 100 0 100 0 100];
lineOx1 = [ 10  20  30 10 0 0];
lineOx2 = [ 10 -20  30 10 0 0];
lineOx3 = [ 10  20 -30 10 0 0];
lineOy1 = [ 10  20  30 0 10 0];
lineOy2 = [-10  20  30 0 10 0];
lineOy3 = [ 10  20 -30 0 10 0];
lineOz1 = [ 10  20  30 0 0 10];
lineOz2 = [-10  20  30 0 0 10];
lineOz3 = [ 10 -20  30 0 0 10];
lines   = [...
    lineOx1;lineOx2;lineOx3;...
    lineOy1;lineOy2;lineOy3;
    lineOz1;lineOz2;lineOz3];
clipped = clipLine3d(lines, box);

testCase.assertEqual(9, size(clipped, 1));
testCase.assertEqual(6, size(clipped, 2));
testCase.assertEqual(6*6, sum(isnan(clipped(:))));

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function header = read_spmeeg_header(filename)

% read_spmeeg_header() - import SPM5 and SPM8 meeg datasets
%
% Usage:
%   >> header = read_spmeeg_header(filename);
%
% Inputs:
%   filename - [string] file name
%
% Outputs:
%   header   - FILEIO toolbox type structure
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
% Vladimir Litvak

if nargin < 1
  help read_spmeeg_header;
  return;
end

try
    D = load(filename, 'D');
    D = D.D;
catch
    ft_error('File not found or wrong format.');
end

header = [];

if isfield(D, 'Radc') % SPM5
    header.Fs          = D.Radc;
    header.nChans      = D.Nchannels;
    header.nSamples    = D.Nsamples;

    if isfield(D, 'events') && isfield(D.events, 'start')
        header.nSamplesPre = D.events.start;
    else
        header.nSamplesPre = 0;
    end

    header.nTrials     = D.Nevents;
    header.label       = D.channels.name';

elseif all(isfield(D, {'type', 'Nsamples', 'Fsample', 'timeOnset'})) % SPM8

    header.Fs          = D.Fsample;
    header.nChans      = numel(D.channels);
    header.nSamples    = D.Nsamples;
    header.nSamplesPre = -D.timeOnset*D.Fsample;
    header.nTrials     = numel(D.trials);
    header.label       = {D.channels.label}';
    header.chantype    = lower({D.channels.type}');
    header.chanunit    = {D.channels.units}';

    if isfield(D.sensors, 'eeg') && ~isempty(D.sensors.eeg)
        header.elec = D.sensors.eeg;
    end

    if isfield(D.sensors, 'meg') && ~isempty(D.sensors.meg)
        header.grad = D.sensors.meg;
    end

else
    ft_error('Cannot recognize an SPM EEG header format');
end

header.orig = D;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [wheelslam_rms, wheelins_rms] = wheelslam_func(paras)

% clear
% close all
format longG

datapath         = paras.datapath;
rollSeqdimension = paras.rollSeqdimension; %Dimension of the roll sequence measurement
NPARTICLES       = paras.NPARTICLES;
gridlen          = paras.gridlen;
NEFFECTIVE       = paras.NEFFECTIVE;   %Threshold of the effective particles
sigmaDis_scale   = paras.sigmaDis_scale; % unit: m
sigmaPhi         = paras.sigmaPhi;     % unit: radians
initheading      = paras.initheading; % Initial heading of the vehicle
conRevisitNumThr = paras.conRevisitNumThr;% The number threshold for determining the revisit of the particle
corrCoefThr      = paras.corrCoefThr;
corrCoefNumThr   = paras.corrCoefNumThr;
gt_path          = paras.gt_path;

freq             = 10;
odom             = read_bin(datapath, 4)';% time, dis(equal), headingInc, roll

%Import odometry info from Wheel-INS, including timestamp, dis_increment,
%heading increment and roll angle

revisit_tThr   = 30;      %Minimal time interval for revisit identification. Sometimes car stop
initmapsize    = 101;
progrssThr     = -100;
start_t        = odom(1,1);
odomsize       = floor(size(odom, 1)*paras.odomdata_scale);

particles      = init_particles(NPARTICLES, initmapsize, odomsize, start_t, initheading, corrCoefNumThr);
res            = zeros(odomsize, 4);
res(1, :)      = [start_t particles(1).xv'];

pre_refstate      = [start_t particles(1).xv'];
wheelinsDR        = zeros(odomsize, 4);
wheelinsDR(1,:)   = pre_refstate;
revisitParticles  = 0;

% This is used to set the randanstream of the particles, so as to
% make the results the same everytime.
% for i = 1:NPARTICLES
%     particles(i).randstream = RandStream.create('mlfg6331_64','Seed', i);
% end
%
% resample_randstream = RandStream.create('mlfg6331_64','Seed', idx);
%---------------------------------------------------------------------

for i = 2:odomsize

    % Current measurement of roll sequence.
    % curRollmeas = odom(i, 5:size(odom, 2));
    curRollmeas     = odom(i, 4);
    ref_state       = state_predict_ref(odom(i, :), pre_refstate);% Get the DR results from the odometry file for comparision.
    wheelinsDR(i,:) = ref_state;
    pre_refstate    = ref_state;

    t_s             = clock;
    revisitNum      = 1;
    curodom = odom(i, :);

    for j = 1:NPARTICLES

        particles(j)                 = state_predict(particles(j), curodom, sigmaDis_scale, sigmaPhi); %State prediction for every particle.
        [particles(j), atRow, atCol] = extendMap(particles(j), gridlen); %Extend (or not) the grid map maintained by every particle.

        if (particles(j).conRevisitNum ~= conRevisitNumThr)
            particles(j).conRevisitNum = particles(j).conRevisitNum + 1;
        else
            % sliding window for loop closure detection.

            particles(j).conCorrCoef(1:(conRevisitNumThr-1)) = particles(j).conCorrCoef(2:conRevisitNumThr);
            particles(j).conCorrCoef(conRevisitNumThr)       = -1;

        end

        particles(j).conCorrCoef(particles(j).conRevisitNum)    = -1;

        if (atRow == particles(j).preIdxinMap(1) && atCol == particles(j).preIdxinMap(2)...
                && particles(j).lastTrajRevisitIdx < i-1)

            % If one step is too short (still in the same grid) and if last step is not a revisit
            numforgrid = particles(j).mapattr(particles(j).gridmap(atRow, atCol)).Num;
            particles(j).mapattr(particles(j).gridmap(atRow, atCol)).Value...
            = (particles(j).mapattr(particles(j).gridmap(atRow, atCol)).Value ...
            * (numforgrid/(numforgrid+1)) + curRollmeas./(numforgrid +1));

            % Average the roll measurement in the same grid.
            particles(j).mapattr(particles(j).gridmap(atRow, atCol)).Num = numforgrid + 1;

        else

            innerGridsIndex  = [];
            revisitGrids     = [];
            allRollSeq       = [];
            allheading       = [];

            curallheading = [allheading;particles(j).xv(3)];
            curAllRollSeq = [allRollSeq;curRollmeas];
            curGrid       = [atRow, atCol];
            curAllGrids   = [innerGridsIndex; curGrid];
            curAllGrids   = unique(curAllGrids, 'rows','stable');

            if (~isempty(curAllGrids))

                %Update (or not) the grid map maintained by every particle
                %and get the revisit grids after revisit time interval check.
                %Check if current grid is a new grid.
                [particles(j), revisitGrids] = updateGridmap(particles(j), curAllGrids, curAllRollSeq, curallheading, revisit_tThr, i);

            end

            if (~isempty(revisitGrids))

                for ii = 1:size(revisitGrids, 1)
                    %Threshold for the heading difference between current entance and last one
                    headingdiff = abs(getheadingdiff(particles(j).xv(3), particles(j).mapattr(particles(j).gridmap(revisitGrids(ii,1), revisitGrids(ii,2))).enterheading));
                    if(headingdiff > (pi/6))
                        break;
                    end

                    particles(j).mapattr(particles(j).gridmap(revisitGrids(ii,1), revisitGrids(ii,2))).revisit = ...
                    particles(j).mapattr(particles(j).gridmap(revisitGrids(ii,1), revisitGrids(ii,2))).revisit + 1;
                    particles(j).mapattr(particles(j).gridmap(revisitGrids(ii,1), revisitGrids(ii,2))).Num = ...
                    particles(j).mapattr(particles(j).gridmap(revisitGrids(ii,1), revisitGrids(ii,2))).Num + 1;

                    %Get the historical index of the revisited grid in the
                    %roll sequence.
                    revisitIdx = particles(j).mapattr(particles(j).gridmap(revisitGrids(ii,1), revisitGrids(ii,2))).idx;
                    if( revisitIdx < rollSeqdimension)
                        break;
                    else
                        particles(j).lastTrajRevisitIdx = i; %where the particle detected a revisit by the traj.

                        currollseq = odom(i-rollSeqdimension+1:i,4);% Current roll sequence from the revisited grid
                        maprollseq = odom(revisitIdx-rollSeqdimension+1:revisitIdx,4); %Historical roll sequence
                        corrMat    = corrcoef(currollseq, maprollseq);
                        corrvalue  = corrMat(1,2); % Get the correlation value
                        particles(j).conCorrCoef(particles(j).conRevisitNum) = corrvalue;
                        corrCoefNumtmp = length(particles(j).conCorrCoef(particles(j).conCorrCoef> corrCoefThr));

                        if(particles(j).conRevisitNum == conRevisitNumThr && ...
                           corrCoefNumtmp > corrCoefNumThr &&...
                           corrvalue > corrCoefThr)
                           %criteria 1. Full of the correlation window;
                           %2.Enough value in the window larger than the
                           %threshold
                           %3.Current correlation value larger than the
                           %threshold;

                            tmpSeq = particles(j).conCorrCoef((particles(j).conCorrCoef > corrCoefThr ));

                            particles(j).w = particles(j).w * exp(rms(tmpSeq)* (corrCoefNumtmp/conRevisitNumThr)); %Update the weight of the particles reported a convinced loop closure
                            particles(j).mapattr(particles(j).gridmap(revisitGrids(ii,1), revisitGrids(ii,2))).last_vt = particles(j).t;
                            particles(j).totalrevisitNum = particles(j).totalrevisitNum + 1;
                            particles(j).revisitpos(particles(j).totalrevisitNum, :) = [particles(j).t particles(j).xv(1) particles(j).xv(2)];

                        end

                        % count the revisited particle only once
                        if (ii == 1)
                            revisitParticles(revisitNum) = j;
                            revisitNum = revisitNum + 1;
                        end

                    end

                end

            end

        end

        particles(j).preIdxinMap = [atRow, atCol];
        particles(j).prepos      = particles(j).xv;
        particles(j).traj(i,:)   = (particles(j).xv)';

    end

    % Check if resample is needed, if yes, do it. Generally this function
    % needn't to be modified, we can set the parameters "NEFFECTIVE"
    particles    = resampleParticles(particles, NEFFECTIVE);%, resample_randstream
    %particles    = resampleParticles(particles, NEFFECTIVE, resample_randstream);

    %save system output
    w            = [particles.w];
    xv           = [particles.xv];
    if (mean(abs(xv(3,:))) > abs(mean(abs(xv(3,:))) - pi))
        xv(3, :)     = zero_to_2pi(xv(3, :));
    else
        xv(3, :)     = pi_to_pi(xv(3, :));
    end

    ii       = find(w == max(w),1);
    w        = w/sum(w);
    xvmean   = [mean(sum(w.* xv(1, :))) mean(sum(w.* xv(2, :))) mean(sum(w.* xv(3, :)))];
    res(i,:) = [particles(1).t xvmean];

    revisitParticles = 0;

    runProgress  = floor((i-1)/(odomsize) * 1000);
    if(runProgress > progrssThr)
        clc;
        progrssThr = runProgress;
        disp(['******Wheel SLAM RUNNING:' num2str(0.1*runProgress) '%******']);
    end

end

% figure,
% plot3(res(:,3), res(:,2), res(:,1)-res(1,1), 'LineWidth', 1.5),grid on;
% title('Weighted mean traj.');hold on;
% plot(res(1,3), res(1,2), 'go');hold on;
% plot(res(end,3), res(end,2), 'ro');
% set(gca,'fontsize',10,'fontname','Times');

% plotMap(particles(ii),i, path);

% write_bin(resfile, res);
% write_bin(wheelinsDRpath, wheelinsDR);

particleres = [res(:, 1) particles(ii).traj];% trajectory of the particle with highest weight
[wheelslam_err, wheelslam_rms, wheelins_err, wheelins_rms, wheelslam_err_p, wheelslam_rms_p] =...
    ResCompare(res, particleres, wheelinsDR, gt_path, freq);

%% plot the total revisit num of all the particles
% revisitNums           = [particles.totalrevisitNum];
% figure,
% set(gcf,'unit','normalized','position',[0.05,0.05,0.64,0.48]);
% plot(revisitNums,'*'); set(gca,'fontsize',10,'fontname','Times');
% title('Total Revisit Num');
% figure,
% plot3(particles(ii).revisitpos(:,3), particles(ii).revisitpos(:,2), particles(ii).revisitpos(:,1),'*');
% title(['Particle =' num2str(ii) 'idx = ' num2str(idx)]);

% save particle data
% save([path 'particle_' num2str(ii) '_' s_curt], 'particles');

% format longG
% fprintf(fp_res,'%s\n', s_curt);
% fprintf(fp_res,'Particles = %d NEFFECTIVE = %d\n', NPARTICLES, NEFFECTIVE);
% fprintf(fp_res,'Gridlen = %.1f mode = %d\n', gridlen, mode);
% fprintf(fp_res,'Disstd    = %.5f  Phistd  = %.5f \n', sigmaDis_scale, sigmaPhi);
% fprintf(fp_res,'RollSampleDis = %.1f    RollSeqDim = %d\n', paras.sampleDis, rollSeqdimension);
% fprintf(fp_res,'RollSeqWindow = %d  CorrCoefThr = %.1f  CorrCoefNumThr = %d\n',conRevisitNumThr,corrCoefThr,corrCoefNumThr);
% fprintf(fp_res,'Start_t = %f  TotalData_t = %f  Process_t = %s s\n',start_t, totaldata_t, processtime);
% fprintf(fp_res,'%s\n', 'RMS        [North     East]');
% fprintf(fp_res,'WheelSLAM    %.5f  %.5f\n', wheelslam_rms(1), wheelslam_rms(2));
% fprintf(fp_res,'WheelSLAM_p  %.5f  %.5f\n', wheelslam_rms_p(1), wheelslam_rms_p(2));
% fprintf(fp_res,'WheelINS     %.5f  %.5f\n', wheelins_rms(1), wheelins_rms(2));
% fprintf(fp_res,'Heading RMSE:  Wheel-SLAM:%.5f Wheel-INS:%.5f\n\n', wheelslam_rms(3), wheelins_rms(3));

end

function p = init_particles(np, mapsize, odomsize, start_t, initheading, corrCoefNumThr)

    %Grid map attribute
    gridattr.Num     = 1; %how many times has this grid been visited
    gridattr.visited = 1; %if this grid has been visited
    gridattr.revisit = 0; %if this grid has been revisited
    gridattr.last_vt = start_t;%last visit time of the grid

    gridattr.Value   = zeros(1, 1);%only single roll angle
    gridattr.idx     = 1;% the index of the current roll measurement in the full odometry file
    initheadingrad   = DEG2RAD(initheading);
    gridattr.enterheading = initheadingrad;%the heading of the vehicle when it entering this grid

    p.conRevisitNum  = 0;%length of the correlation sliding window
    p.conCorrCoef    = zeros(1,corrCoefNumThr);
    p.gridmap        = zeros(mapsize, mapsize);% save the idx in the grid attibute matrix at every grid of the map
    %the coordinates of the original point is (0,0), but it is at the
    %central of the map
    p.originIdx      = [(mapsize+1)/2, (mapsize+1)/2];% the idx of the original point in the grid map (which would be updated with the extension of the map).
    p.w              = 1/np;          %weight
    p.xv             = [0;0;initheadingrad]; %state: N, E, Phi

    p.prepos         = [0;0;initheadingrad]; %last position and heading of the vehicle
    p.preIdxinMap    = p.originIdx; %last idx of the vehicle in the grid map
    p.t              = start_t;

    p.totalrevisitNum = 0; %how many convinced loop closure proposed by this particle
    p.revisitpos      = zeros(1,3);%The positions of the vehicle where the convinced loop closure proposed by this particle

    p.mapattr(1)     = gridattr;% attibute of the first grid in the map. A large grid attribute matrix would be built with the moving of the vehicle
    p.totalGrids     = 1;%
    p.gridmap(p.originIdx(1), p.originIdx(2)) = 1;% grid map maintaind by every particle
    p.traj           = zeros(odomsize,3);% the trajectory of the vehicle estimated by this particle
    p.traj(1,:)      = (p.xv)';
    p.lastTrajRevisitIdx = 0;%The Idx in the odom file where the particle detects a revisit by its traj, (no roll needed)
    p                = repmat(p, [np 1]);

end

function x_vec = zero_to_2pi(x_vec)
    lenth = length(x_vec);
    for i = 1:lenth
        if(x_vec(i)<0)
            x_vec(i) = x_vec(i) + 2*pi;
        end
    end

end

function x = pi_to_pi(x)
    if x > pi
        x= x - 2*pi;
    elseif x < -pi
        x= x + 2*pi;
    end
end

function hdiff = getheadingdiff(heading1, heading2)
    hdiff = heading1 - heading2;
    if hdiff > pi
        hdiff= hdiff - 2*pi;
    elseif hdiff < -pi
        hdiff= hdiff + 2*pi;
    end
end

function [bestIdx, lagestcorrValue] = getBestMacth(currollseq, revisitIdx, searchRad, odom, rollSeqdimension)

    lagestcorrValue = -exp(5);
    for i = (revisitIdx-searchRad):(revisitIdx+searchRad)
        maprollseq = odom(i-rollSeqdimension+1:i, 4);
        corrMat    = corrcoef(currollseq, maprollseq);
        corrvalue  = (corrMat(1,2) + corrMat(2, 1))/2;
        if corrvalue > lagestcorrValue
            lagestcorrValue = corrvalue;
            bestIdx = i;
        end
    end
end

repo: Wheel-SLAM


text:
function [App,Dis] = subFitness(PopObj,P,Center,R)
% Update intersection point solution in each subregion

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    K          = length(R);
    [N,M]      = size(PopObj);

    % Normalize the population
    fmin       = min(PopObj,[],1);
    fmax       = max(PopObj,[],1);
    Obj        = (PopObj-repmat(fmin,N,1))./repmat(fmax-fmin,N,1);

    %% Calculate intersection point in each subregion
    if K == 1
        InterPoint = interPoint(Obj,P);
    else
        InterPoint = ones(N,M);

        % Allocation
        transformation = Allocation(Obj,Center,R);

        for i = 1 : K
            current     = find(transformation == i);
            if ~isempty(current)
                sInterPoint = interPoint(Obj(current,:),P(i,:));
                InterPoint(current,:) = sInterPoint;
            end
        end
    end

    % Calculate the diversity and convergence of intersection points
    App = min(InterPoint-Obj,[],2);
%     App = sqrt(sum(InterPoint.^2,2)) - sqrt(sum(Obj.^2,2));

    % Calculate the diversity of each solution
%     Dis = pdist2(InterPoint,InterPoint);
    Dis = distMax(InterPoint);
    Dis(logical(eye(length(Dis)))) = inf;
end

function InterPoint = interPoint(PopObj,P)
% Calcualte the approximation degree of each solution, and the distances
% between the intersection points of the solutions

    [N,~] = size(PopObj);

    %% Calculate the intersections by gradient descent
    P     = repmat(P,N,1);      % Powers
    r     = ones(N,1);          % Parameters to be optimized
    lamda = zeros(N,1) + 0.002;   % Learning rates
    E     = sum((r.*PopObj).^P,2) - 1;   % errors
    for i = 1 : 1000
        newr = r - lamda.*E.*sum(P.*PopObj.^P.*r.^(P-1),2);
        newE = sum((newr.*PopObj).^P,2) - 1;
        update         = newr > 0 &sum(newE.^2) < sum(E.^2);
        r(update)      = newr(update);
        E(update)      = newE(update);
        lamda(update)  = lamda(update)*1.002;
        lamda(~update) = lamda(~update)/1.002;
    end
    InterPoint = PopObj.*r;
end

function Dis = distMax(X)
% distMax pairwise distance between one set of observations
% Dis = distMax(X) returns a matrix D containing the maximum absolute
%   distance per dimension between each pair of observations in the MX-by-N
%   data matrix X and MX-by-N data matrix X.

%   Example:
%      X = randn(100, 5);
%      D = distMax(X,Y);
%   >>size(D) = 100*100

    if isempty(X)
        error('X must be a non-empty matrix');
    end

    [N,~] = size(X); % nx,p
    Dis = zeros(N,N);
    for i = 1 : N
        for j = i+1 : N
            Dis(i,j) = max(abs(X(i,:)-X(j,:)));
        end
    end
    Dis = Dis + Dis';
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function pass = test_feval(pref)
% Test CHEBFUN3/FEVAL

if ( nargin == 0)
    pref = chebfunpref;
end
tol = 100*pref.cheb3Prefs.chebfun3eps;

seedRNG(42);

f = chebfun3t(@(x,y,z) x, [-1 2 -pi/2 pi -3 1]);
pass(1) = abs(f(0,0,0)) < tol*f.vscale;

pass(2) = abs(f(pi/6,pi/12,-1)-pi/6) < tol*f.vscale;

f = chebfun3t(@(x,y,z) y, [-1 2 -pi/2 pi -3 1]);
pass(3) = abs(f(0,0,0)) < tol;

pass(4) = abs(f(pi/6,pi/12,-1)-pi/12) < tol*f.vscale;

f = chebfun3t(@(x,y,z) z, [-1 2 -pi/2 pi -3 1]);
pass(5) = abs(f(0,0,0)) < tol;

pass(6) = abs(f(pi/6,pi/12,-1)+1) < tol*f.vscale;

% some harder tests.
f = @(x,y,z) cos(x) + sin(x.*y) + sin(z.*x);
g = chebfun3t(f);
pts = 2*rand(3,1) - 1;
pass(7) = abs(f(pts(1),pts(2),pts(3)) - g(pts(1),pts(2),pts(3))) < ...
    tol*g.vscale;

% Are we evaluating on arrays correctly?
r = rand(10,1);
s = rand(10,1);
t = rand(10,1);
[rr, ss, tt]=meshgrid(r,s,t);
pass(8) = max(abs((f(r,s,t) - g(r,s,t)))) < tol*g.vscale;

pass(9) = max(max(max(abs(f(rr,ss,tt) - g(rr,ss,tt))))) < tol*g.vscale;

% Does this work off [-1,1]^2
g = chebfun3t(f,[-pi/6 pi/2 -pi/12 sqrt(3) -3 1]); % strange domain.
r = 0.126986816293506; s = 0.632359246225410; t = 0.351283361405006;
% three fixed random number in domain.
pass(10) = abs(f(r,s,t) - g(r,s,t))<tol*g.vscale;

% Are we evaluating on arrays correctly
pass(11) = abs(f(r,s,t) - g(r,s,t)) < tol*g.vscale;

pass(12) = max(max(max(abs(f(rr,ss,tt) - g(rr,ss,tt))))) < tol*g.vscale;

% vector inputs
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3t(ff);
xx = linspace(-1, 1, 100)';
yy = linspace(-1, 1, 100)';
zz = linspace(-1, 1, 100)';
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(13) = norm(F(:) - FF(:)) < 100*tol;

% random vector inputs
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3t(ff);
xx = rand(100, 1);
yy = rand(100, 1);
zz = rand(100, 1);
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(14) = norm(F - FF) < 100*tol;

% meshgrid
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3t(ff);
[xx, yy, zz] = meshgrid(linspace(-1, 1, 100));
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(15) = norm(F(:) - FF(:)) < 100*tol;

% ndgrid
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3t(ff);
[xx, yy, zz] = ndgrid(linspace(-1, 1, 100));
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(16) = norm(F(:) - FF(:)) < 100*tol;

% random tensor inputs
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3t(ff);
xx = rand(10, 20, 30);
yy = rand(10, 20, 30);
zz = rand(10, 20, 30);
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(17) = norm(F(:) - FF(:)) < 100*tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ vX ] = SolveBpAdmm( mA, vB, paramLambda )
% ----------------------------------------------------------------------------------------------- %
%[ vX ] = SolveBpAdmm( mA, vB, paramLambda )
% Solve Basis Pursuit (Q1) problem using ADMM.
% Input:
%   - mA                -   Input Matirx.
%                           The model matrix.
%                           Structure: Matrix (m X n).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vB                -   Input Vector.
%                           The model known data.
%                           Structure: Vector (m X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - paramLambda       -   Parameter Lambda.
%                           Sets the balance between L1 minimization and
%                           Least Squares minimization.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: (0, inf).
% Output:
%   - vX                -   Output Vector.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  A
% Remarks:
%   1.  A
% Known Issues:
%   1.  A
% TODO:
%   1.  B
% Release Notes:
%   -   1.0.000     03/04/2018
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

numIterations = 350;
primalEps = 1e-6; %<! Stopping Condition

numRows = size(mA, 1);
numCols = size(mA, 2);

vX = zeros([numCols, 1]);
vU = zeros([numCols, 1]);
vV = zeros([numCols, 1]);

mAA = mA.' * mA;
vAb = mA.' * vB;
mI = eye(numCols);

mAAI = mAA + mI;

mL = chol(mAAI, 'lower');
mU = mL.';

sL.LT = true();
sU.UT = true();

for ii = 1:numIterations
    % vX = mAAI \ (vAb + vV - vU);
    vX = linsolve(mU, linsolve(mL, vAb + vV - vU, sL), sU); %<! https://www.mathworks.com/help/matlab/ref/linsolve.html
    vV = SoftThresholding(vX + vU, paramLambda);
    vU = vU + vX - vV;

%     if(sqrt(mean((vX - vV) .^ 2)) < primalEps) %<! Primal Convergence
%         break;
%     end
end

end

function [ vX ] = SoftThresholding( vX, lambdaFactor )

% Soft Thresholding
vX = max(vX - lambdaFactor, 0) + min(vX + lambdaFactor, 0);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function [TRI4, x4, y4, z4]=make_sphere(TRI, x, y, z)
% % Syntax;
% %
% % [TRI3, x3, y3, z3]=make_sphere(TRI, x, y, z);
% %
% % ***********************************************************
% %
% % Description
% %
% % Program makes the dome of an even frequency icosahedron
% % into a full spherical icosahedron.
% %
% % ***********************************************************
% %
% % Input Variables
% %
% % x, y,and z are row vector of the x, y, and z coordinates of the
% % icosahedron nodes of the dome.
% %
% % TRI is the triangularization of the icosahedron nodes.
% %
% % ***********************************************************
% %
% % Output Variables
% %
% % x4,y4,and z4 are row vector of the x4, y4, and z4 coordinates of the
% % full spherical icosahedron nodes.
% %
% % TRI4 is the triangularization of the spherical icosahedron nodes.
% %
% % ***********************************************************
% %
% Example
%
% f=2; % 2 freuqency similar to soccer ball
%
% r=4; % 4 meter radius
%
% [x2, y2, z2, it]=icosahedron_nodes(f, r);
%
% [x, y, z]=build_icosa(x2, y2, z2);
%
% [rho, theta, phi]=spherical_angle_ed(x, y, z);
%
% [x10, y10, z10]=splat(rho, theta, phi);
%
% TRI = delaunay(x10, y10);
%
% [TRI, x, y, z]=make_sphere(TRI, x, y, z);
%
% figure(1);
% h=trisurf(TRI,x,y,z, 'FaceColor', [1 0 0], 'EdgeColor', 0*[1 1 1], ...
% 'LineWidth', 1 );
% title('Red Soccer Ball');
% axis equal
%
% % ***********************************************************
% %
% % This program was written by Edward L. Zechmann
% %
% % date 7   March   2008
% %
% % modified 11   March   2008  added examples and comments
% %
% % ***********************************************************
% %
% % Feel free to modify this code.
% %

nn=length(x);
r=sqrt((x(1, 1).^2+y(1, 1).^2+z(1, 1).^2));

% Calculate the frequency from the length of x
% Assuming that x,y,and z are for a dome.
% length(x)=5*f^2+5/2*f+1;
% using the quadratic equation
f=round(0.25*(sqrt(1/5*(16*length(x)-11))-1));

% For a full sphere the freuqency can be caculated from teh number
% of nodes by the formulas
% length(x)=10*f^2+2;
% Using the quadratic equation
% f=round(sqrt((length(x)-2)/10));

ix1=find( abs(z) < 0.01*r/(f+1));
ix2=setdiff(1:nn, ix1);

trix2=[];
trix3=[];
m=length(TRI);

for e1=1:m;

   tri1=TRI(e1, :);

    aa=intersect(tri1,ix1);

    if ~isempty(aa)
        trix2=[trix2 e1];
    else
        trix3=[trix3 e1];
    end

end

TRI2=TRI(trix2, :);
TRI3=TRI(trix3, :);
TRI33=zeros(size(TRI3));

for e1=1:length(TRI3);
    for e2=1:3;
        TRI33(e1, e2)=find(TRI3(e1, e2)==ix2);
    end
end

TRI22=zeros(size(TRI2));

for e1=1:length(TRI2);
    for e2=1:3;

        ix3=find(TRI2(e1, e2)==ix1);
        if ~isempty(ix3)
            TRI22(e1, e2)=ix1(ix3);
        else
            ix4=nn+find(TRI2(e1, e2)==ix2);
            TRI22(e1, e2)=ix4;
        end

    end
end

x4=[x x(ix2)];
y4=[y y(ix2)];
z4=[z -z(ix2)];

TRI4=[TRI' [nn+TRI33]' TRI22']';

author: Sable

repo: mcbench-benchmarks


text:
function pass = test_chebmatrix

% TODO: Tests 1 and 2 assume a chebcolloc2 discretization.

pref = cheboppref();
pref.discretization = @chebcolloc2;

%% Building blocks
dom = [-2 -0.5 1 2];
I = operatorBlock.eye(dom);
D = operatorBlock.diff(dom);
Z = operatorBlock.zeros(dom);
x = chebfun('x', dom);
u = sin(x.^2);
U = operatorBlock.mult(u);

D5 = [
  -5.499999999999999   6.828427124746189  -2.000000000000000   1.171572875253810  -0.500000000000000
  -1.707106781186547   0.707106781186547   1.414213562373095  -0.707106781186548   0.292893218813452
   0.500000000000000  -1.414213562373095                   0   1.414213562373095  -0.500000000000000
  -0.292893218813452   0.707106781186548  -1.414213562373095  -0.707106781186547   1.707106781186547
   0.500000000000000  -1.171572875253810   2.000000000000000  -6.828427124746189   5.499999999999999
];

%%
A = [ I,Z; D,U ];
M = matrix(A, [5 5 5], pref);
DD = blkdiag(2/1.5*D5,2/1.5*D5,2/1*D5);
[xx, ww] = chebpts([5 5 5], dom);
UU = diag(u(xx));

err(1) = norm( M - [ eye(15), zeros(15); DD, UU ]);

%% more complicated chebmatrix
A = [ I, x, -3*I;
    functionalBlock.sum(dom), 5, functionalBlock.feval(dom(end),dom);
    D, chebfun(1,dom), U ];
M = matrix(A, [5 5 5], pref);
MM = [ eye(15), xx, -3*eye(15);
    ww, 5, [zeros(1,14) 1];
    DD, ones(15,1), UU ];

err(2) = norm( M - MM );

%%
% application to an appropriate chebmatrix
v = [ exp(x); pi; cos(x) ];
Av = A*v;
err(3) = norm( (v{1} + pi*x -3*v{3}) - Av{1} );
err(4) = norm( sum(v{1})+5*v{2}+feval(v{3},dom(end)) - Av{2} );
err(5) = norm( diff(v{1})+pi+(u).*v{3} - Av{3} );

pass = err < 1e-14;

end


text:
function [vnz,rnz,parent,c,leftmost,p,q] = cs_sqr (A)                       %#ok
%CS_SQR symbolic sparse QR factorization.
%   [vnz,rnz,parent,c,leftmost,p] = cs_sqr(A): symbolic QR of A(p,:).
%   [vnz,rnz,parent,c,leftmost,p,q] = cs_sqr(A) computes the symbolic QR
%   factorization of A(p,q).  The fill-reducing ordering q is found via
%   q = cs_amd(A,3).
%
%   vnz is the number of entries in the matrix of Householder vectors, V.
%   rnz is the number of entries in R.  parent is elimination tree.
%   c(i) is the number of entries in R(i,:).  leftmost(i) = min(find(A(i,q))).
%   p is the row permutation used to ensure R has a symbolically zero-free
%   diagonal (it can be larger than m if A is structurally rank deficient).
%   q is the fill-reducing ordering, if requested.
%
%   Example:
%       Prob = UFget ('HB/ibm32') ; A = Prob.A ;
%       [vnz, rnz, parent, c, leftmost, p, q] = cs_sqr (A) ;
%       cspy (A (p,q)) ;
%
%   See also CS_AMD, CS_QR.
%
%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

error ('cs_sqr mexFunction not found') ;


text:
function [vertices, faces] = removeInvalidBorderFaces(varargin)
%REMOVEINVALIDBORDERFACES Remove faces whose edges are connected to 3, 3, and 1 faces.
%
%   [V2, F2] = removeInvalidBorderFaces(V, F)
%
%   Example
%   removeInvalidBorderFaces
%
%   See also
%     isManifoldMesh, collapseEdgesWithManyFaces

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2019-01-31, using Matlab 9.5.0.944444 (R2018b)
% Copyright 2019-2022 INRA - Cepia Software Platform

vertices = varargin{1};
faces = varargin{2};

% compute edge to vertex array
if nargin == 3
    edges = faces;
    faces = varargin{3};
else
    % compute edge to vertex array
    edges = meshEdges(faces);
end

% compute face to edge indices array
% as a nFaces-by-3 array (each face connected to exactly three edges)
faceEdgeInds = meshFaceEdges(vertices, edges, faces);

% compute number of faces incident each edge
edgeFaces = trimeshEdgeFaces(faces);
edgeFaceDegrees = sum(edgeFaces > 0, 2);

% for each face, concatenate the face degree of each edge
faceEdgeDegrees = zeros(size(faces, 1), 3);
for iFace = 1:size(faces, 1)
    edgeInds = faceEdgeInds{iFace};
    faceEdgeDegrees(iFace, :) = edgeFaceDegrees(edgeInds);
end

% remove faces containing edges connected to 1 face and edges connected to
% 3 faces
inds = sum(faceEdgeDegrees == 1, 2) > 0 & sum(faceEdgeDegrees == 3, 2);
% inds = sum(ismember(faceEdgeDegrees, [1 3 4]), 2) == 3;
faces(inds, :) = [];


text:
function str = cstring(x)
%
%CSTRING Print to a string like c
%  CSTRING(x) takes an array x, of type char, and prints it to
%  a string the way C would.
%  the first null, or non-ascii character is the terminator

i=1;
while (x(i)>0 & x(i)<128 & i<length(x))
  i=i+1;
end

if i==1
  str = char(0);
else
  str = char(x(1:i-1));
end

str = sprintf(str,'%s');

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function x = calcAxisDistribution(S3F,varargin)
% axis distribution function
%
% Syntax
%
%   value = calcAxisDistribution(odf, a)
%   adf = calcAxisDistribution(odf)
%
% Input
%  S3F - orientation or misorientation distribution function, @SO3Fun
%  a   - rotational axis, @vector3d
%
% Output
%  afd   - axis distribution function, @S2Fun
%  value - value of axis distribution function for rotational axis a
%
% See also
% symmetry/calcAxisDistribution

[oR,dcs,nSym] = fundamentalRegion(S3F.CS,S3F.SS,varargin{:});

if nargin == 1 || ~isa(varargin{1},'vector3d')
  adf = @(h) calcAxisDistribution(S3F,h,varargin{:});
  x = S2FunHarmonicSym.quadrature(adf,dcs,'bandwidth',64,varargin{:});
  return
end

h = varargin{1};

maxOmega = oR.maxAngle(project2FundamentalRegion(h,dcs));
res = get_option(varargin,'resolution',2.5*degree);
nOmega = round(max(maxOmega(:))/res);

% define a grid for quadrature
omega = linspace(0,1,nOmega);
omega = maxOmega(:) * omega(:).';
h = repmat(h(:),1,nOmega);
S3G = orientation.byAxisAngle(h,omega,S3F.CS,S3F.SS);

% quadrature weights
weights = sin(omega./2).^2 ./ nOmega;

% eval ODF
f = eval(S3F,S3G,varargin{:});

% sum along axes
x = 2*nSym / pi * sum(f .* weights,2) .* maxOmega(:);

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function P = Gamma2aCCDF (v, a)
% This function calculates the complementary cumulative probability
% distribution function (CDF) for a unit-variance (two-sided) generalized
% gamma probability density function (PDF) with parameter a,
%
%   p(x,a) = b exp(-b|x|) (b|x|)^(a-1) / (2 G(a)),   a >= 0.
%
% where G(a) is the complete gamma function.
%
% The complementary cumulative distribution function is
%                       inf
%   G(v,a) = b/(2 G(a)) Int exp(-b|x|) (b|x|)^(a-1) dx.
%                        v
% To develop this formula, we start from the one-sided gamma PDF
%
%  G(x,a) = 1/G(a) exp(-x) x^(a-1),  for x >= 0,
%         = 0,                       for x < 0.
% We now make the PDF two-sided
%   Gamma2PDF(x,a) = 1/2 GammaPDF(|x|,a).
% This is a zero-mean (symmetric) distribution. Its variance can be
% calculated integrating the product of t^2 and Gamma2PDF(x,a). But
% multiplying the PDF by t^2 is the same as replacing the exponent a-1
% by a+1. Then the variance is
%   Var(Gamma2PDF(x,a) = G(a+2)/G(a) = a(a+1),
% since G(a+1) = a G(a).
%
% Let the standard deviation be denoted as b = sqrt(a(a+1)). Then scaling
% x by b and scaling the result, gives the unit-variance generalized gamma
% density shown above. The complementary cumulative distribution function
% is then (for v >= 0),
%                      inf
%  G(v,a) = b/(2 G(a)) Int exp(-b|x|) (b|x|)^(a-1) dx
%                       v
%                      inf
%         = 1/(2 G(a)) Int exp(-u) u^(a-1) du
%                      bv
%         = 1/2 Gamma1aCCDF(bv,a).

% Moments
% Let p(x,a1) be the generalized gamma distribution with parameter a1.
%   F0(x,a1) = Int(x,Inf) p(x,a1) dx
% The moment integrals are
%   Fn(x,a1) = Int(x,Inf) x^n p(x,a1) dx
%            = G(a2)/(G(a1) b1^n) F0(b1 x/b2, a2)       n even
%            = G(a2)/(G(a1) b1^n) F0(b1 abs(x)/b2, a2), n odd
% where a2=a1+n, b1=sqrt(b1(b1+1)), and b2=sqrt(a2(a2+1)). Note that the
% generalized gamma PDF uses abs(x). So for odd moments (n odd), we have to
% compensate by changing the sign of the argument.

% Gaussian
% The Gaussian complementary cumulative distribution function can be
% computed from the generalized gamma distribution. Consider the moments
% of the Gaussian distribution
%   Fn(u) = 1/sqrt(2pi) Int u^n e^(-u^2/2) du
% Consider u positive. Substitute x = u^2/(2b) with u = sqrt(2bx). Then
% du = b/sqrt(2) 1/sqrt(bx), giving
%   Fn(x) = 1/sqrt(2pi) Int (2bx)^(n/2) e^(-bx) b/sqrt(2) 1/sqrt(bx) dx
%         = 1/sqrt(2pi) 2^(n/2)/sqrt(2) b Int (bx)^((n-1)/2) e^(-bx) dx
% Let a = (n+1)/2,
%   Fn(x) = 2^(n/2) G(a)/sqrt(pi) b/(2G(a)) Int exp(-bx) (bx)^(a-1) dx
%         = 2^(n/2) G(a)/sqrt(pi) G(x,a)  with x>0 and b=sqrt(a(a+1)).
%
% For the integral of the Gaussian, n=0, a=1/2, b=sqrt(3)/2), and
% G(a)=sqrt(pi))
%   F0(x) = G(x,1/2)
% or
%   F0(u) = G(sgn(x) x^2/(2b),1/2).
% The sgn() function compensates for the loss of sign information in
% squaring x.
%
% The Gaussian first moment is given by n=1, a=1, b=sqrt(2), and
% G(a)=sqrt(2),
%   F1(u) = sqrt(2/pi) G(x^2/(2b), 1)
% The sgn() function is not necessary since the integrand is odd.
%
% The Gaussian second moment is given by n=2, a=3/2, b=sqrt(3)/2, and
% G(a)=sqrt(pi)/2,
%   F2(u) = G(sgn(x) x^2/(2b), 3/2)

b = sqrt(a*(a+1));
if (v >= 0)
  P = (1/2) * Gamma1aCCDF(b*v, a);
else
  P = 1 - (1/2) * Gamma1aCCDF(-b*v, a);
end

return


text:
function y = choosvd(n, d)

if n <= 100
    if d / n <= 0.02
        y = 1;
    else
        y = 0;
    end
elseif n <= 200
    if d / n <= 0.06
        y = 1;
    else
        y = 0;
    end
elseif n <= 300
    if d / n <= 0.26
        y = 1;
    else
        y = 0;
    end
elseif n <= 400
    if d / n <= 0.28
        y = 1;
    else
        y = 0;
    end
elseif n <= 500
    if d / n <= 0.34
        y = 1;
    else
        y = 0;
    end
else
    if d / n <= 0.38
        y = 1;
    else
        y = 0;
    end
end

author: andrewssobral

repo: lrslibrary


text:
function [mimgR,mimgG] = regRedGreenChannel(ops)

% numPlanes = length(ops.planesToProcess);

% build file list with red channel

if (isfield(ops, 'SubDirsRed') && ~isempty(ops.SubDirsRed))
    subDirsRed = ops.SubDirsRed;
else
    if (isfield(ops, 'expred') && ~isempty(ops.expred))
        for i = 1:length(ops.expred)
            subDirsRed{i} = sprintf('%d', ops.expred(i));
        end
    else
        warning('could not find red channel info, returning...')
        return;
    end
end
%
% build file list
for i = 1:length(ops.expred)
    fsRED = dir(fullfile(ops.RootDir, subDirsRed{i}, '*.tif'));
    for k = 1:length(fsRED)
        fsRED(k).name = fullfile(ops.RootDir, subDirsRed{i}, fsRED(k).name);
    end
end

root = ops.ResultsSavePath;
fregops =  sprintf('regops_%s_%s.mat', ops.mouse_name, ops.date);
if exist(fullfile(root, fregops), 'file')
    load(fullfile(root, fregops))
else
    ops1 = cell(ops.nplanes, 1);
    for j = 1:ops.nplanes

        fname = sprintf('regops_%s_%s_plane%d.mat', ops.mouse_name, ops.date, j);
        dat = load(fullfile(root, fname));
        ops1{j} = dat.ops;
        ops1{j}.useGPU = ops.useGPU;
    end
end

%
ntf0 = 0;
numPlanes = ops.nplanes;
%iplane0 = 1:1:ops.nplanes;

totFrames=0;
for k = 1:length(fsRED)
    %iplane0 = mod(iplane0-1, numPlanes) + 1;
    startPlane=((mod(totFrames+1,ops.nplanes*2)-1)/2)+1;

    nFr = nFramesTiff(fsRED(k).name);
    totFrames=totFrames+nFr;
    data = loadFramesBuff(fsRED(k).name, 1, nFr, 1, ops.temp_tiff);

    if ~exist('mimgR', 'var')
        [Ly, Lx, ~] = size(data);
        mimgR = zeros(Ly, Lx, ops.nplanes);
    end
    if ~exist('mimgG', 'var')
        [Ly, Lx, ~] = size(data);
        mimgG = zeros(Ly, Lx, ops.nplanes);
    end
    %

    for iPlane=1:ops.nplanes

        idx0=mod((ops.nplanes-startPlane+iPlane)*2,ops.nplanes*2);
        planesG=(idx0+1):(2*ops.nplanes):nFr;
        planesR=(idx0+2):(2*ops.nplanes):nFr;
        dataG0=data(:,:,planesG);
        dataR0=data(:,:,planesR);

        BiDiPhase = ops1{iPlane}.BiDiPhase;
        if abs(BiDiPhase) > 0
            yrange = 2:2:Ly;
            if BiDiPhase>0
                dataG0(yrange,(1+BiDiPhase):Lx,:,:) = dataG0(yrange, 1:(Lx-BiDiPhase),:,:);
                dataR0(yrange,(1+BiDiPhase):Lx,:,:) = dataR0(yrange, 1:(Lx-BiDiPhase),:,:);
            else
                dataG0(yrange,1:Lx+BiDiPhase,:,:)   = dataG0(yrange, 1-BiDiPhase:Lx,:,:);
                dataR0(yrange,1:Lx+BiDiPhase,:,:)   = dataR0(yrange, 1-BiDiPhase:Lx,:,:);
            end
        end

        [ds, ~]  = regoffKriging(dataG0, ops1{iPlane}, 0);
        %[ds, ~]  = registration_offsets(dataG0, ops1{iPlane}, 0);

        if k==1
            ds(1,:) = 0;
        end
        dataR       = ...
            register_movie(dataR0, ops1{iPlane}, ds);
        dataG     = ...
            register_movie(dataG0, ops1{iPlane}, ds);

        mimgR(:,:,iPlane) = mimgR(:,:,iPlane) + mean(dataR, 3);
        mimgG(:,:,iPlane) = mimgG(:,:,iPlane) + mean(dataG, 3);

    end

    ntf0 = ntf0 + 1;

    %iplane0 = iplane0 - nFr/ops.nchannels_red;
    fprintf('processing tiff %d/%d\n',k,length(fsRED))
end

mimgR = mimgR/ntf0;
mimgG = mimgG/ntf0;

author: cortex-lab


text:
function [trackers, indices] = order_trackers(trackers, order)
% order_trackers Change the order of trackers
%
% Change the order of trackers in a cell array of tracker structures.
%
% The functions upports the following ordering options:
%
% - IdentifierAscending
% - IdentifierDescending
%
% Input:
% - trackers: Cell array of tracker structures.
% - order: A string that defines the order of trackers. See the description above for available options.
%
% Output:
% - trackers: A modified cell array of tracker structures.
% - indices: The mapping from the old to the new array.

if nargin < 2
    order = 'identifierascending';
end;

order = lower(order);

switch order
    case 'identifierascending'
        identifiers = cellfun(@(t) t.identifier, trackers, 'UniformOutput', false);
        [~, indices] = sort(identifiers);
    case 'identifierdescending'
        identifiers = cellfun(@(t) t.identifier, trackers, 'UniformOutput', false);
        [~, indices] = sort(identifiers);
        indices = fliplr(indices);
    otherwise
        return;
end;

trackers = trackers(indices);

author: votchallenge

repo: toolkit-legacy


text:
% capillary pressure curve for drainage
% Written by Ali A. Eftekhari
% TBD, probably not required in the new IMPES formulation
function res=dpc_drain(sw, pce, swc, labda, pc_max)
sw0=swc+(1-labda*log(pc_max/pce)+sqrt((-1+labda*log(pc_max/pce))^2+...
      4*swc/(1-swc)))/2*(1-swc);
res = (sw>sw0).*(-1.0/((1-swc)*labda)*pce*((sw-swc)/(1-swc)).^(-1.0/labda-1))+...
    (0.0<=sw).*(sw<=sw0).*(-1.0/((1-swc)*labda)*pce*((sw0-swc)/(1-swc)).^(-1.0/labda-1))+...
    (sw<0)*0.0;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function value = p19_f ( dim_num, point_num, x )

%*****************************************************************************80
%
%% P19_F evaluates the integrand for problem 19.
%
%  Dimension:
%
%    DIM_NUM is arbitrary.
%
%  Region:
%
%    0 <= X(1:DIM_NUM) <= 1
%
%  Integral Parameters:
%
%    Z defaults to (1/3,1/3,...,1/3).
%    You can reset Z by calling P19_R8VEC.
%
%  Integrand:
%
%    f(x) = product ( sqrt ( abs ( x(1:dim_num) - z(1:dim_num) ) ) )
%
%  Exact Integral:
%
%    With Z as given,
%
%      (2/3)**DIM_NUM * ( (2/3)**(3/2) + (1/3)**(3/2) )**DIM_NUM
%
%    or approximately 0.49**DIM_NUM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arnold Krommer, Christoph Ueberhuber,
%    Numerical Integration on Advanced Systems,
%    Springer, 1994,
%    ISBN: 3540584102,
%    LC: QA299.3.K76.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the evaluation points.
%
%    Output, real VALUE(POINT_NUM), the integrand values.
%
  z = [];
  z = p19_r8vec ( 'G', 'Z', dim_num, z );

  value(1:point_num) = 0.0;

  for point = 1 : point_num
    value(point) = prod ( sqrt ( abs ( x(1:dim_num,point) - z(1:dim_num)' ) ) );
  end

  p19_i4 ( 'I', '#', point_num );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function V = tapas_uniqc_spm_create_vol_with_tr(V)
% Create a NIfTI image volume
% FORMAT V = spm_create_vol(V)
% V        - image volume information (see spm_vol.m)
%__________________________________________________________________________
% Copyright (C) 2005-2014 Wellcome Trust Centre for Neuroimaging

% John Ashburner

for i=1:numel(V)
    v = create_vol(V(i));

    f = fieldnames(v);
    for j=1:numel(f)
        V(i).(f{j}) = v.(f{j});
    end
end

%==========================================================================
%-function V = create_vol(V)
%==========================================================================
function V = create_vol(V)

if ~isstruct(V),        error('tapas:uniqc:NotStructure', ...
        'Not a structure.'); end

if ~isfield(V,'fname'), error('tapas:uniqc:MissingField', ...
        'No "fname" field'); end

if ~isfield(V,'dim'),   error('tapas:uniqc:MissingField', ...
        'No "dim" field'); end
if ~all(size(V.dim)==[1 3])
    error('tapas:uniqc:WrongFieldSizeDim', ...
        ['"dim" field is the wrong size (' num2str(size(V.dim)) ').']);
end

if ~isfield(V,'n')
    V.n = [1 1];
else
    V.n = [V.n(:)' 1 1];
    V.n =  V.n(1:2);
end
if V.n(1)>1 && V.n(2)>1
    error('tapas:uniqc:TooManyDimensions', ...
        'Can only do up to 4D data (%s).',V.fname);
end

if ~isfield(V,'dt')
    V.dt = [spm_type('float64') spm_platform('bigend')];
end
dt{1} = spm_type(V.dt(1));
if strcmp(dt{1},'unknown')
    error('tapas:uniqc:UnrecognizedDatatype', ...
        ['"' dt{1} '" is an unrecognised datatype (' num2str(V.dt(1)) ').']);
end
if V.dt(2), dt{2} = 'BE'; else dt{2} = 'LE'; end

if ~isfield(V,'pinfo'), V.pinfo      = [Inf Inf 0]'; end
if size(V.pinfo,1)==2,  V.pinfo(3,:) = 0;            end

V.fname  = deblank(V.fname);
ext      = spm_file(V.fname,'ext');

switch ext
    case {'img'}
        minoff = 0;
    case {'nii'}
        minoff = 352; % or 544 for NIfTI-2
    otherwise
        error('tapas:uniqc:InvalidFileExtension', ...
            ['".' ext '" is not a recognised extension.']);
end
bits   = spm_type(V.dt(1),'bits');
minoff = minoff + ceil(prod(V.dim(1:2))*bits/8)*V.dim(3)*(V.n(1)-1+V.n(2)-1);
V.pinfo(3,1) = max(V.pinfo(3,:),minoff);

if ~isfield(V,'descrip'), V.descrip = '';     end
if ~isfield(V,'private'), V.private = struct; end

dim    = [V.dim(1:3) V.n];
dat    = file_array(V.fname,dim,[dt{1} '-' dt{2}],0,V.pinfo(1),V.pinfo(2));
N      = nifti;
N.dat  = dat;
N.mat  = V.mat;
N.mat0 = V.mat;
N.mat_intent  = 'Aligned';
N.mat0_intent = 'Aligned';
N.descrip = V.descrip;
%% changed by Saskia Bollmann, 25/09/2015
try N.timing = V.private.timing; end %#ok<TRYNC>
%% end

try
    N0  = nifti(V.fname);

    % Just overwrite if both are single volume files.
    tmp = [N0.dat.dim ones(1,5)];
    if prod(tmp(4:end))==1 && prod(dim(4:end))==1
        N0 = [];
    end
catch
    N0  = [];
end

if ~isempty(N0)

    % If the dimensions differ, then there is the potential for things to go badly wrong.
    tmp = [N0.dat.dim ones(1,5)];
    if any(tmp(1:3) ~= dim(1:3))
        warning('tapas:uniqc', ...
            ['Incompatible x,y,z dimensions in file "' V.fname '" [' num2str(tmp(1:3)) ']~=[' num2str(dim(1:3)) '].']);
    end
    if dim(5) > tmp(5) && tmp(4) > 1
        warning('tapas:uniqc', ...
            ['Incompatible 4th and 5th dimensions in file "' V.fname '" (' num2str([tmp(4:5) dim(4:5)]) ').']);
    end
    N.dat.dim = [dim(1:3) max(dim(4:5),tmp(4:5))];

    if ~strcmp(dat.dtype,N0.dat.dtype)
        warning('tapas:uniqc', ...
            ['Incompatible datatype in file "' V.fname '" ' N0.dat.dtype ' ~= ' dat.dtype '.']);
    end
    if single(N.dat.scl_slope) ~= single(N0.dat.scl_slope) && (size(N0.dat,4)>1 || V.n(1)>1)
        warning('tapas:uniqc', ...
            ['Incompatible scalefactor in "' V.fname '" ' num2str(N0.dat.scl_slope) '~=' num2str(N.dat.scl_slope) '.']);
    end
    if single(N.dat.scl_inter) ~= single(N0.dat.scl_inter)
        warning('tapas:uniqc', ...
            ['Incompatible intercept in "' V.fname '" ' num2str(N0.dat.scl_inter) '~=' num2str(N.dat.scl_inter) '.']);
    end
    if single(N.dat.offset) ~= single(N0.dat.offset)
        warning('tapas:uniqc', ...
            ['Incompatible intercept in "' V.fname '" ' num2str(N0.dat.offset) '~=' num2str(N.dat.offset) '.']);
    end

    if V.n(1)==1

        % Ensure volumes 2..N have the original matrix
        nt = size(N.dat,4);
        if nt>1 && sum(sum((N0.mat-V.mat).^2))>1e-8
            M0 = N0.mat;
            if ~isfield(N0.extras,'mat')
                N0.extras.mat = zeros([4 4 nt]);
            else
                if size(N0.extras.mat,4)<nt
                    N0.extras.mat(:,:,nt) = zeros(4);
                end
            end
            for i=2:nt
                if sum(sum(N0.extras.mat(:,:,i).^2))==0
                    N0.extras.mat(:,:,i) = M0;
                end
            end
            N.extras.mat = N0.extras.mat;
        end

        N0.mat = V.mat;
        if strcmp(N0.mat0_intent,'Aligned'), N.mat0 = V.mat; end
        if ~isempty(N.extras) && isstruct(N.extras) && isfield(N.extras,'mat') &&...
                size(N.extras.mat,3)>=1
            N.extras.mat(:,:,V.n(1)) = V.mat;
        end
    else
        if sum(sum((N0.mat-V.mat).^2))>1e-8
            N.extras.mat(:,:,V.n(1)) = V.mat;
        end
    end

    if ~isempty(N0.extras) && isstruct(N0.extras) && isfield(N0.extras,'mat')
        N0.extras.mat(:,:,V.n(1)) = N.mat;
        N.extras                  = N0.extras;
    end
    if sum((V.mat(:)-N0.mat(:)).^2) > 1e-4
        N.extras.mat(:,:,V.n(1)) = V.mat;
    end
end

if isfield(N.extras,'mat')
    M0 = N.mat;
    for i=1:size(N.extras.mat,3)
        if sum((M0-N.extras.mat(:,:,i)).^2) < 1e-8
            N.extras.mat(:,:,i) = 0;
        end
    end
    if sum(N.extras.mat(:).^2) < 1e-8*size(N.extras.mat,3)
        N.extras = [];
        if ~isempty(N0) && ~isempty(N0.extras) && isstruct(N0.extras) && isfield(N0.extras,'mat')
            warning('tapas:uniqc', ...
                'Forcing deletion of MAT-file.');
            spm_unlink(spm_file(N.dat.fname,'ext','mat'));
        end
    end
end

create(N);
V.private = N;

author: translationalneuromodeling

repo: tapas


text:
function varargout = spm_graph_ui(action, varargin)
% Graphical display of adjusted data
% FORMAT [Y,y,beta,Bcov] = spm_graph_ui(xSPM,SPM,hReg)
%
% xSPM   - structure containing SPM, distributional & filtering details
%          about the excursion set
% SPM    - structure containing generic details about the analysis
% hReg   - handle of MIP register or [x y z] coordinates
%
% Y      - fitted   data for the selected voxel
% y      - adjusted data for the selected voxel
% beta   - parameter estimates (ML or MAP)
% Bcov   - Covariance of parameter estimates (ML or conditional)
%
% See spm_getSPM for details.
%__________________________________________________________________________
%
% spm_graph is a Callback function that uses the structures above to:
% (1) send adjusted (y) and fitted data (Y), for the selected voxel, to the
% workspace and (ii) provide graphics for:
%
% a) Contrasts of parameter estimates (e.g. activations) and their
% standard error.
%
% b) Fitted and adjusted responses that can be plotted against time, scan,
% or an indicator variable in the design matrix.
%
% c) (fMRI only).  Evoked responses using the basis functions to give
% impulse responses that would have been seen in the absence of other
% effects. The PSTH (peristimulus-time histogram) option provides a finite
% impulse response (FIR) estimate of the trial-specific evoked response as
% a function of peristimulus time.  This is estimated by refitting a
% convolution model to the selected voxel using an FIR basis set.  This is
% simply a set of small boxes covering successive time bins after trial
% onset.  The width of each bin is usually the TR.  This option provides a
% more time-resolved quantitative characterisation of the evoked
% hemodynamic response.  However, it should not be over-interpreted because
% inference is usually made using a simpler and more efficient basis set
% (e.g., canonical hrf, or canonical plus time derivative).
%
% Getting adjusted data:
% Ensuring the data are adjusted properly can be important (e.g. in
% constructing explanatory variables such as in a psychophysiological
% interaction). To remove or correct for specific effects, specify an
% appropriate F contrast and simply plot the fitted (and adjusted)
% responses after selecting that F contrast. The vectors Y (fitted) and y
% (adjusted) in the workspace will now be corrected for the effects in the
% reduced design matrix (X0) specified in the contrast manager with the
% column indices (iX0) of the confounds in this adjustment.
%
% Plotting data:
% All data and graphics use filtered/whitened data and residuals. In PET
% studies the parameter estimates and the fitted data are often the same
% because the explanatory variables are simply indicator variables taking
% the value of one.  Only contrasts previously defined can be plotted. This
% ensures that the parameters plotted are meaningful even when there is
% collinearity among the design matrix subpartitions.
%
% Selecting contrasts used for PPMs will automatically give plots
% based on conditonal estimates.
%
% The structure     contrast.contrast      = cbeta;
%                   contrast.standarderror = SE;
%                   contrast.interval      = 2*CI;
%
% is assigned in base workspace for plots of contrasts and their error.
%__________________________________________________________________________
% Copyright (C) 1996-2013 Wellcome Trust Centre for Neuroimaging

% Karl Friston, Guillaume Flandin
% $Id: spm_graph_ui.m 6985 2017-01-11 11:51:49Z guillaume $

if nargin && ~ischar(action)
    varargin = {action varargin{:}};
    action   = 'Graph';
end

%==========================================================================
switch lower(action), case 'graph'                                  %-Graph
%==========================================================================
% FORMAT [Y,y,beta,Bcov] = spm_graph_ui('Graph',xSPM,SPM,hReg)

    if nargin < 3, error('Not enough input arguments.'); end
    [xSPM,SPM,hReg] = deal(varargin{:});

    [xG, XYZ] = spm_graph_ui('UI',xSPM,SPM,hReg);
    if isempty(xG), varargout = {[],[],[],[]}; return; end

    G = spm_graph_ui('Compute',SPM,XYZ,xG);

    XYZmm = SPM.xVol.M(1:3,:) * [XYZ(:);1];
    spm_graph_ui('Plot',SPM,XYZmm,xG,G,hReg);

    varargout = { G.Y, G.y, G.beta, G.Bcov };

%==========================================================================
case 'ui'                                                              %-UI
%==========================================================================
% FORMAT [xG, XYZ] = spm_graph_ui('UI',xSPM,SPM,hReg)

    if nargin < 3, error('Not enough input arguments.'); end
    [xSPM,SPM,hReg] = deal(varargin{:});

    %-Delete previous axis and their pagination controls (if any)
    %----------------------------------------------------------------------
    Fgraph = spm_figure('GetWin','Graphics');
    spm_results_ui('Clear',Fgraph,2);

    %-Find nearest voxel [Euclidean distance] in point list & update GUI
    %----------------------------------------------------------------------
    if isempty(xSPM.XYZmm)
        spm('alert!','No suprathreshold voxels!',mfilename,0);
        varargout = {[],[]};
        return
    end

    if numel(hReg) == 1
        XYZmm = spm_XYZreg('GetCoords',hReg);
    else
        XYZmm = hReg;
    end
    [XYZmm,i] = spm_XYZreg('NearestXYZ',XYZmm,xSPM.XYZmm);
    if numel(hReg) == 1, spm_XYZreg('SetCoords',XYZmm,hReg); end
    XYZ       = xSPM.XYZ(:,i);

    %-Find out what to plot
    %======================================================================
    Cplot   = {'Contrast estimates and 90% C.I.',...
               'Fitted responses',...
               'Event-related responses',...
               'Parametric responses',...
               'Volterra Kernels'};

    % ensure options are appropriate
    %----------------------------------------------------------------------
    if ~isfield(SPM,'Sess'), Cplot = Cplot(1:2); end
    i       = spm_input('Plot',-1,'m',Cplot);
    xG.def  = Cplot{i};

    switch xG.def

        % select contrast if
        %------------------------------------------------------------------
        case {'Contrast estimates and 90% C.I.'}

            % determine which contrast
            %--------------------------------------------------------------
            xG.spec.Ic = spm_input('Which contrast?','!+1','m',...
                {SPM.xCon.name,'Specify a temporary contrast...'});
            if xG.spec.Ic > numel(SPM.xCon)
                xG.spec.Ic = spm_input('Contrast weights','!+1','x','',Inf,SPM.xX.X);
            end

        % select contrast if
        %------------------------------------------------------------------
        case {'Fitted responses'}

            % determine which contrast
            %--------------------------------------------------------------
            xG.spec.Ic = spm_input('Which contrast?','!+1','m',{SPM.xCon.name});

        % select session and trial if
        %------------------------------------------------------------------
        case {'Event-related responses','Parametric responses','Volterra Kernels'}

            % get session
            %--------------------------------------------------------------
            s     = length(SPM.Sess);
            if  s > 1
                s = spm_input('which session','+1','n1',1,s);
            end
            xG.spec.Sess = s;

            % effect names
            %--------------------------------------------------------------
            switch xG.def
                case 'Volterra Kernels'
                    u = length(SPM.Sess(s).Fc);
                otherwise
                    u = length(SPM.Sess(s).U);
            end
            Uname = {};
            for i = 1:u
                Uname{i} = SPM.Sess(s).Fc(i).name;
            end
            if isempty(Uname), error('No conditions found.'); end

            % get effect
            %--------------------------------------------------------------
            str   = 'which effect';
            u     = spm_input(str,'+1','m',Uname);
            xG.spec.u = u;

    end

    switch xG.def
        case 'Fitted responses'

            % predicted or adjusted response
            %--------------------------------------------------------------
            str   = 'predicted or adjusted response?';
            xG.spec.predicted = spm_input(str,'!+1','b',{'predicted','adjusted'},[1 0]);

            % get ordinates
            %--------------------------------------------------------------
            Xplot = {'an explanatory variable',...
                     'scan or time',...
                     'a user specified ordinate'};
            Cx    = spm_input('plot against','!+1','m',Xplot);
            if     Cx == 1
                str  = 'Which explanatory variable?';
                xG.spec.x.i = spm_input(str,'!+1','m',SPM.xX.name);
            elseif Cx == 2
                xG.spec.x.scan = true;
            elseif Cx == 3
                xG.spec.x.x = spm_input('enter ordinate','!+1','e','',numel(SPM.xY.VY));
            end

        %-Modeling evoked responses based on Sess
        %==================================================================
        case 'Event-related responses'

            % get plot type
            %--------------------------------------------------------------
            Rplot = {'fitted response and PSTH',...
                     'fitted response and 90% C.I.',...
                     'fitted response and adjusted data'};

            Rplot = Rplot{spm_input('plot in terms of','+1','m',Rplot)};
            xG.spec.Rplot = Rplot;

        %-Modeling evoked responses based on Sess
        %==================================================================
        case 'Parametric responses'

            % return gracefully if no parameters
            %--------------------------------------------------------------
            if ~SPM.Sess(s).U(u).P(1).h
                warning('No parametric modulation in condition "%s".',...
                    SPM.Sess(s).Fc(u).name);
                varargout = {[],[]};
                return
            end

            str   = 'which parameter';
            p     = spm_input(str,'+1','m',{SPM.Sess(s).U(u).P.name});
            xG.spec.p = p;
    end

    varargout = { xG, XYZ };

%==========================================================================
case 'compute'                                                    %-Compute
%==========================================================================
% FORMAT [Y,y,beta,Bcov,G] = spm_graph_ui('Compute',SPM,XYZ,xG)

spm('Pointer','Watch');
[Y,y,beta,Bcov,G] = spm_graph(varargin{:});
spm('Pointer','Arrow');

G.Y = Y; G.y = y; G.beta = beta; G.Bcov = Bcov;
varargout = { G };

%==========================================================================
case 'plot'                                                          %-Plot
%==========================================================================
% FORMAT spm_graph_ui('Plot',SPM,XYZmm,xG,G,hReg)

    [SPM,XYZmm,xG,G,hReg] = deal(varargin{:});

    %-Get Graphics figure handle
    %----------------------------------------------------------------------
    Fgraph = spm_figure('GetWin','Graphics');

    %-Scaled font sizes
    %----------------------------------------------------------------------
    FS     = spm('FontSizes');

    %-Voxel coordinates
    %----------------------------------------------------------------------
    XYZstr = sprintf(' at [%g, %g, %g]',XYZmm);

    %-Colour specification
    %----------------------------------------------------------------------
    Col    = [0 0 0; .8 .8 .8; 1 .5 .5];

    switch xG.def

        %-Plot parameter estimates
        %==================================================================
        case 'Contrast estimates and 90% C.I.'

            % bar chart
            %--------------------------------------------------------------
            figure(Fgraph)
            ax = subplot(2,1,2,'Parent',Fgraph);
            cla(ax)
            hold(ax,'on')

            % estimates
            %--------------------------------------------------------------
            cbeta = G.contrast;
            h     = bar(ax,cbeta);
            set(h,'FaceColor',Col(2,:))

            % standard error
            %--------------------------------------------------------------
            CI    = G.interval / 2;
            for j = 1:length(cbeta)
                line([j j],([CI(j) -CI(j)] + cbeta(j)),...
                    'LineWidth',6,'Color',Col(3,:),'Parent',ax)
            end

            TTLstr = {xG.def};
            if numel(xG.spec.Ic) == 1
                TTLstr = [TTLstr {SPM.xCon(xG.spec.Ic).name}];
            else
                TTLstr = [TTLstr {'user-defined contrast'}];
            end
            if isfield(SPM,'VCbeta')
                TTLstr = [TTLstr {'(conditional estimates)'}];
            end
            title(ax,TTLstr,'FontSize',FS(12))
            xlabel(ax,'contrast','FontSize',FS(12))
            ylabel(ax,['contrast estimate',XYZstr],'FontSize',FS(12))
            set(ax,'XLim',[0.4 (length(cbeta) + 0.6)])
            hold(ax,'off')

            % for compatibility with earlier versions of SPM
            %--------------------------------------------------------------
            assignin('base','contrast',G)

        %-All fitted effects or selected effects
        %==================================================================
        case 'Fitted responses'

            % plot
            %--------------------------------------------------------------
            figure(Fgraph)
            ax = subplot(2,1,2,'Parent',Fgraph);
            cla(ax)
            hold(ax,'on')

            x = G.x;
            y = G.y;
            Y = G.Y;
            [p,q] = sort(x);
            if all(diff(x(q)))
                plot(x(q),Y(q),'LineWidth',4,'Color',Col(2,:));
                plot(x(q),y(q),':','Color',Col(1,:));
                plot(x(q),y(q),'.','MarkerSize',8, 'Color',Col(3,:));
            else
                plot(x(q),Y(q),'.','MarkerSize',16,'Color',Col(1,:));
                plot(x(q),y(q),'.','MarkerSize',8, 'Color',Col(2,:));
                xlim = get(ax,'XLim');
                xlim = [-1 1]*diff(xlim)/4 + xlim;
                set(ax,'XLim',xlim)
            end

            TTLstr = {xG.def SPM.xCon(xG.spec.Ic).name};
            if isfield(SPM,'VCbeta')
                TTLstr = [TTLstr {'(conditional estimates)'}];
            end
            title(ax,TTLstr,'FontSize',FS(12))
            switch char(fieldnames(xG.spec.x))
                case 'i' % an explanatory variable
                    XLAB = SPM.xX.name{xG.spec.x.i};
                case 'scan' % scan or time
                    if isfield(SPM.xY,'RT')
                        XLAB = 'time {seconds}';
                    else
                        XLAB = 'scan number';
                    end
                case 'x' % user specified
                    XLAB = 'ordinate';
            end
            xlabel(ax,XLAB,'FontSize',FS(12))
            ylabel(ax,['response',XYZstr],'FontSize',FS(12))
            legend(ax,'fitted','plus error')
            hold(ax,'off')

        %-Modeling evoked responses based on Sess
        %==================================================================
        case 'Event-related responses'

            % plot
            %--------------------------------------------------------------
            figure(Fgraph)
            ax = subplot(2,1,2,'Parent',Fgraph);
            hold(ax,'on')

            x = G.x;
            y = G.y;
            Y = G.Y;
            s = xG.spec.Sess;
            u = xG.spec.u;
            Rplot = xG.spec.Rplot;
            if strcmp(Rplot,'fitted response and PSTH') && ~isfield(G,'PSTH')
                % Data not available
                Rplot = 'fitted response and 90% C.I.';
            end
            switch Rplot

                case 'fitted response and PSTH'
                    %------------------------------------------------------
                    PST  = G.PST;
                    PSTH = G.PSTH;
                    PCI  = G.PCI;
                    errorbar(ax,PST,PSTH,PCI)
                    plot(ax,PST,PSTH,'LineWidth',4,'Color',Col(2,:))
                    plot(ax,x,Y,'-.','Color',Col(3,:))

                case 'fitted response and 90% C.I.'
                    %------------------------------------------------------
                    CI   = G.CI;
                    plot(ax,x,Y,'Color',Col(2,:),'LineWidth',4)
                    plot(ax,x,Y + CI,'-.',x,Y - CI,'-.','Color',Col(1,:))

                case 'fitted response and adjusted data'
                    %------------------------------------------------------
                    pst  = G.pst;
                    plot(ax,x,Y,'Color',Col(2,:),'LineWidth',4)
                    plot(ax,pst,y,'.','Color',Col(3,:))

            end

            % label
            %--------------------------------------------------------------
            [i,j] = max(Y);
            text(ceil(1.1*x(j)),i,SPM.Sess(s).Fc(u).name,'FontSize',FS(8),'Parent',ax);
            title(ax,Rplot,'FontSize',FS(12))
            xlabel(ax,'peristimulus time {secs}','FontSize',FS(12))
            ylabel(ax,['response',XYZstr],'FontSize',FS(12))
            hold(ax,'off')

        %-Parametric responses
        %==================================================================
        case 'Parametric responses'

            % plot
            %--------------------------------------------------------------
            figure(Fgraph)

            pst = G.pst;
            P   = G.P;
            Y   = G.Y;
            s   = xG.spec.Sess;
            u   = xG.spec.u;
            p   = xG.spec.p;

            ax  = subplot(2,2,3,'Parent',Fgraph);
            surf(ax,pst,P,Y')
            shading(ax,'flat')
            title(ax,SPM.Sess(s).U(u).name{1},'FontSize',FS(12))
            xlabel(ax,'PST {secs}','FontSize',FS(12))
            ylabel(ax,SPM.Sess(s).U(u).P(p).name,'FontSize',FS(12))
            zlabel(ax,['responses',XYZstr],'FontSize',FS(12))
            axis(ax,'square')

            ax = subplot(2,2,4,'Parent',Fgraph);
            [j,i] = max(mean(Y,2));
            plot(ax,P,Y(i,:),'LineWidth',4,'Color',Col(2,:))
            str   = sprintf('response at %0.1fs',i*SPM.xBF.dt);
            title(ax,str,'FontSize',FS(12))
            xlabel(ax,SPM.Sess(s).U(u).P(p).name,'FontSize',FS(12))
            axis(ax,'square')
            grid(ax,'on')

        %-Volterra Kernels
        %==================================================================
        case 'Volterra Kernels'

            pst = G.pst;
            Y   = G.Y;
            s   = xG.spec.Sess;
            u   = xG.spec.u;

            % second order kernel
            %--------------------------------------------------------------
            if u > length(SPM.Sess(s).U)

                % plot
                %----------------------------------------------------------
                figure(Fgraph)

                ax = subplot(2,2,3,'Parent',Fgraph);
                imagesc(pst,pst,Y,'Parent',ax)
                axis(ax,'xy')
                axis(ax,'image')
                title(ax,'2nd order Kernel','FontSize',FS(12))
                xlabel(ax,'peristimulus time {secs}','FontSize',FS(12))
                ylabel(ax,'peristimulus time {secs}','FontSize',FS(12))

                ax = subplot(2,2,4,'Parent',Fgraph);
                plot(ax,pst,Y)
                axis(ax,'square')
                grid(ax,'on')
                title(ax,SPM.Sess(s).Fc(u).name,'FontSize',FS(12))
                xlabel(ax,'peristimulus time {secs}','FontSize',FS(12))

            % first  order kernel
            %--------------------------------------------------------------
            else

                % plot
                %----------------------------------------------------------
                figure(Fgraph)
                ax = subplot(2,1,2,'Parent',Fgraph);
                plot(ax,pst,Y)
                grid(ax,'on')
                axis(ax,'square')
                title(ax,{'1st order Volterra Kernel' SPM.Sess(s).Fc(u).name},...
                    'FontSize',FS(12))
                xlabel(ax,'peristimulus time {secs}','FontSize',FS(12))
                ylabel(ax,['impulse response',XYZstr],'FontSize',FS(12))

            end

    end

    %-Turn hold button off - this will alert the user to press it again
    %----------------------------------------------------------------------
    try, set(findobj('Tag','holdButton'),'Value',0); end

    %-Call Plot UI
    %----------------------------------------------------------------------
    spm_results_ui('PlotUi',ax);

    %-Setup registry
    %----------------------------------------------------------------------
    hPB = findobj('Tag','plotButton');
    if ~isempty(hPB)
        setappdata(hPB,'Reg',struct('xG',xG,'Fgraph',Fgraph));
        spm_XYZreg('Add2Reg',hReg,hPB,'spm_graph_ui');
        set(ax,'DeleteFcn',[get(ax,'DeleteFcn') ';' ...
            'spm_XYZreg(''Del2Reg'',hReg,findobj(''Tag'',''plotButton''));']);
    end

%==========================================================================
case 'setcoords'                                        %-Coordinate change
%==========================================================================
% FORMAT spm_graph_ui('SetCoords',XYZmm,h,hReg)
    if nargin<3, error('Not enough input arguments.'); end
    XYZmm = varargin{1};
    h     = varargin{2};
    hReg  = varargin{3};
    UD    = getappdata(h,'Reg');

    SPM   = evalin('base','SPM');
    XYZ   = SPM.xVol.iM(1:3,:) * [XYZmm(:);1];
    G     = spm_graph_ui('Compute',SPM,XYZ,UD.xG);

    spm_results_ui('Clear',UD.Fgraph,2);

    spm_graph_ui('Plot',SPM,XYZmm,UD.xG,G,hReg);

%==========================================================================
otherwise                                           %-Unknown action string
%==========================================================================
    error('Unknown action string')

end

author: spm


text:
function [ gx, gy ] = g05_xy ( gn )

%*****************************************************************************80
%
%% G05_XY returns the grid points for grid 5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer GN, the grid size.
%
%    Output, real GX(GN,1), GY(GN,1), the grid coordinates.
%
  gx = [ ...
    0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000, ...
    0.125,  0.125,  0.125,  0.125,  0.125,  0.125,  0.125,  0.125,  0.125, ...
    0.250,  0.250,  0.250,  0.250,  0.250,  0.250,  0.250,  0.250,  0.250, ...
    0.375,  0.375,  0.375,  0.375,  0.375,  0.375,  0.375,  0.375,  0.375, ...
    0.500,  0.500,  0.500,  0.500,  0.500,  0.500,  0.500,  0.500,  0.500, ...
    0.625,  0.625,  0.625,  0.625,  0.625,  0.625,  0.625,  0.625,  0.625, ...
    0.750,  0.750,  0.750,  0.750,  0.750,  0.750,  0.750,  0.750,  0.750, ...
    0.875,  0.875,  0.875,  0.875,  0.875,  0.875,  0.875,  0.875,  0.875, ...
    1.000,  1.000,  1.000,  1.000,  1.000,  1.000,  1.000,  1.000,  1.000 ]';

  gy = [ ...
    0.000,  0.125,  0.250,  0.375,  0.500,  0.625,  0.750,  0.875,  1.000, ...
    0.000,  0.125,  0.250,  0.375,  0.500,  0.625,  0.750,  0.875,  1.000, ...
    0.000,  0.125,  0.250,  0.375,  0.500,  0.625,  0.750,  0.875,  1.000, ...
    0.000,  0.125,  0.250,  0.375,  0.500,  0.625,  0.750,  0.875,  1.000, ...
    0.000,  0.125,  0.250,  0.375,  0.500,  0.625,  0.750,  0.875,  1.000, ...
    0.000,  0.125,  0.250,  0.375,  0.500,  0.625,  0.750,  0.875,  1.000, ...
    0.000,  0.125,  0.250,  0.375,  0.500,  0.625,  0.750,  0.875,  1.000, ...
    0.000,  0.125,  0.250,  0.375,  0.500,  0.625,  0.750,  0.875,  1.000, ...
    0.000,  0.125,  0.250,  0.375,  0.500,  0.625,  0.750,  0.875,  1.000 ]';

  return
end


text:
function s = bear_settings_WGP2016_test(excelPath)

s = BEARsettings(2, 'ExcelFile', excelPath);

s.frequency=3;
s.startdate='2014m5';
s.enddate='2018m12';
s.varendo='hicp gdp app 10y stockindex';
s.varexo='';
s.lags=2;
s.const=1;
s.results_path = fullfile(fileparts(mfilename('fullpath')),'results');
s.results_sub='results_test_data_WGP2016_temp';
s.results=1;
s.plot=0;
s.workspace=1;

s.favar.FAVAR=0; % augment VAR model with factors (1=yes, 0=no)

s.prior=41;
s.ar=0.8; % this sets all AR coefficients to the same prior value (if PriorExcel is equal to 0)
s.PriorExcel=0; % set to 1 if you want individual priors, 0 for default
s.priorsexogenous=0; % set to 1 if you want individual priors, 0 for default
s.lambda1=1000; % quasi-flat (diffuse) normal-wishart prior here, in the spirit of Uhlig (2005)
s.lambda2=0.5;
s.lambda3=1;
s.lambda4=1;
s.lambda5=0.001;
s.lambda6=1;
s.lambda7=0.1;
s.lambda8=1;
s.It=5000;
s.Bu=1000;
s.hogs=0;
s.bex=0;
s.scoeff=0;
s.iobs=0;
s.lrp=0;
s.priorf=100;

s.IRF=1;
s.IRFperiods=36;
s.F=0;
s.FEVD=0;
s.HD=0;
s.HDall=0;%if we want to plot the entire decomposition, all contributions (includes deterministic part)HDall
s.CF=0;
s.IRFt=4;

s.strctident.MM=0; % option for Median model (0=no (standard), 1=yes)
s.strctident.CorrelShock=''; % exact labelname of the shock defined in one of the "...res values" excel sheets, otherwise if the shock is not identified yet name it 'CorrelShock'
s.strctident.CorrelInstrument=''; % provide the IV variable in excel sheet "IV"
s.Feval=0;
s.CFt=1;
s.Fstartdate='2014q1';
s.Fenddate='2016q4';
s.Fendsmpl=0;
s.hstep=1;
s.window_size=0;
s.evaluation_size=0.5;
s.cband=0.95;
s.IRFband=0.68;
s.Fband=0.68;
s.FEVDband=0.68;
s.HDband=0.68;

author: european-central-bank

repo: BEAR-toolbox


text:
function J=mtfftc(data,tapers,nfft,Fs)
% Multi-taper fourier transform - continuous data
%
% Usage:
% J=mtfftc(data,tapers,nfft,Fs) - all arguments required
% Input:
%       data (in form samples x channels/trials or a single vector)
%       tapers (precalculated tapers from dpss)
%       nfft (length of padded data)
%       Fs   (sampling frequency)
%
% Output:
%       J (fft in form frequency index x taper index x channels/trials)
if nargin < 4; error('Need all input arguments'); end;
data=change_row_to_column(data);
[NC,C]=size(data); % size of data
[NK K]=size(tapers); % size of tapers
if NK~=NC; error('length of tapers is incompatible with length of data'); end;
tapers=tapers(:,:,ones(1,C)); % add channel indices to tapers
data=data(:,:,ones(1,K)); % add taper indices to data
data=permute(data,[1 3 2]); % reshape data to get dimensions to match those of tapers
data_proj=data.*tapers; % product of data with tapers
J=fft(data_proj,nfft)/Fs;   % fft of projected data

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [p1ext p2ext degen] = extline(p1, p2, imgwidth, imgheight)
% p1: [x y]
% p2: [x y]

dir = p2 - p1;

% intersection with top
% p1(1) + alpha*dir(2) = 1
% int1 = p1 + alpha*dir
alpha = (1 - p1(2)) / dir(2);
intp{1} = p1 + alpha*dir;
intp{1}(2) = 1; % numerical precision issues....
intp{1} = intp{1}(:)'; % [2x1]

% intersection with bottom
% p1(1) + alpha*dir(2) = imgheight
% int2 = p1 + alpha*dir
alpha = (imgheight - p1(2)) / dir(2);
intp{2} = p1 + alpha*dir;
intp{2}(2) = imgheight; % numerical precision issues....
intp{2} = intp{2}(:)'; % [2x1]

% intersection with left
% p1(2) + alpha*dir(1) = 1
% int3 = p1 + alpha*dir
alpha = (1 - p1(1)) / dir(1);
intp{3} = p1 + alpha*dir;
intp{3}(1) = 1; % numerical precision issues....
intp{3} = intp{3}(:)'; % [2x1]

% intersection with right
% p1(2) + alpha*dir(2) = imgwidth
% int4 = p1 + alpha*dir
alpha = (imgwidth - p1(1)) / dir(1);
intp{4} = p1 + alpha*dir;
intp{4}(1) = imgwidth; % numerical precision issues....
intp{4} = intp{4}(:)'; % [2x1]

b(1) = isinimage(intp{1}, imgwidth, imgheight);
b(2) = isinimage(intp{2}, imgwidth, imgheight);
b(3) = isinimage(intp{3}, imgwidth, imgheight);
b(4) = isinimage(intp{4}, imgwidth, imgheight);

% should touch at least 2 edges
% greater than 2 when intersects exactly at corner of image
% assert(b(1)+b(2)+b(3)+b(4) >= 2);
if b(1)+b(2)+b(3)+b(4) < 2
    degen = 1;
    p1ext = [0 0];
    p2ext = [0 0];
    return;
else
    degen = 0;
end

% b = find(b);
% p1ext = intp{b(1)};
% p2ext = intp{b(2)};

p = cat(1, intp{logical(b)});
% p = unique(p, 'rows');
% assert(size(p,1)==2);
p1ext = p(1,:);
p2ext = p(2,:);

%%
function flag = isinimage(p, imgwidth, imgheight)
if p(1)>=1 && p(1)<=imgwidth && p(2)>=1 && p(2)<=imgheight
    flag = 1;
else
    flag = 0;
end

author: zouchuhang

repo: LayoutNet


text:
clear all; close all; clc

x=0.2:0.1:5;
y=1./(x);

x2=0.2:0.1:5; n=length(x2);
y2=1./(x2) + 0.5*randn(1,n);

y3=y2+2*rand(5,1)+1

fill([0.5 1.4 1.4 0.5],[0.4 0.4 2 2],[0.8 0.8 0.8])
hold on

plot(x,y,'k','Linewidth',[2])
axis([0.2 4 0 5.5])
set(gca,'Xtick',[],'Ytick',[])

plot(x2,y3,'o','Linewidth',[1],'MarkerEdgeColor','k','MarkerFaceColor',[0 1 0.2],'MarkerSize',8)
plot(x2,y2,'o','Linewidth',[1],'MarkerEdgeColor','k','MarkerFaceColor',[0.9 0 1],'MarkerSize',8)

author: dynamicslab


text:
function [Int, Noise] = regEstFilIntGrad(inp, PbyPflag, lpf);
% regEstFilIntGrad - Estimates the intensity gradient, using local mean
%
%    [Int, Noise] = regEstFilIntGrad(inp, <PbyPflag>, <lpf>);
%
% Inputs:
%  inp - input inplanes affected by the intensity gradient
%  PbyPflag - operates plane by plane if activated (default 0)
%  lpf - low pass filter (applied separably to x,y,z) used to
%        compute the local mean
%
% Outputs:
%  Int   - Estimated intensity
%  Noise - Estimated power-spatial distribution of the noise, as the
%          local variance
%
% Oscar Nestares - 5/99
%

% default low-pass filter
if ~exist('lpf')
  lpf = conv([1 4 6 4 1]/16, [1 4 6 4 1]/16);
end

lpfZ = lpf;
if exist('PbyPflag')
   if PbyPflag
      lpfZ = 1;
   end
end

B = (length(lpf)-1)/2;

% adding border to the original inplanes
for k=1:size(inp,3)
   inpB(:,:,k) = regPutBorde(inp(:,:,k), B, B, 2);
end

% estimates the intensity as the local mean
Int = regConvXYZsep(inpB, 'repeat', lpf, lpf, lpfZ);

% estimates the noise as the mean local variance
for k=1:size(Int,3)     % adding border to the estimated intensity
   IntB(:,:,k) = regPutBorde(Int(:,:,k), B, B, 2);
end
Noise = regConvXYZsep((inpB-IntB).^2, 'repeat', lpf, lpf, lpfZ);


text:
function WindowButtonDownFcnCallback(a, b)

global h configFile selections current files imgId;

buttonPressed = get(h, 'SelectionType');

if(strcmp(buttonPressed, 'normal'))
  %Left click
  selections = [selections; current];
elseif(strcmp(buttonPressed, 'alt'))
  %Right click
  imgId = imgId + 1;
  if(imgId > numel(files))
    imgId = 1;
  end
  InteractiveCenterDemo(configFile);
end

end

author: Cloud-CV

repo: object-proposals