text:
function [lfp] = bz_GetLFP(varargin)
% bz_GetLFP - Get local field potentials.
%
%  Load local field potentials from disk. No longer dependent on
%  FMAT/SetCurrentSession.
%
%  USAGE
%
%    [lfp] = bz_GetLFP(channels,<options>)
%
%  INPUTS
%
%    channels(required) -must be first input, numeric
%                        list of channels to load (use keyword 'all' for all)
%                        channID is 0-indexing, a la neuroscope
%  Name-value paired inputs:
%    'basepath'           - folder in which .lfp file will be found (default
%                           is pwd)
%                           folder should follow buzcode standard:
%                           whateverPath/baseName
%                           and contain file baseName.lfp
%    'basename'           -base file name to load
%    'intervals'          -list of time intervals [0 10; 20 30] to read from
%                           the LFP file (default is [0 inf])
%    'downsample'         -factor to downsample the LFP (i.e. 'downsample',5
%                           will load a 1250Hz .lfp file at 250Hz)
%    'noPrompts'          -logical (default) to supress any user prompts
%
%  OUTPUT
%
%    lfp             struct of lfp data. Can be a single struct or an array
%                    of structs for different intervals.  lfp(1), lfp(2),
%                    etc for intervals(1,:), intervals(2,:), etc
%    .data           [Nt x Nd] matrix of the LFP data
%    .timestamps     [Nt x 1] vector of timestamps to match LFP data
%    .interval       [1 x 2] vector of start/stop times of LFP interval
%    .channels       [Nd X 1] vector of channel ID's
%    .samplingRate   LFP sampling rate [default = 1250]
%    .duration       duration, in seconds, of LFP interval
%
%
%  EXAMPLES
%
%    % channel ID 5 (= # 6), from 0 to 120 seconds
%    lfp = bz_GetLFP(5,'restrict',[0 120]);
%    % same, plus from 240.2 to 265.23 seconds
%    lfp = bz_GetLFP(5,'restrict',[0 120;240.2 265.23]);
%    % multiple channels
%    lfp = bz_GetLFP([1 2 3 4 10 17],'restrict',[0 120]);
%    % channel # 3 (= ID 2), from 0 to 120 seconds
%    lfp = bz_GetLFP(3,'restrict',[0 120],'select','number');

% Copyright (C) 2004-2011 by Michaël Zugaro
% editied by David Tingley, 2017
%
% NOTES
% -'select' option has been removed, it allowed switching between 0 and 1
%   indexing.  This should no longer be necessary with .lfp.mat structs
%
% TODO
% add saveMat input
% expand channel selection options (i.e. region or spikegroup)
% add forcereload
%% Parse the inputs!

channelsValidation = @(x) isnumeric(x) || strcmp(x,'all');

% parse args
p = inputParser;
addRequired(p,'channels',channelsValidation)
addParameter(p,'basename','',@isstr)
addParameter(p,'intervals',[],@isnumeric)
addParameter(p,'restrict',[],@isnumeric)
addParameter(p,'basepath',pwd,@isstr);
addParameter(p,'downsample',1,@isnumeric);
addParameter(p,'saveMat',false,@islogical);
addParameter(p,'forceReload',false,@islogical);
addParameter(p,'noPrompts',false,@islogical);

parse(p,varargin{:})
basename = p.Results.basename;
channels = p.Results.channels;
downsamplefactor = p.Results.downsample;
basepath = p.Results.basepath;
noPrompts = p.Results.noPrompts;

% doing this so you can use either 'intervals' or 'restrict' as parameters to do the same thing
intervals = p.Results.intervals;
restrict = p.Results.restrict;
if isempty(intervals) && isempty(restrict) % both empty
    intervals = [0 Inf];
elseif isempty(intervals) && ~isempty(restrict) % intervals empty, restrict isn't
    intervals = restrict;
end

%% let's check that there is an appropriate LFP file
if isempty(basename)
   %disp('No basename given, so we look for a *lfp/*eeg file...')
   d = dir([basepath filesep '*lfp']);
   if length(d) > 1 % we assume one .lfp file or this should break
       error('there is more than one .lfp file in this directory?');
   elseif length(d) == 0
       d = dir([basepath filesep '*eeg']);
       if isempty(d)
           error('could not find an lfp/eeg file..')
       end
   end
   lfp.Filename = d.name;
   basename = strsplit(lfp.Filename,'.');
   if length(basename) > 2
       base = [];
       for i=1:length(basename)-1
          base = [base basename{i} '.'];
       end
       basename = base(1:end-1);  % this is an fugly hack to make things work with Kenji's naming system...
   else
       basename = basename{1};
   end

else
   d = dir([basepath filesep basename '.lfp']);
   if length(d) > 1 % we assume one .lfp file or this should break
       error('there is more than one .lfp file in this directory?');
   elseif length(d) == 0
       d = dir([basepath filesep basename '.eeg']);
       if isempty(d)
           error('could not find an lfp/eeg file..')
       end
   end
   lfp.Filename = d.name;
end

%% things we can parse from sessionInfo or xml file

sessionInfo = bz_getSessionInfo(basepath, 'noPrompts', noPrompts);

try
    samplingRate = sessionInfo.lfpSampleRate;
catch
    samplingRate = sessionInfo.rates.lfp; % old ugliness we need to get rid of
end
samplingRateLFP = samplingRate./downsamplefactor;

if mod(samplingRateLFP,1)~=0
    error('samplingRate/downsamplefactor must be an integer')
end
%% Channel load options
%Right now this assumes that all means channels 0:nunchannels-1 (neuroscope
%indexing), we could also add options for this to be select region or spike
%group from the xml...
if strcmp(channels,'all')
    channels = sessionInfo.channels;
else
    %Put in something here to collapse into X-Y for consecutive channels...
    display(['Loading Channels ',num2str(channels),' (0-indexing, a la Neuroscope)'])
end

%% get the data
disp('loading LFP file...')
nIntervals = size(intervals,1);
% returns lfp/bz format
for i = 1:nIntervals
    lfp(i).duration = (intervals(i,2)-intervals(i,1));
    lfp(i).interval = [intervals(i,1) intervals(i,2)];

    % Load data and put into struct
    % we assume 0-indexing like neuroscope, but bz_LoadBinary uses 1-indexing to
    % load....
    lfp(i).data = bz_LoadBinary([basepath filesep lfp.Filename],...
        'duration',double(lfp(i).duration),...
                  'frequency',samplingRate,'nchannels',sessionInfo.nChannels,...
                  'start',double(lfp(i).interval(1)),'channels',channels+1,...
                  'downsample',downsamplefactor);
    lfp(i).timestamps = [lfp(i).interval(1):(1/samplingRateLFP):...
                        (lfp(i).interval(1)+(length(lfp(i).data)-1)/...
                        samplingRateLFP)]';
    lfp(i).channels = channels;
    lfp(i).samplingRate = samplingRateLFP;
    % check if duration is inf, and reset to actual duration...
    if lfp(i).interval(2) == inf
        lfp(i).interval(2) = length(lfp(i).timestamps)/lfp(i).samplingRate;
        lfp(i).duration = (lfp(i).interval(i,2)-lfp(i).interval(i,1));
    end

    if isfield(sessionInfo,'region') && isfield(sessionInfo,'channels')
        [~,~,regionidx] = intersect(lfp(i).channels,sessionInfo.channels,'stable');
        lfp(i).region = sessionInfo.region(regionidx); % match region order to channel order..
    end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
classdef TestTonemapDurand
    %TestTonemapDurand

    methods (Static)
        function test_1
            hdr = cv.imread(fullfile(mexopencv.root(),'test','memorial.hdr'), 'Flags',-1);
            hdr = cv.resize(hdr, 0.5, 0.5);

            obj = cv.TonemapDurand();
            obj.Gamma = 2.2;
            obj.Contrast = 4.0;
            obj.Saturation = 1.0;
            obj.SigmaSpace = 2.0;
            obj.SigmaColor = 2.0;

            ldr = obj.process(hdr);
            validateattributes(ldr, {'single'}, ...
                {'size',size(hdr), 'nonnegative'});
            %TODO: output HDR is not always in [0,1] range
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [merged_sample, new_sample, merged_sample_id, new_sample_id, distance_matrix, gram_matrix, prior_weights] = ...
    update_sample_space_model(samplesf, new_train_sample, distance_matrix, gram_matrix, prior_weights,...
    num_training_samples,params)

% Updates the sample space model
% There are 4 possible cases
% 1: Memory is not full. In this case the new train sample is placed in the
% next empty slot.
%
% 2: Memory is full, one of the samples is outdated. In this case, the new
% train sample replaces the outdated sample.
%
% 3: Memory is full, the min distance of the new train sample to the existing samples
% is less than the min distance amongst any of the existing samples. The
% new train sample is merged with the nearest existing sample
%
% 4: Memory is full, the min distance of the new train sample to the existing samples
% is more than the min distance amongst any of the existing samples. The
% closest existing samples are merged. New train sample is placed in the
% free slot

num_feature_blocks = numel(new_train_sample);

% Find the inner product of the new sample with existing samples
gram_vector = find_gram_vector(samplesf, new_train_sample, num_training_samples, params);

% Find the distance of the new sample with the existing samples
% Note: Since getting the 'exact' distance between the samples is not important,
% the distance computation is done by using only the half spectrum for
% efficiency. In practice the error incurred by this is negligible. Also
% since we wannt to merge samples that are similar, and finding the best
% match is not important, small error in the distance computation doesn't
% matter
new_train_sample_norm =  zeros(1,'like', params.data_type);

for k = 1:num_feature_blocks
    new_train_sample_norm = new_train_sample_norm + real(2*(new_train_sample{k}(:)' * new_train_sample{k}(:)));
end

dist_vector = max(new_train_sample_norm + diag(gram_matrix) - 2*gram_vector,0);
dist_vector(num_training_samples+1:end) = inf;

merged_sample = [];
new_sample = [];
merged_sample_id = -1;
new_sample_id = -1;

% Check if we have filled the memory
if num_training_samples == params.nSamples

    % Check if any sample weight is too low
    [min_sample_weight, min_sample_id] = min(prior_weights);

    if min_sample_weight < params.minimum_sample_weight
        % If any prior weight is less than the minimum allowed weight,
        % replace that sample with the new sample

        % Update distance matrix and the gram matrix
        [distance_matrix, gram_matrix] = update_distance_matrix(distance_matrix, gram_matrix, gram_vector, new_train_sample_norm, min_sample_id, -1, 0, 1);

        % Normalise the prior weights so that the new sample gets weight as
        % the learning rate
        prior_weights(min_sample_id) = 0;
        prior_weights = prior_weights*(1 - params.learning_rate)/sum(prior_weights);
        prior_weights(min_sample_id) = params.learning_rate;

        % Set the new sample and new sample position in the samplesf
        new_sample_id = min_sample_id;
        new_sample = new_train_sample;
    else
        % If no sample has low enough prior weight, then we either merge
        % the new sample with an existing sample, or merge two of the
        % existing samples and insert the new sample in the vacated
        % position

        % Find sample closest to the new sample
        [new_sample_min_dist, closest_sample_to_new_sample] = min(dist_vector);

        % Find the closest pair amongst existing samples
        [existing_samples_min_dist, closest_existing_sample_pair] = min(distance_matrix(:));
        [closest_existing_sample1,closest_existing_sample2] = ind2sub(size(distance_matrix),closest_existing_sample_pair);

        if closest_existing_sample1 == closest_existing_sample2
            error('Score matrix diagonal filled wrongly');
        end

        if new_sample_min_dist < existing_samples_min_dist
            % If the min distance of the new sample to the existing samples
            % is less than the min distance amongst any of the existing
            % samples, we merge the new sample with the nearest existing
            % sample

            % Renormalize prior weights
            prior_weights = prior_weights*(1 - params.learning_rate);

            % Set the position of the merged sample
            merged_sample_id = closest_sample_to_new_sample;

            % Extract the existing sample to merge
            existing_sample_to_merge = cell(1, 1, num_feature_blocks);

            for k = 1:num_feature_blocks
                existing_sample_to_merge{k} = samplesf{k}(merged_sample_id,:,:,:);
            end

            % Merge the new_train_sample with existing sample
            merged_sample = merge_samples(existing_sample_to_merge, new_train_sample, prior_weights(merged_sample_id), params.learning_rate, params.sample_merge_type);

            % Update distance matrix and the gram matrix
            [distance_matrix, gram_matrix] = update_distance_matrix(distance_matrix, gram_matrix, gram_vector, new_train_sample_norm, merged_sample_id, -1, prior_weights(merged_sample_id), params.learning_rate);

            % Update the prior weight of the merged sample
            prior_weights(closest_sample_to_new_sample) = prior_weights(closest_sample_to_new_sample) + params.learning_rate;
        else
            % If the min distance amongst any of the existing
            % samples is less than the min distance of the new sample to the existing samples,
            % we merge the nearest existing samples and insert the new
            % sample in the vacated position

            % Renormalize prior weights
            prior_weights = prior_weights*(1 - params.learning_rate);

            % Ensure that the sample with higher prior weight is assigned
            % id1. This is only relevant if the sample_update_criteria is
            % 'replace'
            if prior_weights(closest_existing_sample2) > prior_weights(closest_existing_sample1)
                temp = closest_existing_sample1;
                closest_existing_sample1 = closest_existing_sample2;
                closest_existing_sample2 = temp;
            end

            sample_to_merge1 = cell(1, 1, num_feature_blocks);
            sample_to_merge2 = cell(1, 1, num_feature_blocks);

            % Extract the old sample
            for k = 1:num_feature_blocks
                sample_to_merge1{k} = samplesf{k}(closest_existing_sample1,:,:,:);
                sample_to_merge2{k} = samplesf{k}(closest_existing_sample2,:,:,:);
            end

            % Merge the existing closest samples
            merged_sample = merge_samples(sample_to_merge1, sample_to_merge2, prior_weights(closest_existing_sample1), prior_weights(closest_existing_sample2), params.sample_merge_type);

            % Update distance matrix and the gram matrix
            [distance_matrix, gram_matrix] = update_distance_matrix(distance_matrix, gram_matrix, gram_vector, new_train_sample_norm, closest_existing_sample1,closest_existing_sample2, ...
                prior_weights(closest_existing_sample1),prior_weights(closest_existing_sample2));

            % Update prior weights for the merged sample and the new sample
            prior_weights(closest_existing_sample1) = prior_weights(closest_existing_sample1) + prior_weights(closest_existing_sample2);
            prior_weights(closest_existing_sample2) = params.learning_rate;

            % Set the merged sample position and new sample position
            merged_sample_id = closest_existing_sample1;
            new_sample_id = closest_existing_sample2;

            new_sample = new_train_sample;
        end
    end
else
    % If the memory is not full, insert the new sample in the next empty
    % location

    sample_position = num_training_samples + 1;

    % Update the distance matrix and the gram matrix
    [distance_matrix, gram_matrix] = update_distance_matrix(distance_matrix, gram_matrix, gram_vector, new_train_sample_norm, sample_position, -1, 0, 1);

    % Update the prior weight
    if sample_position == 1
        prior_weights(sample_position) = 1;
    else
        prior_weights = prior_weights*(1 - params.learning_rate);
        prior_weights(sample_position) = params.learning_rate;
    end

    new_sample_id = sample_position;
    new_sample = new_train_sample;
end

% Ensure that prior weights always sum to 1
if (abs(1 - sum(prior_weights)) > 1e-5)
    error('Weights not properly updated');
end
end

author: martin-danelljan

repo: ECO


text:
function convnfft_install
% function convnfft_install
% Installation by building the C-mex file needed for convnfft
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% History
%  Original: 16/Sept/2009

arch=computer('arch');
mexopts = {'-O' '-v' ['-' arch]};
% 64-bit platform
if ~isempty(strfind(computer(),'64'))
    mexopts(end+1) = {'-largeArrayDims'};
end

% invoke MEX compilation tool
mex(mexopts{:},'inplaceprod.c');

author: thomaskuestner


text:
function A1m1 = A1m1f(epsi)
%A1M1F  Evaluate A_1 - 1
%
%   A1M1 = A1M1F(EPSI) evaluates A_1 - 1 using Eq. (17).  EPSI and A1M1 are
%   K x 1 arrays.

  eps2 = epsi.^2;
  t = eps2.*(eps2.*(eps2+4)+64)/256;
  A1m1 = (t + epsi) ./ (1 - epsi);
end

author: Sable

repo: mcbench-benchmarks


text:
function [periods,quiescence] = QuietPeriods(v,velocity,duration,brief)

%QuietPeriods - Find periods of immobility.
%
%  Find periods of immobility, i.e. periods of sufficient duration
%  where instantaneous linear velocity remains low. Brief movements
%  can be ignored.
%
%  USAGE
%
%    [periods,quiescence] = QuietPeriods(v,velocity,duration,brief)
%
%    v              linear velocity samples [t v]
%    velocity       maximum velocity of a quiet period
%    duration       minimum duration of a quiet period
%    brief          optional maximum duration of a 'brief' movement
%
%  OUTPUT
%
%    periods        list of [start stop] pairs
%    quiescence     list of [t s] pairs, where s is 1 if the animal
%                   is quiet at time t (and 0 otherwise)
%
%  SEE
%
%    See also BrainStates, PlotIntervals.

% Copyright (C) 2008-2011 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

if nargin < 3,
	error('Incorrect number of parameters (type ''help QuietPeriods'' for details).');
end
if nargin == 3,
	brief = 0;
end

% Determine beginning/end of quiet periods
below = v(:,2) < velocity;
crossings = diff(below); % yields -1 for upward crossings, and 1 for downward crossings
start = find(crossings == 1);
stop = find(crossings == -1);

% The previous code would ignore quiet periods beginning at the first sample, or ending at the last sample; correct for this
if below(1),
	start = [1;start];
end
if below(end),
	stop = [stop;length(below)];
end

% Determine durations of movements, and discard brief ones
durations = v(start(2:end),1) - v(stop(1:end-1),1);
ignore = find(durations <= brief);
start(ignore+1) = [];
stop(ignore) = [];

% Keep only long enough periods
durations = v(stop,1)-v(start,1);
discard = durations < duration;
start(discard) = [];
stop(discard) = [];

% Outputs
periods = [v(start,1) v(stop,1)];
quiescence = zeros(size(v,1),2);
quiescence(:,1) = v(:,1);
for i = 1:length(start),
	quiescence(start(i):stop(i),2) = 1;
end


text:
function dist = sll_distanceeuclcv(data,condid,varargin)

%preliminary code, undocumented; please do not distribute
%data is 3D (variables, timepoints,observations)
%condid is a cell of strings denoting the condition id for each observation
%
%temporalgen output: dsvm.d: (test time x train time x conditions)
%multiclass output: dsvm.d: (actual condition x predicted condition x conditions)
%
% example:
%  d = sll_distanceeuclcv(data,condid,'numpermutation',10,'verbose',2);
%
% Author: Dimitrios Pantazis

%% parse inputs

numpermutation  = sll_inputparser(varargin,'numpermutation',100, @(x) isscalar(x) && x>0 && x == round(x) );
kfold           = sll_inputparser(varargin,'kfold',5, @(x) isscalar(x) && x>=2 && x == round(x));
method          = sll_inputparser(varargin,'method','pairwise',{'pairwise','temporalgen'});
whiten          = sll_inputparser(varargin,'whiten',true);
verbose         = sll_inputparser(varargin,'verbose',false);

%% check inputs and convert string labels to numbers

[condlabel,condidval,conditiongen,data2,condlabel2,condidval2] = sll_checkdatacondid(data,condid,varargin{:});

%% initialize variables

numcond = length(condlabel); %number of conditions
rng('shuffle'); %seed the random number generator based on the current time
d = 0; %store decoding results
d2 = 0; %store decoding results (if cross-condition generalized)

%% compute cross-validated Euclidean distance
for p = 1:numpermutation

    %% verbose permutation

    if verbose & ~rem(p,verbose)
        disp(['Permutation: ' num2str(p) ' out of ' num2str(numpermutation)]);
    end

    %% assign data to k folds

    foldid = sll_createfolds(condidval,kfold);
    for c = 1:numcond
        traindata{c} = mean(data(:,:,condidval==c & foldid<kfold),3);
        testdata{c} = mean(data(:,:,condidval==c & foldid==kfold),3);
        if conditiongen
            n = nnz(condidval==c & foldid==kfold); %we need that many elements from data2
            cndx = find(condidval2==c); %available elements in condition c
            testdata2{c} = mean(data2(:,:,cndx(randperm(length(cndx),n))),3);
        end
    end

    %% multivariate noise normalization (whiten data)

    if whiten
        %compute whitening matrix using only training data
        down = max(round(size(data,2)/100),1); %downsample the data to make computation fast
        W = sll_whitenmatrix(data(:,1:down:end,foldid<kfold),condidval(foldid<kfold));
        %apply whitening matrix
        for c = 1:numcond
            traindata{c} = W*traindata{c}; %multiply data with whitener
            testdata{c} = W*testdata{c}; %multiply data with whitener
            if conditiongen
                testdata2{c} = W*testdata2{c};
            end
        end
    end

    %% Euclidean cross-validated distance

    d = d + distanceeuclcv(traindata,testdata,method);
    if conditiongen
        d2 = d2 + distanceeuclcv(traindata,testdata2,method);
    end

end

%% normalize output

d = d/numpermutation; %divide by number of permutations
d2 = d2/numpermutation; %divide by number of permutations

%% parse output

dist.d                 = d;
dist.condlabel         = condlabel;
dist.numpermutation    = numpermutation;
dist.kfold             = kfold;
dist.whiten            = whiten;
dist.method            = method;
if conditiongen
    dist.d2            = d2;
    dist.condlabel2    = condlabel2;
end

%% local function; computes Euclidean cross-validated distance
function d = distanceeuclcv(traindata,testdata,method)

%initial variables
ncond = length(traindata);
ntimes = size(traindata{1},2);
[x,y] = find(tril(ones(ncond,ncond),-1)); %indices of lower triangle

switch method

    case 'pairwise'

        %pairwise Euclidean cross-validated distance
        d = zeros(ntimes,ncond*(ncond-1)/2,'single');
        for i = 1:length(x)
            d(:,i) = sum( (traindata{x(i)}-traindata{y(i)}).*(testdata{x(i)}-testdata{y(i)}))'; %cross-validated Euclidean distance
        end

    case 'temporalgen'

        %pairwise Euclidean cross-valideated with temporal generalization
        d = zeros(ntimes,ntimes,ncond*(ncond-1)/2,'single');
        for i = 1:length(x)
            d(:,:,i) = (traindata{x(i)}-traindata{y(i)})' * (testdata{x(i)}-testdata{y(i)}); %cross-validated Euclidean distance
        end

end

author: brainstorm-tools


text:
nProj = size(uproj,2);

uBase = zeros(1e4, nProj);
nS = zeros(1e4, 1);
ncurr = 1;

for ibatch = 1:numel(indBatch)
    % merge in with existing templates
    uS = uproj(indBatch{ibatch}, :);
    [nSnew, iNonMatch] = merge_spikes0(uBase(1:ncurr,:), nS(1:ncurr), uS, ops.crit);
    nS(1:ncurr) = nSnew;
    %
    % reduce non-matches
    [uNew, nSadd] = reduce_clusters0(uS(iNonMatch,:), ops.crit);

    % add new spikes to list
    uBase(ncurr + [1:size(uNew,1)], :) = uNew;
    nS(ncurr + [1:size(uNew,1)]) = nSadd;

    ncurr = ncurr + size(uNew,1);

    if ncurr>1e4
        break;
    end
end
%
nS = nS(1:ncurr);
uBase = uBase(1:ncurr, :);

[~, itsort] = sort(nS, 'descend');

%% initialize U
% compute covariance matrix
sigDrift = 15;
sigShift = 20;
chanDists= bsxfun(@minus, rez.yc, rez.yc').^2 + bsxfun(@minus, rez.xc, rez.xc').^2;
iCovChans = my_inv(exp(-chanDists/(2*sigDrift^2)), 1e-6);

Nfilt = ops.Nfilt;
lam = ops.lam(1) * ones(Nfilt, 1, 'single');

U = gpuArray(uBase(itsort(1:Nfilt), :))';
mu = sum(U.^2,1)'.^.5;
U = normc(U);
%
deltay = zeros(ops.Nchan, numel(indBatch));

for i = 1:10
    % resample spatial masks up and down
    Uup   = shift_data(reshape(U, ops.Nchan, []),  sigShift, rez.yc, rez.xc, iCovChans, sigDrift, rez.Wrot);
    Uup   = reshape(Uup, size(U));
    Udown = shift_data(reshape(U, ops.Nchan, []), -sigShift, rez.yc, rez.xc, iCovChans, sigDrift, rez.Wrot);
    Udown = reshape(Udown, size(U));

    mu = repmat(mu, 3, 1);
    lam = repmat(lam, 3, 1);

    dWU = zeros(Nfilt, nProj, 'single');
    nToT = gpuArray.zeros(Nfilt, 1, 'single');
    Cost = gpuArray(single(0));
    %
    for ibatch = 1:numel(indBatch)
        % find clusters
        clips = gpuArray(uproj(indBatch{ibatch}, :))';
        nSpikes = size(clips,2);
        clips = reshape(clips, ops.Nchan, []);

        % resample clips by the delta y
        clips = shift_data(clips, deltay(:, ibatch), rez.yc, rez.xc, iCovChans, sigDrift, rez.Wrot);
        clips = reshape(clips, size(U,1), [])';

        ci = clips * [Udown U Uup];

        ci = bsxfun(@plus, ci, (mu .* lam)');
        cf = bsxfun(@rdivide, ci.^2, 1 + lam');
        cf = bsxfun(@minus, cf, (mu.^2.*lam)');

        cf = reshape(cf,[], Nfilt, 3);
        [Mmax, imax] = max(cf, [], 2);
        [~, i3max] = max(Mmax, [], 3);
%         keyboard;

        % determine added correction to delta y

        % determine cluster assignment for this iteration
        [max_cf, id] = max(cf(:,:,2), [], 2);
        id = gather_try(id);
        L = gpuArray.zeros(Nfilt, nSpikes, 'single');
        L(id' + [0:Nfilt:(Nfilt*nSpikes-1)]) = 1;
        dWU = dWU + L * clips;
        nToT = nToT + sum(L, 2);
        Cost = Cost + mean(max_cf);
    end
    dWU  = bsxfun(@rdivide, dWU, nToT);

    U = dWU';
    mu = sum(U.^2,1)'.^.5;
    U = normc(U);
    Cost = Cost/size(inds,2);

    % smooth out corrections

    %     disp(Cost)

    %     plot(sort(log(1+nToT)))
    %     drawnow
end
%%
Nchan = ops.Nchan;
Nfilt = ops.Nfilt;
wPCA = ops.wPCA(:,1:3);
Urec = reshape(U, Nchan, size(wPCA,2), Nfilt);

nt0 = 61;
Urec= permute(Urec, [2 1 3]);
Wrec = reshape(wPCA * Urec(:,:), nt0, Nchan, Nfilt);

Wrec = gather_try(Wrec);
Nrank = 3;
W = zeros(nt0, Nfilt, Nrank, 'single');
U = zeros(Nchan, Nfilt, Nrank, 'single');
for j = 1:Nfilt
    [w sv u] = svd(Wrec(:,:,j));
    w = w * sv;

    Sv = diag(sv);
    W(:,j,:) = w(:, 1:Nrank)/sum(Sv(1:ops.Nrank).^2).^.5;
    U(:,j,:) = u(:, 1:Nrank);
end

Uinit = U;
Winit = W;
mu = gather_try(single(mu));
muinit = mu;

WUinit = zeros(nt0, Nchan, Nfilt);
for j = 1:Nfilt
    WUinit(:,:,j) = muinit(j)  * Wrec(:,:,j);
end
WUinit = single(WUinit);
%%

author: cortex-lab

repo: KiloSort

save_path: github-repos/MATLAB/cortex-lab-KiloSort


text:
function [ add, mul, sub ] = setfld2 ( dummy )

%*****************************************************************************80
%
%% SETFLD2 sets up arithmetic tables for the finite field of order 2.
%
%  Discussion:
%
%    SETFLD2 sets up addition, multiplication, and subtraction tables
%    for the finite field of order 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2003
%
%  Author:
%
%    Original FORTRAN77 version by Paul Bratley, Bennett Fox, Harald Niederreiter.
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, integer DUMMY, a dummy argument.
%
%    Output, integer ADD(2,2), MUL(2,2), SUB(2,2), the addition,
%    multiplication, and subtraction tables, mod 2.
%
  q = 2;

  p = 2;

  for i = 0 : 1
    for j = 0 : 1
      add(i+1,j+1) = mod ( i + j, p );
    end
  end

  for i = 0 : 1
    for j = 0 : 1
      mul(i+1,j+1) = mod ( i * j, p );
    end
  end

  for i = 0 : 1
    for j = 0 : 1
      sub(add(i+1,j+1)+1, i+1) = j;
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [FV,C] = eeg_interp_sph_spline(Zi,Ei)

% eeg_interp_sph_spline - Spherical Spline Interpolation of Potential
%
% Useage: [FV,C] = eeg_interp_sph_spline(Zi,Ei)
%
% where:    Zi is Nelec x 1, an EEG/ERP measurement at time t
%           Ei is Nelec x 3, [X Y Z] electrode positions.
%           The origin of Ei is assumed (0,0,0).
%
% FV => interpolated spherical surface (see sphere_tri)
%
% FV.faces    => triangulation of FV.vertices
% FV.vertices => cartesian coordinates (Nx3)
% FV.Cdata    => spherical spline potential at FV.vertices
%
% C => interpolation coefficients of Ei (includes co = C(1))
%
% Notes:    This function calculates the spherical spline of
%           Perrin et al (1989).  Electroenceph. & Clin.
%             Neurophysiology, 72: 184-187. Corrigenda (1990),
%             Electroenceph. & Clin. Neurophysiology, 76: 565.
%             (see comments in the .m file for details).

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:51 $

% Licence:  GNU GPL, no implied or express warranties
% History:  08/2001 Darren.Weber_at_radiology.ucsf.edu, with
%                   mathematical advice from
%                   Dr. Murk Bottema (Flinders University of SA)
%           10/2003 Darren.Weber_at_radiology.ucsf.edu, with
%                   mathematical advice and LegendreP function from
%                   Dr. Tom.Ferree_at_radiology.ucsf.edu
%                   revised, tested & initial verification complete
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check for correct size & orientation of Ei & Zi
[e1,e2] = size(Ei);
[v1,v2] = size(Zi);
if e1 < e2, Ei = Ei'; [e1,e2] = size(Ei); end
if v1 < v2, Zi = Zi'; [v1,v2] = size(Zi); end
if ~and(isequal(e1,v1),and(isequal(e2,3),isequal(v2,1)))
  error('...Ei must be Nx3 & Zi must be Nx1');
end
nElectrodes = e1; % The number of electrodes
clear e1 e2 v1 v2;

% -------------------------------------------------------------------------
% estimate spherical radius of the electrodes and
% obtain spherical projections of Ei
[r,x,y,z] = elec_sphere_project(Ei(:,1),Ei(:,2),Ei(:,3));
%Ei = [ x y z ]; clear x y z;

% create spherical interpolation surface
FV = sphere_tri('ico',4,r);

% -------------------------------------------------------------------------
% Calculate the cosines, if Ei is Nx3, COS is NxN matrix
% We use (Ei,Ei) here because it gives the cosines between
% each electrode and every other electrode.  This is required
% here because we solve a linear system of equations below
% that will find the interpolated value at a given electrode
% location, which must be equal to the measured
% potential at that location.
EiCOS = elec_cosines(Ei,Ei);

% create zeros on the diagonal elements
% [not sure why this works, but it does.]
for i = 1:length(EiCOS), EiCOS(i,i) = 0; end

% -------------------------------------------------------------------------
% Calculate g(x), nElectrodes x nElectrodes
Gx = eeg_interp_sph_spline_g(EiCOS);

% -------------------------------------------------------------------------
% calculate the spherical interpolation coefficients
C = eeg_interp_sph_spline_c(Zi,Gx); clear Gx

% -------------------------------------------------------------------------
% cosines between electrodes and interpolation points
FvCOS = elec_cosines(Ei,FV.vertices);

% Calculate g(x), nElectrodes x NinterpolationPoints
Gx = eeg_interp_sph_spline_g(FvCOS);

eegversion = '$Revision: 1.1 $';
fprintf('EEG_INTERP_SPH_SPLINE [v %s]\n',eegversion(11:15)); tic

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now that Ci is solved, we can obtain interpolated potentials at Ej (Eq. 1)
% U(Ej) = c(0) + ( for i=1:n, sum = (sum + (c(i) * g(cos(Ei,Ej)))) )
% U(Ej) = c(0) + sum( Ci * g(x) )

% Solve Eq 1. (where FV.Cdata = U)
Co = C(1);
Ci = C(2:end);
FV.Cdata = Co + ( Ci' * Gx );

t=toc; fprintf('...done (%6.2f sec)\n',t);

return

author: PatternRecognition

repo: OpenBMI


text:
function mex_files_to_compile_map = PTKGetMexFilesToCompile(reporting)
    % PTKGetMexFilesToCompile. Returns a list of mex files used by PTK
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    root_dir = PTKDirectories.GetSourceDirectory;
    mex_dir = PTKDirectories.GetMexSourceDirectory;

    % Populate list with known mex files
    mex_files_to_compile = CoreCompiledFileInfo.empty(0);
    mex_files_to_compile(end + 1) = CoreCompiledFileInfo(2, 'PTKFastEigenvalues', 'cpp', mex_dir, [], []);
    mex_files_to_compile(end + 1) = CoreCompiledFileInfo(1, 'PTKFastIsSimplePoint', 'cpp', mex_dir, [], []);
    mex_files_to_compile(end + 1) = CoreCompiledFileInfo(2, 'PTKWatershedFromStartingPoints', 'cpp', mex_dir, [], []);
    mex_files_to_compile(end + 1) = CoreCompiledFileInfo(3, 'PTKWatershedMeyerFromStartingPoints', 'cpp', mex_dir, [], []);
    mex_files_to_compile(end + 1) = CoreCompiledFileInfo(2, 'PTKSmoothedRegionGrowingFromBorderedImage', 'cpp', mex_dir, [], []);

    mex_files_to_compile(end + 1) = CoreCompiledFileInfo(3, 'mba_surface_interpolation', 'cpp', fullfile(root_dir, 'External', 'gerardus', 'matlab', 'PointsToolbox'), ...
        {['-I' fullfile(root_dir, 'External')], ['-I' fullfile(root_dir, 'External', 'mba', 'include')]}, ...
        {fullfile(root_dir, 'External', 'mba', 'src', 'MBA.cpp'), fullfile(root_dir, 'External', 'mba', 'src', 'UCBsplines.cpp'), ...
        fullfile(root_dir, 'External', 'mba', 'src', 'UCBsplineSurface.cpp'), fullfile(root_dir, 'External', 'mba', 'src', 'MBAdata.cpp')});

    % Transfer to a map
    mex_files_to_compile_map = containers.Map;
    for mex_file = mex_files_to_compile
        mex_files_to_compile_map(mex_file.Name) = mex_file;
    end

    % Now add unknown .cpp files
    file_list = dir(fullfile(mex_dir, '*.cpp'));
    for index = 1 : numel(file_list)
        filename = file_list(index).name;
        [~, name_part, ~] = fileparts(filename);
        if ~mex_files_to_compile_map.isKey(name_part)
            mex_files_to_compile_map(name_part) = CoreCompiledFileInfo(0, name_part, 'cpp', mex_dir, [], []);
            reporting.ShowWarning('PTKGetMexFilesToCompile:MexFileNotInList', ['The file ' filename ' was found in the mex folder but has not been added to PTKGetMexFilesToCompile.m. Mex files should be added to PTKGetMexFilesToCompile.m so they are correctly versioned.'], []);
        else
        end
    end

end

author: tomdoel

repo: pulmonarytoolkit


text:
function X=NIG(th,k,s,ts,J)

T=length(ts);
DXs=zeros(J,1);
for t=1:T
    Dt=ts(1)-0;
    if t>1
        Dt=ts(t)-ts(t-1);
    end
    l=1/k*(Dt^2);
    m=Dt;
    DS=IG(l,m,J);
    N=randn(J,1);

    DX=s*N.*sqrt(DS)+th*DS;
    DXs=[DXs DX];
end
X=cumsum(DXs,2);


text:
% Calculate the number of star motifs of given (subgraph) size.
% Note 1: Easily extendible to return the actual stars as k-tuples of nodes.
% Note 2: Star of size 1 is the trivial case of a single node.
%
% INPUTs: adjacency list {} (1xn), k - star motif size
% OUTPUTs: number of stars with k nodes (k-1 spokes)
%
% GB: last updated, Oct 5, 2012

function num = numStarMotifs(adjL,k)

num = 0;

for i=1:length(adjL)
  if length(adjL{i})>=(k-1); num = num + nchoosek(length(adjL{i}),k-1); end
end

% ALTERNATIVE
% function num = numStarMotifs(adj,k)
%
% % INPUTs: adjacency matrix, k - star motif size
% % OUTPUTs: number of stars with k nodes (k-1 spokes)
%
% [deg,~,~]=degrees(adj);
%
% num=0;
%
% for i=1:length(deg)
%     if deg(i)>=(k-1); num=num+nchoosek(deg(i),k-1); end
% end

author: aeolianine

repo: octave-networks-toolbox


text:
% Load, modify and export a fig file as an eps file.

clear all;
addpath('../lib');

%% lets plot 3 cycles of 50Hz AC voltage
f = 50;  % frequency
Vm = 10; % peak
phi = 0; % phase

% generate the signal
t = [0:0.0001:3/f];
th = 2*pi*f*t;
v = Vm*sin(th+phi);

% plot it
plt = Plot(t*1E3, v);

plt.Title = 'Voltage as a function of time'; % plot title
plt.XLabel = 'Time, t (ms)'; % xlabel
plt.YLabel = 'Voltage, V (V)'; %ylabel

plt.export('plotSimple1.png');

author: masumhabib

repo: PlotPub

save_path: github-repos/MATLAB/masumhabib-PlotPub


text:
function [hh] = plotEmissionParams( varargin )

MAX_K   = 12;
MIN_THR = 0.01;

% ================================================ Process User Input
if isstruct( varargin{1} )
    if isfield( varargin{1}, 'ThetaM' )
        theta = varargin{1}.ThetaM.theta;
    elseif isfield( varargin{1}, 'theta' ) % || isprop(varargin{1},'theta')
        theta = varargin{1}.theta;
    else
        theta = varargin{1};
    end
    if isfield( varargin{1}, 'stateSeq' )
        stateSeq = varargin{1}.stateSeq;
    end

    for ll = 2:length( varargin )
       curArg = varargin{ll};
       if ishandle(curArg)
           figH = curArg;
       elseif isobject(curArg);
           data = curArg;
       end
    end

elseif isobject( varargin{1} )
    if isprop( varargin{1}, 'theta' )
        theta = varargin{1}.theta;
    end
else
    jobID = varargin{1};
    taskID = varargin{2};
    DATA = loadSamplerOutput( jobID, taskID, {'iters', 'Psi'} );
    data = loadSamplerInfo( jobID, taskID, {'data'} );

    if length( varargin ) >= 3
        queryIter = varargin{3};
        [~, idx] = min( abs( queryIter - DATA.iters.Psi ) );
    else
        idx = length( DATA.Psi );
    end
    Psi = DATA.Psi( idx );
    theta = Psi.theta;
    stateSeq = Psi.stateSeq;

    if length( varargin ) >= 4
        objIDs = varargin{4};
    end
end

K = length( theta );

if ~exist('figH','var')
    figH = gca;
end

if exist( 'stateSeq', 'var' )
    Zall = horzcat( stateSeq(:).z );
    N = length( Zall );
    [Zcounts, sortIDs] = sort( histc( Zall, 1:K ), 'descend' );
    sortIDs = sortIDs( Zcounts > 0 );
    K = length( sortIDs );
    K = min( K, MAX_K );
    theta = theta( sortIDs(1:K) );
else

    K = min( K, MAX_K );
    theta = theta( 1:K );
end

fNames = fieldnames(theta);
if fNames{1} == 'mu'
    obsType = 'Gaussian';
elseif fNames{1} == 'A'
    obsType = 'AR-Gaussian';
elseif fNames{1} == 'logp'
    obsType = 'Multinomial';
end

switch obsType
    case 'Multinomial'

        Px = exp( vertcat( theta(:).logp ) );
        set( gcf, 'Units', 'normalized', 'Position', [0 0.5 0.25 0.5] );
        imagesc( figH, Px, [0 0.01] );
        set( gca, 'Position', [0.1 0.1 0.8 0.06*K] );
    case 'Gaussian'

        if exist( 'data', 'var' )
            X = data.Xdata;
            plot(figH, X(1,:), X(2,:), 'k.' );
        end

        for kk = 1:K
            Mu(kk,:) = theta(kk).mu;
            invSigma(:,:,kk) = theta(kk).invSigma;
        end
        plotGauss( Mu, invSigma, 1, jet( MAX_K ), figH );

        B = 1.6;
        axis( [-B B -B B] );

    case 'AR-Gaussian'
        D = size( theta(1).A, 1);
        for kk = 1:K
            Mu(kk,:) = [ mean( diag( theta(kk).A ) )  zeros( 1, D-1) ];
            invSigma(:,:,kk) = theta(kk).invSigma;
        end
        plotGauss( Mu, 1000*invSigma, 1, jet(MAX_K), figH );

    case 'Bernoulli'
        error( 'TO DO' );
    case 'Multinomial'
        error( 'TO DO' );
end

hold off;
drawnow;

author: michaelchughes

repo: NPBayesHMM


text:
function [ H, Lh ] = homographyHarker( DataA, DataB, LA, LB)
%
% Purpose : Computes the general projective transformation between two sets
% of 2D data using a linear algorithm (the homography).
%
% Uses (syntax) :
%   H = homography( DataA, DataB )
%
% Input Parameters :
%   DataA, DataB := 2D homogeneous data sets in matrix form (3xn)
%
% Return Parameters :
%   H := the 3x3 homography
%
% Description and algorithms:
%   The algorithm is based on the Direct Linear Transform (DLT) method
%   outlined in Hartley et al.  The method uses orthogonal projections of
%   matrices, such that the vanishing line is treated as the principal
%   component of the reduction.  In this manner, the statistical behaviour
%   of the errors in variables are treated uniformly, see Harker and
%   O'Leary 2005.
%
% References :
%   Harker, M., O'Leary, P., Computation of Homographies, to appear in
%   Proceedings of the British Machine Vision Conference 2005, Oxford,
%   England.
%   Hartley, R., Zisserman, A., Multiple View Geometry in Computer Vision,
%   Cambridge University Press, Cambridge, 2001
%
% Cite this as :
%
% Author : Matthew Harker
% Date : July 25, 2005
% Version : 1.0
%--------------------------------------------------------------------------
% (c) 2005, O'Leary, Harker, University of Leoben, Leoben, Austria
% email: automation@unileoben.ac.at, url: automation.unileoben.ac.at
%--------------------------------------------------------------------------
% History:
%   Date:           Comment:
%   July 25, 2005   Original Version 1.0
%--------------------------------------------------------------------------
%
% Check input parameters:
%
if nargin == 2
    %
    if ~is2DData( DataA ) | ~is2DData( DataB )
        %
        error('Input does not represent 2D data') ;
        %
    end
    %
    [mA,nA] = size( DataA ) ;
    [mB,nB] = size( DataB ) ;
    %
    if (nA ~= nB)
        %
        error('Data sets must be the same size') ;
        %
    end
    %
    %
elseif nargin == 4
    %
    if ~is2dcovdata( DataA, LA ) | ~is2dcovdata( DataB, LB )
        %
        error('Input data does not match') ;
        %
    end
    %
    [mA,nA] = size( DataA ) ;
    [mB,nB] = size( DataB ) ;
    %
    if (nA ~= nB)
        %
        error('Data sets must be the same size') ;
        %
    end
    %
    %
    error('Error propagation not implemented as of yet') ;
    %
else
    %
    error('Incorrect number of input parameters') ;
    %
end
%
% Normalize the input data:
%
if nargin == 2
    %
    [DataA,TA,TAi] = normalizeData( DataA ) ;
    [DataB,TB,TBi] = normalizeData( DataB ) ;
    %
elseif nargin == 4
    %
    [DataA,TA,TAi,LA] = normalizeData( DataA, LA ) ;
    [DataB,TB,TBi,LB] = normalizeData( DataB, LB ) ;
    %
end
%
% Construct the orthogonalized design matrix :
%
C1 = -DataB(1,:) .* DataA(1,:) ;
C2 = -DataB(1,:) .* DataA(2,:) ;
%
C3 = -DataB(2,:) .* DataA(1,:) ;
C4 = -DataB(2,:) .* DataA(2,:) ;
%
mC1 = mean( C1 ) ;
mC2 = mean( C2 ) ;
mC3 = mean( C3 ) ;
mC4 = mean( C4 ) ;
%
Mx = [ C1' - mC1, C2' - mC2, -DataB(1,:)' ] ;
My = [ C3' - mC3, C4' - mC4, -DataB(2,:)' ] ;
%
Pp = pinv(DataA(1:2,:)') ;
%
Bx = Pp * Mx ;
By = Pp * My ;
%
D = [ Mx - DataA(1:2,:)'*Bx ;...
      My - DataA(1:2,:)'*By ] ;
%
% Find v_min and backsubstitute :
%
[U,S,V] = svd( D, 0 ) ;
%
h789 = V(:,end) ;
h12 = -Bx * h789 ;
h45 = -By * h789 ;
h3 = -[mC1, mC2] * h789(1:2) ;
h6 = -[mC3, mC4] * h789(1:2) ;
%
% Reshape vector h to matrix H, and transform :
%
H = reshape( [h12; h3; h45; h6; h789], 3, 3)' ;
%
H = TB * H * TAi ;
%
% Compute the covariance of the output :
%
%
%
% END
%

function [DataN, T, Ti, LN] = normalizeData( Data, L ) ;
%
% Purpose : Computes a set of data corresponding to the input data with its
% centroid subtracted, and scaled such that the root-mean-square distance
% to the origin is sqrt(2).  The transformation T, carries the scaled data
% back to its original form.  Optionally, the first order estimation of
% covariance matrices are computed.
%
% Uses (syntax) :
%   [DataN, T, LN] = normalizeData( Data, L )
%   [DataN, T] = normalizeData( Data )
%
% Input Parameters :
%   Data := a 3xn matrix of homogeneous points.
%   L    := is a 3x3 covariance matrix (all points have identical covariance), or
%           a 3x3xn array of n covariance matrices.
%
% Return Parameters :
%   DataN := mean-free data scaled s.t. d_RMS = sqrt(2)
%   T     := transformation to bring DataN to the Affine coordinates
%            corresponding to Data (NOTE: T*DataN is in affine coords).
%   LN    := the covariance of the scaled normalized data (size is
%            generally 2x2xn, due to the normalization)
%
% Description and algorithms:
%
% References :
%   Clarke, J.C., Modelling Uncertainty: A Primer, Dept. of Engineering
%   Science, Oxford University, Technical Report.
%
% Cite this as :
%
% Author : Matthew Harker
% Date : July 7, 2005
% Version :
%
% (c) 2005, Institute for Automation, University of Leoben, Leoben, Austria
% email: automation@unileoben.ac.at, url: automation.unileoben.ac.at
%
% History:
%   Date:           Comment:
%                   Original Version 1.0
%--------------------------------------------------------------------------
%
% Check input arguments :
%
if nargin == 1
    %
    if ~is2DData( Data )
        %
        error('Input does not represent 2D data') ;
        %
    end
    %
elseif nargin == 2
    %
    if ~is2dcovdata( Data, L )
        %
        error('Input data and covariances do not match') ;
        %
    end
    %
else
    %
    error('Not enough input arguments') ;
    %
end
%
s = Data(1,:) ;
t = Data(2,:) ;
u = Data(3,:) ;
%
x = s./u ;
y = t./u ;
%
xm = mean( x ) ;
ym = mean( y ) ;
%
xh = x - xm ;
yh = y - ym ;
%
n = length( xh ) ;
%
kappa = sum( xh.^2 + yh.^2 ) ;
%
beta = sqrt( 2*n / kappa ) ;
%
xn = beta * xh ;
yn = beta * yh ;
%
DataN = [ xn; yn; ones(size(xn)) ] ;
%
T = [ 1/beta,   0   , xm ;...
        0   , 1/beta, ym ;...
        0   ,   0   ,  1 ] ;
%
Ti = [ beta ,  0   , -beta * xm ; ...
        0   , beta , -beta * ym ; ...
        0   ,  0   ,       1    ] ;
%
%
%
if nargin == 1
    %
    % Just the data is needed.
    %
    if nargout == 4
        warning('No covariance matrix for the input data given.') ;
        LN = [] ;
    end
    %
    %
    %
elseif nargin == 2
    %
    % There is a 3x3 homogeneous covariance matrix given
    %
    if nargout == 4
        %
        % The 2x2 covariance matrix of the Euclidean output data is required
        %
        dim = length( size( L ) ) ;
        %
        % dim == 2 --> all covariance matrices are assumed the same
        % dim == 3 --> each point has a given covariance
        %
        % DEHOMOGENIZE
        %
        for k = 1:n
            %
            %
            % Jacobian for dehomogenizing:
            Jd = ( 1/u(k)^2 ) * [ u(k),  0  , -s(k) ;...
                                   0  , u(k), -t(k) ] ;
            %
            %
            if dim == 2
                LN(:,:,k) = Jd * L * Jd' ;
            else
                LN(:,:,k) = Jd * L(:,:,k) * Jd' ;
            end
            %
            % NOTE: If all data has the same covariance, after
            % dehomogenizing they will generally all have different
            % covariance matrices (as implemented here).  Only if the data
            % is already affine, and the covariance matrix is Euclidean
            % will the "normalized" points all have the same covariance.
            %
        end
        %
        % SUBTRACT MEAN, SCALE:
        %
        % Jacobian for subtracting mean:
        Jm = [ (1-(1/n)),   0       , -1 ,  0 ;...
                  0     , (1-(1/n)) ,  0 , -1 ] ;
        %
        % Covariance of the mean value:
        Lm = (1/n^2) * [ sum(LN(1,1,:)), sum(LN(1,2,:)) ;...
                         sum(LN(2,1,:)), sum(LN(2,2,:)) ] ;
        %
        %
        Lkappa = 0 ;
        %
        for k = 1:n
            %
            % Subtract Mean Covariance
            %
            LN(:,:,k) = Jm * [ LN(:,:,k) , zeros(2,2) ; zeros(2,2) , Lm ] * Jm' ;
            %
            % Compute Covariance of 'kappa':
            %
            Jk = 2 * [ xh(k), yh(k) ] ;
            Lkappa = Lkappa + Jk*LN(:,:,k)*Jk' ;
            %
        end
        %
        %
        for k = 1:n
            %
            % Jacobian for scaling the data:
            %
            a = kappa^(-3/2) ;
            b = kappa^(-1/2) ;
            %
            Js = sqrt(2*n) * [ (b - a*xh(k)^2) , -a*xh(k)*yh(k)  , -(a*xh(k))/2 ;...
                                -a*xh(k)*yh(k) , (b - a*yh(k)^2) , -(a*yh(k))/2 ] ;
            %
            LN(:,:,k) = Js * [ LN(:,:,k), [0;0] ; [0,0] , Lkappa ] * Js' ;
            %
            %
            %
        end
        %
        %
        %
    elseif nargout == 3
        %
        warning('No covariance matrix has been returned') ;
        %
    end
    %
    %
end
%
% END
%

author: tobycollins

repo: IPPE

save_path: github-repos/MATLAB/tobycollins-IPPE


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = INVERSEKINEMATIC_C8(robot, T)
%   Solves the inverse kinematic problem for the KUKA KR5 ARC robot
%   where:
%   robot stores the robot parameters.
%   T is an homogeneous transform that specifies the position/orientation
%   of the end effector.
%
%   A call to Q=INVERSEKINEMATIC__KUKA_KR5_ARCHW returns 8 possible solutions, thus,
%   Q is a 6x8 matrix where each column stores 6 feasible joint values.
%
%
%   Example code:
%
%   robot=load_robot('EPSON', 'x');
%   q = [0 0 0 0 0 0];
%   T = directkinematic(robot, q);
%   %Call the inversekinematic for this robot
%   qinv = inversekinematic(robot, T);
%   check that all of them are feasible solutions!
%   and every Ti equals T
%   for i=1:8,
%        Ti = directkinematic(robot, qinv(:,i))
%   end
%	See also DIRECTKINEMATIC.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function q = inversekinematic_C8(robot, T)

%initialize q,
%eight possible solutions are generally feasible
q=zeros(6,8);

% %Evaluate the parameters
% theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
L6=abs(d(6));

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Compute the position of the wrist, being W the Z component of the end effector's system
W = T(1:3,3);

% Pm: wrist position
Pm = [-Py Px Pz]' - L6*W;

%first joint, two possible solutions admited:
% if q(1) is a solution, then q(1) + pi is also a solution
q1=atan2(Pm(2), Pm(1));

%solve for q2
q2_1=solve_for_theta2(robot, [q1 0 0 0 0 0 0], Pm);
%the other possible solution is q1 + pi
q2_2=solve_for_theta2(robot, [q1+pi 0 0 0 0 0 0], Pm);

%solve for q3
q3_1=solve_for_theta3(robot, [q1 0 0 0 0 0 0], Pm);
%solver for q3 for both cases
q3_2=solve_for_theta3(robot, [q1+pi 0 0 0 0 0 0], Pm);

%the next matrix doubles each column. For each two columns, two different
%configurations for theta4, theta5 and theta6 will be computed. These
%configurations are generally referred as wrist up and wrist down solution
q = [q1         q1         q1        q1       q1+pi   q1+pi   q1+pi   q1+pi;
     q2_1(1)    q2_1(1)    q2_1(2)   q2_1(2)  q2_2(1) q2_2(1) q2_2(2) q2_2(2);
     q3_1(1)    q3_1(1)    q3_1(2)   q3_1(2)  q3_2(1) q3_2(1) q3_2(2) q3_2(2);
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0];

%leave only the real part of the solutions
q=real(q);

%Note that in this robot, the joint q3 has a non-simmetrical range. In this
%case, the joint ranges from 60 deg to -219 deg, thus, the typical normalizing
%step is avoided in this angle (the next line is commented). When solving
%for the orientation, the solutions are normalized to the [-pi, pi] range
%only for the theta4, theta5 and theta6 joints.

%normalize q to [-pi, pi]
q(1,:) = normalize(q(1,:));
q(2,:) = normalize(q(2,:));
% solve for the last three joints
% for any of the possible combinations (theta1, theta2, theta3)
for i=1:2:size(q,2),
    qtemp = solve_spherical_wrist(robot, q(:,i), T, 1,'geometric'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i)=qtemp;

    qtemp = solve_spherical_wrist(robot, q(:,i), T, -1, 'geometric'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i+1)=qtemp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for second joint theta2, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q2 = solve_for_theta2(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry
L2=abs(a(2));
L3=abs(d(4));
A2 = abs(a(3));

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%The inverse kinematic problem can be solved as in the IRB 140 (for example)

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

beta = atan2(-p1(2), p1(1));
gamma = real(acos((L2^2+r^2-L4^2)/(2*r*L2)));

%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta3
q2(1) = pi/2 - beta - gamma; %elbow up
q2(2) = pi/2 - beta + gamma; %elbow down

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for third joint theta3, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q3 = solve_for_theta3(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry
L2=abs(a(2));
L3=abs(d(4));

A2 = abs(a(3));

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%the angle phi is fixed
phi=acos((A2^2+L4^2-L3^2)/(2*A2*L4));

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

beta = real(acos((L2^2 + L4^2 - r^2)/(2*L2*L4)));

%return two possible solutions
%elbow up and elbow down solutions
%the order here is important
q3(1) = pi - phi - beta;
q3(2) = pi - phi + beta;

repo: ARTE


text:
function X = normrow(X)
if nargin < 1
    error('Not enough input arguments.'); end

C = size(X,2);
if C == 1
    X = X ./ abs(X);
else
    X = sqrt(ones./(sum(X.*X,2)'))'*ones(1,C).*X;
end

repo: REST


text:
function [ grid_index, grid_base ] = levels_index ( dim_num, level_max, rule, ...
  point_num )

%*****************************************************************************80
%
%% LEVELS_INDEX indexes a sparse grid.
%
%  Discussion:
%
%    The sparse grid is the logical sum of product grids with total LEVEL
%    between LEVEL_MIN and LEVEL_MAX.
%
%    The necessary dimensions of GRID_INDEX can be determined by
%    calling LEVELS_INDEX_SIZE first.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Input, integer RULE, the index of the rule.
%    1, "CC", Clenshaw Curtis Closed Fully Nested rule.
%    2, "F1", Fejer 1 Open Fully Nested rule.
%    3, "F2", Fejer 2 Open Fully Nested rule.
%    4, "GP", Gauss Patterson Open Fully Nested rule.
%    5, "GL", Gauss Legendre Open Weakly Nested rule.
%    6, "GH", Gauss Hermite Open Weakly Nested rule.
%    7, "LG", Gauss Laguerre Open Non Nested rule.
%
%    Input, integer POINT_NUM, the total number of points
%    in the grids.
%
%    Output, integer GRID_INDEX(DIM_NUM,POINT_NUM), a list of
%    point indices, representing a subset of the product grid of level
%    LEVEL_MAX, representing (exactly once) each point that will show up in a
%    sparse grid of level LEVEL_MAX.
%
%    Output, integer GRID_BASE(DIM_NUM,POINT_NUM), a list of
%    the orders of the rules associated with each point and dimension.
%
  if ( rule == 1 )
    [ grid_index, grid_base ] = levels_index_cfn ( dim_num, level_max, ...
      point_num );
  elseif ( 2 <= rule & rule <= 4 )
    [ grid_index, grid_base ] = levels_index_ofn ( dim_num, level_max, ...
      point_num );
  elseif ( 5 <= rule & rule <= 6 )
    [ grid_index, grid_base ] = levels_index_own ( dim_num, level_max, ...
      point_num );
  elseif ( 7 == rule )
    [ grid_index, grid_base ] = levels_index_onn ( dim_num, level_max, ...
      point_num );
  else
    grid_index = [];
    grid_base = [];
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LEVELS_INDEX - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized rule number = %d\n', rule );
    error ( 'LEVELS_INDEX - Fatal error!' );
  end

  return
end


text:
function varargout = hikmeans(varargin)
% VL_HIKMEANS  Hierachical integer K-means
%   [TREE,ASGN] = VL_HIKMEANS(DATA,K,NLEAVES) applies integer K-menas
%   recursively to cluster the data DATA, returing a structure TREE
%   representing the clusters and a vector ASGN with the data to
%   cluster assignments. The depth of the recursive partition is
%   computed so that at least NLEAVES are generated.
%
%   VL_HIKMEANS() is built on top of VL_IKMEANS() and requires the
%   data to be of class UINT8.
%
%   TREE is a structure representing the hierarchical clusters.  Each
%   node of the tree is also a structure with fields:
%
%   DEPTH::
%     Depth of the tree (only at the root node)
%
%   CENTERS::
%     K cluster centers
%
%   SUB::
%     Array of K node structures representing subtrees
%     (this field is missing at leaves).
%
%   ASGN is a matrix with one column per datum and height equal to the
%   depth of the tree. Each column encodes the branch of the tree that
%   correspond to each datum.
%
%   Example::
%     ASGN(:,7) = [1 5 3] means that the tree as depth equal to 3 and
%     that the datum X(:,7) corresponds to the branch
%     ROOT->SUB(1)->SUB(5)->SUB(3).
%
%   See also: VL_HIKMEANSPUSH(), VL_HIKMEANSHIST(), VL_IKMEANS(), VL_HELP().
[varargout{1:nargout}] = vl_hikmeans(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
%CODEGENERATOR.GENMFUNCORIOLIS Generate M-functions for Coriolis matrix
%
% cGen.genmfuncoriolis() generates a robot-specific M-function to compute
% the Coriolis matrix.
%
% Notes::
% - Is called by CodeGenerator.gencoriolis if cGen has active flag genmfun
% - The Coriolis matrix is stored row by row to avoid memory issues.
% - The generated M-function recombines the individual M-functions for each row.
% - Access to generated function is provided via subclass of SerialLink
%   whose class definition is stored in cGen.robjpath.
%
% Author::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.CodeGenerator, CodeGenerator.gencoriolis, CodeGenerator.geninertia.

% Copyright (C) 2012-2014, by Joern Malzahn
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

function [ ] = genmfuncoriolis( CGen )

%% Does robot class exist?
if ~exist(fullfile(CGen.robjpath,[CGen.getrobfname,'.m']),'file')
    CGen.logmsg([datestr(now),'\tCreating ',CGen.getrobfname,' m-constructor ']);
    CGen.createmconstructor;
    CGen.logmsg('\t%s\n',' done!');
end

%%
CGen.logmsg([datestr(now),'\tGenerating m-function for the Coriolis matrix row' ]);

[q, qd] = CGen.rob.gencoords;
nJoints = CGen.rob.n;

for kJoints = 1:nJoints
    CGen.logmsg(' %s ',num2str(kJoints));
    symname = ['coriolis_row_',num2str(kJoints)];
    fname = fullfile(CGen.sympath,[symname,'.mat']);

    if exist(fname,'file')
        tmpStruct = load(fname);
    else
        error ('genmfuncoriolis:SymbolicsNotFound','Save symbolic expressions to disk first!')
    end

    funfilename = fullfile(CGen.robjpath,[symname,'.m']);

    matlabFunction(tmpStruct.(symname),'file',funfilename,...              % generate function m-file
        'outputs', {'Crow'},...
        'vars', {'rob',q,qd});
    hStruct = createHeaderStructRow(CGen.rob,kJoints,symname);                 % replace autogenerated function header
    replaceheader(CGen,hStruct,funfilename);

end
CGen.logmsg('\t%s\n',' done!');

CGen.logmsg([datestr(now),'\tGenerating full Coriolis matrix m-function: ']);

    funfilename = fullfile(CGen.robjpath,'coriolis.m');
    hStruct = createHeaderStructFull(CGen.rob,funfilename);

    fid = fopen(funfilename,'w+');

    fprintf(fid, '%s\n', ['function C = coriolis(rob,q,qd)']);                 % Function definition
    fprintf(fid, '%s\n',constructheaderstring(CGen,hStruct));                   % Header

    fprintf(fid, '%s \n', 'C = zeros(length(q));');                        % Code
    for iJoints = 1:nJoints
        funCall = ['C(',num2str(iJoints),',:) = ','rob.coriolis_row_',num2str(iJoints),'(q,qd);'];
        fprintf(fid, '%s \n', funCall);
    end

    fclose(fid);

    CGen.logmsg('\t%s\n',' done!');
end

function hStruct = createHeaderStructRow(rob,curJointIdx,fName)
[~,hStruct.funName] = fileparts(fName);
hStruct.shortDescription = ['Computation of the robot specific Coriolis matrix row for joint ', num2str(curJointIdx), ' of ',num2str(rob.n),'.'];
hStruct.calls = {['Crow = ',hStruct.funName,'(rob,q,qd)'],...
    ['Crow = rob.',hStruct.funName,'(q,qd)']};
hStruct.detailedDescription = {'Given a full set of joint variables and their first order temporal derivatives this function computes the',...
                               ['Coriolis matrix row number ', num2str(curJointIdx),' of ',num2str(rob.n),' for ',rob.name,'.']};
hStruct.inputs = { ['rob: robot object of ', rob.name, ' specific class'],...
                   ['q:  ',int2str(rob.n),'-element vector of generalized'],...
                   '     coordinates';
                   ['qd:  ',int2str(rob.n),'-element vector of generalized'],...
                   '     velocities', ...
                   'Angles have to be given in radians!'};
hStruct.outputs = {['Crow:  [1x',int2str(rob.n),'] row of the robot Coriolis matrix']};
hStruct.references = {'1) Robot Modeling and Control - Spong, Hutchinson, Vidyasagar',...
    '2) Modelling and Control of Robot Manipulators - Sciavicco, Siciliano',...
    '3) Introduction to Robotics, Mechanics and Control - Craig',...
    '4) Modeling, Identification & Control of Robots - Khalil & Dombre'};
hStruct.authors = {'This is an autogenerated function.',...
    'Code generator written by:',...
    'Joern Malzahn',...
    '2012 RST, Technische Universitaet Dortmund, Germany',...
    'http://www.rst.e-technik.tu-dortmund.de'};
hStruct.seeAlso = {'coriolis'};
end

function hStruct = createHeaderStructFull(rob,fName)
[~,hStruct.funName] = fileparts(fName);
hStruct.shortDescription = ['Coriolis matrix for the ',rob.name,' arm'];
hStruct.calls = {['Crow = ',hStruct.funName,'(rob,q,qd)'],...
    ['Crow = rob.',hStruct.funName,'(q,qd)']};
hStruct.detailedDescription = {'Given a full set of joint variables and their first order temporal derivatives the function computes the',...
                               'Coriolis matrix of the robot.'};
hStruct.inputs = { ['rob: robot object of ', rob.name, ' specific class'],...
                   ['q:  ',int2str(rob.n),'-element vector of generalized'],...
                   '     coordinates';
                   ['qd:  ',int2str(rob.n),'-element vector of generalized'],...
                   '     velocities', ...
                   'Angles have to be given in radians!'};
hStruct.outputs = {['C:  [',int2str(rob.n),'x',int2str(rob.n),'] Coriolis matrix']};
hStruct.references = {'1) Robot Modeling and Control - Spong, Hutchinson, Vidyasagar',...
    '2) Modelling and Control of Robot Manipulators - Sciavicco, Siciliano',...
    '3) Introduction to Robotics, Mechanics and Control - Craig',...
    '4) Modeling, Identification & Control of Robots - Khalil & Dombre'};
hStruct.authors = {'This is an autogenerated function!',...
    'Code generator written by:',...
    'Joern Malzahn',...
    '2012 RST, Technische Universitaet Dortmund, Germany',...
    'http://www.rst.e-technik.tu-dortmund.de'};
hStruct.seeAlso = {'inertia'};
end

author: petercorke

repo: robotics-toolbox-matlab


text:
% J. Shi and J. Malik,
% "Normalized Cuts and Image Segmentation",
% In Proc. IEEE Conf. Computer Vision and Pattern Recognition,
% pages 731-737, 1997.

% Asad Ali
% GIK Institute of Engineering Sciences & Technology, Pakistan
% Email: asad_82@yahoo.com

% CONCEPT: Introduced the use of 2nd generalized eigenvector for clustering
clear all;
close all;

% generate the data
data = GenerateData(2);
figure,plot(data(:,1), data(:,2),'r+'),title('Original Data Points'); grid on;shg

affinity = CalculateAffinity(data);
figure,imshow(affinity,[]), title('Affinity Matrix');

% compute the degree matrix
for i=1:size(affinity,1)
    D(i,i) = sum(affinity(i,:));
end

% compute the unnormalized graph laplacian matrix
L = D - affinity;

[eigVectors,eigValues] = eig(L);

% plot the eigen vector corresponding to the 2nd largest eigen value
figure,plot(eigVectors(:,2),'r*'),title('2nd Largest Eigenvector');

% threshold the eigen vectors
[xx1,yy1,val1] = find(eigVectors(:,2) > 0.15);
[xx2,yy2,val2] = find(eigVectors(:,2) > 0 & eigVectors(:,2) < 0.15);
[xx3,yy3,val3] = find(eigVectors(:,2) < 0);

figure,
hold on;
plot(data(xx1,1),data(xx1,2),'m*');
plot(data(xx2,1),data(xx2,2),'b*');
plot(data(xx3,1),data(xx3,2),'g*');
hold off;
title('Clustering Results using 2nd Generalized Eigen Vector');
grid on;shg


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

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

% Last Modified by GUIDE v2.5 06-Feb-2008 23:04:44

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

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

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

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

% Update handles structure
guidata(hObject, handles);

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

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

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


text:
function [pred,model] = mdlTrans_lapridge(Xl,Yl,Xu,param)
%Laplacian ridge regression
%
%	This is a model-level semi-supervised (transductive) learning
% algorithm. It may be used as a baseline to other domain adaptation
% methods.
% Application scope:
%	+ partially labeled data
%	+ label type: regression

% Xl:	labeled X
% Yl:	labels of X
% Xu:	unlabeled X
%
% param: Struct of hyper-parameters, please see the first cell of this
%	program ("default parameters") for details. You can set parameter p to
%	x by setting param.p = x. For parameters that are not set, default
%	values will be used.
%
% pred : predicted labels for Xu
%
% ref: M. Belkin, P. Niyogi, and V. Sindhwani, "Manifold regularization: A
% geometric framework for learning from labeled and unlabeled examples,"
% J. Mach. Learn. Res., 2006.
% Copyright 2016 Ke YAN, Tsinghua Univ. http://yanke23.com , xjed09@gmail.com

addpath(genpath('lapsvmp_v02'))

%% default parameters
% please find the details and other parameters in functions calckernel and laplacian
t = 2;
g = .001;
gamma_I = 1;
gamma_A = 1e-5;
knn = 5;

defParam

%% kernels
if t==0, Kernel = 'linear'; KernelParam = 0;
elseif t==1, Kernel = 'poly'; KernelParam = g;
elseif t==2, Kernel = 'rbf'; KernelParam = sqrt(1/2/g);
else error('wrong ker'); end

%% make options
options = make_options('gamma_I',gamma_I,'gamma_A',gamma_A,'NN',knn,...
	'Kernel',Kernel,'KernelParam',KernelParam);
options.LaplacianNormalize = 0;
if t==1,options.polyDegree = 2;end

%% sort data
nTrain = size(Xl,1);
nTest = size(Xu,1);
Xl = [ones(nTrain,1),Xl];
Xu = [ones(nTest,1),Xu];

%% compute
X = [Xl;Xu];
K = calckernel(options,X,X);
L = laplacian(options,X);

J = diag([ones(1,nTrain),zeros(1,nTest)]);
I = eye(nTrain+nTest);
Yl = [Yl;zeros(nTest,1)];
alpha = (J*K + gamma_A*nTrain*I + gamma_I*nTrain/(nTrain+nTest)^2*L*K) \ Yl;

%% predict
out = K*alpha;
pred = out(nTrain+1:nTrain+nTest);
model.K = K;
model.alpha = alpha;

end

author: viggin

repo: domain-adaptation-toolbox


text:
A=magic(4)
sum(A)
sum(A')
sum(diag(A))
sum(diag(fliplr(A)))


text:
function deepLearningExamples(type)
%  deepLearningExamples(demo)
%------------------------------------------------------------------------------
%  A set of demos for the MEDAL toolbox
%------------------------------------------------------------------------------
% <demo> is a string indicating the demonstration type. Options for <demo> are:
%
%    'rbm'          -- Train a binary-input RBM on the MNIST dataset and use
%                      it to classify test cases
%
%    'grbm'         -- Train a gaussian-input RBM on a toy Gaussian Mixture
%                      dataset and use it to classify test cases
%
%    'mlnn'         -- Train single and ulti-layer Neural Networks on the
%                      MNIST dataset and use them to classify test cases
%
%    'dbn'          -- Train a Deep Belief Network on MNIST and use the
%                      top layer to classify test cases.
%
%    'ae'           -- Train sparse and denoising autoencoders on MNIST digits
%
%    'dae'          -- Train a deep autoencoder on MNIST and finetune
%                      the network to classify test cases.
%
%    'crbm'         -- Train a Convolutional RBM on the MNIST dataset.
%
%    'mlcnn'        -- Train a mult-layer Convolutional Neural Network
%                      on MNIST and use it to classify test cases
%
%    'drbm'         -- Train a dynamic RBM on a toy spatiotemporal dataset
%
%    'mcrbm'        -- Train a Mean-Covariance RBM on a color image patches
%
%    'all'          -- Run all demos in sequence
%------------------------------------------------------------------------------
% DES
% stan_s_bury@berkeley.edu

if notDefined('type'), type = 'all'; end

switch lower(type)

case {'binary rbm','rbm'}
% BERNOULLI-BERNOULLI RESTRICTED BOLZMANN MACHINES
clear all; close all; clc
fprintf('\nRunning Binary RBM demo for Classification (demoBinaryRBM_MNIST.m)\n')
demoBinaryRBM_MNIST
fprintf('\nRBM demo finished.\n')

case {'gaussian rbm','grbm'}
% GAUSSIAN-BERNOULLI RESTRICTED BOLZMANN MACHINE ON GMM DATASET
clear all; close all; clc
fprintf('\nRunning Gaussian-Bernoulli RBM demo (demoGaussianRBM_Classifer.m)\n')
demoGaussianRBM_GMM
fprintf('\nGaussian RBM demo finished.\n')

case {'mlnn'}
% MULTI-LAYER NEURAL NETWORK
clear all; close all; clc
fprintf('\nRunning Multi-layer Neural Net demo (demoMLNN.m)\n')
demoMLNN
fprintf('\nNeural Net demo finished\n')

case {'dbn'}
% DEEP BELIEF NETWORK FOR CLASSIFICATION
clear all; close all; clc
fprintf('\nRunning DBN for Classification demo (demoDBN.m)\n')
demoDBN
fprintf('\nDBN demo finished\n')

case {'autoencoder','ae'}
% DENOISING AUTOENCODER
clear all; close all; clc
fprintf('\nRunning Denoising and Sparse Autoencoder demo (demoAE.m)\n')
demoAE;
fprintf('\nAutoencoder demo finished.\n')

case {'deep autoencoder','dae'}
% DENOISING AUTOENCODER
clear all; close all; clc
fprintf('\nRunning Deep Autoencoder demo (demoDAE.m)\n')
demoDAE;
fprintf('\nDeep Autoencoder demo finished.\n')

case {'crbm','convolutional rbm'}
% CONVOLUTIONAL RBM
clear all; close all; clc
fprintf('\nRunning Convolutional RBM demo (demoBinaryCRBM)\n')
demoBinaryCRBM;
fprintf('\nConvolutional RBM demo finished.\n')

case {'mlcnn','cnn','convolutional neural net'}
% CONVOLUTIONAL NEURAL NETWORK
clear all; close all; clc
fprintf('\nRunning Convolutional Neural Net demo (demoMLCNN.m)\n')
demoMLCNN
fprintf('\nConvolutional Neural Net demo finished.\n')

case {'drbm','dynamic rbm'}
% DYNAMIC/CONDITIONAL RBM
clear all; close all; clc
fprintf('\nRunning Dynamic RBM demo (demoDRBM.m)\n')
demoDRBM
fprintf('\nDynamic RBM Demo finished.\n')

case {'mcrbm','mean covariance rbm'}
% MEAN-COVARIANCE RBM
clear all; close all; clc
fprintf('\nRunning Mean-Covariance RBM demo (demoMCRBM.m)\n')
demoMCRBM
fprintf('\nMean-Covariance RBM Demofinished.\n')

case 'all'
demos = {'rbm','grbm','mlnn','dbn','ae','dae','crbm','mlcnn','drbm','mcrbm'};
	for iD = 1:numel(demos)
		deepLearningExamples(demos{iD});
		fprintf('\nPress any key to resume');
		pause
	end
	clc; close all; drawnow; clear all;
end

author: dustinstansbury

repo: medal


text:
function pass = test_dot(pref)
% Test DOT
if ( nargin == 0 )
    pref = chebfunpref;
end
tol = 50*pref.cheb3Prefs.chebfun3eps;

% Check definition:
F = chebfun3v(@(x,y,z) cos(x), @(x,y,z) sin(y), @(x,y,z) exp(z));
G = chebfun3v(@(x,y,z) x, @(x,y,z) y, @(x,y,z) z);
dotF1 = dot(F, G);
dotF2 = F' * G;
pass(1) = norm(dotF1 - dotF2) < tol;

% Check definition again:
F = chebfun3v(@(x,y,z) cos(x), @(x,y,z) sin(y), @(x,y,z) x.*y);
G = chebfun3v(@(x,y,z) x, @(x,y,z) y, @(x,y,z) z);
dotF1 = dot(F, G);
dotF2 = F' * G;
pass(2) = norm(dotF1 - dotF2) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef KalmanPredictorX < PredictorX
% KalmanPredictorX class
%
% Summary of KalmanPredictorX:
% This is a class implementation of a standard Kalman Predictor.
%
% KalmanPredictorX Methods:
%   + KalmanPredictorX  - Constructor method
%   + predict - Performs full KF prediction step (both state and measurement)
%   + predictState - Performs KF state prediction step
%   + predictMeasurement - Preforms KF measurement prediction step
%
% (+) denotes puplic properties/methods
%
% See also DynamicModelX, ObservationModelX and ControlModelX template classes

    methods (Static)

        function [xPred, PPred, yPred, S, Pxy] = predict(x,P,F,Q,H,R,u,B,O)
        % predict Perform the discrete-time KF state and measurement
        % prediction steps, under the assumption of additive process noise.
        %
        % Parameters
        % ----------
        % x: column vector
        %   The (xDim x 1) state estimate at the previous time-step.
        % P: matrix
        %   The (xDim x xDim) state covariance matrix at the previous
        %   time-step.
        % F: matrix
        %   An (xDim x xDim) state transition matrix.
        % Q: matrix
        %   The (xDim x xDim) process noise covariance matrix.
        % H: matrix
        %   A (xDim x yDim) measurement matrix.
        % R: matrix
        %   The (yDim x yDim) measurement noise covariance matrix.
        % u: column vector, optional
        %   A optional (xDim x 1) control input.
        %   If omitted, no control input is used.
        % B: matrix, optional
        %   An optional (xDim x xDim) control gain matrix.
        %   If omitted, B is assumed to be 1.
        % O: matrix, optional
        %   An optional (xDim x xDim) control noise covariance
        %   matrix. If omitted, Q is assumed to be 0.
        %
        % Returns
        % -------
        % xPred: column vector
        %   The (xDim x 1) predicted state estimate.
        % PPred: matrix
        %   The (xDim x xDim) predicted state covariance matrix.
        % yPred: column vector
        %   The (yDim x 1) predicted measurement estimate.
        % Pxy: matrix
        %   The (xDim x yDim) cross-covariance matrix.
        % S: matrix
        %   The (yDim x yDim) innovation covariance matrix.
        %
        %October 2017 Lyudmil Vladimirov, University of Liverpool.

            switch(nargin)
                case(6)
                    u = 0;
                    B = 0;
                    O = 0;
                case(7)
                    B = 1;
                    O = 0;
                case(8)
                    O = 0;
            end

           [xPred, PPred] = KalmanPredictorX.predictState(x,P,F,Q,u,B,O);
           [yPred, S, Pxy] = KalmanPredictorX.predictMeasurement(xPred,PPred,H,R);
        end

        function [xPred, PPred] = predictState(x,P,F,Q,u,B,Qu)
        % predictState Perform the discrete-time KF state prediction
        % step, under the assumption of additive process noise.
        %
        % Parameters
        % ----------
        % x: column vector
        %   The (xDim x 1) state estimate at the previous time-step.
        % P: matrix
        %   The (xDim x xDim) state covariance matrix at the previous
        %   time-step.
        % F: matrix
        %   An (xDim x xDim) state transition matrix.
        % Q: matrix
        %   The (xDim x xDim) process noise covariance matrix.
        % u: column vector, optional
        %   An optional (xDim x 1) control input.
        %   If omitted, no control input is used.
        % B: matrix, optional
        %   An optional (xDim x xDim) control gain matrix.
        %   If omitted, B is assumed to be 1.
        % O: matrix, optional
        %   An optional (xDim x xDim) control noise covariance
        %   matrix. If omitted, Q is assumed to be 0.
        %
        % Returns
        % -------
        % xPred: column vector
        %   The (xDim x 1) predicted state estimate.
        % PPred: matrix
        %   The (xDim x xDim) predicted state covariance matrix.
        %
        %October 2017 Lyudmil Vladimirov, University of Liverpool.

            switch(nargin)
                case(4)
                    u  = 0;
                    B  = 0;
                    Qu = 0;
                case(5)
                    B  = 1;
                    Qu = 0;
                case(6)
                    Qu = 0;
            end

            % Compute predicted state mean and covariance
            xPred = F*x + B*u;
            PPred =F*P*F' + Q + B*Qu*B';
        end

        function [yPred, S, Pxy] = predictMeasurement(xPred,PPred,H,R)
        % predictMeasurement Perform the discrete-time KF observation prediction
        % step, under the assumption of additive process noise.
        %
        % Parameters
        % ----------
        % xPred: column vector
        %   The (xDim x 1) predicted state estimate at the current
        %   time-step.
        % PPred: matrix
        %   The (xDim x xDim) predicted state covariance matrix at
        %   the current time-step.
        % H: matrix
        %   An (xDim x yDim) measurement matrix.
        % R: matrix
        %   The (yDim x yDim) measurement noise covariance matrix.
        %
        % Returns
        % -------
        % yPred: column vector
        %   The (yDim x 1) predicted measurement estimate.
        % Pxy: matrix
        %   The (xDim x yDim) cross-covariance matrix.
        % S: matrix
        %   The (yDim x yDim) innovation covariance matrix.
        %
        %October 2017 Lyudmil Vladimirov, University of Liverpool.

            % Compute predicted measurement mean and covariance
            yPred   = H*xPred;
            Pxy     = PPred*H';
            S       = H*PPred*H' + R;
        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2019.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtOprHelmholtzBWneu.m                        |
%|    #    |   VERSION    : 0.61                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 05.09.2019                                    |
%| ( === ) |   SYNOPSIS   : Solve neumann scatering problem with          |
%|  `---'  |                Brackage-Werner formulation                   |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Parameters
N   = 1e3
tol = 1e-3
typ = 'P1'
gss = 3
X0  = [0 0 -1]

% Spherical mesh
sphere = mshSphere(N,1);
sigma  = dom(sphere,gss);
figure
plot(sphere)
axis equal

% Radiative mesh
square     = mshSquare(5*N,[5 5]);
square.vtx = [square.vtx(:,1) zeros(size(square.vtx,1),1) square.vtx(:,2)];
hold on
plot(square)

% Frequency adjusted to maximum esge size
stp = sphere.stp;
k   = 1/stp(2)
f   = (k*340)/(2*pi)

% Incident wave
PW = @(X) exp(1i*k*X*X0');
gradxPW{1} = @(X) 1i*k*X0(1) .* PW(X);
gradxPW{2} = @(X) 1i*k*X0(2) .* PW(X);
gradxPW{3} = @(X) 1i*k*X0(3) .* PW(X);

% Incident wave representation
plot(sphere,real(PW(sphere.vtx)))
plot(square,real(PW(square.vtx)))
title('Incident wave')
xlabel('X');   ylabel('Y');   zlabel('Z');
hold off
view(0,10)

%%% PREPARE OPERATORS
disp('~~~~~~~~~~~~~ PREPARE OPERATORS ~~~~~~~~~~~~~')

% Green kernel function --> G(x,y) = exp(ik|x-y|)/|x-y|
Gxy      = '[exp(ikr)/r]';
gradxGxy = {'gradx[exp(ikr)/r]1','gradx[exp(ikr)/r]2','gradx[exp(ikr)/r]3'};
gradyGxy = {'grady[exp(ikr)/r]1','grady[exp(ikr)/r]2','grady[exp(ikr)/r]3'};

% Finite elements
u = fem(sphere,typ);
v = fem(sphere,typ);

% Coupling coeff
beta = 1i*k;

% Number of pool
Nlab = length(Composite);

% Domain decomposition for u
[Ilab,sigmaLab,uLab] = femSubdivide(sigma,u,Nlab,10);
drawnow

% Parallel loop for Full matrix
tic
spmd
    % Initialize composite
    M = cell(1,numlabs);
    P = cell(1,numlabs);

    % Normal loop
    for j = 1:Nlab
        % Mass matrix
        Id = integral(sigmaLab{labindex},uLab{labindex},uLab{j});

        % Hypersingular
        H  = oprIntegral('H',k,sigmaLab{labindex},uLab{labindex},sigmaLab{j},uLab{j},tol);
        Hr = 1/(4*pi).*(k^2 ...
            * regularize(sigmaLab{labindex},sigmaLab{j},ntimes(uLab{labindex}),'[1/r]',ntimes(uLab{j})) ...
            - regularize(sigmaLab{labindex},sigmaLab{j},nxgrad(uLab{labindex}),'[1/r]',nxgrad(uLab{j})));

        % Double layer transpose : switch indices for block transposition
        Dt  = oprIntegral('Dt',k,sigmaLab{labindex},uLab{labindex},sigmaLab{j},uLab{j},tol);
        Dtr = 1/(4*pi).*regularize(sigmaLab{j},sigmaLab{labindex},uLab{j},'grady[1/r]',ntimes(uLab{labindex})).';

        % Neumann Brackage-Werner : [1i*k*beta*(-Id/2 + Dt) - H]mu
        M{j} = beta.*(-0.5*Id + (Dt+Dtr)) - (H+Hr);
        P{j} = sparse(M{j},Id);
    end
end
toc

% Define LHS
LHS = @(V) spmdProduct(Ilab,M,V);

% Finite element incident wave trace --> \int_Sx psi(x)' pw(x) dx
RHS = - integral(sigma,ntimes(u),gradxPW);

tic
LHS(RHS);
toc

%%% SOLVE LINEAR PROBLEM
disp('~~~~~~~~~~~~~ SOLVE LINEAR PROBLEM ~~~~~~~~~~~~~')

% Block matrix
tic
P = sparse(bmm(Ilab,Ilab,P));
toc

% Factorization for preconditionning
tic
[L,U] = ilu(P);
toc

% Solve linear system : [1i*k*beta*S - (Id/2 + D)] = P0
tic
mu = mgcr(LHS,RHS,[],tol,100,L,U);
toc

% Jump for derivative
lambda = beta * mu;

%%% INFINITE SOLUTION
disp('~~~~~~~~~~~~~ INFINITE RADIATION ~~~~~~~~~~~~~')

% Plane waves direction
theta = 2*pi/1e3 .* (1:1e3)';
nu    = [sin(theta),zeros(size(theta)),cos(theta)];

% Green kernel function
Ginf      = '[exp(-ikxy)]';
gradxGinf = {'gradx[exp(-ikxy)]1','gradx[exp(-ikxy)]2','gradx[exp(-ikxy)]3'};

% Finite element infinite operators
Sinf = 1/(4*pi) .* integral(nu,sigma,Ginf,k,v,tol);
Dinf = 1/(4*pi) .* integral(nu,sigma,gradxGinf,k,ntimes(v),tol);

% Finite element radiation
sol = Sinf*lambda - Dinf*mu;

% Analytical solution
ref = sphereHelmholtz('inf','neu',1,k,nu);
norm(ref-sol,2)/norm(ref,2)
norm(ref-sol,'inf')/norm(ref,'inf')

% Graphical representation
figure
plot(theta,log(abs(sol)),'b',theta,log(abs(ref)),'--r')

%%% DOMAIN SOLUTION
disp('~~~~~~~~~~~~~ RADIATION ~~~~~~~~~~~~~')

% Finite element mass matrix --> \int_Sx psi(x)' psi(x) dx
Id = integral(sigma,u,v);

% Finite element boundary operator --> \int_Sx \int_Sy psi(x)' G(x,y) psi(y) dx dy
tic
Sbnd = 1/(4*pi) .* (integral(sigma,sigma,u,Gxy,k,v,tol) + ...
    regularize(sigma,sigma,u,'[1/r]',v));
toc

% Finite element boundary operator --> \int_Sx \int_Sy psi(x)' dny G(x,y) psi(y) dx dy
tic
Dbnd = 1/(4*pi) .* (integral(sigma,sigma,u,gradyGxy,k,ntimes(v),tol) + ...
    regularize(sigma,sigma,u,'grady[1/r]',ntimes(v)));
toc

% Boundary solution
Psca = Id\(Sbnd*lambda - (0.5*Id*mu + Dbnd*mu));
Pinc = PW(u.dof);
Pbnd = Pinc + Psca;

% Finite element radiative operator --> \int_Sy G(x,y) psi(y) dy
tic
Sdom = 1/(4*pi) .* (integral(square.vtx,sigma,Gxy,k,v,tol) + ...
    regularize(square.vtx,sigma,'[1/r]',v));
toc

% Finite element radiative operator --> \int_Sx \int_Sy psi(x)' grady(G(x,y)) ny.psi(y) dx dy
tic
Ddom = 1/(4*pi) .* ( integral(square.vtx,sigma,gradyGxy,k,ntimes(v),tol) + ...
    regularize(square.vtx,sigma,'grady[1/r]',ntimes(v)) );
toc

% Domain solution
Psca = Sdom*lambda - Ddom*mu;
Pinc = PW(square.vtx);
Pdom = Pinc + Psca;

% Annulation sphere interieure
r             = sqrt(sum(square.vtx.^2,2));
Pdom(r<=1.01) = Pinc(r<=1.01);

% Graphical representation
figure
plot(sphere,abs(Pbnd))
axis equal;
hold on
plot(square,abs(Pdom))
title('Total field solution')
colorbar
hold off
view(0,10)

%%% ANAYTICAL SOLUTIONS FOR COMPARISONS
% Analytical solution
Pbnd = sphereHelmholtz('dom','neu',1,k,1.001*sphere.vtx) + PW(sphere.vtx);
Pdom = sphereHelmholtz('dom','neu',1,k,square.vtx) + PW(square.vtx);

% Solution representation
figure
plot(sphere,abs(Pbnd))
axis equal;
hold on
plot(square,abs(Pdom))
title('Analytical solution')
colorbar
hold off
view(0,10)

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function varargout = ode113(varargin)
%ODE113   Solve stiff differential equations and DAEs. Output a CHEBFUN.
%
% This syntax is deprecated. Please use chebfun.ode113(...) instead.

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

warning('CHEBFUN:DOMAIN:ode113:deprecated', ...
    ['Usage of ODE113 via the @DOMAIN class is deprecated and may be ', ...
    'removed from future releases. Please use chebfun.ode113(...) instead.']);

varargin = domain.toDouble(varargin{:});

varargout{1:nargout} = chebfun.ode113(varargin{:});

end


text:
function [w,h]=addbg(win, hin, bg)
% [w,h]=addbg(win, hin, bg)
% add the background component

sizew=size(win);
sizeh=size(hin);
w_bg = ones(sizew(1), 1)/sizew(1); %background component
w = [win, w_bg];

h_bg = bg*size(w,1)*ones(1, sizeh(2)); %background component
h=[hin; h_bg];

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function isBipartite=graphIsBipartite(A)
%%GRAPHISBIPARTITE Given a nondirectional graph represented by a symmetric
%                  adjacency matrix A, determine whether the graph is
%                  bipartite.
%
%INPUTS: A An nXn symmetric adjacency matrix. If element i is adjacent to
%          element j, then A(i,j) is nonzero. Otherwise, it is zero.
%          Diagonal elements are typically zero (if not, then the graph is
%          not bipartite).
%
%OUTPUTS: isBipartite This is true if the graph is bipartite and false
%                     otherwise.
%
%This function implements Algorithm B in Chapter 7 of [1], modified for an
%adjacency matrix. When very large, space graphs are used, then the
%adjacency matrix representation is not ideal.
%
%EXAMPLE:
%We consider four graphs. Two are bipartite and two are not.
%A1 and A2 are:
%A1:     A2:
%1--2    1--2
% \/      \/
% /\      /\
%3--4    3--6
%         \/
%         5
%and graphs A3 and A4 are:
%A3:        A4:
%1-----2    1-----2
%|\   /|    |\   /|
%| 5-6 |    | 5-6 |
%| | | |    | |\| |
%| 7-8 |    | 7-8 |
%|/   \|    |/   \|
%3-----4    3-----4
% A1=[0,1,0,1;
%     1,0,1,0;
%     0,1,0,1;
%     1,0,1,0];
% A2=[0,1,0,1,0;
%     1,0,1,0,0;
%     0,1,0,1,1;
%     1,0,1,0,1;
%     0,0,1,1,0];
% A3=[0,1,1,0,1,0,0,0;
%     1,0,0,1,0,0,0,0;
%     1,0,0,1,0,0,0,0;
%     0,1,1,0,0,0,0,0;
%     1,0,0,0,0,1,1,0;
%     0,0,0,0,1,0,0,1;
%     0,0,0,0,1,0,0,1;
%     0,0,0,0,0,1,1,0];
% A4=[0,1,1,0,1,0,0,0;
%     1,0,0,1,0,0,0,0;
%     1,0,0,1,0,0,0,0;
%     0,1,1,0,0,0,0,0;
%     1,0,0,0,0,1,1,1;
%     0,0,0,0,1,0,0,1;
%     0,0,0,0,1,0,0,1;
%     0,0,0,0,1,1,1,0];
% isBi1=graphIsBipartite(A1)
% isBi2=graphIsBipartite(A2)
% isBi3=graphIsBipartite(A3)
% isBi4=graphIsBipartite(A4)
%One will find that A1 and A3 are bipartite and A2 and A4 are not.
%
%REFERENCES:
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4A: Combinatorial
%    Algorithms, Part I, Boston: Addison-Wesley, 2011.
%
%October 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n=size(A,1);

%Allocate space.
link=zeros(n,1);

%Step B1
color=-1*ones(n,1);%In [1], it is indexed from 0, not 1; we index from 1.
w=n+1;%Only need to consider the ones originating an edge.

isBipartite=true;
while(w~=1)
    w=w-1;

    %Step B3
    if(color(w)>=0)
        continue;
    end
    color(w)=0;
    link(w)=0;
    s=w;

    while(1)
        %Step B4
        u=s;
        s=link(s);

        aList=find(A(u,:));

        for a=1:length(aList)
            v=aList(a);

            if(color(v)<0)
                color(v)=1-color(u);
                link(v)=s;
                s=v;
            elseif(color(v)==color(u))
                isBipartite=false;
                return;
            end
        end
        if(s==0)
            break;
        end
    end
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function n = length( x )
s = size( x );
if any( s == 0 ),
   n = 0;
else
   n = max( s );
end

% 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:
function test_bug1618

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/bug1618/bug1618.dat');

h = ft_read_header(filename);

disp(h)
disp(h.chanunit)
disp(h.chantype)

X = ft_read_data(filename);
assert(~any(isnan(X(:))), 'NaN values are not expected for this dataset!');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function path = voc2007_devkit()
    path = './datasets/VOCdevkit2007';
end

author: ShaoqingRen


text:
function a = givens_llt ( n )

%*****************************************************************************80
%
%% GIVENS_LLT returns the Cholesky factor of the GIVENS matrix.
%
%  Example:
%
%    N = 5
%
%    1    0        0        0       0
%    1  sqrt(2)    0        0       0
%    1  sqrt(2)  sqrt(2)    0       0
%    1  sqrt(2)  sqrt(2)  sqrt(2)   0
%    1  sqrt(2)  sqrt(2)  sqrt(2) sqrt(2)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  a(1:n,1) = 1.0;

  for i = 1 : n
    a(i,2:i) = sqrt ( 2.0 );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function turn = angle_turn_2d ( p1, p2, p3 )

%*****************************************************************************80
%
%% ANGLE_TURN_2D computes a turning angle in 2D.
%
%  Discussion:
%
%    This routine is most useful when considering the vertices of a
%    polygonal shape.  We wish to distinguish between angles that "turn
%    in" to the shape, (between 0 and 180 degrees) and angles that
%    "turn out" (between 180 and 360 degrees), as we traverse the boundary.
%
%    If we compute the interior angle and subtract 180 degrees, we get the
%    supplementary angle, which has the nice property that it is
%    negative for "in" angles and positive for "out" angles, and is zero if
%    the three points actually lie along a line.
%
%    Assuming P1, P2 and P3 define an angle, the TURN can be
%    defined to be either:
%
%    * the supplementary angle to the angle formed by P1=P2=P3, or
%
%    * the angle between the vector ( P3-P2) and the vector -(P1-P2),
%      where -(P1-P2) can be understood as the vector that continues
%      through P2 from the direction P1.
%
%    The turning will be zero if P1, P2 and P3 lie along a straight line.
%
%    It will be a positive angle if the turn from the previous direction
%    is counterclockwise, and negative if it is clockwise.
%
%    The turn is given in radians, and will lie between -PI and PI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2), P2(2), P3(2), the points that form
%    the angle.
%
%    Output, real TURN, the turn angle, between -PI and PI.
%
  p(1) = ( p3(1) - p2(1) ) * ( p1(1) - p2(1) ) ...
       + ( p3(2) - p2(2) ) * ( p1(2) - p2(2) );

  p(2) = ( p3(1) - p2(1) ) * ( p1(2) - p2(2) ) ...
       - ( p3(2) - p2(2) ) * ( p1(1) - p2(1) );

  if ( p(1) == 0.0 & p(2) == 0.0 )
    turn = 0.0;
  else
    turn = pi - r8_atan ( p(2), p(1) );
  end

  return
end


text:
function pts = DIR_get_all_points_for_the_case(caseN, basepth)
    cn = num2str(caseN);
    path = [basepth, ...
            '4dCT/Case', cn, 'Pack/Case', cn, 'Pack/ExtremePhases/'];
    pts = struct();
    pts.extreme.b = read_DIR_points_file([path, 'Case', cn, ...
        '_300_T00_xyz.txt']);
    pts.extreme.e = read_DIR_points_file([path, 'Case', cn, ...
        '_300_T50_xyz.txt']);

    path2 = [basepth, ...
            '4dCT/Case', cn, 'Pack/Case', cn, 'Pack/Sampled4D/'];
    for i = 1 : 6
        ni = num2str(i - 1);
        pts.smp{i}.pts = read_DIR_points_file([path2, 'case', cn, ...
            '_4D-75_T', ni, '0.txt']);
    end
end

repo: pTVreg