text:
function spikes = bz_GetSpikes(varargin)
% bz_getSpikes - Get spike timestamps.
%       if loading from clu/res/fet/spk files - must be formatted as:
%       baseName.clu.shankNum
%
% USAGE
%
%    spikes = bz_getSpikes(varargin)
%
% INPUTS
%
%    spikeGroups     -vector subset of shank IDs to load (Default: all)
%    region          -string region ID to load neurons from specific region
%                     (requires sessionInfo file or units->structures in xml)
%    UID             -vector subset of UID's to load
%    basepath        -path to recording (where .dat/.clu/etc files are)
%    getWaveforms    -logical (default=true) to load mean of raw waveform data
%    forceReload     -logical (default=false) to force loading from
%                     res/clu/spk files
%    onlyLoad        -[shankID cluID] pairs to EXCLUSIVELY LOAD from
%                       clu/res/fet to spikes.cellinfo.mat file
%    saveMat         -logical (default=false) to save in buzcode format
%    noPrompts       -logical (default=false) to supress any user prompts
%
% OUTPUTS
%
%    spikes - cellinfo struct with the following fields
%          .sessionName    -name of recording file
%          .UID            -unique identifier for each neuron in a recording
%          .times          -cell array of timestamps (seconds) for each neuron
%          .spindices      -sorted vector of [spiketime UID], useful for
%                           input to some functions and plotting rasters
%          .region         -region ID for each neuron (especially important large scale, high density probes)
%          .shankID        -shank ID that each neuron was recorded on
%          .maxWaveformCh  -channel # with largest amplitude spike for each neuron
%          .rawWaveform    -average waveform on maxWaveformCh (from raw .dat)
%          .cluID          -cluster ID, NOT UNIQUE ACROSS SHANKS
%          .numcells       -number of cells/UIDs
%
% NOTES
%
% This function can be used in several ways to load spiking data.
% Specifically, it loads spiketimes for individual neurons and other
% sessionInfodata that describes each neuron.  Spiketimes can be loaded using the
% UID(1-N), the shank the neuron was on, or the region it was recorded in.
% The default behavior is to load all spikes in a recording. The .shankID
% and .cluID fields can be used to reconstruct the 'units' variable often
% used in FMAToolbox.
% units = [spikes.shankID spikes.cluID];
%
%
% first usage recommendation:
%
%   spikes = bz_getSpikes('saveMat',true); Loads and saves all spiking data
%                                          into buzcode format .cellinfo. struct
% other examples:
%
%   spikes = bz_getSpikes('spikeGroups',1:5); first five shanks
%
%   spikes = bz_getSpikes('region','CA1'); cells tagged as recorded in CA1
%
%   spikes = bz_getSpikes('UID',[1:20]); first twenty neurons
%
%
% written by David Tingley, 2017
%% Deal With Inputs
spikeGroupsValidation = @(x) assert(isnumeric(x) || strcmp(x,'all'),...
    'spikeGroups must be numeric or "all"');

p = inputParser;
addParameter(p,'spikeGroups','all',spikeGroupsValidation);
addParameter(p,'region','',@isstr); % won't work without sessionInfodata
addParameter(p,'UID',[],@isvector);
addParameter(p,'basepath',pwd,@isstr);
addParameter(p,'getWaveforms',true)
addParameter(p,'forceReload',false,@islogical);
addParameter(p,'saveMat',false,@islogical);
addParameter(p,'noPrompts',false,@islogical);
addParameter(p,'onlyLoad',[]);

parse(p,varargin{:})

spikeGroups = p.Results.spikeGroups;
region = p.Results.region;
UID = p.Results.UID;
basepath = p.Results.basepath;
getWaveforms = p.Results.getWaveforms;
forceReload = p.Results.forceReload;
saveMat = p.Results.saveMat;
noPrompts = p.Results.noPrompts;
onlyLoad = p.Results.onlyLoad;

[sessionInfo] = bz_getSessionInfo(basepath, 'noPrompts', noPrompts);
baseName = bz_BasenameFromBasepath(basepath);

spikes.samplingRate = sessionInfo.rates.wideband;
nChannels = sessionInfo.nChannels;

cellinfofile = [basepath filesep sessionInfo.FileName '.spikes.cellinfo.mat'];
%% if the cellinfo file exist and we don't want to re-load files
if exist(cellinfofile,'file') && forceReload == false
    disp('loading spikes from cellinfo file..')
    load(cellinfofile)
    %Check that the spikes structure fits cellinfo requirements
    [iscellinfo] = bz_isCellInfo(spikes);
    switch iscellinfo
        case false
            warning(['The spikes structure in baseName.spikes.cellinfo.mat ',...
                'does not fit buzcode standards. Sad.'])
    end

    %If regions have been added since creation... add them
    if ~isfield(spikes,'region') & isfield(sessionInfo,'region')
        if ~isfield(spikes,'numcells')
            spikes.numcells = length(spikes.UID);
        end
        for cc = 1:spikes.numcells
            spikes.region{cc} = sessionInfo.region{spikes.maxWaveformCh(cc)==sessionInfo.channels};
        end

        if saveMat
            save(cellinfofile,'spikes')
        end
    end

else % do the below then filter by inputs... (Load from clu/res/fet)

    if ~noPrompts & saveMat == 0 %Inform the user that they should save a file for later
        savebutton = questdlg(['Would you like to save your spikes in ',...
            sessionInfo.FileName,'.spikes.cellinfo.mat?  ',...
            'This will save significant load time later.']);
        if strcmp(savebutton,'Yes'); saveMat = true; end
    end

disp('loading spikes from clu/res/spk files..')
% find res/clu/fet/spk files here
cluFiles = dir([basepath filesep '*.clu*']);
resFiles = dir([basepath filesep '*.res.*']);
if any(getWaveforms)
    spkFiles = dir([basepath filesep '*.spk*']);
end

% remove *temp*, *autosave*, and *.clu.str files/directories
tempFiles = zeros(length(cluFiles),1);
for i = 1:length(cluFiles)
    dummy = strsplit(cluFiles(i).name, '.'); % Check whether the component after the last dot is a number or not. If not, exclude the file/dir.
    if ~isempty(findstr('temp',cluFiles(i).name)) | ~isempty(findstr('autosave',cluFiles(i).name)) | ...
            isempty(str2num(dummy{length(dummy)})) | find(contains(dummy, 'clu')) ~= length(dummy)-1  | ...
            ~strcmp(dummy{1},baseName)
        tempFiles(i) = 1;
    end
end
cluFiles(tempFiles==1)=[];
tempFiles = zeros(length(resFiles),1);
for i = 1:length(resFiles)
    dummy = strsplit(resFiles(i).name, '.');
    if ~isempty(findstr('temp',resFiles(i).name)) | ~isempty(findstr('autosave',resFiles(i).name)) | ...
            ~strcmp(dummy{1},baseName)
        tempFiles(i) = 1;
    end
end
if any(getWaveforms)
    resFiles(tempFiles==1)=[];
    tempFiles = zeros(length(spkFiles),1);
    for i = 1:length(spkFiles)
        dummy = strsplit(spkFiles(i).name, '.');
        if ~isempty(findstr('temp',spkFiles(i).name)) | ~isempty(findstr('autosave',spkFiles(i).name)) | ...
            ~strcmp(dummy{1},baseName)
            tempFiles(i) = 1;
        end
    end
    spkFiles(tempFiles==1)=[];
end

if isempty(cluFiles)
    disp('no clu files found...')
    spikes = [];
    return
end

% ensures we load in sequential order (forces compatibility with FMAT
% ordering)
for i = 1:length(cluFiles)
    temp = strsplit(cluFiles(i).name,'.');
    shanks(i) = str2num(temp{length(temp)});
end
[shanks ind] = sort(shanks);
cluFiles = cluFiles(ind); %Bug here if there are any files x.clu.x that are not your desired clus
resFiles = resFiles(ind);
if any(getWaveforms)
    spkFiles = spkFiles(ind);
end

% check if there are matching #'s of files
if length(cluFiles) ~= length(resFiles) & length(cluFiles) ~= length(spkFiles)
    error('found an incorrect number of res/clu/spk files...')
end

% use the .clu files to get spike ID's and generate UID and spikeGroup
% use the .res files to get spike times
count = 1;

if isempty(sessionInfo.spikeGroups.groups)
    sessionInfo.spikeGroups = sessionInfo.AnatGrps;
end
for i=1:length(cluFiles)
    disp(['working on ' cluFiles(i).name])

    temp = strsplit(cluFiles(i).name,'.');
    shankID = str2num(temp{length(temp)}); %shankID is the spikegroup number
    clu = load(fullfile(basepath,cluFiles(i).name));
    clu = clu(2:end); % toss the first sample to match res/spk files
    res = load(fullfile(basepath,resFiles(i).name));
    spkGrpChans = sessionInfo.spikeGroups.groups{shankID}; % we'll eventually want to replace these two lines

    if any(getWaveforms) && sum(clu)>0 %bug fix if no clusters
        nSamples = sessionInfo.spikeGroups.nSamples(shankID);
        % load waveforms
        chansPerSpikeGrp = length(sessionInfo.spikeGroups.groups{shankID});
        fid = fopen(fullfile(basepath,spkFiles(i).name),'r');
        wav = fread(fid,[1 inf],'int16=>int16');
        try %bug in some spk files... wrong number of samples?
            wav = reshape(wav,chansPerSpikeGrp,nSamples,[]);
        catch
            if strcmp(getWaveforms,'force')
                wav = nan(chansPerSpikeGrp,nSamples,length(clu));
                display([spkFiles(i).name,' error.'])
            else
            error(['something is wrong with ',spkFiles(i).name,...
                ' Use ''getWaveforms'', false to skip waveforms or ',...
                '''getWaveforms'', ''force'' to write nans on bad shanks.'])
            end
        end
        wav = permute(wav,[3 1 2]);
    end

    cells  = unique(clu);
    % remove MUA and NOISE clusters...
    cells(cells==0) = [];
    cells(cells==1) = [];  % consider adding MUA as another input argument...?

    for c = 1:length(cells)
       spikes.UID(count) = count; % this only works if all shanks are loaded... how do we optimize this?
       ind = find(clu == cells(c));
       spikes.times{count} = res(ind) ./ spikes.samplingRate;
       spikes.shankID(count) = shankID;
       spikes.cluID(count) = cells(c);

       %Waveforms
       if any(getWaveforms)
           wvforms = squeeze(mean(wav(ind,:,:)))-mean(mean(mean(wav(ind,:,:)))); % mean subtract to account for slower (theta) trends
           if prod(size(wvforms))==length(wvforms)%in single-channel groups wvforms will squeeze too much and will have amplitude on D1 rather than D2
               wvforms = wvforms';%fix here
           end
           for t = 1:size(wvforms,1)
              [a(t) b(t)] = max(abs(wvforms(t,:)));
           end
           [aa bb] = max(a,[],2);
           spikes.rawWaveform{count} = wvforms(bb,:);
           spikes.maxWaveformCh(count) = spkGrpChans(bb);
           %Regions (needs waveform peak)
           if isfield(sessionInfo,'region') %if there is regions field in your metadata
                spikes.region{count} = sessionInfo.region{find(spkGrpChans(bb)==sessionInfo.channels)};
           elseif isfield(sessionInfo,'Units') %if no regions, but unit region from xml via Loadparamteres
                %Find the xml Unit that matches group/cluster
                unitnum = cellfun(@(X,Y) X==spikes.shankID(count) && Y==spikes.cluID(count),...
                    {sessionInfo.Units(:).spikegroup},{sessionInfo.Units(:).cluster});
                if sum(unitnum) == 0
                    display(['xml Missing Unit - spikegroup: ',...
                        num2str(spikes.shankID(count)),' cluster: ',...
                        num2str(spikes.cluID(count))])
                    spikes.region{count} = 'missingxml';
                else %possible future bug: two xml units with same group/clu...
                    spikes.region{count} = sessionInfo.Units(unitnum).structure;
                end
           end
           clear a aa b bb
       end

       count = count + 1;
    end
end

if ~isempty(onlyLoad)
    toRemove = true(size(spikes.UID));
    for cc = 1:size(onlyLoad,1)
        whichUID = ismember(spikes.shankID,onlyLoad(cc,1)) & ismember(spikes.cluID,onlyLoad(cc,2));
        toRemove(whichUID) = false;
        if ~any(whichUID)
            display(['No unit with shankID:',num2str(onlyLoad(cc,1)),...
                ' cluID:',num2str(onlyLoad(cc,2))])
        end
    end
    spikes = removeCells(toRemove,spikes,getWaveforms);
end

spikes.sessionName = sessionInfo.FileName;
end

%% save to buzcode format (before exclusions)
if saveMat
    save(cellinfofile,'spikes')
end

%% EXCLUSIONS %%

%filter by spikeGroups input
if ~strcmp(spikeGroups,'all')
    [toRemove] = ~ismember(spikes.shankID,spikeGroups);
    spikes = removeCells(toRemove,spikes,getWaveforms);
end

%filter by region input
if ~isempty(region)
    if ~isfield(spikes,'region') %if no region information in metadata
        error(['You selected to load cells from region "',region,...
            '", but there is no region information in your sessionInfo'])
    end

  toRemove = ~ismember(spikes.region,region);
    if sum(toRemove)==length(spikes.UID) %if no cells from selected region
        warning(['You selected to load cells from region "',region,...
            '", but none of your cells are from that region'])
    end

    spikes = removeCells(toRemove,spikes,getWaveforms);
end

%filter by UID input
if ~isempty(UID)
	[toRemove] = ~ismember(spikes.UID,UID);
    spikes = removeCells(toRemove,spikes,getWaveforms);
end

%% Generate spindices matrics
spikes.numcells = length(spikes.UID);
for cc = 1:spikes.numcells
    groups{cc}=spikes.UID(cc).*ones(size(spikes.times{cc}));
end
if spikes.numcells>0
    alltimes = cat(1,spikes.times{:}); groups = cat(1,groups{:}); %from cell to array
    [alltimes,sortidx] = sort(alltimes); groups = groups(sortidx); %sort both
    spikes.spindices = [alltimes groups];
end

%% Check if any cells made it through selection
if isempty(spikes.times) | spikes.numcells == 0
    spikes = [];
end

end

%%
function spikes = removeCells(toRemove,spikes,getWaveforms)
%Function to remove cells from the structure. toRemove is the INDEX of
%the UID in spikes.UID

    spikes.UID(toRemove) = [];
    spikes.times(toRemove) = [];
    spikes.region(toRemove) = [];
    spikes.cluID(toRemove) = [];
    spikes.shankID(toRemove) = [];

    if any(getWaveforms)
        spikes.rawWaveform(toRemove) = [];
        spikes.maxWaveformCh(toRemove) = [];
    end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [dataout] = ft_denoise_dssp(cfg, datain)

% FT_DENOISE_DSSP implements a dual signal subspace projection algorithm
% to suppress interference outside a predefined source region of
% interest. It is based on: Sekihara et al. J. Neural Eng. 2016 13(3), and
% Sekihara et al. J. Neural Eng. 2018 15(3).
%
% Use as
%   dataout = ft_denoise_dssp(cfg, datain)
% where cfg is a configuration structure that contains
%   cfg.channel          = Nx1 cell-array with selection of channels (default = 'all'), see FT_CHANNELSELECTION for details
%   cfg.trials           = 'all' or a selection given as a 1xN vector (default = 'all')
%   cfg.pertrial         = 'no', or 'yes', compute the temporal projection per trial (default = 'no')
%   cfg.sourcemodel      = structure, source model with precomputed leadfields, see FT_PREPARE_LEADFIELD
%   cfg.demean           = 'yes', or 'no', demean the data per epoch (default = 'yes')
%   cfg.dssp             = structure with parameters that determine the behavior of the algorithm
%   cfg.dssp.n_space     = 'all', or scalar. Number of dimensions for the
%                          initial spatial projection.
%   cfg.dssp.n_in        = 'all', or scalar. Number of dimensions of the
%                          subspace describing the field inside the ROI.
%   cfg.dssp.n_out       = 'all', or scalar. Number of dimensions of the
%                          subspace describing the field outside the ROI.
%   cfg.dssp.n_intersect = scalar (default = 0.9). Number of dimensions (if
%                          value is an integer>=1), or threshold for the
%                          included eigenvalues (if value<1), determining
%                          the dimensionality of the intersection.
%
% See also FT_DENOISE_PCA, FT_DENOISE_SYNTHETIC, FT_DENOISE_TSR

% Copyright (C) 2018-2023, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar    datain
ft_preamble provenance datain

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  % do not continue function execution in case the outputfile is present and the user indicated to keep it
  return
end

% check the input data
datain = ft_checkdata(datain, 'datatype', {'raw'}); % FIXME how about timelock and freq?

% ensure the external cellfunction toolbox is on the path
ft_hastoolbox('cellfunction', 1);

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels', 'trial'}); % prevent accidental typos, see issue 1729

% set the defaults
cfg.trials            = ft_getopt(cfg, 'trials',  'all', 1);
cfg.channel           = ft_getopt(cfg, 'channel', 'all');
cfg.pertrial          = ft_getopt(cfg, 'pertrial', 'no');
cfg.sourcemodel       = ft_getopt(cfg, 'sourcemodel');
cfg.demean            = ft_getopt(cfg, 'demean', 'yes');
cfg.dssp              = ft_getopt(cfg, 'dssp');         % sub-structure to hold the parameters
cfg.dssp.n_space      = ft_getopt(cfg.dssp, 'n_space', 'interactive'); % number of spatial components to retain from the Gram matrix
cfg.dssp.n_in         = ft_getopt(cfg.dssp, 'n_in',    'interactive'); % dimensionality of the Bin subspace to be used for the computation of the intersection
cfg.dssp.n_out        = ft_getopt(cfg.dssp, 'n_out',   'interactive'); % dimensionality of the Bout subspace to be used for the computation of the intersection
cfg.dssp.n_intersect  = ft_getopt(cfg.dssp, 'n_intersect', 'interactive'); % dimensionality of the intersection
cfg.output            = ft_getopt(cfg, 'output', 'original');

pertrial = istrue(cfg.pertrial);

% select channels and trials of interest, by default this will select all channels and trials
tmpcfg = keepfields(cfg, {'trials', 'channel', 'tolerance', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
datain = ft_selectdata(tmpcfg, datain);
% restore the provenance information
[cfg, datain] = rollback_provenance(cfg, datain);

if istrue(cfg.demean)
  ft_info('demeaning the time series');
  tmpcfg = [];
  tmpcfg.demean = 'yes';
  datain = ft_preprocessing(tmpcfg, datain);
  % restore the provenance information
  [cfg, datain] = rollback_provenance(cfg, datain);
end

% match the input data's channels with the labels in the leadfield
sourcemodel = cfg.sourcemodel;
if ~isfield(sourcemodel, 'leadfield')
  ft_error('cfg.sourcemodel needs to contain leadfields');
end
[indx1, indx2] = match_str(datain.label, sourcemodel.label);
if ~isequal(indx1(:),(1:numel(datain.label))')
  ft_error('unsupported mismatch between data channels and leadfields');
end
if islogical(sourcemodel.inside)
  inside = find(sourcemodel.inside);
else
  inside = sourcemodel.inside;
end
for k = inside(:)'
  sourcemodel.leadfield{k} = sourcemodel.leadfield{k}(indx2,:);
end

% compute the Gram-matrix of the supplied forward model
lf = cat(2, sourcemodel.leadfield{:});
G  = lf*lf';

%dat     = cat(2,datain.trial{:});
[Bclean, subspace] = dssp(datain.trial, G, cfg.dssp.n_in, cfg.dssp.n_out, cfg.dssp.n_space, cfg.dssp.n_intersect, pertrial);
% datAe   = datain.trial*cellfun(@transpose, Ae, 'UniformOutput', false); % the projection is a right multiplication
% with a matrix (eye(size(Ae,1))-Ae*Ae'), since Ae*Ae' can become quite
% sizeable, it's computed slightly differently here.

% put some diagnostic information in the output cfg.
cfg.dssp.subspace = subspace;

% replace the input cfg values
cfg.dssp.n_space = subspace.S(1).n;
cfg.dssp.n_in    = subspace.Sin(1).n;
cfg.dssp.n_out   = subspace.Sout(1).n;
cfg.dssp.n_intersect = subspace.T(1).n;

% compute the cleaned data and put in a cell-array
switch cfg.output
  case 'original'
    trial = Bclean;
  case 'complement'
    trial = datain.trial-Bclean;
  otherwise
    ft_error(sprintf('cfg.output = ''%s'' is not implemented',cfg.output));
end

% create the output argument
dataout       = keepfields(datain, {'label', 'time', 'fsample', 'trialinfo', 'sampleinfo', 'grad', 'elec', 'opto'}); % grad can be kept and does not need to be balanced, since the cleaned data is a mixture over time, not space.
dataout.trial = trial;

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   datain
ft_postamble provenance dataout
ft_postamble history    dataout
ft_postamble savevar    dataout

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions for the computation of the projection matrix
% kindly provided by Kensuke, and adjusted a bit by Jan-Mathijs
function [Bclean, subspace] = dssp(B, G, Nin, Nout, Nspace, Nintersect, pertrial)

% Nc: number of sensors
% Nt: number of time points
% inputs
% B(Nc,Nt):  interference overlapped sensor data
% G(Nc,Nc): Gram matrix of voxel lead field
% Nout and Nin: dimensions of the two row spaces
% recom_Nspace: recommended value for the dimension of the pseudo-signal subspace
% outputs
% Bclean(Nc,Nt): cleaned sensor data
% Nintersect: dimension of the intersection
% Nspace: dimension of the pseudo-signal subspace
%  ------------------------------------------------------------
%  programmed by K. Sekihara,  Signal Analysis Inc.
%  All right reserved by Signal Analysis Inc.
% -------------------------------------------------------------
%
% The code below is modified by Jan-Mathijs, no functional changes
% merely cosmetics, added the possibility to run the temporal subspace per
% trial

% eigen decomposition of the Gram matrix, matrix describing the spatial components of the defined 'in' compartment
fprintf('Computing the spatial subspace projection\n');
fprintf('Eigenvalue decomposition of the Gram matrix\n');
[Uspace,S] = eig(G);
Sspace     = abs(diag(S));

[Sspace, iorder] = sort(-Sspace);
Sspace           = -Sspace;
Uspace(:,:)      = Uspace(:,iorder);
Nspace           = getN(Nspace, Sspace, 'spatial');

% spatial subspace projector
Us   = Uspace(:,1:Nspace);

%USUS = Us*Us';
% % Bin and Bout creation
%Bin  =                  USUS  * B;
%Bout = (eye(size(USUS))-USUS) * B;

% computationally more efficient than the above
fprintf('Applying the spatial subspace projector\n');
Bin  = Us*(Us'*B);
Bout = B - Bin;

fprintf('Computing the subspace projector based on signal correlations\n');
[Ae, subspace] = CSP01(Bin, Bout, Nin, Nout, Nintersect, pertrial);

% add the first spatial subspace projection information as well
subspace.S.U = Uspace;
subspace.S.S = Sspace;
subspace.S.n = Nspace;
subspace.trial = Ae;

fprintf('Applying the subspace projector\n');
%Bclean    = B - (B*Ae)*Ae'; % avoid computation of Ae*Ae'
Bclean = B - (B*cellfun(@transpose, Ae, 'UniformOutput', false))*Ae;

function [Ae, subspace] = CSP01(Bin, Bout, Nin, Nout, Nintersect, pertrial)
%
% interference rejection by removing the common temporal subspace of the two subspaces
% K. Sekihara,  March 28, 2012
% Golub and Van Loan, Matrix computations, The Johns Hopkins University Press, 1996
%
%  Nc: number of channels
%  Nt: number of time points
% inputs
%  Bout(1:Nc,1:Nt): interference data
%  Bin(1:Nc,1:Nt): signal plus interference data
%  Nout: dimension of the interference subspace
%  Nin: dimension of the signal plus interference subspace
%  Nintersect: dimension of the intersection of the two subspaces
% outputs
% Ae = matrix from which the projector onto the intersection can
%      be obtained:
% subspace: struct containing information about the different subspace
%      projections
%  ------------------------------------------------------------
%  programmed by K. Sekihara,  Signal Analysis Inc.
%  All right reserved by Signal Analysis Inc.
% -------------------------------------------------------------
%

if ~pertrial
  % compute the projection across trials
  trllist = 1:numel(Bout);
else
  % compute the projection per trial
  trllist = (1:numel(Bout))';
end

Ae = cell(size(Bin));
for k = 1:size(trllist,1)
  indx = trllist(k,:); % this is either a scalar, or a vector
  [Uout,Sout,Vout] = svd(cat(2, Bout{indx}),'econ');
  [Uin, Sin, Vin]  = svd(cat(2, Bin{indx}), 'econ');
  Sout = diag(Sout);
  Sin  = diag(Sin);

  Nout = getN(Nout, Sout, 'outside');
  Nin  = getN(Nin,  Sin,  'inside');

  % compute unit-norm orthogonal time courses
  Qout = diag(1./Sout(1:Nout))*Uout(:,1:Nout)'*Bout(indx); % keep it in cell representation
  Qin  = diag(1./Sin(1:Nin)  )* Uin(:,1:Nin)' *Bin(indx);
  C    = Qin * cellfun(@transpose, Qout, 'UniformOutput', false);
  C    = sum(cat(3, C{:}), 3);

  % store the subspace information that is used in the next step
  subspace.Sin(k).U  = Uin;
  subspace.Sin(k).S  = Sin;
  subspace.Sin(k).n  = Nin;
  subspace.Sout(k).U = Uout;
  subspace.Sout(k).S = Sout;
  subspace.Sout(k).n = Nout;

  % covariance matrix of unit-norm 'components' -> how does this relate to
  % multivariate decomp? This is I guess equivalent mathematically
  [U,S] = svd(C);
  S     = diag(S);
  Nintersect = getN(Nintersect, S, 'intersection');

  Ae(indx) = U(:, 1:Nintersect)'*Qin;

  % keep the subspace information
  subspace.T(k).U = U;
  subspace.T(k).S = S;
  subspace.T(k).C = C; clear C;
  subspace.T(k).n = Nintersect;
end

function N = getN(N, S, name)

ttext = sprintf('enter the dimension for the %s field: ', name);
if isempty(N)
  N  = input(ttext);
elseif ischar(N) && isequal(N, 'interactive') && ~any(strcmp(name, {'outside' 'intersection'}))
  figure, plot(log10(S),'-o'); drawnow
  N = input(ttext);
elseif ischar(N) && isequal(N, 'interactive') && any(strcmp(name, {'outside' 'intersection'}))
  figure, plot(S, '-o'); drawnow
  N = input(ttext);
elseif ischar(N) && isequal(N, 'all')
  N = find(S./S(1)>1e5*eps, 1, 'last');
elseif isnumeric(N) && N<1
  N = find(S<=N, 1, 'last');
end
fprintf('Using %d dimensions for the %s field\n', N, name);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%DEMO_BAYESIANOPTIMIZATION3  A demonstration program for Bayesian
%                            optimization with constraints
%
% The set of BO demos
%  Part 1: see demo_bayesoptimization1
%  One dimensional example
%
%  Part 2: see demo_bayesoptimization3
%  Two dimensional example
%
%  Part 3: this file
%  Two dimensional example with constraints
%  * The implementation of constraints follows Gelbart et al. (2014)
%
%  References:
%    Jones, D., Schonlau, M., & Welch, W. (1998). Efficient global
%    optimization of expensive black-box functions. Journal of Global
%    Optimization, 13(4), 455-492. doi:10.1023/a:1008306431147
%
%    Michael A. Gelbart, Jasper Snoek, and Ryan P. Adams
%    (2014). Bayesian Optimization with Unknown Constraints.
%    http://arxiv.org/pdf/1403.5607v1.pdf
%
%    Snoek, J., Larochelle, H, Adams, R. P. (2012). Practical Bayesian
%    Optimization of Machine Learning Algorithms. NIPS 25
%
%  Copyright (c) 2015 Jarno Vanhatalo

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

%%  Part 3:
%  Two dimensional example with constraints
% For testing purposes:
stack = dbstack;
if (~isempty(stack) && (strcmp(stack(end).name, 'runtestset') || strcmp(stack(end).name, 'runtests'))) test = 1; else test = 0; end;

rng(3)
% Construct function handles to objective function (fx) and two constraint
% functions (fxc, fxc2)
fx = @(x) -log( (mvnpdf([x(:,1) x(:,2)],[3.5 2.5], [1 0.3; 0.3 1]) + 0.3*mvnpdf([x(:,1) x(:,2)],[7 8], [3 0.5; 0.5 4])).*...
    mvnpdf([x(:,1) x(:,2)],[5 5], [100 0; 0 100])) ./15 -1;
fxc = @(x) ((x(:,1)-5) .^2 + (x(:,2)-5).^2 -1)/30;
fxc2 = @(x) ( (x(:,1)-5) .^2)./30 - 0.5;

% The upper and lower limits for the constraints
const = [0 0.8 ; -10 0.1];

% Help variables for visualization
lb=0;
ub=10;
[X,Y] = meshgrid(linspace(lb,ub,100),linspace(lb,ub,100));
xl = [X(:) Y(:)];
Z = reshape(fx(xl),100,100);
Zc1 = fxc(xl); Zc1(Zc1<const(1,1) | Zc1>const(1,2)) = nan;
Zc1(~isnan(Zc1))=1; Zc1 = reshape(Zc1,100,100);
Zc2 = fxc2(xl); Zc2(Zc2<const(2,1) | Zc2>const(2,2)) = nan;
Zc2(~isnan(Zc2))=1; Zc2 = reshape(Zc2,100,100);

% ----- conduct Bayesian optimization -----

% construct GP models for the objective function and constraint functions
cfc = gpcf_constant('constSigma2',10,'constSigma2_prior', prior_fixed);
cfse = gpcf_sexp('lengthScale',[1 1]);
cfl = gpcf_linear('coeffSigma2', 10);
cfl2 = gpcf_squared('coeffSigma2', 10, 'interactions', 'on');
lik = lik_gaussian('sigma2', 0.001, 'sigma2_prior', prior_fixed);
% GP model for objective function
%gp1 = gp_set('cf', {cfc, cfl, cfl2, cfse}, 'lik', lik);cfl, cfl2,
gp1 = gp_set('cf', {cfc, cfse}, 'lik', lik);
% GP models for constraint functions
gpc1 = {gp_set('cf', {cfc, cfse}, 'lik', lik, 'jitterSigma2', 1e-6),...
    gp_set('cf', {cfc, cfse}, 'lik', lik, 'jitterSigma2', 1e-6)};

% Set the options for optimizer of the acquisition function
optimf = @fmincon;
optdefault=struct('GradObj','on','LargeScale','on','Algorithm','interior-point','TolFun',1e-9,'TolX',1e-6, 'Display', 'iter');
opt=optimset(optdefault);
lb=[0 0];     % lower bound of the input space
ub=[10 10];   % upper bound of the input space

% draw initial points
% we assume that at first we don't have any observation from the objective
% function but only observations from constraint functions. Hence x and y
% are initialized to zero
x = [];
y = [];
xc1 = 10*rand(2,2);
yc1 = fxc(xc1);
xc2 = 10*rand(5,2);
yc2 = fxc2(xc2);

figure, % figure for visualization
i1 = 1;
maxiter = 25;
improv = inf;   % improvement between two successive query points
while i1 < maxiter && improv>1e-6
%while i1 < maxiter
    % Train the GP models and calculate variables that are needed when
    %   calculating the Expected improvement (Acquisition function)
    % Objective function
    if ~isempty(x)
        gp = gp_optim(gp1,x,y);
        [K, C] = gp_trcov(gp,x);
        invC = inv(C);
        a = C\y;
        fmin = min( fx(x) );
    else
        a=[];
        x=[];
        invC=[];
        fmin=[];
        gp=gp1;
    end
    % constrain function 1
    gpct = gp_optim(gpc1{1},xc1,yc1);
    [~, Cct] = gp_trcov(gpct,xc1);
    const1.gpc = gpct;
    const1.invCc = inv(Cct);
    const1.ac = Cct\yc1;
    const1.const = const(1,:);
    const1.xc = xc1;
    % constrain function 2
    gpct = gp_optim(gpc1{2},xc2,yc2);
    [~, Cct] = gp_trcov(gpct,xc2);
    const2.gpc = gpct;
    const2.invCc = inv(Cct);
    const2.ac = Cct\yc2;
    const2.const = const(2,:);
    const2.xc = xc2;

    % Calculate EI and the posterior of the functions for visualization
    if ~isempty(x)
        [Ef,Varf] = gp_pred(gp, x, y, xl);
        EI = expectedimprovement_eg(xl, gp, x, a, invC, fmin, const1, const2);
    else
        Ef = zeros(size(xl,1),1);
        Varf = zeros(size(xl,1),1);
        EI = zeros(size(xl,1),1);
    end
    [Efc1] = gp_pred(const1.gpc, xc1, yc1, xl);
    [Efc2] = gp_pred(const2.gpc, xc2, yc2, xl);

    % optimize acquisition function
    %  * Note! Opposite to the standard notation we minimize negative Expected
    %    Improvement since Matlab optimizers seek for functions minimum
    %  * Note! We alternate the acquisition function between Expected
    %    Improvement and expected variance. The latter helps the
    %    optimization so that it does not get stuck in local mode
    % Here we use multiple starting points for the optimization so that we
    % don't crash into suboptimal mode of acquisition function
%     if mod(i1,5)==0  %Do just exploration by finding the maimum variance location
%         fh_eg = @(x_new) expectedvariance_eg(x_new, gp, x, [], invC);
%     else
        fh_eg = @(x_new) expectedimprovement_eg(x_new, gp, x, a, invC, fmin, const1, const2);
%     end
    nstarts = 20;
    xstart = [repmat(lb,nstarts,1) + repmat(ub-lb,nstarts,1).*rand(nstarts,2) ]; %; repmat(x(indbest,:),2,1)+0.1*randn(2,size(x,2))
    for s1=1:nstarts
        x_new(s1,:) = optimf(fh_eg, xstart(s1,:), [], [], [], [], lb, ub, [], opt);
    end
    xnews = x_new;
    EIs = fh_eg(x_new);
    x_new = x_new( find(EIs==min(EIs),1), : );

    % New sample point
    x(end+1,:) = x_new;
    y(end+1,:) = fx(x(end,:));
    xc1(end+1,:) = x_new;
    yc1(end+1,:) = fxc(x(end,:));
    xc2(end+1,:) = x_new;
    yc2(end+1,:) = fxc2(x(end,:));

    % visualize
    clf
    % Plot the objective function
    subplot(2,4,1),hold on, title('Objective, query points')
    pcolor(X,Y,Z),shading flat
    clim = caxis;
    plot(x(1:end-1,1),x(1:end-1,2), 'rx', 'MarkerSize', 10),
    plot(x(end,1),x(end,2), 'ro', 'MarkerSize', 10, 'linewidth', 3)
    % Plot the posterior mean of the GP model
    subplot(2,4,2),hold on, title(sprintf('GP prediction, mean, iter: %d',i1))
    pcolor(X,Y,reshape(Ef,100,100)),shading flat
    caxis(clim)
    % Plot the posterior variance of GP model
    subplot(2,4,6),hold on, title('GP prediction, variance')
    pcolor(X,Y,reshape(Varf,100,100)),shading flat
    plot(xnews(:,1),xnews(:,2), 'ro', 'MarkerSize', 10)
    plot(x(end,1),x(end,2), 'ro', 'MarkerSize', 10, 'linewidth', 3)
    plot(x(1:end-1,1),x(1:end-1,2), 'rx', 'MarkerSize', 10),

    % The expected information
    subplot(2,4,5), hold on, title(sprintf('Expected improvement %.2e', min(EIs)))
    pcolor(X,Y,reshape(EI,100,100)),shading flat
    plot(xnews(:,1),xnews(:,2), 'ro', 'MarkerSize', 10)
    plot(x(end,1),x(end,2), 'ro', 'MarkerSize', 10, 'linewidth', 3)
    plot(x(1:end-1,1),x(1:end-1,2), 'rx'),

    % constraint 1
    subplot(2,4,3), hold on, title(sprintf('constraint 1'))
    pcolor(X,Y,Zc1),shading flat
    plot(xc1(1:end-1,1),xc1(1:end-1,2), 'rx', 'MarkerSize', 10);
    plot(xnews(:,1),xnews(:,2), 'ro', 'MarkerSize', 10);
    plot(xc1(end,1),xc1(end,2), 'ro', 'MarkerSize', 10, 'linewidth', 3);
    % constraint 2
    subplot(2,4,4), hold on, title(sprintf('constraint 2'))
    pcolor(X,Y,Zc2),shading flat
    l1= plot(xc2(1:end-1,1),xc2(1:end-1,2), 'rx', 'MarkerSize', 10);
    l2=plot(xnews(:,1),xnews(:,2), 'ro', 'MarkerSize', 10);
    l3=plot(xc2(end,1),xc2(end,2), 'ro', 'MarkerSize', 10, 'linewidth', 3);
    legend([l1,l2,l3], {'function evaluation points','local modes of acquisition function','The next query point'})

    % prediction of constraint 1
    subplot(2,4,7), hold on, title(sprintf('prediction for const 1'))
    Efc1(Efc1<const(1,1) | Efc1>const(1,2)) = nan;
    Efc1(~isnan(Efc1))=1; Efc1 = reshape(Efc1,100,100);
    pcolor(X,Y,reshape(Efc1,100,100)),shading flat
    plot(xnews(:,1),xnews(:,2), 'ro', 'MarkerSize', 10)
    plot(xc1(end,1),xc1(end,2), 'ro', 'MarkerSize', 10, 'linewidth', 3)
    plot(xc1(1:end-1,1),xc1(1:end-1,2), 'rx'),
    % prediction of constraint 2
    subplot(2,4,8), hold on, title(sprintf('prediction for const 2'))
    Efc2(Efc2<const(2,1) | Efc2>const(2,2)) = nan;
    Efc2(~isnan(Efc2))=1; Efc2 = reshape(Efc2,100,100);
    pcolor(X,Y,reshape(Efc2,100,100)),shading flat
    plot(xc2(1:end-1,1),xc2(1:end-1,2), 'rx', 'MarkerSize', 10),
    plot(xnews(:,1),xnews(:,2), 'ro', 'MarkerSize', 10)
    plot(xc2(end,1),xc2(end,2), 'ro', 'MarkerSize', 10, 'linewidth', 3)

    if length(y)>1
        improv = abs(y(end) - y(end-1));
    end
    i1=i1+1;

    if test == 0
        pause
    end
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function str = showdetail(item)

% function str = showdetail(item)
% Display details for a cfg_choice and all of its options.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: showdetail.m 4867 2012-08-30 13:04:51Z volkmar $

rev = '$Rev: 4867 $'; %#ok

str = showdetail(item.cfg_item);
str{end+1} = 'Set the options listed as a struct:';
% Display detailed help for each cfg_choice value item
tags = tagnames(item, true);
for k = 1:numel(tags)
    str{end+1} = sprintf('.%s', tags{k});
end;

author: spm


text:
%% clear the workspace and select data
clear; clc; close all;

%% choose data
neuron = Sources2D();
nam = './data_2p.tif';          % this demo data is very small, here we just use it as an example
nam = neuron.select_data(nam);  %if nam is [], then select data interactively

%% parameters
% -------------------------    COMPUTATION    -------------------------  %
pars_envs = struct('memory_size_to_use', 8, ...   % GB, memory space you allow to use in MATLAB
    'memory_size_per_patch', 0.3, ...   % GB, space for loading data within one patch
    'patch_dims', [30, 40]);  %GB, patch size

% -------------------------      SPATIAL      -------------------------  %
gSig = 0.5;           % pixel, gaussian width of a gaussian kernel for filtering the data. 0 means no filtering
gSiz = 10;          % pixel, neuron diameter
ssub = 1;           % spatial downsampling factor
with_dendrites = false;   % with dendrites or not
if with_dendrites
    % determine the search locations by dilating the current neuron shapes
    updateA_search_method = 'dilate';  %#ok<UNRCH>
    updateA_bSiz = 20;
    updateA_dist = neuron.options.dist;
else
    % determine the search locations by selecting a round area
    updateA_search_method = 'ellipse'; %#ok<UNRCH>
    updateA_dist = 5;
    updateA_bSiz = neuron.options.dist;
end
spatial_constraints = struct('connected', true, 'circular', false);  % you can include following constraints: 'circular'

% -------------------------      TEMPORAL     -------------------------  %
Fs = 5;             % frame rate
tsub = 1;           % temporal downsampling factor
deconv_options = struct('type', 'ar1', ... % model of the calcium traces. {'ar1', 'ar2'}
    'method', 'foopsi', ... % method for running deconvolution {'foopsi', 'constrained', 'thresholded'}
    'smin', -3, ...         % minimum spike size. When the value is negative, the actual threshold is abs(smin)*noise level
    'optimize_pars', true, ...  % optimize AR coefficients
    'optimize_b', true);% optimize the baseline);
nk = 1;             % detrending the slow fluctuation. usually 1 is fine (no detrending)
                    % when changed, try some integers smaller than total_frame/(Fs*30)
detrend_method = 'local_min';  % compute the local minimum as an estimation of trend.

% -------------------------     BACKGROUND    -------------------------  %
bg_model = 'svd';  % model of the background {'ring', 'svd'(default), 'nmf'}
nb = 1;             % number of background sources for each patch (only be used in SVD and NMF model)
ring_radius = 13;  % when the ring model used, it is the radius of the ring used in the background model.
                    %otherwise, it's just the width of the overlapping area

% -------------------------      MERGING      -------------------------  %
show_merge = false;  % if true, manually verify the merging step
merge_thr = 0.7;     % thresholds for merging neurons; [spatial overlap ratio, temporal correlation of calcium traces, spike correlation]
method_dist = 'max';   % method for computing neuron distances {'mean', 'max'}
dmin = 10;       % minimum distances between two neurons. it is used together with merge_thr
dmin_only = 2;  % merge neurons if their distances are smaller than dmin_only.

% -------------------------  INITIALIZATION   -------------------------  %
K = [];             % maximum number of neurons per patch. when K=[], take as many as possible
min_corr = 0.7;     % minimum local correlation for a seeding pixel
min_pnr =10;       % minimum peak-to-noise ratio for a seeding pixel
min_pixel = 2^2;      % minimum number of nonzero pixels for each neuron
bd = 0;             % number of rows/columns to be ignored in the boundary (mainly for motion corrected data)
frame_range = [];   % when [], uses all frames
save_initialization = false;    % save the initialization procedure as a video.
use_parallel = true;    % use parallel computation for parallel computing
show_init = true;   % show initialization results
choose_params = false; % manually choose parameters
center_psf = false;  % set the value as true when the background fluctuation is large (usually 1p data)
                    % set the value as false when the background fluctuation is small (2p)

% ----------------------   MANUAL INTERVENTION  --------------------  %
with_manual_intervention = false;

% -------------------------  FINAL RESULTS   -------------------------  %
save_demixed = true;    % save the demixed file or not
kt = 1;                 % frame intervals

% -------------------------    UPDATE ALL    -------------------------  %
neuron.updateParams('gSig', gSig, ...       % -------- spatial --------
    'gSiz', gSiz, ...
    'ring_radius', ring_radius, ...
    'ssub', ssub, ...
    'search_method', updateA_search_method, ...
    'bSiz', updateA_bSiz, ...
    'dist', updateA_bSiz, ...
    'spatial_constraints', spatial_constraints, ...
    'tsub', tsub, ...                       % -------- temporal --------
    'deconv_options', deconv_options, ...    '
    'nk', nk, ...
    'detrend_method', detrend_method, ...
    'background_model', bg_model, ...       % -------- background --------
    'nb', nb, ...
    'ring_radius', ring_radius, ...
    'merge_thr', merge_thr, ...             % -------- merging ---------
    'dmin', dmin, ...
    'method_dist', method_dist, ...
    'min_corr', min_corr, ...               % ----- initialization -----
    'min_pnr', min_pnr, ...
    'min_pixel', min_pixel, ...
    'bd', bd, ...
    'center_psf', center_psf);
neuron.Fs = Fs;

%% distribute data and be ready to run source extraction
neuron.getReady(pars_envs);

%% initialize neurons from the video data within a selected temporal range
if choose_params
    % change parameters for optimized initialization
    [gSig, gSiz, ring_radius, min_corr, min_pnr] = neuron.set_parameters();
end

[center, Cn, PNR] = neuron.initComponents_parallel(K, frame_range, save_initialization, use_parallel);
if show_init
    figure;
    imagesc(Cn, [0, 1]); colormap gray;
    hold on;
    plot(center(:, 2), center(:, 1), '.r', 'markersize', 10);
end
neuron_init = neuron.copy();

%% estimation of background components
neuron.update_background_parallel(use_parallel);

%% pick neurons from the residual
[center_res, Cn_res, PNR_res] =neuron.initComponents_residual_parallel(K, save_initialization, use_parallel);
neuron_init_res = neuron.copy();

%% estimate the temporal components
neuron.update_temporal_parallel(use_parallel);

%% estimate the spatial components
neuron.update_spatial_parallel(use_parallel);

%% deal with bad neurons
tags = neuron.tag_neurons_parallel();  % find neurons with fewer nonzero pixels than min_pixel and silent calcium transients
neuron.delete(tags>0);

%% merge neurons
neuron.merge_neurons_dist_corr(show_merge);
neuron.merge_close_neighbors(show_merge, dmin_only);

%% add a manual intervention and run the whole procedure for a second time
neuron.orderROIs('snr');   % order neurons in different ways {'snr', 'decay_time', 'mean', 'circularity'}
if with_manual_intervention
    neuron.viewNeurons([], neuron.C_raw);
end
neuron.update_background_parallel(use_parallel);
neuron.update_temporal_parallel(use_parallel);
neuron.update_spatial_parallel(use_parallel);

K = size(neuron.A,2);
tags = neuron.tag_neurons_parallel();  % find neurons with fewer nonzero pixels than min_pixel and silent calcium transients
neuron.delete(tags>0);
neuron.merge_neurons_dist_corr(show_merge);
neuron.merge_close_neighbors(show_merge, dmin_only);
if K~=size(neuron.A,2)
    neuron.update_temporal_parallel(use_parallel);
    neuron.update_spatial_parallel(use_parallel);
    neuron.update_background_parallel(use_parallel);
end

%% save the workspace for future analysis
neuron.save_workspace();

%% show neuron contours
Coor = neuron.show_contours();

%% create a video for displaying the
amp_ac = 5000;
range_ac = [0, 5000];
range_Y = [0, 10000];
neuron.show_demixed_video(save_demixed, kt, [], amp_ac, range_ac, range_Y);

%% save neurons shapes
neuron.save_neurons();

%% save neurons shapes
neuron.save_neurons();

author: zhoupc


text:
function val=geometricMean(x,dim)
%%GEOMETRICMEAN Compute the geometric mean of a set of samples.
%
%INPUTS: x The set of samples over which the geometric mean is to be found.
%          How the function handles vectors and matrices is determined by
%          the parameter dim.
%      dim An optional parameter specifying over which dimension of x the
%          geometric mean is to be taken. If dim is omitted or an empty
%          matrix is passsed, then if x is a vector, the geometric mean is
%          found over all elements of the vector. If x is a matrix, then
%          the mean is found over the columns of the matrix (resulting in a
%          row vector), and if x is an n-dimensional (n>2) matrix, then the
%          geometric mean is found over the first non-singleton dimension
%          of the matrix.
%
%OUTPUT: val The geometric mean of x taken over the appropriate dimension.
%
%The geometric mean of n items is the nth-root of the product of those
%items. The computation is performed using logarithms to avoid overflows.
%This can mean that one does not obtain the desired complex nth root when
%negative values are complex numbers are the argument as there are multiple
%solutions to the nth root.
%
%October 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(isempty(x))
    val=[];
    return
end

xDims=size(x);
if(nargin<2||isempty(dim))
    if(length(xDims)<3)%If x is a vector or matrix.
        if(xDims(1)==1||xDims(2)==1)%If x is a vector
            x=x(:);
            dim=1;
            n=length(x);
        else%x is a matrix, go over the columns.
            dim=2;
            n=xDims(2);
        end
    else%If x is an n-dimensional matrix.
        dim=find(xDims>1,1);
        if(isempty(dim))
            dim=1;
            n=1;
        else
            n=xDims(dim);
        end
    end
else
    n=xDims(dim);
end

%We could do prod(x,dim).^(1/n);, but that is likely to overflow if there
%are many values. Thus, we take the exponent of the sum of the logarithm.
val=exp((1/n)*sum(log(x),dim));

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function varargout = double_pendulum2(varargin)

% the function simulate double pendulum, using 4th-order Runge-Kutta algorithm for the diffrential equations.
% the diffrential equation are very similar to eq. (12) and (13) here:
%                           http://home1.fvcc.edu/~dhicketh/DiffEqns/spring09projects/LauraStickel/Double%20Pendulum.pdf
% where omega=d(theta)/dt
% Moshe Lindner & Orit Peleg, August 2010 (C)

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @double_pendulum_OpeningFcn, ...
                   'gui_OutputFcn',  @double_pendulum_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 double_pendulum is made visible.
function double_pendulum_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to double_pendulum (see VARARGIN)

% Choose default command line output for double_pendulum
handles.output = hObject;
handles.len1=.5;
handles.len2=.5;
handles.mass1=.5;
handles.mass2=.5;
handles.t1=0;
handles.t2=0;
handles.w1=0;
handles.w2=0;
handles.G=9.8;
pend_plot(handles);

% Update handles structure
guidata(hObject, handles);

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

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

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

% --- Executes on button press in start_stop.
function start_stop_Callback(hObject, eventdata, handles)
% hObject    handle to start_stop (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
set(hObject,'string','Stop');
set(handles.omega1,'enable','off');
set(handles.theta1,'enable','off');
set(handles.omega2,'enable','off');
set(handles.theta2,'enable','off');
set(handles.length_rat,'enable','off');
set(handles.mass_rat,'enable','off');
set(handles.reset,'enable','off');
set(handles.gravity,'enable','off');
i=1;
while get(hObject,'value')==1
    if (get(handles.traj,'value'))==1
        x_1=handles.len1*sin(handles.t1);
        y_1=-handles.len1*cos(handles.t1) ;
        x_2=handles.len2*sin(handles.t2) + x_1;
        y_2=-handles.len2*cos(handles.t2) + y_1;
        handles.trajectory(1,i)=x_2;
        handles.trajectory(2,i)=y_2;
        i=i+1;
        if i>400
            i=1;
        end
    elseif i~=1
        handles.trajectory=[];
        i=1;
    end
    [handles.t1,handles.t2,handles.w1,handles.w2]=db_pendulum(handles.t1,handles.t2,handles.w1,handles.w2,handles.mass1,handles.mass2,handles.len1,handles.len1,handles.G,.01);
    pend_plot(handles);
    set(handles.omega1,'string',num2str(handles.w1));
    set(handles.theta1,'string',num2str(handles.t1));
    set(handles.omega2,'string',num2str(handles.w2));
    set(handles.theta2,'string',num2str(handles.t2));
    guidata(hObject, handles);
end
set(hObject,'string','Start');
set(handles.omega1,'enable','on');
set(handles.theta1,'enable','on');
set(handles.omega2,'enable','on');
set(handles.theta2,'enable','on');
set(handles.length_rat,'enable','on');
set(handles.mass_rat,'enable','on');
set(handles.reset,'enable','on');
set(handles.gravity,'enable','on');
% --- Executes on slider movement.
function length_rat_Callback(hObject, eventdata, handles)
% hObject    handle to length_rat (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.len1=get(hObject,'value');
handles.len2=1-handles.len1;
pend_plot(handles);
guidata(hObject, handles);
% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

% --- Executes during object creation, after setting all properties.
function length_rat_CreateFcn(hObject, eventdata, handles)
% hObject    handle to length_rat (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function mass_rat_Callback(hObject, eventdata, handles)
% hObject    handle to mass_rat (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.mass1=get(hObject,'value');
handles.mass2=1-handles.mass1;
pend_plot(handles);
guidata(hObject, handles);
% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

% --- Executes during object creation, after setting all properties.
function mass_rat_CreateFcn(hObject, eventdata, handles)
% hObject    handle to mass_rat (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

function omega1_Callback(hObject, eventdata, handles)
% hObject    handle to omega1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.w1=str2num(get(hObject,'String'));
energy(handles)
guidata(hObject, handles);
% Hints: get(hObject,'String') returns contents of omega1 as text
%        str2double(get(hObject,'String')) returns contents of omega1 as a double

% --- Executes during object creation, after setting all properties.
function omega1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to omega1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function theta1_Callback(hObject, eventdata, handles)
% hObject    handle to theta1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.t1=str2num(get(hObject,'String'));
pend_plot(handles);
guidata(hObject, handles);
% Hints: get(hObject,'String') returns contents of theta1 as text
%        str2double(get(hObject,'String')) returns contents of theta1 as a double

% --- Executes during object creation, after setting all properties.
function theta1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to text_1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function omega2_Callback(hObject, eventdata, handles)
% hObject    handle to omega2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.w1=str2num(get(hObject,'String'));
energy(handles)
guidata(hObject, handles);
% Hints: get(hObject,'String') returns contents of omega2 as text
%        str2double(get(hObject,'String')) returns contents of omega2 as a
%        double

% --- Executes during object creation, after setting all properties.
function omega2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to omega2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function theta2_Callback(hObject, eventdata, handles)
% hObject    handle to theta2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.t2=str2num(get(hObject,'String'));
pend_plot(handles);
guidata(hObject, handles);
% Hints: get(hObject,'String') returns contents of theta2 as text
%        str2double(get(hObject,'String')) returns contents of theta2 as a double

% --- Executes during object creation, after setting all properties.
function theta2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to theta2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in traj.
function traj_Callback(hObject, eventdata, handles)
% hObject    handle to traj (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.trajectory=[];
guidata(hObject, handles);
% Hint: get(hObject,'Value') returns toggle state of traj

% --- Executes on button press in reset.
function reset_Callback(hObject, eventdata, handles)
% hObject    handle to reset (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.len1=.5;
handles.len2=.5;
handles.mass1=.5;
handles.mass2=.5;
handles.t1=0;
handles.t2=0;
handles.w1=0;
handles.w2=0;
handles.G=9.8;
set(handles.gr,'string','9.8')
set(handles.gravity,'value',9.8)
set(handles.omega1,'string','0')
set(handles.omega2,'string','0')
set(handles.theta1,'string','0')
set(handles.theta2,'string','0')
set(handles.mass_rat,'value',0.5)
set(handles.length_rat,'value',0.5)
pend_plot(handles);
guidata(hObject, handles);

function pend_plot(data);
[xx yy zz]=cylinder(0.1,20);
    x_1=data.len1*sin(data.t1);
    y_1=-data.len1*cos(data.t1) ;
    x_2=data.len2*sin(data.t2) + x_1;
    y_2=-data.len2*cos(data.t2) + y_1;
hold off
plot([0 x_1 x_2],[0 y_1 y_2],'k','linewidth',2);
hold on
if get(data.traj,'value') ==1
    plot(data.trajectory(1,:),data.trajectory(2,:),'r.','markersize',3)
end
fill(x_1+xx(1,:).*(data.mass1)^.3,y_1+yy(1,:).*(data.mass1)^.3,'b')
fill(x_2+xx(1,:).*(data.mass2)^.3,y_2+yy(1,:).*(data.mass2)^.3,'b')
set(gca,'xtick',[],'ytick',[])
axis(1.2*[-1 1 -1 1]);
box on;
drawnow
energy(data)

function [theta_1 theta_2 omega_1 omega_2]=db_pendulum(theta1, theta2, w1, w2, mass1, mass2, d1, d2, G, dt)
l1 = dt * (l(w1));
m1 = dt * (m(w2));
f1 = dt * (f(theta1, theta2, w1, w2, mass1, mass2, d1, d2, G));
g1 = dt * (g(theta1, theta2, w1, w2, mass1, mass2, d1, d2, G));

l2 = dt * (l(w1+(f1/2.)));
m2 = dt * (m(w2+(g1/2.)));
f2 = dt * (f(theta1+(l1/2.),theta2+(m1/2.),w1+(f1/2.),w2+(g1/2.), mass1, mass2, d1, d2, G));
g2 = dt * (g(theta1+(l1/2.),theta2+(m1/2.),w1+(f1/2.),w2+(g1/2.), mass1, mass2, d1, d2, G));

l3 = dt * (l(w1+(f2/2.)));
m3 = dt * (m(w2+(g2/2.)));
f3 = dt * (f(theta1+(l2/2.),theta2+(m2/2.),w1+(f2/2.),w2+(g2/2.), mass1, mass2, d1, d2, G));
g3 = dt * (g(theta1+(l2/2.),theta2+(m2/2.),w1+(f2/2.),w2+(g2/2.), mass1, mass2, d1, d2, G));

l4 = dt * (l(w1+f3));
m4 = dt * (m(w2+g3));
f4 = dt * (f(theta1+(l3),theta2+(m3),w1+(f3),w2+(g3), mass1, mass2, d1, d2, G));
g4 = dt * (g(theta1+(l3),theta2+(m3),w1+(f3),w2+(g3), mass1, mass2, d1, d2, G));
theta_1 = theta1 + (l1+(2.*l2)+(2.*l3)+l4)/6.;
theta_2 = theta2 + (m1+(2.*m2)+(2.*m3)+m4)/6.;
omega_1 = w1 + (f1+(2.*f2)+(2.*f3)+f4)/6.;
omega_2 = w2 + (g1+(2.*g2)+(2.*g3)+g4)/6.;

function f1=f(theta1, theta2, w1, w2, mass1, mass2, d1, d2, G)
f1= (  (( -d2*mass2*w2*w2*sin(theta1- theta2)/(d1*(mass1+mass2)) )  - ...
    (  mass2*w1*w1*cos(theta1-theta2)*sin(theta1- theta2)/(mass1+mass2) )  +...
    (  mass2*G*cos(theta1- theta2)*sin(theta2)/(d1*(mass1+mass2))  )  -...
    (  G*sin(theta1)/d1)...
    ) /  ( 1. - (mass2*cos(theta1- theta2)*cos(theta1- theta2)/(mass1+mass2)))  );
function g1=g(theta1, theta2 ,w1, w2, mass1, mass2, d1, d2, G)
g1=(( (-G*sin(theta2)/d2) + (mass2*w2*w2*cos(theta1-theta2)*sin(theta1-theta2)/(mass1+mass2)) +...
    (G*cos(theta1-theta2)*sin(theta1)/d2) + (d1*w1*w1*sin(theta1-theta2)/d2)  )  /  (1. - (mass2*cos(theta1-theta2)*cos(theta1-theta2)/(mass1+mass2))) );
function L1=l(w1)
L1=w1;
function M1=m(w2)
M1= w2;

% --- Executes on slider movement.
function gravity_Callback(hObject, eventdata, handles)
% hObject    handle to gravity (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.G=get(hObject,'value');
set(handles.gr,'string',num2str(handles.G));
guidata(hObject, handles);

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

% --- Executes during object creation, after setting all properties.
function gravity_CreateFcn(hObject, eventdata, handles)
% hObject    handle to gravity (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

function energy(data)
  Epot1=-data.mass1*data.G*data.len1*(cos(data.t1)-1.000) ;
  Epot2=-data.mass2*data.G*((cos(data.t1)-1.000)*data.len1+(cos(data.t2)-1.000)*data.len2);
  Ekin1=0.5*data.mass1*(data.w1*data.len1).^2;
  Ekin2=0.5*data.mass2*(   (data.w2*data.len2*cos(data.t2) +data.w1*data.len1*cos(data.t1)).^2 + (data.w2*data.len2*sin(data.t2) +data.w1*data.len1*sin(data.t1)).^2 );
  Etot=Ekin1+Ekin2+Epot1+Epot2;

  set(data.U1,'string',num2str(Epot1));
  set(data.U2,'string',num2str(Epot2));
  set(data.T1,'string',num2str(Ekin1));
  set(data.T2,'string',num2str(Ekin2));
  set(data.totE,'string',num2str(Etot));

repo: Aero-matlab