text:
% Compare the speeds of various inference engines on the DBN in Kjaerulff
% "dHugin: A computational system for dynamic time-sliced {B}ayesian networks",
% Intl. J. Forecasting 11:89-111, 1995.
%
% The intra structure is (all arcs point downwards)
%
%  1 -> 2
%   \  /
%     3
%     |
%     4
%    / \
%   5   6
%   \  /
%     7
%     |
%     8
%
% The inter structure is 1->1, 4->4, 8->8

seed = 0;
rand('state', seed);
randn('state', seed);

ss = 8;
intra = zeros(ss);
intra(1,[2 3])=1;
intra(2,3)=1;
intra(3,4)=1;
intra(4,[5 6])=1;
intra([5 6], 7)=1;
intra(7,8)=1;

inter = zeros(ss);
inter(1,1)=1;
inter(4,4)=1;
inter(8,8)=1;

ns = 2*ones(1,ss);
onodes = 2;
bnet = mk_dbn(intra, inter, ns, 'observed', onodes, 'eclass2', (1:ss)+ss);
for i=1:2*ss
  bnet.CPD{i} = tabular_CPD(bnet, i);
end

T = 4;

engine = {};
engine{end+1} = jtree_unrolled_dbn_inf_engine(bnet, T);
engine{end+1} = jtree_dbn_inf_engine(bnet);
engine{end+1} = smoother_engine(jtree_2TBN_inf_engine(bnet));
%engine{end+1} = smoother_engine(hmm_2TBN_inf_engine(bnet)); % observed nodes have children

inf_time = cmp_inference_dbn(bnet, engine, T)
learning_time = cmp_learning_dbn(bnet, engine, T)

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function S = renorm_1st_phi(S,U,Wop,epsilon)
	U0p = U{1};
	U0p.signal{1} = abs(U0p.signal{1});
	S1 = Wop{2}(U0p);
	S1 = S1.signal{1};
	sub_multiplier = 2^(S{2}.meta.resolution(1)/2);
	S{2}.signal = cellfun(@(x)(x./(S1+epsilon*sub_multiplier)),S{2}.signal,'UniformOutput',0);
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [ HypA ] = computeHyperAFromMean(mean, ecc)
%computeHyperAFromMean Summary of this function goes here
%   Detailed explanation goes here

%     M=@(H) ecc*sinh(H) - H - mean;
%     HypA=fzero(M,mean);
    HypA = solveKepler(mean, ecc);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function result = torus_square_5c2 ( func, r1, r2 )

%*****************************************************************************80
%
%% TORUS_SQUARE_5C2 approximates an integral in a "square" torus in 3D.
%
%  Integration region:
%
%    Points (X,Y,Z) such that:
%
%      R1 - R2 <= SQRT ( X**2 + Y**2 ) <= R1 + R2,
%      -R2 <= Z <= R2.
%
%  Discussion:
%
%    A 24 point 5-th degree formula is used, Stroud number TOR3-C2:5-1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 November 2004
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Arthur H Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied
%    function of three variables which is to be integrated,
%    of the form:
%      function value = func ( x, y, z )
%
%    Input, real R1, the primary radius of the torus.
%
%    Input, real R2, one-half the length of a side of the
%    square cross-section.
%
%    Output, real RESULT, the approximate integral of the function.
%
  b1 = 5.0 / 108.0;
  b2 = 4.0 / 108.0;
  quad = 0.0;

  u1 = sqrt ( r1 * r1 + r2 * r2 );

  v = r2 * sqrt ( 0.6 );

  u2 = sqrt ( r1 * r1 - sqrt ( 3.0 ) * r1 * r2 + r2 * r2 );

  u3 = sqrt ( r1 * r1 + sqrt ( 3.0 ) * r1 * r2 + r2 * r2 );

  for i = 1 : 6

    cs = cos ( i * pi / 3.0 );
    sn = sin ( i * pi / 3.0 );

    x = u1 * cs;
    y = u1 * sn;
    z = v;
    quad = quad + b1 * feval ( func, x, y, z );

    z = -v;
    quad = quad + b1 * feval ( func, x, y, z );

    x = u2 * cs;
    y = u2 * sn;
    z = 0.0;
    quad = quad + b2 * feval ( func, x, y, z );

    x = u3 * cs;
    y = u3 * sn;
    z = 0.0;
    quad = quad + b2 * feval ( func, x, y, z );

  end

  volume = torus_square_volume_3d ( r1, r2 );
  result = quad * volume;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function qber = get_wbqber(obj, qbit)

global Ns Nk Nr Vn H hMod hDemod;

V_RF = quantizeRF(obj.V_RF,qbit);
W_RF = quantizeRF(obj.W_RF,qbit);

for k = 1:Nk
    V_equal(:,:,k) = V_RF * obj.V_B(:,:,k);
     V_equal(:,:,k) =  V_equal(:,:,k)/norm( V_equal(:,:,k),'fro');
    W_equal(:,:,k) = W_RF * obj.W_B(:,:,k);
end

data = randi([0 1],Nk*Ns*2,1);

s = reshape(step(hMod,data),Ns, Nk);

u = sqrt(Vn/2).*(randn(Nr,Nk)+1i*randn(Nr,Nk));

r = zeros(Ns,Nk);
for k = 1:Nk
    r(:,k) = W_equal(:,:,k)' * H(:,:,k) * V_equal(:,:,k) * s(:,k) + W_equal(:,:,k)' * u(:,k);
end

de_data = step(hDemod, r(:));
qber = biterr(data,de_data);

author: Zzhaoxingyu

repo: hybrid-beamforming-for-three-scenes


text:
function [R, log_c_binary] = p2p_ht_converse(a, n_vals, eps)
% Compute the converse sum rate (at the symmetrical rate point) from the binary hypothesis testing
% converse for the joint source distribution: [a 1/3*(1-a); 1/3*(1-a) 1/3*(1-a)].
% This is a helper function called by sw_ht_converse to compute initial values for
% the parameters in computing the Slepian-Wolf composite hypothesis testing converse.

% Input:
% a: parameter to define the joint source distribution, should be in range [0.25, 1)
% n_vals: a vector containing the blocklengths
% eps: target error probability
% Output:
% R: converse sum rate, which is equivalent to the point-to-point optimum sum rate computed by p2p_optimum.
% log_c_binary: the initial values of the parameters for the blocklengths specified in n_vals

    % Compute log binomial coefficients
    log_b = binomial_coeff(max(n_vals));

    R = zeros(length(n_vals),1);
    log_c_binary = zeros(length(n_vals),1);

    for i = 1:length(n_vals)
        n = n_vals(i);
        m = 0:n;
        log_Pr = m*log(a)+(n-m)*log(1/3*(1-a));

        % Compute the optimal rate
        k = n;
        log_sum_Pr = log_b{n+1}(k+1)+log_Pr(k+1);
        while exp(log_sum_Pr) < 1 - eps
            k = k - 1;
            log_sum_Pr = my_logsumexp([log_sum_Pr log_b{n+1}(k+1)+(n-k)*log(3)+log_Pr(k+1)]);  % Add up the probabilities of the strings until the total probability >= 1 - eps
        end
        log_c_binary(i) = log_Pr(k+1)/n;
        d = 1 - eps - (exp(log_sum_Pr) - exp(log_b{n+1}(k+1)+(n-k)*log(3)+log_Pr(k+1)));
        log_lambda = log(d) - (log_Pr(k+1) + (n-k)*log(3));
        U1 = my_logsumexp(log_b{n+1}(k+2:n+1)+(n-(k+1:n))*log(3));
        U2 = log_lambda + (n-k)*log(3);
        R(i) = 1/n*my_logsumexp([U1 U2]);
    end
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
% compute recall and viewpoint accuracy
% modified from PASCAL3D+ VDPM/compute_recall_precision_accuracy.m
% added an input param: prediction filename
% changed a few relative path names
function [recall, precision, accuracy, ap, aa] = compute_recall_precision_accuracy_azimuth(cls, vnum_train, vnum_test, prediction_filename, show_curve)

BASE_DIR = fullfile(mfilename('fullpath'),'../');
addpath(fullfile(BASE_DIR, '../'));
global_variables;

% backward compatibility
if nargin < 3
    vnum_test = vnum_train;
end
if nargin < 4
    prediction_filename = sprintf('../../data/pascal3d/VDPM/data/%s_%d_test_flip.mat', cls, vnum_train);
end
if nargin < 5
    show_curve = 1;
end

azimuth_interval = [0 (360/(vnum_test*2)):(360/vnum_test):360-(360/(vnum_test*2))];

% viewpoint annotation path
%path_ann_view = '../Annotations';
path_ann_view = fullfile(g_pascal3d_root_folder, 'Annotations');

% read ids of validation images
addpath(fullfile(g_pascal3d_root_folder, 'VDPM'));
addpath(fullfile(g_pascal3d_root_folder, 'PASCAL/VOCdevkit/VOCcode'));
VOCinit;
%pascal_init;
ids = textread(sprintf(VOCopts.imgsetpath, 'val'), '%s');
M = numel(ids);

% open prediction file
%filename = sprintf('data/%s_%d_test_flip.mat', cls, vnum_train);
filename = prediction_filename;
object = load(filename);
dets_all = object.dets;

% load pre-packaged records (from extract_records.m)
% this step greatly accelerates testing.
%object = load('voc12val_records.mat');
try
    object = load(fullfile(BASE_DIR, 'voc12val_records.mat'));
catch
    run(fullfile(BASE_DIR, 'extract_records.m'));
    object = load(fullfile(BASE_DIR, 'voc12val_records.mat'));
end
voc12val_records = object.voc12val_records;

energy = [];
correct = [];
correct_view = [];
overlap = [];
count = zeros(M,1);
num = zeros(M,1);
num_pr = 0;
for i = 1:M
    fprintf('%s train view %d test view %d: %d/%d\n', cls, vnum_train, vnum_test, i, M);
    % read ground truth bounding box
    rec = voc12val_records{i}; %PASreadrecord(sprintf(VOCopts.annopath, ids{i}));
    clsinds = strmatch(cls, {rec.objects(:).class}, 'exact');
    diff = [rec.objects(clsinds).difficult];
    clsinds(diff == 1) = [];
    n = numel(clsinds);
    bbox = zeros(n, 4);
    for j = 1:n
        bbox(j,:) = rec.objects(clsinds(j)).bbox;
    end
    count(i) = size(bbox, 1);
    det = zeros(count(i), 1);

    % read ground truth viewpoint
    if isempty(clsinds) == 0
        filename = fullfile(path_ann_view, sprintf('%s_pascal/%s.mat', cls, ids{i}));
        object = load(filename);
        record = object.record;
        view_gt = zeros(n, 1);
        for j = 1:n
            if record.objects(clsinds(j)).viewpoint.distance == 0
                azimuth = record.objects(clsinds(j)).viewpoint.azimuth_coarse;
            else
                azimuth = record.objects(clsinds(j)).viewpoint.azimuth;
            end
            view_gt(j) = find_interval(azimuth, azimuth_interval);
        end
    else
        view_gt = [];
    end

    % get predicted bounding box
    dets = dets_all{i};
    num(i) = size(dets, 1);
    % for each predicted bounding box
    for j = 1:num(i)
        num_pr = num_pr + 1;
        energy(num_pr) = dets(j, 6);
        bbox_pr = dets(j, 1:4);
        view_pr = find_interval((dets(j, 5) - 1) * (360 / vnum_train), azimuth_interval);

        % compute box overlap
        if isempty(bbox) == 0
            o = box_overlap(bbox, bbox_pr);
            [maxo, index] = max(o);
            if maxo >= 0.5 && det(index) == 0
                overlap{num_pr} = index;
                correct(num_pr) = 1;
                det(index) = 1;
                % check viewpoint
                if view_pr == view_gt(index)
                    correct_view(num_pr) = 1;
                else
                    correct_view(num_pr) = 0;
                end
            else
                overlap{num_pr} = [];
                correct(num_pr) = 0;
                correct_view(num_pr) = 0;
            end
        else
            overlap{num_pr} = [];
            correct(num_pr) = 0;
            correct_view(num_pr) = 0;
        end
    end
end
overlap = overlap';

[threshold, index] = sort(energy, 'descend');
correct = correct(index);
correct_view = correct_view(index);
n = numel(threshold);
recall = zeros(n,1);
precision = zeros(n,1);
accuracy = zeros(n,1);
num_correct = 0;
num_correct_view = 0;
for i = 1:n
    % compute precision
    num_positive = i;
    num_correct = num_correct + correct(i);
    if num_positive ~= 0
        precision(i) = num_correct / num_positive;
    else
        precision(i) = 0;
    end

    % compute accuracy
    num_correct_view = num_correct_view + correct_view(i);
    if num_correct ~= 0
        accuracy(i) = num_correct_view / num_positive;
    else
        accuracy(i) = 0;
    end

    % compute recall
    recall(i) = num_correct / sum(count);
end

ap = VOCap(recall, precision);
fprintf('AP = %.4f\n', ap);

aa = VOCap(recall, accuracy);
fprintf('AA = %.4f\n', aa);

% draw recall-precision and accuracy curve
if show_curve
figure;
hold on;
plot(recall, precision, 'r', 'LineWidth',3);
plot(recall, accuracy, 'g', 'LineWidth',3);
xlabel('Recall');
ylabel('Precision/Accuracy');
tit = sprintf('Average Precision = %.1f / Average Accuracy = %.1f', 100*ap, 100*aa);
title(tit);
hold off;
end

function ind = find_interval(azimuth, a)

for i = 1:numel(a)
    if azimuth < a(i)
        break;
    end
end
ind = i - 1;
if azimuth > a(end)
    ind = 1;
end

author: ShapeNet

repo: RenderForCNN


text:
function [res, eyec2, cropped, resize_scale] = align_face_WX(img, f5pt, crop_size, ec_mc_y, ec_y)
f5pt = double(f5pt);
ang_tan = (f5pt(1,2)-f5pt(2,2))/(f5pt(1,1)-f5pt(2,1));
ang = atan(ang_tan) / pi * 180;
img_rot = imrotate(img, ang, 'bicubic');
imgh = size(img,1);
imgw = size(img,2);

% eye center
x = (f5pt(1,1)+f5pt(2,1))/2;
y = (f5pt(1,2)+f5pt(2,2))/2;
% x = ffp(1);
% y = ffp(2);

ang = -ang/180*pi;
%{
x0 = x - imgw/2;
y0 = y - imgh/2;
xx = x0*cos(ang) - y0*sin(ang) + size(img_rot,2)/2;
yy = x0*sin(ang) + y0*cos(ang) + size(img_rot,1)/2;
%}
[xx, yy] = transform(x, y, ang, size(img), size(img_rot));
eyec = round([xx yy]);
x = (f5pt(4,1)+f5pt(5,1))/2;
y = (f5pt(4,2)+f5pt(5,2))/2;
[xx, yy] = transform(x, y, ang, size(img), size(img_rot));
mouthc = round([xx yy]);

resize_scale = ec_mc_y/(mouthc(2)-eyec(2));

img_resize = imresize(img_rot, resize_scale);

res = img_resize;
eyec2 = (eyec - [size(img_rot,2)/2 size(img_rot,1)/2]) * resize_scale + [size(img_resize,2)/2 size(img_resize,1)/2];
eyec2 = round(eyec2);
img_crop = zeros(crop_size, crop_size, size(img_rot,3));
% crop_y = eyec2(2) -floor(crop_size*1/3);
crop_y = eyec2(2) - ec_y;
crop_y_end = crop_y + crop_size - 1;
crop_x = eyec2(1)-floor(crop_size/2);
crop_x_end = crop_x + crop_size - 1;

box(1:2) = guard([crop_x crop_x_end], size(img_resize,2));
box(3:4) = guard([crop_y crop_y_end], size(img_resize,1));
img_crop(box(3)-crop_y+1:box(4)-crop_y+1, box(1)-crop_x+1:box(2)-crop_x+1,:) = img_resize(box(3):box(4),box(1):box(2),:);

% img_crop = img_rot(crop_y:crop_y+img_size-1,crop_x:crop_x+img_size-1);
cropped = img_crop;
end

function r = guard(x, N)
x(x<1)=1;
x(x>N)=N;
r = x;
end

function [xx, yy] = transform(x, y, ang, s0, s1)
% x,y position
% ang angle
% s0 size of original image
% s1 size of target image

x0 = x - s0(2)/2;
y0 = y - s0(1)/2;
xx = x0*cos(ang) - y0*sin(ang) + s1(2)/2;
yy = x0*sin(ang) + y0*cos(ang) + s1(1)/2;
end

author: happynear

repo: FaceVerification


text:
function  kern = sdlfmvKernParamInit(kern)

% SDLFMVKERNPARAMINIT SDLFMV kernel initialization
% FORMAT
% DESC
% Initializes the switching dynamical latent force model structure for
% velocities with some initial parameters. The initial parameters are
% passed through an option in kern.
% RETURN kern : the kernel structure with the default parameters placed in.
% ARG kern : the kernel structure which requires initialisation.
%
% SEEALSO : kernCreate, kernParamInit, sdlfmKernParamInit
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

% Actually, this works just as a wrapper function to use the same
% parameters that the SDLFM kernel.

kern = sdlfmKernParamInit(kern);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [layer, para, expTag] = BuildGeneralMaskBF_CHiME4(modelDir, iteration, nCh, poolingType, poolingType2, nPass)
% modelDir =  'SplitMaskBF_LSTM_MVDR_DNN.U8634_mixed_randPair.771-1024-AM0-7_2048-1981.L2_0.LR_3E-2';
% iteration = 4;
% useChannel = 6;
% poolingType = 'Median';

expTag = [modelDir '_itr' num2str(iteration) '_' num2str(nCh) 'ch_' poolingType '_' poolingType2 '_pass' num2str(nPass)];

% switch modelDir
%     case {'SplitMaskBF_LSTM_MVDR_DNN.U8634_mixed_randPair.771-1024-AM0-7_2048-1981.L2_0.LR_3E-2',...    % it4
%             'Split0_2chMaskBF_LSTM_MVDR_DNN.U8634_mixed_randPair.771-1024-AM0-7_2048-1981.L2_0.LR_3E-3',...
%             'MaskBF_LSTM_MVDR_DNN.U8634_mixed_randPair.771-1024-AM0-7_2048-1981.L2_0.LR_3E-3',...
%             'LSTM_Mask1ch.U69953.771-1024-257.L2_3E-4.LR_3E-2'}    % it2
        modelfile = dir(['nnet/' modelDir '/nnet.itr' num2str(iteration) '.*']);
        modelfile = ['nnet/' modelDir '/' modelfile.name];
% end

dnn = load(modelfile);
para = dnn.para;
layer = dnn.layer;
CMN_idx = ReturnLayerIdxByName(layer, 'CMN');
if length(CMN_idx)==1   % this is a MSE trained mask estimator
    para.topology = SetDefaultValue(para.topology, 'nChMask', 1);
    layerCE = genNetworkMaskBF_CE(para.topology);
    scm_idx = GetScmLayer(layerCE);
    layer = [layer(1:end-2) layerCE(scm_idx:end)];
    CMN_idx = ReturnLayerIdxByName(layer, 'CMN');
    layer = layer(1:CMN_idx(end));
    layer = InitMelLayer(layer, para);
else
    layer = layer(1:CMN_idx(end));
end

para.useGPU = 0;                            % do not use GPU at test time
para.IO = RemoveIOStream(para.IO, [2]);     % as we only have one input stream at test time, strip off the second stream settings
para.IO.inputFeature(1) = 1;
para = ParseOptions2(para);
para.topology.freqBin = (0:1/para.topology.fft_len:0.5)*2*pi;
para.topology.nFreqBin = length(para.topology.freqBin);
para.topology.nCh = nCh;

% revise the first 3 layers (STFT) according to the number of channels
layer{1}.dim =  [1 1] * nCh;
layer{2}.dim =  [para.topology.nFreqBin 1] * nCh;
layer{3}.dim =  [1 1] * para.topology.nFreqBin * nCh;
layer{3}.W = layer{3}.W(1) * eye(layer{3}.dim(1));
layer{3}.b = zeros(layer{3}.dim(1),1);
STFT_layer_idx = 3;

% decide how many layers to feed into LSTM
if ~strcmpi(poolingType, 'none')    % if we use pooling, we will process each channels individually.
    ExtractDims_idx = ReturnLayerIdxByName(layer, 'ExtractDims');
    layer{ExtractDims_idx}.name = 'reshape';    % reshape speech into a tensor of nFreqBin x nFrame x nCh
    layer{ExtractDims_idx}.sourceDims = para.topology.nFreqBin;
    layer{ExtractDims_idx}.targetDims = [para.topology.nFreqBin nCh];
    layer{ExtractDims_idx}.dim = [1 nCh] * para.topology.nFreqBin;
    permuteLayer.name = 'permute';
    permuteLayer.permute_order = [1 3 2];
    permuteLayer.prev = -1;
    permuteLayer.dim = [1 1]*layer{ExtractDims_idx}.dim(1);
    layer = [layer(1:ExtractDims_idx) permuteLayer layer(ExtractDims_idx+1:end)];
end

% add pooling layers
if ~strcmpi(poolingType, 'none')
    lstm_idx =  ReturnLayerIdxByName(layer, 'LSTM');
    pool_layer.name = lower(poolingType);
    pool_layer.prev = -1;
    pool_layer.pool_idx = 3;     % pool over the third dimension, i.e. overal all channels
    pool_layer.dim = [1 1] * para.topology.nFreqBin;

    [scm_layer, split] = GetScmLayer(layer);
    if ~split
        layer2 = [layer(1:lstm_idx+2) pool_layer layer(lstm_idx+3:end)];
        scm_layer = ReturnLayerIdxByName(layer, 'SpatialCovMask');
        layer2{scm_layer}.prev = [-1 STFT_layer_idx-scm_layer];
    else
        layer2 = [layer(1:lstm_idx+2) pool_layer layer(lstm_idx+3:lstm_idx+4) pool_layer layer(lstm_idx+5:end)];
        scm_layer = ReturnLayerIdxByName(layer2, 'SpatialCovSplitMask');
        layer2{scm_layer}.prev = [-4 -1 STFT_layer_idx-scm_layer];
        pool_idx = ReturnLayerIdxByName(layer2, poolingType);
        layer2{pool_idx(1)+1}.prev = -4;
        layer2{pool_idx(2)+1}.prev = [-4 -1 STFT_layer_idx-pool_idx(2)-1];
    end
    [layer2, scm_idx, split] = HandleSTFT(layer2, STFT_layer_idx, 1, para);
else
    layer2 = layer;
end

bf_idx = ReturnLayerIdxByName(layer2, 'beamforming');
layer2{bf_idx}.freqBin = para.topology.freqBin;
layer2{bf_idx}.prev = [-1 STFT_layer_idx-bf_idx];

% add multiple passes of mask estimation
if nPass>1
    power_idx = ReturnLayerIdxByName(layer, 'power');
    layerBF_withoutPooling = layer(power_idx(1)+1:power_idx(2));
    power_idx = ReturnLayerIdxByName(layer2, 'power');
    layerBF = layer2(1:power_idx(2));
    layerFE = layer2(power_idx(2)+1:end);

    layerCascaded = layerBF;
    for i=1:nPass-1
        layerCascaded = [layerCascaded layerBF_withoutPooling];
    end
    layerCascaded = [layerCascaded layerFE];

    [layerCascaded, scm_idx, split] = HandleSTFT(layerCascaded, STFT_layer_idx, nPass, para);

    if strcmpi(poolingType2, 'mean')   % pool between the first pass mask and second pass mask
        % find the first and second pass masks
        add_layer.name = 'add';
        add_layer.prev = [];
        add_layer.dim = [1 1] * para.topology.nFreqBin;

        for pass = 2:nPass
            [scm_idx, split] = GetScmLayer(layerCascaded);
            prev_pass_mask_idx = layerCascaded{scm_idx(1)}.prev(1:2)+scm_idx(1);
            curr_pass_mask_idx = layerCascaded{scm_idx(pass)}.prev(1:2)+scm_idx(pass);

            break_idx = curr_pass_mask_idx(2);
            layerCascaded = [layerCascaded(1:break_idx) add_layer add_layer layerCascaded(break_idx+1:end)];
            add_idx = break_idx+[1 2];
            for i=1:2
                layerCascaded{add_idx(i)}.prev = [prev_pass_mask_idx(i) curr_pass_mask_idx(i)] - add_idx(i);
            end
            layerCascaded{break_idx+3}.prev(1) = -2;
        end

        [layerCascaded, scm_idx, split] = HandleSTFT(layerCascaded, STFT_layer_idx, nPass, para);
    end
    layer = layerCascaded;
else
    layer = layer2;
end

para.out_layer_idx = [STFT_layer_idx ReturnLayerIdxByName(layer, 'beamforming')];
[scm_idx, split] = GetScmLayer(layer);
mask_idx = [];
for i=1:length(scm_idx)
    if split
        mask_idx = [mask_idx layer{scm_idx(i)}.prev(1:2)+scm_idx(i)];
    else
        mask_idx = [mask_idx layer{scm_idx(i)}.prev(1)+scm_idx(i)];
    end
end
para.out_layer_idx = [para.out_layer_idx ReturnLayerIdxByName(layer, 'MVDR_spatialCov') mask_idx length(layer)];
layer = FinishLayer(layer);
end

author: singaxiong

repo: SignalGraph


text:
function [tanu, tanv] = calctangent(RDip)

%% Based on calcrads.m, only difference is that RDip is alread
%% with respect to the sphere origin in calctangent.m
%% MODIFIED 13th JAN 2005 MATT BROOKES

x=RDip(1);
y=RDip(2);
z=RDip(3);
r=sqrt(x*x+y*y+z*z);

if (x==0) && (y==0)
  tanu(1)=1.0; tanu(2)=0; tanu(3)=0;
  tanv(1)=0; tanv(2)=1.0; tanv(3)=0;
else
  RZXY= -(r-z)*x*y;
  X2Y2= 1/(x*x+y*y);

  tanu(1)= (z*x*x + r*y*y) * X2Y2/r;
  tanu(2)= RZXY * X2Y2/r;
  tanu(3)= -x/r;

  tanv(1)= RZXY * X2Y2/r;
  tanv(2)= (z*y*y + r*x*x) * X2Y2/r;
  tanv(3)= -y/r;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% [y_norm,coeff] = NormalizationNoise(y,qmfilter)

% This function estimates the noise level and normalizes the signal so as to
% have signal to noise ratio = 1. The output signal is scaled so that the
% median absolute deviation (MAD) of the wavelet coefficients at the
% finest level is 1.
% sigma_estimated = K*MAD; in a normal distribution K = 1.4826
%
% INPUTS:
% y:   	       signal to normalize
% qmfilter:    quadrature mirror filter
%
% OUTPUTS:
% y_norm:       normalized signal
% coeff:       1/sigma_estimated
%
% LOG:
% 10/17/2012 S. Brigadoi

function [y_norm,coeff] = NormalizationNoise(y,qmf)

    c = cconv(y,qmf,length(y)); % circular convolution (final length = length(y))
	y_downsampled = dyaddown(c); % downsample by 2

	medianAbsDev = mad(y_downsampled);

	if medianAbsDev ~= 0
		y_norm =  (1/1.4826).*y./medianAbsDev;
                coeff = 1/(1.4826*medianAbsDev);
	else
		y_norm = y;
                coeff = 1;
	end

author: BUNPC


text:
function err = getH1error3Q1(node,elem,Du,uh,K,option)
%% GETH1ERRORQ1 H1 norm of the approximation error on hex mesh.
%
%  err = getH1error3Q1(node,elem,@Du,uh) computes the H1 norm of the
%  error between the exact solution Du and finite element approximation
%  uh on a hex mesh described by node and elem.
%
%  The input parameter Du is a function handle uh is a column array with
%  length N representing a piecewise linear function on the mesh.
%
% Author: Huayi Wei < huayiwei1984@gmail.com>
% Modified by Long Chen. The original computation is not quite right. The
% diffusion coefficient cannot be computed separately.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Set up
if ~exist('option','var'), option = []; end
[N,  dim] = size(node);
[NT, nv] = size(elem);
if ~isfield('option','quadOrder')
   if exist('K','var') && ~isempty(K) && ~isnumeric(K)
        quadOrder = 3;
   else
        quadOrder = 2; % number of Guass pts in one direction
   end
else
   quadOrder = option.quadorder;
end

%% Mesh size and Volume
hx = node(elem(:,2),1) - node(elem(:,1),1);
hy = node(elem(:,4),2) - node(elem(:,1),2);
hz = node(elem(:,5),3) - node(elem(:,1),3);
volume = abs(hx.*hy.*hz);

%% Compute error
[GaussPt, weight] = quadptscube(quadOrder);
lambdax = GaussPt(:,1);
lambday = GaussPt(:,2);
lambdaz = GaussPt(:,3);

nQuad = size(GaussPt,1);
err = zeros(NT,1);
for p = 1:nQuad
    % Compute gradient of basis.
    if isfield(option,'isoparametric') && option.isoparametric
        [phi, Dphi, J] = hexbasis(node,elem,GaussPt(p,:));
    else
        Dphi = zeros(NT,3,8);
        phi2d = kron([lambday(p) 1-lambday(p)],[lambdax(p); 1-lambdax(p)]);
        phi = transpose(kron([lambdaz(p) 1-lambdaz(p)],phi2d([1 2 4 3])));
        Dx2d = kron([lambday(p) 1-lambday(p)],1./hx*[-1 1]);
        Dphi(:,1,:) = kron([lambdaz(p) 1-lambdaz(p)],Dx2d(:,[1 2 4 3]));
        Dy2d = kron(1./hy*[-1 1],[lambdax(p) 1-lambdax(p)]);
        Dphi(:,2,:) = kron([lambdaz(p) 1-lambdaz(p)],Dy2d(:,[1 2 4 3]));
        Dphi(:,3,:) = kron(1./hz*[-1 1],phi2d([1 2 4 3]));
        J = volume;
    end
%     [phi, Dphi, J] = hexbasis(node,elem, GaussPt(p,:));
    pxyz = zeros(NT, dim);
    for i = 1:dim
        xi = node(:,i);
        pxyz(:,i) = xi(elem)*phi;
    end
    Duh = zeros(NT,dim);
    for i = 1:dim
        Dphi_i(:,1:nv) = Dphi(:,i,:);
        Duh(:,i) = sum(uh(elem).*Dphi_i,2);
    end
    if exist('K','var') && ~isempty(K) && ~isnumeric(K) % K is a function
        err = err + weight(p)*K(pxyz).*sum((Du(pxyz)-Duh).^2,2);
    else
        err = err + weight(p)*sum((Du(pxyz)-Duh).^2,2);
    end
end
if exist('K','var') && ~isempty(K) && isnumeric(K) && size(K,1) == NT
    err = K.*err;    % K is piecewise constant
end
err = J.*err;
err(isnan(err)) = 0; % singular values are excluded
err = sqrt(sum(err));

repo: ifem


text:
function [ a, b ] = p03_lim ( )

%*****************************************************************************80
%
%% P03_LIM returns the integration limits for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 1.0;

  return
end


text:
function P = ComputeImageFactor (img, imgModel)
% This function computes the singleton OCR factor values for a single
% image.
%
% Input:
%   img: The 16x8 matrix of the image
%   imgModel: The provided, trained image model
%
% Output:
%   P: A K-by-1 array of the factor values for each of the K possible
%     character assignments to the given image
%
% Copyright (C) Daphne Koller, Stanford University, 2012

X = img(:);
N = length(X);
K = imgModel.K;

theta = reshape(imgModel.params(1:N*(K-1)), K-1, N);
bias  = reshape(imgModel.params((1+N*(K-1)):end), K-1, 1);

W = [ bsxfun(@plus, theta * X, bias) ; 0 ];
W = bsxfun(@minus, W, max(W));
W = exp(W);

P=bsxfun(@rdivide, W, sum(W));

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function F = vector_fusion_method(img1_or, img2_or, L, unit, isZCA, de_level, norm, is_overlap, stride)
% fusion method - lrr parts and salient parts

[t1,t2] = size(img1_or);
% calculate image patch number, h_num and w_num
h = t1;
w = t2;
s = stride;
if is_overlap == 1
    % resize input image
    h_dis = ceil((t1-unit)/s)*(s) - (t1-unit);
    w_dis = ceil((t2-unit)/s)*(s) - (t2-unit);
    h = h_dis + t1;
    w = w_dis + t2;
    h_num = ceil((h-unit)/s)+1;
    w_num = ceil((w-unit)/s)+1;
end
img1 = zeros(h, w);
img2 = zeros(h, w);
img1(1:t1, 1:t2) = img1_or(:,:);
img2(1:t1, 1:t2) = img2_or(:,:);
if h_dis ~= 0
    img1(t1:t1+h_dis, 1:t2) = img1_or(t1-h_dis:end, 1:end);
    img2(t1:t1+h_dis, 1:t2) = img2_or(t1-h_dis:end, 1:end);
end
if w_dis ~= 0
    img1(:, t2:t2+w_dis) = img1(:, t2-w_dis:t2);
    img2(:, t2:t2+w_dis) = img2(:, t2-w_dis:t2);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% is_overlap = 0;
[I_b1,I_d1_deep, show_matrix1, count1, I_vector1] = vector_deep_latent(img1, de_level, L, unit, is_overlap, stride, w_num);
[I_b2,I_d2_deep, show_matrix2, count2, I_vector2] = vector_deep_latent(img2, de_level, L, unit, is_overlap, stride, w_num);
count_all = count1;
% fuion for base parts
I_bf = fuison_base_parts(I_b1, I_b2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fuion for salient parts
if is_overlap==1
    I_d_temp = zeros(h, w);
    for i=1:de_level
        temp1 = I_d1_deep(:,:,i);
        temp2 = I_d2_deep(:,:,i);
        [temp_vector] = vector_fuison_detail_parts(temp1, temp2, norm, isZCA);
        c1 = 0;
        for ii=1:s:h-unit+1
            c2 = 0;
            c1 = c1+1;
            for jj=1:s:w-unit+1
                c2 = c2+1;
                temp = temp_vector(:, (c1-1)*(w_num)+c2);
                I_d_temp(ii:(ii+unit-1), jj:(jj+unit-1)) = I_d_temp(ii:(ii+unit-1), jj:(jj+unit-1)) + reshape(temp, [unit unit]);
            end
        end
        % average operation for overlapping position
        I_df_temp(:,:,i) = I_d_temp./count_all;
    end
    I_df = sum(I_df_temp,3);
else
    if de_level>1
        temp1 = sum(I_d1_deep,3);
        temp2 = sum(I_d2_deep,3);
    else
        temp1 = I_d1_deep;
        temp2 = I_d2_deep;
    end

    I_df_vector = vector_fuison_detail_parts(temp1, temp2, norm, isZCA);

    [t1,t2] = size(img1);
    I_df = col2im(I_df_vector,[unit, unit],[t1,t2],'distinct');
    I_df(I_df<0) = 0;
end
% figure;imshow(I_df);

% fused image
F = I_bf + I_df;
F = F(1:t1, 1:t2);
% figure;imshow(F);

end

author: thfylsty


text:
function [ w, xy ] = triangle_unit_o03b ( )

%*****************************************************************************80
%
%% TRIANGLE_UNIT_O03B returns a 3 point quadrature rule for the unit triangle.
%
%  Discussion:
%
%    This rule is precise for monomials through degree 2.
%
%    The integration region is:
%
%      0 <= X
%      0 <= Y
%      X + Y <= 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carlos Felippa,
%    A compendium of FEM integration formulas for symbolic work,
%    Engineering Computation,
%    Volume 21, Number 8, 2004, pages 867-890.
%
%  Parameters:
%
%    Output, real W(3), the weights.
%
%    Output, real XY(2,3), the abscissas.
%
  w(1:3,1) = [ ...
    0.33333333333333333333, ...
    0.33333333333333333333, ...
    0.33333333333333333333 ];

  xy(1:2,1:3) = [ ...
    0.0,  0.5; ...
    0.5,  0.0; ...
    0.5,  0.5 ]';

  return
end


text:
clear all; close all;
I=imread('rice.png');
I=im2double(I);
h=[0,1,0; 1, -4, 1; 0, 1, 0];
J=conv2(I, h, 'same');
K=I-J;
figure;
subplot(121);
imshow(I);
subplot(122);
imshow(K);

author: HuangCongQing


text:
function [Torr] = bar2Torr(bar)
% Convert pressure from bar to torr.
% Chad Greene 2012
Torr = bar*750.062;

author: Sable

repo: mcbench-benchmarks


text:
function buzsakilab_database_submit(bz_datasets,database_name)
% Update existing dataset in the Buzsakilab metadata-database
%
% v1.3
%
% INPUTS
% bz_datasets :
% structure with subfields named after the mySQL columns. The subfields has
% to have the right format to make updates to the database.
%
% database_name :
% Name of database to submit to. Default 'buzsakid_submitdataset'
%
% By Peter Petersen
% petersen.peter@gmail.com

if nargin == 1
    database_name = 'buzsakid_submitdataset';
end

fprintf(['\nUpdating dataset(s) in database: ', database_name, ' from server. \n'])

% Make connection to database
bz_mySQL = bz_database_credentials; % loading credentials

conn = database(database_name,bz_mySQL.userid,bz_mySQL.password,'Vendor','MYSQL','Server',bz_mySQL.address,'PortNumber',3306);
setdbprefs('NullNumberWrite','NaN');
setdbprefs('DataReturnFormat','table');
curs = exec(conn,['SELECT * FROM ' database_name '.' bz_mySQL.table]);
curs = fetch(curs);

if ~strcmp(curs.Message,'Invalid connection.')
    bz_database = table2struct(curs.Data);
    % attributes = attr(curs);
    close(curs);

    for j = 1:length(bz_datasets)
        bz_dataset = bz_datasets(j);
        id = find(bz_dataset.Id==[bz_database.Id]);
        if ~isempty(id) && strcmp(bz_dataset.Session,bz_database(id).Session) & strcmp(bz_dataset.Animal,bz_database(id).Animal) & strcmp(bz_dataset.Investigator,bz_database(id).Investigator)

            whereclause = ['WHERE Id = ' num2str(bz_dataset.Id)];
            disp(['Updating dataset ' num2str(bz_dataset.Id) ': ' bz_dataset.Investigator, ', ', bz_dataset.Animal,', ', bz_dataset.Session])
            bz_dataset = rmfield(bz_dataset, 'Id');

            if any(structfun(@isempty, bz_dataset))
                idx = find(structfun(@isempty, bz_dataset));
                names = fieldnames(bz_dataset);
                for field = 1:length(idx)
                    % warning([names{idx(field)} ' was empty, filling with null.'])
                    bz_dataset = setfield(bz_dataset,names{idx(field)},'null');
                end
            end

            % Removing fields with NaN and null values
            data_table = struct2table(bz_dataset);
            temp = fieldnames(bz_dataset);
            temp2 = find(ismissing(data_table,{NaN,'null'}));
            bz_dataset = rmfield(bz_dataset,temp(temp2));
            % Data to insert
            colnames = fieldnames(bz_dataset)'; %
            data_table = struct2table(bz_dataset);

            % data_table = buzsakilab_database_verify_fieldtypes_and_dataformat(data_table,attributes)
            update(conn,bz_mySQL.table,colnames,data_table, whereclause);
        else
            warning(['Dataset with id=' num2str(bz_dataset.Id) ' not consistent with Investigator, Animal and Sessions in database with same id'])
        end
    end
else
    disp('Connection failed. mySQL server not available')
end
% Close connection to database
close(conn)
% Clear variables
clear conn curs

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function ans = multiple(x, m)

  % multiple : checks whether the first number is
  % divisible through the second number.
  % INPUTS: The inputs must be integers.
  % OUTPUT: Is a boolean value 0 or 1.
  % THROWS: Throws a assertion error if the inputs invalid.

  assert((x == round(x) && m == round(m)), "The inputs must be integers")

  if m == 0
    ans = (x == 0);
  else
    ans = (mod(x,m) == 0);
  end

endfunction

author: TheAlgorithms

repo: MATLAB-Octave


text:
 classdef LLR_TrackInitiatorX < TrackInitiatorX
% LLR_TrackInitiatorX Abstract class
%
% Summary of LLR_TrackInitiatorX:
% This is a simple M out of N track initiator
%
% LLR_TrackInitiatorX Properties:
%   TODO!!
%
% LLR_TrackInitiatorX Methods:
%   + LLR_TrackInitiatorX - Constructor method
%
% (+) denotes public properties/methods
%
% March 2018 Lyudmil Vladimirov, University of Liverpool.

    properties
       TagGenerator
       TrackList
       MeasurementList
       AssocWeightsMatrix
       DataAssociator
       InitFilter
       CustomConfirmConditionFcn
       CustomDeleteConditionFcn
    end

    methods
        function this = LLR_TrackInitiatorX(varargin)
        % LLR_TrackInitiatorX Constructor method
        %
        % DESCRIPTION:
        % * LLR_TrackInitiatorX()
        %   returns a LLR_TrackInitiatorX object handle configured with the provided
        %   PHD filter instance PHDFilter and probability of gating ProbOfGating

            if(nargin==0)
                error('Not enough input arguments');
            end

            % First check to see if a structure was received
            if(nargin==1)
                if(isstruct(varargin{1}))
                    config = varargin{1};
                    this.TagGenerator = config.TagGenerator;
                    this.InitFilter = config.InitFilter;
                    this.DataAssociator = config.DataAssociator;
                    if(isfield(config,"CustomConfirmConditionFcn"))
                        this.CustomConfirmConditionFcn = config.CustomConfirmConditionFcn;
                    else
                       this.CustomConfirmConditionFcn=0;
                    end
                    if(isfield(config,"CustomDeleteConditionFcn"))
                        this.CustomDeleteConditionFcn = config.CustomDeleteConditionFcn;
                    else
                       this.CustomDeleteConditionFcn=0;
                    end
                end
                return;
            end

            % Otherwise, fall back to input parser
            parser = inputParser;
            parser.KeepUnmatched = true;
            parser.parse(varargin{:});
            config = parser.Unmatched;
            this.TagGenerator = config.TagGenerator;
            this.InitFilter = config.InitFilter;
            this.DataAssociator = config.DataAssociator;
            if(isfield(config,"CustomConfirmConditionFcn"))
                this.CustomConfirmConditionFcn = config.CustomConfirmConditionFcn;
            else
               this.CustomConfirmConditionFcn=0;
            end
            if(isfield(config,"CustomDeleteConditionFcn"))
                this.CustomDeleteConditionFcn = config.CustomDeleteConditionFcn;
            else
               this.CustomDeleteConditionFcn=0;
            end
        end

        function [ConfirmedTrackList, TentativeTrackList, deletedTracks] = initiateTracks(this, ConfirmedTrackList, MeasurementList, AssocWeightsMatrix, LikelihoodMatrix)
        % INITIATETRACKS Perform track initiation
        %
        % DESCRIPTION:
        % * [TrackList] = initiateTracks(this,TrackList,MeasurementList)
        %   utilises the provided lists of tracks (TrackList) and measurements
        %   (MeasurementList)
        %
        % INPUTS:
        % * TrackList - (1-by-NumTracks) array of TrackX objects
        % * MeasurementList - (NumObsDims-by-NumMeas) matrix of
        %   measurements
        %
        % See also EllipsoidalGaterX/EllipsoidalGaterX

            THD = -8.9; %-5.9;
            %% Confirmed tracks
            numConfirmedTracks = numel(ConfirmedTrackList);
            deletedTracks = {};
            if numConfirmedTracks>0
                for t=1:numConfirmedTracks
                    % Update LPR
                    Pd = this.DataAssociator.DetectionModel.pdf(ConfirmedTrackList{t}.Filter.StatePrediction.Mean);
                    Pg = this.DataAssociator.Gater.GatingProbability;
                    validInds = find(AssocWeightsMatrix(t,1:end));
                    DL = 0;
                    for j=1:numel(validInds)
                        measInd = validInds(j);
                        if measInd == 1
                            DL = DL + AssocWeightsMatrix(t,1)*(1-Pd*Pg);
%                             ConfirmedTrackList{t}.LPR = ConfirmedTrackList{t}.LPR + AssocWeightsMatrix(t,1)*log((1-Pd*Pg));
                        else
                            lambda = this.DataAssociator.ClutterModel.pdf(ConfirmedTrackList{t}.Filter.MeasurementPrediction.Mean);
                            DL = DL + AssocWeightsMatrix(t,measInd)*Pd*LikelihoodMatrix(t,measInd-1)/lambda;
%                             ConfirmedTrackList{t}.LPR = ConfirmedTrackList{t}.LPR + AssocWeightsMatrix(t,measInd)*log(Pd*LikelihoodMatrix(t,measInd-1)/lambda);
                        end
                    end
                    ConfirmedTrackList{t}.LPR = ConfirmedTrackList{t}.LPR + log(DL);
                    % Update detection history
                    validInds = find(AssocWeightsMatrix(t,2:end));
                    validMeas = MeasurementListX(MeasurementList.getMeasurements(validInds));
                    ConfirmedTrackList{t}.Associations(end+1) = validMeas;
                    % Update maximum LPR
                    if (ConfirmedTrackList{t}.LPR>ConfirmedTrackList{t}.LPR_max)
                        ConfirmedTrackList{t}.LPR_max = ConfirmedTrackList{t}.LPR;
                    end
                end

                % Delete Tracks
                del_tracks = 0;
                del_flag = 0;
                for t = 1:numConfirmedTracks
                    % Check condition
                    del_cond = isnan(ConfirmedTrackList{t}.LPR) || ConfirmedTrackList{t}.LPR < ConfirmedTrackList{t}.LPR_max+THD;
                    if isa(this.CustomDeleteConditionFcn, 'function_handle')
                        c_del_cond = this.CustomDeleteConditionFcn(this,ConfirmedTrackList{t});
                        del_cond = del_cond || c_del_cond;
                    end
                    if(del_cond)
                        deletedTracks{end+1} = ConfirmedTrackList{t};
                        ConfirmedTrackList{t} = [];
                        del_tracks = del_tracks + 1;
                        del_flag = 1;
                    end
                end
                if(del_flag)
                    ConfirmedTrackList = ConfirmedTrackList(~cellfun('isempty',ConfirmedTrackList));
                    numConfirmedTracks = numConfirmedTracks - del_tracks;
                end
            end

            %% Tentative tracks

            % Get indices of all measurements where there exist no possible association
            numMeas = MeasurementList.NumMeasurements;
            rhi = zeros(1, numMeas);
            for measInd = 1:numMeas
                rhi_tmp = 1;
                if(AssocWeightsMatrix>-1) % Check if beta exists
                    for trackInd = 1:numConfirmedTracks
                        rhi_tmp = rhi_tmp*(1-AssocWeightsMatrix(trackInd,measInd+1));
                    end
                end
                rhi(measInd) = rhi_tmp;
            end
            UnassocMeasInd = rhi==1; %any(sum(AssocWeightsMatrix,1)==0,1);

            this.DataAssociator.MeasurementList = MeasurementListX(MeasurementList.Measurements(UnassocMeasInd)); % Only use unassociated measurements
            this.DataAssociator.predictTracks();
            this.DataAssociator.associate();
            this.DataAssociator.updateTracks();

            % Update Detection and Miss counters for all tracks
            for t=1:numel(this.DataAssociator.TrackList)
                % Update counts
                Pd = this.DataAssociator.DetectionModel.pdf(this.DataAssociator.TrackList{t}.Filter.StatePrediction.Mean);
                Pg = this.DataAssociator.Gater.GatingProbability;
                validInds = find(this.DataAssociator.AssocWeightsMatrix(t,1:end));
                DL = 0;
                for j=1:numel(validInds)
                    measInd = validInds(j);
                    if measInd == 1
                        DL = DL + this.DataAssociator.AssocWeightsMatrix(t,1)*(1-Pd*Pg);
%                         this.DataAssociator.TrackList{t}.LPR = this.DataAssociator.TrackList{t}.LPR + this.DataAssociator.AssocWeightsMatrix(t,1)*log((1-Pd*Pg));
                    else
                        lambda = this.DataAssociator.ClutterModel.pdf(this.DataAssociator.TrackList{t}.Filter.MeasurementPrediction.Mean);
                        DL = DL + this.DataAssociator.AssocWeightsMatrix(t,measInd)*Pd*this.DataAssociator.LikelihoodMatrix(t,measInd-1)/lambda;
%                         this.DataAssociator.TrackList{t}.LPR = this.DataAssociator.TrackList{t}.LPR + this.DataAssociator.AssocWeightsMatrix(t,measInd)*log(Pd*this.DataAssociator.LikelihoodMatrix(t,measInd-1)/lambda);
                    end
                end
                this.DataAssociator.TrackList{t}.LPR = this.DataAssociator.TrackList{t}.LPR + log(DL);
                % Update detection history
                validInds = find(this.DataAssociator.AssocWeightsMatrix(t,2:end));
                validMeas = MeasurementListX(this.DataAssociator.MeasurementList.getMeasurements(validInds));
                this.DataAssociator.TrackList{t}.Associations(end+1) = validMeas;
                % Update maximum LPR
                if (this.DataAssociator.TrackList{t}.LPR>this.DataAssociator.TrackList{t}.LPR_max)
                    this.DataAssociator.TrackList{t}.LPR_max = this.DataAssociator.TrackList{t}.LPR;
                end
            end

            % Delete or Confirm Tentative Tracks
            del_tracks = 0;
            del_flag = 0;
            for t = 1:numel(this.DataAssociator.TrackList)
                P_TM = 0.1; %0.1;       % Track miss probability (from Blackman & Popoli)
                P_FC = 10^-4; %10^-6;    % False confirm probability (from Blackman & Popoli)

%                 % High clutter region
%                 xk = this.DataAssociator.TrackList{t}.State.Mean;
%                 x = [-469.710602736631,-2572.92295458513,-2516.31234127506,-1320.74163107899,-835.307363807233, -469.710602736631];
%                 y = [-4915.41679013513,-4065.50122858976,-2858.60918927224,-1380.19032924249,-1261.08679763585, -4915.41679013513];
%                 b = find(inpolygon(xk(1,:),xk(3,:),x,y));
%                 if numel(b)>0
%                     P_FC = 10^-10;
%                 end
                % High and low thresholds for confirmation and deletion
                %  of tentative tracks
                gamma_low = log(P_TM/(1-P_FC));
                gamma_high = log((1-P_TM)/P_FC);

                P = this.DataAssociator.TrackList{t}.State.Covar;

                % Check condition
                del_cond = isnan(this.DataAssociator.TrackList{t}.LPR) || this.DataAssociator.TrackList{t}.LPR < gamma_low;
                if isa(this.CustomDeleteConditionFcn, 'function_handle')
                    c_del_cond = this.CustomDeleteConditionFcn(this,this.DataAssociator.TrackList{t});
                    del_cond = del_cond || c_del_cond;
                end
                if(del_cond)
                    this.DataAssociator.TrackList{t} = [];
                    del_tracks = del_tracks + 1;
                    del_flag = 1;
                elseif(this.DataAssociator.TrackList{t}.LPR > gamma_high)
%                     disp("Confirming new track");
                    ConfirmedTrackList{end+1} = this.DataAssociator.TrackList{t};
                    this.DataAssociator.TrackList{t} = [];
                    del_tracks = del_tracks + 1;
                    del_flag = 1;
                end
            end
            if(del_flag)
                this.DataAssociator.TrackList = this.DataAssociator.TrackList(~cellfun('isempty',this.DataAssociator.TrackList));
            end

            % Initiate new tracks
            UnassocMeasInd = find(~any(this.DataAssociator.ValidationMatrix,1));
            for i=1:numel(UnassocMeasInd)

                measurement = this.DataAssociator.MeasurementList.Vectors(:,UnassocMeasInd(i));
                if(~isempty(this.DataAssociator.MeasurementList.Tags))
                    tag = this.DataAssociator.MeasurementList.Tags(UnassocMeasInd(i));
                else
                    tag = TagX(i);
                end
                stateMean = this.InitFilter.Model.Measurement.finv(measurement);
                dist = GaussianDistributionX(stateMean,...
                                            this.InitFilter.StatePrior.Covar);
                if(~isempty(MeasurementList.Timestamp))
                    statePrior = GaussianStateX(dist,...
                                                MeasurementList.Timestamp);
                else
                    statePrior = GaussianStateX(dist);
                end
                track_config = struct('Filter', copy(this.InitFilter), ...
                                      'Tag', this.TagGenerator.generate(), ...
                                      'LPR', 1,...
                                      'LPR_max', 1,...
                                      'ExistenceProbability', 0.5,...
                                      'Associations', MeasurementListX(measurement,MeasurementList.Timestamp,tag));
                track_config.Filter.initialise('Model',this.InitFilter.Model,...
                                               'StatePrior', statePrior);
                track_config.Trajectory = track_config.Filter.StatePosterior;
                newTrack = TrackX(track_config);

                this.DataAssociator.TrackList{end+1} = newTrack;
            end
            TentativeTrackList = this.DataAssociator.TrackList;
        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function [similarity indSim] = SSSimSize(a, b, blobStruct)
% function similarity = SSSimSize(a, b, blobStruct)
%
% Calculate size similarity

similarity = (blobStruct.imSize - blobStruct.size(a) - blobStruct.size(b)) ...
           ./ blobStruct.imSize;

indSim = similarity;

author: zouchuhang

repo: LayoutNet


text:
function gmsh_mesh2d_write ( gmsh_filename, m, node_num, node_x, ...
  element_order, element_num, element_node )

%*****************************************************************************80
%
%% GMSH_MESH2D_WRITE writes 2D mesh data as a Gmsh mesh file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 October 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Christophe Geuzaine, Jean-Francois Remacle,
%    Gmsh: a three-dimensional finite element mesh generator with
%    built-in pre- and post-processing facilities,
%    International Journal for Numerical Methods in Engineering,
%    Volume 79, Number 11, pages 1309-1331, 2009.
%
%  Parameters:
%
%    Input, string GMSH_FILENAME, the name of the Gmsh file.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_X(M,NODE_NUM), the node coordinates.
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM), the nodes
%    that make up each element.
%

%
%  Enforce 1-based indexing.
%
  element_node = mesh_base_one ( node_num, element_order, element_num, ...
    element_node );
%
%  Open the file.
%
  gmsh = fopen ( gmsh_filename, 'wt' );

  if ( gmsh < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GMSH_MESH2D_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'GMSH_MESH2D_WRITE - Error!' );
  end
%
%  Write the data.
%
  fprintf ( gmsh, '$MeshFormat\n' );
  fprintf ( gmsh, '  2.2 0 8\n' );
  fprintf ( gmsh, '$EndMeshFormat\n' );

  fprintf ( gmsh, '$Nodes\n' );
  fprintf ( gmsh, '  %d\n', node_num );
  for node = 1 : node_num
    fprintf ( gmsh, '  %d %g %g 0.0\n', node, node_x(1:m,node) );
  end
  fprintf ( gmsh, '$EndNodes\n' );

  if ( element_order == 3 )
    element_type = 2;
  elseif ( element_order == 6 )
    element_type = 9;
  end

  tag_num = 2;
  tag1 = 0;
  fprintf ( gmsh, '$Elements\n' );
  fprintf ( gmsh, '  %d\n', element_num );
  for element = 1 : element_num
    fprintf ( gmsh, '  %d  %d  %d  %d  %d', ...
      element, element_type, tag_num, tag1, element );
    for vertex = 1 : element_order
      fprintf ( gmsh, '  %d', element_node(vertex,element) );
    end
    fprintf ( gmsh, '\n' );
  end
  fprintf ( gmsh, '$EndElements\n' );

  fclose ( gmsh );

  return
end

author: johannesgerer


text:
function A=triLower(A)
%%TRILOWER Zero the upper-triangular potion of the matrix and multiply the
%          diagonal terms by 0.5
%
%INPUTS: A A square matrix.
%
%OUTPUTS:  A The matrix A with the upper-triangular portion zeroed and the
%            diagonal elements divided by 2.
%
%This function plays a role in turning the derivative of the covariance
%matrix into the derivative of the lower-triangular square root of the
%covariance matrix. It is described in Equation (80) in [1].
%
%REFERENCES
%[1] D. F. Crouse, "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numRows=size(A,1);
numCol=size(A,2);

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function varargout = ldivide(varargin)
%.\   Pointwise DISKFUN left array divide.
%   F.\G if G is a DISKFUN and F is a double this returns (1/F)*G
%
%   F.\G if G is a double and F is a DISKFUN this returns G\F, but this
%   does not work if F becomes numerically close to zero.
%
%   F.\G we do not allow F and G to both be DISKFUN objects.
%
%   F.\G is the same as the command ldivide(F,G)

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

[varargout{1:nargout}] = ldivide@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
clear all;
OpenBMI('C:\Users\Administrator\Desktop\BCI_Toolbox\git_OpenBMI\OpenBMI') % Edit the variable BMI if necessary
global BMI;
BMI.EEG_DIR=['C:\Users\Administrator\Desktop\BCI_Toolbox\git_OpenBMI\DemoData'];

%% DATA LOAD MODULE
file=fullfile(BMI.EEG_DIR, '\calibration_motorimageryVPkg');
marker={'1','left';'2','right';'3','foot';'4','rest'};
[EEG.data, EEG.marker, EEG.info]=Load_EEG(file,{'device','brainVision';'marker', marker;'fs', [100]});

field={'x','t','fs','y_dec','y_logic','y_class','class', 'chan'};
CNT=opt_eegStruct({EEG.data, EEG.marker, EEG.info}, field);
CNT=prep_selectClass(CNT,{'class',{'right', 'left'}});

%% PRE-PROCESSING MODULE
CNT=prep_filter(CNT, {'frequency', [7 13]});
SMT=prep_segmentation(CNT, {'interval', [750 3500]});

aCNT=prep_addTrials(CNT, CNT);

aveCNT=prep_average(SMT);

prep_rejectArtifactMAxMin(SMT, 70);

author: PatternRecognition

repo: OpenBMI


text:
function plugin_reactivate(foldername);

    % get plugin path
    % ---------------
    fullpluginfolder = fullfile(fileparts(which('eeglab.m')), 'deactivatedplugins', foldername);
    if ~exist(fullpluginfolder)
        error([ 'Could not find folder ' foldername 'in deactivatedplugins folder' ]);
    end;

    disp([ 'Moving plugin ' foldername ' to plugins folder' ]);
    fulldeactivatedpluginfolder = fullfile(fileparts(which('eeglab.m')), 'plugins');
    try
        movefile(fullpluginfolder, fulldeactivatedpluginfolder);
    catch
        eeglab_error;
    end;


text:
function [Vref, dVref] = opf_vref_fcn(x, mpc, refs, mpopt)
%OPF_VREF_FCN  Evaluates voltage angle reference and their gradients.
%   [Vref, dVref] = OPF_VREF_FCN(X, mpc, ref, MPOPT)
%
%   Computes the voltage angle reference using real and imaginary part of complex voltage for
%   AC optimal power flow. Computes constraint vectors and their gradients.
%
%   Inputs:
%     X : optimization vector
%     MPC : MATPOWER case struct
%     REFS : reference vector
%     MPOPT : MATPOWER options struct
%
%   Outputs:
%     VREF  : vector of voltage angle reference
%     DVREF : (optional) angle reference gradients
%
%   Examples:
%       Vref = opf_vref_fcn(x, mpc, refs, mpopt);
%       [Vref, dVref] = opf_vref_fcn(x, mpc, refs, mpopt);
%
%   See also OPF_VREF_HESS

%   MATPOWER
%   Copyright (c) 2018, Power Systems Engineering Research Center (PSERC)
%   by Baljinnyam Sereeter, Delft University of Technology
%   and Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% define named indices into data matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;

%% unpack data
[Vr, Vi] = deal(x{:});
n = length(refs);
nb = length(Vr);

%% compute voltage angle reference
Vref = angle(Vr(refs) + 1j* Vi(refs)) - mpc.bus(refs, VA) * pi/180;

if nargout > 1
    %% compute partials of voltage angle reference w.r.t Vr and Vi
    Vm2 = Vr(refs).^2 + Vi(refs).^2;
    dVa_dVr = sparse(1:n, refs, -Vi(refs) ./ Vm2, n, nb);
    dVa_dVi = sparse(1:n, refs,  Vr(refs) ./ Vm2, n, nb);
    dVref = [dVa_dVr dVa_dVi];          %% Vref w.r.t Vr, Vi
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%
%  REALIZATION OF BASIC MATRIX OPERATIONS BY USER-SUPPLIED
%  FUNCTIONS 'munu_*'
%
%
%  This demo program shows how the user-supplied functions 'munu_*' work.
%  In this particular case, we consider a generalized dynamical system
%  with symmetric matrices M and N, but we will use non-real shift
%  parameters. For this reason, 'munu_*' is used instead of 'msns_*'.

% -----------------------------------------------------------------------
% Generate test problem
% -----------------------------------------------------------------------
%
% As test example, we use an FEM-semidiscretized problem, which leads to
% a generalized system where M (the mass matrix) and N (the negative
% stiffness matrix) are sparse, symmetric, and definite.

load rail821   % load the matrices M N

disp('Problem dimensions:')

n = size(M,1)   % problem order

t = 3; j = sqrt(-1);            % generate complex matrix X0, which
X0 = randn(n,t)+j*randn(n,t);   % contains much fewer columns than rows

% -----------------------------------------------------------------------
% Preprocessing
% -----------------------------------------------------------------------

[M0,ML,MU,N0,dummy,dummy,prm,iprm] = munu_pre(M,N,[],[]);

                        % munu_pre realizes a preprocessing:
                        % - The columns and rows of M and N are
                        %   simultaneously reordered to reduce the
                        %   bandwidth. The result is M0 and N0. prm and
                        %   iprm are the corresponding permutation and
                        %   inverse permutation.
                        % - A Cholesky factorization of M is computed,
                        %   so that the implicit system matrix A0 is
                        %   A0 = inv(ML)*N0*inv(MU).
                        % - Since we consider only the matrix A0 but not
                        %   a dynamical system, we use [] as 2nd and 3rd
                        %   input parameter. dummy = [] is returned.

figure(1), hold off, clf
spy(M)
title('Before preprocessing: nonzero pattern of M. That of N is the same.')

figure(2), hold off, clf
spy(M0)
title('After preprocessing: nonzero pattern of M_0. That of N_0 is the same.')

disp('Verification (test_1, test_2, ... should be small):')

% -----------------------------------------------------------------------
% Multiplication of matrix A0 with X0
% -----------------------------------------------------------------------

munu_m_i(M0,ML,MU,N0)   % initialization and generation of data needed
                        % for matrix multiplications with A0

Y0 = munu_m('N',X0);   % compute Y0 = A0*X0
T0 = ML\(N0*(MU\X0));
test_1 = norm(Y0-T0,'fro')

% -----------------------------------------------------------------------
% Solution of system of linear equations with A0
% -----------------------------------------------------------------------

munu_l_i;   % initialization for solving systems with A0

Y0 = munu_l('N',X0);   % solve A0*Y0 = X0
test_2 = norm(ML\(N0*(MU\Y0))-X0,'fro')

% -----------------------------------------------------------------------
% Solve shifted systems of linear equations, i.e.
% solve (A0+p(i)*I)*Y0 = X0.
% -----------------------------------------------------------------------

disp('Shift parameters:')
p = [ -1; -2+3*j; -2-3*j ]

munu_s_i(p)   % initialization for solution of shifted systems of linear
              % equations with system matrix A0+p(i)*I  (i = 1,...,3)

Y0 = munu_s('N',X0,1);
test_3 = norm(ML\(N0*(MU\Y0))+p(1)*Y0-X0,'fro')

Y0 = munu_s('N',X0,2);
test_4 = norm(ML\(N0*(MU\Y0))+p(2)*Y0-X0,'fro')

Y0 = munu_s('N',X0,3);
test_5 = norm(ML\(N0*(MU\Y0))+p(3)*Y0-X0,'fro')

% -----------------------------------------------------------------------
% Postprocessing
% -----------------------------------------------------------------------
%
% There is no postprocessing.

% -----------------------------------------------------------------------
% Destroy global data structures (clear "hidden" global variables)
% -----------------------------------------------------------------------

munu_m_d;   % clear global variables initialized by munu_m_i
munu_l_d;   % clear global variables initialized by munu_l_i
munu_s_d(p);   % clear global variables initialized by munu_s_i


text:
clear;
clc
close all;

%% 已知
Vn = [0 0 1];
Qa2n = [0.981 0.010 -0.191 0.011]';
Qb2n = [0.936 -0.259 -0.233 -0.046]';

Va = [0.375 0.021 0.876]';
Vb = [0.465  -0.447 0.712]';

%% 四元数转成矩阵
Ca2n = ch_q2m(Qa2n);
Cb2n = ch_q2m(Qb2n);

%%计算 Ca2b, Qa2b
Ca2b = Cb2n' * Ca2n;
Qa2b = ch_m2q(Ca2b);

%% 转换并显示结果
fprintf("b2c的四元数为:%.3f %.3f %.3f %.3f\n", Qa2b(1), Qa2b(2), Qa2b(3), Qa2b(4));
Vb_ = Ca2b*Va;
fprintf("通过Ca2b 将Va转成Vb: %.3f %.3f %.3f, 约等于: %.3f %.3f %.3f \n", Vb_(1), Vb_(2), Vb_(3), Vb(1), Vb(2), Vb(3));

Qa2b = ch_qmul(ch_qconj(Qb2n), Qa2n);
Vb_ = ch_qmulv(Qa2b, Va);
 fprintf("通过Qa2b 将Va转成Vb: %.3f %.3f %.3f, 约等于: %.3f %.3f %.3f \n", Vb_(1), Vb_(2), Vb_(3), Vb(1), Vb(2), Vb(3));

 Qb2a = ch_qconj(Qa2b);
 Va_ = ch_qmulv(Qb2a, Vb);
 fprintf("通过Qb2a 将Vb转成Va: %.3f %.3f %.3f, 约等于: %.3f %.3f %.3f \n", Va_(1), Va_(2), Va_(3), Va(1), Va(2), Va(3));

author: yandld


text:
function rmReconCovPlots(vw,numOfStimuli,ROIs,method,threshold)

    % open view struct
    if ~exist('vw','var') || isempty(vw)
        vw = initHiddenGray;
    end

    % Define some stuff
    % Number of different stimuli used
    %numOfStimuli = 6;

    % All the ROIs that need to be loaded
    if ~exist('ROIs','var') || isempty(ROIs)
        ROIs = {'V1'};
    end

    if ~exist('method','var') || isempty(method)
        method = 'clippedweightavg';
    end

    if ~exist('threshold','var') || isempty(threshold)
	threshold = 0.15;
    end

    % Load the retinotopy model
    vw = rmSelect(vw);
    vw = rmLoadDefault(vw);

    fprintf('[%s]: Retinotopy model loaded \n',mfilename);

    % Load VE values from GLM fit
    load([viewGet(vw,'homedir') '/Recon/VE/VE.mat']);

    vw.dualVE = VE;

    % Loop through the ROIs
    for roi_index = 1 : numel(ROIs)

        fprintf('[%s]: ROI: %s \n',mfilename,ROIs{roi_index});

        % Load the ROI
        vw = loadROI(vw,ROIs{roi_index},[],[],[],1);

        % Load the t-values of the stimuli and create the cov plots
        for stim_index = 1 : numOfStimuli
            vw = loadParameterMap(vw,[viewGet(vw,'homedir') '/Gray/' viewGet(vw,'datatype') '/B-values-stim' num2str(stim_index) '.mat']);

                if ~exist( [viewGet(vw,'homedir') '/Gray/' viewGet(vw,'datatype') '/methods/' method], 'dir');
                    mkdir([viewGet(vw,'homedir') '/Gray/' viewGet(vw,'datatype') '/methods/' method]);
                end

                [CovPlot fig all_models weight] = rmPlotCoverage(vw,...
                                        'fieldRange',7.5,...            % Define the field of view
                                        'method',method,...             % Use this method for combining pRFs
                                        'newfig',-1,...                 % Just give us the data
                                        'normalizeRange',0,...          % Don't normalize the data
                                        'cothresh',threshold,...        % Coherence threshold
                                        'removeOutliers',true,...            % Remove odd sized pRFs
                                        'weight','parameter map');      % Use parameter map for weighting

                 % Save the results

                 save_dir = [ viewGet(vw,'homedir') '/' viewGet(vw,'type') '/' viewGet(vw,'datatype') '/methods/' method];
                 save_name = [ 'CovPlot_' ROIs{roi_index} '_stim' num2str(stim_index) '-' method '.mat' ];

                 disp(sprintf('[%s]: Saving result for stimulus #%d, to %s',mfilename,stim_index,[save_dir '/' save_name]));

                 save([save_dir '/' save_name],'CovPlot');

                 % Save the t-values for this ROI & stimulus (useful
                 % later on)

                 if ~exist([ viewGet(vw,'homedir') '/Recon/ROItValues/' viewGet(vw,'datatype') '/' ROIs{roi_index}],'dir')
                     mkdir([viewGet(vw,'homedir') '/Recon/ROItValues/' viewGet(vw,'datatype') '/' ROIs{roi_index}]);
                 end

                 save([viewGet(vw,'homedir') '/Recon/ROItValues/' viewGet(vw,'datatype') '/' ROIs{roi_index} '/tvalues-stim' num2str(stim_index) '.mat'],'weight');

                    if ~exist( [viewGet(vw,'homedir') '/Recon/pRF models'], 'dir');
                        mkdir([viewGet(vw,'homedir') '/Recon/pRF models']);
                    end

                     save_dir = [ viewGet(vw,'homedir') '/Recon/pRF models'];

                     save_name = [ 'pRF_models-' ROIs{roi_index} '-stim' num2str(stim_index) '-' method '.mat' ];

                     save([save_dir '/' save_name],'all_models');

                     fprintf('[%s]: Saved pRF models to %s \n',mfilename,[save_dir '/' save_name]);

                 clear CovPlot fig all_models weight;
        end
    end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [x, funVal, ValueL]=mc_sgLeastR(A, y, z, opts)
%
%%
% Function mc_cgLassoLeast:
%      Least Squares Loss for Multi-class (task) Learning
%             via the sparse group Lasso penalty
%
%% Problem
%
%  min  1/2 || A x - y||^2 + lambda_1 * \|x\|_1 + \lambda_2 *\|x\|_{q,1}
%
%  q=2, infty
%
%  When y are of binary values (1 and -1), this problem is the
%  multi-class classificaiton problem, which learns the projection
%  matrix x, by sharing information across k classification tasks (with the
%  composite group Lasso penalty)
%
%  The elementes in y can be any real value. In this case, this
%  problem can be regarded as the multi-task learning problem, sharing the
%  same data matrix.
%
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size mxk)
%  z -        cgLasso penalty [lambda1, lambda2]
%  opts-      Optional inputs (default value: opts=[])
%
%% Output parameters:
%  x-         Solution (of size n x k)
%  funVal-    Function value during iterations
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified 1 December 2009.
%
%%
%% Related papers
%
% [1] Jun Liu and Jieping Ye, Moreau-Yosida Regularization for
%     Grouped Tree Structure Learning, NIPS 2010
%
%% Related functions:
%
%  sll_opts
%
%%

%% Verify and initialize the parameters
%%
if (nargin <3)
    error('\n Inputs: A, y and z should be specified!\n');
end

[m,n]=size(A);

k=size(y,2);

if (size(y,1) ~=m)
    error('\n Check the length of y!\n');
end

if (length(z)~=2)
    error('\n z should contain two nonnegative parameters!\n');
end

lambda1=z(1);
lambda2=z(2);

opts=sll_opts(opts); % run sll_opts to set default values (flags)

%% Detailed initialization
%%

%% Normalization

% Please refer to sll_opts for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly
%     This implicit normalization is suggested for the sparse matrix
%                                    but not for the dense matrix
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

% Initialize q
if (~isfield(opts,'q'))
    opts.q=2;
else
    if (opts.q~=0 && opts.q~=2)
        error('\n opts.q should be either 0 or 2');
        % 0 stands for inf
        % 2 stands for 2
    end
end

%% Starting point initialization

% compute AT y
if (opts.nFlag==0)
    ATy =A'* y;
elseif (opts.nFlag==1)
    ATy= A'* y -  mu' * sum(y, 1);    ATy=ATy./repmat(nu, 1, k);
else
    invNu=y./repmat(nu, 1, k);        ATy=A'*invNu-mu' * sum(invNu, 1);
end

% process the regularization parameter
if (opts.rFlag~=0) % lambda1 and lambda2 are the scaling factor lying in [0,1]
    if ( lambda1<0 || lambda1>1 ||  lambda2<0 || lambda2>1)
        error('\n opts.rFlag=1, and lambda1 and lambda2 should be in [0,1]');
    end

    temp=abs(ATy);
    lambda1_max=max(max(abs(temp)));
    lambda1=lambda1*lambda1_max;

    temp=max(temp-lambda1,0);

    if opts.q==0
        lambda2_max=max( sum(temp,2) );
    else
        lambda2_max=sqrt( max( sum(temp.^2,2) ) );
    end

    lambda2=lambda2*lambda2_max;

    %[lambda1_max,lambda2_max, lambda1,lambda2]
end

% initialize a starting point
if opts.init==2
    x=zeros(n,k);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if ( size(x,1)~=n || size(x,2)~=k )
            error('\n Check the input .x0');
        end
    else
        x=ATy;  % if .x0 is not specified, we use ratio*ATy,
        % where ratio is a positive value
    end
end

% compute Ax= A * x
if (opts.nFlag==0)
    Ax=A* x;
elseif (opts.nFlag==1)
    invNu=x./repmat(nu, 1, k);  mu_invNu=mu * invNu;
    Ax=A*invNu -repmat(mu_invNu, m, 1);
else
    Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./repmat(nu, 1, k);
end

if (opts.init==0) % If .init=0, we set x=ratio*x by "initFactor"
    % Please refer to the function initFactor for detail

    x_norm=0;
    for i=1:n
        x_norm=x_norm+ norm( x(i,:), opts.q );
    end

    if x_norm>=1e-6
        ratio=initFactor(x_norm, Ax, y, lambda,'mcLeastR');
        x=ratio*x;    Ax=ratio*Ax;
    end
end

% restart the program for better efficiency
%  this is a newly added function
if (~isfield(opts,'rStartNum'))
    opts.rStartNum=opts.maxIter;
else
    if (opts.rStartNum<=0)
        opts.rStartNum=opts.maxIter;
    end
end

%% The main program

bFlag=0; % this flag tests whether the gradient step only changes a little

% The Armijo Goldstein line search scheme + accelearted gradient descent
if (opts.lFlag==0)

    L=1;
    % We assume that the maximum eigenvalue of A'A is over 1

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,k);

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute the gradient (g) at s
        As=Ax + beta* (Ax-Axp);

        % compute AT As
        if (opts.nFlag==0)
            ATAs=A'*As;
        elseif (opts.nFlag==1)
            ATAs= A'* As -  mu' * sum(As, 1);  ATAs=ATAs./repmat(nu, 1, k);
        else
            invNu=As./repmat(nu, 1, k);        ATAs=A'*invNu-mu' * sum(invNu, 1);
        end

        % obtain the gradient g
        g=ATAs-ATy;

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the projection
            v=s-g/L;

            % projection
            [x,normx]=epsgLasso(v, n, k, lambda1 / L, lambda2 / L, opts.q);

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute Ax= A * x
            if (opts.nFlag==0)
                Ax=A* x;
            elseif (opts.nFlag==1)
                invNu=x./repmat(nu, 1, k);  mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./repmat(nu, 1, k);
            end

            Av=Ax -As;
            r_sum=norm(v,'fro')^2; l_sum=norm(Av,'fro')^2;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is ||Av||_2^2 <= L * ||v||_2^2
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                % fprintf('\n L=%5.6f',L);
            end
        end

        ValueL(iterStep)=L;

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        xxp=x-xp;   Axy=Ax-y;
        funVal(iterStep)=norm(Axy,'fro')^2/2;

       if (opts.q==2)
           funVal(iterStep)=funVal(iterStep)+...
               lambda1* normx(1) +...
               lambda2* normx(2);
       else
           funVal(iterStep)=funVal(iterStep)+...
               lambda1* sum(abs( x(:) ) ) +...
               lambda2* sum( max(abs(x), [], 2) );
       end

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=norm(xp,'fro');    norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end

        % restart the program every opts.rStartNum
        if (~mod(iterStep, opts.rStartNum))
            alphap=0; alpha=1;
            xp=x; Axp=Ax; xxp=zeros(n,k); L =1;
        end
    end
else
    error('\n The current code is only applicable to the case of opts.lFlag=0');
end

author: thomaskuestner


text:
addpath('./libqp/matlab/');

% Example: Training two-class SVM classifier with L2-regularized bias term.

% load training and testing data
load('riply_dataset','trn_X','trn_y','tst_X','tst_y');

% ensure that the labels are +1/-1
trn_y(find(trn_y~=1)) = -1;
tst_y(find(tst_y~=1)) = -1;

% input arguments
% lambda = 1e-2;        % regularization parameter
lambda = 10^-0;        % regularization parameter
options.verb = 1;     % will display progress info

% prepare training examples for BMRM
data.X0 = 1;        % added constant feature -> biased classification rule
data.X = trn_X;     % inputs
data.y = trn_y;     % labels

% call BMRM solver
[W,stat]= bmrm(data,@risk_svm,lambda,options);

if data.X0 == 0
    W0 = 0;
else
    W0 = W(end);
    W = W(1:end-1);
end

% compute training and testing error
ypred = sign(W'*trn_X + W0);
trn_error = sum(ypred(:)~=trn_y(:))/length(trn_y)

ypred = sign(W'*tst_X + W0);
tst_error = sum(ypred(:)~=tst_y(:))/length(tst_y)

% visualize bounary of the trained liner classification rule
%figure;
%ppatterns(trn_X,trn_y);
%pline(W,W0);

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function merge_Callback(hObject, eventdata, handles)

cd(handles.data.squeakfolder);
[detectionFilename, detectionFilepath] = uigetfile([handles.data.settings.detectionfolder '/*.mat'],'Select Detection File(s) for Merging','MultiSelect', 'on');
if isnumeric(detectionFilename); return; end

[audiofile, audiopath] = uigetfile({
    '*.wav;*.ogg;*.flac;*.UVD;*.au;*.aiff;*.aif;*.aifc;*.mp3;*.m4a;*.mp4' 'Audio File'
    '*.wav' 'WAVE'
    '*.flac' 'FLAC'
    '*.ogg' 'OGG'
    '*.UVD' 'Ultravox File'
    '*.aiff;*.aif', 'AIFF'
    '*.aifc', 'AIFC'
    '*.mp3', 'MP3 (it''s probably a bad idea to record in MP3'
    '*.m4a;*.mp4' 'MPEG-4 AAC'
    }, 'Select Corresponding Audio File',handles.data.settings.audiofolder);

if isnumeric(audiofile); return; end

hc = waitbar(0,'Merging Output Structures');

cd(handles.data.squeakfolder);
detectionFilename = cellstr(detectionFilename);

AllBoxes = [];
AllScores = [];
AllClass = [];
AllAccept = [];

for j = 1:length(detectionFilename)
    Calls = loadCallfile(fullfile(detectionFilepath, detectionFilename{j}),handles);

    AllBoxes = [AllBoxes; Calls.Box];
    AllScores = [AllScores; Calls.Score];
    AllClass = [AllClass; Calls.Type];
    AllAccept = [AllAccept; Calls.Accept];
end

% Audio info
audio_info = audioinfo([audiopath audiofile]);
if audio_info.NumChannels > 1
    warning('Audio file contains more than one channel. Use channel 1...')
end

clear('Calls')
%% Merge overlapping boxes
waitbar(.5,hc,'Writing Output Structure');
Calls = merge_boxes(AllBoxes, AllScores .* AllAccept, AllClass, audio_info, 1, 0, 0);

[FileName, PathName] = uiputfile(fullfile(handles.data.settings.detectionfolder, '*.mat'), 'Save Merged Detections');
waitbar(1/2, hc, 'Saving...');
save(fullfile(PathName, FileName),'Calls','-v7.3');
update_folders(hObject, eventdata, handles);
close(hc);

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak


text:
function vectors = maps2vectors(maps)
%MAPS2VECTORS Converts maps in the output of a cnn to vectors.
%   VECTORS = MAPS2VECTORS(MAPS) coverts the MAPS at the output of a
%   convolutional neural net (CNN) and converts them to vectors for
%   input into a fully-connected neural net (FCN). MAPS are in the form
%   of a cell array each cell of which is an array of class double of
%   size M-by-M-by-NumImages where M is the spatial size of the (square)
%   maps, and numImages is the number of images submitted to the system
%   as one minibatch. The number of cells in the cell array is the
%   number of output maps. VECTORS is an array of size
%   (M*M*NumOutputMaps) x (NumImages). Thus, each column of VECTORS is
%   formed by concatenating in the vertical direction all the output
%   maps corresponding to ONE input image. See Fig. 14.FF in DIPUM3E for
%   further details.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% The number of cells in equals the number of output maps in the cnn.
numOutputMaps = numel(maps);

numImages = size(maps{1},3);

% Spatial (M-by-M) size.
M = size(maps{1},1);

% Construct the vectors by reshaping each feature map into a column.
% This is done for all images at the same time.
vectors = [];
for j = 1:numOutputMaps
   vectors = cat(1,vectors,reshape(maps{j},[M*M,numImages]));
end

author: dipum

repo: dipum-toolbox

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


text:
function value = mono_between_enum ( d, n1, n2 )

%*****************************************************************************80
%
%% MONO_BETWEEN_ENUM enumerates monomials in D dimensions of degrees in a range.
%
%  Discussion:
%
%    For D = 3, we have the following table:
%
%     N2 0  1  2  3  4  5  6   7   8
%   N1 +----------------------------
%    0 | 1  4 10 20 35 56 84 120 165
%    1 | 0  3  9 19 34 55 83 119 164
%    2 | 0  0  6 16 31 52 80 116 161
%    3 | 0  0  0 10 25 46 74 110 155
%    4 | 0  0  0  0 15 36 64 100 145
%    5 | 0  0  0  0  0 21 49  85 130
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer D, the spatial dimension.
%
%    Input, integer N1, N2, the minimum and maximum degrees.
%    0 <= N1 <= N2.
%
%    Output, integer VALUE, the number of monomials in D variables,
%    of total degree between N1 and N2 inclusive.
%
  n1 = max ( n1, 0 );

  if ( n2 < n1 )
    value = 0;
    return
  end

  if ( n1 == 0 )
    value = nchoosek ( n2 + d, n2 );
  else if ( n1 == n2 )
    value = nchoosek ( n2 + d - 1, n2 );
  else
    n0 = n1 - 1;
    value = nchoosek ( n2 + d, n2 ) - nchoosek ( n0 + d, n0 );
  end

  return
end

repo: jburkardt-m


text:
function [V,F] = readOBJ_sequence(prefix)
  % [V,F] = readOBJ_sequence(prefix)
  %
  % Inputs:
  %   prefix  glob style paths (e.g., 'mymesh.00*.obj')
  % Outputs:
  %   V  #files cell array of mesh vertex positions
  %   F  #files cell array of mesh faces indices into corresponding vertex
  %     positions
  %   ** Whenever possible, cell arrays are reduced into arrays: i.e., if all
  %   elements of V have the same matrix size, then a #V by 3 by #files matrix
  %   is output instead (similarly for F).
  %   ** Whenever possible, duplicated information across files is removed:
  %   i.e., if all elements of F are the same for each file, then a single #F by
  %   3 matrix is returned.
  %
  [status,res] = system(sprintf('printf "%%s\n" %s',prefix));
  files = strsplit(res,'\n');
  files = files(~cellfun(@isempty,files));
  V = cell(1,numel(files));
  F = cell(1,numel(files));
  parfor fi = 1:numel(files)
    progressbar(fi,numel(files));
    [V{fi},F{fi}] = load_mesh(files{fi});
  end
  V = convert_to_array_if_possible(V);
  F = convert_to_array_if_possible(F);
  function V = convert_to_array_if_possible(V)
    if numel(V) > 0 && ...
      all(cell2mat(cellfun(@size,V','UniformOutput',false))==size(V{1}),[1 2])
      V = cell2mat(permute(V,[1 3 2]));
      if size(V,3) > 1 &&  ...
        V(1,1,1)==V(1,1,2) && ...
        all(V(:,:,1)==V(:,:,2),[1 2]) && ...
        all(V(:,:,1)==V(:,:,3:end),[1 2 3])
        V = V(:,:,1);
      end
    end
  end
end

author: alecjacobson

repo: gptoolbox


text:
function recf = geodet3 (lat, long, alt)

% convert geodetic coordinates to ecf position vector

% input

%  lat  = geodetic latitude (radians)
%         (+north, -south; -pi/2 <= lat <= +pi/2)
%  long = geodetic longitude (radians)
%  alt  = geodetic altitude (kilometers)

% output

%  recf = ecf position vector (kilometers)

% global constants

%  req  = equatorial radius (kilometers)
%  flat = flattening factor (non-dimensional)

% Orbital Mechanics with Matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global req flat

% eccentricity squared

esqr = 2.0 * flat - flat * flat;

n = req / sqrt(1 - esqr * sin(lat) * sin(lat));

% ecf position vector components

recf(1) = (n + alt) * cos(lat) * cos(long);

recf(2) = (n + alt) * cos(lat) * sin(long);

recf(3) = (n * (1 - esqr) + alt) * sin(lat);

author: Sable


text:
%HOMWARP Warp image by an homography
%
% OUT = HOMWARP(H, IM, OPTIONS) is a warp of the image IM obtained by
% applying the homography H to the coordinates of every input pixel.
%
% [OUT,OFFS] = HOMWARP(H, IM, OPTIONS) as above but OFFS is the offset of the
% warped tile OUT with respect to the origin of IM.
%
% Options::
% 'full'            output image contains all the warped pixels, but its
%                   position with respect to the input image is given by the
%                   second return value OFFS.
% 'extrapval',V     set unmapped pixels to this value (default NaN)
% 'roi',R           output image contains the specified ROI in the input image
% 'scale',S         scale the output by this factor
% 'dimension',D     ensure output image is DxD
% 'size',S          size of output image S=[W,H]
% 'coords',{U,V}    coordinate matrices for im, each same size as im.
%
% Notes::
% - The edges of the resulting output image will in general not be
%   be vertical and horizontal lines.
%
% See also HOMOGRAPHY, ITRIM, INTERP2.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function [w,foffs] = homwarp(H, im, varargin)

    % the result of a warp can have negative pixel coordinates or very
    % large coordinates, well away from the origin.

    opt.full = false;
    opt.extrapval = NaN;
    opt.size = [];
    opt.roi = [];
    opt.scale = 1;
    opt.dimension = [];
    opt.coords = [];

    opt = tb_optparse(opt, varargin);

    [w,h] = isize(im);

    if opt.roi
        % opt.box is specified in standard ROI format
        l = opt.roi(1,1); t = opt.roi(2,1);
        r = opt.roi(1,2); b = opt.roi(2,2);
        box = [l r r l; t t b b];
    else
        % bounding box in the input image is the full extent
        box = [1 w w 1; 1 1 h h];
    end

    % map the box vertices in input image to vertices in output image
    Hbox = homtrans(H, box);

    % determine the extent of the image after warping
    xmin = min(Hbox(1,:)); xmax = max(Hbox(1,:));
    ymin = min(Hbox(2,:)); ymax = max(Hbox(2,:));

    % we want the pixel coordinates to map to positive values, determine the minimum
    if opt.full
        offs = floor([xmin, ymin]);

        % and prepend a translational homography that translates the output image
        H = [1 0 -offs(1); 0 1 -offs(2); 0 0 1] * H;
    end

    sz = round([xmax-xmin+1, ymax-ymin+1]);

    % we can specify the maxmimum dimension of the resulting image
    if ~isempty(opt.dimension)
        s = opt.dimension / max(sz);
        H = diag([s s 1]) * H;

        Hbox = homtrans(H, box);

        % determine the extent
        xmin = min(Hbox(1,:)); xmax = max(Hbox(1,:));
        ymin = min(Hbox(2,:)); ymax = max(Hbox(2,:));

        % we want the pixel coordinates to map to positive values, determine the minimum
        offs = floor([xmin, ymin]);

        % and prepend a translational homography that translates the output image
        H = [1 0 -offs(1); 0 1 -offs(2); 0 0 1] * H;
        sz = round([xmax-xmin+1, ymax-ymin+1]);
    end

    if isempty(opt.coords)
        [Ui,Vi] = imeshgrid(im);
    else
        Ui = opt.coords{1};
        Vi = opt.coords{2};
    end

    % determine size of the output image
    if ~isempty(opt.size)
        [Uo,Vo] = imeshgrid(opt.size);
    else
        if opt.full
            [Uo,Vo] = imeshgrid(sz);
        else
            [Uo,Vo] = imeshgrid(im);
        end
    end

    % warp the coordinates of the output pixels
    UV = homtrans(inv(H), [Uo(:) Vo(:)]');
    U = reshape(UV(1,:), size(Uo));
    V = reshape(UV(2,:), size(Vo));

    % interpolate for each color plane
    for p=1:size(im,3)
        imh(:,:,p) = interp2(Ui, Vi, idouble(im(:,:,p)), U, V, 'linear', opt.extrapval);
    end

    if nargout > 0
        w = imh;
    else
        idisp(imh);
    end

    if nargout > 1 && opt.full
        foffs = offs;
    end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function l=gmmlpdf(x,m,v,w)
%GMMPDF calculated the pdf of a mixture of gaussians p=(x,m,v,w)
%
% Inputs: n data values, k mixtures, p parameters
%
%     X(n,p)   Input data vectors, one per row.
%     M(k,p)   mixture means, one row per mixture.
%     V(k,p)   mixture variances, one row per mixture (singlton dimensions will be replicated as required)
%              or else V(p,p,k) for full mixture covariance matrixes
%     W(k,1)   mixture weights, one per mixture. The weights will be normalized by their sum. [default: all equal]
%
% Outputs: (Note that M, V and W are omitted if L==0)
%
%     L(n,1)   log PDF values

%  Bugs/Suggestions
%     (1) Sort out full covariance maatrices
%     (2) Improve plotting
%     (3) Add an extra arument for plotting control

%      Copyright (C) Mike Brookes 2000-2006
%      Version: $Id: gmmlpdf.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[n,p]=size(x);
l=[];           % in case n=0
if nargout>0 || n>0
    x2=x.^2;            % need x^2 for variance calculation
    k=size(m,1);        % number of mixtures
    if size(m,2)~=p
        error('x and m must have the same number of columns');
    end
    if nargin<4
        if nargin<3
            v=1;
        end
        w=ones(k,1);
    end
    w=w/sum(w);         % normalize the weights
    sv=size(v);
    if length(sv)>2 || k==1 && p>1 && sv(1)==p     % full covariance matrices
        error('full covariance matrices not yet implemented');
    else                            % diagonal (or constant) covariance matrices
        if sv(1)==1
            v=v(ones(k,1),:);
        end
        if sv(2)==1
            v=v(:,ones(1,p));
        end

        % If data size is large then do calculations in chunks

        memsize=voicebox('memsize');
        nb=min(n,max(1,floor(memsize/(8*p*k))));    % chunk size for testing data points
        nl=ceil(n/nb);                  % number of chunks

        im=repmat(1:k,1,nb); im=im(:);

        lpx=zeros(n,1);             % log probability of each data point
        wk=ones(k,1);
        wnb=ones(1,nb);
        vi=v.^(-1);                 % calculate quantities that depend on the variances
        vm=sqrt(prod(vi,2)).*w;
        vi=-0.5*vi;

        % first do partial chunk

        jx=n-(nl-1)*nb;                % size of first chunk
        ii=1:jx;
        kk=repmat(ii,k,1);
        km=repmat(1:k,1,jx);
        py=reshape(sum((x(kk(:),:)-m(km(:),:)).^2.*vi(km(:),:),2),k,jx);
        mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
        px=exp(py-mx(wk,:)).*vm(:,ones(1,jx));  % find normalized probability of each mixture for each datapoint
        lpx(ii)=log(sum(px,1))+mx;
        ix=jx+1;

        for il=2:nl
            jx=jx+nb;        % increment upper limit
            ii=ix:jx;
            kk=repmat(ii,k,1);
            py=reshape(sum((x(kk(:),:)-m(im,:)).^2.*vi(im,:),2),k,nb);
            mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
            px=exp(py-mx(wk,:)).*vm(:,wnb);  % find normalized probability of each mixture for each datapoint
            lpx(ii)=log(sum(px,1))+mx;
            ix=jx+1;
        end
        l=lpx-0.5*p*log(2*pi);   % log of total probability of each data point
    end
end
if nargout==0                        % attempt to plot the result
    if p==1                            % one dimensional data
        plot(x,l);
    end
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function model = kbrParamInit(model)

% KBRPARAMINIT KBR model parameter initialisation.
% FORMAT
% DESC initialises the kernel based regression
%  model structure with some default parameters.
% ARG model : the model structure which requires initialisation.
% RETURN model : the model structure with the default parameters placed in.
%
% SEEALSO : kbrCreate, modelCreate, modelParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2007

% MLTOOLS

model.A = randn(model.numData, model.outputDim)/sqrt(model.numData+1);
model.bias = randn(1, model.outputDim)/sqrt(model.numData+1);

author: SheffieldML

repo: GPmat