text:
% Mseq Toolbox
% Vesion 1.0 10-28-2001
% written by Giedrius Buracas,
% SNL-B, Salk Institute
%
%
% Purpose:
%
%  Allows to evaluate estimation efficiency of m-sequence and
%  rando-sequence based experimental designs for event-related
%  fMRI experiments
%
% Sequence generation
%
%  balancedRnd		Generates random sequences, with either
%			overlaping or non-overlaping events
%			Equal numbers of events of each type is assumed
%
%  mseq			Generates binary, ternary, or five level
%			m-sequences. A unique feature of this code is that
%			many m-sequences of given parameters
%			can be generated
%
%  cycorr		generates a cyclical autocorrelation
%			function of a given sequence. This is
%			the fastest way to test whether a given sequence
%			is an m-sequence
%
%  m2bin		converts a binary m-sequence [-1,1] to
%			a sequence of [0,1]
%
%  bin2m		converts a sequence of [0,1] to [-1,1]
%
% Calculation of estimation efficiency
%
%  makeEventMtrx	generates an event matrix from a matrix of
%			[0,1] whose each column gives timing for
%			each event type
%
%  efficiencyOFexpDesigns2 script that compares estimation
%			efficiency of random and m-seuence-based
%			experimental designs. In this case
%			simultaneously occuring events are permitted
%			(overlaping events)
%
%  effOFexpDesignsNoOverlap script that compares estimation
%			efficiency of random and m-seuence-based
%			experimental designs. In this case
%			simultaneously occuring events are not permitted
%			(non-overlaping events)
%
%
%  efficiencyOFexpDesignsCorrMtrx script that compares estimation
%			efficiency of random and m-seuence-based
%			experimental designs.
%			Efficiency is calculated with and without
%			fMRI noise
%
%
%
%

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function net = remove_layers_from_block(net, type)
  [names, ~] = find_layers_from_block(net, type);
  for i = 1:numel(names)
    layer = net.layers(net.getLayerIndex(names{i}));
    net.removeLayer(names{i});
    net.renameVar(layer.outputs{1}, layer.inputs{1}, 'quiet', true);
  end
end

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet


text:
function [bitstream,char_table]=make_bits(msg)
% make a test transmission from msg
% Copyright 2002-2010 The MathWorks, Inc.
[vctable char_table]=load_alpha;   % varicode table index to 1+abs('c')
bitstream = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];   % idle bits
for k=1:length(msg)
    index=1+abs(msg(k));        % get ascii index (interesting use of the abs() function)
    vc=char(vctable(index));    % look up varicode bits
    kof=length(bitstream);
    bitstream([kof+(1:2)])=[0,0]; % 2 consecutive zeros is a character break
    for i=1:length(vc)
        bitstream(kof+2+i)=str2num(vc(i));
    end;
end;

author: Sable

repo: mcbench-benchmarks


text:
function [assignments] = FMC_interpret(AllSals, numClusters, AllPs, A_D, beta)
%% Interpretation of FMC data and grain assignment

%% run interpret 3D

A_D = A_D|A_D';

[i,~] = find(A_D);
cs = [0 full(cumsum(sum(A_D)))];

NNlist = cell(size(A_D,1),1);
for k = 1:size(A_D,1)
   NNlist{k} = i(cs(k)+1:cs(k+1))';
end

%% sort by saliency
% sort all course nodes for all s by saliency from lowest to highest
[tmp,elements] = sort(AllSals,'ascend'); %elements is list of indices representing the sorted AllSals list
% all points assigned to own node Nx2 [assignment,confidence]
% pointData is euler1 euler2 euler3 x y z confidence
numPoints = numClusters(1);
assignmentsX=zeros(numPoints,1);  % assignment of each poing
assignmentsY=zeros(numPoints,1);  % probability of that assignment

%% s values of each node
% vector sAssign has s value of all nodes
numS = length(numClusters); %number of iterations of s
sAssign=ones(numClusters(1),1); %initiating sAssign
lengthsub=zeros(1,numS);

%1 spot in sAssign for each node, increasing number for each saliency set
for curS=2:numS
    sAssign=[sAssign; curS*ones(numClusters(curS),1)];
end

% build assignments
for curS = numS-4:numS
    vdisp(['Considering S = ', int2str(curS)]);

    %find nodes for current s
    elemS = elements(sAssign(elements) == curS)';   %nodes of curS
    sLen = numClusters(curS);
    numNodes = length(elemS);

    % for each s from current s to s=1, iteratively find which nodes
    % are inside of the node by multiplying by Ps
    backPs = AllPs{curS};
    for iterS=curS-1:-1:2
        backPs=AllPs{iterS}*backPs;
    end

    % find probability of points assigned to each node
    allUs = [backPs assignmentsY];

    %assign clusters to points based on the strongest probability of any
    %point in the cluster
    [Prob, NodeIndex] = max(allUs,[],2);
    %assign clusters to points based on the strongest probability of any
%     %point in the cluster
    allUs = allUs >= beta;

    for point = find((NodeIndex <= numNodes) & (Prob < 1))'
      nodePoints = allUs(:,NodeIndex(point));
      assignmentsX(nodePoints) = NodeIndex(point);
      assignmentsY(nodePoints) = Prob(point);
    end

end

%% Break up non-continuous (Michigan) clusters
rLimit = 15;   %Recursion limit (MATLAB tends to crash around 840)
%N.B. This runs considerably faster with a low recursion limit due to
%MATLAB's huge overhead. Better to make more calls than fewer, deeper calls
% set(0,'RecursionLimit',rLimit)

vdisp('Breaking non-continuous clusters')

flag = zeros(length(assignmentsX),1);
assignmentsN = zeros(length(assignmentsX),1);
contPoints = [];
for point = 2:length(flag)
    if assignmentsN(point) == 0
        oldClust = assignmentsX(point);
        newClust = point;
        if (flag(point)>1)     %if the point has been reserved, continue the same grouping
            newClust = flag(point);
        end
        assignmentsN(point) = newClust;
        rDepth = 1;

        [assignmentsN, flag, contPoints] = clusterBreaker(assignmentsN, point, assignmentsX, flag, NNlist, oldClust, newClust, contPoints, rDepth, rLimit-5);

        while isempty(contPoints) == 0
            pointC = contPoints(1);
            contPoints = contPoints(2:end);
            assignmentsN(pointC) = newClust;
            rDepth = 1;
            [assignmentsN, flag, contPoints] = clusterBreaker(assignmentsN, pointC, assignmentsX, flag, NNlist, oldClust, newClust, contPoints, rDepth, rLimit-5);
        end
    end
end

singlePoints = find(flag == -1 & sum(A_D,2)>1);
for pt = singlePoints'
assignmentsN(pt) = assignmentsN(NNlist{pt}(2));
end

assignmentsX = assignmentsN;
% set(0,'RecursionLimit',500)

assignments = [assignmentsX assignmentsY];

end

author: mtex-toolbox

repo: mtex

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


text:
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% TEMPLATE FOR THE COMPILEMODEL FCN (TO BE USED WITH THE STATIC GUI)

% Compile the Simulink model through Matlab command line.

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Compiling the Simulink model...')

pause(1.5)

try
    SIMULINK_MODEL_NAME([], [], [], 'compile')
    SIMULINK_MODEL_NAME([], [], [], 'term')

catch ME

    errorMessages = ME;
end

clc

disp('Compilation done.')

% warning about Simulink timing
warning('The model will anyways start with FEW SECONDS OF DELAY after pressing the ''Start Model'' button.')

author: robotology

repo: whole-body-controllers


text:
function [ ra, dec ] = xyz_to_radec ( p )

%*****************************************************************************80
%
%% XYZ_TO_RADEC converts (X,Y,Z) to right ascension/declination coordinates.
%
%  Discussion:
%
%    Given an XYZ point, compute its distance R from the origin, and
%    regard it as lying on a sphere of radius R, whose axis is the Z
%    axis.
%
%    The right ascension of the point is the "longitude", measured in hours,
%    between 0 and 24, with the X axis having right ascension 0, and the
%    Y axis having right ascension 6.
%
%    Declination measures the angle from the equator towards the north pole,
%    and ranges from -90 (South Pole) to 90 (North Pole).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P(3,1), the coordinates of a point in 3D.
%
%    Output, real RA, DEC, the corresponding right ascension
%    and declination.
%
  p_norm = sqrt ( sum ( p(1:3,1).^2 )  );

  phi = r8_asin ( p(3,1) / p_norm );

  if ( cos ( phi ) == 0.0 )
    theta = 0.0;
  else
    theta = r8_atan ( p(2,1), p(1,1) );
  end

  dec = radians_to_degrees ( phi );
  ra = radians_to_degrees ( theta ) / 15.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = styblinskitang(X)
% Styblinski-Tang function
%
%   STYBLINSKYTANG([x1, x2, ..., xn]) returns the value of the
%   Styblinski-Tang at the specified points. All [xi] may be vectors.
%   The search domain is
%
%               -5 < x_i < 5
%
%   The global minimum is
%
%               f(x1, x2, ..., xn) =
%               f(-2.903534, -2.903534, ..., -2.903534) = -39.16599 * n

% Author: Rody P.S. Oldenhuis
% Delft University of Technology
% E-mail: oldenhuis@dds.nl
% Last edited 20/Jul/2009

    % if no input is given, return dimensions, bounds and minimum
    if (nargin == 0)
        varargout{1} = inf;  % # dims
        varargout{2} = -5; % LB
        varargout{3} = +5; % LB
        varargout{4} = -2.903534018185960e+000; % solution
        varargout{5} = -3.916616570377142e+001; % *(#dims) .. function value at solution

    % otherwise, output function value
    else

        % keep all values in the search domain
        X(X < -5) = inf;  X(X > 5) = inf;

        % NOTE: orientation can not be determined automatically
        % defaults to columnsums...
        varargout{1} = sum(X.^4 - 16*X.^2 + 5*X, 1)/2;

    end

end


text:
function [B,twom] = multicatbipartite_f(A,gamma,omega)
% MULTICATBIPARTITE_F  returns multilayer Barber modularity matrix for unordered undirected bipartite networks, function handle version
%
% Version: 2.2.0
% Date: Thu 11 Jul 2019 12:25:42 CEST
%
% MULTICATBIPARTITE_F [B,twomu] = MULTICATBIPARTITE_F(A,gamma,omega)
%
% Input: A: Cell array of MxN adjacency matrices for each layer of a
%           multilayer undirected bipartite network
%        gamma: resolution parameter
%        omega: interlayer coupling strength
%
% Output: B: function handle where B(i) returns the ith column of the
%            [(M+N)xT]x[(M+N)xT] flattened modularity tensor for the
%            multilayer bipartite network with uniform ordinal coupling (T
%            is the number of layers of the network)
%         twomu: normalisation constant
%
% Usage: [B,twomu]=multicatbipartite_f(A,gamma,omega);
%        [S,Q]=genlouvain(B); % see iterated_genlouvain.m and
%          postprocess_categorical_multilayer.m for how to improve output
%          multilayer partition
%        Q=Q/twom;
%        S=reshape(S,M+N,T);
%
%  [B,twom] = MULTICATBIPARTITE_F(A,GAMMA, OMEGA) with A a cell array of
%   matrices of equal size each representing an undirected bipartite network
%   "layer" computes the multilayer Barber modularity matrix using the quality
%   function described in Mucha et al. 2010, with intralayer resolution
%   parameter GAMMA, and with interlayer coupling OMEGA connecting all-to-all
%   layers. Once the mulilayer modularity matrix is computed,
%   optimization can be performed by the generalized Louvain code GENLOUVAIN
%   or ITERATED_GENLOUVAIN. The  output B can be used with other heuristics,
%   provided the same mapping is used to go from the multilayer tensor to
%   the multilayer flattened matrix. That is, the node-layer tuple (i,s)
%   is mapped to i + (s-1)*(M+N). [Note that we can define a mapping between
%   a multilayer partition S_m stored as an (M+N) by T matrix and the
%   corresponding flattened partition S stored as an MNT by 1 vector. In
%   particular S_m = reshape(S,M+N,T) and S = S_m(:). Note that nodes i=1:M
%   correspond to the first class (i.e. the rows of A) and nodes i=M+1:M+N
%   correspond to the second class (i.e. the columns of A) of the bipartite
%   network.]
%
%   Notes:
%     The matrices in the cell array A are assumed to be of equal size.
%     This assumption is not checked here.
%
%     For smaller systems, it is potentially more efficient (and easier) to
%     directly use the sparse quality/modularity matrix B in MULTICATBIPARTITE.
%
%     This code serves as a template and can be modified for situations
%     with other wrinkles (e.g., different intralayer null models,
%     different numbers of nodes from layer-to-layer, or systems which are
%     both multiplex and longitudinal).  That is, this code is only a
%     starting point; it is by no means exhaustive.
%
%     By using this code, the user implicitly acknowledges that the authors
%     accept no liability associated with that use.  (What are you doing
%     with it anyway that might cause there to be a potential liability?!?)
%
% References:
%       Barber, M. Modularity and community detection in bipartite networks.
%           Phys. Rev. E 76, 066102 (2007).
%
%       Mucha, P. J., Richardson, T., Macon, K., Porter, M. A. & Onnela, J.-P.
%           Community structure in time-dependent, multiscale, and multiplex networks.
%           Science 328, 876-878 (2010).

if nargin<2||isempty(gamma)
    gamma=1;
end

if nargin<3
	omega=1;
end

[m,n]=size(A{1});
N=m+n;
T=length(A);

if length(gamma)==1
    gamma=repmat(gamma,T,1);
end

k=zeros(m,T);
d=zeros(T,n);
mm=zeros(T,1);

twom=0;
for j=1:T
    twom = twom + sum(sum(A{j}));
    k(:,j)=sum(A{j},2);
    d(j,:)=sum(A{j});
    mm(j)=sum(k(:,j));
end

%interslice connections
all2all= N*[(-T+1):-1,1:(T-1)];
C=omega*spdiags(ones(N*T,2*T-2),all2all,N*T,N*T);

%bipartite modularity matrix
    function modi=modf(i)

        s=ceil(i/(N+eps));
        if mm(s)~=0
            ii=i-(s-1)*N;
            if ii<=m
                indx=(m+1:N)+(s-1)*N;
                v=A{s}(ii,:)-gamma(s)*k(ii,s)*d(s,:)/mm(s);

                modi=sparse(indx,1,v,N*T,1,n+2);
            else
                indx=(1:m)+(s-1)*N;
                v=A{s}(:,ii-m)-gamma(s)*k(:,s)*d(s,ii-m)/mm(s);

                modi=sparse(indx,1,v,N*T,1,m+2);
            end

            modi=modi+C(:,i);
        else
            modi=C(:,i);

        end
    end

B=@modf;
twom=2*twom+2*N*(T-1)*T*omega;
end

author: GenLouvain

repo: GenLouvain


text:
function [cell_metrics] = bz_CellMetricsSimple(spikes,sr)
% [cell_metrics] = bz_CellMetricsSimple(spikes,sr)
%   Calculate waveform and ACG metrics to use for unit classification
%
% USAGE
%
%   [cell_metrics] = bz_CellMetricsSimple(spikes)
%
% INPUTS
%
%   spikes = buzcode spikes structure (i.e. from bz_LoadPhy)
%   sr = sampling rate
%
% OUTPUT
%
%   cell_metrics
%       .filtWaveform
%       .rawWaveform
%           .PeaktoTrough
%           .TroughtoPeak
%           .AB_ratio
%           .trough
%       .refractory
%       .burstRoyer
%       .burstiness
%       .FR

% Antonio FR, 11/2018

if ~isnumeric(sr)
    sr = 30000;
end

%% waveform metrics
   % reshape data
   for i = 1:length(spikes.rawWaveform)
       rawWaveforms(i,:) = spikes.rawWaveform{i};
       filtWaveforms(i,:) = spikes.filtWaveform{i};
   end

   wave = [];t_before = [];t_after = [];peakA = [];peakB = [];trough = [];
   for m = 1:size(rawWaveforms,1)
        wave = interp1([1:size(filtWaveforms,2)],zscore(filtWaveforms(m,:)),[1:0.5:size(filtWaveforms,2),size(filtWaveforms,2)],'spline');
        midPoint = round((length(wave)/2));
        [MIN2,I2] = min(wave(midPoint-10:midPoint+10));
        [MAX3,I3] = max(wave(1:midPoint));
        [MAX4,I4] = max(wave((I2+midPoint-11):end));
        t_before(m) = (I2+midPoint-11)-I3;
        t_after(m) = I4;
        peakA(m) = MAX3;
        peakB(m) = MAX4;
        trough(m) = MIN2;
        clear MIN2 MAX3 MAX4 I2 I3 I4
   end
    cell_metrics.filtWaveform.PeaktoTrough = (t_before/(sr*2))';
    cell_metrics.filtWaveform.TroughtoPeak = (t_after/(sr*2))';
    cell_metrics.filtWaveform.AB_ratio = ((peakB-peakA)./(peakA+peakB))';
    cell_metrics.filtWaveform.trough = (trough)';

   wave = [];t_before = [];t_after = [];peakA = [];peakB = [];trough = [];
   for m = 1:size(rawWaveforms,1)
        wave = interp1([1:size(rawWaveforms,2)],zscore(rawWaveforms(m,:)),[1:0.5:size(rawWaveforms,2),size(rawWaveforms,2)],'spline');
        midPoint = round((length(wave)/2));
        [MIN2,I2] = min(wave(midPoint-10:midPoint+10));
        [MAX3,I3] = max(wave(1:midPoint));
        [MAX4,I4] = max(wave((I2+midPoint-11):end));
        t_before(m) = (I2+midPoint-11)-I3;
        t_after(m) = I4;
        peakA(m) = MAX3;
        peakB(m) = MAX4;
        trough(m) = MIN2;
        clear MIN2 MAX3 MAX4 I2 I3 I4
   end
    cell_metrics.rawWaveform.PeaktoTrough = (t_before/(sr*2))';
    cell_metrics.rawWaveform.TroughtoPeak = (t_after/(sr*2))';
    cell_metrics.rawWaveform.AB_ratio = ((peakB-peakA)./(peakA+peakB))';
    cell_metrics.rawWaveform.trough = (trough)';

%%  ACG metrics
for i = 1:length(spikes.times)
    [ccg,time] = CCG(spikes.times{i},ones(length(spikes.times{i}),1),'binSize',0.0005,'duration',0.100); %100ms wide CCG with 0.5ms bins
    ccg0 = find(time == 0);ccg10 = find(time == 0.01);ccg20 = find(time == 0.02);
    ccg40 = find(time == 0.04);ccg50 = find(time == 0.05);

    % Refractory period and burst index Royer 2012
    CCGmax = find(ccg(ccg0:ccg20)== max(ccg(ccg0:ccg20)));
    CCGmax = CCGmax+ccg0-1;
    if ~isempty(find(diff(ccg(ccg0:CCGmax(1)))>std(diff(ccg(ccg0:CCGmax(1)))),1))% refractory period according to Royer 2012
        cell_metrics.refractory(i,1) = find(diff(ccg(ccg0:CCGmax(1)))>std(diff(ccg(ccg0:CCGmax(1)))),1);
    else
        cell_metrics.refractory(i,1) = NaN;
    end

    CCGbase=mean(ccg(ccg40:ccg50));clear CCGmax; % mean CCG 40-50 ms
    CCGmax=find(ccg(ccg0:ccg10)==max(ccg(ccg0:ccg10)));
    CCGmax = CCGmax+ccg0-1; % only from 0 to 10ms

    if ccg(CCGmax)>CCGbase
       cell_metrics.burstRoyer(i,1)=(ccg(CCGmax(1))-CCGbase)/ccg(CCGmax(1));
    else
       cell_metrics.burstRoyer(i,1)=(ccg(CCGmax(1))-CCGbase)/(CCGbase);
    end

    % Burstiness
     binssum = 6; %for isis below 6ms
     binSize = 1e-3; nBins = 100; binEdges = (0:nBins) * binSize;
     t = ((0:nBins-1) + 0.5)' * binSize;
     ISIs = diff(spikes.times{i});
     hISI = histc(ISIs, binEdges);
     cell_metrics.burstiness(i,1) = sum(hISI(1:binssum))/sum(hISI);

    % Sam's ACG fit
    rsquare = [];   offset = 106;
    g = fittype('c*exp(-x/a)-d*exp(-x/b)');
    for j = 1:size(ccg,2)
        x = [1:95]';
        y = ccg(x+offset,j)/max(ccg(x+offset,j));

        [f0,gof] = fit(x/2,y,g,'StartPoint',[25, 1, 1, 1.5],'Lower',[1,0.1,-Inf, -Inf],'Upper',[100, 10,2, 10]);
    %     f0 = fit(x/2,y,'exp2','StartPoint',[1, -0.015, -1, -1]);
        fit_params(:,j) = coeffvalues(f0);
        xx = linspace(0,48,100);
        rsquare(j) = gof.rsquare;
        %figure, plot(x/2,y,'-o',x/2,f0(x/2),'r-'); title(['Unit ', num2str(j), ' r^2: ' num2str(rsquare(j))])
    end

    ACG.tau_decay = fit_params(1,:);
    ACG.tau_rise = fit_params(2,:);
    ACG.c = fit_params(3,:);
    ACG.d = fit_params(4,:);
    ACG.fit_rsquare = rsquare;

    clear s ccg CCGmax refT CCGbase ISIs hISI ;
end

%% FR
for i = 1:length(spikes.times)
    cell_metrics.FR(i,1) = length(spikes.times{i})/(spikes.times{i}(end)-spikes.times{i}(1));
end

%% Preliminary pyr/int classification
    % avg FR > 6Hz or narrow waveform (trought to peak < 5 ms) = INT
for i = 1:length(spikes.times)
    if cell_metrics.FR(i,1) > 6 || cell_metrics.filtWaveform.TroughtoPeak(i) < 5e-4
        cell_metrics.putativeClass(i,1) = 2; % int
    else
        cell_metrics.putativeClass(i,1) = 1; % pyr
    end
end

end

%%

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [coef, C] = SS_to_AR(F, Q, k, diagonal)
%
% Extract the parameters of a vector autoregresssive process of order k from the state-space form.
% [coef, C] = SS_to_AR(F, Q, k, diagonal)

if nargin<4, diagonal = 0; end

s = length(Q) / k;
bs = s*ones(1,k);
coef = zeros(s,s,k);
for i=1:k
  if diagonal
    coef(:,:,i) = diag(diag(F(block(1,bs), block(i,bs))));
  else
    coef(:,:,i) = F(block(1,bs), block(i,bs));
  end
end
C = Q(block(1,bs), block(1,bs));
if diagonal
  C = diag(diag(C));
end
%C = sqrt(Q(block(1,bs), block(1,bs))); % since cov(1,1) of full vector = C C'

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function create_cache(cache_name)
  % CREATE_CACHE Save the workspace of the caller function using the provided
  % cache_name, possibly overwriting it.
  %
  % create_cache(cache_name)
  %
  % Inputs:
  %   cache_name  path to .mat cache file
  %
  % See also: find_cache

  if isunix
    % get a list of current variables in this scope, this is the input "state"
    variables = evalin('caller','who');
    variable_names = sprintf('^%s$|',variables{:});
    save_cmd = sprintf('save(''%s'',''-regexp'',''%s'');',cache_name,variable_names);
    evalin('caller',save_cmd);
  end
end

author: alecjacobson

repo: gptoolbox


text:
function physio = tapas_physio_job2physio(job)
% Converts job from SPM batch editor to physio-structure
%
%   physio = tapas_physio_job2physio(job)
%
% IN
%
% OUT
%   physio  physio input structure, as use by
%           tapas_physio_main_create_regressors
%
% EXAMPLE
%   physio = tapas_physio_job2physio(job)
%
%   See also tapas_physio_cfg_matlabbatch tapas_physio_main_create_regressors

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

physio                      = tapas_physio_new();

%% Use existing properties that are cfg_choices in job to overwrite
% properties of physio and set corresponding method

physio = tapas_physio_update_from_job(physio, job, ...
    {'preproc.cardiac.posthoc_cpulse_select', ...
    'preproc.cardiac.initial_cpulse_select', 'scan_timing.sync'}, ...
    {'preproc.cardiac.posthoc_cpulse_select', ...
    'preproc.cardiac.initial_cpulse_select', 'scan_timing.sync'}, ...
    true, ...
    'method');

%% Take over model substructs as is
modelArray =  ...
    {'movement', 'retroicor', 'rvt', 'hrv', ...
    'noise_rois', 'other'};

physio = tapas_physio_update_from_job(physio, job, ...
    strcat('model.', modelArray), strcat('model.', modelArray), ...
    true, 'include');

%% Convert yes => true (=1) and no => false (=0)
nModels = numel(modelArray);
for iModel = 1:nModels
    physio.model.(modelArray{iModel}).include = strcmpi(...
        physio.model.(modelArray{iModel}).include, 'yes');
end

%% Take over yes/no substructs as is, yes/no will become 'include' property
yesNoArray =  ...
    {'preproc.cardiac.filter'};

physio = tapas_physio_update_from_job(physio, job, ...
    yesNoArray, yesNoArray,  true, 'include');

%% Convert yes => true (=1) and no => false (=0)
nChoices = numel(yesNoArray);
for iChoice = 1:nChoices
    try
        eval(sprintf(['physio.%s.include = strcmpi(' ...
            'physio.%s.include, ''yes'');'], yesNoArray{iChoice}, ...
            yesNoArray{iChoice}));
    catch err
        tapas_physio_log(sprintf('No property %s defined in job (error: %s)', ...
            yesNoArray{iChoice}, err.message), [], 1);
    end
end

%% Use existing properties in job to overwrite properties of physio
physio = tapas_physio_update_from_job(physio, job, ...
    {'preproc.cardiac.modality', 'preproc.respiratory', 'scan_timing.sqpar', ...
    'log_files', 'verbose', 'save_dir',...
    'model.orthogonalise', 'model.censor_unreliable_recording_intervals', ...
    'model.output_multiple_regressors', ...
    'model.output_physio'}, ...
    {'preproc.cardiac.modality', 'preproc.respiratory', 'scan_timing.sqpar', ...
    'log_files', 'verbose', 'save_dir',...
    'model.orthogonalise', 'model.censor_unreliable_recording_intervals', ...
    'model.output_multiple_regressors', ...
    'model.output_physio'}, ...
    false);

author: translationalneuromodeling

repo: tapas


text:
function datain(imagefig,varargins)
global M
global Nsamples
set(gcf,'WindowButtonMotionFcn',[]);
set(gcf,'windowbuttondownfcn',{@track});

%%% This function performes the FFT and spectrum analysis of the
%%% Input data. It is not related dirctely to the filtering process,
%%% But performs post-processing.

ini=Nsamples;
temp=get(gcf,'userdata');
temp=temp(ini:size(temp,1),:);
len=length(temp(:,3));                  % Determines the length of x
fs=1/(temp(end,3)-temp(ini,3))*len                % mean sample rate

Pad = 2 * fs * len;                  % Setting the length of ZeroPadded x
Pad = floor(Pad);
w=0:(Pad-1);
half_w=floor(length(w)/2);
win=(kaiser(len,5))';   % Creates a Kaiser window
for k=1:2
    temp(:,k)=temp(:,k).*win';    % Passes x through the Kaiser window
end
FFX=(abs(fft(temp(:,1) ,Pad)));
FFY=(abs(fft(temp(:,2) ,Pad)));

figure
subplot(2,1,1), plot(fs / Pad * w(1:half_w), 20*log10(FFX(1:half_w)));
grid on, ylabel(['[dB]']), title('X axis motion spectrum');

subplot(2,1,2), plot(fs / Pad * w(1:half_w), 20*log10(FFY(1:half_w)));
grid on, ylabel(['[dB]']), xlabel(['[Hz]']), title('Y axis motion spectrum');


text:
function hdr = empty_hdr
% Create an empty NIFTI-1 header
% FORMAT hdr = empty_hdr
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

org = niftistruc;
for i=1:length(org),
    hdr.(org(i).label) = org(i).def;
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [kJ] = kWh2kJ(kWh)
% Convert energy or work from kilowatt-hours to kilojoules.
% Chad A. Greene 2012
kJ = kWh*3600;


text:
classdef EuclideanLoss < dagnn.Loss

  methods
    function outputs = forward(obj, inputs, params)
      outputs{1} = vl_nneuclideanloss(inputs{1}, inputs{2}, ...
                                      'instanceWeights', inputs{3}) ;
      n = obj.numAveraged ;
      m = n + size(inputs{1}, 4) ;
      obj.average = (n * obj.average + gather(outputs{1})) / m ;
      obj.numAveraged = m ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      derInputs{1} = vl_nneuclideanloss(inputs{1}, inputs{2}, derOutputs{1}, ...
                                        'instanceWeights', inputs{3}) ;
      derInputs{2} = [] ;
      derInputs{3} = [] ;
      derParams = {} ;
    end

    function obj = EuclideanLoss(varargin)
      obj.load(varargin) ;
    end
  end
end

repo: MFT


text:
function isargxs(varargin)
%ISARGXS throws an error if not a [1x3] and/or [1x6] vector
%
%   Usage: isargxs(args)
%
%   Input parameters:
%       args        - list of args
%
%   See also: isargscalar, isargvector, isargposition

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

%% ===== Checking for vector =============================================
for ii = 1:nargin
    if ~isvector(varargin{ii}) || ...
            ( ~all(size(varargin{ii})==[1 3]) && ...
        ~all(size(varargin{ii})==[1 6]) )
        error('%s needs to be a [1x3] or [1x6] vector.',inputname(ii));
    end
end

author: sfstoolbox

repo: sfs-matlab


text:
function [CrowdDis]=Crowding(Pop)
%Harmonic average distance of each solution in the decision space
%Return: the crowding distance of each individual

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Wenhua Li

    [N, ~]=size(Pop);
    K=N-1;
    Z = min(Pop,[],1);
    Zmax = max(Pop,[],1);
    pop=(Pop-repmat(Z,N,1))./repmat(Zmax-Z,N,1);
    distance=pdist2(pop,pop);
    [value,~]=sort(distance,2);
    CrowdDis=K./sum(1./value(:,2:N),2);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [spikeTimes, clusterIDs, amplitudes, templates, templateFeatures, ...
    templateFeatureInds, pcFeatures, pcFeatureInds] = rezToPhy(rez, savePath)
% pull out results from kilosort's rez to either return to workspace or to
% save in the appropriate format for the phy GUI to run on. If you provide
% a savePath it should be a folder, and you will need to have npy-matlab
% available (https://github.com/kwikteam/npy-matlab)
%
% spikeTimes will be in samples, not seconds

outputs = {'amplitudes.npy', 'channel_map.npy', 'channel_positions.npy', 'pc_features.npy', ...
           'pc_feature_ind.npy', 'similar_templates.npy', 'spike_clusters.npy', 'spike_templates.npy', ...
           'spike_times.npy', 'templates.npy', 'templates_ind.npy', 'template_features.npy', ...
           'template_feature_ind.npy', 'whitening_mat.npy', 'whitening_mat_inv.npy'};

fs = dir(fullfile(savePath, '*.npy'));
for i = 1:length(fs)
    fname = fs(i).name;
    % don't delete .npy files which have nothing to do with us
    if find(strcmp(fname, outputs))
        delete(fullfile(savePath, fname));
    end
end
if exist(fullfile(savePath, '.phy'), 'dir')
    rmdir(fullfile(savePath, '.phy'), 's');
end

spikeTimes = uint64(rez.st3(:,1));
% [spikeTimes, ii] = sort(spikeTimes);
spikeTemplates = uint32(rez.st3(:,2));
if size(rez.st3,2)>4
    spikeClusters = uint32(1+rez.st3(:,5));
end
amplitudes = rez.st3(:,3);

Nchan = rez.ops.Nchan;

% try
%     load(rez.ops.chanMap);
% catch
%    chanMap0ind  = [0:Nchan-1]';
%    connected    = ones(Nchan, 1);
%    xcoords      = ones(Nchan, 1);
%    ycoords      = (1:Nchan)';
% end
% chanMap0 = chanMap(connected>1e-6);

connected   = rez.connected(:);
xcoords     = rez.xcoords(:);
ycoords     = rez.ycoords(:);
chanMap     = rez.ops.chanMap(:);
chanMap0ind = chanMap - 1;

nt0 = size(rez.W,1);
U = rez.U;
W = rez.W;

% for i = 1:length(chanMap0)
%     chanMap0(i) = chanMap0(i) - sum(chanMap0(i) > chanMap(connected<1e-6));
% end
% [~, invchanMap0] = sort(chanMap0);

templates = zeros(Nchan, nt0, rez.ops.Nfilt, 'single');
for iNN = 1:rez.ops.Nfilt
   templates(:,:,iNN) = squeeze(U(:,iNN,:)) * squeeze(W(:,iNN,:))';
end
templates = permute(templates, [3 2 1]); % now it's nTemplates x nSamples x nChannels
templatesInds = repmat([0:size(templates,3)-1], size(templates,1), 1); % we include all channels so this is trivial

templateFeatures = rez.cProj;
templateFeatureInds = uint32(rez.iNeigh);
pcFeatures = rez.cProjPC;
pcFeatureInds = uint32(rez.iNeighPC);

if ~isempty(savePath)

    writeNPY(spikeTimes, fullfile(savePath, 'spike_times.npy'));
    writeNPY(uint32(spikeTemplates-1), fullfile(savePath, 'spike_templates.npy')); % -1 for zero indexing
    if size(rez.st3,2)>4
        writeNPY(int32(spikeClusters-1), fullfile(savePath, 'spike_clusters.npy')); % -1 for zero indexing
    else
        writeNPY(int32(spikeTemplates-1), fullfile(savePath, 'spike_clusters.npy')); % -1 for zero indexing
    end
    writeNPY(amplitudes, fullfile(savePath, 'amplitudes.npy'));
    writeNPY(templates, fullfile(savePath, 'templates.npy'));
    writeNPY(templatesInds, fullfile(savePath, 'templates_ind.npy'));

%     Fs = rez.ops.fs;
    conn        = logical(connected);
    chanMap0ind = int32(chanMap0ind);

    writeNPY(chanMap0ind(conn), fullfile(savePath, 'channel_map.npy'));
    %writeNPY(connected, fullfile(savePath, 'connected.npy'));
%     writeNPY(Fs, fullfile(savePath, 'Fs.npy'));
    writeNPY([xcoords(conn) ycoords(conn)], fullfile(savePath, 'channel_positions.npy'));

    writeNPY(templateFeatures, fullfile(savePath, 'template_features.npy'));
    writeNPY(templateFeatureInds'-1, fullfile(savePath, 'template_feature_ind.npy'));% -1 for zero indexing
    writeNPY(pcFeatures, fullfile(savePath, 'pc_features.npy'));
    writeNPY(pcFeatureInds'-1, fullfile(savePath, 'pc_feature_ind.npy'));% -1 for zero indexing

    whiteningMatrix = rez.Wrot/200;
    whiteningMatrixInv = whiteningMatrix^-1;
    writeNPY(whiteningMatrix, fullfile(savePath, 'whitening_mat.npy'));
    writeNPY(whiteningMatrixInv, fullfile(savePath, 'whitening_mat_inv.npy'));

    if isfield(rez, 'simScore')
        similarTemplates = rez.simScore;
        writeNPY(similarTemplates, fullfile(savePath, 'similar_templates.npy'));
    end

     %make params file
    if ~exist(fullfile(savePath,'params.py'),'file')
        fid = fopen(fullfile(savePath,'params.py'), 'w');

        [~, fname, ext] = fileparts(rez.ops.fbinary);

        fprintf(fid,['dat_path = ''',fname ext '''\n']);
        fprintf(fid,'n_channels_dat = %i\n',rez.ops.NchanTOT);
        fprintf(fid,'dtype = ''int16''\n');
        fprintf(fid,'offset = 0\n');
        if mod(rez.ops.fs,1)
            fprintf(fid,'sample_rate = %i\n',rez.ops.fs);
        else
            fprintf(fid,'sample_rate = %i.\n',rez.ops.fs);
        end
        fprintf(fid,'hp_filtered = False');
        fclose(fid);
    end
end

author: cortex-lab

repo: KiloSort

save_path: github-repos/MATLAB/cortex-lab-KiloSort


text:
function options = rbfcreate(x, y, varargin)
%RBFCREATE Creates an RBF interpolation
%   OPTIONS = RBFSET(X, Y, 'NAME1',VALUE1,'NAME2',VALUE2,...) creates an
%   radial base function interpolation
%
%   RBFCREATE with no input arguments displays all property names and their
%   possible values.
%
%RBFCREATE PROPERTIES
%

%
% Alex Chirokov, alex.chirokov@gmail.com
% 16 Feb 2006
tic;
% Print out possible values of properties.
if (nargin == 0) & (nargout == 0)
  fprintf('               x: [ dim by n matrix of coordinates for the nodes ]\n');
  fprintf('               y: [   1 by n vector of values at nodes ]\n');
  fprintf('     RBFFunction: [ gaussian  | thinplate | cubic | multiquadrics | {linear} ]\n');
  fprintf('     RBFConstant: [ positive scalar     ]\n');
  fprintf('       RBFSmooth: [ positive scalar {0} ]\n');
  fprintf('           Stats: [ on | {off} ]\n');
  fprintf('\n');
  return;
end
Names = [
    'RBFFunction      '
    'RBFConstant      '
    'RBFSmooth        '
    'Stats            '
];
[m,n] = size(Names);
names = lower(Names);

options = [];
for j = 1:m
  options.(deblank(Names(j,:))) = [];
end

%**************************************************************************
%Check input arrays
%**************************************************************************
[nXDim nXCount]=size(x);
[nYDim nYCount]=size(y);

if (nXCount~=nYCount)
  error(sprintf('x and y should have the same number of rows'));
end;

if (nYDim~=1)
  error(sprintf('y should be n by 1 vector'));
end;

options.('x')           = x;
options.('y')           = y;
%**************************************************************************
%Default values
%**************************************************************************
options.('RBFFunction') = 'linear';
options.('RBFConstant') = (prod(max(x')-min(x'))/nXCount)^(1/nXDim); %approx. average distance between the nodes
options.('RBFSmooth')   = 0;
options.('Stats')       = 'off';

%**************************************************************************
% Argument parsing code: similar to ODESET.m
%**************************************************************************

i = 1;
% A finite state machine to parse name-value pairs.
if rem(nargin-2,2) ~= 0
  error('Arguments must occur in name-value pairs.');
end
expectval = 0;                          % start expecting a name, not a value
while i <= nargin-2
  arg = varargin{i};

  if ~expectval
    if ~isstr(arg)
      error(sprintf('Expected argument %d to be a string property name.', i));
    end

    lowArg = lower(arg);
    j = strmatch(lowArg,names);
    if isempty(j)                       % if no matches
      error(sprintf('Unrecognized property name ''%s''.', arg));
    elseif length(j) > 1                % if more than one match
      % Check for any exact matches (in case any names are subsets of others)
      k = strmatch(lowArg,names,'exact');
      if length(k) == 1
        j = k;
      else
        msg = sprintf('Ambiguous property name ''%s'' ', arg);
        msg = [msg '(' deblank(Names(j(1),:))];
        for k = j(2:length(j))'
          msg = [msg ', ' deblank(Names(k,:))];
        end
        msg = sprintf('%s).', msg);
        error(msg);
      end
    end
    expectval = 1;                      % we expect a value next

  else
    options.(deblank(Names(j,:))) = arg;
    expectval = 0;
  end
  i = i + 1;
end

if expectval
  error(sprintf('Expected value for property ''%s''.', arg));
end

%**************************************************************************
% Creating RBF Interpolatin
%**************************************************************************

switch lower(options.('RBFFunction'))
      case 'linear'
        options.('rbfphi')   = @rbfphi_linear;
      case 'cubic'
        options.('rbfphi')   = @rbfphi_cubic;
      case 'multiquadric'
        options.('rbfphi')   = @rbfphi_multiquadrics;
      case 'thinplate'
        options.('rbfphi')   = @rbfphi_thinplate;
      case 'gaussian'
        options.('rbfphi')   = @rbfphi_gaussian;
    otherwise
        options.('rbfphi')   = @rbfphi_linear;
end

phi       = options.('rbfphi');

A=rbfAssemble(x, phi, options.('RBFConstant'), options.('RBFSmooth'));

b=[y'; zeros(nXDim+1, 1)];

%inverse
rbfcoeff=A\b;

%SVD
% [U,S,V] = svd(A);
%
% for i=1:1:nXCount+1
%     if (S(i,i)>0) S(i,i)=1/S(i,i); end;
% end;
% rbfcoeff = V*S'*U*b;

options.('rbfcoeff') = rbfcoeff;

if (strcmp(options.('Stats'),'on'))
    fprintf('%d point RBF interpolation was created in %e sec\n', length(y), toc);
    fprintf('\n');
end;

function [A]=rbfAssemble(x, phi, const, smooth)
[dim n]=size(x);
A=zeros(n,n);
for i=1:n
    for j=1:i
        r=norm(x(:,i)-x(:,j));
        temp=feval(phi,r, const);
        A(i,j)=temp;
        A(j,i)=temp;
    end
    A(i,i) = A(i,i) - smooth;
end
% Polynomial part
P=[ones(n,1) x'];
A = [ A      P
      P' zeros(dim+1,dim+1)];

%**************************************************************************
% Radial Base Functions
%**************************************************************************
function u=rbfphi_linear(r, const)
u=r;

function u=rbfphi_cubic(r, const)
u=r.*r.*r;

function u=rbfphi_gaussian(r, const)
u=exp(-0.5*r.*r/(const*const));

function u=rbfphi_multiquadrics(r, const)
u=sqrt(1+r.*r/(const*const));

function u=rbfphi_thinplate(r, const)
u=r.*r.*log(r+1);

repo: REST


text:
function inspect_issue2138

% WALLTIME 00:10:00
% MEM 2gb

%%

bidsroot = dccnpath('/home/common/matlab/fieldtrip/data/test/issue2138');

% this is the original
filename1 = fullfile(bidsroot, 'sub-001/ses-01/eeg/sub-001_ses-01_task-WorkingMemory_run-1_eeg.set');

% this is the one with enriched HED tags
filename2 = fullfile(bidsroot, 'sub-002/ses-01/eeg/sub-002_ses-01_task-WorkingMemory_run-1_eeg.set');

%%
% just to check that this works

hdr1 = ft_read_header(filename1, 'readbids', false);
event1 = ft_read_event(filename1, 'readbids', false);

hdr2 = ft_read_header(filename2, 'readbids', false);
event2 = ft_read_event(filename2, 'readbids', false);

%%
% default is readbids=true, which causes the sidecars to be read

hdr1 = ft_read_header(filename1);
event1 = ft_read_event(filename1);

hdr2 = ft_read_header(filename2);
event2 = ft_read_event(filename2);

%%
% show the events

cfg = [];
cfg.dataset = filename1;
cfg.trialdef.eventtype = '?';
ft_definetrial(cfg);

cfg = [];
cfg.dataset = filename2;
cfg.trialdef.eventtype = '?';
ft_definetrial(cfg);

%%
% this is some random ERP analysis using "normal" events

cfg = [];
cfg.dataset = filename1;
cfg.trialdef.eventtype = 'show_letter';
cfg.trialdef.eventvalue = {'A' 'B' 'C' 'D' 'E' 'F' 'G' 'H' 'I' 'J' 'K' 'L' 'M' 'N' 'O' 'P' 'Q' 'R' 'S' 'T' 'U' 'V' 'W' 'X' 'Y' 'Z'};
cfg.trialdef.prestim = 0.3;
cfg.trialdef.poststim = 0.7;
cfg = ft_definetrial(cfg);

cfg.baselinewindow = [-inf 0];
cfg.demean = 'yes';
rawdata1 = ft_preprocessing(cfg);

% make a selective average
cfg = [];
cfg.trials = ismember(rawdata1.trialinfo.eventvalue, {'A' 'B' 'C'});
erp1 = ft_timelockanalysis(cfg, rawdata1);

assert(sum(erp1.cfg.trials)==9);

%%
% this is exactly the same analysis for a dataset that also contains HED tags
% but it does not use those HED tags

cfg = [];
cfg.dataset = filename2;
cfg.trialdef.eventtype = 'show_letter';
cfg.trialdef.eventvalue = {'A' 'B' 'C' 'D' 'E' 'F' 'G' 'H' 'I' 'J' 'K' 'L' 'M' 'N' 'O' 'P' 'Q' 'R' 'S' 'T' 'U' 'V' 'W' 'X' 'Y' 'Z'};
cfg.trialdef.prestim = 0.3;
cfg.trialdef.poststim = 0.7;
cfg = ft_definetrial(cfg);

cfg.baselinewindow = [-inf 0];
cfg.demean = 'yes';
rawdata2 = ft_preprocessing(cfg);

% make a selective average
cfg = [];
cfg.trials = ismember(rawdata2.trialinfo.eventvalue, {'A' 'B' 'C'});
erp2 = ft_timelockanalysis(cfg, rawdata2);

assert(sum(erp2.cfg.trials)==9);

%%
% this is the same analysis, now using the HED tags

cfg = [];
cfg.dataset = filename2;
cfg.trialdef.eventtype = 'HED';
cfg.trialdef.eventvalue = 'Sensory-event,Visual-presentation,(Character,(Center-of,Computer-screen)),Target,*';
cfg.trialdef.prestim = 0.3;
cfg.trialdef.poststim = 0.7;
cfg = ft_definetrial(cfg);

cfg.baselinewindow = [-inf 0];
cfg.demean = 'yes';
rawdata3 = ft_preprocessing(cfg);

% make a selective average
cfg = [];
cfg.trials = ~cellfun(@isempty, regexp(rawdata3.trialinfo.eventvalue, '(Character,Label/[ABC])'));
erp3 = ft_timelockanalysis(cfg, rawdata3);

assert(sum(erp3.cfg.trials)==9);

%%
% now using a trialfun that is specialized in HED tags

cfg = [];
cfg.dataset = filename2; % this has complete tags
cfg.trialfun = 'ft_trialfun_hed';
cfg.trialdef.regexp = '.*,Target,\(Character,Label/[ABC]\),.*';
cfg.trialdef.prestim = 0.3;
cfg.trialdef.poststim = 0.7;
cfg = ft_definetrial(cfg);

cfg.baselinewindow = [-inf 0];
cfg.demean = 'yes';
rawdata4 = ft_preprocessing(cfg);

% make an average, it would be possible to make a selection using rawdata4.trialinfo
cfg = [];
erp4 = ft_timelockanalysis(cfg, rawdata4);

%%
% assemble the HED tags when not present in the file

cfg = [];
cfg.dataset = filename1; % HED tags need to be assembled
cfg.trialfun = 'ft_trialfun_hed';
cfg.trialdef.regexp = '.*,Target,\(Character,Label/[ABC]\),.*';
cfg.trialdef.prestim = 0.3;
cfg.trialdef.poststim = 0.7;
cfg = ft_definetrial(cfg);

cfg.baselinewindow = [-inf 0];
cfg.demean = 'yes';
rawdata5 = ft_preprocessing(cfg);

% make an average, it would be possible to make a selection using rawdata4.trialinfo
cfg = [];
erp5 = ft_timelockanalysis(cfg, rawdata5);

%%
% just to demonstrate that reading+filtering and segmenting can also be done the other way around

cfg = [];
cfg.dataset = filename2;
cfg.hpfreq = 0.2;
cfg.hpfilter = 'yes';
rawdata6_continuous = ft_preprocessing(cfg);

% determine the selection of epochs
cfg = [];
cfg.dataset = filename2;
cfg.trialfun = 'ft_trialfun_hed';
cfg.trialdef.regexp = '.*,Target,\(Character,Label/[ABC]\),.*';
cfg.trialdef.prestim = 0.3;
cfg.trialdef.poststim = 0.7;
cfg = ft_definetrial(cfg);

% make the selection of epochs
rawdata6_epoched = ft_redefinetrial(cfg, rawdata6_continuous);

% make an average
cfg = [];
erp5 = ft_timelockanalysis(cfg, rawdata5_epoched);

%%
% all ERPs should be the same

assert(isequal(erp1.avg, erp2.avg));
assert(isequal(erp1.avg, erp3.avg));
assert(isequal(erp1.avg, erp4.avg));
assert(isequal(erp1.avg, erp5.avg));
% assert(isequal(erp1.avg, erp6.avg)); % this has been preprocessed differently


text:
function F = lmi(X,handlestring,dummy,noprune,symmetryKnown)

superiorto('sdpvar')

persistent temp
persistent F0

fast = 0;
noprune = 0;
symmetryKnown = 0;
switch nargin
    case 0
        F.clauses={};
        F.LMIid    = [];
        F.savedata = [];
        F.dualized = 0;
        F = class(F,'lmi');
        F0 = F;
        return
    case 1
        handlestring = '';
        % ok probably, check later
    case 2
        if ~isa(handlestring,'char')
            error('The handle must be a string')
        end
    case 4
        handlestring = '';
        noprune = 1;
    case 5
        handlestring = '';
        symmetryKnown = 1;

    otherwise
        error('Wrong number of arguments')
end

if isempty(F0)
    F.clauses  = {};
    F.LMIid    = [];
    F.savedata = [];
    F.dualized = 0;
    try
        F = class(F,'lmi');
    catch
        disp('Failure in creating SET object. Restart MATLAB and try again, or type clear classes')
    end
    F0 = F;
else
    F = F0;
end

% Return empty LMI in degenerate case
if isempty(X)
    % A bit inconsitently, we create a constraint, which is empty. We count
    % the number of ocnstraints by looking at LMIid, which still is zero
    F.clauses{1}.data=[];
    F.clauses{1}.type = [];
    F.clauses{1}.symbolic=[];
    F.clauses{1}.handle=[];
    F.clauses{1}.strict = [];
    F.clauses{1}.cut = [];
    F.clauses{1}.expanded = [];
    F.clauses{1}.lift = [];
    F.clauses{1}.schurfun  = [];
    F.clauses{1}.schurdata = [];
    F.clauses{1}.jointprobabilistic = [];
    F.clauses{1}.confidencelevel = [];
    F.clauses{1}.extra = [];
    return
end

% 0 : Nothing (i.e error)
% 1 : LMI
% 2 : Elementwise
% 3 : Equality
% 4 : SOCC
% 5 : Rotated Lorentz
% 9 : KYP constraint
%10 : Eigenvalue constraint
%11 : Sum of square constraint
%12 : Logic CNF
%15 : Deterministic uncertain variable
%16 : Random uncertain variable
%20 : Power cone
%21 : EXP
%22 : Stacked EXPCONE
%30 : User generated Schur
%40 : Generalized KYP
%50 : SOS2
%51 : SOS1
%52 : semivar
%53 : semiintvar
%54 : Stacked SOCP
%55 : Complementary
%56 : Meta constraint to be expanded (implies, iff)
%58 : Stacked power cone
%60 : Chance constraint

switch class(X)
    case 'lmi'
        F = X;
        return
    case 'char'
        try
            % Parse (check for old notation, take care of ||...|| notation)
            sdpvarExpr = parseLMI(X);
            % Old notation
            X = strrep(X,'.>','>');
            X = strrep(X,'.<','<');
        catch
            error(lasterr)
        end
        % Evaluate
        try
            Fi=evalin('caller',sdpvarExpr);
            switch class(Fi)
                case 'sdpvar'
                    Fi = {Fi};
                    strict = 0;
                case 'constraint'
                    [Fi,strict] = getlist(Fi);
                case 'lmi'
                    F = Fi;
                    return;
                otherwise
                    error('The string does not define an SDPVAR object')
            end

        catch
            error(lasterr)
        end

    case 'sdpvar'
        Fi = {X};
        strict = 0;
        X='Numeric value';

    case 'logic'
        Fi = {X};
        strict = 0;
        X = 'Numeric value';

    case {'struct','constraint'}
        [Fi,strict,LMIIdentifiers,tags] = getlist(X);
        if isempty(handlestring) || length(handlestring)==0
            handlestring = tags{1};
        end
        X='Numeric value';

    otherwise
        error('The first argument should be an sdpvar object or string.')
end

% Return empty LMI in degenerate case
if all(cellfun('isempty',Fi))
    return
end

TypeofConstraint = zeros(length(Fi),1)-1;

% support for SOS constraints placed in vector
if length(Fi) == 1
    if isa(Fi{1},'sdpvar')% is(Fi{1},'sos')
        if gethackflag(Fi{1})==1
            % Expand to a set of SOS constraints
            if ~issymmetric(Fi{1})
                p = Fi{1}(:);
                for i = 1:length(p)
                    Fi{i} = p(i);
                end
                TypeofConstraint = zeros(length(Fi),1)-1;
                strict = zeros(length(Fi),1)-1;
            end
        end
    end
end

i = 1;
while i <= length(Fi)
    thisFi = Fi{i};
    if ~isempty(thisFi)
        if isa(thisFi,'logic')
            TypeofConstraint(i) = 12;
        else
            TypeofConstraint(i) = gethackflag(Fi{i});
        end
        % User just sent a sdpvar object, interpret as X>0
        if (TypeofConstraint(i)==0)
            TypeofConstraint(i) = 1;
        end

        if (TypeofConstraint(i)==1) && (symmetryKnown == 0)
            [n,m]=size(thisFi);
            if (n~=m) || ((TypeofConstraint(i) == 1) && (n*m==1)) || ~ishermitian(thisFi)
                TypeofConstraint(i) = 2;
            end
        end

        if TypeofConstraint(i) == 2

            % remove constraint of the type set(0 >= 0)
            B = getbase(thisFi(:));
            if ~noprune
                Bv = B(:,2:end);
                notused = find((~any(Bv,2)) & (B(:,1)>=0));
                if ~isempty(notused)
                    used = setdiff(1:size(Bv,1),notused);
                    thisFi = thisFi(used);
                end
            end
        end

        % save cleaned version
        Fi{i} = thisFi;

        switch TypeofConstraint(i)
            case {1,2,3,4,5,7,8,9,10,11,12,13,15,16,20,21,22,30,40,50,51,52,53,54,57,58}
                i = i + 1;
            otherwise
                error('Error in argument in LMI. Please report bug');
        end
    end
end

% Special case for x<y<z<...
% Merge into one constraint, gives nicer displays
if ~exist('LMIIdentifiers','var')
   LMIIdentifiers = yalmip('lmiid');
end

if all(TypeofConstraint == 2) && all(strict==strict(1))
    if length(Fi)>1
        vecF = [];
        sizes = zeros(length(Fi),1);
        % Speed up concatenation of merging
        for i = 1:length(Fi)
            sizes(i) = prod(size(Fi{i}));
        end
        if all(sizes == 1)
            vecF = [Fi{:}]';
        else
            for i = 1:length(Fi)
                fi = Fi{i};
                if sizes(i) > 1
                    fi = reshape(fi,prod(size(fi)),1);
                end
                vecF = [vecF;fi];
            end
        end
    else
        vecF = Fi{1};
        if prod(size(vecF))>1
            vecF = vecF(:);
        end
    end
    if isempty(temp)
        temp.data=vecF;
        temp.type = 2;
        temp.symbolic=X;
        temp.handle=handlestring;
        temp.strict = strict(1);
        temp.cut = 0;
        temp.expanded = 0;
        temp.lift = 0;
        temp.schurfun  = '';
        temp.schurdata = [];
        temp.jointprobabilistic = [];
        temp.confidencelevel = [];
        temp.extra = [];
    else
        temp.data=vecF;
        temp.symbolic=X;
        temp.handle=handlestring;
        temp.strict = strict(1);
    end
    F.clauses{1}{1} = temp;

%     F.clauses{1}{1}.data=vecF;
%     F.clauses{1}{1}.type = 2;
%     F.clauses{1}{1}.symbolic=X;
%     F.clauses{1}{1}.handle=handlestring;
%     F.clauses{1}{1}.strict = strict(1);
%     F.clauses{1}{1}.cut = 0;
%     F.clauses{1}{1}.expanded = 0;
%     F.clauses{1}{1}.lift = 0;
%     F.clauses{1}{1}.schurfun  = '';
%     F.clauses{1}{1}.schurdata = [];
%     F.clauses{1}{1}.jointprobabilistic = [];
%     F.clauses{1}{1}.confidencelevel = [];
%     F.clauses{1}{1}.extra = [];
    F.LMIid = [F.LMIid LMIIdentifiers(1)];
else
    for i = 1:length(Fi)
        switch TypeofConstraint(i)
            case {1,2,3,4,5,7,8,9,10,11,12,13,15,16,20,21,22,30,40,50,51,52,53,54,57,58}
                F.clauses{1}{i}.data=Fi{i};
                F.clauses{1}{i}.type = TypeofConstraint(i);
                F.clauses{1}{i}.symbolic=X;
                F.clauses{1}{i}.handle=handlestring;
                F.clauses{1}{i}.strict = strict(i);
                F.clauses{1}{i}.cut = 0;
                F.clauses{1}{i}.expanded = 0;
                F.clauses{1}{i}.lift = 0;
                F.clauses{1}{i}.schurfun  = '';
                F.clauses{1}{i}.schurdata = [];
                F.clauses{1}{i}.jointprobabilistic = [];
                F.clauses{1}{i}.confidencelevel = [];
                F.clauses{1}{i}.extra = [];
                F.LMIid = [F.LMIid LMIIdentifiers(i)];
                i = i + 1;
            otherwise
                error('Error in argument in LMI. Please report bug');
        end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [diff_dose] = dicomrt_dosediff(doseone,dosetwo,method,dose_xmesh,dose_ymesh,dose_zmesh,VOI,voi2use)
% dicomrt_dosediff(doseone,dosetwo,method,dose_xmesh,dose_ymesh,dose_zmesh,VOI,voi2use)
%
% Calculate dose difference between two 3D matrices
%
% doseone and dosetwo can be rtplan and/or monte carlo 3D dose distributions.
% method is an OPTIONAL rameter which specify the way doses are normalised
%
% 1. method=a          matrices are normalised to the value a expressed in Gy
% 2. method=[x, y, z]  matrices are independently normalised to the dose value at point (x,y,z) (in cm)
% 3. method=dmean      matrices are independently normalised to the mean dose value in VOI voi2use (key insensitive)
% 4. method=0          matrices are not normalised (default)
%
% dose_xmesh,dose_ymesh,dose_zmesh are OPTIONAL x-y-z coordinates of the center of the matrix voxels
% VOI is a cell array which contain the patients VOIs (OPTIONAL to use with option 3)
% voi2use is a vector pointing to the number of VOI (OPTIONAL to use with option 3)
%
% Examples:
%
% C=dicomrt_dosediff(A,B,0,dose_xmesh,dose_ymesh,dose_zmesh,VOI,0)
% Store in C the dose difference between B and A. C=(B-A).
%
% C=dicomrt_dosediff(A,B,60,method,dose_xmesh,dose_ymesh,dose_zmesh,VOI,voi2use)
% Store in C the dose difference between Bnorm and Anorm. C=(Bnorm-Anorm), where
% Bnorm and Anorm are normalised to 60Gy (100%).
%
% C=dicomrt_dosediff(A,B,[10.5 -18 7],method,dose_xmesh,dose_ymesh,dose_zmesh)
% Store in C the dose difference between Bnorm and Anorm. C=(Bnorm-Anorm), where
% Bnorm and Anorm are normalised to the respective dose values at dnorm=(10.5 -18 7).
%
% C=dicomrt_dosediff(A,B,'Dmean',method,dose_xmesh,dose_ymesh,dose_zmesh,VOI,3)
% Store in C the dose difference between Bnorm and Anorm. C=(Bnorm-Anorm), where
% Bnorm and Anorm are normalised to the respective mean dose values in VOI number 3
%
% See also: dicomrt_doseratio, dicomrt_loadmcdose
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check number of argument and set-up some parameters and variables
error(nargchk(2,8,nargin))

% Define parameter
voilookup=1;
filter=0;
dosenorm=0;

% Check case and set-up some parameters and variables
[doseone_dose_temp,type_doseone_dose,labeld1]=dicomrt_checkinput(doseone,1);
[dosetwo_dose_temp,type_dosetwo_dose,labeld1]=dicomrt_checkinput(dosetwo,1);
doseone_dose=dicomrt_varfilter(doseone_dose_temp);
dosetwo_dose=dicomrt_varfilter(dosetwo_dose_temp);

if exist('VOI')==1 & exist('voi2use')==1
    if strcmpi(type_doseone_dose,'mc')==1
        doseone_dose_temp=dicomrt_mask(VOI,doseone_dose_temp,dose_xmesh,dose_ymesh,dose_zmesh,voilookup,filter,'y');
    end
    if strcmpi(type_dosetwo_dose,'mc')==1
        dosetwo_dose_temp=dicomrt_mask(VOI,dosetwo_dose_temp,dose_xmesh,dose_ymesh,dose_zmesh,voilookup,filter,'y');
    end
end

% Create label for dose diff
%if type_doseone_dose=='mc' & type_dosetwo_dose=='rtplan'
%    dose_diff_label='rtplan -mc'
%elseif type_doseone_dose=='mc' & type_dosetwo_dose=='mc'
%    dose_diff_label='mc2 - mc1'
%elseif type_doseone_dose=='rtplan' & type_dosetwo_dose=='mc'
%    dose_diff_label='mc - rtplan'
%else
%    dose_diff_label='rtplan2 - rtplan1'
%end

if exist('method')==1
    if isnumeric(method)==1 & length(method)==1 & method~=0
        % normalize to a specific dose level passed through "method"
        doseone_dose=doseone_dose./method.*100;
        dosetwo_dose=dosetwo_dose./method.*100;
        diff_dose=dosetwo_dose-doseone_dose;
    elseif isnumeric(method)==1 & length(method)==1 & method==0
        % no normalisation is carried out
        diff_dose=dosetwo_dose-doseone_dose;
    elseif isnumeric(method)==1 & length(method)==3 & ...
            (exist('dose_xmesh')~=1 | exist('dose_ymesh')~=1 | exist('dose_zmesh')~=1)
        error('dicomrt_dosediff: This normalisation method requires mesh to be provided. Exit now!');
    elseif isnumeric(method)==1 & length(method)==3 & ...
            (exist('dose_xmesh')==1 | exist('dose_ymesh')==1 | exist('dose_zmesh')==1)
        % normalize to a specific point in 3D passed through "method"
        locx=dicomrt_findpointVECT(dose_xmesh,method(1));
        locy=dicomrt_findpointVECT(dose_ymesh,method(2));
        locz=dicomrt_findpointVECT(dose_zmesh,method(3));
        doseone_dose=doseone_dose./doseone_dose(locy,locx,locz).*100;
        dosetwo_dose=dosetwo_dose./dosetwo_dose(locy,locx,locz).*100;
        diff_dose=dosetwo_dose-doseone_dose;
    elseif ischar(method)==1 & strcmpi(method,'dmean')==1 & ...
            (exist('VOI')~=1 | exist('voi2use')~=1)
        error('dicomrt_dosediff: This normalisation method requires VOI and voi2use to be provided. Exit now!');
    elseif ischar(method)==1 & strcmpi(method,'dmean')==1 & ...
            (exist('VOI')==1 | exist('voi2use')==1)
        % normalize to dmean
        dmean_one=dicomrt_MDcal(doseone,0,dose_xmesh,dose_ymesh,dose_zmesh,VOI,voi2use);
        dmean_two=dicomrt_MDcal(dosetwo,0,dose_xmesh,dose_ymesh,dose_zmesh,VOI,voi2use);
        doseone_dose=doseone_dose./dmean_one.*100;
        dosetwo_dose=dosetwo_dose./dmean_two.*100;
        diff_dose=dosetwo_dose-doseone_dose;
    end
else
    % no normalisation is carried out
    diff_dose=dosetwo_dose-doseone_dose;
end

% Restore original variable format
[diff_dose]=dicomrt_restorevarformat(dosetwo_dose_temp,diff_dose);

% Label Plan and update time of creation
if iscell(diff_dose)==1
    diff_dose{1,1}{1}.RTPlanLabel=[diff_dose{1,1}{1}.RTPlanLabel,'-DDIFF'];
    diff_dose{1,1}{1}.RTPlanDate=date;
    time=fix(clock);
    creationtime=[num2str(time(4)),':',num2str(time(5))];
    diff_dose{1,1}{1}.RTPlanTime=creationtime;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% Animation Object member functions
% See also
%
% Animation             Animation class (wth 3D data, projection, camera info ...)
% addNoise              Add noise to the animation
% centerW               Generate camera position from R and t in an Animation
% computeError          Compute several SFM errors
% display               Display the different elements contained in the object
% generateCamFromRt     Generate camera position from R and t in an Animation
% generateP             Compute the projection matrices using K, R and t
% generateS             Generate S tranformed by R and t
% generateSFromLSBasis  Generate 3D data from the basis of an Animation
% generateW             Generate projected data
% sampleFrame           Sample an Animation by reducing the number of frames
% setFirstPRtToId       Sets the first P or (R,t) of an Animation to eye(3,4)
% subsasgn              Set an element in the Animation object
% subsref               Reference an element in the Animation object

author: vrabaud


text:
if 0
	s = sparse([1 2],[1 2],[2 3]);
	%awf_sparse(int32([1 2]),int32([1 2]),[2 3])
	setnonzeros(s,[4 5])
	return
end

n = 1000;
s = rand(n,n);
[i,j,v] = find(s);
[m,n] = size(s);
tic, for iter = 1:10, s = sparse(i,j,v,m,n); end; t1=toc;
fprintf('time for sparse = %g\n', t1);
%tic, for iter = 1:10, s = awf_sparse(int32(i),int32(j),v); end; t1=toc;
%fprintf('time for awf_sparse = %g\n', t1);
tic, for iter = 1:10, s = setnonzeros(s,v); end; t2=toc;
fprintf('time for setnonzeros = %g (%g times faster)\n', t2, t1/t2);

author: andrewssobral

repo: lrslibrary


text:
%   ****** See nozzle.m for instructions ******
%  Program to extrapolate the data points from nozzle design and make a
%  uniform grid spacing in the x-direction
%  Change nothing... simply run this script

%  Find the minimum spacing given by the method of characteristics
%  and set as the dx value
dx = 0;
for i=1: size(noz,1)-1
    len = noz(i+1,1) - noz(i,1);
    if (len < dx || i == 1)
        dx = len;
    end
end

%  Explicitly give the dx value here
dx = max(noz(:,1))/ceil(max(noz(:,1))/dx);
n = max(noz(:,1))/dx;

% len = max(noz(:,1));
% n = 50;   %  Note # of points is actually n+1
% dx = len/n

%  Pick m points in y as some factor of x points
yfactor = .8;
m = ceil(yfactor*n);

%  Make uniform x-distribution of points
xmax = 0;
i = 1;
while xmax < max(noz(:,1))
    xmax = dx*(i-1);
    x(i,1:m) = xmax;
    i = i+1;
end

%  Make the y-points and extrapolate linearly from closest points to fit
%  the nozzle geometry

%  Initialize and assign last value
y(1:size(x,1),1:size(x,2)) = 0;
y(1,size(y,2)) = noz(1,size(noz,2));
y(size(y,1),size(y,2)) = noz(size(noz,1),size(noz,2));
for i = 1 : size(x,1)-1

    j = 1;
    while x(i,1) >= noz(j,1)
        x1 = noz(j,1);
        x2 = noz(j+1,1);
        y1 = noz(j,2);
        y2 = noz(j+1,2);
        j = j + 1;
    end

    slope = (y2 - y1)/(x2 - x1);
    y(i,size(y,2)) = y1 + slope*(x(i,1)-x1);

    %Fill in mesh
    dy = y(i,size(y,2))/(size(y,2)-1);
    for k = 1: size(y,2)
        y(i,k) = dy*(k-1);
    end

end

%Fill in mesh
dy = y(size(y,1),size(y,2))/(size(y,2)-1);
for k = 1: size(y,2)
    y(size(y,1),k) = dy*(k-1);
end

'Mesh done'

author: Sable

repo: mcbench-benchmarks


text:
% Assigns the given value to the named property of an object or structure.
% This function can deal with nested properties.
%
% Input arguments:
% obj:
%    the structure, handle or value object the value should be assigned to
% name:
%    a property name with dot (.) separating property names at
%    different hierarchy levels
% value:
%    the value to assign to the property at the deepest hierarchy
%    level
%
% Output arguments:
% obj:
%    the updated object or structure, optional for handle objects
%
% Example:
%    obj = struct('surface', struct('nested', 10));
%    obj = nestedassign(obj, 'surface.nested', 23);
%    disp(obj.surface.nested);  % prints 23
%
% See also: nestedfetch

% Copyright 2010 Levente Hunyadi
function obj = nestedassign(obj, name, value)
    if ~iscell(name)
        nameparts = strsplit2(name, '.');
    else
        nameparts = name;
    end
    obj = nestedassign_recurse(obj, nameparts, value);
end

function obj = nestedassign_recurse(obj, name, value)
% Assigns the given value to the named property of an object.
%
% Input arguments:
% obj:
%    the handle or value object the value should be assigned to
% name:
%    a cell array of the composite property name
% value:
%    the value to assign to the property at the deepest hierarchy
%    level
    if numel(name) > 1
        obj.(name{1}) = nestedassign_recurse(obj.(name{1}), name(2:end), value);
    else
        obj.(name{1}) = value;
    end
end


text:
function analyse_dev_lin_fusion(key,scores_obj_array,scr_lin,prior)
% Analyses the dev scores for all input systems and their fusion
% so that their relative strengths can be compared.  It calculates
% the percentage of target trials calculated by each system and the
% actual dcf, min dcf and prbep values for each system.
% Inputs:
%   key: The Key for the dev trials.
%   scores_obj_array: An array of score objects (one object for
%     each system) that was fused.
%   scr_lin: An object of type Scores containing the result of the
%     linear fusion.
%   prior: The effective target prior.

assert(nargin==4)
assert(isa(key,'Key'))
assert(key.validate())

stacked_scores = stackScores(key,scores_obj_array);
trialmask = key.to_ndx().trialmask;
dev_lin_fusion = scr_lin.scoremat(trialmask(:))';
dev_scores = [stacked_scores;dev_lin_fusion];

true_prop = sum(key.tar(:))/(sum(key.tar(:))+sum(key.non(:)));
logprint(Logger.Info,'True Dev target proportion is %g%%. System estimates are:\n',100*true_prop);
logprint(Logger.Info,'%s',sprintfmatrix(100*estimate_target_proportions(dev_scores,true_prop)'));

numsystems = size(stacked_scores,1);
mf = size(dev_scores,1);
logprint(Logger.Info,'\n\nAnalysing Dev: %i subsystems and %i fusions:\n',numsystems,mf-numsystems);
E = evalAll(dev_scores,key,prior);
logprint(Logger.Info,'\nDev all: dcf*100, mindcf*100, prbep:\n%s',sprintfmatrix(E,1:mf));
end

author: nesl


text:
function [C,CompositeObjectiveValue, SingleObjectiveValues]=ITERevaluate_X0...
    (Problem,Constraints,Spec,Objectives,askflag,printflag)

if nargin < 5; askflag = true; end
if nargin < 6; printflag = true; end
if askflag
choice = txtmenu('What do you want to evaluate?',...
    'Current design point','Lower bound','Upper bound');
% choice = 0;
switch choice
    case 0
        x0 = Problem.x0(Problem.activeX);
    case 1
        x0 = Problem.lb(Problem.activeX);
    case 2
        x0 = Problem.ub(Problem.activeX);
end
else
    x0 = Problem.x0(Problem.activeX);
end

[C] = EVALnonlcons(x0,Problem,Constraints,Spec);

if printflag
jnk=Constraints.conLabels(Constraints.active_cons);
fprintf('%30s\t\t%s\n','CONSTRAINT','MARGIN (%)');
for ii = 1:nnz(Constraints.active_cons)
    fprintf('%30s\t\t%-3.1f\n',jnk{ii},-100*C(ii));
end
disp(' ')
end

[CompositeObjectiveValue, SingleObjectiveValues] = EVALobjective...
    (x0,Problem,Objectives,Spec,0,'single');

if printflag
jnk=Objectives.ObjLabels;

fprintf('%15s\t\t%10s\t\t%s\n','OBJECTIVE','WEIGHTING','VALUE');
for ii = 1:Objectives.nObj
    if Objectives.weightings(ii)
     fprintf('%15s\t\t%9.2f\t\t%-3.2d\n',jnk{ii},...
         Objectives.weightings(ii),SingleObjectiveValues(ii));
    end
end

disp(' ')
fprintf('Composite Objective Value:   %d\n',CompositeObjectiveValue')
end


text:
function [V, S, varexp, w, Yhat] = plssquash(X, Y, varargin)
% Decomposes data X into K components that are ordered in their
% covariance with Y, designed to predict orthogonal parts of Y
%
% :Usage:
% ::
%
%     [V, S, varexp, w, Yhat] = plssquash(X, Y, varargin)
%
% :Optional Inputs:
%   - case {'noplot'}, turn off plotting
%   - case 'ndims', save only first ndims (K) vectors
%
% :Outputs:
%
%   **V:**
%        'eigenvetors', or weights, on data (columns)
%
%   **S:**
%        score matrix, N x K
%
%   **varexp:**
%        sqrt(r-square) with first k components predicting Y
%
%   **w:**
%        V*b, integrated weights.  for predicting new data, pred = X*w
%
%   **Yhat:**
%        X*V*b, or S*b
%
% ..
%    Tor Wager, 9/12/09
% ..

doplot = 1;
ndims = length(Y) - 1;

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            case {'noplot'}, doplot = 0;
            case 'ndims', ndims = varargin{i + 1};

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

% univariate covariance-based weights
% imperfect prediction
X = X - repmat(mean(X), size(X, 1), 1); %scale(X);
Y = scale(Y);

%
% neither V nor S appears to have orthogonal columns, though many are
% later columns, as df in Y is approached, appear to be highly colinear
% this may be because all variance in Y is basically explained...

clear V  % V is voxel weights for each component 1:K
clear S
rY = Y; % Y values to successively predict; intialize to Y

for i = 1:ndims  %ceil(length(Y)./2)

    % V = voxel weights, based on univariate relationship with rY
    V(:,i) = (X' * rY); %.^ .5;
    V(:,i) = V(:,i) ./ norm(V(:,i));

    % score matrix, N x K
    % these will ultimately be predictors for Y
    % chosen to maximize predictive power with few orthogonal components
    S = X*V;

    rY = Y - S * pinv(S) * Y;

    b = pinv(S) * Y;
    Yhat = S * b;

    rsq(i) = 1 - var(rY) ./ var(Y);

    if doplot
        create_figure('Fit'); plot(Yhat, Y, 'kx'); refline;

        rsq_adj(i) = 1 - (1-rsq(i)) * ((size(Y,1)-1) ./ (size(Y,1)-i-1));

        title(sprintf('cum. r-squared = %3.1f, adj = %3.1f', 100*rsq(i), 100*rsq_adj(i)));
        drawnow; pause(.3)
    end

end

varexp = sqrt(rsq);

b = pinv(S) * Y;
Yhat = S * b;
w = V * b;  % final weight vector

%prediction = mean(training) + data*S*b
%Yhat = X * V * b;
%figure; plot(Yhat, Y, 'kx')

end

author: canlab

repo: CanlabCore


text:
function dig = fif_read_headpoints(fid, tree)
% FIF_READ_HEADPOINTS: Read digitized head points from a FIF file.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@

% FIFF Constants
global FIFF;
if isempty(FIFF)
   FIFF = fiff_define_constants();
end

% Get DIG points
isotrak = fiff_dir_tree_find(tree, FIFF.FIFFB_ISOTRAK);
dig=struct('kind',{},'ident',{},'r',{},'coord_frame',{});
if length(isotrak) == 1
    p = 0;
    for k = 1:isotrak.nent
        kind = isotrak.dir(k).kind;
        pos  = isotrak.dir(k).pos;
        if kind == FIFF.FIFF_DIG_POINT
            p = p + 1;
            tag = fiff_read_tag(fid,pos);
            dig(p) = tag.data;
	    dig(p).coord_frame = FIFF.FIFFV_COORD_HEAD;
        end
    end
end

author: brainstorm-tools


text:
function view = saveGlmMap(view, fieldName, condNum, mapName, scan)
%
% view = saveGlmMap(view, fieldName, condNum, [mapName], [scan]);
%
% Load a GLM model for a mrVista view (saved for the specified scan),
% and save the specifield field name as a parameter map.
%
% To use this, you should have run 'applyGlm' on the view, to save
% the GLM data files.
%
% ras, 1/03/06
if ieNotDefined('view'), view = getSelectedInplane; end
if ieNotDefined('mapName'), mapName = fieldName; end
if ieNotDefined('scan'), scan = getCurScan(view); end

% initialize map as zeros
map = cell(1, numScans(view));
mapSize = dataSize(view);
map{scan} = zeros(mapSize);

h = mrvWaitbar(0, 'Saving GLM Map...');

for slice = 1:numSlices(view)
   % load the results of the GLM
   model = loadGlmSlice(view, slice, scan);

   % find the indices in the map which correspond to this slice
   [X Y Z] = meshgrid(1:mapSize(1), 1:mapSize(2), slice);
   coords = [X(:) Y(:) Z(:)]';
   ind = sub2ind(mapSize, coords(1,:), coords(2,:), coords(3,:));

   % get the data from the appropriate field
   if ~isfield(model, fieldName)
       error('field not found')
   elseif isequal(lower(fieldName), 'residual')
       % take the mean residual from each voxel
       vals = mean(model.(fieldName), 1);
   else
       vals = model.(fieldName)(:, condNum, :);
   end

   % map the values from the specified field to the appropriate
   % place in the map:
   map{scan}(ind) = vals;

   mrvWaitbar(slice/numSlices(view), h);
end

close(h);

% set the absolute value of the parameter as the 'co' field, so
% you can threshold by that as well:
co = cell(1, numScans(view));
co{scan} = abs(map{scan});

% save the map
mapPath = fullfile(dataDir(view), mapName);
save(mapPath, 'map', 'co', 'mapName');

% set as the active map in the view, so you can view the
% results right away:
view = setParameterMap(view, map, mapName);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [p] = mesh_write(p)

% mesh_write - a switch yard for several mesh writing functions
%
% Usage: [p] = mesh_write(p)
%
% p is a parameter structure (see eeg_toolbox_defaults for
% more details). In this function, it should contain mesh data
% in p.mesh.data and the following string fields:
%
%       p.mesh.path - the directory location of the file to load
%       p.mesh.file - the name of the file to load
%       p.mesh.type - the file formats (case insensitive) are:
%
%                    'emse'
%                    'brainstorm'
%                    'ascii' or 'freesurfer' for *.asc file
%                    'fs_surf' for freesurfer binary surface
%                    'fs_curv' for freesurfer binary curvature
%                    'fs_overlay' for freesurfer binary overlay (.w)
%
% If you have a file format you would like supported, please email
% the author with an example and description of the format.  The
% file formats supported here are described in more detail at
% their respective websites:
%
%       EMSE:       http://www.sourcesignal.com
%       BrainStorm: http://neuroimage.usc.edu/brainstorm/
%       FreeSurfer: http://surfer.nmr.mgh.harvard.edu/
%
% See also, mesh_open, mesh_write_emse,
%           mesh_write_brainstorm, mesh_write_freesurfer
%           mesh_write_fs_surf, mesh_write_fs_curv,
%           mesh_write_fs_overlay
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:58 $

% Licence:  GNU GPL, no express or implied warranties
% History:  02/2002 Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('p','var'),
    msg = sprintf('\nMESH_WRITE...no input p struct.\n');
    error(msg);
elseif isempty(p),
    msg = sprintf('\nMESH_WRITE...input p struct is empty.\n');
    error(msg);
elseif isempty(p.mesh.data),
    msg = sprintf('\nMESH_WRITE...input p struct has no mesh data.\n');
    error(msg);
end

type = lower(p.mesh.type);

switch type,

case 'emse',

    mesh_write_emse(p);

case 'brainstorm',

    mesh_write_brainstorm(p);

case {'ascii','freesurfer'}, % for *.asc file

    mesh_write_freesurfer(p);

case 'fs_surf', % for freesurfer binary surface

    mesh_write_fs_surf(p);

case 'fs_curv', % for freesurfer binary curvature

    mesh_write_fs_curv(p);

case 'fs_overlay', % for freesurfer binary overlay (.w)

    mesh_write_fs_overlay(p);

otherwise,
    fprintf('...mesh format: %s\n', p.mesh.type);
    fprintf('...sorry, cannot write this data format at present.\n');
    return;
end

return

author: PatternRecognition

repo: OpenBMI


text:
function a = r8cvv_rset ( mn, a, m, roff, i, ai )

%*****************************************************************************80
%
%% R8CVV_RSET sets row I from an R8CVV.
%
%  Discussion:
%
%    An R8CVV is a "vector of vectors" of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 November 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer MN, the size of the cell array.
%
%    Input, real A(MN), the cell array.
%
%    Input, integer M, the number of rows in the array.
%
%    Input, integer ROFF(M+1), the row offsets.
%
%    Input, integer I, the row.
%    1 <= I <= M.
%
%    Input, real AI(NR(I)), the new value of A(I,*).
%
%    Output, real A(MN), the cell array.
%
  k1 = roff(i) + 1;
  k2 = roff(i+1);
  nv = k2 + 1 - k1;
  a(k1:k2) = ai(1:nv);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [A xy] = wheel_graph(n)
% WHEEL_GRAPH Construct a wheel graph of order n
%
% The wheel graph is a cycle graph of order n-1 along with an additional
% vertex that connects all the remaining vertices.  (Run the example and it
% will be extremely clear if you are still confused.)
%
% [A xy] = wheel_graph(n) returns the adjacency matrix for the wheel graph
% of order n.  The matrix xy stores two-dimensional coordinates for each
% vertex.
%
% Example:
%   [A xy] = wheel_graph(10);
%   gplot(A,xy);
%
% See also CYCLE_GRAPH, STAR_GRAPH

% David Gleich
% Copyright, Stanford University, 2007-2008

%% History
%  2007-09-29: Changed output to double, fixed output for i=0,1
%%

if n>1
    i = 1:(n-1);
    j = [i(2:end) i(1)];
    i = [i i];
    j = [j n*ones(1,n-1)];
    A = sparse(i,j,1,n,n);
    A = A|A';
    A = double(A);
else
    i=[];
    A = sparse(n,n);
end

i = i(1:end/2);
if n>0
    xy = [ [cos(2*pi*(i./(n-1))) 0]' [sin(2*pi*(i./(n-1))) 0]' ];
else
    xy = zeros(0,2);
end


text:
% -------------------------------------------------------------------------------------------------------------------------
function net = remove_layers_from_prefix(net, prefix)
% -------------------------------------------------------------------------------------------------------------------------
    L = net.layers;
    num_layers = numel(L);
    to_remove = {};
    for i = 1:num_layers
        if strfind(L(i).name, prefix)
            to_remove{end+1} = L(i).name;
        end
    end
    net.removeLayer(to_remove);
end

author: bertinetto

repo: siamese-fc


text:
function bc = specific_bc(xbd,ybd)
%collide_bc   Reference problem 5.4 boundary condition
%   bc = specific_bc(xbd,ybd);
%   input
%          xbd          x boundary coordinate vector
%          ybd          y boundary coordinate vector
%
%   specifies streamfunction associated with colliding flow
%   IFISS function: DJS; 6 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
bc=-(xbd.^5)+5*xbd.*ybd.^4;
return


text:
function [XTrain XTest yTrain yTest] = splitData(X, y)

% Splits the data into training and testing

% Useful Variables
totalEgs = size(X,1);
numTrain = 0.8*totalEgs; % 80% data reserved for training
numTest = 0.2*totalEgs; % 20% reserved for testing
numLabels = unique(y); % how many unique labels present
TrainPerLable = numTrain/length(numLabels); % how many egs per lable in training
TestPerLable = numTest/length(numLabels);  % how many egs per lable in test

for i=1:length(numLabels)
    temp = (y==numLabels(i,1)); % find out where each label is present
    idx = find(temp==1); % get the index of that row where the lable is present
    Xtemp = X(idx(:,1),:); % get values of X stored at that particular index
    idxTemp = (randperm(size(Xtemp,1)))'; % randomly choose rows to put into train
    XTrainTemp = Xtemp(idxTemp(1:TrainPerLable,1),:);
    yTrainTemp = (ones(TrainPerLable,1))*i;
    XTestTemp = Xtemp(idxTemp(TrainPerLable+1:end),:); % select remaining rows for testing
    yTestTemp = (ones(TestPerLable,1))*i;

    if i==1 % set up first input to train and test
        XTrain = XTrainTemp;
        yTrain = yTrainTemp;
        XTest = XTestTemp;
        yTest = yTestTemp;
    else % keep adding new egs to the training and testing sets
        XTrain = [XTrain; XTrainTemp];
        yTrain = [yTrain; yTrainTemp];
        XTest = [XTest; XTestTemp];
        yTest = [yTest; yTestTemp];
    end
end


text:
classdef nme_bus_acp < mp.nme_bus & mp.form_acp

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

%     properties
%         name = 'bus';
%     end

    methods
        function obj = add_vvars(obj, nm, dm, idx)
            dme = obj.data_model_element(dm);
            nb = obj.nk;

            %% prepare angle bounds for ref buses
            va_lb = -Inf(nb, 1);
            va_ub =  Inf(nb, 1);
            k = find(dme.type == mp.NODE_TYPE.REF);
            va_lb(k) = dme.va_start(k);
            va_ub(k) = dme.va_start(k);

            nm.add_var('va', 'Va', nb, dme.va_start, va_lb, va_ub);
            nm.add_var('vm', 'Vm', nb, dme.vm_start, dme.vm_lb, dme.vm_ub);
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%% RES = factorial(NUM)
%
% Factorial function that works on matrices (matlab's does not).

% EPS, 11/02

function res = factorial(num)

res = ones(size(num));

ind = find(num > 0);
if ( ~isempty(ind) )
  subNum = num(ind);
  res(ind) = subNum .* factorial(subNum-1);
end

repo: SelfExSR


text:
classdef io
  % a class for input and output functions

  methods (Static)
    function im_data = load_image(im_file)
      % im_data = load_image(im_file)
      %   load an image from disk into Caffe-supported data format
      %   switch channels from RGB to BGR, make width the fastest dimension
      %   and convert to single
      %   returns im_data in W x H x C. For colored images, C = 3 in BGR
      %   channels, and for grayscale images, C = 1
      CHECK(ischar(im_file), 'im_file must be a string');
      CHECK_FILE_EXIST(im_file);
      im_data = imread(im_file);
      % permute channels from RGB to BGR for colored images
      if size(im_data, 3) == 3
        im_data = im_data(:, :, [3, 2, 1]);
      end
      % flip width and height to make width the fastest dimension
      im_data = permute(im_data, [2, 1, 3]);
      % convert from uint8 to single
      im_data = single(im_data);
    end
    function mean_data = read_mean(mean_proto_file)
      % mean_data = read_mean(mean_proto_file)
      %   read image mean data from binaryproto file
      %   returns mean_data in W x H x C with BGR channels
      CHECK(ischar(mean_proto_file), 'mean_proto_file must be a string');
      CHECK_FILE_EXIST(mean_proto_file);
      mean_data = caffe_('read_mean', mean_proto_file);
    end
    function write_mean(mean_data, mean_proto_file)
      % write_mean(mean_data, mean_proto_file)
      %   write image mean data to binaryproto file
      %   mean_data should be W x H x C with BGR channels
      CHECK(ischar(mean_proto_file), 'mean_proto_file must be a string');
      CHECK(isa(mean_data, 'single'), 'mean_data must be a SINGLE matrix');
      caffe_('write_mean', mean_data, mean_proto_file);
    end
  end
end

author: aimerykong

repo: deepImageAestheticsAnalysis


text:
% collapseDim() - collapse a range of values across a chosen dimension
%                 using the chosen method (e.g., average, peak, integral).
%
% Inputs:
%   C           -   [Nchs x Nchs x DIM3 x Dim4 x ... ] connectivity matrix.
%
% Optional:
%   'dim'         -   the dimension to collapse {def: 3}
%   'range'       -   [min-idx max-idx] index range to collapse.
%                     If range is a single value, than the layer for that
%                     index only is returned. {def: collapse entire dim}
%   'method'      -   ['mean' | 'net' | 'peak']. 'net' applies the trapz()
%                     integral ('dx' option must be set) {def: 'mean'}
%                 -   if method=@fcn, then function fcn() will be applied.
%                     This function must take a matrix as first argument
%                     and a dimension across which to evaluate as the
%                     second arg.
%   'dx'          -   [real] Required if 'method'='net'. This is the
%                     spacing increment (in same units as dataset)for trapz
%                     e.g., if integrating across frequency with freq.
%                     spacing of 0.5Hz, then dx=0.5;
%   'peakWidth'   -   [real] The width (in points) over which to identify
%                     a local peak. e.g., if peakWidth is 10, C(n,n,f) is
%                     only a peak if C(n,n,f:f-5)<C(n,n,f)>C(n,n,f:f+5)
%   'numPeaks'    -   [real] maximum number of peaks to return [if
%                     method='peaks' or 'peaks2d']
% Output:
%   C       -   collapsed [Nchs x Nchs x ... x 1 x ...] matrix with
%               dimension 'dim' reduced to a singleton
%   peaks   -   [Nchs x Nchs x ... x (numdims(C)-1) x Npeaks]
%               matrix containing indicies of peaks along the collapsed
%               dimension.
%
% Example:  collapse C(:,:,10:50) using integral with unit spacing
%           C2 = collapseDim(C,'dim',3,'range',[10 50],'method','net','dx',1)
%
% SEE ALSO
% localmax(), CAT
%
% Copyright (c) 2008 by Tim Mullen
function [C peaks] = collapseDim(C,varargin)

if nargin<2
    help collapseFreqs;
    return;
end

peaks = [];

range = [];
method = 'mean';
peakWidth = 20;
dim = 3;
dx = 0;
mph = -Inf;
numPeaks = 1;
keepedges = 1;
extent = 0;

for i=1:2:length(varargin)-1
    switch varargin{i}
        case 'range',       range = varargin{i+1};
        case 'method',      method = varargin{i+1};
        case 'dx',          dx = varargin{i+1};
        case 'dim',         dim = varargin{i+1};
        case 'peakWidth',   peakWidth = varargin{i+1};
        case 'numPeaks',    numPeaks = varargin{i+1};
        case 'minpeakheight', mph = varargin{i+1};
        case 'keepedges',   keepedges = varargin{i+1};
        case 'extent',      extent = varargin{i+1};
    end
end

if size(C,dim)==1
    return;
end

maxrange = size(C,dim);

if dim>length(size(C))
    error('The matrix does not have this many dimensions!');
end

% collapse over full range?
if isempty(range)
    range = [1 maxrange];
end

if length(range)>2
    error('range must be singleton or a 2-element [min max] vector');
end

% handle when single freq. input
if length(range)==1,
    C = getRange(C,dim,[range range]);
    return;
end

% range checking
if range(1)<1 || range(2)>maxrange
    error(sprintf(['range out of bounds\n' ...
                   'range must be between %d and %d'],...
                   1,maxrange));
   return;
end

switch lower(method)
    case 'mean'
        C = nan_mean(getRange(C,dim,range),dim);
    case 'net'
        if ~dx, error('''dx'' must be specified for ''net'' method'); end
        C = trapz(getRange(C,dim,range),dim);
        C = C.*dx;
    case 'peak'     % 1D peak detection
        [C peaks] = findConnPeaks(getRange(C,dim,range),dim, mph,peakWidth,numPeaks,keepedges,extent);
    case 'peak2d'   % 2D peak detection
%         C = getRange(C,dim,range);
        newC = zeros([size(C,1) size(C,2)],'single');
        peaks.freqs = zeros(size(C,1),size(C,2),numPeaks);
        peaks.times = zeros(size(C,1),size(C,2),numPeaks);
        if issymmetric(C)    % only estimate peaks for upper triangle (faster)
            for i=1:size(C,1)
                for j=i:size(C,2)
                    [pkval peaks.freqs(i,j,:) peaks.times(i,j,:)] = findpeaksND(squeeze(C(i,j,:,:)),peakWidth, true,numPeaks);   %max([size(C,3) size(C,4)])
                    newC(i,j) = pkval; %max(pkval);  % only return the maximum peak
                    newC(j,i) = newC(i,j);   % symmetric, so fill in lower triangle too
                end
            end
        else
            for i=1:size(C,1)
                for j=1:size(C,2)  % all channels

                    % new stuff
                    [pk] = imregionalmax(squeeze(C(i,j,:,:)));
                    if all(pk(:))
                        continue;  % no peak here
                    else
                        cc = squeeze(C(i,j,:,:));
                        [rows cols] = size(cc);
                        if extent
                            [ii jj] = find(pk);
                            for pp=1:length(ii)  % for each peak
                                % check if it's extent-neighborhood is
                                % greater than mph

                                if all(cc(setdiff_bc(max(1,ii-extent):min(rows,ii+extent),ii),setdiff_bc(max(1,jj-extent):min(cols,jj+extent),jj))<=mph)
                                    pk(ii,jj)=0;
                                end
                            end
                            if ~any(pk(:)), continue; end
                        end

                        cc = cc.*pk;
                        [newC(i,j) ii] = max(cc(:));
                        [peaks.freqs(i,j) peaks.times(i,j)] = ind2sub(size(cc),ii);

                    end

                end
            end
        end
        C = newC;

    case 'maxmag'
        % 1-D max magnitude (max of absval) along dimension dim
        [C peaks] = max(abs(getRange(C,dim,range)),[],dim);
    case 'max'
        % 1-D max along dimension dim
        [C peaks] = max(getRange(C,dim,range),[],dim);
    case 'min'
        % 1-D min along dimension dim
        [C peaks] = min(getRange(C,dim,range),[],dim);
    case {'getrange' 'shrinkonly'}
        C = getRange(C,dim,range);
    otherwise
        try
            C = feval(method,getRange(C,dim,range),dim);
        catch
            error('could not evaluate function');
        end
end

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

% extract a range of values from a specified dim
function C = getRange(C,dim,range)

numdims = length(size(C));

st = '';
for i=1:numdims
    st = fastif(dim==i,[st ' range(1):range(2),'],[st ':,']);
end
st(end) = [];

eval(['C = C( ' st ');']);

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

function [Cout pidxs] = findConnPeaks(C,dim,mph,mpd,npeaks,keepedges,extent)
    % C is [nchs nchs dim3 dim4]
    % dim is the dimension to find 1-d peaks along (this dimension will be
    % collapsed to a singleton)
    % output is same size as C with singleton dimension dim
    % mph is the minimum peak heights
    % mpd is the minimum distance between peaks
    % npeaks is the max number of peaks to return
    % keepedges determines whether or not to count a max at the series edge
    %   as a peak
    % extent determines the number of pixels adjecent to a peak that must
    % be greater than mph

    if nargin<6
        keepedges = 0;
    end

    if nargin<7
        extent = 0;
    end

    warning off

    verb = 0;

    if nargin<5
        npeaks = 10;    % number of top peaks to store in pidxs
    end

    % permute C so we are finding peaks over last (4th) dimension
    permutation = [1 2 fastif(dim==3,[4 3],[3 4])];
    C=permute(C,permutation);

    [rows cols nd3 nd4] = size(C);
    pidxs = zeros(rows,cols,nd3,npeaks);
    Cout = zeros(rows,cols,nd3);

    if nd3>10
        fprintf('warning: searching for peaks over %d series -- this may take a while...\n',nd3);
        verb=1;
    end
    if verb, h=waitbar(0,'searching for peaks...'); cur=1; end
    for kk=1:nd3
        for ch1=1:rows
            for ch2=1:cols
                series = [squeeze(C(ch1,ch2,kk,:))];
                if keepedges==0
                    [peaks locs] = findpeaks(series,length(series),extent,0,mph);
                else
                    [peaks locs] = findpeaks([0 ; series ; 0],length(series),extent,0,mph);
                end
                locs = round(locs);

                if ~isempty(locs)
                    pidxs(ch1,ch2,kk,1:min(length(locs),npeaks)) = ...          % keep indices of npeaks largest peak
                        locs(1:min(length(locs),npeaks))-keepedges;             % -1 to counteract zero-padding
                    Cout(ch1,ch2,kk)=peaks(1);                                  % store largest peak value
                else
                    pidxs(ch1,ch2,kk,:)=nan;
                    Cout(ch1,ch2,kk)=nan;
                end

                if verb, waitbar(cur/(nd3*rows*cols),h); cur=cur+1; end

            end
        end

        if ~mod(kk,10), fprintf('.'); end
    end

    if verb, close(h); end
    warning on;

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

function sym = issymmetric(C)
% check if C is symmetric along the first two dimensions

sym = 1;

sz = size(C);

for i=1:sz(1)
    for j=1:sz(2)
        if ~isequal(C(i,j,:,:,:,:),C(j,i,:,:,:,:))
            sym = 0;
            return;
        end
    end
end