text:
function varargout = lambda_min(varargin)
% LAMBDA_MIN Returns largest eigenvalue of Hermitian matrix.
%
% r = LAMBDA_MIN(X)
%
% See also LAMBDA_MAX

switch class(varargin{1})

    case 'double' % What is the numerical value of this argument (needed for displays etc)
        error(nargchk(1,1,nargin));
        X = varargin{1};
        [n,m] = size(X);
        if ~ishermitian(X)
            error('LAMBDA_MIN can only be applied on Hermitian matrices');
        else
            varargout{1} = min(real(eig(X)));
        end

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.
        error(nargchk(1,1,nargin));
        X = varargin{1};
        [n,m] = size(X);
        if ~ishermitian(X)
            error('LAMBDA_MIN can only be Hermitian matrices');
        else
            varargout{1} = yalmip('define',mfilename,varargin{:});
        end

    case 'char' % YALMIP sends 'model' when it wants the epigraph or hypograph
        if isequal(varargin{1},'graph')
            t = varargin{2}; % Second arg is the extended operator variable
            X = varargin{3}; % Third arg and above are the args user used when defining t.
            varargout{1} = (t*eye(size(X,1)) <= X);
            varargout{2} = struct('convexity','concave','monotonicity','none','definiteness','none','model','graph');
            varargout{3} = X;
        else
            varargout{1} = [];
            varargout{2} = [];
            varargout{3} = [];
        end
    otherwise
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function plaza = clear_boundary(plaza)

[a,b] = size(plaza);
for i = 1:b
    if plaza(a,i) > 0
        plaza(a,i) = 0;
    end
end

author: Eurus-Holmes


text:
%Copyright 2013 The MathWorks, Inc
function [sCol,sRx,sRD,sMFilt,sTVG,threshold] = setupRx(nint,nf,pfa,maxrange,range_gates,sWav,sAnt,fc)
%% Collector
sCol = phased.Collector('Sensor',sAnt,'Wavefront','Plane',...
                        'OperatingFrequency',fc);

%% Receiver Preamp
sRx = phased.ReceiverPreamp('Gain',20,'NoiseBandwidth',sWav.SweepBandwidth,...
                            'NoiseFigure',nf,'EnableInputPort',true,...
                            'SeedSource','Property','Seed',2007);

%% Range Doppler Estimator
sRD = phased.RangeDopplerResponse(...
    'DopplerWindow','Chebyshev','DopplerSidelobeAttenuation',60,...
    'SampleRate',sWav.SampleRate,'DopplerOutput','Speed','OperatingFrequency',fc);

%% Matched Filter
match_sig = getMatchedFilter(sWav);
sMFilt    = phased.MatchedFilter('Coefficients',match_sig);

%% Time Varying Gain
lambda   = physconst('LightSpeed')/fc;
rng_loss = 2*fspl(range_gates,lambda);   % Factor 2 for round trip
ref_loss = 2*fspl(maxrange,lambda);      % Reference loss for maximum range

sTVG = phased.TimeVaryingGain('RangeLoss',rng_loss,'ReferenceLoss',ref_loss);

%% Detection Threshold
npower    = noisepow(sWav.SweepBandwidth,nf,sRx.ReferenceTemperature);
sG        = phased.ArrayGain('SensorArray',sAnt);
gain      = step(sG,fc,0);
threshold = npower * norm(match_sig)^2 * db2pow(gain)...
                   * db2pow(npwgnthresh(pfa,nint,'noncoherent'));

author: Sable

repo: mcbench-benchmarks


text:
function varargout = singleSignTest(varargin)
%SINGLESIGNTEST   Heuristic check to see if a DISKFUN changes sign.
%   SINGLESIGNTEST(F) returns 1 if the values of F on a tensor grid are of
%   the same sign.
%
%   [OUT, WZERO] = SINGLESIGNTEST( F ), returns WZERO = 1 if a zero has
%   been found.
%
%   The algorithm works by sampling F on a tensor-grid and checking if
%   those values are of the same sign. This command is mainly for internal
%   use in DISKFUN commands.
%
% See also DISKFUN/ABS, DISKFUN/SQRT, and DISKFUN/LOG.

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

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

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function failed_ft_headmodel_bemcp

% MEM 12gb
% WALLTIME 03:00:00
% DEPENDENCY ft_headmodel_localspheres ft_prepare_localspheres

% to test actual numerical output of bemcp, rather than simply successful
% running and correct inputs (which is what test_ft_prepare_headmodel tests).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% method 1: segment MRI then compute bnd from it

% read in the mri, this is already (non?)linearly aligned with MNI
mri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/template/headmodel/standard_mri.mat'));

cfg = [];
cfg.output = {'brain', 'skull', 'scalp'};
segmentedmri_new = ft_volumesegment(cfg, mri);

cfg = [];
cfg.tissue = {'brain', 'skull', 'scalp'};
cfg.numvertices = [1500 1000 500]; % these numbers later match the precomputed bnd
mesh_new = ft_prepare_mesh(cfg, segmentedmri_new);

% I'm choosing 'mm' as units since that is what both the MRI and the
% standard_bem are in already. Whether that is optimal for 'bemcp' is not
% clear to me.

% Create simple concentricspheres to compare to BEMCP
cfg = [];
cfg.method = 'concentricspheres';
vol_new_cs = ft_prepare_headmodel(cfg, ft_convert_units(mesh_new, 'mm'));

cfg = [];
cfg.method = 'bemcp';
vol_new_bemcp = ft_prepare_headmodel(cfg, ft_convert_units(mesh_new, 'mm'));

try
  cfg = [];
  cfg.method = 'dipoli';
  vol_new_dipoli = ft_prepare_headmodel(cfg, ft_convert_units(mesh_new, 'mm'));
catch
  vol_new_dipoli = [];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% method 2: compute bnd from segmented MRI

% read in the already-segmented mri (in case its segmentation is better than above)
segmentedmri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/template/headmodel/standard_seg.mat'));

cfg = [];
cfg.tissue = [3 2 1]; % brain, skull, scalp
cfg.numvertices = [3000 2000 1000];
mesh_seg = ft_prepare_mesh(cfg, segmentedmri);

figure;
ft_plot_mesh(mesh_new(1), 'facecolor', [0.2 0.2 0.2], 'facealpha', 0.3, 'edgecolor', [1 1 1], 'edgealpha', 0.05);
hold on;
ft_plot_mesh(mesh_seg(3), 'facecolor', [0.4 0.6 0.4], 'facealpha', 0.3, 'edgecolor', 'none');
% aside from possible different order, they should align spatially

% Create simple concentricspheres to compare to BEMCP
cfg = [];
cfg.method = 'concentricspheres';
vol_seg_cs = ft_prepare_headmodel(cfg, ft_convert_units(mesh_seg, 'mm'));

cfg = [];
cfg.method = 'bemcp';
vol_seg_bemcp = ft_prepare_headmodel(cfg, ft_convert_units(mesh_seg, 'mm'));

try
  cfg = [];
  cfg.method = 'dipoli';
  vol_seg_dipoli = ft_prepare_headmodel(cfg, ft_convert_units(mesh_new, 'mm'));
catch
  vol_seg_dipoli = [];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% method 3: use existing bnd in standard_bem

% presumably the mesh that was used to create this 'standard' output from
% dipoli should work (irrespective of segmentation above).
load(dccnpath('/home/common/matlab/fieldtrip/template/headmodel/standard_bem.mat'));
vol_exist_dipoli = vol;
mesh_exist = vol.bnd;
mesh_exist = ft_convert_units(mesh_exist, 'mm');
clear vol

figure;
ft_plot_mesh(mesh_new(3), 'facecolor', [0.2 0.2 0.2], 'facealpha', 0.3, 'edgecolor', [1 1 1], 'edgealpha', 0.05);
hold on;
ft_plot_mesh(mesh_exist(3), 'edgecolor', 'none', 'facecolor', [0.4 0.6 0.4]);
% aside from possible different order, they should align spatially

% Create simple concentricspheres to compare to BEMCP
cfg = [];
cfg.method = 'concentricspheres';
vol_exist_cs = ft_prepare_headmodel(cfg, mesh_exist);

cfg = [];
cfg.method = 'bemcp';
vol_exist_bemcp = ft_prepare_headmodel(cfg, mesh_exist);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Check coregistration and view results

load standard_sourcemodel3d8mm;
sourcemodel = ft_convert_units(sourcemodel, 'mm');

elec = ft_read_sens(dccnpath('/home/common/matlab/fieldtrip/template/electrode/standard_1005.elc'));
elec = ft_convert_units(elec, 'mm');

% check for general spatial alignment
figure; ft_plot_mesh(vol_new_bemcp.bnd(3), 'facecolor', [0.2 0.2 0.2], 'facealpha', 0.3, 'edgecolor', [1 1 1], 'edgealpha', 0.05);
hold on; ft_plot_sens(elec);
hold on; ft_plot_mesh(sourcemodel.pos(sourcemodel.inside, :));

figure; ft_plot_mesh(vol_seg_bemcp.bnd(3), 'facecolor', [0.2 0.2 0.2], 'facealpha', 0.3, 'edgecolor', [1 1 1], 'edgealpha', 0.05);
hold on; ft_plot_sens(elec);
hold on; ft_plot_mesh(sourcemodel.pos(sourcemodel.inside, :));

figure; ft_plot_mesh(vol_exist_bemcp.bnd(3), 'facecolor', [0.2 0.2 0.2], 'facealpha', 0.3, 'edgecolor', [1 1 1], 'edgealpha', 0.05);
hold on; ft_plot_sens(elec);
hold on; ft_plot_mesh(sourcemodel.pos(sourcemodel.inside, :));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% View leadfields

% make sure vol*.mat of bemcp not NaN
% see bug 1954
if any(isnan(vol_new_bemcp.mat(:)))
  warning('NaN in vol_new_bemcp.mat')
end

close all
clear volnames vol
volnames = whos('vol*');

% these would need to be recomputed given the headmodel geometry
% sourcemodel = rmfield(sourcemodel, 'inside');
% sourcemodel = rmfield(sourcemodel, 'outside');

% or even faster
sourcemodel = [];
sourcemodel.pos = [-20 0 50];
sourcemodel.unit = 'mm';

for ll = 1:length(volnames)
  vol = eval(volnames(ll).name);

  if ~isempty(vol)
    % vol can be empty in case it failed to be created above, e.g. one of
    % the dipoli examples, for one reason or another
    cfg = [];
    cfg.sourcemodel = sourcemodel;
    cfg.elec = elec;
    cfg.headmodel = vol;

    lf = ft_prepare_leadfield(cfg);

    % Hack: to make LF appear as if it were an ERP
    tlock = [];
    tlock.time = [1 2 3];
    tlock.avg = lf.leadfield{dsearchn(sourcemodel.pos, [-20 0 50])};
    tlock.label = lf.cfg.channel;
    tlock.dimord = 'chan_time';

    cfg = [];
    cfg.layout = 'elec1010.lay';
    cfg.xlim = [0.9 1.1];
    figure;
    ft_topoplotER(cfg, tlock);
    title(volnames(ll).name)
    cfg.xlim = [1.9 2.1];
    figure;
    ft_topoplotER(cfg, tlock);
    title(volnames(ll).name)
    cfg.xlim = [2.9 3.1];
    figure;
    ft_topoplotER(cfg, tlock);
    title(volnames(ll).name)

    disp(volnames(ll).name)
  end
end
% For me, sensible patterns come from the vol*cs and vol*dipoli, but not from the vol*bemcp

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = Sim_MonteCarlo_Diffusion_GUI(varargin)
% ## Purpose
% This Monte Carlo simulator for 2D diffusion is able to generate synthetic diffusion signal from any 2D axon packing.
%
% ## Approach
% First, the [axonpacking](https://github.com/neuropoly/axonpacking) software is used to generate dense packing of circles.
% Then, the Monte Carlo simulator moves water molecules around in the packing and apply the phase shift induced by the diffusion gradients.
%
% ## Interface
% First step is to generate an axon packing. Pre-packed examples are
% proposed.
%
% Then, when clicking on `Update` button, the Monte Carlo simulation runs in this packing.
% _Note: blues dots correspond to extra-axonal water molecules and red dots to intra-axonal water_
% _Note: The `MPG strength` and `Signal Strength` plots correspond to the last bvalue in the protocol. But Signal Strength is computed for all bvalues from the protocol._
%
% Once the simulation is finished, synthetic data are plotted and fitted by the charmed model.
% The user is able to save this synthetic signal in a `.mat` file.
%
% Authors: Yasuhiko Tachibana, Tanguy Duval, Tom Mingasson, 2019

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

% Last Modified by GUIDE v2.5 26-Aug-2019 15:37:40

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

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

% --- Executes just before Sim_SimMCdiff is made visible.
function Sim_SimMCdiff_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = hObject;
handles.Model = varargin{1};
setappdata(0,'Model',handles.Model);

if ~isfield(handles,'opened')
    if ~ismac
        set(findobj(hObject,'Type','uicontrol'),'FontSize',7);
    end % everything is bigger on windows or linux

    savedPacks = dir(fullfile(fileparts(which('Sim_MonteCarlo_Diffusion_GUI.m')),'savedPacks','*.mat'));
    set(handles.preset_packing,'String',{savedPacks.name});
    [axons, packing] = loadPreset(handles);
    handles.axonpacking.axons = axons;
    handles.axonpacking.packing = packing;

    % Statistics from the packing
    k=1;
    [FVF, FR, MVF, AVF] = compute_statistics(axons.d{k}/2, axons.Delta{k}, packing.final_positions{k}, [], axons.g_ratio{k});
    set(handles.tableVolumes,'Data',[FVF; FR; MVF; AVF])

    % set values
    slider_Naxons_Callback(handles.slider_Naxons, [], handles)
    slider_dvar_Callback(handles.slider_dvar, [], handles)
    slider_dmean_Callback(handles.slider_dmean, [], handles)
    slider_gap_Callback(handles.slider_gap, [], handles)
    slider_trans_Callback(handles.slider_trans, [], handles)
    slider_numelparticle_Callback(handles.slider_numelparticle, [], handles)
    slider_Dcoef_Callback(handles.slider_Dcoef, [], handles)

    % clear axe
    handles.opened = 1;
end
% Update handles structure
guidata(hObject, handles);

% UIWAIT makes Sim_SimMCdiff wait for user response (see UIRESUME)
% uiwait(handles.Simu);

% --- Outputs from this function are returned to the command line.
function varargout = Sim_SimMCdiff_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure

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

% --- Executes on button press in SimMCdiffUpdate.
function SimMCdiffUpdate_Callback(hObject, eventdata, handles)
set(findobj('Name','SimMCdiff'),'pointer', 'watch'); drawnow;
MonteCarloSim(handles, handles.axonpacking.axons, handles.axonpacking.packing)
set(findobj('Name','SimMCdiff'),'pointer', 'arrow'); drawnow;

% GETAPPDATA
function varargout = GetAppData(varargin)
for k=1:nargin; varargout{k} = getappdata(0, varargin{k}); end

%SETAPPDATA
function SetAppData(varargin)
for k=1:nargin; setappdata(0, inputname(k), varargin{k}); end

% RMAPPDATA
function RmAppData(varargin)
for k=1:nargin; rmappdata(0, varargin{k}); end

% --------------------------------------------------------------------
function helpbutton_ClickedCallback(hObject, eventdata, handles)
doc Sim_MonteCarlo_Diffusion_GUI

function slider_gap_Callback(hObject, eventdata, handles)
%update_axonSetup(handles);
set(handles.text_gap,'String', ['Gap between axons: ' num2str(get(hObject,'Value')) 'um'])
function slider_dvar_Callback(hObject, eventdata, handles)
update_axonSetup(handles);
set(handles.text_dvar,'String', ['Diameter variance: ' num2str(get(hObject,'Value')) 'um'])
function slider_dmean_Callback(hObject, eventdata, handles)
update_axonSetup(handles);
set(handles.text_dmean,'String', ['mean diameter: ' num2str(get(hObject,'Value')) 'um'])
function slider_Naxons_Callback(hObject, eventdata, handles)
update_axonSetup(handles);
set(handles.text_Naxons,'String', ['# axons: ' num2str(round(get(hObject,'Value')))])

function [d, x0, side, axons] = update_axonSetup(handles)
k=1;
axons.N{k}      = round(get(handles.slider_Naxons,'Value'));
axons.d_mean{k} = round(get(handles.slider_dmean,'Value')*10)/10;
set(handles.slider_dmean,'Value',axons.d_mean{k})
axons.d_var{k}  = round(get(handles.slider_dvar,'Value')*10)/10;
set(handles.slider_dvar,'Value',axons.d_var{k});
axons.Delta{k}  = get(handles.slider_gap,'Value');
axons.threshold_high{k} = 20;
axons.threshold_low{k}  = .1;
axons
[d, x0, side] = axons_setup(axons,'gamma', k, handles.axes_axonDist);
if ~ismac, set(handles.axes_axonDist,'FontSize',7); end

function run_pack_Callback(hObject, eventdata, handles)
k=1;
[d, x0, side, axons] = update_axonSetup(handles);

axons.d{k} = d;
axons.g_ratio{k} = compute_gratio(d);

% packing process of the axons
iter_max = 10000;
iter_fvf = iter_max/10;
[final_positions, final_overlap, fvf_historic] = process_packing(x0, d/2, axons.Delta{k}, side, iter_max, iter_fvf);

% store packing results
% main results
packing.initial_positions{k}    = reshape(x0,2,length(x0)/2);
packing.final_positions{k}      = final_positions;
% secondary results
packing.final_overlap{k}        = final_overlap;
packing.FVF_historic{k}         = fvf_historic;
packing.iter_max{k}             = iter_max;

% Statistics from the packing
[FVF, FR, MVF, AVF] = compute_statistics(axons.d{k}/2, axons.Delta{k}, packing.final_positions{k}, side, axons.g_ratio{k});
set(handles.tableVolumes,'Data',[FVF; FR; MVF; AVF])

if ishandle(201), close(201); end
plotPacking(handles,axons,packing)
handles.axonpacking.axons = axons;
handles.axonpacking.packing = packing;

guidata(hObject, handles);

function save_pack_Callback(hObject, eventdata, handles)
axons   = handles.axonpacking.axons;
packing = handles.axonpacking.packing;
fname = sprintf('pack_d%.1fvar%.1fgap%.1f.mat',axons.d_mean{1},axons.d_var{1},axons.Delta{1});
fname = fullfile(fileparts(which('Sim_MonteCarlo_Diffusion_GUI.m')),'savedPacks',fname);
[fname, path] = uiputfile('*.mat','Save Packing as...',fname);
if fname
    save(fullfile(path,fname),'axons','packing')
    savedPacks = get(handles.preset_packing,'String');
    set(handles.preset_packing,'String',unique([savedPacks; {fname}]));
    savedPacks = get(handles.preset_packing,'String');
    set(handles.preset_packing,'Value',length(savedPacks));
end

function preset_packing_Callback(hObject, eventdata, handles)
[axons, packing] = loadPreset(handles);
handles.axonpacking.axons = axons;
handles.axonpacking.packing = packing;
k=1;
[FVF, FR, MVF, AVF] = compute_statistics(axons.d{k}/2, axons.Delta{k}, packing.final_positions{k}, [], axons.g_ratio{k});
set(handles.tableVolumes,'Data',[FVF; FR; MVF; AVF])

guidata(hObject, handles);

function     [axons, packing] = loadPreset(handles)
file = get(handles.preset_packing,'String');
load(file{get(handles.preset_packing,'Value')})
axons_setup(axons,'gamma', 1, handles.axes_axonDist);
if ~ismac, set(handles.axes_axonDist,'FontSize',7); end
plotPacking(handles,axons,packing)

function plotPacking(handles,axons,packing)
% plot disks
numelobj = size(packing.final_positions{1},2);
cen = packing.final_positions{1}'; cen=cen(:,[2 1]); cen=cen-repmat([mean(cen(:,1)) mean(cen(:,2))],[size(cen,1) 1]); % center position of the cells
cen = [cen, zeros(numelobj,1)];

% Radius of the "cells"
R = [axons.d{1}.*axons.g_ratio{1} axons.d{1}]/2; % internal / external diameter in um
R = R(:,2); %use external radius only.

t = linspace(0,2*pi);

axes(handles.axes_axonPack)
hold off
for k =1:size(R,1)
    plot(R(k,1)*cos(t)+cen(k,1), R(k)*sin(t)+cen(k,2),'b','linewidth',1);
    hold on
end
hold off
axis equal tight
xlabel x(\mum)
ylabel y(\mum)

if ~ismac, set(handles.axes_axonPack,'FontSize',7); end

function MonteCarloSim(handles, axons, packing)

% Read updated Model
Model_new = getappdata(0,'Model');
if ~isempty(Model_new) && strcmp(class(Model_new),class(handles.Model))
    handles.Model = Model_new;
end

axes(handles.uipanel12);

% Read parameters
numelparticle = round(get(handles.slider_numelparticle,'Value'));

trans_mean = get(handles.slider_trans,'Value'); % mean probability of penetrating the cell walls [0-1]

D = get(handles.slider_Dcoef,'Value'); % Diffusion coefficient

handles.Model.Sim_MonteCarlo_Diffusion(numelparticle, trans_mean, D, packing, axons);

function slider_trans_Callback(hObject, eventdata, handles)
set(handles.text_trans,'String', ['Permeability: ' num2str(get(hObject,'Value'))])
function slider_numelparticle_Callback(hObject, eventdata, handles)
set(hObject,'Value',round(get(hObject,'Value')))
set(handles.text_numelparticle,'String', ['Number of particles: ' num2str(get(hObject,'Value'))])
function slider_bv_Callback(hObject, eventdata, handles)
set(handles.text_bv,'String', ['bvalue: ' num2str(get(hObject,'Value')) 'sec/mm2'])
function slider_Dcoef_Callback(hObject, eventdata, handles)
D = get(hObject,'Value');
set(handles.text_Dcoef,'String', sprintf('Diffusion coefficient: \n\tD = %.2gx 10-3 mm2/sec',D))
set(handles.text_Step,'String', sprintf('steptime = %.1g [ms]\nstepflight = %.1g [um]\n(stepflight^2 = 4*D*steptime)',0.5,sqrt(4*D*.5)))

function preset_packing_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function slider_gap_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end
function slider_dvar_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end
function slider_dmean_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end
function slider_Naxons_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end
function slider_trans_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end
function slider_numelparticle_CreateFcn(hObject, eventdata, handles)
function slider_bv_CreateFcn(hObject, eventdata, handles)
function slider9_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end
function slider_Dcoef_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function xdot=Duffing(t,x)
global Gamma;
xdot(1)=x(2);
xdot(2)=x(1)-0.1*x(2)-(x(1))^3+Gamma*cos(1.25*t);
xdot=[xdot(1);xdot(2)];


text:
function [orig, data] = load_curry_data_file(datafile)

[PathName, FileName, Extension] = fileparts(datafile);

if isempty(PathName)
  % the dataset is in the current working directory, refer to this with a .
  % to ensure proper path detection later on
  PathName = '.';
end

FileName = [FileName, Extension];
PathName = [PathName filesep]; % add the / or \

% loads Curry dat/dap/rs3 files into MATLAB and displays their waveforms
%{
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ask user to specify file
Title = 'Open Curry Data File';

[FileName,PathName,FilterIndex] = uigetfile({'*.cdt;*.dat' 'All Curry Data Files';'*.cdt' 'Curry 8 Data Files'; '*.dat' 'Curry 7 Data Files'},Title,pwd);

if ( FilterIndex == 0 )
    return;
end

if ( FileName == 0 )
    msgbox('Invalid File Name',Title,'error');
    return;
end
%}

Length = size(FileName,2);
DataFile = [PathName,FileName];
Extension = FileName(Length-2:Length);

if ( strcmpi ( Extension,'dat' ) )
    BaseName = [PathName,FileName(1:Length-4)];
    ParameterFile = [BaseName,'.dap'];
    LabelFile = [BaseName,'.rs3'];
    EventFile = [BaseName,'.cef'];
    EventFile2 = [BaseName,'.ceo'];
elseif ( strcmpi ( Extension,'cdt' ) )
    ParameterFile = [PathName,FileName,'.dpa'];
    LabelFile = [PathName,FileName,'.dpa'];
    EventFile = [PathName,FileName,'.cef'];
    EventFile2 = [PathName,FileName,'.ceo'];
else
    %msgbox('Unsupported File Name (choose a .cdt or .dat file)',Title,'error');
    error('Unsupported File Name (choose a .cdt or .dat file)');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% open parameter file
fid = fopen_or_error(ParameterFile,'rt');
cell = textscan(fid,'%s','whitespace','','endofline',char(167));
fclose(fid);
cont = cell2mat(cell{1});

% check for compressed file format
ctok = 'DataGuid';
ix = strfind(cont,ctok);
if ~isempty ( ix )
    text = sscanf(cont(ix+numel(ctok):end),' = %s');
    if (strcmp(text, '{2912E8D8-F5C8-4E25-A8E7-A1385967DA09}') == 1)
        msgbox('Unsupported data format (compressed). Use Curry to convert this file to Raw Float format.',Title, 'error')
        return;
    end
end

% read parameters from parameter file
% tokens (second line is for Curry 6 notation)
tok = { 'NumSamples'; 'NumChannels'; 'NumTrials'; 'SampleFreqHz';  'TriggerOffsetUsec';  'DataFormat'; 'DataSampOrder';   'SampleTimeUsec';
        'NUM_SAMPLES';'NUM_CHANNELS';'NUM_TRIALS';'SAMPLE_FREQ_HZ';'TRIGGER_OFFSET_USEC';'DATA_FORMAT';'DATA_SAMP_ORDER'; 'SAMPLE_TIME_USEC'};

% scan in cell 1 for keywords - all keywords must exist!
nt = size(tok,1);
a = zeros(nt,1);
for i = 1:nt
     ctok = tok{i,1};
     ix = strfind(cont,ctok);
     if ~isempty ( ix )
         text = sscanf(cont(ix+numel(ctok):end),' = %s');     % skip =
         if strcmp ( text,'ASCII' ) || strcmp ( text,'CHAN' ) % test for alphanumeric values
             a(i) = 1;
         else
             c = sscanf(text,'%f');         % try to read a number
             if ~isempty ( c )
                 a(i) = c;                  % assign if it was a number
             end
         end
     end
end

% derived variables. numbers (1) (2) etc are the token numbers
nSamples    = a(1)+a(1+nt/2);
nChannels   = a(2)+a(2+nt/2);
nTrials     = a(3)+a(3+nt/2);
fFrequency  = a(4)+a(4+nt/2);
fOffsetUsec = a(5)+a(5+nt/2);
nASCII      = a(6)+a(6+nt/2);
nMultiplex  = a(7)+a(7+nt/2);
fSampleTime = a(8)+a(8+nt/2);

if ( fFrequency == 0 && fSampleTime ~= 0 )
    fFrequency = 1000000 / fSampleTime;
end

% try to guess number of samples based on datafile size
if nSamples < 0
    if nASCII == 1
        msgbox('Number of samples cannot be guessed from ASCII data file. Use Curry to convert this file to Raw Float format.',Title,'error');
        return;
    else
        fileInfo = dir(DataFile);
        fileSize = fileInfo.bytes;
        nSamples = fileSize / (4 * nChannels * nTrials);
    end
end

%Search for Impedance Values
tixstar = strfind(cont,'IMPEDANCE_VALUES START_LIST');
tixstop = strfind(cont,'IMPEDANCE_VALUES END_LIST');

impedancelist = [];
impedancematrix = [];

if (~isempty(tixstar)) && (~isempty(tixstop))
    text = cont(tixstar:tixstop-1);
    tcell = textscan(text,'%s');
    tcell = tcell{1,1};
    for tcC = 1:size(tcell,1)
       tcell{tcC} = str2num(tcell{tcC}); % data was read in as strings - force to numbers
       if ~isempty(tcell{tcC}) % skip if it is not a number
           impedancelist(end+1) = tcell{tcC};
       end
    end

    % Curry records last 10 impedances
    impedancematrix = reshape(impedancelist,[(size(impedancelist,2)/10),10])';
    impedancematrix(impedancematrix == -1) = NaN; % screen for missing
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% open label file
fid = fopen_or_error(LabelFile,'rt');
cell = textscan(fid,'%s','whitespace','','endofline',char(167));
fclose(fid);
cont = cell2mat(cell{1});

% read labels from label file
% initialize labels
labels = num2cell(1:nChannels);

for i = 1:nChannels
    text = sprintf('EEG%d',i);
    labels(i) = cellstr(text);
end

% scan in cell 1 for LABELS (occurs four times per channel group)
ix = strfind(cont,[char(10),'LABELS']);
nt = size(ix,2);
nc = 0;

for i = 4:4:nt                                                      % loop over channel groups
    newlines = ix(i-1) + strfind(cont(ix(i-1)+1:ix(i)),char(10));   % newline
    last = nChannels - nc;
    for j = 1:min(last,size(newlines,2)-1)                          % loop over labels
        text = cont(newlines(j)+1:newlines(j+1)-1);
        if isempty(strfind(text,'END_LIST'))
            nc = nc + 1;
            labels(nc) = cellstr(text);
        else
            break
        end
    end
end

% read sensor locations from label file
% initialize sensor locations
sensorpos = zeros(3,0);

% scan in cell 1 for SENSORS (occurs four times per channel group)
ix = strfind(cont,[char(10),'SENSORS']);
nt = size(ix,2);
nc = 0;

for i = 4:4:nt                                                      % loop over channel groups
    newlines = ix(i-1) + strfind(cont(ix(i-1)+1:ix(i)),char(10));   % newline
    last = nChannels - nc;
    for j = 1:min(last,size(newlines,2)-1)                          % loop over labels
        text = cont(newlines(j)+1:newlines(j+1)-1);
        if isempty(strfind(text,'END_LIST'))
            nc = nc + 1;
            tcell = textscan(text,'%f');
            posx = tcell{1}(1);
            posy = tcell{1}(2);
            posz = tcell{1}(3);
            sensorpos = cat ( 2, sensorpos, [ posx; posy; posz ] );
        else
            break
        end
    end
end

%Search for Epoch Labels
tixstar = strfind(cont,'EPOCH_LABELS START_LIST');
tixstop = strfind(cont,'EPOCH_LABELS END_LIST');

epochlabelslist = [];

if (~isempty(tixstar)) && (~isempty(tixstop))
    text = cont(tixstar:tixstop-1);
    tcell = textscan(text,'%s', 'delimiter','\n','whitespace','', 'headerlines', 1);
    epochlabelslist = tcell{1,1};
end

%Search for Epoch Information
tixstar = strfind(cont,'EPOCH_INFORMATION START_LIST');
tixstop = strfind(cont,'EPOCH_INFORMATION END_LIST');

epochinformationlist = [];

if (~isempty(tixstar)) && (~isempty(tixstop))
    text = cont(tixstar:tixstop-1);
    tcell = textscan(text,'%d%d%d%d%d%d%d', 'delimiter','\n','headerlines', 1);
    epochinformationlist = cell2mat(tcell);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read events from event file
% initialize events
ne = 0;                                                             % number of events
events = zeros(4,0);
annotations = cellstr('empty');

% find appropriate file
fid = fopen(EventFile,'rt');

if fid < 0
    fid = fopen(EventFile2,'rt');
end

if fid >= 0
    cell = textscan(fid,'%s','whitespace','','endofline',char(167));
    fclose(fid);
    cont = cell2mat(cell{1});

    % scan in cell 1 for NUMBER_LIST (occurs five times)
    ix = strfind(cont,'NUMBER_LIST');

    newlines = ix(4) - 1 + strfind(cont(ix(4):ix(5)),char(10));     % newline
    last = size(newlines,2)-1;
    for j = 1:last                                                  % loop over labels
        text = cont(newlines(j)+1:newlines(j+1)-1);
        tcell = textscan(text,'%d');
        sample = tcell{1}(1);                                       % access more content using different columns
        type = tcell{1}(3);
        startsample = tcell{1}(5);
        endsample = tcell{1}(6);
        ne = ne + 1;
        events = cat ( 2, events, [ sample; type; startsample; endsample ] );
    end

    % scan in cell 1 for REMARK_LIST (occurs five times)
    ix = strfind(cont,'REMARK_LIST');
    na = 0;

    newlines = ix(4) - 1 + strfind(cont(ix(4):ix(5)),char(10));     % newline
    last = size(newlines,2)-1;
    for j = 1:last                                                  % loop over labels
        text = cont(newlines(j)+1:newlines(j+1)-1);
        na = na + 1;
        annotations(na) = cellstr(text);
    end
end

if nargout>1

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % read data file
    if nASCII == 1
        fid = fopen_or_error(DataFile,'rt');
        cell = textscan(fid,'%f',nChannels*nSamples*nTrials);
        fclose(fid);
        data = reshape([cell{1}],nChannels,nSamples*nTrials);
    else
        fid = fopen_or_error(DataFile,'rb');
        data = fread(fid,[nChannels,nSamples*nTrials],'float32');
        fclose(fid);
    end

    if nSamples*nTrials ~= size(data,2)
        msgbox('Inconsistent number of samples. File may be displayed incompletely.', 'WARNING', 'warn')
        nSamples = size(data,2)/nTrials;
    end
    % transpose?
    if nMultiplex == 1
        data = data';
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% time axis
time = linspace(fOffsetUsec/1000,fOffsetUsec/1000+(nSamples*nTrials-1)*1000/fFrequency,nSamples*nTrials);

%{
% simple plot
subplot(2,1,1);
plot(time,data);

% stacked plot
subplot(2,1,2);
range = max([abs(min(min(data))) abs(max(max(data)))]);
shift = linspace((nChannels-1)*range*0.3,0,nChannels);
plot(time,data+repmat(shift,nSamples*nTrials,1)');
set(gca,'ytick',flip(shift),'yticklabel',flip(labels),'GridLineStyle',':','XGrid','on','YGrid','off');
ylim([min(min(data+repmat(shift,nSamples*nTrials,1)')) max(max(data+repmat(shift,nSamples*nTrials,1)'))]);
%}

orig = [];
orig.nSamples    = nSamples;
orig.nChannels   = nChannels;
orig.nTrials     = nTrials;
orig.fFrequency  = fFrequency;
orig.fOffsetUsec = fOffsetUsec;
orig.nASCII      = nASCII;
orig.nMultiplex  = nMultiplex;
orig.fSampleTime = fSampleTime;
orig.events      = events;
orig.labels      = labels;
orig.time        = time;
orig.impedancematrix = impedancematrix;
orig.sensorpos       = sensorpos;
orig.epochlabelslist = epochlabelslist;
orig.epochinformationlist = epochinformationlist;
orig.ne          = ne;
orig.events      = events;
orig.annotations = annotations;


text:
function [nu, g] = orderedNoiseUpdateParams(noise, mu, varsigma, y, index)

% ORDEREDNOISEUPDATEPARAMS Update parameters for ordered categorical noise model.

% NOISE

% NOISE

[g, dlnZ_dvs] = orderedNoiseGradVals(noise, mu(index, :), ...
                                            varsigma(index, :), ...
                                            y(index, :));

nu = g.*g - 2*dlnZ_dvs;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%%
% Test for power diagrams

addpath('../toolbox/');
addpath('power_bounded/');
addpath('power_diagrams/');

rep = '../results/semi-discrete/';
[~,~] = mkdir(rep);

% number of sites
N = 200;
% x and y coordinate of the sites
xy = rand(N, 2);
sigma = .1;
xy = sigma*randn(N, 2)+1/2;

% weights
w = ones(N,1);

% the bounding box in clockwise order
bb = [0,0; 0,1; 1,1; 1,0];
% get the power diagram
[V,C] = power_bounded(xy(:,1),xy(:,2), w, bb);
% [PD, PDinf] = powerDiagramWrapper(xy, w);
% draw the resulted power diagram
plot_power(xy,V,C,bb);

q = .01;
axis([-q 1+q -q 1+q]);

% target histogram
Atgt = ones(N,1)/N;

% optimization run
tau = .0004 * N; % gradient step size
niter = 200;
iter_disp = round( [1:8 niter/16 niter/8 niter/4 niter/2 niter] );
iter_disp = 1:2:niter;
err = []; kdisp = 0;
for i=1:niter
    progressbar(i,niter);
    [V,C] = power_bounded(xy(:,1),xy(:,2), w, bb);
    A = area_power(xy,V,C,bb);
    if intersect(i,iter_disp)
        kdisp = kdisp+1;
        clf; plot_power(xy,V,C,bb); axis([-q 1+q -q 1+q]);
        saveas(gcf, [rep 'iteration-' num2str(kdisp), '.eps'], 'epsc');
        drawnow;
    end
    A = A/sum(A); % be sure to be normalized ...
    wnew = w - tau * (A(:)-Atgt(:));
    err(i) = norm(w-wnew);
    w = wnew;
end

% error decay
plot(err); axis tight;

% plot final OT
clf; plot_power(xy,V,C,bb, 2);
saveas(gcf, [rep 'matching.eps'], 'epsc');
axis([-q 1+q -q 1+q]);

author: optimaltransport


text:
function coast = ma_createNBodyCoast(name, coastType, coastToValue, revs, refBody, vars, soiSkipIds, lineColor, lineStyle, lineWidth, massLoss, forceModel, maxPropTime)
%ma_createNBodyCoast Summary of this function goes here
%   Detailed explanation goes here

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Valid coast type strings:
    % goto_ut        - goes to a specific Universal Time
    % goto_dt        - goes to a specific delta-time from the input state
    % goto_tru       - goes to a specific true anomaly
    % goto_apo       - goes to the next apogee
    % goto_peri      - goes to the next perigee
    % goto_asc_node  - goes to ascending node
    % goto_desc_node - goes to the descending node (asc node + pi)
    % goto_soi_trans - goes to the next SOI transition
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    coast = struct();
    coast.name          = name;
    coast.type          = 'NBodyCoast';
    coast.coastType     = coastType;
    coast.coastToValue  = coastToValue;
    coast.revs          = revs;
    coast.refBody       = refBody;
    coast.id            = rand(1);
    coast.vars          = vars;
    coast.soiSkipIds    = soiSkipIds;
    coast.forceModel    = forceModel;
    coast.lineColor     = lineColor;
    coast.lineStyle     = lineStyle;
    coast.lineWidth     = lineWidth;
    coast.massloss      = massLoss;
    coast.maxPropTime   = maxPropTime;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function X = positive_definite_karcher_mean(A)
% Computes a Karcher mean of a collection of positive definite matrices.
%
% function X = positive_definite_karcher_mean(A)
%
% Input:  A 3D matrix A of size nxnxm such that each slice A(:,:,k) is a
%         positive definite matrix of size nxn.
%
% Output: A positive definite matrix X of size nxn which is a Karcher mean
%         of the m matrices in A, that is, X minimizes the sum of squared
%         Riemannian distances to the matrices in A:
%            f(X) = sum_k=1^m .5*dist^2(X, A(:, :, k))
%         The distance is defined by the natural metric on the set of
%         positive definite matrices: dist(X,Y) = norm(logm(X\Y), 'fro').
%
% This simple example is not the best way to compute Karcher means. Its
% purpose it to serve as base code to explore other algorithms. In
% particular, in the presence of large noise, this algorithm seems to not
% be able to reach points with a very small gradient norm. This may be
% caused by insufficient accuracy in the gradient computation.

% This file is part of Manopt and is copyrighted. See the license file.
%
% Main author: Nicolas Boumal, Sept. 3, 2013
% Contributors:
%
% Change log:
%

    % Generate some random data to test the function if none is given.
    if ~exist('A', 'var') || isempty(A)
        n = 5;
        m = 10;
        A = zeros(n, n, m);
        ref = diag(max(.1, 1+.1*randn(n, 1)));
        for i = 1 : m
            noise = 0.01*randn(n);
            noise = (noise + noise')/2;
            [V, D] = eig(ref + noise);
            A(:, :, i) = V*diag(max(.01, diag(D)))*V';
        end
    end

    % Retrieve the size of the problem:
    % There are m matrices of size nxn to average.
    n = size(A, 1);
    m = size(A, 3);
    assert(n == size(A, 2), ...
           ['The slices of A must be square, i.e., the ' ...
	        'first and second dimensions of A must be equal.']);

    % Our search space is the set of positive definite matrices of size n.
    % Notice that this is the only place we specify on which manifold we
    % wish to compute Karcher means. Replacing this factory for another
    % geometry will yield code to compute Karcher means on that other
    % manifold, provided that manifold is equipped with a dist function and
    % a logarithmic map log.
    M = sympositivedefinitefactory(n);

    % Define a problem structure, specifying the manifold M, the cost
    % function and its gradient.
    problem.M = M;
    problem.cost = @cost;
    problem.grad = @grad;

    % Explicitly pick an approximate Hessian for the trust-region method
    problem.approxhess = approxhessianFD(problem, struct('stepsize', 1e-4));

    % The functions below make many redundant computations. This
    % performance hit can be alleviated by using the caching system. We go
    % for a simple implementation here, as a tutorial example.

    % Cost function
    function f = cost(X)
        f = 0;
        for k = 1 : m
            f = f + M.dist(X, A(:, :, k))^2;
        end
        f = f/(2*m);
    end

    % Riemannian gradient of the cost function
    function g = grad(X)
        g = M.zerovec(X);
        for k = 1 : m
            % Update g in a linear combination of the form
            % g = g - [something]/m.
            g = M.lincomb(X, 1, g, -1/m, M.log(X, A(:, :, k)));
        end
    end

    % Execute some checks on the derivatives for early debugging.
    % These things can be commented out of course.
    % The slopes should agree on part of the plot at least. In this case,
    % it is sometimes necessary to inspect the plot visually to make the
    % call, but it is indeed correct.
    % checkgradient(problem);
    % pause;

    % Execute this if you want to force using a proper parallel vector
    % transport. This is not necessary. If you omit this, the default
    % vector transport is the identity map, which is (of course) cheaper
    % and seems to perform well in practice.
    % M.transp = M.paralleltransp;

    % Issue a call to a solver. Default options are selected.
    % Our initial guess is the first data point.
    X = trustregions(problem, A(:, :, 1));

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
classdef tattention < matlab.unittest.TestCase
    % tattention   Tests for transformer.layer.attention

    % Copyright 2020 The Mathworks, Inc.

    % The attention function has 3 purposes:
    % 1. Apply the pre and post attention fully connected layers.
    % 2. Rearrange data format for applying multi-head attention.
    % 3. Preserve the key and value matrices, and append to them when using
    % values from the past. This is done pre-attention.

    properties(Constant,Access=private)
        attention = @transformer.layer.attention
        multiheadAttention = @transformer.layer.multiheadAttention
    end

    properties(TestParameter)
        NumQueries = {1,2}
        NumObs = {1,3}
    end

    methods(Test)
        function checkSingleHeadNoFullyConnected(test,NumQueries,NumObs)
            % Verify that multiheadAttention is used.
            % Do this by setting the fully connected weights to identity
            % and the bias to 0.
            % With a single head we do not need split and merge heads.
            latentDim = 10;
            hyperParams.NumHeads = 1;
            past = [];
            % Set up the fc weights to be identity matrices and biases to
            % be 0.
            weights = test.prepareWeightsStructWithIdentityFC(hyperParams.NumHeads,latentDim);
            % Call attention on an arbitrary input.
            x = dlarray(rand(latentDim*3,NumQueries,NumObs));
            yAct = test.attention(x,past,weights,hyperParams);
            % Verify against multiheadAttention by splitting the arbitrary
            % input into query, key and value.
            [q,k,v] = iSplitQKV(x,hyperParams.NumHeads,latentDim);
            yExp = test.multiheadAttention(q,k,v);
            % Note that for numObs > 1 we need to merge heads. Since we
            % have a single head here, no merger needs to take place -- we
            % just need to permute the observation dimension from 4 to 3
            yExp = permute(yExp, [1 2 4 3]);
            test.verifyEqual(yAct,yExp);
        end

        function checkMultiHeadNoFullyConnected(test,NumQueries,NumObs)
            % Verify the multiple head case. Multi head attention is
            % achieved by simply creating a "head dimension" and pushing it
            % to the back so that it is treated as a "batch"
            % dimension for matrix multiplication.
            latentDim = 10;
            hyperParams.NumHeads = 2;
            past = [];
            % Set up the fc weights to be identity matrices and biases to
            % be 0.
            weights = test.prepareWeightsStructWithIdentityFC(hyperParams.NumHeads,latentDim);
            % Call attention on an arbitrary input.
            x = dlarray(rand(latentDim*hyperParams.NumHeads*3,NumQueries,NumObs));
            yAct = test.attention(x,past,weights,hyperParams);
            [Q,K,V] = iSplitQKV(x,hyperParams.NumHeads,latentDim);
            % Verify against the multiheadAttention function.
            yExp = test.multiheadAttention(Q,K,V);
            yExp = iMergeHeads(yExp);
            test.verifyEqual(yAct,yExp);
        end

        function checkPastPresentCaching(test,NumQueries,NumObs)
            % Verify the 2nd input and output of attention - the pasts
            % passed in are the keys and values for the previous time step.
            % These are concatenated to the keys and values for the current
            % time step before the multiheadAttention call, and these
            % concatenated key and values are passed out as a second
            % output.
            latentDim = 10;
            hyperParams.NumHeads = 2;
            % Set up a fake past by making an initial call to attention.
            past = [];
            % Set up the fc weights to be identity matrices and biases to
            % be 0.
            weights = test.prepareWeightsStructWithIdentityFC(hyperParams.NumHeads,latentDim);
            % Call attention on an arbitrary input.
            x = dlarray(rand(latentDim*hyperParams.NumHeads*3,NumQueries,NumObs));
            [~,past] = test.attention(x,past,weights,hyperParams);
            % Verify the expected value of past - it is the key and values
            % concatenated on the 4th dimension.
            [~,K,V] = iSplitQKV(x,hyperParams.NumHeads,latentDim);
            test.verifyEqual(past,cat(5,K,V));
            % Now verify second call to attention is possible with the first
            % past as input - and verify the value of the attention output.
            [yAct,present] = test.attention(x,past,weights,hyperParams);
            [Q,K,V] = iSplitQKV(x,hyperParams.NumHeads,latentDim);
            % Verify the correct value for present.
            pastK = past(:,:,:,:,1);
            pastV = past(:,:,:,:,2);
            test.verifyEqual(extractdata(present),extractdata(cat(5,cat(2,pastK,K),cat(2,pastV,V))),'AbsTol',1e-5);
            % To compute the expected value, concatenate the pasts
            K = cat(2,K,pastK);
            V = cat(2,V,pastV);
            yExp = test.multiheadAttention(Q,K,V);
            yExp = iMergeHeads(yExp);
            test.verifyEqual(extractdata(yAct),extractdata(yExp),'AbsTol',1e-5);
        end

        function checkInputOutputFC(test,NumQueries,NumObs)
            % The tests above ensure multiheadAttention is used by
            % attention. Now verify the input and output FC operations are
            % used.
            latentDim = 10;
            hyperParams.NumHeads = 2;
            % Set up a fake past by making an initial call to attention.
            past = [];
            params = test.prepareWeightsStructUsingGenerators(hyperParams.NumHeads,latentDim,@randn,@randn);
            % Call attention on an arbitrary input.
            x = dlarray(rand(latentDim*hyperParams.NumHeads*3,NumQueries,NumObs));
            yAct = test.attention(x,past,params,hyperParams);
            % Verify this matches the full attention implementation
            z = fullyconnect(x,params.attn_c_attn_w_0,params.attn_c_attn_b_0,'DataFormat','CBT');
            [Q,K,V] = iSplitQKV(z,hyperParams.NumHeads,latentDim);
            z = test.multiheadAttention(Q,K,V);
            z = iMergeHeads(z);
            yExp = fullyconnect(z,params.attn_c_proj_w_0,params.attn_c_proj_b_0,'DataFormat','CBT');
            test.verifyEqual(extractdata(yAct),extractdata(yExp),'AbsTol',1e-5);
        end

        function defaultIsMaksed(test)
            % Verify the default for the 'CausalMask' NVP is true
            latentDim = 10;
            hyperParams.NumHeads = 2;
            past = [];
            numQueries = 2;
            % Set up the fc weights to be identity matrices and biases to
            % be 0.
            weights = test.prepareWeightsStructWithIdentityFC(hyperParams.NumHeads,latentDim);
            % Call attention on an arbitrary input.
            x = dlarray(rand(latentDim*hyperParams.NumHeads*3,numQueries));
            default = test.attention(x,past,weights,hyperParams);
            masked = test.attention(x,past,weights,hyperParams,'CausalMask',true);
            test.verifyEqual(extractdata(default),extractdata(masked));
        end

        function canTurnOffMask(test)
            % Verify the 'CausalMask' NVP can be set to false - the expectation
            % is this simply sets 'Masked' to false for the
            % multiheadAttention call.
            latentDim = 10;
            hyperParams.NumHeads = 2;
            past = [];
            numQueries = 2;
            % Set up the fc weights to be identity matrices and biases to
            % be 0.
            weights = test.prepareWeightsStructWithIdentityFC(hyperParams.NumHeads,latentDim);
            % Call attention on an arbitrary input.
            x = dlarray(rand(latentDim*hyperParams.NumHeads*3,numQueries));
            yAct = test.attention(x,past,weights,hyperParams,'CausalMask',false);
            [Q,K,V] = iSplitQKV(x,hyperParams.NumHeads,latentDim);
            % Verify against the multiheadAttention function with 'Masked'
            % set to false.
            yExp = test.multiheadAttention(Q,K,V,'CausalMask',false);
            yExp = iMergeHeads(yExp);
            test.verifyEqual(yAct,yExp);
        end

        function canDropout(test)
            % Verify dropout can be applied to attention - the expectation
            % is that this matches simply setting the Dropout NVP in the
            % multiheadAttention call
            latentDim = 10;
            hyperParams.NumHeads = 2;
            past = [];
            numQueries = 2;
            % Set up the fc weights to be identity matrices and biases to
            % be 0.
            weights = test.prepareWeightsStructWithIdentityFC(hyperParams.NumHeads,latentDim);
            % Call attention on an arbitrary input.
            x = dlarray(rand(latentDim*hyperParams.NumHeads*3,numQueries));
            p = 0.5;
            % Set global rng between non-deterministic calls.
            rng(0);
            yAct = test.attention(x,past,weights,hyperParams,'Dropout',p);
            [Q,K,V] = iSplitQKV(x,hyperParams.NumHeads,latentDim);
            % Verify against the multiheadAttention function with 'Dropout'
            % set to p
            % Set global rng between non-deterministic calls.
            rng(0);
            yExp = test.multiheadAttention(Q,K,V,'Dropout',p);
            yExp = iMergeHeads(yExp);
            test.verifyEqual(yAct,yExp);
        end

        function defaultIsNoDropout(test)
            % Verify the default Dropout is 0.
            latentDim = 10;
            hyperParams.NumHeads = 2;
            numQueries = 2;
            % Set up a fake past by making an initial call to attention.
            past = [];
            params = test.prepareWeightsStructUsingGenerators(hyperParams.NumHeads,latentDim,@randn,@randn);
            % Call attention on an arbitrary input.
            x = dlarray(rand(latentDim*hyperParams.NumHeads*3,numQueries));
            yDefault = test.attention(x,past,params,hyperParams);
            yNoDropout = test.attention(x,past,params,hyperParams,'Dropout',0);
            test.verifyEqual(yDefault,yNoDropout);
        end
    end

    methods(Access=private)
        function s = prepareWeightsStruct(~,W1,b1,W2,b2)
            % Prepare a struct compatible with the weights input of
            % attention. These are for the fully connected layers.
            s = struct(...
                'attn_c_attn_w_0',W1,...
                'attn_c_attn_b_0',b1,...
                'attn_c_proj_w_0',W2,...
                'attn_c_proj_b_0',b2);
        end

        function s = prepareWeightsStructUsingGenerators(test,nHeads,latentDim,weightGenerator,biasGenerator)
            % Use function handles weightGenerator and biasGenerator to
            % create weight and bias values for the fc layers in attention.
            % The expectation is these function handles map a size vector
            % to an array of that size.

            % The query, key and value vectors are flattened into 1
            % dimension for the input.
            % There are nHeads of each, and each is a latentDim
            % dimensional vector
            inputDim = 3*nHeads*latentDim;
            W1 = dlarray(weightGenerator([inputDim,inputDim]));
            b1 = dlarray(biasGenerator([inputDim,1]));

            % For the output we don't need the 3 - the multiheadAttention
            % has taken an expectation of the value vectors under some
            % probability distribution.
            outputDim = nHeads*latentDim;
            W2 = dlarray(weightGenerator([outputDim,outputDim]));
            b2 = dlarray(biasGenerator([outputDim,1]));
            s = test.prepareWeightsStruct(W1,b1,W2,b2);
        end

        function s = prepareWeightsStructWithIdentityFC(test,nHeads,latentDim)
            % Prepare a struct compatible with the weights input of
            % attention such that the fc layers are identity operations.
            s = test.prepareWeightsStructUsingGenerators(nHeads,latentDim,@eye,@zeros);
        end
    end
end

function [Q,K,V] = iSplitQKV(X, nHeads, latentDim)
% Split an input array X into the corresponding Q, K and V arrays.
splitSize = latentDim*nHeads;
Q = iSplitHeads(X(1:splitSize,:,:),splitSize,nHeads);
K = iSplitHeads(X((splitSize+1):2*splitSize,:,:),splitSize,nHeads);
V = iSplitHeads(X((2*splitSize+1):3*splitSize,:,:),splitSize,nHeads);
end

function X = iSplitHeads(X, splitSize, numHeads)
X = reshape(X, splitSize/numHeads, numHeads, [], size(X,3));   % Split states
X = permute(X,[1 3 2 4]);
end

function X = iMergeHeads(X)
X = permute(X, [1 3 2 4]);
X = reshape(X, size(X,1)*size(X,2), [], size(X,4));            % Merge states
end

author: matlab-deep-learning

repo: transformer-models


text:
function problem1 ( )

%*****************************************************************************80
%
%% PROBLEM1 is the main program for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2009
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PROBLEM1\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  A test problem for FD1D_HEAT_STEADY.\n' );

  n = 11;

  a = 0.0;
  b = 1.0;

  ua = 0.0;
  ub = 1.0;

  [ x, u ] = fd1d_heat_steady ( n, a, b, ua, ub, @k1, @f1 );

  plot ( x, u )
  title ( 'Problem 1: Uniform K(X)' );
  xlabel ( 'A <= X <= B' );
  ylabel ( 'Temperature' );

  filename = 'problem1_nodes.txt';
  r8mat_write ( filename, 1, n, x );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  X data written to "%s".\n', filename );
  filename = 'problem1_values.txt';
  r8mat_write ( filename, 1, n, u' );
  fprintf ( 1, '  U data written to "%s".\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PROBLEM1\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );

  timestamp ( );

  return
end
function value = k1 ( x )

%*****************************************************************************80
%
%% K1 evaluates the heat transfer coefficient K(X).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the position.
%
%    Output, real VALUE, the value of K(X).
%
  value = 1.0;

  return
end
function value = f1 ( x )

%*****************************************************************************80
%
%% F1 evaluates the right hand side function F(X).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the position.
%
%    Output, real VALUE, the value of F(X).
%
  value = 0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [kdynpcm2] = kPa2dynpcm2(kPa)
% Convert units of pressure from kilopascals to dynes per centimeter squared.
% Chad A Greene 2012
dynpcm2 = kPa*10000;


text:
function iOut = iVals(trialOrder,ignoreNulls);
% iOut = iVals(trialOrder,ignoreNulls);
% iVals: Calculate i-value (intervening trials) for fMR-A Expts (but maybe more generally useful).
%
% Here's the idea: trialOrder is an array of numbers, generally with values that are repeated.
% The i-value of each element counts the number of intervening entries between repeated values.
% If the value for a given element in trialOrder has not occurred before in trialOrder, i = -1.
% If the value has occurred immediately before the current element, i = 0;
% If the value last occurred two elements before, i = 1 (one element between repetitions), and so on.
%
% So if trialOrder is:
%
% 1 1 2 1 3 3 1 1 1 2
%
% nOut will be:
%
% -1 0 -1 1 -1 -1 2 0 0 6
%
% If ignoreNulls is set to 1, iOut will ignore intervening elements with a value of 0 in counting
% (e.g., 1 0 0 1 will have an iOut of -1 -1 -1 0 instead of -1 -1 0 2).
%
% 1/03 by ras.
% 10/03 ras: updated to work with matrices. If a matrix is passed, it will run through eac
% row and get the ivals.
if ~exist('ignoreNulls','var')	ignoreNulls = 0;		end

if size(trialOrder,1) > 1 & size(trialOrder,2) > 1
	iOut = [];
	for i = 1:size(trialOrder,1)
		iOut = [iOut; iVals(trialOrder(i,:))];
	end
	return
end

if ignoreNulls
	nonNulls = find(trialOrder~=0);
	nulls = find(trialOrder==0);
	fullOrder = trialOrder;
	trialOrder = trialOrder(nonNulls);
end

for i = 1:length(trialOrder)
	currVal = trialOrder(i);
	ind = find(trialOrder==currVal);
	if ind(1)==i
		iOut(i) = -1;
	else
		loc = find(ind==i);
		lastInstance = ind(loc-1);
		iOut(i) = i - lastInstance - 1;
	end
end

if ignoreNulls
	tmp = iOut;
	iOut = zeros(1,length(fullOrder));
	iOut(nonNulls) = tmp;
	iOut(nulls) = -1;
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function d=v_distchar(ar1,ar2,mode)
%V_DISTCHAR calculates the cosh spectral distance between AR coefficients D=(AR1,AR2,MODE)
%
% Inputs: AR1,AR2     AR coefficient sets to be compared. Each row contains a set of coefficients.
%                     AR1 and AR2 must have the same number of columns.
%
%         MODE        Character string selecting the following options:
%                         'x'  Calculate the full distance matrix from every row of AR1 to every row of AR2
%                         'd'  Calculate only the distance between corresponding rows of AR1 and AR2
%                              The default is 'd' if AR1 and AR2 have the same number of rows otherwise 'x'.
%
% Output: D           If MODE='d' then D is a column vector with the same number of rows as the shorter of AR1 and AR2.
%                     If MODE='x' then D is a matrix with the same number of rows as AR1 and the same number of columns as AR2'.
%
% The COSH spectral distance is the average over +ve and -ve frequency of
%
%                     cosh(log(p1/p2))-1   =   (p1-p2)^2/(2p1*p2)   =   (p1/p2 + p2/p1)/2 - 1
%
% Where p1 and p2 are the power spectra corresponding to the AR coefficient sets AR1 and AR2.
% The COSH distance is a symmetrical version of the Itakura-Saito distance: v_distchar(x,y)=(v_distisar(x,y)+v_distisar(y,x))/2

% Since the power spectrum is the fourier transform of the autocorrelation, we can calculate
% the average value of p1/p2 by taking the 0'th order term of the convolution of the autocorrelation
% functions associated with p1 and 1/p2. Since 1/p2 corresponds to an FIR filter, this convolution is
% a finite sum even though the autocorrelation function of p1 is infinite in extent.

% The Cosh distance can also be calculated directly from the power spectra; providing np is large
% enough, the values of d0 and d1 in the following will be very similar:
%
%         np=255; d0=v_distchar(ar1,ar2); d1=v_distchpf(v_lpcar2pf(ar1,np),v_lpcar2pf(ar2,np))
%

% Ref: A.H.Gray Jr and J.D.Markel, "Distance measures for speech processing", IEEE ASSP-24(5): 380-391, Oct 1976
%      L. Rabiner abd B-H Juang, "Fundamentals of Speech Recognition", Section 4.5, Prentice-Hall 1993, ISBN 0-13-015157-2

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: v_distchar.m 10865 2018-09-21 17:22:45Z 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf1,p1]=size(ar1);
nf2=size(ar2,1);
p2=p1+1;
m1=zeros(nf1,2*p1);
m2=zeros(nf2,2*p1);
m1(:,1:p1)=v_lpcar2rr(ar1);
m1(:,p2:end)=v_lpcar2ra(ar1);
m1(:,1)=m1(:,1)*0.5;
m1(:,p2)=m1(:,p1+1)*0.5;
m2(:,p2:end)=v_lpcar2rr(ar2);
m2(:,1:p1)=v_lpcar2ra(ar2);

if nargin<3 | isempty(mode) mode='0'; end
if any(mode=='d') | (mode~='x' & nf1==nf2)
   nx=min(nf1,nf2);
   d=sum(m1(1:nx,:).*m2(1:nx,:),2)-1;
else
   d=m1*m2'-1;
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [ChipsOut, Scrambler] = PacketBuilder(DataBits, G, Gs);
%
% PACKET BUILDER 	This function build the IS-95 forward Channel Data Packets.
% 						Block Diagram:
%						InData -> [Viterbi Encoder] -> [Interleaver] -> [Scrambler]
%
% 						9.6 KBps for 20 msec packet --> 192 bits (including 8 zeros tail) 184 data bits netto
% 						9.6*2 = 19.2 KBps (384 bits) after Viterbi incoding
% 						19.2 KBps * 64 = 1.2288 Mbps after Walsh and Spreading
%
% 						Inputs: DataBits - data to transmit (binary form)
%                         G - Viterbi Encoder generation polynom
%                         Gs - Long sequence generation polynom (scrambler)
%
% 						Outputs: ChipsOut - The result chip sequnece entered to modulator (binary form)
%                          Scrambler - Scrambler sequence%

global Zs

K = size(G, 2); 		% memory of Viterbi Polynoms
L = size(G, 1);      % number of chips per each data bit

N = 64*L*(length(DataBits)+K-1);		% Number of chips (9.6 Kbps -> 1.288 Mbps)

%====================== PACKET BUILDER =================
%-------- Create a Data Block (184 data bits) concatenated by 8 zeros
chips = VitEnc(G, [DataBits; zeros(K-1,1)]);	% Viterbi Encoder

%-------- Interleaver ---------
INTERL = reshape(chips, 24, 16);			% IN-> columns, OUT-> rows
chips = reshape(INTERL', length(chips), 1);  % Rate = 19.2 KBps

%----------- Scrambler---------
% Rate = 19.2 KBps
[LongSeq Zs] = PNGen(Gs, Zs, N);  % Long sequence Generation (at rate 1.2288 Mbps)
Scrambler = LongSeq(1:64:end); 	 % Decimation of Long Sequence

ChipsOut = xor(chips, Scrambler);


text:
%% load param
param.net_def_dir = 'model/ResNet_3b_s16_f2r';
param.init_model = fullfile(pwd, 'output/ResNet_3b_s16/tot_wometa_1epoch'); %fullfile(pwd, 'output/ResNet-101L/1000000');
param.solverfile = 'solver_test.prototxt';
param.gpu_id = 0;
mins = 3;
maxs = 14;
%% init caffe solver
caffe.reset_all;
caffe_solver = caffe.get_solver(fullfile(param.net_def_dir, param.solverfile), param.gpu_id);
if ~isempty(param.init_model)
    assert(exist(param.init_model)==2, 'Cannot find caffemodel.');
    caffe_solver.use_caffemodel(param.init_model);
end
caffe_solver.set_phase('test');
num = length(featmap_trans);
cons = 0;
hit = 0;
sums = 0;
missed = cell(num, 1);
recall = 0;
for i = 1 : num
    ti = tic();
    drawnow;
    fprintf('Step 3 rpn:  %d/%d...', i, num);
    parsed(i) = detect_all_by_featmap( featmap_trans{i}, unique(min(scale{i}, 5)), param, caffe_solver, mins);
    toc(ti);
end
fprintf('Modify scale...')
for i = 1 : num
    s = max(size(imread(fullfile(param.test_root, list{i}))))/param.max_img;
    parsed(i).box(:,1:4) = parsed(i).box(:,1:4) * s;
    parsed(i).point = parsed(i).point * s;
end
fprintf('Done\n');

author: liuyuisanai

repo: RSA-for-object-detection


text:
function y = loadPcwConst(type, nSamples)
%LOADPCWCONST Creates a piecewise constant signal or image

% written by M. Storath
% $Date: 2014-09-02 14:56:23 +0200 (Di, 02. Sep 2014) $	$Revision: 105 $

if not(exist('nSamples', 'var'))
    nSamples = 2^8;
end

x = linspace(-1, 1, nSamples)';
y = zeros(size(x));

switch type
    case 'rect'
        y = (-0.5 < x) & ( x < 0.5);
    case {'step', 'heaviside'}
        y = x >= 0;
    case 'jumps'
        % two jumps
        y = (-0.3 < x) & ( x < -0.1) + 2 * ((0.1 < x) & ( x < 0.4));
    case 'equidistant'
        n = 8; % number of partitions
        step = nSamples / n;
        heights = [3, 1, 7, 6, 5, 0,6,4];
        for i = 1:n
            iv = ((i-1) * step + 1) : (i * step);
            y(iv) = heights(i) / max(heights);
        end

    case 'sampleDec'
        heights = [3, 1, 7, 6, 5, 0,6,6];
        jumps = diff(heights);
        jumpDist = [30, 40, 32, 32, 32, 50, 60];
        jumpIdx = cumsum(jumpDist);
        y(jumpIdx) = jumps;
        y = cumsum(y);
        y = mat2gray(y);

    case 'sampleDec2'
        heights = [3, 1, 7, 6, 5, 0,6,6];
        jumps = diff(heights);
        jumpDist = [30, 40, 32, 32, 32, 50, 40];
        jumpIdx = cumsum(jumpDist);
        y(jumpIdx) = jumps;
        y = cumsum(y);
        y = mat2gray(y);

    case 'equidistant2'
        n = 8; % number of partitions
        step = nSamples / n;
        heights = [3, 1, 7, 6, 5, 4,6,4];
        for i = 1:n
            iv = ((i-1) * step + 1) : (i * step);
            y(iv) = heights(i) / max(heights);
        end

    case 'sample1'
        steps = [2, 4, 3, 3, 2, 5, 3, 3, 3];
        steps = (cumsum(steps ./ sum(steps)) - 0.5) * 2;
        heights = [3, 1, 7, 5.5, 4.5, 3.5,6,0];
        n= numel(heights);
        for i = 1:n
            idx = find((steps(i) < x) & (x <= steps(i+1)));
            y( idx ) = heights(i) / max(heights);
        end

    case 'sample2'
        steps = [2, 4, 3, 3, 2, 5, 2, 3, 3,6,4,2,8,3,4, 4]; %15
        steps = (cumsum(steps ./ sum(steps)) - 0.5) * 2;
        heights = [3, 1, 7, 6, 5, 4,6,0, 2, 4, 3, 1.5, 7, 9, 2];
        n= numel(heights);
        for i = 1:n
            idx = find((steps(i) < x) & (x <= steps(i+1)));
            y( idx ) = heights(i) / max(heights);
        end

   case 'sample3'
        steps = [2, 4, 3, 6, 4, 5, 4, 3, 3];
        steps = (cumsum(steps ./ sum(steps)) - 0.5) * 2;
        heights = [3, 1, 7, 6, 5, 4,6,0];
        n= numel(heights);
        for i = 1:n
            idx = find((steps(i) < x) & (x <= steps(i+1)));
            y( idx ) = heights(i) / max(heights);
        end

   case 'sample4'
        %steps = [3, 2, 5, 10, 4, 3, 3];
        %steps = [ -Inf, (cumsum(steps ./ sum(steps)) - 0.5) * 2];
        %heights = [4, 2,3, 3.1, 3.05, 1, 1];
        steps = [2, 4, 2, 2, 1, 1];
        steps = [ -Inf, (cumsum(steps ./ sum(steps)) - 0.5) * 2];
        heights = [0.0, 0.1, 0.05, 1, 0.5, 0];
        n= numel(heights);
        for i = 1:n
            idx = find((steps(i) < x) & (x <= steps(i+1)));
            y( idx ) = heights(i) / max(abs(heights));
        end

   case 'sample5'
        steps = [3, 5, 3, 4, 6, 4, 3, 5, 3];
        steps = (cumsum(steps ./ sum(steps)) - 0.5) * 2;
        heights = [3, 7, 2, 4, 5, 2,6,0];
        n= numel(heights);
        for i = 1:n
            idx = find((steps(i) < x) & (x <= steps(i+1)));
            y( idx ) = heights(i) / max(heights);
        end

    case 'geo1'
        y = double(imread('geo1.png'))/255;

    case 'geo2'
        y = double(imread('geo2.png'))/255;

    case 'geo3'
        y = double(imread('geo3.png'))/255;

    case 'geo4'
        y = double(imread('geo4.png'))/255;

    case 'geo5'
        y = double(imread('geo5.png'))/255;

    case 'geo6'
        y = double(imread('geo6.png'))/255;

    case 'geo7'
        y = double(imread('geo7.png'))/255;

    case 'geo8'
        y = double(imread('geo8.png'))/255;

    case 'overlay'
        I1 = double(imread('rectangle.png'));
        I2 = double(imread('octo.png'));
        y =  (I1 + 2*I2)/ (3*255);

    otherwise
        error('This option does not exist.')
end

% cast to double
y = double(y);

end

author: mstorath

repo: Pottslab

save_path: github-repos/MATLAB/mstorath-Pottslab


text:
function phi = correct_azimuth(phi)
%CORRECT_AZIMUTH ensures azimuth angle between -pi and +pi-eps
%
%   Usage: phi = correct_azimuth(phi)
%
%   Input parameters:
%       phi     - azimuth / rad. Can be a single value or a matrix.
%
%   Output paramteres:
%       phi     - angle between -pi and +pi-eps / rad
%
%   See also: correct_elevation, get_ir

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax);

%% ===== Computation ====================================================
% Ensure -2pi <= phi <= 2pi
phi = rem(phi,2*pi);
% Ensure -pi <= phi < pi
phi(phi<-pi) = phi(phi<-pi) + 2*pi;
phi(phi>=pi) = phi(phi>=pi) - 2*pi;

author: sfstoolbox

repo: sfs-matlab


text:
depthfileNamesFile = fopen('sunrgbd_train_depth_fileNames.txt','w');
res_file = fopen('sunrgbd_train_res_fileNames.txt','w');
K_file = fopen('sunrgbd_train_K_fileNames.txt','w');
T_file = fopen('sunrgbd_train_T_fileNames.txt','w');

for i = 1:5285

    i

    k = i+5050;

    fName = mat2str(cell2mat(alltrain(i)));
    fName = strcat('/scratch/SUNRGBD_13',fName(18:end-1),'/depth_bfx/');

    X = dir(fName);

    depth = X(3).name;
    depth = strcat(fName,depth);

    fprintf(depthfileNamesFile,'%s\n',depth);

    Id = double(imread(depth))/10000.0;
%     I = SUNRGBD2Dseg(k).seglabel;

    K = SUNRGBDMeta(k).K;
    Rtilt = SUNRGBDMeta(k).Rtilt;

%     size(Rtilt)

%     X = [1 0 0;
%          0 0 1;
%          0 -1 0];

% %      Rtilt = X * Rtilt * X';

%     Rtilt = X*Rtilt*X';

    fprintf(res_file,'%d %d\n',size(Id,2), size(Id,1));

    fprintf(K_file,'%f %f %f\n',K(1,1),K(1,2),K(1,3));
    fprintf(K_file,'%f %f %f\n',K(2,1),K(2,2),K(2,3));
    fprintf(K_file,'%f %f %f\n\n',K(3,1),K(3,2),K(3,3));

    val_zero = 0.0;

    fprintf(T_file,'%f %f %f %f\n',Rtilt(1,1),Rtilt(1,2),Rtilt(1,3),val_zero);
    fprintf(T_file,'%f %f %f %f\n',Rtilt(2,1),Rtilt(2,2),Rtilt(2,3),val_zero);
    fprintf(T_file,'%f %f %f %f\n\n',Rtilt(3,1),Rtilt(3,2),Rtilt(3,3),val_zero);

%     normals = compute_normals(Id,size(Id,2),size(Id,1),K(1,1), K(2,2), K(1,3), K(2,3), 1);
%
%     [U, V] = meshgrid(1:size(Id,2),1:size(Id,1));
%
%     X = Id.*(U - K(1,3))/K(1,1);
%     Y = Id.*(V - K(2,3))/K(2,2);
%     Z = Id;
%
%     XDash = Rtilt(1,1)*X + Rtilt(1,2)*Y + Rtilt(1,3)*Z;
%     YDash = Rtilt(2,1)*X + Rtilt(2,2)*Y + Rtilt(2,3)*Z;
%     ZDash = Rtilt(3,1)*X + Rtilt(3,2)*Y + Rtilt(3,3)*Z;
%
%     YDash = (YDash - min(min(YDash)));
%
%     normal_x= normals(:,:,1);
%     normal_y= normals(:,:,2);
%     normal_z= normals(:,:,3);
%
%     new_normals_x = Rtilt(1,1)*normal_x + Rtilt(1,2)*normal_y + Rtilt(1,3)*normal_z;
%     new_normals_y = Rtilt(2,1)*normal_x + Rtilt(2,2)*normal_y + Rtilt(2,3)*normal_z;
%     new_normals_z = Rtilt(3,1)*normal_x + Rtilt(3,2)*normal_y + Rtilt(3,3)*normal_z;
%
%     new_normals(:,:,1) = new_normals_x;
%     new_normals(:,:,2) = new_normals_y;
%     new_normals(:,:,3) = new_normals_z;
%
%     angle = acosd(new_normals(:,:,2))/180.0;
%
%     size(normals);
%
%     file1 = fopen(sprintf('train-dha-%04d.bin',i),'wb');
%
%     Id    = imresize(Id,   [224 224],'bilinear');
%     YDash = imresize(YDash,[224 224],'bilinear');
%     angle = imresize(angle,[224 224],'bilinear');
%
%     YDash(isnan(YDash))=0;
%     angle(isnan(angle))=0;
%
%     for i = 1:224
%         for j = 1:224
%             fwrite(file1,single(Id(i,j)),'float');
%             fwrite(file1,single(YDash(i,j)),'float');
%             fwrite(file1,single(angle(i,j)),'float');
%         end
%     end
%
%
%     fclose(file1);

%   [sum(isnan(I(:))), sum(isnan(YDash(:))), sum(isnan(angle(:)))]

end

%     fclose(depthfileNamesFile);
%     fclose(res_file);
%     fclose(T_file);
%     fclose(K_file);

author: ankurhanda

repo: sunrgbd-meta-data


text:
function this = load_data(this)
%ONE_LINE_DESCRIPTION
%
%   Y = MrRoi()
%   Y.load_data(inputs)
%
% This is a method of class MrRoi.
%
% IN
%
% OUT
%
% EXAMPLE
%   load_data
%
%   See also MrRoi

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2014-11-21
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.
%

% get filename
filename = fullfile(this.parameters.save.path, this.parameters.save.fileName);
% load data
data = struct2cell(load(filename, 'data'));
this.data = data{:};

author: translationalneuromodeling

repo: tapas


text:
%function tfrscalt
%TFRSCALT Unit test for the time-frequency representation TFRSCALO.

%	O. Lemoine - June 1996.

% We test each property of the corresponding TFR :

N=128;

% Covariance by translation in time
t1=60; t2=70; f=0.3; W=0;
sig1=amgauss(N,t1).*fmconst(N,f,t1);
sig2=amgauss(N,t2).*fmconst(N,f,t2);
tfr1=tfrscalo(sig1,1:N,W,0.1,0.4,128);
tfr2=tfrscalo(sig2,1:N,W,0.1,0.4,128);
[tr,tc]=size(tfr1);
nu=round(f*(tc-1)*2)+1;
tfr=tfr1-tfr2(:,modulo((1:tc)-t1+t2,tc));
if any(any(abs(tfr)>sqrt(eps))),
 error('tfrscalo test 1 failed');
end

% Covariance by dilation
t=N/2; f=0.2; T=2*sqrt(N); a=2; W=8;
sig1=amgauss(N,t,T).*fmconst(N,f,t);
sig2=amgauss(a*N,a*t,T*a).*fmconst(a*N,f/a,a*t);
[tfr1,t1,f1]=tfrscalo(sig1,1:N  ,W,0.01,0.49,N);
[tfr2,t2,f2]=tfrscalo(sig2,1:a*N,W,0.01,0.49,N);
Max1=max(max(tfr1)); Max2=max(max(tfr2));
[I1,J1]=find(tfr1==Max1); [I2,J2]=find(tfr2==Max2);
if abs(f1(I1)-a*f2(I2))>1e-2 | J2~=a*J1,
 error('tfrscalo test 2 failed');
end

% Reality of the TFR
sig=noisecg(N); W=5;
tfr=tfrscalo(sig,1:N,W,0.01,0.5,N);
if sum(any(abs(imag(tfr))>sqrt(eps)))~=0,
 error('tfrscalo test 3 failed');
end

% Energy conservation
sig=fmsin(N,.1,.4); W=6; Nf=2*N ;
[tfr,t,f]=tfrscalo(sig,1:N,W,0.01,0.49,Nf);
Es=norm(sig)^2/Nf;
Etfr=integ2d(tfr,t,f)/N;
if abs(Es-Etfr)>sqrt(eps),
 error('tfrscalo test 4 failed');
end

% Positivity
if any(any(tfr<0)),
 error('tfrscalo test 5 failed');
end

% Same energy in the time-scale plane for 2 gaussian atoms at different scales
sig=amgauss(256).*(fmconst(256,.15)+fmconst(256,.35));
[tfr,t,f]=tfrscalo(sig,1:256,12,.01,.49,512);
int1=integ2d(tfr(1:430,:),t,f(1:430));
int2=integ2d(tfr(431:512,:),t,f(431:512));
if abs(int1-int2)>1e-4,
 error('tfrscalo test 6 failed');
end

N=127;

% Covariance by dilation
t=ceil(N/2); f=0.2; T=2*sqrt(N); a=2; W=8;
sig1=amgauss(N,t,T).*fmconst(N,f,t);
sig2=amgauss(a*N,a*t,T*a).*fmconst(a*N,f/a,a*t);
[tfr1,t1,f1]=tfrscalo(sig1,1:N  ,W,0.01,0.49,N);
[tfr2,t2,f2]=tfrscalo(sig2,1:a*N,W,0.01,0.49,N);
Max1=max(max(tfr1)); Max2=max(max(tfr2));
[I1,J1]=find(tfr1==Max1); [I2,J2]=find(tfr2==Max2);
if abs(f1(I1(1))-a*f2(I2))>1e-2 | J2~=a*J1(1),
 error('tfrscalo test 7 failed');
end

% Reality of the TFR
sig=noisecg(N); W=5;
tfr=tfrscalo(sig,1:N,W,0.01,0.5,N);
if sum(any(abs(imag(tfr))>sqrt(eps)))~=0,
 error('tfrscalo test 8 failed');
end

% Energy conservation
sig=fmsin(N,.1,.4); W=6; Nf=2*N+1 ;
[tfr,t,f]=tfrscalo(sig,1:N,W,0.01,0.49,Nf);
SP = fft(hilbert(real(sig)));
indmin = 1+round(0.01*(N-2));
indmax = 1+round(0.49*(N-2));
SPana = SP(indmin:indmax);
Es=SPana'*SPana/Nf;
Etfr=integ2d(tfr,t,f);
if abs(Es-Etfr)>1e-2,
 error('tfrscalo test 9 failed');
end

% Positivity
if any(any(tfr<0)),
 error('tfrscalo test 10 failed');
end

% Same energy in the time-scale plane for 2 gaussian atoms at different scales
sig=amgauss(N).*(fmconst(N,.15)+fmconst(N,.35));
[tfr,t,f]=tfrscalo(sig,1:N,12,.01,.49,2*N+1);
int1=integ2d(tfr(1:210,:),t,f(1:210));
int2=integ2d(tfr(211:255,:),t,f(211:255));
if abs(int1-int2)>5e-4,
 error('tfrscalo test 11 failed');
end

%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

author: HeLiangHIT


text:
%CALCOPTICALFLOWFARNEBACK  Computes a dense optical flow using the Gunnar Farneback's algorithm
%
%     flow = cv.calcOpticalFlowFarneback(prevImg, nextImg)
%     flow = cv.calcOpticalFlowFarneback(prevImg, nextImg, 'OptionName',optionValue, ...)
%
% ## Input
% * __prevImg__ First 8-bit single-channel input image.
% * __nextImg__ Second input image of the same size and the same type as
%   `prevImg`.
%
% ## Output
% * __flow__ Computed flow image that has the same size as `prevImg` and
%   `single` type (2-channels). Flow for `(x,y)` is stored in the third
%   dimension.
%
% ## Options
% * __InitialFlow__ Initial flow approximation. Not set by default.
% * __PyrScale__ Parameter specifying the image scale (`<1`) to build pyramids
%   for each image. `PyrScale=0.5` means a classical pyramid, where each next
%   layer is twice smaller than the previous one. default 0.5.
% * __Levels__ Number of pyramid layers including the initial image.
%   `Levels=1` means that no extra layers are created and only the original
%   images are used. default 5.
% * __WinSize__ Averaging window size. Larger values increase the algorithm
%   robustness to image noise and give more chances for fast motion detection,
%   but yield more blurred motion field. default 13.
% * __Iterations__ Number of iterations the algorithm does at each pyramid
%   level. default 10.
% * __PolyN__ Size of the pixel neighborhood used to find polynomial expansion
%   in each pixel. Larger values mean that the image will be approximated with
%   smoother surfaces, yielding more robust algorithm and more blurred motion
%   field. Typically, `PolyN` is 5 or 7. default 5.
% * __PolySigma__ Standard deviation of the Gaussian that is used to smooth
%   derivatives used as a basis for the polynomial expansion. For `PolyN=5`,
%   you can set `PolySigma = 1.1`. For `PolyN=7`, a good value would be
%   `PolySigma = 1.5`. default 1.1.
% * __Gaussian__ Use the Gaussian `WinSize x WinSize` filter instead of a box
%   filter of the same size for optical flow estimation. Usually, this option
%   gives z more accurate flow than with a box filter, at the cost of lower
%   speed. Normally, `WinSize` for a Gaussian window should be set to a larger
%   value to achieve the same level of robustness. default false.
%
% The function finds an optical flow for each `prevImg` pixel using the
% [Farneback2003] alorithm so that:
%
%     prevImg(y,x) ~ nextImg(y + flow(y,x,2), x + flow(y,x,1))
%
% ## References
% [Farneback2003]:
% > Gunnar Farneback. "Two-frame motion estimation based on polynomial
% > expansion". In Image Analysis, pages 363-370. Springer, 2003.
%
% See also: cv.FarnebackOpticalFlow, cv.calcOpticalFlowPyrLK, opticalFlowLK,
%  opticalFlowHS, opticalFlowFarneback, vision.OpticalFlow, vision.BlockMatcher
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function sS = som_vs2to1(sS)

%SOM_VS2TO1 Convert version 2 struct to version 1.
%
% sSold = som_vs2to1(sSnew)
%
%  sMold = som_vs2to1(sMnew);
%  sDold = som_vs2to1(sDnew);
%
%  Input and output arguments:
%   sSnew   (struct) a SOM Toolbox version 2 struct
%   sSold   (struct) a SOM Toolbox version 1 struct
%
% For more help, try 'type som_vs2to1' or check out online documentation.
% See also  SOM_SET, SOM_VS1TO2.

%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% som_vs2to1
%
% PURPOSE
%
% Converts SOM Toolbox version 2 structs to version 1 structs.
%
% SYNTAX
%
%  sS1 = som_vs2to1(sS2)
%
% DESCRIPTION
%
% This function is offered to allow the change of new map and data structs
% to old ones. There are quite a lot of changes between the versions,
% especially in the map struct, and this function makes it possible to
% use the old functions with new structs.
%
% Note that part of the information is lost in the conversion. Especially,
% training history is lost, and the normalization is, except in the simplest
% cases (like all have 'range' or 'var' normalization) screwed up.
%
% REQUIRED INPUT ARGUMENTS
%
%  sS2       (struct) som SOM Toolbox version 2.0 struct (map, data,
%                     training or normalization struct)
%
% OUTPUT ARGUMENTS
%
%  sS1       (struct) the corresponding SOM Toolbox version 2.0 struct
%
% EXAMPLES
%
%  sM = som_vs2to1(sMnew);
%  sD = som_vs2to1(sDnew);
%  sT = som_vs2to1(sMnew.trainhist(1));
%
% SEE ALSO
%
%  som_set          Set values and create SOM Toolbox structs.
%  som_vs1to2       Transform structs from 1.0 version to 2.0.

% Copyright (c) 1999-2000 by the SOM toolbox programming team.
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta juuso 101199

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% check arguments

error(nargchk(1, 1, nargin));   % check no. of input arguments is correct

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% set field values

switch sS.type,
 case 'som_map',
  msize = sS.topol.msize;
  [munits dim] = size(sS.codebook);

  % topology
  if strcmp(sS.topol.shape,'sheet'), shape = 'rect';
  else shape = sS.shape;
  end

  % labels
  labels = cell(munits,1);
  nl = size(sS.labels,2);
  for i=1:munits,
    labels{i} = cell(nl,1);
    for j=1:nl, labels{i}{j} = sS.labels{i,j}; end
  end

  % trainhist
  tl = length(sS.trainhist);
  if tl==0 || strcmp(sS.trainhist(1).algorithm,'lininit'),
    init_type = 'linear';
  else
    init_type = 'random';
  end
  if tl>1,
    for i=2:tl,
      train_seq{i-1} = som_vs2to1(sS.trainhist(i));
    end
    train_type = sS.trainhist(tl).algorithm;
  else
    train_seq = [];
    train_type = 'batch';
  end
  if tl>0, data_name = sS.trainhist(tl).data_name; else data_name = ''; end

  % component normalizations
  sN = convert_normalizations(sS.comp_norm);
  if strcmp(sN.name,'som_hist_norm'),
    sS.codebook = redo_hist_norm(sS.codebook,sS.comp_norm,sN);
  end

  % map
  sSnew = struct('init_type', 'linear', 'train_type', 'batch', 'lattice' ,...
		 'hexa', 'shape', 'rect', 'neigh', 'gaussian', 'msize', msize, ...
		 'train_sequence', [], 'codebook', [], 'labels', [], ...
		 'mask', [], 'data_name', 'unnamed', 'normalization', [], ...
		 'comp_names', [], 'name', 'unnamed');
  sSnew.init_type = init_type;
  sSnew.train_type = train_type;
  sSnew.lattice = sS.topol.lattice;
  sSnew.shape = shape;
  sSnew.neigh = sS.neigh;
  sSnew.msize = sS.topol.msize;
  sSnew.train_sequence = train_seq;
  sSnew.codebook = reshape(sS.codebook,[sS.topol.msize dim]);
  sSnew.labels = labels;
  sSnew.mask = sS.mask;
  sSnew.data_name = data_name;
  sSnew.normalization = sN;
  sSnew.comp_names = sS.comp_names;
  sSnew.name = sS.name;

 case 'som_data',
  [dlen dim] = size(sS.data);

  % component normalizations
  sN = convert_normalizations(sS.comp_norm);
  if strcmp(sN.name,'som_hist_norm'),
    sS.codebook = redo_hist_norm(sS.codebook,sS.comp_norm,sN);
  end

  % data
  sSnew = struct('data', [], 'name', '', 'labels' , [], 'comp_names', ...
		 [], 'normalization', []);
  sSnew.data = sS.data;
  sSnew.name = sS.name;
  sSnew.labels = sS.labels;
  sSnew.comp_names = sS.comp_names;
  sSnew.normalization = sN;

 case 'som_norm',
  sSnew = struct('name','som_var_norm','inv_params',[]);

  switch sS.method,
   case 'var',   sSnew.name = 'som_var_norm';
   case 'range', sSnew.name = 'som_lin_norm';
   case 'histD', sSnew.name = 'som_hist_norm';
   otherwise,
    warning(['Method ' method ' does not exist in version 1.'])
  end

  if strcmp(sS.status,'done'),
    switch sS.method,
     case 'var',
      sSnew.inv_params = zeros(2,1);
      sSnew.inv_params(1) = sS.params(1);
      sSnew.inv_params(2) = sS.params(2);
     case 'range',
      sSnew.inv_params = zeros(2,1);
      sSnew.inv_params(1) = sS.params(1);
      sSnew.inv_params(2) = sS.params(2) + sS.params(1);
     case 'histD',
      bins = length(sS.params);
      sSnew.inv_params = zeros(bins+1,1) + Inf;
      sSnew.inv_params(1:bins,i) = sS.params;
      sSnew.inv_params(end,i) = bins;
    end
  end

 case 'som_train',
  sSnew = struct('algorithm', sS.algorithm, 'radius_ini', ...
		 sS.radius_ini, 'radius_fin', sS.radius_fin, 'alpha_ini', ...
		 sS.alpha_ini, 'alpha_type', sS.alpha_type, 'trainlen', sS.trainlen, ...
		 'qerror', NaN, 'time', sS.time);

 case 'som_topol',
  disp('Version 1 of SOM Toolbox did not have topology structure.\n');

 otherwise,

  error('Unrecognized struct.');
end

sS = sSnew;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% subfunctions

function sN = convert_normalizations(cnorm)

  dim = length(cnorm);
  sN = struct('name','som_var_norm','inv_params',[]);

  % check that there is exactly one normalization per component
  % and that their status and method is the same
  ok = 1;
  nof = zeros(dim,1);
  for i=1:dim, nof(i) = length(cnorm{i}); end
  if any(nof>1), ok=0;
  elseif any(nof==1) && any(nof==0), ok=0;
  elseif any(nof>0),
    status = cnorm{1}.status;
    method = cnorm{1}.method;
    for i=2:dim,
      if ~strcmp(cnorm{i}.status,status) || ~strcmp(cnorm{i}.method,method),
	ok = 0;
      end
    end
  elseif all(nof==0),
    return;
  end
  if ~ok,
    warning(['Normalization could not be converted. All variables can' ...
	     ' only be normalized with a single, and same, method.']);
    return;
  end

  % method name
  switch method,
   case 'var', sN.name = 'som_var_norm';
   case 'range', sN.name = 'som_lin_norm';
   case 'histD', sN.name = 'som_hist_norm';
   otherwise,
    warning(['Normalization could not be converted. Method ' method ...
	     'does not exist in version 1.']);
    return;
  end

  % if not done, inv_params is empty
  if ~strcmp(status,'done'), return; end

  % ok, make the conversion
  switch method,
   case 'var',
    sN.inv_params = zeros(2,dim);
    for i=1:dim,
      sN.inv_params(1,i) = cnorm{i}.params(1);
      sN.inv_params(2,i) = cnorm{i}.params(2);
    end
   case 'range',
    sN.inv_params = zeros(2,dim);
    for i=1:dim,
      sN.inv_params(1,i) = cnorm{i}.params(1);
      sN.inv_params(2,i) = cnorm{i}.params(2) + cnorm{i}.params(1);
    end
   case 'histD',
    bins = zeros(dim,1);
    for i=1:dim, bins(i) = length(cnorm{i}.params); end
    m = max(bins);
    sN.inv_params = zeros(m+1,dim) + Inf;
    for i=1:dim,
      sN.inv_params(1:bins(i),i) = cnorm{i}.params;
      if bins(i)<m, sN.inv_params(bins(i)+1,i) = NaN; end
      sN.inv_params(end,i) = bins(i);
    end
  end

function D = redo_hist_norm(D,cnorm,sN)

  dim = size(D,2);

  % first - undo the new way
  for i=1:dim,
    bins = length(cnorm{i}.params);
    D(:,i) = round(D(:,i)*(bins-1)+1);
    inds = find(~isnan(D(:,i)) & ~isinf(D(:,i)));
    D(inds,i) = cnorm{i}.params(D(inds,i));
  end
  % then - redo the old way
  n_bins = sN.inv_params(size(sN.inv_params,1),:);
  for j = 1:dim,
    for i = 1:size(D, 1)
      if ~isnan(D(i, j)),
	[d ind] = min(abs(D(i, j) - sN.inv_params(1:n_bins(j), j)));
	if (D(i, j) - sN.inv_params(ind, j)) > 0 && ind < n_bins(j),
	  D(i, j) = ind + 1;
	else
	  D(i, j) = ind;
	end
      end
    end
  end
  D = D * sparse(diag(1 ./ n_bins));

author: ilarinieminen

repo: SOM-Toolbox


text:
classdef VademecumTxiRhoPlotter < VademecumPlotter

    properties (SetAccess = private, GetAccess = public)
        feasibleIndex
    end

    properties (Access = protected)
       name = 'TxiRho';
    end

    properties (Access = private)
        mxT
        myT
        chi
    end

    methods (Access = public)

        function obj = VademecumTxiRhoPlotter(d)
            obj.init(d);
            obj.computeTxiMxMyVariables();
            obj.computeFeasibleIndex();
            obj.feasibleIndex = 1:length(obj.mxV)*length(obj.myV);
            ind = obj.feasibleIndex;
            obj.xV = obj.chi(ind);
            obj.yV = obj.volume(ind);
        end

        function plot(obj)
            obj.plotMxMy();
            obj.plotHomogenizedTensor();
            obj.plotHomogenizedTensorIsotropy();
            obj.plotAmplificatorTensor();
        end

    end

    methods (Access = private)

        function computeTxiMxMyVariables(obj)
            for i = 1:length(obj.mxV)
                for j = 1:length(obj.myV)
                    mx = obj.mxV(i);
                    my = obj.myV(j);
                    obj.chi(i,j) = atan(mx/my);
                    obj.mxT(i,j) = mx;
                    obj.myT(i,j) = my;
                end
            end
        end

        function computeFeasibleIndex(obj)
            d.mx = obj.mxV;
            d.my = obj.myV;
            d.chi = obj.chi;
            d.rho = obj.volume;
            fC = FeasibleIndexComputer(d);
            obj.feasibleIndex = fC.index;
        end

        function plotMxMy(obj)
            obj.printWidthVariable(obj.mxT,'m_1');
            obj.printWidthVariable(obj.myT,'m_2');
        end

        function printWidthVariable(obj,val,name)
            obj.fileName    = name;
            obj.titleName   = name;
            obj.value2print = val;
            obj.plotFigure();
            obj.printFigure();
        end

    end

    methods (Access = protected)

        function plotFigure(obj)
            x(:,1) = obj.xV;
            y(:,1) = obj.yV;
            z(:,1) = obj.value2print(obj.feasibleIndex);

            s.fileName = fullfile(obj.outPutPath,[obj.fileName,'XiRho']);
            s.title    = obj.titleName;
            s.axisAdder = XiRhoAxisAdder();
            p  = SuperEllipseExponentContourPlotter(s);
            p.plot(x,y,z);
%
%             ncolors = 50;
%             tri = delaunay(x,y);
%             obj.fig = figure;
%             tricontour(tri,x,y,z,ncolors)
%             colorbar
%             hold on
%             plot(x,y,'+');
%             ylim([0 1])
%             xlabel('$\xi$','Interpreter','latex');
%             ylabel('\rho');
%             tN = obj.titleName;
%             title(['$',tN,'$'],'interpreter','latex')
%             set(gca,'xtick',[0:pi/8:pi/2]) % where to set the tick marks
%             set(gca,'xticklabels',{'0','\pi/8','\pi/4','3\pi/8','\pi/2'})
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan