text:
close all
clear
clc
graphics_toolkit qt

h.user = 'username';
h.ii = 1;
h.measurements = {};
h.measurements(1, :) = {'Measurement', 'Frequency', 'Amplitude', 'Phase'};
h.freqs = [177, 297, 500, 841, 1414, 2378, 4000, 6727, 11314];

h.gui = figure('Position', [10 10 2000 1200],
       'NumberTitle', 'off',
       'Name', 'Comb-filter phase measurement GUI',
       'toolbar', 'none',
       'menubar', 'none');

% Selection of the pure tone frequency and adaption of the figure axes
function set_freq(obj)
  h = guidata(obj);
  string = get(obj, 'string');
  h.f = str2num(string);
  set(h.dot,'visible','off')
  h.dot = plot(0, 0, 'o', 'color', 'r', 'markersize', 1);
  guidata(obj, h);
end

% Callback for click in the figure
function clickcallback(obj)
  h = guidata(obj);
  axesHandle  = get(obj, 'Parent');
  coordinates = get(gca, 'CurrentPoint');
  coordinates = coordinates(1, 1:2);
  [x_max, y_max] = size(get(obj, 'cdata'));
  y = coordinates(2);
  coordinates(2) = abs((coordinates(2) - 128));
  h.x = coordinates(1);
  h.y = coordinates(2)./128;
  guidata(obj, h);
  coord2param(obj);
  h = guidata(obj);
  h.fs = 48000;
  delay = 1+h.phi/(2*pi).*h.fs./h.f;
  b = zeros(1, 100);
  b(round(delay)) = h.a;
  system(['echo mha.transducers.mhachain.injector.B = [' sprintf('%.15f ',b) '] | nc -w 1 127.0.0.1 33337']);
  hold all
  delete(h.dot);
  h.dot = plot(h.x, y, 'o', 'markerfacecolor', 'r', 'markeredgecolor', 'r', 'markersize', 8);
  drawnow;
  guidata(obj, h);
end

% Convert input values to amplitude (a) and phase (phi)
function coord2param(obj)
  h = guidata(obj);
  h.a = h.y;
  h.phi = (5/2*pi)/160.*h.x;
  guidata(obj, h);
  amplification_rule(obj);
  h = guidata(obj);
  guidata(obj, h);
end

% Scale amplitude
function amplification_rule(obj)
  h = guidata(obj);
  switch h.f
    case 177
      h.a = 10.^((h.a*30-30)/20); % 0 dB to -30 dB
    case 297
      h.a = 10.^((h.a*30-30)/20); % 0 dB to -30 dB
    case 500
      h.a = 10.^((h.a*30-30)/20); % 0 dB to -30 dB
    case 841
      h.a = 10.^((h.a*30-35)/20); % -5 dB to -35 dB
    case 1414
      h.a = 10.^((h.a*35-45)/20); % -10 dB to -40 dB
    case 2378
      h.a = 10.^((h.a*35-50)/20); % -15 dB to -45 dB
    case 4000
      h.a = 10.^((h.a*30-30)/20); %  -5 dB to -35 dB  (0 dB to -30 dB)
    case 6727
      h.a = 10.^((h.a*30-35)/20); %  -5 dB to -35 dB
    case 11314
      h.a = 10.^((h.a*30-35)/20); %  -5 dB to -35 dB
  endswitch
  guidata(obj, h);
end

% Saving of the data
function finish(obj)
  h = guidata(obj);
  h.measurements(end+1, :) = {h.ii, h.f, h.a, h.phi};
  progress = [num2str(h.ii), ' out of 9 complete'];
  set(h.measurement_complete, 'string', progress);
  if h.ii == 9 % # of measurements
    set(h.measurement_complete, 'string', 'Measurement complete!');
    filename = num2str([h.user, '.csv']);
    cell2csv(filename, h.measurements)
  end
  h.ii = h.ii + 1;
  guidata(obj, h);
end

% Toggle between amplification and no amplification
function play_reference(obj)
  h = guidata(obj);
  b = zeros(1, 100);
  string = get(obj, 'string');
  switch string
    case 'Amp off'
      system(['echo mha.transducers.mhachain.injector.B = [' sprintf('%.15f ',b) '] | nc -w 1 127.0.0.1 33337']);
    case 'Amp on'
      delay = 1+h.phi/(2*pi).*h.fs./h.f;
      b(round(delay)) = h.a;
      system(['echo mha.transducers.mhachain.injector.B = [' sprintf('%.15f ',b) '] | nc -w 1 127.0.0.1 33337']);
  endswitch
  guidata(obj, h);
end

% Keyboard event detection
function keyPress(obj, e)
  h = guidata(obj);
  switch e.Key
    case '1'
      play_reference(h.reference, [], []);
    case '2'
      play_reference(h.adjusted, [], []);
  endswitch
  guidata(obj);
end

h.finish = uicontrol('style', 'pushbutton',
  'units', 'normalized',
  'string', 'Save',
  'fontsize', 18,
  'fontweight', 'bold',
  'callback', @finish,
  'position', [0.75 0.1 0.2 0.1]);

h.amplification_buttons = uibuttongroup('position', [0.75 0.35 0.2 0.1]);
h.amp_off = uicontrol(h.amplification_buttons, 'style', 'radiobutton',
  'units', 'normalized',
  'string', 'Amp off',
  'fontsize', 10,
  'callback', @play_reference,
  'position', [0.1 0.4 0.3 0.2]);
h.amp_on = uicontrol(h.amplification_buttons, 'style', 'radiobutton',
  'units', 'normalized',
  'string', 'Amp on',
  'fontsize', 10,
  'callback', @play_reference,
  'position', [0.6 0.4 0.3 0.2]);

h.freq_buttons = uibuttongroup('position', [0.75 0.475 0.09 0.3]);
h.freq_177 = uicontrol(h.freq_buttons, 'style', 'radiobutton',
  'units', 'normalized',
  'selected', 'on',
  'string', '177',
  'fontsize',10,
  'fontweight','bold',
  'callback', @set_freq,
  'position', [0.05 0.85 0.9 0.05]);
h.freq_297 = uicontrol(h.freq_buttons, 'style', 'radiobutton',
  'units', 'normalized',
  'string', '297',
  'fontsize',10,
  'fontweight','bold',
  'callback', @set_freq,
  'position', [0.05 0.75 0.9 0.05]);
h.freq_500 = uicontrol(h.freq_buttons, 'style', 'radiobutton',
  'units', 'normalized',
  'string', '500',
  'fontsize',10,
  'fontweight','bold',
  'callback', @set_freq,
  'position', [0.05 0.65 0.9 0.05]);
h.freq_841 = uicontrol(h.freq_buttons, 'style', 'radiobutton',
  'units', 'normalized',
  'string', '841',
  'fontsize', 10,
  'fontweight','bold',
  'callback', @set_freq,
  'position', [0.05 0.55 0.9 0.05]);
h.freq_1414 = uicontrol(h.freq_buttons, 'style', 'radiobutton',
  'units', 'normalized',
  'string', '1414',
  'fontsize',10,
  'fontweight','bold',
  'callback', @set_freq,
  'position', [0.05 0.45 0.9 0.05]);
h.freq_2378 = uicontrol(h.freq_buttons, 'style', 'radiobutton',
  'units', 'normalized',
  'string', '2378',
  'fontsize',10,
  'fontweight','bold',
  'callback', @set_freq,
  'position', [0.05 0.35 0.9 0.05]);
h.freq_4000 = uicontrol(h.freq_buttons, 'style', 'radiobutton',
  'units', 'normalized',
  'string', '4000',
  'fontsize',10,
  'fontweight','bold',
  'callback', @set_freq,
  'position', [0.05 0.25 0.9 0.05]);
h.freq_6727 = uicontrol(h.freq_buttons, 'style', 'radiobutton',
  'units', 'normalized',
  'string', '6727',
  'fontsize',10,
  'fontweight','bold',
  'callback', @set_freq,
  'position', [0.05 0.15 0.9 0.05]);
h.freq_11314 = uicontrol(h.freq_buttons, 'style', 'radiobutton',
  'units', 'normalized',
  'string', '11314',
  'fontsize',10,
  'fontweight','bold',
  'callback', @set_freq,
  'position', [0.05 0.05 0.9 0.05]);

h.measurement_complete = uicontrol('style', 'text',
  'units', 'normalized',
  'string', '',
  'fontsize', 11,
  'fontweight', 'bold',
  'position', [0.75 0.225 0.2 0.1]);

img = imread('grid_extended.png');
min_x = 0;
max_x = 160;
min_y = 0;
max_y = 128;

fig_pos_click = {[0.05 0.1 0.65 0.8]};
h.ax_click = axes('position', fig_pos_click{});
box('on')
h_image = imagesc([min_x max_x], [min_y, max_y], img);

set(gca, 'xtick', [0 32 64 96 128 160]);
set(gca, 'xticklabel', {'0', '^{\pi}/_{2}', '\pi', '^{3\pi}/_{2}', '2\pi', '^{5\pi}/_{2}'});
set(gca, 'ytick', [0 64/3 128/3 64 256/3 320/3 128]);
set(gca, 'yticklabel', {'+15', '+10', '+5', '0', '-5', '-10', '-15'});
set(gca, 'linewidth', 3, 'fontsize', 20);
hold all

h.dot = plot(0, 0, 'o', 'color', 'r', 'markersize', 1);

set(h_image,'ButtonDownFcn',@clickcallback);
set(gcf, 'KeyPressFcn', @keyPress);
guidata(gcf, h);

author: m-r-s

repo: hearingaid-prototype


text:
function img = discreteReuleauxRevol(varargin)
%DISCRETEREULEAUXREVOL Discretize the revolution of a Reuleaux triangle
%
%   A Reuleaux triangle is formed by 3 circle arcs, and has the property to
%   have a constant width whatever the  orientation.
%   In 3 dimension, the surface of revolution obtained by rotating a
%   Reuleaux triangle around one of its axis of symetry has also a constant
%   width. The resulting shape looks like a hazelnut.
%
%   IMG = discreteReuleauxRevol(LX, LY, LZ, REULEAUX)
%   REULEAUX is the domaine enclosed by the revolution of a Reuleaux
%   triangle around one of its axis. REULEAUX is defined by:
%   [XC YC ZC R THETA PHI]
%   where (XC, YC, ZC) is the geometric centroid the triangle, R is the
%   radius of the circle arcs, and THETA and PHI define the orientation of
%   the Reuleaux Apex: THETA is the colatitude, between 0 and 180 degrees,
%   and PHI is the azimut, between 0 and 360 degrees.
%
%   IMG = discreteReuleauxRevol(LX, LY, LZ, REULEAUX, D)
%   Use the dilation of a Reuleaux shape by a ball of radius D. The
%   constant diameter of the resulting shape equals R + 2*D.
%
%   IMG = discreteReuleauxRevol(DIM, ...)
%   DIM is the size of image, with the format [x0 dx x1;y0 dy y1;z0 dz z1]
%
%   Examples
%   % Display oriented reuleaux triangle revolution
%     reuleaux = [50 50 50  60  30 40];
%     img = discreteReuleauxRevol(1:100, 1:100, 1:100, reuleaux);
%     figure();
%     isosurface(img, .5);
%     axis equal;
%
%   % The same, with a little bit dilation
%     reuleaux = [50 50 50  60  30 40]
%     img2 = discreteReuleauxRevol(1:100, 1:100, 1:100, reuleaux, 10);
%     figure();
%     isosurface(img2, .5);
%     axis equal;
%
%   % Alternative syntax
%     img = discreteReuleauxRevol([1 1 100;1 1 100;1 1 100], reuleaux);
%
%   See Also
%   imShapes, discreteEllipsoid, discreteTorus, discreteBall
%

%
% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2006-02-27
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

%   HISTORY
%   04/01/2007: concatenate transforms before transforming points
%   04/03/2009: use meshgrid
%   30/04/2009: update transforms
%   29/05/2009: use more possibilities for specifying grid

% compute coordinate of image voxels
[lx, ly, lz, varargin] = parseGridArgs3d(varargin{:});
[x, y, z] = meshgrid(lx, ly, lz);

% parameters of Revolved Reuleaux triangle
reul    = varargin{1};
center  = reul(1:3);
R       = reul(4);
theta   = deg2rad(reul(5));
phi     = deg2rad(reul(6));
varargin(1) = [];

% extract dilation factor if present
dil = 0;
if ~isempty(varargin)
    dil = varargin{1};
end

% top point of triangle
x1  = 0;
y1  = 0;
z1  = 0;

% the height of the triangle
h   = R * sqrt(3) / 2;
dil2 = dil * sqrt(3) / 2;

% create empty image
img = false(size(x));

% transforms voxels according to shape orientation
trans = composeTransforms3d(...
    createTranslation3d(-center),...
    createRotationOz(-phi),...
    createRotationOy(-theta), ...
    createTranslation3d([0 0 h*2/3]));
[x, y, z] = transformPoint3d(x, y, z, trans);

% distance of the points to the revolution axis
rho = hypot(x-x1, y-y1);

% points above the top of the triangle
ind = find((z <= z1) & (z >= z1-dil));
dist = rho(ind).^2 + (z(ind)-z1).^2;
img(ind) = dist <= dil ^ 2;

% points above the basis of the triangle and below the top
ind  = find(z >= z1-dil2 & z <= (z1+h));
dist = (rho(ind)+R/2).^2 + (z(ind)-z1-h).^2;
img(ind) = dist < (R+dil)^2;

% points around the 2 lower corners
dist = (rho-R/2).^2 + (z-z1-h).^2;
img(dist < dil^2) = 1;
ind  = find(z > z1+h & z <= z1+h+dil2);
dist = rho(ind) ;
img(ind(dist < R/2)) = 1;

% points below the basis of the triangle
ind  = find(z >= z1+h+dil2);
dist = rho(ind).^2 + (z(ind)-z1).^2;
img(ind) = dist < (R+dil)^2;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
% readneuronedata() -  Read NeurOne Binary files
%
% Usage:  >> data = readneuronedata(dataFiles,nChannels,chans)
%
% ===================================================================
% Inputs:
%         dataFiles      - A cell structure containing full paths for data
%                          files to be read. This is input argument is
%                          required.
%         nChannels      - Another required input argument:
%                          the total number of channels.
%         chans          - A vector containing the channel numbers to be
%                          read. They have to be arranged in an ascending
%                          order. This input argument is optional.
% ====================================================================
% Output:
%         data           - A matrix where all numerical data is stored.
%                          The data is arranged so that the data for each
%                          channel is stored in a row according
%                          to its number.
%
% ========================================================================
% NOTE:
% This file is part of the NeurOne data import plugin for EEGLAB.
% ========================================================================
%
% Current version: 1.0.3.4 (2016-06-17)
% Author: Mega Electronics

function data = readneuronedata(dataFiles,nChannels,chans)
%% Parse input arguments
p = inputParser;
p.addRequired('dataFiles', @iscellstr);
p.addRequired('nChannels', @isnumeric);
p.addOptional('chans', 0, @isnumeric);

p.parse(dataFiles, nChannels, chans);
arglist=p.Results;

%% Prepare reading data
nDataFiles = numel(arglist.dataFiles);
chans = arglist.chans;
nChannels = arglist.nChannels;

if chans==0
    chans=1:nChannels;
end

% Get all file sizes
fSize=zeros(1,nDataFiles);
tmp={};

for k=1:nDataFiles
    tmp{k,1}=dir(dataFiles{k,1});
    fSize(1,k)=tmp{k,1}.bytes;
end

fSizes=cumsum(fSize);

% Total size of the data
totalSize=sum(fSize);
% Total number of data points (per channel)
dataPntsTotal=(totalSize/4)/nChannels;
% Number of data points in each binary file (per channel)
dataPnts=(fSize./4)./nChannels;

%% Read binary data

fprintf('Allocating memory...\n')
data=zeros(numel(chans),dataPntsTotal);

% Option 1: Read all channels at once (faster method)
if numel(chans)==nChannels
    fprintf('Reading data...');
    data=data(:);
    readidx=[0 (fSizes/4)];
    for n=1:nDataFiles
        fid = fopen([dataFiles{n}], 'rb');
        data((1+readidx(n)):readidx(n+1),1)=fread(fid, ...
            fSize(n)/4,'int32');
        fclose(fid);
    end
    data=reshape(data,nChannels,dataPntsTotal);
    fprintf('%s','Done');
else
    % Option 2: Read only specific channels to save memory
    fprintf('Reading data... 0%%');
    readidx=[0 (fSizes/4)./nChannels];

    % Read channels one by one
    for k=1:numel(chans)
        for n=1:nDataFiles

            fid = fopen([dataFiles{n}], 'rb');
            fseek(fid, 4*(chans(k) - 1), 'bof');
            data(k,(1+readidx(n)):readidx(n+1)) = fread(fid, ...
                dataPnts(n),'int32',4*(nChannels-1))';
            fclose(fid);
        end

    % Print progress to the command window
    ii = floor(k/numel(chans)*100);
        if ii<10
            ii = [num2str(ii) '%'];
            fprintf(1,'\b\b%s',ii);
        elseif ii==100;
            ii = 'Done';
            fprintf(1,'\b\b\b%s',ii);
        else
            ii = [num2str(ii) '%'];
            fprintf(1,'\b\b\b%s',ii);

        end
    end
end

fprintf('\n');

%% Convert data from nanovolts to microvolts
data=data./1000;

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [dim lf] = readBESAlft(filename, chanid)

% readBESAlft reads the leadfield from the BESA leadfield file format.
% The leadfield matrix and a vector with the leadfield dimensions are
% returned.
%
% Parameters:
%     [filename]
%         In the case that the current folder is not the folder containing
%         the file it should be the full path including the name of the
%         lft file else only the name of the file should be specified.
%
%     [chanid]
%         Additional parameter specifying particular channel for which
%         leadfield information should be returned.
%
% Return:
%     [dim]
%         Leadfield dimensions {<number sensors>, <number source space
%         nodes>, <number source directions>}. Please note, in case that
%         the additional parameter is specified, the number of sensors will
%         be set to 1.
%
%     [lf]
%         Leadfield matrix with <number sensors> rows and <number source
%         space nodes> x <number source directions> columns. Each row
%         contains first the potentials for sources at all nodes in x-dir,
%         then for sources at all nodes in y-dir, and finally for all
%         nodes in z-dir.
%
% Copyright (C) 2013, BESA GmbH
%
% File name: readBESAlft.m
%
% Author: Benjamin Lanfer/Robert Spangler
% Created: 2013-11-27

% Check additional parameter
if exist('chanid', 'var')
    % check if chanid is scalar
    if isscalar(chanid)
        % get leadfield info from specific channels
    else
        % error; only scalar values allowed for chanid
        error('Parameter chanid is not a scalar value!');
    end;
else
    % get leadfield info all channels
    chanid = [];
end;

% Try to open file.
FileID = fopen(filename, 'rb');
if(FileID < 0)
	fprintf('Error opening file.\n');
	return
end

% Read version number (int32)
[VersionNumber NrReadElements] = fread(FileID, 1, 'int32');
if(NrReadElements ~= 1)
	fprintf('Could not read number of elements.\n');
	return
end

% Check version number
ExpectedVersionNumber = 1;
if(VersionNumber ~= ExpectedVersionNumber)
	fprintf('Wrong version number. Expected: %d, read %d.\n', ExpectedVersionNumber, ...
		VersionNumber);
	return
end

% Read number of sensors (int32)
[NumberSensors NrReadElements] = fread(FileID, 1, 'int32');
if(NrReadElements ~= 1)
	fprintf('Could not read number of sensors.\n');
	return
end

% In case the leadfield for one specific channel has to be returned, check
% if chanid is in correct range.
if ~isempty(chanid)
    if ((chanid > NumberSensors) ||...
        (chanid < 1))
        % error; chanid out of range
        error('Parameter chanid is out of range!');
    end;
end;

% Read number of source space nodes (int32)
[NumberSourceSpaceNodes NrReadElements] = fread(FileID, 1, 'int32');
if(NrReadElements ~= 1)
	fprintf('Could not read number of source space nodes.\n');
	return
end

% Read number of source directions per node (int32)
[NumberDirections NrReadElements] = fread(FileID, 1, 'int32');
if(NrReadElements ~= 1)
	fprintf('Could not read number of source directions.\n');
	return
end

% Read maximum LF values for each source (source nodes x directions) (float32)
NumberColumns = NumberDirections * NumberSourceSpaceNodes;
[MaxVals NrReadElements] = fread(FileID, NumberColumns, 'float32');
if(NrReadElements ~= NumberColumns)
	fprintf('Could not read maximum leadfield values from file.\n');
	return
end

if ~isempty(chanid)
    % Get leadfield for specific channel only
    lf = [];
    NumberSensors = 1;

    % Go to correct position in file where leadfield info for current
    % channel is stored
    NumSkippedBytes = (chanid-1)*NumberSourceSpaceNodes*NumberDirections*2;
    status = fseek(FileID, NumSkippedBytes, 0);
    if status ~= 0
        error('Unable to move to specified position in file.');
    end;

    % Read compactly stored LF values for channel #chanid (int16)
    NumberLFValues = NumberColumns * 1;
    [CompactLFVals NrReadElements] = fread(FileID, NumberLFValues, 'int16');
    if(NrReadElements ~= NumberLFValues)
        fprintf('Could not read leadfield values from file.\n');
        return
    end

else
    % Get leadfield for all channels

    % Read compactly stored LF values (int16)
    NumberLFValues = NumberColumns * NumberSensors;
    [CompactLFVals NrReadElements] = fread(FileID, NumberLFValues, 'int16');
    if(NrReadElements ~= NumberLFValues)
        fprintf('Could not read leadfield values from file.\n');
        return
    end
end

% Reshape matrix with compact LF values.
CompactLFVals = reshape(CompactLFVals, NumberColumns, NumberSensors);
CompactLFVals = CompactLFVals';

% Compute factors for converting compactly stored LF values
% to float values.
ConversionFactors = MaxVals / 32767.0;

% Clear superfluous parameters
clear('MaxVals');

% Convert compact LF values to float.
%lf = CompactLFVals * repmat(ConversionFactors, 1, NumberColumns);

% NOTE RS:
% * operator in line %lf = CompactLFVals * repmat(...) is a real
% matrix multiplication, not a scalar multiplication. Is this intended?

% NOTE RS:
% Use alternative conversion of compact LF values to float LF
% values due to the extrem memory consumption of the repmat function:
for iCol=1:NumberColumns
    CompactLFVals(:,iCol) = CompactLFVals(:,iCol)*ConversionFactors(iCol);
end
lf = CompactLFVals;

% Copy dimensions.
dim = [NumberSensors; NumberSourceSpaceNodes; NumberDirections];

% Close file.
fclose(FileID);
end


text:
filename='Gripping_quad_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'PROJECTED GRADIENT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.75;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function read_bv
% read_bv - read data from an eeg-File
%
% SYNOPSIS
%    data = read_bv(file, HDR, OPT);
%
% ARGUMENTS
%                file - Name of EEG file (.eeg) is appended)
%                HDR  - Information about the file (read from the *.vhdr header file)
%                   .fs      - Sampling rate
%                   .nChans  - Number of channels
%                   .nPoints - Number of data points in the file (optional)
%                   .scale   - Scaling factors for each channel
%                   .endian  - Byte ordering: 'l' little or 'b' big
%                OPT  - Struct with following fields
%                   .chanidx         - Indices of the channels that are to be read
%                   .fs              - Down sample to this sampling rate
%                   .filt_b          - Filter coefficients of IIR filter
%                                      applied to raw data (b part)
%                                      (optional)
%                   .filt_a          - Filter coefficients of IIR filter
%                                      applied to raw data (a part)
%                                      (optional)
%                   .filt_subsample  - Filter coefficients of FIR filter
%                                      used for sub sampling (optional)
%                   .data            - A matrix where the data is stored
%                                      (optional)
%                   .dataPos         - The position in the matrix
%                                      [dataStart dataEnd fileStart
%                                      fileEnd](optional)
%
%      The filter parts of the OPT structure are optional fields.
%      The default for the filt_subsample is a filter which takes the last
%      value of filtered block e.g. [0 ... 0 1]
%
%      With opt.data and opt.dataPos read_bv can write directly to a
%      matrix. dataPos is an optional value for opt.data where you can set
%      the position of the read data. dataStart is the position in data
%      where the first read datasample is stored.
%
%       Please note, that the fields chanidx and dataPos used as c indices
%       starting at 0.
% RETURNS
%          data: [nChans, len] the actual data
%
% DESCRIPTION
%    Open a file and read the eeg data. The data is filtered with an IIR
%    filter and an FIR filter.
%
%    For the eeg file we assume that it was written with the following
%    settings: DataFormat = BINARY
%              DataOrientation = MULTIPLEXED
%              BinaryFormat = INT_16
%
% COMPILE WITH
%    mex read_bv.c
%
% AUTHOR
%    Max Sagebaum
%
%    2008/04/15 - Max Sagebaum
%                   - file created
% (c) 2005 Fraunhofer FIRST

author: bbci


text:
function [net, stats] = cnn_train(net, imdb, getBatch, varargin)
%CNN_TRAIN  An example implementation of SGD for training CNNs
%    CNN_TRAIN() is an example learner implementing stochastic
%    gradient descent with momentum to train a CNN. It can be used
%    with different datasets and tasks by providing a suitable
%    getBatch function.
%
%    The function automatically restarts after each training epoch by
%    checkpointing.
%
%    The function supports training on CPU or on one or more GPUs
%    (specify the list of GPU IDs in the `gpus` option).

% Copyright (C) 2014-16 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.expDir = fullfile('data','exp') ;
opts.continue = true ;
opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.train = [] ;
opts.val = [] ;
opts.gpus = [] ;
opts.prefetch = false ;
opts.numEpochs = 300 ;
opts.learningRate = 0.001 ;
opts.weightDecay = 0.0005 ;
opts.momentum = 0.9 ;
opts.saveMomentum = true ;
opts.nesterovUpdate = false ;
opts.randomSeed = 0 ;
opts.memoryMapFile = fullfile(tempdir, 'matconvnet.bin') ;
opts.profile = false ;
opts.parameterServer.method = 'mmap' ;
opts.parameterServer.prefix = 'mcn' ;

opts.conserveMemory = true ;
opts.backPropDepth = +inf ;
opts.sync = false ;
opts.cudnn = true ;
opts.errorFunction = 'multiclass' ;
opts.errorLabels = {} ;
opts.plotDiagnostics = false ;
opts.plotStatistics = true;
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isnan(opts.train), opts.train = [] ; end
if isnan(opts.val), opts.val = [] ; end

% -------------------------------------------------------------------------
%                                                            Initialization
% -------------------------------------------------------------------------

net = vl_simplenn_tidy(net); % fill in some eventually missing values
net.layers{end-1}.precious = 1; % do not remove predictions, used for error
vl_simplenn_display(net, 'batchSize', opts.batchSize) ;

evaluateMode = isempty(opts.train) ;
if ~evaluateMode
  for i=1:numel(net.layers)
    J = numel(net.layers{i}.weights) ;
    if ~isfield(net.layers{i}, 'learningRate')
      net.layers{i}.learningRate = ones(1, J) ;
    end
    if ~isfield(net.layers{i}, 'weightDecay')
      net.layers{i}.weightDecay = ones(1, J) ;
    end
  end
end

% setup error calculation function
hasError = true ;
if isstr(opts.errorFunction)
  switch opts.errorFunction
    case 'none'
      opts.errorFunction = @error_none ;
      hasError = false ;
    case 'multiclass'
      opts.errorFunction = @error_multiclass ;
      if isempty(opts.errorLabels), opts.errorLabels = {'top1err', 'top5err'} ; end
    case 'binary'
      opts.errorFunction = @error_binary ;
      if isempty(opts.errorLabels), opts.errorLabels = {'binerr'} ; end
    otherwise
      error('Unknown error function ''%s''.', opts.errorFunction) ;
  end
end

state.getBatch = getBatch ;
stats = [] ;

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
  fprintf('%s: resuming by loading epoch %d\n', mfilename, start) ;
  [net, state, stats] = loadState(modelPath(start)) ;
else
  state = [] ;
end

for epoch=start+1:opts.numEpochs

  % Set the random seed based on the epoch and opts.randomSeed.
  % This is important for reproducibility, including when training
  % is restarted from a checkpoint.

  rng(epoch + opts.randomSeed) ;
  prepareGPUs(opts, epoch == start+1) ;

  % Train for one epoch.
  params = opts ;
  params.epoch = epoch ;
  params.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
  params.train = opts.train(randperm(numel(opts.train))) ; % shuffle
  params.val = opts.val(randperm(numel(opts.val))) ;
  params.imdb = imdb ;
  params.getBatch = getBatch ;

  if numel(params.gpus) <= 1
    [net, state] = processEpoch(net, state, params, 'train') ;
    [net, state] = processEpoch(net, state, params, 'val') ;
    if ~evaluateMode
      saveState(modelPath(epoch), net, state) ;
    end
    lastStats = state.stats ;
  else
    spmd
      [net, state] = processEpoch(net, state, params, 'train') ;
      [net, state] = processEpoch(net, state, params, 'val') ;
      if labindex == 1 && ~evaluateMode
        saveState(modelPath(epoch), net, state) ;
      end
      lastStats = state.stats ;
    end
    lastStats = accumulateStats(lastStats) ;
  end

  stats.train(epoch) = lastStats.train ;
  stats.val(epoch) = lastStats.val ;
  clear lastStats ;
  saveStats(modelPath(epoch), stats) ;

  if params.plotStatistics
    switchFigure(1) ; clf ;
    plots = setdiff(...
      cat(2,...
      fieldnames(stats.train)', ...
      fieldnames(stats.val)'), {'num', 'time'}) ;
    for p = plots
      p = char(p) ;
      values = zeros(0, epoch) ;
      leg = {} ;
      for f = {'train', 'val'}
        f = char(f) ;
        if isfield(stats.(f), p)
          tmp = [stats.(f).(p)] ;
          values(end+1,:) = tmp(1,:)' ;
          leg{end+1} = f ;
        end
      end
      subplot(1,numel(plots),find(strcmp(p,plots))) ;
      plot(1:epoch, values','o-') ;
      xlabel('epoch') ;
      title(p) ;
      legend(leg{:}) ;
      grid on ;
    end
    drawnow ;
    print(1, modelFigPath, '-dpdf') ;
  end
end

% With multiple GPUs, return one copy
if isa(net, 'Composite'), net = net{1} ; end

% -------------------------------------------------------------------------
function err = error_multiclass(params, labels, res)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
[~,predictions] = sort(predictions, 3, 'descend') ;

% be resilient to badly formatted labels
if numel(labels) == size(predictions, 4)
  labels = reshape(labels,1,1,1,[]) ;
end

% skip null labels
mass = single(labels(:,:,1,:) > 0) ;
if size(labels,3) == 2
  % if there is a second channel in labels, used it as weights
  mass = mass .* labels(:,:,2,:) ;
  labels(:,:,2,:) = [] ;
end

m = min(5, size(predictions,3)) ;

error = ~bsxfun(@eq, predictions, labels) ;
err(1,1) = sum(sum(sum(mass .* error(:,:,1,:)))) ;
err(2,1) = sum(sum(sum(mass .* min(error(:,:,1:m,:),[],3)))) ;

% -------------------------------------------------------------------------
function err = error_binary(params, labels, res)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
error = bsxfun(@times, predictions, labels) < 0 ;
err = sum(error(:)) ;

% -------------------------------------------------------------------------
function err = error_none(params, labels, res)
% -------------------------------------------------------------------------
err = zeros(0,1) ;

% -------------------------------------------------------------------------
function [net, state] = processEpoch(net, state, params, mode)
% -------------------------------------------------------------------------
% Note that net is not strictly needed as an output argument as net
% is a handle class. However, this fixes some aliasing issue in the
% spmd caller.

% initialize with momentum 0
if isempty(state) || isempty(state.momentum)
  for i = 1:numel(net.layers)
    for j = 1:numel(net.layers{i}.weights)
      state.momentum{i}{j} = 0 ;
    end
  end
end

% move CNN  to GPU as needed
numGpus = numel(params.gpus) ;
if numGpus >= 1
  net = vl_simplenn_move(net, 'gpu') ;
  for i = 1:numel(state.momentum)
    for j = 1:numel(state.momentum{i})
      state.momentum{i}{j} = gpuArray(state.momentum{i}{j}) ;
    end
  end
end
if numGpus > 1
  parserv = ParameterServer(params.parameterServer) ;
  vl_simplenn_start_parserv(net, parserv) ;
else
  parserv = [] ;
end

% profile
if params.profile
  if numGpus <= 1
    profile clear ;
    profile on ;
  else
    mpiprofile reset ;
    mpiprofile on ;
  end
end

subset = params.(mode) ;
num = 0 ;
stats.num = 0 ; % return something even if subset = []
stats.time = 0 ;
adjustTime = 0 ;
res = [] ;
error = [] ;

start = tic ;
for t=1:params.batchSize:numel(subset)
  fprintf('%s: epoch %02d: %3d/%3d:', mode, params.epoch, ...
          fix((t-1)/params.batchSize)+1, ceil(numel(subset)/params.batchSize)) ;
  batchSize = min(params.batchSize, numel(subset) - t + 1) ;

  for s=1:params.numSubBatches
    % get this image batch and prefetch the next
    batchStart = t + (labindex-1) + (s-1) * numlabs ;
    batchEnd = min(t+params.batchSize-1, numel(subset)) ;
    batch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
    num = num + numel(batch) ;
    if numel(batch) == 0, continue ; end

    [im, labels] = params.getBatch(params.imdb, batch) ;

    if params.prefetch
      if s == params.numSubBatches
        batchStart = t + (labindex-1) + params.batchSize ;
        batchEnd = min(t+2*params.batchSize-1, numel(subset)) ;
      else
        batchStart = batchStart + numlabs ;
      end
      nextBatch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
      params.getBatch(params.imdb, nextBatch) ;
    end

    if numGpus >= 1
      im = gpuArray(im) ;
    end

    if strcmp(mode, 'train')
      dzdy = 1 ;
      evalMode = 'normal' ;
    else
      dzdy = [] ;
      evalMode = 'test' ;
    end
    net.layers{end}.class = labels ;
    res = vl_simplenn(net, im, dzdy, res, ...
                      'accumulate', s ~= 1, ...
                      'mode', evalMode, ...
                      'conserveMemory', params.conserveMemory, ...
                      'backPropDepth', params.backPropDepth, ...
                      'sync', params.sync, ...
                      'cudnn', params.cudnn, ...
                      'parameterServer', parserv, ...
                      'holdOn', s < params.numSubBatches) ;

    % accumulate errors
    error = sum([error, [...
      sum(double(gather(res(end).x))) ;
      reshape(params.errorFunction(params, labels, res),[],1) ; ]],2) ;
  end

  % accumulate gradient
  if strcmp(mode, 'train')
    if ~isempty(parserv), parserv.sync() ; end
    [net, res, state] = accumulateGradients(net, res, state, params, batchSize, parserv) ;
  end

  % get statistics
  time = toc(start) + adjustTime ;
  batchTime = time - stats.time ;
  stats = extractStats(net, params, error / num) ;
  stats.num = num ;
  stats.time = time ;
  currentSpeed = batchSize / batchTime ;
  averageSpeed = (t + batchSize - 1) / time ;
  if t == 3*params.batchSize + 1
    % compensate for the first three iterations, which are outliers
    adjustTime = 4*batchTime - time ;
    stats.time = time + adjustTime ;
  end

  fprintf(' %.1f (%.1f) Hz', averageSpeed, currentSpeed) ;
  for f = setdiff(fieldnames(stats)', {'num', 'time'})
    f = char(f) ;
    fprintf(' %s: %.3f', f, stats.(f)) ;
  end
  fprintf('\n') ;

  % collect diagnostic statistics
  if strcmp(mode, 'train') && params.plotDiagnostics
    switchFigure(2) ; clf ;
    diagn = [res.stats] ;
    diagnvar = horzcat(diagn.variation) ;
    diagnpow = horzcat(diagn.power) ;
    subplot(2,2,1) ; barh(diagnvar) ;
    set(gca,'TickLabelInterpreter', 'none', ...
      'YTick', 1:numel(diagnvar), ...
      'YTickLabel',horzcat(diagn.label), ...
      'YDir', 'reverse', ...
      'XScale', 'log', ...
      'XLim', [1e-5 1], ...
      'XTick', 10.^(-5:1)) ;
    grid on ;
    subplot(2,2,2) ; barh(sqrt(diagnpow)) ;
    set(gca,'TickLabelInterpreter', 'none', ...
      'YTick', 1:numel(diagnpow), ...
      'YTickLabel',{diagn.powerLabel}, ...
      'YDir', 'reverse', ...
      'XScale', 'log', ...
      'XLim', [1e-5 1e5], ...
      'XTick', 10.^(-5:5)) ;
    grid on ;
    subplot(2,2,3); plot(squeeze(res(end-1).x)) ;
    drawnow ;
  end
end

% Save back to state.
state.stats.(mode) = stats ;
if params.profile
  if numGpus <= 1
    state.prof.(mode) = profile('info') ;
    profile off ;
  else
    state.prof.(mode) = mpiprofile('info');
    mpiprofile off ;
  end
end
if ~params.saveMomentum
  state.momentum = [] ;
else
  for i = 1:numel(state.momentum)
    for j = 1:numel(state.momentum{i})
      state.momentum{i}{j} = gather(state.momentum{i}{j}) ;
    end
  end
end

net = vl_simplenn_move(net, 'cpu') ;

% -------------------------------------------------------------------------
function [net, res, state] = accumulateGradients(net, res, state, params, batchSize, parserv)
% -------------------------------------------------------------------------
numGpus = numel(params.gpus) ;
otherGpus = setdiff(1:numGpus, labindex) ;

for l=numel(net.layers):-1:1
  for j=numel(res(l).dzdw):-1:1

    if ~isempty(parserv)
      tag = sprintf('l%d_%d',l,j) ;
      parDer = parserv.pull(tag) ;
    else
      parDer = res(l).dzdw{j}  ;
    end

    if j == 3 && strcmp(net.layers{l}.type, 'bnorm')
      % special case for learning bnorm moments
      thisLR = net.layers{l}.learningRate(j) ;
      net.layers{l}.weights{j} = vl_taccum(...
        1 - thisLR, ...
        net.layers{l}.weights{j}, ...
        thisLR / batchSize, ...
        parDer) ;
    else
      % Standard gradient training.
      thisDecay = params.weightDecay * net.layers{l}.weightDecay(j) ;
      thisLR = params.learningRate * net.layers{l}.learningRate(j) ;

      if thisLR>0 || thisDecay>0
        % Normalize gradient and incorporate weight decay.
        parDer = vl_taccum(1/batchSize, parDer, ...
                           thisDecay, net.layers{l}.weights{j}) ;

        % Update momentum.
        state.momentum{l}{j} = vl_taccum(...
          params.momentum, state.momentum{l}{j}, ...
          -1, parDer) ;

        % Nesterov update (aka one step ahead).
        if params.nesterovUpdate
          delta = vl_taccum(...
            params.momentum, state.momentum{l}{j}, ...
            -1, parDer) ;
        else
          delta = state.momentum{l}{j} ;
        end

        % Update parameters.
        net.layers{l}.weights{j} = vl_taccum(...
          1, net.layers{l}.weights{j}, ...
          thisLR, delta) ;
      end
    end

    % if requested, collect some useful stats for debugging
    if params.plotDiagnostics
      variation = [] ;
      label = '' ;
      switch net.layers{l}.type
        case {'conv','convt'}
          variation = thisLR * mean(abs(state.momentum{l}{j}(:))) ;
          power = mean(res(l+1).x(:).^2) ;
          if j == 1 % fiters
            base = mean(net.layers{l}.weights{j}(:).^2) ;
            label = 'filters' ;
          else % biases
            base = sqrt(power) ;%mean(abs(res(l+1).x(:))) ;
            label = 'biases' ;
          end
          variation = variation / base ;
          label = sprintf('%s_%s', net.layers{l}.name, label) ;
      end
      res(l).stats.variation(j) = variation ;
      res(l).stats.power = power ;
      res(l).stats.powerLabel = net.layers{l}.name ;
      res(l).stats.label{j} = label ;
    end
  end
end

% -------------------------------------------------------------------------
function stats = accumulateStats(stats_)
% -------------------------------------------------------------------------

for s = {'train', 'val'}
  s = char(s) ;
  total = 0 ;

  % initialize stats stucture with same fields and same order as
  % stats_{1}
  stats__ = stats_{1} ;
  names = fieldnames(stats__.(s))' ;
  values = zeros(1, numel(names)) ;
  fields = cat(1, names, num2cell(values)) ;
  stats.(s) = struct(fields{:}) ;

  for g = 1:numel(stats_)
    stats__ = stats_{g} ;
    num__ = stats__.(s).num ;
    total = total + num__ ;

    for f = setdiff(fieldnames(stats__.(s))', 'num')
      f = char(f) ;
      stats.(s).(f) = stats.(s).(f) + stats__.(s).(f) * num__ ;

      if g == numel(stats_)
        stats.(s).(f) = stats.(s).(f) / total ;
      end
    end
  end
  stats.(s).num = total ;
end

% -------------------------------------------------------------------------
function stats = extractStats(net, params, errors)
% -------------------------------------------------------------------------
stats.objective = errors(1) ;
for i = 1:numel(params.errorLabels)
  stats.(params.errorLabels{i}) = errors(i+1) ;
end

% -------------------------------------------------------------------------
function saveState(fileName, net, state)
% -------------------------------------------------------------------------
save(fileName, 'net', 'state') ;

% -------------------------------------------------------------------------
function saveStats(fileName, stats)
% -------------------------------------------------------------------------
if exist(fileName)
  save(fileName, 'stats', '-append') ;
else
  save(fileName, 'stats') ;
end

% -------------------------------------------------------------------------
function [net, state, stats] = loadState(fileName)
% -------------------------------------------------------------------------
load(fileName, 'net', 'state', 'stats') ;
net = vl_simplenn_tidy(net) ;
if isempty(whos('stats'))
  error('Epoch ''%s'' was only partially saved. Delete this file and try again.', ...
        fileName) ;
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

% -------------------------------------------------------------------------
function switchFigure(n)
% -------------------------------------------------------------------------
if get(0,'CurrentFigure') ~= n
  try
    set(0,'CurrentFigure',n) ;
  catch
    figure(n) ;
  end
end

% -------------------------------------------------------------------------
function clearMex()
% -------------------------------------------------------------------------
%clear vl_tmove vl_imreadjpeg ;
disp('Clearing mex files') ;
clear mex ;
clear vl_tmove vl_imreadjpeg ;

% -------------------------------------------------------------------------
function prepareGPUs(params, cold)
% -------------------------------------------------------------------------
numGpus = numel(params.gpus) ;
if numGpus > 1
  % check parallel pool integrity as it could have timed out
  pool = gcp('nocreate') ;
  if ~isempty(pool) && pool.NumWorkers ~= numGpus
    delete(pool) ;
  end
  pool = gcp('nocreate') ;
  if isempty(pool)
    parpool('local', numGpus) ;
    cold = true ;
  end
end
if numGpus >= 1 && cold
  fprintf('%s: resetting GPU\n', mfilename) ;
  clearMex() ;
  if numGpus == 1
    disp(gpuDevice(params.gpus)) ;
  else
    spmd
      clearMex() ;
      disp(gpuDevice(params.gpus(labindex))) ;
    end
  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
%--- help for generic/setup_priors ---
%
%  INTERNAL FUNCTION: Format the priors for the estimation of models in RISE.
%
%  ::
%
%    obj=setup_priors(obj,MyPriors)
%    obj=setup_priors(obj,MyPriors,error_control)
%
%  Args:
%
%     obj (rise | dsge | svar | rfvar): model object. This function is not
%       meant to be called directly by the user. But it shows how to write the
%       priors when those are not written inside the model file.
%
%     MyPriors (struct): Each field of the structure is one of the
%       following alternatives
%
%       - Maximum likelihood or uniform priors::
%
%           P.pname={start_value,lower_bound,upper_bound};
%
%       - Bayesian prior using mean and standard deviation::
%
%           P.pname={start_value,prior_mean,prior_stdev,'distribution'};
%
%       - Same as above except for adding a hard lower bound::
%
%           P.pname={start_value,prior_mean,prior_stdev,'distribution',lower_bound};
%
%       - Same as above except for adding a hard upper bound::
%
%           P.pname={start_value,prior_mean,prior_stdev,'distribution',lower_bound,upper_bound};
%
%       - Bayesian prior using quantiles of the distribution::
%
%           P.pname={start_value,lower_quantile,lower_quantile,'distribution(prob)'};
%
%       - Same as above except for adding a hard lower bound::
%
%           P.pname={start_value,lower_quantile,lower_quantile,'distribution(prob)',lower_bound};
%
%       - Same as above except for adding a hard upper bound::
%
%           P.pname={start_value,lower_quantile,lower_quantile,'distribution(prob)',lower_bound,upper_bound};
%
%       - Dirichlet
%
%           Denote by dirichlet_j, the jth dirichlet distribution, j=1,2,... we
%           have P.dirichlet_j={sd_ii,pname1,m1,pname2,m2,...,pnamen,mn} where
%
%             - sd_ii: is the standard deviation of the diagonal element (a
%               parameter never listed by RISE) of the dirichlet distribution.
%             - pname1, pname2,...,pnamen: are the names of the off diagonal
%               parameters
%             - m1, m2,...,mn: are the means of each parameters
%
%     error_control (empty | cell): element constructed by RISE for
%       controling the syntax used in the model file.
%
%  Returns:
%     :
%
%     - **obj** [rise|dsge|svar|rfvar]: updated model object.
%
%  Note:
%
%     - This function is also indirectly used for svar and rfvar objects.
%     - Possible distributions include:
%
%        - beta
%        - cauchy
%        - gamma
%        - inv_gamma
%        - laplace
%        - left_triang
%        - logistic
%        - lognormal
%        - normal
%        - pareto
%        - right_triang
%        - uniform
%        - weibull
%
%  See also:
%     - rise_generic/estimate
%     - dsge/estimate
%
%

author: jmaih


text:
% Profiler extension for Normalized Leaky Kernel Affine Projection
% Algorithm
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

classdef nlkapa_profiler < nlkapa

    properties (GetAccess = 'public', SetAccess = 'private')
        elapsed = 0; % elapsed time
    end

    methods

        function kaf = nlkapa_profiler(parameters) % constructor
            if nargin<1, parameters = struct(); end
            kaf = kaf@nlkapa(parameters);
        end

        function flops = lastflops(kaf) % flops for last iteration
            % numbers of operations
            m1 = 1;
            m2 = 1;
            m3 = 1;
            m4 = 1;

            floptions = struct(...
                'sum', m1, ...
                'mult', m2, ...
                'div', m3, ...
                sprintf('%s_kernel',kaf.kerneltype), [m4,1,size(kaf.dict,2)]);

            flops = nan; % dummy. kflops(floptions);
        end

        %% flops breakdown

        % [space for remarks on number of operations used above]

        %%

        function train_profiled(kaf,x,y)
            t1 = tic;
            kaf.train(x,y);
            t2 = toc(t1);
            kaf.elapsed = kaf.elapsed + t2;
        end

        function bytes = lastbytes(kaf) % bytes used in last iteration
            m = size(kaf.dict,1);
            m2 = 1;
            bytes = nan; % dummy. (m2 + m*size(kaf.dict,2)); % 8 bytes for double precision
            % [list variables]
        end

    end
end

repo: kafbox


text:
function kernel = MimNormalisedGaussianKernel(voxel_size_mm, filter_size_mm, minimum_grid_size_mm)
    % MimNormalisedGaussianKernel.
    %
    %
    %     The input and output images are of class PTKImage.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    if nargin < 3
        minimum_grid_size_mm = [];
    end

    if numel(minimum_grid_size_mm) == 1
        minimum_grid_size_mm = repmat(minimum_grid_size_mm, [1, 3]);
    end

    sigma_mm = filter_size_mm;

    epsilon = 1e-3;
    sigma_voxels = sigma_mm./voxel_size_mm;
    grid_size = 2*(ceil((sigma_voxels).*sqrt(-2*log(sqrt(2*pi).*(sigma_voxels)*epsilon)))) + 1;

    if ~isempty(minimum_grid_size_mm)
        minimum_grid_size = 2*(ceil((minimum_grid_size_mm./voxel_size_mm)/2));
        grid_size = max(grid_size, minimum_grid_size);
    end

    grid_size_i = grid_size(1);
    grid_size_j = grid_size(2);
    grid_size_k = grid_size(3);

    center_i = grid_size_i/2 + 0.5;
    center_j = grid_size_j/2 + 0.5;
    center_k = grid_size_k/2 + 0.5;

    n_i = 1 : grid_size_i;
    n_j = 1 : grid_size_j;
    n_k = 1 : grid_size_k;

    sigmai = sigma_voxels(1);
    sigmaj = sigma_voxels(2);
    sigmak = sigma_voxels(3);

    keri = zeros(1, 1, grid_size_i, 'single');
    kerj = zeros(1, 1, grid_size_j, 'single');
    kerk = zeros(1, 1, grid_size_k, 'single');

    keri(1,1,:) = (1/((2*pi*sigmai.^2).^(1/2))) * exp(-((n_i - center_i).^2)/(2*sigmai.^2));
    kerj(1,1,:) = (1/((2*pi*sigmaj.^2).^(1/2))) * exp(-((n_j - center_j).^2)/(2*sigmaj.^2));
    kerk(1,1,:) = (1/((2*pi*sigmak.^2).^(1/2))) * exp(-((n_k - center_k).^2)/(2*sigmak.^2));

    % Normalise
    keri = keri./sum(keri);
    kerj = kerj./sum(kerj);
    kerk = kerk./sum(kerk);

    ker1 = repmat(shiftdim(keri, 2), [1, grid_size_j, grid_size_k]);
    ker2 = repmat(shiftdim(kerj, 1), [grid_size_i, 1, grid_size_k]);
    ker3 = repmat(kerk, [grid_size_i, grid_size_j, 1]);
    kernel = ker1.*ker2.*ker3;

    kernel = kernel/max(kernel(:));
end

author: tomdoel

repo: pulmonarytoolkit


text:
function mr=rotqr2mr(qr)
%ROTQR2MR converts a matrix of real quaternion vectors to quaternion matrices
% Inputs:
%
%     QR(4m,n)   mxn matrix of real quaternion vectors (each 4x1)
%
% Outputs:
%
%     MR(4m,4n)   mxn matrix of real quaternion matrices (each 4x4)
%
% In matrix form, quaternions can be multiplied and added using normal matrix
% arithmetic. Each element of an mxn matrix of quaternions is itself a 4x4 block
% so the total dimension of MR is 4m x 4n.

%
%      Copyright (C) Mike Brookes 2000-2006
%      Version: $Id: rotqr2mr.m 1615 2012-03-15 09:10:51Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent a b c
if isempty(a)
    a=[1 2 3 3 1 2];    % destination row of +ve entries (from 0)
    b=[1 2 3 2 3 1];    % destination col of +ve entries (from 0)
    c=[0 0 0 1 2 3];    % source row of +ve entries (from 0)
end
[m,n]=size(qr);
mr=repmat(qr,4,1);
mn=m*n;
j=repmat(4*m*(0:n-1),m/4,1);
i=repmat((1:4:m)',n,1)+j(:);
ni=length(i);
i6=repmat(i,1,6);
mr(i6+repmat(a+m*b,ni,1))=mr(i6+repmat(c,ni,1));
mr(i6+repmat(c+m*b,ni,1))=-mr(i6+repmat(a,ni,1));
mr=reshape(mr,m,4*n);

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function [pos vel] = IntKalman(z)
%
%
persistent A H Q R
persistent x P
persistent firstRun

if isempty(firstRun)
  firstRun = 1;

  dt = 0.1;

  A = [ 1 dt;
        0 1  ];
  H = [0 1];

  Q = [ 1 0;
        0 3 ];
  R = 10;

  x = [ 0 20 ]';
  P = 5*eye(2);
end

xp = A*x;
Pp = A*P*A' + Q;

K = Pp*H'*inv(H*Pp*H' + R);

x = xp + K*(z - H*xp);
P = Pp - K*H*Pp;

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

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
% SIS模型演示程序
clear;clc;
ts = 0 : 150;

x0=[0.2];   % 最初感染率，百万分之一

% lambda   日接触率
%  miu     日治愈率，即每天治愈病人数占病人总数的比例
lambda_miu = [0.27, 0.1;...
              0.25, 0.1;...
              0.20, 0.1;...
              0.15, 0.1;...
              0.08, 0.1];
infective_matrix = zeros(length(ts),length(lambda_miu));

for i = 1 : length(lambda_miu)
  lambda = lambda_miu(i,1);
  miu = lambda_miu(i,2);
  [t,x] = ode45(@(t,y) ill_sis(t,y,lambda,miu), ts, x0);
  infective_matrix(:,i) = x(:,1);
end

plot(t, infective_matrix);
legend(num2str(lambda_miu(1,:)), num2str(lambda_miu(2,:)), num2str(lambda_miu(3,:)), num2str(lambda_miu(4,:)), num2str(lambda_miu(5,:)));
title('SIS模型');
grid;

repo: NumericalAnalysis


text:
function chi = imEuler3dDensity(img, varargin)
% Compute Euler density in a 3D image.
%
%   CHI_V = imEuler3dDensity(IMG)
%   Compute Euler number estimate in a 3D image, and normalize by the
%   observed volume. This function is well suited for estimating
%   topological properties of a random material observed through a sampling
%   window.
%
%   CHI_V = imEuler3dDensity(IMG, CONN)
%   Specifies the connectivity to use to define whether two voxels are
%   neihbors or not. Can be either 6 (the default) or 26.
%
%   Example
%   imEuler3dDensity
%
%   See also
%     imEuler3d, imEuler3dEstimate, imSurfaceAreaDensity, imMeanBreadthDensity

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2010-07-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRAE - Cepia Software Platform.

% check image dimension
if ndims(img) ~= 3
    error('first argument should be a 3D image');
end

% in case of a label image, return a vector with a set of results
if ~islogical(img)
    labels = unique(img);
    labels(labels==0) = [];
    epcd = zeros(length(labels), 1);
    for i = 1:length(labels)
        epcd(i) = imEuler3dDensity(img==labels(i), varargin{:});
    end
    return;
end

% Process user input arguments
delta = [1 1 1];
conn = 6;
while ~isempty(varargin)
    var = varargin{1};
    if ~isnumeric(var)
        error('option should be numeric');
    end

    % option is either connectivity or resolution
    if isscalar(var)
        conn = var;
    else
        delta = var;
    end
    varargin(1) = [];
end

% Euler-Poincare Characteristic of each component in image
chi = imEuler3dEstimate(img, conn);

% total volume of image
obsVolume = prod(size(img) - 1) * prod(delta);

% compute area density
chi = chi / obsVolume;


text:
% Demonstration of the basic operation of a kernel adaptive filtering
% algorithm with fixed budget. This script generates several images and
% saves them to disk.
%
% KRLS-T is chosen as the example algorithm. Inspired by
% http://www.tsc.uc3m.es/~miguel/MLG/adjuntos/slidesKRLST.pdf
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

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

%% PARAMETERS

n = 10; % number of training data
n_star = 201; % number of points to evaluate

ell = 2; % length scale
sigvar = 1; % signal variance
noisevar = .02; % noise variance

budget = 5; % KRLS-T budget

savefig = 0;
manual = 0;

%% GENERATE DATA

fprintf('Generating data...\n');
covfunc = {'kafbox_covSum', {'kafbox_covSEiso','kafbox_covNoise'}};
loghyper = [log(ell); log(sigvar); log(noisevar)];

x = (1:n)' + 0.2*randn(n,1);
y = chol(feval(covfunc{:}, loghyper, x))'*randn(n,1);      % Cholesky decomp.

x_star = linspace(0,n+2,n_star)';

covfunc_no_noise = {'kafbox_covSEiso'};
loghyper_no_noise = [log(ell); log(sigvar)];

%% RUN ALGORITHM

kaf = krlst_split(struct('kerneltype','gauss','kernelpar',ell,...
    'sn2',noisevar/sigvar,'M',budget+1));

fprintf('Running KRLS-T\n');
t1 = tic;
for i=1:n
    if ~mod(i,floor(n/10)), fprintf('.'); end
    Qold = kaf.Q;
    kaf.train(x(i),y(i));

    for stage = 1:2
        switch stage
            case 1 % stage 1: after including a sample
                stage_string = 'after inclusion';
                fname_i = sprintf('fig/kafbox_sample_run_%02da',i);
            case 2 % stage 2: after pruning
                if length(kaf.mu) > budget
                    kaf = kaf.prune(Qold,1);
                    stage_string = 'after pruning  ';
                    fname_i = sprintf('fig/kafbox_sample_run_%02db',i);
                end
        end

        y_est = kaf.evaluate(x(1:i));
        [y_star_est,y_star_var] = kaf.evaluate(x_star);

        y_star_var = y_star_var/kaf.s02*sigvar; % fix scale

        x_dict = kaf.dict;
        y_dict = kaf.evaluate(x_dict);

        papersize = [8 4];

        figure(1); clf
        set(gcf,'Position',[500 400 papersize(1)*100 papersize(2)*100])
        hold all
        f = [y_star_est+2*sqrt(y_star_var); flipud(y_star_est-2*sqrt(y_star_var))];
        p0 = fill([x_star; flipud(x_star)], f, [7 7 7]/8,'EdgeColor','none');
        p1 = plot(x_star,y_star_est,'r','LineWidth',1.5);
        p2 = plot(x_dict,y_dict,'ro','LineWidth',2,'MarkerSize',8);
        p3 = plot(x(1:i),y(1:i),'bx','LineWidth',2,'MarkerSize',8);

        set(gca,'box','on');
        set(gca, 'FontSize', 14);
        set(gca, 'FontName', 'Helvetica');

        axis([0 n+2 floor(min(y))-1 ceil(max(y))+1])
        legend([p3,p1,p2],'data','prediction','dictionary bases')

        title(sprintf('KRLS-T iteration %d: %s',i,stage_string));

        set(gca, 'Position', get(gca, 'OuterPosition') - ...
            get(gca, 'TightInset') * [-1 0 1 0; 0 -1 0 1; 0 0 1 0; 0 0 0 1]);

        set(gcf, 'PaperUnits', 'inches');
        set(gcf, 'PaperSize', papersize);
        set(gcf, 'PaperPositionMode', 'manual');
        set(gcf, 'PaperPosition', [0 0 papersize(1) papersize(2)]);
        set(gcf, 'renderer', 'painters');

        drawnow
        if manual
            keyboard; %#ok<UNRCH>
        end

        if savefig
            if ~exist('fig','file') %#ok<UNRCH>
                mkdir('fig');
            end
            % write individual png
            print('-dpng','-r100',fname_i);
            print('-dpdf','-r100',fname_i);
            % write animated gif
            fname = 'fig/kafbox_sample_run.gif';
            frame = getframe(1);
            im = frame2im(frame);
            [imind,cm] = rgb2ind(im,256);
            if i==1 && stage==1
                imwrite(imind,cm,fname,'gif','Loopcount',inf,'DelayTime',1.5);
            else
                if stage==2 && i==n && kaf.pruned
                    imwrite(imind,cm,fname,'gif','WriteMode','append','DelayTime',2.5);
                elseif ~(stage==2 && ~kaf.pruned)
                    imwrite(imind,cm,fname,'gif','WriteMode','append','DelayTime',1.5);
                end
            end
        end
    end

end
fprintf('\n')


text:
classdef ZDT3 < PROBLEM
% <multi> <real> <large/none> <expensive/none>
% Benchmark MOP proposed by Zitzler, Deb, and Thiele

%------------------------------- Reference --------------------------------
% E. Zitzler, K. Deb, and L. Thiele, Comparison of multiobjective
% evolutionary algorithms: Empirical results, Evolutionary computation,
% 2000, 8(2): 173-195.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj(:,1) = PopDec(:,1);
            g = 1 + 9*mean(PopDec(:,2:end),2);
            h = 1 - (PopObj(:,1)./g).^0.5 - PopObj(:,1)./g.*sin(10*pi*PopObj(:,1));
            PopObj(:,2) = g.*h;
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - R(:,1).^0.5 - R(:,1).*sin(10*pi*R(:,1));
            R      = R(NDSort(R,1)==1,:);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R(:,1) = linspace(0,1,100)';
            R(:,2) = 1 - R(:,1).^0.5 - R(:,1).*sin(10*pi*R(:,1));
            R(NDSort(R,1)>1,:) = nan;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% AUTHOR :   CHANDAN KUMAR
% SUBJECT:   COMPUTER GRAPHICS AND SOLID MODELLING.
% TITLE  :   LINE CLIPPING ALGORITHM (CYRUS BECK )
% DATE   :   31/08/09
function lineclip()
clc, clear all;
x = input('Give coord[ x0 y0 x1 y1]: ');
v = input('Give Viewport coord[vx0 vy0 vx1 vy1]: ');
d = [(x(3)-x(1)) (x(4)-x(2)) -(x(3)-x(1))  -(x(4)-x(2))];
for i = 1:4
    if (i==1)||(i==3)
        j =1;
    else
        j =2;
    end
    t(i) = (v(i)-x(j))/d(j);
end
tl = min(t(1),t(3));
ti = max(t(2),t(4));

z(1) = x(1) +d(1)*ti ;   % x coord of clipped line
z(3) = x(1) +d(1)*tl ;   % end X coord of clipped line
z(2) = x(2) +d(2)*ti ;
z(4) = x(2) +d(2)*tl ;

plot([z(1) z(3)] ,[z(2) z(4)],'r-','LineWidth',2);   % Plots Clipped line
hold on,grid on
plot([v(1) v(1)] ,[v(2) v(4)],'b-','LineWidth',2);   % Plots Viewport
plot([x(1) z(3)],[x(2) z(4)],'k-');                  % Plots Unclipped line
legend('Clipped line','Viewport','Unclipped line',3)
plot([z(1) x(3)],[z(2) x(4)],'k-');                  % Plots Unclipped line
plot([v(1) v(3)] ,[v(2) v(2)],'b-','LineWidth',2);   % Plots Viewport
plot([v(3) v(3)] ,[v(2) v(4)],'b-','LineWidth',2);   % Plots Viewport
plot([v(1) v(3)] ,[v(4) v(4)],'b-','LineWidth',2);   % Plots Viewport
axis([v(1)-1 v(3)+1 v(2)-1 v(4)+1]);
title('Line clipping by Cyrus Beck algorithm')

author: Sable

repo: mcbench-benchmarks


text:
function title = p43_title ( )

%*****************************************************************************80
%
%% P43_TITLE returns a title for problem 43.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 January 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The Himmelblau function.';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
close all;
%%
load('ExperimentInformation.mat');
load('NAT.mat');
%% paremater settings
SpeedBinning=1.5; %cm/s
MinSpeed=2.5;
MaxSpeed=16;
SpeedRange=[MinSpeed:SpeedBinning:MaxSpeed]; %cm0
MinSpanTime=10;%second
ShuffleNum=200;
Shuffling_mininterval=30;
SpeedSmooth=round(0.25*ExperimentInformation.Trackingframerate);
MinEventCount=100;

%%
SpeedTunning=cell(ExperimentInformation.Session,ExperimentInformation.TotalCell); % tuning curve
Speed_KStest=zeros(ExperimentInformation.Session,ExperimentInformation.TotalCell,ShuffleNum+3); % K-S value for shuffles (1 to ShuffleNum, 95th, 5th, and observed data)
Speed_Spearman=zeros(ExperimentInformation.Session,ExperimentInformation.TotalCell,ShuffleNum+3); % Spearman for shuffles (1 to ShuffleNum, 95th, 5th, and observed data)

%% K-S analysis and Spearmean test
% generate speed tuning curve with finner binning

for j=1:1:ExperimentInformation.Session
    k=1;
    for i=1:1:ExperimentInformation.TotalCell

        Event_filtered=intersect(find(NAT{1,j}(:,4*i+12)>0),find(NAT{1,j}(:,6)==1));% filter out the frames with speed valid
        Event_filtered=intersect(Event_filtered,find(NAT{1,j}(:,5)>MinSpeed));% filter out the frames with speed threadholds
       if length(Event_filtered)>MinEventCount && ExperimentInformation.CellSNR(i)>3 && ~ismember(i,ExperimentInformation.RepeatCell)
       SelectedFrame_raw=find(~isnan(NAT{1,j}(:,12+4*i))); %all the position are select
       SpeedTrain=NAT{1,1}(SelectedFrame_raw,5);
       SpeedTrain_smooth=general.smoothGauss(SpeedTrain,SpeedSmooth);
       SelectedFrame_filtered=logical((SpeedTrain_smooth>=MinSpeed).*(SpeedTrain_smooth<=MaxSpeed+SpeedBinning));
       SpeedTrain_filtered=SpeedTrain_smooth(SelectedFrame_filtered);
       EventTrain_original=NAT{1,j}(SelectedFrame_raw,4*i+12);
       EventTrain_smooth=general.smoothGauss(EventTrain_original,SpeedSmooth);
       EventTrain_filtered=EventTrain_smooth(SelectedFrame_filtered);
       SpeedTunning{j,i}=SpatialTuning_BNT.SpeedTuningCalcultation(SpeedTrain_filtered,EventTrain_filtered,...
                         ExperimentInformation.Trackingframerate./double(ExperimentInformation.ImagingPlane),...
                         SpeedRange,MinSpanTime,...
                         MaxSpeed+SpeedBinning);
       SpeedTunning{j,i}.Rate=smoothdata(SpeedTunning{k,1}.Rate,'lowess',3);
       SpeedCount=round(SpeedTunning{j,i}.Rate./max(SpeedTunning{j,i}.Rate)*100);
       KS_SpeedBinVector{k,1}=zeros(sum(SpeedCount(~isnan(SpeedCount))),1);
       m=1;
              for p=1:1:length(SpeedCount)
                  if ~isnan(SpeedCount(p))
                    for q=1:1:SpeedCount(p)
                            KS_SpeedBinVector{k,1}(m,1)=SpeedTunning{k,1}.SpeedRange(p);
                            m=m+1;
                    end
                  else
                  end
              end
      for n=1:1:ShuffleNum
          xmin=Shuffling_mininterval*ExperimentInformation.Trackingframerate;
          xmax=size(SelectedFrame_filtered,1)-xmin;
          ShiftFrame=round(xmin+rand(1,1)*(xmax-xmin));
          SpikeTrain_raw_shuffled=circshift(EventTrain_filtered,ShiftFrame);
          Speedtuning_shuffled=SpatialTuning_BNT.SpeedTuningCalcultation(SpeedTrain_filtered,SpikeTrain_raw_shuffled,...
                                ExperimentInformation.Trackingframerate./double(ExperimentInformation.ImagingPlane),...
                                SpeedTunning{k,1}.SpeedRange,MinSpanTime,17.5);
          RateSmoothed=smoothdata(Speedtuning_shuffled.Rate,'lowess',3);
          SpeedCount=round(RateSmoothed./max(SpeedTunning{k,1}.Rate)*100);

          % calculate the KS value for shuffles
          KS_SpeedBinVector_shuffle=zeros(sum(SpeedCount(~isnan(SpeedCount))),1);
          m=1;
          for p=1:1:length(SpeedCount)
              if ~isnan(SpeedCount(p))
                  for q=1:1:SpeedCount(p)
                      KS_SpeedBinVector_shuffle(m,1)=Speedtuning_shuffled.SpeedRange(p);
                      m=m+1;
                  end
              else
              end
          end
          [~,~,Speed_KStest(j,i,n)]=kstest2(KS_SpeedBinVector_shuffle,SpeedTunning{k,1}.SpeedRange);
          % calculate the spearman correlation
          speedScore_tem=SpatialTuning_BNT.speedScore(SpeedTrain_filtered,SpikeTrain_raw_shuffled(:,2),...
                         ExperimentInformation.Trackingframerate./double(ExperimentInformation.ImagingPlane));
          Speed_Spearman(j,i,n)=speedScore_tem(1);

      end
      Speed_KStest(j,i,ShuffleNum+1)=prctile(Speed_KStest(j,i,1:ShuffleNum),95);
      Speed_KStest(j,i,ShuffleNum+2)=prctile(Speed_KStest(j,i,1:ShuffleNum),5);
      [~,~,Speed_KStest(k,ShuffleNum+3)]=kstest2(KS_SpeedBinVector{k,1},SpeedTunning{k,1}.SpeedRange);

      Speed_Spearman(j,i,ShuffleNum+1)=prctile(Speed_Spearman(j,i,1:ShuffleNum),95);
      Speed_Spearman(j,i,ShuffleNum+2)=prctile(Speed_Spearman(j,i,1:ShuffleNum),5);
      speedScore_tem=SpatialTuning_BNT.speedScore(SpeedTrain_filtered,EventTrain_filtered(:,2),...
                                                  ExperimentInformation.Trackingframerate./double(ExperimentInformation.ImagingPlane));
      Speed_Spearman(j,i,ShuffleNum+3)=speedScore_tem(1);
       else
       end
    end
end

%% identify speed cells
IsPosSpeedCell=cell(ExperimentInformation.Session,1);
IsNagSpeedCell=cell(ExperimentInformation.Session,1);
IsNonlinearSpeedCell=cell(ExperimentInformation.Session,1);
for j=1:1:ExperimentInformation.Session
    g=1;
    k=1;
    h=1;
    for i=1:1:ExperimentInformation.TotalCell
            if Speed_KStest(j,i,Shuffling+3,j)>Speed_KStest(j,i,Shuffling+1,j)
                if Speed_Spearman(j,i,ShuffleNum+3)>Speed_Spearman(j,i,ShuffleNum+1)
                    IsPosSpeedCell{j,1}(g)=i;
                    g=g+1;
                elseif Speed_Spearman(j,i,ShuffleNum+3)<Speed_Spearman(j,i,ShuffleNum+2)
                    IsNagSpeedCell{j,1}(k)=i;
                    k=k+1;
                else
                    IsNonlinearSpeedCell{j,1}(h)=i;
                    h=h+1;
                end
            else
            end
    end
end

author: kavli-ntnu


text:
% Fig. 6.21   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%

clear all;
close all;

num=1;
den=[1 2 1];
rlocus(num,den);
grid;
title('Fig. 6.21: Root locus of G(s)=1/(s+2)^2 vs K');


text:
function vecCross = crossMat(vec)
%
% Computes the cross-product matrix of a 3x1 vector
%
    if(size(vec,1) ~= 3 || size(vec,2) ~= 1)
        error('Input vector must be 3x1');
    end

    vecCross = [ 0, 		-vec(3), 	vec(2);
				 vec(3),	0,			-vec(1);
				 -vec(2),	vec(1),		0        ];
end

repo: MSCKF


text:
function PrepareTrainData()

global param;

sceneFolder = param.trainingScenes;
outputFolder = param.trainingData;

[sceneNames, scenePaths, numScenes] = GetFolderContent(sceneFolder, [], true);

MakeDir(outputFolder);

for i = 1 : numScenes

    count = fprintf('Started working on scene %d of %d', i, numScenes);

    %%% reading input data
    curExpo = ReadExpoTimes(scenePaths{i});
    [curImgsLDR, curLabel] = ReadImages(scenePaths{i});

    %%% processing data
    [inputs, label] = ComputeTestExamples(curImgsLDR, curExpo, curLabel);

    %%% writing data
%     label=[];
%     inputs = ((curImgsLDR{1}));
%     inputs = cat(3, inputs, (curImgsLDR{2}));
%     inputs = cat(3, inputs, (curImgsLDR{3}));
    WriteTestExamples(inputs, label, [outputFolder, '/', sceneNames{i}, '.h5']);
%     WriteTestExamples(inputs, label, [outputFolder, '/', 'TrainSequence', num2str(i) '.h5']);

    fprintf(repmat('\b', [1, count]));
end

fprintf('Done\n\n');

author: qingsenyangit

repo: AHDRNet


text:
function midi=matrix2midi(M,ticks_per_quarter_note,timesig)
% midi=matrix2midi(M,ticks_per_quarter_note)
%
% generates a midi matlab structure from a matrix
%  specifying a list of notes.  The structure output
%  can then be used by writemidi.m
%
% M: input matrix:
%   1     2    3  4   5  6
%  [track chan nn vel t1 t2] (any more cols ignored...)
%
% optional arguments:
% - ticks_per_quarter_note: integer (default 300)
% - timesig: a vector of len 4 (default [4,2,24,8])
%

% Copyright (c) 2009 Ken Schutte
% more info at: http://www.kenschutte.com/midi

% TODO options:
%  - note-off vs vel=0
%  - tempo, ticks, etc

if nargin < 2
  ticks_per_quarter_note = 300;
end

if nargin < 3
  timesig = [4,2,24,8];
end

tracks = unique(M(:,1));
Ntracks = length(tracks);

% start building 'midi' struct

if (Ntracks==1)
  midi.format = 0;
else
  midi.format = 1;
end

midi.ticks_per_quarter_note = ticks_per_quarter_note;

tempo = 500000;   % could be set by user, etc...
% (microsec per quarter note)

for i=1:Ntracks

  trM = M(tracks(i)==M(:,1),:);

  note_events_onoff = [];
  note_events_n = [];
  note_events_ticktime = [];

  % gather all the notes:
  for j=1:size(trM,1)
    % note on event:
    note_events_onoff(end+1)    = 1;
    note_events_n(end+1)        = j;
    note_events_ticktime(end+1) = 1e6 * trM(j,5) * ticks_per_quarter_note / tempo;

    % note off event:
    note_events_onoff(end+1)    = 0;
    note_events_n(end+1)        = j;
    note_events_ticktime(end+1) = 1e6 * trM(j,6) * ticks_per_quarter_note / tempo;
  end

  msgCtr = 1;

  % set tempo...
  midi.track(i).messages(msgCtr).deltatime = 0;
  midi.track(i).messages(msgCtr).type = 81;
  midi.track(i).messages(msgCtr).midimeta = 0;
  midi.track(i).messages(msgCtr).data = encode_int(tempo,3);
  midi.track(i).messages(msgCtr).chan = [];
  msgCtr = msgCtr + 1;

  % set time sig...
  midi.track(i).messages(msgCtr).deltatime = 0;
  midi.track(i).messages(msgCtr).type = 88;
  midi.track(i).messages(msgCtr).midimeta = 0;
  midi.track(i).messages(msgCtr).data = timesig(:);
  midi.track(i).messages(msgCtr).chan = [];
  msgCtr = msgCtr + 1;

  [junk,ord] = sort(note_events_ticktime);

  prevtick = 0;
  for j=1:length(ord)

    n = note_events_n(ord(j));
    cumticks = note_events_ticktime(ord(j));

    midi.track(i).messages(msgCtr).deltatime = cumticks - prevtick;
    midi.track(i).messages(msgCtr).midimeta = 1;
    midi.track(i).messages(msgCtr).chan = trM(n,2);
    midi.track(i).messages(msgCtr).used_running_mode = 0;

    if (note_events_onoff(ord(j))==1)
      % note on:
      midi.track(i).messages(msgCtr).type = 144;
      midi.track(i).messages(msgCtr).data = [trM(n,3); trM(n,4)];
    else
      %-- note off msg:
      %midi.track(i).messages(msgCtr).type = 128;
      %midi.track(i).messages(msgCtr).data = [trM(n,3); trM(n,4)];
      %-- note on vel=0:
      midi.track(i).messages(msgCtr).type = 144;
      midi.track(i).messages(msgCtr).data = [trM(n,3); 0];
    end
    msgCtr = msgCtr + 1;

    prevtick = cumticks;
  end

  % end of track:
  midi.track(i).messages(msgCtr).deltatime = 0;
  midi.track(i).messages(msgCtr).type = 47;
  midi.track(i).messages(msgCtr).midimeta = 0;
  midi.track(i).messages(msgCtr).data = [];
  midi.track(i).messages(msgCtr).chan = [];
  msgCtr = msgCtr + 1;

end

% return a _column_ vector
% (copied from writemidi.m)
function A=encode_int(val,Nbytes)

A = zeros(Nbytes,1);  %ensure col vector (diff from writemidi.m...)
for i=1:Nbytes
  A(i) = bitand(bitshift(val, -8*(Nbytes-i)), 255);
end

author: kts

repo: matlab-midi

save_path: github-repos/MATLAB/kts-matlab-midi


% 0.05 700

% 0.1 1200

% 0.2 2200

%0.02

totalNum = 447026;

loopCompress(0.1, '/home/yinhuan/compression_exps/park/weightVector.txt', '/home/yinhuan/compression_exps/park/visMatrix/', 50, totalNum, 350, '/home/yinhuan/compression_exps/park/gurobi_compress_0.02/');

ratio = salientNumCnt( '/home/yinhuan/compression_exps/park/gurobi_compress_0.02/', totalNum ) / totalNum


text:
function pass = test_subsassgn(pref)

%% Setup.
if ( nargin == 0 )
    pref = chebfunpref();
end

%% See #2280

V = chebmatrix({1,2});

V(1) = pi;
pass(1) = V{1} == pi;

V{1} = -pi;
pass(2) = V{1} == -pi;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
M = 5;
S1 = 4;
S2 = 4;
S = S1 + S2;
K = 2;
Y = zeros(M, S);
Y(1, 1:4) = [1 2 3 4];
Y(2, 1:4) = [3 3 2 6];
Y(3, 5:8) = [3 3 6 6];
Y(4, 5:8) = [4 5 6 7];
Y
Y = spx.norm.normalize_l2(Y)
C = spx.fast.omp_spr(Y, K, 0)

Y * C
abs(Y - Y * C)

author: indigits

repo: sparse-plex


text:
function [mask,xmesh_new,ymesh_new] = dicomrt_trimZ(slice,voiref,inpmatrix,dose_xmesh,dose_ymesh,voi,voiselect)
% dicomrt_trimZ(slice,voiref,inpmatrix,dose_xmesh,dose_ymesh,voi,voiselect)
%
% Trim a matrix to min and max boundaries in the selected voi along X and Y.
%
% slice is the number of the section (relative to voiref) to be trimmed.
%   Slice is used to locate the matrix along Z and to match matrix with the corresponding contour in voiselect.
% inpmatrix is a 2D matrix.
% dose_xmesh, dose_ymesh, are the coordinates of the center of the pixel for eval and ref.
% voi and voiselect are the vois' cell array and the # of the voi to be used for the gamma calculation respectively.
%   They have to be specified together. Both matrices will be masked and reduced in size accordingly with the selected
%   voi's dimensions.
%
% NOTE: Use dicomrt_mask if you want to zeroes all the voxel data outside the voi of interest, without altering the
%       dimensions of the dose matrix.
%       This function leave unchanged matrix values outside the voi of interest: no mask is performed.
%       This functions changes the X, and Y dimensions of the input matrix to fit as close as possible
%       to the max and min contour position in x and y direction.
%
% Example:
%
% [trimdose,new_xmesh,new_ymesh]=dicomrt_trimZ(15,1,dose(:,:,15),dose_xmesh,dose_ymesh,VOI,4);
%
% calculates the position of the maximum and minimum contour for VOI number 4 in x and y and cut from "dose" all
% the voxels positioned outside these boundaries. Returns the new matrix in trimdose, and the new xmesh, ymesh
% in new_xmesh and new_ymesh. Here slice number 15 is defined in voiref 1.
%
% See also dicomrt_loadmcdose, dicomrt_dosediff
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check case and set-up some parameters and variables
[matrix_temp,type,label,PatientPosition]=dicomrt_checkinput(inpmatrix,1);
matrix=dicomrt_varfilter(matrix_temp);
[voi_temp,type,label]=dicomrt_checkinput(voi);
voi=dicomrt_varfilter(voi_temp);

[locate_voi_min_x,locate_voi_max_x,locate_voi_min_y,locate_voi_max_y] = ...
    dicomrt_voiboundariesZ(slice,voiref,dose_xmesh,dose_ymesh,voi_temp,voiselect,PatientPosition);

mask=matrix(locate_voi_min_y:locate_voi_max_y,locate_voi_min_x:locate_voi_max_x);

xmesh_new=dose_xmesh(locate_voi_min_x:locate_voi_max_x);
ymesh_new=dose_ymesh(locate_voi_min_y:locate_voi_max_y);

% Restore original variable format
[mask]=dicomrt_restorevarformat(matrix_temp,mask);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [TorF, vstr, rdate] = have_feature_intlinprog()
%HAVE_FEATURE_INTLINPROG  Detect availability/version info for INTLINPROG
%
%   Feature detection function implementing 'intlinprog' tag for HAVE_FEATURE
%   to detect availability/version of INTLINPROG, MILP solver from MATLAB
%   Optimization Toolbox 7.0 (R2014a) and later.
%
%   See also HAVE_FEATURE, MIQPS_MASTER, INTLINPROG.

%   MP-Opt-Model
%   Copyright (c) 2004-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

TorF = 0;
vstr = '';
rdate = '';
v = have_feature('optim', 'all');
if v.av && have_feature('matlab')   %% ignore Octave version
    TorF = exist('intlinprog', 'file') == 2;
    if TorF
        vstr = v.vstr;
        rdate = v.date;
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% Mathematics Q2595199
% https://math.stackexchange.com/questions/2595199
% Proximal Mapping of Least Squares with L1 and L2 Mixed Norm Regularization (Elastic Net)
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     13/03/2018
%   *   First release.

%% General Parameters

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

paramLam1 = 1;
paramLam2 = 2;

numElements = 8;

numIterations   = 1000;
stepSize        = 0.0075;

%% Generate Data

vB = 10 * randn([numElements, 1]);

hObjFun = @(vX) (0.5 * sum((vX - vB) .^ 2)) + (paramLam1 * norm(vX, 1)) + (paramLam2 * norm(vX, 2));
hL2SubGrad = @(vX) vX ./ max(norm(vX, 2), 1e-9);

hSoftThresholdL1 = @(vX, paramLambda) sign(vX) .* max(abs(vX) - paramLambda, 0);
hSoftThresholdL2 = @(vX, paramLambda) vX .* (1 - (paramLambda / (max(norm(vX, 2), paramLambda))));

%% Solution by CVX

cvx_begin('quiet')
    cvx_precision('best');
    variable vX(numElements)
    minimize( (0.5 * square_pos(norm(vX - vB, 2))) + (paramLam1 * norm(vX, 1)) + (paramLam2 * norm(vX, 2)) );
cvx_end

disp([' ']);
disp(['CVX Solution Summary']);
disp(['The CVX Solver Status - ', cvx_status]);
disp(['The Optimal Value Is Given By - ', num2str(cvx_optval)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Analytic Solution

vX = hSoftThresholdL2(hSoftThresholdL1(vB, paramLam1), paramLam2);
analyticObjVal = hObjFun(vX);

disp([' ']);
disp(['Analytic Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(analyticObjVal)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Solution by Sub Gradient Method

vObjValSgm = zeros([numIterations, 1]);

vX = zeros([numElements, 1]);
vObjValSgm(1) = hObjFun(vX);

for ii = 2:numIterations
    vG = (vX - vB) + (paramLam1 * sign(vX)) + (paramLam2 * hL2SubGrad(vX));
    vX = vX - (stepSize * vG);

    vObjValSgm(ii) = hObjFun(vX);
end

disp([' ']);
disp(['Sub Gradient Method Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(vObjValSgm(numIterations))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Solution by https://doi.org/10.1007/s10957-012-0245-9 (A Primal Dual Splitting Method for Convex Optimization Involving Lipschitzian, Proximable and Linear Composite Terms)
% Implementing Algorithm 3.2
% The Matrix L is identity
% The Gradient of F is vX - vB

vObjValSplit = zeros([numIterations, 1]);

paramSigma  = 0.05;
paramTau    = 0.05;
paramPhi    = 0.5;

hSoftThresholdL1 = @(vX, paramLambda) sign(vX) .* max(abs(vX) - paramLambda, 0);
hSoftThresholdL2 = @(vX, paramLambda) vX .* (1 - (paramLambda / (max(norm(vX, 2), paramLambda))));

hProxG = @(vX, paramTau) hSoftThresholdL1(vX, paramTau * paramLam1); %<! Soft Thresholding L1
hProxH = @(vX, paramSigma) hSoftThresholdL2(vX, paramSigma * paramLam2); %<! Soft Thresholding L2

hProxHConj = @(vX, paramSigma) vX - (paramSigma * hProxH(vX / paramSigma, (1 / paramSigma))); %<! Moreau’s Identity

vX = zeros([numElements, 1]);
vY = zeros([numElements, 1]);
vXX = vX;
vYY = vY;

vObjValSplit(1) = hObjFun(vX);

for ii = 2:numIterations
    vYY = hProxHConj(vY + (paramSigma * vX), paramSigma * paramLam2);
    vXX = hProxG(vX - (paramTau * (vX - vB + (2 * vYY) - vY)), paramTau);

    vX = (paramPhi * vXX) + (1 - paramPhi) * vX;
    vY = (paramPhi * vYY) + (1 - paramPhi) * vY;

    vObjValSplit(ii) = hObjFun(vX);
end

disp([' ']);
disp(['Split Method Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(vObjValSplit(numIterations))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Solution by ADMM (3 Blocks)
% Using Scaled Form (See Distributed Optimization and Statistical Learning
% via the Alternating Direction Method of Multipliers Pg. 15).

vObjValAdmm = zeros([numIterations, 1]);

paramRho = 0.25;

mA = [eye(numElements); eye(numElements)];
mB = [eye(numElements); zeros(numElements)];
mC = [zeros(numElements); eye(numElements)];

mI = eye(numElements);
mII = (1 / paramRho) * mI;
vBB = (1 / paramRho) * vB;

mAA = mA.' * mA;
mAB = mA.' * mB;
mAC = mA.' * mC;

mAAInv = inv(mAA + mII);

vX = zeros([numElements, 1]);
vY = zeros([numElements, 1]);
vZ = zeros([numElements, 1]);
vU = zeros([2 * numElements, 1]);

vObjValAdmm(1) = hObjFun(vX);

for ii = 2:numIterations

    vX = mAAInv * (vBB + (mAB * vY) + (mAC * vZ) - (mA.' * vU));
    % vX = (mAA + mII) \ (vBB + (mAB * vY) + (mAC * vZ) - (mA.' * vU));
    vY = SolveLsL1Prox(mB, ((mA * vX) + vU - (mC * vZ)), (paramLam1 / paramRho), numIterations);
    vZ = SolveLsL2Prox(mC, ((mA * vX) + vU - (mB * vY)), (paramLam2 / paramRho), numIterations);

    vU = vU + (mA * vX) - (mB * vY) - (mC * vZ);

    vObjValAdmm(ii) = hObjFun(vX);
end

disp([' ']);
disp(['ADMM 3 Blocks Method Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(vObjValAdmm(numIterations))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Display Reesults

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);
hAxes       = axes();
hLineSeries = plot([1:numIterations], [cvx_optval * ones([numIterations, 1]), analyticObjVal * ones([numIterations, 1]), vObjValSgm, vObjValSplit, vObjValAdmm]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(hLineSeries(2:end), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', ['Least Squares with Mixed Norm Regularization'], ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', 'Iteration Index', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', 'Objective Function Value', ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['CVX'], ['Analytic Solution'], ['Sub Gradient Method'], ['Split Method'], ['ADMM - 3 Blocks']});
set(hAxes, 'LooseInset', [0.07, 0.07, 0.07, 0.07]);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function integral = circle01_monomial_integral ( e )

%*****************************************************************************80
%
%% CIRCLE01_MONOMIAL_INTEGRAL: integrals on circumference of unit circle in 2D.
%
%  Discussion:
%
%    The integration region is
%
%      X^2 + Y^2 = 1.
%
%    The monomial is F(X,Y) = X^E(1) * Y^E(2).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Academic Press, 1984, page 263.
%
%  Parameters:
%
%    Input, integer E(2), the exponents of X and Y in the
%    monomial.  Each exponent must be nonnegative.
%
%    Output, real INTEGRAL, the integral.
%
  if ( any ( e(1:2) < 0 ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CIRCLE01_MONOMIAL_INTEGRAL - Fatal error!\n' );
    fprintf ( 1, '  All exponents must be nonnegative.\n' );
    error ( 'CIRCLE01_MONOMIAL_INTEGRAL - Fatal error!' );
  end

  if ( any ( mod ( e(1:2), 2 ) == 1 ) )

    integral = 0.0;

  else

    integral = 2.0;

    for i = 1 : 2
      integral = integral * gamma ( 0.5 * ( e(i) + 1 ) );
    end

    integral = integral / gamma ( 0.5 * sum ( e(1:2) + 1 ) );

  end

  return
end


text:
clear; close all; clc;

% vector field

[x,y]=meshgrid(-2:0.3:2,-2:.3:2);

u=x.^2;
v=y;

quiver(x,y,u,v)

grid on;
xlabel('x');
ylabel('y');
set(gca,'fontsize',15);

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
%DISTANCEMAPS_LEAF  One-line description here, please.
%
%   output = distanceMaps_leaf(input)
%
%   Example
%   distanceMaps_leaf
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
% Created: 2021-05-14,    using Matlab 9.8.0.1323502 (R2020a)
% Copyright 2021 INRAE.

% read data
img = imread('IJ-leaf-bin.tif');
imwrite(255 - img, 'ijLeaf-binInv.png', 'png');

% compute distance map
distMap = imDistanceMap(img);

% save RGB display
rgb = double2rgb(distMap, [1 1 1; parula]);
imwrite(rgb, 'ijLeaf-bin-distMap-rgb.png', 'png');

% compute thickness map
thMap = imThicknessMap(img);

% save RGB display
rgb = double2rgb(thMap, [1 1 1; parula]);
imwrite(rgb, 'ijLeaf-bin-thickMap-rgb.png', 'png');

author: mattools

repo: matImage


text:
function [planC] = dicomrt_dicomrt2cerr(plan,xmesh,ymesh,zmesh,ct,ct_xmesh,ct_ymesh,ct_zmesh,voi,optName,MRCell,PETCell,SPECTCell)
% dicomrt_dicomrt2CERR(dose,xmesh,ymesh,zmesh,ct,ct_xmesh,ct_ymesh,ct_zmesh,voi)
%
% Convert DICOMRT-TOOLBOX data in CERR format .
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)
%
% LM DK merging doses and beam information incase there are more than one

% Check number of argument and set-up some parameters and variables

error(nargchk(9,13,nargin))

%Create CERR data Container & templates for certain cells.
planInitC       = initializeCERR;
indexS          = planInitC{end};
headerInitS     = planInitC{indexS.header};
commentInitS    = planInitC{indexS.comment};
scanInitS       = planInitC{indexS.scan};
structureInitS  = planInitC{indexS.structures};
beamsInitS      = planInitC{indexS.beams};
doseInitS       = planInitC{indexS.dose};

tagMapS = initTagMapS;
tags = dicomrt_d2c_rtogtags;

% Initialize planC
planC = planInitC;
nCells = length(planInitC);
tags.nimages = 0;

% NOTE that conversion sequence is important and must be retained
%

% Convert header
disp('(+) Converting DICOM headers');
if isempty(plan)==0
    headerInitS=dicomrt_d2c_header(headerInitS,plan);
    disp('(+) Fetching data from RTPLAN');
elseif isempty(plan)==1 & isempty(ct)==0
    headerInitS=dicomrt_d2c_header(headerInitS,ct);
    disp('(+) Fetching data from CT scan');
elseif ~isempty(voi)
    headerInitS=dicomrt_d2c_header(headerInitS,voi);
    disp('(+) Fetching data from RTSTRUCT');
end
disp('(=) Conversion DICOM headers completed.');

% Convert comments
% commentInitS

% Convert scans
if isempty(ct)==0
    disp('(+) Converting DICOM scans');
    [tmp_scan,tags]=dicomrt_d2c_scan(scanInitS,indexS,ct,ct_xmesh,ct_ymesh,ct_zmesh,tags);
    disp('(=) Conversion DICOM scans completed.');
else
    disp('(=) CT scan not available.');
    tmp_scan = scanInitS;
end

% Validating and converting structures
% Validation
if isempty(voi)==0
    voitype=dicomrt_checkvoitype(voi);
    if isequal(voitype,'3D')
        voi=dicomrt_3dto2dVOI(voi);
    end

    if isempty(ct)==0
        disp('(+) Fitting DICOM structures to CT scans to comply with RTOG format');
        voi=dicomrt_fitvoi2ct(voi,ct_zmesh,ct);
        disp('(=) Fitting DICOM structures completed.');

        disp('(+) Validating DICOM structures');
        [voi]=dicomrt_validatevoi(voi,ct_xmesh,ct_ymesh,ct_zmesh);
        [newvoi]=dicomrt_closevoi(voi);
        disp('(=) Validation DICOM structures completed.');

        % Conversion
        disp('(+) Converting DICOM structures');
        %         names = fieldnames(structureInitS);
        %         if length(structureInitS) == 0 & ~isempty(names)
        %             structureInitS(1).(names{1}) = deal([]);
        %         end
        [tmp_str,tags]=dicomrt_d2c_voi(structureInitS,indexS,newvoi,ct_xmesh,ct_ymesh,ct_zmesh,tags,tmp_scan.scanUID);
        disp('(=) Conversion DICOM structures completed.');

    else
        % Conversion
        disp('(+) Converting DICOM structures');
        %         names = fieldnames(structureInitS);
        %         if length(structureInitS) == 0 & ~isempty(names)
        %             structureInitS(1).(names{1}) = deal([]);
        %         end
        [tmp_str,tags]=dicomrt_d2c_voi(structureInitS,indexS,voi,ct_xmesh,ct_ymesh,ct_zmesh,tags,tmp_scan.scanUID);
        disp('(=) Conversion DICOM structures completed.');
    end
else
    disp('(=) RTSTRUCT not available.');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DK
global dosetype;

if isempty(plan)==0
    if strcmpi(dosetype,'fraction')
        try
            beams=fieldnames(plan{1,1}{1}{1}.BeamSequence);
        catch
            warning('Beam Sequence not present in RTPLAN !');
        end
    else
        try
            try
                beams=fieldnames(plan{2,1}.BeamSequence);
            catch
                beams=fieldnames(plan{1}{1}.BeamSequence);
            end
        catch
            warning('Beam Sequence not present in RTPLAN !');
        end
    end
end

if exist('beams')==1
    % Convert beam geometry
    disp('(+) Converting DICOM beam geometry');
    %     names = fieldnames(beamsInitS);
    %     if length(beamsInitS) == 0 & ~isempty(names)
    %         beamsInitS(1).(names{1}) = deal([]);
    %     end
    if strcmpi(dosetype,'fraction')
        planSiz = length(plan);
        if planSiz == 1
            [tmp_beam,tags]=dicomrt_d2c_beam(beamsInitS,indexS,plan{1},tags);
        else
            for i = 1:planSiz(end)
                if i == 1
                [tmp_beam,tags]=dicomrt_d2c_beam(beamsInitS,indexS,plan{i},tags);
                else
                [tmp_beam1,tags]=dicomrt_d2c_beam(beamsInitS,indexS,plan{i},tags);
                tmp_beam = dissimilarInsert(tmp_beam,tmp_beam1); % insert new dose intp planC
                end
            end
        end
    else
        [tmp_beam,tags]=dicomrt_d2c_beam(beamsInitS,indexS,plan,tags);
        disp('(=) Conversion DICOM beam geometry completed.');
    end
else
    warning('Beam Sequence not present in RTPLAN !');
end

% Convert digital films
nfilms=0;
if isempty(plan)==0 & strcmpi(dosetype,'fraction')
    disp('(+) Converting DICOM dose fraction');
    planSiz = size(plan);
    for i = 1:planSiz(end)
        if i == 1
            [tmp_dose,tags]=dicomrt_d2c_dose(doseInitS,indexS,plan{i},xmesh{i},ymesh{i},zmesh{i},tags,tmp_scan.scanUID);
        else
            disp('(+) Converting DICOM dose fraction');
            [tmp_dose1,tags]=dicomrt_d2c_dose(doseInitS,indexS,plan{i},xmesh{i},ymesh{i},zmesh{i},tags,tmp_scan.scanUID);
            tmp_dose = dissimilarInsert(tmp_dose,tmp_dose1); % insert new dose intp planC
        end
    end
elseif isempty(plan)==0 & ~strcmpi(dosetype,'fraction') & isempty(plan{2,1})==0
    % Convert dose
    disp('(+) Converting DICOM dose');
    %     names = fieldnames(doseInitS);
    %     if length(doseInitS) == 0 & ~isempty(names)
    %         doseInitS(1).(names{1}) = deal([]);
    %     end
    [tmp_dose,tags]=dicomrt_d2c_dose(doseInitS,indexS,plan,xmesh,ymesh,zmesh,tags,tmp_scan.scanUID);
    disp('(=) Conversion DICOM dose completed.');
else
    disp('(=) RTDOSE not available');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end DK
% Convert DVH
% Convert seed geometry
% Convert MR
if ~isempty(MRCell)
    tmp_MR = dicomrt_d2c_MR(initializeMRInfo,MRCell);
end

% Convert PET
if ~isempty(PETCell)
    tmp_PET  = dicomrt_d2c_PET(initializePETInfo, PETCell);
end

if ~isempty(SPECTCell)
    tmp_SPECT = dicomrt_d2c_SPECT(initializeSPECTInfo, SPECTCell);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%Populate planC fields %%%%%%%%%%%%%%%%%%%%%%%%%%
% Store header info
planC{indexS.header} = headerInitS;

if isempty(plan)==0 & exist('beams')==1 & exist('tmp_dose')==1
    % Store beams
    planC{indexS.beams}=tmp_beam;
    % Store dose
    planC{indexS.dose}=tmp_dose;
elseif isempty(plan)==0 & exist('beams')==1 & exist('tmp_dose')==0
    % Store beams only
    planC{indexS.beams}=tmp_beam;
elseif isempty(plan)==0 & exist('beams')==0 &  exist('tmp_dose')==1
    % Store dose only
    planC{indexS.dose}=tmp_dose;
end

if isempty(ct)==0
    % Store scan
    planC{indexS.scan}=tmp_scan;
end

if isempty(voi)==0
    % Store structures
    planC{indexS.structures}=tmp_str;
end

% Store options and indexS
pathStr = getCERRPath;
optName = [pathStr 'CERROptions.m'];
optS = opts4Exe(optName);
planC{indexS.CERROptions} = optS;
planC{indexS.indexS} = indexS;

% Post processing actions

if isempty(ct)==0
    % Determine voxel thicknesses
    planC=dicomrt_d2c_setVoxelThicknesses(planC,indexS);
end

if isempty(ct)==0 & isempty(voi)==0
    % Get structure scan segments
    planC =  getRasterSegs(planC, optS);

    % Get any dose surface points
    planC =  getDSHPoints(planC,optS);
end

if isempty(ct)==0
    % Create uniformized datasets
    flagVOI = isempty(voi);
    planC=dicomrt_d2c_uniformizescan(planC,optS,flagVOI);
end

if ~isempty(MRCell)
    planC{indexS.scan}= dissimilarInsert(planC{indexS.scan},tmp_MR);
end

if ~isempty(PETCell)
    planC{indexS.scan} = dissimilarInsert(planC{indexS.scan},tmp_PET);
end

if ~isempty(SPECTCell)
    planC{indexS.scan} = dissimilarInsert(planC{indexS.scan},tmp_SPECT);
end

% Set UID's for Scan, Dose and Structures
planC = guessPlanUID(planC);


text:
function replay_scatter_abund(frames_scatter, frames_abund, fps)
%REPLAY_SCATTER_ABUND Summary of this function goes here
%   Detailed explanation goes here
if nargin < 3
    fps = 5;
end

concat = 0;
if length(frames_scatter) == length(frames_abund)
    scatter1 = frames_scatter(1);
    abund1 = frames_abund(1);
    if isempty(scatter1.colormap) && isempty(abund1.colormap)
        concat = 1;
    end
end

if ~concat
    scatter_fh = figure('name','Replay evolution of the Gaussians');
    movie(scatter_fh,frames_scatter,1,fps);
    abund_fh = figure('name','Replay evolution of the abundances');
    movie(abund_fh,frames_abund,1,fps);
    return;
end

N = length(frames_scatter);
disp(['Replay ',num2str(N), ' frames']);
for i = 1:N
    scatter_data = frames_scatter(i).cdata;
    abund_data = frames_abund(i).cdata;
    new_data = cat(2,scatter_data,abund_data);
    new_frames(i).cdata = new_data;
    new_frames(i).colormap = [];
end

fh = figure('name','Replay evolution of the Gaussians and the abundances');
pos = get(fh,'position');
pos(1) = pos(1) - pos(3)/2;
pos(3) = pos(3) * 2;
set(fh,'position',pos);
movie(fh,new_frames,1,fps);

author: zhouyuanzxcv

repo: Hyperspectral


text:
#!/usr/bin/octave -q
% Copyright 2018 Marc René Schädler
%
% This file is part of the mobile hearing aid prototype project
% The the mobile hearing aid prototype project 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.
%
% The mobile hearing aid prototype project 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 the mobile hearing aid prototype project. If not, see http://www.gnu.org/licenses/.

close all
clear
clc

fs = 48000;

[signal, fs] = audioread('/dev/shm/feedback.wav');
assert(fs==48000);

[feedback, range] = estimate_feedback(signal, fs);

feedback1 = feedback(:,1);
feedback2 = feedback(:,2);
range1 = range(1,:);
range2 = range(2,:);

% Write coefficients and ranges to files
fp = fopen('abhang/src/configuration/feedback1.bin','wb');
fwrite(fp,single(feedback1),'single');
fclose(fp);
fp = fopen('abhang/src/configuration/feedback2.bin','wb');
fwrite(fp,single(feedback2),'single');
fclose(fp);
fp = fopen('abhang/src/configuration/range1.bin','wb');
fwrite(fp,int32(range1),'int32');
fclose(fp);
fp = fopen('abhang/src/configuration/range2.bin','wb');
fwrite(fp,int32(range2),'int32');
fclose(fp);

[~, maxidx1] = max(abs(feedback1));
[~, maxidx2] = max(abs(feedback2));
latency1 = (maxidx1-1).*1000./fs;
latency2 = (maxidx2-1).*1000./fs;

printf("latency1 = %.2fms\nlatency2 = %.2fms\n",latency1,latency2);
printf("range1 = [%i %i]\nrange2 = [%i %i]\n",range1(1),range1(2),range2(1),range2(2));

author: m-r-s

repo: hearingaid-prototype