text:
clc;
close all;
clear;

mDir = fileparts(mfilename('fullpath'));
cd(mDir);
addpath(genpath('../'));

%% Params:

exp = 'ALL_4_SEGS'; % Select which configuration to use: 'JUST_LAB_SEG', 'ALL_4_SEGS'
calcInLog=false; %false (for linear vus), true (for logarithmic vus)

%Change these paths:
originalVOCDirectory = TOFILL; %Top directory of VOC 2007 dataset
parsedVOCDir = TOFILL; %Directory where parsed VOC07 will be stored

%VOC params:

imgListDir= 'complete_test_set.txt'; %File with ids of the 4952 images of the test set

params.includeBoxesWithClasses={'aeroplane', 'bicycle','boat','bottle','bus','chair','diningtable','horse','motorbike','person','pottedplant','sofa','train','tvmonitor'};
params.omitImagesWithClasses={'bird','car','cat','cow','dog','sheep'};
params.considerDifficult=true;

nMaximumWindows = 10000;
evalParams.nWindows = unique(round(logspace(0, log10(nMaximumWindows), 100)));
evalParams.ious = 0 : 0.05 : 1;

if(strcmp(exp, 'ALL_4_SEGS'))
  configFile = '../config/rp_4segs.mat';
  proposalsDir = [mDir '/tmp_prop_dir'];
elseif(strcmp(exp, 'JUST_LAB_SEG'))
  configFile = '../config/rp.mat';
  proposalsDir = [mDir '/tmp_prop_dir_lab'];
else
  assert(false);
end
resultDirs={proposalsDir};

%% Parse VOC2007 data:

AdaptVOC2007Data(originalVOCDirectory,imgListDir,parsedVOCDir,params);

%% Compute proposals and individual detection rates:

ComputeProposals(parsedVOCDir, proposalsDir, configFile, evalParams);

%% Aggregate detection rates:

AggregateDRs(proposalsDir, parsedVOCDir, evalParams);

%% Compute Volume Under Surface:

nWindowsRange=[0 10000];
iouRange=[0.5 1];

nResults=numel(resultDirs);
drs=[];
for k=1:nResults
  data=load([resultDirs{k} '/aggregated_drs/drs.mat']);

  [cropped,legendStr]=CropEvaluations({data}, nWindowsRange, iouRange,{''});

  data=cropped{1};

  if(exist('ious','var'))
    assert( all(ious==data.ious));
  else
    ious=data.ious;
    nIous=length(ious);
  end

  if(exist('nWindows','var'))
    assert(all(nWindows==data.nWindows));
  else
    nWindows=data.nWindows;
    nnWindows=length(nWindows);
  end

  assert( ~exist('iNWindows','var') && ~exist('iIous','var'));

  drs=cat(3,drs,data.drs);
end

if(calcInLog)
  nWindows=log10(nWindows);
end

maxDrs=max(drs,[],3);
volumes=[];
maxVolume=trapz(ious, trapz(nWindows, maxDrs));
for k=1:nResults
  assert(all(all(maxDrs>=drs(:,:,k))));
  volumes=[ volumes, trapz(ious, trapz(nWindows, drs(:,:,k)))];
end

%Worse point computation
worsePoint=[];
for k=1:nResults
  worsePoint=[worsePoint, -max(max(maxDrs-drs(:,:,k)))];
end

perVolumes = volumes./((nWindows(end)-nWindows(1))*(ious(end)-ious(1)));

assert(all(perVolumes)>=0 && all(perVolumes)<=1);

volumeRatios=volumes/maxVolume;
[~, ids]=sort(volumeRatios,'descend');

volumeRatios=volumeRatios(ids);
worsePoint=worsePoint(ids);
perVolumes=perVolumes(ids);
assert(issorted(flipdim(perVolumes,2)));
for k=1:nResults
  fName=resultDirs{ids(k)};
  if(fName(end)=='/')
    fName=fName(1:(end-1));
  end
  [~,name,~]=fileparts(fName);
  disp('############################################################################');
  disp([ 'Volume Under Surface (VUS):' num2str(perVolumes(k))]);
end

%% Display DR curve:

selIou = 0.5;
semilogx(data.nWindows, data.drs(:, find(data.ious==selIou)), 'LineWidth', 3, 'Color', 'r');
ylim([0, 1]);
ylabel('Detection Rate');
xlabel('Number of proposals')
axis square;


text:
function [ value, ifault ] = prncst ( st, idf, d )

%*****************************************************************************80
%
%% PRNCST computes the lower tail of noncentral T distribution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 January 2008
%
%  Author:
%
%    Original FORTRAN77 version by BE Cooper.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    BE Cooper,
%    Algorithm AS 5:
%    The Integral of the Non-Central T-Distribution,
%    Applied Statistics,
%    Volume 17, Number 2, 1968, page 193.
%
%  Parameters:
%
%    Input, real ST, the argument.
%
%    Input, integer IDF, the number of degrees of freedom.
%
%    Input, real D, the noncentrality parameter.
%
%    Output, real PRNCST, the value of the lower tail of
%    the noncentral T distribution.
%
%    Output, integer IFAULT, error flag.
%    0, no error occurred.
%    nonzero, an error occurred.
%
%  Local Parameters:
%
%    Local, real G1, 1.0 / sqrt(2.0 * pi)
%
%    Local, real G2, 1.0 / (2.0 * pi)
%
%    Local, real G3, sqrt(2.0 * pi)
%
  emin = 12.5;
  g1 = 0.3989422804;
  g2 = 0.1591549431;
  g3 = 2.5066282746;

  f = idf;
%
%  For very large IDF, use the normal approximation.
%
  if ( 100 < idf )

    ifault = 1;

    a = sqrt ( 0.5 * f ) * exp ( alngam ( 0.5 * ( f - 1.0 ) ) ...
    - alngam ( 0.5 * f ) ) * d;

    value = alnorm ( ( st - a ) / sqrt ( f * ( 1.0 + d * d ) ...
    / ( f - 2.0 ) - a * a ), 0 );

    return
  end

  ifault = 0;
  ioe = mod ( idf, 2 );
  a = st / sqrt ( f );
  b = f / ( f + st * st );
  rb = sqrt ( b );
  da = d * a;
  drb = d * rb;

  if ( idf == 1 )
    value = alnorm ( drb, 1 ) + 2.0 * tfn ( drb, a );
    return
  end

  sum = 0.0;

  if ( abs ( drb ) < emin )
    fmkm2 = a * rb * exp ( - 0.5 * drb * drb ) * alnorm ( a * drb, 0 ) * g1;
  else
    fmkm2 = 0.0;
  end

  fmkm1 = b * da * fmkm2;

  if ( abs ( d ) < emin )
    fmkm1 = fmkm1 + b * a * g2 * exp ( - 0.5 * d * d );
  end

  if ( ioe == 0 )
    sum = fmkm2;
  else
    sum = fmkm1;
  end

  ak = 1.0;
  fk = 2.0;

  for k = 2 : 2 : idf - 2

    fkm1 = fk - 1.0;
    fmkm2 = b * ( da * ak * fmkm1 + fmkm2 ) * fkm1 / fk;
    ak = 1.0 / ( ak * fkm1 );
    fmkm1 = b * ( da * ak * fmkm2 + fmkm1 ) * fk / ( fk + 1.0 );

    if ( ioe == 0 )
      sum = sum + fmkm2;
    else
      sum = sum + fmkm1;
    end

    ak = 1.0 / ( ak * fk );
    fk = fk + 2.0;

  end

  if ( ioe == 0 )
    value = alnorm ( d, 1 ) + sum * g3;
  else
    value = alnorm ( drb, 1 ) + 2.0 * ( sum + tfn ( drb, a ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ y, symm ] = cvx_s_symmetric_ut( m, n, symm )
%CVX_S_SYMMETRIC_UT Symmetric matrices (upper triangle storage).
if m ~= n,
    cvx_throw( 'Symmetric structure requires square matrices.' );
end
symm = false;
nsq = n * n;
ntr = 0.5 * ( nsq + n );
c   = 0 : n - 1;
c   = c( ones( 1, n ), : );
r   = c';
mn  = min( r, c );
mx  = max( r, c );
y   = mn + 0.5 * mx .* ( mx + 1 ) + 1;
y   = sparse( y( : ), 1 : nsq, 1, ntr, nsq );

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

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function C = cluster_rotation_classes(vecs, k, vis)

styles = get_vis_styles();

[idx, C]= kmeans(vecs,k);

if vis
    figure(1);
    hold on;

    % do PCA to display stuff
    if size(vecs, 2) > 2
        [coeff,score,latent,tsquared,explained,mu] = pca(vecs);
        explained
        vecs = score;
    end

    for i = 1:k
        iset = (idx == i);
        if size(vecs, 2) > 2
            scatter3(vecs(iset, 1), vecs(iset, 2), vecs(iset, 3));
        else
            scatter(vecs(iset, 1), vecs(iset, 2), styles{i});
        end
    end
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function d = xcorrrow(d,c,index)

%XCORRROW Cross correlation for one or more traces.
% D = XCORR1XR(D,C,STYLE) This private function performs cross correlations
% for a subset of traces. It differs from the XCORR1XR algorithm in that it
% fills both halves of the correlation matrix. While this is unnecessary
% when correlating all waveforms, it is necessary when doing just a subset
% of traces. This should be of little concern to most users. All steps are
% included in this function. That is, no calls to the Matlab built-in xcorr
% are used. This code is a little kludgy since it mixes both the original d
% structure (predates the internal use of waveform objects) with the later
% wave-based correlation objects. This is the reason it reads in both D and
% C arguments. STYLE denotes whether or not polynomial interpolation should
% be used to refine cross correlations to subsample precision.

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% PREP NECESSARY TERMS
[M,N] = size(d.w);
pretrig = 86400*(d.trig-d.start);   % time between trace start and trigger
l = (1./d.Fs)*[-M+1:M-1]';        % lag vector
% next two lines are equivalent ways to get normalization coefficients
wcoeff = 1./sqrt(sum(d.w.*d.w));
%for i = 1:size(d.w,2), wcoeff(i) = 1./norm(d.w(:,i)); end;

% CREATE MATRICES IF NEEDED
nn = get(c,'traces');
if (size(c.C,1) == 0)
    c.C = nan(nn);
end
if (size(c.L,1) == 0)
    c.L = nan(nn);
end

% do not overwrite matrices if they are only being added to
if size(c.C,1)==0
    d.C = eye(length(d.trig),'single');
    eraseC = 1;
else
    d.C = c.C;
    eraseC = 0;
end

if size(c.L,1)==0
    d.L = zeros(length(d.trig),'single');
    eraseL = 1;
else
    d.L = c.L;
    eraseL = 0;
end

% GET FFT OF TRACES
X = fft(d.w,2^nextpow2(2*M-1));
Xc = conj(X);
[MX,NX] = size(X);

% LOOP THROUGH ROWS OF SIMILARITY MATRIX
for n = index
    cols = 1:N;
    % multiply fourier series and transform back to time domain
    %CC = (X(:,n) * ones(1,length(cols))) .* Xc(:,cols);
    CC = repmat(X(:,n),1,length(cols)) .* Xc(:,cols);
    corr = ifft(CC);
    corr = corr([end-M+2:end,1:M],:);

    % USE POLYNOMIAL INTERPOLATION
    [maxtest,indx1] = max(corr(2:end-1,:));
    [mm,nn] = size(corr);
    indx2 = (indx1+1) + mm*[0:nn-1];     % convert to matrix index
    lag = repmat( l , 1 , size(corr,2) );  % replace size statement with nn
    lagM   = lag([ indx2-1 ; indx2 ; indx2+1 ]) + repmat(pretrig(cols)'-pretrig(n)',3,1);
    corrM  = corr([ indx2-1 ; indx2 ; indx2+1 ]);
    for z = 1:numel(cols)
        p = polyfit( lagM(:,z) , corrM(:,z) , 2 );
        Ltmp = -0.5*p(2)/p(1);
        if abs(Ltmp)<eps('single')
            Ltmp = 0;
        end
        d.L(n,cols(z))  = Ltmp;
        d.C(n,cols(z)) = polyval( p , d.L(n,cols(z)) ) .* wcoeff(n) .* wcoeff(cols(z));
    end

end

% FILL IN OTHER HALF OF MATRIX
d.C(:,index) = d.C(index,:)';

if eraseL
    d.L = [];
else
    d.L(:,index) = -1* d.L(index,:)';
end

if eraseC
    d.C = [];
else
    d.C(:,index) = d.C(index,:)';
end

author: geoscience-community-codes

repo: GISMO


text:
function failed_tutorial_natmeg_beamforming

% WALLTIME 00:20:00
% MEM 8gb
% DEPENDENCY

% this script executes the MATLAB content from
% http://www.fieldtriptoolbox.org/tutorial/natmeg/beamforming
%
% it corresponds to the wiki version of 7 October 2014

clear all
close all

cd(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/natmeg'));

load timefrequency/data_clean_MEG_responselocked.mat;
load dipolefitting/headmodel_meg.mat;

% Select time window of interest
cfg = [];
cfg.toilim = [0.35 0.85];
data_timewindow = ft_redefinetrial(cfg,data_clean_MEG_responselocked);

% Freqanalysis for beamformer
cfg = [];
cfg.channel      = {'MEG*2', 'MEG*3'};
cfg.method       = 'mtmfft';
cfg.taper        = 'dpss';
cfg.output       = 'powandcsd';
cfg.keeptrials   = 'no';
cfg.foi          = 18;
cfg.tapsmofrq    = 4;

% for common filter over conditions
powcsd_all      = ft_freqanalysis(cfg, data_timewindow);

% for conditions
cfg.trials       = find(data_timewindow.trialinfo(:,1) == 256);
powcsd_left      = ft_freqanalysis(cfg, data_timewindow);
cfg.trials       = find(data_timewindow.trialinfo(:,1) == 4096);
powcsd_right     = ft_freqanalysis(cfg, data_timewindow);

% Create leadfield grid
cfg                 = [];
cfg.channel         = {'MEG*2', 'MEG*3'};
cfg.grad            = powcsd_all.grad;
cfg.headmodel       = headmodel_meg;
cfg.dics.reducerank = 2; % default for MEG is 2, for EEG is 3
cfg.sourcemodel.resolution = 0.5;   % use a 3-D grid with a 0.5 cm resolution
cfg.sourcemodel.unit       = 'cm';
cfg.sourcemodel.tight      = 'yes';
[grid] = ft_prepare_leadfield(cfg);

cfg              = [];
cfg.channel      = {'MEG*2', 'MEG*3'};
cfg.method       = 'dics';
cfg.frequency    = 18;
cfg.sourcemodel         = grid;
cfg.headmodel    = headmodel_meg;
cfg.senstype     = 'MEG'; % Must me 'MEG', although we only kept MEG channels, information on EEG channels is still present in data
cfg.dics.keepfilter   = 'yes'; % We wish to use the calculated filter later on
cfg.dics.projectnoise = 'yes';
cfg.dics.lambda  = '5%';
source_all = ft_sourceanalysis(cfg, powcsd_all);

cfg              = [];
cfg.channel      = {'MEG*2', 'MEG*3'};
cfg.method       = 'dics';
cfg.frequency    = 18;
cfg.sourcemodel         = grid;
cfg.sourcemodel.filter  = source_all.avg.filter;
cfg.headmodel    = headmodel_meg;
cfg.senstype     ='MEG';

source_left = ft_sourceanalysis(cfg, powcsd_left);
source_right = ft_sourceanalysis(cfg, powcsd_right);

load dipolefitting/mri_realigned2.mat;

mri_resliced = ft_volumereslice([], mri_realigned2);

cfg            = [];
cfg.parameter = 'avg.pow';
source_left_int  = ft_sourceinterpolate(cfg, source_left, mri_resliced);
source_right_int  = ft_sourceinterpolate(cfg, source_right, mri_resliced);

source_diff_int  = source_left_int;
source_diff_int.avg.pow  = (source_left_int.avg.pow - source_right_int.avg.pow) ./ (source_left_int.avg.pow + source_right_int.avg.pow);

% plot
cfg = [];
cfg.method        = 'ortho';
cfg.funparameter  = 'avg.pow';
cfg.funcolorlim   = 'maxabs';
cfg.opacitylim    = [0 1e-4];
cfg.opacitymap    = 'rampup';

ft_sourceplot(cfg, source_left_int);

cfg.location = [35 -13 76];
ft_sourceplot(cfg, source_diff_int);

load dipolefitting/headmodel_eeg.mat
load timefrequency/data_clean_EEG_responselocked.mat

% select time window
cfg = [];
cfg.toilim = [0.35 0.85];
data_timewindow = ft_redefinetrial(cfg,data_clean_EEG_responselocked);

% Freqanalysis for beamformer
cfg = [];
cfg.method       = 'mtmfft';
cfg.taper        = 'dpss';
cfg.output       = 'powandcsd';
cfg.keeptrials   = 'no';
cfg.foi          = 18;
cfg.tapsmofrq    = 4;

% for common filter over conditions and full duration
powcsd_all      = ft_freqanalysis(cfg, data_timewindow);

% for conditions
cfg.trials       = find(data_timewindow.trialinfo(:,1) == 256);
powcsd_left      = ft_freqanalysis(cfg, data_timewindow);
cfg.trials       = find(data_timewindow.trialinfo(:,1) == 4096);
powcsd_right     = ft_freqanalysis(cfg, data_timewindow);

% common grid/filter
cfg                 = [];
cfg.elec            = powcsd_all.elec;
cfg.headmodel       = headmodel_eeg;
cfg.reducerank      = 3; % default is 3 for EEG, 2 for MEG
cfg.sourcemodel.resolution = 0.5;   % use a 3-D grid with a 0.5 cm resolution
cfg.sourcemodel.unit       = 'cm';
cfg.sourcemodel.tight      = 'yes';
[grid] = ft_prepare_leadfield(cfg);

% beamform common filter
cfg              = [];
cfg.method       = 'dics';
cfg.frequency    = 18;
cfg.sourcemodel         = grid;
cfg.headmodel    = headmodel_eeg;
cfg.senstype     = 'EEG'; % Remember this must be specified as either EEG, or MEG
cfg.dics.keepfilter   = 'yes';
cfg.dics.lambda       = '15%';
source_all = ft_sourceanalysis(cfg, powcsd_all);

% beamform conditions
cfg              = [];
cfg.method       = 'dics';
cfg.frequency    = 18;
cfg.sourcemodel         = grid;
cfg.sourcemodel.filter  = source_all.avg.filter; % Use the common filter
cfg.headmodel    = headmodel_eeg;
cfg.senstype     = 'EEG';

source_left = ft_sourceanalysis(cfg, powcsd_left);
source_right = ft_sourceanalysis(cfg, powcsd_right);

cfg            = [];
cfg.parameter = 'avg.pow';
source_left_int  = ft_sourceinterpolate(cfg, source_left, mri_resliced);
source_right_int  = ft_sourceinterpolate(cfg, source_right, mri_resliced);

source_diff_int  = source_left_int;
source_diff_int.avg.pow  = (source_left_int.avg.pow - source_right_int.avg.pow) ./ (source_left_int.avg.pow + source_right_int.avg.pow);

cfg = [];
cfg.method        = 'ortho';
cfg.funparameter  = 'avg.pow';
cfg.funcolorlim   = 'maxabs';

ft_sourceplot(cfg, source_left_int);

cfg.location = [-19.5 -18.5 70.5];
ft_sourceplot(cfg, source_diff_int);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
clear;close all;
setenv('LC_ALL','C');

%% settings
folder = './Train_291/';
size_input = 31;
size_label = 31;
stride = 21;
savepath = ['./train_291_' num2str(size_input) '_x234.h5'];

% scale augmentation
scale_2 = 2;
scale_3 = 3;
scale_4 = 4;

%% initialization
data = zeros(size_input, size_input, 1, 1);
label = zeros(size_label, size_label,1, 1);
padding = abs(size_input - size_label)/2;
count = 0;

%% generate data
filepaths = dir(fullfile(folder,'*.bmp'));

for i = 1 : length(filepaths)
    disp(['i = ' num2str(i)]);
    image = imread(fullfile(folder,filepaths(i).name));
    image = rgb2ycbcr(image);
    image_luminance = image(:,:,1); % lumination
    image = im2double(image_luminance);

    % scale 2
    im_label_2 = modcrop(image, scale_2);
    [hei_2,wid_2] = size(im_label_2);
    im_input_2 = imresize(imresize(im_label_2,1/scale_2,'bicubic'),[hei_2,wid_2],'bicubic');
%     imshow([im_label_2 im_input_2],[]);

    for x = 1 : stride : hei_2-size_input+1
        for y = 1 :stride : wid_2-size_input+1

            subim_input = im_input_2(x : x+size_input-1, y : y+size_input-1);
            subim_label = im_label_2(x+padding : x+padding+size_label-1, y+padding : y+padding+size_label-1);

            count=count+1;
            data_temp = subim_input;
            label_temp = subim_label;

            data(:, :, :, count) = single(data_temp);
            label(:, :, :, count) = single(label_temp);
        end
    end

    % scale 3
    im_label_3 = modcrop(image, scale_3);
    [hei_3,wid_3] = size(im_label_3);
    im_input_3 = imresize(imresize(im_label_3,1/scale_3,'bicubic'),[hei_3,wid_3],'bicubic');
%     imshow([im_label_3 im_input_3],[]);

    for x = 1 : stride : hei_3-size_input+1
        for y = 1 :stride : wid_3-size_input+1

            subim_input = im_input_3(x : x+size_input-1, y : y+size_input-1);
            subim_label = im_label_3(x+padding : x+padding+size_label-1, y+padding : y+padding+size_label-1);

            count=count+1;
            data_temp = subim_input;
            label_temp = subim_label;

            data(:, :, :, count) = single(data_temp);
            label(:, :, :, count) = single(label_temp);
        end
    end

    % scale 4
    im_label_4 = modcrop(image, scale_4);
    [hei_4,wid_4] = size(im_label_4);
    im_input_4 = imresize(imresize(im_label_4,1/scale_4,'bicubic'),[hei_4,wid_4],'bicubic');
%     imshow([im_label_4 im_input_4],[]);

    for x = 1 : stride : hei_4-size_input+1
        for y = 1 :stride : wid_4-size_input+1

            subim_input = im_input_4(x : x+size_input-1, y : y+size_input-1);
            subim_label = im_label_4(x+padding : x+padding+size_label-1, y+padding : y+padding+size_label-1);

            count=count+1;
            data_temp = subim_input;
            label_temp = subim_label;

            data(:, :, :, count) = single(data_temp);
            label(:, :, :, count) = single(label_temp);
        end
    end

end

order = randperm(count);
data = data(:, :, :, order);
label = label(:, :, :, order);

%% writing to HDF5
chunksz = 64;
created_flag = false;
totalct = 0;

for batchno = 1:floor(count/chunksz)
    last_read=(batchno-1)*chunksz;
    batchdata = data(:,:,:,last_read+1:last_read+chunksz);
    batchlabs = label(:,:,:,last_read+1:last_read+chunksz);

    startloc = struct('dat',[1,1,1,totalct+1], 'lab', [1,1,1,totalct+1]);
    curr_dat_sz = store2hdf5(savepath, batchdata, batchlabs, ~created_flag, startloc, chunksz);
    created_flag = true;
    totalct = curr_dat_sz(end);
end
h5disp(savepath);

author: tyshiwo


text:
function La = YeePattanaikLuminanceAdaptation(img, maxLayers)
%
%
%       La = YeePattanaikLuminanceAdaptation(img, maxLayers)
%
%
%       Input:
%           -img: HDR image
%           -maxLayers: the number of layers in [16,96]
%
%       Output:
%           -La: luminance adaptation
%
%     Copyright (C) 2010-2020  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%     The paper describing this technique is:
%     "Segmentation and Adaptive Assimilation for Detail-Preserving Display of High-Dynamic Range Images"
% 	  by Hector Yee, Sumanta N. Pattanaik
%     in Elsevier The Visual Computer 2003
%

%is it a three color channels image?
check13Color(img);

checkNegative(img);

%check parameters
if(~exist('maxLayers', 'var'))
    maxLayers = 32;
end

if((maxLayers < 16) || (maxLayers > 96))
    maxLayers = 32;
end

maxLayers = round(maxLayers);

%these could be parameters
bin_size1 = 0.5;
bin_size2 = 2.0;

%compute luminance channel
L = lum(img);

%compute the adaptation
L_log = log10(L + 1e-6);
minL_Log = min(L_log(:));

La = zeros(size(L));

for i=0:(maxLayers - 1)
    bin_size = bin_size1 + (bin_size2 - bin_size1) * i / (maxLayers - 1);
    category = round((L_log - minL_Log) / bin_size) + 1;

    %compute layers
    [imgLabel, ~, ~] = computeConnectedComponents(category, 8);
    labels = unique(imgLabel);

    for j=1:length(labels) %adaptation group
        indx = find(imgLabel == labels(j));
        La(indx) = La(indx) + mean(L_log(indx));
    end
end

La = 10.^(La / maxLayers);
La(La < 0.0) = 0.0;

end

author: banterle


text:
classdef Camera_Rotatory < Camera_Abstract

    properties (Access = protected)
        currentView = [90 30]
    end

    properties (Access = private)
        period = 200
        speed
    end

    methods (Access = public)

        function obj = Camera_Rotatory(axes)
            obj@Camera_Abstract(axes);
        end

        function updateView(obj)
            obj.applyAngularMotion();
            obj.checkView();
            obj.refresh();
        end

    end

    methods (Access = private)

        function applyAngularMotion(obj)
            obj.currentView  = obj.currentView + [obj.speed 0];
        end

        function checkView(obj)
            az = obj.currentView(1);
            el = obj.currentView(2);

            az = obj.keepInRange(az,0,360);
            el = obj.keepInRange(el,-90,90);

            obj.currentView(1) = az;
            obj.currentView(2) = el;

        end

    end

    methods (Access = private, Static)

        function x = keepInRange(x,xMin,xMax)
            X = abs(xMin-xMax);
            while x > xMax
                x = x - X;
            end

            while x < xMin
                x = x + X;
            end
        end

    end

    methods

        function w = get.speed(obj)
            w = 360/obj.period;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [Abar, Bbar, Cbar] = genCompMat(phi, gamma, lambda, Np, Nc)
% genCompMat Generate composite matrices for FTCOCP
%
% Author : Ajinkya Khade, askhade@ncsu.edu

Abar = cell(Np, 1);
Bbar = cell(Np, Nc);

b0 = zeros(size(gamma));

for i = 1:Np
    Abar{i} = phi^i;

    for j = 1:Nc
        if i >= j
            Bbar{i,j} = (phi^(i-j))*gamma;
        else
            Bbar{i,j} = b0;
        end
    end
end

Abar = cell2mat(Abar);
Bbar = cell2mat(Bbar);
Cbar = kron(eye(Np),lambda);

end

author: ajinkya-khade


text:
function ac = lpcweight(ar,c)
%  lpcweight --> LPC based perceptual weighting filter.
%
%    ac = lpcweight(ar,c)
%
%    The function takes the LP coefficients, ar = [1 -a(1) ... -a(M)],
%    and the parameter, c, as inputs, and returns the coefficients of
%    the filter function A(z/c) in the vector ac.

% Linear predictor order.
M = length(ar);

% The i'th coefficient of A(z/c) is given by ar(i)*c^(i-1).
ac = ar;
ci = c;
for (i=2:M)
  ac(i) = ar(i)*ci;
  ci = ci*c;
end

author: Sable

repo: mcbench-benchmarks


text:
%sim_spoil.m
%Jamie Near, 2014.
%
% USAGE:
% d_out = sim_spoil(d_in,H,angle)
%
% DESCRIPTION:
% This function simulates the effect of a rotation about the z-axis.
%
% INPUTS:
% d_in      = input density matrix structure.
% H         = Hamiltonian operator structure.
% angle     = Spoil angle (degrees).
%
% OUTPUTS:
% d_out     = output density matrix following z-rotation.

function d_out = sim_spoil(d_in,H,angle)

for m=1:length(H)
    %make vector of spoil angles;
    angle=angle*ones(length(H(m).shifts),1);

    %Make spoiler hamiltonian;
    spoil=zeros(2^H(m).nspins,2^H(m).nspins);
    for n=1:H(m).nspins
        spoil=spoil+((angle(n)*pi/180)*H(m).Iz(:,:,n));
    end

    %Do matrix multiplication;
    p=expm(1i*spoil);
    d_out{m} = p' * d_in{m} * p;
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function ROIid = getROIdbId(view)
% function ROIid = getROIdbId(view)
% User should select a ROI from a list (a list of all
% ROIs available in the DB).
% Or maybe we must show ROIs for the current session only?..

OKflag = openDbConnection(view);
if(~OKflag)
  ROIid = 0;
  return;
end;

[IdList,NameList,SessionIdList] = ...
  mysql('SELECT rois.id,rois.ROIname,rois.sessionid FROM rois');
[KnownSessionsId,KnownSessionsName] = ...
  mysql('SELECT sessions.id, sessions.sessionCode FROM sessions,rois WHERE sessions.id=rois.sessionid');
mysql('close');

for(Q=1:length(IdList))
  if(SessionIdList(Q)~=0)
    NameList{Q} = ['Session [' KnownSessionsName{max(find(KnownSessionsId==SessionIdList(Q)))} ']: ' NameList{Q}];
  else
    NameList{Q} = ['Unknown session: ' NameList{Q}];
  end
end

[ROIid,OK] = selector(IdList,NameList,'Choose a ROI');

if(OK==0)
  ROIid = 0;
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function omega = angle_outer(q1,q2,varargin)
% calcualtes the rotational angle between all rotations q1 and q2
%
% Syntax
%   omega = angle_outer(q1,q2)
%
% Input
%  q1, q2 - @quaternion
%
% Output
%  omega  - double [q1 x q2]

if nargin == 1, q2 = q1; end

omega = real(2*acos(abs(dot_outer(q1,q2,varargin{:}))));

author: mtex-toolbox

repo: mtex

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


text:
%io_loadspec_sdat.m
%Jamie Near, McGill University 2014.
%Georg Oeltzschner, Johns Hopkins University 2018.
%
% USAGE:
% out=io_loadspec_sdat(filename,subspecs);
%
% DESCRIPTION:
% Reads in Philips MRS data (.spar and .sdat files) using code adapted from
% PhilipsRead.m, provided as part of the Gannet software package by Richard
% Edden (gabamrs.com).
%
% io_loadspec_sdat outputs the data in structure format, with fields
% corresponding to time scale, fids, frequency scale, spectra, and header
% fields containing information about the acquisition.  The resulting
% matlab structure can be operated on by the other functions in this MRS
% toolbox. ALSO:  This code is not currently smart enough to parse out all
% of the relevant information from the header file, such as the number of
% subspectra.  So for now, these details must be passed to the function as
% input arguments.  Help implementing these improvements are most welcome!!
%
% INPUTS:
% filename   = filename of Philips sdat file to be loaded.
% subspecs   = number of subspectra in the data (from spectral editing, ISIS, etc.)
%
% OUTPUTS:
% out        = Input dataset in FID-A structure format.

function out = io_loadspec_sdat(filename,subspecs)

% Read in the data and header information
[data, header] = philipsLoad(filename);

% Determine the dimensions of the data
data_size       = size(data);
dims.t          = find(data_size == header.samples);
dims.averages   = find(data_size == header.rows);
dims.coils      = 0; % SDAT is already coil-combined
% Now arrange in the standard order (samples-avgs-subspecs):
data = permute(data ,[dims.t dims.averages]);
dims.t = 1;
dims.averages = 2;
dims.extras = 0;

% We have no way of actually knowing the number of sub-spectra (e.g. for
% MEGA-PRESS or HERMES acquisitions, so we will split the averages
% according to the 'subspecs' input.
% Initialize fids array:
fids = squeeze(zeros(header.samples, header.rows/subspecs, subspecs));
if subspecs == 2
    %Split the subspectra out of the "averages" dimension:
    fids(:,:,1) = data(:,[1:2:end]);
    fids(:,:,2) = data(:,[2:2:end]);
elseif subspecs == 4
    fids(:,:,1) = data(:,[1:4:end]);
    fids(:,:,2) = data(:,[2:4:end]);
    fids(:,:,3) = data(:,[3:4:end]);
    fids(:,:,4) = data(:,[4:4:end]);
else
    fids = data;
end

if subspecs > 1
    dims.subSpecs = 3;
else
    dims.subSpecs = 0;
end

sz = size(fids);

% Fill in the header information
txfrq = header.synthesizer_frequency; % transmitter frequency [Hz]
Bo = txfrq/42577000; % B0 [T]
averages = size(fids,2)*size(fids,3); % number of rows in the file
rawAverages = averages;
spectralwidth = header.sample_frequency; % bandwidth [Hz]
dwelltime = 1/spectralwidth; % dwelltime [s]
te = header.echo_time; % echo time [ms]
tr = header.repetition_time; % repetition time [ms]
sequence = header.scan_id; % sequence type
Ncoils = 1; % SDAT already coil-combined
rawSubspecs = subspecs;
date = header.scan_date; % scan date [yyyy.mm.dd hh:mm:ss]

% Produce specs
specs = fftshift(ifft(fids,[],dims.t),dims.t);
% Calculate t and ppm arrays using the calculated parameters:
f = [(-spectralwidth/2)+(spectralwidth/(2*sz(1))):spectralwidth/(sz(1)):(spectralwidth/2)-(spectralwidth/(2*sz(1)))];
ppm = -f/(Bo*42.577);
ppm = ppm+4.65;
t = [0:dwelltime:(sz(1)-1)*dwelltime];

%FILLING IN DATA STRUCTURE
out.fids=fids;
out.specs=specs;
out.sz=sz;
out.ppm=ppm;
out.t=t;
out.spectralwidth=spectralwidth;
out.dwelltime=dwelltime;
out.txfrq=txfrq;
out.date=date;
out.dims=dims;
out.Bo=Bo;
out.averages=averages;
out.rawAverages=rawAverages;
out.subspecs=subspecs;
out.rawSubspecs=rawSubspecs;
out.seq=sequence;
out.te=te;
out.tr=tr;
out.pointsToLeftshift=0;

%FILLING IN THE FLAGS
out.flags.writtentostruct=1;
out.flags.gotparams=1;
out.flags.leftshifted=0;
out.flags.filtered=0;
out.flags.zeropadded=0;
out.flags.freqcorrected=0;
out.flags.phasecorrected=0;
out.flags.averaged=0;
out.flags.addedrcvrs=1;
out.flags.subtracted=0;
out.flags.writtentotext=0;
out.flags.downsampled=0;
if out.dims.subSpecs==0
    out.flags.isFourSteps=0;
else
    out.flags.isFourSteps=(out.sz(out.dims.subSpecs)==4);
end

%DONE


text:
function [D,G,B] = autoGen_acrobotDynamics(q1,q2,dq1,dq2,m1,m2,g,l1,l2)
%AUTOGEN_ACROBOTDYNAMICS
%    [D,G,B] = AUTOGEN_ACROBOTDYNAMICS(Q1,Q2,DQ1,DQ2,M1,M2,G,L1,L2)

%    This function was generated by the Symbolic Math Toolbox version 6.2.
%    12-Jun-2015 16:56:47

t2 = cos(q1);
t3 = l1.^2;
t4 = sin(q1);
t5 = cos(q2);
t6 = l1.*t2;
t7 = l2.*t5;
t8 = t6+t7;
t9 = sin(q2);
t10 = l1.*t4;
t11 = l2.*t9;
t12 = t10+t11;
t13 = l2.^2;
D = reshape([-m1.*t2.^2.*t3-m1.*t3.*t4.^2-l1.*m2.*t2.*t8-l1.*m2.*t4.*t12,-l1.*l2.*m2.*t2.*t5-l1.*l2.*m2.*t4.*t9,-l2.*m2.*t5.*t8-l2.*m2.*t9.*t12,-m2.*t5.^2.*t13-m2.*t9.^2.*t13],[2, 2]);
if nargout > 1
    t14 = dq1.^2;
    t15 = dq2.^2;
    t16 = l1.*t2.*t14;
    t17 = l2.*t5.*t15;
    t18 = t16+t17;
    t19 = l1.*t4.*t14;
    t20 = l2.*t9.*t15;
    t21 = t19+t20;
    G = [-g.*m2.*t12+m2.*t8.*t21-m2.*t12.*t18-g.*l1.*m1.*t4;-g.*l2.*m2.*t9+l2.*m2.*t5.*t21-l2.*m2.*t9.*t18];
end
if nargout > 2
    B = [0.0;-1.0];
end

author: MatthewPeterKelly

repo: dscTutorials


text:
function u = mldivide(L, f, varargin)
%\         Solve a linear system (same as LINOP.LINSOLVE).
%   Important: A CHEBOPPREF object PREFS has to be passed. When this method
%   is called via CHEBOP/MLDIVIDE, PREFS is inherited from the CHEBOP level.
%
%  Copyright 2017 by The University of Oxford and The Chebfun Developers.
%  See http://www.chebfun.org/ for Chebfun information.
u = linsolve(L, f, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function T = conj(T)
% conjugate of a tensor
%
% Input
%  T - @tensor
%
% Output
%  T - @tensor
%

switch T.rank

  case 1

  case 2

    T.M = conj(T.M);

  case 3

  case 4

    % check for symmetry

    % convert to a matrix
    M = tensor42(T.M);

    % invert the matrix
    M = conj(M);

    % make some corrections
    w = 1./(1+((1:6)>3));
    w = w.' * w;
    M = M .* w;

    % convert to back a 4 rank tensor
    T.M = tensor24(M);
end


text:
function [h,b] = montage(varargin)
%MONTAGE Display multiple image frames as rectangular montage.
%   MONTAGE displays all the frames of a multiframe image array in a single
%   image object, arranging the frames so that they roughly form a square.
%
%   MONTAGE(I) displays the K frames of the intensity image array I. I is
%   M-by-N-by-1-by-K.
%
%   MONTAGE(BW) displays the K frames of the binary image array BW. BW is
%   M-by-N-by-1-by-K.
%
%   MONTAGE(X,MAP) displays the K frames of the indexed image array X, using the
%   colormap MAP for all frames. X is M-by-N-by-1-by-K.
%
%   MONTAGE(RGB) displays the K frames of the truecolor image array RGB. RGB is
%   M-by-N-by-3-by-K.
%
%   H = MONTAGE(...) returns the handle to the image object.
%
%   Class support
%   -------------
%   An intensity image can be logical, uint8, uint16, int16, single, or double.
%   An indexed image can be logical, uint8, uint16, single, or double. The map
%   must be double. A truecolor image can be uint8, uint16, single, or
%   double. The output is a handle to the graphics objects produced by
%   this function.
%
%   Example
%   -------
%       load mri
%       montage(D,map)
%
%   See also IMMOVIE.

%   Copyright 1993-2005 The MathWorks, Inc.
%   $Revision: 1.1.8.4 $  $Date: 2005/11/15 01:03:29 $

[a,cm,aspectRatio] = parse_inputs(varargin{:});

[nRows, nCols, nBands, nFrames] = size(a);

% Estimate nMontageColumns and nMontageRows given the desired ratio of
% Columns to Rows to be one (square montage).
nMontageCols = sqrt(aspectRatio * nRows * nFrames / nCols);

% Make sure montage rows and columns are integers. The order in the adjustment
% matters because the montage image is created horizontally across columns.
nMontageCols = ceil(nMontageCols);
nMontageRows = ceil(nFrames / nMontageCols);

% Create the montage image.
b = a(1,1); % to inherit type
b(1,1) = 0; % from a
b = repmat(b, [nMontageRows*nRows, nMontageCols*nCols, nBands, 1]);

rows = 1 : nRows;
cols = 1 : nCols;

for i = 0:nMontageRows-1
  for j = 0:nMontageCols-1,
    k = j + i * nMontageCols + 1;
    if k <= nFrames
      b(rows + i * nRows, cols + j * nCols, :) = a(:,:,:,k);
    else
      break;
    end
  end
end

if isempty(cm)
  hh = imshow(b);
else
  hh = imshow(b,cm);
end

if nargout > 0
    h = hh;
end

%--------------------------------------------------------------
%Parse Inputs Function

function [I,map,aspectRatio] = parse_inputs(varargin)

% initialize variables
map = [];
aspectRatio = 1;

iptchecknargin(1,4,nargin,mfilename);
iptcheckinput(varargin{1},{'uint8' 'double' 'uint16' 'logical' 'single' ...
                    'int16'},{},mfilename, 'I, BW, or RGB',1);
I = varargin{1};

if nargin>=2
  if isa(I,'int16')
    eid = sprintf('Images:%s:invalidIndexedImage',mfilename);
    msg1 = 'An indexed image can be uint8, uint16, double, single, or ';
    msg2 = 'logical.';
    error(eid,'%s %s',msg1, msg2);
  end
  if nargin >= 3 & isstr(varargin{2}),
    if strcmp(varargin{2},'aspectratio'),
      aspectRatio = varargin{3};
    end;
  else,
    map = varargin{2};
    iptcheckinput(map,{'double'},{},mfilename,'MAP',1);
    if ((size(map,1) == 1) && (prod(map) == numel(I)))
      % MONTAGE(D,[M N P]) OBSOLETE
      eid = sprintf('Images:%s:obsoleteSyntax',mfilename);
      msg1 = 'MONTAGE(D,[M N P]) is an obsolete syntax.';
      msg2 = 'Use multidimensional arrays to represent multiframe images.';
      error(eid,'%s\n%s',msg1,msg2);
    end
  end
end;

author: kristinbranson

repo: JAABA


text:
function tuned = TuneArtefactTimes(x,times,varargin)

%TuneArtefactTimes - Fine-tune artefact (e.g. stimulation) times.
%
%  Fine tuning can be done using either the LFP (faster), or the wideband
%  data (slower but more accurate). The latter is especially useful in cases
%  where the artefacts are no longer visible in the LFP signal due to low-
%  pass filtering.
%
%  USAGE
%
%    % Using LFP (pass the actual samples)
%    times = TuneArtefactTimes(lfp,artefacts,<options>)
%
%    % Using wideband data (pass the channel ID because the data would not fit
%    % in memory, so it will have to be read piecewise from disk)
%    times = TuneArtefactTimes(channel,artefacts,<options>)
%
%    lfp            local field potentials
%    channel        recording channel (wideband data)
%    artefacts      approximate artefact times
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'method'      artefacts can be detected using either the derivative
%                   ('slope', default) or the amplitude ('amplitude') of the
%                   signal
%     'durations'   interval (in s) around each approximate artefact time
%                   where the actual artefact should be sought
%                   (default = [-0.0075 0.0075])
%     'verbose'     display information about ongoing processing
%                   (default = 'off')
%    =========================================================================
%
%  SEE
%
%    See also RemoveArtefacts.

% 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.

% Default values
durations = [-0.0075 0.0075];
verbose = 'off';
nChunks = 100;
method = 'slope';

% Check number of parameters
if nargin < 2 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help TuneArtefactTimes'' for details).');
end

% LFP or wideband?
if isiscalar(x),
	source = 'wideband';
	channel = x;
else
	source = 'lfp';
	lfp = x;
	if size(lfp,2) ~= 2,
		error('Parameter ''lfp'' is not a Nx2 matrix (type ''help TuneArtefactTimes'' for details).');
	end
end

% Check parameter sizes
if ~isdvector(times),
	error('Parameter ''times'' is not a vector (type ''help TuneArtefactTimes'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+1) ' is not a property (type ''help TuneArtefactTimes'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'durations',
			durations = varargin{i+1};
			if ~isdvector(durations,'#2','>'),
			error('Incorrect value for property ''durations'' (type ''help TuneArtefactTimes'' for details).');
			end
		case 'method',
			method = lower(varargin{i+1});
			if ~isstring_FMAT(method,'slope','amplitude'),
				error('Incorrect value for property ''method'' (type ''help Sync'' for details).');
			end
		case 'verbose',
			verbose = lower(varargin{i+1});
			if ~isstring_FMAT(verbose,'on','off'),
				error('Incorrect value for property ''verbose'' (type ''help Sync'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help TuneArtefactTimes'' for details).']);

	end
end

% Discard duplicates
cleaned = unique(times);
if length(times) ~= length(cleaned),
	warning('Data contains duplicates - these will be discarded.');
end
times = cleaned;

if strcmp(source,'lfp'),
	% Process LFP all at once
	% Time intervals where the artefacts should be sought
	intervals = [times+durations(1) times+durations(2)];
	tuned = process(lfp,times,intervals,verbose);
else
	% Split data in chunks
	start = times(1) + durations(1);
	stop = times(end) + durations(2);
	totalTime = stop-start;
	dt = totalTime/nChunks;
	tuned = [];
	for i = 1:nChunks,
		% Time intervals where the artefacts should be sought
		t = Restrict(times,[start start+dt]);
		intervals = [t+durations(1) t+durations(2)];
		data = GetWidebandData(channel,'intervals',[start start+dt]);
		tuned = [tuned;process(data,t,intervals,method,verbose)];
		start = start + dt;
	end
end

% Actual processing

function tuned = process(signal,times,intervals,method,verbose)

% Differentiate signal
if strcmp(method,'slope'),
	signal(:,2) = [0;diff(signal(:,2))];
end

% Select the samples within these intervals
[in,i,j] = InIntervals(signal(:,1),intervals,'verbose',verbose);
% Transform this into a matrix where each line contains the data for one artefact, and find minima
m = full(sparse(i(i~=0),j(j~=0),signal(in,2)));
infs = ~logical(full(sparse(i(i~=0),j(j~=0),1)));
m(infs) = inf;
[unused,where] = min(m,[],2);
% Build the same matrix for data timestamps
t = full(sparse(i(i~=0),j(j~=0),signal(in,1)));
% Times
s1 = size(m,1);
s2 = size(m,2);
where = logical(full(sparse(1:s1,where,1,s1,s2)));
tuned = sort(t(where));

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [yd3] = m32yd3(m3)
% Convert volume from cubic meters to cubic yards.
% Chad Greene 2012
yd3 = m3*1.3079506193;


text:
%% (Internal) Calculate the performance data of a QRS detection task.
% This function is used by the ECGtask related with QRS detection to
% calculate performance of detectors.
%
% Example
%
%   payload_out = CalculatePerformanceECGtaskQRSdet(payload_out, ECG_annotations, ECG_header);
%
% Arguments:
%     +payload_out: Internal payload with all detections
%
%     +ECG_annotations: Gold standard to calculate performances.
%
%     +ECG_header: ECG header information.
%
%     +ECG_start_offset: Offset from ECG start.
%
% Output:
%     + payload: payloads generated by ECGtask_QRS_Detection.Process result
%        of concatenating plA and plB.
%
% See also ECGtask_QRS_Detection
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 06/10/2015
% Birthdate  : 06/10/2015
% Copyright 2008-2015
function payload_out = CalculatePerformanceECGtaskQRSdet(payload_out, ECG_annotations, ECG_header, ECG_start_offset)

    if( isempty(payload_out) || isempty(ECG_annotations) )
        return
    end

    AnnNames = payload_out.series_quality.AnnNames(:,1);
    cant_lead_name = size(AnnNames,1);
    payload_out.series_performance.conf_mat = zeros(2,2,cant_lead_name);
    payload_out.series_performance.error = nan(cant_lead_name,2);

    if(isempty(ECG_annotations))
        disp_string_framed(2, sprintf('Trusted references not found for %s', ECG_header.recname) );
    else
        % offset refs, produced anns were already shifted
        %ECG_annotations.time = ECG_annotations.time + ECG_start_offset - 1;
        ECG_annotations.time = ECG_annotations.time(ECG_annotations.time >= ECG_start_offset & ...
            ECG_annotations.time <= ECG_header.nsamp + ECG_start_offset);

        payload_out.series_performance.conf_mat_details = cell(cant_lead_name,4);

        for kk = 1:cant_lead_name
            this_lead_det = payload_out.(AnnNames{kk}).time;
            [payload_out.series_performance.conf_mat(:,:,kk), ...
                TP_gs_idx, ...
                TP_det_idx, ...
                FN_idx, ...
                FP_idx ] = bxb(ECG_annotations, this_lead_det, ECG_header );

            % details about the confusion matrix
            payload_out.series_performance.conf_mat_details(kk,:) = { TP_gs_idx, TP_det_idx, FN_idx, FP_idx };

            % meadian/mad of the error wrt the gold standard
            if( ~isempty(TP_det_idx) )
                aux_val = (colvec(this_lead_det(TP_det_idx)) - colvec(ECG_annotations.time(TP_gs_idx)));
                payload_out.series_performance.error(kk,:) = [ median(aux_val) mad(aux_val) ] /ECG_header.freq;
            end
        end
    end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function test_tutorial_clusterpermutationtimelock(dataset, datadir)

% MEM 3gb
% WALLTIME 00:10:00
% DEPENDENCY ft_timelockanalysis ft_multiplotER ft_singleplotER ft_timelockstatistics ft_topoplotER ft_clusterplot ft_megplanar ft_combineplanar

if nargin<1 || isempty(dataset)
  dataset = dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.ds');
end
if nargin<2
  datadir = dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/cluster_permutation_timelock');
end

cfg                         = [];
cfg.dataset                 = dataset;
cfg.trialfun                = 'ft_trialfun_general'; % this is the default
cfg.trialdef.eventtype      = 'backpanel trigger';
cfg.trialdef.eventvalue     = [3 5 9]; % the values of the stimulus trigger for the three conditions
% 3 = fully incongruent (FIC), 5 = initially congruent (IC), 9 = fully congruent (FC)
cfg.trialdef.prestim        = 1; % in seconds
cfg.trialdef.poststim       = 2; % in seconds

cfg = ft_definetrial(cfg);

% remove the trials that have artifacts from the trl
cfg.trl([2, 5, 6, 8, 9, 10, 12, 39, 43, 46, 49, 52, 58, 84, 102, 107, 114, 115, 116, 119, 121, 123, 126, 127, 128, 133, 137, 143, 144, 147, 149, 158, 181, 229, 230, 233, 241, 243, 245, 250, 254, 260],:) = [];

% preprocess the data
cfg.channel    = {'MEG', '-MLP31', '-MLO12'};        % read all MEG channels except MLP31 and MLO12
cfg.demean     = 'yes';
cfg.baselinewindow  = [-0.2 0];
cfg.lpfilter   = 'yes';                              % apply lowpass filter
cfg.lpfreq     = 35;                                 % lowpass at 35 Hz.

data_all = ft_preprocessing(cfg);

cfg        = [];
cfg.trials = data_all.trialinfo == 3;
dataFIC_LP = ft_redefinetrial(cfg, data_all);

cfg = [];
cfg.trials = data_all.trialinfo == 9;
dataFC_LP = ft_redefinetrial(cfg, data_all);

%save dataFIC_LP dataFIC_LP
%save dataFC_LP dataFC_LP

cfg = [];
cfg.keeptrials = 'yes';
timelockFIC    = ft_timelockanalysis(cfg, dataFIC_LP);
timelockFC     = ft_timelockanalysis(cfg, dataFC_LP);

%% ## Permutation test

cfg                  = [];
cfg.method           = 'montecarlo'; % use the Monte Carlo Method to calculate the significance probability
cfg.statistic        = 'indepsamplesT'; % use the independent samples T-statistic as a measure to
                                   % evaluate the effect at the sample level
cfg.correctm         = 'cluster';
cfg.clusteralpha     = 0.05;       % alpha level of the sample-specific test statistic that
                                   % will be used for thresholding
cfg.clusterstatistic = 'maxsum';   % test statistic that will be evaluated under the
                                   % permutation distribution.
cfg.minnbchan        = 2;          % minimum number of neighborhood channels that is
                                   % required for a selected sample to be included
                                   % in the clustering algorithm (default=0).
% cfg.neighbours     = neighbours; % see below
cfg.tail             = 0;          % -1, 1 or 0 (default = 0); one-sided or two-sided test
cfg.clustertail      = 0;
cfg.alpha            = 0.025;      % alpha level of the permutation test
cfg.numrandomization = 100;        % number of draws from the permutation distribution

n_fc  = size(timelockFC.trial, 1);
n_fic = size(timelockFIC.trial, 1);

cfg.design           = [ones(1,n_fic), ones(1,n_fc)*2]; % design matrix
cfg.ivar             = 1; % number or list with indices indicating the independent variable(s)

cfg_neighb        = [];
cfg_neighb.method = 'distance';
neighbours        = ft_prepare_neighbours(cfg_neighb, dataFC_LP);

cfg.neighbours    = neighbours;  % the neighbours specify for each sensor with
                                 % which other sensors it can form clusters
cfg.channel       = {'MEG'};     % cell-array with selected channel labels
cfg.latency       = [0 1];       % time interval over which the experimental
                                 % conditions must be compared (in seconds)
[stat] = ft_timelockstatistics(cfg, timelockFIC, timelockFC);

% Save the output to disk:
%save stat_ERF_axial_FICvsFC stat;

%% ## Plotting of the results
%
% To plot the results of the permutation test, we use the plotting function **[ft_topoplotER](https://github.com/fieldtrip/fieldtrip/blob/release/ft_topoplotER.m)**. In doing so, we will plot a topography of the difference between the two experimental conditions (FIC and FC). On top of that, and for each timestep of interest, we will highlight the sensors which are members of significant clusters. First, however, we must calculate the difference between conditions using **[ft_math](https://github.com/fieldtrip/fieldtrip/blob/release/ft_math.m)**.
%
cfg    = [];
avgFIC = ft_timelockanalysis(cfg, dataFIC_LP);
avgFC  = ft_timelockanalysis(cfg, dataFC_LP);

% Then take the difference of the averages using ft_math
cfg           = [];
cfg.operation = 'subtract';
cfg.parameter = 'avg';
raweffectFICvsFC = ft_math(cfg, avgFIC, avgFC);

% We then construct a boolean matrix indicating whether a channel/time point belongs to a cluster that we deem interesting to inspect. This matrix has size [Number_of_MEG_channels x Number_of_time_samples], like stat.posclusterslabelmat. We'll make two such matrices: one for positive clusters (named pos), and one for negative (neg). All (channel,time)-pairs belonging to the large clusters whose probability of occurrence is sufficiently low in relation to the associated randomization distribution of clusterstats will be coded in the new boolean matrix as 1, and all those that don't will be coded as 0.
%
%
pos_cluster_pvals = [stat.posclusters(:).prob];

% Then, find which clusters are deemed interesting to visualize, here we use a cutoff criterion based on the
% cluster-associated p-value, and take a 5% two-sided cutoff (i.e. 0.025 for the positive and negative clusters,
% respectively
pos_clust = find(pos_cluster_pvals < 0.025);
pos       = ismember(stat.posclusterslabelmat, pos_clust);

% and now for the negative clusters...
neg_cluster_pvals = [stat.negclusters(:).prob];
neg_clust         = find(neg_cluster_pvals < 0.025);
neg               = ismember(stat.negclusterslabelmat, neg_clust);

% Alternatively, we can manually select which clusters we want to plot. If we only want to see the extext of the first (i.e. most significant) positive and negative clusters, for instance, we can do so as follows:
%
pos = stat.posclusterslabelmat == 1; % or == 2, or 3, etc.
neg = stat.negclusterslabelmat == 1;

% To plot a sequence of twenty topographic plots equally spaced between 0 and 1 second, we define the vector j of time steps. These time intervals correspond to the samples m in stat and in the variables pos and neg. m and j must, therefore, have the same length.
%
% To be sure that your sample-based time windows align with your time windows in seconds, check the following:
%
timestep      = 0.05; % timestep between time windows for each subplot (in seconds)
sampling_rate = dataFC_LP.fsample; % Data has a temporal resolution of 300 Hz
sample_count  = length(stat.time);
% number of temporal samples in the statistics object
j = [0:timestep:1]; % Temporal endpoints (in seconds) of the ERP average computed in each subplot
m = [1:timestep*sampling_rate:sample_count]; % temporal endpoints in M/EEG samples

% To plot the data use the following for-loop:
%
% First ensure the channels to have the same order in the average and in the statistical output.
% This might not be the case, because ft_math might shuffle the order
[i1,i2] = match_str(raweffectFICvsFC.label, stat.label);

for k = 1:20
   subplot(4,5,k);
   cfg = [];
   cfg.xlim = [j(k) j(k+1)];   % time interval of the subplot
   cfg.zlim = [-2.5e-13 2.5e-13];
   % If a channel is in a to-be-plotted cluster, then
   % the element of pos_int with an index equal to that channel
   % number will be set to 1 (otherwise 0).

   % Next, check which channels are in the clusters over the
   % entire time interval of interest.
   pos_int = zeros(numel(raweffectFICvsFC.label),1);
   neg_int = zeros(numel(raweffectFICvsFC.label),1);
   pos_int(i1) = all(pos(i2, m(k):m(k+1)), 2);
   neg_int(i1) = all(neg(i2, m(k):m(k+1)), 2);

   cfg.highlight   = 'on';
   % Get the index of the to-be-highlighted channel
   cfg.highlightchannel = find(pos_int | neg_int);
   cfg.comment     = 'xlim';
   cfg.commentpos  = 'title';
   cfg.layout      = 'CTF151_helmet.mat';
   cfg.interactive = 'no';
   cfg.figure      = 'gca';
   ft_topoplotER(cfg, raweffectFICvsFC);
end

%% ## Using planar gradient data
%
% To perform the permutation test using synthetic planar gradient data, the data must first be converted using the functions **[ft_megplanar](https://github.com/fieldtrip/fieldtrip/blob/release/ft_megplanar.m)** and **[ft_combineplanar](https://github.com/fieldtrip/fieldtrip/blob/release/ft_combineplanar.m)**. These functions were described in the tutorial on event-related fields. After running these functions, the statistical analysis using **[ft_timelockstatistics ](https://github.com/fieldtrip/fieldtrip/blob/release/ft_timelockstatistics.m)** involves the same configuration options as for ordinary event-related averages (no synthetic planar gradients). There is only one additional step, which is needed to add the gradiometer structure to one of the planar gradient data sets.
%
cfg = [];
cfg.planarmethod   = 'sincos';
cfg.neighbours     = neighbours; % also here, neighbouring sensors needs to be defined
timelockFIC_planar = ft_megplanar(cfg, timelockFIC);
timelockFC_planar  = ft_megplanar(cfg, timelockFC);

timelockFIC_planar_cmb = ft_combineplanar(cfg, timelockFIC_planar);
timelockFC_planar_cmb  = ft_combineplanar(cfg, timelockFC_planar);

timelockFIC_planar_cmb.grad = timelockFIC.grad;  % add the gradiometer structure
timelockFC_planar_cmb.grad  = timelockFC.grad;

% Having calculated synthetic planar gradient data, one can use the same configuration parameters as used for the analysis of the original data.
%
cfg                  = [];
cfg.channel          = {'MEG'};
cfg.latency          = [0 1];
cfg.neighbours       = neighbours;
cfg.method           = 'montecarlo';
cfg.statistic        = 'indepsamplesT';
cfg.correctm         = 'cluster';
cfg.clusteralpha     = 0.05;
cfg.clusterstatistic = 'maxsum';
cfg.minnbchan        = 2;
cfg.tail             = 0;
cfg.clustertail      = 0;
cfg.alpha            = 0.025;
cfg.numrandomization = 100;

n_fc  = size(timelockFC_planar_cmb.trial, 1);
n_fic = size(timelockFIC_planar_cmb.trial, 1);

cfg.design           = [ones(1,n_fic), ones(1,n_fc)*2]; % design matrix
cfg.ivar             = 1; % number or list with indices indicating the independent variable(s)

cfg.ivar = 1;

[stat] = ft_timelockstatistics(cfg, timelockFIC_planar_cmb, timelockFC_planar_cmb);

%save stat_ERF_planar_FICvsFC stat

% The output can also be obtained from [ftp://ftp.fieldtriptoolbox.org/pub/fieldtrip/tutorial/cluster_permutation_timelock/stat_ERF_planar_FICvsFC.mat](ftp://ftp.fieldtriptoolbox.org/pub/fieldtrip/tutorial/cluster_permutation_timelock/stat_ERF_planar_FICvsFC.mat). If you need to reload the statistics output, use:
%
%load stat_ERF_planar_FICvsFC

% We now calculate the raw effect in the average with planar gradient data using the following configuration:
%
cfg = [];
cfg.keeptrials = 'no';   % now only the average, not the single trials
avgFIC_planar  = ft_timelockanalysis(cfg, timelockFIC_planar);
avgFC_planar   = ft_timelockanalysis(cfg, timelockFC_planar);
cfg = [];
avgFIC_planar_cmb = ft_combineplanar(cfg, avgFIC_planar);
avgFC_planar_cmb  = ft_combineplanar(cfg, avgFC_planar);

% subtract avgFC from avgFIC
cfg = [];
cfg.operation    = 'subtract';
cfg.parameter    = 'avg';
raweffectFICvsFC = ft_math(cfg, avgFIC_planar_cmb, avgFC_planar_cmb);

% Using the following configuration for **[ft_topoplotER](https://github.com/fieldtrip/fieldtrip/blob/release/ft_topoplotER.m)** we can plot the raw effect and highlight the channels contributing to the largest cluster
%
figure;
timestep      = 0.05; %(in seconds)
sampling_rate = dataFC_LP.fsample;
sample_count  = length(stat.time);
j = [0:timestep:1]; % Temporal endpoints (in seconds) of the ERP average computed in each subplot
m = [1:timestep*sampling_rate:sample_count]; % temporal endpoints in M/EEG samples

pos_cluster_pvals = [stat.posclusters(:).prob];
pos_clust         = find(pos_cluster_pvals < 0.025);
pos               = ismember(stat.posclusterslabelmat, pos_clust);

% Remember to do the same for negative clusters if you want them!

% First ensure the channels to have the same order in the average and in the statistical output.
% This might not be the case, because ft_math might shuffle the order
[i1,i2] = match_str(raweffectFICvsFC.label, stat.label);

for k = 1:20
   subplot(4,5,k);
   cfg = [];
   cfg.xlim =[j(k) j(k+1)];
   cfg.zlim = [-1.0e-13 1.0e-13];
   pos_int = zeros(numel(raweffectFICvsFC.label),1);
   pos_int(i1) = all(pos(i2, m(k):m(k+1)), 2);
   cfg.highlight = 'on';
   cfg.highlightchannel = find(pos_int);
   cfg.comment = 'xlim';
   cfg.commentpos = 'title';
   cfg.layout = 'CTF151_helmet.mat';
   cfg.figure = 'gca';
   ft_topoplotER(cfg, raweffectFICvsFC);
end

%% # Within-subjects experiments
load(fullfile(datadir, 'ERF_orig.mat'));

% The configuration looks as follows:
cfg         = [];
cfg.channel = {'MEG'};
cfg.latency = [0 1];

cfg.method           = 'montecarlo';
cfg.statistic        = 'depsamplesT';
cfg.correctm         = 'cluster';
cfg.clusteralpha     = 0.05;
cfg.clusterstatistic = 'maxsum';
cfg.minnbchan        = 2;
cfg.neighbours       = neighbours;  % same as defined for the between-trials experiment
cfg.tail             = 0;
cfg.clustertail      = 0;
cfg.alpha            = 0.025;
cfg.numrandomization = 500;

Nsubj  = 10;
design = zeros(2, Nsubj*2);
design(1,:) = [1:Nsubj 1:Nsubj];
design(2,:) = [ones(1,Nsubj) ones(1,Nsubj)*2];

cfg.design = design;
cfg.uvar   = 1;
cfg.ivar   = 2;

[stat] = ft_timelockstatistics(cfg, allsubjFIC{:}, allsubjFC{:});

%%save stat_ERF_planar_FICvsFC_GA stat

% The output can also be obtained from [stat_ERF_planar_FICvsFC_GA.mat](ftp://ftp.fieldtriptoolbox.org/pub/fieldtrip/tutorial/cluster_permutation_timelock/stat_ERF_planar_FICvsFC_GA.mat). If you need to reload the statistics output, use:
%load stat_ERF_planar_FICvsFC_GA

%% ## Plotting the results
% load individual subject data
load(fullfile(datadir, 'ERF_orig.mat'));

% calculate the grand average for each condition
cfg = [];
cfg.channel   = 'all';
cfg.latency   = 'all';
cfg.parameter = 'avg';
GA_FIC        = ft_timelockgrandaverage(cfg, allsubjFIC{:});
GA_FC         = ft_timelockgrandaverage(cfg, allsubjFC{:});

cfg = [];
cfg.operation = 'subtract';
cfg.parameter = 'avg';
GA_FICvsFC    = ft_math(cfg, GA_FIC, GA_FC);

figure;
% define parameters for plotting
timestep      = 0.05; %(in seconds)
sampling_rate = dataFIC_LP.fsample;
sample_count  = length(stat.time);
j = [0:timestep:1];   % Temporal endpoints (in seconds) of the ERP average computed in each subplot
m = [1:timestep*sampling_rate:sample_count];  % temporal endpoints in M/EEG samples

% get relevant values
pos_cluster_pvals = [stat.posclusters(:).prob];
pos_clust = find(pos_cluster_pvals < 0.025);
pos       = ismember(stat.posclusterslabelmat, pos_clust);

% First ensure the channels to have the same order in the average and in the statistical output.
% This might not be the case, because ft_math might shuffle the order
[i1,i2] = match_str(GA_FICvsFC.label, stat.label);

% plot
for k = 1:20
   subplot(4,5,k);
   cfg            = [];
   cfg.xlim       = [j(k) j(k+1)];
   cfg.zlim       = [-5e-14 5e-14];
   pos_int        = zeros(numel(GA_FICvsFC.label),1);
   pos_int(i1)    = all(pos(i2, m(k):m(k+1)), 2);
   cfg.highlight  = 'on';
   cfg.highlightchannel = find(pos_int);
   cfg.comment    = 'xlim';
   cfg.commentpos = 'title';
   cfg.layout     = 'CTF151_helmet.mat';
   cfg.figure     = 'gca';
   ft_topoplotER(cfg, GA_FICvsFC);
end


text:
%compute dsumspinerelanglehc

function [data,units]=compute_dsumspinerelanglehc(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
dsumspinerelanglehc=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    dsumspinerelanglehc{1,i}=(bsxfun(@minus,trx(larva).sumspinerelanglehc(2:end),trx(larva).sumspinerelanglehc(1:end-1)))./trx(larva).dt;
end

units=parseunits('rad/s');
data=dsumspinerelanglehc;


text:
classdef quaternion
%
% The class quaternion realizes the internal representation of rotations
% and orientations in MTEX. An important difference is that an unit
% quaternion *q* and its antipodal *-q* are considered as different while
% the corresponding rotations are equal.
%
% Syntax
%   q = quaternion(a,b,c,d)
%   q = quaternion.rand
%   q = quaternion.id
%   q = quaternion.nan
%
% Input
%  a, b, c, d - double
%
% Class Properties
%  a, b, c, d - quaternion coefficients
%
% See also
% rotation.rotation orientation.orientation

  properties
    a % real part
    b % * i
    c % * j
    d % * k
  end

  methods
    function q = quaternion(varargin)

      if nargin == 0, return;end

      if isa(varargin{1},'quaternion')   % copy constructor

        if nargin == 1
          q.a = varargin{1}.a;
          q.b = varargin{1}.b;
          q.c = varargin{1}.c;
          q.d = varargin{1}.d;
        else
          q.a = varargin{1}.a(varargin{2:end});
          q.b = varargin{1}.b(varargin{2:end});
          q.c = varargin{1}.c(varargin{2:end});
          q.d = varargin{1}.d(varargin{2:end});
        end
      elseif isa(varargin{1},'vector3d')

        q.a = zeros(size(varargin{1}));
        [q.b,q.c,q.d] = double(varargin{1});

      elseif isnumeric(varargin{1})

        switch nargin

          case 1

            D = varargin{1};

            q.a = D(1,:); q.b = D(2,:); q.c = D(3,:); q.d = D(4,:);

            s = size(D);
            s = [1 s(2:ndims(D))];

            q = reshape(q,s);

          case 2

            if length(varargin{1}) ~= 1
              q.a = varargin{1};
            else
              q.a = repmat(varargin{1},size(varargin{2}));
            end

            [q.b,q.c,q.d] = double(varargin{2});

          case 4

            q.a = varargin{1};
            q.b = varargin{2};
            q.c = varargin{3};
            q.d = varargin{4};

        end
      end

    end

    function d = quat_dot(g1,g2)
      d = g1.a .* g2.a + g1.b .* g2.b + g1.c .* g2.c + g1.d .* g2.d;
    end

    function d = quat_dot_outer(g1,g2,varargin)

      if ~isempty(g1) && ~isempty(g2)

        q1 = [g1.a(:) g1.b(:) g1.c(:) g1.d(:)];
        q2 = [g2.a(:) g2.b(:) g2.c(:) g2.d(:)];

        d = q1 * q2';

      else
        d = [];
      end
    end

    function n = numArgumentsFromSubscript(varargin)
      n = 0;
    end

  end

  methods (Static = true)

    function q = nan(varargin)
      a = nan(varargin{:});
      q = quaternion(a,a,a,a);
    end

    function q = id(varargin)
      a = ones(varargin{:});
      b = zeros(varargin{:});
      q = quaternion(a,b,b,b);
    end

    function q = rand(varargin)

      isNum = cellfun(@isnumeric,varargin);
      [~,last] = find(~isNum,1,'first');

      s = varargin;
      if ~isempty(last), s = s(1:last-1); end

      if check_option(varargin,'maxAngle')

        omega = linspace(0,get_option(varargin,'maxAngle'),1e4);

        v = vector3d.rand(s{:});

        id = discretesample(sin(omega).^2,prod([s{:}]));
        omega = reshape(omega(id),size(v));

        q = axis2quat(v,omega);

      else
        if length(s) < 2, s = [s 1]; end

        alpha = 2*pi*rand(s{:});
        beta  = acos(2*(rand(s{:})-0.5));
        gamma = 2*pi*rand(s{:});

        q = euler2quat(alpha,beta,gamma);
      end
    end

  end

end


text:
function out = getPolCoeffs(T,a,b,wf,N,q0)
% -----------------------------------------------------------------------
% The function computes coefficeints of the 5-th order polynomail
% from a and b coefficients of the finite furier series and
% initial position of the robot.
% Inputes:
%   T - period of motion
%   a - sine coeffs in finite fourier series
%   b - cosine coeffs in finite fourier series
%   wf - fundamental frequency
%   N - number of harmonics
%   q0 - initial position of the robot
% Output:
%   out - coefficients of the fifth order polynomail that guarantees
%         that initial and final constraints on position, velocities
%         and accelerations are satisfied
% -----------------------------------------------------------------------
qh0 = -sum(b./((1:N)*wf),2);
qdh0 = sum(a,2);
q2dh0 = sum(b.*(1:1:N)*wf,2);

[qhT,qdhT,q2dhT] = fourier_series_traj(T,zeros(6,1),a,b,wf,N);

I_6x6 = eye(6);
O_6x6 = zeros(6);

Ac = [I_6x6, O_6x6, O_6x6, O_6x6, O_6x6, O_6x6;
      O_6x6, I_6x6, O_6x6, O_6x6, O_6x6, O_6x6;
      O_6x6, O_6x6, 2*I_6x6, O_6x6, O_6x6, O_6x6;
      I_6x6, T*I_6x6, T^2*I_6x6, T^3*I_6x6, T^4*I_6x6, T^5*I_6x6;
      O_6x6, I_6x6, 2*T*I_6x6, 3*T^2*I_6x6, 4*T^3*I_6x6, 5*T^4*I_6x6;
      O_6x6, O_6x6, 2*I_6x6, 6*T*I_6x6, 12*T^2*I_6x6, 20*T^3*I_6x6];

% q0 = deg2rad([0    -90   0     -90   0     0 ]');  % !!!!!!!!!!!
c = Ac\[q0-qh0; -qdh0; -q2dh0; q0-qhT; -qdhT; -q2dhT];
out = reshape(c,[6,6]);

author: shamilmamedov


text:
classdef dme_load3p_opf < mp.dme_load3p & mp.dme_shared_opf

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

%     properties
%     end     %% properties

%     methods
%     end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function im = read_tiny_gist_binary(words_in,ind)
%
%%% Routine to read in tiny images from giant binary file holding all
% 384-dim gist descriptors.
%
%%% There are two modes of operation:
%
% 1. raw indexing mode
% --------------------
% Inputs:
% 1. words_in - 1 x nImages vector, each element being an integer
% in the range 1 to 73,777,893 (or whatever the total number of images in the dataset is).
% Outputs:
% 1. im - 384 x nImages single data of Gist descriptors.
%
% ++++++++++++++++++++++++++++++++++++++++++++++++++++
%
% 2. word based indexing.
% -----------------------
% Inputs:
% 1. word_in - 1 x nWords cell array, each entry holding one noun
% 2. ind - 1 x nWords cell array, each holding indices of images
% for the corresponding noun
% Outputs:
% 1. im - 384 x nImages single data of Gist descriptors.
%
% Requires: read_tiny_binary_gist_core.mex file
%
%%%% Hints regarding operation:
% The access time per image drops as the total number of images to be
% gathered increases. e.g. for 500,000 images load time is
% 1ms/image. For 100,000 images it is 2ms/image. So try and pass in big
% blocks of indices.

% path to giant binary file (currently ~114Gb)
binary_fname = '/mit/tiny/data/tinygist80million.bin';
% indexing info .mat file
binary_data_fname = '/mit/tiny/data/tiny_index.mat';

% total # imgs
TOTAL_NUM_IMGS = 79302017;

if (nargin==1)
    RAW_INDEX_MODE = 1;
else
    RAW_INDEX_MODE = 0;
end

if (nargin==2)
    if iscell(ind)
        % do nothing
    else
        % convert to cell arrays
        words_in = {words_in};
        ind = {ind};
    end

    nWords = length(words_in);
end

if RAW_INDEX_MODE
    %% just go and call binary with raw indices
    %% check that all indices are >0 and <TOTAL_NUM_IMGS
    if ((sum(words_in<1)==0) & (sum(words_in>TOTAL_NUM_IMGS)==0))

        % now sort for speed
        [sorted_words,sort_ind] = sort(words_in);

        %% use MEX for sppedy reading in
        im = read_tiny_binary_gist_core(binary_fname, uint64(sorted_words));

        % sort back
        im(:,sort_ind) = im;
    else
        error('indices <0 or greater than total number of images');
    end

else % words_in and indices passed in

    % load up index data
    load(binary_data_fname);

    % compute hashes for all words_in
    hashes = compute_hash_function(words_in,hash_key);

    img_offsets = [];

    for a=1:nWords

        % find word using hash
        match_ind = find(hashes(a) == hash);

        if any(match_ind)
            % check all ind are >=1 and <=total number of images

            if ((sum(ind{a}<1)==0) & (sum(ind{a}>num_imgs(match_ind))==0))
                word_offset = offset(match_ind);
                img_offsets = [img_offsets , uint64(double(word_offset) + [ind{a}])];
            else
                error('indices <0 or greater than number per image');
            end
        else
           error('word not found');
        end

        % now sort for speed
        [sorted_offsets,sort_ind] = sort(img_offsets);

        im = read_tiny_binary_gist_core(binary_fname, uint64(sorted_offsets));

        % sort back
        im(:,sort_ind) = im;
    end
end

author: huashiyiqike

repo: LSTM-MATLAB


text:
function [hhist,hax,hfig,hleg,hxlabel,hylabel,frac,frac_outside,edges,centers_plot] = ...
  HistogramPerFrameFeature(data,perframefn,varargin)

[hax,hfig,...
  minprctile,maxprctile,minv,maxv,...
  edges,nbins,...
  binmode,...
  labelcolors,unknowncolor,...
  frac,...
  frac_outside,...
  centers_plot] = ...
  myparse(varargin,'axes',[],'figure',[],...
  'minprctile',[],'maxprctile',[],...
  'minv',[],'maxv',[],...
  'edges',[],...
  'nbins',100,...
  'binmode','linear',...
  'labelcolors',[],...
  'unknowncolor',[1,1,1],...
  'frac',[],...
  'frac_outside',[],...
  'centers_plot',[]);

nbehaviors = data.nbehaviors;

if isempty(hax),
  if isempty(hfig),
    hfig = figure;
  end
  hax = gca;
end

holdstate = ishold(hax);
hold(hax,'off');

docompute = isempty(frac) || isempty(frac_outside) || ...
  isempty(edges) || isempty(centers_plot);

if docompute,

  perframeidx = find(strcmpi(perframefn,data.allperframefns),1);
  if isempty(perframeidx),
    error('Unknown perframe field %s',perframefn);
  end

  % use current fly to choose hist edges
  if isempty(edges),
    x = data.perframedata{perframeidx};
    x = x(~isnan(x));
    if isempty(minv),
      if isempty(minprctile),
        minv = min(x);
      else
        minv = prctile(x,minprctile);
      end
    end
    if isempty(maxv),
      if isempty(minprctile),
        maxv = max(x);
      else
        maxv = prctile(x,maxprctile);
      end
    end
    [edges,centers] = SelectHistEdges(nbins,[minv,maxv],binmode);
  else
    centers = (edges(1:end-1)+edges(2:end))/2;
    nbins = numel(centers);
  end

  counts = zeros(nbehaviors+1,nbins);
  counts_outside = zeros(nbehaviors+1,2);
  n = zeros(nbehaviors+1,1);

  for expi = 1:data.nexps,

    % load per-frame data for this experiment
    perframedir = data.GetFile('perframedir',expi);
    file = fullfile(perframedir,[perframefn,'.mat']);
    if ~exist(file,'file'),
      warning('Per-frame data file %s does not exist',file);
      continue
    end
    perframedata = load(file);

    % TODO: extend to multiple flies
    for fly = 1:data.nflies_per_exp(expi),

      x = perframedata.data{fly}(:)';

      % load labels for this fly and experiment
      if expi == data.expi && all(fly == data.flies),
        % use current labelidx
        labelidx = data.labelidx.vals;
      else
        labelidx = data.GetLabelIdx(expi,fly);
        labelidx = labelidx.vals;
      end
      noff = numel(labelidx) - numel(x);
      noff_start = floor(noff/2);
      noff_end = noff - noff_start;
      labelidx = labelidx(1+noff_start:end-noff_end);

      % histogram
      for i = 0:nbehaviors,
        idx = ~isnan(x) & labelidx == i;
        if ~any(idx), continue; end
        n(i+1) = n(i+1) + nnz(idx);
        counts_curr = histc(x(idx),edges);
        counts(i+1,:) = counts(i+1,:) + counts_curr(1:end-1);
        counts_outside(i+1,1) = counts_outside(i+1,1) + nnz(x(idx)<edges(1));
        counts_outside(i+1,2) = counts_outside(i+1,1) + nnz(x(idx)>edges(end));
      end
    end
  end

  frac = bsxfun(@rdivide,counts,n);
  frac_outside = bsxfun(@rdivide,counts_outside,n);

  centers_plot = [2*centers(1)-centers(2),centers,2*centers(end)-centers(end-1)];

end

hhist = zeros(1,nbehaviors+1);

if size(labelcolors,1) < nbehaviors,
  labelcolors(size(labelcolors,1)+1:nbehaviors,:) = jet(nbehaviors-size(labelcolors,1));
end

for i = 1:nbehaviors+1,
  hhist(i) = plot(hax,centers_plot,...
    [frac_outside(i,1),frac(i,:),frac_outside(i,2)],...
    '-','linewidth',2);
  if i == 1,
    set(hhist(i),'Color',unknowncolor);
    hold(hax,'on');
  else
    set(hhist(i),'Color',labelcolors(i-1,:));
  end
end

axisalmosttight(1/20,hax);
xtick = get(hax,'XTick');
xticklabel = get(hax,'XTickLabel');
xticklabel = cellstr(xticklabel);
i = find(xtick >= centers_plot(5),1);
if ~isempty(i),
  xtick = xtick(i:end);
  xticklabel = xticklabel(i:end);
end
i = find(xtick <= centers_plot(end-4),1,'last');
if ~isempty(i),
  xtick = xtick(1:i);
  xticklabel = xticklabel(1:i);
end
xtick = xtick(xtick>centers_plot(1));
xtick = xtick(xtick<centers_plot(end));
if centers_plot(1)==centers_plot(end)
  xtick = [centers_plot(1)];
else
  xtick = [centers_plot(1),xtick,centers_plot(end)];
end
xticklabel = [{sprintf('<%.1f',edges(1))};xticklabel;{sprintf('>=%.1f',edges(end))}];
set(hax,'XTick',xtick,'XTickLabel',xticklabel);

s = [{'Unknown'},data.labelnames];
hleg = legend(hax,hhist(n>0),s(n>0));
hxlabel = xlabel(hax,perframefn,'Interpreter','none');
hylabel = ylabel(hax,'Fraction of frames');

if ~holdstate,
  hold(hax,'off');
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Model/simulation parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sim=0;
OS_RATE = 8;
SNR = 1;
fc = 10e3/20e6; % sample rate is 20 MHz, top is 10 kHz offset
muFOC = floor(.01*2^12)/2^12;
muTOC = floor(.01*8*2^12)/2^12;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize LUTs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
make_srrc_lut;
make_train_lut;
make_trig_lut;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (sim)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Emulate microprocessor packet creation
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % data payload creation
    messageASCII = 'hello world!';
    message = double(unicode2native(messageASCII));
    % add on length of message to the front with four bytes
    msgLength = length(message);
    messageWithNumBytes =[ ...
        mod(msgLength,2^8) ...
        mod(floor(msgLength/2^8),2^8) ...
        mod(floor(msgLength/2^16),2^8) ...
        1 ... % message ID
        message];
    % add two bytes at the end, which is a CRC
    messageWithCRC = CreateAppend16BitCRC(messageWithNumBytes);
    ml = length(messageWithCRC);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % FPGA radio transmit core
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    data_in = 0;
    empty_in = 1;
    tx_en_in = 0;
    numBytesFromFifo = 0;
    num_samp = ml*8*2*2*3;
    x = zeros(1,num_samp);
    CORE_LATENCY = 4;
    data_buf = zeros(1,CORE_LATENCY);
    empty_buf = ones(1,CORE_LATENCY);
    clear_buf = zeros(1,CORE_LATENCY);
    tx_en_buf = zeros(1,CORE_LATENCY);
    for i1 = 1:num_samp
        % first thing the processor does is clear the internal tx fifo
        if i1 == 1
            clear_fifo_in = 1;
        else
            clear_fifo_in = 0;
        end

        if i1 == 5 % wait a little bit then begin to load the fifo
            empty_in = 0;
            numBytesFromFifo = 0;
        end

        data_buf = [data_buf(2:end) data_in];
        empty_buf = [empty_buf(2:end) empty_in];
        clear_buf = [clear_buf(2:end) clear_fifo_in];
        tx_en_buf = [tx_en_buf(2:end) tx_en_in];
        [i_out, q_out, re_byte_out, tx_done_out, d1, d2, d3] = ...
            qpsk_tx(data_buf(1),empty_buf(1),clear_buf(1),tx_en_buf(1));
        x_out = complex(i_out,q_out)/2^11;
        x(i1) = x_out;

        %%% Emulate read FIFO AXI interface
        if re_byte_out == 1 && numBytesFromFifo < length(messageWithCRC)
            data_in = messageWithCRC(numBytesFromFifo+1);
            numBytesFromFifo = numBytesFromFifo + 1;
        end
        % processor loaded all bytes into FIFO so begin transmitting
        if numBytesFromFifo == length(messageWithCRC)
            empty_in = 1;
            tx_en_in = 1;
        end
    end

    index = find(abs(x) > sum(SRRC))+24*8; % constant is pad bits
    offset = index(1)+6+length(TB_i)*OS_RATE;
    idx = offset:8:(offset+8*ml*4-1);
    y = x(idx); % four symbos per byte of data
    sc = zeros(1,18*8);
    sc(1:2:end) = real(y);
    sc(2:2:end) = imag(y);
    sh = sign(sc);
    sb = (sh+1)/2;
    d = zeros(1,ml);
    for i1 = 1:ml
        si = sb(1+(i1-1)*8:i1*8);
        d(i1) = bi2de(si);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Emulate channel
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % pad on either side with zeros
    p = complex(zeros(1,100),zeros(1,100));
    xp = [p x p]; % pad

    % Apply frequency offset and receive/over-the-air AWGN
    y = xp.*exp(1i*2*pi*fc*(0:length(xp)-1));
    rC = y/max(abs(y))*.1*2^11; % this controls receive gain
    r = awgn(rC,SNR,0,1);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load Chipscope samples
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (~sim)
    fid = fopen('rx.prn');
    M = textscan(fid,'%d %d %d %d','Headerlines',1);
    fclose(fid);
    is = double(M{3});
    qs = double(M{4});
    r = complex(is,qs);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main receiver core
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r_out = zeros(1,length(r));
s_f = zeros(1,length(r));
s_t = zeros(1,length(r));
s_c = zeros(1,length(r));
f_est = zeros(1,length(r));
t_est = zeros(1,length(r));
for i1 = 1:length(r)+200
    if i1 > length(r)
        r_in = 0;
    else
        r_in = r(i1);
    end
    i_in = round(real(r_in));
    q_in = round(imag(r_in));
    [r_out(i1), s_f(i1), s_c(i1), s_t(i1), t_est(i1), f_est(i1)] = ...
        qpsk_rx(i_in, q_in, floor(muFOC*2^12), floor(muTOC*2^12));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
subplot(2,3,1)
scatter(real(r),imag(r))
title('Pre FOC Signal');
subplot(2,3,4)
plot(real(r_out));
title('Pre FOC Signal (real part)');
subplot(2,3,2)
scatter(real(s_c),imag(s_c))
title('Post FOC Signal');
subplot(2,3,5)
plot(real(s_c));
title('Post FOC Signal(real part)');
subplot(2,3,3)
scatter(real(s_t),imag(s_t))
title('Post TOC Signal');
subplot(2,3,6)
plot(real(s_t));
title('Post TOC Signal(real part)');

figure(3)
plot(t_est);
title('Timing Estimate');


text:
% Fig. 5.18   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%
clf
 np=1;
dp=[1 .2 6.6*6.6+.1*.1 0 0];
nc=[1 1];
dc=[1 12];
nol=conv(np,nc);
dol=conv(dp,dc);
rlocus(nol,dol)
title(' Fig. 5.18 Root locus for noncollocated system')
axis([-15, 5, -7.5, 7.5])
z=0:.1:.9;
 wn=2:2:14;
 sgrid(z, wn)
 hold off


text:
function [pos vel alt] = RadarUKF(z, dt)
%
%
persistent Q R
persistent x P
persistent n m
persistent firstRun

if isempty(firstRun)
  Q = [ 0.01  0      0;
        0  0.01  0;
        0  0      0.01 ];

  R = 100;

  x = [0 90 1100]';
  P = 100*eye(3);

  n = 3;
  m = 1;

  firstRun = 1;
end

[Xi W] = SigmaPoints(x, P, 0);

fXi = zeros(n, 2*n+1);
for k = 1:2*n+1
  fXi(:, k) = fx(Xi(:,k), dt);
end

[xp Pp] = UT(fXi, W, Q);
%norm(xp - fx(x, dt))

hXi = zeros(m, 2*n+1);
for k = 1:2*n+1
  hXi(:, k) = hx(fXi(:,k));
end

[zp Pz] = UT(hXi, W, R);
%norm(zp - hx(xp))

Pxz = zeros(n, m);
for k = 1:2*n+1
  Pxz = Pxz + W(k)*(fXi(:, k) - xp)*(hXi(:, k) - zp)';
end

K = Pxz*inv(Pz);

x = xp + K*(z - zp);
P = Pp - K*Pz*K';

pos = x(1);
vel = x(2);
alt = x(3);

%------------------------------
function xp = fx(x, dt)
%
%
A = eye(3) + dt*[ 0 1 0;
                  0 0 0;
                  0 0 0 ];

xp = A*x;

%------------------------------
function yp = hx(x)
%
%
yp = sqrt(x(1)^2 + x(3)^2);

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
function prob = survProbGenHazFun(t,b,hazardFun)
% survProbGeneralHazardFunction: Computes survival probability using the
% hazard function 'hazardFun'

h = @(t)hazardFun(t,b);

H = zeros(size(t));
H(1) = quad(h,0,t(1));
for i=2:length(H)
   H(i) = H(i-1) + quad(h,t(i-1),t(i));
end

prob = exp(-H);


text:
function [] = test_lasso()

    %clc;
    clear;
    close all;

    rng('default')

    %% Set algorithms
    %algorithms = {'PG-BKT', 'PG-TFOCS-BKT', 'PG-FIX', 'SUBG-DEC', 'SUBG-BKT', 'APG-TFOCS-BKT', 'FISTA', 'CD-LASSO', 'ADMM-LASSO'};
    algorithms = {'SMOOTH-FIX', 'SMOOTH-BKT', 'PG-WOLFE', 'PG-TFOCS-BKT', 'PG-FIX', 'SUBG-DEC', 'SUBG-BKT', 'APG-BKT', 'APG-TFOCS-BKT', 'FISTA', 'ADMM-LASSO'};
    %algorithms = {'SUBG-DEC', 'SUBG-BKT', 'SUBG-TFOCS-BKT'};
    %algorithms = {'SUBG-DEC', 'SUBG-BKT', 'SUBG-TFOCS-BKT', 'SMOOTH-FIX', 'FISTA', 'ADMM-LASSO'};
    algorithms = {'SUBG-DEC', 'SUBG-BKT'};

    %% prepare dataset
    if 1
        n = 1280;
        d = 100;
        k = 15;                                     % cardinality of nonzero elements
        [A,~] = qr(randn(n,d),0);
        A = A';
        p = randperm(n);
        p = p(1:k);                                 % select location of k nonzeros
        x0 = zeros(n,1);
        x0(p) = randn(k,1);
        b = A*x0 + .02*randn(d, 1);                 % add random noise
        lambda_max = norm( A'*b, 'inf' );
        lambda = 0.1*lambda_max;
    else
        n = 500;
        d = 100;
        A = randn(d,n);
        b = randn(d,1);
        lambda = 5;
    end

    %% initialize
    w_init = rand(n,1);
    w_list = cell(length(algorithms),1);
    info_list = cell(length(algorithms),1);

    %% perform algorithms
    for alg_idx=1:length(algorithms)

        %% define problem definitions
        clear problem;
        problem = lasso(A, b, lambda, 'prox_reg');

        fprintf('\n\n### [%02d] %s ###\n\n', alg_idx, algorithms{alg_idx});

        clear options;
        % general options for optimization algorithms
        options.w_init = w_init;
        options.tol_gnorm = 1e-10;
        options.max_epoch = 100;
        options.max_iter = options.max_epoch;
        options.verbose = true;
        options.f_opt = 0;
        options.store_w = false;

        switch algorithms{alg_idx}
            case {'SMOOTH-FIX'}
                smooth_mu = 0.1;
                clear problem;
                problem = lasso(A, b, lambda, 'smooth', smooth_mu);

                options.step_alg = 'fix';
                options.step_init = 1/(problem.L+lambda/smooth_mu);
                [w_list{alg_idx}, info_list{alg_idx}] = smoothing_gd(problem, options);

            case {'SMOOTH-BKT'}
                smooth_mu = 0.1;
                clear problem;
                problem = lasso(A, b, lambda, 'smooth', smooth_mu);

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = smoothing_gd(problem, options);

            case {'SUBG-DEC'}

                options.step_alg = 'decay-7';
                options.step_init = 1/problem.L;
                [w_list{alg_idx}, info_list{alg_idx}] = subg(problem, options);

            case {'SUBG-BKT'}

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = subg(problem, options);

            case {'SUBG-TFOCS-BKT'}

                options.step_alg = 'tfocs_backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = subg(problem, options);

            case {'PG-BKT'}

                options.step_alg = 'backtracking';
                %options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'PG-FIX'}

                options.step_alg = 'fix';
                options.step_alg = 'fix';
                options.step_init = 1/problem.L;
                %options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'PG-WOLFE'}

                options.step_alg = 'strong_wolfe';
                options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'PG-TFOCS-BKT'}

                options.step_alg = 'tfocs_backtracking';
                options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'APG-BKT'}

                options.step_alg = 'backtracking';
                options.step_init_alg = 'bb_init';
                %[w_list{alg_idx}, info_list{alg_idx}] = sd_nesterov(problem, options);
                [w_list{alg_idx}, info_list{alg_idx}] = ag(problem, options);

            case {'APG-TFOCS-BKT'}

                options.step_alg = 'tfocs_backtracking';
                options.step_init_alg = 'bb_init';
                %[w_list{alg_idx}, info_list{alg_idx}] = sd_nesterov(problem, options);
                [w_list{alg_idx}, info_list{alg_idx}] = ag(problem, options);

            case {'FISTA'}

                options.sub_mode  = 'FISTA';
                [w_list{alg_idx}, info_list{alg_idx}] = ista(problem, options);

            case {'ADMM-LASSO'}

                options.rho = 0.1;
                [w_list{alg_idx}, info_list{alg_idx}] = admm_lasso(problem, options);

            case {'CD-LASSO'}

                options.sub_mode = 'lasso';
                [w_list{alg_idx}, info_list{alg_idx}] = cd_lasso_elasticnet(problem, options);

            case {'L-BFGS-BKT'}

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'L-BFGS-WOLFE'}

                options.step_alg = 'strong_wolfe';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'P-Newton-CHOLESKY'}

                options.sub_mode = 'CHOLESKY';
                options.step_init_alg = 'bb_init';
                options.step_alg = 'tfocs_backtracking';
                %options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            otherwise
                warn_str = [algorithms{alg_idx}, ' is not supported.'];
                warning(warn_str);
                w_list{alg_idx} = '';
                info_list{alg_idx} = '';
        end

    end

    fprintf('\n\n');

    %% plot all

    % display iter vs cost
    display_graph('iter','cost', algorithms, w_list, info_list);
    % display time vs cost
    display_graph('time','cost', algorithms, w_list, info_list);
    % display iter vs. l1 norm, i.e. the toral number of non-zero elements
    display_graph('iter','reg', algorithms, w_list, info_list);

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function copyScan(sessionDir,dataType,scan)
%
% function copyScan(sessionDir,dataType,scan)
%
% Copy a scan from a different session.
% Creates dataType 'Copied' if it doesn't already exist.
%
% sessionDir: path string to the other session
% dataType, scan: specifies where to find the tSeries within that session
%
% Output is a new set of tSeries files in a new Scan subdirectory
% under the Copied dataType directory.
%
% djh, 2/28/2001

mrGlobals

% Make sure dataType 'Copied' exists
if ~existDataType('Copied')
    addDataType('Copied');
end

% Open hidden gray
hiddenView = initHiddenGray;
hiddenView = selectDataType(hiddenView,existDataType('Copied'));

% update dataTYPES.scanParams with the new scan,
% copying scanParams, blockedAnalysisParams, & eventAnalysisParams
% from the other session.
newScanNum = numScans(hiddenView)+1;
ndataType = hiddenView.curDataType;
% load dataTYPES from other session
s = load(fullfile(sessionDir,'mrSESSION.mat'));
% get scanParams, blockedAnalysisParams, & eventAnalysisParams from other session
ndataTypeOtherSession = existDataType(dataType,s.dataTYPES);
scanParams = s.dataTYPES(ndataTypeOtherSession).scanParams(scan);
% blast slices & cropSize because inplane view of new scan is bogus
scanParams.slices = [];
scanParams.cropSize = [];
blockedAnalysisParams = s.dataTYPES(ndataTypeOtherSession).blockedAnalysisParams(scan);
eventAnalysisParams = s.dataTYPES(ndataTypeOtherSession).eventAnalysisParams(scan);
dataTYPES(ndataType).scanParams(newScanNum) = scanParams;
dataTYPES(ndataType).blockedAnalysisParams(newScanNum) = blockedAnalysisParams;
dataTYPES(ndataType).eventAnalysisParams(newScanNum) = eventAnalysisParams;
dataTYPES(ndataType).scanParams(newScanNum).annotation = ['Copied from ',sessionDir,', dataType: ',dataType,', scan: ',num2str(scan)];

% Load gray coords from other session.
% Compute intersection of their coords with my gray coords.
grayCoordsFile = fullfile(sessionDir,'Gray','coords.mat');
s = load(grayCoordsFile);
[c,ia,ib] = intersectCols(hiddenView.coords,s.coords);
% Load the tSeries from the other session
tSeriesPath = fullfile(sessionDir,'Gray',dataType,'TSeries',['Scan',num2str(scan)],'tSeries1.mat');
s = load(tSeriesPath);

tSeries = NaN*ones([size(s.tSeries,1) size(hiddenView.coords,2)]);
tSeries(:,ia) = s.tSeries(:,ib);
%This should be able to be used the same as before as the gray view has not
%yet been updated
savetSeries(tSeries,hiddenView,newScanNum,1);
saveSession

% Loop through the open views, switch their curDataType appropriately,
% and update the dataType popups
INPLANE = resetDataTypes(INPLANE,ndataType);
VOLUME = resetDataTypes(VOLUME,ndataType);
FLAT = resetDataTypes(FLAT,ndataType);

return;

function viewList=resetDataTypes(viewList,ndataType)
% Loops through the views, updating the dataType to reflect the additional scan.
% The call to selectDataType updates the dataType popup.
for s=1:length(viewList)
    if ~isempty(viewList{s})
        if viewList{s}.curDataType == ndataType
            viewList{s} = selectDataType(viewList{s},1);
            viewList{s} = selectDataType(viewList{s},ndataType);
        end
    end
end
return

% Debug/test

clear all; close all;
sessionDirs = {'e:/mri/test','j:/djhRetinotopy/030701a'};
newSession = 'e:/mri/foo';
mergeSessions(sessionDirs,newSession);
chdir(newSession)
mrLoadRet
copyScan(sessionDirs{1},'Original',1);
copyScan(sessionDirs{1},'Original',2);
copyScan(sessionDirs{2},'Averages',1);
copyScan(sessionDirs{2},'Averages',2);


text:
function a = mirt_dctn(a)
%MIRT_DCTN N-D (multidimensional) discrete cosine transform.
%   B = MIRT_DCTN(A) returns the discrete cosine transform of A.
%   The array B is the same size as A and contains the
%   discrete cosine transform coefficients.
%
%   This function works much faster than Matlab standard
%   dct (for 1D case) or dct2 (for 2D case), and also allows the ND input.
%   The function takes advantage of 1) FFT 2) fast permutation through indicies
%   3) persistent precomputation (coefficients and indicies are precomputed during
%   the first run).
%
%   This transform can be inverted using MIRT_IDCTN.
%   Author: Andriy Myronenko, see www.bme.ogi.edu/~myron
%   revised on May 06, 2010
%
%   See also MIRT_IDCT.

persistent siz ww ind isreala;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Check input
if (nargin == 0) || isempty(a) ,
    error('Insufficient input');
end

isreala=isreal(a);
ndim=ndims(a);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Check for the row vector
transpose=0;
if (ndim==2) && (size(a,1)==1)
    transpose=1; a=a';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Check if the variable size has changed and we need to
%%% precompute weights and indicies

precompute=0;
if  ~exist('siz','var'),
    precompute=1;
elseif abs(numel(siz)-ndims(a))>0
    precompute=1;
elseif sum(abs(siz-size(a)),2)>0,
    precompute=1;
elseif isreala~=isreal(a),
    precompute=1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Precompute weights and indicies
if precompute,
    siz=size(a);
    ndim=ndims(a);

    for i=1:ndim,
        n=siz(i);

        ww{i} = 2*exp(((-1i*pi)/(2*n))*(0:n-1)')/sqrt(2*n);
        ww{i}(1) = ww{i}(1) / sqrt(2);
        ind{i}=bsxfun(@plus, [(1:2:n) fliplr(2:2:n)]', 0:n:n*(prod(siz)/n-1));
        if (siz(i)==1), break; end;
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Actual multidimensional DCT. Handle 1D and 2D cases
%%% separately, because .' is much faster than shiftdim.

% check for 1D or 2D cases
if ndim==2
    if (min(siz)==1),
        a=dct(a,ww{1},ind{1});if transpose, a=a'; end;  % 1D case
    else a = dct(dct(a,ww{1},ind{1}).',ww{2},ind{2}).'; % 2D case
    end
else
    % ND case (3D and higher)
    for i=1:ndim
        a=reshape(dct(reshape(a,siz(1),[]),ww{i},ind{i}), siz); % run dct along vectors (1st dimension)
        siz=[siz(2:end) siz(1)];                   % circular shift of size array by 1
        a=shiftdim(a,1);                           % circular shift dimensions by 1
    end
end

function a=dct(a,ww,ind)
%DCT  Discrete cosine transform 1D (operates along first dimension)

isreala=isreal(a);
k=1; if ~isreala, ia = imag(a); a = real(a); k=2; end; % check if complex

% k=1 (if 'a' is real) and 2 (if 'a' is complex)
for i=1:k
     a=a(ind);     % rearange
     a = fft(a);  % ifft
     a = real(bsxfun(@times,  ww, a)); % multiply weights

    % check if the data is not real
    if ~isreala,
        if i==1, ra = a; a = ia; clear ia;  % proceed to imaginary part
        else a = complex(ra,a); end;        % finalize the output
    end

end