text:
function a = angle(v1,v2,varargin)
% angle between two vectors
%
% Syntax
%   omega = angle(v1,v2)
%   omega = angle(v1,v2,'antipodal')
%   omega = angle(v1,v2,N)
%
% Input
%  v1, v2 - @vector3d
%       N - @vector3d normal to plane
%
% Output
%  omega - double
%
% Options
%  antipodal  - include <VectorsAxes.html antipodal symmetry>

if nargin == 3 && isa(varargin{1},'vector3d')

  N = normalize(varargin{1});
  v1 = normalize(v1 - dot(v1,N,'noSymmetry') .* N);
  v2 = normalize(v2 - dot(v2,N,'noSymmetry') .* N);

  a  = angle(v1,v2);

  ind = dot(N,cross(v1,v2),'noSymmetry')<0;
  a(ind) = 2*pi - a(ind);

else

  a = dot(v1.normalize,v2.normalize,varargin{:});

  a = real(acos(a));

end

author: mtex-toolbox

repo: mtex

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


text:
function x = r8mat_read_bin ( input_filename, m, n )

%*****************************************************************************80
%
%% R8MAT_READ_BIN reads an R8MAT from a binary file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the input filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of data items.
%
%    Output, real X(M,N), the data.
%
  input_unit = fopen ( input_filename, 'rb' );

  x = fread ( input_unit, [ m, n ], 'double' );

  fclose ( input_unit );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function result = MTrick_enterFunc(alpha,beta,numCluster,maxIter,inputPath)
    %% ======================================================================
    %%STEP 1: load the data
    fprintf('start load the data...\n');
%     TrainData = load('MTrick/Train.data');
%     TrainData = spconvert(TrainData);
%     TrainLabel = load('MTrick/Train.label');
%     TrainLabel = TrainLabel';
%     TestData = load('MTrick/Test.data');
%     TestData = spconvert(TestData);
%     TestLabel = load('MTrick/Test.label');
%     TestLabel = TestLabel';
    [TrainData, TestData, TrainLabel] = MTrick_loadData(inputPath);
    TrainData = double(TrainData);
    TestData = double(TestData);
    TrainLabel = double(TrainLabel);
%
%     TrainData = sparse(TrainData);
%     TestData = sparse(TestData);
%
%     mex -largeArrayDims mex_Pw_d.c
%     mex -largeArrayDims mex_EMstep.c
%     mex -largeArrayDims mex_logL.c
    TrainData = sparse(TrainData);
    TestData = sparse(TestData);

    result = MTrick1(TrainData,TrainLabel,TestData,alpha,beta,numCluster,maxIter);
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function [m,im,iu] = unique(m,varargin)
% disjoint list of Miller indices%
%
% Syntax
%   u = unique(m) % find disjoined elements of the vector v
%   u = unique(m,'tolerance',0.01) % use tolerance 0.01
%   [u,im,iu] = unique(m,varargin)]
%
% Input
%  m   - @Miller
%  tol - double (default 1e-7)
%
% Output
%  u - @Miller
%  im - index such that u = m(im)
%  iu - index such that m = u(iu)
%
% Flags
%  stable     - prevent sorting
%  noSymmetry - ignore symmetry
%
% See also
% unique
%

if check_option(varargin,'noSymmetry')

  [~,im,iu] = unique@vector3d(m,varargin{:});

else
  v = vector3d(symmetrise(m,varargin{:}));

  [~,~,iu] = unique(v,varargin{:});

  [~,im,iu] = unique(min(reshape(iu,size(v)),[],1));

end

m = m.subSet(im);


text:
function parser = m2tInputParser()
%MATLAB2TIKZINPUTPARSER   Input parsing for matlab2tikz..
%   This implementation exists because Octave is lacking one.

%   Copyright (c) 2008--2014 Nico Schlömer
%   All rights reserved.
%
%   Redistribution and use in source and binary forms, with or without
%   modification, are permitted provided that the following conditions are
%   met:
%
%       * Redistributions of source code must retain the above copyright
%         notice, this list of conditions and the following disclaimer.
%       * Redistributions in binary form must reproduce the above copyright
%         notice, this list of conditions and the following disclaimer in
%         the documentation and/or other materials provided with the distribution
%
%   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%   AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%   IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%   ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%   LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%   CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%   SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%   INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%   CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%   ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%   POSSIBILITY OF SUCH DAMAGE.
% =========================================================================
  % Initialize the structure.
  parser = struct ();
  % Public Properties
  parser.Results = {};
  % Enabel/disable parameters case sensitivity.
  parser.CaseSensitive = false;
  % Keep parameters not defined by the constructor.
  parser.KeepUnmatched = false;
  % Enable/disable warning for parameters not defined by the constructor.
  parser.WarnUnmatched = true;
  % Enable/disable passing arguments in a structure.
  parser.StructExpand = true;
  % Names of parameters defined in input parser constructor.
  parser.Parameters = {};
  % Names of parameters not defined in the constructor.
  parser.Unmatched = struct ();
  % Names of parameters using default values.
  parser.UsingDefaults = {};
  % Names of deprecated parameters and their alternatives
  parser.DeprecatedParameters = struct();

  % Handles for functions that act on the object.
  parser.addRequired    = @addRequired;
  parser.addOptional    = @addOptional;
  parser.addParamValue  = @addParamValue;
  parser.deprecateParam = @deprecateParam;
  parser.parse          = @parse;

  % Initialize the parser plan
  parser.plan = {};
end
% =========================================================================
function p = parser_plan (q, arg_type, name, default, validator)
  p = q;
  plan = p.plan;
  if (isempty (plan))
    plan = struct ();
    n = 1;
  else
    n = numel (plan) + 1;
  end
  plan(n).type      = arg_type;
  plan(n).name      = name;
  plan(n).default   = default;
  plan(n).validator = validator;
  p.plan = plan;
end
% =========================================================================
function p = addRequired (p, name, validator)
  p = parser_plan (p, 'required', name, [], validator);
end
% =========================================================================
function p = addOptional (p, name, default, validator)
  p = parser_plan (p, 'optional', name, default, validator);
end
% =========================================================================
function p = addParamValue (p, name, default, validator)
  p = parser_plan (p, 'paramvalue', name, default, validator);
end
% =========================================================================
function p = deprecateParam (p, name, alternatives)
  if isempty(alternatives)
      alternatives = {};
  elseif ischar(alternatives)
      alternatives = {alternatives}; % make cellstr
  elseif ~iscellstr(alternatives)
      error('m2tInputParser:BadAlternatives',...
            'Alternatives for a deprecated parameter must be a char or cellstr');
  end
  p.DeprecatedParameters.(name) = alternatives;
end
% =========================================================================
function p = parse (p, varargin)
  plan = p.plan;
  results = p.Results;
  using_defaults = {};
  if (p.CaseSensitive)
    name_cmp = @strcmp;
  else
    name_cmp = @strcmpi;
  end
  if (p.StructExpand)
    k = find (cellfun (@isstruct, varargin));
    for m = numel(k):-1:1
      n = k(m);
      s = varargin{n};
      c = [fieldnames(s).'; struct2cell(s).'];
      c = c(:).';
      if (n > 1 && n < numel (varargin))
        varargin = horzcat (varargin(1:n-1), c, varargin(n+1:end));
      elseif (numel (varargin) == 1)
        varargin = c;
      elseif (n == 1);
        varargin = horzcat (c, varargin(n+1:end));
      else % n == numel (varargin)
        varargin = horzcat (varargin(1:n-1), c);
      end
    end
  end
  if (isempty (results))
    results = struct ();
  end
  type = {plan.type};
  n = find( strcmp( type, 'paramvalue' ) );
  m = setdiff (1:numel( plan ), n );
  plan = plan ([n,m]);
  for n = 1 : numel (plan)
    found = false;
    results.(plan(n).name) = plan(n).default;
    if (~ isempty (varargin))
      switch plan(n).type
      case 'required'
        found = true;
        if (strcmpi (varargin{1}, plan(n).name))
          varargin(1) = [];
        end
        value = varargin{1};
        varargin(1) = [];
      case 'optional'
        m = find (cellfun (@ischar, varargin));
        k = find (name_cmp (plan(n).name, varargin(m)));
        if (isempty (k) && validate_arg (plan(n).validator, varargin{1}))
          found = true;
          value = varargin{1};
          varargin(1) = [];
        elseif (~ isempty (k))
          m = m(k);
          found = true;
          value = varargin{max(m)+1};
          varargin(union(m,m+1)) = [];
        end
      case 'paramvalue'
        m = find( cellfun (@ischar, varargin) );
        k = find (name_cmp (plan(n).name, varargin(m)));
        if (~ isempty (k))
          found = true;
          m = m(k);
          value = varargin{max(m)+1};
          varargin(union(m,m+1)) = [];
        end
      otherwise
        error( sprintf ('%s:parse', mfilename), ...
               'parse (%s): Invalid argument type.', mfilename ...
             )
      end
    end
    if (found)
      if (validate_arg (plan(n).validator, value))
        results.(plan(n).name) = value;
      else
        error( sprintf ('%s:invalidinput', mfilename), ...
               '%s: Input argument ''%s'' has invalid value.\n', mfilename, plan(n).name ...
             );
      end
      p.Parameters = union (p.Parameters, {plan(n).name});
    elseif (strcmp (plan(n).type, 'required'))
      error( sprintf ('%s:missinginput', mfilename), ...
             '%s: input ''%s'' is missing.\n', mfilename, plan(n).name ...
           );
    else
      using_defaults = union (using_defaults, {plan(n).name});
    end
  end

  if ~isempty(varargin)
    % Include properties that do not match specified properties
      for n = 1:2:numel(varargin)
        if ischar(varargin{n})
          if p.KeepUnmatched
            results.(varargin{n}) = varargin{n+1};
          end
          if p.WarnUnmatched
            warning(sprintf('%s:unmatchedArgument',mfilename), ...
                    'Ignoring unknown argument "%s"', varargin{n});
          end
          p.Unmatched.(varargin{n}) = varargin{n+1};
        else
          error (sprintf ('%s:invalidinput', mfilename), ...
                 '%s: invalid input', mfilename)
        end
      end
  end

  % Store the results of the parsing
  p.Results = results;
  p.UsingDefaults = using_defaults;

  warnForDeprecatedParameters(p);
end
% =========================================================================
function result = validate_arg (validator, arg)
  try
    result = validator (arg);
  catch %#ok
    result = false;
  end
end
% =========================================================================
function warnForDeprecatedParameters(p)
  usedDeprecatedParameters = intersect(p.Parameters, fieldnames(p.DeprecatedParameters));

  for iParam = 1:numel(usedDeprecatedParameters)
      oldParameter = usedDeprecatedParameters{iParam};
      alternatives = p.DeprecatedParameters.(oldParameter);

      switch numel(alternatives)
          case 0
              replacements = '';
          case 1
              replacements = ['''' alternatives{1} ''''];
          otherwise
              replacements = deblank(sprintf('''%s'' and ',alternatives{:}));
              replacements = regexprep(replacements,' and$','');
      end
      if ~isempty(replacements)
          replacements = sprintf('From now on, please use %s to control the output.\n',replacements);
      end

      message = ['\n===============================================================================\n', ...
                  'You are using the deprecated parameter ''%s''.\n', ...
                  '%s', ...
                  '==============================================================================='];
      warning('matlab2tikz:deprecatedParameter', ...
                message, oldParameter, replacements);

  end
end

author: thomas-koehler

repo: SupER


text:
function [fx] = f_FitzHughNagumo_calcium(Xt,Theta,ut,inF)

% FitzHugh-Nagumo membrane potential evolution function (+convolution)
% function [fx] = FitzHughNagumo(Xt,Theta,ut,inF)
% IN:
%   - Xt: system's states, ie:
%       Xt(1): calcium imaging prediction
%       Xt(2): membrane depolarization (mV)
%       Xt(3): proxy for ion channel opening probabilities
%   - Theta: evolution parameters (see bellow)
%   - ut: input current
%   - inF: [optional]
% OUT:
%   - fx: the evolution function evaluated at Xt

deltat = inF.delta_t;

try
    a = inF.a.^-1;
catch
    a = 1;
end
a = a.*exp(Theta(1));
K1 = (1/3)*exp(Theta(2));
K2 = 0.08*exp(Theta(3));
K3 = 0.7*exp(Theta(4));
K4 = 0.8*exp(Theta(5));
tmp = [ -a*Xt(1)+ exp(Xt(2))
        Xt(2) - K1*Xt(2).^3 - Xt(3) + ut
        K2*(Xt(2) - K3 - K4*Xt(3))];

fx = Xt + deltat.*tmp;

% A = [ -a  exp(Xt(2))
%        0  -b          ];
% dF_dX = eye(2) + deltat.*A';
%
% dF_dTheta = [ -deltat.*Xt(1)*a  0
%               0                 -deltat.*Xt(2)*b];

author: MBB-team

repo: VBA-toolbox


text:
function [dmodedr, dmodeds] = GradSimplex2DP(a,b,id,jd)

% function [dmodedr, dmodeds] = GradSimplex2DP(a,b,id,jd)
% Purpose: Return the derivatives of the modal basis (id,jd) on the 2D simplex at (a,b).

fa = JacobiP(a, 0, 0, id);     dfa = GradJacobiP(a, 0, 0, id);
gb = JacobiP(b, 2*id+1,0, jd); dgb = GradJacobiP(b, 2*id+1,0, jd);

% r-derivative
% d/dr = da/dr d/da + db/dr d/db = (2/(1-s)) d/da = (2/(1-b)) d/da
dmodedr = dfa.*gb;
if(id>0)
  dmodedr = dmodedr.*((0.5*(1-b)).^(id-1));
end

% s-derivative
% d/ds = ((1+a)/2)/((1-b)/2) d/da + d/db
dmodeds = dfa.*(gb.*(0.5*(1+a)));
if(id>0)
 dmodeds = dmodeds.*((0.5*(1-b)).^(id-1));
end

tmp = dgb.*((0.5*(1-b)).^id);
if(id>0)
  tmp = tmp-0.5*id*gb.*((0.5*(1-b)).^(id-1));
end
dmodeds = dmodeds+fa.*tmp;

% Normalize
dmodedr = 2^(id+0.5)*dmodedr; dmodeds = 2^(id+0.5)*dmodeds;
return;

repo: Aero-matlab


text:
function PCs = pca(this, cutoff, pcDimension)
% Computes principal component analysis (PCA) of image time series
%
%   Y = MrImage()
%   PCs = pca(this, cutoff, pcDimension)
%
% This is a method of class MrImage.
%
% Spatial PCA (default):
% The principal component analysis report representative spatial
% distributions ("principal components") of fluctuation patterns that explain most variance in the
% data (over all time points).
% The corresponding temporal evolution of weights (or scores or
% projections) describes how each of these spatial patterns fluctuates over
% time.
%
% Temporal PCA (?)
% The principal component analysis reports the representative time series
% ("principal components") that explain most of the variance in the data
% (pooled over all voxels).
% The corresponding spatial maps of weights (or scores or projections) give
% the relative contribution of this particular time series to the variance
% in each voxel
%
% IN
%   cutoff      value determining number of principal components extracted
%               0 < cutoff < 1      interpreted as relative amount of
%                                   variance explained;
%                                   nPCs will be determined as the number
%                                   of components that explain at least
%                                   cutoff*100 % of the variance in the
%                                   data
%               cutoff = 1,2,3...   interpreted as number of PCs extracted
%   pcDimension 'spatial' (default) or 'temporal'
%               determines which dimension i.e. spatial image or time
%               (volumes) will be principal component, and consequently,
%               which other one will be the projection dimension
%               Technically, the non-PC dimension is the one considered to
%               "generate the variance", in that the covariance matrix
%               entries would relate PC vector components, co-varying over
%               the non-PC dimension,
%               e.g.    spatial PCA: How do two voxels co-vary over time
%                       temporal PCA: How do two time-points co-vary over
%                       space
%               Typically, the 1st projection of the temporal PCA looks
%               like the mean, and for the other components, it usually
%               holds:
%               PC_spatial(n) approximately equal to Proj_temporal(n+1)
% OUT
%   PCs         cell(nPCs,1) of MrImages
%               Each element is a 4D image, the nth volume is computed as
%                   PCs{k}_n = PC{k}*Projection_n,
%               i.e. for
%                   spatial PCA: principal component * nth time point pf
%                   projection
%                   temporal PCA: nth element of principal component *
%                   projection vector of all voxels
%
% EXAMPLE
%   pca
%
%   See also MrImage

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

% if number of components specified explicity, no additional variance
% threshold is needed

if nargin < 2
    cutoff = 1;
end

isSpatialPca = nargin < 3 || strcmpi(pcDimension, 'spatial');

hasVarianceThreshold = cutoff<1;

if hasVarianceThreshold
    nComponents = 1;
else
    nComponents = cutoff;
    cutoff = 0;
end

% created reshaped data matrix for PCA

% data matrix X: [nVoxels, nVolumes] for 4D, [nVoxel2D, nSlices] for 3D...
applicationDimension = find(this.geometry.nVoxels>1, 1, 'last');
X = reshape(this.data, [], this.geometry.nVoxels(applicationDimension));

% [nVolumes, nVoxels] for spatial PCA...
if isSpatialPca
    X = X';
end

% remove invalid data for PCA
X(isinf(X)) = 0;
X(isnan(X)) = 0;

% iteratively increase number of components, if variance threshold given,
% otherwise compute once with number of components specified
doPca = 1;
while doPca

    % Explanation for temporal PCA
    % COEFF = [nVolumes, nPCs]  principal components (PCs) ordered by variance
    %                           explained
    % SCORE = [nVoxel, nPCs]    loads of each component in each voxel, i.e.
    %                           specific contribution of each component in
    %                           a voxel's variance
    % LATENT = [nPCs, 1]        eigenvalues of data covariance matrix,
    %                           stating how much variance was explained by
    %                           each PC overall
    % TSQUARED = [nVoxels,1]    Hotelling's T-Squared test whether PC
    %                           explained significant variance in a voxel
    % EXPLAINED = [nPCs, 1]     relative amount of variance explained (in
    %                           percent) by each component
    % MU = [1, nVolumes]        mean of all time series
    [COEFF, SCORE, LATENT, TSQUARED, EXPLAINED, MU] = ...
        pca(X, 'NumComponents', nComponents);

    explainedVariance   = sum(EXPLAINED(1:nComponents))/100;
    doPca               = hasVarianceThreshold && ...
        (explainedVariance < cutoff);

    if doPca
        nComponents = nComponents + 1;

        % somehow, pca also gives out EXPLAINED variance also for more than one
        % component, jump to 1st achieved cutoff-threshold
        nComponentsTemp = find(cumsum(EXPLAINED) >= cutoff*100, 1, 'first');
        if ~isempty(nComponentsTemp)
            nComponents = nComponentsTemp;
        end
    end
end

% create 4D PCs i.e. 3D Pcs which are co-varied along the volume
% dimension with the computed projections

PCs = cell(nComponents,1);
for c = 1:nComponents
    PCs{c} = this.copyobj;
    PCs{c}.rois = {};
    PCs{c}.name = sprintf(['PC %d (explained variance: %5.2f %%, ', ...
        'cumulative %5.2f %%)) - %s'], ...
        c, EXPLAINED(c), sum(EXPLAINED(1:c)), this.name);

    % multiply the PC with each weight entry to generate whole time series
    Y = kron(SCORE(:,c)', COEFF(:,c));

    if ~isSpatialPca
        Y = Y';
    end

    PCs{c}.data = reshape(Y, this.geometry.nVoxels);
end

author: translationalneuromodeling

repo: tapas


text:
function spec = cfg_findspec(cellspec)

% function spec = cfg_findspec(cellspec)
% Create a find specification. cellspec should contain a cell array of
% cells, each of them containing name/value pairs that will be combined
% into a struct suitable for e.g. @cfg_item/match and @cfg_item/list.
% These methods will be used to e.g. select items in a configuration tree or
% to match dependencies and input items.
%
% Name/value pairs within a cell will be OR concatenated, while cells
% will be AND concatenated.
%
% A cellspec
%  {{'field1','val1','field2','val2'},{'field3','val3'}}
%
% matches an item if
%  (item.field1==val1 || item.field2==val2) && item.field3==val3
%
% If the field name is 'class', an item matches, if its class name is equal to
% spec.value.
%
% For class specific matching rules, see the help for the
% resp. @cfg_.../match method.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_findspec.m 3944 2010-06-23 08:53:40Z volkmar $

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

if nargin == 0 || isempty(cellspec)
    spec = {};
    return;
end;
spec = cell(size(cellspec));
for k = 1:numel(cellspec)
    spec{k} = struct('name',{}, 'value',{});
    for l = 1:2:numel(cellspec{k})
        spec{k}(end+1).name = cellspec{k}{l};
        spec{k}(end).value = cellspec{k}{l+1};
    end;
end;

author: spm


text:
filename='Cantileverbeam_Hexahedra_Linear_Structured_fixedcorners';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'holes';
cost = {'compliance'};
weights = [1];
constraint = {'volumeConstraint'};
optimizer = 'HAMILTON-JACOBI';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';
constraint_case = 'INEQUALITY';

nsteps = 10;
Vfrac_final = 0.15;
Perimeter_target=3.5;
optimality_final =1e-4;
constr_final =1e-4;

BCscale_factor = 0.3;
HJiter0 = 1;
e2 = 100;
N_holes = [7 3 3];
R_holes = 0.9;
phase_holes = [pi/2 0 0];

Vfrac_initial = 0.3;
optimality_initial = 5e-2;
constr_initial = 5e-2;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;

% maxiter = 1;
TOL.nu_plus = 0.3;
TOL.nu_minus = 0.3;

plotting = 1;
printing = 0;
monitoring = 1;
monitoring_interval = 1;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [A, b, x, ProbInfo] = PRblurdefocus(varargin)
% PRblurdefocus Image deblurring problem with a defocus point spread function
%
% [A, b, x, ProbInfo] = PRblurdefocus
% [A, b, x, ProbInfo] = PRblurdefocus(n)
% [A, b, x, ProbInfo] = PRblurdefocus(n, options)
% [A, b, x, ProbInfo] = PRblurdefocus(options)
%
% Input:
%  n      -  size of the image. Can be either a scalar n (in this case
%            the size is n x n) or a vector [nrow, ncol] (in this case the
%            size is (nrow x ncol).
%            Default: n = 256.
%  options - Structure containing the following optional fields:
%    trueImage  : test image of size n, of type numeric, 2-D only,
%                 or character string indicating
%                 'pattern1'  : geometrical image
%                 'pattern2'  : geometrical image
%                 'ppower'    : random image with patterns of nonzero pixels
%                 'smooth'    : very smooth image
%                 'dot2'      : two small Gaussian shaped dots, e.g., a
%                               binary star
%                 'dotk'      : n/2 small Gaussian shaped dots, e.g., stars
%                               (placement is random, reset using rng(0))
%                 'satellite' : satellite test image
%                 'hst'       : image of the Hubble space telescope
%                 Default: 'hst'.
%                 This image is then stored in the output vector x.
%    BlurLevel  : If choosing one of the built-in PSFs, this sets the
%                 severity of the blur to one of the following:
%                 'mild'
%                 'medium'
%                 'severe'
%                 Default is 'medium'
%    BC         : Specify boundary condition:
%                 'zero'
%                 'periodic'
%                 'reflective' (or 'neumann' or 'reflexive')
%                 Default: 'reflective'
%                 Note that in this case an extended (or padded) test image
%                 is blurred using 'zero' boundary conditions, and then the
%                 central subimage of size n is extracted from the exact and
%                 the blurred image. No inverse crime is committed,
%                 i.e., A*x ~= b.
%    CommitCrime: To get an exact system Ax = b (i.e., commit the inverse
%                 crime), set this to:
%                 'on'
%                 Default is 'off' (do not commit the inverse crime).
%
% Output:
%  A        - blurring matrix (psfMatrix class)
%  b        - blurred vector (i.e., blurred image with stacked columns)
%  x        - image vector, i.e., exact (unknown) image with stacked columns
%  ProbInfo - structure whose fields contain information about problem:
%               problemType : kind of test problem generated
%                              (in this case: 'deblurring')
%                xType       : solution type (in this case 'image2D')
%                bType       : data type (in this case 'image2D')
%                xSize       : size of image x
%                bSize       : size of image b
%                psf         : point spread function
%
% See also: PRblur, PRblurgauss, PRblurmotion, PRblurrotation,
% PRblurshake, PRblurspeckle, PRdiffusion, PRinvinterp2, PRnmr,
% PRseismic, PRspherical, PRtomo, PRnoise, PRshowb, PRshowx, fspecial

% Silvia Gazzola, University of Bath
% Per Christian Hansen, Technical University of Denmark
% James G. Nagy, Emory University
% April, 2018.

% This file is part of the IR Tools package and is distributed under the
% 3-Clause BSD License. A separate license file should be provided as part
% of the package.

% Set default values for options.
defaultopt = struct('trueImage', 'hst', 'BlurLevel', 'medium', ...
    'BC', 'reflective', 'CommitCrime', 'off');

% If input is 'defaults,' return the default options in X
if nargin == 1 && nargout <= 1 && strcmp(varargin,'defaults')
    A = defaultopt;
    return;
end

% Check for acceptable number of optional input arguments
switch length(varargin)
    case 0
        n = []; options = [];
    case 1
        if isa(varargin{1}, 'double')
            n = varargin{1}; options = [];
        else
            n = []; options = varargin{1};
        end
    case 2
        if isa(varargin{1}, 'double')
            n = varargin{1}; options = varargin{2};
        else
            n = varargin{2}; options = varargin{1};
        end
    otherwise
        error('Too many input parameters')
end

if isempty(options)
    options = defaultopt;
end

options = PRset(defaultopt, options);
options = PRset(options, 'PSF', 'defocus');
[A, b, x, ProbInfo] = PRblur(n, options);

repo: IRtools


text:
function c = min_errorbar_scale(stderrs,significance)
% c = min_errorbar_scale(stderrs,significance)
% returns the minimum scale factor c such that any pair of non-overlapping
% error bars represents statistical significance exceeding the given level.
% stderrs is a vector.
% significance is a scalar.
%
% See "Judging significance from error bars" by Tom Minka (2002)
% http://research.microsoft.com/~minka/papers/minka-errorbars.pdf
%
% Examples:
%   min_errorbar_scale(ones(10,1))  % returns 1.16
%   min_errorbar_scale(0:10)        % returns 1.64

% Written by Tom Minka

% Algorithm:
% We want c*(stderr(i)+stderr(j)) >= z*sqrt(stderr(i)^2+stderr(j)^2)
% for all pairs (i,j).  Therefore
% c = max_(i,j)  z*sqrt(stderr(i)^2+stderr(j)^2)/(stderr(i) + stderr(j));

if nargin < 2
  significance = 0.95;
end

z = erfinv(2*significance - 1)*sqrt(2);

n = length(stderrs);
ratio = zeros(n,1);
for i = 1:length(stderrs)
  exact = sqrt(stderrs(i).^2 + stderrs.^2);
  approx = stderrs(i) + stderrs;
  ratio(i) = max(exact./approx);
end
c = z*max(ratio);

author: andrewssobral

repo: lrslibrary


text:
function [ objectIndex ] = GetScenario_random(varargin)

fprintf('[SCENARIO]\tGetting the random distribution of objects scenario.\n');

% DEFAULT AGENT/OBSTACLE CONFIGURATION IN THIS EXAMPLE
defaultConfig = struct(...
    'file','scenario.mat',...
    'objects',[],...
    'velocity',18,...
    'positionGain',10,...
    'velocityGain',1,...
    'poseGain',pi,...
    'plot',false,....
    'noiseFactor',0);

% Instanciate the scenario builder
SBinstance = scenarioBuilder();
% Parse user inputs
[inputConfig] = SBinstance.configurationParser(defaultConfig,varargin);

% Get the random object configuration
objectConfig = SBinstance.random(...
    'objects',numel(inputConfig.objects),...
    'positionGain',inputConfig.positionGain,...
    'velocityGain',inputConfig.velocityGain,...
    'poseGain',inputConfig.poseGain);

for index = 1:numel(inputConfig.objects)
    % Assign object
    objectIndex{index} = inputConfig.objects{index};
    % APPLY GLOBAL STATE VARIABLES
    objectIndex{index}.SetGLOBAL('position',objectConfig.positions(:,index) + inputConfig.noiseFactor*randn(3,1));
    objectIndex{index}.SetGLOBAL('velocity',objectConfig.velocities(:,index) + inputConfig.noiseFactor*randn(3,1));
    objectIndex{index}.SetGLOBAL('quaternion',objectConfig.quaternions(:,index));
end
% PLOT THE SCENE
if inputConfig.plot
    SBinstance.plotObjectIndex(objectIndex);
end
% CLEAR THE REMAINING VARIABLES
clearvars -except objectIndex
end

repo: OpenMAS


text:
function [post, a] = gtmpost(net, data)
%GTMPOST Latent space responsibility for data in a GTM.
%
%	Description
%	 POST = GTMPOST(NET, DATA) takes a GTM structure NET, and computes
%	the  responsibility at each latent space sample point NET.X for each
%	data point in DATA.
%
%	[POST, A] = GTMPOST(NET, DATA) also returns the activations A of the
%	GMM NET.GMMNET as computed by GMMPOST.
%
%	See also
%	GTM, GTMEM, GTMLMEAN, GMLMODE, GMMPROB
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Check for consistency
errstring = consist(net, 'gtm', data);
if ~isempty(errstring)
  error(errstring);
end

net.gmmnet.centres = rbffwd(net.rbfnet, net.X);

[post, a] = gmmpost(net.gmmnet, data);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function U = rhs_computation(M,sphere_translational_velocity)
%+========================================================================+
%|                                                                        |
%|           This function uses the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT :                                                            |
%| PROPERTY  :                                                            |
%| LICENCE   :                                                            |
%| CONTACT   :                                                            |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : rhs_computation.m                             |
%|    #    |   VERSION    : 0.10                                          |
%|   _#_   |   AUTHOR(S)  : Luca Berti                                    |
%|  ( # )  |   CREATION   : 25.12.2018                                    |
%|  / 0 \  |   LAST MODIF :                                               |
%| ( === ) |   SYNOPSIS   :                                               |
%|  `---'  |                                                              |
%+========================================================================+
% COMPUTE THE RIGHT-HAND SIDE OF THE STOKES SINGLE LAYER PROBLEM
% ONLY P1 CASE TESTED FOR THE MOMENT

N_unknown = size(M,1)/3;

U = M*[sphere_translational_velocity(1)*ones(N_unknown,1);
     sphere_translational_velocity(2)*ones(N_unknown,1);
     sphere_translational_velocity(3)*ones(N_unknown,1)];

end


text:
function [C,top, bottom, left, right] = imtrim(im,varargin)
  % IMTRIM auto-crop an image like Photoshop's Edit>Trim feature, as of yet only
  % grayscale imagse are supported
  %
  % C = imtrim(im)
  % [C,top,bottom,left,right] = imtrim(im,'ParamName',ParamValue, ...)
  %
  % Inputs:
  %   im  h by w by nc image
  %   Optional:
  %     'Location'  followed by {'NorthWest'} or 'SouthEast' specifying where
  %       to find trim color
  %     'Dims' followed by 2 long logic whether to use each dimension: [1 0]
  %       means only crop in h and not in w
  %     'Threshold' followed by distance threshold {0.1}
  %     'Map'  followed by colormap {[]}
  % Outputs:
  %   C  sh by wh cropped image
  %   top  index of top pixel in im
  %   bottom  index of bottom pixel in im
  %   left  index of left pixel in im
  %   righ  index of right pixel in im
  %
  %   Copyright Alec Jacobson, 2010
  %
  % Example:
  %  % Load in a one-way animation
  %  [X,M] = imread('bunny.gif');
  %  % Trim to animations extent
  %  Y = imtrim(X,'Map',M,'Threshold',0);
  %  writeGIF('output.gif',Y,M,0);
  %

  location = 'NorthWest';
  threshold = [];
  dims = [1 1];
  map = [];

  % parse optional input parameters
  v = 1;
  while v < numel(varargin)
    switch varargin{v}
    case 'Map'
      assert(v+1<=numel(varargin));
      v = v+1;
      map = varargin{v};
    case 'Location'
      assert(v+1<=numel(varargin));
      v = v+1;
      location = varargin{v};
    case 'Threshold'
      assert(v+1<=numel(varargin));
      v = v+1;
      threshold = varargin{v};
    case 'Dims'
      assert(v+1<=numel(varargin));
      v = v+1;
      dims = varargin{v};
    otherwise
      error('Unsupported parameter: %s',varargin{v});
    end
    v = v+1;
  end

  if isempty(threshold)
    % default threshold parameter, works equivalently with Photoshop's
    % hardcoded parameter
    threshold = 0.1;
  end

  im_in = im;

  if ~isempty(map)
    nf = size(im,4);
    im = cell2mat(permute(arrayfun(@(C) ind2rgb(im(:,:,:,C),map(:,:,min(end,C))),1:size(im,4),'UniformOutput',false),[1 4 3 2]));
  end

  if isempty(map) && ~isfloat(im)
    warning('converting input to double');
    im = im2double(im);
  end

  % gather corner value to which the image is compared
  if(strcmp(location,'NorthWest'))
    corner_value = im(1,1,:,1);
  elseif(strcmp(location,'SouthEast'))
    corner_value = im(1,1,:,1);
  else
    error([location ' is not a valid location']);
  end

  % get difference of image with corner value
  %difference = abs(im - corner_value)>0.1;
  % should work for any number of channels
  %difference = sqrt(sum((im - corner_value).^2,3)) > ...
  %  sqrt(threshold^2*size(im,3));
  difference = sqrt(sum(bsxfun(@minus,im,corner_value).^2,3)) > ...
    sqrt(threshold^2*size(im,3));
  difference = max(difference,[],4);

  if dims(2)
    [left_i,left] = ind2sub(size(difference),find(difference,1));
    [right_i,right] = ind2sub(size(difference),find(difference,1,'last'));
  else
    left = 1;
    right = size(im,2);
  end
  if dims(1)
    [top_j,top] = ind2sub(size(difference'),find(difference',1));
    [bottom_j,bottom] = ind2sub(size(difference'),find(difference',1,'last'));
  else
    top = 1;
    bottom = size(im,1);
  end
  C = im_in(top:bottom,left:right,:,:);
end

author: alecjacobson

repo: gptoolbox


text:
function doseCompare(command,varargin)
% doseCompare
% Used to show dose associated with scan in a compare mode.
%
% Written DK
%
% Usage: doseCompare(command,varargin)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC stateS
indexS = planC{end};
hCSVA = stateS.handle.CERRSliceViewerAxis;
switch upper(command)
    case 'INIT'
        if stateS.layout == 7
            return
        end
        scanNum = getAxisInfo(hCSVA,'scanSets');

        if isempty(scanNum)
            scanNum = stateS.scanSet;
        end

        % doseNum = length(planC{indexS.dose});%Get number of dose
        doseNum = getScanAssociatedDose(scanNum,'all');

        if length(doseNum) < 2
            warndlg('Dose Comparison tool requires 2 or more doses associated to currnet scan');
            return
        elseif length(doseNum) > 4
            newAxis = 3;
        else
            newAxis = length(doseNum)-1;
        end
        stateS.doseCompare.newAxis = newAxis;

        if stateS.layout == 6
            scanCompare('exit');
        end

        if length( stateS.handle.CERRAxis)>4
            delete(stateS.handle.CERRAxis(4:end));
            stateS.handle.CERRAxisLabel1(5:end)=[];
            stateS.handle.CERRAxisLabel2(5:end)=[];
            stateS.handle.CERRAxis(5:end)=[];
        end

        %Calculating positions for each viewport
        leftMarginWidth = 195; bottomMarginHeight = 70;
        pos = get(stateS.handle.CERRSliceViewer, 'position');
        figureWidth = pos(3); figureHeight = pos(4);
        wid = (figureWidth-leftMarginWidth-70-10)/5;
        hig = (figureHeight-bottomMarginHeight-20)/2;

        stateS.oldDoseSet = stateS.doseSet;
        stateS.doseSet = doseNum(1);

        %change the layout after checking the number of dose in the plan
        stateS.layout = 7;
        stateS.Oldlayout = 7;
        doseSetsLast = getAxisInfo(hCSVA,'doseSets');
        scanSetsLast = getAxisInfo(hCSVA,'scanSets');
        setAxisInfo(hCSVA,'doseSets',doseNum(1),'doseSelectMode', 'manual','doseSetsLast',doseSetsLast,'scanSetsLast',scanSetsLast);
        %delete any new duplicate axis.
        if length(stateS.handle.CERRAxis) > 4
            delete(stateS.handle.CERRAxis(5:end));
            stateS.handle.CERRAxis = stateS.handle.CERRAxis(1:4);
        end

        for i = 1:newAxis % create linked axis to the transverse axis
            stateS.handle.CERRAxis(end+1) = axes('parent',stateS.handle.CERRSliceViewer, 'units', 'pixels', 'position', [1 1 1 1],...
                'color', [0 0 0], 'xTickLabel', [], 'yTickLabel', [], 'xTick', [], 'yTick', [], 'buttondownfcn',...
                'sliceCallBack(''axisClicked'')', 'nextplot', 'add', 'yDir', 'reverse', 'linewidth', 2,'visible','off','Tag','doseCompareAxes');

            stateS.handle.CERRAxisLabel1(end+1) = text('parent', stateS.handle.CERRAxis(end), 'string', '', 'position', [.02 .98 0],...
                'color', [1 0 0], 'units', 'normalized', 'visible', 'off', 'horizontalAlignment', 'left', 'verticalAlignment', 'top');

            stateS.handle.CERRAxisLabel2(end+1) = text('parent', stateS.handle.CERRAxis(end), 'string', '', 'position', [.90 .98 0],...
                'color', [1 0 0], 'units', 'normalized', 'visible', 'off', 'horizontalAlignment', 'left', 'verticalAlignment', 'top');

            AI(i).miscHandles = [stateS.handle.CERRAxisLabel1(end) stateS.handle.CERRAxisLabel2(end)];
        end

        axisInfo = get(hCSVA, 'userdata');
        axisInfo.scanObj(1:end) = [];
        axisInfo.doseObj(1:end) = [];
        axisInfo.structureGroup(1:end) = [];
        axisInfo.miscHandles = [];
        axisInfo.coord       = {'Linked', hCSVA};
        axisInfo.view        = {'Linked', hCSVA};
        axisInfo.xRange      = {'Linked', hCSVA};
        axisInfo.yRange      = {'Linked', hCSVA};

        for i = 1:newAxis
            axisInfo.miscHandles = AI(i).miscHandles;
            set(stateS.handle.CERRAxis(4+i), 'userdata', axisInfo);
            setAxisInfo(stateS.handle.CERRAxis(4+i),'doseSets',doseNum(i+1),'doseSelectMode', 'manual');
            CERRAxisMenu(stateS.handle.CERRAxis(4+i));
            set(stateS.handle.CERRAxis(4+i),'visible','on');
        end

        if stateS.MLVersion >= 8.4
            set(stateS.handle.CERRAxis,'ClippingStyle','rectangle')
        end

        stateS.handle.doseColorbar.Compare = axes('units', 'pixels', 'position', [leftMarginWidth+60+wid*4+20 bottomMarginHeight+30 50 hig-40],...
            'xTickLabel', [], 'yTickLabel', [], 'xTick', [], 'yTick', [], 'Tag', 'ColorbarCompare', 'visible', 'off');
        stateS.doseSetChanged = 1;
        CERRRefresh
        sliceCallBack('RESIZE');

    case 'COMPAREMODE'
        hAxis = gca;
        stateS.doseSetChanged = 1;
        refDose= getAxisInfo(hCSVA,'doseSets');
        tarDose= getAxisInfo(hAxis,'doseSets');

        doseNumOld = length(planC{indexS.dose});

        switch upper(varargin{1})
            case 'DEFAULT'
                doseNum = getAssociatedDose(planC{indexS.dose}(getAxisInfo(hAxis,'doseSetsLast')).doseUID);
                setappdata(hAxis,'compareMode',[]);
            otherwise
                setappdata(hAxis,'compareMode',upper(varargin{1}));
                [doseNum flag] = chkDose(refDose,tarDose);
                if ~flag & isempty(doseNum)%if not calculated
                    DoseAdditionSubtraction(planC{indexS.dose}(refDose),planC{indexS.dose}(tarDose),'Subtract',0);
                    doseNum = length(planC{indexS.dose});
                    if doseNumOld == doseNum
                        return
                    end
                end
        end
        setAxisInfo(hAxis,'doseSets',doseNum);
        CERRRefresh

    case 'EXIT'
        stateS.Oldlayout = [];
        stateS.doseSet = stateS.oldDoseSet;
        stateS.doseSetChanged = 1;
        delete(stateS.handle.doseColorbar.Compare)

        hDel = (findobj('Tag','doseCompareAxes'));
        for i = 1:length(hDel)
            hDel_ind = find(stateS.handle.CERRAxis == hDel(i));
            delete(stateS.handle.CERRAxis(hDel_ind));
            stateS.handle.CERRAxisLabel1(hDel_ind)=[];
            stateS.handle.CERRAxisLabel2(hDel_ind)=[];
            stateS.handle.CERRAxis(hDel_ind)=[];
        end
        for i = 1:length(stateS.handle.CERRAxis)
            setAxisInfo(stateS.handle.CERRAxis(i),'doseSets',getAxisInfo(stateS.handle.CERRAxis(i),'doseSetsLast'),...
                'scanSets',getAxisInfo(stateS.handle.CERRAxis(i),'scanSetsLast'));
            setappdata(stateS.handle.CERRAxis(i),'compareMode',[]);
        end
        sliceCallBack('resize');
        CERRRefresh
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ back, sigma ] = image_background( im, make_plot )
% IMAGE_BRACKGROUND Automatically find the pixel-value of the background,
% for images that are predominantly black
%
% [ back, sigma ] = image_background( im )
% [ back, sigma ] = image_background( im, make_plot )
%
% Automatically find the pixel-value of the background, for images that are
% predominantly black.
%
% In that case, the histogram of the image has a peak in the low values,
% of Gaussian shape which correspond to the background pixels.
% image_background returns the mean and the sigma of this peak
%
% if 'make_plot' is set and true, an histogram and detected values is displayed
%
% TODO: automatically detect when the method fails

% F. Nedelec, 2005 - April 2008
%
% Ramon Casero <rcasero@gmail.com>: Minor edits

%compatibility with tiffread
if ( isfield( im, 'data') )
    filename = im.filename;
    im = im(1).data;
else
    filename = inputname(1);
end

% add a plot:
if nargin < 2
    make_plot = 0;
end

h = image_histogram( im );
%remove under-exposed pixels, where there might be a peak:
h(1) = 0;

%% calculate a Gaussian filter:

filt_sig = 20;
filt = exp( -( (-2*filt_sig:2*filt_sig) / filt_sig ) .^ 2 );
filt = filt ./ sum( filt );

% convolve by the filter to get a smooth profile:

hs = double( conv( h, filt ) );

%crop to the same size as h
hs = hs( 2*filt_sig+(1:numel(h)) );

%% find the maximum in the smoothed histogram

%[ hsmax, upper ] = max(hs);

% iteratively search for the first maximum
%back = 1;
%while  back < upper  &&  hs(back+1) >= hs(back)
%    back = back + 1;
%end

[ hsmax, back ] = max(hs);
upper = back;

if sum( h(1:back) ) < numel(im) / 5
    fprintf(['  Background detection might have failed for ', filename, '\n']);
end

%% Make a figure to debug things
if make_plot

    x = (0:size(h)-1)';
    figure('Name','image_background', 'Position', [100 150 800 300]);
    axes('Position', [0.05 0.1 0.9 0.8] );
    hold on;
    xlim([0 size(h,1)]);
    set(gca, 'ytick', [])

    plot(x, h, 'g.' );
    plot(x, hs, 'k:', 'linewidth', 1);

    plot([back, back], ylim, 'b-', 'linewidth', 2);
    plot([upper, upper], ylim, 'k:');

    text(back, hsmax, sprintf(' %i', back), 'FontSize', 18 );
    title('Histogram of pixel values');

end

%% Gaussian fit of the distribution of dark pixels

    function f = gauss(x, a)
        f = exp( -(x-a(1)).^2 / a(2) );
    end

    function err = gauss_err(a)
        f = gauss(rx, a);
        s = sum(ry.*f) / sum(f.*f);
        err = sum( abs(f*s-ry) );   %robust fitting
        %err = norm( f*s-ry );      %least-square fitting
    end

if nargout > 1

    % Gaussian fit to obtain the variance of the black pixel values

    rwd = fix(back*0.8);
    ri = max(back-rwd, 1) : min(back+rwd, length(hs));
    rh = h(ri);

    rx = (ri-1)';
    ry = rh / sum(rh);
    esp = sum(rx.*ry);
    var = sum(rx.*rx.*ry);
    sigma = sqrt(var-esp*esp);

    if nargin > 1
        fprintf('Histogram  : mean %.2f, sigma %.2f\n', back, sigma);
    end

    [ pamG, pval, flag ] = fminsearch(@gauss_err, [ esp, 2*(var-esp*esp) ]);

    %if the fit is successful, we change our estimate of sigma:
    if flag == 1
        f     = gauss(rx, pamG);
        fitG  = f * ( sum(rh) * sum(ry.*f) / sum(f.*f) );

        backG = pamG(1);
        sigma = sqrt(pamG(2)/2);

        if nargin > 1
            fprintf('Gauss   fit: mean %.2f, sigma %.2f\n', backG, sigma);
        end
    end

    if make_plot
        %make an inset figure:
        axes('Position', [0.5 0.25 0.4 0.62]);
        plot( rx, h(1+rx), 'g.');
        hold on;
        plot( rx, hs(1+rx), 'b--');
        if exist('fitG', 'var')
            plot( rx, fitG, 'k-');
        end
        plot([back, back], ylim, 'k-');
        if exist('sigma', 'var')
            w = sqrt(2) * sigma;
            plot([back-w, back-w], ylim, 'k:');
            plot([back+w, back+w], ylim, 'k:');
        end
    end
end

%% Gamma fit (could be appropriate for low photon counts)

    function g = gamma_dis(x, a)
        %a is of dimension 2: a = { k, theta }
        %gamma density function
        g = exp( (a(1)-1)*log(x) - x/a(2) - a(1)*log(a(2)) ) / gamma(a(1));
    end

    function err = gamma_err(a)
        f = gamma_dis(rx, a);
        s = sum(ry.*f) / sum(f.*f);
        err = sum( abs(f*s-ry) );   %robust fitting
        %err = norm( f*s-ry );      %least-square fitting
    end

if nargout > 1

    [ pamP, pval, flag ] = fminsearch(@gamma_err, [ 5, esp/5 ]);

    if flag == 1

        f      = gamma_dis(rx, pamP);
        fitP   = f * ( sum(rh) * sum(ry.*f) / sum(f.*f) );

        %backP  = pamP(1) * pamP(2);  %mean-value
        sigmaP = sqrt(pamP(1)) * pamP(2); %sigma
        [v, id] = max(f);
        backP  = rx(id);  %position of peak

        if nargin > 1
            fprintf('Gamma   fit: mean %.2f, sigma %.2f\n', backP, sigmaP);
        end

        if make_plot
            plot( rx, fitP, 'r-');
        end

    end

end

%% Poisson fit

    function p = poisson(x, a)
        n = x ./ a(2);
        % discrete poisson distribution over n,
        p = exp( n .* ( 1 - log( n ./ a(1) ) ) + 0.5*log(n) );
    end

    function err = poisson_err(a)
        f = poisson(rx, a);
        s = sum(ry.*f) / sum(f.*f);
        err = sum( abs(f*s-ry) );   %robust fitting
        %err = norm( f*s-ry );      %least-square fitting
    end

if nargout > 1

    [ pamS, pval, flag ] = fminsearch(@poisson_err, [ esp / 10, 10 ]);

    if flag == 1

        grain  = pamS(2);
        backS  = pamS(2) * pamS(1);      %mean
        sigmaS = pamS(2) * sqrt(pamS(1));

        f      = poisson(rx, pamS);
        fitS   = f * ( sum(rh) * sum(ry.*f) / sum(f.*f) );

        if nargin > 1
            fprintf('Poisson fit: mean %.2f, sigma %.2f, grain %.2f\n', backS, sigmaS, grain);
        end

        if make_plot
            plot( rx, fitS, 'b-');
        end

    end

end

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function h = circle(x,y,r)
hold on
th = 0:pi/50:2*pi;
xunit = r * cos(th) + x;
yunit = r * sin(th) + y;
h = plot(xunit, yunit);
hold off

author: deng-haoyang

repo: ParNMPC


text:
function this = blank(this, fnamedat)
% Creates a blank datafile matching in the header in dimensions
% Will not erase existing datafile it it's there
% FORMAT this = blank(this)
%   Will create the datafile using fname and path
% FORMAT this = blank(this, fnamedat)
%   Will create the datafile using the provided name and path
% _________________________________________________________________________
% Copyright (C) 2011-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: blank.m 5648 2013-09-23 12:40:40Z vladimir $

if nargin == 1
    [p, f] = fileparts(fullfile(this));
    fnamedat = fullfile(p, [f '.dat']);
else
   [p, f, x] = fileparts(fnamedat);
   if isempty(p)
       p = path(this);
   end
   if isempty(x)
       x = '.dat';
   end
   fnamedat = fullfile(p, [f x]);
end

if exist(fnamedat, 'file')
    error('Data file exists. Use rmdata to delete.')
end

if isempty(this)
   error('All header dimensions should be >0');
end

this.data = file_array(fnamedat, size(this), 'float32-le');

initialise(this.data);

this = check(this);


text:
classdef GemPanel < GemUserInterfaceObject
    % GemPanel GEM class for a panel
    %
    %     Licence
    %     -------
    %     Part of GEM. https://github.com/tomdoel/gem
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = protected)
        BorderAxes
        InnerPosition
    end

    properties
        BackgroundColour
    end

    properties (SetObservable)
        BorderColour
    end

    properties (Dependent)
        LeftBorder
        RightBorder
        TopBorder
        BottomBorder
    end

    methods
        function obj = GemPanel(parent_handle)
            obj = obj@GemUserInterfaceObject(parent_handle);
            obj.BackgroundColour = obj.StyleSheet.BackgroundColour;
            obj.AddPostSetListener(obj, 'BorderColour', @obj.BorderColourChangedCallback);
        end

        function CreateGuiComponent(obj, position)
            obj.GraphicalComponentHandle = uipanel('Parent', obj.Parent.GetContainerHandle, 'BorderType', 'none', 'Units', 'pixels', ...
                'BackgroundColor', obj.BackgroundColour, 'ForegroundColor', obj.StyleSheet.TextPrimaryColour, 'ResizeFcn', '', 'Position', position);
        end

        function Resize(obj, position)
            Resize@GemUserInterfaceObject(obj, position);

            if ~isempty(obj.BorderAxes)
                obj.BorderAxes.Resize([1, 1, position(3), position(4)]);
            end

            inner_position = [1, 1, position(3), position(4)];

            if obj.LeftBorder
                inner_position(1) = inner_position(1) + 1;
                inner_position(3) = inner_position(3) - 1;
            end
            if obj.RightBorder
                inner_position(3) = inner_position(3) - 1;
            end
            if obj.TopBorder
                % FIXME: -2 instead of -1 to ensure top line is visible
                inner_position(4) = inner_position(4) - 2;
            end
            if obj.BottomBorder
                inner_position(2) = inner_position(2) + 1;
                inner_position(4) = inner_position(4) - 1;
            end

            obj.InnerPosition = inner_position;
        end

        function inner_position = GetInnerPosition(obj, inner_position)
            % Gets the position minus any borders
            if obj.LeftBorder
                inner_position(1) = inner_position(1) + 1;
                inner_position(3) = inner_position(3) - 1;
            end
            if obj.RightBorder
                inner_position(3) = inner_position(3) - 1;
            end
            if obj.TopBorder
                % FIXME: -2 instead of -1 to ensure top line is visible
                inner_position(4) = inner_position(4) - 2;
            end
            if obj.BottomBorder
                inner_position(2) = inner_position(2) + 1;
                inner_position(4) = inner_position(4) - 1;
            end
        end

        function set.LeftBorder(obj, border)
            if border
                obj.lazyCreateBorderAxes();
            end
            if ~isempty(obj.BorderAxes)
                obj.BorderAxes.LeftLine = border;
            end
        end

        function border = get.LeftBorder(obj)
            if isempty(obj.BorderAxes)
                border = false;
            else
                border = obj.BorderAxes.LeftLine;
            end
        end
        function set.RightBorder(obj, border)
            if border
                obj.lazyCreateBorderAxes();
            end
            if ~isempty(obj.BorderAxes)
                obj.BorderAxes.RightLine = border;
            end
        end

        function border = get.RightBorder(obj)
            if isempty(obj.BorderAxes)
                border = false;
            else
                border = obj.BorderAxes.RightLine;
            end
        end
        function set.TopBorder(obj, border)
            if border
                obj.lazyCreateBorderAxes();
            end
            if ~isempty(obj.BorderAxes)
                obj.BorderAxes.TopLine = border;
            end
        end

        function border = get.TopBorder(obj)
            if isempty(obj.BorderAxes)
                border = false;
            else
                border = obj.BorderAxes.TopLine;
            end
        end
        function set.BottomBorder(obj, border)
            if border
                obj.lazyCreateBorderAxes();
            end
            if ~isempty(obj.BorderAxes)
                obj.BorderAxes.BottomLine = border;
            end
        end

        function border = get.BottomBorder(obj)
            if isempty(obj.BorderAxes)
                border = false;
            else
                border = obj.BorderAxes.BottomLine;
            end
        end
    end

    methods (Access = protected)
        function BorderColourChangedCallback(obj, ~, ~, ~)
            if ~isempty(obj.BorderAxes)
                obj.BorderAxes.LineColour = obj.BorderColour;
            end
        end

        function lazyCreateBorderAxes(obj)
            if isempty(obj.BorderAxes)
                obj.BorderAxes = GemLineAxes(obj);
                if ~isempty(obj.BorderColour)
                    obj.BorderAxes.LineColour = obj.BorderColour;
                end
                obj.AddChild(obj.BorderAxes);
            end
        end

    end

end

author: tomdoel

repo: pulmonarytoolkit


text:
function varargout = loadxml(filename,varargin)
%LOADXML Load workspace variables from disk (XML file).
%  LOADXML FILENAME retrieves all variables from a file given a full
%  pathname or a MATLABPATH relative partial pathname (see PARTIALPATH).
%  If FILENAME has no extension LOAD looks for FILENAME and FILENAME.xml
%  and treats it as an XML file.
%
%  LOAD, by itself, uses the XML file named 'matlab.xml'. It is an error
%  if 'matlab.xml' is not found.
%
%  LOAD FILENAME X loads only X.
%  LOAD FILENAME X Y Z ... loads just the specified variables.  The
%  wildcard '*' loads variables that match a pattern.
%  Requested variables from FILENAME are created in the workspace.
%
%  S = LOAD(...) returns the contents of FILENAME in variable S. S is
%  a struct containing fields matching the variables retrieved.
%
%  Use the functional form of LOAD, such as LOAD('filename'), when the
%  file name is stored in a string, when an output argument is requested,
%  or if FILENAME contains spaces.
%
%  See also LOAD, XML2MAT, XMLTREE.

%  Copyright 2003 Guillaume Flandin.
%  $Revision: 6530 $  $Date: 2003/07/10 13:50 $

%  $Id: loadxml.m 6530 2015-08-21 14:43:52Z guillaume $

if nargin == 0
    filename = 'matlab.xml';
    fprintf('\nLoading from: %s\n\n',filename);
end

if ~ischar(filename)
    error('[LOADXML] Argument must contain a string.');
end

if ~exist(filename,'file')
    filename = [filename '.xml'];
    if ~exist(filename,'file')
        error(sprintf(...
        '[LOADXML] Unable to read file %s: file does not exist',filename));
    end
end

if nargout > 1,
    error('[LOADXML] Too many output arguments.');
end

t = xmltree(filename);

uid = children(t,root(t));

if nargout == 1
    % varargout{1} = struct([]); % Matlab 6.0 and above
end

flagfirstvar = 1;
for i=1:length(uid)
    if strcmp(get(t,uid(i),'type'),'element')
        vname = get(t,uid(i),'name');
        % TODO % No need to parse the whole tree
        if isempty(varargin) | ismember(varargin,vname)
            v = xml_create_var(t,uid(i));
            if nargout == 1
                if flagfirstvar
                    varargout{1} = struct(vname,v);
                    flagfirstvar = 0;
                else
                    varargout{1} = setfield(varargout{1},vname,v);
                end
            else
                assignin('caller',vname,v);
            end
        end
    end
end

%=======================================================================
function v = xml_create_var(t,uid)
    type = attributes(t,'get',uid,'type');
    sz   = str2num(attributes(t,'get',uid,'size'));

    switch type
        case 'double'
            v = str2num(get(t,children(t,uid),'value'));
            if ~isempty(sz)
                v = reshape(v,sz);
            end
        case 'sparse'
            u = children(t,uid);
            for k=1:length(u)
                if strcmp(get(t,u(k),'name'),'row')
                    i = str2num(get(t,children(t,u(k)),'value'));
                elseif strcmp(get(t,u(k),'name'),'col')
                    j = str2num(get(t,children(t,u(k)),'value'));
                elseif strcmp(get(t,u(k),'name'),'val')
                    s = str2num(get(t,children(t,u(k)),'value'));
                end
            end
           v = sparse(i,j,s,sz(1),sz(2));
        case 'struct'
            u = children(t,uid);
            v = []; % works with Matlab < 6.0
            for i=1:length(u)
                s(1).type = '()';
                s(1).subs = {str2num(attributes(t,'get',u(i),'index'))};
                s(2).type = '.';
                s(2).subs = get(t,u(i),'name');
                v = subsasgn(v,s,xml_create_var(t,u(i)));
            end
            if isempty(u),
                v = struct([]); % Need Matlab 6.0 and above
            end
        case 'cell'
            v = cell(sz);
            u = children(t,uid);
            for i=1:length(u)
                v{str2num(attributes(t,'get',u(i),'index'))} = ...
                    xml_create_var(t,u(i));
            end
        case 'char'
            if isempty(children(t,uid))
                v = '';
            else
                v = get(t,children(t,uid),'value');
            end
            try % this can fail if blank spaces are lost or entity escaping
                if ~isempty(sz)
                    if sz(1) > 1
                        v = reshape(v,fliplr(sz))'; % row-wise order
                    else
                        v = reshape(v,sz);
                    end
                end
            end
        case  {'int8','uint8','int16','uint16','int32','uint32'}
            % TODO % Handle integer formats
            warning(sprintf('%s matrices not handled.',type));
            v = 0;
        otherwise
            try,
                v = feval(class(v),get(t,uid,'value'));
            catch,
                warning(sprintf(...
                '[LOADXML] Cannot convert from XML to %s.',type));
            end
    end


text:
function G = ipopt_callback_dg(x,model)%returnStructOnly,model)

global latest_x_g
global latest_G
global latest_g

x = x(:);
if isequal(x,latest_x_g)
    % Jacobian was computed already in the call for constraints
    G = latest_G;
else
    % Compute the nonlinear terms in the constraints
    [g,geq,dg,dgeq] = fmincon_con_liftlayer(x,model);

    % Append with linear terms
    if isempty(dg)
        G = dgeq';
    else
        G = [dg';dgeq'];
    end

    if ~isempty(model.A)
        G = [G;model.A];
    end
    if ~isempty(model.Aeq)
        G = [G;model.Aeq];
    end
    G = sparse(G);
end
if model.dense
    G = full(G);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%% Enumerated
% Abstract superclass for enumerated quantities.

%%% Description
% |Enumerated| is the superclass of classes representing enumerated quantities.
% Each subclass of |Enumerated| allows only instances predefined by
% |enumeration|.  For more details about |enumeration|, see
% <matlab:web(['jar:file://',fullfile(matlabroot,'help','techdoc','help.jar!','matlab_oop','br4imrp.html')],'-helpbrowser')
% Enumerations>.
%
% |Enumerated| supports conversion of its instance into a string by |char()|,
% and into an integer by |int()|.
%
% Each concrete subclass of |Enumerated| allows comparison between its
% instances, i.e., supports relational operators <, <=, ==, >=, and >.  Also,
% the result consistent with this comparison is returned by |sort()| on a vector
% of these instances.

%%% Concrete Subclasses
% * <Axis.html |Axis|>: Catesion axes
% * <BC.html |BC|>: boundary conditions
% * <Dir.html |Dir|>: horizontal and vertical directions
% * <GT.html |GT|>: grid kinds
% * <PhysQ.html |PhysQ|>: physical quantities
% * <PML.html |PML|>: PML kinds
% * <Sign.html |Sign|>: negative and positive signs

%%% Methods (static)
% Below, |Enumerated| must be replaced by a concrete subclass.
%
% * |Enumerated.elems()|: returns a row vector of all instances
% * |Enumerated.count()|: returns the number of all instances

%%% Methods (public)
% Below, |enumerated| is an instance of a concrete subclass of |Enumerated|.
%
% * |char(enumerated)|: name of |enumerated|
% * |int(enumerated)|: integer representation of |enumerated|

%%% Example
% Below, the methods of |Enumerated| are demonstrated using a concrete subclass
% <Axis.html |Axis|>.
%
%   % Test conversion to integers and strings.
%   fprintf('Total # of Axis objects: %d\n', Axis.count);
%   for w = Axis.elems
%       fprintf('Name of Axis object %d: %s\n', int(w), char(w));
%   end
%
%   % Test sort().
%   unsorted = [Axis.z, Axis.x, Axis.y];
%   fprintf('\nBefore sort: ');
%   for w = unsorted
%       fprintf('%s  ', char(w));
%   end
%
%   fprintf('\nAfter sort: ');
%   sorted = sort(unsorted);
%   for w = sorted
%       fprintf('%s  ', char(w));
%   end
%   fprintf('\n');

classdef Enumerated < handle
	properties (SetAccess = immutable, GetAccess = private)
		name
	end

	methods (Abstract, Static)
		elems = elems(ind)
		count = count()
	end

%   (Deleted because it uses Axis rather than Enumerated, and also because
%   loading does not generate warnings without loadobj(); only saveobj() seems
%   necessary.)
% 	methods(Static)
% 		function obj = loadobj(S)
% 			% If saveobj() and loadobj() are not implemented, a warning is
% 			% issued when loading a saved Enumerated object.
% 			obj = Axis(char(S));
% 		end
% 	end

	methods
		function this = Enumerated(name)
			chkarg(ischar(name), '"name" should be string.');
			this.name = name;
		end

		function S = saveobj(this)
			% If saveobj() and loadobj() are not implemented, a warning is
			% issued when loading a saved Enumerated object.
			S.name = char(this);
		end

		function name = char(this)
			% This function returned name = this.name originally, but it is
			% modified to handle arguments given as arrays.
% 			name = this.name;
			name = cell2mat(arrayfun(@(x) x.name, this, 'UniformOutput', false));
		end

		function n = int(this)
			% This function returned n = find(this.elems==this) originally, but
			% it is modified to handle arguments given as arrays.
% 			[~, n] = ismember(this, this.elems, 'legacy');
% 			[~, n] = ismember(this, this.elems);
			n = arrayfun(@(x) find(this(1).elems==x), this);
		end

		function ind = subsindex(this)
			ind = int(this) - 1;
		end

		function c = plus(this, another)
				c = int(this) + double(another);
		end

		function [sorted, ind] = sort(this, varargin)
			n = length(this);
			nums = NaN(1, n);
			for i = 1:n
				nums(i) = int(this(i));
			end
			[~, ind] = sort(nums, varargin{:});
			sorted = this(ind);
		end

		function truth = lt(this, another)
			cn = class(this);  % class name
			chkarg(istypesizeof(another, cn), 'cannot compare %s with %s.', cn, class(another));
			truth = int(this) < int(another);
		end

		function truth = le(this, another)
			cn = class(this);  % class name
			chkarg(istypesizeof(another, cn), 'cannot compare %s with %s.', cn, class(another));
			truth = int(this) <= int(another);
		end

		function truth = gt(this, another)
			cn = class(this);  % class name
			chkarg(istypesizeof(another, cn), 'cannot compare %s with %s.', cn, class(another));
			truth = int(this) > int(another);
		end

		function truth = ge(this, another)
			cn = class(this);  % class name
			chkarg(istypesizeof(another, cn), 'cannot compare %s with %s.', cn, class(another));
			truth = int(this) >= int(another);
		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
 function [basis, Hw] = deconv_z(psf, fun, t)
% psf: 3-tap symmetric psf: [a _b_ a]
% H(z) = 1 / (az + b + a/z)

if numel(psf) == 1
	basis = fun(t) / psf;
	Hw = @(t) psf * ones(size(t));
return
end

if length(psf) ~= 3, error 'not done', end
a = psf(1);
b = psf(2);
c = b / a / 2;
p = -c + sign(c) * sqrt(c^2 - 1); % pole

scale = 1/a * 1/(p - 1/p);
basis = fun(t);
for n=1:9
	basis = basis + p^n * (fun(t-n) + fun(t+n));
end
basis = scale * basis;

Hw = @(om) 1 ./ (b + 2 * a * cos(om));

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function varargout=upolyhedron(varargin)
%
% For calling details please see v_upolyhedron.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_upolyhedron.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: upolyhedron.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_upolyhedron(varargin{:});
else
    v_upolyhedron(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function engine = ff_inf_engine(bnet)
% FF_INF_ENGINE Factored frontier inference engine for DBNs
% engine = ff_inf_engine(bnet)
%
% The model must be topologically isomorphic to an HMM.
% In addition, each hidden node is assumed to have at most one observed child,
% and each observed child is assumed to have exactly one hidden parent.
%
% For details of this algorithm, see
%  "The Factored Frontier Algorithm for Approximate Inference in DBNs",
%   Kevin Murphy and Yair Weiss, UAI 2001.
%
% THIS IS HIGHLY EXPERIMENTAL CODE!

ss = length(bnet.intra);
onodes = bnet.observed;
hnodes = mysetdiff(1:ss, onodes);

[persistent_nodes, transient_nodes] = partition_dbn_nodes(bnet.intra, bnet.inter);
assert(isequal(onodes, transient_nodes));
assert(isequal(hnodes, persistent_nodes));

engine.onodes = onodes;
engine.hnodes = hnodes;
engine.marginals = [];
engine.fwd = [];
engine.back = [];
engine.CPDpot = [];
engine.filter = [];

obschild = zeros(1,ss);
for i=engine.hnodes(:)'
  %ocs = myintersect(children(bnet.dag, i), onodes);
  ocs = children(bnet.intra, i);
  assert(length(ocs) <= 1);
  if length(ocs)==1
    obschild(i) = ocs(1);
  end
end
engine.obschild = obschild;

engine = class(engine, 'ff_inf_engine', inf_engine(bnet));


text:
function asa299_test01 ( )

%*****************************************************************************80
%
%% TEST01 tests SIMPLEX_LATTICE_POINT_NEXT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 January 2007
%
%  Author:
%
%    John Burkardt
%
  n = 4;
  t = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  SIMPLEX_LATTICE_POINT_NEXT generates lattice points\n' );
  fprintf ( 1, '  in the simplex\n' );
  fprintf ( 1, '    0 <= X\n' );
  fprintf ( 1, '    sum ( X(1:N) ) <= T\n' );
  fprintf ( 1, '  Here N = %d\n', n );
  fprintf ( 1, '  and T =  %d\n', t );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     Index        X(1)      X(2)      X(3)      X(4)\n' );
  fprintf ( 1, '\n' );

  more = 0;
  x = [];

  i = 0;

  while ( 1 )

    [ more, x ] = simplex_lattice_point_next ( n, t, more, x );

    i = i + 1;

    fprintf ( 1, '  %8d  ', i );
    for j = 1 : n
      fprintf ( 1, '  %8d', x(j) );
    end
    fprintf ( 1, '\n' );

    if ( ~more )
      break
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function obj=reportTestOutcome(obj, test_outcome)
% reports the test outcome to the stream of the current object
%
% obj=reportTestOutcome(obj, test_outcome)
%
% Inputs:
%   obj             MOxUnitTestReport instance
%   test_outcome    MOxUnitTestReport instance with the outcome to report.
%                   This function will print the outcome to the stream used
%                   when obj was instantiated. As a side effect,
%                   test_outcome is added to the report
%
% Output:
%   obj             MOxUnitTestReport instance with the test_outcome added

    obj=addTestOutcome(obj,test_outcome);

    verbosity=obj.verbosity;

    outcome_str=getProgressStr(test_outcome, verbosity);
    stream=obj.stream;

    postfix=''; % the default

    max_line_length=60;
    if mod(countTestOutcomes(obj),max_line_length)==0 && verbosity==1
        % when showing single characters, print a newline every sixty
        % characters
        postfix=sprintf('\n');
    end

    fprintf(stream,'%s%s',outcome_str,postfix);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
clc;
clear;
close all;
restoredefaultpath;

sdpnalpath  = '../../SDPNAL+v1.0';
pgdpath     = '../STRIDE';
utilspath   = '../utils';

addpath(genpath(pwd));
addpath(genpath(utilspath))
%% construct space of n1 x n2 hankel matrices (n1 <= n2)
n1 = 20;
n2 = 20;
[S,k,Scell] = hankel_struct(n1,n2);
% generate random hankel matrix
u1 = randn(k,1);
U1 = applyAffineMapCell(Scell,u1);

%% Convert the nearest hankel matrix problem to SDP
addpath(genpath(pgdpath));
SDP     = nearest_hankel_sdp(S,u1);
fprintf('SDP size: n = %d, m = %d.\n\n\n',SDP.n,SDP.m);
rmpath(genpath(pgdpath));

%% Optimistic initialization using SLRA
s.m             = n1;
s.n             = n2;
tol             = 0;
[zhtls,Hu]      = htls(s.m-1,u1,'linear'); % good initialization
opts.Rini       = zhtls';
opts.maxiter    = 5000;
opts.tol        = tol;
opts.epsrel     = tol;
opts.epsabs     = tol;
opts.epsgrad    = tol;
opts.method     = 'll'; % 'll': Levenberg–Marquardt, 'qb': BFGS
[uslra, info]   = slra(u1, s, s.m-1, opts);
zslra           = info.Rh'; zslra = zslra/norm(zslra);
Uslra           = applyAffineMapCell(Scell,uslra);
ztUnorm         = norm(zslra'*Uslra); % measure of rank deficientness
fprintf('SLRA: norm(zt*U) = %3.2e.\n',ztUnorm);
% lift to SDP solution
xtld            = kron([uslra;1],zslra);
X0              = {xtld * xtld'};

%% STRIDE
addpath(genpath(pgdpath))
pgdopts.pgdStepSize     = 10;
pgdopts.maxiterSGS      = 300;
pgdopts.maxiterLBFGS    = 1000;
pgdopts.SDPNALpath      = sdpnalpath;
pgdopts.tolADMM         = 5e-5;
pgdopts.phase1          = 1;
pgdopts.rrOpt           = 1:3;
pgdopts.lbfgsmemory     = 10;
pgdopts.tolPGD          = 1e-8;
pgdopts.rrFunName       = 'rr_hankel';
rrPar.m = n1; rrPar.n = n2; rrPar.k = k; rrPar.theta = u1; rrPar.S = S;
pgdopts.rrPar           = rrPar;
[outPGD,Xopt,yopt,Sopt] = PGDSDP(SDP.blk, SDP.At, SDP.b, SDP.C, X0, pgdopts);
f_sdp                   = outPGD.pobj;
time_pgd                = outPGD.totaltime;

[z,u,f_est]     = recover_solution(Xopt,SDP.C,n1,k);
eta             = abs(f_est-f_sdp) / (1+abs(f_est)+abs(f_sdp)); % relative suboptimality gap
U               = applyAffineMapCell(Scell,u);
ztUnorm         = norm(z'*U); % measure of rank deficientness
fprintf('norm(zt*U) = %3.2e, eta = %3.2e.\n',ztUnorm,eta);

rmpath(genpath(pgdpath))
fprintf('\n\n\n\n\n')

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function value = i4_log_r8 ( x, b )

%*****************************************************************************80
%
%% I4_LOG_R8 returns the integer part of the logarithm base ABS(B) of ABS(X).
%
%  Example:
%
%    If B is greater than 1, and X is positive:
%
%    if 1/B^2  <  X <= 1/B   I4_LOG_R8(X) = -1,
%    if 1/B    <  X <= 1     I4_LOG_R8(X) = 0,
%    if 1      <= X <  B,    I4_LOG_R8(X) = 0,
%    if B      <= X <  B^2   I4_LOG_R8(X) = 1,
%    if B^2    <= X <  B^3   I4_LOG_R8(X) = 2.
%
%    For positive I4_LOG_R8(X), it should be true that
%
%      ABS(B)^I4_LOG_R8(X) <= ABS(X) < ABS(B)^(I4_LOG_R8(X)+1).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X, the number whose logarithm base B is desired.
%    If X is 0, then I4_LOG_R8 is returned as -HUGE().
%
%    Input, real B, the absolute value of the base of the
%    logarithms.  B must not be -1, 0, or 1.
%
%    Output, integer VALUE, the integer part of the logarithm
%    base abs(B) of X.
%
  i4_huge = 2147483647;

  x = floor ( x );

  if ( x == 0 )
    value = - i4_huge;
    return
  end

  b = abs ( b );
  value = 0;

  if ( b == 1.0 )
    return
  end

  if ( b == 0.0 )
    return
  end

  x = abs ( x );

  if ( b < 1.0 )
    value_sign = -1;
    b = 1.0 / b;
  else
    value_sign = +1;
  end

  if ( 1.0 <= x && x < b )
    value = value_sign * value;
    return
  end

  while ( b < x )
    x = x / b;
    value = value + 1;
  end

  while ( x * b <= 1.0 )
    x = x * b;
    value = value - 1;
  end
%
%  If the absolute value of the base was less than 1, we inverted
%  earlier.  Now negate the logarithm to account for that.
%
  value = value_sign * value;

  return
end


text:
function [ x, w ] = j_quadrature_rule ( n, alpha, beta )

%*****************************************************************************80
%
%% J_QUADRATURE_RULE: Gauss-Jacobi quadrature based on J(n,a,b,x).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2012
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer N, the order.
%
%    Input, real ALPHA, BETA, the parameters.
%    -1 < ALPHA, BETA.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%
  ab = alpha + beta;
  abi = 2.0 + ab;
%
%  Define the zero-th moment.
%
  zemu = 2.0^( ab + 1.0 ) * gamma ( alpha + 1.0 ) ...
    * gamma ( beta + 1.0 ) / gamma ( abi );
%
%  Define the Jacobi matrix.
%
  x = zeros ( n, 1 );
  bj = zeros ( n, 1 );

  x(1) = ( beta - alpha ) / abi;
  bj(1) = 4.0 * ( 1.0 + alpha ) * ( 1.0 + beta ) ...
    / ( ( abi + 1.0 ) * abi * abi );
  a2b2 = beta * beta - alpha * alpha;

  for i = 2 : n
    abi = 2.0 * i + ab;
    x(i) = a2b2 / ( ( abi - 2.0 ) * abi );
    abi = abi^2;
    bj(i) = 4.0 * i * ( i + alpha ) * ( i + beta ) * ( i + ab ) ...
      / ( ( abi - 1.0 ) * abi );
  end
  bj(1:n) =  sqrt ( bj(1:n) );

  w = zeros ( n, 1 );
  w(1) = sqrt ( zemu );
%
%  Diagonalize the Jacobi matrix.
%
  [ x, w ] = imtqlx ( n, x, bj, w );

  w(1:n) = w(1:n).^2;

  return
end


text:
function y = prtRvUtilDiscretePdf(X,values,probabilities)
%y = discretepdf(X,values,probabilities);
%	Return the PDF of a discrete distribution with values and
%	probabilities.

y = zeros(size(X));
for i = 1:size(X(:),1);
    cI = find(X(i) == values);
    if isempty(cI)
        y(i) = 0;
    else
        y(i) = probabilities(cI);
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [dat, w] = ft_preproc_denoise(dat, refdat, hilbertflag)

% FT_PREPROC_DENOISE performs a regression of the matrix dat onto
% refdat, and subtracts the projected data. This is for the
% purpose of removing signals generated by coils during continuous
% head motion tracking, for example.
%
% Use as
%   [dat] = ft_preproc_denoise(dat, refdat, hilbertflag)
% where
%   dat         data matrix (Nchan1 X Ntime)
%   refdat      data matrix (Nchan2 X Ntime)
%   hilbertflag specifying to regress out the real and imaginary parts of
%                 the hilbert transformed signal. Only meaningful for narrow
%                 band reference data
%
% See also PREPROC

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

if nargin<3
  hilbertflag = 0;
end

% preprocessing fails on channels that contain NaN
if any(isnan(dat(:))) || any(isnan(refdat(:)))
  ft_warning('FieldTrip:dataContainsNaN', 'data contains NaN values');
end

n1 = size(dat,2);
n2 = size(refdat,2);
m1 = mean(dat,2);
m2 = mean(refdat,2);

%remove mean
refdat  = refdat-m2(:,ones(n2,1));
tmpdat  = dat-m1(:,ones(n1,1));

%do hilbert transformation
if hilbertflag>0
  hrefdat = hilbert(refdat')';
  refdat  = [real(hrefdat);imag(hrefdat)];
end

c12 = tmpdat*refdat'; %covariance between signals and references
c1  = refdat*refdat'; %covariance between references and references
w   = (pinv(c1)*c12')'; %regression weights

%subtract
dat = dat-w*refdat;

author: spm


text:
function imgOut = KimKautzConsistentTMO(img, Ld_max, Ld_min, KK_c1, KK_c2)
%
%
%      imgOut = KimKautzConsistentTMO(img, Ld_max, Ld_min, KK_c1, KK_c2)
%
%
%       Input:
%           -img: input HDR image
%           -Ld_max: max luminance of the LDR monitor in cd/m^2
%           -Ld_min: max luminance of the LDR monitor in cd/m^2
%           -KK_c1: this parameter adjusts the shape of Gaussian fall-off
%           within the width of tis characteristic curve. It influcences
%           the resulting brightness and local details of the tone-mapped
%           image. A good value is 3.0 (tradeoff between compression and
%           lost details)
%           -KK_c2: the ratio between the dynamic range (in log10) of an
%           8-bit imag (2.4) and the dynamic range (in log10) of the
%           LDR monitor for visualization
%
%       Output:
%           -imgOut: output tone mapped image in linear domain
%
%     Copyright (C) 2013-15  Francesco Banterle
%
%     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/>.
%
%     The paper describing this technique is:
%     "Consistent Tone Reproduction"
% 	  by Min H. Kim, Jan Kautz
%     in CGIM '08 Proceedings of the Tenth IASTED
%     International Conference on Computer Graphics and Imaging  2008
%

checkNegative(img);

if(~exist('Ld_max', 'var'))
    Ld_max = 300; %300 cd/m^2
end

if(~exist('Ld_min', 'var'))
    Ld_min = 0.3; %0.3 cd/m^2
end

if(~exist('KK_c1', 'var'))
    KK_c1 = 3.0; %as in the original paper
end

if(~exist('KK_c2', 'var'))
    KK_c2 = 0.5; %as in the original paper
end

L = lum(img);

L_log = log(L + 1e-6);

mu = mean(L_log(:));

maxL = max(L_log(:));
minL = min(L_log(:));

maxLd = log(Ld_max);
minLd = log(Ld_min);

k1 = (maxLd - minLd) / (maxL - minL);

d0 = maxL - minL;
sigma = d0 / KK_c1;

sigma_sq_2 = (sigma^2) * 2;
w = exp(-(L_log - mu).^2 / sigma_sq_2);
k2 = (1 - k1) * w + k1;

Ld = exp(KK_c2 * k2 .* (L_log - mu) + mu);

%robust min and max
minLd = MaxQuart(Ld, 0.01);
maxLd = MaxQuart(Ld, 0.99);

Ld(Ld < minLd) = minLd;
Ld(Ld > maxLd) = maxLd;

Ld = (Ld - minLd) / (maxLd - minLd);

%change luminance
imgOut = ChangeLuminance(img, L, Ld);
imgOut = RemoveSpecials(imgOut);

end

author: banterle


text:
function [y,deriv] = sum_ai_f_of_w_i(w,a,f,b)
% This is an MV2DF. See MV2DF_API_DEFINITION.readme.
%
% Does y = sum_i a_i f(w_i) + b, where f is non-linear.
%
%Notes:
%
%  f is a function handle, with behaviour as demonstrated in the test code
%  of this function.
%
%  b is optional, defaults to 0 if omitted

if nargin==0
    test_this();
    return;
end

if ~exist('b','var')
    b = 0;
end

if isempty(w)
    y = @(w)sum_ai_f_of_w_i(w,a,f,b);
    return;
end

if isa(w,'function_handle')
    outer = sum_ai_f_of_w_i([],a,f,b);
    y = compose_mv(outer,w,[]);
    return;
end

ntot = length(a);
nz = find(a~=0);
a = a(nz);

if nargin==1
    y = f(w(nz));
else
    [y,dfdw,f2] = f(w(nz));
    deriv = @(Dy) deriv_this(Dy,dfdw.*a,f2,a,nz,ntot);
end
y = y(:);
y = a.'*y + b;

function [g,hess,linear] = deriv_this(Dy,g0,f2,a,nz,ntot)
g = zeros(ntot,1);
g(nz) = Dy*g0(:);
hess = @(d) hess_this(d,g0,f2,Dy,a,nz,ntot);
linear = false;

function [h,Jd] = hess_this(d,g0,f2,Dy,a,nz,ntot)
d = d(nz);
hnz = f2();
hnz = hnz(:).*d(:);
h = zeros(ntot,1);
h(nz) = Dy*(hnz.*a);
if nargout>1
    Jd = g0.'*d(:);
end

function [y,ddx,f2] = test_f(x)
y = log(x);
if nargout>1
    ddx = 1./x;
    f2 = @() -1./(x.^2);
end

function test_this()
n = 10;
a = randn(n,1);
a = bsxfun(@max,a,0);
b = 5;
f = sum_ai_f_of_w_i([],a,@(x)test_f(x),b);

w = 1+rand(n,1);
test_MV2DF(f,w);

author: nesl


text:
function [est, x, k] = pnorm(A, p, tol, prnt)
%PNORM   Estimate of matrix p-norm (1 <= p <= inf).
%        [EST, x, k] = PNORM(A, p, TOL) estimates the Holder p-norm of a
%        matrix A, using the p-norm power method with a specially
%        chosen starting vector.
%        TOL is a relative convergence tolerance (default 1E-4).
%        Returned are the norm estimate EST (which is a lower bound for the
%        exact p-norm), the corresponding approximate maximizing vector x,
%        and the number of power method iterations k.
%        A nonzero fourth input argument causes trace output to the screen.
%        If A is a vector, this routine simply returns NORM(A, p).
%
%        See also NORM, NORMEST, NORMEST1.

%        Note: The estimate is exact for p = 1, but is not always exact for
%        p = 2 or p = inf.  Code could be added to treat p = 2 and p = inf
%        separately.
%
%        Calls DUAL.
%
%        Reference:
%        N. J. Higham, Estimating the matrix p-norm, Numer. Math.,
%             62 (1992), pp. 539-555.
%        N. J. Higham, Accuracy and Stability of Numerical Algorithms,
%           Second edition, Society for Industrial and Applied Mathematics,
%           Philadelphia, PA, 2002; sec. 15.2.

if nargin < 2, error('Must specify norm via second parameter.'), end
[m,n] = size(A);
if min(m,n) == 1, est = norm(A,p); return, end

if nargin < 4, prnt = 0; end
if nargin < 3 | isempty(tol), tol = 1e-4; end

% Stage I.  Use Algorithm OSE to get starting vector x for power method.
% Form y = B*x, at each stage choosing x(k) = c and scaling previous
% x(k+1:n) by s, where norm([c s],p)=1.

sm = 9;  % Number of samples.
y = zeros(m,1); x = zeros(n,1);

for k=1:n

    if k == 1
       c = 1; s = 0;
    else
       W = [A(:,k) y];

       if p == 2   % Special case.  Solve exactly for 2-norm.
          [U,S,V] = svd(full(W));
          c = V(1,1); s = V(2,1);

       else

          fopt = 0;
          for th=linspace(0,pi,sm)
              c1 = cos(th); s1 = sin(th);
              nrm = norm([c1 s1],p);
              c1 = c1/nrm; s1 = s1/nrm;   % [c1 s1] has unit p-norm.
              f = norm( W*[c1 s1]', p );
              if f > fopt
                 fopt = f;
                 c = c1; s = s1;
              end
          end

       end
    end

    x(k) = c;
    y = x(k)*A(:,k) + s*y;
    if k > 1, x(1:k-1) = s*x(1:k-1); end

end

est = norm(y,p);
if prnt, fprintf('Alg OSE: %9.4e\n', est), end

% Stage II.  Apply Algorithm PM (the power method).

q = dual(p);
k = 1;

while 1

    y = A*x;
    est_old = est;
    est = norm(y,p);

    z = A' * dual(y,p);

    if prnt
        fprintf('%2.0f: norm(y) = %9.4e,  norm(z) = %9.4e', ...
                 k, norm(y,p), norm(z,q))
        fprintf('  rel_incr(est) = %9.4e\n', (est-est_old)/est)
    end

    if ( norm(z,q) <= z'*x | abs(est-est_old)/est <= tol ) & k > 1
       return
    end

    x = dual(z,q);
    k = k + 1;

end

author: opencobra

repo: cobratoolbox


text:
function selected = atlasSelectAtlas
%
%     selected = atlasSelectAtlas
%
% Author: BW, AAB,
% Purpose:
%    List dialogue to let the user select from the existing atlases (if selected >=1)
%    indicate the user wants a new atlas (selected == 0)
%    or the user decided to cancel, selected = [].
%

global dataTYPES;

atlasTypeNum = existDataType('Atlases',[],0);
str{1} = 'New Atlas';
if all(atlasTypeNum > 0)   % only if we have at least one atlas
	for ii=2:(length(atlasTypeNum)+1);
        str{ii} = dataTYPES(atlasTypeNum(ii-1)).name;
    end
end

[s,valid] = listdlg('PromptString','Select an Atlas:',...
    'SelectionMode','single',...
    'ListString',str);

%
if valid
    s = s - 1;
    if s > 0,  selected = atlasTypeNum(s);
    else       selected = 0;
    end
else
    selected = [];
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% sound2gyro.m
% evaluate effect of sound signal on phone internal gyro
%  publish('sound2gyro.m','doc');

Fs = 200;                   % 200 sps default
[filenam, pathnam] = uigetfile('*.txt', 'get measurements data file');
fid = fopen(filenam,'r');
% 0.0
% 0.000000 -0.044288
% 5046272.000000 -0.044288
% L2:L212 =[txyz] (tinnsec)
% repeat n times
n=20;
f_in = zeros(n,1);
txyz = zeros(200,4,n);
frewind(fid),
i=1;
f_in(i) = str2double(fgetl(fid)); % read a line, 1st is noise freq
while ~isnan(f_in(i))
    linenums = 4;
    cntlines = 1;
    while linenums == 4
        % # of nums per line
        % counting for sampnum
        % loop to read one noise sequence
        tt = fgetl(fid);
        ttnum = str2num(tt);
        linenums = length(ttnum);      %
        if linenums ==4,
            txyz(cntlines,:,i) = ttnum;
            cntlines = cntlines+1;
        end;
    end;
    % get next line chars
    % make numbers
    sampnum = cntlines - 1;             % last was 1 or nan
    i = i+1;
    if isempty(ttnum), f_in(i)=nan; else f_in(i) = ttnum; end;
end

f_in(i:end)=[];             % clean the rest
txyz(:,:,i:end) = [];
n = i-1;                    % last one was NAN
% set  # of columns for plots
if n>5
    pcols = 2;
else
    pcols = 1;
end;
prows = n/pcols;
fclose(fid);

filenam = 'RAFDOCS-#7775044Magn.TXT';

if findstr(filenam,'77750')     % .. 28 | 44
    savetxyz = txyz;
    resamp;     % returns newtxyz, f_in, Fs, nlen
    len = min(nlen);
    txyz = zeros(len,4, length(f_in));
    for kf = 1:n,
        txyz(:,:,kf) = newtxyz{kf}(1:len,:);
    end
    xyz = txyz(:,2:4,:);
    time = squeeze(txyz(:,1,:));
    t = time;
else   %% verify, update ->  len, n
    len = 199;
    % clean
    % orig in nsec, Fs
    % clean
    t = repmat((0:len-1)'/Fs,1,n);
end;

xyzstr = 'XYZ';

xyz = txyz(1:len,2:4,:);      % 199x3xn
% x10 = squeeze(txyz(1:len,2,:));
% y10 = squeeze(txyz(1:len,3,:));
% z10 = squeeze(txyz(1:len,4,:));

for d = 1:3            % dimensions
    figure
    mn = min(min(squeeze(xyz(:,d,:))));
    mx = max(max(squeeze(xyz(:,d,:))));
    for i=1:n
        subplot(prows, pcols,i)
        curr = xyz(:,d,i);
        plot(t(:,i),curr,'.-')
        axis([0 t(end,i) mn mx]);
        title([ xyzstr(d) '-meas. for freq = ' num2str(f_in(i)) ' Hz'])
        xlabel('time [s]'), grid on
    end;
    if pcols>1,
        set(gcf,'posi',[80 80 750 1000])
    end;
    emarkpl( filenam)
end;

figure;
plot(diff(t(:,:)))
xlabel ('sample #')
ylabel('diff(time) [sec]')
title('Sample time interval variation')
dtser = diff(t);
[Val, Bin]=hist(dtser(:),100);      % Val, Bin
line(Val*len/max(Val),Bin,'linew',3,'colo','k')
axis tight
tmu = mean(dtser(:));
tsi = std (dtser(:));
legend(['\mu= ' num2str(tmu*1e3) ', \sigma= ' num2str(tsi*1e3) ' [ms]'])
emarkpl( filenam);

means = zeros(3,n);
sigms = means;
s1=24;              % 1st sample
if findstr(filenam,'7775028')
    s1 = 10;
end;
xyzCln = cell(3,1); % each cell is X, Y or Z, no transient
linecol='bgrcmyk';          % colors order
for d = 1:3        % dim
    xyzCln{d} = squeeze(xyz(s1:len,d,:));      %176xn samples, dimension d, all n freqs
    mn = min(xyzCln{d}(:)); mx = max(xyzCln{d}(:));
    aa = [0 1 mn mx];   % for all subplot axes
    figure
    for i=1:n
        subplot(prows, pcols,i)
        plot(t(s1:len,i),xyzCln{d}(:,i),'.-')
        axis(aa); grid on
        title([ xyzstr(d) '-meas. for freq = ' num2str(f_in(i)) ' Hz'])
        %sqd = squeeze(txyz(s1:len,d+1,i));
        means(d,i) = mean(xyzCln{d}(:,i));
        sigms(d,i) = std (xyzCln{d}(:,i));
        xlabel(['[\mu \sigma]=[' num2str([means(d,i) sigms(d,i)]) ']'])
    end
    if pcols>1
        set(gcf,'posi',[80 80 750 1000])
    end;
    emarkpl( filenam)
end;

subplot(211)
%semilogx(f_in, means','.-')
plot(f_in, means','.-')
title('mean')
legend('X', 'Y', 'Z',0)
axis tight, grid on
subplot(212)
%semilogx(f_in, sigms','.-')
plot(f_in, sigms','.-')
title('sigma')
legend('X', 'Y', 'Z',0)
xlabel ('disturb. freq')
axis tight, grid on
emarkpl( filenam)

figure
for d=1:3
    subplot(3,1,d)
    stem3(f_in,t(s1:end,1), xyzCln{d},'.')
    view(-26,34)
    mesh(audiof,t(s1:end,1), squeeze(xyz(s1:end,d,:)))
    contour3(audiof,t(s1:end,1), squeeze(xyz(s1:end,d,:)))
    xlabel('freq[Hz]'), ylabel('time[sec]'), title(xyzstr(d))
    axis tight
end;

if pcols>1,
    set(gcf,'posi',[80 80 750 1000])
end
emarkpl( filenam);

Hs=spectrum.welch;
% EstimationMethod: 'Welch'
%    SegmentLength: 64
% OverlapPercent: 50
%     WindowName: 'Hamming'
% SamplingFlag: 'symmetric'
% psd(Hs,,'Fs',Fs)

axnum = zeros(3,1);
for ax = 1:3,        % dim x, y or z
    axnum(ax) = figure    ;
    for i=1:n
        subplot(prows, pcols,i)
        psd(Hs,xyzCln{d}(:,i),'Fs',Fs);
        %  set(get(gca,'chi'),'colo','m'),
        axis tight,
        aa=axis;
        ine (rem(f_in(i),Fs/2)*[1;1],aa(3:4),'color', 'r')
        he=legend('PSD',['Noise freq=' num2str(f_in(i)) ' Hz'],0);
        set(he,'fontsi',8, 'box', 'off');
    end;
    xlabel(['Frequency (Hz) Axis: ' xyzstr(ax)]);
end;
for ax = 1:3,
    figure(axnum(ax))
    if pcols>1,
        set(gcf,'posi',[80 80 750 1000])
    end
    emarkpl( filenam)
end;

Hs.SegmentLength = len-s1+1; % 176

axnum = zeros(3,1);
for ax = 1:3        % dim x, y or z
    % t_frqin = squeeze(txyz(s1:len,ax+1,:)); dimentions x freq.in
    axnum(ax) = figure    ;
    for i=1:n
        subplot(prows, pcols,i)
        psd(Hs,xyzCln{d}(:,i),'Fs',Fs);
        axis tight, aa=axis;
        line (rem(f_in(i),Fs/2)*[1;1],aa(3:4),'color', 'r')
        he=legend('PSD',['Noise freq=' num2str(f_in(i)) ' Hz']);
        set(he,'fontsi',8, 'box', 'off');
    end;
    xlabel(['Frequency (Hz) Axis: ' xyzstr(ax)]);
end;
for ax = 1:3
    figure(axnum(ax))
    if pcols>1,
        set(gcf,'posi',[80 80 750 1000])
    end
    emarkpl( filenam)
end;

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone