%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: various distances and Multi-Level Parametric Image Registration
%
%   - data                 PETCT, Omega=(0,140)x(0,151), level=4:7, m=[128,128]
%   - viewer               viewImage2D
%   - interpolation        splineInter
%   - distance             {'SSD','NCC','MI','NGF'}
%   - transformation       affine2D
% see also E7_PETCT_MLPIR_ext
%==============================================================================

function E7(example)

if nargin == 0, example = 'PETCT';  end;
imagePath  = fullfile(FAIRpath,'temp',example);
switch example
  case 'HNSP'
    setup2DHNSPData;
    theta = 1;
  case 'MRIhead',
    setup3DMRIData;
    theta = 1;
  case 'PETCT',
    setup2DPETCTData;
    theta = 0;
  otherwise, return;
end
level = 6; omega = ML{level}.omega; m = ML{level}.m;

DM = {'SSD','NCC','MImex','NGFdot'};
viewImage('set',viewPara{:})
filename = fullfile(FAIRpath,'temp',sprintf('%s-%s-%s.mat',mfilename,'rotation',example));

% 2e1 MRIhead
imgModel('reset','imgModel','splineInter','regularizer','moments','theta',theta);
[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega);
xc    = getCellCenteredGrid(omega,m);
xc    = reshape(xc,[],2);
Tin   = imgModel(T,omega,xc);
Rin   = imgModel(R,omega,xc);

%%  plots
k = figureh(1); clf; set(k,'color','w');
subplot(1,2,1); viewImage(Tin,omega,m);
subplot(1,2,2); viewImage(Rin,omega,m);
FAIRpause(1);

center = (omega(2:2:end)-omega(1:2:end))'/2;
trafo('reset','trafo','rotation2D','c',center);
trafo('w0');
if ~exist(filename,'file'),
  w = pi/4*linspace(-1,1,101)';
  save(filename,'w','DM');
else
  clear DM*
  load(filename)
end;

edge = 100;

for k=1:length(DM),
  variable = ['DM',DM{k}];
  var = whos('-file',filename);
  j = find(strcmp({var(:).name},variable)==1);

  if isempty(j),
    fprintf('============== %s ====================\n\n',variable)
    disp([variable,'=dm;']);
    dm = zeros(size(w));
    for j=1:length(w),
      Y  = trafo(w(j),xc(:));
      Tc = imgModel(T,omega,Y);
      dm(j) = feval(DM{k},Tc,Rin,omega,m,'edge',edge);
      if j== 1,
        figureh(3);
        ph = viewImage(Tc,omega,m);
      else
        set(ph,'cdata',reshape(Tc,m)'); drawnow
        FAIRpause(1/100)
      end;
      title(sprintf('%d/%d:%d/%d',k,length(DM),j,numel(w)));
    end;
    eval([variable,'=dm;']);
    save(filename,'-append',variable);
  end;
end;

%%
load(filename)

filename = fullfile(FAIRpath,'temp',...
  sprintf('%s-%s-%s.mat',mfilename,'rotation',example));

for k=1:4,
  variable = ['DM',DM{k}];
  eval(['dm=',variable,';']);
  dm = dm(1:length(w));
  [wOpt,j] = min(dm);
  figureh(k); clf; set(k,'color','w');
  ph = plot(w,dm,'-',w(j),dm(j),'*');
  set(ph,'linewidth',2,'color','k','markersize',20);
  set(gca,'fontsize',30);
  a = max(dm)-0.2*(max(dm)-min(dm));
  th = text(w(j),a,['$w^*=',num2str(w(j)),'$']);
  set(th,'fontsize',30,'interpreter','latex','horizontalalignment','center');
end;

%% finish
filename = fullfile(FAIRpath,'temp',...
  sprintf('%s-%s-%s.mat',mfilename,'translation',example));

trafo('reset','trafo','translation2D');
if ~exist(filename,'file'),
  [w1,w2] = ndgrid(0.1*(omega(2)-omega(1))*linspace(-1,1,21),...
                   0.2*(omega(4)-omega(3))*linspace(-1,1,21));
  save(filename,'w1','w2','DM');
else
  clear DM*
  load(filename)
end;

for k=1:length(DM),
  variable = ['DM',DM{k}];
  disp([variable,'=dm;']);

  var = whos('-file',filename);
  j = find(strcmp({var(:).name},variable)==1);

  if isempty(j),
    fprintf('============== %s ====================\n\n',variable)
    dm = zeros(size(w1));
    for j=1:numel(w1),
      Y = trafo([w1(j);w2(j)],xc(:));
      Tc = imgModel(T,omega,Y);
      dm(j) = feval(DM{k},Tc,Rin,omega,m,'edge',edge);
      if j== 1,
        figureh(3);
        ph = viewImage(Tc,omega,m);
      else
        set(ph,'cdata',reshape(Tc,m)'); drawnow
        pause(1/100)
      end;
      title(sprintf('%d/%d:%d/%d',k,length(DM),j,numel(w1)));
    end;
    eval([variable,'=dm;']);
    save(filename,'-append',variable);
  end;
end;

load(filename)

%figureh(1); close(1); figureh(1); clf; set(1,'position',position(800),'color','w');
%figureh(2); close(2); figureh(2); clf; set(2,'position',position(800),'color','w');
shift = [3000,0.5,-1,0.08];
Ztick = {[0:1000:4000],[0:0.1:0.3],[-1:0:1],[0:0.02:0.06]};

for k=1:4,

  variable = ['DM',DM{k}];
  eval(['dm=',variable,';']);

  fig = figureh(10+k);
  if ~isnumeric(fig), fig = fig.Number; end;
  close(fig); figureh(fig); set(fig,'color','w');
  ph=mesh(w1,w2,dm-shift(k)); grid off;
  set(gca,'fontsize',30);
  view(-135,25);
  set(ph,'linewidth',2)

  fig = figureh(20+k);
  if ~isnumeric(fig), fig = fig.Number; end;
  close(fig); figureh(fig); clf; set(fig,'color','w');
  contour(w1,w2,dm,10,'linewidth',2);
  set(gca,'fontsize',30);
end;
%==============================================================================


text:
function [ap, cmc] = compute_AP(good_image, junk_image, index)

cmc = zeros(length(index), 1);
ngood = length(good_image);

old_recall = 0;
old_precision = 1.0;
ap = 0;
intersect_size = 0;
j = 0;
good_now = 0;
njunk = 0;
for n = 1:length(index)
    flag = 0;
    if ~isempty(find(good_image == index(n), 1))
        cmc(n-njunk:end) = 1;
        flag = 1; % good image
        good_now = good_now+1;
    end
    if ~isempty(find(junk_image == index(n), 1))
        njunk = njunk + 1;
        continue; % junk image
    end

    if flag == 1%good
        intersect_size = intersect_size + 1;
    end
    recall = intersect_size/ngood;
    precision = intersect_size/(j + 1);
    ap = ap + (recall - old_recall)*((old_precision+precision)/2);
    old_recall = recall;
    old_precision = precision;
    j = j+1;

    if good_now == ngood
        return;
    end
end

end

author: JDAI-CV

repo: VeRidataset

save_path: github-repos/MATLAB/JDAI-CV-VeRidataset


text:
function [positionDecodingGLM_gaussian] = bz_positionDecodingGLM_gauss(varargin)
% USAGE
%   [] = bz_positionDecodingGLM_gauss()
%
% INPUTS
%
%   spikes
%
%   behavior
%
%   lfp
%
%   smoothingRange
%
%   plotting   - cell UID to plot
%
%   saveMat
%
% OUTPUTS
%
%   positionDecodingGLM_gaussian     cellinfo struct with the following fields
%                      .UID
%                      .region
%                      .sessionName
%                      .results{1:nCells}
%                      .tau
%                      .condition
%
%
%
%
% SEE
%
% written by david tingley, 2017

%% parse inputs
p = inputParser();
addRequired(p,'spikes',@isstruct);
addRequired(p,'behavior',@isstruct);
addRequired(p,'lfp',@isstruct);
addRequired(p,'smoothingRange',@isvector);

addParameter(p,'plotting',true,@isnumeric);
addParameter(p,'saveMat',[],@islogical);
parse(p,varargin{:})

spikes = p.Results.spikes;
behavior = p.Results.behavior;
lfp = p.Results.lfp;
smoothingRange = p.Results.smoothingRange;
plotting = p.Results.plotting;
saveMat = p.Results.saveMat;

%% set up data format and cellinfo struct that will be returned
disp('initializing data structs and formatting input data...')

positionDecodingGLM_gaussian.UID = spikes.UID;
positionDecodingGLM_gaussian.region = spikes.region;
positionDecodingGLM_gaussian.sessionName = spikes.sessionName; % session name

conditions = unique(behavior.events.trialConditions);
nCells = length(spikes.times);
positionSamplingRate = behavior.samplingRate;

% find a better way to get spike phase relationship...
[firingMaps] = bz_firingMap1D(spikes,behavior,lfp,5);

% iterate through conditions and compile spike trains and spike-phase
% trains
for cond = conditions
    trials = find(behavior.events.trialConditions==cond);
    intervals = behavior.events.trialIntervals(trials,:);
    for t = 1:length(trials)
        trial = trials(t);
        spk_trains{cond}{t} = zeros(nCells,ceil((intervals(t,2)-intervals(t,1))*1000)); % assumes intervals are in seconds, rounds to nearest millisecond
        phase_trains{cond}{t} = zeros(nCells,ceil((intervals(t,2)-intervals(t,1))*1000));
        for cell = 1:nCells
            if ~isempty(firingMaps.phaseMaps{cond}{cell})
                f = find(firingMaps.phaseMaps{cond}{cell}(:,2)==t);
                if ~isempty(f)
                for s=1:length(f)
                    phase_trains{cond}{t}(cell,ceil(firingMaps.phaseMaps{cond}{cell}(f(s),5)*1000)) = ...
                        firingMaps.phaseMaps{cond}{cell}(f(s),end);
                end
                end
            end
            sp = find(InIntervals(spikes.times{cell},intervals(t,:)));
            if ~isempty(sp)
                spks = Restrict(ceil((spikes.times{cell}(sp)-intervals(t,1))*1000+.0000001),[1 size(spk_trains{cond}{t},2)]);
                spk_trains{cond}{t}(cell,spks)=1;
            end
        end

%         position{cond}{t} = interp1(1:length(behavior.events.trials{trial}.x)...
%             ,behavior.events.trials{trial}.mapping,1:positionSamplingRate/1000:length(...
%             behavior.events.trials{trial}.x));
        nBins = size(spk_trains{cond}{t},2);
        nPos = length(behavior.events.trials{trial}.x);
        position{cond}{t} = interp1(1:length(behavior.events.trials{trial}.x)...
            ,behavior.events.trials{trial}.mapping,1:(nPos-1)/nBins:length(...
            behavior.events.trials{trial}.x)); % the -1 gaurantees the length to be longer than the above spk/phase trains
        position{cond}{t} = position{cond}{t}(1:length(spk_trains{cond}{t}));

    end
end

% collapse across trials..
for cond = conditions
%    phase_trains{cond} = cell2mat(phase_trains{cond});
%    spk_trains{cond} = cell2mat(spk_trains{cond});
% this was a bad idea because it lead to smoothing across trial boundaries,
% fixed below by smoothing first, then concatenating trials...
   position{cond} = cell2mat(position{cond});
end

%% set up GLM
for cell=1:nCells  % initialize table for data
    positionDecodingGLM_gaussian.results{cell} = table;
end
disp('running models...')
for cond = conditions
    warning off
    for wind = smoothingRange
       for cell = 80%1:nCells
            %smoothing
            phase_trains_smooth=[];
            cos_phase_trains_smooth=[];
            sin_phase_trains_smooth=[];
            rates_trains_smooth = [];
            for t = 1:length(phase_trains{cond})
%                 phase_trains_smooth=[phase_trains_smooth;...
%                     circ_smoothTS(phase_trains{cond}{t}(cell,:),wind,'method','mean','exclude',0)];
                phase_trains_smooth=[phase_trains_smooth;...
                     gaussfilt(1:length(phase_trains{cond}{t}(cell,:)),phase_trains{cond}{t}(cell,:),wind)'];
                cos_phase_trains_smooth=[cos_phase_trains_smooth;...
                     gaussfilt(1:length(phase_trains{cond}{t}(cell,:)),cos(phase_trains{cond}{t}(cell,:)),wind)'];
                sin_phase_trains_smooth=[sin_phase_trains_smooth;...
                     gaussfilt(1:length(phase_trains{cond}{t}(cell,:)),sin(phase_trains{cond}{t}(cell,:)),wind)'];
                rates_trains_smooth = [rates_trains_smooth; ...
                                       gaussfilt(1:length(phase_trains{cond}{t}(cell,:)),spk_trains{cond}{t}(cell,:),wind)'];
            end
            % phase coding
            [b dev stats] = glmfit(phase_trains_smooth',position{cond}','normal');
            yfit = glmval(b,phase_trains_smooth,'identity');
            struct.mse_phase = mean((yfit-position{cond}').^2); % mean squared error for rate code
            struct.mse_phase_pval = stats.p';

            % phase coding
            [b dev stats] = glmfit(cos_phase_trains_smooth',position{cond}','normal');
            yfit = glmval(b,cos_phase_trains_smooth,'identity');
            struct.mse_phase_cos = mean((yfit-position{cond}').^2); % mean squared error for rate code
            struct.mse_phase_cos_pval = stats.p';

            % phase coding
            [b dev stats] = glmfit(sin_phase_trains_smooth',position{cond}','normal');
            yfit = glmval(b,sin_phase_trains_smooth,'identity');
            struct.mse_phase_sin = mean((yfit-position{cond}').^2); % mean squared error for rate code
            struct.mse_phase_sin_pval = stats.p';

            % phase coding all
            [b dev stats] = glmfit([sin_phase_trains_smooth,cos_phase_trains_smooth],position{cond}','normal');
            yfit = glmval(b,[sin_phase_trains_smooth,cos_phase_trains_smooth],'identity');
            struct.mse_phase_all = mean((yfit-position{cond}').^2); % mean squared error for rate code
            struct.mse_phase_all_pval = stats.p';

            % rate coding
            [b dev stats] = glmfit(rates_trains_smooth',position{cond}','normal');
            yfit = glmval(b,rates_trains_smooth,'identity');
            struct.mse_rate = mean((yfit-position{cond}').^2);  % mean squared error for rate code
            struct.mse_rate_pval = stats.p';

            % extra variables to save
            struct.dfe = stats.dfe;
            struct.tau = wind;
            struct.condition = cond;
            positionDecodingGLM_gaussian.results{cell} = ...
                [positionDecodingGLM_gaussian.results{cell};struct2table(struct)];

            if ~isempty(plotting) & cell == plotting  % can only display one neuron for now..
%                 figure(cond);
                subplot(2,2,1)
                title('GLM decoding of pos, r-rate, g-phase')
                rows = positionDecodingGLM_gaussian.results{plotting}.condition == cond;
                plot(positionDecodingGLM_gaussian.results{plotting}.tau(rows),...
                    positionDecodingGLM_gaussian.results{plotting}.mse_rate(rows),'r')
                hold on
                plot(positionDecodingGLM_gaussian.results{plotting}.tau(rows),...
                    positionDecodingGLM_gaussian.results{plotting}.mse_phase_all(rows),'g')
                hold off
                subplot(2,2,2)
                plot(positionDecodingGLM_gaussian.results{plotting}.tau(rows),...
                    positionDecodingGLM_gaussian.results{plotting}.mse_phase_cos(rows),'g')
                subplot(2,2,3)
                plot(positionDecodingGLM_gaussian.results{plotting}.tau(rows),...
                    positionDecodingGLM_gaussian.results{plotting}.mse_phase_sin(rows),'g')
                subplot(2,2,4)
                plot(positionDecodingGLM_gaussian.results{plotting}.tau(rows),...
                    positionDecodingGLM_gaussian.results{plotting}.mse_phase(rows),'g')
                pause(.0001)
            end
       end
       disp(['finished with wind: ' num2str(wind) ' out of ' num2str(smoothingRange(end)) ' total'])
    end
    positionDecodingGLM_gaussian.dateRun = date;  % this can take a very long time so lets save each loop...
    if saveMat
        save([spikes.sessionName '.positionDecodingGLM_gaussian.cellinfo.mat'],'positionDecodingGLM_gaussian')
    end
    disp(['finished with condition: ' num2str(cond) ' out of ' num2str(length(conditions)) ' total'])
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function K = Abrarov2010CPC(x,y,tau,N)
% This is an approximation of the Voigt function within the Humlicek
% regions 3 and 4. The approximation is one given by S.M. Abrarov et. al.
% "High-accurace approximation of the complex probability function by
% Fourier expansion of exponential multiplier" (2010).
% x = sqrt(ln(2))*(nu - nu0)/alphaD
% y = sqrt(ln(2))*alphaL/alphaD
% where 'ln' denotes the natural log, nu the wavenumber, nu0 the wavenumber
% at center, alphaD and alphaL the Doppler and Lorentzian half-width at
% half-maximum.
% Suggested values for N & tau are 23 and 12 respectively.

if nargin == 2
    tau = 12;
    N = 23;
end

K = zeros(size(x));
a = zeros(1,N);
for c = 1:length(x)
    summation = 0;
    for n = 0:N
        a(n+1) = 2*sqrt(pi)/tau*exp(-(n*pi/tau)^2);
        first = ((1i*n*pi*tau+tau^2*y)*(1-exp(-(1i*n*pi+tau*y))*cos(tau*x(c))) + exp(-(1i*n*pi+tau*y))*tau^2*x(c)*sin(tau*x(c)))/(tau^2*x(c)^2 - (n*pi - 1i*tau*y)^2);
        second = ((1i*n*pi*tau-tau^2*y)*(1-exp(1i*n*pi-tau*y)*cos(tau*x(c))) - exp(1i*n*pi-tau*y)*tau^2*x(c)*sin(tau*x(c)))/(tau^2*x(c)^2 - (n*pi + 1i*tau*y)^2);
        summation = summation + a(n+1)*(first - second);
    end
    third = (y-exp(-(tau*y))*(y*cos(tau*x(c))-x(c)*sin(tau*x(c))))/(2*sqrt(pi)*(x(c)^2+y^2));
    K(c) = summation/(2*sqrt(pi)) - a(1)*third;
end

author: Sable

repo: mcbench-benchmarks


text:
function evaluatePCKh(predidxs,bSave)
% implementation of PCKh measure,
% as defined in [Andriluka et al., CVPR'14]

fprintf('evaluatePCKh()\n');

if (nargin < 2)
    bSave = true;
end

range = 0:0.01:0.5;

tableDir = './latex'; if (~exist(tableDir,'dir')), mkdir(tableDir); end
plotsDir = './plots'; if (~exist(plotsDir,'dir')), mkdir(plotsDir); end
tableTex = cell(length(predidxs)+1,1);

% load ground truth
p = getExpParams(-1);
load([p.gtDir '/annolist_dataset_v12'],'annolist');  % do not have this file. so the author does not release ground truth?
load([p.gtDir '/mpii_human_pose_v1_u12'],'RELEASE');

annolist_test = annolist(RELEASE.img_train == 0);
% evaluate on the "single person" subset only
single_person_test = RELEASE.single_person(RELEASE.img_train == 0);
% convert to annotation list with a single pose per entry
[annolist_test_flat, single_person_test_flat] = flatten_annolist(annolist_test,single_person_test);
% represent ground truth as a matrix 2x14xN_images
gt = annolist2matrix(annolist_test_flat(single_person_test_flat == 1));
% compute head size
headSize = getHeadSizeAll(annolist_test_flat(single_person_test_flat == 1));

pckAll = zeros(length(range),16,length(predidxs));

for i = 1:length(predidxs);

    % load predictions
    p = getExpParams(predidxs(i));
    load([fileparts(p.predFilename) '/pred_keypoints_mpii' ],'pred');
    assert(length(annolist_test) == length(pred));
    pred_flat = flatten_annolist(pred,single_person_test);
    pred = annolist2matrix(pred_flat(single_person_test_flat == 1));
    % only gt is allowed to have NaN
    pred(isnan(pred)) = inf;

    % compute distance to ground truth joints
    dist = getDistPCKh(pred,gt,headSize);

    % compute PCKh
    pck = computePCK(dist,range);

    % plot results
    [row, header] = genTablePCK(pck(end,:),p.name);
    tableTex{1} = header;
    tableTex{i+1} = row;

    pckAll(:,:,i) = pck;

    auc = area_under_curve(scale01(range),pck(:,end));
    fprintf('%s, AUC: %1.1f\n',p.name,auc);
end

if (bSave)
    fid = fopen([tableDir '/pckh.tex'],'wt');assert(fid ~= -1);
    for i=1:length(tableTex),fprintf(fid,'%s\n',tableTex{i}); end; fclose(fid);
end

% plot curves
plotCurve(squeeze(pckAll(:,end,:)),range,predidxs,'PCKh total, MPII',[plotsDir '/pckh-total-mpii'],bSave,range(1:5:end));
plotCurve(squeeze(mean(pckAll(:,[1 6],:),2)),range,predidxs,'PCKh ankle, MPII',[plotsDir '/pckh-ankle-mpii'],bSave,range(1:5:end));
plotCurve(squeeze(mean(pckAll(:,[2 5],:),2)),range,predidxs,'PCKh knee, MPII',[plotsDir '/pckh-knee-mpii'],bSave,range(1:5:end));
plotCurve(squeeze(mean(pckAll(:,[3 4],:),2)),range,predidxs,'PCKh hip, MPII',[plotsDir '/pckh-hip-mpii'],bSave,range(1:5:end));
plotCurve(squeeze(mean(pckAll(:,[7 12],:),2)),range,predidxs,'PCKh wrist, MPII',[plotsDir '/pckh-wrist-mpii'],bSave,range(1:5:end));
plotCurve(squeeze(mean(pckAll(:,[8 11],:),2)),range,predidxs,'PCKh elbow, MPII',[plotsDir '/pckh-elbow-mpii'],bSave,range(1:5:end));
plotCurve(squeeze(mean(pckAll(:,[9 10],:),2)),range,predidxs,'PCKh shoulder, MPII',[plotsDir '/pckh-shoulder-mpii'],bSave,range(1:5:end));
plotCurve(squeeze(mean(pckAll(:,[13 14],:),2)),range,predidxs,'PCKh head, MPII',[plotsDir '/pckh-head-mpii'],bSave,range(1:5:end));

end

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
function [g1, g2] = ggwhiteXggwhiteKernGradient(ggwhiteKern1, ggwhiteKern2, ...
    x, x2, covGrad)

% GGWHITEXGGWHITEKERNGRADIENT Compute a cross gradient between two GG WHITE kernels.
% FORMAT
% DESC computes cross gradient of parameters of a cross kernel between two
%	gg white kernels for the multiple output kernel.
% RETURN g1 : gradient of the parameters of the first kernel, for ordering
%	   see ggwhiteKernExtractParam.
% RETURN g2 : gradient of the parameters of the second kernel, for ordering
%	   see ggwhiteKernExtractParam.
% ARG ggwhiteKern1 : the kernel structure associated with the first GG white
%	   kernel.
% ARG ggwhiteKern2 : the kernel structure associated with the second GG white
%	   kernel.
% ARG x : inputs for which kernel is to be computed.
% ARG covgrad : gradient of the objective function with respect to the
%	   elements of the cross kernel matrix.
%
% FORMAT
% DESC computes cross kernel terms between two GG WHITE kernels for the multiple
%	output kernel.
% RETURN g1 : gradient of the parameters of the first kernel, for ordering
%	   see ggwhiteKernExtractParam.
% RETURN g2 : gradient of the parameters of the second kernel, for ordering
%	   see ggwhiteKernExtractParam.
% ARG ggwhiteKern1 : the kernel structure associated with the first GG WHITE
%	   kernel.
% ARG ggwhiteKern2 : the kernel structure associated with the second GG WHITE
%	   kernel.
% ARG x : row inputs for which kernel is to be computed.
% ARG x2 : column inputs for which kernel is to be computed.
% ARG covgrad : gradient of the objective function with respect to the
%	   elements of the cross kernel matrix.
%
% SEEALSO : multiKernParamInit, multiKernCompute, ggwhiteKernParamInit,
% ggwhiteKernExtractParam
%
% COPYRIGHT : Mauricio A. Alvarez and Neil D. Lawrence, 2008
%
% MODIFICATIONS : Mauricio A. Alvarez, 2009.

% KERN

if nargin < 5
    covGrad = x2;
    x2 = x;
end

[K, Pinv, Lqrinv, Lsrinv, Kbase, factorNoise, factorVar1, factorVar2, dist] = ...
    ggwhiteXggwhiteKernCompute(ggwhiteKern1, ggwhiteKern2, x, x2);
P = 1./Pinv;

if ggwhiteKern1.isArd
    matGradLqr = zeros(ggwhiteKern1.inputDimension,1);
    matGradLsr = zeros(ggwhiteKern1.inputDimension,1);
    for i=1:ggwhiteKern1.inputDimension
        X = repmat(x(:,i),1, size(x2,1));
        X2 = repmat(x2(:,i)',size(x,1),1);
        X_X2 = (X - X2).*(X - X2);
        matGradLqr(i) = sum(sum(0.5*covGrad.*K.*...
            (-0.5*Lqrinv(i) + Lqrinv(i)*P(i)*Lqrinv(i) - Lqrinv(i)*P(i)*X_X2*P(i)*Lqrinv(i) )));
        matGradLsr(i) = sum(sum(0.5*covGrad.*K.*...
            (-0.5*Lsrinv(i) + Lsrinv(i)*P(i)*Lsrinv(i) - Lsrinv(i)*P(i)*X_X2*P(i)*Lsrinv(i) )));
    end
else
    matGradLqr = sum(sum(0.5*covGrad.*K.*( -0.5*ggwhiteKern1.inputDimension*Lqrinv...
        + ((P*Lqrinv)^2)*(ggwhiteKern1.inputDimension*Pinv - dist))    ));
    matGradLsr = sum(sum(0.5*covGrad.*K.*( -0.5*ggwhiteKern1.inputDimension*Lsrinv...
        + ((P*Lsrinv)^2)*(ggwhiteKern1.inputDimension*Pinv - dist))   ));
end
grad_sigma2Noise = factorNoise*sum(sum(covGrad.*Kbase));
grad_variance1 = factorVar1*sum(sum(covGrad.*Kbase));
grad_variance2 = factorVar2*sum(sum(covGrad.*Kbase));
% only pass the gradient with respect to the inverse width to one
% of the gradient vectors ... otherwise it is counted twice.
g1 = [matGradLqr(:)' grad_sigma2Noise grad_variance1];
%g1 = [matGradLqr(:)' 0 grad_variance1];
g2 = [matGradLsr(:)' 0 grad_variance2];

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function aa=lpcrf2aa(rf)
%LPCRF2AA Convert reflection coefficients to area function AA=(RF)
%The areas are normalised so that aa(p+2)=1: the effective area of the free air beyond the lips.
% aa(1) is the area of the glottis. This will be zero if rf(:,1)=1.

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: lpcrf2aa.m,v 1.4 2007/05/04 07:01:39 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

aa = fliplr(cumprod([ones(1,size(rf,1)); fliplr((1-rf)./(1+rf)).']).');

author: decouples


text:
% example_2 - inversion of a fractional F(s)
clear, close all
[t1,ft1]=INVLAP('1/(sqrt(s)*s)',0.01,5,200,6,40,20);
[t2,ft2]=INVLAP('(20.5+3.7343*s^1.15)/(21.5+3.7343*s^1.15+0.8*s^2.2+0.5*s^0.9)/s',0.01,5,200);
figure(4)
set(4,'color','white')
subplot(2,1,1)
plot(t1,ft1), grid on, zoom on
xlabel('t [s]'), ylabel('f(t)')
%title('
subplot(2,1,2)
plot(t2,ft2), grid on, zoom
xlabel('t [s]'), ylabel('f(t)')
title('step response of a fractional control system')


text:
classdef CustomSequence < IncrementalSequence

    properties (Access = public)
        factor
    end

    methods (Access = public)

        function obj = CustomSequence(x0,x1,nSteps,initialValue,finalValue)
            obj@IncrementalSequence(x0,x1,nSteps,initialValue,finalValue);
        end

    end

    methods (Access = protected)

        function generateAlphaSequence(obj)
            if obj.nSteps < 2
                x = x1;
            else
                iSteps = 0:obj.nSteps-1;
                x = 1-(1-iSteps/(obj.nSteps-1)).^(obj.factor);
                x = (x1-x0)*x + x0;
            end
            obj.alpha = x;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef TensorRepresentation < RepresentationDescriptor

    properties (Access = protected)
        dim
    end

    methods (Access = public)

        function d = getDimension(obj)
            d = obj.dim;
        end
    end

    methods (Access = protected)

        function obj = TensorRepresentation()
            obj.loadDimensionValue()
        end

        function loadRepresentationVariable(obj)
            obj.representation = 'tensor';
        end
    end

    methods (Access = private)

        function loadDimensionValue(obj)
            tensSize = obj.getTensorSize();
            obj.dim = tensSize(1);
        end

    end
end


text:
  function ob = Fatrix(dim, arg, varargin)
%|function ob = Fatrix(dim, arg, handles, options)
%|
%| WARNING: any new objects should use fatrix2 instead of this!
%|
%| Construct Fatrix object, which is a 'fake matrix' or 'function-based matrix'
%| object, designed to generalize matrices by representing any linear operator.
%| The caller can provide a variety of (overloaded) "methods" for this object,
%| particularly multiplication and transposed multiplication.
%| Using this "meta object" allows the object designer to focus on the key
%| methods, rather than reinventing matrix methods for each new object.
%| Basic methods include things like "size" and "disp" etc., whereas more
%| advanced methods included ob(:,1) or [ob1; ob2] among many others.
%| See Gblur.m for examples.
%|
%|	ob * x		multiplication
%|	ob' * x		transposed multiplication
%|
%| It is up to the caller to provide function handles for both!
%|
%| in
%|	dim [2]		Fatrix "dimensions"
%|	arg		arguments passed to handle functions
%|				(usually a cell or struct)
%|
%| handles (all optional): 'name1', handle1, 'name2', handle2, ...
%|	'forw'		forw(arg, x)	ob * x
%|	'back'		back(arg, x)	ob' * x
%|	'gram'		build_gram(ob, W, reuse), build G' * W * G
%|	'free'		free(arg)
%|	'abs'		ob = abs(ob): absolute value operation
%|	'ufun'		out = ufun(ob, varargin) : user-defined function
%|	'block_setup'	ob = block_setup(ob, varargin)
%|	'blockify_data'	cell_data = blockify_data(ob, array_data, varargin)
%|	'mtimes_block'	mtimes_block(arg, is_transpose, x, iblock, nblock)
%|	'power'		ob = power(ob, sup)	ob.^sup
%|
%| options
%|	'caller', string	name of calling routine ("meta class")
%|				default: name determined automatically.
%|	'cascade_after', thing		ob * x -> thing * forw(arg, x)
%|	'cascade_before', thing		ob * x -> forw(arg, thing * x)
%|		these "things" can also be function handles (see mtimes_block.m)
%|
%| out
%|	ob		Fatrix object
%|
%| Copyright 2004-6-29, Jeff Fessler, University of Michigan

% create default object, as required by Mathworks
ob.caller = '';
ob.arg = {};
ob.dim = [];
ob.is_transpose = false;
ob.is_subref = false;
ob.nblock = [];
ob.iblock = [];

% trick: default to some simple anonymous functions that will work if
% "arg" is anything that knows how to multiply, e.g., a matrix.
ob.handle_back = @(M,y) M' * y;
ob.handle_forw = @(M,x) M * x;
ob.handle_power = @(M,p) M .^ p;

ob.handle_ufun = [];
ob.handle_abs = [];
ob.handle_free = [];
ob.handle_gram = [];
ob.handle_mtimes_block = [];
ob.handle_block_setup = [];
ob.handle_blockify_data = [];
ob.cascade_after = [];
ob.cascade_before = [];

if nargin == 0 % required by Mathworks
	if nargout == 0, help(mfilename), end
	ob = class(ob, 'Fatrix');
return
end

if nargin < 2
	ir_usage
end

ob.arg = arg;
ob.dim = dim;

ob = vararg_pair(ob, varargin, 'subs', { ...
	'ufun', 'handle_ufun';
	'abs', 'handle_abs';
	'free', 'handle_free';
	'back', 'handle_back';
	'forw', 'handle_forw';
	'gram', 'handle_gram';
	'power', 'handle_power';
	'block_setup', 'handle_block_setup';
	'blockify_data', 'handle_blockify_data';
	'mtimes_block', 'handle_mtimes_block'});

if isempty(ob.caller)
	ob.caller = caller_name;
end

ob = class(ob, mfilename);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function x = gather_try(x)
try
    x = gather(x);
catch
end
end

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function [myhandle, x_circle, y_circle] = circles(radii,centers,varargin)
% Plot multiple circles as a single line object...very fast, very efficient!
%
% circles(radii,centers) plots circles of radius values specified by RADII,
%     at centers specified by CENTERS.
% circles(radii,centers,'PropertyName',propertyvalue,...)
%         Will accept as inputs any valid properties
%         accepted by PLOT.
% [h,x_circle,y_circle] = circles(...)
%         Also provides handle to (single-plot) circle
%         object, plotted x-values, and plotted y-values.
%
% NOTE: x- and y- values for each circle will be terminated
%       by NaNs in the output x_circle, y_circle vectors!
%
% INPUTS:
%     r        A vector of n radii
%     centers  An n x 2 array of x,y center coordinates
%
%     Optional Parameter-Value Pairs:
%        Resolution: Stepsize for plotting from theta = 0 to
%                    theta = 360 (Default = 1)
%        Any valid PV-pair accepted by PLOT
%
% EXAMPLES
%
%%% Example 1: Plot a single circle:
%
% r = 10; c = [3,5];
% circles(r,c,'color','g','linewidth',2)
%
%
%%% Example 2: Plot 200 random-sized and centered circles
%%% (Note: this is MUCH faster than drawing them one-by-one,
%%% but the speed gain is at the cost of getting a single
%%% line-object handle rather than a vector of them.)
%
% tic;
% n = 1000;
% r = rand(n,1)*100;
% c = rand(n,2)*300;
% [h,xs,ys] = circles(r,c,...
%         'color','r','linewidth',1,'resolution',1);
% axis equal
% t_allAtOnce = toc
%
%
%%% Example 3: (FOR COMPARISON ONLY!!!)
%%% (Using the r,c generated above:)
%
% tic;
% for ii = 1:numel(r)
%   circles(r(ii),c(ii,:),'color','r','linewidth',1,'resolution',1);
% end
% t_oneAtATime = toc
% fprintf('Plotting %d circles was %0.2f x faster as a single object!\n',n,t_oneAtATime/t_allAtOnce)
%
% Written by Brett Shoelson, PhD 8/3/98.
% brett.shoelson@mathworks.com
%
% MODIFICATIONS:
% Vectorized 11/20/01
% Arg for resolution 03/08/08
% 12/06/11 Modified as CIRCLES to accept multiple radii and centers,
%   and to vectorize all as a single line object.
%
% Copyright 2012 The MathWorks, Inc.

if nargin < 2
    error('CIRCLES requires a minimum of 2 input arguments.');
end

if rem(nargin,2) ~= 0
    % Odd number of inputs!
    error('CIRCLES: Parameter-Values must be entered in valid pairs')
end

[PVs,resolution] = parsePVs(varargin);

theta=0:resolution:360;

x_circle = bsxfun(@times,radii,cos(theta*pi/180));
x_circle = bsxfun(@plus,x_circle,centers(:,1));
x_circle = cat(2,x_circle,nan(size(x_circle,1),1));
x_circle =  x_circle';
x_circle = x_circle(:);

y_circle = bsxfun(@times,radii,sin(theta*pi/180));
y_circle = bsxfun(@plus,y_circle,centers(:,2));
y_circle = cat(2,y_circle,nan(size(y_circle,1),1));
y_circle =  y_circle';
y_circle = y_circle(:);

oldhold = ishold;
hold on;
myhandle = plot(x_circle,y_circle);

% Apply PV pairs
for ii = 1:2:numel(PVs)
    set(myhandle,PVs{ii},PVs{ii+1})
end

if ~oldhold,hold off;end
if nargout < 3
    clear y_circle
end
if nargout < 2
    clear x_circle
end
if nargout < 1
    clear myhandle
end

function [PVs,resolution] = parsePVs(PVs)
resolution = 1;
for ii = 1:2:numel(PVs)
    if strcmpi(PVs{ii},'resolution')
        resolution = PVs{ii+1};
        PVs = PVs(setdiff(1:numel(PVs),[ii,ii+1]));
        return
    end
end


text:
function res = adjD(y)

res = zeros(size(y,1),size(y,2));

%y1 = ones(imsize)*y(1)/sqrt(prod(imsize));
%yx = (reshape(y(2:prod(imsize)+1), imsize(1), imsize(2)));
%yy = (reshape(y(prod(imsize)+2:end), imsize(1), imsize(2)));

res = adjDx(y(:,:,1)) + adjDy(y(:,:,2));

return;

function res = adjDy(x)
res = x(:,[1,1:end-1]) - x;
res(:,1) = -x(:,1);
res(:,end) = x(:,end-1);

function res = adjDx(x)
res = x([1,1:end-1],:) - x;
res(1,:) = -x(1,:);
res(end,:) = x(end-1,:);

author: thomaskuestner


text:
function [x,resvec,state] = pcg_ccot(A,b,opts,M1,M2,ip,x0,state)

% This is a modified version of Matlab's pcg function, that performs
% preconditioned conjugate gradient.

% tol = opts.tol;
maxit = opts.maxit;

if ~isfield(opts, 'init_forget_factor')
    opts.init_forget_factor = 1;
end

if opts.debug
    n2b = sqrt(ip(b,b)); % Norm of rhs vector, b
end

existM1 = ((nargin >= 4) && ~isempty(M1));
existM2 = ((nargin >= 5) && ~isempty(M2));

x = x0;

% Load the CG state
p = [];
rho = 1;
r_prev = [];
load_state = nargin > 7 && ~isempty(state) && opts.init_forget_factor > 0;
if load_state
    if isfield(state, 'p')
        p = state.p;
    end
    if isfield(state, 'rho') && ~isempty(state.rho)
        rho = state.rho / opts.init_forget_factor;
    end
    if isfield(state, 'r_prev') && ~opts.CG_use_FR
        r_prev = state.r_prev;
    end
end

% Set up for the method
state.flag = 1;

% r = cellfun(@minus, b, iterapp('mtimes',afun,atype,afcnstr,x,varargin{:}), 'uniformoutput', false);
r = cellfun(@minus, b, A(x), 'uniformoutput', false);

if opts.debug
    normr = sqrt(ip(r,r));                   % Norm of residual
    normr_act = normr;
end

if opts.debug
    resvec = zeros(maxit+1,1);         % Preallocate vector for norm of residuals
    resvec(1,:) = normr;               % resvec(1) = norm(b-A*x0)
else
    resvec = [];
    relres = [];
end

% loop over maxit iterations (unless convergence or failure)

for ii = 1 : maxit
    if existM1
        y = M1(r);
    else % no preconditioner
        y = r;
    end

    if existM2
        z = M2(y);
    else % no preconditioner
        z = y;
    end

    rho1 = rho;
    rho = ip(r, z);
    if ((rho == 0) || isinf(rho))
        state.flag = 4;
        break
    end

    if (ii == 1 && isempty(p))
        p = z;
    else
        if opts.CG_use_FR
            % Use Fletcher-Reeves
            beta = rho / rho1;
        else
            % Use Polak-Ribiere
            rho2 = ip(r_prev, z);
            beta = (rho - rho2) / rho1;
        end
        if ((beta == 0) || isinf(beta))
            state.flag = 4;
            break
        end
        beta = max(0, beta);
        p = cellfun(@(z,p) z + beta * p, z, p, 'uniformoutput', false);
    end

    q = A(p);
    pq = ip(p, q);
    if ((pq <= 0) || isinf(pq))
        state.flag = 4;
        break
    else
        if opts.CG_standard_alpha
            alpha = rho / pq;
        else
            alpha = ip(p, r) / pq;
        end
    end
    if isinf(alpha)
        state.flag = 4;
        break
    end

    % Save old r if not using FR formula for beta
    if ~opts.CG_use_FR
        r_prev = r;
    end

    % form new iterate
    x = cellfun(@(x,p) x + alpha * p, x, p, 'uniformoutput', false);

    if ii < maxit || opts.debug
        r = cellfun(@(r,q) r - alpha * q, r, q, 'uniformoutput', false);
    end

    if opts.debug
        normr = sqrt(ip(r,r));
        normr_act = normr;
        resvec(ii+1,1) = normr;
    end
end

iter = ii;
if opts.debug
    relres = normr_act / n2b;
end

% truncate the zeros from resvec
if opts.debug
    if ((state.flag <= 1) || (state.flag == 3))
        resvec = resvec(1:ii+1,:);
    else
        resvec = resvec(1:ii,:);
    end
end

% Save the state
if nargout > 2
    state.p = p;
    state.rho = rho;
    if ~opts.CG_use_FR
        state.r_prev = r_prev;
    end
end

author: martin-danelljan

repo: ECO


text:
function sigma = massey_berlekamp_M2(n,k,t,S,field)

%http://www.ee.ucla.edu/~matache/rsc/node8.html#SECTION00051000000000000000

%Step 2: Initialize variables
kk = 0;

for i = 1:n
    Kappa(1,i) = -Inf;
end
Kappa(1,1) = 0;

%Kappa

LAMBDA = 0;
Tau = [-inf 0];

done = 0;

%Step 3:
while (done ~= 1)
    %disp('K');

    kk = kk + 1;

    %disp('S(kk)');
    %S(kk)

    %disp('LAMBDA')
    %LAMBDA

    sum = -Inf;
    for i = 1:LAMBDA
        %Kappa(kk,i+1)
        %S(kk-i)
        sum = gfadd(sum,gfmul(Kappa(kk,i+1),S(kk-i),field),field);
    end

    %disp('Delta - sum')
    %sum

    delta(kk) = gfadd(S(kk),sum,field);

    %disp('delta');
    %delta

    %Step 4:
    if (delta(kk) == -Inf)
        for i = 1:n
            Kappa(kk+1,i) = Kappa(kk,i);
        end
    end

    if (delta(kk) ~= -Inf)

        for i = 1:n
            Kappa_i(i) = Kappa(kk-1+1,i);
        end

        Kappa_k = gfadd(Kappa_i,gfconv(delta(kk),Tau,field),field);

        while length(Kappa_k) < n
            Kappa_k = [Kappa_k -Inf];
        end

        for i = 1:length(Kappa_k)
            Kappa(kk+1,i) = Kappa_k(i);
        end

        %Step 7:
        if (2*LAMBDA < kk)
            LAMBDA = kk - LAMBDA;

            for i = 1:n
                Kappa_k(i) = Kappa(kk+1-1,i);
            end

            Tau = gfconv(Kappa_k,gfdiv(0,delta(kk),field),field);
        end
    end

    %Step 8:
    Tau = gfconv([-Inf 0],Tau,field);

    %step 9:
    if kk >= 2*t
        done = 1;
    end

    %Kappa
    %LAMBDA
    %Tau

end

for i = 1:n
    sigma(i) = Kappa(kk+1,i);
end


text:
%  This callback routine has to be an m-file, not a function
%  because it modifies the workspace.  We know which function to
%  execute by examining the label of the current menu item.
%
%  Keith Rogers 11/30/93

SavedFile = get(get(gcm2,'Parent'),'UserData');

if (strcmp(get(gcm2,'Label'),'Load'))
    [filename,pathname] = uigetfile('*.mat',...
        'Choose a MATLAB Data file',50,50);
    if (filename ~= 0)
        SavedFile = [pathname filename];
        clear filename pathname;
        eval(['load ' SavedFile ';']);
        set(get(gcm2,'Parent'),'UserData',SavedFile);
    end
elseif (strcmp(get(gcm2,'Label'),'Save'))
    if (strcmp(SavedFile,''))
        [filename,pathname] = uiputfile('*.mat','Data Filename',...
            50,50);
        if (filename ~= 0)
            SavedFile = [pathname filename];
            clear filename pathname;
            eval(['save ' SavedFile ';']);
            set(get(gcm2,'Parent'),'UserData',SavedFile);
        end
    else
        eval(['save ' SavedFile ';']);
    end
elseif (strcmp(get(gcm2,'Label'),'Save As'))
    if (strcmp(SavedFile,''))
        [filename,pathname] = uiputfile('*.mat','Data Filename',...
            50,50);
    else
        [filename,pathname] = uiputfile(SavedFile,'Data Filename',...
            50,50);
    end
    if (filename ~= 0)
        SavedFile = [pathname filename];
        clear filename pathname;
        eval(['save ' SavedFile ';']);
        set(get(gcm2,'Parent'),'UserData',SavedFile);
    end
elseif (strcmp(get(gcm2,'Label'),'New Figure'))
    figure;
    makemenus;
end

author: CelsoReyes


text:
function [A,b,Aeq,beq,lb,ub,report] = separateBounds( A,b,Aeq,beq,tol )
%separateBounds - given a quadruplet {A, b, Aeq, beq} of linear constraint
%matrices and vectors corresponding to the region A*x<=b, Aeq*x=beq, the routine
%will convert to a sextuplet {A, b, Aeq, beq, lb, ub} of the kind used
%in the Optimization Toolbox.
%
%      [A,b,Aeq,beq,lb,ub,report] = separateBounds( A,b,Aeq,beq )
%
%In other words, the code will look through the rows of A and Aeq for
%constraints corresponding to simple upper and lower bounds. It will then
%separate these constraints from the rest, expressing them instead using
%vectors lb, ub. The "report" output argument is a structure containing
%some relevnt stats,
%
%    report.infeasibleBounds: indices, i, of any lb(i)>ub(i)
%    report.inequalitiesRemoved: indices of rows removed from A,b
%    report.equalitiesRemoved: indices of rows removed from Aeq,beq
%
%EXAMPLE:
%
%     A =  [1     1     1
%          -1     0     0
%           0    -1     0
%           1     0     0
%           0     1     0] ;
%
%      b=[1 0 0 1 1].';
%
%     Aeq=[0 0 1];
%     beq=1;
%
%     >>[A,b,Aeq,beq,lb,ub] = separateBounds( A,b,Aeq,beq )
%
%         A =
%
%              1     1     1
%
%
%         b =
%
%              1
%
%
%         Aeq =
%
%            Empty matrix: 0-by-3
%
%
%         beq =
%
%            Empty matrix: 1-by-0
%
%
%         lb =
%
%              0
%              0
%              1
%
%
%         ub =
%
%              1
%              1
%              1
%
%By default, the code will consider a constraint A(i,:)<=b(i) to correspond
%to a pure bound if a row A(i,:) or similarly Aeq(i,:) contains no more
%than one non-zero element. In certain situations, however, the matrices A
%and Aeq are the output of floating point calculations and such rows will
%contain non-zero elements representing numerical noise. In such instances,
%one can call the code with a tolerance parameter
%
%      [A,b,Aeq,beq,lb,ub,report] = separateBounds( A,b,Aeq,beq,tol )
%
%where 0<=tol<=1. The criterion for deciding whether a constraint is a pure
%bound is then,
%
%        max(abs(A(:,i)))>=tol*sum(abs(A(i,:)))
%
%The default beavior corresponds to tol=1.

  %%%%begin parsing

          if ~exist('Aeq','var'), Aeq=[]; end
          if ~exist('beq','var'),  beq=[]; end
          if ~exist('tol','var')||isempty(tol),
              tol=1;
          elseif tol<0 || tol>1
              error 'Tolerance parameter must satisfy 0<=tol<=1';
          end

          Na=nan; Nb=nan;
          if xor(isempty(A),isempty(b))
            error 'If A is empty then b must also be empty and vice versa'
          elseif ~isempty(A)
              Na=size(A,2);
          elseif isempty(A)
              A=[];b=[];
          end

          if xor(isempty(Aeq),isempty(beq))
              error 'If Aeq is empty then beq must also be empty and vice versa'
          elseif ~isempty(Aeq)
              Nb=size(Aeq,2);
          elseif isempty(Aeq)
                Aeq=[]; beq=[];
          end

          if ~isnan(Na) && ~isnan(Nb)

              if Na~=Nb
               error 'If both A and Aeq are both non-empty, they must have same number of columns'
              end

          elseif isnan(Na) && isnan(Nb)

              lb=[]; ub=[];
              return

          end

          N=max(Na,Nb);

          b=b(:);  beq=beq(:); %henceforth we are sure b-data are column vectors

           if size(A,1)~=length(b)
               error 'Incompatible inequality matrix data sizes: size(A,1) ~= length(b)'
           end

           if size(Aeq,1)~=length(beq)
               error 'Incompatible equality matrix data sizes: size(Aeq,1) ~= length(beq)'
           end

  %%%%end parsing

   [ls1,lv1,us1,uv1,A1,b1,rows1] = extract(A,b,tol);
   [ls2,lv2,us2,uv2,A2,b2,rows2] = extract(Aeq,beq,tol);
   [ls3,lv3,us3,uv3] = deal(us2,uv2,ls2,lv2);

   A=A1;
   b=b1;
   Aeq=A2;
   beq=b2;
   lb=accumarray([ls1;ls2;ls3]  ,  [lv1;lv2;lv3]  , [N,1], @max , -inf);
   ub=accumarray([us1;us2;us3]  ,  [uv1;uv2;uv3]  , [N,1], @min , +inf);

   report.infeasibleBounds=find(ub<lb).';
   report.inequalitiesRemoved=rows1.';
   report.equalitiesRemoved=rows2.';

   if nargout<7 &&  ~isempty(report.infeasibleBounds)

       warnstr=[sprintf('Infeasible combinations of bounds ( LB(i)>UB(i) ) were detected.\n\n'),...
               sprintf( 'To disable this warning, call with 7 output argument syntax:\n\n'),...
               sprintf( '          [A,b,Aeq,beq,lb,ub,report] = separateBounds(...)')];

       warning('separateBounds:infeasible', warnstr);

   end

 function [lsubs,lvals,usubs,uvals,A,b,rows] = extract(A,b,tol)

      lsubs=[]; usubs=[];
      lvals=[];  uvals=[];
      rows=[];

     if isempty(A), return; end

    b=b(:);

    absA=abs(A);

    [maxA,idx]=max(absA,[],2);
    sumA=sum(absA,2);

    rows =find(  maxA>=tol*sumA & sumA>0 );

    subs=idx(rows);

      ind=sub2ind(size(A),rows, subs);
      num=b(rows);
      den=A(ind);

    vals=num./den;

    lidx=(den<0);
    uidx=(den>0);

    lsubs=subs(lidx);
    lvals=vals(lidx);

    usubs=subs(uidx);
    uvals=vals(uidx);

   zrows=(sumA==0);
   rowsInfeas=find(zrows & b<0);
   rowsFeas=find(zrows & b>=0);

   if ~isempty(rowsInfeas) %row of all zeros

       N=size(A,2);
       lsubs=(1:N).';
       lvals=inf(N,1);
       usubs=lsubs;
       uvals=-lvals;
       rows=[rows;rowsInfeas];
   end

   if ~isempty(rowsFeas) %row of all zeros
       rows=[rows;rowsFeas];
   end

    A(rows,:)=[];
    b(rows)=[];


text:
% Wrapper for vl_nnscale2D block
% inputs{1} :   X     : 1 x 2 x n x b
% inputs{2} :   S     : 1 x 1 x 1 x b
% outputs{1}:   y     : 1 x 2 x n x b

classdef scale2D < dagnn.Layer
    methods
        function outputs = forward(~, inputs, ~)

            useGPU = isa(inputs{1}, 'gpuArray');
            if useGPU
                outputs{1} = gpuArray( vl_nnscale2D(gather(inputs{1}), gather(inputs{2})) );
            else
                outputs{1} = vl_nnscale2D(inputs{1}, inputs{2});
            end

        end

        function [derInputs, derParams] = backward(~, inputs, ~, derOutputs)

            useGPU = isa(inputs{1}, 'gpuArray');
            if useGPU
                [y,dsdy] =  vl_nnscale2D(gather(inputs{1}), gather(inputs{2}), gather(derOutputs{1}));
                derInputs = {gpuArray(y),gpuArray(dsdy)};
            else
                [y,dsdy] =  vl_nnscale2D(inputs{1}, inputs{2}, derOutputs{1});
                derInputs = {y,dsdy};
            end

            derParams = {};
        end

        function outputSizes = getOutputSizes(~, inputSizes)
            outputSizes = inputSizes{1};
        end

        function obj = scale2D(varargin)
            obj.load(varargin);
        end
    end
end

author: anilbas


text:
function [varargout] = likGaussHe(hyp, y, s, mu, s2, inf, i)

% likGauss - Gaussian likelihood function for regression. The expression for the
% likelihood is
%   likGauss(t) = exp(-(t-y)^2/2*sn^2) / sqrt(2*pi*sn^2),
% where y is the mean and sn is the standard deviation.
%
% The hyperparameters are:
%
% hyp = [  log(sn)  ]
%
% Several modes are provided, for computing likelihoods, derivatives and moments
% respectively, see likFunctions.m for the details. In general, care is taken
% to avoid numerical issues when the arguments are extreme.
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2015-07-13.
%                                      File automatically generated using noweb.
%
% See also LIKFUNCTIONS.M.

if nargin<3, varargout = {'1'}; return; end   % report number of hyperparameters

if isempty(s)
    sn2 = exp(2*hyp);
else
    sn2_base = exp(2*hyp);
    sn2 = sn2_base + s.^2;
end

if nargin<6                              % prediction mode if inf is not present
  if isempty(y);  y = zeros(size(mu)); end
  if isempty(s);  s = zeros(size(mu)); end
  s2zero = 1; if nargin>4&&numel(s2)>0&&norm(s2)>eps, s2zero = 0; end  % s2==0 ?
  if s2zero                                                    % log probability
    lp = -(y-mu).^2./sn2/2-log(2*pi*sn2)/2; s2 = 0;
  else
    lp = likGaussHe(hyp, y, s, mu, s2, 'infEP');                         % prediction
  end
  ymu = {}; ys2 = {};
  if nargout>1
    ymu = mu;                                                   % first y moment
    if nargout>2
      ys2 = s2 + sn2;                                          % second y moment
    end
  end
  varargout = {lp,ymu,ys2};
else
  switch inf
  case 'infLaplace'
    if nargin<6                                             % no derivative mode
      if isempty(y), y=0; end
      ymmu = y-mu; dlp = {}; d2lp = {}; d3lp = {};
      lp = -ymmu.^2/(2*sn2) - log(2*pi*sn2)/2;
      if nargout>1
        dlp = ymmu/sn2;                      % dlp, derivative of log likelihood
        if nargout>2                    % d2lp, 2nd derivative of log likelihood
          d2lp = -ones(size(ymmu))/sn2;
          if nargout>3                  % d3lp, 3rd derivative of log likelihood
            d3lp = zeros(size(ymmu));
          end
        end
      end
      varargout = {lp,dlp,d2lp,d3lp};
    else                                                       % derivative mode
      lp_dhyp = (y-mu).^2/sn2 - 1;  % derivative of log likelihood w.r.t. hypers
      dlp_dhyp = 2*(mu-y)/sn2;                               % first derivative,
      d2lp_dhyp = 2*ones(size(mu))/sn2;   % and also of the second mu derivative
      varargout = {lp_dhyp,dlp_dhyp,d2lp_dhyp};
    end

  case 'infEP'
    if nargin<7                                             % no derivative mode
      lZ = -(y-mu).^2./(sn2+s2)/2 - log(2*pi*(sn2+s2))/2;    % log part function
      dlZ = {}; d2lZ = {};
      if nargout>1
        dlZ  = (y-mu)./(sn2+s2);                    % 1st derivative w.r.t. mean
        if nargout>2
          d2lZ = -1./(sn2+s2);                      % 2nd derivative w.r.t. mean
        end
      end
      varargout = {lZ,dlZ,d2lZ};
    else                                                       % derivative mode
      dlZhyp = ((y-mu).^2./(sn2+s2)-1) ./ (1+s2./sn2);   % deriv. w.r.t. hyp.lik
      varargout = {dlZhyp};
    end

  case 'infVB'
    % variational lower site bound
    % t(s) = exp(-(y-s)^2/2sn2)/sqrt(2*pi*sn2)
    % the bound has the form: (b+z/ga)*f - f.^2/(2*ga) - h(ga)/2
    n = numel(s2); b = zeros(n,1); y = y.*ones(n,1); z = y;
    varargout = {b,z};
  end
end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function spm_MDP_trust
% Demo of active inference for trust games
%__________________________________________________________________________
%
% This routine uses the Markov decision process formulation of active
% inference (with variational Bayes) to model a simple trust game. In trust
% games, one plays an opponent who can either cooperate or defect. The
% payoff contingencies depend upon the joint choices of you and your
% opponent, which in turn depend upon your inferences about the nature of
% the opponent (pro-social or non-social). This example illustrates single
% round games with a special focus on Bayesian belief updating between
% games. This is illustrated in terms of evidence accumulation about
% the nature of the opponent by using the posterior marginal distributions
% following one game as the prior distribution over beliefs about the
% opponent in the next. This accumulation is shown in the final figures.
%
% In this example, there are nine states. The first is a starting state
% and the subsequent eight states model the four combinations of
% cooperation and defection (between you and your opponent) under the
% prior beliefs that the opponent is either pro-social or non-social.
% Initially, these prior beliefs are uninformative but are subsequently
% informed through experience. prior beliefs about behaviour are based on
% relative entropy or KL divergence in the usual way - which requires the
% specification of utility functions over states based upon standard payoff
% tables in these sorts of games. It is interesting to see how precision
% or confidence in beliefs about choices, fluctuates with beliefs about
% the nature of one's opponent.
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MDP_trust.m 6044 2014-06-14 10:22:46Z karl $

% set up and preliminaries
%==========================================================================
rng('default')
%
% Payoffs (reward):
% _________________________________________________________________________
%                                 Trustee
% _________________________________________________________________________
% Investor              Cooperate (fT high) Defect (fT low)
% _________________________________________________________________________
% Cooperate            A (e.g.=26)         C (e.g.= 10)
% (fI high)            a (e.g.=26)         b (e.g.= 42)
%
% Defect               B (e.g.=21)         D (e.g.= 18)
% (fI low)             c (e.g.= 7)         d (e.g.= 10)
% _________________________________________________________________________

U{1} = [26 10;
        21 18]/8;                   % self payoff (utility)

U{2} = [26 42;
         7 10]/8;                   % other payoff (utility)

% initial state - encoding the actual type of trustee
%--------------------------------------------------------------------------
S    = [0 1]';                      % indicator - [prosocial nonsocial]
S    = kron(S,[1 0 0 0 0]');

% prior beliefs about initial state
%--------------------------------------------------------------------------
k    = [1 1];
p    = spm_softmax(k(:));
k    = log(p);
D    = kron(p,[1 0 0 0 0]');

% investor's payoffs or prior beliefs (softmax(utility))
%--------------------------------------------------------------------------
a    = 1/2;
pp   = [0; spm_softmax(spm_vec((1 - a)*U{1} + a*U{2}))];
pn   = [0; spm_softmax(spm_vec(        U{1}         ))];

C    = [pp*p(1); pn*p(2)];

% investor's belief (based on a prosocial and nonsocial trustee)
%--------------------------------------------------------------------------
cp   = spm_softmax(((1 - a)*U{2}(1,:) + a*U{1}(1,:))');
dp   = spm_softmax(((1 - a)*U{2}(2,:) + a*U{1}(2,:))');
cn   = spm_softmax((        U{2}(1,:)              )');
dn   = spm_softmax((        U{2}(2,:)              )');

% transition probabilities (B{1} - (c)ooperate; B{2} - (d)efect)
%--------------------------------------------------------------------------
B{1} = ...                        % cooperate:
   [0     0 0 0 0  0 0 0 0 0;     % start - prosocial trustee
    cp(1) 1 0 0 0  0 0 0 0 0;     % cc - prosocial - state 2
    0     0 1 0 0  0 0 0 0 0;     % dc - prosocial - state 3
    cp(2) 0 0 1 0  0 0 0 0 0;     % cd - prosocial - state 4
    0     0 0 0 1  0 0 0 0 0;     % dd - prosocial - state 5the

    0 0 0 0 0  0     0 0 0 0;     % cc - nonsocial
    0 0 0 0 0  cn(1) 1 0 0 0;     % cc - nonsocial
    0 0 0 0 0  0     0 1 0 0;     % dc - nonsocial
    0 0 0 0 0  cn(2) 0 0 1 0;     % cd - nonsocial
    0 0 0 0 0  0     0 0 0 1];    % dd - nonsocial

B{2} = ...                        % defect:
   [0     0 0 0 0  0 0 0 0 0;     % start - nonsocial trustee
    0     1 0 0 0  0 0 0 0 0;     % ...
    dp(1) 0 1 0 0  0 0 0 0 0;
    0     0 0 1 0  0 0 0 0 0;
    dp(2) 0 0 0 1  0 0 0 0 0;

    0 0 0 0 0  0     0 0 0 0;
    0 0 0 0 0  0     1 0 0 0;
    0 0 0 0 0  dn(1) 0 1 0 0;
    0 0 0 0 0  0     0 0 1 0;
    0 0 0 0 0  dn(2) 0 0 0 1];

% observation probabilities
%--------------------------------------------------------------------------
A    = kron([1 1],speye(5,5));

% allowable policies - (of depth T); here, simply defect will cooperate
%--------------------------------------------------------------------------
V    = [1 2;
        1 1];

% MDP Structure
%==========================================================================
MDP.N = 8;                          % number of variational iterations
MDP.T = 2;                          % process depth (one-shot game)
MDP.S = S;                          % true initial state
MDP.A = A;                          % observation model
MDP.B = B;                          % transition probabilities (priors)
MDP.C = C;                          % terminal cost probabilities (priors)
MDP.D = D;                          % initial state probabilities (priors)
MDP.V = V;                          % allowable policies

MDP.alpha = 2;                      % gamma hyperparameters
MDP.beta  = 1/2;

% Solve - an example game
%==========================================================================
spm_figure('GetWin','Figure 1'); clf
MDP.plot = gcf;
MDP      = spm_MDP_game(MDP);

% now iterate repeated games accumulating posterior beliefs
%==========================================================================
MDP.plot = 0;
NG       = 64;

for i = 1:NG

    % solve and marginalise over posterior beliefs about hidden states
    %----------------------------------------------------------------------
    MDP    = spm_MDP_game(MDP);
    Q(:,i) = spm_softmax(k);
    O(:,i) = MDP.O(:,end);
    W(:,i) = MDP.W(:,end);
    P(:,i) = MDP.P(:,1);

    % update prior beliefs about initial state (pro-social or
    % non-social) and associated utility functions
    %----------------------------------------------------------------------
    a      = find(MDP.U(:,1));
    p      = MDP.O(:,2)'*MDP.A*MDP.B{a}(:,[1 6]);
    p      = p(:)/sum(p);
    k      = k + log(p);
    p      = spm_softmax(k);
    MDP.D  = kron(p,[1 0 0 0 0]');
    MDP.C  = [pp*p(1); pn*p(2)];

end

% graphics
%==========================================================================
spm_figure('GetWin','Figure 2'); clf

% posterior beliefs about hidden states (prosocial versus nonsocial)
%--------------------------------------------------------------------------
subplot(2,2,1)
plot(1:NG,Q)
title('Beliefs about other','FontSize',16)
xlabel('Number of games','FontSize',12)
ylabel('True and posterior expectations','FontSize',12)
spm_axis tight, axis square
legend({'prosocial','nonsocial'})

subplot(2,2,2)
plot(1:NG,P)
title('Beliefs about control','FontSize',16)
xlabel('Number of games','FontSize',12)
ylabel('True and posterior expectations','FontSize',12)
spm_axis tight, axis square
legend({'cooperate','defect'})

subplot(2,2,3)
imagesc(O)
title('Outcomes','FontSize',16)
xlabel('Number of games','FontSize',12)
ylabel('Observed outcome','FontSize',12)
axis square

subplot(2,2,4)
plot(W)
title('Precision','FontSize',16)
xlabel('Number of games','FontSize',12)
ylabel('Expected precision','FontSize',12)
axis square

author: spm


text:
classdef IMMOEA_F10 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing IM-MOEA

%------------------------------- Reference --------------------------------
% R. Cheng, Y. Jin, K. Narukawa, and B. Sendhoff, A multiobjective
% evolutionary algorithm using Gaussian process-based inverse modeling,
% IEEE Transactions on Evolutionary Computation, 2015, 19(6): 838-856.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = [1,zeros(1,obj.D-1)+10];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            t = X(:,2:obj.D).^(1./(1+3*repmat(2:obj.D,size(X,1),1)/obj.D)) - repmat(X(:,1),1,obj.D-1);
            g = 1 + 10*(obj.D-1) + sum(t.^2-10*cos(2*pi*t),2);
            PopObj(:,1) = X(:,1);
            PopObj(:,2) = g.*(1-sqrt(PopObj(:,1)./g));
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - sqrt(R(:,1));
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
classdef PostProcessingConfiguratorUIElements < handle
    properties (Access = public)
        algorithmsList;
        algorithmsPropertiesTable;
        annotationsPanel;
        annotationMappingList;
        annotationMappingTable;
    end

    methods (Access = public)
        function obj = PostProcessingConfiguratorUIElements(algorithmsList,...
                algorithmsPropertiesTable, annotationsPanel, annotationMappingList,...
            annotationMappingTable)
            obj.algorithmsList = algorithmsList;
            obj.annotationsPanel = annotationsPanel;
            obj.algorithmsPropertiesTable = algorithmsPropertiesTable;
            obj.annotationMappingList = annotationMappingList;
            obj.annotationMappingTable = annotationMappingTable;
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function [y_splined, x_spline, s_weights, y_splined_ext] = splinerMat(x, y, dxs, reg_factor, x_ext)
% SYNTAX:
%   [y_splined, x_spline, s_weights, y_splined_ext] = splinerMat(x, y, dxs, reg_factor, <x_ext>)
%
% EXAMPLE:
%   [y_splined, x_spline, s_weights, y_splined_ext] = splinerMat(x, y, 4, 0, x_ext);
%   [y_splined, x_spline, s_weights, y_splined_ext] = splinerMat(x, [y y_var], 4, 0, x_ext);
%
% INPUT:
%   x            [n x 1] observation time
%   y            [n x 1] observation value
%                [n x 2] observation value and variances (in this case the variances of the data are taken into account)
%   dxs          [1 x 1] spline base size
%   reg_factor   [1 x 1] regularization factor on the first derivative
%   x_ext        [m x 1] points in which to compute the interpolation
%
% OUTPUT:
%   y_splined     [n x 1] interpolated observation (on the observation epochs)
%   x_spline      [o x 1] center of the (o) splines
%   s_weights     [o x 1] weights of the splines
%   y_splined_ext [m x 1] spline interpolated in x_ext positions
%
% DESCRIPTION:
%   Interpolate with cubic splines a given dataset
%

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:        Andrea Gatti
%  Contributors:      Andrea Gatti, Giulio Tagliaferro
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

    if (nargin < 4)
        reg_factor = 0;
    end
    if isempty(x)
        x = 1:size(y, 1);
    end

    inan = isnan(y);
    if (size(y,2) == 2)
        inan = inan(:,1) | inan(:,2);
    end

    x(inan) = [];
    y(inan) = [];

    [x, id] = sort(x);
    y = y(id,:);
    if ~isempty(y)
        if ((nargin == 3) || (reg_factor == 0))
            if (size(y,2) == 2)
                [y_splined, x_spline, s_weights] = spliner_v51(x,y(:,1),y(:,2),dxs);
            else
                [y_splined, x_spline, s_weights] = spliner_v5(x,y,dxs);
            end
        else
            if (size(y,2) == 2)
                [y_splined, x_spline, s_weights] = spliner_v51R(x,y(:,1),y(:,2),dxs, reg_factor);
            else
                [y_splined, x_spline, s_weights] = spliner_v5R(x,y,dxs, reg_factor);
            end
        end
    else
        y_splined = y;
        x_spline = [];
        s_weights = [];
    end

    % Interpolation => using spline to predict in different coordinates
    % (not present in the C version)
    if (nargin == 5)
        if ~isempty(x_spline)
            mask = (isnan(s_weights));
            if (length(mask) > 2)
                mask = mask | [mask(2:end); 0] | [0; mask(1:end-1)];
            end
            s_weights = interp1(x_spline(~mask),s_weights(~mask),x_spline);
            if (size(x_ext,1)==1)
                x_ext = x_ext';
            end
            y_splined_ext = zeros(size(x_ext,1),1);
            for s = 1:length(x_spline)
                tau = round((x_ext-repmat(x_spline(s),length(x_ext),1))/dxs *1e13)/1e13; % 1e13 rounding necessary to avoid numerical problems
                y_splined_ext = y_splined_ext + s_weights(s)*cubicSpline(tau);
            end
        else
            y_splined_ext = nan(numel(x_ext),1);
        end
    else
        y_splined_ext = y_splined;
    end
    tmp = nan(numel(inan),1);
    tmp(~inan) = y_splined;
    y_splined = tmp;
end

% No Regularization + variances
function [y_splined, x_spline, s_weights] = spliner_v51(x, y, y_var, dxs)
    nObs = length(x);

    % size of the intervall to interpolate
    x_span = x(nObs) - x(1);
    x0 = x(1);
    x = x - x0;

    % compute the number of splines needed for the interpolation
    n_splines = ceil(x_span/dxs) + 3;

    % compute spline centers
    x_spline = zeros(n_splines,1);
    s_weights = [];
    s_center = x(1) - (((n_splines-3)*dxs-x_span)/2) - dxs;
    for i = 1:n_splines
        x_spline(i) = s_center+(i-1)*dxs;
    end

    % init A matrix
    A = zeros(nObs, 4);
    skips = zeros(n_splines,1);
    N = sparse(n_splines,n_splines);
    TN = zeros(n_splines, 1);

    cur_spline = 1;          % first spline whose domain intersect the observation
    %tau = 0;                % normalized distance between the observation and the center of the cur_spline
    i = 1;                  % index of the first observation
    first_obs = i;          % first observation used in the current A matrix
    first_spline = 1;       % first spline used in the current A matrix
    n_skip = 0;              % number of spline to "skip" because ain't intersecting an observation
    used_obs = 0;            % number of observation used in building the N matrix
    y_splined = zeros(length(y),1);           % output
    skips(1) = 0;
    while (i <= nObs)
        % Compute the distance between the current observation and the current spline
        tau = round((x(i)-x_spline(cur_spline))/dxs *1e13)/1e13; % 1e13 rounding necessary to avoid numerical problems
        if (tau <= 2)
            % fill the design matrix
            A(i-first_obs+1,:) = cubicSpline([tau (tau-1) (tau-2) (tau-3)]);
            used_obs = used_obs+1;
            n_skip = 0;
            i = i+1;
        else
            cur_local_spline = cur_spline-first_spline+1;
            skips(cur_spline+1) = i-first_obs;
            % This block of the A matrix is completed
            % Computing N
            n_skip = n_skip +1;
            if (n_skip < 4)
                if (n_skip == 1)
                    A2 = A((skips(cur_spline)+1):i-first_obs,:);
                    iQ = sparse(diag(1./y_var(first_obs+skips(cur_spline):i-1)));
                    N(cur_local_spline:cur_local_spline+3,cur_local_spline:cur_local_spline+3) = sparse(N(cur_local_spline:cur_local_spline+3,cur_local_spline:cur_local_spline+3) + A2'*iQ*A2);

                    % Computing TN
                    TN(cur_local_spline:cur_local_spline+3) = TN(cur_local_spline:cur_local_spline+3) + A2' * iQ * y(first_obs+skips(cur_spline):i-1);
                end
                cur_spline = cur_spline +1;
            else
                % If I skip more than 3 times the spline solutions are independent,
                % I can start solving my filtering for the first i points

                s_par = [];
                if (used_obs < size(N,2))
                    fprintf('WARNING: Regularization is needed observations are less than splines.\n         Adding 1e-9 on the normal matrix diagonal\n');
                    R = sparse(eye(cur_local_spline)*1e-9);
                    s_par = (N(1:cur_local_spline,1:cur_local_spline)+R)\TN(1:cur_local_spline);
                else
                    s_par = (N(1:cur_local_spline,1:cur_local_spline))\TN(1:cur_local_spline);
                end
                s_weights = [s_weights; s_par];

                for s = first_spline:cur_spline-3
                    y_splined(skips(s)+first_obs:skips(s+1)+first_obs-1) = y_splined(skips(s)+first_obs:skips(s+1)+first_obs-1) + A((skips(s):skips(s+1)-1)+1,:) * s_par(s-first_spline+1:s+3-first_spline+1);
                end

                % find the next spline whose domain intersect the next observation
                tau = (x(i)-x_spline(cur_spline))/dxs;
                while (tau > 2)
                    cur_spline = cur_spline+1;
                    tau = (x(i)-x_spline(cur_spline))/dxs;
                    if (tau > 2),
                        s_weights(cur_spline) = nan;
                    end
                end
                first_spline = cur_spline;
                first_obs = i;
                skips(cur_spline) = 0;

                A = zeros(nObs-i, 4);
                N = sparse(n_splines-cur_spline+1, n_splines-cur_spline+1);
                TN = zeros(n_splines-cur_spline+1,1);
                used_obs = 0;
            end
        end
    end

    skips(cur_spline+1) = i-first_obs;
    if (n_skip == 0)
        A2 = A(skips(cur_spline)+1:i-first_obs,:);
        iQ = sparse(diag(1./y_var(first_obs+skips(cur_spline):i-1)));
        N(cur_spline-first_spline+1:cur_spline+3-first_spline+1,cur_spline-first_spline+1:cur_spline+3-first_spline+1) = sparse(N(cur_spline-first_spline+1:cur_spline+3-first_spline+1,cur_spline-first_spline+1:cur_spline+3-first_spline+1) + A2'*iQ*A2);

        % Computing TN
        TN(cur_spline-first_spline+1:cur_spline+3-first_spline+1) = TN(cur_spline-first_spline+1:cur_spline+3-first_spline+1) + A2' * iQ * y(first_obs+skips(cur_spline):i-1);
    end

    % find the interpolation for the last subset of observations
    s_par = [];
    if (used_obs < size(N,2))
        fprintf('WARNING: Regularization is needed observations are less than splines.\n         Adding 1e-9 on the normal matrix diagonal\n');
        R = speye(size(N,2), size(N,2))*1e-9;
        s_par = (N+R)\TN;
    else
        s_par = N\TN;
    end
    s_weights = [s_weights; s_par];

    for s = first_spline:cur_spline
        if ((skips(s)<skips(s+1)))
            y_splined(skips(s)+first_obs:skips(s+1)+first_obs-1) = y_splined(skips(s)+first_obs:skips(s+1)+first_obs-1) + A((skips(s):skips(s+1)-1)+1,:) * s_par(s-first_spline+1:s+3-first_spline+1);
        end
    end

    x_spline = x_spline + x0;
end

% Regularization + variances
function [y_splined, x_spline, s_weights] = spliner_v51R(x, y, y_var, dxs, reg_factor)
    nObs = length(x);

    y_var(y_var < reg_factor/2) = reg_factor/2;
    % size of the intervall to interpolate
    x_span = x(nObs) - x(1);
    x0 = x(1);
    x = x - x0;

    % compute the number of splines needed for the interpolation
    n_splines = ceil(x_span/dxs) + 3;

    % compute spline centers
    x_spline = zeros(n_splines,1);
    s_weights = [];
    s_center = x(1) - (((n_splines-3)*dxs-x_span)/2) - dxs;
    for i = 1:n_splines
        x_spline(i) = s_center+(i-1)*dxs;
    end

    % init A matrix
    A = zeros(nObs, 4);
    skips = zeros(n_splines,1);
    N = sparse(n_splines,n_splines);
    TN = zeros(n_splines, 1);

    cur_spline = 1;          % first spline whose domain intersect the observation
    tau = 0;                % normalized distance between the observation and the center of the cur_spline
    i = 1;                  % index of the first observation
    first_obs = i;          % first observation used in the current A matrix
    first_spline = 1;       % first spline used in the current A matrix
    n_skip = 0;              % number of spline to "skip" because ain't intersecting an observation
    used_obs = 0;            % number of observation used in building the N matrix
    y_splined = zeros(length(y),1);           % output
    skips(1) = 0;
    while (i <= nObs)
        % Compute the distance between the current observation and the current spline
        tau = round((x(i)-x_spline(cur_spline))/dxs *1e13)/1e13; % 1e13 rounding necessary to avoid numerical problems
        if (tau <= 2)
            % fill the design matrix
            A(i-first_obs+1,:) = cubicSpline([tau (tau-1) (tau-2) (tau-3)]);
            used_obs = used_obs+1;
            n_skip = 0;
            i = i+1;
        else
            cur_local_spline = cur_spline-first_spline+1;
            skips(cur_spline+1) = i-first_obs;
            % This block of the A matrix is completed
            % Computing N
            n_skip = n_skip +1;
            if (n_skip < 4)
                if (n_skip == 1)
                    A2 = A((skips(cur_spline)+1):i-first_obs,:);
                    iQ = sparse(diag(1./y_var(first_obs + skips(cur_spline):i-1)));
                    N(cur_local_spline:cur_local_spline+3,cur_local_spline:cur_local_spline+3) =   sparse(N(cur_local_spline:cur_local_spline+3,cur_local_spline:cur_local_spline+3) + A2'*iQ*A2);

                    % Computing TN
                    TN(cur_local_spline:cur_local_spline+3) = TN(cur_local_spline:cur_local_spline+3) + A2' * iQ * y(first_obs+skips(cur_spline):i-1);
                end
                cur_spline = cur_spline +1;
            else
                % find the next spline whose domain intersect the next observation
                tau = (x(i)-x_spline(cur_spline))/dxs;
                while (tau > 2)
                    cur_spline = cur_spline+1;
                    tau = (x(i)-x_spline(cur_spline))/dxs;
                end
                skips(cur_spline) = skips(cur_spline-1);
                used_obs = -1e10;
            end
        end
    end

    skips(cur_spline+1) = i-first_obs;
    if (n_skip == 0)
        A2 = A(skips(cur_spline)+1:i-first_obs,:);
        iQ = sparse(diag(1./y_var(first_obs+skips(cur_spline):i-1)));
        N(cur_spline-first_spline+1:cur_spline+3-first_spline+1,cur_spline-first_spline+1:cur_spline+3-first_spline+1) = sparse(N(cur_spline-first_spline+1:cur_spline+3-first_spline+1,cur_spline-first_spline+1:cur_spline+3-first_spline+1) + A2'*iQ*A2);

        % Computing TN
        TN(cur_spline-first_spline+1:cur_spline+3-first_spline+1) = TN(cur_spline-first_spline+1:cur_spline+3-first_spline+1) + A2' * iQ * y(first_obs+skips(cur_spline):i-1);
    end

    % find the interpolation for the last subset of observations
    s_par = [];
    if (size(N,2)>2)
        if reg_factor == 0
            fprintf('WARNING: Regularization is needed observations are less than splines.\n         Adding 1e-9 on the normal matrix diagonal\n');
            reg_factor = 1e-9;
        end
        R = sparse(eye(size(N,2))-diag(ones(size(N,2)-1,1),1)-diag(ones(size(N,2)-1,1),-1) + diag([0; ones(size(N,2)-2,1); 0]));
        R = R*reg_factor;
        s_par = (N+R)\TN;
    else
        if (used_obs < size(N,2))
            if reg_factor == 0
                fprintf('WARNING: Regularization is needed observations are less than splines.\n         Adding 1e-9 on the normal matrix diagonal\n');
                reg_factor = 1e-9;
            end
            R = sparse(eye(size(N,2))*reg_factor);
            s_par = (N+R)\TN;
        else
            s_par = N\TN;
        end
    end
    s_weights = [s_weights; s_par];

    for s = first_spline:cur_spline
        if ((skips(s)<skips(s+1)))
            y_splined(skips(s)+first_obs:skips(s+1)+first_obs-1) = y_splined(skips(s)+first_obs:skips(s+1)+first_obs-1) + A((skips(s):skips(s+1)-1)+1,:) * s_par(s-first_spline+1:s+3-first_spline+1);
        end
    end

    x_spline = x_spline + x0;
end

% No Regularization - no variances
function [y_splined, x_spline, s_weights] = spliner_v5(x, y, dxs)
    nObs = length(x);

    % size of the intervall to interpolate
    x_span = x(nObs) - x(1);
    x0 = x(1);
    x = x - x0;

    % compute the number of splines needed for the interpolation
    n_splines = ceil(x_span/dxs) + 3;

    % compute spline centers
    x_spline = zeros(n_splines,1);
    s_center = x(1) - (((n_splines-3)*dxs-x_span)/2) - dxs;
    for i = 1:n_splines
        x_spline(i) = s_center+(i-1)*dxs;
    end

    % init A matrix
    A = zeros(nObs, 4);
    A_idx = zeros(nObs, 4);

    tau   = round(rem(x',dxs)/dxs*1e13)/1e13;  % 1e13 rounding necessary to avoid numerical problems
    idx   = floor((x')/dxs)+1;
    A     = cubicSpline4Col(tau);
    A_idx = [idx(:) idx(:)+1 idx(:)+2 idx(:)+3];
    n_par = max(A_idx(:,4));
    n_obs = numel(x);
    rows  = repmat((1:n_obs)',1,4);

    A = sparse(rows, A_idx, A, n_obs, n_par);

    idx_null = sum(A~=0) == 0;
    A(:,idx_null) = [];

    N = A'*A;
    B = A'*y;

    x = N\B;

    s_weights = nan(numel(idx_null),1);
    s_weights(~idx_null) = x;
    y_splined = A*x;

    x_spline = x_spline + x0;
end

% Regularization - no variances
function [y_splined, x_spline, s_weights] = spliner_v5R(x, y, dxs, reg_factor)
    nObs = length(x);

    % size of the intervall to interpolate
    x_span = x(nObs) - x(1);
    x0 = x(1);
    x = x - x0;

    % compute the number of splines needed for the interpolation
    n_splines = ceil((x_span+eps(x_span))/dxs) + 3;

    % compute spline centers
    x_spline = zeros(n_splines,1);
    s_weights = [];
    s_center = x(1) - (((n_splines-3)*dxs-x_span)/2) - dxs;
    for i = 1:n_splines
        x_spline(i) = s_center+(i-1)*dxs;
    end

    % init A matrix
    A = zeros(nObs, 4);
    skips = zeros(n_splines,1);
    N = sparse(n_splines,n_splines);
    TN = zeros(n_splines, 1);

    cur_spline = 1;          % first spline whose domain intersect the observation
    tau = 0;                % normalized distance between the observation and the center of the cur_spline
    i = 1;                  % index of the first observation
    first_obs = i;          % first observation used in the current A matrix
    first_spline = 1;       % first spline used in the current A matrix
    n_skip = 0;              % number of spline to "skip" because ain't intersecting an observation
    used_obs = 0;            % number of observation used in building the N matrix
    y_splined = zeros(length(y),1);           % output
    skips(1) = 0;
    while (i <= nObs)
        % Compute the distance between the current observation and the current spline
        tau = round((x(i) - x_spline(cur_spline))/dxs * 1e13) / 1e13; % 1e13 rounding necessary to avoid numerical problems
        if (tau <= 2)
            % fill the design matrix
            A(i-first_obs+1,:) = cubicSpline([tau (tau-1) (tau-2) (tau-3)]);
            used_obs = used_obs+1;
            n_skip = 0;
            i = i+1;
        else
            cur_local_spline = cur_spline-first_spline+1;
            skips(cur_spline+1) = i-first_obs;
            % This block of the A matrix is completed
            % Computing N
            n_skip = n_skip +1;
            if (n_skip < 4)
                if (n_skip == 1)
                    A2 = A((skips(cur_spline)+1):i-first_obs,:);
                    N(cur_local_spline:cur_local_spline+3,cur_local_spline:cur_local_spline+3) =   sparse(N(cur_local_spline:cur_local_spline+3,cur_local_spline:cur_local_spline+3) + A2'*A2);

                    % Computing TN
                    TN(cur_local_spline:cur_local_spline+3) = TN(cur_local_spline:cur_local_spline+3) + A2' * y(first_obs+skips(cur_spline):i-1);
                end
                cur_spline = cur_spline +1;
            else
                % find the next spline whose domain intersect the next observation
                tau = (x(i)-x_spline(cur_spline))/dxs;
                while (tau > 2)
                    cur_spline = cur_spline+1;
                    tau = (x(i)-x_spline(cur_spline))/dxs;
                end
                skips(cur_spline) = skips(cur_spline-1);
                used_obs = -1e10;
            end
        end
    end

    skips(cur_spline+1) = i-first_obs;
    if (n_skip == 0)
        A2 = A(skips(cur_spline)+1:i-first_obs,:);
        N(cur_spline-first_spline+1:cur_spline+3-first_spline+1,cur_spline-first_spline+1:cur_spline+3-first_spline+1) = N(cur_spline-first_spline+1:cur_spline+3-first_spline+1,cur_spline-first_spline+1:cur_spline+3-first_spline+1) + A2'*A2;

        % Computing TN
        TN(cur_spline-first_spline+1:cur_spline+3-first_spline+1) = TN(cur_spline-first_spline+1:cur_spline+3-first_spline+1) + A2' * y(first_obs+skips(cur_spline):i-1);
    end

    % find the interpolation for the last subset of observations
    s_par = [];
    if (size(N,2)>2)
        %fprintf('WARNING: Regularization is needed observations are less than splines.\n         Adding 1e-9 on the normal matrix diagonal\n');
        R = (speye(size(N,2), size(N,2)) - spdiags(ones(size(N,2), 1), 1, size(N,2), size(N,2)) - spdiags(ones(size(N,2), 1), -1, size(N,2), size(N,2)) + spdiags([0; ones(size(N,2) - 2, 1); 0], 0, size(N,2), size(N,2))) * reg_factor;
        s_par = (N+R)\TN;
    else
        if (used_obs < size(N,2))
            %fprintf('WARNING: Regularization is needed observations are less than splines.\n         Adding 1e-9 on the normal matrix diagonal\n');
            R = speye(size(N,2)) * reg_factor;
            s_par = (N+R)\TN;
        else
            s_par = N\TN;
        end
    end
    s_weights = [s_weights; s_par];

    for s = first_spline:cur_spline
        if (skips(s) == 0 && s > 1)
            skips(s) = skips(s-1);
        end
        if ((skips(s)<skips(s+1)))
            y_splined(skips(s)+first_obs:skips(s+1)+first_obs-1) = y_splined(skips(s)+first_obs:skips(s+1)+first_obs-1) + A((skips(s):skips(s+1)-1)+1,:) * s_par(s-first_spline+1:s+3-first_spline+1);
        end
    end

    x_spline = x_spline + x0;
end

% SYNTAX:
%   [y] = cubicSpline(t)
%
% EXAMPLE:
%   [y] = cubicSpline(1)
%
% INPUT:
%   t = normalized value from [-2 to 2] of the distance between the
%       observation and the center of the cubic spline
%
% OUTPUT:
%   y = value of the spline in the given point
%
% DESCRIPTION:
%   Get the value of the cubic spline with a base of 4 at the given t normalized
%   value
%
% USEFUL VALUES:
%   2/3 = cubicSpline(0);
%   1/6 = cubicSpline(-1);
%   1/6 = cubicSpline(1);
%
% by Andrea Gatti
%
function [y] = cubicSpline(t)

% FAST VECTORIAL IMPLEMENTATION
y = zeros(size(t));
pos = (t > -2) + (t > -1) - (t > 1) - (t > 2);

% pos = 0    => spline = 0;
% pos = 1    => spline = (-abs(t)+2)^3/6;
% pos = 2    => spline = ((-abs(t)+2)^3 - 4*(-abs(t)+1)^3)/6;

p1 = pos==1;
p2 = pos==2;
t = abs(t);
y(p1) = (2-t(p1)).^3/6;
y(p2) = ((2-t(p2)).^3 - 4*(1-t(p2)).^3)/6;

%  PLAIN IMPLEMENTATION
%
% 	if ((t < -2) || (t > 2))
% 		y = 0;
%     else
%         if (t < 0)
%             if (t <= -1)
%                 tmp2=(t+2);
%                 y = (tmp2*tmp2*tmp2 / 6);
%             else
%                 tmp1=(t+1);
%                 tmp2=(t+2);
%                 y = ((tmp2*tmp2*tmp2 - 4*tmp1*tmp1*tmp1) / 6);
%             end
%         else
%             if (t < 1)
%                 tmp1=(1-t);
%                 tmp2=(2-t);
%                 y = ((tmp2*tmp2*tmp2 - 4*tmp1*tmp1*tmp1) / 6);
%             else
%                 tmp2=(2-t);
%                 y = (tmp2*tmp2*tmp2 / 6);
%             end
%         end
%     end
end

function [val] = cubicSpline4Col(t)
    % Compute matrix entry for cubic spline
    %
    % INPUT
    %   t -> 0 : 1
    %   order -> 1,3
    %
    % SYNTAX:
    %  Core_Utils.cubicSplic(t)
    val = zeros(numel(t),4);
    val(:,1) = (1 - t).^3/6;
    val(:,2) = ((2-t).^3 - 4*(1-t).^3)/6;
    val(:,3) = ((1+t).^3 - 4*(t).^3)/6;
    val(:,4) = (t).^3/6;
end

author: goGPS-Project


text:
function rescaleROIs(sourceRes, targetRes, sourceDir, targetDir)
% rescaleROIs(sourceRes, targetRes, [sourceDir], [targetDir])
%
% Purpose: Convert ROI from one spatial resolution to another.
%
% written by JW 04.02.08
%
%   inputs:
%
%       sourceRes: 3d resolution of ROIs in source dir, in mm (e.g, [1 1 1])
%       targetDir: 3d resolution of new ROIs, in mm (e.g, [0.7 0.7 0.7])
%       sourceDir: a source directory containing ROIs; default = pwd
%       targetDir: a destination directory; default = sourceDir

if ieNotDefined('sourceDir'), sourceDir = pwd; end
if ieNotDefined('targetDir'), targetDir = sourceDir; end

scaleFactor = diag(sourceRes/diag(targetRes));

%% find all .mat files in the source directory
currentDir = pwd;
cd(sourceDir);

w = dir('*.mat');
fileList = {w.name};

for i = 1:length(fileList)
    load (fileList{i})
    ROI.coords = single(round(scaleFactor * ROI.coords));
    savePath = fullfile(targetDir, [ROI.name '_resampled.mat']);
    save (savePath, 'ROI')
end

cd(currentDir);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function value = i4_bset ( i4, pos )

%*****************************************************************************80
%
%% I4_BSET returns a copy of an I4 in which the POS-th bit is set to 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Military Standard 1753,
%    FORTRAN, DoD Supplement To American National Standard X3.9-1978,
%    9 November 1978.
%
%  Parameters:
%
%    Input, integer I4, the integer to be tested.
%
%    Input, integer POS, the bit position, between 0 and 31.
%
%    Output, integer VALUE, a copy of I4, but with the POS-th bit
%    set to 1.
%
  i4_huge = 2147483647;

  value = round ( i4 );

  if ( pos < 0 )

  elseif ( pos < 31 )

    add = 1;

    if ( 0 <= i4 )
      j = i4;
    else
      j = ( i4_huge + i4 ) + 1;
    end

    for k = 1 : pos
      j = floor ( j / 2 );
      add = add * 2;
    end

    if ( mod ( j, 2 ) == 0 )
      value = i4 + add;
    end

  elseif ( pos == 31 )

    if ( 0 < i4 )
      value = - ( i4_huge - i4 ) - 1;
    end

  elseif ( 31 < pos )

    value = i4;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function dist = getNormGTJointDistMulti(keypointsAll,gt_annolist,jidx,bUseHeadSize,bHeadSizeFromRect,jidxGT,nrects)

if (nargin < 5)
    bHeadSizeFromRect = true;
end

if (nargin < 6)
    jidxGT = jidx;
end

scCorr = 1.33;% WAF scale correction factor when using head annopoints

dist = nan(nrects,1);
n = 0;
for imgidx = 1:length(gt_annolist)

    rect = gt_annolist(imgidx).annorect;
    for ridx = 1:length(rect)
        n = n + 1;
        if (bUseHeadSize)
            if (bHeadSizeFromRect)
                refDist = util_get_head_size(rect(ridx));
            else
                if (isfield(rect(ridx), 'annopoints'))
                    points = rect(ridx).annopoints.point;
                    p1 = util_get_annopoint_by_id(points,8);
                    p2 = util_get_annopoint_by_id(points,9);
                    refDist = scCorr*norm([p1.x p1.y] - [p2.x, p2.y]);
                else
                    assert(false);
                end
            end
        else
            refDist = util_get_torso_size(rect(ridx));
            if (isnan(refDist))
                continue;
            end
        end

        if (isfield(rect(ridx), 'annopoints'))
            points = rect(ridx).annopoints.point;
            pointsAll = [];
            for ppidx = 1:length(points)
                pointsAll = [pointsAll; points(ppidx).x points(ppidx).y];
            end
            bbox = [min(pointsAll(:,1)) min(pointsAll(:,2)) max(pointsAll(:,1)) max(pointsAll(:,2))];
            d = refDist*0.1;
            bbox = bbox + [-d -d d d];
            p = util_get_annopoint_by_id(rect(ridx).annopoints.point, jidxGT);
            if (~isempty(p))
                det = keypointsAll(imgidx).det{jidx+1};
                idxs = det(:,1) >= bbox(1) & det(:,2) >= bbox(2) & det(:,1) <= bbox(3) & det(:,2) <= bbox(4);
                detRect = det(idxs,:);
                if (~isempty(detRect))
                    [val,id] = max(detRect(:,3));
                    points_det = detRect(id,1:2);
                else
                    points_det = [inf inf];
                end
                dist(n) = norm([p.x p.y] - points_det)/refDist;
                bVis = false;
                if (bVis)
                    figure(100); clf; imagesc(imread(gt_annolist(imgidx).image.name)); axis equal; hold on;
                    vis_bbox(bbox,'g');
                    plot(points_det(:,1),points_det(:,2),'ro','MarkerSize',6,'MarkerFaceColor','r');
                    plot(p.x,p.y,'go','MarkerSize',6,'MarkerFaceColor','g');
                    text(bbox(1),bbox(2),sprintf('%1.2f',dist(n)),'FontSize',10,'color','k','BackgroundColor','g',...
                        'verticalalignment','top','horizontalalignment','left');
                end
            end
        end
    end
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [x] = boundProject(x,LB,UB)
% function [x] = boundProject(x,LB,UB)
%   Computes projection of x onto constraints LB <= x <= UB

x(x < LB) = LB(x < LB);
if nargin > 2
x(x > UB) = UB(x > UB);
end

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function tests = test_savestack(varargin)
%TEST_SAVESTACK  Test case for the file savestack
%
%   Test case for the file savestack

%   Example
%   test_savestack
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-09-27,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

tests = functiontests(localfunctions);

function test_SaveGrayStack(testCase)
% Test call of function without argument

D = load ('mri');
D = squeeze(D.D);

fName = 'testSaveGrayStack.tif';
if exist(fName, 'file')
    delete(fName);
end

savestack(D, fName);

info = imfinfo(fName);
assertEqual(testCase, 27, length(info));

delete(fName);

function test_SaveGraySlices(testCase)

D = load ('mri');
D = squeeze(D.D);

genericName = 'testSaveGraySlices*.tif';
patternName = 'testSaveGraySlices###.tif';

list = dir(genericName);
for i = 1:length(list)
    delete(list(i).name);
end

savestack(D, patternName);

list = dir(genericName);
assertEqual(testCase, 27, length(list));

list = dir(genericName);
for i = 1:length(list)
    delete(list(i).name);
end

function test_SaveRGBStack(testCase)
% Test call of function without argument

% create a 3D color image
mri = load ('mri');
mri = squeeze(mri.D);
bin = imclose(mri>0, ones(3, 3));
bnd = imdilate(imBoundary(bin), ones(3, 3));
ovr = imOverlay(mri*3, bnd);

fName = 'testSaveRGBStack.tif';
if exist(fName, 'file')
    delete(fName);
end

savestack(ovr, fName);

info = imfinfo(fName);
assertEqual(testCase, 27, length(info));

delete(fName);

function test_SaveRGBSlices(testCase)
% Test call of function without argument

% create a 3D color image
mri = load ('mri');
mri = squeeze(mri.D);
bin = imclose(mri > 0, ones(3, 3));
bnd = imdilate(imBoundary(bin), ones(3, 3));
ovr = imOverlay(mri*3, bnd);

genericName = 'testSaveRGBSlices*.tif';
patternName = 'testSaveRGBSlices###.tif';

list = dir(genericName);
for i = 1:length(list)
    delete(list(i).name);
end

savestack(ovr, patternName);

list = dir(genericName);
assertEqual(testCase, 27, length(list));

list = dir(genericName);
for i = 1:length(list)
    delete(list(i).name);
end

function test_SaveGrayStackMap(testCase)

D = load ('mri');
D = squeeze(D.D);

fName = 'testSaveGrayStackMap.tif';
if exist(fName, 'file')
    delete(fName);
end
map = jet(256);

savestack(D, map, fName);

info = imfinfo(fName);
assertEqual(testCase, 27, length(info));
assertFalse(testCase, isempty(info(1).Colormap));

delete(fName);

function test_SaveGraySlicesMap(testCase)

D = load ('mri');
D = squeeze(D.D);
map = jet(256);

genericName = 'testSaveGraySlicesMap*.tif';
patternName = 'testSaveGraySlicesMap###.tif';

list = dir(genericName);
for i = 1:length(list)
    delete(list(i).name);
end

savestack(D, map, patternName);

list = dir(genericName);
assertEqual(testCase, 27, length(list));
info = imfinfo(list(1).name);
assertFalse(testCase, isempty(info.Colormap));

list = dir(genericName);
for i = 1:length(list)
    delete(list(i).name);
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
classdef YawAngleTermCondOptimVar < AbstractOptimizationVariable
    %YawAngleTermCondOptimVar Summary of this class goes here
    %   Detailed explanation goes here

    properties
        varObj(1,1) YawTermCondition = YawTermCondition(0);

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        useTf(1,1) = false;
    end

    methods
        function obj = YawAngleTermCondOptimVar(varObj)
            obj.varObj = varObj;
            obj.varObj.optVar = obj;

            obj.id = rand();
        end

        function x = getXsForVariable(obj)
            x = [];

            if(obj.useTf)
                x = obj.varObj.targetYawAngle;
            end
        end

        function [lb, ub] = getBndsForVariable(obj)
            lb = obj.lb(obj.useTf);
            ub = obj.ub(obj.useTf);
        end

        function [lb, ub] = getAllBndsForVariable(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function setBndsForVariable(obj, lb, ub)
            obj.lb = lb;
            obj.ub = ub;
        end

        function useTf = getUseTfForVariable(obj)
            useTf = obj.useTf;
        end

        function setUseTfForVariable(obj, useTf)
            obj.useTf = useTf;
        end

        function updateObjWithVarValue(obj, x)
            obj.varObj.targetYawAngle = x;
        end

        function nameStrs = getStrNamesOfVars(obj, evtNum, varLocType)
            nameStrs = {sprintf('Event %i Yaw Angle Termination Condition', evtNum)};
        end

        function varsStoredInRad = getVarsStoredInRad(obj)
            varsStoredInRad = true;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%
% Set the default parameters for a specific experiment (protocol/subject)
%
% Parameters
% ----------
% protocol : string
%   Folder for a specific acquisition protocol, possible containing different subjects
% subject : string
%   Folder of the specific subject
%
function AMICO_SetSubject( protocol, subject )

	global AMICO_data_path
	global CONFIG

	CONFIG = [];

	CONFIG.protocol			= protocol;
	CONFIG.subject			= subject;

	CONFIG.DATA_path		= fullfile( AMICO_data_path, CONFIG.protocol, CONFIG.subject );
	CONFIG.dwiFilename		= fullfile( CONFIG.DATA_path, 'DWI.nii' );
	CONFIG.maskFilename     = fullfile( CONFIG.DATA_path, 'roi_mask.nii' );
	CONFIG.dim              = [];
	CONFIG.pixdim           = [];
	CONFIG.OUTPUT_path      = [];

	CONFIG.schemeFilename	= fullfile( CONFIG.DATA_path, 'DWI.scheme' );
	CONFIG.scheme			= [];
    CONFIG.b0_thr           = 1;

	CONFIG.OPTIMIZATION     = [];

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% A real linear filter, usually lowpass, Abstract Method
%
% Copyright (c) 2018 Department of Computer Science,
%                    University of Toronto, Canada,
%                    Vector Institute, Canada
%
% License
% This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Yingxue Wang <yingxue@cs.toronto.edu>
%  Sean Robertson <sdrobert@cs.toronto.edu>
%

classdef WindowFunction < handle
    methods (Abstract)
        window = get_impulse_response(obj, width)
    end
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep