text:
function [estimate] = ft_inverse_pcc(sourcemodel, sens, headmodel, dat, C, varargin)

% FT_INVERSE_PCC implements a linearly-constrained miminum variance beamformer that
% allows for post-hoc computation of canonical or partial coherence or correlation.
% Moreover, if cortico-cortical interactions are computed, the spatial filters are
% computed with a paired dipole as sourcemodel, thus suppressing the distortive
% effect of correlated activity from the seed region.
%
% Use as
%   [estimate] = ft_inverse_pcc(sourcemodel, sens, headmodel, dat, cov, ...)
% where
%   sourcemodel is the input source model, see FT_PREPARE_SOURCEMODEL
%   sens        is the gradiometer or electrode definition, see FT_DATATYPE_SENS
%   headmodel   is the volume conductor definition, see FT_PREPARE_HEADMODEL
%   dat         is the data matrix with the ERP or ERF
%   cov         is the data covariance or cross-spectral density matrix
% and
%   estimate    contains the estimated source parameters
%
% Additional input arguments should be specified as key-value pairs and can include
%   'refchan'
%   'refdip'
%   'supchan'
%   'supdip'
%   'feedback'
%   'keepcsd'
%   'keepfilter'
%   'keepleadfield'
%   'keepmom'
%   'lambda'
%   'projectnoise'
%   'realfilter'
%   'fixedori'
%
% These options influence the forward computation of the leadfield
%   'reducerank'      = 'no' or number  (default = 3 for EEG, 2 for MEG)
%   'backproject'     = 'yes' or 'no', in the case of a rank reduction this parameter determines whether the result will be backprojected onto the original subspace (default = 'yes')
%   'normalize'       = 'no', 'yes' or 'column' (default = 'no')
%   'normalizeparam'  = parameter for depth normalization (default = 0.5)
%   'weight'          = number or Nx1 vector, weight for each dipole position to compensate for the size of the corresponding patch (default = 1)
%
% See also FT_SOURCEANALYSIS, FT_PREPARE_HEADMODEL, FT_PREPARE_SOURCEMODEL

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

if mod(nargin-5,2)
  % the first 5 arguments are fixed, the other arguments should come in pairs
  ft_error('invalid number of optional arguments');
end

% get the optional input arguments, or use defaults
refchan        = ft_getopt(varargin, 'refchan');
refdip         = ft_getopt(varargin, 'refdip');
supchan        = ft_getopt(varargin, 'supchan');
supdip         = ft_getopt(varargin, 'supdip');
feedback       = ft_getopt(varargin, 'feedback',      'text');
keepcsd        = ft_getopt(varargin, 'keepcsd',       'yes');
keepfilter     = ft_getopt(varargin, 'keepfilter',    'no');
keepleadfield  = ft_getopt(varargin, 'keepleadfield', 'no');
keepmom        = ft_getopt(varargin, 'keepmom',       'yes');
projectnoise   = ft_getopt(varargin, 'projectnoise',  'yes');
realfilter     = ft_getopt(varargin, 'realfilter',    'yes');
fixedori       = ft_getopt(varargin, 'fixedori',      'no');
weightnorm     = ft_getopt(varargin, 'weightnorm',    'no');

% construct the low-level options for the covariance matrix inversion as key-value pairs, these are passed to FT_INV
invopt = {};
invopt = ft_setopt(invopt, 'lambda',    ft_getopt(varargin, 'lambda', 0));
invopt = ft_setopt(invopt, 'kappa',     ft_getopt(varargin, 'kappa'));
invopt = ft_setopt(invopt, 'tolerance', ft_getopt(varargin, 'tol'));
invopt = ft_setopt(invopt, 'method',    ft_getopt(varargin, 'invmethod'));

% construct the low-level options for the leadfield computation as key-value pairs, these are passed to FT_COMPUTE_LEADFIELD
leadfieldopt = {};
leadfieldopt = ft_setopt(leadfieldopt, 'reducerank',     ft_getopt(varargin, 'reducerank'));
leadfieldopt = ft_setopt(leadfieldopt, 'backproject',    ft_getopt(varargin, 'backproject'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalize',      ft_getopt(varargin, 'normalize'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', ft_getopt(varargin, 'normalizeparam'));
leadfieldopt = ft_setopt(leadfieldopt, 'weight',         ft_getopt(varargin, 'weight'));

% convert the yes/no arguments to the corresponding logical values
fixedori       = istrue(fixedori);
keepcsd        = istrue(keepcsd);  % see below
keepfilter     = istrue(keepfilter);
keepleadfield  = istrue(keepleadfield);
keepmom        = istrue(keepmom);
projectnoise   = istrue(projectnoise);
realfilter     = istrue(realfilter);

if ~isequal(weightnorm, 'no')
  ft_error('weight normalization is not supported');
end

if ~keepcsd
  ft_warning('setting keepcsd to ''yes'', overruling specified option');
  keepcsd = true; % this needs to be true by definition
end

% flags to avoid calling isfield repeatedly in the loop over grid positions (saves a lot of time)
hasmom        = isfield(sourcemodel, 'mom');
hasleadfield  = isfield(sourcemodel, 'leadfield');
hasfilter     = isfield(sourcemodel, 'filter');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find the dipole positions that are inside/outside the brain
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isfield(sourcemodel, 'inside')
  if hasfilter
    sourcemodel.inside = ~cellfun(@isempty, sourcemodel.filter);
  elseif hasleadfield
    sourcemodel.inside = ~cellfun(@isempty, sourcemodel.leadfield);
  else
    sourcemodel.inside = ft_inside_headmodel(sourcemodel.pos, headmodel);
  end
end

% convert to logical representation
sourcemodel = fixinside(sourcemodel);

% keep the original details on inside and outside positions
originside = sourcemodel.inside;
origpos    = sourcemodel.pos;

% select only the dipole positions inside the brain for scanning
sourcemodel.pos    = sourcemodel.pos(originside,:);
sourcemodel.inside = true(size(sourcemodel.pos,1),1);

if hasmom
  sourcemodel.mom = sourcemodel.mom(:, originside);
end

if hasfilter
  ft_info('using precomputed filters\n');
  sourcemodel.filter = sourcemodel.filter(originside);
elseif hasleadfield
  ft_info('using precomputed leadfields\n');
  sourcemodel.leadfield = sourcemodel.leadfield(originside);
else
  ft_info('computing forward model on the fly\n');
end

if ~isempty(refdip)
  rf = ft_compute_leadfield(refdip, sens, headmodel, leadfieldopt{:});
else
  rf = [];
end

if ~isempty(supdip)
  sf = ft_compute_leadfield(supdip, sens, headmodel, leadfieldopt{:});
else
  sf = [];
end

% sanity check
if (~isempty(rf) || ~isempty(sf)) && hasfilter
  ft_error('precomputed filters cannot be used in combination with a refdip or supdip')
end

megchan  = setdiff(1:size(C,1), [refchan supchan]);
Nrefchan = length(refchan);
Nsupchan = length(supchan);
Nmegchan = length(megchan);
Cmeg     = C(megchan,megchan);   %  the filter uses the csd between all MEG channels

% it is difficult to give a quantitative estimate of lambda, therefore also
% support relative (percentage) measure that can be specified as string (e.g. '10%')
% the converted value needs to be passed on to ft_inv
lambda = ft_getopt(invopt, 'lambda');
if ~isempty(lambda) && ischar(lambda) && lambda(end)=='%'
  ratio  = sscanf(lambda, '%f%%');
  ratio  = ratio/100;
  lambda = ratio * trace(C)/size(C,1);
  invopt = ft_setopt(invopt, 'lambda', lambda);
end

if projectnoise
  % estimate the noise level in the covariance matrix by the smallest singular (non-zero) value
  noise = svd(Cmeg);
  noise = noise(rank(Cmeg));
  % estimated noise floor is equal to or higher than a numeric lambda
  noise = max(noise, lambda);
end

if realfilter
  % construct the filter only on the real part of the CSD matrix, i.e. filter is real
  invCmeg = ft_inv(real(Cmeg), invopt{:});
else
  % construct the filter on the complex CSD matrix, i.e. filter contains imaginary component as well
  % this results in a phase rotation of the channel data if the filter is applied to the data
  invCmeg = ft_inv(Cmeg, invopt{:});
end

ft_progress('init', feedback, 'scanning grid');
for i=1:size(sourcemodel.pos,1)
  ft_progress(i/size(sourcemodel.pos,1), 'scanning grid %d/%d\n', i, size(sourcemodel.pos,1));

  if hasfilter
    % precomputed filter is provided, the leadfield is not needed
  elseif hasleadfield && hasmom && size(sourcemodel.mom, 1)==size(sourcemodel.leadfield{i}, 2)
    % reuse the leadfield that was previously computed and project
    lf = sourcemodel.leadfield{i} * sourcemodel.mom(:,i);
  elseif  hasleadfield &&  hasmom
    % reuse the leadfield that was previously computed but don't project
    lf = sourcemodel.leadfield{i};
  elseif  hasleadfield && ~hasmom
    % reuse the leadfield that was previously computed
    lf = sourcemodel.leadfield{i};
  elseif ~hasleadfield &&  hasmom
    % compute the leadfield for a fixed dipole orientation
    lf = ft_compute_leadfield(sourcemodel.pos(i,:), sens, headmodel, leadfieldopt{:}) * sourcemodel.mom(:,i);
  else
    % compute the leadfield
    lf = ft_compute_leadfield(sourcemodel.pos(i,:), sens, headmodel, leadfieldopt{:});
  end

  % determine the dimensionality of the source (usually 3, or 6 for a dipole pair)
  if hasfilter
    sourcedim = size(sourcemodel.filter{i},1);
  else
    % concatenate scandip, refdip and supdip
    lfa = [lf rf sf];
    sourcedim = size(lfa,2);
  end

  if fixedori
    if isempty(refdip) && isempty(supdip) && isempty(refchan) && isempty(supchan)
      % compute the leadfield for the optimal dipole orientation
      % subsequently the leadfield for only that dipole orientation will
      % be used for the final filter computation
      if hasfilter && size(sourcemodel.filter{i},1)==1
        % nothing to do
        ft_warning('Ignoring ''fixedori''. The fixedori option is supported only if there is ONE dipole for location.')
      else
        if hasfilter && size(sourcemodel.filter{i},1)~=1
          filt = sourcemodel.filter{i};
        else
          filt = pinv(lfa' * invCmeg * lfa) * lfa' * invCmeg;
        end
        [u, s, v] = svd(real(filt * Cmeg * ctranspose(filt)));
        maxpowori = u(:,1);
        if numel(s)>1
          eta = s(1,1)./s(2,2);
        else
          eta = nan;
        end
        lfa  = lfa * maxpowori;
        estimate.ori{i} = maxpowori;
        estimate.eta(i) = eta;
        % update the number of dipole components
        sourcedim = size(lfa,2);
      end
    else
      ft_warning('Ignoring ''fixedori''. The fixedori option is supported only if there is ONE dipole for location.')
    end
  end

  if hasfilter
    % use the provided filter
    filt = sourcemodel.filter{i};
  else
    % construct the spatial filter
    filt = pinv(lfa' * invCmeg * lfa) * lfa' * invCmeg;                % use PINV/SVD to cover rank deficient leadfield
  end

  % concatenate the source filters with the channel filters
  filtn = zeros(sourcedim+Nrefchan+Nsupchan, Nmegchan+Nrefchan+Nsupchan);
  % this part of the filter relates to the sources
  filtn(1:sourcedim,megchan) = filt;
  % this part of the filter relates to the channels
  filtn((sourcedim+1):end,setdiff(1:(Nmegchan+Nrefchan+Nsupchan), megchan)) = eye(Nrefchan+Nsupchan);
  filt = filtn;
  clear filtn

  if keepcsd
    estimate.csd{i,1} = filt * C * ctranspose(filt);
  end
  if projectnoise
    estimate.noisecsd{i,1} = noise * (filt * ctranspose(filt));
  end
  if keepmom && ~isempty(dat)
    estimate.mom{i,1} = filt * dat;
  end
  if keepfilter
    estimate.filter{i,1} = filt;
  end
  if keepleadfield && needleadfield
    estimate.leadfield{i,1} = lf;
  end

  % remember how all components in the output csd should be interpreted
  %scandiplabel = repmat({'scandip'}, 1, size(lf, 2));    % based on last leadfield
  scandiplabel = repmat({'scandip'}, 1, size(filt, 1)-size(rf, 2)-size(sf, 2)-Nrefchan-Nsupchan); % robust if lf does not exist
  refdiplabel  = repmat({'refdip'},  1, size(rf, 2));
  supdiplabel  = repmat({'supdip'},  1, size(sf, 2));
  refchanlabel = repmat({'refchan'}, 1, Nrefchan);
  supchanlabel = repmat({'supchan'}, 1, Nsupchan);
  % concatenate all the labels
  estimate.csdlabel{i,1} = [scandiplabel refdiplabel supdiplabel refchanlabel supchanlabel];

end % for each dipole position
ft_progress('close');

% wrap it all up, prepare the complete output
estimate.inside  = originside;
estimate.pos     = origpos;

fnames_cell   = {'leadfield' 'filter' 'ori' 'csd' 'noisecsd' 'subspace' 'mom' 'csdlabel'};
for k = 1:numel(fnames_cell)
  if isfield(estimate, fnames_cell{k})
    estimate.(fnames_cell{k})( originside) = estimate.(fnames_cell{k});
    estimate.(fnames_cell{k})(~originside) = {[]};
  end
end

fnames_scalar = {'pow' 'noise' 'eta' 'coh'};
for k = 1:numel(fnames_scalar)
  if isfield(estimate, fnames_scalar{k})
    estimate.(fnames_scalar{k})( originside) = estimate.(fnames_scalar{k});
    estimate.(fnames_scalar{k})(~originside) = nan;
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to compute the pseudo inverse. This is the same as the
% standard MATLAB function, except that the default tolerance is twice as
% high.
%   Copyright 1984-2004 The MathWorks, Inc.
%   $Revision$  $Date: 2009/01/07 13:12:03 $
%   default tolerance increased by factor 2 (Robert Oostenveld, 7 Feb 2004)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function X = pinv(A,varargin)
[m,n] = size(A);
if n > m
  X = pinv(A',varargin{:})';
else
  [U,S,V] = svd(A,0);
  if m > 1, s = diag(S);
  elseif m == 1, s = S(1);
  else s = 0;
  end
  if nargin == 2
    tol = varargin{1};
  else
    tol = 10 * max(m,n) * max(s) * eps;
  end
  r = sum(s > tol);
  if (r == 0)
    X = zeros(size(A'),class(A));
  else
    s = diag(ones(r,1)./s(1:r));
    X = V(:,1:r)*s*U(:,1:r)';
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function p = presolveSOS(p)
if ~isempty(p.sosgroups)
    for i = 1:length(p.sosgroups)
        for j = 1:length(p.sosgroups{i})
            % Probe what happens if we activate this binary
            k = p.sosgroups{i}(j);
            other = setdiff(p.sosgroups{i},k);
            if p.ub(k) && ~any(p.lb(p.sosgroups{i}))
                % Still free to be activated
                p_ = p;
                p_.ub(other) = 0;
                p_.lb(k) = 1;
                p_ = smashFixed(p_);
                if ~isempty(fixedInfeasibleEquality(p_))
                    p.ub(k) = 0;
                end
                % Try de-activate too
                p_ = p;
                p_.ub(k) = 0;
                p_ = smashFixed(p_);
                if ~isempty(fixedInfeasibleEquality(p_))
                    p.lb(k) = 1;
                end
            end
        end
    end
    % Prune in case some are presolved completely
    for i = 1:length(p.sosgroups)
        if nnz(p.ub(p.sosgroups{i}))==1
            keep(i) = 0;
        else
            keep(i) = 1;
        end
    end
    if ~all(keep)
        p.sosgroups = {p.sosgroups{find(keep)}};
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function scanInfoS = initializeScanInfo
%"initializeScanInfo"
%   Return a skeleton scanInfo struct of dimension 1x0.  Only the fields
%   exist.
%
%JRA 06/28/06
%AI  12/28/16 Added scale slope, intercept
%Usage:
%   scanInfoS = initializeScanInfo
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

if nargout > 1
   error('Invalid number of outputs requested from initializeCERR.  Check for the new single output prototype!');
end

scanInfoS =  struct(...
'imageNumber'              ,   '', ...
'imageType'                ,   '', ...
'caseNumber'               ,   '', ...
'patientName'              ,   '', ...
'patientID'                ,   '', ...
'patientBirthDate'         ,   '', ...
'scanType'                 ,   '', ...
'CTOffset'                 ,   '', ...
'rescaleSlope'             ,   '', ...
'rescaleIntercept'         ,   '', ...
'scaleSlope'               ,   '', ...
'scaleIntercept'           ,   '', ...
'realWorldValueSlope'      ,   '', ...
'realWorldValueIntercept'  ,   '', ...
'realWorldMeasurCodeMeaning',  '', ...
'philipsImageUnits'        ,   '', ...
'grid1Units'               ,   '', ...
'grid2Units'               ,   '', ...
'numberRepresentation'     ,   '', ...
'bytesPerPixel'            ,   '', ...
'numberOfDimensions'       ,   '', ...
'sizeOfDimension1'         ,   '', ...
'sizeOfDimension2'         ,   '', ...
'zValue'                   ,   '', ...
'xOffset'                  ,   '', ...
'yOffset'                  ,   '', ...
'CTAir'                    ,   '', ...
'CTWater'                  ,   '', ...
'sliceThickness'           ,   '', ...
'siteOfInterest'           ,   '', ...
'unitNumber'               ,   '', ...
'seriesDescription'          ,   '', ...
'scannerType'              ,   '', ...
'manufacturer'             ,   '', ...
'scanFileName'             ,   '', ...
'headInOut'                ,   '', ...
'positionInScan'           ,   '', ...
'patientAttitude'          ,   '', ...
'bValue'                   ,   '', ...
'acquisitionDate'          ,   '', ...
'acquisitionTime'          ,   '', ...
'patientWeight'            ,   '', ...
'patientSize'              ,   '', ...
'patientBmi'               ,   '', ...
'patientSex'               ,   '', ...
'radiopharmaInfoS'         ,   '', ...
'injectionTime'            ,   '', ...
'injectionDate'            ,   '', ...
'injectedDose'             ,   '', ...
'halfLife'                 ,   '', ...
'imageUnits'               ,   '', ...
'suvType'                  ,   '', ...
'petCountSource'           ,   '', ...
'petSeriesType'            ,   '', ...
'petActivityConcentrationScaleFactor', '', ...
'petNumSlices'             ,   '', ...
'petIsDecayCorrected'      ,   '', ...
'petPrimarySourceOfCounts' ,   '', ...
'petDecayCorrectionDateTime',  '', ...
'decayCorrection'          ,   '', ...
'correctedImage'           ,   '', ...
'seriesDate'               ,   '', ...
'seriesTime'               ,   '', ...
'studyDate'                ,   '', ...
'studyTime'                ,   '', ...
'tapeOfOrigin'             ,   '', ...
'studyNumberOfOrigin'      ,   '', ...
'scanID'                   ,   '', ...
'scanNumber'               ,   '', ...
'scanDate'                 ,   '', ...
'CTScale'                  ,   '', ...
'distrustAbove'            ,   '', ...
'imageSource'              ,   '', ...
'transferProtocol'         ,   '', ...
'LRflippedToMatchPACS'     ,   '', ...
'APflippedToMatchPACS'     ,   '', ...
'SIflippedToMatchPACS'     ,   '', ...
'studyInstanceUID'         ,   '', ...
'seriesInstanceUID'        ,   '', ...
'sopInstanceUID'           ,   '', ...
'sopClassUID'              ,   '', ...
'frameOfReferenceUID'      ,   '', ...
'patientPosition'          ,   '', ...
'imageOrientationPatient'  ,   '', ...
'imagePositionPatient'     ,   '', ...
'windowCenter'             ,   '', ...
'windowWidth'              ,   '', ...
'temporalPositionIndex'    ,   '', ...
'frameAcquisitionDuration' ,   '', ...
'frameReferenceDateTime'   ,   '', ...
'DICOMHeaders'             ,   '');

scanInfoS(1) = [];

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function fem1d_pack_test ( )

%*****************************************************************************80
%
%% FEM1D_PACK_TEST tests the FEM1D_PACK library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2011
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_PACK_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the FEM1D_PACK library.\n' );

  fem1d_pack_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_PACK_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function fixed = p00_fixed_points ( test, m, fixed_num );

%*****************************************************************************80
%
%% P00_FIXED_POINTS returns the fixed points in a problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TEST, the index of the test problem
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer FIXED_NUM, the number of fixed points.
%
%    Output, real FIXED(M,FIXED_NUM), the
%    coordinates of the fixed points.
%
  if ( test == 1 )
    fixed = p01_fixed_points ( m, fixed_num );
  elseif ( test == 2 )
    fixed = p02_fixed_points ( m, fixed_num );
  elseif ( test == 3 )
    fixed = p03_fixed_points ( m, fixed_num );
  elseif ( test == 4 )
    fixed = p04_fixed_points ( m, fixed_num );
  elseif ( test == 5 )
    fixed = p05_fixed_points ( m, fixed_num );
  elseif ( test == 6 )
    fixed = p06_fixed_points ( m, fixed_num );
  elseif ( test == 7 )
    fixed = p07_fixed_points ( m, fixed_num );
  elseif ( test == 8 )
    fixed = p08_fixed_points ( m, fixed_num );
  elseif ( test == 9 )
    fixed = p09_fixed_points ( m, fixed_num );
  elseif ( test == 10 )
    fixed = p10_fixed_points ( m, fixed_num );
  elseif ( test == 11 )
    fixed = p11_fixed_points ( m, fixed_num );
  elseif ( test == 12 )
    fixed = p12_fixed_points ( m, fixed_num );
  elseif ( test == 13 )
    fixed = p13_fixed_points ( m, fixed_num );
  elseif ( test == 14 )
    fixed = p14_fixed_points ( m, fixed_num );
  elseif ( test == 15 )
    fixed = p15_fixed_points ( m, fixed_num );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_FIXED_POINTS - Fatal error!\n' );
    fprintf ( 1, '  Input value of TEST is out of range.\n' );
    error ( 'P00_FIXED_POINTS - Fatal error!' );
  end

  return
end


text:
% function to load and view the registered binaries output from Suite2P

% db file
make_db_MP024;
% which experiment to view
iexp = 5;

db = db(iexp);
ops.ResultsSavePath = 'D:/DATA/F/';
CharSubDirs = '';
for k = 1:length(db.expts)
    CharSubDirs    = [CharSubDirs num2str(db.expts(k)) '_'];
end
CharSubDirs = CharSubDirs(1:end-1);

ops.ResultsSavePath = sprintf('%s//%s//%s//%s//', ops.ResultsSavePath, db.mouse_name, db.date, ...
    CharSubDirs);
regpath = sprintf('%s/regops_%s_%s.mat', ops.ResultsSavePath, db.mouse_name, db.date);
load(regpath);

%%
iplane = 5; % which plane of the recording
ops = ops1{iplane};
fid  = fopen(ops.RegFile, 'r'); % opens the registered binary
[Ly Lx] = size(ops.mimg1);    % size of binary in x,y

clf;
nt0 = 0;
nimgbatch = 2000; % how many frames you can load at a time (will depend on RAM)
while 1
  data = fread(fid, Ly*Lx*nimgbatch, '*int16');
  if isempty(data)
    break;
  end
  data = reshape(data, Ly, Lx, []);
  data = data(ops.yrange, ops.xrange, :);
  NT   = size(data,3);
  for j = 1:NT
    imagesc(data(:,:,j),[0 2000]);%,[1000 3000]);
    title(sprintf('frame %d', j+nt0));
    axis square;
    drawnow;

    pause(0.25);
  end
  nt0 = nt0 + NT;
end
fclose all;

author: cortex-lab


text:
function  [data_hull,sources] = get_reduced_graph_weights_sources(data,hull);

data_hull.hull = hull;
switch data.dag_type
	case 'grid + input_space'

		affinity_hull = cell(1,size(hull,1));
		weights_hull = zeros(1,size(hull,1))';
		for i=1:size(hull,1)
			affinity_hull{i} = find(all(hull - repmat(hull(i,:),size(hull,1),1) >= 0,2));
			temp = 1;
			if all(hull(i,:)==1)

				weights_hull(i) = data.weight0;
			else
				for j=1:data.p, temp = temp * data.weightA^(hull(i,j)-1); end
				weights_hull(i) = temp;
			end
		end

		% prepare data
		X_hull = [];
		groups_hull = cell(1,size(hull,1));
		groupdescs_hull = cell(1,size(hull,1));

		ai1=1;
		for i1=1:size(hull,1)
			Xloc = get_data_reduced(hull(i1,:),data);
			X_hull = [ X_hull, Xloc ];
			groups_hull{i1} = ai1:ai1+size(Xloc,2)-1;   % records the indices of that group
			ai1 = ai1 + size(Xloc,2);
		end

		for i1=1:size(hull,1)
			groupdescs_hull{i1}=[];
			for i2=1:length(affinity_hull{i1})
				groupdescs_hull{i1} = [ groupdescs_hull{i1}, groups_hull{affinity_hull{i1}(i2)} ];
			end
		end

		qq =data.qs;

		% sources = find_sources_complement_grid_fast_c(hull,data.qs);
        % this function finds the sources of the complement of the active
        % set defined by the hull

        if nargout>1

        		sources = compute_sources(hull,data.qs);
% 		sources_old =find_sources_complement_grid_fast_int_c( int32(hull'),int32(data.qs));

        end
		data_hull.X = X_hull;
		data_hull.weights = weights_hull;
		data_hull.affinity = affinity_hull;
		data_hull.groups = groups_hull;
		data_hull.groupdescs = groupdescs_hull;

	case 'mkl + input_space'

		affinity_hull = cell(1,size(hull,1));
		weights_hull = zeros(1,size(hull,1))';
		for i=1:size(hull,1)
			affinity_hull{i} = find(all(hull - repmat(hull(i,:),size(hull,1),1) >= 0,2));
			temp = 1;
			if all(hull(i,:)==1)

				weights_hull(i) = data.weight0;
			else

				weights_hull(i) = data.weightA;
			end
		end

		% prepare data
		X_hull = [];
		groups_hull = cell(1,size(hull,1));
		groupdescs_hull = cell(1,size(hull,1));

		included_kernels = [];
		ai1=1;
		for i1=1:size(hull,1)
			Xloc = get_data_reduced(hull(i1,:),data);
			X_hull = [ X_hull, Xloc ];
			groups_hull{i1} = ai1:ai1+size(Xloc,2)-1;   % records the indices of that group
			ai1 = ai1 + size(Xloc,2);
		end

		for i1=1:size(hull,1)
			groupdescs_hull{i1}=[];
			for i2=1:length(affinity_hull{i1})
				groupdescs_hull{i1} = [ groupdescs_hull{i1}, groups_hull{affinity_hull{i1}(i2)} ];
			end
		end

		%sources = find_sources_complement_mkl(hull,data.p,data.qs);
		included_kernels = find(max(hull,[],1)==2);
		temp = 1:data.p;
		temp(included_kernels) = [];
		sources = [];
		for i=1:length(temp)
			hloc = ones(1,data.p);

			hloc(temp(i))=2;
			sources =  [ sources; hloc];
		end

		data_hull.X = X_hull;
		data_hull.weights = weights_hull;
		data_hull.affinity = affinity_hull;
		data_hull.groups = groups_hull;
		data_hull.groupdescs = groupdescs_hull;

	case 'bimkl + input_space'

		affinity_hull = cell(1,size(hull,1));
		weights_hull = zeros(1,size(hull,1))';
		for i=1:size(hull,1)
			temp = 1;
			if all(hull(i,:)==1)
				affinity_hull{i} =  1:size(hull,1);
				weights_hull(i) = data.weight0;
			else
				affinity_hull{i} =i;
				for j=1:data.p, temp = temp * data.weightA^(hull(i,j)-1); end
				weights_hull(i) = temp;
			end
		end

		% prepare data
		X_hull = [];
		groups_hull = cell(1,size(hull,1));
		groupdescs_hull = cell(1,size(hull,1));

		ai1=1;
		for i1=1:size(hull,1)
			Xloc = get_data_reduced(hull(i1,:),data);
			X_hull = [ X_hull, Xloc ];
			groups_hull{i1} = ai1:ai1+size(Xloc,2)-1;   % records the indices of that group
			ai1 = ai1 + size(Xloc,2);
		end

		for i1=1:size(hull,1)
			groupdescs_hull{i1}=[];
			for i2=1:length(affinity_hull{i1})
				groupdescs_hull{i1} = [ groupdescs_hull{i1}, groups_hull{affinity_hull{i1}(i2)} ];
			end
		end

		qq =data.qs;
		save temp hull qq

		p = data.p;

		included_kernels_one = zeros(p,1);
		included_kernels_two = zeros(p,p);
		for i=1:size(hull,1);
			tt = find(hull(i,:)==2);
			if length(tt)==1
				included_kernels_one(tt)=1;
			elseif length(tt)==2
				included_kernels_one(tt(1),tt(2)) = 1;
				included_kernels_one(tt(2),tt(1)) = 1;
			end
		end

		sources = [];
		for i=1:p
			if (included_kernels_one(i)==0),
				hloc = ones(1,p);
				hloc(i)=2;
				sources =  [ sources; hloc];
			end
		end
		for i=2:p
			for j=1:i-1
				if (included_kernels_two(i,j)==0),
					hloc = ones(1,p);
					hloc(i)=2; hloc(j)=2;
					sources =  [ sources; hloc];
				end
			end
		end

		data_hull.X = X_hull;
		data_hull.weights = weights_hull;
		data_hull.affinity = affinity_hull;
		data_hull.groups = groups_hull;
		data_hull.groupdescs = groupdescs_hull;

end
data_hull.dag_type = data.dag_type;

% check that the total dimension of the input space is smaller than the
% number of observations. If not -> store kernel matrices
switch data.dag_type
	case {'grid + input_space', 'mkl + input_space', 'bimkl + input_space'}

		if size(data_hull.X,2) > size(data_hull.X,1)
			n = size(data_hull.X,1);
			data_hull.kernels = zeros(n*(n+1)/2,size(hull,1),'single');
			for i=1:size(hull,1)
				data_hull.kernels(:,i) = symmetric_vectorize_single( single(data_hull.X(:,data_hull.groups{i})*data_hull.X(:,data_hull.groups{i})' ) );
			end
			switch data.dag_type
				case 'grid + input_space'
					data_hull.dag_type = 'grid + kernels';
				case 'mkl + input_space'
					data_hull.dag_type = 'mkl + kernels';
				case 'bimkl + input_space'
					data_hull.dag_type = 'bimkl + kernels';
			end

		end

end

repo: REST


text:
function A = meanOne(hyp, x, i)

% One mean function. The mean function does not have any parameters.
%
% m(x) = 1
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2010-08-04.
%
% See also MEANFUNCTIONS.M.

if nargin<2, A = '0'; return; end             % report number of hyperparameters
if nargin==2
  A = ones(size(x,1),1);                                         % evaluate mean
else
  A = zeros(size(x,1),1);                                           % derivative
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [ varargout ] = Beam_Propagation_Method_GUI ( varargin )

%% GUI for Beam Propagation Method (BPM)
%
%%
%    Author:   Ido Schwartz (IS)
%    April 2011 (ver. 1.0.0)
%
%% build GUI
fig = build_gui;
h   = guidata ( fig );

switch nargout      %----------------------- check nargout
    case 0
        % do nothing
    case 1                                  % return a handle to GUI
        varargout(1) = {fig};
    otherwise
        error ('Invalid number of output arguments');
end

e = exist ('data_bpm_tmp.mat','file');
if e==2
    delete('data_bpm_tmp.mat');
end

end

%% uifunctions
function fig = build_gui

setappdata(0,'UseNativeSystemDialogs',false);

fig = figure;

set (fig, 'Name', ' Beam Propagation Method GUI','NumberTitle', 'off', 'Visible','off');
set (fig, 'deleteFcn', @close_figure_callback ); %%%%%%%%%%%%%%%
set (fig, 'Color', get(0,'defaultUicontrolBackgroundColor') );

back_color = get (fig,'Color');
edit_box_back_color = [1 1 1];      % white

%% ------------------------------------ draw tools panel
tools_panel = uipanel('units', 'normalized', 'position' , [.005 .005 .99 .11] ,...  %% draw panel
    'parent' , fig,    'tag' , ' tools_panel'  , 'TitlePosition', 'righttop', ...
    'title' , 'tools');

uicontrol ('style', 'pushbutton', 'units', 'normalized', 'position', [0.003 0.1 0.12 0.8],...
    'HorizontalAlignment','center', 'FontWeight', 'bold','fontsize' , 12 ,...
    'string', 'simulation panel'  ,'tag', 'graphs_pushbutton', 'callback', @graphs_pushbutton_callback, 'parent' , tools_panel );

uicontrol ('style', 'pushbutton', 'units', 'normalized', 'position', [0.126 0.1 0.12 0.8],...
    'HorizontalAlignment','center', 'FontWeight', 'bold','fontsize' , 12 ,...
    'string', 'set parameters'  ,'tag', 'setparameters_pushbutton', 'callback', @setparameters_pushbutton_callback, 'parent' , tools_panel );

uicontrol ('style', 'pushbutton', 'units', 'normalized', 'position', [0.877 0.1 0.12 0.8],...
    'HorizontalAlignment','center', 'FontWeight', 'bold','fontsize' , 12 ,...
    'string', 'Exit'  ,'tag', 'exit_pushbutton', 'callback', @exit_pushbutton_callback, 'parent' , tools_panel );

%% ------------------------------------ draw graphs panel
graph_panel = uipanel('units', 'normalized', 'position' , [.005 .11 .99 .885] ,...  %% draw panel
    'parent' , fig,    'tag' , 'g_panel'  , 'TitlePosition', 'righttop', ...
    'title' , 'Graphs');

uicontrol ('style', 'pushbutton', 'units', 'normalized', 'position', [0.003 0.003 0.12 0.07],...
    'HorizontalAlignment','center', 'FontWeight', 'bold','fontsize' , 12 ,...
    'string', 'start'  ,'tag', 'start_pushbutton', 'callback', @start_pushbutton_callback, 'parent' , graph_panel );

uicontrol ('style', 'pushbutton', 'units', 'normalized', 'position', [0.126 0.003 0.12 0.07],...
    'HorizontalAlignment','center', 'FontWeight', 'bold','fontsize' , 12 ,'enable','off',...
    'string', 'stop'  ,'tag', 'stop_pushbutton', 'callback', @stop_pushbutton_callback, 'parent' , graph_panel );

uicontrol ('style', 'pushbutton', 'units', 'normalized', 'position', [0.875 0.003 0.12 0.07],...
    'HorizontalAlignment','center', 'FontWeight', 'bold','fontsize' , 12 ,'enable','on',...
    'string', 'save'  ,'tag', 'save_pushbutton', 'callback', @save_pushbutton_callback, 'parent' , graph_panel );

axes('parent',graph_panel,'tag','ax1','Position', [.04 .125 .58 .83]'); %define the axes space
axes('parent',graph_panel,'tag','ax2','Position',[.665 .4 .3 .5]); %define the axes space

%% ------------------------------------ draw Set Parameters panel
parameters_panel = uipanel('units', 'normalized', 'position' , [.005 .11 .99 .885] ,...  %% draw panel
    'parent' , fig,    'tag' , 'p_panel'  , 'TitlePosition', 'righttop', ...
    'title' , 'Set Parameters');
% --------------------------------------------------------- Field parameters
field_panel = uipanel('units', 'normalized', 'position' , [.005 .005 .3 .99] ,...  %% draw panel
    'parent' , parameters_panel,    'tag' , 'field_panel'  , 'TitlePosition', 'righttop', ...
    'title' , 'Field Initial Parameters');

uicontrol ('style', 'text', 'units', 'normalized', 'position', [0.05 0.9 0.9 0.1],...
    'HorizontalAlignment','left', 'FontWeight', 'bold','fontsize' , 12 ,...
    'string', 'F(r,z=0)=A*exp(-0.5*(r/w)^2)'  ,'tag', 'field_text', 'parent' , field_panel );

uicontrol ('style', 'text', 'units', 'normalized', 'position', [0.05 0.83 0.11 0.07],...
    'HorizontalAlignment','left', 'FontWeight', 'normal','fontsize' , 12 ,...
    'string', 'A = '  ,'tag', 'A_text', 'parent' , field_panel );

uicontrol ('style', 'edit', 'units', 'normalized', 'position', [0.165 0.83 0.15 0.07],...
    'HorizontalAlignment','center', 'FontWeight', 'normal','fontsize' , 12 ,...
    'string', '1' , 'backgroundcolor',[1 1 1],'tag', 'A_edit','enable','on', 'parent' , field_panel );

uicontrol ('style', 'text', 'units', 'normalized', 'position', [0.05 0.75 0.11 0.07],...
    'HorizontalAlignment','left', 'FontWeight', 'normal','fontsize' , 12 ,...
    'string', 'w = '  ,'tag', 'w_text', 'parent' , field_panel );

uicontrol ('style', 'edit', 'units', 'normalized', 'position', [0.165 0.75 0.15 0.07],...
    'HorizontalAlignment','center', 'FontWeight', 'normal','fontsize' , 12 ,...
    'string', '10' ,'backgroundcolor',[1 1 1],'tag', 'w_edit','enable','on', 'parent' , field_panel );

uicontrol ('style', 'text', 'units', 'normalized', 'position', [0.33 0.75 0.22 0.07],...
    'HorizontalAlignment','left', 'FontWeight', 'normal','fontsize' , 12 ,...
    'string', 'micrometer '  ,'tag', 'w_units_text', 'parent' , field_panel );

uicontrol ('style', 'text', 'units', 'normalized', 'position', [0.05 0.67 0.11 0.07],...
    'HorizontalAlignment','left', 'FontWeight', 'normal','fontsize' , 12 ,...
    'string', 'wave length'  ,'tag', 'lambda_text', 'parent' , field_panel );

uicontrol ('style', 'edit', 'units', 'normalized', 'position', [0.165 0.67 0.15 0.07],...
    'HorizontalAlignment','center', 'FontWeight', 'normal','fontsize' , 12 ,...
    'string', '0.5' ,'backgroundcolor',[1 1 1],'tag', 'lambda_edit','enable','on', 'parent' , field_panel );

uicontrol ('style', 'popupmenu', 'units', 'normalized', 'position', [0.33 0.67 0.3 0.07],...
    'HorizontalAlignment','left', 'FontWeight', 'normal','fontsize' , 12,'backgroundcolor',[1 1 1] ,...
    'string', {'micrometer','nanometer'} ,'tag', 'lambda_popup', 'callback', @lambda_popupmenu_callback, 'parent' , field_panel );

% --------------------------------------------------------- Medium Parameters ---
medium_panel = uipanel('units', 'normalized', 'position' , [.305 .005 .3 .99] ,...  %% draw panel
    'parent' , parameters_panel,    'tag' , 'medium_panel'  , 'TitlePosition', 'righttop', ...
    'title' , 'Medium Parameters');

tex=axes('parent',medium_panel,'tag','tex','Position', [.005 .98 .0001 .0001],'XTick',[],'YTick',[]); %define the axes space
text( 1,1, '\bfn=n_0+n_1exp(-\alphar^m/w^m)', 'parent' , tex ,'FontSize',12);

uicontrol ('style', 'popupmenu', 'units', 'normalized', 'position', [0.005 0.87 0.5 0.07],...
    'HorizontalAlignment','left', 'FontWeight', 'normal','fontsize' , 12 ,'backgroundcolor',[1 1 1],...
    'string', {'uniform','linear','quadratic','customize'} ,'tag', 'medium_popup', 'callback', @medium_popupmenu_callback, 'parent' , medium_panel );

tex=axes('parent',medium_panel,'tag','tex','Position', [.005 .8 .0001 .0001],'XTick',[],'YTick',[]); %define the axes space
text( 1,1, 'n_0 =', 'parent' , tex ,'FontSize',12);

uicontrol ('style', 'edit', 'units', 'normalized', 'position', [0.15 0.775 0.15 0.06],...
    'HorizontalAlignment','center', 'FontWeight', 'normal','fontsize' , 12 ,...
    'string', '1' ,'backgroundcolor',[1 1 1],'tag', 'n0_edit','enable','on', 'parent' , medium_panel );

tex=axes('parent',medium_panel,'tag','tex','Position', [.005 .73 .0001 .0001],'XTick',[],'YTick',[]); %define the axes space
text( 1,1, 'n_1 =', 'parent' , tex ,'FontSize',12);

uicontrol ('style', 'edit', 'units', 'normalized', 'position', [0.15 0.705 0.15 0.06],...
    'HorizontalAlignment','center', 'FontWeight', 'normal','fontsize' , 12 ,...
    'string', '0' ,'backgroundcolor',[1 1 1],'tag', 'n1_edit','enable','off', 'parent' , medium_panel );

tex=axes('parent',medium_panel,'tag','tex','Position', [.35 .74 .0001 .0001],'XTick',[],'YTick',[]); %define the axes space
text( 1,1, '*10^-^4', 'parent' , tex ,'FontSize',12);

tex=axes('parent',medium_panel,'tag','tex','Position', [.005 .66 .0001 .0001],'XTick',[],'YTick',[]); %define the axes space
text( 1,1, '\alpha =', 'parent' , tex ,'FontSize',12);

uicontrol ('style', 'edit', 'units', 'normalized', 'position', [0.15 0.635 0.15 0.06],...
    'HorizontalAlignment','center', 'FontWeight', 'normal','fontsize' , 12 ,...
    'string', '0.5' ,'backgroundcolor',[1 1 1],'tag', 'alpha_edit','enable','off', 'parent' , medium_panel );

tex=axes('parent',medium_panel,'tag','tex','Position', [.005 .59 .0001 .0001],'XTick',[],'YTick',[]); %define the axes space
text( 1,1, 'm =', 'parent' , tex ,'FontSize',12);

uicontrol ('style', 'edit', 'units', 'normalized', 'position', [0.15 0.565 0.15 0.06],...
    'HorizontalAlignment','center', 'FontWeight', 'normal','fontsize' , 12 ,...
    'string', '0' ,'backgroundcolor',[1 1 1],'tag', 'm_edit','enable','off', 'parent' , medium_panel );

uicontrol ('style', 'text', 'units', 'normalized', 'position', [0.001 0.47 0.8 0.07],...
    'HorizontalAlignment','left', 'FontWeight', 'bold','fontsize' , 12 ,...
    'string', 'max value of r = 1000*(wave length)'  ,'tag', 'r_max_text', 'parent' , medium_panel );

uicontrol ('style', 'text', 'units', 'normalized', 'position', [0.001 0.415 0.8 0.07],...
    'HorizontalAlignment','left', 'FontWeight', 'bold','fontsize' , 12 ,...
    'string', 'size of propagation step:'  ,'tag', 'dz_size_text', 'parent' , medium_panel );

uicontrol ('style', 'text', 'units', 'normalized', 'position', [0.05 0.365 0.15 0.07],...
    'HorizontalAlignment','left', 'FontWeight', 'bold','fontsize' , 12 ,...
    'string', 'dz ='  ,'tag', 'dz_text', 'parent' , medium_panel );

 uicontrol ('style', 'edit', 'units', 'normalized', 'position', [0.16 0.38 0.15 0.06],...
     'HorizontalAlignment','center', 'FontWeight', 'normal','fontsize' , 12 ,...
     'string', '5' ,'backgroundcolor',[1 1 1],'tag', 'dz_edit','enable','on', 'parent' , medium_panel );

 uicontrol ('style', 'text', 'units', 'normalized', 'position', [0.35 0.365 0.3 0.07],...
     'HorizontalAlignment','left', 'FontWeight', 'bold','fontsize' , 12 ,...
     'string', 'micrometer'  ,'tag', 'dz_unit_text', 'parent' , medium_panel );

 uicontrol ('style', 'text', 'units', 'normalized', 'position', [0.001 0.285 0.8 0.07],...
     'HorizontalAlignment','left', 'FontWeight', 'bold','fontsize' , 12 ,...
     'string', 'propagation distance:'  ,'tag', 'z_distance_text', 'parent' , medium_panel );

 uicontrol ('style', 'text', 'units', 'normalized', 'position', [0.05 0.235 0.15 0.07],...
     'HorizontalAlignment','left', 'FontWeight', 'bold','fontsize' , 12 ,...
     'string', 'z ='  ,'tag', 'z_text', 'parent' , medium_panel );

uicontrol ('style', 'edit', 'units', 'normalized', 'position', [0.16 0.25 0.15 0.06],...
    'HorizontalAlignment','center', 'FontWeight', 'normal','fontsize' , 12 ,...
    'string', '1000' ,'backgroundcolor',[1 1 1],'tag', 'steps_edit','enable','on', 'parent' , medium_panel );

uicontrol ('style', 'popupmenu', 'units', 'normalized', 'position', [0.35 0.237 0.28 0.07],...
    'HorizontalAlignment','left', 'FontWeight', 'normal','fontsize' , 12 ,'backgroundcolor',[1 1 1],...
    'string', {'dz steps','micrometer'} ,'tag', 'steps_popup', 'callback', @steps_popupmenu_callback, 'parent' , medium_panel );

%
%% ------------------------------------ prepare gui_handles
h = guihandles;                                                                            % get handles of uicontrols
h.fig = fig;                                                                                % figures handle
h.object = [];                                                                                % instrument object
h.back_color = back_color;                                                 % figure's background color
h.edit_box_back_color = edit_box_back_color;   % edit box background color
h.start_itr = 1;
h.Output = [];

h.dir_path = fileparts(mfilename('fullpath')); % path of the gui mfile
%
%% ------------------------------------ show figure
% update guidata
guidata ( fig, h );

% show figure
hide_object(parameters_panel);
show_object(graph_panel);
movegui(fig,'center');
set(gcf,'Visible','on');

% focus on uicontrol
focus_on(h.start_pushbutton);
end
%
%% ============================================== focus_on uicontrol ===
function focus_on ( hObject )
% focus on uicontrol only if figure is visible
h = guidata (hObject);
visible_fig = findobj(h.fig,'flat','Visible','on');
if ~isempty(visible_fig)
    uicontrol(hObject);
    drawnow expose;
end
end
%
%% ============================================== hide_object ===
function hide_object (hObject)
set ( hObject, 'Visible', 'off' );
end
%
%% ============================================== show_object ===
function show_object (hObject)
set ( hObject, 'Visible', 'on' );
drawnow expose;
end
%
%% ============================================== callbacks ===
%
%% ============================================== exit_pushbutton_callback ===
function exit_pushbutton_callback (hObject,eventdata)
e = exist ('data_bpm_tmp.mat','file');
if e==2
    delete('data_bpm_tmp.mat');
end
close;                                      % close GUI and all open ports
                                            % (see close_figure_callback)
end
%
%% ============================================== close_figure_callback ===
function close_figure_callback (hObject,eventdata)
stop_pushbutton_callback (hObject,eventdata);
end
%
%% ============================================== stop_pushbutton_callback ===
function stop_pushbutton_callback (hObject,eventdata)
h = guidata(hObject);
set(findobj('tag','stop_pushbutton'),'enable','off');
set(findobj('tag','start_pushbutton'),'enable','on','BackGroundColor',h.back_color,'string','start');
set(findobj('tag','setparameters_pushbutton'),'enable','on');
set(findobj('tag','graphs_pushbutton'),'enable','on');
end
%
%% ============================================== start_pushbutton_callback ===
function flag = start_pushbutton_callback (hObject,flag)
h = guidata(hObject);
get(findobj('tag','start_pushbutton'),'string');

set(findobj('tag','stop_pushbutton'),'enable','on');
set(findobj('tag','setparameters_pushbutton'),'enable','off');
set(findobj('tag','graphs_pushbutton'),'enable','off');
set(findobj('tag','start_pushbutton'),'enable','off','string','running','BackGroundColor',[0 1 0]);
% --------------------------------------------------------------------------

[A w Lam n0 deltaN alpha m dz z]=get_data_from_gui(h);
% Propagation of gaussian using BPM

NR          =     2^13;                % Numerical Resolution

%************Propagating Parameters ********%
steps=round(z/dz);                                          % max value of z-axis
%************Medium Parameters **********%
x_max=1000*Lam;                                        % maximum value of x-axis
%........................................................................%
dx = 2*x_max/(NR-1);
x = -x_max : dx : x_max;
dkx = pi/(x_max)*(NR-1)/NR;
kx = ( (-NR/2):(NR/2-1) )*dkx;
k0  = 2*pi/Lam;
%........................................................................%

%************Inpout Field  Parameters****************%

fz=A*exp(-0.5*(abs(x)/w).^2);        % initial field distribution

N=deltaN*exp(-alpha*(abs(x)/w).^m);

    start_itr = 1;
    Output = zeros(NR,steps);
    Output(:,1) = abs(fz).^2;

%%%%%%%%%%%%%%%%%%%%%%%%
%************PROPAGATION *****************%
%%%%%%%%%%%%%%%%%%%%%%%%
    Fz=ifftshift(ifft(ifftshift(fz)));
    Fz_d=Fz.*exp(1i*(kx.^2)./(2*k0*(n0+N))*dz/2);
    fz=fftshift(fft(fftshift(Fz_d)));

for itr=start_itr:steps-1

    fz_N=fz.*exp(-1i*N*k0*dz);
    Fz=ifftshift(ifft(ifftshift(fz_N)));
    D=exp(1i*(kx.^2)./(2*k0*(n0+N))*dz);
    Fz_D=Fz.*D;
    fz=fftshift(fft(fftshift(Fz_D)));

    plot(h.ax1,x,abs(fz).^2);
    axis(h.ax1,[-0.5*10^-4 0.5*10^-4 0 4*A]);
    set(get(h.ax1,'XLabel'),'String','\bfr [m]');
    set(get(h.ax1,'YLabel'),'String','\bfI [r]');
    set(get(h.ax1,'Title'),'String','\bfIntensity');

    Output(:,itr)=abs(fz).^2;

    imagesc((0:itr)*dz,x,Output(:,1:itr),'parent',h.ax2);
    colormap('hot');
    ylim(h.ax2,[-0.5*10^-4 0.5*10^-4]);
    set(get(h.ax2,'XLabel'),'String','\bfz [m]');
    set(get(h.ax2,'YLabel'),'String','\bfr [m]');
    set(get(h.ax2,'Title'),'String','\bfIntensity Map');
    drawnow;

    val1 = get(findobj('tag','stop_pushbutton'),'enable');
    if strcmp(val1,'off')
        break
    end

end

Fz=ifftshift(ifft(ifftshift(fz)));
Fz_d=Fz.*exp(1i*(kx.^2)./(2*k0*(n0+N))*dz/2);
fz=fftshift(fft(fftshift(Fz_d)));

Output(:,itr+1)=abs(fz).^2;
save('data_bpm_tmp.mat','Output');

plot(h.ax1,x,abs(fz).^2);
axis(h.ax1,[-0.5*10^-4 0.5*10^-4 0 4*A]);
set(get(h.ax1,'XLabel'),'String','\bfr [m]');
set(get(h.ax1,'YLabel'),'String','\bfI [r]');
set(get(h.ax1,'Title'),'String','\bfIntensity');

imagesc((0:itr)*dz,x,Output(:,1:itr+1),'parent',h.ax2);
colormap('hot');
ylim(h.ax2,[-0.5*10^-4 0.5*10^-4]);
set(get(h.ax2,'XLabel'),'String','\bfz [m]');
set(get(h.ax2,'YLabel'),'String','\bfr [m]');
set(get(h.ax2,'Title'),'String','\bfIntensity Map');

% ---------------------------------------------
    set(findobj('tag','start_pushbutton'),'enable','on','string','running','BackGroundColor',h.back_color,'string','start');
    set(findobj('tag','stop_pushbutton'),'enable','off');
    set(findobj('tag','setparameters_pushbutton'),'enable','on');
    set(findobj('tag','graphs_pushbutton'),'enable','on');
end
%
%% ============================================== save_pushbutton_callback ===
function save_pushbutton_callback (hObject,eventdata)
h = guidata(hObject);
e = exist ('data_bpm_tmp.mat','file');
if e==2
    load('data_bpm_tmp.mat');
else
    Output = [];
end
uisave('Output','BPM Simulation.mat');
end
%
%% ============================================== graphs_pushbutton_callback ===
function graphs_pushbutton_callback (hObject,eventdata)
h=guidata(hObject);
hide_object(h.p_panel);
show_object(h.g_panel);
end
%
%% ============================================== setparameters_pushbutton_callback ===
function setparameters_pushbutton_callback (hObject,eventdata)
h=guidata(hObject);
hide_object(h.g_panel);
show_object(h.p_panel);

% focus on uicontrol
focus_on(h.start_pushbutton);
end
%
%% ============================================== lambda_popupmenu_callback ===
function lambda_popupmenu_callback (hObject,eventdata)
h=guidata(hObject);
val=get(h.lambda_popup,'value');
lambda_edit=findobj('tag','lambda_edit');
lambda=get(lambda_edit,'string');

switch val
    case 1
        new_lambda=num2str(str2double(lambda)/1000);
        set(lambda_edit,'string',new_lambda);
    case 2
        new_lambda=num2str(str2double(lambda)*1000);
        set(lambda_edit,'string',new_lambda);
end
end
%
%% ============================================== medium_popupmenu_callback ===
function medium_popupmenu_callback (hObject,eventdata)
h=guidata(hObject);
val=get(h.medium_popup,'value');
n0=findobj('tag','n0_edit');
n1=findobj('tag','n1_edit');
alpha=findobj('tag','alpha_edit');
m=findobj('tag','m_edit');

switch val
    case 1
        set(n0,'enable','on');
        set(n1,'enable','off','string','0');
        set(alpha,'enable','off','string','0.5');
        set(m,'enable','off','string','0');
    case 2
        set(n0,'enable','on');
        set(n1,'enable','on','string','10');
        set(alpha,'enable','on');
        set(m,'enable','off','string','1');
    case 3
        set(n0,'enable','on');
        set(n1,'enable','on','string','10');
        set(alpha,'enable','on');
        set(m,'enable','off','string','2');
    case 4
        set(n0,'enable','on');
        set(n1,'enable','on');
        set(alpha,'enable','on');
        set(m,'enable','on');
end

end
%
%% ============================================== steps_popupmenu_callback ===
function steps_popupmenu_callback (hObject,eventdata)
h=guidata(hObject);
val=get(h.steps_popup,'value');
steps_edit=findobj('tag','steps_edit');
steps=get(steps_edit,'string');
dz=str2double(get(findobj('tag','dz_edit'),'string'));

switch val
    case 1
        new_steps=num2str(round(str2double(steps)/dz));
        set(steps_edit,'string',new_steps);
    case 2
        new_steps=num2str(round(str2double(steps)*dz));
        set(steps_edit,'string',new_steps);
end

end
%
%% ============================================== get_data_from_gui ===
function [A w lambda n0 n1 alpha m dz z]=get_data_from_gui(h)
A_text=get(h.A_edit,'string');
A=str2double(A_text);

w_text=get(h.w_edit,'string');
w=str2double(w_text);
w=w*10^(-6);

lambda_text=get(h.lambda_edit,'string');
lambda=str2double(lambda_text);
w_val=get(h.lambda_popup,'value');
switch w_val
    case 1
        lambda=lambda*10^(-6);
    case 2
        lambda=lambda*10^(-9);
end

n0_text=get(h.n0_edit,'string');
n0=str2double(n0_text);

n1_text=get(h.n1_edit,'string');
n1=str2double(n1_text);
n1=n1*10^(-4);

alpha_text=get(h.alpha_edit,'string');
alpha=str2double(alpha_text);

m_text=get(h.m_edit,'string');
m=str2double(m_text);

dz_text=get(h.dz_edit,'string');
dz=str2double(dz_text);
dz=dz*10^(-6);

z_text=get(h.steps_edit,'string');
z=str2double(z_text);
z_val=get(h.steps_popup,'value');
switch z_val
    case 1
        z=z*dz;
    case 2
        z=z*10^(-6);
end
end

author: Sable

repo: mcbench-benchmarks


text:
function linplus_test22 ( )

%*****************************************************************************80
%
%% TEST22 tests R8GB_FA, R8GB_SL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = m;
  ml = 1;
  mu = 2;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST22\n' );
  fprintf ( 1, '  For a general banded matrix,\n' );
  fprintf ( 1, '  R8GB_FA computes the PLU factors.\n' );
  fprintf ( 1, '  R8GB_SL solves a factored linear system.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of matrix rows M    = %d\n', m );
  fprintf ( 1, '  Number of matrix columns N = %d\n', n );
  fprintf ( 1, '  Lower bandwidth ML         = %d\n', ml );
  fprintf ( 1, '  Upper bandwidth MU         = %d\n', mu );
%
%  Set the matrix.
%
  [ a, seed ] = r8gb_random ( m, n, ml, mu, seed );

  r8gb_print ( m, n, ml, mu, a, '  The banded matrix:' );
%
%  Set the desired solution.
%
  x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
  b = r8gb_mxv ( m, n, ml, mu, a, x );
%
%  Factor the matrix.
%
  [ alu, pivot, info ] = r8gb_fa ( n, ml, mu, a );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST022 - Fatal error!\n' );
    fprintf ( 1, '  R8GB_FA declares the matrix is singular!\n' );
    fprintf ( 1, '  The value of INFO is %d\n', info );
    return;
  end
%
%  Solve the linear system.
%
  job = 0;
  x = r8gb_sl ( n, ml, mu, alu, pivot, b, job );

  r8vec_print ( n, x, '  Solution to A*x=b:' );
%
%  Set the desired solution.
%
  x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
  job = 1;
  b = r8gb_ml ( n, ml, mu, alu, pivot, x, job );

  r8vec_print ( n, b, '  Right hand side of transposed system:' );
%
%  Solve the linear system.
%
  job = 1;
  x = r8gb_sl ( n, ml, mu, alu, pivot, b, job );

  r8vec_print ( n, x, '  Solution to A''x=b:' );

  return
end


text:
function varargout = horzcat(varargin)
%HORZCAT   Horizontal concatenation of SPHEREFUNV objects.
%   This is not allowed and returns an error.  This function exists so that the
%   error message is meaningful to a SPHEREFUN user.

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

error('SPHEREFUN:SPHEREFUNV:horzcat:notSupported', ...
    'Horizontal concatenation of SPHEREFUNV objects is not supported.')

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function lines3d(varargin)
%LINES3D Description of functions operating on 3D lines.
%
%   A 3D Line is represented by a 1-by-6 row vector containing a 3D point
%   (its origin) and a 3D vector (its direction):
%   LINE = [X0 Y0 Z0 DX DY DZ];
%
%   See also
%   createLine3d, distancePointLine3d, isPointOnLine3d, linePosition3d
%   intersectLinePlane, distanceLines3d, parallelLine3d, projPointOnLine3d
%   clipLine3d, fitLine3d, drawLine3d, transformLine3d
%   edgeToLine3d, lineToEdge3d
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2008-10-13, using Matlab 7.4.0.287 (R2007a)
% Copyright 2008-2022 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas

help('lines3d');

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
classdef NomadDirectionTypeEnum < matlab.mixin.SetGet
    %NomadDirectionTypeEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        Ortho1('ortho 1', 'Ortho 1')
        Ortho2('ortho 2', 'Ortho 2')
        OrthoN1('ortho n+1', 'Ortho N+1')
        OrthoN1Quad('ortho n+1 quad', 'Ortho N+1 Quad')
        OrthoN1Neg('ortho n+1 neg', 'Ortho N+1 Neg')
        Ortho('ortho', 'Ortho')
        Ortho2N('ortho 2n', 'Ortho 2N')

        LT1('lt 1', 'LT 1')
        LT2('lt 2', 'LT 2')
        LT2N('lt 2n', 'LT 2N')
        LTN1('lt n+1', 'LT N+1')

        GPSBin('gps binary', 'GPS Binary')
        GPS2NStatic('gps 2n static', 'GPS 2N Static')
        GPS2NRand('gps 2n rand', 'GPS 2N Random')
        GPSN1StaticUniform('gps n+1 static uniform', 'GPS N+1 Static Uniform')
        GPSN1Static('gps n+1 static', 'GPS N+1 Static')
        GPSN1RandUniform('gps n+1 rand uniform', 'GPS N+1 Random Uniform')
        GPSN1Rand('gps n+1 rand', 'GPS N+1 Random')
        GPS1Static('gps 1 static', 'GPS 1 Static')
    end

    properties
        optionStr char = ''
        name char = '';
    end

    methods
        function obj = NomadDirectionTypeEnum(optionStr, name)
            obj.optionStr = optionStr;
            obj.name = name;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('NomadDirectionTypeEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('NomadDirectionTypeEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('NomadDirectionTypeEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%io_readlcmcoord_getBackground.m
%Jamie Near, McGill University 2014.
%
% USAGE:
% out = io_readlcmcoord_getBackground(filename,part)
%
% DESCRIPTION:
% Reads a LCModel .coord file and extracts the desired part.
%
% INPUTS:
% filename   = filename of the LCModel .coord file.
% part       = Which part of the .coord file to extract - 'bg' extracts the
%             LCModel baseline signal, 'sp' extracts the spectrum, and
%             'fit' extracts the fit.
%
% OUTPUTS:
% out        = Desired background component in simplified FID-A structure format.

function [out]=io_readlcmcoord_getBackground(filename,part)

fid=fopen(filename);
linenum=1;
line=fgets(fid)

%READ PPM AXIS
ppmstr='points on ppm-axis ='
ppm_index=findstr(line,ppmstr);
while isempty(ppm_index)
    line=fgets(fid);
    ppm_index=findstr(line,ppmstr);
end

Npts=str2num(line(2:5))

line=fgets(fid);
while linenum<=(Npts/10)
    [A,count, errmsg, nextindex] = sscanf(line, '%f', inf);
    A
    ppm((linenum-1)*10+1:linenum*10,1)=A;
    linenum=linenum+1;
    line=fgets(fid);
end

%GET BACKGROUND

linenum=1;
switch part
    case 'bg'
        bgstr='NY background values follow'
    case 'sp'
        bgstr='NY phased data points follow'
    case 'fit'
        bgstr='NY points of the fit to the data follow'
    otherwise
        error('ERROR:  part not found');
end
bg_index=findstr(line,bgstr);
while isempty(bg_index)
    line=fgets(fid);
    bg_index=findstr(line,bgstr);
end

line=fgets(fid);
while linenum<=(Npts/10)
    [A,count, errmsg, nextindex] = sscanf(line, '%f', inf);
    A;
    bg((linenum-1)*10+1:linenum*10,1)=A;
    linenum=linenum+1;
    line=fgets(fid);
end

out.ppm=ppm;
out.specs=bg;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
classdef RemeshingTests < handle & matlab.unittest.TestCase

    methods (Test, TestTags = {'Remesh'})

        function testRemeshMesh(obj)
            mC = obj.createCoarseMesh();
            mF = mC.remesh(1);
            s = load('test_RemeshMesh');
            err(1) = norm(mF.coord(:)  - s.meshFine.coord(:));
            err(2) = norm(mF.connec(:) - s.meshFine.connec(:));
            tol = 1e-6;
            obj.verifyLessThanOrEqual(norm(err), tol)
        end

        function testRemeshP1ContinousFunction(obj)
            m = obj.createCoarseMesh();
            f = obj.createP1ContinousFunction(m);
            for i = 1:2
                mF = m.remesh(1);
                f = f.refine(m,mF);
                m    = mF;
            end
            s = load('test_RemeshP1Function');
            err = norm(s.fValues(:) - f.fValues(:));
            tol = 1e-6;
            obj.verifyLessThanOrEqual(err, tol)
        end

        function testRemeshP1DiscontinousFunction(obj)
            m = obj.createCoarseDiscontinousMesh();
            f = obj.createP1DiscontinousFunction(m);
            for i = 1:2
                mF = m.remesh(1); % mF is CONTINUOUS, m is DISCONTINUOUS
                f = f.refine(m,mF); %fNew = fOld.refine(mF); -> fOld has m, fNew has mF
                m = mF.createDiscontinuousMesh(); % care: discont/cont -> interpolation through continuous
                % ideally: work only with discontinuous meshes
            end
            s = load('test_RemeshP1DiscFunction');
            err = norm(s.fValues(:) - f.fValues(:));
            tol = 1e-6;
            obj.verifyLessThanOrEqual(err, tol)
        end

    end

    methods (Access = private)

        function m = createCoarseMesh(obj)
            s.coord = [1 0; 0 1; 0 0; 1 1];
            s.connec = [3 1 2; 1 4 2];
            m = Mesh(s);
        end

        function mD = createCoarseDiscontinousMesh(obj)
            m   = obj.createCoarseMesh();
            mD = m.createDiscontinuousMesh();
        end

        function fC = createP1ContinousFunction(obj,m)
            f         = obj.createFunctionToRemesh();
            s.mesh    = m;
            s.fValues = f(m.coord);
            fC        = P1Function(s);
        end

        function fC = createP1DiscontinousFunction(obj,m)
            f         = obj.createFunctionToRemesh();
            s.fValues = f(m.computeBaricenter()');
            s.mesh    = m;
            f0 = P0Function(s);
            fC = f0.project('P1D');
        end

        function f = createFunctionToRemesh(obj)
            f = @(x) x(:,1).^2+x(:,2).^2;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function imgBlend = sc_poisson_blend(imgTrg, imgSrc, holeMask)
% SC_POISSON_BLEND: Blend the source and target image using a simple
% discrete Poisson solver
%
% Input:
%   - imgTrg: original image with hole
%   - imgSrc: synthesized image
%   - holeMask: specify the hole pixels
% Output:
%   - imgBlend: blended image

[imgH, imgW, nCh] = size(imgSrc);

% Initialize the reconstructed image
imgRecon = zeros(imgH, imgW, nCh);

% Independently process each channel
for ch = 1: nCh
    % Source and target image
    imgS = imgSrc(:,:,ch);
    imgT = imgTrg(:,:,ch);

    % prepare discrete Poisson equation
    [A, b] = solvePoisson(holeMask, imgS, imgT);

    % solve Poisson equation
    x = A\b;
    imgRecon(:,:,ch) = reshape(x, [imgH, imgW]);
end

% Combined with the known region in the target
holeMaskC = cat(3, holeMask, holeMask, holeMask);
imgBlend = holeMaskC.*imgRecon + ~holeMaskC.*imgTrg;

end

function [A, b] = solvePoisson(holeMask, imgS, imgT)

% Prepare the linear system of equations for Poisson blending

[imgH, imgW] = size(holeMask);
N = imgH*imgW;

% Number of unknown variables
numUnknownPix = sum(holeMask(:));

maxNumInd = 8*numUnknownPix; % Max number of indices for constructing the sparse matrix
maxNumEqn = 4*numUnknownPix; % Max number of equations (4-neighbor)

% 4-neighbors: dx and dy
dx = [1, 0, -1,  0];    dy = [0, 1,  0, -1];

% Initialize (I, J, S), for sparse matrix A where A(I(k), J(k)) = S(k)
I = zeros(maxNumInd, 1);
J = zeros(maxNumInd, 1);
S = zeros(maxNumInd, 1);

% Initialize b
b = zeros(maxNumEqn, 1);

% Precompute unkonwn pixel position
[pi, pj] = find(holeMask == 1);
pind = sub2ind([imgH, imgW], pi, pj);

% Precompute the 4-neighbor of the unkonwn pixel positions
qi = bsxfun(@plus, pi, dy);
qj = bsxfun(@plus, pj, dx);

% Handling cases at image borders
validN = (qi >= 1) & (qi <= imgH) & (qj >= 1) & (qj <= imgW);
qind = zeros(size(validN));
qind(validN) = sub2ind([imgH, imgW], qi(validN), qj(validN));

c = 1; % index counter for matrix A
e = 1; % equation counter

% Set up the matrix A and the vector b
for k = 1: numUnknownPix
    pind_cur = pind(k);
    for n = 1: 4 % 4-neighbor
        if(validN(k, n)) % if the neighbor pixel q lies in the image
            qind_cur = qind(k, n);
            if(holeMask(qind_cur))
                % A(e, pind_cur) = 1;
                I(c) = e;   J(c) = pind_cur;    S(c) = 1;   c = c + 1;
                % A(e, qind_cur) = 1;
                I(c) = e;   J(c) = qind_cur;    S(c) = -1;  c = c + 1;

                % gradient constraint
                b(e) = imgS(pind_cur) - imgS(qind_cur);
                e = e + 1;
            else
                % A(e, pind_cur) = 1;
                I(c) = e;   J(c) = pind_cur;    S(c) = 1;   c = c + 1;

                % boundary constraint
                b(e) = imgS(pind_cur) - imgS(qind_cur) + imgT(qind_cur);
                e = e + 1;
            end
        end
    end
end

% Clean up unused entries
nEqn = e - 1;
b = b(1:e-1);

nInd = c - 1;
I = I(1:nInd);  J = J(1:nInd);  S = S(1:nInd);

% Construct the sparse matrix A
A = sparse(I, J, S, nEqn, N);

end

repo: StructCompletion


text:
% demo for mixed-effects analysis (2-levels hierachical model)
% This script simulates and analyses data under two different hierarchical
% models, namely (i) H0: group-mean=0, and (ii) H1: group-mean~=0.
% The effect of H1 and H0 onto both the group-mean parameters estimates
% (given H1) and the log-Bayes factors (comparisons of H1 vs H0) is
% assessed using Monte-Carlo simulations.
% The trick for performing a mixed-effects analysis using the structure of
% the generative model of VBA is to use hidden states, where "initial
% conditions" (x_0) serve as the group mean, x_1 are subject-dependant
% effects, and the precision (alpha) of the transition density
% p(x_1|x_0,alpha) captures the between-subject variance. This trick
% however, can be used here because we deal with a static model (a
% within-subject GLM).

clear variables
close all

% Choose basic settings for simulations
ns = 32; % number of subjects
n = 256; % number of observations per subject
f_fname = @g_RFX;
g_fname = @g_RFX;
alpha = 1e1;
sigma = 1e0;
theta = [];
phi = [];
x0 = zeros(ns,1);

X2 = [ones(ns,1),zeros(ns,ns-1)]; % 2d-level design matrix
u = [];

% Build options structure for temporal integration of SDE
inF.ns = ns;
inF.X = X2;
options.inF = inF;

% Build priors for model inversion
priors.muX0 = zeros(ns,1);
priors.SigmaX0 = 0*eye(ns);
priors.a_alpha = 1e0;
priors.b_alpha = 1e0;
priors.a_sigma = 1e0;
priors.b_sigma = 1e0;

% Build options and dim stuctures for model inversion
options.priors = priors;
options.DisplayWin = 0;
dim.n_theta = 0;
dim.n_phi = 0;
dim.n = ns;

Nmcmc = 32;
p = cell(2,2,Nmcmc);
o = cell(2,2,Nmcmc);
F = zeros(2,2,Nmcmc);
X0 = zeros(2,Nmcmc);
for ii=1:Nmcmc
    % simulate data with and without 2nd-level effect
    X1 = randn(n,ns); % 1st-level design matrix
    options.inG.X = X1;
    for i=1:2
        x0(1) = 2-i;
        [y,x,x0,eta,e] = VBA_simulate (1,f_fname,g_fname,theta,phi,u,alpha,sigma,options,x0);
        % Invert model with and without 2nd-level effect
        for j=1:2
            options.priors.SigmaX0(1,1) = 2-j; % group mean effect
            [p{i,j,ii},o{i,j,ii}] = VBA_NLStateSpaceModel(y,u,f_fname,g_fname,dim,options);
            F(i,j,ii) = o{i,j,ii}.F;
            if j==1
                % extract estimaetd group-mean
                X0(i,ii) = p{i,j,ii}.muX0(1);
            end
        end
    end
end

hf = figure('color',[1 1 1]);
pos = get(hf,'position');
set(hf,'position',pos.*[1 1 1.5 1]);
models = {'H1','H0'};

ha = subplot(1,2,1,'parent',hf);
mx0 = mean(X0,2);
vx0 = var(X0,[],2)./Nmcmc;
plotUncertainTimeSeries(mx0,vx0,[],ha);
set(ha,'xlim',[0,3],'xtick',[1,2],'xticklabels',models)
xlabel(ha,'type of simulated data')
ylabel(ha,['E[X0|y,H1]'])
title(ha,'estimated group effect (under H1)')
box(ha,'off')

ha = subplot(1,2,2,'parent',hf);
dF = F(:,1,:) - F(:,2,:);
mdF = mean(dF,3);
vdF = var(dF,[],3)./Nmcmc;
plotUncertainTimeSeries(mdF,vdF,[],ha);
set(ha,'xlim',[0,3],'xtick',[1,2],'xticklabels',models)
xlabel(ha,'type of simulated data')
ylabel(ha,['log p(y|',models{1},') - log p(y|',models{2},')'])
title(ha,'evidence for a group effect')
box(ha,'off')

VBA_getSubplots ();

author: MBB-team

repo: VBA-toolbox


text:
%VideoCamera Abstract class to read from local video camera
%
% A concrete subclass of ImageSource that acquires images from a local
% camera using the MATLAB Image Acquisition Toolbox (imaq).  This Toolbox
% provides a multiplatform interface to a range of cameras, and this
% class provides a simple wrapper.
%
% This class is not intended to be used directly, instead use the factory
% method Video which will return an instance of this class if the Image
% Acquisition Toolbox is installed, for example
%
%         vid = VideoCamera();
%
% Methods::
% grab    Aquire and return the next image
% size    Size of image
% close   Close the image source
% char    Convert the object parameters to human readable string
%
% See also VideoCamera, ImageSource, AxisWebCamera, Movie.
function vid = VideoCamera(varargin)

    opt.simple = false;

    [opt,args] = tb_optparse(opt, varargin);

    % this function looks like a class, the only way to implement the
    % factory design pattern

    if exist('imaqhwinfo', 'file') && ~opt.simple
        % we have the Mathworks Image Acquisition Toolbox
        if length(args) == 1 && strcmp(args{1}, '?')
            VideoCamera_IAT.list();
        else
            vid = VideoCamera_IAT(args{:});
        end
    elseif exist('framegrabber', 'file') == 3
        % we have the MVTB framegrabber MEX interface, for either
        % MacOS, Linux or Windows
        if length(args) == 1 && strcmp(args{1}, '?')
            VideoCamera_fg.list();
        else
            vid = VideoCamera_fg(args{:});
        end
    else
        error('no video capture capability on this computer');
    end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function varargout = normalization(PopObj)

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

global normalization_str
if strcmp(normalization_str, 'normalize')
    [N, ~] = size(PopObj);
    a = intercepts(PopObj);
    PopObj = (PopObj - repmat(min(PopObj, [], 1), size(PopObj, 1), 1)) ./ repmat(a,N,1);
end
varargout{1} = PopObj;
if nargout == 2
    varargout{2} = a;
end
end

function a = intercepts(PopObj)
[N, M] = size(PopObj);
%% Find the extreme points
[~, Choosed(1:M)] = min(PopObj, [], 1);
L2NormABO = zeros(N, M);
for i = 1 : M
    L2NormABO(:, i) = sum(PopObj(:, [1: i - 1, i + 1: M]) .^ 2, 2);
end
[~, Choosed(M + 1: 2 * M)] = min(L2NormABO, [], 1);
[~, Extreme] = max(PopObj(Choosed, :), [], 1);
Extreme = unique(Choosed(Extreme));
%% Calculate the intercepts
if length(Extreme) < M
    a = max(PopObj, [], 1);
else
    lastwarn('');
    Hyperplane = mldivide(PopObj(Extreme,:), ones(M, 1));
    [~, msgid] = lastwarn();
    if strcmp(msgid, 'MATLAB:nearlySingularMatrix')
        % error('error in normalize');
    end
    a = 1 ./ Hyperplane';
end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function varargout = size(T,dim)
% overloads size

s = size(T.M);
s = s(T.rank+1:end);

if nargin == 1
  if length(s) <= 1
    s = [s,ones(1,2-length(s))];
  end
elseif dim > numel(s)
  s = 1;
else
  s = s(dim);
end

if nargout > 1
  varargout = num2cell(s);
else
  varargout{1} = s;
end

author: mtex-toolbox

repo: mtex

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


text:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% file: test_prop_diff_ci.m
%
% use: run "test_prop_diff_ci" and then examine "test_prop_diff_ci_results.csv" in the work folder and
%      compare it to the "test_prop_diff_ci_results_reference.csv" distributed with the code.
%      Failure to run or significant differences in results may indicate installation problems.
%
% change log:
%   030114 tdr created from test_prop_ci.m

function test_prop_diff_ci

methods = [1 3 4 5 6];
alphas = [1e-4 1e-2 0.9];
ns = [1 30 1e5];

fid = fopen('test_prop_diff_ci_results.csv','w');
fprintf(fid,'%s\n','delta_p_hat, Lower CI Bound, Upper CI Bound, x1, n1, x2, n2, Desired alpha, Method, Length, Lower Tail, Upper Tail, Actual alpha, Delta alpha, Run Time');

warning off;

for method = methods
   for alpha = alphas
      for n1 = ns
         if n1 < 5
            x1s = [0];
         else
            x1s = [0 1 round(n1/3)];
         end;
         for n2 = ns
            if n2 < 5
                x2s = [n2];
            else
                x2s = [0 1 round(n2/2) n2];
            end;
            for x1 = x1s
                for x2 = x2s
                    ci = prop_diff_ci(x1,n1,x2,n2,alpha,method,1);
                    fprintf(fid,'%8.6g, %8.6g, %8.6g, %d, %d, %d, %d, %8.6g, %d, %8.6g, %8.6g, %8.6g, %8.6g, %8.6g, %8.2g\n',ci);
                end;
			end;
		 end;
     end;
   end;
end;

warning on;
fclose(fid);


text:
function vn = normalize(v)
%NORMALIZE Normalize a vector.
%
%   V2 = normalize(V);
%   Returns the normalization of vector V, such that ||V|| = 1. V can be
%   either a row or a column vector.
%
%   When V is a MxN array, normalization is performed for each row of the
%   array.
%
%   See also
%   vectors2d, normalizeVector, vectorNorm
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2004-11-29
% Copyright 2004-2022 INRA - TPV URPOI - BIA IMASTE

dim = size(v);

if dim(1)==1 || dim(2)==1
    vn = v/sqrt(sum(v.*v));
else
    vn = v./repmat(sqrt(sum(v.*v, 2)), [1 dim(2)]);
end


text:
function [reverseStr] = cmdWaitbar(i,nItems,reverseStr,varargin)
	% Puts a txt waitbar into the cmd window, less intrusive than pop-up waitbar.
	% Biafra Ahanonu
	% started: 2014.01.14
	% thanks to:
		% http://www.mathworks.com/matlabcentral/newsreader/view_thread/32291
		% http://stackoverflow.com/questions/11050205/text-progress-bar-in-matlab
	% inputs
			%
	% outputs
			%

	% changelog
			% 2014.02.14 [16:35:55] now is mostly
	% TODO
			% Should reverseStr be made global so function is entirely self-contained? - NO, globals are evil.
			% change so waitbarOn = 0 can short circuit getOptions to save speed execution time
	% example
		% before loop = reverseStr
		% reverseStr = cmdWaitbar(i,nItems,reverseStr,'inputStr','loading hdf5','waitbarOn',options.waitbarOn,'displayEvery',50);

	[reverseStr] = ciapkg.view.cmdWaitbar(i,nItems,reverseStr,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function y = vl_nnreshape(x, shape, varargin)
% VL_NNRESHAPE Feature reshaping
%   Y = VL_NNRESHAPE(X, SHAPE) reshpaes the input data X to have
%   the dimensions specified by SHAPE. X is a SINGLE array of
%   dimension H x W x D x N where (H,W) are the height and width of
%   the map stack, D is the image depth (number of feature channels)
%   and N the number of of images in the stack. SHAPE is a 1 x 3 cell
%   array, the contents of which are passed in order to the MATLAB
%   reshape function. As a consequence, `[]` an be used to specify a
%   dimension which should be computed from the other two. The batch size
%   (the fourth dimension of the input) is left unchanged by this
%   reshaping operation.
%
%   Alternatively, SHAPE can be specified in the "caffe style", as an
%   array, using `0` to indicate that a dimension should be preserved and
%   `-1` to indicate a dimension that should be computed from the others
%   (this serves the same role as [] in the standard MATLAB convention).
%
%   Example:
%       Inputs: X with shape [100 100 3 5] and SHAPE = { 100 3 [] }
%       will produce an output Y with shape [100 3 100 5]
%
%   DZDX = VL_NNRESHAPE(X, SHAPE, DZDY) computes the derivatives of the
%   block projected onto DZDY. DZDX and DZDY have the same dimensions
%   as X and Y respectively.
%
% Copyright (C) 2017 Samuel Albanie and Andrea Vedaldi.
% Licensed under The MIT License [see LICENSE.md for details]

  [~, dzdy] = vl_argparsepos(struct(), varargin) ;

  if isnumeric(shape) % apply caffe style conventions if needed
    shape_ = num2cell(shape) ;
    if numel(shape_) == 2, shape_{3} = [] ; end
    k = find(shape == -1) ; if k, shape_{k} = [] ; end
    k = find(shape == 0) ;
    if k, rep = arrayfun(@(i) {size(x,i)}, k) ; shape_(k) = rep ; end
    shape = shape_ ;
  end

  batchSize = size(x, 4);

  if isempty(dzdy)
    y = reshape(x, shape{1}, shape{2}, shape{3}, batchSize) ;
  else
    y = reshape(dzdy{1}, size(x)) ;
  end

repo: MFT


text:
function data = read_tck(filename)

% READ_TCK reads tractography information from an mrtrix-generated .tck
% file. Requires the matlab functions from mrtrix.

tmp  = read_mrtrix_tracks(filename);
npos = cellfun(@numel, tmp.data)./3;
maxnpos = max(npos);
lines   = zeros(2*maxnpos, numel(npos))+nan;
offset  = 0;
for k = 1:numel(npos)
  indx1 = 1:npos(k)*2;
  indx2 = [1:npos(k) npos(k):-1:1]+offset;
  lines(indx1,k) = indx2;
  offset = max(indx2);
end

data.pos  = cat(1, tmp.data{:});
data.line = lines.';
data.hdr  = rmfield(tmp, 'data');

author: fieldtrip

repo: fieldtrip


text:
function varargout = imregrotate(varargin)
%ROIROTATE Rotate & crop region of interest avoiding 'black corners'.
% 	B = ROIROTATE(A,ROI,ANGLE,METHOD) rotates the region of intererest ROI
% 	found in the image A by ANGLE degrees in a counter-clockwise direction,
% 	using the specified interpolation method and gives the rotated
% 	region as output in B cropped out of the original image.
%
% 	Usefulness
% 	----------
% 	ROIROTATE was designed to avoid the 'black corners' typical of
%   MATLAB's imrotate function, which are apearing from its behaviour
%   of padding by zeros the space outside the image borders prior to
%   the rotation. This is possible if the image to be rotated is
%   part of a bigger image - which is the definition of a 'region
%   of interest' -, because there is enough information to fill
%   the corners.
%
% 	Description
% 	-----------
% 	Apart of the first entry which is based on the IMCROP syntax,
% 	all other are identical to the IMROTATE function.
%
%   ROI is a four-element vector with the form [xmin ymin width height];
%   these values are specified in spatial coordinates.
%
% 	xmin - starting pixel to be cut, x-th of the original image
% 	ymin - same on y-axis
% 	width - number of pixels to be cut minus one (width)
% 	height - same on height
%
% 	If you want to point the 5 pixel wide & 7 pixel heigh NW corner
% 	of an image, replace R with [1,1,4,6].
%
%   Because ROI is specified in terms of spatial coordinates,
%   the WIDTH and HEIGHT of ROI do not always correspond exactly
%   with the size of the output image. For example, suppose ROI
%   is [20 20 40 30], using the default spatial coordinate
%   system. The upper left corner of the specified rectangle is
%   the center of the pixel (20,20) and the lower right corner is
%   the center of the pixel (50,60). The resulting output image
%   is 31-by-41, not 30-by-40, because the output image includes
%   all pixels in the input that are completely or partially
%   enclosed by the rectangle.
%
%         Note for IMCROP:
%         To cut the NW pixel only:
%             NWcornerPixel = imcrop(I, [1,1,0,0]);
%
%         To cut the first top stripe of pixels:
%             topStripe = imcrop(I, [1,1,size(I,2)-1,0]);
%
% 	METHOD is a string that can have one of these values:
%
%        'nearest'  (default) nearest neighbor interpolation
%        'bilinear' bilinear interpolation
%        'bicubic'  bicubic interpolation
%
% 	If you omit the METHOD argument, ROIROTATE uses the default
% 	method of 'nearest'.
%
%   To rotate the image clockwise, specify a negative angle.
%
% 	Class Support
% 	-------------
% 	The input image can be of class uint8, uint16, or double.
%   The output image is of the same class as the input image.
%
% 	Example
% 	-------
%        I = imread('circuit.tif');
%        J = imrotate(imcrop(I,[60 80 100 90]),-21,'bilinear','crop');
%        [K x y w h] = roirotate(I,[60 80 100 90],-21,'bilinear');
%        figure, subplot(1,3,1), imshow(I), title('roi')
%        rectangle('Position',[60 80 100 90],'EdgeColor','red');
%        rectangle('Position',[x y w h],'EdgeColor','green');
%        subplot(1,3,2), imshow(J), title('imrotate')
%        subplot(1,3,3), imshow(K), title('roirotate')
%
% 	See also IMROTATE, IMCROP, IMRESIZE, IMTRANSFORM, TFORMARRAY.

%   Vlad Atanasiu - 06/05/2002
%   Revisions: 3  > 01/11/2008

% check image input data
[A,dim,x,y,width,height,ang,method] = parseInputs(varargin{:});

% calculate padding borders
padX = abs(ceil(height*sin(ang*pi/180)/2));
padY = abs(ceil(width*sin(ang*pi/180)/2));

% padding zeros if extended region of interest coordinates outside image
if x-padX < 1   % left
    A = [zeros(size(A,1),padX-x+1,dim) A];
end
if y-padY < 1   % top
    A = [zeros(padY-y+1,size(A,2),dim);...
            A];
end
if x+width+padX > size(A,2)   % right
    A = [A zeros(size(A,1),x+width+padX-size(A,2),dim)];
end
if y+height+padY > size(A,1)   % bottom
    A = [A;...
            zeros(y+height+padY-size(A,1),size(A,2),dim)];
end

% cut extended region
ROI = imcrop(A,[x-padX y-padY width+padX*2 height+padY*2]);

% rotate extended region
ROI = imrotate(ROI,ang,method,'crop');

% cut region
varargout{1} = imcrop(ROI,[padX padY width height]);

% gives extended region coordinates
if nargout == 5
    varargout{2} = x-padX;
    varargout{3} = y-padY;
    varargout{4} = width+padX*2;
    varargout{5} = height+padY*2;
end

% ------------------------------------------------

function [A,dim,x,y,width,height,ang,method] = parseInputs(varargin)
% Input checkings

% Defaults
method = 'n';

% Check number of input parameters
error(nargchk(3,4,nargin));
switch nargin
case 3,             % regrotate(A,reg,ang,)
    A = varargin{1};
	x = varargin{2}(1);
	y = varargin{2}(2);
	width = varargin{2}(3);
	height = varargin{2}(4);
    ang = varargin{3};
case 4,             % regrotate(A,reg,ang,method)
    A = varargin{1};
	x = varargin{2}(1);
	y = varargin{2}(2);
	width = varargin{2}(3);
	height = varargin{2}(4);
    ang = varargin{3};
    method = varargin{4};
otherwise,
    error('Invalid number of input arguments.');
end

% Check validity of the input parameters
if ischar(method),
    strings = {'nearest','bilinear','bicubic'};
    idx = strmatch(lower(method),strings);
    if isempty(idx),
        error(sprintf('Unknown interpolation method: %s',method));
    elseif length(idx)>1,
        error(sprintf('Ambiguous interpolation method: %s',method));
    else
        method = strings{idx};
    end
else
    error('Interpolation method have to be a string.');
end

if size(varargin{2},2) ~= 4
    error('Region of interest should be defined by 4 arguments.');
end

% transforms the format of A into an image if it is a filename
if (isstr(A))
    A = imread(A);
end

% gets the dimensions of A
dim = ndims(A);
if dim == 2, dim=1; end

% modifies if needed the region of interest coordinates if outside image
if x < 1, x = 1; end
if y < 1, y = 1; end
if x > size(A,2)-1, x = size(A,2)-1; end
if y > size(A,1)-1, y = size(A,1)-1; end
if x+width > size(A,2)-1, width = size(A,2)-x-1; end
if y+height > size(A,1)-1, height = size(A,1)-y-1; end


text:
function pass = test_separableFormat( pref )
% Test code for computing separable format of PDOs.
%
% Alex Townsend, September 2014.

% Constant coefficient, Laplace:
N = chebop2(@(x,y,u) diff(u, 2, 1) + diff(u, 2, 2), [-1 1 -1 1]);
[U, S, V] = chebop2.separableFormat( N );
pass(1)  = checkCoeffs(N.coeffs, U, S, V);

% Constant coefficient, helmholtz:
N = chebop2(@(x,y,u) diff(u, 2, 1) + diff(u, 2, 2) + 100*u, [-1 1 -1 1]);
[U, S, V] = chebop2.separableFormat( N );
pass(2)  = checkCoeffs(N.coeffs, U, S, V);

% Rank 3 PDO:
N = chebop2(@(x,y,u) diff(u, 2, 1) + diff(u, 2, 2) + (10*x + y).*u, [-1 1 -1 1]);
[U, S, V] = chebop2.separableFormat( N );
pass(3)  = checkCoeffs(N.coeffs, U, S, V);

% Highly oscillatory variable coefficient:
N = chebop2(@(x,y,u) diff(u, 2, 1) + diff(u, 2, 2) + cos(100*x).*u, [-1 1 -1 1]);
[U, S, V] = chebop2.separableFormat( N );
pass(4)  = checkCoeffs(N.coeffs, U, S, V);

% Helmholtz with gravity:
N = chebop2(@(x,y,u) diff(u, 2, 1) + diff(u, 2, 2) + (10 + y).*u, [-1 1 -1 1]);
[U, S, V] = chebop2.separableFormat( N );
pass(5)  = checkCoeffs(N.coeffs, U, S, V);

% On non-standard domain:
N = chebop2(@(x,y,u) diff(u,2,1) + diff(u,2,2) + (10+y).*u, ...
    [-1 pi/3 -sqrt(3)/2 1.1]);
[U, S, V] = chebop2.separableFormat( N );
pass(6)  = checkCoeffs(N.coeffs, U, S, V);

% Variable coefficients on higher derivatives:
N = chebop2(@(x,y,u) (1+x.^2).*diff(u,2,1) + diff(u,2,2) + u, ...
    [-1 pi/3 -sqrt(3)/2 1.1]);
[U, S, V] = chebop2.separableFormat( N );
pass(7)  = checkCoeffs(N.coeffs, U, S, V);

% Heat equation:
N = chebop2(@(x,y,u) diff(u,1,1) - diff(u, 2, 2), [-1 1 0 10]);
[U, S, V] = chebop2.separableFormat( N );
pass(8) = checkCoeffs(N.coeffs, U, S, V);

% Wave equation:
N = chebop2(@(x,y,u) diff(u, 1, 2) - x.^2.*diff(u, 2, 2) + y.*u, [-1 1 0 10]);
[U, S, V] = chebop2.separableFormat( N );
pass(9) = checkCoeffs(N.coeffs, U, S, V);

% complex scalar:
N = chebop2(@(x,y,u) 1i*diff(u, 1, 1));
[U, S, V] = chebop2.separableFormat( N );
pass(10) = checkCoeffs(N.coeffs, U, S, V);

% Complex heat equation?
N = chebop2(@(x,y,u) 1i*diff(u,1,1) + (1+1i)*diff(u,1,2));
[U, S, V] = chebop2.separableFormat( N );
pass(11) = checkCoeffs(N.coeffs, U, S, V);

% Time-independent Schrodinger (complex coefficients):
hb = 0.0256;
N = chebop2(@(x,y,u) 1i*hb*diff(u, 1, 1) + hb^2*diff(u, 2, 2) - x.^2.*u, ...
    [-1 1 0 10]);
[U, S, V] = chebop2.separableFormat( N );
pass(12) = checkCoeffs(N.coeffs, U, S, V);

end

function bol = checkCoeffs( A, cellU, S, cellV )
% CheckCoeffs  debugging script.
% Alex Townsend, September 2014.
for jj = 1:size(cellU, 1)
    for kk = 1:size(cellV, 1)
        a = cellU{jj,1} * S(1,1) * cellV{kk,1}.';
        for r = 2:size(cellU,2)
            a = a + cellU{jj,r} * S(r,r) * cellV{kk,r}.';
        end
        store{jj,kk} = a;
        err(jj,kk) = norm( a - A{jj,kk} );
    end
end
if ( norm( err ) < 1e-8 )
    bol = 1;
else
    bol = err;
end
% debug plot:
% sp = 1;
% M = size(A, 1);
% N = size(A, 2);
% for jj = 1:M
%     for kk = 1:N
%         subplot(M,N,sp), plot(store{jj,kk} )
%         sp = sp + 1;
%     end
% end
end


text:
% pop_readegi() - load a EGI EEG file (pop out window if no arguments).
%
% Usage:
%   >> EEG = pop_readegi;             % a window pops up
%   >> EEG = pop_readegi( filename );
%   >> EEG = pop_readegi( filename, datachunks, forceversion, fileloc);
%
% Inputs:
%   filename       - EGI file name
%   datachunks     - desired frame numbers (see readegi() help)
%                    option available from the command line only
%   forceversion   - [integer] force reading a specfic file version
%   fileloc        - [string] channel location file name. Default is
%                    'auto' (autodetection)
%
% Outputs:
%   EEG            - EEGLAB data structure
%
% Author: Arnaud Delorme, CNL / Salk Institute, 12 Nov 2002
%
% See also: eeglab(), readegi(), readegihdr()

% Copyright (C) 12 Nov 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [EEG, command] = pop_readegi(filename, datachunks, forceversion, fileloc);

EEG = [];
command = '';
if nargin < 4, fileloc = 'auto'; end;
disp('Warning: This function can only import continuous files or');
disp('         epoch files with only one length for data epochs');

if nargin < 1
    % ask user
    [filename, filepath] = uigetfile('*.RAW;*.raw', ...
        'Choose an EGI RAW file -- pop_readegi()');
    drawnow;
    if filename == 0 return; end;
    filename = [filepath filename];

    fid = fopen(filename, 'rb', 'b');
    if fid == -1, error('Cannot open file'); end
    head = readegihdr(fid); % read EGI file header
    fclose(fid);
    if head.segments ~= 0
        fileloc = '';
        floc = { 'GSN-HydroCel-32.sfp' 'GSN65v2_0.sfp' 'GSN129.sfp' 'GSN-HydroCel-257.sfp'};
        switch head.nchan
            case { 32 33 }, fileloc = floc;
            case { 64 65 }, fileloc =  {floc{2} floc{3:4} floc{1}};
            case { 128 129 }, fileloc = {floc{3} floc{4} floc{1:2}};
            case { 256 257 }, fileloc = {floc{4} floc{1:3}};
        end;
        uilist = { { 'style' 'text' 'string' sprintf('Segment/frame number (default: 1:%d)', head.segments) } ...
                   { 'style' 'edit' 'string' '' } ...
                   { 'style' 'text' 'string' 'Channel location file (in eeglab/sample_locs)' } ...
                   { 'style' 'popupmenu' 'string' fileloc } ...
                   { } ...
                   { 'style' 'text' 'string' ...
                   [ 'Note: Choosing the correct electrode location file for your data is critical.' 10 ...
        'Note that in some cases none of the channel location files listed will correspond' 10 ...
        'to your montage as EGI has different versions of caps and is creating new ones' 10 ...
        'constantly. Remember to check your montage in the channel editor and import. Channel' 10 ...
        'location files are stored in the sample_locs sub-folder of the EEGLAB distribution.' ] } };
        uigeometry = { [2 1] [2 1] [1] [1] };
        uigeomvert = [1 1 1.5 3];
        result = inputgui('uilist', uilist, 'geometry', uigeometry, 'geomvert', uigeomvert);

%         promptstr    = { sprintf('Segment/frame number (default: 1:%d)', head.segments) 'Channel location file (in eeglab/sample_locs)' };
%         inistr       = { '' fileloc(res{2})};
%         result       = inputdlg2( promptstr, 'Import EGI file -- pop_readegi()', 1,  inistr, 'pop_readegi');
        if length(result) == 0 return; end;
        datachunks   = eval( [ '['  result{1} ']' ] );
        fileloc      = char(fileloc(result{2}));
    else
        datachunks   = [];
        disp('Only one segment, cannot read portion of the file');
    end;
end;

% load data
% ----------
EEG = eeg_emptyset;
if exist('datachunks') && exist('forceversion') && ~isempty(forceversion)
    [Head EEG.data Eventdata SegCatIndex] = readegi( filename, datachunks,forceversion);
elseif exist('forceversion') && ~isempty(forceversion)
    [Head EEG.data Eventdata SegCatIndex] = readegi( filename,[],forceversion);
elseif exist('datachunks')
    [Head EEG.data Eventdata SegCatIndex] = readegi( filename, datachunks);
    forceversion = [];
else
    [Head EEG.data Eventdata SegCatIndex] = readegi( filename);
    forceversion = [];
end
if ~isempty(Eventdata) & size(Eventdata,2) == size(EEG.data,2)
    EEG.data(end+1:end+size(Eventdata,1),:) = Eventdata;
end;
EEG.comments        = [ 'Original file: ' filename ];
EEG.setname 		= 'EGI file';
EEG.nbchan          = size(EEG.data,1);
EEG.srate           = Head.samp_rate;
EEG.trials          = Head.segments;
EEG.pnts            = Head.segsamps;
EEG.xmin            = 0;

% importing the events
% --------------------
EEG = eeg_checkset(EEG);
if ~isempty(Eventdata)
    orinbchans = EEG.nbchan;
    for index = size(Eventdata,1):-1:1
        EEG = pop_chanevent( EEG, orinbchans-size(Eventdata,1)+index, 'edge', 'leading', ...
            'delevent', 'off', 'typename', Head.eventcode(index,:), ...
            'nbtype', 1, 'delchan', 'on');
        Head.eventcode(end,:) = [];
    end;

    % renaming event codes
    % --------------------
    try,
        tmpevent = EEG.event;
        alltypes = { tmpevent.type };
        if isstr(alltypes{1})
            indepoc = strmatch('epoc', lower(alltypes), 'exact');
            indtim  = strmatch('tim0', lower(alltypes), 'exact');

            % if epoc but no tim0 then epoc represent pauses in recording
            if isempty(indtim) & ~isempty(indepoc)
                for index = indepoc
                    EEG.event(index).type = 'boundary';
                end;
            end;
            % other wise if both non-empty data epochs
            if ~isempty(indtim) & ~isempty(indepoc)
                if rem(size(EEG.data,2) / (length(indepoc)+1),1) == 0
                    EEG.event(index) = []; % remove epoch events
                    EEG.trials       = length(indepoc)+1;
                else
                    disp('Warning: data epochs detected but wrong data size');
                end;
            end;
        end;
    catch, disp('Warning: event renaming failed'); end;
end;

% adding segment category indices
% -------------------------------
if ~isempty(SegCatIndex) && EEG.trials > 1
    try
        if ~isempty(EEG.event)
            for index = 1:length(EEG.event)
                EEG.event(index).category = Head.catname{SegCatIndex(EEG.event(index).epoch)};
            end;
        else % create time-locking events
            for trial = 1:EEG.trials
                EEG.event(trial).epoch    = trial;
                EEG.event(trial).type     = 'TLE';
                EEG.event(trial).latency  = -EEG.xmin*EEG.srate+1+(trial-1)*EEG.pnts;
                EEG.event(trial).category = Head.catname{SegCatIndex(trial)};
            end;
        end;
    catch,
        disp('Warning: error while importing trial categories');
        EEG.event = rmfield(EEG.event, 'category');
    end;
end;
EEG = eeg_checkset(EEG, 'makeur');
EEG = eeg_checkset(EEG, 'eventconsistency');

% importing channel locations
% ---------------------------
if nargin < 1
    warndlg2( [ 'EEGLAB will now import a default electrode location file' 10 ...
        'for your data. Note that this might not correspond' 10 ...
        'to your montage as EGI has different versions of caps.' 10 ...
        'Check your montage in the channel editor and import' 10 ...
        'the correct location file if necessary.' ]);
end;
if all(EEG.data(end,1:10) == 0)
    disp('Deleting empty data reference channel (reference channel location is retained)');
    EEG.data(end,:)   = [];
    EEG.nbchan        = size(EEG.data,1);
    EEG = eeg_checkset(EEG);
end;
if ~isempty(fileloc)
    if strcmpi(fileloc, 'auto')
        EEG = readegilocs(EEG);
    else
        EEG = readegilocs(EEG, fileloc);
    end;
end;

if nargin < 1
    command = sprintf('EEG = pop_readegi(''%s'', %s);', filename, vararg2str({datachunks forceversion fileloc }) );
end;

return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function f = fxy5 ( n, x, y )

%*****************************************************************************80
%
%% FXY5 is the fourth 2D example.
%
%  Discussion:
%
%    This is example 3.1 in the reference.
%
%    It is known as the discontinuous medium wave function.
%
%    Here, we are computing the second component of the solution, U(X,Y).
%
%    It should be plotted on (x,y) in [-1,0]x[0,0.1].
%
%    The second variable y actually represents time.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Rick Archibald, Anne Gelb, Jungho Yoon,
%    Determining the locations and discontinuities in the derivatives
%    of functions,
%    Applied Numerical Mathematics,
%    Volume 58, 2008, pages 577-592.
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real X(N), Y(N), the arguments.
%
%    Output, real F(N,1), the function values.
%

%
%  Destroy all row vectors!
%
  x = x ( : );
  y = y ( : );

  cl = 0.87879;
  cr = 1.0;
  omega = 12.0;
  rhol = 0.55556;
  rhor = 1.0;

  f = zeros ( n, 1 );

  i = find ( x(1:n) <= -0.5 );
  j = find ( -0.5 < x(1:n) );

  f(i) = sin ( pi * omega * ( y(i) - ( x(i) + 0.5 ) / cl ) ) ...
    + ( rhol * cl - rhor * cr ) / ( rhol * cl + rhor * cr ) ...
    / ( rhol * cl ) ...
    * sin ( pi * omega * ( y(i) + ( x(i) + 0.5 ) / cl ) );

  f(j) = 2.0 / ( rhol * cl + rhor * cr ) ...
    * sin ( pi * omega * ( y(j) - ( x(j) + 0.5 ) / cl ) );

  return
end

repo: jburkardt-m