text:
function [snr_mean, segsnr_mean]= comp_SNR(cleanFile, enhdFile);
%
%   Segmental Signal-to-Noise Ratio Objective Speech Quality Measure
%
%     This function implements the segmental signal-to-noise ratio
%     as defined in [1, p. 45] (see Equation 2.12).
%
%   Usage:  [SNRovl, SNRseg]=comp_snr(cleanFile.wav, enhancedFile.wav)
%
%         cleanFile.wav - clean input file in .wav format
%         enhancedFile  - enhanced output file in .wav format
%         SNRovl        - overall SNR (dB)
%         SNRseg        - segmental SNR (dB)
%
%     This function returns 2 parameters.  The first item is the
%     overall SNR for the two speech signals.  The second value
%     is the segmental signal-to-noise ratio (1 seg-snr per
%     frame of input).  The segmental SNR is clamped to range
%     between 35dB and -10dB (see suggestions in [2]).
%
%   Example call:  [SNRovl,SNRseg]=comp_SNR('sp04.wav','enhanced.wav')
%
%  References:
%
%     [1] S. R. Quackenbush, T. P. Barnwell, and M. A. Clements,
%	    Objective Measures of Speech Quality.  Prentice Hall
%	    Advanced Reference Series, Englewood Cliffs, NJ, 1988,
%	    ISBN: 0-13-629056-6.
%
%     [2] P. E. Papamichalis, Practical Approaches to Speech
%	    Coding, Prentice-Hall, Englewood Cliffs, NJ, 1987.
%	    ISBN: 0-13-689019-9. (see pages 179-181).
%
%  Authors: Bryan L. Pellom and John H. L. Hansen (July 1998)
%  Modified by: Philipos C. Loizou  (Oct 2006)
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $
%-------------------------------------------------------------------------

if nargin ~=2
    fprintf('USAGE: [snr_mean, segsnr_mean]= comp_SNR(cleanFile, enhdFile) \n');
    return;
end

[data1, Srate1, Nbits1]= wavread(cleanFile);
[data2, Srate2, Nbits2]= wavread(enhdFile);
if (( Srate1~= Srate2) | ( Nbits1~= Nbits2))
    error( 'The two files do not match!\n');
end

len= min( length( data1), length( data2));
data1= data1( 1: len);
data2= data2( 1: len);

[snr_dist, segsnr_dist]= snr( data1, data2,Srate1);

snr_mean= snr_dist;
segsnr_mean= mean( segsnr_dist);

% =========================================================================
function [overall_snr, segmental_snr] = snr(clean_speech, processed_speech,sample_rate)

% ----------------------------------------------------------------------
% Check the length of the clean and processed speech.  Must be the same.
% ----------------------------------------------------------------------

clean_length      = length(clean_speech);
processed_length  = length(processed_speech);

if (clean_length ~= processed_length)
  disp('Error: Both Speech Files must be same length.');
  return
end

% ----------------------------------------------------------------------
% Scale both clean speech and processed speech to have same dynamic
% range.  Also remove DC component from each signal
% ----------------------------------------------------------------------

%clean_speech     = clean_speech     - mean(clean_speech);
%processed_speech = processed_speech - mean(processed_speech);

%processed_speech = processed_speech.*(max(abs(clean_speech))/ max(abs(processed_speech)));

overall_snr = 10* log10( sum(clean_speech.^2)/sum((clean_speech-processed_speech).^2));

% ----------------------------------------------------------------------
% Global Variables
% ----------------------------------------------------------------------

winlength   = round(30*sample_rate/1000); %240;		   % window length in samples for 30-msecs
skiprate    = floor(winlength/4); %60;		   % window skip in samples
MIN_SNR     = -10;		   % minimum SNR in dB
MAX_SNR     =  35;		   % maximum SNR in dB

% ----------------------------------------------------------------------
% For each frame of input speech, calculate the Segmental SNR
% ----------------------------------------------------------------------

num_frames = clean_length/skiprate-(winlength/skiprate); % number of frames
start      = 1;					% starting sample
window     = 0.5*(1 - cos(2*pi*(1:winlength)'/(winlength+1)));

for frame_count = 1: num_frames

   % ----------------------------------------------------------
   % (1) Get the Frames for the test and reference speech.
   %     Multiply by Hanning Window.
   % ----------------------------------------------------------

   clean_frame = clean_speech(start:start+winlength-1);
   processed_frame = processed_speech(start:start+winlength-1);
   clean_frame = clean_frame.*window;
   processed_frame = processed_frame.*window;

   % ----------------------------------------------------------
   % (2) Compute the Segmental SNR
   % ----------------------------------------------------------

   signal_energy = sum(clean_frame.^2);
   noise_energy  = sum((clean_frame-processed_frame).^2);
   segmental_snr(frame_count) = 10*log10(signal_energy/(noise_energy+eps)+eps);
   segmental_snr(frame_count) = max(segmental_snr(frame_count),MIN_SNR);
   segmental_snr(frame_count) = min(segmental_snr(frame_count),MAX_SNR);

   start = start + skiprate;

end

author: vipchengrui

repo: traditional-speech-enhancement


text:
% In this example, we consider three disjoint subspaces.
% Any two of them span the third one.
% We take random points from each subspace
% We apply sparse subspace clustering
% via orthogonal matching pursuit algorithm
% to cluster them.
% We finally verify the labels assigned by the clustering
% algorithm with the true labels and measure the
% quality of clustering using f-score.

clear all; close all; clc; rng('default');

shuffle_indices = false;
perform_ssc_omp = true;

% Ambient space dimensions
D = 50;
% number of subspaces = number of clusters
ns = 3;
% subspace dimension
K = 4;
% Angle between subspaces A-B and B-C.
theta = 30; % in degree
% Let us form the subspaces
[A, B, C] = spx.la.spaces.three_disjoint_spaces_at_angle(K, deg2rad(theta));
% Put them together
X = [A B C];
% Put them to bigger dimension
X = spx.la.spaces.k_dim_to_n_dim(X, D);
% Perform a random orthonormal transformation
O = orth(randn(D));
X = O * X;
% Split them again
A = X(:, 1:K);
B = X(:, K + (1:K));
C = X(:, 2*K + (1:K));
fprintf('Spaces after moving to n dimensions and an orthonormal transformation');
spx.la.spaces.describe_three_spaces(A, B, C);
% Number of points per subspace
Ng = 16*K;
% coefficients for Ng vectors chosen randomly in subspace A
coeffs_a = randn(K,Ng);
% avoid small values
coeffs_a = 2*sign(coeffs_a) + coeffs_a;
% coefficients for Ng vectors chosen randomly in subspace B
coeffs_b = randn(K,Ng);
% avoid small values
coeffs_b = 2*sign(coeffs_b) + coeffs_b;
% coefficients for Ng vectors chosen randomly in subspace C
coeffs_c = randn(K,Ng);
% avoid small values
coeffs_c = 2*sign(coeffs_c) + coeffs_c;
% Actual vectors from the three subspaces
XA = A * coeffs_a;
XB = B * coeffs_b;
XC = C * coeffs_c;
% Prepare the overall set of signals
X = [XA XB XC];
% ground through clustering data
true_labels = [1*ones(Ng,1) ; 2*ones(Ng,1) ; 3*ones(Ng,1)];
% All signals are expected to  have a K-sparse representation
if shuffle_indices
    % We also shuffle the signals
    shuffled_indices = randperm(3*Ng);
    X = X(:, shuffled_indices);
    true_labels = true_labels(shuffled_indices);
end
fprintf('Number of signals: %d\n', size(X, 2));

mf = spx.graphics.Figures();

if perform_ssc_omp
    tstart = tic;
    % Application of Sparse subspace clustering
    solver = spx.cluster.ssc.SSC_OMP(X, K, ns);
    result = solver.solve();
    elapsed_time = toc(tstart);
    fprintf('Sparse subspace clustering time spent: %.2f seconds\n', elapsed_time);

    cluster_labels = result.Labels;
    combined_labels = [true_labels cluster_labels]';
    disp(combined_labels);

    % Time to compare the clustering
    comparer = spx.cluster.ClusterComparison(true_labels, cluster_labels);
    result = comparer.fMeasure();
    fprintf('Sparse subspace clustering results:\n');
    comparer.printF1MeasureResult(result);

    mf.new_figure('SSC-OMP Coefficients');
    colormap gray;
    imagesc(abs(solver.Representation));
end

author: indigits

repo: sparse-plex


text:
function [R, G, B] = Lab2RGB(L, a, b)
% function [R, G, B] = Lab2RGB(L, a, b)
% Lab2RGB takes matrices corresponding to L, a, and b in CIELab space
% and transforms them into RGB.  This transform is based on ITU-R
% Recommendation  BT.709 using the D65 white point reference.
% and the error in transforming RGB -> Lab -> RGB is approximately
% 10^-5.  By Mark Ruzon from C code by Yossi Rubner, 23 September 1997.
% Updated for MATLAB 5 28 January 1998.
% Fixed a bug in conversion back to uint8 9 September 1999.

if (nargin == 1)
  b = L(:,:,3);
  a = L(:,:,2);
  L = L(:,:,1);
end

% Thresholds
T1 = 0.008856;
T2 = 0.206893;

[M, N] = size(L);
s = M * N;
L = reshape(L, 1, s);
a = reshape(a, 1, s);
b = reshape(b, 1, s);

% Compute Y
fY = ((L + 16) / 116) .^ 3;
YT = fY > T1;
fY = (~YT) .* (L / 903.3) + YT .* fY;
Y = fY;

% Alter fY slightly for further calculations
fY = YT .* (fY .^ (1/3)) + (~YT) .* (7.787 .* fY + 16/116);

% Compute X
fX = a / 500 + fY;
XT = fX > T2;
X = (XT .* (fX .^ 3) + (~XT) .* ((fX - 16/116) / 7.787));

% Compute Z
fZ = fY - b / 200;
ZT = fZ > T2;
Z = (ZT .* (fZ .^ 3) + (~ZT) .* ((fZ - 16/116) / 7.787));

X = X * 0.950456;
Z = Z * 1.088754;

MAT = [ 3.240479 -1.537150 -0.498535;
       -0.969256  1.875992  0.041556;
        0.055648 -0.204043  1.057311];

RGB = max(min(MAT * [X; Y; Z], 1), 0);

R = reshape(RGB(1,:), M, N) * 255;
G = reshape(RGB(2,:), M, N) * 255;
B = reshape(RGB(3,:), M, N) * 255;

if ((nargout == 1) | (nargout == 0))
  R = uint8(round(cat(3,R,G,B)));
end

author: Cloud-CV

repo: object-proposals


text:
function [p] = dot(tt1,tt2,chunk_start,chunk_end,do_qr)
%Dot  product of two TT tensors
%   [PR]=DOT(TT1,TT2) -- dot product of two TT-tensors
%
%   [PR]=DOT(TT1,TT2, DO_QR) if DO_QR==true is specified, perform the
%   left-to-right QRs of TT1,TT2
%   before the scalar product. It increases the  accuracy in some cases.
%
% In general, returns a 4D tensor of sizes
% r0(tt1), r0(tt2), rd(tt1), rd(tt2)
% If r0(tt1) = r0(tt2) = 1 it returns a matrix of size rd(tt1) x rd(tt2)
%
%   [PR]=DOT(TT1,TT2,CHUNK_START,CHUNK_END[,DO_QR]) If TT2.d>TT1.d, returns a
%   tt_tensor obtained from TT2 by projecting its chunk from CHUNK_START to
%   CHUNK_END onto TT1. Convenient for extracting slices.
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

if (nargin<3)
    do_qr = false;
end;
if (nargin==3)
    do_qr = chunk_start;
end;
if (nargin>3) % dot is only applied to a chunk of tt2
    if (nargin<5)
        do_qr = false;
    end;
    if (tt2.d>tt1.d)
        % chunk_start and chunk_end refer to tt1
        tt2_chunk1 = chunk(tt2, chunk_start, chunk_end);
        D = dot(tt1, tt2_chunk1, do_qr);
        D = reshape(D, tt2.r(chunk_start), tt2.r(chunk_end+1));
        p = [];
        if (chunk_start>1)
            tt2_chunk1 = chunk(tt2, 1, chunk_start-1);
            p = tt2_chunk1*D;
        end;
        if (chunk_end<tt2.d)
            tt2_chunk2 = chunk(tt2, chunk_end+1, tt2.d);
            if (isempty(p))
                p = D*tt2_chunk2;
            else
                p = tkron(p, tt2_chunk2);
            end;
        end;
    else
        error('chunky dot is defined only if tt2.d>tt1.d');
    end;
    return;
end;

if (do_qr)
    [tt1,rv1]=qr(tt1, 'lr');
    [tt2,rv2]=qr(tt2, 'lr');
end;

d=tt1.d;
r1=tt1.r; r2=tt2.r; ps1=tt1.ps; ps2=tt2.ps;
n=tt1.n;
core1=tt1.core; core2=tt2.core;

%ps is always r1(i-1)xr; but if there is a hanging thing? what to do?
%That means, I define a dot product of two "hanging" tensors as a matrix...
%brr.... And if it is hanging on the right?
%
% p=ones(r1(1),r2(1)); % Over r1(1) and r2(1) there is not summation blin.
% %So, if we just sum over n(1) separatedly and leave a strange QxR1(I)xR2(I)
% %matrix...
% for i=1:d
%   cr1=core1(ps1(i):ps1(i+1)-1);
%   cr2=core2(ps2(i):ps2(i+1)-1);
%   p=reshape(p,[r1(i),r2(i)]);
%   cr2=reshape(cr2,[r2(i),numel(cr2)/r2(i)]);
%   p=p*cr2; %p is Q*r1(i)xn(i)xr2(i+1);
%   cr1=reshape(cr1,[r1(i)*n(i),numel(cr1)/(r1(i)*n(i))]);
%   p=reshape(p,[r1(i)*n(i),numel(p)/(r1(i)*n(i))]);
%   p=cr1'*p;
% end

% If the first indices are not ones
p=eye(r1(1)*r2(1));
p = reshape(p, r1(1)*r2(1)*r1(1), r2(1));

for i=1:d
    cr1=core1(ps1(i):ps1(i+1)-1);
    cr2=core2(ps2(i):ps2(i+1)-1);
    cr2=reshape(cr2,[r2(i), n(i)*r2(i+1)]);

    p = p*cr2; % size r11*r21*r1-, n*r2+
    p = reshape(p,r1(1)*r2(1), r1(i)*n(i), r2(i+1));
    p = permute(p, [1, 3, 2]);
    p = reshape(p, r1(1)*r2(1)*r2(i+1), r1(i)*n(i));

    cr1=reshape(cr1,[r1(i)*n(i), r1(i+1)]);

    p = p*conj(cr1); % size r11*r12*r2+, r1+
    p = reshape(p, r1(1)*r2(1), r2(i+1), r1(i+1));
    p = permute(p, [1, 3, 2]);
    p = reshape(p, r1(1)*r2(1)*r1(i+1), r2(i+1));
end;

if (do_qr)
    r2old = size(rv2, 2);
    r1old = size(rv1,2);
    p = p*rv2;
    p = reshape(p, r1(1)*r2(1), r1(d+1), r2old);
    p = permute(p, [1, 3, 2]);
    p = reshape(p, r1(1)*r2(1)*r2old, r1(d+1));
    p = p*conj(rv1);
    p = reshape(p, r1(1), r2(1), r2old, r1old);
    p = permute(p, [1,2,4,3]);
    if ( r1(1)*r2(1) == 1 ) %Save the rabbits
        p=reshape(p, r1old, r2old);
    end
else
    p = reshape(p, r1(1), r2(1), r1(d+1), r2(d+1));
    if ( r1(1)*r2(1) == 1 ) %Save the rabbits
        p=reshape(p, r1(d+1), r2(d+1));
    end
end;
end

author: oseledets

repo: TT-Toolbox


text:
%% Example 9.18: Hermite expansion of Benes SDE
%
% Copyright:
%   2018 - Simo Särkkä and Arno Solin
%
% License:
%   This software is provided under the MIT License. See the accompanying
%   LICENSE file for details.

%% Plot the Hermite expansion example

    x0 = 1/2;
    trans_dens = @(xx,t) 1./sqrt(2*pi*t).*cosh(xx)./cosh(x0).*exp(-0.5*t).*exp(-1./(2*t).*(xx-x0).^2);

    %
    % % Symbolic solution (requires the Symbolic Math Toolbox)
    % syms x;
    % f = matlabFunction(tanh(x));
    % df = matlabFunction(diff(tanh(x),1));
    % d2f = matlabFunction(diff(tanh(x),2));
    % d3f = matlabFunction(diff(tanh(x),3));
    % d4f = matlabFunction(diff(tanh(x),4));
    % d5f = matlabFunction(diff(tanh(x),5));
    %

    f = @(x) tanh(x);
    df = @(x) 1 - tanh(x)^2;
    d2f = @(x) 2*tanh(x)*(tanh(x)^2 - 1);
    d3f = @(x) -2*(tanh(x)^2 - 1)^2 - 4*tanh(x)^2*(tanh(x)^2 - 1);
    d4f = @(x) 16*tanh(x)*(tanh(x)^2 - 1)^2 + 8*tanh(x)^3*(tanh(x)^2 - 1);
    d5f = @(x) -16*(tanh(x)^2 - 1)^3 - 16*tanh(x)^4*(tanh(x)^2 - 1) - 88*tanh(x)^2*(tanh(x)^2 - 1)^2;

    xx = -15:0.1:15;

    figure(1); clf; hold on

      % Solve at t=2
      t = 2;
      herm_x = herm_exp(f(x0),df(x0),d2f(x0),d3f(x0),d4f(x0),d5f(x0),t,xx,x0);

      % Exact density
      fill(xx,trans_dens(xx,t),1,'FaceColor',[.7 .7 .7],'EdgeColor',[.7 .7 .7])

      % Series approximation
      plot(xx,herm_x,'-k','LineWidth',1)

      ylim([-0.21 0.25]);
      legend('Exact density','Series approximation','Location','south');
      xlabel('$x$');
      ylabel('$p(x)$');

    figure(2); clf; hold on

      % Solve at t=5
      t = 5;
      herm_x = herm_exp(f(x0),df(x0),d2f(x0),d3f(x0),d4f(x0),d5f(x0),t,xx,x0);

      % Exact density
      fill(xx,trans_dens(xx,t),1,'FaceColor',[.7 .7 .7],'EdgeColor',[.7 .7 .7])

      % Series approximation
      plot(xx,herm_x,'-k','LineWidth',1)

      ylim([-0.21 0.25]);
      xlabel('$x$');
      ylabel('$p(x)$');

author: AaltoML

repo: SDE

save_path: github-repos/MATLAB/AaltoML-SDE


text:
function Y = spe(X, no_dims, varargin)
%SPE Perform the Stochastic Proximity Embedding algorithm
%
%   Y = spe(X, no_dims, 'Global')
%   Y = spe(X, no_dims, 'Local', k)
%
% Perform the Stochastic Proximity Embedding algorithm. The SPE algorithm
% can be compared to the MDS algorithm, although it is much more efficient.
% The total number of updates necessary after the random initialization is
% usually less than n^2 (where n is the number of datapoints). X is the set
% of datapoints on which SPE has to be applied, and no_dims the number of
% dimensions in the embedding space. If the method is used with the 'Global'
% switch, a minimization of the traditional MDS raw stress function is
% performed (the default). Execution with the 'Local' switch only retains
% Euclidean distances in the local neighborhood of the datapoints. The size
% of the neighborhood is defined by the variable k (default = 12).
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    if nargin <= 2
        variant = 'Global';
    else
        variant = varargin{1};
    end
    if strcmp(variant, 'Local')
        if nargin > 3, k = varargin{2};
        else k = 12; end
    end
    if ~strcmp(variant, 'Global') && ~strcmp(variant, 'Local')
        error('Unknown parameter.');
    end
    if exist('k', 'var') && ischar(k)
        error('Adaptive neighborhood selection is not yet supported in SPE.');
    end

    % Initialize parameters
    lambda = 1;                                         % initial learning parameter
    s = 100;                                            % number of updates per iteration
    max_iter = 20000 + round(.04 * size(X, 1) ^ 2);     % number of iterations
    tol = 1e-5;                                         % regularlization parameter
    n = size(X, 1);                                     % number of datapoints
    if strcmp(variant, 'Local')
        max_iter = max_iter * 3;
    end

    % Compute proximity matrix in original space
    if strcmp(variant, 'Global')
        R = L2_distance(X', X');
        R = R / max(max(R)) * sqrt(2);
    else
        [R, n_ind] = find_nn(X, k);
    end

    % Initialize datapoints randomly
    Y = rand(n, no_dims);

    % Perform SPE
    for i=1:max_iter
        if rem(i, 10000) == 0
            disp(['Iteration ' num2str(i) ' of ' num2str(max_iter) '...']);
        end

        % Select points that should be updated
        J = randperm(n);
        ind1 = J(1:s);
        if strcmp(variant, 'Global')
            ind2 = J(s+1:2*s);
        else
            ind2 = double(n_ind(ind1,:))';
            J = round(rand(1, size(ind2, 2)) * (k - 1)) + 1;
            J = J + ((0:length(J) - 1) * k);
            ind2 = ind2(J);
        end

        % Compute distances between points in embedded space
        D = sqrt(sum((Y(ind1,:) - Y(ind2,:)) .^ 2, 2));

        % Get corresponding distances in real space
        Rt = R((ind1 - 1) * size(R, 1) + ind2)';

        % Update locations of points
        Y(ind1,:) = Y(ind1,:) + lambda * (1/2) * repmat(((Rt - D) ./ (D + tol)), 1, no_dims) .* (Y(ind1,:) - Y(ind2,:));
        Y(ind2,:) = Y(ind2,:) + lambda * (1/2) * repmat(((Rt - D) ./ (D + tol)), 1, no_dims) .* (Y(ind2,:) - Y(ind1,:));

        % Update learning parameter
        lambda = lambda - (lambda / max_iter);
    end

repo: cluster


text:
function dim = size(ob)
%function dim = size(ob)
%       "size" method for "mridft" class

dim = [size(ob.t,1) size(ob.we,1)];

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% Fig. 5.26  Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%script for Figure 5.26
n=[1 5.4];
d1=conv([1 1 0],[1 7 49]);
d=conv(d1,[1 20]);
pzmap(n,d)
axis([-20 0 -7.5 7.5])
hold on
r=roots([1 7 49]);
plot(r,'*')
title('Fig.5.26 Construction for placing a specific point')
z=0:.1:.9;
 wn=2:2:19;
 sgrid(z, wn)
 hold off

author: Sable

repo: mcbench-benchmarks


text:
% Note that for PCA-SIFT, the descriptor dimension changes and you have to set
% kDescDim = 128 in the colmap-tools source code.

load(fullfile(fileparts(mfilename('fullpath')), '../data/pca-sift.mat'));

pool = gcp('nocreate');
if isempty(pool)
    pool = parpool(maxNumCompThreads());
end

parfor i = 1:num_images
    fprintf('Computing features for %s [%d/%d]', ...
            image_names{i}, i, num_images);

    if exist(keypoint_paths{i}, 'file') ...
            && exist(descriptor_paths{i}, 'file')
        fprintf(' -> skipping, already exist\n');
        continue;
    end

    tic;

    % Read the image for keypoint detection, patch extraction and
    % descriptor computation.
    image = imread(image_paths{i});
    if ismatrix(image)
        image = single(image);
    else
        image = single(rgb2gray(image));
    end

    % Read the pre-computed SIFT keypoints.
    keypoints = read_keypoints(keypoint_paths{i});

    % Compute the descriptors for the detected keypoints.
    if size(keypoints, 1) == 0
        descriptors = zeros(0, 80);
    else
        % Extract the descriptors from the patches.
        [~, descriptors] = vl_covdet(image, 'Frames', keypoints', ...
                                     'Descriptor', 'SIFT');
        % Perform PCA-SIFT projection and extract top 80 principal components.
        descriptors = pca_sift_eigvecs * double(descriptors);
        descriptors = single(data.descriptors(1:80,:))';
    end

    % Make sure that each keypoint has one descriptor.
    assert(size(keypoints, 1) == size(descriptors, 1));

    % Write the descriptors to disk for matching.
    write_descriptors(descriptor_paths{i}, descriptors);

    fprintf(' in %.3fs\n', toc);
end

author: ahojnnes

repo: local-feature-evaluation


text:
function [x_hat,PSNR,estFin,estHist] = VAMP(y,iters,n,M_func,Mt_func,PSNR_func,U_func,Ut_func,d,UtM_func,MtU_func)
%
% [x_hat,PSNR] = VAMP(y,iters,n,M_func,Mt_func,PSNR_func)
%
% This function implements VAMP for sparse recovery
%
% Inputs:
%       y       : the measurements
%       iters   : the number of iterations
%       n       : signal size
%       M_func  : function handle that projects onto A. Or a matrix A.
%       Mt_func : function handle that projects onto A'. Or no entry
%       PSNR_func : optional function handle to evaluate PSNR
%
% Output:s
%       x_hat   : the recovered signal.
%       PSNR    : the PSNR trajectory.

vampOpt = VampSlmOpt3;
vampOpt.nitMax = iters;

if (nargin>4)&&(~isempty(Mt_func)) % function handles
    vampOpt.Ah = Mt_func;
    vampOpt.N = n;
end
if (nargin>5)&&(~isempty(PSNR_func)) % error trajectory
    vampOpt.fxnErr = @(x2) PSNR_func(x2);
end
if (nargin>6) % external eigendecomposition [U,D]=eig(M_func(Mt_func(I)))
    vampOpt.U = @(x) U_func(x); % function handle to U
    vampOpt.Uh = @(x) Ut_func(x); % function handle to U'
    vampOpt.d = d; % eigenvalues d=diag(D)
end
if (nargin>9) % external SVD [U,S,V]=svd(M_func) with U matching eigenvectors above
    vampOpt.UhA = @(x) UtM_func(x); % function handle to U'*M = S*V'
    vampOpt.AhU = @(x) MtU_func(x); % function handle to M'*U = V*S'
end

% prepare for VAMP
alf = 1.0;
debias = false; % automatically debias the soft thresholder?
denoi = SoftThreshDMMEstimIn(alf,'debias',debias);

% run VAMP
if nargout==4
  [x_hat,estFin,estHist] = VampSlmEst3(denoi,y,M_func,vampOpt);
else
  [x_hat,estFin] = VampSlmEst3(denoi,y,M_func,vampOpt);
end
PSNR = estFin.err;

author: ricedsp


text:
function [nv_stage] = nvars(N, nx, ni, nu, np)

% N control interval which each have states, integrator vars,
% controls, parameters, and timesteps.
% Ends with a single state.
nv_stage = N*nx + N*ni + N*nu + N*np + N + nx + np;

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
%  Core QP Solver. Use qpng.m instead.
%
% This routine solves the following optimization problem:
%
% min_x .5x' H x + q' x
% s.t.  Aeq x =  beq
%       Ain x <= bin
%
% note that x0 is a feasible starting point.
%
% (C) Nicolo Giorgetti, 2006.
function [x, lam, k, status]=qpsolng(H, q, Aeq, beq, Ain, bin, x0)

nmax=200; % max number of iterations
tol=1e-6; % tolerance

neq=length(beq); % number of eqs
nin=length(bin); % number of ineqs

x=x0;
n=size(x);
Wact=[];
nact=0; % number of rows in Wact (only active inequalities)

status=1; % problem feasible
lam=zeros(neq+nin,1);
k=[];

for k=1:nmax
   % Construct KKT
   K=H;
   r=-q-H*x;
   if neq > 0
      % Add equality constraints
      A=Aeq;
   end
   if nin > 0
      % Add active inequality constraints
      for j=1:nact
         i=Wact(j);
         if j+neq==1
            A=Ain(i,:);
         else
            A=[A; Ain(i,:)];
         end
      end
   end
   nneq=neq+nact;
   if nneq>0
      K=[K, A'; A, zeros(nneq,nneq)];
      r=[r; zeros(nneq,1)];
   end

   y=K\r; %%%%%%%% Check this: we should use pinv instead. Possible numerical problems
   p=y(1:n);

   if norm(p)<tol
      % check optimality or add to work set
      if nact==0
         x=x+p;
         status=1;
         return % successfully
      else
         lam=y(n+neq+1:n+neq+nact);
         [lmin,arg]=min(lam);
         if lmin >= 0
            x=x+p;
            status=1;
            return; % successfully
         else
            % remove constraints from W
            nact=nact-1;
            for j=arg:nact
               Wact(j)=Wact(j+1);
            end
         end
      end
   else % x is not the minimizer in W
      count=nin+1;
      val(1:count)=1.1;
      for j=1:nin
         if Ain(j,1:n)*p > max(tol^2,1e-15)
            val(j)=(bin(j)-Ain(j,1:n)*x)/(Ain(j,1:n)*p);
         end
      end
      val(count)=1;

      [alpha,ind]=min(val); % this is the allowed step size

      x=x+alpha*p;

      if ind < count % there are blocking constraints, add one to W
         nact=nact+1;
         Wact(nact)=ind;
      end
   end
end

eq_infeas=0;
in_infeas=0;

if neq>0
   eq_infeas = (norm(Aeq*x-beq) > rtol*(1+norm(beq)));
end
if nin>0
   in_infeas = (any(Ain*x-bin > -rtol*(1+norm(bin))));
end

if eq_infeas | in_infeas
   status=2; % Max iterations reached, no solution found.
else
   status=3; % Max iterations reached but a feasible solution found.
end

return;

author: opencobra

repo: cobratoolbox


text:
function dirs = getSubDirs(root)
% getSubDirs - gets recursively all subdirectories of a given directory.

dirs = {root};

% get all file
file = dir(root);

% extract directories
file = file([file.isdir] & ~strncmp({file.name},'.',1) ...
  & ~strncmp({file.name},'private',7));

% search subdirectories
for k=1:length(file)
  dirs = [dirs; getSubDirs(fullfile(root,file(k).name))];
end

author: mtex-toolbox

repo: mtex

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


text:
function optionsFile = getEEGOptionsFile()
% This function will get the location of the EEGLAB options file.
%
% Output parameters:
%  optionsFile               The path to the file storing the
%                            EEGLAB options. An empty string will be
%                            returned if there is no options file found.

% Copyright (C) 2015  Kay Robbins with contributions from Nima
% Bigdely-Shamlo, Christian Kothe, Tim Mullen, Jeremy Cockfield,
% and Cassidy Matousek
%
% 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 should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  US
optionsFile = '';
eeglab_options;
if ~isempty(EEGOPTION_PATH)
    homefolder = EEGOPTION_PATH;
elseif ispc
    if ~exist('evalc') %#ok<EXIST>
        eval('evalc = @(x)(eval(x));');
    end
    homefolder = deblank(evalc('!echo %USERPROFILE%'));
else
    homefolder = '~';
end
location = fullfile(homefolder, 'eeg_options.m');
if exist(location, 'file')
    optionsFile = location;
end
end % getEEGOptionsFile

author: VisLab

repo: EEG-Clean-Tools


text:
function UpdateSliceOMatic(X,hSlico)
% This should get the data you want:

% Copyright 2004-2010 The MathWorks, Inc.

if nargin<2
  hSlico=gcf;
else                %%RAB%%
  figure(hSlico)    %%RAB%%
end

D = getappdata(hSlico,'sliceomatic');
D.data = X;
setappdata(hSlico,'sliceomatic',D)

% Get axes children objects that are surfaces
ch = findobj('Type','Surface');

% If no slice then no need to update data in the image

if isempty(ch)
    return
end

% Loop thorugh slices
for i=1:length(ch)
    mySlice = ch(i);

    % Get the X, Y, Z Data
    XData = get(mySlice,'XData');
    YData = get(mySlice,'YData');
    ZData = get(mySlice,'ZData');

    % Find out which one doesn't change by finding the differene in X,Y,Z values
    diffxyz = [all(diff(XData(:))==0) all(diff(YData(:))==0) all(diff(ZData(:))==0)];
    [nodiff, i] = max(diffxyz);

    % Make updates to appropriate slice
    switch i
        case 1
            slice = round(XData(1));
            set(mySlice,'CData',squeeze(X(:,slice,:)))
        case 2
            slice = round(YData(1));
            set(mySlice,'CData',squeeze(X(slice,:,:)))
        case 3
            slice = round(ZData(1));
            set(mySlice,'CData',squeeze(X(:,:,slice)))
    end

end

% % The subroutine is attached. It takes the 3D image matrix as its input. So
% % if using sliceomatic you would do the following:
% %
% % 1. Read Images
% % 2. Manipulate
% % 3. Call sliceomatic(ImageData)
% % 4. Make some changes (i.e. ImageData(ImageData<100)=0
% % 5. call mySubRoutine(ImageData)
% %
% % Everything should be updated and it runs fast
% %
% % --L--
% %
% % -----Original Message-----
% % From: Robert Bemis
% % Sent: Friday, February 20, 2004 2:26 PM
% % To: Laurens Schalekamp
% % Subject: movtext
% %
% % Please enjoy file: movtext
% % courtesy of rbemis@mathworks.com


text:
function [L_,U_] = glover_sherali_raw(a,a0,u)

L_ = -inf(length(a),1);
U_ = inf(length(a),1);
% Map to sorted
[a_, loc] = sort(a,'descend');
n = length(a_);

%From here, we only work with binary and in sorted
for j = 1:n
    s = a_(1:j);
    SN(j) = sum(s(s>0));
end
%SN = cumsum(max(a_,0));
try
    SN_u_plus_1 = sum(a_(1:u+1));
catch
    1
end
SN_u_minus_1 = sum(a_(1:u-1));

fixed_at_one = find(a_ > SN_u_plus_1-a0);
if ~isempty(fixed_at_one)
    L_(fixed_at_one) = 1;
    U_(fixed_at_one) = 1;
end

jhat = u+1:n;
fixed_at_zero = jhat(find(a_(jhat) < a0-SN_u_minus_1));
if ~isempty(fixed_at_zero)
    U_(fixed_at_zero) = 0;
    L_(fixed_at_zero) = 0;
end

L = L_(loc);
U = U_(loc);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%................................................................

function stiffness=formStiffness2D(GDof,numberElements,...
    elementNodes,numberNodes,nodeCoordinates,C,thickness)

% compute stiffness matrix
% for plane stress Q8 elements

stiffness=zeros(GDof);

% 3 by 3 quadrature
[gaussWeights,gaussLocations]=gauss2d('3x3');

for e=1:numberElements
  numNodePerElement = length(elementNodes(e,:));
  numEDOF = 2*numNodePerElement;
  elementDof=zeros(1,numEDOF);
  for i = 1:numNodePerElement
      elementDof(2*i-1)=2*elementNodes(e,i)-1;
      elementDof(2*i)=2*elementNodes(e,i);
  end

  % cycle for Gauss point
  for q=1:size(gaussWeights,1)
    GaussPoint=gaussLocations(q,:);
    xi=GaussPoint(1);
    eta=GaussPoint(2);

% shape functions and derivatives
    [shapeFunction,naturalDerivatives]=shapeFunctionQ8(xi,eta);

% Jacobian matrix, inverse of Jacobian,
% derivatives w.r.t. x,y
    [Jacob,invJacobian,XYderivatives]=...
        Jacobian(nodeCoordinates(elementNodes(e,:),:),naturalDerivatives);

%  B matrix
    B=zeros(3,numEDOF);
    B(1,1:2:numEDOF)       = XYderivatives(:,1)';
    B(2,2:2:numEDOF)  = XYderivatives(:,2)';
    B(3,1:2:numEDOF)       = XYderivatives(:,2)';
    B(3,2:2:numEDOF)  = XYderivatives(:,1)';

% stiffness matrix
    stiffness(elementDof,elementDof)=...
        stiffness(elementDof,elementDof)+...
        B'*C*thickness*B*gaussWeights(q)*det(Jacob);
  end
end

repo: Aero-matlab


text:
clear, clc;

% This is an example for running the function sgLeastR
%
%  Problem:
%
%  min  1/2 || A x - y||^2 + z_1 \|x\|_1 + z_2 * sum_j w_j ||x_{G_j}||
%
%  G_j's are the non-overlapping groups
%
%    we apply L1 for each element
%    and the L2 for the non-overlapping group
%
%  The group information is contained in
%  opts.ind, which is a 3 x nodes matrix, where nodes denotes the number of
%  nodes of the tree.
%  opts.ind(1,:) contains the starting index
%  opts.ind(2,:) contains the ending index
%  opts.ind(3,:) contains the corresponding weight (w_j)
%
%% Related papers
%
% [1] Jun Liu and Jieping Ye, Moreau-Yosida Regularization for
%     Grouped Tree Structure Learning, NIPS 2010

cd ..
cd ..

root=cd;
addpath(genpath([root '/SLEP']));
                     % add the functions in the folder SLEP to the path

% change to the original folder
cd Examples/sgLasso;

m=50;  n=100;       % The data matrix is of size m x n

%randNum=10;           % a random number

% ---------------------- generate random data ----------------------
%randn('state',(randNum-1)*3+1);
A=randn(m,n);        % the data matrix

%randn('state',(randNum-1)*3+2);
xOrin=full(sprandn(n, 1,1));
xOrin(1:50,:)=0;

%randn('state',(randNum-1)*3+3);
noise=randn(m,1);
y=A*xOrin +...
    noise*0.01;      % the response

%----------------------- Set optional items -----------------------
opts=[];

% Starting point
opts.init=2;        % starting from a zero point

% Termination
opts.tFlag=5;       % run .maxIter iterations
opts.maxIter=100;   % maximum number of iterations

% regularization
opts.rFlag=1;       % use ratio

% Normalization
opts.nFlag=0;       % without normalization

% Group Property
opts.ind=[ [1, 20, sqrt(20)]', [21, 40, sqrt(20)]',...
    [41, 50, sqrt(10)]', [51, 70, sqrt(20)]', [71,100, sqrt(30)]'];

%----------------------- Run the code -----------------------
z=[0.1,0.1];
tic;
[x1, funVal1, ValueL]= sgLeastR(A, y, z, opts);
toc;

a=xOrin~=0;
b=x1~=0;
c=a.*b;

author: zzstefan

repo: BrainNetClass


text:
function mtype = mwrite (filename, A, Z, comments_filename)		    %#ok
%MWRITE write a matrix to a file in Matrix Market form.
%
%   Example:
%        mtype = mwrite (filename, A, Z, comments_filename)
%
% A can be sparse or full.
%
% If present and non-empty, A and Z must have the same dimension.  Z contains
% the explicit zero entries in the matrix (which MATLAB drops).  The entries
% of Z appear as explicit zeros in the output file.  Z is optional.  If it is
% an empty matrix it is ignored.  Z must be sparse or empty, if present.
% It is ignored if A is full.
%
% filename is the name of the output file.  comments_filename is the file
% whose contents are include after the Matrix Market header and before the
% first data line.  Ignored if an empty string or not present.
%
% See also mread.

% Copyright 2006-2007, Timothy A. Davis

error ('mwrite mexFunction not found') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function bernoulli_poly2_test ( )

%*****************************************************************************80
%
%% BERNOULLI_POLY2_TEST tests BERNOULLI_POLY2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  n = 15;
  x = 0.2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BERNOULLI_POLY2_TEST\n' );
  fprintf ( 1, '  BERNOULLI_POLY2 evaluates Bernoulli polynomials.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  X = %f\n', x );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I          BX\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    bx = bernoulli_poly2 ( i, x );
    fprintf ( 1, '  %2d  %14f\n', i, bx )
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function images = findAllImagesInFolders(folderName,fileType,frontConstraint)
%finds all images within 'folderName' (recursively) whose names end in
%'fileType' and start with 'frontConstraint

    if nargin==1
        fileType = '.tiff';
    end

    if nargin < 3 || isempty(frontConstraint) == 1
        frontConstraint = '';
    end

    if folderName(end) ~= '/'
        folderName = strcat(folderName, '/');
    end

    [~,temp] = unix(['find ' folderName ' -name ' frontConstraint '*' fileType]);
    images = regexp(temp,'\n','split')';
    imageLengths = returnCellLengths(images);
    images = images(imageLengths > length(fileType));

author: gordonberman

repo: MotionMapper


text:
% -------------------------------------------------------------------------
% This file is copied from matconvnet-fcn repository (has minor modification)
% (https://github.com/vlfeat/matconvnet-fcn)
% -------------------------------------------------------------------------
classdef SegmentationAccuracy < dagnn.Loss

    properties (Transient)
        pixelAccuracy = 0
        meanAccuracy = 0
        meanIntersectionUnion = 0
        confusion = 0
    end
    properties
        numClasses = 0;
    end

    methods
        function outputs = forward(obj, inputs, params)
            [~,predictions] = max(inputs{1}, [], 3) ;
            predictions = gather(predictions) ;
            labels = gather(inputs{2}) ;

            % compute statistics only on accumulated pixels
            ok = labels > 0 ;
            numPixels = sum(ok(:)) ;
            obj.confusion = obj.confusion + accumarray([labels(ok),predictions(ok)],1,[obj.numClasses obj.numClasses]) ;

            % compute various statistics of the confusion matrix
            pos = sum(obj.confusion,2) ;
            res = sum(obj.confusion,1)' ;
            tp = diag(obj.confusion) ;

            obj.pixelAccuracy = sum(tp) / max(1,sum(obj.confusion(:))) ;
            obj.meanAccuracy = mean(tp ./ max(1, pos)) ;
            obj.meanIntersectionUnion = mean(tp ./ max(1, pos + res - tp)) ;

            obj.average = [obj.pixelAccuracy ; obj.meanAccuracy ; obj.meanIntersectionUnion] ;
            obj.numAveraged = obj.numAveraged + numPixels ;
            outputs{1} = obj.average ;
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
            derInputs{1} = [] ;
            derInputs{2} = [] ;
            derParams = {} ;
        end

        function reset(obj)
            obj.confusion = 0 ;
            obj.pixelAccuracy = 0 ;
            obj.meanAccuracy = 0 ;
            obj.meanIntersectionUnion = 0 ;
            obj.average = [0;0;0] ;
            obj.numAveraged = 0 ;
        end

        function str = toString(obj)
            str = sprintf('acc:%.2f, mAcc:%.2f, mIU:%.2f', ...
                obj.pixelAccuracy, obj.meanAccuracy, obj.meanIntersectionUnion) ;
        end

        function obj = SegmentationAccuracy(varargin)
            obj.load(varargin) ;
            obj.numClasses = obj.numClasses;
        end
    end
end

author: aimerykong


text:
% Fig. 6.43   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clear all;
%close all;
clf

num=86*conv([1 1],[1 2 43.25]);
p=conv([1 2 82],[1 2 101]);
den=conv([1 0 0],p);
w=logspace(-1,2,1000);
[mag,phas]=bode(num,den,w);
figure(1)
loglog(w,mag,w,ones(size(w)));
axis([.1 100 .01 10])
xlabel('\omega (rad/sec)');
ylabel('magnitude');
title('Fig. 6.43 : Bode plot for Example 6.12 (a) magnitude');
bodegrid;
%pause;
figure(2)
semilogx(w,phas,w,-180*ones(size(w)));
xlabel('\omega (rad/sec)');
ylabel('phase (deg)');
title('Fig. 6.43 : Bode plot for Example 6.13 (b) phase');
bodegrid;

% actually, the GMs and PMs are slightly different for the system here
% than that described in the text, the differences arose from
% roundoff because the initial design was done by hand.  The ideas
% are unchanged.


text:
clc, clear
x1=[6.683, 6.681, 6.676, 6.678, 6.679, 6.672];
x2=[6.661, 6.661, 6.667, 6.667, 6.664];
[h,p,ci,st]=ttest2(x1,x2,'Alpha',0.1)

author: Eurus-Holmes


text:
% @Author: aaronmishkin
% @Date:   2018-06-07T19:03:42-07:00
% @Email:  amishkin@cs.ubc.ca
% @Last modified by:   aaronmishkin
% @Last modified time: 2018-07-26T13:09:11-07:00

function [y, X, y_te, X_te] = get_data_log_reg(name, seed)
%name: name of a dataset
%seed: seed used for creating a train/test set
%y (y_te): labels for a train (test) set
%X (X_te): features for a train (test) set
%Note: for classification problems, we use 0-1 encoding in labels.
switch name
case 'murphy_synth'
    setSeed(seed)
    N=30;
    D=2;
    mu1=[ones(N,1) 5*ones(N,1)];
    mu2=[-5*ones(N,1) 1*ones(N,1)];
    class1_std = 1;
    class2_std = 1.1;
    X = [class1_std*randn(N,2)+mu1;2*class2_std*randn(N,2)+mu2];
    y = [ones(N,1);zeros(N,1)];
    y = 2*y-1;
    % Set the test set to the be same as the training set.
    X_te = X; y_te = y;
case 'synth'
  setSeed(seed);
  N = 5000;
  D = 2;
  s2 = .01;
  X = randn(N,D);
  D = D + 1;
  w = [0.1; -1; +1];
  eta = [ones(N,1) X]*w + s2*randn(N,1);
  y = sign(eta);
  [X, y, X_te, y_te] = split_data(y, X, 0.5);

case {'a2a','a3a','a4a','a5a','a6a','a7a'}
    setSeed(seed);
  load(name);
  [N,D] = size(X);
  X = [ones(N,1) X];
  [N_te,D] = size(X_te);
  X_te = [ones(N_te,1) X_te];
  y = (y+1)/2;
  y_te = (y_te+1)/2;
  assert ( length(unique(abs(y))) == 2)
  assert ( length(unique(abs(2*y-1))) == 1)

case {'svmguide3'}
    setSeed(seed);
  load(name);
  [N,D] = size(X);
  X = [ones(N,1) X];
  [N_te,D] = size(X_te);
  X_te = [ones(N_te,1) X_te];
  y = (y+1)/2;
  y_te = (y_te+1)/2;
  assert ( length(unique(abs(y))) == 2)
  assert ( length(unique(abs(2*y-1))) == 1)

case 'svmguide1'
    setSeed(seed);
  load(name);
  [N,D] = size(X);
  X = [ones(N,1) X];
  [N_te,D] = size(X_te);
  X_te = [ones(N_te,1) X_te];
  assert ( length(unique(abs(y))) == 2)
  assert ( length(unique(abs(2*y-1))) == 1)

case 'a1a'
    setSeed(seed);
  load(name);
  [N,D] = size(X);
  X = [ones(N,1) X zeros(N,4)];
  [N_te,D] = size(X_te);
  X_te = [ones(N_te,1) X_te];
  y = (y+1)/2;
  y_te = (y_te+1)/2;
  assert ( length(unique(abs(y))) == 2)
  assert ( length(unique(abs(2*y-1))) == 1)

case {'colon-cancer'}
  load(name);
  [N,D] = size(X);
  X = [ones(N,1) X];
  y = (y+1)/2;
  setSeed(seed);
  [X, y, X_te, y_te] = split_data(y, X, 0.5);

case {'duke'}
  load(name);
  [N,D] = size(X);
  X = [ones(N,1) full(X)];
  y = (y+1)/2;
  setSeed(seed);
  [X, y, X_te, y_te] = split_data(y, X, 0.5);

case {'leukemia'}
    setSeed(seed);
  load(name);
  [N,D] = size(X);
  X = [ones(N,1) X];
  [N_te,D] = size(X_te);
  X_te = [ones(N_te,1) X_te];
  y = (y+1)/2;
  y_te = (y_te+1)/2;
  assert ( length(unique(abs(y))) == 2)
  assert ( length(unique(abs(2*y-1))) == 1)

case {'gisette_scale'}
  load(name);
  X = [X; X_te];
  y = [y; y_te];
  [N,D] = size(X);
  X = [ones(N,1) full(X)];
  y = (y+1)/2;
  unique(y)
  setSeed(seed);
  [X, y, X_te, y_te] = split_data(y, X, 0.5);

case {'covtype_binary_scale'}
  load(name);
  [N,D] = size(X);
  X = [ones(N,1) full(X)];
  y = y-1;
  setSeed(seed);
  [X, y, X_te, y_te] = split_data(y, X, 0.5);
  assert ( length(unique(abs(y))) == 2)
  assert ( length(unique(abs(2*y-1))) == 1)

case {'SUSY'}
  load('SUSY.amat','-mat');
  [N,D] = size(X);
  X = [ones(N,1) full(X)];
  setSeed(seed);
  [X, y, X_te, y_te] = split_data(y, X, 0.5);
  assert ( length(unique(abs(y))) == 2)
  assert ( length(unique(abs(2*y-1))) == 1)

case {'australian_scale', 'diabetes_scale'}
  load(name);
  [N,D] = size(X);
  X = [ones(N,1) full(X)];
  y = (y+1)/2;
  setSeed(seed);
  [X, y, X_te, y_te] = split_data(y, X, 0.5);
  assert ( length(unique(abs(y))) == 2)
  assert ( length(unique(abs(2*y-1))) == 1)

case {'breast_cancer_scale'}
    setSeed(seed);
  load(name);
  [N,D] = size(X);
  X = [ones(N,1) full(X)];
  y = (y-2)/2;
  [X, y, X_te, y_te] = split_data(y, X, 0.5);
  assert ( length(unique(abs(y))) == 2)
  assert ( length(unique(abs(2*y-1))) == 1)

case 'usps_3vs5'
  setSeed(seed);
  load('usps_resampled');
  y = ([train_labels test_labels] + 1)/2; % 1540 obs
  X = ([train_patterns test_patterns]);
  y = sum(bsxfun(@times, y, [0:9]'));
  idx = find(or((y==3), (y==5)));
  y = y(idx);
  y = (y==5);
  X = X(:,idx);
  X = X'; % 1540x256
  [N,D] = size(X);
  y = y(:); % in 0/1 encoding
  X = [ones(N,1) full(X)];
  [X, y, X_te, y_te] = split_data(y, X, 0.5);
%
%   y = 2*y - 1;
%   y_te = 2*y_te - 1;
%

otherwise
  error('no such name');
end
end

function [XTr, yTr, XTe, yTe] = split_data(y, X, prop)

  N = size(y,1);
	idx = randperm(N);
  Ntr = floor(prop * N);
	idxTr = idx(1:Ntr);
	idxTe = idx(Ntr+1:end);
  XTr = X(idxTr,:);
  yTr = y(idxTr);
  XTe = X(idxTe,:);
  yTe = y(idxTe);
end

author: emtiyaz

repo: vadam

save_path: github-repos/MATLAB/emtiyaz-vadam


text:
function title = p00_title ( problem )

%*****************************************************************************80
%
%% P00_TITLE returns a title for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the number of the problem.
%
%    Output, string TITLE, a title for the problem.
%
  if ( problem == 1 )
    title = p01_title ( );
  elseif ( problem == 2 )
    title = p02_title ( );
  elseif ( problem == 3 )
    title = p03_title ( );
  elseif ( problem == 4 )
    title = p04_title ( );
  elseif ( problem == 5 )
    title = p05_title ( );
  elseif ( problem == 6 )
    title = p06_title ( );
  elseif ( problem == 7 )
    title = p07_title ( );
  elseif ( problem == 8 )
    title = p08_title ( );
  elseif ( problem == 9 )
    title = p09_title ( );
  elseif ( problem == 10 )
    title = p10_title ( );
  elseif ( problem == 11 )
    title = p11_title ( );
  elseif ( problem == 12 )
    title = p12_title ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_TITLE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of PROBLEM =  %d\n', problem );
    error ( 'P00_TITLE - Fatal error!' );
  end

  return
end


text:
%This Script saves all the variables from the current workspace into a
%single structure array.
%Author: Andres Gonzalez
%Year:2012
%Version:0.2

%#########################################################################
%Summary:
%This function allows to save all the variables from the 'caller' workspace into a struct array
%Description:
%Sometimes you need to save the variables from your base workspace, but using "save" function will have them all stored individually so if you reload them into a new workspace it could be a mess, and some variables could be overwritten.
%With this function, you can save all of them into a struct array, and so they'll be nicely packaged and ready to be saved to a .mat file that, when reloaded, will be easy to identify.
%Tags:
%save, struct, base, workspace, structure array

%Example:
% a='LALALA'
% b=[1:12:258]
% c={'cell1', 'cell2', 'cell3'}
% d=768
% e=true(3)
% theworkspace=ws2struct();
% theworkspace =
%
%     a: 'LALALA'
%     b: [1x22 double]
%     c: {'cell1'  'cell2'  'cell3'}
%     d: 768
%     e: [3x3 logical]

function WStruct=ws2struct()

WSVARS = evalin('caller', 'who');
for wscon=1:size(WSVARS,1)
    thisvar=evalin('caller', WSVARS{wscon});
    THEWORKSPACE.(WSVARS{wscon})=thisvar;
end

WStruct=THEWORKSPACE;

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function sc = patchCorr(imId,pix,ptch,method)

% PATCHCORR  Correlation score of a patch in an image
%   PATCHCORR(IMID,PIX,PATCH,METHOD) computes the correlation score
%   between PATCH and the corresponding patch centered at pixel
%   PIX in the _global_ image Image{IMID}, using METHOD.
%
%   PATCH is a structure containing:
%     I:   Matrix with pixellic values of the patch
%     SI:  Sum of all pixellic values
%     SII: Sum of all squared pixellic values
%
%   METHOD is one of the following strings 'zncc', 'ssd', 'census'
%
%   See also ZNCC, SSD, CENSUS

%   (c) 2005 Joan Sola

% patch size
pSze = size(ptch.I);

% pixel centered image patch
iPatch = pix2patch(imId,pix,pSze(2),pSze(1));

if nargin<4
    method = 'zncc';
end

switch lower(method)
    case 'zncc'
        sc = zncc(...
            ptch.I,iPatch.I,...
            ptch.SI,ptch.SII,...
            iPatch.SI,iPatch.SII);
    case 'ssd'
        sc = ssd(ptch,iPatch);
    case 'census'
        sc = census(ptch,iPatch);
    otherwise
        error('Unknown correlation method.')
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

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

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function asb_vrmfunc(hBlk)
%ASB_VRMFUNC  Helper function for On-Board Camera model

% Copyright 2018 The MathWorks, Inc.

vrmfunc('FnInit', hBlk);
vrmfunc('FnSetupVideo', hBlk);
ud = get_param(hBlk, 'UserData');
if ~strcmp(get(ud.VideoFigure, 'Viewpoint'),'Quadcopter Camera')
    set(ud.VideoFigure, 'Viewpoint', 'Quadcopter Camera');
end
end

author: Firmament-Autopilot

repo: FMT-Model


text:
function fx2 = p02_fx2 ( x )

%*****************************************************************************80
%
%% P02_FX2 evaluates the second derivative of the function for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the abscissa.
%
%    Output, real FX2, the second derivative of the function at X.
%
  fx2 = - exp ( - x );

  return
end


text:
function Control = setupController(P)

if P.Sim.springDoublePendulum
    %Then set everything to zero except for the proportional gains on the
    %leg length actuators. This turns them both into springs, which makes
    %the system conserve energy. This is a handle error check.

    %The Foot One actuator operates on the absolute angle of leg one
    Control.footOne.nomPos = 0;
    Control.footOne.nomVel = 0;
    Control.footOne.Kp = 0;
    Control.footOne.Kd = 0;

    %The Foot Two actuator operates on the absolute angle of leg two
    Control.footTwo.nomPos = 0;
    Control.footTwo.nomVel = 0;
    Control.footTwo.Kp = 0;
    Control.footTwo.Kd = 0;

    %The Leg One actuator operates on the length of leg one
    Control.legOne.nomPos = 1; %(m)
    Control.legOne.nomVel = 0;
    Control.legOne.Kp = 35;
    Control.legOne.Kd = 0;

    %The Leg Two actuator operates on the length of leg two
    Control.legTwo.nomPos = 1; %(m)
    Control.legTwo.nomVel = 0;
    Control.legTwo.Kp = 35;
    Control.legTwo.Kd = 0;

    %The Hip actuator operates on the angle between the legs: th2-th1
    Control.hip.nomPos = 0;
    Control.hip.nomVel = 0;
    Control.hip.Kp = 0;
    Control.hip.Kd = 0;
else
    %Allow system to be fully controlled. This will not necessarily
    %conserve energy.

    %The Foot One actuator operates on the absolute angle of leg one
    Control.footOne.nomPos = 0;
    Control.footOne.nomVel = 0;
    Control.footOne.Kp = 0;
    Control.footOne.Kd = 0;

    %The Foot Two actuator operates on the absolute angle of leg two
    Control.footTwo.nomPos = 0;
    Control.footTwo.nomVel = 0;
    Control.footTwo.Kp = 0;
    Control.footTwo.Kd = 0;

    %The Leg One actuator operates on the length of leg one
    Control.legOne.nomPos = 1; %(m)
    Control.legOne.nomVel = 0;
    Control.legOne.Kp = 25;
    Control.legOne.Kd = 10;

    %The Leg Two actuator operates on the length of leg two
    Control.legTwo.nomPos = 1; %(m)
    Control.legTwo.nomVel = 0;
    Control.legTwo.Kp = 25;
    Control.legTwo.Kd = 10;

    %The Hip actuator operates on the angle between the legs: th2-th1
    Control.hip.nomPos = 0;
    Control.hip.nomVel = 0;
    Control.hip.Kp = 100;
    Control.hip.Kd = 80;
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [ i_choose, w_choose ] = subset_sum_task ( weights, target, range )

%*****************************************************************************80
%
%% SUBSET_SUM_TASK seeks a subset of a set that has a given sum.
%
%  Discussion:
%
%    This function tries to compute a target value as the sum of
%    a selected subset of a given set of weights.
%
%    This function works by brute force, that is, it tries every
%    possible subset to see if it sums to the desired value.
%
%    Given N weights, every possible selection can be described by
%    one of the N-digit binary numbers from 0 to 2^N-1.
%
%    This function includes a range, which allows the user to
%    control which subsets are to be checked.  Thus, if there are
%    N weights, specifying a range of [ 0, 2^N-1] indicates that
%    all subsets should be checked.  On the other hand, this full
%    range could be broken down into smaller subranges, each of
%    which could be checked independently.
%
%    It is possible that, in the given range, there may be multiple
%    solutions of the problem.  This function will only return
%    one such solution, if found.
%
%  Example:
%
%    [ c, wv ] = subset_sum_task ( [1 2 4 8 16 32], 22, [0 2^6 - 1] )
%
%    c = 2 3 5
%    wv = 2 4 16
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 May 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer WEIGHTS(N), a set of weights.
%
%    Input, integer TARGET, the target value.
%
%    Input, integer RANGE(2), the lower and upper limits to be searched.
%
%    Output, integer I_CHOOSE(K), the indices of the chosen weights.
%
%    Output, integer W_CHOOSE(K), the values of the chosen weights.
%

%
%  Initialize the output to empty matrices.
%
  i_choose = [];
  w_choose = [];
%
%  Iterate through all possible combinations in the provided range.
%
  n = length ( weights );

  for i = range(1) : range(2)
%
%  Find the combination of weights represented by the current attempt.
%
    choose = find ( bitget ( i, 1:n ) );
%
%  If the sum of those weights matches the target, return combination.
%
    if ( sum ( weights(choose) ) == target )
      i_choose = choose;
      w_choose = weights(choose);
      return
    end

  end

  return
end


text:
function y_buffer = BF_MakeBuffer(y,bufferSize)
% BF_MakeBuffer    Make a buffered version of a time series
%
% y_buffer contains segments (rows) of length bufferSize (columns) corresponding
% to consecutive, non-overlapping segments of the series of length bufferSize

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

N = length(y);

numBuffers = floor(N/bufferSize);

% May need trimming:
y_buffer = y(1:numBuffers*bufferSize);

% Then reshape:
y_buffer = reshape(y_buffer,bufferSize,numBuffers)';

% (Each buffer is a contiguous subsequence of the time series; a row of y_buffer)

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function x = cvx_setdual( x, y )
x = cvx_setdual( cvx( x ), y );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
classdef (Abstract) SolverApplication < handle
    % SolverApplication defines an abstract interface class for NLP solvers
    %
    %
    % @author Ayonga Hereid @date 2016-10-21
    %
    % Copyright (c) 2016, AMBER Lab
    % All right reserved.
    %
    % Redistribution and use in source and binary forms, with or without
    % modification, are permitted only in compliance with the BSD 3-Clause
    % license, see
    % http://www.opensource.org/licenses/bsd-license.php

    properties

        % A structure for NLP solver options
        %
        % @type struct
        Options

    end

    methods (Abstract)

        % An abstract method that initialize the solver to be ready to
        % solve the given NLP problem
        initialize(obj, nlp);

        % An abstract method that run the NLP problem
        optimize(obj, nlp);
    end

    methods
        function obj = SolverApplication()
            % The default constructor function

            obj.Options = struct;
        end

    end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
%This Matlab script can be used to reproduce Figures 4.21-22 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%Number of BSs
L = 16;

%Number of UEs per BS
K = 10;

%Number of BS antennas
M = 100;

%Define the pilot reuse factor
f = 1;

%Select the number of setups with random UE locations
nbrOfSetups = 10;

%Select the number of channel realizations per setup
nbrOfRealizations = 1000;

%% Propagation parameters

%Communication bandwidth
B = 20e6;

%Total uplink transmit power per UE (mW)
p = 100;

%Total downlink transmit power per UE (mW)
rho = 100;

%Noise figure at the BS (in dB)
noiseFigure = 7;

%Compute noise power
noiseVariancedBm = -174 + 10*log10(B) + noiseFigure;

%Use the approximation of the Gaussian local scattering model
accuracy = 2;

%Angular standard deviation in the local scattering model (in degrees)
ASDdeg = 10;

%Prepare to save simulation results
sumchannelGains_MR = zeros(K,L,K,L,nbrOfSetups);
sumchannelGains_RZF = zeros(K,L,K,L,nbrOfSetups);
sumchannelGains_MMMSE = zeros(K,L,K,L,nbrOfSetups);

%% Go through all setups
for n = 1:nbrOfSetups

    %Output simulation progress
    disp([num2str(n) ' setups out of ' num2str(nbrOfSetups)]);

    %Compute channel statistics for one setup
    [R,channelGaindB] = functionExampleSetup(L,K,M,accuracy,ASDdeg);

    %Compute the normalized average channel gain, where the normalization
    %is based on the noise power
    channelGainOverNoise = channelGaindB - noiseVariancedBm;

    %Generate channel realizations with estimates and estimation
    %error correlation matrices
    [Hhat,C,tau_p,~,H] = functionChannelEstimates(R,channelGainOverNoise,nbrOfRealizations,M,K,L,p,f);

    %Compute channel gain terms between all users
    [channelGains_MR,channelGains_RZF,channelGains_MMMSE] = functionComputeULDLPowerLevels(H,Hhat,C,nbrOfRealizations,M,K,L,p);

    %Delete large matrices
    clear Hhat C H R;

    %Save results
    sumchannelGains_MR(:,:,:,:,n) = channelGains_MR;
    sumchannelGains_RZF(:,:,:,:,n) = channelGains_RZF;
    sumchannelGains_MMMSE(:,:,:,:,n) = channelGains_MMMSE;

end

%Prepare to compute signal and interference gains in UL and DL
signalGainsUL_MR = zeros(K,L,nbrOfSetups);
interferenceGainsUL_MR = zeros(K,L,nbrOfSetups);

signalGainsDL_MR = zeros(K,L,nbrOfSetups);
interferenceGainsDL_MR = zeros(K,L,nbrOfSetups);

signalGainsUL_MMMSE = zeros(K,L,nbrOfSetups);
interferenceGainsUL_MMMSE = zeros(K,L,nbrOfSetups);

signalGainsDL_MMMSE = zeros(K,L,nbrOfSetups);
interferenceGainsDL_MMMSE = zeros(K,L,nbrOfSetups);

%Go through all cells
for j = 1:L

    %Go through all UEs
    for k = 1:K

        %Extract the interference gains with MR in the UL
        signalGainsUL_MR(k,j,:) = reshape(p*sumchannelGains_MR(k,j,k,j,:),[1 1 nbrOfSetups]);
        interferenceGainsUL_MR(k,j,:) = reshape(p*sum(sum(sumchannelGains_MR(k,j,:,:,:),3),4),[1 1 nbrOfSetups])-signalGainsUL_MR(k,j,:);

        %Extract the interference gains with MR in the DL
        signalGainsDL_MR(k,j,:) = reshape(rho*sumchannelGains_MR(k,j,k,j,:),[1 1 nbrOfSetups]);
        interferenceGainsDL_MR(k,j,:) = reshape(rho*sum(sum(sumchannelGains_MR(:,:,k,j,:),1),2),[1 1 nbrOfSetups])-signalGainsDL_MR(k,j,:);

        %Extract the interference gains with M-MMSE in the UL
        signalGainsUL_MMMSE(k,j,:) = reshape(p*sumchannelGains_MMMSE(k,j,k,j,:),[1 1 nbrOfSetups]);
        interferenceGainsUL_MMMSE(k,j,:) = reshape(p*sum(sum(sumchannelGains_MMMSE(k,j,:,:,:),3),4),[1 1 nbrOfSetups])-signalGainsUL_MMMSE(k,j,:);

        %Extract the interference gains with M-MMSE in the DL
        signalGainsDL_MMMSE(k,j,:) = reshape(rho*sumchannelGains_MMMSE(k,j,k,j,:),[1 1 nbrOfSetups]);
        interferenceGainsDL_MMMSE(k,j,:) = reshape(rho*sum(sum(sumchannelGains_MMMSE(:,:,k,j,:),1),2),[1 1 nbrOfSetups])-signalGainsDL_MMMSE(k,j,:);

    end

end

%Plot simulation results in the UL with MR
figure; hold on; box on;

plot(10*log10(signalGainsUL_MR(:)),10*log10(interferenceGainsUL_MR(:)),'k*','LineWidth',1);
xlabel('Signal power over noise floor [dB]');
ylabel('Interference power over noise floor [dB]');
axis([-10 70 0 60]);

%Plot simulation results in the UL with M-MMSE
figure; hold on; box on;

plot(10*log10(signalGainsUL_MMMSE(:)),10*log10(interferenceGainsUL_MMMSE(:)),'ks','LineWidth',1);
xlabel('Signal power over noise floor [dB]');
ylabel('Interference power over noise floor [dB]');
axis([-10 70 0 60]);

%Plot simulation results in the DL with MR
figure; hold on; box on;

plot(10*log10(signalGainsDL_MR(:)),10*log10(interferenceGainsDL_MR(:)),'k*','LineWidth',1);
xlabel('Signal power over noise floor [dB]');
ylabel('Interference power over noise floor [dB]');
axis([-10 70 0 60]);

%Plot simulation results in the DL with M-MMSE
figure; hold on; box on;

plot(10*log10(signalGainsDL_MMMSE(:)),10*log10(interferenceGainsDL_MMMSE(:)),'ks','LineWidth',1);
xlabel('Signal power over noise floor [dB]');
ylabel('Interference power over noise floor [dB]');
axis([-10 70 0 60]);

author: emilbjornson

repo: massivemimobook


text:
function imResult = blendMode_Screen(A, B, offsetW, offsetH)
%% Screen blending mode: the values of the pixels in the two layers are
%   inverted, multiplied, and then inverted again. This yields the opposite
%   effect to multiply. The result is a brighter picture.
%
% Input:
%       A       -       Base Image
%       B       -       Top Image
%   offsetW     -   move picture B horizontally in respect to the top-left
%                   corner of picture A. Default value = 1.
%   offsetH     -   move picture B vertically in respect to the top-left
%                   corner of picture A. Default value = 1.
%
% Output:
%       imResult    -   Result of the blending, having the same size of the
%                       Base Image A.
%

%% Check Input
a = size(A);
b = size(B);
blendMode_checkInput(nargin, a, b, func2str(@blendMode_Screen));

if nargin < 3
    offsetW = 1;
    offsetH = 1;
end

if nargin < 4
    offsetH = 1;
end

%% Implementation
imResult = A;

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    [A, B] = blendMode_ResizeImages(A, B, a, b, offsetW, offsetH);
end

C = 1 - (1 - A) .* (1 - B);

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    imResult = blendMode_CreateResult(imResult, C, offsetW, offsetH);
else
    imResult = C;
end

end

author: Sable

repo: mcbench-benchmarks


text:
function ftrs = wholeMovieFtrs(moviename,trackfilename,blockSize)

tracks = load_tracks(trackfilename);
[readfcn,nframes,fid,headerinfo] = get_readframe_fcn(moviename);

if fid>0,
  fclose(fid);
end

nblocks = ceil(nframes/blockSize);

parfor ndx = 1:nblocks
  [readfcn,nframes,fid,headerinfo] = get_readframe_fcn(moviename);
  % For parfor files has to be opened for every loop.
  blockStart = (ndx-1)*blockSize + 1;
  blockEnd   = min(ndx*blockSize,nframes);
  ftrs{ndx} = genFeatures(readfcn,headerinfo,blockStart,blockEnd,tracks);
  if fid>0,
    fclose(fid);
  end
  fprintf('.');
  if mod(ndx,50)==0, fprintf('\n'); end
end

author: kristinbranson

repo: JAABA


text:
function AutoClusterTest(hObject,~)
%hfig = getParentFigure(hObject);
hfig = hObject;% pass it hFig for now

cIX = getappdata(hfig,'cIX');
gIX = getappdata(hfig,'gIX');
M = getappdata(hfig,'M');
M_0 = getappdata(hfig,'M_0');
thres_size = 10;
thres_split = getappdata(hfig,'thres_split');
thres_stimlock = 1.0;
thres_merge = getappdata(hfig,'thres_merge');
thres_silh = 0.4;

isWkmeans = getappdata(hfig,'isWkmeans');

%% kmeans
if isWkmeans,
    numK = 20;
    disp(['kmeans k = ' num2str(numK)]);
    tic
    rng('default');% default = 0, but can try different seeds if doesn't converge
    if numel(M)*numK < 10^7 && numK~=1,
        disp('Replicates = 5');
        gIX = kmeans(M,numK,'distance','correlation','Replicates',5);
    elseif numel(M)*numK < 10^8 && numK~=1,
        disp('Replicates = 3');
        gIX = kmeans(M,numK,'distance','correlation','Replicates',3);
    else
        gIX = kmeans(M,numK,'distance','correlation');
    end
    toc
    SaveCluster_Direct(hfig,cIX,gIX,'k=20');
end
[gIX, numU] = SqueezeGroupIX(gIX);

%% pushbutton_iter_split(hObject,~);
disp('iter. split all...');
I_rest = [];
iter = 1;
gIX_last = gIX;
I_clean_last = cIX;
cIX = [];
gIX = [];
for i = 1:numU,
    disp(['i = ' num2str(i)]);
    ix = gIX_last == i;
    IX = I_clean_last(ix);
    M_s = M_0(IX,:);
    [I_rest,cIX,gIX,numU] = CleanClus(M_s,IX,I_rest,cIX,gIX,numU,1-thres_split,thres_size);
end
[gIX, ~] = SqueezeGroupIX(gIX);
if isempty(gIX),
    errordlg('nothing to display!');
    return;
end
SaveCluster_Direct(hfig,cIX,gIX,['clean_round' num2str(iter)]);
SaveCluster_Direct(hfig,I_rest,ones(length(I_rest),1),['rest_round' num2str(iter)]);

[gIX, numU] = Merge_direct(thres_merge,M_0,cIX,gIX);

%% rank by stim-lock
disp('stim-lock');
UpdateIndices(hfig,cIX,gIX,numU);
[gIX,rankscore] = RankByStimLock_Direct(hfig,gIX,numU);
disp('ranking complete');
% and threshold
IX = find(rankscore<thres_stimlock);
ix = ismember(gIX,IX);
gIX = gIX(ix);
cIX = cIX(ix);
UpdateIndices(hfig,cIX,gIX);

%% Regression with the centroid of each cluster
[cIX,gIX,~] = AllCentroidRegression_direct(hfig);
disp('auto-reg-clus complete');

[gIX, numU] = Merge_direct(thres_merge,M_0,cIX,gIX);
SaveCluster_Direct(hfig,cIX,gIX,'clean_round2');

%% Silhouette
disp('silhouette analysis');
gIX_last = gIX;
for i = 1:numU,
    disp(['i = ' num2str(i)]);
    IX = find(gIX_last == i);
    cIX_2 = cIX(IX);
    M_s = M_0(cIX_2,:);
    % try k-means with k=2, see whether to keep
    gIX_ = kmeans(M_s,2,'distance','correlation');
    silh = silhouette(M_s,gIX_,'correlation');
    if mean(silh)>thres_silh,
        % keep the k-means k=2 subsplit
        disp('split');
        gIX(IX) = gIX_ + numU; % reassign (much larger) gIX
    end
end
[gIX, ~] = SqueezeGroupIX(gIX);

%% rank by stim-lock ?? bug?
% disp('stim-lock');
% M = M_0(cIX,:);
% [gIX,rankscore] = RankByStimLock_Direct(hfig,cIX,gIX,M,numU);
% disp('ranking complete');
% % and threshold
% IX = find(rankscore<thres_stimlock);
% ix = ismember(gIX,IX);
% gIX = gIX(ix);
% cIX = cIX(ix);
%
% [gIX, ~] = Merge_direct(thres_merge,M_0,cIX,gIX);

% size threshold
thres_size = getappdata(hfig,'thres_size');
[cIX, gIX, numU] = ThresSize(cIX,gIX,thres_size);

%% update GUI
if isempty(gIX),
    errordlg('nothing to display!');
    return;
end
UpdateIndices(hfig,cIX,gIX,numU);
RefreshFigure(hfig);

SaveCluster_Direct(hfig,cIX,gIX,'clean_round3');
beep;

end

%% Functions
function pushbutton_thressize_Callback(hObject,~)
hfig = getParentFigure(hObject);
cIX = getappdata(hfig,'cIX');
gIX = getappdata(hfig,'gIX');
thres_size = getappdata(hfig,'thres_size');
[cIX, gIX, numU] = ThresSize(cIX,gIX,thres_size);
UpdateIndices(hfig,cIX,gIX,numU);
RefreshFigure(hfig);
end

function edit_sizethres_Callback(hObject,~)
str = get(hObject,'String');
if ~isempty(str),
    temp = textscan(str,'%f');
    thres_size = temp{:};
end
hfig = getParentFigure(hObject);
setappdata(hfig,'thres_size',thres_size);
end

function [cIX, gIX, numU] = ThresSize(cIX,gIX,thres_size)
U = unique(gIX);
numU = length(U);
for i=1:numU,
    if length(find(gIX==U(i)))<thres_size,
        cIX(gIX==U(i)) = [];
        gIX(gIX==U(i)) = [];
    end
end
[gIX, numU] = SqueezeGroupIX(gIX);
end

function [gIX, numU] = Merge_direct(thres_merge,M_0,cIX,gIX)
M = M_0(cIX,:);
[gIX, numU] = HierClus(M,gIX);
U = unique(gIX);
M = M_0(cIX,:);
[C,D] = FindCentroid_Direct(gIX,M);
i = 1;
while i<numU,
    c = corr(C(i,:)',C(i+1,:)');
    if c > thres_merge,
        IX = find(gIX == U(i+1));
        gIX(IX)=U(i); %#ok<*FNDSB>
        U = unique(gIX);
        numU = length(U);

        IX = find(gIX == U(i));
        M_s = M(IX,:);
        [~,C1,~,D1] = kmeans(M_s,1,'distance','correlation');
        C(i,:) = C1;
        D(i) = mean(D1);
        C(i+1,:) = [];
        D(i+1) = [];
    else
        i = i+1;
    end
end
[gIX, numU] = HierClus(M,gIX);
disp('merging complete');
end

function checkbox_wkmeans_Callback(hObject,~)
hfig = getParentFigure(hObject);
setappdata(hfig,'isWkmeans',get(hObject,'Value'));
end
function pushbutton_merge_Callback(hObject,~)
% disp('merging...');
hfig = getParentFigure(hObject);
cIX = getappdata(hfig,'cIX');
gIX = getappdata(hfig,'gIX');
M = getappdata(hfig,'M');
U = unique(gIX);
numU = length(U);
[C,D] = FindCentroid(hfig);

thres_merge = getappdata(hfig,'thres_merge');

i = 1;
while i<numU,
    c = corr(C(i,:)',C(i+1,:)');
    if c > thres_merge,
        IX = find(gIX == U(i+1));
        gIX(IX)=U(i);
        U = unique(gIX);
        numU = length(U);

        IX = find(gIX == U(i));
        M_s = M(IX,:);
        [~,C1,~,D1] = kmeans(M_s,1,'distance','correlation');
        C(i,:) = C1;
        D(i) = mean(D1);
        C(i+1,:) = [];
        D(i+1) = [];
    else
        i = i+1;
    end
end

if numU>1,
    [gIX, numU] = HierClus(M,gIX);
end

UpdateIndices(hfig,cIX,gIX,numU);
RefreshFigure(hfig);
disp('merging complete');
end

function edit_mergethres_Callback(hObject,~)
str = get(hObject,'String');
temp = textscan(str,'%f',1);
hfig = getParentFigure(hObject);
setappdata(hfig,'thres_merge',temp{:});
end

function pushbutton_iter_split(hObject,~)
hfig = getParentFigure(hObject);
cIX = getappdata(hfig,'cIX');
gIX = getappdata(hfig,'gIX');
numU = getappdata(hfig,'numK');
thres_split = getappdata(hfig,'thres_split');
M_0 = getappdata(hfig,'M_0');

disp('iter. split all, beep when done...');
thres_size = 10;
thres_H = thres_split;
% thres_H = [0.2;0.15;0.1;0.05]; % could have more rounds...

% initialization
I_rest = [];
% loop
tic
for round = 1:length(thres_H),
    disp(['round ' num2str(round) ', numU = ' num2str(numU)]);
    dthres = 1-thres_H(round);
    gIX_last = gIX;
    I_clean_last = cIX;
    cIX = [];
    gIX = [];
    for i = 1:numU,
        disp(['i = ' num2str(i)]);
        ix = find(gIX_last == i);
        %         IX = ix;
        IX = I_clean_last(ix);
        M_s = M_0(IX,:);
        [I_rest,cIX,gIX,numU] = CleanClus(M_s,IX,I_rest,cIX,gIX,numU,dthres,thres_size);
        %         cIX = I_clean_last(I_clean);
    end

    [gIX, numU] = SqueezeGroupIX(gIX);
    SaveCluster_Direct(hfig,cIX,gIX,['clean_round' num2str(round)]);

    SaveCluster_Direct(hfig,I_rest,ones(length(I_rest),1),['rest_round' num2str(round)]);
end
toc
beep
end

function [I_rest,I_clean,gIX_clean,numU] = CleanClus(M_s,IX,I_rest,I_clean,gIX_clean,numU,dthres,thres_size)
I_clean_s = [];

% find numK_s for kmeans
kmax = min(round(size(M_s,1)/thres_size),30);
% try numK_s = 1
numK_s = 1;
rng('default');
[gIX_s,~,~,D] = kmeans(M_s,numK_s,'distance','correlation');
Dist = min(D,[],2);
if mean(Dist)>dthres,
    % try numK_s = kmax
    numK_s = kmax;
    rng('default');
    [gIX_s,~,~,D] = kmeans(M_s,numK_s,'distance','correlation');
    Dist = min(D,[],2);
    if mean(Dist)<dthres, % find in between value for numK_s
        numK_s = 2;
        while 1,
            % kmeans-cluster by numK_s
            rng('default');
            [gIX_s,~,~,D] = kmeans(M_s,numK_s,'distance','correlation');
            Dist = min(D,[],2);
            if mean(Dist)<dthres,
                break;
            end

            if numK_s < kmax,
                numK_s = numK_s+1;
                disp(['numK_s = ' num2str(numK_s)]);
            else % numK_s = kmax;
                break;
            end
        end
    else disp(['numK_s = ' num2str(numK_s)]);
    end
end
% have numK_s that makes mean(Dist) < thres, or numK_s = kmax

for i = 1:numK_s,
    IX_s = find(gIX_s == i);
    if length(IX_s)>thres_size, % cluster big enough to start
        dst = Dist(IX_s);
        if mean(dst) < dthres,
            I_clean_s = [I_clean_s; IX_s]; %#ok<*AGROW>
            gIX_clean = [gIX_clean; gIX_s(IX_s)+double(numU)];
        else
            ix = find(dst<dthres); % clean
            if length(ix)>=thres_size, % clean cluster still big enough
                I = IX_s(ix);
                I_clean_s = [I_clean_s; I];
                gIX_clean = [gIX_clean; gIX_s(I)+double(numU)];
            end
        end
    end
end

numU = numU + numK_s;
I_rest_s = setdiff(1:size(M_s,1),I_clean_s);
I_rest = [I_rest; IX(I_rest_s)];
I_clean = [I_clean; IX(I_clean_s)];
end

function edit_splitthres_Callback(hObject,~)
str = get(hObject,'String');
temp = textscan(str,'%f',1);
hfig = getParentFigure(hObject);
setappdata(hfig,'thres_split',temp{:});
end

%% Internal functions

function UpdateClusGroupID(hfig,clusgroupID,new_clusgroupID,norefresh) %#ok<INUSD>
% save/update old Cluster into ClusGroup before exiting,
% as Cluster is the variable handled in hfig but not saved elsewhere
ClusGroup = getappdata(hfig,'ClusGroup');
Cluster = getappdata(hfig,'Cluster');
i_fish = getappdata(hfig,'i_fish');

% update into workspace
ClusGroup{clusgroupID} = Cluster;
setappdata(hfig,'ClusGroup',ClusGroup);
global VAR;
VAR(i_fish).ClusGroup = CurrentClusGroup(hfig);

% load new 'Cluster'
Cluster = ClusGroup{new_clusgroupID};
setappdata(hfig,'Cluster',Cluster);
setappdata(hfig,'clusgroupID',new_clusgroupID);

% update GUI: hclusgroupmenu
global hclusgroupmenu hclusgroupname;
if ishandle(hclusgroupmenu),
    menu = MakeNumberedMenu(VAR(i_fish).ClusGroupName);
    set(hclusgroupmenu,'String',menu,'Value',new_clusgroupID+1);
    set(hclusgroupname,'String',VAR(i_fish).ClusGroupName(new_clusgroupID));
end

if ~exist('norefresh','var'),
    if numel(Cluster) == 0, % i.e. for newly created ClusGroup
        SaveCluster(hfig,'new');
    else % load this ClusGroup
        clusID = 1;
        UpdateClusID(hfig,clusID);
    end
end
end

function ClusGroup = CurrentClusGroup(hfig)
ClusGroup = getappdata(hfig,'ClusGroup');
Cluster = getappdata(hfig,'Cluster');
clusgroupID = getappdata(hfig,'clusgroupID');
ClusGroup{clusgroupID} = Cluster;
setappdata(hfig,'ClusGroup',ClusGroup);
end

function [Cluster,clusID] = SaveCluster(hfig,state,clusheader,name)
cIX = getappdata(hfig,'cIX');
gIX = getappdata(hfig,'gIX');
absIX = getappdata(hfig,'absIX');
Cluster = getappdata(hfig,'Cluster');

if strcmp(state,'current'),
    clusID = getappdata(hfig,'clusID');
else %if strcmp(state,'new'),
    clusID = numel(Cluster)+1;
    if ~exist('name','var'),
        name = getappdata(hfig,'newclusname');
    end
    if ~exist('clusheader','var'),
        clusheader = getappdata(hfig,'clusheader');
    end
    Cluster(clusID).name = [clusheader name];
end

cIX_abs = absIX(cIX);

Cluster(clusID).cIX_abs = cIX_abs;
Cluster(clusID).gIX = gIX;
Cluster(clusID).numel = length(cIX_abs);
U = unique(gIX);
numU = length(U);
Cluster(clusID).numK = numU;

setappdata(hfig,'Cluster',Cluster);
UpdateClusID(hfig,clusID);
disp('cluster saved');
end

function [Cluster,clusID] = SaveCluster_Direct(hfig,cIX,gIX,name) %,clusheader,name)
new_clusgroupID = 1;
clusgroupID = getappdata(hfig,'clusgroupID');
UpdateClusGroupID(hfig,clusgroupID,new_clusgroupID);

if ~exist('cIX','var'),
    cIX = getappdata(hfig,'cIX');
end
if ~exist('gIX','var'),
    gIX = getappdata(hfig,'gIX');
end

Cluster = getappdata(hfig,'Cluster');
absIX = getappdata(hfig,'absIX');
cIX_abs = absIX(cIX);

clusID = numel(Cluster)+1;
% if ~exist('name','var'),
%     name = getappdata(hfig,'newclusname');
% end
% if ~exist('clusheader','var'),
%     clusheader = getappdata(hfig,'clusheader');
% end
Cluster(clusID).name = name; %[clusheader name];
Cluster(clusID).cIX_abs = cIX_abs;
Cluster(clusID).gIX = gIX;
Cluster(clusID).numel = length(cIX);
Cluster(clusID).numK = length(unique(gIX));

setappdata(hfig,'Cluster',Cluster);
UpdateClusID(hfig,clusID);
disp('cluster saved');
end

function UpdateClusID(hfig,clusID)
Cluster = getappdata(hfig,'Cluster');

% save
setappdata(hfig,'clusID',clusID);
% update GUI
global hclusname hclusmenu;
set(hclusname,'String',Cluster(clusID).name);
menu = MakeNumberedMenu({Cluster.name});
set(hclusmenu,'String', menu,'Value',clusID+1);
numK = Cluster(clusID).numK;
gIX = Cluster(clusID).gIX;

% convert absolute index to index used for this dataset
absIX = getappdata(hfig,'absIX');
cIX_abs = Cluster(clusID).cIX_abs;
[~,cIX] = ismember(cIX_abs,absIX);

if ~isempty(find(cIX==0,1)),
    errordlg('cell index out of bound for currently loaded dataset');
    IX = cIX==0;
    cIX(IX) = [];
    gIX(IX) = [];
end

UpdateIndices(hfig,cIX,gIX,numK);
RefreshFigure(hfig);
end

function menu = MakeNumberedMenu(name) % e.g. name = {Cluster.name} (note {})
menu = [{'(choose)'},name];
for j=2:length(menu),menu(j)={[num2str(j-1) ': ' menu{j}]};end
end

function [gIX, numU] = HierClus(M,gIX,isplotfig) %#ok<INUSD>
[gIX, numU] = SqueezeGroupIX(gIX);
[C,~] = FindCentroid_Direct(gIX,M);
D = pdist(C,'correlation');
tree = linkage(C,'average','correlation');
leafOrder = optimalleaforder(tree,D);

if numU>1,
    if exist('isplotfig','var'),
        figure('Position',[100 100 600 600]);
        %             subplot(1,3,1);
        %             CORR = corr(C');
        %             CorrPlot(CORR);
        %
        %             subplot(1,3,2);
        dendrogram(tree,numU,'orientation','right','reorder',leafOrder);
        set(gca,'YDir','reverse');
        set(gca,'XTick',[]);

        %             subplot(1,3,3);
        %             C2 = C(leafOrder,:);
        %             CORR2 = corr(C2');
        %             CorrPlot(CORR2);
    end
    % sort for uniform colorscale
    temp = zeros(size(gIX));
    for i = 1:numU,
        temp(gIX==leafOrder(i)) = i; % = T(i) for clusters segmented from tree
    end
    gIX = temp;
end
end

function gIX = HierClusDirect(C,gIX,numU)
D = pdist(C,'correlation');
tree = linkage(C,'average','correlation');
leafOrder = optimalleaforder(tree,D);

% sort for uniform colorscale
temp = zeros(size(gIX));
for i = 1:numU,
    temp(gIX==leafOrder(i)) = i; % = T(i) for clusters segmented from tree
end
gIX = temp;
end

function [gIX, numK] = SqueezeGroupIX(gIX)
U = unique(gIX);
numK = length(U);
for i = 1:numK,
    old = U(i);
    gIX(gIX==old) = i;
end
end

% frequently used, updates cell-index,group-index,cluster-number. set-operations included in here.
function UpdateIndices(hfig,cIX,gIX,numK)
global hback hopID;
if ~exist('gIX','var'),
    gIX = getappdata(hfig,'gIX');
end
if ~exist('cIX','var'),
    cIX = getappdata(hfig,'cIX');
end

% update cache
bC = getappdata(hfig,'bCache');
cIX_last = getappdata(hfig,'cIX');
gIX_last = getappdata(hfig,'gIX');
if ~(isequal(cIX_last,cIX) && isequal(gIX_last,gIX)),
    bC.cIX = [cIX_last,bC.cIX];
    bC.gIX = [gIX_last,bC.gIX];
    bC.numK = [getappdata(hfig,'numK'),bC.numK];
    set(hback,'enable','on');
    if length(bC.cIX)>20,
        bC.cIX(end) = [];
        bC.gIX(end) = [];
        bC.numK(end) = [];
    end
end

% set operations, if applicable
opID = getappdata(hfig,'opID');
if opID ~= 0,
    switch opID,
        case 1,
            disp('union');
            [~,ia,ib] = union(cIX_last,cIX,'stable');
            IX = vertcat(cIX_last(ia),cIX(ib));% IX;
        case 2,
            disp('intersect');
            [IX,ia,~] = intersect(cIX_last,cIX);
            ib = [];
        case 3,
            disp('setdiff');
            [IX,ia] = setdiff(cIX_last,cIX);
            ib = [];
        case 4,
            disp('rev setdiff');
            % swap sequence, then same as opID==3
            temp = cIX;
            cIX = cIX_last;
            cIX_last = temp;
            temp = gIX;
            gIX = gIX_last;
            gIX_last = temp;
            [IX,ia] = setdiff(cIX_last,cIX);
            ib = [];
        case 5,
            disp('setxor');
            [IX,ia,ib] = setxor(cIX_last,cIX);
        case 6,
            disp('smartUnion');
            CIX = vertcat(cIX_last,cIX);
            GIX = [gIX_last;gIX+max(gIX_last)]; % gIX to match
            M_0 = getappdata(hfig,'M_0');
            [cIX,gIX,numK] = SmartUnique(CIX,GIX,M_0(CIX,:));
    end
    if opID<6,
        if ~isempty(IX),
            cIX = IX;
            gIX = vertcat(gIX_last(ia),gIX(ib)+max(gIX_last(ia)));
            numK = length(unique(gIX));
            %         [gIX, numK] = SqueezeGroupIX(gIX);
        else
            errordlg('operation result is empty set!')
            waitforbuttonpress;
        end
    end
    set(hopID,'Value',1,'BackgroundColor',[1,1,1]); % reset
    setappdata(hfig,'opID',0);
end

setappdata(hfig,'bCache',bC);
setappdata(hfig,'cIX',cIX);
setappdata(hfig,'gIX',gIX);

M = GetTimeIndexedData(hfig);
setappdata(hfig,'M',M);

if exist('numK','var'),
    setappdata(hfig,'numK',double(numK));
end

%% Resets: reset flags the NEXT time this function is called (so they only apply to this particular plot)
% handle rankID: >=2 means write numbers as text next to colorbar
% first UpdateIndices sets rankID to 100, second sets back to 0
rankID = getappdata(hfig,'rankID');
if rankID>=2,
    if rankID==100,
        setappdata(hfig,'rankID',0);
    else
        setappdata(hfig,'rankID',100);
    end
end

% toggle 'isWeighAlpha'
isWeighAlpha = getappdata(hfig,'isWeighAlpha');
if isWeighAlpha == 1,
    setappdata(hfig,'isWeighAlpha',100);
elseif isWeighAlpha == 100,
    setappdata(hfig,'isWeighAlpha',0);
end

% FindCentroid reset:
setappdata(hfig,'Centroids',[]);
setappdata(hfig,'D_ctrd',[]);
end

% frequently used, 2 plotting functions are outside ('DrawTimeSeries.m' and 'DrawCellsOnAnatProj.m')
function RefreshFigure(hfig)
watchon; drawnow;
isPopout = 0; % with down-sampling in plots

% clean-up canvas
allAxesInFigure = findall(hfig,'type','axes');
if ~isempty(allAxesInFigure)
    delete(allAxesInFigure);
end

figure(hfig);
h1 = axes('Position',[0.05, 0.04, 0.55, 0.83]);
h2 = axes('Position',[0.63, 0.04, 0.35, 0.83]);

isCentroid = getappdata(hfig,'isCentroid');
isRefAnat = getappdata(hfig,'isRefAnat');

isPlotLines = 0; %getappdata(hfig,'isPlotLines');
isPlotBehavior = 1; %getappdata(hfig,'isPlotBehavior');

% double-check if cIX is valid
cIX = getappdata(hfig,'cIX');
if isempty(cIX),
    errordlg('empty set!');
    % GO BACK to the last step (presumably not empty)
    pushbutton_back_Callback(h1); % using h1 instaed of the usual 'hObject'
    return;
end

% left subplot
axes(h1);
DrawTimeSeries(hfig,h1,isPopout,isCentroid,isPlotLines,isPlotBehavior);

% right subplot
axes(h2);
DrawCellsOnAnatProj(hfig,isRefAnat,isPopout);
watchoff;
end

function [C,D] = FindCentroid_Direct(gIX,M)
U = unique(gIX);
numU = length(U);
C = zeros(numU,size(M,2));
D = zeros(numU,1);
for i=1:numU,
    IX = find(gIX == U(i));
    if length(IX)==1,
        C(i,:) = M(IX,:);
        D(i) = 1;
    else
        M_s = M(IX,:);
        [~,C1,~,D1] = kmeans(M_s,1,'distance','correlation');
        C(i,:) = C1;
        D(i) = mean(D1);
    end
end
end

function UpdateTimeIndex(hfig,isSkipcIX) %#ok<INUSD>
% input params
isAvr = getappdata(hfig,'isAvr');
isRawtime = getappdata(hfig,'isRawtime');
stimrange = getappdata(hfig,'stimrange');
% load
timelists = getappdata(hfig,'timelists');
periods = getappdata(hfig,'periods');
fishset = getappdata(hfig,'fishset');

if fishset == 1,
    if isAvr,
        tIX = 1:periods;
    else
        tIX = timelists{1};
    end

else % fishset>1,
    if isAvr,
        tIX = [];
        for i = 1:length(stimrange),
            ix = stimrange(i);
            i_start = sum(periods(1:ix-1)); % if ix-1<1, sum = 0
            tIX = horzcat(tIX,(i_start+1:i_start+periods(ix)));
            %             tIX = vertcat(tIX,(i_start+1:i_start+periods(ix))');
        end
    else % full range
        if ~isRawtime,
            tIX = cat(2, timelists{stimrange});
        else
            tIX = sort(cat(2, timelists{stimrange}));
        end
    end
end

setappdata(hfig,'tIX',tIX);

% set Matrices to hold time-series
M_0 = GetTimeIndexedData(hfig,'isAllCells');
setappdata(hfig,'M_0',M_0);
if ~exist('isSkipcIX','var'),
    cIX = getappdata(hfig,'cIX');
    setappdata(hfig,'M',M_0(cIX,:));
end
end

function [M,behavior,stim] = GetTimeIndexedData(hfig,isAllCells) %#ok<INUSD>
%{
% naming convention used:
M = GetTimeIndexedData(hfig);
M_0 = GetTimeIndexedData(hfig,'isAllCells');
%}

isZscore = getappdata(hfig,'isZscore');
% main data input
if ~isZscore,
    cellResp = getappdata(hfig,'CellResp');
    cellRespAvr = getappdata(hfig,'CellRespAvr');
else
    cellResp = getappdata(hfig,'CellRespZ');
    cellRespAvr = getappdata(hfig,'CellRespAvrZ');
end
Behavior_full = getappdata(hfig,'Behavior_full');
BehaviorAvr = getappdata(hfig,'BehaviorAvr');
stim_full = getappdata(hfig,'stim_full');
stimAvr = getappdata(hfig,'stimAvr');
% other params
isAvr = getappdata(hfig,'isAvr');
cIX = getappdata(hfig,'cIX');
tIX = getappdata(hfig,'tIX');

%% set data
if isAvr,
    if exist('isAllCells','var'),
        M = cellRespAvr(:,tIX);
    else
        M = cellRespAvr(cIX,tIX);
    end
    behavior = BehaviorAvr(:,tIX);
    stim = stimAvr(:,tIX);
else
    if exist('isAllCells','var'),
        M = cellResp(:,tIX);
    else
        M = cellResp(cIX,tIX);
    end
    behavior = Behavior_full(:,tIX);
    stim = stim_full(:,tIX);
end

setappdata(hfig,'behavior',behavior);
setappdata(hfig,'stim',stim);
end

function closefigure_Callback(hfig,~)
global EXPORT_autorecover;
EXPORT_autorecover = getappdata(hfig);
end

function fig = getParentFigure(fig)
% if the object is a figure or figure descendent, return the figure. Otherwise return [].
while ~isempty(fig) && ~strcmp('figure', get(fig,'type'))
    fig = get(fig,'parent');
end
end

function runscript(flag_script,var_script)
switch flag_script
    case 'push_cIX_gIX'
        UpdateIndices(var_script{:});
        RefreshFigure(var_script{1});
end
end

author: xiuyechen

repo: FishExplorer


text:
function katcherr(err_mode,err_typ)
% KATCHERR Catch error function.

% Author: Craig Borghesani
% Date: 11/15/94
% Revised:
% Copyright (c) 1999, Prentice-Hall

f = gcf;
ui_han = get(f,'userdata');
stat_bar = get(ui_han(43),'userdata');

if err_mode == 1,
 set(stat_bar,'string','An Error has occurred in TERMIMPL');
elseif err_mode == 2,
 if err_typ == 1, typ_str = 'Gain';
 elseif err_typ == 2, typ_str = 'Integrator';
 elseif err_typ == 3, typ_str = 'Differentiator';
 elseif err_typ == 4, typ_str = 'Real Pole';
 elseif err_typ == 5, typ_str = 'Real Zero';
 elseif err_typ == 6, typ_str = 'Complex Pole';
 elseif err_typ == 7, typ_str = 'Complex Zero';
 elseif err_typ == 8, typ_str = 'Lead/Lag';
 elseif err_typ == 9, typ_str = 'Num/Den';
 elseif err_typ == 10, typ_str = 'Delay';
 elseif err_typ == 11, typ_str = 'Design Point'; end
 if err_typ ~= 9,
  set(stat_bar,'string',['Incorrect type of input for ',typ_str,' element.  Please re-enter.']);
 else
  set(stat_bar,'string',['Coefficients must be in brackets. i.e. [1,2,3]']);
 end
elseif err_mode == 3,
 if err_typ == 1, typ_str = 'Proportional Gain';
 elseif err_typ == 2, typ_str = 'Integral Gain';
 elseif err_typ == 3, typ_str = 'Derivative Gain';
 elseif err_typ == 4, typ_str = 'Design Point'; end
 set(stat_bar,'string',['Incorrect type of input for ',typ_str,' element.  Please re-enter.']);
elseif err_mode == 4,
 if err_typ == 1,
  typ_str = 'Gain';
  opt_str = 'Please re-enter.';
 elseif err_typ == 4,
  typ_str = 'Real Pole';
  opt_str = 'Please use an Integrator.';
 elseif err_typ == 5,
  typ_str = 'Real Zero';
  opt_str = 'Please use a Differentiator.';
 elseif any(err_typ == [6,7,11]),
  typ_str = 'wn';
  if err_typ == 6,
   opt_str = 'Please use an Integrator.';
  elseif err_typ == 7,
   opt_str = 'Please use a Differentiator.';
  else
   opt_str = 'Please re-enter.';
  end
 end
 set(stat_bar,'string',[typ_str,' cannot be zero  ',opt_str]);
elseif err_mode == 5,
 set(stat_bar,'string','Phase value must be < 88 degrees or > -88 degrees');
elseif err_mode == 6,
 set(stat_bar,'string','Corner frequency cannot be zero');
end


text:
% irreg_demo
%
% demo of an irregularly sampled problem
%
% Author Piet M.T. Broersen, September 2008
%
% The three statements :
%           ******************************************************
%           *  [asel bsel sellog]=ARMAsel_mis(ti,xi,ARmax,Tr,w); *
%           *  psdsel=arma2psd(asel,bsel);                       *
%           *  corsel=arma2cor(asel,bsel,50);                    *
%           ******************************************************
%  compute the selected estimated spectrum and autocorrelation function.
%  Sellog gives additional information and plots spectra of all estimated
%  AR models
%
%  ti : observation times [s]
%  xi : observation amplitudes
%  ARmax : higherst candidate AR order
%  Tr : equidistant resampling distance [s]
%  w ; slot width as fraction of Tr [s]
%
%  If the mean sampling frequency f0 Hz is high in comparison with the highest
%  frequency in the desired spectrum, it is advisable to use
%  nearest neighbor resampling and armasel.
%  See the example in simple_irreg_demo.
%  This gives very accurate spectra until f0/20 Hz and
%  about 50 % error due to resampling at f0/2pi.
%  The error becomes very large for higher frequencies and
%  armasel_irreg is advised for spectra higher than f0/2pi Hz.

clear all, close all, clc, echo off

load irreg_data
% row vectors of times ti and amplitudes xi of irregularly sampled observations

N=length(ti)

figure
plot(ti,xi,'p',ti,xi,':r')
title([int2str(N),' irregular benchmark data.'])
xlabel('\rightarrow time axis [s]')
axis tight

disp('Mean time between observations')
T0=(ti(N) - ti(1))/(N-1)

disp(' ')
disp('************************************************************************')
disp('Some warning messages from the OPTIM toolbox cannot be suppressed easily')
disp('It is not necessary to provide gradient information')
disp('Warnings generated by the MATLAB OPTIM routine fminunc can be ignored')
disp('************************************************************************')
disp(' ')
disp('Input values for ARMAsel_irreg')
ARmax=3
w=1/2            %*Tres   slot width, fraction of Tres
Tr=1/2000
disp(' ')
disp('[air bir sellogir] = ARMAsel_irreg(ti,xi,ARmax,Tr,w)')

[air bir sellogir] = ARMAsel_irreg(ti,xi,ARmax,Tr,w)

disp('************************************************************************')
disp(' ')
disp('New input for w')
w=1/8
disp(' ')
disp('[air2 bir2 sellogir2] = ARMAsel_irreg(ti,xi,ARmax,Tr,w)')
[air2 bir2 sellogir2] = ARMAsel_irreg(ti,xi,ARmax,Tr,w)

disp('************************************************************************')
disp(' ')
disp('New input for w and Tr')
w=1/2            %*Tres   slot width, fraction of Tres
Tr2=1/4000
disp(' ')
disp('[air3 bir3 sellogir3] = ARMAsel_irreg(ti,xi,ARmax,Tr,w)')
[air3 bir3 sellogir3] = ARMAsel_irreg(ti,xi,ARmax,Tr2,w)

[h_ir f_ir] = arma2psd(air,bir,1000,Tr);
h_ir2 = arma2psd(air2,bir2,1000,Tr);
[h_ir3 f_ir3] = arma2psd(sellogir3.AR_sel_corrected,1,1000,Tr2);

disp(' ')
disp('  rc AR_sel for w = 1/2 and for w = 1/8, Tr = 0.0005 s')
[dum rc_sel]=ar2arset(sellogir.ar);
[dum rc_sel2]=ar2arset(sellogir2.ar);
[dum rc_sel3]=ar2arset(sellogir3.AR_sel_corrected);
disp(rc_sel)
disp(rc_sel2)
disp('  rc AR_sel-corrected for w = 1/2, Tr = 0.00025 s')
disp(rc_sel3)

figure
loglog(f_ir,h_ir,f_ir,h_ir2,f_ir3,h_ir3)
legend('\it{Tr}\rm = 0.0005 s,  \itw\rm = 1/2 * \it{Tr}\rm s', ...
       '\it{Tr}\rm = 0.0005 s,  \itw\rm = 1/8 * \it{Tr}\rm s', ...
       '\it{Tr}\rm = 0.00025 s, \itw\rm = 1/2 * \it{Tr}\rm s',3)
title(['PSD of ',int2str(N),' irregular benchmark data.'])
xlabel('\rightarrow frequency [Hz]')
ylabel('\rightarrow Logarithm of power spectral density')
axis tight