text:
classdef ConvTranspose < dagnn.Layer
  properties
    size = [0 0 0 0]
    hasBias = true
    upsample = [1 1]
    crop = [0 0 0 0]
    numGroups = 1
    opts = {'cuDNN'}
  end

  methods
    function outputs = forward(obj, inputs, params)
      if ~obj.hasBias, params{2} = [] ; end
      outputs{1} = vl_nnconvt(...
        inputs{1}, params{1}, params{2}, ...
        'upsample', obj.upsample, ...
        'crop', obj.crop, ...
        'numGroups', obj.numGroups, ...
        obj.opts{:}) ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      if ~obj.hasBias, params{2} = [] ; end
      [derInputs{1}, derParams{1}, derParams{2}] = vl_nnconvt(...
        inputs{1}, params{1}, params{2}, derOutputs{1}, ...
        'upsample', obj.upsample, ...
        'crop', obj.crop, ...
        'numGroups', obj.numGroups, ...
        obj.opts{:}) ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      outputSizes{1} = [...
        obj.upsample(1) * (inputSizes{1}(1) - 1) + obj.size(1) - obj.crop(1) - obj.crop(2), ...
        obj.upsample(2) * (inputSizes{1}(2) - 1) + obj.size(2) - obj.crop(3) - obj.crop(4), ...
        obj.size(3), ...
        inputSizes{1}(4)] ;
    end

    function rfs = getReceptiveFields(obj)
      rfs.size = (obj.size(1:2) - 1) ./ obj.upsample + 1 ;
      rfs.stride = 1 ./ [obj.upsample] ;
      rfs.offset = (2*obj.crop([1 3]) - obj.size(1:2) + 1) ...
        ./ (2*obj.upsample) + 1 ;
    end

    function params = initParams(obj)
      % todo: test this initialization method
      sc = sqrt(2 / prod(obj.size([1 2 4]))) ;
      params{1} = randn(obj.size,'single') * sc ;
      if obj.hasBias
        params{2} = zeros(obj.size(3),1,'single') * sc ;
      end
    end

    function set.size(obj, ksize)
      % make sure that ksize has 4 dimensions
      ksize = [ksize(:)' 1 1 1 1] ;
      obj.size = ksize(1:4) ;
    end

    function set.crop(obj, crop)
      if numel(crop) == 1
        obj.crop = [crop crop crop crop] ;
      elseif numel(crop) == 2
        obj.crop = crop([1 1 2 2]) ;
      else
        obj.crop = crop ;
      end
    end

    function set.upsample(obj, upsample)
      if numel(upsample) == 1
        obj.upsample = [upsample upsample] ;
      else
        obj.upsample = upsample ;
      end
    end

    function obj = ConvTranspose(varargin)
      obj.load(varargin) ;
      % normalize field by implicitly calling setters defined in
      % dagnn.Filter and here
      obj.size = obj.size ;
      obj.upsample = obj.upsample ;
      obj.crop = obj.crop ;
    end
  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


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

%script to plot the phase plane for a bang-bang example
figure(2)
hold off
clf
N=-1;
x = -100;
for k=7:11;
xdot=2.02*k;
sim('bang');
plot(xbang(:,2),xdotbang(:,2),'b-',xbang(:,2),-xdotbang(:,2),'b-');
xlabel('x_1');
ylabel('x_2');
text(-120,19,'u=-1');
hold on;
end;
N=1;
x = 100;
for k=7:11;
xdot=-2.02*k;
sim('bang')
plot(xbang(:,2),xdotbang(:,2),'r-',xbang(:,2),-xdotbang(:,2),'r-');
hold on;
end
text(100,-17,'u=+1');
title('Switching curves for 1/s^2 plant')
grid on

author: Sable

repo: mcbench-benchmarks


text:
function test_tutorial_plotting

% WALLTIME 00:15:00
% MEM 4gb
% DEPENDENCY ft_multiplotER ft_singleplotER ft_topoplotER ft_singleplotTFR ft_multiplotTRF ft_megplanar ft_combineplanar ft_volumereslice

%%
% *Plotting data at the channel and source level*
%
% _*Introduction*_
%
% To visualize your data you can use FieldTrip's high-level plotting functions,
% which are optimised for the FT data structures. Like other high-level functions
% they take a cfg as first input argument, followed by the data. These functions
% take care of all bookkeeping and some of the functions allow you to interact
% with the data by clicking in the figure.
%
% Alternatively, you can use the low-level FieldTrip plotting functions.
% These are located in the plotting sub-directory and are named ft_plot_xxx. You
% can find them by typing in the command window ?help ft_plot_? and then press
% the Tab key.
%
% Of course you can use the standard MATLAB functions ('plot, 'image', 'images',
% 'patch').
%
% _*Background*_
%
% The goal of the plotting functions in FieldTrip is to ease the visualization
% of complex data structures with multiple dimensions and with data that is not
% trivial to interpret spatially. This is accomplished with high-level functions
% (e.g. ft_topoplotER or ft_multiplotTFR) and low-level functions (with prefixes
% 'ft_plot_' and 'ft_select_', e.g. ft_plot_matrix or ft_select_box). For more
% simple data, such as a set of reaction times of the subject, we expect you to
% use the standard MATLAB plotting functions.
%
% The high-level functions take care of the data bookkeeping and call the
% low-level function. If you want to make more complex figures or tweak all options,
% you can bypass the high-level functions and call the low-level functions instead.
% This is especially useful if there is no data selection and bookkeeping involved,
% e.g. when you want to plot multiple geometrical objects (like sensors, source
% model, head surface, etc).
%
% To determine which high-level functions are suitable for you depends on
% the type of data you have: sensor or source space data. In this tutorial we
% assume that you already have the data from the event related averaging tutorial,
% the time-frequency representations of power tutorial and the applying beamforming
% techniques in the frequency domain tutorial, and we will demonstrate plotting
% at both the sensor and source level.
%
% Let us start by loading some data.

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

load avgFC.mat
load GA_FC.mat
load TFRhann.mat
load statERF.mat
load statTFR.mat
load sourceDiff.mat

%%
% _*Plotting data at the channel level*_
%
% Data at the channel level has a value for each sensor (MEG) or electrode
% (EEG). Those values can change over time and/or over frequency.
%
% *Singleplot functions*
%
% With ft_singleplotER you can make a plot using the avgFC data from the
% ERF tutorial by the following code:
%
% The ft_singleplotER function first selects the data to be plotted, in this
% case channel MLC24, from -0.2 to 1.0 seconds. Subsequently this selected data
% is plotted with the MATLAB plot function.

cfg = [];
cfg.xlim = [-0.2 1.0];
cfg.ylim = [-1e-13 3e-13];
cfg.channel = 'MLC24';
figure; ft_singleplotER(cfg,avgFC);
%%
% You could make the same plot by the following code:

selected_data = avgFC.avg(9,241:601); % MLC24 is the 9th channel, -0.2 to 1.0 is sample 241 to 601
selected_time = avgFC.time(241:601);
figure;
plot(selected_time, selected_data)
xlim([-0.2 1.0])
ylim([-1e-13 3e-13])
%%
% The advantage of ft_singleplotER is that it does the bookeeping for you.
% E.g. you can plot the mean over all MLC channels, which are MEG channels over
% left-central regions:

cfg = [];
cfg.xlim = [-0.2 1.0];
cfg.ylim = [-1e-13 3e-13];
cfg.channel = 'MLC*';
figure; ft_singleplotER(cfg,avgFC);

%%
% and change the line width or color:

cfg = [];
cfg.xlim = [-0.2 1.0];
cfg.ylim = [-1e-13 3e-13];
cfg.channel = 'MLC24';
cfg.linewidth = 10;
cfg.graphcolor = 'g';
figure; ft_singleplotER(cfg,avgFC);

%%
% The ft_singleplotER function might not look so impressive compared to
% standard plot; that is because the ERF data representation is quite simple.
%
% The advantage of using the standard MATLAB plot function is that you can
% easily find all documentation for it ("doc plot") and tweak the low level figure
% characteristics.
%
% The FieldTrip plotting functions have a lot of built-in intelligence to
% make plotting of multidimensional data easier. It is for instance possible to
% do baseline correction before plotting, by specifying the baseline type and
% time limits. In the plotting functions either the FieldTrip function ft_timelockbaseline
% or ft_freqbaseline is called. If you specify multiple channels in cfg.channel
% both singleplot functions will plot the mean over these channels. In the plotting
% functions the FieldTrip function ft_channelselection is called, which makes
% it straightforward to plot for instance the mean TFR (download here, see time-frequency
% analysis tutorial)of the left central channels.

cfg = [];
cfg.baseline = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.zlim = [-1.5e-27 1.5e-27];
cfg.channel = 'MLC24'; % top figure
figure; ft_singleplotTFR(cfg, TFRhann);

%%
% or wiht the average over a group of channels:

cfg.channel  = 'MLC*' % bottom figure
figure; ft_singleplotTFR(cfg,TFRhann);
%%
% You can again make very similar figures using standard MATLAB plotting
% functions.

chanindx = find(strcmp(TFRhann.label, 'MLC24'));
figure; imagesc(squeeze(TFRhann.powspctrm(chanindx,:,:)));
%%
% The default MATLAB image is shown upside down compared to usual MEG/EEG
% conventions. Along the axes you can see the number of time and frequency bins,
% but not the acqual time and frequency. This can be fixed of course.

figure; imagesc(TFRhann.time, TFRhann. freq, squeeze(TFRhann.powspctrm(chanindx,:,:)));
axis xy % flip vertically
colorbar
%%
% Still, the power is more difficult to interpret due to it being larger
% for low frequencies and smaller for high frequencies. This 1/f characteristic
% is removed from the figure by baseline correction. See also ft_freqbaseline.
%
% *Multiplot functions*
%
% The multiplot functions work similarly to the singleplot functions, again
% first by selecting the data and subsequently using the MATLAB functions plot
% and imagesc. But instead of one plot, multiple plots are made; one for each
% channel. These plots are arranged according to a specified layout in one pair
% of axes. In the subsequent figures you can see these axes that are normally
% set to ?off?.

cfg = [];
cfg.layout = 'CTF151.lay';
figure; ft_multiplotER(cfg, avgFC);
axis on % this shows the actual MATLAB axes
%%
% Normally the axes of the figure are not visible, only the ?axis? of each
% channel, but remember these are not real axes on which you can use MATLAB axis
% commands, the are just lines drawn by the function. Of course you can set the
% limits of the channel ?axis? by the cfg structure (cfg.xlim, cfg.ylim).

cfg = [];
cfg.baseline = [-0.5 -0.1];
cfg.zlim = [-3e-27 3e-27];
cfg.baselinetype = 'absolute';
cfg.layout = 'CTF151.lay';
figure; ft_multiplotTFR(cfg,TFRhann);
%%
% In the upper left corner you see the "COMMENT" element and in the upper
% right cornier you see the "SCALE" element. The scale contains the average of
% all channels.
%
% The location of all channels and the COMMENT and SCALE placeholder can
% be visualized with

cfg = [];
cfg.layout = 'CTF151.lay';
cfg.layout = ft_prepare_layout(cfg);
figure; ft_plot_layout(cfg.layout);
%%
% Note that the layout contains all 151 MEG channels; the one channel missing
% in the data is automatically excluded from the plots.
%
% All channels are squeezed in a circle and the nose and ears are indicated
% at the top and sides. The actual 3-D layout of the MEG channels in the CTF system
% is more complex than the circular arrangement suggested by the flat 2-D layout.
% This is a generic challenge for projections onto a 2-D plane and also applies
% to cartography. We can plot the bottom coil of the CTF gradiometer channels
% with

figure;
ft_plot_sens(avgFC.grad, 'chantype', 'meggrad')
ft_plot_axes(avgFC.grad)
view(-90, 90) % you should rotate the figure to get a good 3-D feel of the channel locations.
%%
% Although many tutorials elsewhere on the Fieldtrip wiki are using the
% standard CTF151 layout, we do have another one that better reflects the helmet
% shape:

cfg = [];
cfg.layout = 'CTF151_helmet.mat'; % note that it is a binary *.mat file, not an ascii *.lay file
cfg.layout = ft_prepare_layout(cfg);
figure; ft_plot_layout(cfg.layout);
%%
% The layout is determined by the layout file. Read more on layout files
% [[tutorial:layout|here]], and in the [[faq:what_is_the_format_of_the_layout_file_which_is_used_for_plotting|frequently
% asked questions]].

%%
% For multiplotting planar gradient data from the Elekta Neuromag system
% it is especially relevant to understand the layout files. The Neuromag system
% has two planar gradiometers and one axial magnetometer at each sensor location.
% You do not want to plot those on top of each other. Hence the Neuromag layout
% files contain two (for the old 122 channel system) or three (for the recent
% 306 channel systems) separate subplots for each channel location. Those two
% (or three) subplots hold the data for the two planar gradients and for the magnetometer
% signal.

cfg = [];
cfg.layout = 'neuromag306all.lay';
cfg.layout = ft_prepare_layout(cfg);
figure; ft_plot_layout(cfg.layout);
%%
% You should zoom in on the figure to see the triplets for the three channels
% ar each sensor locations. There are also template layouts available that only
% contain the magnetometer channels at the correct location (i.e. not shifted)
% and that contain the combined planar channels. The latter one can be used after
% applying ft_combineplanar on your data.
%
% *==== Topoplot functions ====*
%
%
%
% Ft_topoplotER and ft_topoplotTFR plot the topographic distribution of 2-Dimensional
% or 3-Dimensional datatypes as a 2-D circular view (looking down at the top of
% the head). The arrangement of the channels is again specified in the layout
% (see above in multiplot functions). The ft_topoplotER and ft_topoplotTFR functions
% first again select the data to be plotted from the 2D or 3D input data and subsequently
% plot the selected data using low-level Fieldtrip functions. Using one value
% for each channel and the x and y coordinates, the values between points are
% interpolated and plotted.

cfg = [];
cfg.xlim = [0.3 0.5];
cfg.zlim = [0 6e-14];
cfg.layout = 'CTF151.lay';
cfg.parameter = 'individual'; % the default (avg) is not present
figure; ft_topoplotER(cfg,GA_FC); colorbar
%%
%

cfg = [];
cfg.xlim = [0.9 1.3];
cfg.ylim = [15 20];
cfg.zlim = [-1e-27 1e-27];
cfg.baseline = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.layout = 'CTF151.lay';
figure; ft_topoplotTFR(cfg,TFRhann); colorbar

%%
% You can compare this to the more realistic helmet-style layout, which
% better acknowledges the fact that the MEG channels extend below the ears and
% that the helmet has sort of an opening at the face.

cfg = [];
cfg.xlim = [0.9 1.3];
cfg.ylim = [15 20];
cfg.zlim = [-1e-27 1e-27]; % this is in T^2
cfg.baseline = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.layout = 'CTF151_helmet';
figure; ft_topoplotTFR(cfg,TFRhann); colorbar
%%
%   The baseline correction computes the power difference compared to the
% pre-stimuls baseline. We can also look at the power ratio

cfg = [];
cfg.xlim = [0.9 1.3];
cfg.ylim = [15 20];
cfg.zlim = [0.6 1.4]; % the value 1 means 100%, so this ranges from 60% to 140% of the baseline power
cfg.baseline = [-0.5 -0.1];
cfg.baselinetype = 'relative';
cfg.layout = 'CTF151_helmet';
figure; ft_topoplotTFR(cfg,TFRhann); colorbar
%%
% The interpretation becomes more clear with

cfg = [];
cfg.xlim = [0.9 1.3];
cfg.ylim = [15 20];
cfg.zlim = [0.6 1.4]; % the value 1 means 100%, so this ranges from 60% to 140% of the baseline power
cfg.baseline = [-0.5 -0.1];
cfg.baselinetype = 'relative';
cfg.layout = 'CTF151_helmet';
figure; ft_topoplotTFR(cfg,TFRhann); colorbar

% so far it was the same as above, now change the colormap
ft_hastoolbox('brewermap', 1);         % ensure this toolbox is on the path
colormap(flipud(brewermap(64,'RdBu'))) % change the colormap
%%
% See <http://colorbrewer2.org http://colorbrewer2.org> for more details.
%
% In the help of ft_topoplotER and ft_topoplotTFR you can find many cfg options.
% For instance by specifying the cfg.xlim as a vector the functions makes selections
% of multiple time-windows and plots them as a series of subplots.

% for the multiple plots also
cfg = [];
cfg.xlim = [-0.4:0.2:1.4];
cfg.ylim = [15 20];
cfg.zlim = [-1e-27 1e-27];
cfg.baseline = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.layout = 'CTF151_helmet.mat';
cfg.comment = 'xlim';
cfg.commentpos = 'title';
figure; ft_topoplotTFR(cfg,TFRhann);
%%
% In ft_topoplotER and ft_topoplotTFR, you can specify many options to fully
% control the appearance of the picture. Subsequently you can use the MATLAB print
% function to write the figure to a file. Preferred file formats are EPS for vector
% drawings that can be edited in Adobe Illustrator or in Canvas (using ?print
% -depsc?) or PNG for bitmaps (using ?print -dpng?).
%
% To make the EPS-files optimally suitable for Adobe Illustrator, use the
% command ?print -depsc -adobecs -painter?.
%
% Since MATLAB uses the 'painter' renderer to export in Illustrator format,
% this method allows one to export quite complex figures that otherwise would
% be exported as bitmaps. Note, however, that the 'painter' renderer also has
% certain limitations compared to the z-buffer and openGL renderers. (See also
% MATLAB help on selecting a renderer).
%
% Some examples of what you can do:

% options for data selection (used with any plotting function):
cfg = [];
cfg.xlim = [0.9 1.3];
cfg.ylim = [15 20];
cfg.zlim = [-1e-27 1e-27];
cfg.baseline = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.layout = 'CTF151.lay';
%%
% Options specific to the topoplot

cfg.gridscale = 300;
cfg.style = 'straight';
cfg.marker = 'labels';
figure; ft_topoplotTFR(cfg,TFRhann);
%%
%

cfg.gridscale = 300;
cfg.contournum = 10;
cfg.colormap = gray(10);
figure; ft_topoplotTFR(cfg,TFRhann);
%%
%

cfg.gridscale = 300;
cfg.contournum = 4;
cfg.colormap = spring(4);
cfg.markersymbol  = '.';
cfg.markersize = 12;
cfg.markercolor = [0 0.69 0.94];
figure; ft_topoplotTFR(cfg,TFRhann);
%%
%
%
% *==== Interactive mode ====*
%
% During the phase of the analysis in which you do a lot of data inspection,
% you can use the interactive modus to go from one plot to the other. You can
% for instance select a certain frequency and time range in a singleplot, to get
% the average over that range plotted in a topoplot. Or select a group of channels
% in a topoplot or a multiplot and get the average over those channels for the
% whole time and frequency range in a single plot.
%
%

%interactive
cfg = [];
cfg.baseline = [-0.5 -0.1];
cfg.zlim = [-3e-27 3e-27];
cfg.baselinetype = 'absolute';
cfg.layout = 'CTF151.lay';
cfg.interactive = 'yes';
figure; ft_multiplotTFR(cfg,TFRhann);
%%
%
%
% You should select channels using the mouse and click, then you get singleplot
% average over those channels. You should select a time-frequency window and click,
% then you get topoplot average over that time-frequency window.
%
%
%
% *==== Plotting clusters ====*
%
% Ft_clusterplot plots a series of topoplots that highlight the clusters
% from cluster-based permutation testing. The output "stat" is 2D data from **[[reference:ft_timelockstatistics]]**
% or **[[reference:ft_freqstatistics]]** with cfg.correctm='cluster'.
%
% The function automatically finds the clusters in the data which are smaller
% than the pre-specified threshold (cfg.alpha) and plots a series of topoplots
% with the data in "stat" field (are for instance t-values) and the sensors which
% are part of the cluster highlighted.
%
% Stat should be 2D, i.e. channels-by-time, or channels-by frequency. You
% cannot visualize channels-by-frequency-by-time, that case requires either averaging
% over time, or averaging over frequency.
%
% *== Clusters in timelocked data ==*

% clusterplot
cfg = [];
cfg.zlim = [-6 6]; % T-values
cfg.alpha = 0.05;
cfg.layout = 'CTF151.lay';
ft_clusterplot(cfg,statERF);
%%
%
%
% *== Clusters in time-frequency data ==*

% clusterplot
cfg = [];
cfg.zlim = [-5 5];
cfg.alpha = 0.05;
cfg.layout = 'CTF151.lay';
ft_clusterplot(cfg,statTFR);
%%
%
%
% *==== Plotting channel-level connectivity ====*
%
% This is possible using ft_connectivityplot,  ft_multiplotCC and ft_topoplotCC.
%
%
%
% *==== Plotting Principal or Independent Components (PCA/ICA)====*
%
% To plot PCA, ICA or other decompositions that result from ft_componentanalysis
% you can use ft_topoplotIC for the topographies and ft_databrowser for the topographies
% combined with the time series. For a viewer that displays the power spectrum,
% topography and variance over time of each component, see [[https://github.com/fieldtrip/fieldtrip/blob/master/contrib/misc/ft_icabrowser.m]].
%
%
%
% *===== Plotting data at the source level =====*
%
% With the ft_sourceplot function you can plot functional source reconstructed
% data. Data structures can be source estimates from ft_sourceanalysis or ft_sourcegrandaverage
% or statistical values from ft_sourcestatistics. The following sections will
% present the available options for source-level plotting, following the same
% structure as previous sections.
%
% In visualising source reconstructed data, you should consider that  there
% are two principled ways of representing the spatial dimension of source reconstructed
% data:
%
%   - On a regular, 3-dimensional grid, i.e. volumetric data like an MRI
% scan
%
%   - On a surface geometry, i.e. a cortical sheet.
%
% Besides the three dimensions that describe the spatial or geometrical aspect
% of source reconstructed brain activity, the data can have additional dimensions
% for time and/or frequency.
%
%
%
% *==== Volumetric data  ====*
%
% Source level data is considered volumetric if the locations at which activity
% is estimated are spaced in a regular 3-dimensional grid, like voxels in an MRI.
% FieldTrip offers multiple plotting options for volumetric data that can be specified
% in the cfg for ft_sourceplot|ft_sourceplot.
%
% You can (1) make multiple 2D axials slices throughout the brain on which
% the functional data is plotted, (2) create three slices in each of the three
% orthorgonal directions (axial, sagittal and coronal) on which you can click
% around to navigate the whole the brain (3) project the functional data onto
% a surface and make a 3-D rendering of that.
%
% Let us first provide the basic code to use one of these plotting methods,
% using the data from Subject01 from the beamformer tutorial.  Later on, we will
% give more details on the configuration options for ft_sourceplot.
%
% *=== Individual anatomical MRI (prior to spatial normalization) ===*
%
%

% load MRI, reslice it and interpolate the low-res functional data onto the high-res anatomy
mri = ft_read_mri('Subject01.mri');
mri = ft_volumereslice([], mri);

cfg            = [];
cfg.downsample = 2;
cfg.parameter  = 'pow';
sourceDiffInt  = ft_sourceinterpolate(cfg, sourceDiff , mri);
%%
%

% spatially normalize the anatomy and functional data to MNI coordinates
cfg = [];
cfg.nonlinear = 'no';
sourceDiffIntNorm = ft_volumenormalise(cfg, sourceDiffInt);
%%
%

% plot multiple 2D axial slices
cfg = [];
cfg.method        = 'slice';
cfg.funparameter  = 'pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0.0 1.2];
cfg.opacitylim    = [0.0 1.2];
cfg.opacitymap    = 'rampup';
ft_sourceplot(cfg, sourceDiffInt);

% idem, in MNI (or SPM) coordinates
ft_sourceplot(cfg, sourceDiffIntNorm);
%%
%

% plot 3 orthogonal slices
cfg = [];
cfg.method        = 'ortho';
cfg.funparameter  = 'pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0.0 1.2];
cfg.opacitylim    = [0.0 1.2];
cfg.opacitymap    = 'rampup';
ft_sourceplot(cfg, sourceDiffInt);

% idem, in MNI (or SPM) coordinates
ft_sourceplot(cfg, sourceDiffIntNorm);
%%
%
%
% The three essential cfg parameters are:
%
%   * cfg.anaparameter the anatomy parameter, specifying the anatomy to be
% plotted
%
%   * cfg.funparameter the functional parameter, specifying the functional
% data to be plotted
%
%   * cfg.maskparameter the mask parameter, specifying the parameter to be
% used to mask the functional data
%
% All data that is used in the figure (anatomical, functional, mask) must
% be interpolated and represented onto the same grid. See ft_sourceinterpolate
% and ft_volumenormalise.
%
% *=== the anatomical parameter ===*
%
% The anatomy can be read with **[[:reference:ft_read_mri]]**. The functional
% data can be interpolated onto the anatomy by **[[reference:ft_sourceinterpolate]]**.
% The anatomy is scaled between 0 and 1 and plotted in gray scale.
%
% *=== the functional parameter ===*
%
% The functional data is plotted in color optionally on top of the anatomy.
% The colors used can be determined by cfg.colormap (see MATLAB function COLORMAP).
% How the functional values are assigned to the colormap is determined by cfg.colorlim.
% It makes sense to plot for instance source data as functional parameter, but
% also statistical values (for instance T-values).
%
% *=== the masking parameter ===*
%
% You can control the opacity of the functional data by the mask parameter.
% Which values are plotted opaque and which transparent is determined by cfg.opacitymap
% and cfg.opacitylim (see MATLAB function ALPHA and ALPHAMAP). The opacity map
% determines the degree of opacity of the functional data going from opaque to
% transparent.  There are multiple ways to determine your opacity scale, as a
% user you can determine the opacity values for each and every single voxel (and
% as such, region of interest).  As such, the opacity limits determine how the
% opacity map is assigned to the values of the mask parameter.
%
% === Example 1: Plotting only positive values ===
%
% Imagine that your functional data has values ranging from -3 to 3. Here
% we plot only the positive values  (zeromax), using the scale whereby the strongest
% values are opaque, and the values close to zero are transparent:
%
% |  cfg.maskparameter = cfg.funparameter|
%
% |  cfg.colorlim      = [0 3] (or'zeromax')|
%
% |  cfg.opacitymap    = 'rampup'|
%
% |  cfg.opacitylim    = [0 3] (or 'zeromax')|
%
%
%
% === Example 2: Plotting high absolute values ===
%
% Suppose the functional data is the same as in example 1, but now we only
% wants to plot the high negative values and high positive values (use "maxabs"
% setting). We set these high absolute values to opaque, and the values around
% zero to transparent:
%
% |  cfg.maskparameter = cfg.funparameter|
%
% |  cfg.colorlim      = [-3 3] (or'maxabs')|
%
% |  cfg.opacitymap    = 'vdown'|
%
% |  cfg.opacitylim    = [-3 3] (or'maxabs')|
%
%
%
% === Example 3: Masking voxels outside values of interest ===
%
% Here, we make a field in the data with an opacity value for each voxel,
% and apply that as your mask. For instance if you only want to plot the values
% between 2 and 2.5 you can specify:
%
% |  data.mask         = (data.fun>2 & data.fun<2.5)|
%
% |  cfg.maskparameter = 'mask'|
%
%
%
% ==== Plotting on a brain surface ====
%
% === Scalar data per vertex ===
%
% The representation of source activity on a surface results from source
% estimation using a cortical sheet as the source model. The cortical sheet is
% represented as a triangulated surface and the activity is assigned to each of
% the vertices, i.e. corner points of the triangles. Besides doing the source
% estimation on a cortical sheet, it is also possible to interpolate the volumetric
% data (i.e. estimated on a regular 3-D grid) onto a cortical sheet for visualization.
%
% Scalar data (e.g., time-averaged activity, frequency-specific power estimates,
% statistics, etc.) can be plotted using the ft_plot_mesh function. Alternatively,
% volumetric data can also be rendered on a surface by projecting it to a surface
% geometry, using ft_sourceplot. An example of the latter is given below, where
% we use the same data as in the preceding section.
%
%
%
% === Project volumetric data to an MNI white-matter surface surface ===
%
%

cfg = [];
cfg.method         = 'surface';
cfg.funparameter   = 'pow';
cfg.maskparameter  = cfg.funparameter;
cfg.funcolorlim    = [0.0 1.2];
cfg.funcolormap    = 'jet';
cfg.opacitylim     = [0.0 1.2];
cfg.opacitymap     = 'rampup';
cfg.projmethod     = 'nearest';
cfg.surffile       = 'surface_white_both.mat'; % Cortical sheet from canonical MNI brain
cfg.surfdownsample = 10;  % downsample to speed up processing
ft_sourceplot(cfg, sourceDiffIntNorm);

view ([90 0])             % rotate the object in the view
%%
% If you enable the camera toolbar (menu option "view", "camera toolbar")
% you have more options for controlling the 3-dimensional rendering, e.g. change
% perspective, change the position of the light (used for relection and shadows).
%
%
%
% As with the channel-level multiplots, you can change the colormap that
% maps the functional values on the colors.

ft_sourceplot(cfg, sourceDiffIntNorm);
colormap jet

ft_sourceplot(cfg, sourceDiffIntNorm);
colormap hot

ft_sourceplot(cfg, sourceDiffIntNorm);
colormap cool

%%
% If you right-click with the mouse on the colobar in the figure, you can
% also also do this interactively.
%
% To  interactively explore higher-dimensional data (such as TFR data) on
% the surface, ft_sourceplot cannot be used. You can make selections of a specific
% time and/or frequency and make a series of sublplots.
%
% *=== Higher dimensional source-level data ===*
%
%
%
% Since source-level function data requires three dimensions for "space"
% and uses the color dimension for "strength", there is no graphical dimension
% in which the temporal evolution or the spectral distribution can be visualized
% in full detail. The ft_sourceplot function allows you with method=ortho to use
% teh lower right corner of the figure to show the timecourse or spectrum of the
% functional activity at a specific location (the one you click). The ft_sourcemovie
% function allows you to use time (as in a movie) to explore the changes of cortical
% activity over time and/or frequency.
%
%
%
% *==== Using external tools ====*
%
%
%
% Although MATLAB is a very flexible development and analysis environment,
% it is not super-fast in visualisation. Hence external visualisation tools are
% sometimes more useful for exploring your data. Volumetric and surface based
% data can be exported to standard file formats using **[[:reference:ft_sourcewrite]]**.
% Subsequently, you can use external tools such as
%
%
%
%   * [[http://www.mccauslandcenter.sc.edu/mricro/mricron/index.html|MRIcron]]
%
%   * [[http://www.openwalnut.org|OpenWallnut]]
%
%   * [[http://www.humanconnectome.org/software/connectome-workbench.html|Connectome
% Workbench]]
%
%   * https://github.com/piermorel/gramm
%
%
%
% *=====Suggested further reading=====*
%
% Plotting channel-level data in a 2-dimensional representation on your flat
% computer screen or on paper requires that the 3-dimensional channel positions
% are mapped or projected onto the 2-dimensional plane. The tutorial on [[layout|specifying
% the channel layout for plotting]] explains how this mapping is constructed.
%
%

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  MATLAB Code for                                                  %
%                                                                   %
%  Multi-Objective Particle Swarm Optimization (MOPSO)              %
%  Version 1.0 - Feb. 2011                                          %
%                                                                   %
%  According to:                                                    %
%  Carlos A. Coello Coello et al.,                                  %
%  "Handling Multiple Objectives with Particle Swarm Optimization," %
%  IEEE Transactions on Evolutionary Computation, Vol. 8, No. 3,    %
%  pp. 256-279, June 2004.                                          %
%                                                                   %
%  Developed Using MATLAB R2009b (Version 7.9)                      %
%                                                                   %
%  Programmed By: S. Mostapha Kalami Heris                          %
%                                                                   %
%         e-Mail: sm.kalami@gmail.com                               %
%                 kalami@ee.kntu.ac.ir                              %
%                                                                   %
%       Homepage: http://www.kalami.ir                              %
%                                                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function rep=DeleteFromRep(rep,EXTRA,gamma)

    if nargin<3
        gamma=1;
    end

    for k=1:EXTRA
        [occ_cell_index occ_cell_member_count]=GetOccupiedCells(rep);

        p=occ_cell_member_count.^gamma;
        p=p/sum(p);

        selected_cell_index=occ_cell_index(RouletteWheelSelection(p));

        GridIndices=[rep.GridIndex];

        selected_cell_members=find(GridIndices==selected_cell_index);

        n=numel(selected_cell_members);

        selected_memebr_index=randi([1 n]);

        j=selected_cell_members(selected_memebr_index);

        rep=[rep(1:j-1); rep(j+1:end)];
    end

end

author: sfvsfv

repo: Mathematical-modeling


text:
% prettyMap 2.0
%--------------------------------------------------------------------------
%
% plot with projection an area of the world
%
% POSSIBLE SINTAXES:
%   prettyMap(map);
%
%   prettyMap(map, shape);
%   prettyMap(map, projection);
%   prettyMap(map, lineCol);
%
%   prettyMap(map, phiGrid, lambdaGrid);
%   prettyMap(map, shape, projection);
%   prettyMap(map, shape, lineCol);
%   prettyMap(map, projection, shape);
%   prettyMap(map, projection, lineCol);
%
%   prettyMap(map, phiGrid, lambdaGrid, shape);
%   prettyMap(map, phiGrid, lambdaGrid, projection);
%   prettyMap(map, phiGrid, lambdaGrid, lineCol);
%
%   prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax)
%   prettyMap(map, phiMin, phiMax, shape, projection);
%   prettyMap(map, phiMin, phiMax, shape, lineCol);
%   prettyMap(map, phiMin, phiMax, projection, shape);
%   prettyMap(map, phiMin, phiMax, projection, lineCol);
%
%   prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, shepe);
%   prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, projection);
%   prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, lineCol);
%
%   prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax)
%   prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection);
%   prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, shape, lineCol);
%   prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape);
%   prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, projection, lineCol);
%
%   prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shepe);
%   prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection);
%   prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, lineCol);
%
%   prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection);
%   prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shape, lineCol);
%   prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape);
%   prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection, lineCol);
%
%   prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape, lineCol);
%   prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection, lineCol);
%
% EXAMPLE:
%   prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, 'Miller Cylindrical');
%
% INPUT:
%   map             matrix containing data of the whole world to be shown
%   phiGrid         array [degree]
%   lambdaGrid      array [degree]
%   phiMin          minimum latitude    [degree]
%   phiMax          maximum latitude    [degree]
%   lambdaMin       minimum longitude   [degree]
%   lambdaMax       maximum longitude   [degree]
%   projection      type of projection to be used "standard" is the default
%   shape           shapefile to load as coast (or country) contour
%                       - fill          fill coasts coarse
%                       - coast         only coasts coarse
%                       - 50m           1:50000000 scale country contours
%                       - 30m           1:30000000 scale country contours
%                       - 10m           1:10000000 scale country contours
%   lineCol         [1 1 1] array of RGB component to draw the contour lines
%
% DEFAULT VALUES:
%    projection = 'Lambert'
%
% AVAILABLE PROJECTION:
%    * Lambert
%      Stereographic
%      Orthographic
%      Azimuthal Equal-area
%      Azimuthal Equidistant
%      Gnomonic
%      Satellite
%      Albers Equal-Area Conic
%      Lambert Conformal Conic
%      Mercator
%    * Miller Cylindrical
%    * Equidistant Cylindrical (world map)
%      Oblique Mercator
%      Transverse Mercator
%      Sinusoidal
%      Gall-Peters
%      Hammer-Aitoff
%      Mollweide
%      Robinson
%    * UTM
%
% SEE ALSO:
%   mapPlot, mapPlot3D, quiver
%
% REQUIREMENTS:
%   M_Map: http://www.eos.ubc.ca/~rich/map.html
%   shape files with contours
%
% VERSION: 2.1
%
% CREDITS:
%   http://www.eos.ubc.ca/~rich/map.html
%
%   Andrea Gatti
%   DIIAR - Politecnico di Milano
%   2013-12-19
%
function h = prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape, lineCol)

%shape = 'coast';
%shape = 'fill';
%shape = '10m';
%shape = '30m';
%shape = '50m';

% lineCol = [0 0 0];
limitsOk = false;

% Manage opening a new figure;
tohold = false;
if length(findall(0,'Type','figure'))>=1
    if ishold
        clf;
        tohold = true;
    else
        figure;
    end
end

switch (nargin)
    case 1                                                                % prettyMap(map);
        shape = 'coast';
        lineCol = [0 0 0];
        projection = 'Miller Cylindrical';
        phiMin = 90;
        phiMax = -90;
        lambdaMin = -180;
        lambdaMax = 180;

        deltaPhi = (phiMax-phiMin)/size(map,1);
        deltaLambda = (lambdaMax-lambdaMin)/size(map,2);

        phiGrid    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
        lambdaGrid = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
    case 2
        shape = 'coast';
        lineCol = [0 0 0];
        projection = 'Miller Cylindrical';
        if (ischar(phiGrid))
            if (sum(strcmp(phiGrid,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}]))) % prettyMap(map, shape);
                shape = phiGrid;
            else                                                          % prettyMap(map, projection);
                projection = phiGrid;
            end
        elseif (length(phiGrid) == 3)                                     % prettyMap(map, lineCol);
            lineCol = phiGrid;
        end

        phiMin = 90;
        phiMax = -90;
        lambdaMin = -180;
        lambdaMax = 180;

        deltaPhi = (phiMax-phiMin)/size(map,1);
        deltaLambda = (lambdaMax-lambdaMin)/size(map,2);

        phiGrid    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
        lambdaGrid = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
    case 3
        shape = 'coast';
        lineCol = [0 0 0];
        if (ischar(phiGrid))
            projection = 'Miller Cylindrical';
            if (sum(strcmp(phiGrid,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))
                shape = phiGrid;
                if (ischar(lambdaGrid))
                    projection = lambdaGrid;                              % prettyMap(map, shape, projection);
                else
                    lineCol = lambdaGrid;                                 % prettyMap(map, shape, lineCol);
                end
            else
                projection = phiGrid;
                if (ischar(lambdaGrid))
                    shape = lambdaGrid;                                   % prettyMap(map, projection, shape);
                else
                    lineCol = lambdaGrid;                                 % prettyMap(map, projection, lineCol);
                end
            end
            phiMin = 90;
            phiMax = -90;
            lambdaMin = -180;
            lambdaMax = 180;

            deltaPhi = (phiMax-phiMin)/size(map,1);
            deltaLambda = (lambdaMax-lambdaMin)/size(map,2);

            phiGrid    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
            lambdaGrid = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
        else                                                              % prettyMap(map, phiGrid, lambdaGrid);
            projection = 'Miller Cylindrical';
            phiMin = max(phiGrid);
            phiMax = min(phiGrid);
            lambdaMin = min(lambdaGrid);
            lambdaMax = max(lambdaGrid);
        end
    case 4
        shape = 'coast';
        lineCol = [0 0 0];
        projection = 'Miller Cylindrical';
        if (ischar(phiMin))
            if (sum(strcmp(phiMin,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))  % prettyMap(map, phiGrid, lambdaGrid, shape);
                shape = phiMin;
            else                                                          % prettyMap(map, phiGrid, lambdaGrid, projection);
                projection = phiMin;
            end
        elseif (length(phiMin) == 3)                                      % prettyMap(map, phiGrid, lambdaGrid, lineCol);
            lineCol = phiMin;
        end

        phiMin = max(phiGrid);
        phiMax = min(phiGrid);
        lambdaMin = min(lambdaGrid);
        lambdaMax = max(lambdaGrid);
    case 5
        shape = 'coast';
        lineCol = [0 0 0];
        projection = 'Miller Cylindrical';
        if (ischar(phiMin))
            if (sum(strcmp(phiMin,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))
                shape = phiMin;
                if (ischar(phiMax))
                    projection = phiMax;                                  % prettyMap(map, phiMin, phiMax, shape, projection);
                else
                    lineCol = phiMax;                                     % prettyMap(map, phiMin, phiMax, shape, lineCol);
                end
            else
                projection = phiMin;
                if (ischar(phiMax))
                    shape = phiMax;                                       % prettyMap(map, phiMin, phiMax, projection, shape);
                else
                    lineCol = phiMax;                                     % prettyMap(map, phiMin, phiMax, projection, lineCol);
                end
            end
            phiMin = max(phiGrid);
            phiMax = min(phiGrid);
            lambdaMin = min(lambdaGrid);
            lambdaMax = max(lambdaGrid);
        else                                                             %  prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax);
            limitsOk = true;
            lambdaMin = phiMin;
            lambdaMax = phiMax;
            phiMin = phiGrid;
            phiMax = lambdaGrid;

            if (phiMin < phiMax)
                tmp = phiMin;
                phiMin = phiMax;
                phiMax = tmp;
            end

            deltaPhi = (phiMax-phiMin)/size(map,1);
            deltaLambda = (lambdaMax-lambdaMin)/size(map,2);

            phiGrid    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
            lambdaGrid = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
         end
    case 6
        shape = 'coast';
        lineCol = [0 0 0];
        limitsOk = true;
        projection = 'Lambert';
        if (ischar(lambdaMin))
            if (sum(strcmp(lambdaMin,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))
                shape = lambdaMin;                                        % prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, shape);
            else
                projection = lambdaMin;                                   % prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, projection);
            end
        elseif (length(lambdaMin) == 3)
            lineCol = lambdaMin;                                          % prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, lineCol);
        end

        lambdaMax = phiMax;
        lambdaMin = phiMin;
        phiMin = phiGrid;
        phiMax = lambdaGrid;

        if (phiMin < phiMax)
            tmp = phiMin;
            phiMin = phiMax;
            phiMax = tmp;
        end

        deltaPhi = (phiMax-phiMin)/size(map,1);
        deltaLambda = (lambdaMax-lambdaMin)/size(map,2);

        phiGrid    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
        lambdaGrid = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
    case 7
        shape = 'coast';
        lineCol = [0 0 0];
        limitsOk = true;
        projection = 'Lambert';
        if (ischar(lambdaMin))
            if (sum(strcmp(lambdaMin,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))
                shape = lambdaMin;
                if (ischar(lambdaMax))
                    projection = lambdaMax;                               % prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection);
                else
                    lineCol = lambdaMax;                                  % prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, shape, lineCol);
                end
            else
                projection = lambdaMin;
                if (ischar(lambdaMax))
                    shape = lambdaMax;                                    % prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape);
                else
                    lineCol = lambdaMax;                                  % prettyMap(map, phiMin, phiMax, lambdaMin, lambdaMax, projection, lineCol);
                end
            end

            lambdaMin = phiMin;
            lambdaMax = phiMax;
            phiMin = phiGrid;
            phiMax = lambdaGrid;

            if (phiMin < phiMax)
                tmp = phiMin;
                phiMin = phiMax;
                phiMax = tmp;
            end

            deltaPhi = (phiMax-phiMin)/size(map,1);
            deltaLambda = (lambdaMax-lambdaMin)/size(map,2);

            phiGrid    = (phiMin + deltaPhi/2 : deltaPhi : phiMax - deltaPhi/2)';
            lambdaGrid = (lambdaMin + deltaLambda/2 :  deltaLambda :  lambdaMax  - deltaLambda/2)';
        else
            projection = 'lambert';                                       % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax);
        end
    case 8
        shape = 'coast';
        lineCol = [0 0 0];
        limitsOk = true;
        if (ischar(projection))
            if (sum(strcmp(projection,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))
                shape = projection;                                       % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shape);
                projection = 'Lambert';
            else
                                                                          % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection);
            end
        elseif (length(projection) == 3)
            lineCol = projection;                                         % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, lineCol);
            projection = 'Lambert';
        end
    case 9
        lineCol = [0 0 0];
        limitsOk = true;
        if (ischar(projection))
            if (sum(strcmp(projection,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))
                tmp = shape;
                shape = projection;
                if (ischar(tmp))
                    projection = tmp;                                     % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection);
                else
                    lineCol = tmp;                                        % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shape, lineCol);
                    projection = 'UTM';
                end
            else
                if (ischar(shape))
                                                                          % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape);
                else
                    lineCol = shape;                                      % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection, lineCol);
                    shape = 'coast';
                end
            end
        end
    case 10                                                               % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, projection, shape, lineCol)
        limitsOk = true;
       if (sum(strcmp(projection,[{'none'},{'coast'},{'fill'},{'10m'},{'30m'},{'50m'}])))   % prettyMap(map, phiGrid, lambdaGrid, phiMin, phiMax, lambdaMin, lambdaMax, shape, projection, lineCol)
            tmp = shape;
            shape = projection;
            projection = tmp;
       end
end

if (phiMin < phiMax)
    tmp = phiMin;
    phiMin = phiMax;
    phiMax = tmp;
end

lambdaGrid = sort(lambdaGrid);

[val idMax] = max(diff(lambdaGrid));
if sum(diff(lambdaGrid) == val) == 1
    lambdaGrid(1:idMax) = lambdaGrid(1:idMax)+360;
    if ~limitsOk
        lambdaMax = lambdaGrid(idMax);
        lambdaMin = lambdaGrid(idMax+1);
    end
end
lambdaGrid = sort(lambdaGrid);

if(lambdaMax<lambdaMin)
    lambdaMax = lambdaMax+360;
end

% setup the projection
if (strcmpi(projection,[{'lambert'}]) && abs(phiMax==-phiMin))
    projection='Miller Cylindrical';
end

if (sum(strcmpi(projection,[{'lambert'},{'UTM'},{'Sinusoidal'},{'Transverse Mercator'},{'Mollweid'},{'Oblique Mercator'},{'Miller Cylindrical'}])))
    m_proj(projection,'long',[lambdaMin lambdaMax],'lat',[phiMax phiMin]);
else
    m_proj(projection);
end

% Printing projection
fprintf('Using projection: %s\n', projection);

if sum(diff(lambdaGrid)<-200)
    lambdaGrid(lambdaGrid<0)=lambdaGrid(lambdaGrid<0)+360;
end
% plot the map
h = m_pcolor(lambdaGrid, phiGrid, map);

% set the light
shading flat;

% read shapefile
if (~strcmp(shape,'none'))
	if (~strcmp(shape,'coast')) && (~strcmp(shape,'fill'))
    	if (strcmp(shape,'10m'))
        	M=m_shaperead('countries_10m');
	    elseif (strcmp(shape,'30m'))
    	    M=m_shaperead('countries_30m');
	    else
    	    M=m_shaperead('countries_50m');
	    end
    	[xMin,yMin] = m_ll2xy(lambdaMin,phiMin);
	    [xMax,yMax] = m_ll2xy(lambdaMax,phiMax);
    	for k=1:length(M.ncst)
        	lamC = M.ncst{k}(:,1);
	        ids = lamC < lambdaMin;
   	    	lamC(ids) = lamC(ids) + 360;
        	phiC = M.ncst{k}(:,2);
        	[x,y] = m_ll2xy(lamC,phiC);
        	if sum(~isnan(x))>1
            	x(find(abs(diff(x))>=abs(xMax-xMin)*0.90)+1) = nan; % Remove lines that occupy more than th 90% of the plot
	            line(x,y,'color', lineCol);
    	    end
        end
    else
        if (strcmp(shape,'coast'))
        	m_coast('line','color', lineCol);
        else
            m_coast('patch',lineCol);
        end
	end
end

m_grid('box','fancy','tickdir','in');
colorbar;

if tohold
    hold on;
end

author: goGPS-Project


text:
function [D,S,C] = signed_distance_direction(P,V,F)
  % SIGNED_DISTANCE_DIRECTION Compute a direction which decreases signed
  % distance with respect to an input mesh (V,F) at a set of points P.
  %
  % [D,S,C] = signed_distance_direction(P,V,F)
  %
  % Inputs:
  %   P  #P by dim list of query points
  %   V  #V by dim list of mesh vertex positions
  %   F  #F by dim+1 list of mesh indices into V
  % Outputs:
  %   D  #P by dim list of normalized directions
  %   S  #P signed distances
  %   C  #P by dim list of closest points
  %

  dim = size(F,2);
  switch dim
  case 2
    % Facets are really edges
    E = F;
    % O(n*m) too slow...
    [T,sqrD] = project_to_lines(P,V(E(:,1),:),V(E(:,2),:),'Segments',true);
    % snap to line segment
    [~,J] = min(sqrD,[],2);
    T = T(sub2ind(size(T),1:size(P,1),J'))';
    C = V(E(J,1),:) + bsxfun(@times,T,(V(E(J,2),:)-V(E(J,1),:)));
    s = -2*(winding_number(V,E,P))+1;
    vec = C-P;
    % signed distance direction
    D = bsxfun(@times,s,normalizerow(vec));
    % signed distance
    S = s.*normrow(D);
  case 3
    [S,I,C,N] = signed_distance(P,V,F,'SignedDistanceType','pseudonormal');
    D = normalizerow(C-P);
    D = bsxfun(@times,sign(S),D);
    min_dist = 1e-5;
    too_close = abs(S) < min_dist;
    % Normals always point outside regardless of the eval point.
    D(too_close,:) = -N(too_close,:);
  end

  %% Find closest points
  %[sqrD,I,C] = point_mesh_squared_distance(P,V,F);
  %% Compute barycentric coordinates on closest faces
  %B = barycentric_coordinates(C,V(F(I,1),:),V(F(I,2),:),V(F(I,3),:));
  %% Direction to closest point
  %D = normalizerow(C-P);
  %% Determine which are too close to trust
  %min_sqr_dist = 1e-10;
  %too_close = sqrD < min_sqr_dist;
  %% Determine if closest point is on vertex, edge, or face AND too clost to
  %% trust direction.
  %epsilon = 1e-15;
  %on_face = (sum(B<=epsilon,2)==0) & too_close;
  %on_edge = (sum(B<=epsilon,2)==1) & too_close;
  %on_vertex = (sum(B<=epsilon,2)==2) & too_close;
  %% Determine which vertex or edge of that face
  %[~,which_vertex] = find(B(on_vertex,:)>epsilon);
  %[~,which_edge] = find(B(on_edge,:)<=epsilon);
  %% Compute normals for faces, vertices and edges
  %% Q: Are area normals OK?
  %% H: [Bærentzen 2002/2005] suggests angle weighting.
  %N_face = normalizerow(normals(V,F));
  %N_vertex = per_vertex_normals(V,F);
  %% map to corners
  %N_vertex = N_vertex(F,:);
  %[N_edge,E,EMAP] = per_edge_normals(V,F);
  %% map to directed edges
  %N_edge = N_edge(EMAP,:);
  %% This is an expensive way to find out if inside/outside a closed mesh
  %w = winding_number(V,F,P);
  %% Flip sign for interior points
  %S = (1-2*w);
  %D = bsxfun(@times,S,D);
  %% Use inverse normals for those that are too close
  %if any(on_vertex)
  %  D(on_vertex,:) = ...
  %    -N_vertex(sub2ind(size(F),I(on_vertex),which_vertex),:);
  %end
  %if any(on_edge)
  %  D(on_edge,:) = -N_edge(sub2ind(size(F),I(on_edge),which_edge),:);
  %end
  %if any(on_face)
  %  D(on_face,:) = -N_face(I(on_face),:);
  %end
end

author: alecjacobson

repo: gptoolbox


text:
function [J, grad] = costFunction(theta, X, y)
%COSTFUNCTION Compute cost and gradient for logistic regression
%   J = COSTFUNCTION(theta, X, y) computes the cost of using theta as the
%   parameter for logistic regression and the gradient of the cost
%   w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
% J = 0;
% grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Note: grad should have the same dimensions as theta
%

temp1 = -1 * (y .* log(sigmoid(X * theta)));
temp2 = (1 - y) .* log(1 - sigmoid(X * theta));

J = sum(temp1 - temp2) / m;

grad = (X' * (sigmoid(X * theta) - y)) * (1/m);

% =============================================================

end

repo: MachineLearning


text:
function Do_clustering(input, varargin)

% PROGRAM Do_clustering.
% Does clustering on all files in Files.txt
% Runs after Get_spikes, which prepares '*_spikes.mat' files.
%
% function Do_clustering(input, par_input)
% Saves spikes, spike times (in ms), coefficients used (inspk), used
% parameters, random spikes selected for clustering (ipermut) and
% results (cluster_class).
%
%input must be:
%               A .txt file with the names of the '*_spikes.mat' files to use.
%               A matlab cell with the names of the '*_spikes.mat' files to use.
%               A vector with channel numbers. In this case the function
%                   will proccess all the '*_spikes.mat' files located in
%                   the folder with those channel numbers (e.g., CSC1_spikes.mat or NSX4_spikes.mat)
%               'all', in this case the functions will process all the
%                   '*_spikes.mat' files in the folder.
%               A string with a filename ended in '_spikes.mat'
% optional argument 'par' and the next input must be a struct with some of
%       the SPC parameters (the detection parameters are taken from the
%       XXX_spikes file or the set_parameters file. All the parameters included
%       in the structure par will overwrite the parameters loaded from set_parameters()
% optional argument 'parallel' : true for use parallel computing.
% optional argument 'make_times': true for computing the sorting from the XXX_spikes files.
% optional argument 'make_plots': true for plotting the results based on the save XXX_times files.
% optional argument 'resolution': resolution string used by the plots (default: '-r150').
% optional argument 'save_spikes': false for disable saving all the spikes
% waveforms on the times* file.
% See also
% Get_spikes

% Example
% param.min_clus = 60;
% param.max_spk = 50000;
% par.maxtemp = 0.251;
% Do_clustering(1:16,'parallel',true,'par',param)
%
% If the times were created but something went wrong with the plots. Then, the
% code can be called again just to create the plots:
% Do_clustering(1:16,'make_times',false,'par',param)
%
%
% Similarly, for don't create the plots:
% Do_clustering(1:16,'make_plots',false,'par',param)
%

min_spikes4SPC = 16; % if are less that this number of spikes, clustering won't be made.

%default config
par_input = struct;
parallel = false;
make_times = true;
make_plots = true;
resolution = '-r150';
save_spikes = true;
%search for optional inputs
nvar = length(varargin);
for v = 1:nvar
    if strcmp(varargin{v},'par')
        if (nvar>=v+1) && isstruct(varargin{v+1})
            par_input = varargin{v+1};
        else
            error('Error in ''par'' optional input.')
        end
    elseif strcmp(varargin{v},'parallel')
        if (nvar>=v+1) && islogical(varargin{v+1})
            parallel = varargin{v+1};
        else
            error('Error in ''parallel'' optional input.')
        end
    elseif strcmp(varargin{v},'make_times')
        if (nvar>=v+1) && islogical(varargin{v+1})
            make_times = varargin{v+1};
        else
            error('Error in ''make_times'' optional input.')
        end
    elseif strcmp(varargin{v},'make_plots')
        if (nvar>=v+1) && islogical(varargin{v+1})
            make_plots = varargin{v+1};
        else
            error('Error in ''make_plots'' optional input.')
        end
    elseif strcmp(varargin{v},'resolution')
        if (nvar>=v+1) && ischar(varargin{v+1})
            resolution = varargin{v+1};
        else
            error('Error in ''resolution'' optional input.')
        end
    elseif strcmp(varargin{v},'save_spikes')
        if (nvar>=v+1) && islogical(varargin{v+1})
            save_spikes = varargin{v+1};
        else
            error('Error in ''save_spikes'' optional input.')
        end

    end
end

run_par_for = parallel;

% get a cell of filenames from the input
if isnumeric(input) || any(strcmp(input,'all'))  %cases for numeric or 'all' input
    filenames = {};
    dirnames = dir();
    dirnames = {dirnames.name};

    for i = 1:length(dirnames)
        fname = dirnames{i};

        if length(fname) < 12
            continue
        end
        if ~ strcmp(fname(end-10:end),'_spikes.mat')
            continue
        end
        if strcmp(input,'all')
            filenames = [filenames {fname}];
        else
            aux = regexp(fname, '\d+', 'match');
            if ~isempty(aux) && ismember(str2num(aux{1}),input)
                filenames = [filenames {fname}];
            end
        end
    end

elseif ischar(input) && length(input) > 4
    if  strcmp (input(end-3:end),'.txt')   %case for .txt input
        filenames =  textread(input,'%s');
    else
        filenames = {input};               %case for cell input
    end

elseif iscellstr(input)
    filenames = input;
else
    ME = MException('MyComponent:noValidInput', 'Invalid input arguments');
    throw(ME)
end

tic
par_file = set_parameters();
if make_times
% open parallel pool, if parallel input is true
    if parallel == true
        if exist('matlabpool','file')
            try
                matlabpool('open');
            catch
                parallel = false;
            end
        else
            poolobj = gcp('nocreate'); % If no pool, do not create new one.
            if isempty(poolobj)
                parpool
            else
                parallel = false;
            end
        end
    end

    initial_date = now;
    Nfiles = length(filenames);
    if run_par_for == true
        parfor fnum = 1:Nfiles
            filename = filenames{fnum};
            do_clustering_single(filename,min_spikes4SPC, par_file, par_input,fnum,save_spikes);
            disp(sprintf('%d of %d ''times'' files finished.',count_new_times(initial_date, filenames),Nfiles))
        end
    else
        for fnum = 1:length(filenames)
            filename = filenames{fnum};
            do_clustering_single(filename,min_spikes4SPC, par_file, par_input,fnum,save_spikes);
            disp(sprintf('%d of %d ''times'' files finished.',count_new_times(initial_date, filenames),Nfiles))
        end
    end
    if parallel == true
        if exist('matlabpool','file')
            matlabpool('close')
        else
            poolobj = gcp('nocreate');
            delete(poolobj);
        end
    end

	log_name = 'spc_log.txt';
	f = fopen(log_name, 'w');
	for fnum = 1:length(filenames)
        filename = filenames{fnum};
        log_name = [filename 'spc_log.txt'];
        if exist(log_name, 'file')
			fi = fopen(log_name,'r');
			result = fread(fi);
			fwrite(f,result);
			fclose(fi);
			delete(log_name);
		end
    end
	fclose(f);

	tocaux = toc;
    disp(['Computations Done (' num2str(tocaux,'%2.2f') 's).'])
end

if make_plots
    disp('Creating figures...')
    numfigs = length(filenames);
    curr_fig = figure('Visible','Off');
    curr_fig2 = figure('Visible','Off');
    set(curr_fig, 'PaperUnits', 'inches', 'PaperType', 'A4', 'PaperPositionMode', 'auto','units','normalized','outerposition',[0 0 1 1],'RendererMode','manual','Renderer','painters')
    set(curr_fig2, 'PaperUnits', 'inches', 'PaperType', 'A4', 'PaperPositionMode', 'auto','units','normalized','outerposition',[0 0 1 1],'RendererMode','manual','Renderer','painters')
    if isfield(get(curr_fig),'GraphicsSmoothing')
        set(curr_fig,'GraphicsSmoothing','off');
        set(curr_fig2,'GraphicsSmoothing','off');
    end

    for fnum = 1:numfigs
        clf(curr_fig2)
        clf(curr_fig)
        set(0, 'CurrentFigure', curr_fig)
        filename = filenames{fnum};
        par = struct;
        par.filename = filename;

        par.cont_segment = true;  %maybe true and save the sample in spikes

        data_handler = readInData(par);
        par = data_handler.update_par(par);

        par = update_parameters(par,par_file,'batch_plot');
        par = update_parameters(par,par_input,'batch_plot');

        if ~data_handler.with_wc_spikes       			%data should have spikes
            continue
        end
        filename = data_handler.nick_name;

        file_pos_names = {'','a','b','c','d','e','f'};
        for i=1:length(file_pos_names)
            new_file_name = ['fig2print_' filename file_pos_names{i} '.png'];
            if exist(new_file_name, 'file')==2
                delete(new_file_name);
            end
        end

        subplot(3,1,1)
        if par.cont_segment && data_handler.with_psegment
            box off; hold on
            %these lines are for plotting continuous data
            [xd_sub, sr_sub] = data_handler.get_signal_sample();
            lx = length(xd_sub);
            plot((1:lx)/sr_sub,xd_sub)
            noise_std_detect = median(abs(xd_sub))/0.6745;
            xlim([0 lx/sr_sub])
            thr = par.stdmin * noise_std_detect;
            plotmax = 15 * noise_std_detect; %thr*par.stdmax/par.stdmin;

            if strcmp(par.detection,'pos')
                line([0 length(xd_sub)/sr_sub],[thr thr],'color','r')
                ylim([-plotmax/2 plotmax])
            elseif strcmp(par.detection,'neg')
                line([0 length(xd_sub)/sr_sub],[-thr -thr],'color','r')
                ylim([-plotmax plotmax/2])
            else
                line([0 length(xd_sub)/sr_sub],[thr thr],'color','r')
                line([0 length(xd_sub)/sr_sub],[-thr -thr],'color','r')
                ylim([-plotmax plotmax])
            end
        end
        title([pwd '/' filename],'Interpreter','none','Fontsize',14)

        if ~data_handler.with_spc
            print2file = par_file.print2file;
            if isfield(par_input,'print2file')
                print2file = par_input.print2file;
            end
            if print2file
                print(curr_fig,'-dpng',['fig2print_' filename '.png'],resolution);
            else
                print(curr_fig)
            end
            clear print2file
            fprintf('%d figs Done. ',fnum);
            continue
        end

        % LOAD SPIKES
        [clu, tree, spikes, index, inspk, ipermut, classes, forced,temp] = data_handler.load_results();
        nspk = size(spikes,1);
        Mclasses = max(classes);
        auto_sort_info = [];
        if data_handler.with_gui_status
            [original_classes, current_temp, auto_sort_info] = data_handler.get_gui_status();
            temperature = par.mintemp+current_temp*par.tempstep;
            org_class = zeros(1,Mclasses);
            for i = 1:Mclasses
                org_class(i) = original_classes(find(classes==i,1,'first'));
            end
        else
            [temp] = find_temp_old(tree,par);
            temp = ones(max(classes),1)*temp;
            org_class = 1:max(classes);
            temperature = par.mintemp+temp*par.tempstep;
        end

        %PLOTS
        clus_pop = [];
        ylimit = [];
        subplot(3,5,11)

        color = [[0.0 0.0 1.0];[1.0 0.0 0.0];[0.0 0.5 0.0];[0.620690 0.0 0.0];[0.413793 0.0 0.758621];[0.965517 0.517241 0.034483];
        [0.448276 0.379310 0.241379];[1.0 0.103448 0.724138];[0.545 0.545 0.545];[0.586207 0.827586 0.310345];
        [0.965517 0.620690 0.862069];[0.620690 0.758621 1.]];
        maxc = size(color,1);

        hold on
        num_temp = floor((par.maxtemp -par.mintemp)/par.tempstep);     % total number of temperatures
        switch par.temp_plot
                case 'lin'
                   if ~isempty(auto_sort_info)
                        [xp, yp] = find(auto_sort_info.peaks);
                        ptemps = par.mintemp+(xp)*par.tempstep;
                        psize = tree(sub2ind(size(tree), xp,yp+4));
                        plot(ptemps,psize,'xk','MarkerSize',7,'LineWidth',0.9);
                        area(par.mintemp+par.tempstep.*[auto_sort_info.elbow,num_temp],max(ylim).*[1 1],'LineStyle','none','FaceColor',[0.9 0.9 0.9]);
                    end
                    plot([par.mintemp par.maxtemp-par.tempstep], ...
                    [par.min_clus par.min_clus],'k:',...
                    par.mintemp+(1:num_temp)*par.tempstep, ...
                    tree(1:num_temp,5:size(tree,2)),[temperature temperature],[1 tree(1,5)],'k:')
                    for i=1:Mclasses
                        tree_clus = tree(temp(i),4+org_class(i));
                        tree_temp = tree(temp(i)+1,2);
                        plot(tree_temp,tree_clus,'.','color',color(mod(i-1,maxc)+1,:),'MarkerSize',20);
                    end
                case 'log'
                    set(gca,'yscale','log');
                    if ~isempty(auto_sort_info)
                        [xp, yp] = find(auto_sort_info.peaks);
                        ptemps = par.mintemp+(xp)*par.tempstep;
                        psize = tree(sub2ind(size(tree), xp,yp+4));
                        semilogy(ptemps,psize,'xk','MarkerSize',7,'LineWidth',0.9);
                        area(par.mintemp+par.tempstep.*[auto_sort_info.elbow,num_temp],max(ylim).*[1 1],'LineStyle','none','FaceColor',[0.9 0.9 0.9],'basevalue',1);
                    end
                    semilogy([par.mintemp par.maxtemp-par.tempstep], ...
                    [par.min_clus par.min_clus],'k:',...
                    par.mintemp+(1:num_temp)*par.tempstep, ...
                    tree(1:num_temp,5:size(tree,2)),[temperature temperature],[1 tree(1,5)],'k:')

                    for i=1:Mclasses
                        tree_clus = tree(temp(i),4+org_class(i));
                        tree_temp = tree(temp(i)+1,2);
                        semilogy(tree_temp,tree_clus,'.','color',color(mod(i-1,maxc)+1,:),'MarkerSize',20);
                    end
        end
        xlim([par.mintemp, par.maxtemp])
        subplot(3,5,6)
        hold on

        class0 = find(classes==0);
            max_spikes=min(length(class0),par.max_spikes_plot);
            plot(spikes(class0(1:max_spikes),:)','k');
            xlim([1 size(spikes,2)]);
        subplot(3,5,10);
            hold on
            plot(spikes(class0(1:max_spikes),:)','k');
            plot(mean(spikes(class0,:),1),'c','linewidth',2)
            xlim([1 size(spikes,2)]);
            title(['Cluster 0: # ' num2str(length(class0))],'Fontweight','bold')
        subplot(3,5,15)
            xa=diff(index(class0));
            [n,c]=hist(xa,0:1:100);
            bar(c(1:end-1),n(1:end-1))
            xlim([0 100])
            xlabel('ISI (ms)');
            title([num2str(nnz(xa<3)) ' in < 3ms']);

        for i = 1:max(classes)
            set(0, 'CurrentFigure', curr_fig)
            class = find(classes==i);
            subplot(3,5,6);
            max_spikes=min(length(class),par.max_spikes_plot);
            plot(spikes(class(1:max_spikes),:)','color',color(mod(i-1,maxc)+1,:));
            xlim([1 size(spikes,2)]);

            if i<=3
                subplot(3,5,6+i);
                hold on
                plot(spikes(class(1:max_spikes),:)','color',color(mod(i-1,maxc)+1,:));
                plot(mean(spikes(class,:),1),'k','linewidth',2)
                xlim([1 size(spikes,2)]);
                title(['Cluster ' num2str(i) ': # ' num2str(length(class)) ' (' num2str(nnz(classes(:)==i & ~forced(:))) ')'],'Fontweight','bold')
                ylimit = [ylimit;ylim];
                subplot(3,5,11+i)
                xa=diff(index(class));
                [n,c]=hist(xa,0:1:100);
                bar(c(1:end-1),n(1:end-1))
                xlim([0 100])
                xlabel('ISI (ms)');
                title([num2str(nnz(xa<3)) ' in < 3ms']);

            elseif i<=8
            	  set(0, 'CurrentFigure', curr_fig2)
                subplot(3,5,2+i);
                hold on
                plot(spikes(class(1:max_spikes),:)','color',color(mod(i-1,maxc)+1,:));
                plot(mean(spikes(class,:),1),'k','linewidth',2)
                xlim([1 size(spikes,2)]);
                title(['Cluster ' num2str(i) ': # ' num2str(length(class)) ' (' num2str(nnz(classes(:)==i & ~forced(:))) ')'],'Fontweight','bold')
                ylimit = [ylimit;ylim];
                subplot(3,5,7+i)
                xa=diff(index(class));
                [n,c]=hist(xa,0:1:100);
                bar(c(1:end-1),n(1:end-1))
                xlim([0 100])
                xlabel('ISI (ms)');
                title([num2str(nnz(xa<3)) ' in < 3ms']);
            end
        end
        numclus = max(classes);
        % Rescale spike's axis
        if ~isempty(ylimit)
            ymin = min(ylimit(:,1));
            ymax = max(ylimit(:,2));
            set(0, 'CurrentFigure', curr_fig)
            for i = 1:min(3,numclus)
               subplot(3,5,6+i); ylim([ymin ymax]);
            end
            set(0, 'CurrentFigure', curr_fig2)
            for i = 4:min(8,numclus)
               subplot(3,5,2+i); ylim([ymin ymax]);
            end

        end

        features_name = par.features;

        outfileclus='cluster_results.txt';
        fout=fopen(outfileclus,'at+');
        if isfield(par,'stdmin')
            stdmin = par.stdmin;
        else
            stdmin = NaN;
        end
        fprintf(fout,'%s\t %s\t %g\t %d\t %g\t', char(filename), features_name, temperature, numclus, stdmin);
        for ii=0:numclus
            fprintf(fout,'%d\t',nnz(classes==ii));
        end
        fclose(fout);

        if par.print2file
            print(curr_fig,'-dpng',['fig2print_' filename '.png'],resolution);
            if numclus>3
                print(curr_fig2,'-dpng',['fig2print_' filename 'a.png'],resolution);
            end
        else
            print(curr_fig)
            if numclus>3
                print(curr_fig2);
            end
        end
        fprintf('%d ',fnum);
    end
    close(curr_fig)
    close(curr_fig2)
    disp(' ')
end

toc

end

function do_clustering_single(filename,min_spikes4SPC, par_file, par_input,fnum,save_spikes)

    par = struct;
    par = update_parameters(par,par_file,'clus');
    par = update_parameters(par,par_file,'batch_plot');
    par = update_parameters(par,par_input,'clus');
    par = update_parameters(par,par_input,'batch_plot');

    par.filename = filename;
    par.reset_results = true;

    data_handler = readInData(par);
    par = data_handler.par;
    check_WC_params(par)
%     if isfield(par,'channels') && ~isnan(par.channels)
%         par.max_inputs = par.max_inputs * par.channels;
%     end

    par.fname_in = ['tmp_data_wc' num2str(fnum)];                       % temporary filename used as input for SPC
    par.fname = ['data_' data_handler.nick_name];
    par.nick_name = data_handler.nick_name;
    par.fnamespc = ['data_wc' num2str(fnum)];

    if data_handler.with_spikes            			%data have some time of _spikes files
        [spikes, index] = data_handler.load_spikes();
    else
        warning('MyComponent:noValidInput', 'File: %s doesn''t include spikes', filename);
        return
    end

    % LOAD SPIKES
    nspk = size(spikes,1);
    naux = min(par.max_spk,size(spikes,1));

    if nspk < min_spikes4SPC
        warning('MyComponent:noValidInput', 'Not enough spikes in the file');
        return
    end

    % CALCULATES INPUTS TO THE CLUSTERING ALGORITHM.
    inspk = wave_features(spikes,par);     %takes wavelet coefficients.
    par.inputs = size(inspk,2);                       % number of inputs to the clustering

	if par.permut == 'n'
        % GOES FOR TEMPLATE MATCHING IF TOO MANY SPIKES.
        if size(spikes,1)> par.max_spk;
            % take first 'par.max_spk' spikes as an input for SPC
            inspk_aux = inspk(1:naux,:);
        else
            inspk_aux = inspk;
        end
	else
        % GOES FOR TEMPLATE MATCHING IF TOO MANY SPIKES.
        if size(spikes,1)> par.max_spk;
            % random selection of spikes for SPC
            ipermut = randperm(length(inspk));
            ipermut(naux+1:end) = [];
            inspk_aux = inspk(ipermut,:);
        else
            ipermut = randperm(size(inspk,1));
            inspk_aux = inspk(ipermut,:);
        end
	end
    %INTERACTION WITH SPC
    save(par.fname_in,'inspk_aux','-ascii');
    try
        [clu, tree] = run_cluster(par,true);
		if exist([par.fnamespc '.dg_01.lab'],'file')
			movefile([par.fnamespc '.dg_01.lab'], [par.fname '.dg_01.lab'], 'f');
			movefile([par.fnamespc '.dg_01'], [par.fname '.dg_01'], 'f');
		end
    catch
        warning('MyComponent:ERROR_SPC', 'Error in SPC');
        return
    end

    [clust_num temp auto_sort] = find_temp(tree,clu,par);

    if par.permut == 'y'
        clu_aux = zeros(size(clu,1),2 + size(spikes,1)) -1;  %when update classes from clu, not selected elements go to cluster 0
        clu_aux(:,ipermut+2) = clu(:,(1:length(ipermut))+2);
        clu_aux(:,1:2) = clu(:,1:2);
        clu = clu_aux;
        clear clu_aux
    end

    classes = zeros(1,size(clu,2)-2);
    for c =1: length(clust_num)
        aux = clu(temp(c),3:end) +1 == clust_num(c);
        classes(aux) = c;
    end

    if par.permut == 'n'
        classes = [classes zeros(1,max(size(spikes,1)-par.max_spk,0))];
    end

    Temp = [];
    % Classes should be consecutive numbers
    classes_names = nonzeros(sort(unique(classes)));
    for i= 1:length(classes_names)
       c = classes_names(i);
       if c~= i
           classes(classes == c) = i;
       end
       Temp(i) = temp(i);
    end

    % IF TEMPLATE MATCHING WAS DONE, THEN FORCE
    if (size(spikes,1)> par.max_spk || ...
            (par.force_auto))
        f_in  = spikes(classes~=0,:);
        f_out = spikes(classes==0,:);
        class_in = classes(classes~=0);
        class_out = force_membership_wc(f_in, class_in, f_out, par);
        forced = classes==0;
        classes(classes==0) = class_out;
        forced(classes==0) = 0;
    else
        forced = zeros(1, size(spikes,1));
    end

    gui_status = struct();
    gui_status.current_temp =  max(temp);
    gui_status.auto_sort_info = auto_sort;
    gui_status.original_classes = zeros(size(classes));

    for i=1:max(classes)
        gui_status.original_classes(classes==i) = clust_num(i);
    end

    current_par = par;
    par = struct;
    par = update_parameters(par, current_par, 'relevant');
    par = update_parameters(par,current_par,'batch_plot');

    par.sorting_date = datestr(now);
    cluster_class = zeros(nspk,2);
    cluster_class(:,2)= index';
    cluster_class(:,1)= classes';

    vars = {'cluster_class', 'par','inspk','forced','Temp','gui_status'};
    if exist('ipermut','var')
        vars{end+1} = 'ipermut';
    end
    if save_spikes
        vars{end+1} = 'spikes';
    else
        spikes_file = filename;
        vars{end+1} = 'spikes_file';
    end

    try
      save(['times_' data_handler.nick_name],vars{:});
    catch
      save(['times_' data_handler.nick_name],vars{:},'-v7.3');
    end

end

function counter = count_new_times(initial_date, filenames)
counter = 0;
for i = 1:length(filenames)
    fname = filenames{i};
    FileInfo = dir(['times_' fname(1:end-11) '.mat']);
    if length(FileInfo)==1 && (FileInfo.datenum > initial_date)
        counter = counter + 1;
    end
end
end

author: csn-le


text:
% SB2_USEROPTIONS  User option specification for SPARSEBAYES
%
% OPTIONS = SB2_USEROPTIONS(parameter1, value1, parameter2, value2,...)
%
% OUTPUT ARGUMENTS:
%
%	OPTIONS		An options structure to pass to SPARSEBAYES
%
% INPUT ARGUMENTS:
%
%	Optional number of parameter-value pairs to specify the following:
%
%	ITERATIONS	Number of interations to run for.
%
%	TIME		Time limit to run for, expressed as a space-separated
%				string. e.g. '1.5 hours', '30 minutes', '1 second'.
%
%	DIAGNOSTICLEVEL	Integer [0,4] or string to determine the verbosity of
%					diagnostic output.
%					0 or 'ZERO' or 'NONE'	No output
%					1 or 'LOW'				Low level of output
%					2 or 'MEDIUM'			etc...
%					3 or 'HIGH'
%					4 or 'ULTRA'
%
%	DIAGNOSTICFILE	Filename to write diagnostics to file instead of
%					the default stdout.
%
%	MONITOR		Integer number: diagnostic information is output
%				every MONITOR iterations.
%
%	FIXEDNOISE	True/false whether the Gaussian noise is to be fixed
%				(default: false.
%
%	FREEBASIS	Indices of basis vectors considered "free" and not
%				constrained by the Bayesian prior (e.g. the "bias").
%
%	CALLBACK	External function to call each iteration of the algorithm
%				(string). Intended to facilitate graphical demos etc.
%
%	CALLBACKDATA	Arbitrary additional data to pass to the CALLBACK
%					function.
%
% EXAMPLE:
%
%	OPTIONS = SB2_UserOptions('diagnosticLevel','medium',...
%				  'monitor',25,...
%				  'diagnosticFile', 'logfile.txt');
%
% NOTES:
%
% Each option (field of OPTIONS) is given a default value in
% SB2_USEROPTIONS. Any supplied property-value pairs over-ride those
% defaults.
%

%
% Copyright 2009, Vector Anomaly Ltd
%
% This file is part of the SPARSEBAYES library for Matlab (V2.0).
%
% SPARSEBAYES 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.
%
% SPARSEBAYES 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 SPARSEBAYES in the accompanying file "licence.txt"; if not, write to
% the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
% MA 02110-1301 USA
%
% Contact the author: m a i l [at] m i k e t i p p i n g . c o m
%
function OPTIONS = SB2_UserOptions(varargin)

% Ensure arguments are supplied in pairs
%
if rem(nargin,2)
  error('Arguments to SB2_UserOptions should be (property, value) pairs')
end
% Any options specified?
numSettings	= nargin/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Set defaults
%
% Assume we will infer the noise in the Gaussian case
%
OPTIONS.fixedNoise		= false;
%
% Option to allow subset of the basis (e.g. bias) to be unregularised
%
OPTIONS.freeBasis		= [];
%
% Option to set max iterations to run for
%
OPTIONS.iterations		= 10000;
%
% Option to set max time to run for
%
OPTIONS.time			= 10000; % seconds
%
% Set options for monitoring and recording the algorithm's progress
%
OPTIONS.monitor			= 0;
OPTIONS.diagnosticLevel	= 0;
OPTIONS.diagnosticFID	= 1; % stdout
OPTIONS.diagnosticFile_	= [];
%
% Option to call a function during each iteration (to create demos etc)
%
OPTIONS.callback		= false;
OPTIONS.callbackFunc	= [];
OPTIONS.callbackData	= {};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Parse string/variable pairs

for n=1:numSettings
  property_	= varargin{(n-1)*2+1};
  value		= varargin{(n-1)*2+2};
  switch upper(property_)
    %
   case 'FIXEDNOISE',
    OPTIONS.fixedNoise	= value;
    %
   case 'FREEBASIS',
    OPTIONS.freeBasis	= value;
    %
   case 'ITERATIONS',
    OPTIONS.iterations	= value;
	%
   case 'TIME',
    OPTIONS.time		= timeInSeconds(value);
    %
   case 'MONITOR',
    OPTIONS.monitor		= value;
    %
   case 'DIAGNOSTICLEVEL',
    OPTIONS.diagnosticLevel	= value;
    MAX_LEVEL	= 4;
    if ischar(value)
      switch upper(value)
       case {'ZERO','NONE'},
	OPTIONS.diagnosticLevel	= 0;
       case 'LOW',
	OPTIONS.diagnosticLevel	= 1;
       case 'MEDIUM',
	OPTIONS.diagnosticLevel	= 2;
       case 'HIGH',
	OPTIONS.diagnosticLevel	= 3;
       case 'ULTRA',
	OPTIONS.diagnosticLevel	= 4;
       otherwise,
	error('Unrecognised textual diagnostic level: ''%s''\n', value)
      end
    elseif isnumeric(value)
      if value>=0 & value<=MAX_LEVEL
	OPTIONS.diagnosticLevel		= value;
      else
	error(['Supplied level should be integer in [0,%d], '...
	       'or one of ZERO/LOW/MEDIUM/HIGH/ULTRA'], MAX_LEVEL);

      end
    end
    %
   case 'DIAGNOSTICFILE',
    OPTIONS.diagnosticFile_	= value;
    OPTIONS.diagnosticFID	= -1; % It will be opened later
    %
   case 'CALLBACK',
    OPTIONS.callback		= true;
    OPTIONS.callbackFunc	= value;
    if exist(OPTIONS.callbackFunc)~=2
      error('Callback function ''%s'' does not appear to exist\n', value)
    end
    %
   case 'CALLBACKDATA',
    OPTIONS.callbackData	= value;
   otherwise,
    error('Unrecognised user option: ''%s''\n', property_)
    %
  end
end

%%
%% Support function: parse time specification
%%
function s = timeInSeconds(value_)
%
[v_ r_]				= strtok(value_,' ');
v					= str2num(v_);
r_(isspace(r_))		= [];
switch upper(r_)
 case {'SECONDS', 'SECOND'}
  s	= v;
 case {'MINUTES', 'MINUTE'}
  s = 60*v;
 case {'HOURS', 'HOUR'}
  s = 3600*v;
 otherwise,
  error('Badly formed time string: ''%s''', value_)
end
% Returns time in seconds
%

repo: REST


text:
function demo_blockproc_denoising(source,varargin) %RUNASSCRIPT
%DEMO_BLOCKPROC_DENOISING Variable coefficients thresholding
%   Usage: demo_blockproc_denoising('gspi.wav')
%
%   For additional help call |demo_blockproc_denoising| without arguments.
%
%   The present demo allows you to set the coefficient threshold during the
%   playback using the control panel.
%

if demo_blockproc_header(mfilename,nargin)
   return;
end

% Control pannel (Java object)
% Each entry determines one parameter to be changed during the main loop
% execution.
p = blockpanel({
               {'GdB','Gain',-20,20,0,21},...
               {'Thr','Treshold',0,0.1,0,1000}
               });

% Number of frequency channels
M = 1000;

% Setup blocktream
try
    fs=block(source,varargin{:},'loadind',p);
catch
    % Close the windows if initialization fails
    blockdone(p);
    err = lasterror;
    error(err.message);
end

% Buffer length (30 ms)
bufLen = floor(30e-3*fs);

% Window length in ms
winLenms = 20; %floor(fs*winLenms/1e3)
[F,Fdual] = framepair('dgtreal',{'hann',floor(fs*winLenms/1e3)},'dual',40,M);
% Or using fwt
%[F,Fdual] = framepair('fwt','ana:symorth3','dual',7);
[Fa,Fs] = blockframepairaccel(F,Fdual, bufLen,'segola');

flag = 1;
%Loop until end of the stream (flag) and until panel is opened
while flag && p.flag

  % Obtain parameters from the control panel
  gain = 10^(p.getParam('GdB')/20); % dB -> val
  thres = p.getParam('Thr');
  %bufLen = floor(p.getParam('bufLen'));

  % Read block of length bufLen
  [f,flag] = blockread(bufLen);
  % Apply analysis frame
  c = blockana(Fa, f*gain);
  % Plot
  % blockplot(fobj,F,c);
  % Apply thresholding
  c = thresh(c,thres,'soft');
  % Apply synthesis frame
  fhat = real(blocksyn(Fs, c, size(f,1)));
  % Play the block
  %fhat = f;
  blockplay(fhat);
end
blockdone(p);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat