text:
function fig=vgg_gui_H(i1,i2,H)
%
%	fig=vgg_gui_H(i1,i2,H)
%
%
% Visualizes an homography matrix of two views
%
%IN:
%	i1 - Matlab image
%	i2 - Matlab image
%	H - 3x3 Homography matrix. Assumes that image coordiantes are 1..width
%		where pixel centers are at integer locations.
%
%OUT:
%	fig - handle to the figure

% $Id: vgg_gui_H.m,v 1.3 2001/10/30 14:42:41 wexler Exp $
% Yoni, Tue Mar 27 19:31:16 2001

if nargin==3
   action='start';
else
   action=i1;
   ud = get(gcf, 'UserData');
end

if strcmp(action,'start'),
   if nargin ~= 3
      error('Must give 3 arguments... read the docs.\n');
   end

   h0 = figure('Color',[0.8 0.8 0.8], ...
	       'NumberTitle','off', ...
	       'Name','Play With Homography Matrix', ...
	       'ButtonDownFcn', 'disp(''Click on images'')',...
	       'WindowButtonUpFcn', 'vgg_gui_H(''none'');',...
	       'WindowButtonMotionFcn', 'vgg_gui_H(''move'')', ...
	       'Pointer', 'crosshair', ...
	       'DoubleBuffer', 'on',...
	       'Units','normalized');
   m0=uimenu('Label', '&Color');
   uimenu(m0, 'Label', '&Red', 'ForegroundColor', [1 0 0], ...
	  'Accelerator', 'r', 'Callback', 'vgg_gui_H(''cr'');');
   uimenu(m0, 'Label', '&Green', 'ForegroundColor', [0 1 0], ...
	  'Accelerator', 'g', 'Callback', 'vgg_gui_H(''cg'');');
   uimenu(m0, 'Label', '&Blue', 'ForegroundColor', [0 0 1], ...
	  'Accelerator', 'b', 'Callback', 'vgg_gui_H(''cb'');');
   m1=uimenu('Label', '&Size');
   uimenu(m1, 'Label', '&Increase', 'Callback', 'vgg_gui_H(''s+'');', 'Accelerator', '+');
   uimenu(m1, 'Label', '&Decrease', 'Callback', 'vgg_gui_H(''s-'');', 'Accelerator', '-');
   uimenu(m1, 'Label', '&1', 'Callback', 'vgg_gui_H(''s1'');', 'Accelerator', '1');
   uimenu(m1, 'Label', '&2', 'Callback', 'vgg_gui_H(''s2'');', 'Accelerator', '2');
   uimenu(m1, 'Label', '&3', 'Callback', 'vgg_gui_H(''s3'');', 'Accelerator', '3');
   uimenu(m1, 'Label', '&4', 'Callback', 'vgg_gui_H(''s4'');', 'Accelerator', '4');
   uimenu(m1, 'Label', '&5', 'Callback', 'vgg_gui_H(''s5'');', 'Accelerator', '5');
   uimenu(m1, 'Label', '&6', 'Callback', 'vgg_gui_H(''s6'');', 'Accelerator', '6');
   uimenu(m1, 'Label', '&7', 'Callback', 'vgg_gui_H(''s7'');', 'Accelerator', '7');
   uimenu(m1, 'Label', '&8', 'Callback', 'vgg_gui_H(''s8'');', 'Accelerator', '8');
   uimenu(m1, 'Label', '&9', 'Callback', 'vgg_gui_H(''s9'');', 'Accelerator', '9');

   ah1 = axes('Parent', h0, ...
	      'Position',[0 0 .5 1]);
   h1=imshow(i1); hold on; title('Image 1');
   set(h1, 'ButtonDownFcn','vgg_gui_H(''b1'');');

   ah2 = axes('Parent',h0, ...
	      'Position',[.5 0 .5 1], ...
	      'Tag','Axes2');
   h2=imshow(i2); hold on; title('Image 2');
   set(h2, 'ButtonDownFcn','vgg_gui_H(''b2'');');

   point=plot(-1000, -1000,'EraseMode','xor');
   point2=plot(-1000, -1000,'EraseMode','xor');

   s1=size(i1); s2=size(i2);
   t(:,:,1)=H;  t(:,:,2)=inv(H);  H=t;

   ud=struct('h0', h0, 'h',[h1 h2], 'ah', [ah1, ah2], ...
	     'sizes', [s1(1:2); s2(1:2)], ...
	     'current', -1, 'color', 'r', 'size', 1, ...
	     'p', point, 'H', H, 'p2', point2 );

   set(h0,'UserData',ud);

   if nargout > 0, fig = h0; end

elseif strcmp(action, 'move')
   if ud.current<0 return; end;
   pt=get(ud.ah(ud.current),'CurrentPoint');
   pt2=ud.H(:,:,ud.current)*pt(1,:)';
   pt2=pt2/pt2(3);

   set(ud.p2, 'XData', pt2(1,1), 'YData', pt2(2,1))
   set(ud.p, 'XData', pt(1,1), 'YData', pt(1,2))

elseif action(1)=='b'
   if action(2)=='1' ud.current=1;
   elseif action(2)=='2' ud.current=2;
   else return;
   end
   pt=get(ud.ah(ud.current),'CurrentPoint');

   p2=ud.H(:,:,ud.current)*pt(1,:)';
   p2=p2/p2(3);

   delete(ud.p);
   delete(ud.p2);
   axes(ud.ah(ud.current));
   ud.p=plot(pt(1,1), pt(1,2), [ud.color '+'], ...
	     'MarkerSize', 8+2*ud.size, 'LineWidth', ud.size,...
	     'EraseMode','xor');
   axes(ud.ah(3-ud.current));
   ud.p2=plot(p2(1,1), p2(2,1), [ud.color '+'], ...
	     'MarkerSize', 8+2*ud.size, 'LineWidth', ud.size,...
	     'EraseMode','xor');

elseif action(1)=='c'
   ud.color=action(2);
   %get(ud.l)
   set(ud.p2, 'Color', ud.color);
   set(ud.p, 'Color', ud.color);

elseif action(1)=='s'
   if action(2)=='+' ud.size=ud.size+1;
   elseif action(2)=='-' ud.size=max(1, ud.size-1);
   else
      ud.size = str2num(action(2));
   end
   set(ud.p, 'LineWidth', ud.size, 'MarkerSize', 8+2*ud.size);
   set(ud.p2, 'LineWidth', ud.size, 'MarkerSize', 8+2*ud.size);

elseif strcmp(action, 'none')
   ud.current = -1;

else
   error(['Unknown command: ' action]);
end

set(ud.h0, 'UserData',ud);

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function pass = test_constructor_inputs(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end

% [TODO]: This test needs to be updated to include more exotic input options.

% Generate a few random points to use as test values.
seedRNG(6178);
xx = 2 * rand(100, 1) - 1;

% Test the vectorise flag:
try
    f = chebfun(@(x) x^2, pref, 'vectorize'); %#ok<NASGU>
    pass(1) = true;
catch
    pass(1) = false;
end

% Test the vectorise flag:
f = chebfun(@(x) sin(x), pref, 5);
pass(2) = ( length(f.funs{1}) == 5 );
% [TODO]: Change this once @CHEBFUN/LENGTH is implemented.

% Test the 'splitting on' flag.
f = chebfun(@(x) abs(x), 'splitting', 'on');
pass(3) = all(size(f.funs) == [1, 2]);

% Test the 'extrapolate' flag.
f = chebfun(@(x) sign(x), [0, 1], 'extrapolate', 'on');
pass(4) = get(f, 'ishappy');

% Test construction from an array of FUN objects:
f = chebfun({0, 1, 2, 3, 4, 5}, [0, 1, 2, 3, 4, 5, 6]);
g = chebfun(f.funs);
pass(5) = all(size(g.funs) == [1, 6]) && all(g.domain == 0:6);

% Test array-valued construction from an array of FUN objects:
f = chebfun({[0, 1], [2, 3], [4, 5]}, [0, 1, 2, 3]);
g = chebfun(f.funs);
pass(6) = all(size(g.funs) == [1, 3]) && all(g.domain == 0:3);

% Test fixed-length construction:
f = chebfun(@sin, 10);
pass(7) = numel(f.funs) == 1 && size(f.funs{1},1) == 10;

% Test equispaced construction:
f = chebfun(rand(10,3), 'equi');
pass(8) = numel(f.funs) == 1 && size(f.funs{1}, 1) > 10;

% Test construction from coefficients.
f = chebfun([1 ; 2 ; 3], 'coeffs');
pass(9) = isequal(f.funs{1}.onefun.coeffs, [1 ; 2 ; 3]);

% Test 'chebkind' and 'kind' flags.
f1 = chebfun(@(x) x, 'chebkind', '1st');
f3 = chebfun(@(x) x, 'chebkind', 1);
pass(10) = isa(f1.funs{1}.onefun, 'chebtech1') && ...
    isa(f3.funs{1}.onefun, 'chebtech1');

% Test construction from numeric string.
f = chebfun('1');
pass(11) = all(feval(f, linspace(-1, 1, 10)) == 1);

% Test 'trunc', flag.
f = chebfun(@abs, 'trunc', 10, 'splitting', 'on');
c = get(f, 'coeffs');
pass(12) = abs(-4/63/pi - c(9)) < 10*eps;

% Test construction from cells of strings:
f = chebfun({'x','x-1'}, [0 1 2]);
pass(13) = norm(feval(f, [.5, 1.5]) - .5) < eps;

% Test construction from a chebfun:
x = chebfun('x', [0, 5]);
f = chebfun(x);
pass(13) = all(domain(f) == [0 5] );

f = chebfun(x, [0 2]);
pass(14) = all(domain(f) == [0 2] );

% Test construction from a piecewise chebfun:
f = chebfun(chebfun({'x','x'}, [-1 0 1]));
x = [-.5, .5];
pass(15) = numel(f.funs) == 1  && norm(feval(f, x) - x) < eps;

% Test 'minSamples' flag.
f_op = @(x) -x - x.^2 + exp(-(50*(x - .5)).^4);
f1 = chebfun(f_op, 'minSamples', 17);
err1 = norm(feval(f1, xx) - f_op(xx), inf);
f2 = chebfun(f_op, 'minSamples', 33);
err2 = norm(feval(f2, xx) - f_op(xx), inf);
pass(16) = (err1 > 1e-3) && (err2 < 1e2*vscale(f2)*eps);

% Test support for "legacy" preferences.
f_op = @(x) sin(200*x);

f = chebfun(f_op, 'resampling', 'on');
pass(17) = ishappy(f);

warnstate = warning('off','CHEBFUN:CHEBFUN:constructor:notResolved');
f = chebfun(f_op, 'maxdegree', 129, 'tech', @chebtech2);
pass(18) = ~ishappy(f) && (length(f) == 129);
warning(warnstate);

f = chebfun(f_op, 'splitting', 'on', 'splitdegree', 65);
pass(19) = ishappy(f) && all(cellfun(@(fk) length(fk) <= 65, f.funs));

warnstate = warning('off','CHEBFUN:CHEBFUN:constructor:funNotResolved');
f = chebfun(f_op, 'splitting', 'on', 'splitLength', 65, 'splitMaxLength', 200);
pass(20) = ~ishappy(f) && (length(f) <= 65*4);
warning(warnstate);

% Test construction with a mixture of preference object and keyword inputs.
p = pref;
p.tech = @chebtech1;
f = chebfun(@(x) 1./x, [0 1], 'exps', [-1 0], p);
pass(21) = get(f, 'ishappy') && isa(f.funs{1}.onefun.smoothPart, 'chebtech1');
f = chebfun(@(x) 1./x, [0 1], p, 'exps', [-1 0]);
pass(22) = get(f, 'ishappy') && isa(f.funs{1}.onefun.smoothPart, 'chebtech1');

p = struct();
p.tech = @chebtech1;
f = chebfun(@(x) 1./x, [0 1], 'exps', [-1 0], p);
pass(23) = get(f, 'ishappy') && isa(f.funs{1}.onefun.smoothPart, 'chebtech1');
f = chebfun(@(x) 1./x, [0 1], p, 'exps', [-1 0]);
pass(24) = get(f, 'ishappy') && isa(f.funs{1}.onefun.smoothPart, 'chebtech1');

% Test the vectorise flag on an array-valued function:
try
    f = chebfun(@(x) [x^2 sin(x)*cos(x)], pref, 'vectorize'); %#ok<NASGU>
    pass(25) = true;
catch
    pass(25) = false;
end

% Test the 'eps' preference.
f = chebfun(@sin);
g = chebfun(@sin, 'eps', 1e-6);
pass(26) = length(g) < length(f);

% Test the use of 'coeffs' and 'chebkind' simultaneously.
c = (1:5).';
try
    chebfun(c, 'coeffs', 'chebkind', 1);
catch ME
    pass(27) = strcmpi(ME.message, [' ''coeffs'' and ''chebkind'' ' ...
        'should not be specified simultaneously.']);
end

try
    chebfun(c, 'chebkind', 2, 'coeffs');
catch ME
    pass(28) = strcmpi(ME.message, [' ''coeffs'' and ''chebkind'' ' ...
        'should not be specified simultaneously.']);
end

% Test construction from a chebfun row (transposed chebfun). The result
% should be a chebfun column (not transposed).
f = chebfun(@(x) cos(pi*x))';
g = chebfun(f);
pass(29) = g.isTransposed == 0;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function SR = superresolve_wnuisr(slidingWindows, magFactor)

numberOfFrames = size(slidingWindows.frames,3);

% Set PSF kernels for the different binning factors.
psf{1} = fspecial('gaussian', 5, 0.8); % PSF for binning factor 2
psf{2} = fspecial('gaussian', 7, 1.2); % PSF for binning factor 3
psf{3} = fspecial('gaussian', 9, 1.6); % PSF for binning factor 4

% Specific preprocessing steps for the non-uniform interpolation based multi-frame methods.
weighting_struct = struct('mec_weight_flag',4,'dis_weight_flag',1,'qps_weight_flag',0,'mec_confMapN',[],...
    'mec_thr_scaler',2,'mec_weight_exp',2,'mec_weighting',[],'dis_rho',0.7,'dis_scaler',magFactor*10,'dis_weighting',[],...
    'qps_map',[],'qps_rho',0.7,'qps_weighting',[],'alpha',1,'beta',1,'gamma',1,'whi_mode',2);
warped_meshXN = zeros(size(slidingWindows.referenceFrame,1),size(slidingWindows.referenceFrame,2),numberOfFrames-1);
warped_meshYN = zeros(size(slidingWindows.referenceFrame,1),size(slidingWindows.referenceFrame,2),numberOfFrames-1);

lr_seq_sorted = slidingWindows.frames;
lr_seq_tmp = lr_seq_sorted(:,:,(numberOfFrames+1)/2);
lr_seq_sorted(:,:,(numberOfFrames+1)/2) = lr_seq_sorted(:,:,1);
lr_seq_sorted(:,:,1) = lr_seq_tmp;

lr_seq_sorted = permute(lr_seq_sorted,[1 2 4 3]);

permutation_vec = 1:numberOfFrames-1;
permutation_vec((numberOfFrames+1)/2:end) = permutation_vec((numberOfFrames+1)/2:end) + 1;
permutation_vec(1:(numberOfFrames-1)/2)   = circshift(permutation_vec(1:(numberOfFrames-1)/2),[0 -1]);

for itera = 1:numberOfFrames-1
    [warped_meshXN(:,:,itera),warped_meshYN(:,:,itera)] = lmsSR_warpImgGridLocalTrans(size(slidingWindows.referenceFrame,1),size(slidingWindows.referenceFrame,2),slidingWindows.flowToReference{permutation_vec(itera)}.mvs_xc);

    % Storing the MEC-SSD-Confidence into the Weighting Struct
    weighting_struct.mec_confMapN(:,:,itera) = slidingWindows.mec_confMap{permutation_vec(itera)};
end

% Actual SR call
[SR,~,~,~] = lmsSR_generateSRusingNUIv4weighted(lr_seq_sorted,warped_meshXN,warped_meshYN,[magFactor, magFactor],psf{magFactor-1},2,weighting_struct);

author: thomas-koehler

repo: SupER


text:
function vw = wangAtlasToROIs(vw, wangAtlasNifti)
%
% Converts Wang Atlas from freesurfer MGZ volume to vistasoft ROIs
%
% Example
%   fsID = 'ernie';
%

hideGrayROIs = viewGet(vw, 'hide gray ROIs');
vw = viewSet(vw, 'hide Gray ROIs', true);

Wang_ROI_Names = {...
    'V1v' 'V1d' 'V2v' 'V2d' 'V3v' 'V3d' 'hV4' 'VO1' 'VO2' 'PHC1' 'PHC2' ...
    'TO2' 'TO1' 'LO2' 'LO1' 'V3B' 'V3A' 'IPS0' 'IPS1' 'IPS2' 'IPS3' 'IPS4' ...
    'IPS5' 'SPL1' 'FEF'};

wangComments = 'ROI made based on maximum probabilty map from Wang et al, 2014 ''Probabilistic Maps of Visual Topography in Human Cortex''';

% count the current number of ROIs
nROIs = length(viewGet(vw, 'ROIs'));

% create new ROIs from atlas
vw = nifti2ROI(vw, wangAtlasNifti);

% how many new ROIs did we get?
numNewROIs = length(viewGet(vw, 'ROIs')) - nROIs;
assert(numNewROIs>0)

for jj = 1:numNewROIs
    this_roi       = viewGet(vw, 'ROI name', nROIs + jj);
    this_roi_num   = str2double(this_roi(regexp(this_roi, '[0-9]')));
    this_roi_name  = sprintf('WangAtlas_%s', Wang_ROI_Names{this_roi_num});
    vw             = viewSet(vw, 'ROI Name', this_roi_name, nROIs + jj);
    comments       = viewGet(vw, 'ROI comments', nROIs + jj);
    comments       = sprintf('%s\n%s', comments, wangComments);
    vw             = viewSet(vw, 'ROI comments', comments, nROIs + jj);
end

vw = viewSet(vw, 'hide Gray ROIs', hideGrayROIs);

return
%
%
% colors = jet(25);
%
% colors = colors(randperm(25),:)
%
% for ii = 1:length(vw.ROIs); vw.ROIs(ii).color = colors(ii,:); end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function verifyVarsEqual(testCase, name, varargin)
%verifyVarsEqual  Verify the equality of demo results within tolerance.
%
%  Description
%    verifyVarsEqual(testCase, name, vars, funcs, OPTIONS)
%    compares the elements in the results of a demo into the expected ones
%    using relative tolerance to each element individually and/or to the
%    range of the expected array and/or using absolute tolerance. The data
%    is loaded from the 'realValues' and 'testValues' folders. Used by the
%    gpstuff tests. Can be used with xUnit Test Framework package by Steve
%    Eddins or with the Matlab's built-in Unit Testing Framework (as of
%    version 2013b).
%
%    Parameters:
%      testCase
%        The testCase object, if using the integrated unit test framework,
%        or empty array, if using the xUnit package.
%      name
%        The name of the demo.
%      vars (optional)
%        Cell array of strings defining the numerical array variables that
%        are compared. Giving string 'same' (default) looks the names of
%        the saved variables in the file 'realValues/<nameOfTheDemo>.mat'
%        and compares them.
%      funcs (optional)
%        Cell array of function handles that are applied into the
%        corresponding actual and expected arrays before comparison. Cell
%        containing an empty array indicates that no function is used for
%        that variable. Empty array in general indicates that no function
%        is applied for any variable. Providing a single function handle
%        applies the function to all variables.
%
%    OPTIONS is optional parameter-value pair
%      RelTolElement - tolerance relative to the magnitude of each element
%                      (default 0.05)
%      RelTolRange   - tolerance relative to the range of the elements
%                      (default 0.01)
%      AbsTol        - absolute tolerance (default 0)
%
%   See also
%     TEST_*
%
% Copyright (c) 2014 Tuomas Sivula

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

% Parse inputs
ip = inputParser;
ip.FunctionName = 'verifyVarsEqual';
ip.addRequired('testCase',@(x) isempty(x) ...
    || isa(testCase, 'matlab.unittest.FunctionTestCase'))
ip.addRequired('name', @ischar)
ip.addOptional('vars', 'same', @(x) iscellstr(x) ...
    || (ischar(x) && strcmp(x,'same')))
ip.addOptional('funcs', [], @(x) isempty(x) || ...
    isa(x,'function_handle') || (iscell(x) && all( ...
    cellfun(@(y)isa(y,'function_handle')||isempty(y),x) )))
ip.addParamValue('RelTolElement', 0, ...
  @(x) isnumeric(x) && isscalar(x) && isreal(x) && x >= 0)
ip.addParamValue('RelTolRange', 0, ...
  @(x) isnumeric(x) && isscalar(x) && isreal(x) && x >= 0)
ip.addParamValue('AbsTol', 0, ...
  @(x) isnumeric(x) && isscalar(x) && isreal(x) && x >= 0)
ip.parse(testCase,name,varargin{:});
testCase = ip.Results.testCase;
name = ip.Results.name;
vars = ip.Results.vars;
funcs = ip.Results.funcs;
tol_el = ip.Results.RelTolElement;
tol_range = ip.Results.RelTolRange;
tol_abs = ip.Results.AbsTol;

% Handle parameter vars
if ischar(vars) && strcmp(vars, 'same')
  % Check the variable names from 'realValues/<nameOfTheDemo>.mat'
  vars = whos(matfile(['realValues/' name '.mat']));
  vars = {vars.name};
end

% Validate parameter funcs
if ~isempty(funcs) && iscell(funcs) && length(funcs) ~= length(vars)
  error('Parameter funcs size mismatch')
end

% Find the path to the test folder
path = mfilename('fullpath');
path = path(1:end-length(mfilename));

% For every variable in vars
for i = 1:length(vars)

  % Load the variables
  if isempty(testCase)
    % xUnit package
    Sw = warning('off','MATLAB:load:variableNotFound')
    try
      actual = load([path, 'testValues/' name '.mat'], vars{i});
      if ~isfield(actual, vars{i})
        error(['Variable ' vars{i} ' not found in testValues/' name '.mat'])
      end
      expected = load([path, 'realValues/' name '.mat'], vars{i});
      if ~isfield(expected, vars{i})
        error(['Variable ' vars{i} ' not found in realValues/' name '.mat'])
      end
    catch err
      warning(Sw)
      rethrow(err)
    end
    warning(Sw)
  else
    % Built-in test framework
    actual = assertWarningFree(testCase, ...
      @() load([path, 'testValues/' name '.mat'], vars{i}));
    expected = assertWarningFree(testCase, ...
      @() load([path, 'realValues/' name '.mat'], vars{i}));
  end

  % Error message (if the variable should fail the validation)
  msg = ['The variable ' vars{i} ' failed the verification.'];

  % Apply funcs
  if isa(funcs,'function_handle')
    actual = funcs(actual.(vars{i}));
    expected = funcs(expected.(vars{i}));
    %Add info to the error message
    msg = sprintf('%s\nFunction applied:\n%s', msg, func2str(funcs));
  elseif ~isempty(funcs) && ~isempty(funcs{i})
    actual = funcs{i}(actual.(vars{i}));
    expected = funcs{i}(expected.(vars{i}));
    %Add info to the error message
    msg = sprintf('%s\nFunction applied:\n%s', msg, func2str(funcs{i}));
  else
    actual = actual.(vars{i});
    expected = expected.(vars{i});
  end

  if isempty(testCase)

    % xUnit package (N.B. using assertions instead of verifications)
    if tol_el == 0 && tol_range == 0 && tol_abs == 0
      % No tolerance
      assertEqual(actual, expected, msg)
    else
      % Manual assertion

      % Compare
      diff = abs(actual(:) - expected(:));
      viol_el = 1;
      viol_range = 1;
      viol_abs = 1;
      if tol_el > 0
        viol_el = diff > tol_el*abs(expected(:));
      end
      if tol_range > 0
        viol_range = diff > tol_range ...
          *( max(expected(:)) - min(expected(:)) );
      end
      if tol_abs > 0
        viol_abs = diff > tol_abs;
      end

      if any( viol_el & viol_range & viol_abs )

        % Compose the message
        msg = sprintf( ...
          ['%s\nInput elements are not all equal within tolerances\n' ...
           '(RelTolElement=%g, RelTolRange=%g, AbsTol=%g)\n\n'], ...
          msg, tol_el, tol_range, tol_abs);
        if tol_el > 0 && any(viol_el)
          msg = sprintf( ...
            ['%sThe following indexes fail the elementwise relative ' ...
             'tolerance:\n%s\n\n'], ...
            msg, xunit.utils.arrayToString(find(viol_el)'));
        end
        if tol_range > 0 && any(viol_range)
          msg = sprintf( ...
            ['%sThe following indexes fail the range relative ' ...
             'tolerance:\n%s\n\n'], ...
            msg, xunit.utils.arrayToString(find(viol_range)'));
        end
        if tol_abs > 0 && any(viol_abs)
          msg = sprintf( ...
            ['%sThe following indexes fail the absolute ' ...
             'tolerance:\n%s\n\n'], ...
            msg, xunit.utils.arrayToString(find(viol_abs)'));
        end
        msg = sprintf( ...
          ['%sThe following indexes fail the comparison(s):\n%s\n\n' ...
           'First input:\n%s\n\nSecond input:\n%s'], ...
          msg, ...
          xunit.utils.arrayToString(find(viol_el & viol_range & viol_abs)'), ...
          xunit.utils.arrayToString(actual), ...
          xunit.utils.arrayToString(expected) );

        % Raise error
        throwAsCaller(MException('assertVars:tolExceeded', msg));

      end
    end

  else

    if tol_el == 0 && tol_range == 0 && tol_abs == 0
      % No tolerance
      verifyEqual(testCase, actual, expected, msg)
    else
      % Construct tolerance object
      if tol_el > 0
        tolObj = matlab.unittest.constraints.RelativeTolerance(tol_el);
      end
      if tol_range > 0
        if exist('tolObj', 'var')
          tolObj = tolObj | matlab.unittest.constraints.AbsoluteTolerance( ...
            tol_range*( max(expected(:)) - min(expected(:)) ));
        else
          tolObj = matlab.unittest.constraints.AbsoluteTolerance( ...
            tol_range*( max(expected(:)) - min(expected(:)) ));
        end
      end
      if tol_abs > 0
        if exist('tolObj', 'var')
          tolObj = tolObj | matlab.unittest.constraints.AbsoluteTolerance( ...
            tol_abs);
        else
          tolObj = matlab.unittest.constraints.AbsoluteTolerance(tol_abs);
        end
      end

      % Perform the verification
      testCase.verifyThat(actual, matlab.unittest.constraints.IsEqualTo( ...
        expected, 'Within', tolObj), msg);

    end

  end

  % Clear variables
  clear actual expected tolObj

end

author: gpstuff-dev

repo: gpstuff

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


text:
function by = cvx_getlog( bx )

% WARNING: This assumes that there is exactly one non-zero element per
% column, that all non-zeros are positive, and that any element appearing
% in rows greater than 1 has a non-zero entry in cvx___.logarithm.

global cvx___
nb = size( bx, 2 );
[ rx, cx, vx ] = find( bx );
tt = rx > 1;
ry = rx(tt);
cy = cx(tt);
logs = cvx___.logarithm( ry, 1 );
nl = max([1,max(logs)]);
by = sparse( logs, cy, 1, nl, nb ) + sparse( 1, cx, log(vx), nl, nb );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.tx for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function dfaddbuttons(dffig)
%DFADDBUTTONS Add buttons to the curve fitting plot figure window

%   $Revision: 1.1.6.8 $  $Date: 2004/01/24 09:35:09 $
%   Copyright 2003-2004 The MathWorks, Inc.

% Clear out any old stuff
h0 = findall(dffig,'Type','uicontrol','Style','pushbutton');
if ~isempty(h0), delete(h0); end
p0 = ones(1,4);   % temporary position before adjustment

h0=uicontrol(dffig,'units','pixels','Tag','selectionframe',...
             'Style','frame','Position',p0);
h1=uicontrol(dffig,'units','pixels','Tag','displaytext',...
            'String','Display type:', 'style','text','Position',p0,...
            'HorizontalAlignment','left','FontWeight','bold');

choices = ['Density (PDF)|Cumulative probability (CDF)|Quantile (inverse CDF)|Probability plot|'...
           'Survivor function|Cumulative hazard'];
h2=uicontrol(dffig,'units','pixels','Tag','displaylist',...
            'String',choices, 'Style','pop','BackgroundColor',ones(1,3), ...
            'Callback', @cbkfunction,'Position',p0);
setappdata(h2,'codenames',...
             {'pdf' 'cdf' 'icdf' 'probplot' 'survivor' 'cumhazard'});

h3=uicontrol(dffig,'units','pixels','Tag','typetext',...
            'String','Distribution:', 'style','text','Position',p0,...
            'HorizontalAlignment','left','FontWeight','bold', 'Enable', 'off');

% Figure out which distributions could be used for probability plots,
% and remember their names and properties for later
alldist = dfgetdistributions;
dnames = {alldist.name};
dcodes = {alldist.code};
islocscale = [alldist.islocscale];
choices = dnames(islocscale);
default = strmatch('Normal',choices);
if length(default)~=1
   default = 1;
end
h4=uicontrol(dffig,'units','pixels','Tag','typelist',...
            'String',choices, 'Style','pop','BackgroundColor',ones(1,3),...
            'Value',default,'Enable','off', 'Callback',@cbkfunction,...
             'Position',p0);

% Store information about all distributions that could be used to make
% a probability plot, then separately list distributions that are ok
% given the current data.  Right now they're the same, but if we read
% in negative data we'll prune the "ok" list to omit positive
% distributions.
setappdata(h4,'allfullnames',choices);
setappdata(h4,'allcodenames',dcodes(islocscale));
setappdata(h4,'alldistspec',alldist(islocscale));
setappdata(h4,'okcodenames',dcodes(islocscale));

hvec = [h0 h1 h2 h3 h4];
setappdata(dffig,'selectioncontrols',hvec);

% Define information for the buttons
strings = {xlate('Data...') xlate('New Fit...') xlate('Manage Fits...') xlate('Evaluate...') xlate('Exclude...')};
tips = {xlate('Import, view, rename, plot and delete data') ...
        xlate('Add a fitted distribution') ...
        xlate('Edit, view, plot and rename fits') ...
        xlate('Evaluate fits to compute a table of results') ...
        xlate('Define rules for excluding data from a fit')};
cbacks = {@cbkdata @cbknewfit @cbkmanfit @cbkeval @cbkexclude};

% Add the buttons to the figure
tags = {'dfdata' 'dfnewfit' 'dfmanfit' 'dfevaluate' 'dfexclude'};
n = length(strings);
h = zeros(1,n);
for j=1:length(strings)
   h(j) = uicontrol(dffig,'Units','pixels', ...
                    'Position',[j,1,1,1],...
                    'String',strings{j}, 'TooltipString',tips{j}, ...
                    'Callback',cbacks{j}, 'Tag',tags{j});
end
setappdata(dffig,'buttoncontrols',h);

% ---------------------- callback for Import button
function cbkdata(varargin)
%CBKDATA Callback for Data button

delete(findall(gcbf,'Tag','dfstarthint'));
com.mathworks.toolbox.stats.Data.showData;

% ---------------------- callback for New Fit button
function cbknewfit(varargin)
%CBKNEWFIT Callback for New Fit button

com.mathworks.toolbox.stats.Fitting.getFitting.showNewFit;

% ---------------------- callback for Manage Fit button
function cbkmanfit(varargin)
%CBKMANFIT Callback for Manage Fit button

com.mathworks.toolbox.stats.FitManager.showFitManager;

% ---------------------- callback for Evaluate button
function cbkeval(varargin)
%CBKEVAL Callback for Evaluate button

% Get the current horizontal axis limits of the main figure.  We'll set the
% default points at which to evaluate to a "pretty" colon expression that
% spans those limits, with 10 steps.
xlims = dfgetset('xminmax');

if isempty(xlims)
   xstr = '';
else
   % Start out by choosing a rounding that will give the smaller (in
   % magnitude) of min(x) and max(x) a single sig digit, but at most two sig
   % digits in the larger.  The smaller may round to zero.  If the min and max
   % round to the same thing, use more digits until we get rounded numbers
   % that differ.
   xmag = max(min(abs(xlims)), max(abs(xlims))/10);
   rounder = 10^floor(log10(xmag));
   while true
       xmin = round(xlims(1)./rounder) * rounder;
       xmax = round(xlims(2)./rounder) * rounder;
       if xmin < xmax, break; end
       rounder = rounder/10;
   end

   % Create 10 steps, where the step size will have three more significant
   % digits than the endpoints.
   stepRounder = rounder ./ 1000;
   step = floor((xmax-xmin)./(10*stepRounder)) * stepRounder;

   % Figure out how many digits we need display in order to distinguish the
   % endpoints.  That's the number of digits to the left of the decimal, plus
   % however many we've rounded to on the right.  Allow for at least four so
   % that we'll get, e.g., "4000", and not "4e+03".
   xminDigits = max(round(log10(max(abs(xmin),1))-log10(rounder)),4);
   xmaxDigits = max(round(log10(max(abs(xmax),1))-log10(rounder)),4);
   xstr = sprintf('%0.*g:%g:%0.*g', xminDigits, xmin, step, xmaxDigits, xmax);
end

% Get the current plot type of the main figure, we'll set the default
% function type to evaluate based on that.
ftype = dfgetset('ftype');
if strcmp(ftype,'probplot')
    ftype = 'cdf'; % can't evaluate a prob plot, use cdf instead
% else {'pdf' 'cdf' 'survivor' 'icdf' 'cumhazard' 'hazrate'}
    % otherwise use the current setting
end

com.mathworks.toolbox.stats.Evaluate.showEvaluate(ftype,xstr);

% ---------------------- callback for Exclude button
function cbkexclude(varargin)
%CBKEXCLUDE Callback for Exclude button

com.mathworks.toolbox.stats.Exclude.showExclude;

% ---------------------- callback for display list
function cbkfunction(varargin)
%CBKFUNCTION Callback for setting function to display

% Get the requested function and distribution types
dffig = gcbf;
dfsetplottype(dffig);

author: zouchuhang

repo: LayoutNet


text:
% pop_newtimef() - Returns estimates and plots of event-related (log) spectral
%           perturbation (ERSP) and inter-trial coherence (ITC) changes
%           timelocked to a set of input events in one data channel.
%
% Usage:
%   >> pop_newtimef(EEG, typeplot); % pop_up window
%   >> pop_newtimef(EEG, typeplot, lastcom); % pop_up window
%   >> pop_newtimef(EEG, typeplot, channel); % do not pop-up
%   >> pop_newtimef(EEG, typeproc, num, tlimits,cycles,
%                        'key1',value1,'key2',value2, ... );
%
% Inputs:
%   INEEG    - input EEG dataset
%   typeproc - type of processing. 1 process the raw
%              data and 0 the ICA components
%   num      - component or channel number
%   tlimits  - [mintime maxtime] (ms) sub-epoch time limits
%   cycles   -  >0 -> Number of cycles in each analysis wavelet
%               0 -> Use FFTs (with constant window length)
%
% Optional inputs:
%    See the newtimef() function.
%
% Outputs: same as newtimef(), no outputs are returned when a
%          window pops-up to ask for additional arguments
%
% Getting the ERSP and ITC output values:
% Simply look up the history using the eegh function (type eegh).
% Then copy and paste the pop_newtimef command call and manually add output
% (see the newtimef function for a list of outputs). For instance
% [ersp itc powbase times frequencies] = pop_newtimef( EEG, ....);
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: newtimef(), eeglab()

% Copyright (C) 2002 arno@salk.edu, Arnaud Delorme, CNL / Salk Institute
%
% 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

% 01-25-02 reformated help & license -ad
% 03-08-02 add eeglab option & optimize variable sizes -ad
% 03-10-02 change newtimef call -ad
% 03-18-02 added title -ad & sm
% 04-04-02 added outputs -ad & sm

function varargout = pop_newtimef(EEG, typeproc, num, tlimits, cycles, varargin );

varargout{1} = '';
% display help if not enough arguments
% ------------------------------------
if nargin < 2
	help pop_newtimef;
	return;
end;
lastcom = [];
if nargin < 3
	popup = 1;
else
	popup = isstr(num) | isempty(num);
	if isstr(num)
		lastcom = num;
	end;
end;

% pop up window
% -------------
if popup
	[txt vars] = gethelpvar('newtimef.m');

    g = [1 0.3 0.6 0.34];
	geometry = { g g g g g g g g [1.025 1.27] [1] [1.2 1 1.2]};
    uilist = { ...
               { 'Style', 'text', 'string', fastif(typeproc, 'Channel number', 'Component number'), 'fontweight', 'bold'  } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,3,[],'1') 'tag' 'chan'} {} {} ...
               ...
			   { 'Style', 'text', 'string', 'Sub epoch time limits [min max] (msec)', 'fontweight', 'bold' } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,4,[],[ int2str(EEG.xmin*1000) ' ' int2str(EEG.xmax*1000) ]) 'tag' 'tlimits' } ...
               { 'Style', 'popupmenu', 'string', 'Use 50 time points|Use 100 time points|Use 150 time points|Use 200 time points|Use 300 time points|Use 400 time points' 'tag' 'ntimesout' 'value' 4} { } ...
			   ...
			   { 'Style', 'text', 'string', 'Frequency limits [min max] (Hz) or sequence', 'fontweight', 'bold' } ...
			   { 'Style', 'edit', 'string', '' 'tag' 'freqs'  } ...
               { 'Style', 'popupmenu', 'string', 'Use limits, padding 1|Use limits, padding 2|Use limits, padding 4|Use actual freqs.' 'tag' 'nfreqs' }  ...
               { 'Style', 'checkbox', 'string' 'Log spaced' 'value' 0 'tag' 'freqscale' } ...
			   ...
			   { 'Style', 'text', 'string', 'Baseline limits [min max] (msec) (0->pre-stim.)', 'fontweight', 'bold' } ...
			   { 'Style', 'edit', 'string', '0' 'tag' 'baseline' } ...
			   { 'Style', 'popupmenu',  'string', 'Use divisive baseline|Use standard deviation' 'tag' 'basenorm' } ...
               { 'Style', 'checkbox', 'string' 'No baseline' 'tag' 'nobase' } ...
               ...
               { 'Style', 'text', 'string', 'Wavelet cycles [min max/fact] or sequence', 'fontweight', 'bold' } ...
               { 'Style', 'edit', 'string', getkeyval(lastcom,5,[],'3 0.5') 'tag' 'cycle' } ...
               { 'Style', 'popupmenu', 'string', 'Use limits|Use actual seq.' 'tag' 'ncycles' } ...
               { 'Style', 'checkbox', 'string' 'Use FFT' 'value' 0 'tag' 'fft' } ...
			   ...
			   { 'Style', 'text', 'string', 'ERSP color limits [max] (min=-max)', 'fontweight', 'bold' } ...
               { 'Style', 'edit', 'string', '' 'tag' 'erspmax'} ...
               { 'Style', 'checkbox', 'string' 'see log power (set)' 'tag' 'scale' 'value' 1} {} ...
			   ...
			   { 'Style', 'text', 'string', 'ITC color limits [max]', 'fontweight', 'bold' } ...
               { 'Style', 'edit', 'string', '' 'tag' 'itcmax'} ...
               { 'Style', 'checkbox', 'string' 'plot ITC phase (set)' 'tag' 'plotphase' } {} ...
			   ...
			   { 'Style', 'text', 'string', 'Bootstrap significance level (Ex: 0.01 -> 1%)', 'fontweight', 'bold' } ...
               { 'Style', 'edit', 'string', getkeyval(lastcom,'alpha') 'tag' 'alpha'} ...
               { 'Style', 'checkbox', 'string' 'FDR correct (set)' 'tag' 'fdr' } {} ...
			   ...
			   { 'Style', 'text', 'string', 'Optional newtimef() arguments (see Help)', 'fontweight', 'bold', ...
				 'tooltipstring', 'See newtimef() help via the Help button on the right...' } ...
			   { 'Style', 'edit', 'string', '' 'tag' 'options' } ...
			   {} ...
			   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotersp','present',0), 'string', ...
				 'Plot Event Related Spectral Power', 'tooltipstring', ...
				 'Plot log spectral perturbation image in the upper panel' 'tag' 'plotersp' } ...
			   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotitc','present',0), 'string', ...
				 'Plot Inter Trial Coherence', 'tooltipstring', ...
				 'Plot the inter-trial coherence image in the lower panel' 'tag' 'plotitc' } ...
			   { 'Style', 'checkbox', 'value', 0, 'string', ...
				 'Plot curve at each frequency' 'tag' 'plotcurve' } ...
			 };
      % { 'Style', 'edit', 'string', '''padratio'', 4, ''plotphase'', ''off''' } ...
			   %{ 'Style', 'text', 'string',  '[set] -> Plot ITC phase sign', 'fontweight', 'bold', ...
			%	 'tooltipstring', ['Plot the sign (+/-) of inter-trial coherence phase' 10 ...
			%		'as red (+) or blue (-)'] } ...
			%   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotphase','present',1) } { } ...

	[ tmp1 tmp2 strhalt result ] = inputgui( geometry, uilist, 'pophelp(''pop_newtimef'');', ...
					   fastif(typeproc, 'Plot channel time frequency -- pop_newtimef()', ...
							  'Plot component time frequency -- pop_newtimef()'));
	if length( tmp1 ) == 0 return; end;

	if result.fft,      result.cycle = '0'; end;
	if result.nobase,   result.baseline = 'NaN'; end;

	num	     = eval( [ '[' result.chan    ']' ] );
	tlimits	 = eval( [ '[' result.tlimits ']' ] );
	cycles	 = eval( [ '[' result.cycle   ']' ] );
    freqs    = eval( [ '[' result.freqs   ']' ] );
    %result.ncycles == 2 is ignored

    % add topoplot
    % ------------
    options = [];
    if isfield(EEG.chanlocs, 'theta')
        if ~isfield(EEG, 'chaninfo'), EEG.chaninfo = []; end;
        if typeproc == 1
            options = [options ', ''topovec'', ' int2str(num) ...
                        ', ''elocs'', EEG.chanlocs, ''chaninfo'', EEG.chaninfo' ];
        else
            options = [options ', ''topovec'', EEG.icawinv(:,' int2str(num) ...
                       '), ''elocs'', EEG.chanlocs, ''chaninfo'', EEG.chaninfo' ];
      end;
    end;

	if ~isempty( result.baseline ),  options = [ options ', ''baseline'',[' result.baseline ']' ]; end;
    if ~isempty( result.alpha ),     options = [ options ', ''alpha'',' result.alpha ];   end;
	if ~isempty( result.options ),   options = [ options ',' result.options ];            end;
	if ~isempty( result.freqs ),     options = [ options ', ''freqs'', [' result.freqs ']'   ]; end;
	if ~isempty( result.erspmax ),   options = [ options ', ''erspmax'', [' result.erspmax ']' ]; end;
	if ~isempty( result.itcmax ),    options = [ options ', ''itcmax'','  result.itcmax ];      end;
	if ~result.plotersp,             options = [ options ', ''plotersp'', ''off''' ];     end;
	if ~result.plotitc,              options = [ options ', ''plotitc'' , ''off''' ];     end;
	if result.plotcurve,             options = [ options ', ''plottype'', ''curve''' ];   end;
	if result.fdr,                   options = [ options ', ''mcorrect'', ''fdr''' ];     end;
	if result.freqscale,             options = [ options ', ''freqscale'', ''log''' ];    end;
	if ~result.plotphase,            options = [ options ', ''plotphase'', ''off''' ];    end;
	if ~result.scale,                options = [ options ', ''scale'', ''abs''' ];        end;
    if result.basenorm == 2,         options = [ options ', ''basenorm'', ''on''' ];      end;
    if result.ntimesout == 1,        options = [ options ', ''ntimesout'', 50' ];         end;
    if result.ntimesout == 2,        options = [ options ', ''ntimesout'', 100' ];        end;
    if result.ntimesout == 3,        options = [ options ', ''ntimesout'', 150' ];        end;
    if result.ntimesout == 5,        options = [ options ', ''ntimesout'', 300' ];        end;
    if result.ntimesout == 6,        options = [ options ', ''ntimesout'', 400' ];        end;
    if result.nfreqs == 1,           options = [ options ', ''padratio'', 1' ];           end;
    if result.nfreqs == 2,           options = [ options ', ''padratio'', 2' ];           end;
    if result.nfreqs == 3,           options = [ options ', ''padratio'', 4' ];           end;
    if result.nfreqs == 4,           options = [ options ', ''nfreqs'', ' int2str(length(freqs)) ]; end;

    % add title
    % ---------
	if isempty( findstr(  '''title''', result.options))
        if ~isempty(EEG.chanlocs) & typeproc
            chanlabel = EEG.chanlocs(num).labels;
        else
            chanlabel = int2str(num);
        end;
	end;

	figure; try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end;
else
    options = [ ',' vararg2str(varargin) ];
end;

% compute epoch limits
% --------------------
if isempty(tlimits)
	tlimits = [EEG.xmin, EEG.xmax]*1000;
end;
pointrange1 = round(max((tlimits(1)/1000-EEG.xmin)*EEG.srate, 1));
pointrange2 = round(min((tlimits(2)/1000-EEG.xmin)*EEG.srate, EEG.pnts));
pointrange = [pointrange1:pointrange2];

% call function sample either on raw data or ICA data
% ---------------------------------------------------
if typeproc == 1
	tmpsig = EEG.data(num,pointrange,:);
else
	if ~isempty( EEG.icasphere )
        eeglab_options; % changed from eeglaboptions 3/30/02 -sm
 	    if option_computeica
    		tmpsig = EEG.icaact(num,pointrange,:);
 	    else
            tmpsig = (EEG.icaweights(num,:)*EEG.icasphere)*reshape(EEG.data(:,pointrange,:), EEG.nbchan, EEG.trials*length(pointrange));
        end;
	else
		error('You must run ICA first');
	end;
end;
tmpsig = reshape( tmpsig, length(num), size(tmpsig,2)*size(tmpsig,3));

% outputs
% -------
outstr = '';
if ~popup
    for io = 1:nargout, outstr = [outstr 'varargout{' int2str(io) '},' ]; end;
    if ~isempty(outstr), outstr = [ '[' outstr(1:end-1) '] =' ]; end;
end;

% plot the datas and generate output command
% --------------------------------------------
if length( options ) < 2
    options = '';
end;
if nargin < 4
    varargout{1} = sprintf('figure; pop_newtimef( %s, %d, %d, [%s], [%s] %s);', inputname(1), typeproc, num, ...
			int2str(tlimits), num2str(cycles), options);
end;
com = sprintf('%s newtimef( tmpsig(:, :), length(pointrange), [tlimits(1) tlimits(2)], EEG.srate, cycles %s);', outstr, options);
eval(com)

return;

% get contextual help
% -------------------
function txt = context(var, allvars, alltext);
	loc = strmatch( var, allvars);
	if ~isempty(loc)
		txt= alltext{loc(1)};
	else
		disp([ 'warning: variable ''' var ''' not found']);
		txt = '';
	end;

author: PatternRecognition

repo: OpenBMI


text:
function G = gsp_graph_default_parameters(G)
%GSP_GRAPH_DEFAULT_PARAMETERS load default parameters for graphs
%   Usage: G = gsp_graph_default_parameters(G);
%          G = gsp_graph_default_parameters();
%
%   Input parameters
%       G   : Graph (Optional)
%   Output parameters
%       G   : Graph
%
%   This function will fill a graph with all missing parameters such that
%   it is compabatible with all functions of the GSPBox. If you create a
%   graph manually, you need to set only the weight matrix *W*. If you have
%   some coordonate, you can also set *G.coords*. *G.coords* is a $N$ x $2$
%   or a $N$ x $3$ matrix with each columns beeing the coordonates in each
%   dimention. Finally, we recommend to set the fiel *G.type* with a name
%   that suits your graph.
%
%   Example::
%
%          W = rand(30);
%          W = (W + W')/2;
%          G.W = W - diag(diag(W));
%          G = gsp_graph_default_parameters(G)
%
%   This function can be used to update the weights of your graph. It will
%   recompute the Laplacian operator. Warning this function does not
%   perform any change to the Fourier basis::
%
%          G.W = Wnew;
%          G = gsp_graph_default_parameters(G);
%
%
%   List of parameters of the graph structure
%   -----------------------------------------
%
%   By default, a graph structure in the GSPbox contains the following
%   parameters:
%
%   * *G.W*: Weight matrix (empty by default)
%   * *G.A*: Adacency matrix (constructed with *W*)
%   * *G.N*: Number of nodes (`size(W,1)`)
%   * *G.type*: Type of graph ('unknown' by default)
%   * *G.directed*: 1 if the graph is directed, 0 if not
%   * *G.lap_type*: Laplacian type (default 'combinatorial') See the
%     function |gsp_create_laplacian| for a exhaustive list of the
%     available laplacians.
%   * *G.d*: Degree vector (Computed with *G.W*)
%   * *G.Ne*: Number of edges
%   * *G.coords*: Coordinates of the vertices (default (0,0) )
%   * *G.plotting*: Plotting parameters
%     * *G.plotting.edge_width*: Width of edges (default 1)
%     * *G.plotting.edge_color*: Color of edges (default [255,88,41]/255 )
%     * *G.plotting.edge_style*: Style of edges (default '-')
%     * *G.plotting.vertex_size*: Size of vertex (default 50)
%     * *G.plotting.vertex_color*: Color of vertex (default 'b')
%
%   Remark: There is redudancy between $A$, $W$, $L$. However, the GSPBox
%   is done in matlab and is not suppose yet to scale to graph sufficiently
%   large that your matlab have memory problem. However, this will be most
%   likely change for milestone 1.0.0. If you do have a urgent need to
%   overcome this problem, please contact the devolper team.
%

% Author: Nathanael Perraudin
% Date  : 09.12.2013

if nargin<1
    G=struct;
end

if ~isfield(G,'W') % Weight matrix
    G.W = sparse(0);

end

G.A=sparse(G.W>0);
G.N = size(G.W,1);

% Type of graph
if ~isfield(G,'type')
    G.type='unknown';
end;

% if ~isfield(G,'directed')
    G.directed = gsp_isdirected(G);
% end;

if ~isfield(G,'hypergraph')
    G.hypergraph = 0;
end;

% Create the graph Laplacian
if ~isfield(G,'lap_type')
    G.lap_type='combinatorial';
end

% if ~isfield(G,'L')
    G = gsp_create_laplacian(G);
% end

G.d = full(sum(G.W,2));

% Number of edges
if G.directed
    G.Ne = nnz(G.W);
else
    G.Ne = nnz(G.W)/2;
end

if ~isfield(G,'coords') % Coordonates
    G.coords = [];
end

G = gsp_graph_default_plotting_parameters(G);

end

repo: gspbox


text:
%GETDEFAULTNEWCAMERAMATRIX  Returns the default new camera matrix
%
%     newCameraMatrix = cv.getDefaultNewCameraMatrix(cameraMatrix)
%     [...] = cv.getDefaultNewCameraMatrix(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __cameraMatrix__ Input camera matrix, a 3x3 double matrix.
%
% ## Output
% * __newCameraMatrix__ Output camera matrix, same size and type as input.
%
% ## Options
% * __ImgSize__ Camera view image size in pixels `[w,h]`. Default [0,0]
% * __CenterPrincipalPoint__ Location of the principal point in the new camera
%   matrix. The parameter indicates whether this location should be at the
%   image center or not. default false
%
% The function returns the camera matrix that is either an exact copy of the
% input `cameraMatrix` (when `CenterPrinicipalPoint=false`), or the modified
% one (when `CenterPrincipalPoint=true`).
%
% In the latter case, the new camera matrix will be:
%
%     [ fx,   0,  (ImgSize(1)-1)*0.5 ;
%        0,  fy,  (ImgSize(2)-1)*0.5 ;
%        0,   0,                   1 ]
%
% where `fx` and `fy` are (0,0) and (1,1) elements of `cameraMatrix`,
% respectively.
%
% By default, the undistortion functions in OpenCV (see cv.undistort,
% cv.initUndistortRectifyMap) do not move the principal point. However, when
% you work with stereo, it is important to move the principal points in both
% views to the same y-coordinate (which is required by most of stereo
% correspondence algorithms), and may be to the same x-coordinate too. So, you
% can form the new camera matrix for each view where the principal points are
% located at the center.
%
% See also: cv.undistort, cv.initUndistortRectifyMap,
%  cv.getOptimalNewCameraMatrix
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Script to scale a phase map so that max = pi and min =-pi radians.
% Writes out scaled image prepended with 'sc'.
% Chloe Hutton 25/02/04
% SPM Update - 13/11/06
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Chloe Hutton
% $Id: pm_scale_phase.m 4446 2011-08-30 10:50:29Z guillaume $

V   = spm_vol(spm_select(1,'image','Select phase image to scale'));
vol = spm_read_vols(V);

mn   = min(vol(:));
mx   = max(vol(:));
svol = -pi+(vol-mn)*2*pi/(mx-mn);

% Output image struct
oV = struct(...
    'fname',   spm_file(V.fname,'prefix','sc'),...
    'dim',     V.dim(1:3),...
    'dt',      [4 spm_platform('bigend')],...
    'mat',     V.mat,...
    'descrip', 'Scaled phase');

spm_write_vol(oV,svol);

author: spm


text:
function [ stg_Ld, stg_Rd, att ] = filter_damping_design( stg_L, stg_C, fs, Zmax)
% function [ stg_Ld, stg_Rd, att ] = filter_damping_design( stg_L, stg_C, fs, Zmax)
% Design the optimal damping for a passive filter in power electronics applications.
% Written by Dr. Yoash Levron
%
% This function designs the optimal L-R damping network for a multistage passive filter.
% It compues the optimal damping network, that achieves the best
% attenuation at a desired switching frequency fs, and has an output impedance
% smaller than a specified threshold Zmax, over all frequencies.
% the number of stages may be 1,2 or 3.
% The damping inductors, if placed at all, are chosen to be equal,
% so practical filters can be built using the same inductor, reducing cost.
% The attached photo shows the topology and component names.
%
% Input:
% stg_L - [H] specifies the primary inductors at each stage.
%                    This is a vector with 1, 2 or 3 components.
% stg_C - [F] specifies the primary capacitors at each stage.
%                     It is the same length as stg_L.
% fs - [Hz] desired switching frequency to attenuate.
% Zmax - [ohm] desired maximal output impedance (over all frequencies).
%
% Output:
% stg_Ld - [H] the resulting damping inductors.
% stg_Rd - [ohm] the resulting damping resistors.
% att - [dB] the resulting optimal attenuation at frequency fs

% numeric parameters
NLd = 10; % number of damping inductors to compare
NRd = 20;  % number of damping resistors to compare
min_resistor_factor = 0.2;
max_resistor_factor = 20;

% default outputs
stg_Ld = NaN;
stg_Rd = NaN;
att = NaN;
Zmax_out = NaN;

% check inputs:
if (length(stg_L) == length(stg_C))
    NS = length(stg_L);  % number of stages
else
    disp('filter_damping_design - error : unequal number of inductors and capacitors');
    beep; return
end
if (NS >3)
    disp('filter_damping_design - error : Three stages maximum');
    beep; return
end

Ltot = sum(stg_L);
Ctot = sum(stg_C);
Rnom = (Ltot/Ctot)^0.5;
Ldvec = linspace((max(stg_L)/NLd),max(stg_L),NLd);
Rdvec = logspace( log10(min_resistor_factor*Rnom), ...
    log10(max_resistor_factor*Rnom),  (NRd-1));
Rdvec(NRd) = 10e6;  % add an infinite resistor

ws = 2*pi*fs;

%%%% design for a first order filter %%%%
if (NS ==1)
    L1 =  stg_L(1);
    C1 =  stg_C(1);
    Ldvec = [0 Ldvec];   NLd=NLd+1;  % add a zero to the possible inductor values
    MAT = zeros(NLd,NRd);  % matrix of results
    for ii1 = 1:NLd  % index for damping inductor
        xx = (ii1-1) / NLd;
        done = floor(1000*xx)/10;
        disp(sprintf('%3.1f%% done',done));

        for ii2 = 1:NRd % index for 1st damping resistor
            % determine components
            Ld1 = Ldvec(ii1);
            Rd1 = Rdvec(ii2);

            % compute coefficients
            a1 = (L1 + Ld1)/Rd1;
            d1 = L1;
            d2 = (L1*Ld1)/Rd1;
            b1 = (L1 + Ld1)/Rd1;
            b2 =  C1*L1;
            b3 = (C1*L1*Ld1)/Rd1;

            % compute the resonance frequencies - vector Wr
            Wr = [(1/b2) (b1/b3)].^0.5;
            ind = find(~isinf(Wr));
            Wr = Wr(ind);

            % compute the impedace at resonance
            s = i*Wr;
            Zout = (d1*s + d2*s.^2) ./ (1 + b1*s + b2*s.^2 + b3*s.^3);
            ind = find(abs(Zout) > Zmax);
            output_impedance_too_big = ~isempty(ind);

            % compute the attenuation
            s = i*ws;
            H = (1 + a1*s) / (1 + b1*s + b2*s^2 + b3*s^3);
            H_dB = 20*log10(abs(H));

            MAT(ii1,ii2) = H_dB + 1e10*output_impedance_too_big;
        end
    end
    %%% find the best filter
    [att, ind] = min(MAT(:));
    if (att > 1e9)
        disp('filter_damping_design - error : Design failed. Output impedance too high');
        beep; return;
    end
    [ii1,ii2] = ind2sub(size(MAT),ind);  % indexes of optimal filter
    Ld1 = Ldvec(ii1);
    Rd1 = Rdvec(ii2);
    stg_Ld = [Ld1];
    stg_Rd = [Rd1];
end

%%%% design for a second order filter %%%%
if (NS ==2)
    disp('This may take a few minutes. please wait ...');
    L1 =  stg_L(1);  L2 = stg_L(2);
    C1 =  stg_C(1);  C2 = stg_C(2);
    MAT = zeros(NLd,NRd,NRd,2,2);  % matrix of results
    for ii1 = 1:NLd  % index for damping inductor
        for ii2 = 1:NRd % index for 1st damping resistor
            xx = ((ii1-1)*NRd+(ii2-1)) / (NLd*NRd);
            done = floor(1000*xx)/10;
            disp(sprintf('%3.1f%% done',done));

            for ii3 = 1:NRd % index for 2st damping resistor
                for ii5 = 1:2  % index for the existance of 1st damping inductor
                    for ii6 = 1:2  % index for the existance of 2st damping inductor
                        % determine components
                        Ld1 = Ldvec(ii1) * (ii5-1);
                        Ld2 = Ldvec(ii1) * (ii6-1);
                        Rd1 = Rdvec(ii2);
                        Rd2 = Rdvec(ii3);

                        % compute coefficients
                        a1 = (L1 + Ld1)/Rd1 + (L2 + Ld2)/Rd2;
                        a2 = ((L1 + Ld1)*(L2 + Ld2))/(Rd1*Rd2);
                        d1 = L1 + L2;
                        d2 = (L1*L2 + L1*Ld1 + L2*Ld1)/Rd1 + (L1*L2 + L1*Ld2 + L2*Ld2)/Rd2;
                        d3 = (L1*L2*Ld1 + L1*L2*Ld2 + L1*Ld1*Ld2 + L2*Ld1*Ld2)/(Rd1*Rd2) + C1*L1*L2;
                        d4 = (C1*L1*L2*(Ld1*Rd2 + Ld2*Rd1))/(Rd1*Rd2);
                        d5 = (C1*L1*L2*Ld1*Ld2)/(Rd1*Rd2);
                        b1 = (L1 + Ld1)/Rd1 + (L2 + Ld2)/Rd2;
                        b2 = C1*L1 + C2*L1 + C2*L2 + ((L1 + Ld1)*(L2 + Ld2))/(Rd1*Rd2);
                        b3 = (C2*L1*L2 + C1*L1*Ld1 + C2*L1*Ld1 + C2*L2*Ld1)/Rd1 + ...
                            (C1*L1*L2 + C2*L1*L2 + C1*L1*Ld2 + C2*L1*Ld2 + C2*L2*Ld2)/Rd2;
                        b4 = (C1*L1*L2*Ld1 + C2*L1*L2*Ld1 + C2*L1*L2*Ld2 + C1*L1*Ld1*Ld2 ...
                            + C2*L1*Ld1*Ld2 + C2*L2*Ld1*Ld2)/(Rd1*Rd2) + C1*C2*L1*L2;
                        b5 = (C1*C2*L1*L2*(Ld1*Rd2 + Ld2*Rd1))/(Rd1*Rd2);
                        b6 = (C1*C2*L1*L2*Ld1*Ld2)/(Rd1*Rd2);

                        % compute the resonance frequencies - vector Wr
                        P1 = [-b6, +b4, -b2, +1];
                        P2 = [+b5, -b3, +b1];
                        WrA = (roots(P1)).^0.5;
                        WrB = (roots(P2)).^0.5;
                        ind1 = find(real(WrA) > 1e6*imag(WrA) );
                        ind2 = find(real(WrB) > 1e6*imag(WrB));
                        Wr = [real(WrA(ind1)).' real(WrB(ind2)).'];

                        % compute the impedace at resonance
                        s = i*Wr;
                        Zout = (d1*s + d2*s.^2 + d3*s.^3 + d4*s.^4 + d5*s.^5) ./ ...
                            (1 + b1*s + b2*s.^2 + b3*s.^3 + b4*s.^4 + b5*s.^5 + b6*s.^6);
                        ind = find(abs(Zout) > Zmax);
                        output_impedance_too_big = ~isempty(ind);

                        % compute the attenuation
                        s = i*ws;
                        H = (1 + a1*s + a2*s^2) / ...
                            (1 + b1*s + b2*s^2 + b3*s^3 + b4*s^4 + b5*s^5 + b6*s^6);
                        H_dB = 20*log10(abs(H));

                        MAT(ii1,ii2,ii3,ii5,ii6) = H_dB + 1e10*output_impedance_too_big;
                    end
                end
            end
        end
    end
    %%% find the best filter
    [att, ind] = min(MAT(:));
    if (att > 1e9)
        disp('filter_damping_design - error : Design failed. Output impedance too high');
        beep; return;
    end
    [ii1,ii2,ii3,ii5,ii6] = ind2sub(size(MAT),ind);  % indexes of optimal filter
    Ld1 = Ldvec(ii1) * (ii5-1);       Ld2 = Ldvec(ii1) * (ii6-1);
    Rd1 = Rdvec(ii2);       Rd2 = Rdvec(ii3);
    stg_Ld = [Ld1 Ld2];
    stg_Rd = [Rd1 Rd2];
end

%%%% design for a third order filter %%%%
if (NS == 3)
    disp('This may take a few minutes. please wait ...');
    L1 =  stg_L(1);  L2 = stg_L(2);  L3 = stg_L(3);
    C1 =  stg_C(1);  C2 = stg_C(2);  C3 = stg_C(3);
    MAT = zeros(NLd,NRd,NRd,NRd,2,2,2);  % matrix of results
    for ii1 = 1:NLd  % index for damping inductor
        for ii2 = 1:NRd % index for 1st damping resistor
            xx = ((ii1-1)*NRd+(ii2-1)) / (NLd*NRd);
            done = floor(1000*xx)/10;
            disp(sprintf('%3.1f%% done',done));

            for ii3 = 1:NRd % index for 2st damping resistor
                for ii4 = 1:NRd % index for 3rd damping resistor
                    for ii5 = 1:2  % index for the existance of 1st damping inductor
                        for ii6 = 1:2  % index for the existance of 2st damping inductor
                            for ii7 = 1:2  % index for the existance of 3rd damping inductor
                                % determine components
                                Ld1 = Ldvec(ii1) * (ii5-1);
                                Ld2 = Ldvec(ii1) * (ii6-1);
                                Ld3 = Ldvec(ii1) * (ii7-1);
                                Rd1 = Rdvec(ii2);
                                Rd2 = Rdvec(ii3);
                                Rd3 = Rdvec(ii4);

                                % compute coefficients
                                a1 = (L1 + Ld1)/Rd1 + (L2 + Ld2)/Rd2 + (L3 + Ld3)/Rd3;
                                a2 = ((L1 + Ld1)/(Rd1*Rd3) + (L2 + Ld2)/(Rd2*Rd3))*(L3 + Ld3) +...
                                    ((L1 + Ld1)*(L2 + Ld2))/(Rd1*Rd2);
                                a3 = ((L1 + Ld1)*(L2 + Ld2)*(L3 + Ld3))/(Rd1*Rd2*Rd3);
                                d1 = L1 + L2 + L3;
                                d2 = (L1*L2*Rd1 + L1*L2*Rd2 + L1*L3*Rd2 + L2*L3*Rd1 + ...
                                    L1*Ld1*Rd2 + L1*Ld2*Rd1 + L2*Ld1*Rd2 + L2*Ld2*Rd1 + ...
                                    L3*Ld1*Rd2 + L3*Ld2*Rd1)/(Rd1*Rd2) +...
                                    (L1*L3*Rd1*Rd2 + L2*L3*Rd1*Rd2 + L1*Ld3*Rd1*Rd2 + ...
                                    L2*Ld3*Rd1*Rd2 + L3*Ld3*Rd1*Rd2)/(Rd1*Rd2*Rd3);
                                d3 =  (L1*L2*L3 + L1*L2*Ld1 + L1*L2*Ld2 + L1*L3*Ld2 + ...
                                    L2*L3*Ld1 + L1*Ld1*Ld2 + L2*Ld1*Ld2 + L3*Ld1*Ld2 + ...
                                    C1*L1*L2*Rd1*Rd2 + C1*L1*L3*Rd1*Rd2 + C2*L1*L3*Rd1*Rd2 + ...
                                    C2*L2*L3*Rd1*Rd2)/(Rd1*Rd2) + (L1*L2*L3*Rd1 + L1*L2*L3*Rd2 + ...
                                    L1*L2*Ld3*Rd1 + L1*L3*Ld1*Rd2 + L1*L3*Ld2*Rd1 + L1*L2*Ld3*Rd2 +...
                                    L2*L3*Ld1*Rd2 + L2*L3*Ld2*Rd1 + L1*L3*Ld3*Rd2 + L2*L3*Ld3*Rd1 +...
                                    L1*Ld1*Ld3*Rd2 + L1*Ld2*Ld3*Rd1 + L2*Ld1*Ld3*Rd2 + L2*Ld2*Ld3*Rd1 +...
                                    L3*Ld1*Ld3*Rd2 + L3*Ld2*Ld3*Rd1)/(Rd1*Rd2*Rd3);
                                d4 = (C1*L1*L2*L3*Rd1 + C2*L1*L2*L3*Rd1 + C2*L1*L2*L3*Rd2 + ...
                                    C1*L1*L2*Ld1*Rd2 + C1*L1*L2*Ld2*Rd1 + C1*L1*L3*Ld1*Rd2 + ...
                                    C1*L1*L3*Ld2*Rd1 + C2*L1*L3*Ld1*Rd2 + C2*L1*L3*Ld2*Rd1 + ...
                                    C2*L2*L3*Ld1*Rd2 + C2*L2*L3*Ld2*Rd1)/(Rd1*Rd2) +...
                                    (L1*L2*L3*Ld1 + L1*L2*L3*Ld2 + L1*L2*L3*Ld3 + L1*L2*Ld1*Ld3 +...
                                    L1*L3*Ld1*Ld2 + L1*L2*Ld2*Ld3 + L2*L3*Ld1*Ld2 + L1*L3*Ld2*Ld3 +...
                                    L2*L3*Ld1*Ld3 + L1*Ld1*Ld2*Ld3 + L2*Ld1*Ld2*Ld3 + L3*Ld1*Ld2*Ld3 +...
                                    C1*L1*L2*L3*Rd1*Rd2 + C1*L1*L2*Ld3*Rd1*Rd2 + C1*L1*L3*Ld3*Rd1*Rd2 +...
                                    C2*L1*L3*Ld3*Rd1*Rd2 + C2*L2*L3*Ld3*Rd1*Rd2)/(Rd1*Rd2*Rd3);
                                d5 =  (C1*L1*L2*L3*Ld1 + C2*L1*L2*L3*Ld1 + ...
                                    C2*L1*L2*L3*Ld2 + C1*L1*L2*Ld1*Ld2 + C1*L1*L3*Ld1*Ld2 +...
                                    C2*L1*L3*Ld1*Ld2 + C2*L2*L3*Ld1*Ld2 + ...
                                    C1*C2*L1*L2*L3*Rd1*Rd2)/(Rd1*Rd2) + (C1*L1*L2*L3*Ld1*Rd2 +...
                                    C1*L1*L2*L3*Ld2*Rd1 + C1*L1*L2*L3*Ld3*Rd1 + C2*L1*L2*L3*Ld3*Rd1 +...
                                    C2*L1*L2*L3*Ld3*Rd2 + C1*L1*L2*Ld1*Ld3*Rd2 + C1*L1*L2*Ld2*Ld3*Rd1 +...
                                    C1*L1*L3*Ld1*Ld3*Rd2 + C1*L1*L3*Ld2*Ld3*Rd1 + C2*L1*L3*Ld1*Ld3*Rd2 +...
                                    C2*L1*L3*Ld2*Ld3*Rd1 + C2*L2*L3*Ld1*Ld3*Rd2 + ...
                                    C2*L2*L3*Ld2*Ld3*Rd1)/(Rd1*Rd2*Rd3);
                                d6 = (C1*L1*L2*L3*Ld1*Ld2 + C1*L1*L2*L3*Ld1*Ld3 + ...
                                    C2*L1*L2*L3*Ld1*Ld3 + C2*L1*L2*L3*Ld2*Ld3 + C1*L1*L2*Ld1*Ld2*Ld3 +...
                                    C1*L1*L3*Ld1*Ld2*Ld3 + C2*L1*L3*Ld1*Ld2*Ld3 + C2*L2*L3*Ld1*Ld2*Ld3 +...
                                    C1*C2*L1*L2*L3*Ld3*Rd1*Rd2)/(Rd1*Rd2*Rd3) +...
                                    (C1*C2*L1*L2*L3*(Ld1*Rd2 + Ld2*Rd1))/(Rd1*Rd2);
                                d7 =  (C1*C2*L1*L2*L3*(Ld1*Ld2*Rd3 + Ld1*Ld3*Rd2 +...
                                    Ld2*Ld3*Rd1))/(Rd1*Rd2*Rd3);
                                d8 = (C1*C2*L1*L2*L3*Ld1*Ld2*Ld3)/(Rd1*Rd2*Rd3);
                                b1 = (L1 + Ld1)/Rd1 + (L2*Rd3 + L3*Rd2 + Ld2*Rd3 + Ld3*Rd2)/(Rd2*Rd3);
                                b2 = (L2*L3 + L2*Ld3 + L3*Ld2 + Ld2*Ld3 + C1*L1*Rd2*Rd3 + C2*L1*Rd2*Rd3 +...
                                    C2*L2*Rd2*Rd3 + C3*L1*Rd2*Rd3 + C3*L2*Rd2*Rd3 + ...
                                    C3*L3*Rd2*Rd3)/(Rd2*Rd3) + ((L1 + Ld1)*(L2*Rd3 +...
                                    L3*Rd2 + Ld2*Rd3 + Ld3*Rd2))/(Rd1*Rd2*Rd3);
                                b3 =  (C1*L1*L2*Rd3 + C1*L1*L3*Rd2 + C2*L1*L2*Rd3 + C2*L1*L3*Rd2 +...
                                    C2*L2*L3*Rd2 + C3*L1*L2*Rd3 + C3*L1*L3*Rd2 + C3*L2*L3*Rd2 +...
                                    C3*L2*L3*Rd3 + C1*L1*Ld2*Rd3 + C1*L1*Ld3*Rd2 + C2*L1*Ld2*Rd3 +...
                                    C2*L1*Ld3*Rd2 + C2*L2*Ld2*Rd3 + C2*L2*Ld3*Rd2 + C3*L1*Ld2*Rd3 +...
                                    C3*L1*Ld3*Rd2 + C3*L2*Ld2*Rd3 + C3*L2*Ld3*Rd2 + C3*L3*Ld2*Rd3 +...
                                    C3*L3*Ld3*Rd2)/(Rd2*Rd3) + (L1*L2*L3 + L1*L2*Ld3 + L1*L3*Ld2 +...
                                    L2*L3*Ld1 + L1*Ld2*Ld3 + L2*Ld1*Ld3 + L3*Ld1*Ld2 + Ld1*Ld2*Ld3 +...
                                    C2*L1*L2*Rd2*Rd3 + C3*L1*L2*Rd2*Rd3 + C3*L1*L3*Rd2*Rd3 +...
                                    C1*L1*Ld1*Rd2*Rd3 + C2*L1*Ld1*Rd2*Rd3 + C2*L2*Ld1*Rd2*Rd3 +...
                                    C3*L1*Ld1*Rd2*Rd3 + C3*L2*Ld1*Rd2*Rd3 +...
                                    C3*L3*Ld1*Rd2*Rd3)/(Rd1*Rd2*Rd3);
                                b4 = (C1*L1*L2*L3 + C2*L1*L2*L3 + C3*L1*L2*L3 + C1*L1*L2*Ld3 +...
                                    C1*L1*L3*Ld2 + C2*L1*L2*Ld3 + C2*L1*L3*Ld2 + C2*L2*L3*Ld2 +...
                                    C3*L1*L2*Ld3 + C3*L1*L3*Ld2 + C3*L2*L3*Ld2 + C3*L2*L3*Ld3 +...
                                    C1*L1*Ld2*Ld3 + C2*L1*Ld2*Ld3 + C2*L2*Ld2*Ld3 + C3*L1*Ld2*Ld3 +...
                                    C3*L2*Ld2*Ld3 + C3*L3*Ld2*Ld3 + C1*C2*L1*L2*Rd2*Rd3 +...
                                    C1*C3*L1*L2*Rd2*Rd3 + C1*C3*L1*L3*Rd2*Rd3 + C2*C3*L1*L3*Rd2*Rd3 +...
                                    C2*C3*L2*L3*Rd2*Rd3)/(Rd2*Rd3) + (C2*L1*L2*L3*Rd2 + C3*L1*L2*L3*Rd2 +...
                                    C3*L1*L2*L3*Rd3 + C1*L1*L2*Ld1*Rd3 + C1*L1*L3*Ld1*Rd2 +...
                                    C2*L1*L2*Ld1*Rd3 + C2*L1*L3*Ld1*Rd2 + C2*L1*L2*Ld2*Rd3 + ...
                                    C2*L1*L2*Ld3*Rd2 + C2*L2*L3*Ld1*Rd2 + C3*L1*L2*Ld1*Rd3 +...
                                    C3*L1*L3*Ld1*Rd2 + C3*L1*L2*Ld2*Rd3 + C3*L1*L2*Ld3*Rd2 + ...
                                    C3*L2*L3*Ld1*Rd2 + C3*L1*L3*Ld2*Rd3 + C3*L1*L3*Ld3*Rd2 + ...
                                    C3*L2*L3*Ld1*Rd3 + C1*L1*Ld1*Ld2*Rd3 + C1*L1*Ld1*Ld3*Rd2 + ...
                                    C2*L1*Ld1*Ld2*Rd3 + C2*L1*Ld1*Ld3*Rd2 + C2*L2*Ld1*Ld2*Rd3 + ...
                                    C2*L2*Ld1*Ld3*Rd2 + C3*L1*Ld1*Ld2*Rd3 + C3*L1*Ld1*Ld3*Rd2 + ...
                                    C3*L2*Ld1*Ld2*Rd3 + C3*L2*Ld1*Ld3*Rd2 + C3*L3*Ld1*Ld2*Rd3 + ...
                                    C3*L3*Ld1*Ld3*Rd2)/(Rd1*Rd2*Rd3);
                                b5 = (C1*C2*L1*L2*L3*Rd2 + C1*C3*L1*L2*L3*Rd2 + C1*C3*L1*L2*L3*Rd3 + ...
                                    C2*C3*L1*L2*L3*Rd3 + C1*C2*L1*L2*Ld2*Rd3 + C1*C2*L1*L2*Ld3*Rd2 + ...
                                    C1*C3*L1*L2*Ld2*Rd3 + C1*C3*L1*L2*Ld3*Rd2 + C1*C3*L1*L3*Ld2*Rd3 + ...
                                    C1*C3*L1*L3*Ld3*Rd2 + C2*C3*L1*L3*Ld2*Rd3 + C2*C3*L1*L3*Ld3*Rd2 + ...
                                    C2*C3*L2*L3*Ld2*Rd3 + C2*C3*L2*L3*Ld3*Rd2)/(Rd2*Rd3) + ...
                                    (C1*L1*L2*L3*Ld1 + C2*L1*L2*L3*Ld1 + C2*L1*L2*L3*Ld2 + C3*L1*L2*L3*Ld1 +...
                                    C3*L1*L2*L3*Ld2 + C3*L1*L2*L3*Ld3 + C1*L1*L2*Ld1*Ld3 + C1*L1*L3*Ld1*Ld2 +...
                                    C2*L1*L2*Ld1*Ld3 + C2*L1*L3*Ld1*Ld2 + C2*L1*L2*Ld2*Ld3 + ...
                                    C2*L2*L3*Ld1*Ld2 + C3*L1*L2*Ld1*Ld3 + C3*L1*L3*Ld1*Ld2 + ...
                                    C3*L1*L2*Ld2*Ld3 + C3*L2*L3*Ld1*Ld2 + C3*L1*L3*Ld2*Ld3 + ...
                                    C3*L2*L3*Ld1*Ld3 + C1*L1*Ld1*Ld2*Ld3 + C2*L1*Ld1*Ld2*Ld3 + ...
                                    C2*L2*Ld1*Ld2*Ld3 + C3*L1*Ld1*Ld2*Ld3 + C3*L2*Ld1*Ld2*Ld3 + ...
                                    C3*L3*Ld1*Ld2*Ld3 + C2*C3*L1*L2*L3*Rd2*Rd3 + C1*C2*L1*L2*Ld1*Rd2*Rd3 +...
                                    C1*C3*L1*L2*Ld1*Rd2*Rd3 + C1*C3*L1*L3*Ld1*Rd2*Rd3 + ...
                                    C2*C3*L1*L3*Ld1*Rd2*Rd3 + C2*C3*L2*L3*Ld1*Rd2*Rd3)/(Rd1*Rd2*Rd3);
                                b6 = (C1*C2*L1*L2*L3*Ld2 + C1*C3*L1*L2*L3*Ld2 + C1*C3*L1*L2*L3*Ld3 + ...
                                    C2*C3*L1*L2*L3*Ld3 + C1*C2*L1*L2*Ld2*Ld3 + C1*C3*L1*L2*Ld2*Ld3 + ...
                                    C1*C3*L1*L3*Ld2*Ld3 + C2*C3*L1*L3*Ld2*Ld3 + C2*C3*L2*L3*Ld2*Ld3 + ...
                                    C1*C2*C3*L1*L2*L3*Rd2*Rd3)/(Rd2*Rd3) + (C1*C2*L1*L2*L3*Ld1*Rd2 + ...
                                    C1*C3*L1*L2*L3*Ld1*Rd2 + C1*C3*L1*L2*L3*Ld1*Rd3 + C2*C3*L1*L2*L3*Ld1*Rd3 +...
                                    C2*C3*L1*L2*L3*Ld2*Rd3 + C2*C3*L1*L2*L3*Ld3*Rd2 + C1*C2*L1*L2*Ld1*Ld2*Rd3 +...
                                    C1*C2*L1*L2*Ld1*Ld3*Rd2 + C1*C3*L1*L2*Ld1*Ld2*Rd3 + ...
                                    C1*C3*L1*L2*Ld1*Ld3*Rd2 + C1*C3*L1*L3*Ld1*Ld2*Rd3 + ...
                                    C1*C3*L1*L3*Ld1*Ld3*Rd2 + C2*C3*L1*L3*Ld1*Ld2*Rd3 + ...
                                    C2*C3*L1*L3*Ld1*Ld3*Rd2 + C2*C3*L2*L3*Ld1*Ld2*Rd3 + ...
                                    C2*C3*L2*L3*Ld1*Ld3*Rd2)/(Rd1*Rd2*Rd3);
                                b7 =  (C1*C2*L1*L2*L3*Ld1*Ld2 + C1*C3*L1*L2*L3*Ld1*Ld2 + ...
                                    C1*C3*L1*L2*L3*Ld1*Ld3 + C2*C3*L1*L2*L3*Ld1*Ld3 + ...
                                    C2*C3*L1*L2*L3*Ld2*Ld3 + C1*C2*L1*L2*Ld1*Ld2*Ld3 + ...
                                    C1*C3*L1*L2*Ld1*Ld2*Ld3 + C1*C3*L1*L3*Ld1*Ld2*Ld3 + ...
                                    C2*C3*L1*L3*Ld1*Ld2*Ld3 + C2*C3*L2*L3*Ld1*Ld2*Ld3 + ...
                                    C1*C2*C3*L1*L2*L3*Ld3*Rd1*Rd2)/(Rd1*Rd2*Rd3) + ...
                                    (C1*C2*C3*L1*L2*L3*(Ld1*Rd2 + Ld2*Rd1))/(Rd1*Rd2);
                                b8 = (C1*C2*C3*L1*L2*L3*(Ld1*Ld2*Rd3 + Ld1*Ld3*Rd2 + ...
                                    Ld2*Ld3*Rd1))/(Rd1*Rd2*Rd3);
                                b9 = (C1*C2*C3*L1*L2*L3*Ld1*Ld2*Ld3)/(Rd1*Rd2*Rd3);

                                % compute the resonance frequencies - vector Wr
                                P1 = [+b8, -b6, +b4, -b2, +1];
                                P2 = [+b9, -b7, +b5, -b3, +b1];
                                WrA = (roots(P1)).^0.5;
                                WrB = (roots(P2)).^0.5;
                                ind1 = find(real(WrA) > 1e6*imag(WrA) );
                                ind2 = find(real(WrB) > 1e6*imag(WrB));
                                Wr = [real(WrA(ind1)).' real(WrB(ind2)).'];

                                % compute the impedace at resonance
                                s = i*Wr;
                                Zout = (d1*s + d2*s.^2 + d3*s.^3 + d4*s.^4 + d5*s.^5 + d6*s.^6 + d7*s.^7 + d8*s.^8) ./ ...
                                    (1 + b1*s + b2*s.^2 + b3*s.^3 + b4*s.^4 + b5*s.^5 + b6*s.^6 + b7*s.^7 + b8*s.^8 + b9*s.^9);
                                ind = find(abs(Zout) > Zmax);
                                output_impedance_too_big = ~isempty(ind);

                                % compute the attenuation
                                s = i*ws;
                                H = (1 + a1*s + a2*s^2 + a3*s^3) / ...
                                    (1 + b1*s + b2*s^2 + b3*s^3 + b4*s^4 + b5*s^5 + b6*s^6 + b7*s^7 + b8*s^8 + b9*s^9);
                                H_dB = 20*log10(abs(H));

                                MAT(ii1,ii2,ii3,ii4,ii5,ii6,ii7) = H_dB + 1e10*output_impedance_too_big;
                            end
                        end
                    end
                end
            end
        end
    end
    %%% find the best filter
    [att, ind] = min(MAT(:));
    if (att > 1e9)
        disp('filter_damping_design - error : Design failed. Output impedance too high');
        beep; return;
    end
    [ii1,ii2,ii3,ii4,ii5,ii6,ii7] = ind2sub(size(MAT),ind);  % indexes of optimal filter
    Ld1 = Ldvec(ii1) * (ii5-1);       Ld2 = Ldvec(ii1) * (ii6-1);       Ld3 = Ldvec(ii1) * (ii7-1);
    Rd1 = Rdvec(ii2);       Rd2 = Rdvec(ii3);       Rd3 = Rdvec(ii4);
    stg_Ld = [Ld1 Ld2 Ld3];
    stg_Rd = [Rd1 Rd2 Rd3];
end

% mark infinite damping resistors
ind = find(stg_Rd > 9e6);
stg_Rd(ind) = inf;
stg_Ld(ind)  = 0;
end

%%% Transfer functions of the filters
% %%% 1 stage %%%
% H = (1 + a1*s) / (1 + b1*s + b2*s^2 + b3*s^3);
% Zout = (d1*s + d2*s^2) / (1 + b1*s + b2*s^2 + b3*s^3);

% %%% 2 stages %%%
% H = (1 + a1*s + a2*s^2) / ...
%     (1 + b1*s + b2*s^2 + b3*s^3 + b4*s^4 + b5*s^5 + b6*s^6);
% Zout = (d1*s + d2*s^2 + d3*s^3 + d4*s^4 + d5*s^5) / ...
%     (1 + b1*s + b2*s^2 + b3*s^3 + b4*s^4 + b5*s^5 + b6*s^6);

% %%% 3 stages %%%
% H = (1 + a1*s + a2*s^2 + a3*s^3) / ...
%     (1 + b1*s + b2*s^2 + b3*s^3 + b4*s^4 + b5*s^5 + b6*s^6 + b7*s^7 + b8*s^8 + b9*s^9);
% Zout = (d1*s + d2*s^2 + d3*s^3 + d4*s^4 + d5*s^5 + d6*s^6 + d7*s^7 + d8*s^8) / ...
%     (1 + b1*s + b2*s^2 + b3*s^3 + b4*s^4 + b5*s^5 + b6*s^6 + b7*s^7 + b8*s^8 + b9*s^9);

author: Sable

repo: mcbench-benchmarks


text:
function [C,PI,PC] = dual(V,F)
  % DUAL Construct the dual polygonal mesh of a given triangle mesh
  %
  % [C,PI,PC] = dual(V,F)
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by 3 list of indices into rows of V
  % Outputs:
  %   C  #C by dim list dual vertex positions
  %   PI  #PI stream of polygon indices into rows of C
  %   PC  #V+1 list of cumulative sum of dual face valences
  %
  % See also: polygons_to_triangles, dual_subdivide
  %
  assert(~any(on_boundary(F),'all'));
  [~,C] = circumradius(V,F);

  % vertex indices
  I = F(:);
  % only keep unique
  [I,J] = unique(I);
  % index in faces
  IF = mod(J-1,size(F,1))+1;
  % order in face
  IC = floor((J-1)/size(F,1))+1;

  [Fp,Fi] = triangle_triangle_adjacency(F);

  p1 = @(I) mod(I,3)+1;

  %NF = Fp(sub2ind(size(Fp),IF,p1(IC)));
  %NC = p1(Fi(sub2ind(size(Fi),IF,p1(IC))));
  %clf;
  %hold on;
  %tsurf(F,V,'FaceColor','w',falpha(0.8,1));
  %qvr(V(I,:),C(IF,:)-V(I,:),0,'LineWidth',2);
  %qvr(C(IF,:),C(NF,:)-C(IF,:),0,'LineWidth',2);
  %qvr(C(NF,:),V(F(sub2ind(size(F),NF,NC)),:)-C(NF,:),0,'LineWidth',2);
  %hold off;
  %axis equal;
  %view(3);

  % starting face
  IF0 = IF;
  % ledger of vertex-face pairs in order of observation
  L = [];
  while true
    L = [L;I IF];
    N = sub2ind(size(Fp),IF,p1(IC));
    NF = Fp(N);
    NC = p1(Fi(N));
    %clf;
    %hold on;
    %tsurf(F,V,'FaceColor','w',falpha(0.8,1));
    %qvr(V(I,:),C(IF,:)-V(I,:),0,'LineWidth',2);
    %qvr(C(IF,:),C(NF,:)-C(IF,:),0,'LineWidth',2);
    %qvr(C(NF,:),V(F(sub2ind(size(F),NF,NC)),:)-C(NF,:),0,'LineWidth',2);
    %hold off;
    %axis equal;
    %view(3);
    %pause
    IF = NF;
    IC = NC;

    keep = find(IF ~= IF0);
    IF = IF(keep);
    IC = IC(keep);
    I = I(keep);
    IF0 = IF0(keep);
    if isempty(keep)
      break;
    end
  end

  % stable sort by vertex
  [~,S] = sort(L(:,1));
  L = L(S,:);

  PC = cumsum([0;accumarray(L(:,1),1)]);
  PI = L(:,2);
end

author: alecjacobson

repo: gptoolbox


text:
function fout = gwaitbar(x,whichbar, varargin)
%WAITBAR Display wait bar.
%   H = WAITBAR(X,'message', property, value, property, value, ...)
%   creates and displays a waitbar of fractional length X.  The
%   handle to the waitbar figure is returned in H.
%   X should be between 0 and 1.  Optional arguments property and
%   value allow to set corresponding waitbar figure properties.
%   Property can also be an action keyword 'CreateCancelBtn', in
%   which case a cancel button will be added to the figure, and
%   the passed value string will be executed upon clicking on the
%   cancel button or the close figure button.
%
%   WAITBAR(X) will set the length of the bar in the most recently
%   created waitbar window to the fractional length X.
%
%   WAITBAR(X,H) will set the length of the bar in waitbar H
%   to the fractional length X.
%
%   WAITBAR(X,H,'message') will update the message text in
%   the waitbar figure, in addition to setting the fractional
%   length to X.
%
%   WAITBAR is typically used inside a FOR loop that performs a
%   lengthy computation.
%
%   Example:
%       h = waitbar(0,'Please wait...');
%       for i=1:1000,
%           % computation here %
%           waitbar(i/1000,h)
%       end
%
%   See also DIALOG, MSGBOX.

%   Copyright 1984-2015 The MathWorks, Inc.

% Error in -nojvm mode
if ~usejava('jvm')
    error(message('MATLAB:HandleGraphics:noJVM'));
end

inWaitbarCreate = true;

if (nargin == 0)
    error(message('MATLAB:waitbar:InvalidArguments'));
elseif (nargin > 0)
    % Must be a numeric value
    if ~isnumeric(x) || ~isscalar(x)
        error(message('MATLAB:waitbar:InvalidFirstInput'));
    elseif ((x < 0) || (x > 1))
        % Throw a warning in this case, clamp it down and keep going
        % this is a behavior change and we want to eventually error out in this scenario,
        % but want to do that gradually
        if (x < 0)
            x = 0;
        elseif (x > 1)
            x = 1;
        end
        % This warning will be enabled when callers no longer send in values
        % outside the allowed range
        % warning('MATLAB:waitbar:invalidValue', '%s\n%s%s', ...
        %    'The first argument must be a numeric value between 0 and 1.',...
        %    'Setting the value to: ', num2str(x));
    end

    if (nargin == 1)
        % A waitbar handle is not provided. Look for one.
        f = findobj(allchild(0),'flat','Tag','TMWWaitbar');

        if isempty(f)
            % No waitbar found. Create a new waitbar with a default string
            name = getString(message('MATLAB:uistring:waitbar:WaitbarLabel'));
        else
            % Found waitbar(s). Update the first (current) waitbar with the new value
            inWaitbarCreate = false;
            f = f(1);
            % Get handles from the saved appdata
            handles = extractHandles(f);
        end
    elseif (nargin > 1)
        % A waitbar message or handle to an existing waitbar has been
        % provided
        if ischar(whichbar) || iscellstr(whichbar)
            % The second arg is a message, we are creating a waitbar
            name = whichbar;
        elseif all(ishghandle(whichbar, 'figure'))
            % The second arg is a handle, we are updating the value
            inWaitbarCreate = false;
            % Get handles from the saved appdata
            handles = extractHandles(whichbar);
        else
            error(message('MATLAB:waitbar:InvalidSecondInput'));
        end
    end
end

x = max(0,min(100*x,100)); % Map any value of x to a value between 0 and 100
try
    if inWaitbarCreate
        % waitbar(x,name)  initialize
        createWaitbar(varargin{:});
    else
        % waitbar(x)    update
        updateWaitbar(varargin{:});
    end  % inWaitbarCreate
catch ex
    err = MException('MATLAB:waitbar:InvalidArguments','%s',...
        getString(message('MATLAB:waitbar:ImproperArguments')));
    err = err.addCause(ex);
    close(findobj(allchild(0),'flat','Tag','TMWWaitbar'));
    throw(err);
end
drawnow;

if nargout==1,
    fout = handles.figure;
end

    function createWaitbar(varargin)
        vertMargin = 0;
        if nargin > 0,
            % we have optional arguments: property-value pairs
            if rem (nargin, 2 ) ~= 0
                error(message('MATLAB:waitbar:InvalidOptionalArgsPass'));
            end
        end

        oldRootUnits = get(0,'Units');
        % Restore the Units safely
        c = onCleanup(@()set(0, 'Units', oldRootUnits));

        set(0, 'Units', 'points');
        screenSize = get(0,'ScreenSize');
        delete(c);

        pointsPerPixel = 72/get(0,'ScreenPixelsPerInch');

        width = 360 * pointsPerPixel;
        height = 75 * pointsPerPixel;
        pos = [screenSize(3)/2-width/2 screenSize(4)/2-height/2 width height];

        handles.figure = figure(...
            'Units', 'points', ...
            'BusyAction', 'queue', ...
            'WindowStyle', 'normal', ...
            'Position', pos, ...
            'Resize','off', ...
            'CreateFcn','', ...
            'NumberTitle','off', ...
            'IntegerHandle','off', ...
            'MenuBar', 'none', ...
            'Tag','TMWWaitbar',...
            'Interruptible', 'off', ...
            'DockControls', 'off', ...
            'Visible','off');

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % set figure properties as passed to the fcn
        % pay special attention to the 'cancel' request
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        visValue = 'on';
        cancelBtnCreated = 0;
        if nargin > 0,
            propList = varargin(1:2:end);
            valueList = varargin(2:2:end);

            visibleExist = strmatch('vis',lower(propList));
            if ~isempty(visibleExist)
                visValue = valueList{visibleExist};
            end

            for ii = 1:length( propList )
                try
                    if strcmpi(propList{ii}, 'createcancelbtn' ) && ~cancelBtnCreated
                        cancelBtnHeight = 23 * pointsPerPixel;
                        cancelBtnWidth = 60 * pointsPerPixel;
                        newPos = pos;
                        vertMargin = vertMargin + cancelBtnHeight;
                        newPos(4) = newPos(4)+vertMargin;
                        callbackFcn = [valueList{ii}];
                        set( handles.figure, 'Position', newPos, 'CloseRequestFcn', callbackFcn );
                        cancelButt = uicontrol('Parent',handles.figure, ...
                            'Units','points', ...
                            'Callback',callbackFcn, ...
                            'ButtonDownFcn', callbackFcn, ...
                            'Enable','on', ...
                            'Interruptible','off', ...
                            'Position', [pos(3)-cancelBtnWidth*1.4, 7,  ...
                            cancelBtnWidth, cancelBtnHeight], ...
                            'String',getString(message('MATLAB:uistring:waitbar:Cancel')), ...
                            'Tag','TMWWaitbarCancelButton'); %#ok<NASGU>
                        cancelBtnCreated = 1;
                    else
                        % simply set the prop/value pair of the figure
                        set( handles.figure, propList{ii}, valueList{ii});
                    end
                catch ex
                   warning(message('MATLAB:uistring:waitbar:WarningCouldNotSetPropertyValue',propList{ii}, num2str(valueList{ii})));
                end
            end
        end

        axNorm=[.05 .3 .9 .2];
        axPos=axNorm.*[pos(3:4),pos(3:4)] + [0 vertMargin 0 0];

        handles.axes = axes('Parent', handles.figure, ...
            'XLim',[0 100],...
            'YLim',[0 1],...
            'Box','on', ...
            'Units','Points',...
            'Position',axPos,...
            'XTickMode','manual',...
            'YTickMode','manual',...
            'XTick',[],...
            'YTick',[],...
            'XTickLabelMode','manual',...
            'XTickLabel',[],...
            'YTickLabelMode','manual',...
            'YTickLabel',[],...
            'Visible', 'off');

        handles.axesTitle = get(handles.axes,'Title');
        oldTitleUnits=get(handles.axesTitle,'Units');
        set(handles.axesTitle,...
            'Units',      'points',...
            'FontSize', 'factory',...
            'FontName', 'factory',...
            'FontWeight', 'factory',...
            'FontAngle', 'factory',...
            'FontUnits', 'factory',...
            'String',     name,...
            'Visible',    'on');

        tExtent=get(handles.axesTitle,'Extent');
        set(handles.axesTitle,'Units',oldTitleUnits);

        titleHeight=tExtent(4)+axPos(2)+axPos(4)+5;
        if titleHeight>pos(4)
            pos(4)=titleHeight;
            pos(2)=screenSize(4)/2-pos(4)/2;
            figPosDirty=true;
        else
            figPosDirty=false;
        end

        if tExtent(3)>pos(3); % If extent width is greater than the width of figure. change the figure width.
            pos(3)=min(tExtent(3)*1.10,screenSize(3));
            pos(1)=screenSize(3)/2-pos(3)/2;

            axPos([1,3])=axNorm([1,3])*pos(3);
            set(handles.axes,'Position',axPos);

            figPosDirty=true;
        end

        if figPosDirty
            set(handles.figure,'Position',pos);
            if (cancelBtnCreated)
                cancelButtonPos = get(cancelButt, 'Position');
                cancelButtonPos(1) = pos(3)-cancelButtonPos(3)*1.4;
                set(cancelButt, 'Position', cancelButtonPos);
            end
        end

        uiwaitbar(handles, 'create', x, axPos);

        set(handles.figure,'HandleVisibility','Callback','Visible', visValue);
    end

    function updateWaitbar(varargin)
        uiwaitbar(handles, 'update', x);

        if nargin>0
            % Update waitbar title:
            set(handles.axesTitle,'String',varargin{1});
        end
    end
end

function handles = extractHandles(fighandle)
% Used to get the handles structure from the given waitbar figure and
% also to check if the structure is of the right size.
handles = getappdata(fighandle,'TMWWaitbar_handles');
if isempty(handles) || (length(fieldnames(handles))~=5)
    error(message('MATLAB:waitbar:WaitbarHandlesNotFound'));
end
end
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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 3 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, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function unfreezeColors(h)
% unfreezeColors  Restore colors of a plot to original indexed color. (v2.3)
%
%   Useful if you want to apply a new colormap to plots whose
%       colors were previously frozen with freezeColors.
%
%   Usage:
%       unfreezeColors          unfreezes all objects in current axis,
%       unfreezeColors(axh)     same, but works on axis axh. axh can be vector.
%       unfreezeColors(figh)    same, but for all objects in figure figh.
%
%       Has no effect on objects on which freezeColors was not already called.
%				(Note: if colorbars were frozen using cbfreeze, use cbfreeze('off') to
%       unfreeze them. See freezeColors for information on cbfreeze.)
%
%
%   See also freezeColors, freezeColors_pub.html, cbfreeze.
%
%
%   John Iversen (iversen@nsi.edu) 3/23/05
%

%   Changes:
%   JRI 9/1/06 now restores any object with frozen CData;
%              can unfreeze an entire figure at once.
%   JRI 4/7/10 Change documentation for colorbars

% Free for all uses, but please retain the following:
%
%   Original Author:
%   John Iversen, 2005-10
%   john_iversen@post.harvard.edu

error(nargchk(0,1,nargin,'struct'))

appdatacode = 'JRI__freezeColorsData';

%default: operate on gca
if nargin < 1,
    h = gca;
end

if ~ishandle(h),
     error('JRI:unfreezeColors:invalidHandle',...
            'The argument must be a valid graphics handle to a figure or axis')
end

%if h is a figure, loop on its axes
if strcmp(get(h,'type'),'figure'),
    h = get(h,'children');
end

for h1 = h', %loop on axes

    %process all children, acting only on those with saved CData
    %   ( in appdata JRI__freezeColorsData)
    ch = findobj(h1);

    for hh = ch',

        %some object handles may be invalidated when their parent changes
        %   (e.g. restoring colors of a scattergroup unfortunately changes
        %   the handles of all its children). So, first check to make sure
        %   it's a valid handle
        if ishandle(hh)
            if isappdata(hh,appdatacode),
                ad = getappdata(hh,appdatacode);
                %get oroginal cdata
                %patches have to be handled separately (see note in freezeColors)
                if ~strcmp(get(hh,'type'),'patch'),
                    cdata = get(hh,'CData');
                else
                    cdata = get(hh,'faceVertexCData');
                    cdata = permute(cdata,[1 3 2]);
                end
                indexed = ad{1};
                scalemode = ad{2};

                %size consistency check
                if all(size(indexed) == size(cdata(:,:,1))),
                    %ok, restore indexed cdata
                    if ~strcmp(get(hh,'type'),'patch'),
                        set(hh,'CData',indexed);
                    else
                        set(hh,'faceVertexCData',indexed);
                    end
                    %restore cdatamapping, if needed
                    g = get(hh);
                    if isfield(g,'CDataMapping'),
                        set(hh,'CDataMapping',scalemode);
                    end
                    %clear appdata
                    rmappdata(hh,appdatacode)
                else
                    warning('JRI:unfreezeColors:internalCdataInconsistency',...
                        ['Could not restore indexed data: it is the wrong size. ' ...
                        'Were the axis contents changed since the call to freezeColors?'])
                end

            end %test if has our appdata
        end %test ishandle

    end %loop on children

end %loop on axes

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ma,ASAsellog,ASAcontrol] = sig2ma(sig,cand_order,last)
%SIG2MA MA model identification
%   [MA,SELLOG] = SIG2MA(SIG) estimates moving average models from the
%   data vector SIG and selects a model with optimal predictive
%   qualities. The selected model is returned in the parameter vector MA.
%   The structure SELLOG provides additional information on the selection
%   process.
%
%   SIG2MA(SIG,CAND_ORDER) selects only from candidate models whose
%   orders are entered in CAND_ORDER. CAND_ORDER must either be a row of
%   ascending orders, or a single order (in which case no true order
%   selection is performed).
%
%   Without user intervention, the mean of SIG is subtracted from the
%   data. To control the subtraction of the mean, see the help topics on
%   ASAGLOB_SUBTR_MEAN and ASAGLOB_MEAN_ADJ.
%
%   SIG2MA is an ARMASA main function.
%
%   See also: SIG2AR, SIG2ARMA, ARMASEL.

%   References: P. M. T. Broersen, Autoregressive Model Orders for
%               Durbin's MA and ARMA estimators, IEEE Transactions on
%               Signal Processing, Vol. 48, No. 8, August 2000,
%               pp. 2454-2457.

%Header
%=========================================================================

%Declaration of variables
%------------------------

%Declare and assign values to local variables
%according to the input argument pattern
switch nargin
case 1
   if isa(sig,'struct'), ASAcontrol=sig; sig=[];
   else, ASAcontrol=[];
   end
   cand_order=[];
case 2
   if isa(cand_order,'struct'), ASAcontrol=cand_order; cand_order=[];
   else, ASAcontrol=[];
   end
case 3
   if isa(last,'struct'), ASAcontrol=last;
   else, error(ASAerr(39))
   end
otherwise
   error(ASAerr(1,mfilename))
end

if isequal(nargin,1) & ~isempty(ASAcontrol)
      %ASAcontrol is the only input argument
   ASAcontrol.error_chk = 0;
   ASAcontrol.run = 0;
end

%Declare ASAglob variables
ASAglob = {'ASAglob_subtr_mean';'ASAglob_mean_adj'; ...
      'ASAglob_rc';'ASAglob_ar';'ASAglob_final_f'; ...
      'ASAglob_final_b';'ASAglob_ar_cond'};

%Assign values to ASAglob variables by screening the
%caller workspace
for ASAcounter = 1:length(ASAglob)
   ASAvar = ASAglob{ASAcounter};
   eval(['global ' ASAvar]);
   if evalin('caller',['exist(''' ASAvar ''',''var'')'])
      eval([ASAvar '=evalin(''caller'',ASAvar);']);
   else
      eval([ASAvar '=[];']);
   end
end

%ARMASA-function version information
%-----------------------------------

%This ARMASA-function is characterized by
%its current version,
ASAcontrol.is_version = [2000 12 30 20 0 0];
%and its compatability with versions down to,
ASAcontrol.comp_version = [2000 12 30 20 0 0];

%This function calls other functions of the ARMASA
%toolbox. The versions of these other functions must
%be greater than or equal to:
ASAcontrol.req_version.burg = [2000 12 30 20 0 0];
ASAcontrol.req_version.cic = [2000 12 30 20 0 0];
ASAcontrol.req_version.rc2arset = [2000 12 30 20 0 0];
ASAcontrol.req_version.cov2arset = [2000 12 30 20 0 0];
ASAcontrol.req_version.armafilter = [2000 12 12 14 0 0];
ASAcontrol.req_version.convol = [2000 12 6 12 17 20];
ASAcontrol.req_version.convolrev = [2000 12 6 12 17 20];

%Checks
%------

if ~any(strcmp(fieldnames(ASAcontrol),'error_chk')) | ASAcontrol.error_chk
   %Perform standard error checks
   %Input argument format checks
   ASAcontrol.error_chk = 1;
   if ~isnum(sig)
      error(ASAerr(11,'sig'))
   end
   if ~isavector(sig)
      error([ASAerr(14) ASAerr(15,'sig')])
   elseif size(sig,2)>1
      sig = sig(:);
      warning(ASAwarn(25,{'row';'sig';'column'},ASAcontrol))
   end
   if ~isempty(cand_order)
      if ~isnum(cand_order) | ~isintvector(cand_order) |...
            cand_order(1)<0 | ~isascending(cand_order)
         error(ASAerr(12,{'candidate';'cand_order'}))
      elseif size(cand_order,1)>1
         cand_order = cand_order';
         warning(ASAwarn(25,{'column';'cand_order';'row'},ASAcontrol))
      end
   end

   %Input argument value checks
   if ~isreal(sig)
      error(ASAerr(13))
   end
   if max(cand_order) > length(sig)-1
      error(ASAerr(21))
   end
end

if ~any(strcmp(fieldnames(ASAcontrol),'version_chk')) | ASAcontrol.version_chk
      %Perform version check
   ASAcontrol.version_chk = 1;

   %Make sure the requested version of this function
   %complies with its actual version
   ASAversionchk(ASAcontrol);

   %Make sure the requested versions of the called
   %functions comply with their actual versions
   burg(ASAcontrol);
   cic(ASAcontrol);
   rc2arset(ASAcontrol);
   cov2arset(ASAcontrol);
   armafilter(ASAcontrol);
   convol(ASAcontrol);
   convolrev(ASAcontrol);
end

if ~any(strcmp(fieldnames(ASAcontrol),'run')) | ASAcontrol.run
      %Run the computational kernel
   ASAcontrol.run = 1;
   ASAcontrol.version_chk = 0;
   ASAcontrol.error_chk = 0;
   ASAtime = clock;
   ASAdate = now;

%Main
%================================================================================================

%Initialization of variables
%---------------------------

if isempty(ASAglob_subtr_mean) | ASAglob_subtr_mean
   sig = sig-mean(sig);
   if isempty(ASAglob_mean_adj)
      ASAglob_mean_adj = 1;
   end
elseif isempty(ASAglob_mean_adj)
   ASAglob_mean_adj = 0;
end

n_obs = length(sig);
ar_stack = cell(4,1);
ar_entry = ones(1,4);
rc = [];
ma_sel = 1;

%Combined determination of the maximum candidate MA
%order and the max. candidate sliding AR order
%--------------------------------------------------

def_max_ma_order = min(fix(n_obs/5),fix(80*log10(n_obs)));
if def_max_ma_order > 400;
   def_max_ma_order = 400;
end

if isempty(cand_order)
   cand_order = 0:def_max_ma_order;
end
max_ma_order = cand_order(end);

if max_ma_order <= def_max_ma_order
   max_slid_ar_order = fix(2.5*def_max_ma_order);
else
   max_slid_ar_order = fix(2.5*max_ma_order);
   if max_slid_ar_order > n_obs-1
      max_slid_ar_order = n_obs-1;
   end
end

%Preparations for the estimation procedure
%-----------------------------------------

l_cand_order = length(cand_order);
ma = zeros(1,max_ma_order);
var = sig'*sig/n_obs;
if cand_order(1)==0
   zero_incl = 1;
   gic3 = zeros(1,l_cand_order);
   gic3(1) = log(var)+3/n_obs;
   pe_est = zeros(1,l_cand_order);
   if ASAglob_mean_adj
      pe_est(1) = var*(n_obs+1)/(n_obs-1);
   else
      pe_est(1) = var;
   end
else
   zero_incl = 0;
   cand_order = [0 cand_order];
   gic3 = zeros(1,l_cand_order+1);
   gic3(1) = inf;
   pe_est = zeros(1,l_cand_order+1);
end

if max_ma_order > 0
   %Conditioning AR orders to the previously selected AR model
   if isequal(ASAglob_ar_cond,1) & ~isempty(ASAglob_ar)
      ar = ASAglob_ar;
      sel_ar_order = length(ar)-1;
      if  2*sel_ar_order+max_ma_order < max_slid_ar_order
         max_slid_ar_order = 2*sel_ar_order+max_ma_order;
      end
   end

   %AR model estimation
   l_rc = length(ASAglob_rc);
   if l_rc>1
      rc = ASAglob_rc;
      if (l_rc < max_slid_ar_order+1)
         if isempty(ASAglob_final_f)
            ar_det = rc2arset(ASAglob_rc(1:end-1),ASAcontrol);
            ASAglob_final_f = convol(sig,ar_det,l_rc-1,n_obs,ASAcontrol);
            ASAglob_final_b = convolrev(ar_det,sig,l_rc-1,n_obs,ASAcontrol);
         end
         rc = [ASAglob_rc burg(ASAglob_final_f, ...
               ASAglob_final_b,max_slid_ar_order+1-l_rc,ASAcontrol)];
      end
   else
      rc = burg(sig,max_slid_ar_order,ASAcontrol);
   end

   %AR model order selection
   if ~isequal(ASAglob_ar_cond,1) | isempty(ASAglob_ar)
      rc(1) = 0;
      res = var*cumprod(1-rc(1:max_slid_ar_order+1).^2);
      rc(1) = 1;
      [min_value,sel_location] = min(cic(res,n_obs,ASAcontrol));
      sel_ar_order = sel_location-1;
   end

   min_ma_order = max(1,cand_order(1));

   slid_ar_order = 2*sel_ar_order+min_ma_order;
   if slid_ar_order > max_slid_ar_order
      slid_ar_order = max_slid_ar_order;
   elseif slid_ar_order < 3
      slid_ar_order = min(3,max_slid_ar_order);
   end

   pred_ar_order = min(3*sel_ar_order+min(9,1+fix(n_obs/10)),max_slid_ar_order);

   %Determine a minimum set of AR parameter vectors, as needed for the preparations
   [cand_ar_order,ar_entry] = sort([sel_ar_order pred_ar_order slid_ar_order max_slid_ar_order]);
   equal_entry = zeros(1,4);
   [dummy,redirect] = sort(ar_entry);
   equal_counter = 0;
   for i = 2:4
      if isequal(cand_ar_order(i),cand_ar_order(i-1));
         equal_counter = equal_counter+1;
         equal_entry(i) = i;
         index = find(max(0,redirect-i+equal_counter));
         redirect(index) = redirect(index)-1;
      end
   end
   cand_ar_order(find(equal_entry)) = [];
   ar_entry = redirect;
   ar_stack = rc2arset(rc,cand_ar_order,ASAcontrol);

   ar_pred = ar_stack{ar_entry(2)};
   l_ar_pred = length(ar_pred);
   l_pred_sig = fix(n_obs/2);
   pred_sig_rev = armafilter(zeros(l_pred_sig,1),ar_pred,1,...
      sig(end:-1:1),convolrev(sig,ar_pred,1,l_ar_pred,ASAcontrol),ASAcontrol);
   pred_sig = pred_sig_rev(end:-1:1);

   counter = 2;
   req_counter = 2;
   sel_index = 1;
   ar_slid = zeros(1,max_slid_ar_order+1);
   ar_slid(1:slid_ar_order+1) = ar_stack{ar_entry(3)};

   %Estimation procedure and model order selection
   %----------------------------------------------

   for order = min_ma_order:max_ma_order
      if cand_order(req_counter)==order
         ar_corr = convolrev(ar_slid(1:slid_ar_order+1),order,ASAcontrol);
         ma = cov2arset(ar_corr,ASAcontrol);

         e = armafilter(sig,ma,1,filter(1,ma,pred_sig),pred_sig,ASAcontrol);
         res = e'*e/n_obs;
         gic3_temp = log(res)+3*(order+1)/n_obs;
         gic3(req_counter) = gic3_temp;
         if gic3_temp < gic3(sel_index)
           sel_index = req_counter;
           ma_sel = ma;
         end
         if ASAglob_mean_adj
            pe_est(req_counter) = res*(n_obs+order+1)/(n_obs-order-1);
         else
            pe_est(req_counter) = res*(n_obs+order)/(n_obs-order);
         end
         req_counter = req_counter+1;
      end

      if slid_ar_order < max_slid_ar_order
         slid_ar_order = slid_ar_order+1;
         rc_temp = rc(slid_ar_order+1);
         ar_slid(2:slid_ar_order) = ar_slid(2:slid_ar_order)+rc_temp*ar_slid(slid_ar_order:-1:2);
         ar_slid(slid_ar_order+1) = rc_temp;
      end

      counter = counter+1;
   end
end

%Arranging output arguments
%--------------------------

ma = ma_sel;

if ~zero_incl
   gic3(1) = [];
   pe_est(1) = [];
   cand_order(1) =[];
end

if ~isempty(rc)
   ASAglob_rc = rc;
end

if nargout>1
   ASAsellog.funct_name = mfilename;
   ASAsellog.funct_version = ASAcontrol.is_version;
   ASAsellog.date_time = [datestr(ASAdate,8) 32 datestr(ASAdate,0)];
   ASAsellog.comp_time = etime(clock,ASAtime);
   ASAsellog.ma = ma_sel;
   ASAsellog.ar_sel = ar_stack{ar_entry(1)};
   ASAsellog.mean_adj = ASAglob_mean_adj;
   ASAsellog.cand_order = cand_order;
   ASAsellog.gic3 = gic3;
   ASAsellog.pe_est = pe_est;
end

%Footer
%=====================================================

else %Skip the computational kernel
   %Return ASAcontrol as the first output argument
   if nargout>1
      warning(ASAwarn(9,mfilename,ASAcontrol))
   end
   ma = ASAcontrol;
   ASAcontrol = [];
end

%Program history
%======================================================================
%
% Version                Programmer(s)          E-mail address
% -------                -------------          --------------
% former versions        P.M.T. Broersen        p.m.t.broersen@tudelft.nl
% [2000 12 30 20 0 0]    W. Wunderink           wwunderink01@freeler.nl