text:
% STD_CHANTOPO - plot ERP/spectral/ERSP topoplot at a specific
%                  latency/frequency.
% Usage:
%          >> std_chantopo( data, 'key', 'val', ...)
% Inputs:
%  data  -  [cell array] mean data for each subject group and/or data
%           condition. These arrays are usually returned by function
%           std_erspplot and std_erpplot. For example
%
%           >> data = { [1x64x12] [1x64x12 }; % 2 groups of 12 subjects, 64 channels
%           >> std_chantopo(data, 'chanlocs', 'chanlocfile.txt');
%
% Scalp map plotting option (mandatory):
%  'chanlocs'    - [struct] channel location structure
%
% Other scalp map plotting options:
%  'chanlocs'    - [struct] channel location structure
%  'topoplotopt' - [cell] topoplot options. Default is { 'style', 'both',
%                  'shading', 'interp' }. See topoplot help for details.
%  'ylim'        - [min max] ordinate limits for ERP and spectrum plots
%                  {default: all available data}
%  'caxis'       - [min max] same as above
%
% Optional display parameters:
%  'datatype'    - ['erp'|'spec'] data type {default: 'erp'}
%  'titles'      - [cell array of string] titles for each of the subplots.
%                  { default: none}
%  'subplotpos'  - [addr addc posr posc] perform plotting in existing figure.
%                  Add "addr" rows, "addc" columns and plot the scalp
%                  topographies starting at position (posr,posc).
%
% Statistics options:
%  'groupstats'  - [cell] One p-value array per group {default: {}}
%  'condstats'   - [cell] One p-value array per condition {default: {}}
%  'interstats'  - [cell] Interaction p-value arrays {default: {}}
%  'threshold'   - [NaN|real<<1] Significance threshold. NaN -> plot the
%                  p-values themselves on a different figure. When possible,
%                  significance regions are indicated below the data.
%                  {default: NaN}
%  'binarypval'  - ['on'|'off'] if a threshold is set, show only significant
%                  channels as red dots. Default is 'off'.
%
% Author: Arnaud Delorme, CERCO, CNRS, 2006-
%
% See also: POP_ERSPPARAMS, POP_ERPPARAMS, POP_SPECPARAMS, STATCOND

% Copyright (C) 2006 Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function std_chantopo(data, varargin)

pgroup = [];
pcond  = [];
pinter = [];
if nargin < 2
    help std_chantopo;
    return;
end

opt = finputcheck( varargin, { 'ylim'        'real'    []              [];
                               'titles'      'cell'    []              cell(20,20);
                               'threshold'   'real'    []              NaN;
                               'chanlocs'    'struct'  []              struct('labels', {});
                               'groupstats'  'cell'    []              {};
                               'condstats'   'cell'    []              {};
                               'interstats'  'cell'    []              {};
                               'effect'      'string' { 'main','marginal' }   'marginal';
                               'subplotpos'  'integer' []              [];
                               'topoplotopt' 'cell'    []              { 'style', 'both' };
                               'binarypval'  'string'  { 'on','off' }  'on';
                               'datatype'    'string'  { 'ersp','itc','erp','spec' }    'erp';
                               'caxis'       'real'    []              [] }, 'std_chantopo', 'ignore'); %, 'ignore');
if ischar(opt), error(opt); end
if ~isempty(opt.ylim), opt.caxis = opt.ylim; end
if isnan(opt.threshold), opt.binarypval = 'off'; end
if strcmpi(opt.binarypval, 'on'), opt.ptopoopt = { 'style' 'blank' }; else opt.ptopoopt = opt.topoplotopt; end

% remove empty entries
datapresent = ~cellfun(@isempty, data);
for c = size(data,1):-1:1, if sum(datapresent(c,:)) == 0, data(c,:) = []; opt.titles(c,:) = []; if ~isempty(opt.groupstats), opt.groupstats(c) = []; end; end; end
for g = size(data,2):-1:1, if sum(datapresent(:,g)) == 0, data(:,g) = []; opt.titles(:,g) = []; if ~isempty(opt.condstats ), opt.condstats( g) = []; end; end; end

nc = size(data,1);
ng = size(data,2);
if nc >= ng, opt.transpose = 'on';
else         opt.transpose = 'off';
end

% plotting parameters
% ------------------
if ng > 1 && ~isempty(opt.groupstats), addc = 1; else addc = 0; end
if nc > 1 && ~isempty(opt.condstats ), addr = 1; else addr = 0; end
if ~isempty(opt.subplotpos)
     if strcmpi(opt.transpose, 'on'), opt.subplotpos = opt.subplotpos([2 1 4 3]); end
     addr = opt.subplotpos(1);
     addc = opt.subplotpos(2);
     posr = opt.subplotpos(4);
     posc = opt.subplotpos(3);
else posr = 0;
     posc = 0;
end

% compute significance mask
% -------------------------
pinterplot = {};
if strcmpi(opt.effect, 'marginal') || ng == 1 || nc == 1
    if ~isnan(opt.threshold) && ( ~isempty(opt.groupstats) || ~isempty(opt.condstats) )
        pcondplot  = opt.condstats;
        pgroupplot = opt.groupstats;
        maxplot = 1;
    else
        for ind = 1:length(opt.condstats),  pcondplot{ind}  = -log10(opt.condstats{ind}); end
        for ind = 1:length(opt.groupstats), pgroupplot{ind} = -log10(opt.groupstats{ind}); end
        maxplot = 3;
    end
elseif strcmpi(opt.effect, 'main') && ~isempty(opt.interstats)
    if ~isnan(opt.threshold) && ( ~isempty(opt.groupstats) || ~isempty(opt.condstats) )
        pcondplot  = { opt.interstats{2} };
        pgroupplot = { opt.interstats{1} };
        pinterplot = opt.interstats{3};
        maxplot = 1;
    else
        if ~isempty(opt.interstats{2}), pcondplot  = { -log10(opt.interstats{2}) }; end
        if ~isempty(opt.interstats{1}), pgroupplot = { -log10(opt.interstats{1}) }; end
        if ~isempty(opt.interstats{3}), pinterplot = -log10(opt.interstats{3}); end
        maxplot = 3;
    end
end

% adjust figure size
% ------------------
if isempty(opt.subplotpos)
    fig = figure('color', 'w');
    pos = get(fig, 'position');
    set(fig, 'position', [ pos(1)+15 pos(2)+15 pos(3)/2.5*(nc+addr), pos(4)/2*(ng+addc) ]);
    pos = get(fig, 'position');
    if strcmpi(opt.transpose, 'off'), set(gcf, 'position', [ pos(1) pos(2) pos(4) pos(3)]);
    else                              set(gcf, 'position', pos);
    end
end

% topoplot
% --------
tmpc = [inf -inf];
for c = 1:nc
    for g = 1:ng
        hdl(c,g) = mysubplot(nc+addr, ng+addc, c, g, opt.transpose);
        if ~isempty(data{c,g})
            tmpplot = double(mean(data{c,g},3));
            if ~isreal(tmpplot(1)), tmpplot = abs(tmpplot); end % comes second for processing single trials
            if ~all(isnan(tmpplot))
                if ~isreal(tmpplot), error('This function cannot plot complex values'); end
                topoplot( tmpplot, opt.chanlocs, opt.topoplotopt{:});
                if isempty(opt.caxis)
                    tmpc = [ min(min(tmpplot), tmpc(1)) max(max(tmpplot), tmpc(2)) ];
                else
                    caxis(opt.caxis);
                end
                title(opt.titles{c,g}, 'interpreter', 'none');
            else
                axis off;
            end
        else
            axis off;
        end

        % statistics across groups
        % -------------------------
        if strcmpi(opt.effect, 'marginal') || (strcmpi(opt.effect, 'main') && c == 1)
            if g == ng && ng > 1 && ~isempty(opt.groupstats)
                if strcmpi(opt.effect, 'main') && nc>1, centerc = nc/2-0.5; else centerc = 0; end
                hdl(c,g+1) = mysubplot(nc+addr, ng+addc, c+centerc, ng + 1, opt.transpose);
                pgroupplot{c}(pgroupplot{c} <0) = 0;
                topoplot( pgroupplot{c}, opt.chanlocs, opt.ptopoopt{:});
                title(opt.titles{c,g+1}, 'interpreter', 'none');
                caxis([-maxplot maxplot]);
            end
        end
    end
end

% color scale
% -----------
if isempty(opt.caxis)
    for c = 1:nc
        for g = 1:ng
            axes(hdl(c,g));
            caxis(tmpc);
        end
    end
end

for g = 1:ng
    % statistics across conditions
    % -----------------------------
    if strcmpi(opt.effect, 'marginal') || (strcmpi(opt.effect, 'main') && g == 1)
        if ~isempty(opt.condstats) && nc > 1
            if strcmpi(opt.effect, 'main') && ng>1, centerg = ng/2-0.5; else centerg = 0; end
            hdl(nc+1,g) = mysubplot(nc+addr, ng+addc, nc+addr, g+centerg, opt.transpose);
            pcondplot{g}(pcondplot{g} < 0) = 0;
            topoplot( pcondplot{g}, opt.chanlocs, opt.ptopoopt{:});
            title(opt.titles{nc+1,g}, 'interpreter', 'none');
            caxis([-maxplot maxplot]);
        end
    end
end

% statistics across group and conditions
% ---------------------------------------
if ~isempty(opt.condstats) && ~isempty(opt.groupstats) && ng > 1 && nc > 1 && ~isempty(pinterplot)
    hdl(nc+1,ng+1) = mysubplot(nc+addr, ng+addc, nc+addr, ng+1, opt.transpose);
    pinterplot(pinterplot<0) = 0;
    topoplot( pinterplot, opt.chanlocs, opt.ptopoopt{:});
    title(opt.titles{nc+1,ng+1}, 'interpreter', 'none');
    caxis([-maxplot maxplot]);
end

% color bars
% ----------
if isnan(opt.threshold(1)) && (nc ~= size(hdl,1) || ng ~= size(hdl,2))
    try
        axes(hdl(end,end));
        cbar_signif(ng, maxplot);
    catch
        cbar_signif(ng, maxplot);
    end
end
axes(hdl(nc,ng));
cbar_standard(opt.datatype, ng);

% remove axis labels
% ------------------
for c = 1:size(hdl,1)
    for g = 1:size(hdl,2)
        if g ~= 1 && size(hdl,2) ~=1, ylabel(''); end
        if c ~= size(hdl,1) && size(hdl,1) ~= 1, xlabel(''); end
    end
end

% colorbar for ERSP and scalp plot
% --------------------------------
function cbar_standard(datatype, ng);
    pos = get(gca, 'position');
    tmpc = caxis;
    fact = fastif(ng == 1, 40, 20);
    tmp = axes('position', [ pos(1)+pos(3)+pos(3)/fact pos(2) pos(3)/fact pos(4) ]);
    set(gca, 'unit', 'normalized');
    if strcmpi(datatype, 'itc')
         cbar(tmp, 0, tmpc, 10); ylim([0.5 1]);
    else cbar(tmp, 0, tmpc, 5);
    end

% colorbar for significance
% -------------------------
function cbar_signif(ng, maxplot);
    % Retrieving Defaults
    icadefs;

    pos = get(gca, 'position');
    tmpc = caxis;
    fact = fastif(ng == 1, 40, 20);
    tmp = axes('position', [ pos(1)+pos(3)+pos(3)/fact pos(2) pos(3)/fact pos(4) ]);
    map = colormap(DEFAULT_COLORMAP);
    n = size(map,1);
    cols = [ceil(n/2):n]';
    image([0 1],linspace(0,maxplot,length(cols)),[cols cols]);
    %cbar(tmp, 0, tmpc, 5);
    tick = linspace(0, maxplot, maxplot+1);
    set(gca, 'ytickmode', 'manual', 'YAxisLocation', 'right', 'xtick', [], ...
        'ytick', tick, 'yticklabel', round(10.^-tick*1000)/1000);
    xlabel('');
    colormap(DEFAULT_COLORMAP);

% mysubplot (allow to transpose if necessary)
% -------------------------------------------
function hdl = mysubplot(nr,nc,r,c,subplottype)

    cmargin = 0.2/nc;
    rmargin = 0.2/nr;
    if strcmpi(subplottype, 'transpose') || strcmpi(subplottype, 'on'),   hdl = subplot('position',[(r-1)/nr+rmargin (nc-c)/nc+cmargin 1/nr-2*rmargin 1/nc-2*cmargin]);
    elseif strcmpi(subplottype, 'normal') || strcmpi(subplottype, 'off'), hdl = subplot('position',[(c-1)/nc+cmargin (nr-r)/nr+rmargin 1/nc-2*cmargin 1/nr-2*rmargin]);
    elseif strcmpi(subplottype, 'noplot'), hdl = gca;
    else error('Unknown subplot type');
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function test_issue1198

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_sourceanalysis ft_sourcedescriptives ft_sourceplot

%%
% http://www.fieldtriptoolbox.org/example/testing_bem_created_leadfields/

% Create a spherical volume conductor of radius 100, conductivity 1 and center [0,0,0]
vol   = [];
vol.r = 100;
vol.c = 1;
vol.o = [0 0 0];
vol.unit = 'mm';

% Create a set of electrodes on the upper half of the sphere
[X, Y, Z] = sphere(10);
pos = unique([X(:) Y(:) Z(:)], 'rows');
pos = pos(pos(:,3)>=0,:);

elec = [];
elec.elecpos = vol.r * pos;
elec.label = {};
nelec = size(pos,1);
for ii = 1:nelec
  elec.label{ii} = sprintf('vertex%03d', ii);
end

% Define positions of dipoles (along z axis)
zp = linspace(0,100,50)';
pos = [zeros(size(zp)) zeros(size(zp)) zp];

% Define the corresponding spatial grid
cfg = [];
cfg.inwardshift = 5;
cfg.sourcemodel.pos = pos;
cfg.headmodel = vol;
cfg.elec = elec;
sourcemodel = ft_prepare_sourcemodel(cfg);

% construct a dense triangulated mesh
[X, Y, Z] = sphere(100);
pos = unique([X(:) Y(:) Z(:)], 'rows');
mesh.pos = pos*100;
mesh.tri = convhulln(mesh.pos);

vertic = [50 80 100 200 2000];
for ll=1:length(vertic)
  cfg = [];
  cfg.method = 'dipoli';
  cfg.conductivity  = 1;
  cfg.numvertices = vertic(ll);
  cfg.isolatedsource = false;
  volbem{ll} = ft_prepare_headmodel(cfg, mesh);
end

% calculate BEM leadfield
for ll=1:length(vertic)
  cfg = [];
  cfg.sourcemodel = sourcemodel;
  cfg.headmodel = volbem{ll};
  cfg.elec = elec;
  leadfieldBEM{ll} = ft_prepare_leadfield(cfg);
end

% calculate theoretical single sphere leadfield
cfg = [];
cfg.grid = sourcemodel;
cfg.headmodel = vol;
cfg.elec = elec;
leadfieldSphere = ft_prepare_leadfield(cfg);

%%
% http://www.fieldtriptoolbox.org/example/compute_forward_simulated_data_and_apply_a_beamformer_scan/

% This example script shows you how to create some simulated channel-level
% MEG data with a single dipole at a specified location in the head.
% Subsequently it does a beamformer source reconstruction to localize that
% source.

% create an array with some magnetometers at 12cm distance from the origin
[X, Y, Z] = sphere(10);
pos = unique([X(:) Y(:) Z(:)], 'rows');
pos = pos(pos(:,3)>=0,:);
grad = [];
grad.coilpos = 12*pos;
grad.coilori = pos; % in the outward direction
% grad.tra = eye(length(pos)); % each coils contributes exactly to one channel
for i=1:length(pos)
  grad.label{i} = sprintf('chan%03d', i);
end

% create a spherical volume conductor with 10cm radius
vol.r = 10;
vol.o = [0 0 0];

% note that beamformer scanning will be done with a 1cm grid, so you should
% not put the dipole on a position that will not be covered by a grid
% location later
cfg = [];
cfg.headmodel = vol;
cfg.grad = grad;
cfg.dip.pos = [0 0 4];    % you can vary the location, here the dipole is along the z-axis
cfg.dip.mom = [1 0 0]';   % the dipole points along the x-axis
cfg.relnoise = 10;
cfg.ntrials = 20;
data = ft_dipolesimulation(cfg);

% compute the data covariance matrix, which will capture the activity of
% the simulated dipole
cfg = [];
cfg.covariance = 'yes';
timelock = ft_timelockanalysis(cfg, data);

% do the beamformer source reconstuction on a 1 cm grid
cfg = [];
cfg.headmodel = vol;
cfg.grad = grad;
cfg.resolution = 1;
cfg.method = 'lcmv';
cfg.lcmv.projectnoise = 'yes'; % needed for neural activity index
source = ft_sourceanalysis(cfg, timelock);

% compute the neural activity index, i.e. projected power divided by
% projected noise
cfg = [];
cfg.powmethod = 'none'; % keep the power as estimated from the data covariance, i.e. the induced power
source_nai = ft_sourcedescriptives(cfg, source);

cfg = [];
cfg.method = 'ortho';
cfg.funparameter = 'nai';
cfg.funcolorlim = [1.4 1.5];  % the voxel in the center of the volume conductor messes up the autoscaling
ft_sourceplot(cfg, source_nai);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function prob_test070 ( )

%*****************************************************************************80
%
%% TEST070 tests FACTORIAL_STIRLING, I4_FACTORIAL;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST070\n' );
  fprintf ( 1, '  FACTORIAL_STIRLING computes Stirling''s\n' );
  fprintf ( 1, '    approximate factorial function;\n' );
  fprintf ( 1, '  I4_FACTORIAL evaluates the factorial function;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  N      Stirling     N!\n' );
  fprintf ( 1, '\n' );

  for i = 0 : 20
    value = factorial_stirling ( i );
    fprintf ( 1, '  %6d  %14f  %20d\n', i, value, i4_factorial ( i ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function jed = ymdf_to_jed_eg_civil ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_JED_EG_CIVIL converts an Egyptian Civil YMDF date to a JED.
%
%  Discussion:
%
%    The Egyptian Civil calendar used a year of 365 days.  The year comprised
%    12 months of 30 days, with 5 epagomenal days occurring at the end of
%    the year.  Since the observed year is about 365.25 days long, and no
%    attempt was made to adjust the Egyptian Civil year to the observed year,
%    the calendar dates gradually drifted with respect to the observed dates.
%
%    The epoch or first day of the Egyptian Civil calendar is taken as
%    JED = 1448638.5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Algorithm E,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999, pages 323-324.
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, real JED, the corresponding Julian Ephemeris Date.
%

%
%  Convert the calendar date to a computational date.
%
  y_prime = y + 3968 - floor ( ( 13 - m ) / 13 );
  m_prime = mod ( m + 12, 13 );
  d_prime = d - 1;
%
%  Convert the computational date to a JED.
%
  jed = 365 * y_prime + 30 * m_prime + d_prime - 47 + 1 - 0.5;
  jed = jed + f;

  return
end


text:
function result = mergecellstruct(result)

% MERGECELLSTRUCT is a helper function for FT_TEST

if iscell(result)
  if numel(result)>1
    result = appendstruct(result{:});
  else
    result = result{1};
  end
end


text:
function [x_tplus1,z_tplus1,pseudo_data] = DIT_oneIter(y,z_t,x_t,width,height,denoiser,M_func,Mt_func)
% function x_hat = DIT_oneIter(y,z_t,x_t,width,height,denoiser,M_func,Mt_func)
% This function computes one iteration of the D-IT algorithm
% Input:
%       y       : the measurements
%       z_t     : current residual
%       x_t     : current signal estimate
%       iters   : the number of iterations
%       width   : width of the sampled signal
%       height  : height of the sampeled signal. height=1 for 1D signals
%       denoiser: string that determines which denosier to use. e.g.
%       denoiser='BM3D'
%       M_func  : function handle that projects onto M. Or a matrix M.
%       Mt_func : function handle that projects onto M'. Or no entry
%Output:
%       x_tplus1 : next signal estimate.
%       z_tplus1 : next residual

if nargin==8%function
    M=@(x) M_func(x);
    Mt=@(z) Mt_func(z);
else%Matrix
    M=@(x)M_func*x;
    Mt=@(z)M_func'*z;
end
denoi=@(noisy,sigma_hat) denoise(noisy,sigma_hat,width,height,denoiser);

n=width*height;
m=length(y);

pseudo_data=Mt(z_t)+x_t;
sigma_hat=sqrt(1/m*sum(abs(z_t).^2));
x_tplus1=denoi(pseudo_data,sigma_hat);
z_tplus1=y-M(x_tplus1);
end

author: ricedsp


text:
function [metrics, emphysema_mask] = PTKComputeEmphysemaFromMask(roi_data, mask)
    % PTKComputeEmphysemaFromMask. Computes emphysema percentage and percentile
    % density
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.

    emphysema_threshold_value_hu_1 = -950;
    emphysema_threshold_value_hu_2 = -910;

    emphysema_threshold_value_1 = roi_data.HounsfieldToGreyscale(emphysema_threshold_value_hu_1);
    emphysema_threshold_value_2 = roi_data.HounsfieldToGreyscale(emphysema_threshold_value_hu_2);
    emphysema_threshold_value_percentile = 15;
    lower_threshold = (roi_data.RawImage <= emphysema_threshold_value_1) & (mask.RawImage > 0);
    upper_threshold = (roi_data.RawImage <= emphysema_threshold_value_2) & (mask.RawImage > 0);
    emphysema_mask_raw = uint8(upper_threshold);
    emphysema_mask_raw(lower_threshold) = 3;
    emphysema_mask = mask.BlankCopy;
    emphysema_mask.ChangeRawImage(emphysema_mask_raw);

    number_of_voxels_in_mask = sum(mask.RawImage(:));
    emphysema_voxels_in_mask = sum(lower_threshold(:)); % We only count voxels in the lower threshold
    emphysema_percentage = 100*emphysema_voxels_in_mask/number_of_voxels_in_mask;

    if ~mask.ImageExists
        emphysema_percentile_density = NaN;
        emphysema_percentile_density_hu = NaN;
    else
        emphysema_percentile_density = prctile(roi_data.RawImage(mask.RawImage(:)), emphysema_threshold_value_percentile);
        emphysema_percentile_density_hu = roi_data.GreyscaleToHounsfield(emphysema_percentile_density);
    end

    metrics = PTKMetrics;
    metrics.AddMetric('EmphysemaPercentage', emphysema_percentage, '% of emphysema');
    metrics.AddMetric('EmphysemaPercentileDensityHU', emphysema_percentile_density_hu, 'Emphysema percentile density (HU)');
end

author: tomdoel

repo: pulmonarytoolkit


text:
function test_suite=test_moxunit_fieldtrip_util_parse_dependencies
% #   For MOxUnit_fieldtrip's copyright information and license terms,   #
% #   see the COPYING file distributed with MOxUnit_fieldtrip.           #

    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function varargout=randint(mx)
    n=nargout;
    varargout=arrayfun(@(unused)ceil(rand()*mx),1:n,'UniformOutput',false);

function assert_output(input, expected_output)
    output=moxunit_fieldtrip_util_parse_dependency(input);
    if isempty(output)
        assertTrue(isempty(expected_output));
    else
        assertEqual(output,expected_output,sprintf('input: ''%s''',input));
    end

function test_fieldtrip_util_parse_dependencies_basics()
    empty={};
    input_outputs={'foo',{'foo'};...
                  'foo  ',{'foo'};...               % whitespace
                  'foo bar',{'foo','bar'};...
                  'foo  bar  baz ',{'foo','bar','baz'};...
                  '',empty;...
                  };
    for k=1:size(input_outputs)
        input_output=input_outputs(k,:);
        assert_output(input_output{1},input_output{2});
    end

function moxunit_fieldtrip_util_parse_dependency_exceptions
    aet=@(varargin)assertExceptionThrown(@()...
            moxunit_fieldtrip_util_parse_dependency(varargin{:}),'');
    % inputs that are not string throw an exception
    aet(struct)
    aet([1 2])
    aet(1);


text:
function [slug] = lbm2slug(lbm)
% convert mass from pounds-mass to slugs.
slug = lbm/32.17405;

author: Sable

repo: mcbench-benchmarks


text:
function [elem,bdFlag] = sortelem(elem,bdFlag)
%% SORTELEM sort elem in ascend ordering
%
% [elem,bdFlag] = SORTELEM(elem,bdFlag) sorts the elem such that
% elem(t,1)< elem(t,2)< elem(t,3). A simple sort(elem,2) cannot
% switch bdFlag.
%
% See also  sortelem3
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Dimension check
if size(elem,2) >= 4 % 3D case
    [elem,bdFlag] = sortelem3(elem,bdFlag);
    return
end

%% Step 1: make elem(:,3) to the largest one
[tempvar,idx] = max(elem,[],2); %#ok<ASGLU>
elem(idx==1,1:3) = elem(idx==1,[2 3 1]);
elem(idx==2,1:3) = elem(idx==2,[3 1 2]);
if exist('bdFlag','var') && ~isempty(bdFlag)
    bdFlag(idx==1,1:3) = bdFlag(idx==1,[2 3 1]);
    bdFlag(idx==2,1:3) = bdFlag(idx==2,[3 1 2]);
end

%% Step 2: swtich the first two vertices such that elem(:,1)<elem(:2)
idx = (elem(:,2) < elem(:,1));
elem(idx,[1 2]) = elem(idx,[2 1]);
if exist('bdFlag','var') && ~isempty(bdFlag)
    bdFlag(idx,[1 2]) = bdFlag(idx,[2 1]);
else
    bdFlag = [];
end

repo: ifem


text:
function scale = p22_scale ( neqn )

%*****************************************************************************80
%
%% P22_SCALE returns scale factors for problem p22.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real SCALE(NEQN), the scaling factors.
%
  scale(1:neqn,1) = [ 1.0; 1.0 ];

  return
end


text:
% RECOMBIN.M       (RECOMBINation high-level function)
%
% This function performs recombination between pairs of individuals
% and returns the new individuals after mating. The function handles
% multiple populations and calls the low-level recombination function
% for the actual recombination process.
%
% Syntax:  NewChrom = recombin(REC_F, OldChrom, RecOpt, SUBPOP)
%
% Input parameters:
%    REC_F     - String containing the name of the recombination or
%                crossover function
%    Chrom     - Matrix containing the chromosomes of the old
%                population. Each line corresponds to one individual
%    RecOpt    - (optional) Scalar containing the probability of
%                recombination/crossover occurring between pairs
%                of individuals.
%                if omitted or NaN, 1 is assumed
%    SUBPOP    - (optional) Number of subpopulations
%                if omitted or NaN, 1 subpopulation is assumed
%
% Output parameter:
%    NewChrom  - Matrix containing the chromosomes of the population
%                after recombination in the same format as OldChrom.

%  Author:    Hartmut Pohlheim
%  History:   18.03.94     file created

function NewChrom = recombin(REC_F, Chrom, RecOpt, SUBPOP);

% Check parameter consistency
   if nargin < 2, error('Not enough input parameter'); end

   % Identify the population size (Nind)
   [Nind,Nvar] = size(Chrom);

   if nargin < 4, SUBPOP = 1; end
   if nargin > 3,
      if isempty(SUBPOP), SUBPOP = 1;
      elseif isnan(SUBPOP), SUBPOP = 1;
      elseif length(SUBPOP) ~= 1, error('SUBPOP must be a scalar'); end
   end

   if (Nind/SUBPOP) ~= fix(Nind/SUBPOP), error('Chrom and SUBPOP disagree'); end
   Nind = Nind/SUBPOP;  % Compute number of individuals per subpopulation

   if nargin < 3, RecOpt = 0.7; end
   if nargin > 2,
      if isempty(RecOpt), RecOpt = 0.7;
      elseif isnan(RecOpt), RecOpt = 0.7;
      elseif length(RecOpt) ~= 1, error('RecOpt must be a scalar');
      elseif (RecOpt < 0 | RecOpt > 1), error('RecOpt must be a scalar in [0, 1]'); end
   end

% Select individuals of one subpopulation and call low level function
   NewChrom = [];
   for irun = 1:SUBPOP,
      ChromSub = Chrom((irun-1)*Nind+1:irun*Nind,:);
      NewChromSub = feval(REC_F, ChromSub, RecOpt);
      NewChrom=[NewChrom; NewChromSub];
   end

% End of function

author: vonsylvia


text:
function [ output_maps ] = convolution( input_maps, kernels, biases )
%CONVOLUTION Summary of this function goes here
%   Detailed explanation goes here

    % If MatConvNet is not installed use Matlab (much slower)
    if(exist('vl_nnconv', 'file') == 3)
        output_maps = vl_nnconv(single(input_maps), kernels, biases);
    else
        n_filters = size(kernels, 4);

        kernels2 = kernels(:,:,end:-1:1,:);
        for i=1:n_filters
            for n_in_maps=1:size(kernels,3)
                kernels2(:,:,n_in_maps,i) = fliplr(squeeze(kernels2(:,:,n_in_maps,i)));
                kernels2(:,:,n_in_maps,i) = flipud(squeeze(kernels2(:,:,n_in_maps,i)));
            end
        end
        output_maps = [];
        for i=1:n_filters
            output_maps = cat(3, output_maps, convn(input_maps, kernels2(:,:,:,i), 'valid') + biases(i));
        end
    end
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function point = intersectLinePlane(line, plane, varargin)
%INTERSECTLINEPLANE Intersection point between a 3D line and a plane
%
%   PT = intersectLinePlane(LINE, PLANE)
%   Returns the intersection point of the given line and the given plane.
%   LINE:  [x0 y0 z0 dx dy dz]
%   PLANE: [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2]
%   PT:    [xi yi zi]
%   If LINE and PLANE are parallel, return [NaN NaN NaN].
%   If LINE (or PLANE) is a matrix with 6 (or 9) columns and N rows, result
%   is an array of points with N rows and 3 columns.
%
%   PT = intersectLinePlane(LINE, PLANE, TOL)
%   Specifies the tolerance factor to test if a line is parallel to a
%   plane. Default is 1e-14.
%
%   Example
%     % define horizontal plane through origin
%     plane = [0 0 0   1 0 0   0 1 0];
%     % intersection with a vertical line
%     line = [2 3 4  0 0 1];
%     intersectLinePlane(line, plane)
%     ans =
%        2   3   0
%     % intersection with a line "parallel" to plane
%     line = [2 3 4  1 2 0];
%     intersectLinePlane(line, plane)
%     ans =
%       NaN  NaN  NaN
%
%   See also:
%   lines3d, planes3d, points3d, clipLine3d
%

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 17/02/2005.
%

%   HISTORY
%   24/11/2005 add support for multiple input
%   23/06/2006 correction from Songbai Ji allowing different number of
%       lines or plane if other input has one row
%   14/12/2006 correction for parallel lines and plane normals
%   05/01/2007 fixup for parallel lines and plane normals
%   24/04/2007 rename as 'intersectLinePlane'
%   11/19/2010 Added bsxfun functionality for improved speed (Sven Holcombe)
%   01/02/2011 code cleanup, add option for tolerance, update doc

% extract tolerance if needed
tol = 1e-14;
if nargin > 2
    tol = varargin{1};
end

% unify sizes of data
nLines  = size(line, 1);
nPlanes = size(plane, 1);

% N planes and M lines not allowed
if nLines ~= nPlanes && min(nLines, nPlanes) > 1
    error('MatGeom:geom3d:intersectLinePlane', ...
        'Input must have same number of rows, or one must be 1');
end

% plane normal
n = crossProduct3d(plane(:,4:6), plane(:,7:9));

% difference between origins of plane and line
dp = bsxfun(@minus, plane(:, 1:3), line(:, 1:3));

% dot product of line direction with plane normal
denom = sum(bsxfun(@times, n, line(:,4:6)), 2);

% relative position of intersection point on line (can be inf in case of a
% line parallel to the plane)
t = sum(bsxfun(@times, n, dp),2) ./ denom;

% compute coord of intersection point
point = bsxfun(@plus, line(:,1:3),  bsxfun(@times, [t t t], line(:,4:6)));

% set indices of line and plane which are parallel to NaN
par = abs(denom) < tol;
point(par,:) = NaN;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function gridutil_shiftAxesUp(up, hc)
%GRIDUTIL_SHIFTAXESUP - Shifts up all axes in HC (default: all axes in the
% current figure) by a factor of UP (default: 10%)

if ~exist('up','var') || isempty(up), up=0.1; end
if ~exist('hc','var'), hc= get(gcf, 'children'); end

for hi= 1:length(hc);
  tag= get(hc(hi), 'tag');
  if isempty(tag) || strcmp(tag, 'Colorbar'),
    pos= get(hc(hi), 'position');
    if length(pos)==4,  %% otherwise it might be a uicontextmenu
      pos(2)= 1 - (1-pos(2))*(1-up);
      pos(4)= pos(4)*(1-up);
      set(hc(hi), 'position', pos);
    end
  end
end

author: bbci


text:
function[varargout] = vmedian(varargin)
%VMEDIAN  Median over finite elements along a specified dimension.
%
%   Y=VMEDIAN(X,DIM) takes the median of all finite elements of X along
%   dimension DIM.
%
%   [Y1,Y2,...YN]=VMEDIAN(X1,X2,...XN,DIM) also works.
%
%   VMEDIAN(X1,X2,...XN,DIM); with no output arguments overwrites the
%   original input variables.
%
%   VMEDIAN, like MATLAB's MEDIAN, defines the median over an even
%   number of values to be the average of the two middle elements.
%
%   VMEDIAN uses a fast algorithm which can be several times faster
%   than MEDIAN, but unlike MEDIAN, excludes both INFs and NANs.
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2008--2015 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
  vmedian_test,return
end

dim=varargin{end};

for i=1:length(varargin)-1
   if isreal(varargin{i})
       varargout{i}=vmedian1(varargin{i},dim);
   else
       varargout{i}=vmedian1(real(varargin{i}),dim)+sqrt(-1)*vmedian1(imag(varargin{i}),dim);
   end
end

eval(to_overwrite(nargin-1))

function[med]=vmedian1(data,dim)

numel=sum(isfinite(data),dim);
med=vzeros(size(numel),'nan');

temp=vswap(data,nan,inf);
temp=vswap(temp,-inf,inf);
sorted=sort(temp,dim,'ascend');
sorted=permute(sorted,[1:dim-1 dim+1:ndims(sorted) dim]);
sorted=reshape(sorted,[length(numel(:)),size(sorted,ndims(sorted))]);

ii=(1:length(numel(:)))';
ii=reshape(ii,size(numel));

boolodd=isodd(numel);
booleven=~boolodd&(numel~=0);

%ii(boolodd)

%(numel(boolodd)+1)./2

%looking to fix bug reported by Sandra
%indexodd=sub2ind(size(sorted),ii(boolodd),(numel(boolodd)+1)./2);
%indexeven1=sub2ind(size(sorted),ii(booleven),numel(booleven)./2);
%indexeven2=sub2ind(size(sorted),ii(booleven),numel(booleven)./2+1);

% vsize(sorted,ii,find(boolodd))
% anyany(~isfinite(boolodd))
% minmin(numel(boolodd))
% maxmax(numel(boolodd))
% median(vcolon(numel(boolodd)),1)
% minmin((numel(boolodd)+1)/2)
% maxmax((numel(boolodd)+1)/2)
% median((numel(boolodd)+1)/2,1)

if ~isempty(boolodd)
    indexodd=sub2ind(size(sorted),ii(boolodd),(numel(boolodd)+1)./2);
    med(boolodd)=sorted(indexodd);
end
if ~isempty(booleven)
    indexeven1=sub2ind(size(sorted),ii(booleven),numel(booleven)./2);
    indexeven2=sub2ind(size(sorted),ii(booleven),numel(booleven)./2+1);
    med(booleven)=sorted(indexeven1)./2+sorted(indexeven2)./2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[]=vmedian_test

N=6;
x1=randn(N,N,N,N,N);
bool=false(5,1);
tol=1e-4;

etime1=0;
etime2=0;
for i=1:5
    tic;y1=vmedian(x1,i);etime1=etime1+toc;
    tic;z1=median(x1,i);etime2=etime2+toc;
    bool(i)=aresame(y1,z1,tol);
end
reporttest('VMEDIAN 5-D with no NANs versus MEDIAN', allall(bool))
disp(['VMEDIAN was ' num2str(etime2./etime1) ' faster than MEDIAN.'])

N=30;
x1=randn(N,N,N);
x1(2:7:end)=nan;

tic
y1=vmedian(x1,1);
y2=vmedian(x1,2);
y3=vmedian(x1,3);
etime1=toc;

z1=zeros(size(y1));
z2=zeros(size(y2));
z3=zeros(size(y3));

tic
for i=1:N
    for j=1:N
        z1(1,i,j)=median(x1(isfinite(x1(:,i,j)),i,j),1);
        z2(i,1,j)=median(x1(i,isfinite(x1(i,:,j)),j),2);
        z3(i,j)=median(x1(i,j,isfinite(x1(i,j,:))),3);
    end
end
etime2=toc;
tol=1e-4;
reporttest('VMEDIAN 3-D including NANs', aresame(y1,z1,tol) && aresame(y2,z2,tol)&&aresame(y3,z3,tol))
disp(['VMEDIAN was ' num2str(etime2./etime1) ' faster than loop with MEDIAN.'])

N=30;
x1=round(randn(N,N,N));
x1(2:7:end)=nan;
x1(1:9:end)=x1(1);

tic
y1=vmedian(x1,1);
y2=vmedian(x1,2);
y3=vmedian(x1,3);
etime1=toc;

z1=zeros(size(y1));
z2=zeros(size(y2));
z3=zeros(size(y3));

tic
for i=1:N
    for j=1:N
        z1(1,i,j)=median(x1(isfinite(x1(:,i,j)),i,j),1);
        z2(i,1,j)=median(x1(i,isfinite(x1(i,:,j)),j),2);
        z3(i,j)=median(x1(i,j,isfinite(x1(i,j,:))),3);
    end
end
etime2=toc;
tol=1e-4;
reporttest('VMEDIAN 3-D including NANs and repeated data', aresame(y1,z1,tol) && aresame(y2,z2,tol)&&aresame(y3,z3,tol))
disp(['VMEDIAN was ' num2str(etime2./etime1) ' faster than loop with MEDIAN.'])

reporttest('VMEDIAN NaNs with row median',aresame([1 nan],vmedian([1 nan; nan nan],1)));
reporttest('VMEDIAN NaNs with column median',aresame([1 nan]',vmedian([1 nan; nan nan],2)));

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function varargout=ceil(varargin)
%CEIL (overloaded)

switch class(varargin{1})

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.

        x = varargin{1};

        dim = size(x);
        x = reshape(x,prod(dim),1);
        y = [];
        for i = 1:prod(dim)
            y = [y;yalmip('addextendedvariable',mfilename,extsubsref(x,i))];
        end
        y = reshape(y,dim);
        varargout{1} = y;

    case 'char' % YALMIP send 'graph' when it wants the epigraph or hypograph
        switch varargin{1}
            case {'milp','graph'}
                % Description using epigraphs
                t = varargin{2};
                X = varargin{3};

                c = intvar(1,1);
                F = (x <= c <= x + 1);

                varargout{1} = F;
                varargout{2} = struct('convexity','milp','monotonicity','milp','definiteness','milp');
                varargout{3} = X;

            otherwise
                error('SDPVAR/SORT called with CHAR argument?');
        end
    otherwise
        error('Strange type on first argument in SDPVAR/SORT');
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
%% Analyzing Neural Time Series Data
% Matlab code for Chapter 19
% Mike X Cohen
%
% This code accompanies the book, titled "Analyzing Neural Time Series Data"
% (MIT Press). Using the code without following the book may lead to confusion,
% incorrect data analyses, and misinterpretations of results.
% Mike X Cohen assumes no responsibility for inappropriate or incorrect use of this code.

%% figure 19.1

% define angles
a = [0.2 2*pi-.2];

figure
% plot unit vectors defined by those angles
polar([0 0; a],[0 0; 1 1])
hold on
% plot a unit vector with the average angle
polar([0 mean(a)],[0 1],'r')
% plot the average vector
polar([0 angle(mean(exp(1i*a)))],[0 abs(mean(exp(1i*a)))],'m')

%% figure 19.2

% load sample scalp EEG data
load sampleEEGdata

% center frequency
centerfreq = 12; % in Hz
chan2plot  = 'pz';
times2plot = [ 200 800 ]; % in ms, from stimulus onset

% definte convolution parameters
n_wavelet     = EEG.pnts;
n_data        = EEG.pnts*EEG.trials;
n_convolution = n_wavelet+n_data-1;
n_conv_pow2   = pow2(nextpow2(n_convolution));

% create wavelet. Note that the time vector to create the wavelet uses the
% same number of points as there are in the EEG data. If the EEG data has
% an even number of points (which is the case here--640), the wavelet will
% not have an exact center point. Though not technically incorrect, using
% a wavelet with an even number of points should be avoided when possible.
time    = -EEG.pnts/EEG.srate/2:1/EEG.srate:EEG.pnts/EEG.srate/2-1/EEG.srate;
wavelet = exp(2*1i*pi*centerfreq.*time) .* exp(-time.^2./(2*((4/(2*pi*centerfreq))^2)))/centerfreq;

% get FFT of data
eegfft = fft(reshape(EEG.data(strcmpi(chan2plot,{EEG.chanlocs.labels}),:,:),1,[]),n_conv_pow2);

% convolution
eegconv = ifft(fft(wavelet,n_conv_pow2).*eegfft);
eegconv = eegconv(1:n_convolution);
% Cut the edges off the result of convolution. Because the wavelet here was
% created with the same number of time points as the EEG data, EEG.pnts
% also corresponds to the number of time points in the wavelet.
eegconv = reshape(eegconv(floor((EEG.pnts-1)/2):end-1-ceil((EEG.pnts-1)/2)),EEG.pnts,EEG.trials);

% plot
figure
for subploti=1:2
    subplot(2,2,subploti)
    [junk,idx]=min(abs(EEG.times-times2plot(subploti)));
    polar([zeros(1,EEG.trials); angle(eegconv(idx,:))],[zeros(1,EEG.trials); ones(1,EEG.trials)])
    title([ 'ITPC at ' num2str(times2plot(subploti)) ' ms = ' num2str(round(1000*abs(mean(exp(1i*angle(eegconv(idx,:))))))/1000) ])

    subplot(2,2,subploti+2)
    hist(angle(eegconv(idx,:)),20)
    set(gca,'xlim',[-pi pi]*1.1,'xtick',-pi:pi/2:pi,'xticklabel',{'-pi';'-pi/2';'0';'pi';'pi/2'},'ylim',[0 20])
end

%% Figure 19.3

vectors{1} = [0 pi/3];
vectors{2} = [0 pi/2];
vectors{3} = [0 2*pi/3];
vectors{4} = [0 pi*.9];

figure
for i=1:length(vectors)
    subplot(1,length(vectors),i)

    % plot individual unit vectors
    polar([0 vectors{i}(1)],[0 1],'k')
    hold on
    polar([0 vectors{i}(2)],[0 1],'k')

    % plot mean vector (ignore the math for now,
    % this will be discussed later in the chapter)
    meanvect = mean(exp(1i*vectors{i}));
    polar([0 angle(meanvect)],[0 abs(meanvect)],'r')
    title([ 'Mean vector length: ' num2str(abs(meanvect)) ])
end

%% Figure 19.4

% get FFT of data
chan2plot  = 'pz';
centerfreq = 12; % in Hz
eegfft = fft(reshape(EEG.data(strcmpi(chan2plot,{EEG.chanlocs.labels}),:,:),1,[]),n_conv_pow2);

% ITPC at one frequency band
wavelet = exp(2*1i*pi*centerfreq.*time) .* exp(-time.^2./(2*((4/(2*pi*centerfreq))^2)))/centerfreq;
% convolution
eegconv = ifft(fft(wavelet,n_conv_pow2).*eegfft);
eegconv = eegconv(1:n_convolution);
eegconv = reshape(eegconv(floor((EEG.pnts-1)/2):end-1-ceil((EEG.pnts-1)/2)),EEG.pnts,EEG.trials);

figure
plot(EEG.times,abs(mean(exp(1i*angle(eegconv)),2)))
set(gca,'xlim',[-200 1000])
xlabel('Time (ms)')
ylabel('ITPC')

% TF plot of ITPC
frequencies = logspace(log10(4),log10(40),20);
s = logspace(log10(3),log10(10),length(frequencies))./(2*pi*frequencies);
itpc = zeros(length(frequencies),EEG.pnts);

for fi=1:length(frequencies)
    % create wavelet
    wavelet = exp(2*1i*pi*frequencies(fi).*time) .* exp(-time.^2./(2*(s(fi)^2)))/frequencies(fi);

    % convolution
    eegconv = ifft(fft(wavelet,n_conv_pow2).*eegfft);
    eegconv = eegconv(1:n_convolution);
    eegconv = reshape(eegconv(floor((EEG.pnts-1)/2):end-1-ceil((EEG.pnts-1)/2)),EEG.pnts,EEG.trials);

    % extract ITPC
    itpc(fi,:) = abs(mean(exp(1i*angle(eegconv)),2));
end

figure
contourf(EEG.times,frequencies,itpc,40,'linecolor','none')
set(gca,'clim',[0 .6],'xlim',[-200 1000])
xlabel('Time (ms)'), ylabel('Frequencies (Hz)')

%% figure 19.5

n_trials    = 500;
itpcByNFake = zeros(1,n_trials);

for n=1:n_trials
    for iteri=1:50
        itpcByNFake(n) = itpcByNFake(n) + abs(mean(exp(1i*(rand(1,n)*2*pi-pi))));
    end
end
itpcByNFake = itpcByNFake./iteri;

% Z and P (you will learn more about these formulae in chapter 34)
itpcByNFakeZ = (1:n_trials).*(itpcByNFake.^2);
itpcByNFakeP = exp(sqrt(1+4.*(1:n_trials)+4*( ((1:n_trials).^2) - ((1:n_trials).*itpcByNFake).^2))-(1+2*(1:n_trials))); % used later
itpcFakeCrit = sqrt( -log(.01)./(1:n_trials) );

figure
plot(1:n_trials,itpcByNFake)
hold on
plot(1:n_trials,itpcFakeCrit,'r')
set(gca,'ylim',[0 1])
legend({'ITPC_r_a_w';'ITPC_C_r_i_t'})

%% figure 19.6

% center frequency
centerfreq = 6; % Hz
chan2plot  = 'fcz';

% definte convolution parameters
n_wavelet     = EEG.pnts;
n_data        = EEG.pnts*EEG.trials;
n_convolution = n_wavelet+n_data-1;
n_conv_pow2   = pow2(nextpow2(n_convolution));

% create wavelet
time    = -EEG.pnts/EEG.srate/2:1/EEG.srate:EEG.pnts/EEG.srate/2-1/EEG.srate;
wavelet = exp(2*1i*pi*centerfreq.*time) .* exp(-time.^2./(2*((4/(2*pi*centerfreq))^2)))/centerfreq;

% get FFT of data
eegfft  = fft(reshape(EEG.data(strcmpi(chan2plot,{EEG.chanlocs.labels}),:,:),1,[]),n_conv_pow2);

% convolution
eegconv = ifft(fft(wavelet,n_conv_pow2).*eegfft);
eegconv = eegconv(1:n_convolution);
eegconv = reshape(eegconv(floor((EEG.pnts-1)/2):end-1-ceil((EEG.pnts-1)/2)),EEG.pnts,EEG.trials);

figure

% compute ITPC as function of N
itpcByN = zeros(1,EEG.trials);
for n=1:EEG.trials
    % multiple iterations to select different sets of random trials
    for iteri=1:50
        trials2use  = randsample(EEG.trials,n);
        itpcByN(n)  = itpcByN(n) + mean(abs(mean(exp(1i*angle(eegconv(282:372,trials2use))),2)),1);
    end
end
itpcByN = itpcByN./iteri;

% Z and P
itpcByNZ = (1:EEG.trials).*(itpcByN.^2);
itpcByNP = exp(sqrt(1+4.*(1:EEG.trials)+4*( ((1:EEG.trials).^2) - ((1:EEG.trials).*itpcByN).^2))-(1+2*(1:EEG.trials)));
itpcCrit = sqrt( -log(.01)./(1:EEG.trials) ); % .01 is the p-value in this case

subplot(211)
plot(itpcByN)
hold on
plot(itpcCrit,'r')
xlabel('Number of trials in analysis')
ylabel('ITPC')
set(gca,'ylim',[.2 1])

% (Note that this figure will look different than that in the book
% because trials are randomly selected.)
subplot(212)
plot(EEG.times,abs(mean(exp(1i*angle(eegconv)),2)))
hold on
randomtrials2plot = randperm(EEG.trials); % could also use randsamp if you have the stats toolbox
plot(EEG.times,abs(mean(exp(1i*angle(eegconv(:,randomtrials2plot(1:20)))),2)),':')
xlabel('Time (ms)'), ylabel('ITPC')
legend({'99 trials';'20 trials'})

%% Figure 34.8

% This cell concerns statistical evaluation of ITPC values. It will be
% discussed in depth in chapter 34, but the code is presented here because
% it relies on calculations from the previous two figures.

% p-values under assumption of von Mises distribution
approx_pval_fake = exp(-itpcByNFakeZ);
approx_pval_real = exp(-itpcByNZ);

ncutoff = 10;

figure

subplot(121)
plot(itpcByNFakeP(ncutoff+1:end),approx_pval_fake(ncutoff+1:end),'mo','markersize',8)
hold on
plot(itpcByNFakeP(1:ncutoff),approx_pval_fake(1:ncutoff),'k.')
plot([0 1],[0 1],'k')
set(gca,'xlim',[0 1],'ylim',[0 1],'xtick',0:.25:1,'ytick',0:.25:1)
axis square
ylabel('approximate p')
xlabel('exact p')
title('P-values from fake data')

subplot(122)
plot(itpcByNP(ncutoff+1:end),approx_pval_real(ncutoff+1:end),'mo','markersize',8)
hold on
plot(itpcByNP(1:ncutoff),approx_pval_real(1:ncutoff),'k.')
plot([0.00001 1],[0.00001 1],'k')
set(gca,'xlim',[.0001 1],'ylim',[.0001 1],'xscale','log','yscale','log')
axis square
ylabel('approximate p')
xlabel('exact p')
title('P-values from real data')

%% figure 19.7
% (This figure takes a while to generate. You could also reduce the number
% of iterations to speed it up.)

% Data for this cell are from figure 19.6. If you want to plot
% results from a different channel, change the electrode in 19.6.

% center frequencies
frequencies = 1:40;

niterations = 50;

% initialize
itpcByNandF = zeros(length(frequencies),EEG.trials);

for fi=1:length(frequencies)

    centerfreq = frequencies(fi);

    wavelet = exp(2*1i*pi*centerfreq.*time) .* exp(-time.^2./(2*((4/(2*pi*centerfreq))^2)))/centerfreq;

    % convolution
    eegconv = ifft(fft(wavelet,n_conv_pow2).*eegfft);
    eegconv = eegconv(1:n_convolution);
    eegconv = reshape(eegconv(floor((EEG.pnts-1)/2):end-1-ceil((EEG.pnts-1)/2)),EEG.pnts,EEG.trials);

    for n=1:EEG.trials
        % multiple iterations to select different random sets of trials
        for iteri=1:niterations
            trials2use  = randsample(EEG.trials,n);
            itpcByNandF(fi,n)  = itpcByNandF(fi,n) + mean(abs(mean(exp(1i*angle(eegconv(282:372,trials2use))),2)),1);
        end
    end
end

figure
contourf(1:EEG.trials,frequencies,itpcByNandF./iteri,40,'linecolor','none')
colormap gray
set(gca,'clim',[.1 .5])
xlabel('Trials')
ylabel('Frequency (Hz)')

%% Figure 19.8

figure, set(gcf,'name','Rayleigh''s Z')
subplot(121)
plot(itpcByN)
hold on
plot(itpcByNFake(1:99),'r')
xlabel('Trial count'), ylabel('ITPC')
axis square

subplot(122)
plot(itpcByNZ)
hold on
plot(itpcByNFakeZ(1:99),'r')
xlabel('Trial count'), ylabel('ITPC_Z')
legend({'real data';'fake data'})
axis square

%% Figure 19.9

% pick electrode and frequencies
chan2plot   = 'fcz';
frequencies = 1:40;

% definte convolution parameters
time          = -EEG.pnts/EEG.srate/2:1/EEG.srate:EEG.pnts/EEG.srate/2-1/EEG.srate;
n_wavelet     = EEG.pnts;
n_data        = EEG.pnts*EEG.trials;
n_convolution = n_wavelet+n_data-1;
n_conv_pow2   = pow2(nextpow2(n_convolution));

plotlegends={'without jitter';'with jitter'};

figure

baselinetime = [ -300 -100 ];
baseidx=dsearchn(EEG.times',baselinetime(1)):dsearchn(EEG.times',baselinetime(2));

for simuli=1:2

    % add time jitter (or not)
    tempdat = squeeze(EEG.data(strcmpi(chan2plot,{EEG.chanlocs.labels}),:,:));
    for ti=1:size(tempdat,2)
        timejitter = ceil(rand*10)*(simuli-1); % (when simuli==1, timejitter==0)
        tempdat(:,ti) = tempdat([ timejitter+1:end 1:timejitter ],ti);
    end

    % get FFT of data
    eegfft = fft(reshape(tempdat,1,[]),n_conv_pow2);

    % initialize
    itpc = zeros(length(frequencies),EEG.pnts);
    powr = zeros(length(frequencies),EEG.pnts);

    for fi=1:length(frequencies)

        centerfreq = frequencies(fi);

        wavelet = exp(2*1i*pi*centerfreq.*time) .* exp(-time.^2./(2*((4/(2*pi*centerfreq))^2)))/centerfreq;

        % convolution
        eegconv = ifft(fft(wavelet,n_conv_pow2).*eegfft);
        eegconv = eegconv(1:n_convolution);
        eegconv = reshape(eegconv(floor((EEG.pnts-1)/2):end-1-ceil((EEG.pnts-1)/2)),EEG.pnts,EEG.trials);

        % compute and store itpc and power
        itpc(fi,:) = abs(mean(exp(1i*angle(eegconv)),2));
        powr(fi,:) = mean(abs(eegconv).^2,2);
        powr(fi,:) = 10*log10(powr(fi,:)./mean(powr(fi,baseidx),2));
    end

    subplot(2,2,simuli)
    contourf(EEG.times,frequencies,itpc,40,'linecolor','none')
    set(gca,'clim',[.1 .5],'xlim',[-200 1000])
    xlabel('Time (ms)')
    ylabel('Frequency (Hz)')
    title([ 'ITPC ' plotlegends{simuli} ])

    subplot(2,2,simuli+2)
    contourf(EEG.times,frequencies,powr,40,'linecolor','none')
    set(gca,'clim',[-3 3],'xlim',[-200 1000])
    xlabel('Time (ms)')
    ylabel('Frequency (Hz)')
    title([ 'DB-power ' plotlegends{simuli} ])
end

%% Figure 19.10

% initialize
data4test  = zeros(2,length(time),EEG.trials);
sensor2use = 'p7';

% amplitude modulation (modulate power by 1 Hz sine wave)
time    = -EEG.pnts/EEG.srate/2:1/EEG.srate:EEG.pnts/EEG.srate/2-1/EEG.srate;
amp_mod = (sin(2*pi*1.*time)+2)-1;

for triali=1:EEG.trials
    % each trial is a random channel and trial
    trialdata = EEG.data(strcmpi(sensor2use,{EEG.chanlocs.labels}),:,triali);

    % Uncomment the next line of code for band-pass filtered data.
    % This uses the eegfilt function, which is part of the eeglab toolbox.
    trialdata = eegfilt(double(trialdata),EEG.srate,10,20);

    data4test(1,:,triali) = trialdata.*amp_mod;
    data4test(2,:,triali) = trialdata;
end

% compute ITPC
itpc_mod   = abs(mean(exp(1i*angle(hilbert(data4test(1,:,:)))),3));
itpc_nomod = abs(mean(exp(1i*angle(hilbert(data4test(2,:,:)))),3));

% plot!
figure
subplot(311)
plot(EEG.times,amp_mod)
set(gca,'ylim',[-.2 2.2])
title('Amplitude modulator')

subplot(312)
plot(EEG.times,data4test(1,:,10))
hold on
plot(EEG.times,data4test(2,:,10),'r')
axis tight
title('Example trials')

subplot(313)
plot(EEG.times,itpc_mod)
hold on
plot(EEG.times,itpc_nomod,'r')
legend({'amplitude modulation';'no amp mod'})
set(gca,'ylim',[0 1])
xlabel('Time (ms)'), ylabel('ITPC')
title('ITPC')

figure
plot(amp_mod,itpc_mod,'.')
xlabel('Power modulation'), ylabel('ITPC')

%% Figure 19.11
% (This figure is populated with randomly generated data,
%  and so will look different from the book figure.)

randvects = rand(50,1)*2*pi-pi;
vectormod = (randvects + randn(size(randvects))).^2;
vectormod2 = vectormod-min(vectormod)+1; % make sure no negative values

figure

% ITPC
subplot(221)
polar([zeros(size(randvects)) randvects]',[zeros(size(randvects)) ones(size(randvects))]','k')
title([ 'ITPC = ' num2str(abs(mean(exp(1i*randvects)))) ])

% wITPC
subplot(222)
polar([zeros(size(randvects)) randvects]',[zeros(size(randvects)) vectormod]','k')

witpc = abs(mean(vectormod.*exp(1i*randvects)));
perm_witpc = zeros(1,1000);

for i=1:1000
    perm_witpc(i) = abs(mean(vectormod(randperm(length(vectormod))).*exp(1i*randvects)));
end

witpc_z = (witpc-mean(perm_witpc))./std(perm_witpc);
title([ '_wITPC_z = ' num2str(witpc_z) ])

% example of one permutation
subplot(223)
polar([zeros(size(randvects)) randvects]',[zeros(size(randvects)) vectormod(randperm(length(vectormod)))]','k')
title('One null hypothesis iteration')

% histogram of null-hypothesis WITPC
subplot(224)
[y,x]=hist(perm_witpc,50);
h=bar(x,y,'histc');
set(h,'linestyle','none')
hold on
plot([witpc witpc],get(gca,'ylim')/2,'m')
title('Histogram of null hypothesis wITPC')

%% Figure 19.12

centerfreq  = 6;
channel2use = 'po7';
times2save  = -200:50:1200;

% initialize matrix to store RTs
rts = zeros(size(EEG.epoch));

for ei=1:EEG.trials

    % find which event is time=0, and take the latency of the event thereafter.
    time0event = find(cell2mat(EEG.epoch(ei).eventlatency)==0);

    % use try-catch in case of no response
    try
        rts(ei) = EEG.epoch(ei).eventlatency{time0event+1};
    catch me;
        rts{ei} = NaN;
    end
end

% definte convolution parameters
time          = -EEG.pnts/EEG.srate/2:1/EEG.srate:EEG.pnts/EEG.srate/2-1/EEG.srate;
n_wavelet     = EEG.pnts;
n_data        = EEG.pnts*EEG.trials;
n_convolution = n_wavelet+n_data-1;
n_conv_pow2   = pow2(nextpow2(n_convolution));

% get FFT of data and wavelet
eegfft  = fft(reshape(EEG.data(strcmpi(channel2use,{EEG.chanlocs.labels}),:,:),1,n_data),n_conv_pow2);
wavefft = fft(exp(2*1i*pi*centerfreq.*time) .* exp(-time.^2./(2*((4/(2*pi*centerfreq))^2))),n_conv_pow2);

% convolution
eegconv = ifft(wavefft.*eegfft);
eegconv = eegconv(1:n_convolution);
eegconv = reshape(eegconv(floor((EEG.pnts-1)/2):end-1-ceil((EEG.pnts-1)/2)),EEG.pnts,EEG.trials);

phase_angles = angle(eegconv);

% initialize
itpc    = zeros(size(times2save));
witpc   = zeros(size(times2save));
witpc_z = zeros(size(times2save));

for ti=1:length(times2save)

    % find index for this time point
    [junk,timeidx] = min(abs(EEG.times-times2save(ti)));

    % ITPC is unmodulated phase clustering
    itpc(ti) = abs(mean(exp(1i*phase_angles(timeidx,:))));

    % wITPC is rts modulating the length of phase angles
    witpc(ti) = abs(mean(rts.*exp(1i*phase_angles(timeidx,:))));

    % permutation testing
    perm_witpc = zeros(1,1000);
    for i=1:1000
        perm_witpc(i) = abs(mean(rts(randperm(EEG.trials)).*exp(1i*phase_angles(timeidx,:))));
    end

    witpc_z(ti) = (witpc(ti)-mean(perm_witpc))./std(perm_witpc);
end

figure

% plot ITPC
subplot(311)
h=plotyy(EEG.times,abs(mean(exp(1i*phase_angles),2)),times2save,witpc);
set(h,'xlim',[times2save(1) times2save(end)]);
title([ 'ITPC and wITPC at ' channel2use ])
legend({'ITPC';'wITPC'})

% plot wITPCz
subplot(312)
plot(times2save,witpc_z)
hold on
plot(get(gca,'xlim'),[0 0],'k')
set(gca,'xlim',[times2save(1) times2save(end)])
xlabel('Time (ms)'), ylabel('wITPCz')
title([ 'wITPCz at ' channel2use ])

subplot(325)
plot(itpc,witpc,'.')
xlabel('ITPC'), ylabel('wITPC')
axis square

subplot(326)
plot(itpc,witpc_z,'.')
xlabel('ITPC'), ylabel('wITPCz')
axis square

%% end.

author: mikexcohen

repo: AnalyzingNeuralTimeSeries


text:
function [ha hb hc] = shadedplot(x, y1, y2, varargin)

% SHADEDPLOT draws two lines on a plot and shades the area between those
% lines.
%
% SHADEDPLOT(x, y1, y2)
%   All of the arguments are vectors of the same length, and each y-vector is
%   horizontal (i.e. size(y1) = [1  N]). Vector x contains the x-axis values,
%   and y1:y2 contain the y-axis values.
%
%   Plot y1 and y2 vs x, then shade the area between those two
%   lines. Highlight the edges of that band with lines.
%
%   SHADEDPLOT(x, y1, y2, areacolor, linecolor)
%   The arguments areacolor and linecolor allow the user to set the color
%   of the shaded area and the boundary lines. These arguments must be
%   either text values (see the help for the PLOT function) or a
%   3-element vector with the color values in RGB (see the help for
%   COLORMAP).
%
%   [HA HB HC = SHADEDPLOT(x, y1, y2) returns three handles to the calling
%   function. HA is a vector of handles to areaseries objects (HA(2) is the
%   shaded area), HB is the handle to the first line (x vs y1), and HC is
%   the handle to the second line (x vs y2).
%
%   Example:
%
%     x1 = [1 2 3 4 5 6];
%     y1 = x1;
%     y2 = x1+1;
%     x3 = [1.5 2 2.5 3 3.5 4];
%     y3 = 2*x3;
%     y4 = 4*ones(size(x3));
%     ha = shadedplot(x1, y1, y2, [1 0.7 0.7], 'r'); %first area is red
%     hold on
%     hb = shadedplot(x3, y3, y4, [0.7 0.7 1]); %second area is blue
%     hold off

% plot the shaded area
y = [y1; (y2-y1)]';
ha = area(x, y);
set(ha(1), 'FaceColor', 'none') % this makes the bottom area invisible
set(ha, 'LineStyle', 'none')

% plot the line edges
hold on
hb = plot(x, y1, 'LineWidth', 1);
hc = plot(x, y2, 'LineWidth', 1);
hold off

% set the line and area colors if they are specified
switch length(varargin)
        case 0
        case 1
                set(ha(2), 'FaceColor', varargin{1})
        case 2
                set(ha(2), 'FaceColor', varargin{1})
                set(hb, 'Color', varargin{2})
                set(hc, 'Color', varargin{2})
        otherwise
end

% put the grid on top of the colored area
set(gca, 'Layer', 'top')
grid off

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function imData=read_file_old(path_to_file,sframe,num2read)

% Reads uncompressed multipage .tiff, .hdf5, or .avi files
% Usage:  my_data=read_file('path_to_data_file, start frame, num to read);

% INPUTS:
% path_to_file:     location of file to be read
% sframe:           first frame to read (optional, default: 1)
% num2read:         number of frames to read (optional, default: read the whole file)

% OUTPUT:
% imData:           data in array format

% Written by Eftychios A. Pnevmatikakis, Simons Foundation

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

[~,~,ext] = fileparts(path_to_file);

if strcmpi(ext,'.tiff') || strcmpi(ext,'.tif');

    %get image info
    tiffInfo = imfinfo(path_to_file);
    T = length(tiffInfo);
    if nargin < 3 || isempty(num2read); num2read = T - sframe + 1; end
    num2read = min(num2read,T - sframe + 1);

    Y1 = imread(path_to_file,'Index',sframe,'Info',tiffInfo);
    imData = zeros([size(Y1),num2read],'like',Y1);
    nd = ndims(Y1);
    if nd == 2
        imData(:,:,1) = Y1;
        for t = sframe+1:sframe+num2read-1
            imData(:,:,t-sframe+1) = imread(path_to_file,'Index',t,'Info',tiffInfo);
        end
    elseif nd == 3
        imData(:,:,:,1) = Y1;
        for t = sframe+1:sframe+num2read-1
            imData(:,:,:,t-sframe+1) = imread(path_to_file,'Index',t,'Info',tiffInfo);
        end
    end

elseif strcmpi(ext,'.hdf5') || strcmpi(ext,'.h5');
    info = hdf5info(path_to_file);
    dims = info.GroupHierarchy.Datasets.Dims;
    name = info.GroupHierarchy.Datasets.Name;
    if nargin < 3
        num2read = dims(end)-sframe+1;
    end
    num2read = min(num2read,dims(end)-sframe+1);
    imData = h5read(path_to_file,name,[ones(1,length(dims)-1),sframe],[dims(1:end-1),num2read]);
elseif strcmpi(ext,'.avi')
    v = VideoReader(path_to_file);
    if nargin < 3
        num2read = v.Duration*v.FrameRate-sframe+1;
    end
    Y1 = readFrame(v);
    imData = zeros(v.Height,v.Width,num2read,'like',Y1);
    i = 0;
    while hasFrame(v)
        video = readFrame(v);
        i = i + 1;
        if i >= sframe
            imData(:,:,i-sframe+1) = video;
        end
        if i - sframe + 1 >= num2read
            break;
        end
    end
else
    error('Unknown file extension. Only .tiff, .avi and .hdf5 files are currently supported');
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [nrows ncols] = subplotSz(n)
% [nrows ncols] = subplotSz(n)
%
% Choose the size, rows x cols, of a subplot. Make it as small as possible
% and as close to square as possible.
%
% Example:
%
% [nrows, ncols] = subplotSz(19)

if nargin < 1,
    help subplotSz
    return
end

nrows = round(sqrt(n));
ncols = ceil(n/nrows);

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = acsch(varargin)

switch class(varargin{1})

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        operator = CreateBasicOperator('callback');
        operator.monotonicity = @(xL,xU)decreasing_except_at(xL,xU,0);
        operator.derivative = @(x)(-1./(abs(x).*sqrt(1 + x.^2)));
        operator.singularity = [0 -inf inf];
        operator.inflection = [-inf -1 0 1];

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
clear all
clf

projV=@(x,A)A*inv(A'*A)*A'*x;
V1=[1 0 0]';
V2=[1 1 0]';

% Alternating projections.
x=[0.5 1 1]';
subplot(121)
plot3([0 V1(1)],[0 V1(2)],[0 V1(3)],'-b','LineWidth',2);hold on
plot3([0 V2(1)],[0 V2(2)],[0 V2(3)],'-b','LineWidth',2);
plot3(x(1),x(2),x(3),'.k','MarkerSize',20);
for i=1:5
  xV1=projV(x,V1);
  h=arrow(x,xV1,'Length',5);
  set(h,'LineStyle',':','EdgeColor','r','FaceColor','r');
  text(xV1(1),xV1(2),xV1(3),['P_{V_1}(x)'])
  x=projV(xV1,V2);
  h=arrow(xV1,x,'Length',5);
  text(x(1),x(2),x(3),['P_{V_2}(P_{V_1}(x))'])
  plot3(x(1),x(2),x(3),'.k','MarkerSize',20);
end
axis([0 1 0 1 0 1]);box on;axis equal;axis tight
title('Alternating projections');

% Average Alternating Reflections/DR.
x=[0.5 1 1]';
subplot(122)
plot3([-1 V1(1)],[0 V1(2)],[0 V1(3)],'-b','LineWidth',2);hold on
plot3([-1 V2(1)],[-1 V2(2)],[0 V2(3)],'-b','LineWidth',2);
plot3(x(1),x(2),x(3),'.k','MarkerSize',20);
for i=1:5
  xV1=2*projV(x,V1)-x;
  h=arrow(x,xV1,'Length',5);
  set(h,'LineStyle',':','EdgeColor','r','FaceColor','r');
  text(xV1(1),xV1(2),xV1(3),['R_{V_1}(x)'])
  xV2=2*projV(xV1,V2)-xV1;
  h=arrow(xV1,xV2,'Length',5);
  set(h,'LineStyle',':','EdgeColor','r','FaceColor','r');
  text(xV2(1),xV2(2),xV2(3),['R_{V_2}(R_{V_1}(x))'])
  x=(xV2+x)/2;
  arrow(xV2,x,'Length',5);
  text(x(1),x(2),x(3),['(R_{V_2}(R_{V_1}(x))+x)/2'])
  plot3(x(1),x(2),x(3),'.k','MarkerSize',20);
end
axis([-1 1 -1 1 -1 1]);box on;axis equal;axis tight
title('Average Alternating Reflections/DR');

author: gpeyre

repo: matlab-toolboxes


text:
function outcome=getOutcomeStr(obj, verbosity)
% give string indicating the outcome of the test
%
% outcome=getOutcomeStr(obj, verbosity)
%
% Inputs:
%   obj                 MOxUnitTestOutcome object
%   verbosity           Integer in the range [0,1,2]
%
% Returns:
%   outcome:            Depending on the value of verbosity
%                       0: ''
%                       1: '.'
%                       2: 'passed'
%

    outcome_cell={'','.','passed'};
    outcome=outcome_cell{verbosity+1};

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [C, Q] = procCosp(DATA, WINDOW_TYPE, WINDOW_SIZE, OVERLAP, FREQUENCY_RANGE, PHASE_CORRECTION)
% [C, Q] = procCosp(DATA, WINDOW_TYPE, WINDOW_SIZE, OVERLAP, FREQUENCY_RANGE)
%
% This function computes cospectra and quadrature spectra on multivariate
% data using Welch's method. For a better estimation, each tappering window
% can be phase corrected according to its length and its position in data.
%
% Inputs:
% - DATA            --> matrix with the data (M samples by N channels matrix)
% - WINDOW_TYPE     --> (optionnal) string with the type of the window
% - WINDOW_SIZE     --> (optionnal) size of the window in number of samples
% - OVERLAP         --> (optionnal) percentage of the overlapping window (from 0 to 0.99)
% - FREQUENCY_RANGE --> (optionnal) range of interesting frequencies to keep,
%                                   format must be [sample_rate min_frec max_frec]
% - PHASE_CORRECTION--> (optionnal) correction of phase shift intrinsic to Welch's windows.
%                                   (1: correction, 0: no correction)
%
% Outputs:
% - C        --> cospectral matrix [N x N x f]
% - Q        --> quadrature matrix [N x N x f]
%
% History
% First version: Jonas Chatel-Goldman @ GIPSA-Lab, 01/12/2011

% default UI display
if(nargin < 9)  UIdisplay  = 0;  end
% default correction
if(nargin < 6)  PHASE_CORRECTION = 0; end
% default frequency range
if(nargin < 5)  FREQUENCY_RANGE = []; end
% default overlapping
if(nargin < 4)  OVERLAP = .75;  end
% default windows size
if(nargin < 3)  WINDOW_SIZE = 128;  end
% default windows type
if(nargin < 2)  WINDOW_TYPE = 'hanning';  end

% Adjust the window to the next pow of 2 (for the FFT)
WINDOW_SIZE = 2 ^ nextpow2(WINDOW_SIZE);
% Calculate the number of frequencies
number_freqs = (WINDOW_SIZE / 2); % +1;

s_data = size(DATA);

% UI DISPLAY
if(UIdisplay)  h= waitbar(0,'Processing cospectrum...');  end

%% Create the window from the inputs and preallocates memory
switch lower(WINDOW_TYPE)
    case 'hamming'
        win = hamming(WINDOW_SIZE);
    case 'hann'
        win = hann(WINDOW_SIZE);
    case 'hanning'
        win = hanning(WINDOW_SIZE);
    case 'blackman'
        win = blackman(WINDOW_SIZE);
    case 'barthannwin'
        win = barthannwin(WINDOW_SIZE);
    case 'blackmanharris'
        win = blackmanharris(WINDOW_SIZE);
    case 'bohmanwin'
        win = bohmanwin(WINDOW_SIZE);
    case 'chebwin'
        win = chebwin(WINDOW_SIZE);
    case 'gausswin'
        win = gausswin(WINDOW_SIZE);
    case 'kaiser'
        win = kaiser(WINDOW_SIZE);
    case 'nuttallwin'
        win = nuttallwin(WINDOW_SIZE);
    case 'parzenwin'
        win = parzenwin(WINDOW_SIZE);
    case 'tukeywin'
        win = tukeywin(WINDOW_SIZE);
    case 'rectangular'
        win = ones(WINDOW_SIZE,1);
    case 'flattopwin'
        win = flattopwin(WINDOW_SIZE);
    case 'welch'
         %win = welchwin(WINDOW_SIZE, 0);
        Nd2=(WINDOW_SIZE-1)/2;
        for i=1 : WINDOW_SIZE
            win(i)=1-((i-Nd2)/Nd2)^2;
        end
        win = win';
    otherwise
        disp('WARNING - Unknown window type!');
        disp('Switching to default window type : Hamming Window');
        win = hamming(WINDOW_SIZE);
end

if s_data(1) < WINDOW_SIZE
    error('Cospectra computation: data segment too short for this window length, consider changing frequency resolution');
end

% Calculate the number of windows because of the overlapping
if(OVERLAP == 0)
    number_windows = floor(s_data(1) / WINDOW_SIZE);
else
    nbFullWin = floor(s_data(1)/WINDOW_SIZE);
    number_windows = 1 + (nbFullWin-1)/(1-OVERLAP) + floor((s_data(1)-((nbFullWin)*WINDOW_SIZE))/((1-OVERLAP)*WINDOW_SIZE));
end

% pre-allocation of memory
S = zeros(s_data(2), s_data(2), number_freqs);

%% Loop on all frequencies
for window_ix = 1 : number_windows

    % UI display
    if(UIdisplay)  waitbar(window_ix/number_windows,h); end

    % time markers to select the data
    t1 = floor((window_ix-1) * (1-OVERLAP) * WINDOW_SIZE) +1;   % marker of the beginning of the time window
    t2 = t1 + WINDOW_SIZE -1;                           % marker of the end of the time window

    % select current window and apodize it
    cdata = DATA(t1:t2, :) .* (win*ones(1,s_data(2)));

    % FFT calculation
    fdata = fft(cdata,WINDOW_SIZE) ./ number_windows; % complex data
    if(PHASE_CORRECTION == 1)
        fdata = fdata.*(exp(-sqrt(-1)*t1*( 0:(WINDOW_SIZE-1) )'/WINDOW_SIZE*2*pi)*ones(1,s_data(2)));
    end

    % Complexe cospectrum averaging over windows
    for f = 1 : number_freqs
        S(:,:,f) = S(:,:,f) + (fdata(f,:)' * fdata(f,:));
    end
end

C = real(S);
Q = imag(S);

% Adjust Frequency range to specified range (in case it is a parameter)
if ~isempty(FREQUENCY_RANGE)
    if  FREQUENCY_RANGE(3) <= (FREQUENCY_RANGE(1) / 2)
        Faxis = [0:(1/number_freqs):1-(1/number_freqs)];
        FindexMin = find(Faxis >= FREQUENCY_RANGE(2)/(FREQUENCY_RANGE(1)/2), 1);
        FindexMax = find(Faxis >= FREQUENCY_RANGE(3)/(FREQUENCY_RANGE(1)/2), 1);
        C = C(:, :, FindexMin:FindexMax);
        Q = Q(:, :, FindexMin:FindexMax);

    else
        % Assertion on frequency, to comply with Nyquist Theorem
        error('Frequency range exceeds data sampling rate in consideration of Shannon''s limitation.')
    end
end

% UI display
if(UIdisplay) close(h); end

author: alexandrebarachant

repo: covariancetoolbox


text:
classdef BankAngleTermCondOptimizationVariable < AbstractOptimizationVariable
    %BankAngleTermCondOptimizationVariable Summary of this class goes here
    %   Detailed explanation goes here

    properties
        varObj(1,1) BankAngleTermCondition = BankAngleTermCondition(0);

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        useTf(1,1) = false;
    end

    methods
        function obj = BankAngleTermCondOptimizationVariable(varObj)
            obj.varObj = varObj;
            obj.varObj.optVar = obj;

            obj.id = rand();
        end

        function x = getXsForVariable(obj)
            x = [];

            if(obj.useTf)
                x = obj.varObj.targetBankAngle;
            end
        end

        function [lb, ub] = getBndsForVariable(obj)
            lb = obj.lb(obj.useTf);
            ub = obj.ub(obj.useTf);
        end

        function [lb, ub] = getAllBndsForVariable(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function setBndsForVariable(obj, lb, ub)
            obj.lb = lb;
            obj.ub = ub;
        end

        function useTf = getUseTfForVariable(obj)
            useTf = obj.useTf;
        end

        function setUseTfForVariable(obj, useTf)
            obj.useTf = useTf;
        end

        function updateObjWithVarValue(obj, x)
            obj.varObj.targetBankAngle = x;
        end

        function nameStrs = getStrNamesOfVars(obj, evtNum, varLocType)
            nameStrs = {sprintf('Event %i Bank Angle Termination Condition', evtNum)};
        end

        function varsStoredInRad = getVarsStoredInRad(obj)
            varsStoredInRad = true;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function mesh = ft_defacemesh(cfg, mesh)

% FT_DEFACEMESH allows you to de-identify a scalp surface mesh by erasing specific
% regions, such as the face and ears. The graphical user interface allows you to
% position a box over the anatomical data inside which all vertices will be removed.
% You might have to call this function multiple times when both face and ears need to
% be removed. Following defacing, you should check the result with FT_PLOT_MESH.
%
% Use as
%   mesh = ft_defacevolume(cfg, mesh)
%
% The configuration can contain the following options
%   cfg.translate  = initial position of the center of the box (default = [0 0 0])
%   cfg.scale      = initial size of the box along each dimension (default is automatic)
%   cfg.rotate     = initial rotation of the box (default = [0 0 0])
%   cfg.selection  = which voxels to keep, can be 'inside' or 'outside' (default = 'outside')
%
% See also FT_ANONYMIZEDATA, FT_DEFACEVCOLUME, FT_ANALYSISPIPELINE, FT_PLOT_MESH

% Copyright (C) 2015-2016, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

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

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

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% the actual work is done by FT_DEFACEVOLUME
tmpcfg = cfg;
tmpcfg.showcallinfo = 'no';
mesh = ft_defacevolume(tmpcfg, mesh);
% restore provenance information and put back cfg.callinfo
tmpcallinfo = cfg.showcallinfo;
[cfg, mesh] = rollback_provenance(cfg, mesh);
cfg.showcallinfo = tmpcallinfo;

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

author: spm


text:
function [indClinic,bestCost] = estimateVolClinicalRF(pathRF,nameOutcome,testComb,testCost,nSplit,testSplit,nBoot,seed,output)

% output: optional.

if nargin < 9
    output = 1; % By default
end

startpath = pwd;
tStart = tic;

% INITIALIZATIONS
cd(pathRF)
rng(seed), seeds = ceil(1000000000*rand(2,1));
training = load('training'); training = struct2cell(training); training = training{1};
load('volumeTrain'), Volume = volumeTrain; tableVol = table(Volume);
nComb = numel(testComb);
nCost = numel(testCost);
nTest = nComb *nCost;

% COMPUTATIONS
if strcmp(nameOutcome,'DeathSign')
    outcome = training.outcomes.Death;
else
    outcome = training.outcomes.(nameOutcome);
end
rng(seeds(1)), [trainSets,testSets] = getSplits(outcome,nSplit,testSplit);
maxMetric = 0; count = 0;
for c = 1:nCost
    cost = testCost(c);
    for t = 1:nComb
        tableComb = [tableVol,training.clinical.table(:,testComb{t})];
        cat = [false,logical(training.clinical.categories(testComb{t}))];
        meanAUC = 0; meanSens = 0; meanSpec = 0; count = count + 1;
        if output
            tic, fprintf('\n--> ESTIMATING RF PERFORMANCE USING %u SPLITS FOR TEST %u OF %u: %s ... ',nSplit,count,nTest,nameOutcome)
        end
        for s = 1:nSplit
            Ttrain = tableComb(trainSets(:,s),:); Ttest = tableComb(testSets(:,s),:);
            Ytrain = outcome(trainSets(:,s)); Ytest = outcome(testSets(:,s));
            rng(seeds(2)), [RF] = trainRF_table(Ttrain,Ytrain,cat,nBoot,cost);
            [prob] = predictRF(Ttest,RF);
            [aucSplit,sensSplit,specSplit,~] = calcPerformMetrics(prob,Ytest,0.5);
            meanAUC = meanAUC + aucSplit;
            meanSens = meanSens + sensSplit;
            meanSpec = meanSpec + specSplit;
        end
        meanAUC = meanAUC/nSplit; meanSens = meanSens/nSplit; meanSpec = meanSpec/nSplit;
        metric = 0.5*meanAUC + 0.5*(1-abs(meanSens - meanSpec));
        %metric = meanAUC;
        if metric >= maxMetric
            bestCost = cost;
            indClinic = testComb{t};
            maxMetric = metric;
        end
        if output
            fprintf('DONE!\n'), toc
        end
    end
end
if output
    time = toc(tStart);
    fprintf('\n\n\n-------------------------------------\n')
    fprintf('TOTAL TIME: %f seconds',time)
    fprintf('\n-------------------------------------\n')
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function combo_test23 ( )

%*****************************************************************************80
%
%% COMBO_TEST23 tests PART_SUCCESSOR and PART_SF_CONJUGATE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
  n = 8;

  fprintf ( 1, ' \n' );
  fprintf ( 1, 'COMBO_TEST23\n' );
  fprintf ( 1, '  PART_SUCCESSOR produces partitions of N,\n' );
  fprintf ( 1, '  PART_SF_CONJUGATE produces the conjugate of a partition.\n' );
  fprintf ( 1, ' \n' );
  fprintf ( 1, '  Partitions of N = %d\n', n );
  fprintf ( 1, ' \n' );
%
%  List.
%
  npart = 0;
  t = [];
  rank = -1;

  while ( 1 )

    rank_old = rank;

    [ npart, t, rank ] = part_successor ( n, npart, t, rank );

    if ( rank <= rank_old )
      break
    end

    fprintf ( 1, '  %3d   ', rank );
    for i = 1 : npart
      fprintf ( 1, '%3d', t(i) );
    end
    fprintf ( 1, '\n' );

    [ npartb, b ] = part_sf_conjugate ( n, npart, t );

    fprintf ( 1, '  Con   ', rank );
    for i = 1 : npartb
      fprintf ( 1, '%3d', b(i) );
    end
    fprintf ( 1, '\n' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%IMREAD  Loads an image from a file
%
%     img = cv.imread(filename)
%     img = cv.imread(filename, 'OptionName',optionValue, ...)
%
% ## Input
% * __filename__ Name of a file to be loaded.
%
% ## Output
% * __img__ Loaded image, whose depth and number of channels depend on options.
%
% ## Options
% * __Unchanged__ If set, return the loaded image as is (with alpha channel,
%   otherwise it gets cropped). Both the depth and number of channels are
%   unchanged as determined by the decoder. default false
% * __AnyDepth__ If set, return 16-bit/32-bit image when the input has the
%   corresponding depth, otherwise convert it to 8-bit. default false
% * __AnyColor__ If set, the image is read in any possible color format.
%   default false
% * __Color__ If set, always convert image to the 3 channel BGR color image.
%   default true
% * __Grayscale__ If set, always convert image to the single channel grayscale
%   image. default false
% * __GDAL__ If set, use the gdal driver for loading the image. default false
% * __ReduceScale__ Loads the image reduced by a scale factor (JPEG library
%   natively supports direct image scaling, other formats are resized after
%   loading). One of:
%   * `1`: no scaling (default).
%   * `2`: image scaled by 1/2 factor.
%   * `4`: image scaled by 1/4 factor.
%   * `8`: image scaled by 1/8 factor.
% * __IgnoreOrientation__ If set, do not rotate the image according to EXIF's
%   orientation flag. default false
% * __Flags__ Advanced option to directly set the flag specifying the depth
%   and color type of a loaded image. Note that setting this integer flag
%   overrides all the other flag options. Not set by default:
%   * `>0`: Return a 3-channel color image. Note that in the current
%           implementation the alpha channel, if any, is stripped from the
%           output image. For example, a 4-channel RGBA image is loaded as
%           RGB if `Flags >= 0`.
%   * `=0`: Return a grayscale image
%   * `<0`: Return the loaded image as is (with alpha channel if present).
% * __FlipChannels__ in case the output is color image, flips the color order
%   from OpenCV's BGR/BGRA to MATLAB's RGB/RGBA order. default true
%
% The function cv.imread loads an image from the specified file and returns
% it. If the image cannot be read (because of missing file, improper
% permissions, unsupported or invalid format), the function issues an error.
%
% Currently, the following file formats are supported:
%
% * Windows bitmaps - `*.bmp`, `*.dib` (always supported)
% * JPEG files - `*.jpeg`, `*.jpg`, `*.jpe` (see the Notes section)
% * JPEG 2000 files - `*.jp2` (see the Notes section)
% * Portable Network Graphics - `*.png` (see the Notes section)
% * WebP - `*.webp` (see the Notes section)
% * Portable image format - `*.pbm`, `*.pgm`, `*.ppm`, `*.pxm`, `*.pnm`,
%   `*.pam` (always supported)
% * Sun rasters - `*.sr`, `*.ras` (always supported)
% * TIFF files - `*.tiff`, `*.tif` (see the Notes section)
% * OpenEXR Image files - `*.exr` (see the Notes section)
% * Radiance HDR - `*.hdr`, `*.pic` (always supported)
% * Raster and Vector geospatial data supported by Gdal (see the Notes section)
% * DICOM medical images - `*.dcm` (see the Notes section)
%
% ### Notes
% * The function determines the type of an image by the content, not by the
%   file extension.
% * In the case of color images, the decoded images will have the channels
%   stored in BGR order. If `FlipChannels` is set, the channels are flipped to
%   RGB order.
% * On Microsoft Windows OS and MacOSX, the codecs shipped with an OpenCV
%   image (libjpeg, libpng, libtiff, and libjasper) are used by default. So,
%   OpenCV can always read JPEGs, PNGs, and TIFFs. On MacOSX, there is also an
%   option to use native MacOSX image readers. But beware that currently these
%   native image loaders give images with different pixel values because of
%   the color management embedded into MacOSX.
% * On Linux, BSD flavors and other Unix-like open-source operating systems,
%   OpenCV looks for codecs supplied with an OS image. Install the relevant
%   packages (do not forget the development files, for example, "libjpeg-dev",
%   in Debian and Ubuntu) to get the codec support or turn on the
%   `OPENCV_BUILD_3RDPARTY_LIBS` flag in CMake.
% * In the case you set `WITH_GDAL` flag to true in CMake and `GDAL` option to
%   load the image, then [GDAL driver](http://www.gdal.org) will be used in
%   order to decode the image by supporting the following formats:
%   * [Raster](http://www.gdal.org/formats_list.html),
%   * [Vector](http://www.gdal.org/ogr_formats.html).
% * If EXIF information are embedded in the image file, the EXIF orientation
%   will be taken into account and thus the image will be rotated accordingly
%   except if the option `IgnoreOrientation` is passed.
%
% See also: cv.imwrite, cv.imdecode, imread, imfinfo, imformats, dicomread,
%  hdrread, Tiff
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function det = dgbdi ( abd, lda, n, ml, mu, ipvt )

%*****************************************************************************80
%
%% DGBDI computes the determinant of a band matrix factored by DGBCO or DGBFA.
%
%  Discussion:
%
%    If the inverse is needed, use DGBSL N times.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real ABD(LDA,N), the output from DGBCO or DGBFA.
%
%    Input, integer LDA, the leading dimension of the array ABD.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer ML, MU, the number of diagonals below and above the
%    main diagonal.  0 <= ML < N, 0 <= MU < N.
%
%    Input, integer IPVT(N), the pivot vector from DGBCO or DGBFA.
%
%    Output, real DET(2), the determinant of the original matrix.
%      determinant = DET(1) * 10.0**DET(2)
%    with  1.0 <= abs ( DET(1) ) < 10.0 or DET(1) = 0.0.
%
  ten = 10.0;

  m = ml + mu + 1;
  det(1) = 1.0;
  det(2) = 0.0;

  for i = 1 : n

    if ( ipvt(i) ~= i )
      det(1) = -det(1);
    end

    det(1) = abd(m,i) * det(1);

    if ( det(1) == 0.0 )
      return
    end

    while ( abs ( det(1) ) < 1.0 )
      det(1) = ten * det(1);
      det(2) = det(2) - 1.0;
    end

    while ( ten <= abs ( det(1) ) )
      det(1) = det(1) / ten;
      det(2) = det(2) + 1.0;
    end

  end

  return
end


text:
function test_bug1850

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_prepare_neighbours ft_channelrepair
%
% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1850

load(dccnpath('/home/common/matlab/fieldtrip/data/test/latest/raw/meg/preproc_ctf275.mat'));

cfg = [];
cfg.channel = {'all', '-MRT23', '-MLP57'};
data = ft_selectdata(cfg, data);

cfg = [];
cfg.method = 'template';
cfg.template = 'CTF275_neighb.mat';
neighbours = ft_prepare_neighbours(cfg);

% get the full list of 275 channel names
allchans = {neighbours.label};
missingchans = setdiff(allchans, data.label);

% repair the two channels that were removed
cfg = [];
cfg.missingchannel = missingchans;
cfg.neighbours = neighbours;
cfg.method = 'spline';
data_repaired = ft_channelrepair(cfg,data);

author: fieldtrip

repo: fieldtrip


text:
function chebcoeffs = chebvals2chebcoeffs(chebvals, kind)
%CHEBVALS2CHEBCOEFFS  Convert Chebyshev values to coefficients.
% 	CHEBCOEFFS = CHEBVALS2CHEBCOEFFS(CHEBVALS), converts the column vector
%   CHEBVALS of values on a second-kind Chebyshev grid (i.e, F(CHEBPTS(N)))
%   to a vector CHEBCOEFFS of the Chebyshev coefficients of the series
%       F(X) = C_CHEB(1)*T0(X) + ... + C_CHEB(N)*T{N-1}(X).
%
% 	CHEBVALS2CHEBCOEFFS(CHEBVALS, 1) is similar, but assumes the entries in
% 	CHEBVALS come from evaluating on a first-kind Chebyshev grid, i.e.,
%   F(CHEBPTS(N,1))).
%
% See also CHEBTECH2.VALS2COEFFS, CHEBPTS.

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

% Default to second-kind points
if ( nargin == 1 )
    kind = 2;
end

if ( kind == 1 )
    % This command is a wrapper for chebtech2/vals2coeffs.
    chebcoeffs = chebtech1.vals2coeffs(chebvals);
elseif ( kind == 2 )
    % This command is a wrapper for chebtech1/vals2coeffs.
    chebcoeffs = chebtech2.vals2coeffs(chebvals);
else
    error('CHEBFUN:chebvals2chebcoeffs:kind', ...
        'Invalid Chebyshev kind. Must be 1 or 2.')
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function out = cfg_run_named_input(job)

% Return evaluated input.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_run_named_input.m 1716 2008-05-23 08:18:45Z volkmar $

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

out.input = job.input;


text:
clear, clc;

% This is an example for running the function LeastC
%
%  min  1/2 || A x - y||^2 + 1/2 * rsL2 * ||x||_2^2
%  s.t. ||x||_1 <= z
%
% For detailed description of the function, please refer to the Manual.
%
%% Related papers
%
% [1]  Jun Liu and Jieping Ye, Efficient Euclidean Projections
%      in Linear Time, ICML 2009.
%
% [2]  Jun Liu and Jieping Ye, Sparse Learning with Efficient Euclidean
%      Projections onto the L1 Ball, Technical Report ASU, 2008.
%
% [3]  Jun Liu, Jianhui Chen, and Jieping Ye,
%      Large-Scale Sparse Logistic Regression, KDD, 2009.
%
%% ------------   History --------------------
%
% First version on August 10, 2009.
%
% September 5, 2009: adaptive line search is added
%
% For any problem, please contact Jun Liu (j.liu@asu.edu)

cd ..
cd ..

root=cd;
addpath(genpath([root '/SLEP']));
                     % add the functions in the folder SLEP to the path

% change to the original folder
cd Examples/L1;

m=1000;  n=1000;    % The data matrix is of size m x n

% for reproducibility
randNum=1;

% ---------------------- generate random data ----------------------
randn('state',(randNum-1)*3+1);
A=randn(m,n);       % the data matrix

randn('state',(randNum-1)*3+2);
xOrin=randn(n,1);

randn('state',(randNum-1)*3+3);
noise=randn(m,1);
y=A*xOrin +...
    noise*0.01;     % the response

z=100;               % the radius of the L1 ball

%----------------------- Set optional items -----------------------
opts=[];

% Starting point
opts.init=2;            % starting from a zero point

% Termination criterion
opts.tFlag=5;          % run .maxIter iterations
opts.maxIter=100;      % maximum number of iterations

% Mormalization
opts.nFlag=0;         % without normalization

%opts.rsL2=0.1;        % the two norm regularization

%----------------------- Run the code LeastC -----------------------
fprintf('\n lFlag=0 \n');
opts.lFlag=0;       % Nemirovski's line search
tic;
[x1, funVal1, ValueL1]= LeastC(A, y, z, opts);
toc;

fprintf('\n lFlag=1 \n');
opts.lFlag=1;       % adaptive line search
opts.tFlag=2; opts.tol= funVal1(end);
tic;
[x2, funVal2, ValueL2]= LeastC(A, y, z, opts);
toc;

figure;
plot(funVal1,'-r');
hold on;
plot(funVal2,'--b');
legend('lFlag=0', 'lFlag=1');
xlabel('Iteration (i)');
ylabel('The objective function value');

% --------------------- compute the pathwise solutions ----------------
opts.fName='LeastC';    % set the function name to 'LeastC'
Z=[10, 100, 200, 500];  % set the parameters

% run the function pathSolutionLeast
fprintf('\n Compute the pathwise solutions, please wait...');
X=pathSolutionLeast(A, y, Z, opts);

author: zzstefan

repo: BrainNetClass


text:
% Usage: [xopt, fmin, retcode] = nlopt_minimize_constrained
%                                          (algorithm, f, f_data,
%                                           fc, fc_data, lb, ub,
%                                           xinit, stop)
%
% Minimizes a nonlinear multivariable function f(x, f_data{:}), subject
% to nonlinear constraints described by fc and fc_data (see below), where
% x is a row vector, returning the optimal x found (xopt) along with
% the minimum function value (fmin = f(xopt)) and a return code (retcode).
% A variety of local and global optimization algorithms can be used,
% as specified by the algorithm parameter described below.  lb and ub
% are row vectors giving the upper and lower bounds on x, xinit is
% a row vector giving the initial guess for x, and stop is a struct
% containing termination conditions (see below).
%
% This function is a front-end for the external routine
% nlopt_minimize_constrained in the free NLopt nonlinear-optimization
% library, which is a wrapper around a number of free/open-source
% optimization subroutines.  More details can be found on the NLopt
% web page (ab-initio.mit.edu/nlopt) and also under
% 'man nlopt_minimize_constrained' on Unix.
%
% f should be a handle (@) to a function of the form:
%
%    [val, gradient] = f(x, ...)
%
% where x is a row vector, val is the function value f(x), and gradient
% is a row vector giving the gradient of the function with respect to x.
% The gradient is only used for gradient-based optimization algorithms;
% some of the algorithms (below) are derivative-free and only require
% f to return val (its value).  f can take additional arguments (...)
% which are passed via the argument f_data: f_data is a cell array
% of the additional arguments to pass to f.  (Recall that cell arrays
% are specified by curly brackets { ... }.  For example, pass f_data={}
% for functions that require no additional arguments.)
%
% A few of the algorithms (below) support nonlinear constraints,
% in particular NLOPT_LD_MMA and NLOPT_LN_COBYLA.  These (if any)
% are specified by fc and fc_data.  fc is a cell array of
% function handles, and fc_data is a cell array of cell arrays of the
% corresponding arguments.  Both must have the same length m, the
% number of nonlinear constraints.  That is, fc{i} is a handle
% to a function of the form:
%
%   [val, gradient] = fc(x, ...)
%
% (where the gradient is only used for gradient-based algorithms),
% and the ... arguments are given by fc_data{i}{:}.
%
% If you have no nonlinear constraints, i.e. fc = fc_data = {}, then
% it is equivalent to calling the the nlopt_minimize() function,
% which omits the fc and fc_data arguments.
%
% stop describes the termination criteria, and is a struct with a
% number of optional fields:
%     stop.ftol_rel = fractional tolerance on function value
%     stop.ftol_abs = absolute tolerance on function value
%     stop.xtol_rel = fractional tolerance on x
%     stop.xtol_abs = row vector of absolute tolerances on x components
%     stop.fmin_max = stop when f < fmin_max is found
%     stop.maxeval = maximum number of function evaluations
%     stop.maxtime = maximum run time in seconds
%     stop.verbose = > 0 indicates verbose output
% Minimization stops when any one of these conditions is met; any
% condition that is omitted from stop will be ignored.  WARNING:
% not all algorithms interpret the stopping criteria in exactly the
% same way, and in any case ftol/xtol specify only a crude estimate
% for the accuracy of the minimum function value/x.
%
% The algorithm should be one of the following constants (name and
% interpretation are the same as for the C function).  Names with
% _G*_ are global optimization, and names with _L*_ are local
% optimization.  Names with _*N_ are derivative-free, while names
% with _*D_ are gradient-based algorithms.  Algorithms:
%
% NLOPT_GD_MLSL_LDS, NLOPT_GD_MLSL, NLOPT_GD_STOGO, NLOPT_GD_STOGO_RAND,
% NLOPT_GN_CRS2_LM, NLOPT_GN_DIRECT_L, NLOPT_GN_DIRECT_L_NOSCAL,
% NLOPT_GN_DIRECT_L_RAND, NLOPT_GN_DIRECT_L_RAND_NOSCAL, NLOPT_GN_DIRECT,
% NLOPT_GN_DIRECT_NOSCAL, NLOPT_GN_ISRES, NLOPT_GN_MLSL_LDS, NLOPT_GN_MLSL,
% NLOPT_GN_ORIG_DIRECT_L, NLOPT_GN_ORIG_DIRECT, NLOPT_LD_AUGLAG_EQ,
% NLOPT_LD_AUGLAG, NLOPT_LD_LBFGS, NLOPT_LD_LBFGS_NOCEDAL, NLOPT_LD_MMA,
% NLOPT_LD_TNEWTON, NLOPT_LD_TNEWTON_PRECOND,
% NLOPT_LD_TNEWTON_PRECOND_RESTART, NLOPT_LD_TNEWTON_RESTART,
% NLOPT_LD_VAR1, NLOPT_LD_VAR2, NLOPT_LN_AUGLAG_EQ, NLOPT_LN_AUGLAG,
% NLOPT_LN_BOBYQA, NLOPT_LN_COBYLA, NLOPT_LN_NELDERMEAD,
% NLOPT_LN_NEWUOA_BOUND, NLOPT_LN_NEWUOA, NLOPT_LN_PRAXIS, NLOPT_LN_SBPLX
%
% For more information on individual algorithms, see their individual
% help pages (e.g. "help NLOPT_LN_SBPLX").
function [xopt, fmin, retcode] = nlopt_minimize_constrained(algorithm, f, f_data, fc, fc_data, lb, ub, xinit, stop)

opt = stop;
if (isfield(stop, 'minf_max'))
  opt.stopval = stop.minf_max;
end
opt.algorithm = algorithm;
opt.min_objective = @(x) f(x, f_data{:});
opt.lower_bounds = lb;
opt.upper_bounds = ub;
% opt.local_optimizer = alg;
% for i = 1:length(fc)
%   opt.fc{i} = @(x) fc{i}(x, fc_data{i}{:});
% end
opt.fc = fc;
[xopt, fmin, retcode] = nlopt(opt, xinit);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
addpath /n/fs/vision/www/pvt/

latexFname = '/Users/xj/Documents/museumIJCV/museumCSG/ijcv.tex';
oriPath = '/Users/xj/Documents/museumIJCV/museumCSG/all';
newPath = '/Users/xj/Documents/museumIJCV/museumCSG/';

% this script parse a latex file and copy all figures into a new location
% this is usually used to get a clean version for submitting latex files
% to places like IJCV or ECCV that requires the latex file for final
% submission

str = file2string(latexFname);

pos = strfind(str,'\includegraphics[');

newline = sprintf('\n');

images = {};

for i=1:length(pos)
    ps = strfind(str(pos(i):end),'{');  ps = ps(1);
    pe = strfind(str(pos(i):end),'}');  pe = pe(1);

    nlp = strfind(str(1:pos(i)),newline);
    if isempty(nlp)
        nlp = 0;
    else
        nlp = nlp(end);
    end

    cmp = strfind(str(1:pos(i)),'%');
    if isempty(cmp)
        cmp = 0;
    else
        cmp = cmp(end);
    end

    if nlp>=cmp
        images{end+1} = str(pos(i)+ps:pos(i)+pe-2);
    end
end

for i=1:length(images)
    newF = fullfile(newPath,images{i});
    if ~exist(fileparts(newF),'dir')
        mkdir(fileparts(newF));
    end
    copyfile(fullfile(oriPath,images{i}),newF);
end

author: jianxiongxiao

repo: ProfXkit


text:
function par_estc
% transport parameter estimation with derivatives             Holzbecher January 2006

global xfit cfit T D c0 c1

% Example values for Chlorid in Marmara Sea Sediment Core
T = 3.15e11;  % [s] 10.000 years
D = 1.0e-5;  % [cm*cm/s]
c0 = 0;      % [mmol/l]
c1 = 619;    % [mmol/l]
xmax = 4000; % [cm]

% specify fitting data
xfit = [0 20 40 60 100 120 140 160 215 255 275 300 375 450 525 600 750 1050 1200 1350 1650 1950 2250 2550 2700 3000 3450 3900];
cfit = [619 597 608 615 619 615 621 571 621 618 619 625 577 612 608 612 609 590 582 582 556 494 457 489 487 444 381 371];

x = [0:xmax/400:xmax];
options = optimset('Display','iter','TolFun',1e-9);
v = fzero(@myfun,0.2e-8,options);
display (['Best fit for v = ' num2str(v)]);

h = 1./(2.*sqrt(D*T)); e = diag(eye(size(x,2)));
plot (xfit,cfit,'o',x,c0+0.5*c1*(erfc(h*(x-v*T*e'))+(exp((v/D)*x)).*erfc(h*(x+v*T*e'))),'-');
legend ('given','modelled');
xlabel ('depth [cm]'); ylabel ('chloride concentration [mmol/l]');
text(0.1*xmax,c1*0.65,['sedimentation velocity [cm/a]: ' num2str(v*3.15e7)]);
e = diag(eye(size(xfit,2)));
normc = norm(cfit-c0+0.5*c1*(erfc(h*(xfit-v*T*e'))+(exp((v/D)*xfit)).*erfc(h*(xfit+v*T*e'))));
text(0.1*xmax,c1*0.6,['norm of residuals: ' num2str(normc)]);

function f = myfun(v);
global xfit cfit T D c0 c1

e=diag(eye(size(xfit,2))); h=1./(2.*sqrt(D*T));
arg1 = h*(xfit-v*T*e'); arg2 = h*(xfit+v*T*e'); arg3 = (v/D)*xfit;

% solve advection diffusion equation for c with c(t=0)=c0 and c(x=0)=c1
c = c0 + 0.5*c1*(erfc(arg1)+(exp(arg3).*erfc(arg2)));

% compute derivative of solution due to v
cv = c1*((T*h/sqrt(pi))*(exp(-arg1.*arg1)-exp(arg3).*exp(-arg2.*arg2))+0.5*(xfit/D).*exp(arg3).*erfc(arg2));

% specify function f to vanish
f = 2*(c-cfit)*cv';

author: Sable

repo: mcbench-benchmarks


text:
%
%
% Copyright (c) 2012 University of Crete - Computer Science Department (UOC-CSD)
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

% Estimate the paramaters of a wrapped normal distribution
% http://en.wikipedia.org/wiki/Wrapped_normal_distribution
function [m s] = wrappednormestim(angles, dim)

    if nargin<2; dim=1; end

    expjangle = exp(1j*angles);

    z = mean(expjangle, dim);

    m = angle(z);

    s = sqrt(-2*log(abs(z)));

    if nargout==0
        hx = -pi:0.1:pi;
        hy = hist(angles, hx);
        hy = hy./sum(hy);

        hold off;
        plot(hx, (hy), 'k');
        hold on;
        p = wrappednormpdf(hx, m, s);
        p = p./sum(p);
        plot(hx, (p), 'b');
        ylim([0 max(p)+0.1]);
        title(['m=' num2str(m) ' s=' num2str(s)]);
    end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep