text:
function [imageCha, ssim_map, metrics] = algo_FCSA_WaTMRI_2D_real( obj,input )
% 2D variant of WaTMRI
% based on Huang et al. paper on WaTMRI
%
% input:
% obj           CS reconstruction object (holding all parameters)
% input         struct containing recon parameters and image
%
% output:
% imageCha      reconstructed channel individual image
% ssim_map      structural similarity map
% metrics       evaluation metrics
%
% (c) Marc Fischer, Thomas Kuestner, May 2015
% -------------------------------------------------------------------------

%%
timer_WaT = tic;

%% variables:
% internal flags
flag_wavetree = true;
flag_fast = true;
flag_extendImage = true;

% internal variables:
mue = 0.001;
L = 1;
t_old = 1;
NLTV_struct.kernelratio = 3;
NLTV_struct.windowratio = 6;
NLTV_struct.nThreads = 1;  % mind this option if used with -singleCompThread on BWCluster
itrNLTV = obj.iNINNER - 5;
% chambolle tv:
parsin.MAXITER=100; parsin.tv='iso'; % 'iso' or 'l1'

% from obj:
maxitr = obj.iNINNER;
% n1 = obj.measPara.dim(1);
% n2 = obj.measPara.dim(2);
n1 = input.n1;
n2 = input.n2;
% nSlices = obj.measPara.dim(3);
nCha = obj.measPara.dim(5);
lambdaWave = obj.lambda;
lambdaTV = obj.lambdaTV;
lambdaGroup = obj.lambdaGroup;
NLTV_struct.filterstrength = obj.lambdaNLTV_h; % 0.03 % converted from NLTV h = 0.01 %old: used: 3e-10
lambdaNLTV = obj.lambdaNLTV;
lambdaNLTV_h = obj.lambdaNLTV_h;
regularizerWeights = obj.regularizerWeights;
flagTV = obj.flagTV;
flagTV_iso = obj.flagTV_iso;
flagWave = obj.flagWave;
flagGroup = obj.flagGroup;
flagNLTV = obj.flagNLTV;
flagSBNLTV = obj.flagSBNLTV;
waveletStages = obj.trafo.waveletStages;
waveletFilterName_l1 = obj.trafo.waveletFilterName_l1;
waveletFilterName_l12 = obj.trafo.waveletFilterName_l12;

% from input:
b=input.b;
mask = input.mask;
G_wat = input.G_wat;
Gt_wat = input.Gt_wat;
% groupnorm_index = input.groupnorm_index;
waveS_l1 = input.waveS_l1;
waveS_l12 = input.waveS_l12;
waveS_WaT = input.waveS_WaT;
WaT_extend_y = waveS_WaT(waveletStages+2,1) - waveS_l12(waveletStages+2,1);
WaT_extend_x = waveS_WaT(waveletStages+2,2) - waveS_l12(waveletStages+2,2);
z_WaT = cell(1,nCha);

% im_ref = input.im_ref;
% im_ref_full = zeros(n1,n2);
% for j = 1:nCha
%     im_ref_full = im_ref_full + abs(im_ref{1,j}).^2;
% end;
% im_ref_full = sqrt(im_ref_full);
clear input

% initialize cells/vectors
for j=1:nCha
    FTb{1,j} = real(iFFT2D(b{1,j}));
    % y{1,j} = real(FTb{1,j});
    % y{1,j+nCha} = imag(FTb{1,j});
end;
z = FTb; % starting point
y = z;

x_wave = cell(1,nCha);
x_helper = x_wave;
x_wave_helper = x_wave;
x_tv = x_wave;
x_nltv = x_wave;
for j=1:nCha
    x_nltv{1,j} = 0;
end;
x_g = x_wave;
x_g_helper = x_wave;
x_g_WaT = x_wave;

%% MAD dependent lambdas:
flag_MAD = true;
if flag_MAD
    x_wavedec = cell(1,nCha);
    threshold = zeros(1:nCha);
    for j=1:nCha % 2*nCha
        x_wavedec{1,j} = wavedec2(z{1,j},waveletStages,waveletFilterName_l1); % atm only based on l1-daubechie
        x_wave_fine_scale = size(x_wavedec{1,j},2) - (3*waveS_l1(waveletStages+1,1)*waveS_l1(waveletStages+1,2));
        threshold(j) = mad(x_wavedec{1,j}(x_wave_fine_scale:end),1);
    end;
    clear x_wavedec
else
    threshold(1:nCha) = 1;
end;

threshold_wave(j) = lambdaWave * threshold(j) * 2/L;
threshold_TV(j) = lambdaTV * threshold(j) * 2/L;
threshold_group(j) = lambdaGroup * threshold(j) * 2/L;
threshold_NLTV(j) = lambdaNLTV; % * threshold(j) * 2/L;
threshold_NLTV_h(j) = lambdaNLTV_h; %*threshold(j); % adjust carefully or NLTV won't find a solution. lambdaNLTV_h should be < 0.01

for j = 1:nCha
    x_max(j) = max(max(z{1,1}))*1.0;
    x_min(j) = min(min(z{1,1}))*1.0;
end;

%% initialize metrics:
    itr = 0;
   	metrics.xtime(itr+1)= 0;
%     [metrics, ssim_map{1,1}] = get_metrics_itr( im_ref, im_ref_full, z, itr, maxitr, nCha, n1, n2, metrics, obj.K_1, obj.K_2, obj.W_size, obj.W_sigma );
    ssim_map = [];

%% recon
dispProgress('Proximal Average', 0, maxitr);
for itr = 1:maxitr  % total iter counter

    t_new = (1+sqrt(1+4*t_old^2))/2;
    t_old = t_new;

    y_old = z;   % y_old = y for complex case

%% WaTMRI-Step
    if itr > 1
        if flagGroup
            for j = 1:nCha
                if flag_extendImage
                    z_WaT{1,j} = extend_image(y{1,j}, waveS_WaT, waveletStages, WaT_extend_y, WaT_extend_x);
                    x_g_helper{1,j} = wavedec2(z_WaT{1,j},waveletStages,waveletFilterName_l12);
                else
                    x_g_helper{1,j} = wavedec2(z{1,j},waveletStages,waveletFilterName_l12);
                end;

                if flag_wavetree
                    x_g_helper{2,j} = (G_wat*x_g_helper{1,j}')';
                else
                    x_g_helper{2,j} = zeros(1,size(x_g_helper{1,j},2));
                end;
            end;

                x_g_helper = softthresh_group_2vec(x_g_helper,threshold_group(j),nCha); % threshold_wave(j)/threshold_group(j) acc. to author

            for j = 1:nCha
                x_g_WaT{1,j} = waverec2(x_g_helper{1,j}+(Gt_wat*x_g_helper{2,j}')',waveS_WaT,waveletFilterName_l12);
                x_g{1,j} =  x_g_WaT{1,j}(1:end-WaT_extend_y,1:end-WaT_extend_x);
            end;
        end;
    else
        for j = 1:nCha
            x_g{1,j} = zeros(n1,n2);
            z_WaT{1,j} = zeros(n1,n2);
        end;
    end;

%% landweber step + WaTMRI-part
    for j = 1:nCha
        x_helper{1,j} = real(iFFT2D(FFT2D_mask(z{1,j},mask))) -FTb{1,j} ;
        if itr > 1  % Phi'*G'*G*Phi*z{1,j}
            z_WaT{1,j} = extend_image(z{1,j}, waveS_WaT, waveletStages, WaT_extend_y, WaT_extend_x);
            z_WaT{1,j} = trafo_WaT(z_WaT{1,j},waveletStages,waveletFilterName_l12,G_wat,Gt_wat);
            z_WaT{1,j} = z_WaT{1,j}(1:end-WaT_extend_y,1:end-WaT_extend_x);
        end;
        z{1,j} = z{1,j} - (x_helper{1,j}  + mue*(z_WaT{1,j} - x_g{1,j}))/L; % mue used as tradeoff parameter.
    end;

%% l1-Wavelet
    if flagWave
            for j = 1:nCha % 2*nCha
                    x_wave_helper{1,j} = wavedec2(z{1,j},waveletStages,waveletFilterName_l1);
                    x_wave_helper{1,j} = softthresh_real(x_wave_helper{1,j},threshold_wave(j));
                    x_wave{1,j} = waverec2(x_wave_helper{1,j},waveS_l1,waveletFilterName_l1);
            end;
    end;

%% TV
    if flagTV
        if ~flagTV_iso
            for j = 1:nCha
                    x_tv{1,j} = MTV_2D(z{1,j},threshold_TV(j),n1,n2);
            end;
        else
            for j = 1:nCha
                if (itr==1)
                    [x_tv{1,j}, P]=denoise_TV_One((z{1,j}), threshold_TV(j),-inf,inf,[],parsin);
                else
                    [x_tv{1,j}, P]=denoise_TV_One((z{1,j}), threshold_TV(j),-inf,inf,P,parsin);
                end;
            end;
        end;
    end;

%% NLTV
    if flagNLTV
        if itr >= itrNLTV
            if flagSBNLTV
                for j = 1:nCha
                    x_nltv{1,j} = SB_NLTVfunc_slim_rescale(z{1,j},n1,n2, threshold_NLTV(j), threshold_NLTV_h(j) );
                end;
            else
                for j = 1:nCha
                    % if mod(itr,5) == 0 || itr == 1
                        x_nltv{1,j} = NLMF(z{1,j},NLTV_struct);
                        x_nltv{1,j} = (L.*z{1,j} + 2*threshold_NLTV(j)*x_nltv{1,j})./(L+2*threshold_NLTV(j));
                    % end;
                end;
            end;
        else
            x_nltv{1,j} = 0;
       end;
    end;

%% add prox(.)
    for j = 1:nCha
        y{1,j} = zeros(n1,n2);
        if flagWave y{1,j} = y{1,j} + x_wave{1,j}.*regularizerWeights(1); end;
        if flagTV y{1,j} = y{1,j} + x_tv{1,j}.*regularizerWeights(2); end;
        % if flagGroup y{1,j} = y{1,j} + x_g{1,j}.*regularizerWeights(3); end;
        if flagNLTV y{1,j} = y{1,j} + x_nltv{1,j}.*regularizerWeights(4); end;

        if ~flagWave && ~flagTV && ~flagNLTV
            y{1,j} = z{1,j};
            y{1,j}(y{1,j} > x_max(j)) = x_max;
            y{1,j}(y{1,j} < x_min(j)) = x_min;
        else
            if itr < itrNLTV
                y{1,j} = y{1,j}/(flagTV.*regularizerWeights(2) + flagWave.*regularizerWeights(1)); % flagGroup.*regularizerWeights(3));
            else
                y{1,j} = y{1,j}/(flagTV.*regularizerWeights(2) + flagNLTV.*regularizerWeights(4) + flagWave.*regularizerWeights(1)); % + flagGroup.*regularizerWeights(3));
            end;
        end;

        flag_fast = true;
        if flag_fast
            z{1,j}=y{1,j}+((t_old-1)/t_new).*(y{1,j}-y_old{1,j});
        end;
    end;

%% metrics of current itr:
%     disp(itr);
    dispProgress('Proximal Average', itr/maxitr);

   	metrics.xtime(itr+1)= toc(timer_WaT);
    for j = 1:nCha
        z{1,j} = z{1,j};
    end;
%     [metrics, ssim_map{1,2}] = get_metrics_itr( im_ref, im_ref_full, z, itr, maxitr, nCha, n1, n2, metrics, obj.K_1, obj.K_2, obj.W_size, obj.W_sigma );

end;
dispProgress('Proximal Average', 'Close');

    imageCha = z;
    for j = 1:nCha
        imageCha{1,j} = turn_image( imageCha{1,j} );
    end;
%     for j = 1:nCha+1
%         ssim_map{1,1}{1,j} = turn_image( ssim_map{1,1}{1,j} );
%         ssim_map{1,2}{1,j} = turn_image( ssim_map{1,2}{1,j} );
%     end;

end

author: thomaskuestner


text:
function [Flux, FBAsolution, model] = testPathway(model, MetIn, MetOut, AdditionalMetsInorOut, ObjectiveOption)
% Allows the user to see if given one metabolite `A`,
% downstream metabolite `B` can be made made. Additional sinks can be added
% for co-factors if needed
% `A -->-->-->-->--> B`
%
% USAGE:
%
%    [Flux, FBAsolution, model] = testPathway(model, MetIn, MetOut, AdditionalMetsInorOut, ObjectiveOption)
%
% INPUTS:
%    model:                    COBRA model structure
%    MetIn:                    The input metabolite(s) (`A`)
%    MetOut:                   The output metabolite (`B`)
%
% OPTIONAL INPUTS:
%    AdditionalMetsInorOut:    Additional metabolites for which sinks will be added
%    ObjectiveOption:          Boolean:
%
%                                * 1 = objective will be production of `B` (default)
%                                * 0 = use objective in model
%
% OUTPUTS:
%    Flux:                     The rate of `B` production
%    FBAsolution:              Solution to the FBA problem
%    model:                    COBRA model with sinks in it
%
% .. Author: Nathan Lewis, Feb 16 2009

if ~iscell(MetIn)
    Met = MetIn; clear MetIn; MetIn{1} = Met;
end
if ~iscell(MetOut)
    Met = MetOut; clear MetOut; MetOut{1} = Met;
end
if nargin > 3
    if ~iscell(AdditionalMetsInorOut)
        Met = AdditionalMetsInorOut; clear AdditionalMetsInorOut; AdditionalMetsInorOut{1} = Met;
    end
    % add sink rxns for all AdditionalMetsInorOut
    for i = 1:length(AdditionalMetsInorOut)
        model = addReaction(model,cat(2,'Tempsink_',AdditionalMetsInorOut{i}),{AdditionalMetsInorOut{i} },-1 ,true);
    end
end
if nargin <5,ObjectiveOption=1;end

for i = 1:length(MetIn) % add inputs
    model = addReaction(model,cat(2,'TempInput_',MetIn{i}),{MetIn{i} },1 ,false);
end
[model, rxnExists] = addReaction(model,cat(2,'TempOutput_',MetOut{1}),{MetOut{1} },-1 ,false);
if (ObjectiveOption==1 && isempty(rxnExists))
	model = changeObjective(model,cat(2,'TempOutput_',MetOut{1}));
elseif (ObjectiveOption == 1 && ~isempty(rxnExists))
	model = changeObjective(model,model.rxns(rxnExists));
end
FBAsolution = optimizeCbModel(model,'max');
Flux = FBAsolution.f;
if ~isempty(FBAsolution.x)
%printFluxVector(model,FBAsolution.x);
else display('zero flux in network')
end

author: opencobra

repo: cobratoolbox


text:
function arffwrite(fileName,dataName,attributeName,attributeType,data)
% ARFFWRITE  Writes the file as a arff formatted file.
%
% USAGE:
%       arffWrite(fileName,dataName,attributeName,attributeType,data);
%
% INPUT:
%       fileName:       file name for writing data
%       dataName:       relation name for arff file
%       attributeName:  attribute name for each variable { 1 by nAttr}
%       attributeType:  data type for each attribute { 1 by nAttr}
%       data:      		data for writing arff formatted(nInstan by nAttr)
%
% See also ARFFREAD

% Copyright 2004-2004 by Durga Lal Shrestha.
% eMail: durgals@hotmail.com
% $Date: 2004/06/23
% $Revision: 3.2.0 $ $Date: 2004/08/16 $

% ***********************************************************************
% Check for input data
if nargin < 5,
	error('Too few input arguments!');
end
if nargin > 5,
	error('Too many input arguments!');
end
nAttribute = size(data,2);
nVar = size(attributeName,2);
nVarType= size(attributeType,2);
if nAttribute ~= nVar | nAttribute ~=nVarType
	error('dimensions (column) of data must agree with number of varible name or type!');
end
% first check for heading
format = [];
for i=1:nAttribute-1
   format = [format ' %6.4f'];
end
format = [format ' %6.4f\n'];
%-------------------------------
fid = fopen(fileName,'w');          % open the file if exists otherwise create new
%-------------------------------

% Writing headings in the arff file format.
fprintf(fid,'%s %s\n','@relation',dataName);
for i=1:nAttribute
    fprintf(fid,'%s %s %s\n' ,'@attribute' ,attributeName{i} ,attributeType{i} );
end
fprintf(fid,'%s\n','@data');
%-------------------------------

% Data are space delimeted matrix
fprintf(fid,format,data' );
fclose(fid);

author: Sable

repo: mcbench-benchmarks


text:
function [map, aps, pr, prs] = compute_map (ranks, gnd, kappas)
% COMPUTE_MAP  This function computes the mAP for a given set of returned results.
%
% Usage:
%   map = compute_map (ranks, gnd)
%         computes mean average precsion (map) only
%
%   [map, aps, pr, prs] = compute_map (ranks, gnd, kappas)
%         computes mean average precision (map), average precision (aps) for each query
%         computes mean precision at kappas (pr), precision at kappas (prs) for each query
%
% Notes:
% 1) ranks starts from 1, size(ranks) = db_size X #queries
% 2) The junk results (e.g., the query itself) should be declared in the gnd stuct array
% 3) If there are no positive images for some query, that query is excluded from the evaluation

  if ~exist('kappas'), kappas = 0; end

  nq = numel (gnd);   % number of queries
  % init map and pr
  map = 0;
  aps = zeros (nq, 1);
  pr = zeros(1, numel(kappas));
  prs = zeros (nq, numel(kappas));
  nempty = 0;

  for i = 1:nq
    qgnd = gnd(i).ok;

    if isempty(qgnd) % no positive at all, skip from the average
      aps (i) = nan;
      prs (i, :) = nan;
      nempty = nempty + 1;
      continue;
    end

    if isfield (gnd(i), 'junk')
      qgndj = gnd(i).junk;
    else
      qgndj = [];
    end

    % positions of positive and junk images
    [~, pos] = intersect (ranks (:,i), qgnd);
    [~, junk] = intersect (ranks (:,i), qgndj);

    pos = sort(pos);
    junk = sort(junk);

    k = 0;
    ij = 1;

    if length (junk)
        % decrease positions of positives based on the number of junk images appearing before them
        ip = 1;
        while ip <= numel (pos)

            while ( ij <= length (junk) & pos (ip) > junk (ij) )
                k = k + 1;
                ij = ij + 1;
            end

            pos (ip) = pos (ip) - k;
            ip = ip + 1;
        end
    end

    % compute ap
    ap = score_ap_from_ranks1 (pos, length (qgnd));
    map = map + ap;
    aps (i) = ap;

    % compute precision@k
    for j = 1:numel(kappas)
      kq = min(max(pos), kappas(j));
      prs(i, j) = numel(find(pos <= kq)) ./ kq;
    end
    pr = pr + prs(i, :);

  end

  map = map / (nq-nempty);
  pr = pr / (nq-nempty);

end

% This function computes the AP for a query
function ap = score_ap_from_ranks1 (ranks, nres)

% number of images ranked by the system
nimgranks = length (ranks);
ranks = ranks - 1;

% accumulate trapezoids in PR-plot
ap = 0;

recall_step = 1 / nres;

for j = 1:nimgranks
  rank = ranks(j);

  if rank == 0
    precision_0 = 1.0;
  else
    precision_0 = (j - 1) / rank;
  end

  precision_1 = j / (rank + 1);
  ap = ap + (precision_0 + precision_1) * recall_step / 2;
end

end

author: filipradenovic

repo: cnnimageretrieval


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

%% Models

lambdaV = 5; % Expected number of clutter measurements over entire surveillance region
V = 10^2;     % Volume of surveillance region (10x10 2D-grid)
V_bounds = [0 10 0 10]; % [x_min x_max y_min y_max]
P_D = 0.5;

% Instantiate a Transitionamic model
transition_model = ConstantVelocityX('NumDims',2,'VelocityErrVariance',0.0001);

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

% Instantiate a clutter model
clutter_model = PoissonRateUniformPositionX('ClutterRate',lambdaV,'Limits',[V_bounds(1:2);V_bounds(3:4)]);

% Instantiate birth model
birth_model = DistributionBasedBirthModelX('Distribution', UniformDistributionX([V_bounds(1:2); ...
                                                                                [-0.1 0.1 ];...
                                                                                V_bounds(3:4);...
                                                                                [-0.1 0.1 ]]),...
                                           'BirthIntensity', 0.00005);
% Instantiate detection model
detection_model = ConstantDetectionProbabilityX('DetectionProbability',P_D);

% Compile the State-Space model
model = StateSpaceModelX(transition_model,...
                       measurement_model,...
                       'Clutter',clutter_model,...
                       'Birth', birth_model,...
                       'Detection', detection_model);

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

% Simulate some measurements from ground-truth data
MeasurementScans = dataSim.simulate(TrueTrack);
measurements = [MeasurementScans.Vectors];

%% Initiation

% Setup prior assuming we have some intuition of the target's initial
% position
priorParticles = TrueTrack.Trajectory(1).Vector ...
                 + model.Measurement.finv(model.Measurement.random(4000));
priorWeights = ones(1,4000)/4000;
config.Model = model;
config.StatePrior = ParticleStateX(priorParticles,priorWeights./sum(priorWeights));
config.StatePrior.Metadata.ExistenceProbability = 0.5;
config.BirthScheme = {'Expansion', 5000};
config.SurvivalProbability = 0.99;

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

%% Estimation
% Instantiate a filter objects
filter = BernoulliParticleFilterX(config);
figure;
for t = 2:NumIter

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

    % Perform filtering
    filter.MeasurementList = MeasurementList;
    filter.predict();
    filter.update();

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

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

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
% Fast computation of the vector q defined by
% q_k = \sum_{l = 1}^{N_y} G(X_k - Y_l) f_l, k = 1 .. Nx
% Using the "Efficient Bessel Decomposition"
% Code developped by Martin Averseng
% See also SCSD method by François Alouges and Matthieu Aussal.

addpath(genpath(pwd)); % Add folders of the toolbox to the path.
clear all;%#ok
close all;
clc;

Nx = 10^4;
Ny = 10^4;
fprintf('10^%s x 10^%s cloud of points \n \n',num2str(round(log(Nx)/log(10))),num2str(round(log(Ny)/log(10))))
% Data points
X = uniformCircle([0,0],1,Nx);
Y = X; %uniformCircle([0,0],1,Ny)+0.5; %uniformDisk([0.2,0],1,Ny);
f = rand(size(Y,1),1); % Vector f

% Parameter for rescaling
rMax = rMaxCalc(X,Y);

% Kernel choice:

% G = LogKernel; % G(x) = log(x)

G = Y0Kernel(0.3); % G(x) = Y0(0.3*x) => Bessel decomposition with Robin
% conditions.

% G = 3*Y0Kernel(4); % G(x) = Y0(2*x) => Method of rescaling to a root of Y0

% G = H0Kernel(1000); % G(x) = Y0(1000*x) => Selects frequencies near 0 and 1000

% G = Kernel(@(r)(exp(-r.^2)),@(r)(-2*r.*exp(-r.^2))); % Arbitrary (smooth)
% kernel

% G = Kernel(@(r)(1./r.^2 ),@(r)(-2./r.^3)); % Arbitrary (singular) kernel

% Choice of the cutoff parameter.
lambda = 3;
a = lambda/sqrt(sqrt(Nx*Ny)); %this value is of the order of the optimal
% value for data uniformly distributed in a disk. Choose lambda by trial
% and error to minimize the online time.

tol = 1e-3; % input tolerance

% Offline computations.
gradOpt = true;
timeOff = tic;
[dxGxy,dyGxy,rq,loc] = offline_dEBD(G,X,Y,a,tol);
timeOff = toc(timeOff);
fprintf('Offline computations performed in \n %s seconds \n',num2str(timeOff))
% show the radial quadrature :
% rq.showDer;

% Online procedure.
timeOn = tic;
qx = dxGxy(f);
qy = dyGxy(f);
timeOn = toc(timeOn);
fprintf('Online products computed in \n %s seconds \n',num2str(timeOn))

% Error on first entry of q
dist = sqrt((X(1,1) - Y(:,1)).^2 + (X(1,2) - Y(:,2)).^2);
distInv = 1./dist;
distInv(dist<1e-12) = 0;
Y_X_x = (Y(:,1) - X(1,1));
Y_X_y = (Y(:,2) - X(1,2));
qvalx = sum((G.evalDer(dist).*Y_X_x.*distInv).*f);
qvaly = sum((G.evalDer(dist).*Y_X_y.*distInv).*f);
disp('error on first entry');
disp(abs(qvalx - qx(1))/(norm(qx,1)));
disp(abs(qvaly - qy(1))/(norm(qy,1)));

% Error when f = [1 0 0 ... 0]
dist = sqrt((Y(1,1) - X(:,1)).^2 + (Y(1,2) - X(:,2)).^2);
distInv = 1./dist;
distInv(dist < 1e-12) = 0;
Y_X_x = (Y(1,1) - X(:,1));
Y_X_y = (Y(1,2) - X(:,2));
qvalx = (G.evalDer(dist).*Y_X_x.*distInv);
qvaly = (G.evalDer(dist).*Y_X_y.*distInv);
f = [1; zeros(size(Y,1)-1,1)];

qx = dxGxy(f);
qy = dyGxy(f);

fprintf('Linf error for f = [1 0 0 ... 0] \n (effective error / target accuracy) \n');
fprintf('%s / %s \n\n',num2str(max(abs(qvalx - qx))),num2str(tol))
fprintf('%s / %s \n\n',num2str(max(abs(qvaly - qy))),num2str(tol))

disp('Done');

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
close all
clear

%% topology graph
D = [1,-1,0,0,0,0,0,0,0,-1,0,1;
    -1,0,0,0,0,0,1,-1,0,0,0,0;
    0,1,-1,0,0,0,0,0,1,0,0,0;
    0,0,0,0,0,1,-1,0,0,1,-1,0;
    0,0,1,-1,0,0,0,0,0,0,1,-1;
    0,0,0,0,1,-1,0,0,-1,0,0,0;
    0,0,0,1,-1,0,0,1,0,0,0,0];
L = D*D';
H = D';

%% dimensions
[n,m] = size(D);
d = 2;
% r represents P(r)
r = [2,0;
    1,1;
    1,-1;
    0,1;
    0,-1;
    -1,1;
    -1,-1];
% P represents \bar P(r)
P = [r,ones(n,1)];
% edge set
edge = mod(reshape(find(D~=0),2,m),n);
edge(edge==0) = n;

%% target formation
figure
for i=1:m
    plot(r(edge(:,i),1),r(edge(:,i),2),'k','linewidth',2); hold on
end
for i=1:n
    plot(r(i,1),r(i,2),'.','markersize',50);
end
axis([-2 2 -2 2]);

%% trajectory
figure
via = [0,0;
    5,0;
    10,0;
    10,-5;
    10,-10;
    5,-10;
    0,-10];
for j=1:size(via,1)
    if ~mod(j,2)
        if j==4
            T1 = diag([0.1,1]);
        else
            T1 = diag([1,0.1]);
        end
    else
        T1 = eye(2);
    end
    T2 = rot2(-pi/2*floor((j-1)/2));
    ra(:,:,j) = r*T2'*T1'+via(j,:);
    qvia(j,:) = [vec(T1*T2)',via(j,:)];
    for i=1:m
        plot(ra(edge(:,i),1,j),ra(edge(:,i),2,j),'k','linewidth',2); hold on
    end
end
[qr,dqr,ddqr,tr] = mstraj_(qvia,ones(1,6),0.1,0.2);
A = reshape(qr(1,1:4)',[2,2]);
b = qr(1,5:6)';
xr = r*A'+b';
% for j=1:length(tr)
%     dA = reshape(dqr(j,1:4)',[2,2]);
%     db = dqr(j,5:6)';
%     dxr = r*dA'+db';
%     xr = xr+dxr*0.1;
%     plot(xr(1,1),xr(1,2),'bo');
%     plot(xr(2,1),xr(2,2),'go');
%     plot(xr(3,1),xr(3,2),'mo');
%     drawnow
% end

%% calculate stress matrix
E = [];
for i=1:n
    E = [E;P'*H'*diag(H(:,i))];
end
[U,S,V] = svd(P);
% the order in S: from big to small
U1 = U(:,1:d+1);
U2 = U(:,d+2:end);
% in fact z is exactly the same as omega given in Fig. 3
z = null(E);
% check if all(svd(U2'*H'*diag(z)*H*U2)>0), then
Omega = H'*diag(z)*H;

%% initial position
x = rand(n,d)*4-2;
dx = zeros(size(x));
% figure
for i=1:m
    fig_edge(i) = plot(x(edge(:,i),1),x(edge(:,i),2),'k','linewidth',2); hold on
end
for i=1:n
    fig_node(i) = plot(x(i,1),x(i,2),'.','markersize',50);
end
axis([-2 12 -12 2]);

%% simulation
dt = 0.05;
loop = 0;
video_on = true;
for t=tr(1):dt:tr(end)
    loop = loop+1;
    dq = interp1(tr,dqr,t);
    dA = reshape(dq(1:4)',[2,2]);
    db = dq(5:6)';
    dxr = r*dA'+db';
    % test leader
    alpha = 1;
%     dx = dxr-alpha*(x-xr);
    % followers apply control law (11)
    D = H';
    gamma = diag(Omega);
    for i=4:n
        err_sum = [0,0];
        edge_ind = find(D(i,:)~=0);
        for k=edge_ind
            node_ind = find(D(:,k)~=0);
            j = node_ind(node_ind~=i);
            err_sum = err_sum+z(k)*(x(i,:)-x(j,:)-dx(j,:));
        end
        dx(i,:) = -1/gamma(i)*err_sum;
    end
    % first 3 agents are leaders
    dx(1:3,:) = dxr(1:3,:)-alpha*(x(1:3,:)-xr(1:3,:));
    % update states
    x = x+dt*dx;
    xr = xr+dt*dxr;
    pos_data(:,:,loop) = x;
    % update figure
    for i=1:m
        set(fig_edge(i),'xdata',x(edge(:,i),1),'ydata',x(edge(:,i),2));
    end
    for i=1:n
        set(fig_node(i),'xdata',x(i,1),'ydata',x(i,2));
    end
    % video
    if video_on
        frame(loop) = getframe(gcf);
    end
    drawnow
end
% write video
if video_on
    savevideo('affine_maneuver',frame);
end
figure
t_data = tr(1):dt:tr(end);
xpos_data = squeeze(pos_data(:,1,:));
ypos_data = squeeze(pos_data(:,2,:));
plot3(kron(ones(n,1),t_data)',xpos_data',ypos_data');
xlabel('time/s');ylabel('x/m');zlabel('y/m');grid

repo: paper-simulation


text:
function idct = myidct1(vector)
%Determines the inverse discrete cosine tranform of a vector using an fft.
%
%Given a row vector, returns the inverse discrete cosine transformation of
%the vector. Since the DCT-I is its own inverse, this simply entails
%renormalizing.

    idct = (2/(numel(vector)-1))*mydct1(vector);

end


text:
function f1 = p07_f1 ( x )

%*****************************************************************************80
%
%% P07_F1 evaluates the first derivative for problem 7.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the value of the variable.
%
%    Output, real F1, the first derivative of the
%    objective function.
%
  f1 = ( 1.0 - 2.0 * x * x + cos ( x ) ...
         - 2.0 * x * sin ( x ) ) * exp ( - x * x );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function p=polelague(n)

% p=polegend(n)
% Almacena en las filas de la matriz p los coefs de los polinomios de Legendre

p(1,1)=1;
p(2,1:2)=[-1 1];
for k=2:n
   p(k+1,1:k+1)=((2*(k-2)*[0 p(k,1:k)]+3*[0 p(k,1:k)]-[p(k,1:k) 0]-(k-1).^2*[0 0 p(k-1,1:k-1)]));
end


text:
function mpc = t_case_int
%T_CASE_INT  Case data in internal format used to test EXT2INT and INT2EXT.

%   MATPOWER

%% MATPOWER Case Format : Version 2
mpc.version = '2';

%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = 100;

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [
	1	3	0	0	0	0	1	1	0	345	1	1.1	0.9;
	2	2	0	0	0	0	1	1	0	345	1	1.1	0.9;
	3	2	0	0	0	0	1	1	0	345	1	1.1	0.9;
	4	1	0	0	0	0	1	1	0	345	1	1.1	0.9;
	5	1	90	30	0	0.2	1	1	0	345	1	1.1	0.9;
	6	1	0	0	0	0	1	1	0	345	1	1.1	0.9;
	7	1	100	35	0	0	1	1	0	345	1	1.1	0.9;
	8	1	0	0	0.1	0	1	1	0	345	1	1.1	0.9;
	9	1	125	50	0	0	1	1	0	345	1	1.1	0.9;
];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
mpc.gen = [
	3	85	0	300	-300	1	100	1	270	10	0	0	0	0	0	0	0	0	0	0	0;
	2	163	0	300	-300	1	100	1	300	10	0	0	0	0	0	0	0	0	0	0	0;
	1	0	0	300	-300	1	100	1	250	90	0	0	0	0	0	0	0	0	0	0	0;
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [
	1	4	0	0.0576	0	0	250	250	0	0	1	-360	360;
	4	5	0.017	0.092	0.158	0	250	250	0	0	1	-360	360;
	5	6	0.039	0.17	0.358	150	150	150	0	0	1	-360	360;
	3	6	0	0.0586	0	0	300	300	0	0	1	-360	360;
	6	7	0.0119	0.1008	0.209	40	150	150	0	0	1	-360	360;
	7	8	0.0085	0.072	0.149	250	250	250	0	0	1	-360	360;
	8	2	0	0.0625	0	250	250	250	0	0	1	-360	360;
	8	9	0.032	0.161	0.306	250	250	250	0	0	1	-360	360;
	9	4	0.01	0.085	0.176	250	250	250	0	0	1	-360	360;
];

%%-----  OPF Data  -----%%
%% generator cost data
%	1	startup	shutdown	n	x1	y1	...	xn	yn
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [
	2	0	0	2	15	0	0	0	0	0	0	0;
	1	0	0	4	0	0	100	2500	200	5500	250	7250;
	1	0	0	4	0	0	100	2000	200	4403.5	270	6363.5;
];

%% bus names
mpc.bus_name = {
            'BUS 1';
            'BUS 2';
            'BUS 30';
            'BUS 4';
            'BUS 5';
            'BUS 6';
            'BUS 7';
            'BUS 8';
            'BUS 9';
};

%% generator unit type (see GENTYPES)
mpc.gentype = {
	'ST';
	'GT';
	'HY';
};

%% generator fuel type (see GENFUELS)
mpc.genfuel = {
	'coal';
	'ng';
	'hydro';
};

mpc.A = [
	1	2	3	4	5	7	8	9	10	11	12	13	14	15	17	18	19	20	21	22	24	25	26	28	29	30;
	2	4	6	8	10	14	16	18	20	22	24	26	28	30	34	36	38	40	42	44	48	50	52	56	58	60;
];

mpc.N = [
	30	29	28	27	26	24	23	22	21	20	19	18	17	16	14	13	12	11	10	9	7	6	5	3	2	1;
	60	58	56	54	52	48	46	44	42	40	38	36	34	32	28	26	24	22	20	18	14	12	10	6	4	2;
];

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
classdef TankMassFlowRateConstraint < AbstractConstraint
    %TankMassFlowRateConstraint Summary of this class goes here
    %   Detailed explanation goes here

    properties
        normFact = 1;
        tank LaunchVehicleTank
        event LaunchVehicleEvent
        eventNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;

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

        evalType(1,1) ConstraintEvalTypeEnum = ConstraintEvalTypeEnum.FixedBounds;
        stateCompType(1,1) ConstraintStateComparisonTypeEnum = ConstraintStateComparisonTypeEnum.Equals;
        stateCompEvent LaunchVehicleEvent
        stateCompNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;
    end

    methods
        function obj = TankMassFlowRateConstraint(tank, event, lb, ub)
            obj.tank = tank;
            obj.event = event;
            obj.lb = lb;
            obj.ub = ub;

             obj.id = rand();
        end

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

        function [c, ceq, value, lwrBnd, uprBnd, type, eventNum, valueStateComp] = evalConstraint(obj, stateLog, celBodyData)
            type = obj.getConstraintType();

            switch obj.eventNode
                case ConstraintStateComparisonNodeEnum.FinalState
                    stateLogEntry = stateLog.getLastStateLogForEvent(obj.event);

                case ConstraintStateComparisonNodeEnum.InitialState
                    stateLogEntry = stateLog.getFirstStateLogForEvent(obj.event);

                otherwise
                    error('Unknown event node.');
            end

            value = lvd_TankMassTasks(stateLogEntry, 'tankMDot', obj.tank);

            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                switch obj.stateCompNode
                    case ConstraintStateComparisonNodeEnum.FinalState
                        stateLogEntryStateComp = stateLog.getLastStateLogForEvent(obj.stateCompEvent);

                    case ConstraintStateComparisonNodeEnum.InitialState
                        stateLogEntryStateComp = stateLog.getFirstStateLogForEvent(obj.stateCompEvent);

                    otherwise
                        error('Unknown event node.');
                end

                valueStateComp = lvd_TankMassTasks(stateLogEntryStateComp, 'tankMDot', obj.tank);
            else
                valueStateComp = NaN;
            end

            [c, ceq] = obj.computeCAndCeqValues(value, valueStateComp);

            lwrBnd = obj.lb;
            uprBnd = obj.ub;

            eventNum = obj.event.getEventNum();
        end

        function sF = getScaleFactor(obj)
            sF = obj.normFact;
        end

        function setScaleFactor(obj, sF)
            obj.normFact = sF;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = obj.tank == tank;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesEvent(obj, event)
            tf = obj.event == event;
            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                tf = tf || obj.stateCompEvent == event;
            end
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end

        function tf = usesExtremum(obj, extremum)
            tf = false;
        end

        function tf = usesGroundObj(obj, grdObj)
            tf = false;
        end

        function tf = canUseSparseOutput(obj)
            tf = true;
        end

        function event = getConstraintEvent(obj)
            event = obj.event;
        end

        function type = getConstraintType(obj)
            type = 'Tank Mass Flow Rate';
        end

        function [unit, lbLim, ubLim, usesLbUb, usesCelBody, usesRefSc] = getConstraintStaticDetails(obj)
            unit = 'mT/s';
            lbLim = -Inf;
            ubLim = Inf;
            usesLbUb = true;
            usesCelBody = false;
            usesRefSc = false;
        end

        function addConstraintTf = openEditConstraintUI(obj, lvdData)
            [~, tanks] = lvdData.launchVehicle.getTanksListBoxStr();
            if(numel(tanks) >= 1)
%                 addConstraintTf = lvd_EditTankConstraintGUI(obj, lvdData);

                output = AppDesignerGUIOutput({false});
                lvd_EditTankConstraintGUI_App(obj, lvdData, output);
                addConstraintTf = output.output{1};
            else
                errordlg('There are currently no tanks assigned to the launch vehicle in this scenario.  Add at least one tank first.');

                addConstraintTf = false;
            end
        end
    end

    methods(Static)
        function constraint = getDefaultConstraint(~, ~)
            constraint = TankMassFlowRateConstraint(LaunchVehicleTank.empty(1,0), LaunchVehicleEvent.empty(1,0), 0, 0);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%MTIT		creates a major title in a figure with many axes
%
%		MTIT
%		- creates a major title above all
%		  axes in a figure
%		- preserves the stack order of
%		  the axis handles
%
%SYNTAX
%-------------------------------------------------------------------------------
%		P = MTIT(TXT,[OPT1,...,OPTn])
%		P = MTIT(FH,TXT,[OPT1,...,OPTn])
%
%INPUT
%-------------------------------------------------------------------------------
%    FH	:	a valid figure handle		[def: gcf]
%   TXT	:	title string
%
% OPT	:	argument
% -------------------------------------------
%  xoff	:	+/- displacement along X axis
%  yoff	:	+/- displacement along Y axis
%  zoff	:	+/- displacement along Z axis
%
%		title modifier pair(s)
% -------------------------------------------
%   TPx	:	TVx
%		see: get(text) for possible
%		     parameters/values
%
%OUTPUT
%-------------------------------------------------------------------------------
% par	:	parameter structure
%  .pos :	position of surrounding axis
%   .oh	:	handle of last used axis
%   .ah :	handle of invisible surrounding axis
%   .th :	handle of main title
%
%EXAMPLE
%-------------------------------------------------------------------------------
%	subplot(2,3,[1 3]);		title('PLOT 1');
%	subplot(2,3,4); 		title('PLOT 2');
%	subplot(2,3,5); 		title('PLOT 3');
%	axes('units','inches',...
%	     'color',[0 1 .5],...
%	     'position',[.5 .5 2 2]);	title('PLOT 41');
%	axes('units','inches',...
%	     'color',[0 .5 1],...
%	     'position',[3.5 .5 2 2]);	title('PLOT 42');
%	shg;
%	p=mtit('the BIG title',...
%	     'fontsize',14,'color',[1 0 0],...
%	     'xoff',-.1,'yoff',.025);
% % refine title using its handle <p.th>
%	set(p.th,'edgecolor',.5*[1 1 1]);

% created:
%	us	24-Feb-2003		/ R13
% modified:
%	us	24-Feb-2003		/ CSSM
%	us	06-Apr-2003		/ TMW
%	us	13-Nov-2009 17:38:17

%-------------------------------------------------------------------------------
function	par=mtit(varargin)

		defunit='normalized';
	if	nargout
		par=[];
	end

% check input
	if	nargin < 1
		help(mfilename);
		return;
	end
	if	isempty(get(0,'currentfigure'))
		disp('MTIT> no figure');
		return;
	end

		vl=true(size(varargin));
	if	ischar(varargin{1})
		vl(1)=false;
		figh=gcf;
		txt=varargin{1};
	elseif	any(ishandle(varargin{1}(:)))		&&...
		ischar(varargin{2})
		vl(1:2)=false;
		figh=varargin{1};
		txt=varargin{2};
	else
		error('MTIT> invalid input');
	end
		vin=varargin(vl);
		[off,vout]=get_off(vin{:});

% find surrounding box
		ah=findall(figh,'type','axes');
	if	isempty(ah)
		disp('MTIT> no axis');
		return;
	end
		oah=ah(1);

		ou=get(ah,'units');
		set(ah,'units',defunit);
		ap=get(ah,'position');
	if	iscell(ap)
		ap=cell2mat(get(ah,'position'));
	end
		ap=[	min(ap(:,1)),max(ap(:,1)+ap(:,3)),...
			min(ap(:,2)),max(ap(:,2)+ap(:,4))];
		ap=[	ap(1),ap(3),...
			ap(2)-ap(1),ap(4)-ap(3)];

% create axis...
		xh=axes('position',ap);
% ...and title
		th=title(txt,vout{:});
		tp=get(th,'position');
		set(th,'position',tp+off);
		set(xh,'visible','off','hittest','on');
		set(th,'visible','on');

% reset original units
		ix=find(~strcmpi(ou,defunit));
	if	~isempty(ix)
	for	i=ix(:).'
		set(ah(i),'units',ou{i});
	end
	end

% ...and axis' order
		uistack(xh,'bottom');
		axes(oah);				%#ok

	if	nargout
		par.pos=ap;
		par.oh=oah;
		par.ah=xh;
		par.th=th;
	end
end
%-------------------------------------------------------------------------------
function	[off,vout]=get_off(varargin)

% search for pairs <.off>/<value>

		off=zeros(1,3);
		io=0;
	for	mode={'xoff','yoff','zoff'};
		ix=strcmpi(varargin,mode);
	if	any(ix)
		io=io+1;
		yx=find(ix);
		ix(yx+1)=1;
		off(1,io)=varargin{yx(end)+1};
		varargin=varargin(xor(ix,1));
	end
	end
		vout=varargin;
end
%-------------------------------------------------------------------------------

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function w = ymdf_to_weekday_republican ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_WEEKDAY_REPUBLICAN returns the weekday of a Republican YMDF date.
%
%  Discussion:
%
%    The Republican calendar used a 10 day week.
%    There was a final "month" of 5 or 6 days.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, integer W, the day of the week of the date.
%    The days are numbered from Sunday through Saturday, 1 through 7.
%
  w = i4_wrap ( d, 1, 10 );

  return
end


text:
function [dets, boxes, t] = cascade_detect(pyra, model, thresh)

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2009-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

th = tic();

% gather PCA root filters for convolution
numrootfilters = length(model.rootfilters);
rootfilters = cell(numrootfilters, 1);
for i = 1:numrootfilters
  rootfilters{i} = model.rootfilters{i}.wpca;
end

% compute PCA projection of the feature pyramid
projpyra = project_pyramid(model, pyra);

% stage 0: convolution with PCA root filters is done densely
% before any pruning can be applied

% Precompute location/scale scores
loc_f      = loc_feat(model, pyra.num_levels);
loc_scores = cell(model.numcomponents, 1);
for c = 1:model.numcomponents
  loc_w         = model.loc{c}.w;
  loc_scores{c} = loc_w * loc_f;
end
pyra.loc_scores = loc_scores;

numrootlocs = 0;
nlevels = size(pyra.feat,1);
rootscores = cell(model.numcomponents, nlevels);
s = 0;  % will hold the amount of temp storage needed by cascade()
for i = 1:pyra.num_levels
  s = s + size(pyra.feat{i},1)*size(pyra.feat{i},2);
  if i > model.interval
    scores = fconv_var_dim(projpyra.feat{i}, rootfilters, 1, numrootfilters);
    for c = 1:model.numcomponents
      u = model.components{c}.rootindex;
      v = model.components{c}.offsetindex;
      rootscores{c,i} = scores{u} + model.offsets{v}.w + loc_scores{c}(i);
      numrootlocs = numrootlocs + numel(scores{u});
    end
  end
end
s = s*length(model.partfilters);
model.thresh = thresh;
% run remaining cascade stages and collect object hypotheses
coords = cascade(model, pyra, projpyra, rootscores, numrootlocs, s);

boxes = coords';
dets = boxes(:,[1:4 end-1 end]);
t = toc(th);

author: rbgirshick

repo: voc-dpm

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


text:
function out = times( a,b )
%TIMES (element-wise) for nested cells
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

% determine input types
inType = {class(a), class(b)};
if(strcmp(inType{1},'TRAFO') && strcmp(inType{2},'TRAFO'))
    [outA, idx] = flattenCellMatrix(a.data);
    [outB, idxB] = flattenCellMatrix(b.data);
    meta = a.meta;
    meta_b = b.meta;
    clear 'a' 'b'

    if(~isempty(idx) && ~isempty(idxB) && ~all(cellfun(@(x,y) isequal(x,y), idx, idxB)))
        error('TRAFO::times: Nested cells must have the same size');
    end
    if (iscell(meta) && iscell(meta_b))
        if(~isempty(meta) && ~isempty(meta_b))
            if(iscell(meta{1}) && iscell(meta_b{1}))
                if(~all(cellfun(@(x,y) isequal(x,y), flattenCellMatrix(meta), flattenCellMatrix(meta_b))))
                    error('TRAFO::times: Unequal reconstruction information');
                end
            else
                if(~all(cellfun(@(x,y) isequal(x,y), meta, meta_b)))
                    % happens in kernelMult for kernel .* zpad(W) -> take kernel
                    % meta information (meta) => corrected
                    [ST, ~] = dbstack;
                    if(~strcmp(ST(2,1).name,'kernelMult'))
                        error('TRAFO::times: Unequal reconstruction information');
                    end
                end
            end
        end
    else
        if ~isequal(meta, meta_b)
            error('TRAFO::times: Unequal reconstruction information');
        end
    end

%     try
        out = cellfun(@(x,y) x.*y, outA, outB, 'UniformOutput', false);
%     catch msg
%         out = cell(size(outA));
%         for i=1:length(out)
%             out{i} = outA{i} .* outB{i};
%         end
%     end

elseif(strcmp(inType{1},'TRAFO')) % assume b is double/uint/int
    [out, idx] = flattenCellMatrix(a.data);
    meta = a.meta;
    clear 'a'

    out = cellfun(@(x) x .* b, out, 'UniformOutput', false);

elseif(strcmp(inType{2},'TRAFO'))
    [out, idx] = flattenCellMatrix(b.data);
    meta = b.meta;
    clear 'b'

    out = cellfun(@(x) a .* x, out, 'UniformOutput', false);

else
    error('TRAFO::times: Impossible constellation!');
end

out = reconFlatCellMatrix(out,idx);
out = TRAFO(out,meta); % return TRAFO object again without modifying input

end


text:
function v = logdet(A, op)
%LOGDET Computation of logarithm of determinant of a matrix
%
%   v = logdet(A);
%       computes the logarithm of determinant of A.
%
%       Here, A should be a square matrix of double or single class.
%       If A is singular, it will returns -inf.
%
%       Theoretically, this function should be functionally
%       equivalent to log(det(A)). However, it avoids the
%       overflow/underflow problems that are likely to
%       happen when applying det to large matrices.
%
%       The key idea is based on the mathematical fact that
%       the determinant of a triangular matrix equals the
%       product of its diagonal elements. Hence, the matrix's
%       log-determinant is equal to the sum of their logarithm
%       values. By keeping all computations in log-scale, the
%       problem of underflow/overflow caused by product of
%       many numbers can be effectively circumvented.
%
%       The implementation is based on LU factorization.
%
%   v = logdet(A, 'chol');
%       If A is positive definite, you can tell the function
%       to use Cholesky factorization to accomplish the task
%       using this syntax, which is substantially more efficient
%       for positive definite matrix.
%
%   Remarks
%   -------
%       logarithm of determinant of a matrix widely occurs in the
%       context of multivariate statistics. The log-pdf, entropy,
%       and divergence of Gaussian distribution typically comprises
%       a term in form of log-determinant. This function might be
%       useful there, especially in a high-dimensional space.
%
%       Theoretially, LU, QR can both do the job. However, LU
%       factorization is substantially faster. So, for generic
%       matrix, LU factorization is adopted.
%
%       For positive definite matrices, such as covariance matrices,
%       Cholesky factorization is typically more efficient. And it
%       is STRONGLY RECOMMENDED that you use the chol (2nd syntax above)
%       when you are sure that you are dealing with a positive definite
%       matrix.
%
%   Examples
%   --------
%       % compute the log-determinant of a generic matrix
%       A = rand(1000);
%       v = logdet(A);
%
%       % compute the log-determinant of a positive-definite matrix
%       A = rand(1000);
%       C = A * A';     % this makes C positive definite
%       v = logdet(C, 'chol');
%

%   Copyright 2008, Dahua Lin, MIT
%   Email: dhlin@mit.edu
%
%   This file can be freely modified or distributed for any kind of
%   purposes.
%

%% argument checking

assert(isfloat(A) && ndims(A) == 2 && size(A,1) == size(A,2), ...
    'logdet:invalidarg', ...
    'A should be a square matrix of double or single class.');

if nargin < 2
    use_chol = 0;
else
    assert(strcmpi(op, 'chol'), ...
        'logdet:invalidarg', ...
        'The second argument can only be a string ''chol'' if it is specified.');
    use_chol = 1;
end

%% computation

if use_chol
    v = 2 * sum(log(diag(chol(A))));
else
    [L, U, P] = lu(A);
    du = diag(U);
    c = det(P) * prod(sign(du));
    v = log(c) + sum(log(abs(du)));
end

author: OHBA-analysis

repo: HMM-MAR


function varargout = plot_function(varargin)
% PLOT_FUNCTION M-file for plot_function.fig
%      PLOT_FUNCTION, by itself, creates a new PLOT_FUNCTION or raises the existing
%      singleton*.
%
%      H = PLOT_FUNCTION returns the handle to a new PLOT_FUNCTION or the handle to
%      the existing singleton*.
%
%      PLOT_FUNCTION('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in PLOT_FUNCTION.M with the given input arguments.
%
%      PLOT_FUNCTION('Property','Value',...) creates a new PLOT_FUNCTION or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before plot_function_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to plot_function_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help plot_function

% Last Modified by GUIDE v2.5 20-Mar-2004 17:42:25

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @plot_function_OpeningFcn, ...
                   'gui_OutputFcn',  @plot_function_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before plot_function is made visible.
function plot_function_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to plot_function (see VARARGIN)

% Choose default command line output for plot_function
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes plot_function wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = plot_function_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

function a_Callback(hObject, eventdata, handles)
% hObject    handle to a (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of a as text
%        str2double(get(hObject,'String')) returns contents of a as a double

% --- Executes during object creation, after setting all properties.
function a_CreateFcn(hObject, eventdata, handles)
% hObject    handle to a (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function b_Callback(hObject, eventdata, handles)
% hObject    handle to b (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of b as text
%        str2double(get(hObject,'String')) returns contents of b as a double

% --- Executes during object creation, after setting all properties.
function b_CreateFcn(hObject, eventdata, handles)
% hObject    handle to b (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function c_Callback(hObject, eventdata, handles)
% hObject    handle to c (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of c as text
%        str2double(get(hObject,'String')) returns contents of c as a double

% --- Executes during object creation, after setting all properties.
function c_CreateFcn(hObject, eventdata, handles)
% hObject    handle to c (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function xmin_Callback(hObject, eventdata, handles)
% hObject    handle to xmin (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of xmin as text
%        str2double(get(hObject,'String')) returns contents of xmin as a double

% --- Executes during object creation, after setting all properties.
function xmin_CreateFcn(hObject, eventdata, handles)
% hObject    handle to xmin (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function xmax_Callback(hObject, eventdata, handles)
% hObject    handle to xmax (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of xmax as text
%        str2double(get(hObject,'String')) returns contents of xmax as a double

% --- Executes during object creation, after setting all properties.
function xmax_CreateFcn(hObject, eventdata, handles)
% hObject    handle to xmax (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in popupmenu1.
function popupmenu1_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns popupmenu1 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu1

% --- Executes during object creation, after setting all properties.
function popupmenu1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popupmenu1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in popupmenu2.
function popupmenu2_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns popupmenu2 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu2

% --- Executes during object creation, after setting all properties.
function popupmenu2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popupmenu2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function edit7_Callback(hObject, eventdata, handles)
% hObject    handle to edit7 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit7 as text
%        str2double(get(hObject,'String')) returns contents of edit7 as a double

% --- Executes during object creation, after setting all properties.
function edit7_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit7 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end


text:
function desired_velocity = variable_v(position,velocity)
%{
This function  receives the following inputs:
    1. Position of the aircraft- [x, y, z]
    2. Velocity of the aircraft - [x,y,z]
Using the inputs, the function returns the value of the desired velocity,
cartasian coordinates.
Created by Roni Peer, 11.7.2010
%}
%#eml
height = position(3);                                           % Get the Height.
if height > 500                                                         % Maintain velocity.
    desired_descent_rate = velocity(3);
elseif height > 200                                                % Slow down.
    desired_descent_rate = -10;
elseif height > 100
    desired_descent_rate = -2;
elseif height > 10
    desired_descent_rate = -1;
else
    desired_descent_rate = -0.5;                        % Break.
end
desired_velocity = [0; 0; desired_descent_rate];


text:
function [is_planar ksubgraph EI]=boyer_myrvold_planarity_test(A,varargin)
% BOYER_MYRVOLD_PLANARITY_TEST Test a graph for planarity
%
% is_planar = boyer_myrvold_planaity_test(A) yields 1 if A is a planar
% graph or 0 otherwise.  A planar graph can be drawn on the xy plane with
% no edge crossings.
%
% [is_planar K] = ... identifies a Kuratowski subgraph of A when A is not
% planar.
% [is_planar K EI] = ... computes the planar embedding edge order in
% EI.edge_order.  EI.vp gives the starting and ending entries in EI.edge_order
% for the order of edges from each vertex.  (EI is the compressed
% sparse row representation of the matrix A with the order of the column
% indices permuted to the planar graph embedding order.)
%
% ... = boyer_myrvold_planaity_test(A,...) takes a set of
% key-value pairs or an options structure.  See set_matlab_bgl_options
% for the standard options.
%   No additional options for this function
%
% Example:
%   G = grid_graph(6,5);
%   boyer_myrvold_planarity_test(G) % G is planar
%   K5 = clique_graph(5);
%   boyer_myrvold_planarity_test(K5) % K5 is not planar

% David Gleich
% Copyright, Stanford University, 2008

%% History
%  2007-10-06: Initial coding
%%

[trans check full2sparse] = get_matlab_bgl_options(varargin{:});
if full2sparse && ~issparse(A), A = sparse(A); end

options = struct();
options = merge_options(options,varargin{:});

if nargout <= 1
    is_planar = planar_test_mex(A,0);
else
    if nargout <= 2
        [is_planar ki kj] = planar_test_mex(A,0);
    else
        [is_planar ki kj eip eie] = planar_test_mex(A,0);
        EI = struct('vp',eip,'edge_order',eie);
    end
    if ~isempty(ki)
        ksubgraph = sparse(ki,kj,1,size(A,1),size(A,2));
    else
        ksubgraph = sparse([]);
    end
    ksubgraph = ksubgraph|ksubgraph'; % placed here to get the type right
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [varargout]=spmax(varargin)

% function [maxVal,maxInd]=spmax(A,B,vecdim,nanflag,logicRelevant,nanOut)
% ------------------------------------------------------------------------
%
%
% Kevin Moerman
% ------------------------------------------------------------------------

%%

A=varargin{1};

if nargin>=3
    dimDir=varargin{3};
else
    dimDir=1;
end

switch nargin
    case 5
        logicRelevant=varargin{5};
        varargin=varargin(1:4);
        nanOut=0;
        if isempty(varargin{4})
            varargin=varargin(1:3);
        end
    case 6
        logicRelevant=varargin{5};
        nanOut=varargin{6};
        varargin=varargin(1:4);
        if isempty(varargin{4})
            varargin=varargin(1:3);
        end
    otherwise
        logicRelevant=[];
        nanOut=0;
end

if isempty(logicRelevant)
    logicRelevant=A~=0;
end

logicRelevantRow=full(any(logicRelevant,dimDir));

%%

minOffset=min(A(logicRelevant));

A(logicRelevant)=A(logicRelevant)-minOffset; %Shift to negative to zeros loose
varargin{1}=A;
[maxVal,maxInd]=max(varargin{:});
maxVal(logicRelevantRow)=maxVal(logicRelevantRow)+minOffset; %Shift back
if nanOut==1
    maxVal(~logicRelevantRow)=NaN;
end

%% Collect output
varargout{1}=maxVal;
varargout{2}=maxInd;

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% APM Load Data File
%
% response = csv_load(server,app,filename)
%
% Function csv_load uploads the data file (csv) to the web-server
%   with the following arguments:
%
%   server = address of server
%      app = application name
% filename = data filename
%
% A response is returned indicating whether the file was
%   uploaded successfully to the server
%
function [response] = csv_load(server,app,filename)
   % convert to lowercase and deblank
   app = lower(deblank(app));

   % load model
   fid=fopen(filename,'r');
   tline = [];
   while 1
      aline = fgets(fid);
      if ~ischar(aline), break, end
      % remove any double quote marks
      aline = [strrep(aline,'"',' ')];
      tline = [tline aline];
   end
   fclose(fid);

   % send to server once for every 2000 characters
   ts = size(tline,2);
   block = 2000;
   cycles = ceil(ts/block);
   for i = 1:cycles,
      if i<cycles,
         csv_block = ['csva ' tline((i-1)*block+1:i*block)];
      else
         csv_block = ['csv ' tline((i-1)*block+1:end)];
      end
      response = apm(server,app,csv_block);
   end
   response = 'Successfully loaded CSV file';

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function value = monomial_value ( dim_num, point_num, x, expon )

%*****************************************************************************80
%
%% MONOMIAL_VALUE evaluates a monomial.
%
%  Discussion:
%
%    This routine evaluates a monomial of the form
%
%      product ( 1 <= dim <= dim_num ) x(dim)^expon(dim)
%
%    where the exponents are nonnegative integers.  Note that
%    if the combination 0^0 is encountered, it should be treated
%    as 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer POINT_NUM, the number of points at which the
%    monomial is to be evaluated.
%
%    Input, real X(DIM_NUM,POINT_NUM), the point coordinates.
%
%    Input, integer EXPON(DIM_NUM), the exponents.
%
%    Output, real VALUE(POINT_NUM), the value of the monomial.
%
  value(1,1:point_num) = 1.0;

  for dim = 1 : dim_num
    if ( 0 ~= expon(dim) )
      value(1,1:point_num) = value(1,1:point_num) .* ( x(dim,1:point_num).^expon(dim) );
    end
  end

  return
end


text:
%kkasc2mask 'Import ASCII Data File to Mask Segment'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kasc2mask.pane file
%
% Parameters:
% InputFile: i1 'ASCII Input File', required: 'ASCII input data file containing new data segment'
% InputFile: i2 'Data Object Input', optional: 'data object input to insert the new segment into'
% Integer: so 'Start Offset ', default: 0: 'starting point in ASCII file to start importing at'
% Integer: sf 'Skip Factor ', default: 0: 'number of points in ASCII file to skip between reads'
% Integer: wsize ' ', default: 1: 'width size must be > 0'
% Integer: hsize ' ', default: 1: 'height size must be > 0'
% Integer: dsize ' ', default: 1: 'depth size must be > 0'
% Integer: tsize ' ', default: 1: 'time size must be > 0'
% Integer: esize ' ', default: 1: 'elements size must be > 0'
% OutputFile: o 'Output Object', required: 'output data object'
%
% Example: o = kkasc2mask({i1, i2}, {'i1','';'i2','';'so',0;'sf',0;'wsize',1;'hsize',1;'dsize',1;'tsize',1;'esize',1;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kasc2mask - Import ASCII Data File to Mask Segment
%
%  DESCRIPTION
% This routine will import an ASCII data file and insert it into a data
% object if it is specified.  If there is no source object to insert the
% new segment into, the data object specified as the output file will be
% created.  The new segment will be inserted into this data object in the
% mask segment.  If a source object was specified, then all of its
% data segments and the global object attributes will be copied to the
% destination.
% The index order is specified by setting the -windex, -hindex, -dindex,
% -tindex and -eindex flags.  They can not be set to the same value.
% The start offset is the point to start importing in the
% ASCII data.  All points before this one are ignored.
% The skip factor is the number of points to skip in the ASCII
% data between reads.  If the size of the new segment's lowest
% index is 40 and the skip factor is 5, then this routine will
% read in 40 points skip 5 and read in 40 more, etc. until the
% entire ASCII file read in.
% Complex values in ASCII are represented by number pairs, surrounded
% by parenthesis and separated by a comma.  For example, (8, 3) represents
% the complex number 8+3i.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% kasc2val, kasc2loc, kasc2time, kasc2map, kimportasc
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkasc2mask(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkasc2mask(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'i2', '__input';'so', 0;'sf', 0;'wsize', 1;'hsize', 1;'dsize', 1;'tsize', 1;'esize', 1;'o', '__output'};
maxval={0,1,1,1,2,2,2,2,2,0};
minval={0,1,1,1,2,2,2,2,2,0};
istoggle=[0,1,1,1,1,1,1,1,1,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','Integer','Integer','Integer','Integer','Integer','Integer','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kasc2mask"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function sts = isequaltarget(dep1, dep2)

% function sts = isequaltarget(dep1, dep2)
% Compare source references of two dependencies and return true if both
% point to the same object. If multiple dependencies are given, the
% number and order of dependencies must match.
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

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

rev = '$Rev: 1716 $'; %#ok
sts = numel(dep1) == numel(dep2);
if sts
    for k = 1:numel(dep1)
        sts = sts && isequal(dep1.tgt_exbranch, dep2.tgt_exbranch) && ...
              isequal(dep1.tgt_input, dep2.tgt_input);
        if ~sts
            break;
        end;
    end;
end;

author: spm


text:
function rx = rxSaveBestrotvol(rx,savePath);
%
% rxSaveBestrotvol([rx],[savePath]);
%
% Save a mrAlign3 bestrotvol.mat file.
%
% ras 03/05
if ieNotDefined('rx')
    cfig = findobj('Tag','rxControlFig');
    rx = get(cfig,'UserData');
end

if ieNotDefined('savePath')
    savePath = fullfile(pwd,'bestrotvol.mat');
end

newXform = rx.xform;

A = newXform(1:3,1:3);
b = newXform(1:3,4)';

% account for scale factors
scaleFac(1,:) = 1./rx.rxVoxelSize;
scaleFac(2,:) = 1./rx.volVoxelSize;
rot = diag(1./scaleFac(2,:))*A*diag(scaleFac(1,:));
trans = b ./ scaleFac(2,:);

save(savePath,'rot','trans','scaleFac');
fprintf('Saved %s.\n',savePath);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ya, yd] = dwt_process(x, lf, hf, num)
x = double(x);
ya = x;
yd = [];
for i = 1 : num
    yli = conv(ya,lf);
    yai = downsample_prcoess(yli);
    yhi = conv(ya, hf);
    ydi = downsample_prcoess(yhi);
    ya = yai;
    yd = [yd ydi];
end

author: decouples


text:
function [yn,en,S] = PRAadapt(un,dn,S)

% PRAadapt          Partial Rank Algorithm (See Section 6.3)
%
% Arguments:
% un                Input signal
% dn                Desired signal
% S                 Adptive filter parameters as defined in PRAinit.m
% yn                History of output signal
% en                History of error signal
%
% by Lee, Gan, and Kuo, 2008
% Subband Adaptive Filtering: Theory and Implementation
% Publisher: John Wiley and Sons, Ltd

M = length(S.coeffs);        % Length of FIR filter
mu = S.step;                 % Step size
P = S.order;                 % Projection order
alpha = S.alpha;             % Small constant
AdaptStart = S.AdaptStart;
w = S.coeffs;                % Weight vector of FIR filter
u = zeros(M,1);              % Input signal vector
A = zeros(M,P);              % Projection matrix

ITER = length(un);           % Length of input sequence
yn = zeros(1,ITER);          % Initialize output sequence to zero
en = zeros(1,ITER);          % Initialize error sequence to zero

if isfield(S,'unknownsys')
    b = S.unknownsys;
    norm_b = norm(b);
    eml = zeros(1,ITER);
    ComputeEML = 1;
else
    ComputeEML = 0;
end

for n = 1:ITER
    u = [un(n); u(1:end-1)]; % Input signal vector [u(n),u(n-1),...,u(n-M+1)]'
    A = [u A(:,1:end-1)];    % Projection matrix
    yn(n) = u'*w;            % Output of adaptive filter
    en(n) = dn(n) - yn(n);
    if ComputeEML == 1;
        eml(n) = norm(b-w)/norm_b;        % System error norm (normalized)
    end
    if (mod(n-1,P)==0)&(n >= AdaptStart)
        e = en(n:-1:n-P+1)';
        w = w + mu*A*inv(A'*A + alpha)*e; % Tap-weight adaptation for the next cycle
    end
end

S.coeffs = w;                             % Coefficient values at the final iteration
if ComputeEML == 1;
    S.eml = eml;
end

author: CharlesThaCat

repo: acoustic-interference-cancellation


text:
%% Example of clustering with RSS, AIC and BIC
clear all; close all;

%% Generate Random Clusters from GMM
clear all; close all;
num_samples         = 400;
num_classes         = 4;
dim                 = 2;
[X,labels,gmm]      = ml_clusters_data(num_samples,dim,num_classes);
options.title       = '2D Mixture of Gaussians Dataset';
options.labels      = labels;
[N,D]               = size(X);

if exist('h1','var') && isvalid(h1), delete(h1);end
h1 = ml_plot_data(X,options);
axis equal;

%% Generate Perfect circle
clear all; clc;
nb_samples  = 100;
phi         = linspace(0,2 * pi,nb_samples);
r           = 10;
X           = [r .* cos(phi(:)), r .* sin(phi(:))];
% Plot circle points

plot_options                = [];
plot_options.title          = 'Circle';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1  = ml_plot_data(X,plot_options);
axis square;

%% Find optimal K  for K-means according to metrics
[N,D]                       = size(X);
cluster_options             = [];
cluster_options.method_name = 'kmeans';
repeats                     = 1;
Ks                          = 1:N;
[mus, stds]                 = ml_clustering_optimise(X,Ks,repeats,cluster_options,'Start','plus','Distance','sqeuclidean','MaxIter',500);

RSS_curve = mus(1,:);
AIC_curve = mus(2,:);
BIC_curve = mus(3,:);

%% Plot RSS curve
clc;
if exist('h_rss','var')     && isvalid(h_rss),     delete(h_rss);    end
h_rss = figure;hold on;

plot(RSS_curve,'--', 'LineWidth', 1); hold on;

% Find optimal value on RSS curve
[opt_RSS_vals, opt_Ks_RSS] = ml_curve_opt(RSS_curve,'line');
plot(opt_Ks_RSS, RSS_curve(opt_Ks_RSS) ,'s','Color',[1 0 0]); hold on;
[opt_RSS_vals, opt_Ks_RSS] = ml_curve_opt(RSS_curve,'derivatives');
plot(opt_Ks_RSS, RSS_curve(opt_Ks_RSS) ,'s','Color',[1 0 0]); hold on;
grid on;
box on;
legend('RSS','Location','northeast')
xlabel('K')
title('K-means Evaluation Metrics')

disp(['(RSS) Optimal K: ', num2str(opt_Ks_RSS(1))])

%% Find optimal value on AIC/BIC curve
clc;
if exist('h_rss','var')     && isvalid(h_rss),     delete(h_rss);    end
h_rss = figure;hold on;

plot(RSS_curve,'--', 'LineWidth', 1); hold on;
plot(AIC_curve,'--', 'LineWidth', 1); hold on;
plot(BIC_curve,'--', 'LineWidth', 1); hold on;

% Find optimal value on RSS curve
[opt_RSS_vals, opt_Ks_RSS] = ml_curve_opt(RSS_curve,'derivatives');
plot(opt_Ks_RSS, RSS_curve(opt_Ks_RSS) ,'s','Color',[1 0 0]); hold on;

% Find optimal value on AIC curve
[~, opt_Ks_AIC] = min(AIC_curve);
plot(opt_Ks_AIC, AIC_curve(opt_Ks_AIC) ,'s','Color',[1 0 0]); hold on;

% Find optimal value on BIC curve
[~, opt_Ks_BIC] = min(BIC_curve);
plot(opt_Ks_BIC, BIC_curve(opt_Ks_BIC) ,'s','Color',[1 0 0]); hold on;

grid on;
box on;
legend('RSS', 'AIC', 'BIC','Location','northwest')
xlabel('K')
title('K-means Evaluation Metrics')

disp(['(RSS) Optimal K: ', num2str(opt_Ks_RSS)])
disp(['(AIC) Optimal K: ', num2str(opt_Ks_AIC)])
disp(['(BIC) Optimal K: ', num2str(opt_Ks_BIC)])

%% K-means

cluster_options             = [];
cluster_options.method_name = 'kmeans';
cluster_options.K           = opt_Ks_RSS(1);

result1                     = ml_clustering(X,cluster_options,'Start','plus','Distance','sqeuclidean','MaxIter',500);
result1.lambda              = 2;
[rss,aic,bic]               = ml_clustering_eval(X,result1);

% Plot decision boundary

if exist('hd','var') && isvalid(hd), delete(hd);end
hd = ml_plot_class_boundary(X,result1);
%% Search for optimal \lambda value in AIC(RSS)
clc;
cluster_options             = [];
cluster_options.method_name = 'kmeans';

% Compute Maximum RSS values (when K=1,\alpha)
[N,D]       = size(X);
alpha       = N/4;
repeats     = 1;
K           = 1:alpha;
[mus, stds] = ml_clustering_optimise(X,K,repeats,cluster_options,'Start','plus','Distance','sqeuclidean','MaxIter',500);

% Compute lower/upper bounds of \lambda values according to limits derived
% in slide 9,12,13 of TP2-Metrics-Recap
rss_k1         = mus(1,1);
rss_k2         = mus(1,alpha);
upper_lambda   = (rss_k1 - rss_k2) / (alpha*D);
lower_lambda   = (rss_k1) / (D*N);
lambda_range   = [ceil(lower_lambda):1:ceil(upper_lambda)+1];

% Plot AIC(RSS) with a constrained range of lambda's
clc;
if exist('h_aic','var')     && isvalid(h_aic),     delete(h_aic);    end
h_aic = figure;hold on;
leg_names = [];
aic_mins = [];
for l=1:length(lambda_range)
    cluster_options.lambda    = lambda_range(l); %Magic number
    repeats                   = 1;
    Ks                        = 1:N;
    [mus, stds]               = ml_clustering_optimise(X,Ks,repeats,cluster_options,'Start','plus','Distance','sqeuclidean','MaxIter',500);
    plot(mus(2,:),'-');
    hold on;
    [min_aic, min_k] = min(mus(2,:));
    leg_names = [leg_names; {strcat('\lambda= ',num2str(lambda_range(l)))}];
    aic_mins = [aic_mins; [min_aic, min_k]];
end
legend(leg_names);
grid on;
box on;
xlabel('K');
title('AIC(RSS) for K-means')

if exist('h_bic','var')     && isvalid(h_bic),     delete(h_bic);    end
h_bic = figure;hold on;
plot(mus(3,:),'-');
grid on;
box on;
xlabel('K');
title('BIC(RSS) for K-means')

author: epfl-lasa


text:
function traject = compute_vf_trajectory(points0, vf, options)

% compute_vf_trajectory - compute the trajectories of a vector field
%
%   traject = compute_vf_trajectory(points0, vf, options);
%
%   For each point points0(:,i) in R^2, traject(:,j,i) is the jth point
%   along the trajectory (x(t),y(t)) satisfying
%       d/dt (x,y) = vf(x,y,:)
%
%   This function is very similar to perform_vf_integration.
%
%   The number points along the trajectory is options.niter
%   The timestep is options.dt
%
%   Copyright (c) 2007 Gabriel Peyr?

options.null = 0;
if isfield(options, 'niter')
    niter = options.niter;
else
    niter = 200;
end
if isfield(options, 'dt')
    dt = options.dt;
else
    dt = 0.2;
end

npoints = size(points0,2);
n = size(vf,1);
m = size(vf,2);

traject = zeros(2,niter+1,npoints);
traject(:,1,:) = reshape(points0, [2 1 npoints]);
points = points0;
for i=1:niter
    progressbar(i,niter);
    % evaluate the gradient at the location
    gx = interp2( vf(:,:,1), points(2,:), points(1,:) );
    gy = interp2( vf(:,:,2), points(2,:), points(1,:) );
    points = points + dt * [gx;gy];
    points(1,:) = clamp(points(1,:), 1,n);
    points(2,:) = clamp(points(2,:), 1,m);
    traject(:,i+1,:) = reshape(points, [2 1 npoints]);
end

author: gpeyre

repo: matlab-toolboxes


text:
function all_bbox_proposals = extract_selective_search_boxes_from_dataset(image_db, ss_boxes_dst_file)
%
% This file is part of the code that implements the following ICCV2015 accepted paper:
% title: "Object detection via a multi-region & semantic segmentation-aware CNN model"
% authors: Spyros Gidaris, Nikos Komodakis
% institution: Universite Paris Est, Ecole des Ponts ParisTech
% Technical report: http://arxiv.org/abs/1505.01749
% code: https://github.com/gidariss/mrcnn-object-detection
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2015 Spyros Gidaris
%
% "Object detection via a multi-region & semantic segmentation-aware CNN model"
% Technical report: http://arxiv.org/abs/1505.01749
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

try
    ld = load(ss_boxes_dst_file);
    all_bbox_proposals = ld.all_bbox_proposals;
catch
    chunk_size = 1000;
    num_imgs   = numel(image_db.image_paths);
    num_chunks = ceil(num_imgs/chunk_size);

    ss_boxes_dst_file_in_progress1 = regexprep(ss_boxes_dst_file, '.mat', '_in_progress.mat');
    ss_boxes_dst_file_in_progress2 = regexprep(ss_boxes_dst_file, '.mat', '_in_progress_prev.mat');

    try
        try
            ld = load(ss_boxes_dst_file_in_progress1);
            all_bbox_proposals = ld.all_bbox_proposals;
            first_chunk = ld.chunk + 1;
        catch
            ld = load(ss_boxes_dst_file_in_progress2);
            all_bbox_proposals = ld.all_bbox_proposals;
            first_chunk = ld.chunk + 1;
        end
    catch exception
        fprintf('Exception message %s\n', getReport(exception));
        all_bbox_proposals = cell(num_imgs,1);
        first_chunk = 1;
    end

    total_num_elems = 0;
    total_time = 0;
    count = 0;
    for chunk = first_chunk:num_chunks
        start_idx = (chunk-1) * chunk_size + 1;
        stop_idx  = min(chunk * chunk_size, num_imgs);
        th = tic;
        all_bbox_proposals(start_idx:stop_idx) = extract_selective_search_prposlas(image_db.image_paths(start_idx:stop_idx));
        for i = start_idx:stop_idx
            count = count + 1;
            total_num_elems = total_num_elems + numel(all_bbox_proposals{i});
        end
        elapsed_time = toc(th);
        total_time = total_time + elapsed_time;
        est_rem_time = (total_time / count) * (num_imgs - stop_idx);
        est_num_bytes = (total_num_elems / count) * num_imgs * 4 / (1024*1024*1024);
        fprintf('Extract Selective Search boxes %s %d/%d: ET %.2fmin | ETA %.2fmin | EST. NUM BYTES %.2f giga\n', ...
            image_db.image_set_name, stop_idx, num_imgs, ...
            total_time/60, est_rem_time/60, est_num_bytes);

        if (exist(ss_boxes_dst_file_in_progress1,'file')>0)
            copyfile(ss_boxes_dst_file_in_progress1,ss_boxes_dst_file_in_progress2);
        end

        save(ss_boxes_dst_file_in_progress1, 'all_bbox_proposals', 'chunk', '-v7.3');
    end

    save(ss_boxes_dst_file, 'all_bbox_proposals', '-v7.3');
end
end

function all_box_proposals = extract_selective_search_prposlas(image_paths)
fast_mode = true;
num_imgs = length(image_paths);
all_box_proposals = cell(num_imgs,1);
parfor (i = 1:num_imgs)
%     th = tic;
    img = imread(image_paths{i});
    all_box_proposals{i} = selective_search_boxes(img, fast_mode);
    all_box_proposals{i} = single(all_box_proposals{i}(:,[2 1 4 3]));
%     fprintf(' image %d/%d: elapsed time %.2f\n', i, num_imgs, toc(th))
end
end

author: gidariss

repo: mrcnn-object-detection


text:
function template = PTKCreateTemplateForOriginalImage(original_image, context, reporting)
    % PTKCreateTemplateForOriginalImage. Function for creating a template image
    %     for the original loaded image before any cropping or resizing
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    template = original_image.BlankCopy;
end

author: tomdoel

repo: pulmonarytoolkit


text:
% NUMDIM - estimate a lower bound on the (minimum) number of discrete sources
%            in the data via their second-order statistics.
% Usage:
%   >> num = numdim( data );
%
% Inputs:
%   data   - 2-D data (nchannel x npoints)
%
% Outputs:
%   num    - number of sources (estimated from second order measures)
%
% References:
%   WACKERMANN, J. 1996. Beyond mapping: estimating complexity
%   of multichannel EEG recordings. Acta Neurobiologiae
%   Experimentalis, 56, 197-208.
%   WACKERMANN, J. 1999. Towards a quantitative characterization
%   of functional states of the brain: from non-linear methodology
%   to the global linear description. International Journal of
%   Psychophysiology, 34, 65-80.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 23 January 2003

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

function lambda = numdim( a )

    if nargin < 1
        help numdim;
        return;
    end

% Akaike, Identification toolbox (linear identification)

    a = a';
    b = a'*a/100; % correlation
    [v d] = eig(b);
    %det(d-b); % checking

    l = diag(d);
    l = l/sum(l);
    lambda = real(exp(-sum(l.*log(l))));

    return;

    % testing by duplicating columns
    a = rand(100,5)*2-1;
    a = [a a];
    numdim( a )

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [IdxParent, IdxChildren, Ms] = WaveTreeStructure2D(N, L);
% WaveTreeStructure: parent/children relationships for coefficients
%     of 2-dimensional wavelet transfom
% Invoked after W = @(x) midwt(x,wav);WT = @(x) mdwt(x,wav);
% It finds the relations between the coefficient images obtained by Rice
%     wavelet toolbox
%
% USAGE: [IdxParent, IdxChildren]=WaveTreeStructure2D(C, S)
%
% INPUT:  N: image with size of N x N, scaling/coefficients
%            after 2-d wavelet transform "mdwt"
%         L: denotes total decomposition levels
%         [See instruction of "mdwt" function for definition of C and L]
%
% OUTPUT: IdxParent: N^2 x 1, parent index.
%         IdxParent(i)=j means the parent of the ith coefficient is coefficient j
%
%         IdxChildren: nChildren x (M-Ms(N)), children index.
%         IdxChildren(i,2:end)=[j(1),...,j(nChildren)] means the children of
%         the ith coefficient are coefficients j(1),...,j(nChildren).
%         Since coefficients at level L (finest level) do not have any children,
%         they are not stored.
%
%         The order of coefficients in IdxParent and IdxChildren are the
%         same as the coefficient order in C(:)
%
%         Ms: 1 x N, number of wavelet coefficients for each decomposition
%             level (scaling coeff. are excluded)

%--------------------------------------------------------------------------
% Junzhou Huang, CSE, University of Texas at Arlington, Feb, 2012
%--------------------------------------------------------------------------
for i=1:L
    Ms(i,1)=3*(N/2^(L-i+1))^2;
end
%---------------------------------------------------------------------------%
% INITIALIZATIONS
IdxParent=zeros(N^2,1);
IdxChildren=zeros(N^2/4,5); %Col 1: index i; col 2-5 the 4 children of index i

IJ=zeros(N^2, 2);
%    col 1: R  (row pointer)
%    col 2: C  (column pointer)
[Isub, Jsub]=ind2sub([N, N], [1:N^2]');
IJ(:,1)=Isub;IJ(:,2)=Jsub;

%%%% Parent %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=zeros(N, N);
B(1:N/2^(L-1),1:N/2^(L-1)) = ones(N/2^(L-1), N/2^(L-1));
index0=find(B==1);
IdxParent(index0)=0;

index=find(B==0);
[Isub, Jsub]=ind2sub([N, N], index);

indexPar=sub2ind([N, N], ceil(Isub/2), ceil(Jsub/2));
IdxParent(index)=indexPar;

%%%%%%%% child  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=zeros(N, N);
B(1:N/2^(L-1),1:N/2^(L-1)) = ones(N/2^(L-1), N/2^(L-1));
index0=find(B==1);

IdxChildren(1:length(index0), 1)=index0;
IdxChildren(1:length(index0), 2:5)=0;

B(N/2+1:N, :) = ones(N/2, N);
B(:, N/2+1:N) = ones(N, N/2);

index=find(B==0);
[Isub, Jsub]=ind2sub([N, N], index);
indexChild4=sub2ind([N, N], Isub*2, Jsub*2);
indexChild3=sub2ind([N, N], Isub*2, Jsub*2-1);
indexChild2=sub2ind([N, N], Isub*2-1, Jsub*2);
indexChild1=sub2ind([N, N], Isub*2-1, Jsub*2-1);

IdxChildren(length(index0)+1:end,:)=[index, indexChild1, indexChild2, indexChild3, indexChild4];
IdxChildren2=IdxChildren;
[YY,II] = sort(IdxChildren(:,1));

IdxChildren=zeros(size(IdxChildren));
IdxChildren=IdxChildren2(II,:);

author: thomaskuestner