text:
function metric = metric_euclidean(varargin)
%METRIC_EUCLIDEAN An euclidean metric function
%
%  Description
%    METRIC = METRIC_EUCLIDEAN('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a an euclidean metric function structure in which the
%    named parameters have the specified values. Either
%    'components' or 'deltadist' has to be specified. Any
%    unspecified parameters are set to default values.
%
%    METRIC = METRIC_EUCLIDEAN(METRIC,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a metric function structure with the named parameters
%    altered with the specified values.
%
%    Parameters for Euclidean metric function [default]
%       components        - cell array of vectors specifying which
%                           inputs are grouped together with a same
%                           scaling parameter. For example, the
%                           component specification {[1 2] [3]}
%                           means that distance between 3
%                           dimensional vectors computed as
%                           r = (r_1^2 + r_2^2 )/l_1 + r_3^2/l_2,
%                           where r_i are distance along component
%                           i, and l_1 and l_2 are lengthscales for
%                           corresponding component sets. If
%                           'components' is not specified, but
%                           'deltadist' is specified, then default
%                           is {1 ... length(deltadist)}
%       deltadist         - indicator vector telling which component sets
%                           are handled using the delta distance
%                           (0 if x=x', and 1 otherwise). Default is
%                           false for all component sets.
%       lengthScale       - lengthscales for each input component set
%                           Default is 1 for each set
%       lengthScale_prior - prior for lengthScales [prior_unif]
%
%  See also
%    GP_SET, GPCF_SEXP
%
% Copyright (c) 2008 Jouni Hartikainen
% Copyright (c) 2008 Jarno Vanhatalo
% Copyright (c) 2010 Aki Vehtari

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

  ip=inputParser;
  ip.FunctionName = 'METRIC_EUCLIDEAN';
  ip.addOptional('metric', [], @isstruct);
  ip.addParamValue('components',[], @(x) isempty(x) || iscell(x));
  ip.addParamValue('deltadist',[], @(x) isvector(x) || isempty(x));
  ip.addParamValue('lengthScale',[], @(x) isvector(x) && all(x>0));
  ip.addParamValue('lengthScale_prior',prior_unif, ...
                   @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  metric=ip.Results.metric;

  if isempty(metric)
    % Initialize a Gaussian process
    init=true;
  else
    % Modify a Gaussian process
    if ~isfield(metric,'type') && isequal(metric.type,'metric_euclidean')
      error('First argument does not seem to be a metric structure')
    end
    init=false;
  end

  if init
    % Type
    metric.type = 'metric_euclidean';
  end

  % Components
  if init || ~ismember('components',ip.UsingDefaults)
    metric.components = ip.Results.components;
  end
  % Deltadist
  if init || ~ismember('deltadist',ip.UsingDefaults)
    metric.deltadist = ip.Results.deltadist;
  end
  % Components+Deltadist check and defaults
  if isempty(metric.components) && isempty(metric.deltadist)
    error('Either ''components'' or ''deltadist'' has to be specified')
  elseif isempty(metric.components)
    metric.components=num2cell(1:length(metric.components));
  elseif isempty(metric.deltadist)
    metric.deltadist = false(1,length(metric.components));
  end
  % Lengthscale
  if init || ~ismember('lengthScale',ip.UsingDefaults)
    metric.lengthScale = ip.Results.lengthScale;
    if isempty(metric.lengthScale)
      metric.lengthScale = repmat(1,1,length(metric.components));
    end
  end
  % Prior for lengthscale
  if init || ~ismember('lengthScale_prior',ip.UsingDefaults)
    metric.p=[];
    metric.p.lengthScale = ip.Results.lengthScale_prior;
  end

  if init
    % Set the function handles to the subfunctions
    metric.fh.pak       = @metric_euclidean_pak;
    metric.fh.unpak     = @metric_euclidean_unpak;
    metric.fh.lp        = @metric_euclidean_lp;
    metric.fh.lpg       = @metric_euclidean_lpg;
    metric.fh.dist      = @metric_euclidean_dist;
    metric.fh.distg     = @metric_euclidean_distg;
    metric.fh.ginput    = @metric_euclidean_ginput;
    metric.fh.recappend = @metric_euclidean_recappend;
  end

end

function [w s] = metric_euclidean_pak(metric)
%METRIC_EUCLIDEAN_PAK  Combine GP covariance function
%                      parameters into one vector.
%
%  Description
%    W = METRIC_EUCLIDEAN_PAK(GPCF) takes a covariance function
%    structure GPCF and combines the covariance function
%    parameters and their hyperparameters into a single row
%    vector W and takes a logarithm of the covariance function
%    parameters.
%
%       w = [ log(metric.lengthScale(:))
%             (hyperparameters of metric.lengthScale)]'
%
%  See also
%    METRIC_EUCLIDEAN_UNPAK

  w = []; s = {};
  if ~isempty(metric.p.lengthScale)
    w = log(metric.lengthScale);
    if numel(metric.lengthScale)>1
      s = [s; sprintf('log(metric.lengthScale x %d)',numel(metric.lengthScale))];
    else
      s = [s; 'log(metric.lengthScale)'];
    end
    % Hyperparameters of lengthScale
    [wh sh] = metric.p.lengthScale.fh.pak(metric.p.lengthScale);
    w = [w wh];
    s = [s; sh];
  end

end

function [metric, w] = metric_euclidean_unpak(metric, w)
%METRIC_EUCLIDEAN_UNPAK  Separate metric parameter vector into components
%
%  Description
%    METRIC, W] = METRIC_EUCLIDEAN_UNPAK(METRIC, W) takes a
%    metric structure GPCF and a parameter vector W, and returns
%    a covariance function structure identical to the input,
%    except that the covariance parameters have been set to the
%    values in W. Deletes the values set to GPCF from W and
%    returns the modified W.
%
%    The covariance function parameters are transformed via exp
%    before setting them into the structure.
%
%  See also
%    METRIC_EUCLIDEAN_PAK
%

  if ~isempty(metric.p.lengthScale)
    i2=length(metric.lengthScale);
    i1=1;
    metric.lengthScale = exp(w(i1:i2));
    w = w(i2+1:end);

    % Hyperparameters of lengthScale
    [p, w] = metric.p.lengthScale.fh.unpak(metric.p.lengthScale, w);
    metric.p.lengthScale = p;
  end
end

function lp = metric_euclidean_lp(metric)
%METRIC_EUCLIDEAN_LP  Evaluate the log prior of metric parameters
%
%  Description
%    LP = METRIC_EUCLIDEAN_LP(METRIC) takes a metric structure
%    METRIC and returns log(p(th)), where th collects the
%    parameters.
%
%  See also
%    METRIC_EUCLIDEAN_PAK, METRIC_EUCLIDEAN_UNPAK, METRIC_EUCLIDEAN_G, GP_E
%

% Evaluate the prior contribution to the error. The parameters that
% are sampled are from space W = log(w) where w is all the "real" samples.
% On the other hand errors are evaluated in the W-space so we need take
% into account also the  Jakobian of transformation W -> w = exp(W).
% See Gelman et al. (2013), Bayesian Data Analysis, third edition, p. 21.
  if ~isempty(metric.p.lengthScale)
    lp = metric.p.lengthScale.fh.lp(metric.lengthScale, metric.p.lengthScale) + sum(log(metric.lengthScale));
  else
    lp=0;
  end

end

function lpg = metric_euclidean_lpg(metric)
%METRIC_EUCLIDEAN_LPG  d log(prior)/dth of the metric parameters th
%
%  Description
%    LPG = METRIC_EUCLIDEAN_LPG(METRIC) takes a likelihood
%    structure METRIC and returns d log(p(th))/dth, where th
%    collects the parameters.
%
%  See also
%    METRIC_EUCLIDEAN_PAK, METRIC_EUCLIDEAN_UNPAK, METRIC_EUCLIDEAN, GP_E
%

% Evaluate the prior contribution of gradient with respect to lengthScale
  if ~isempty(metric.p.lengthScale)
    i1=1;
    lll = length(metric.lengthScale);
    lpgs = metric.p.lengthScale.fh.lpg(metric.lengthScale, metric.p.lengthScale);
    lpg(i1:i1-1+lll) = lpgs(1:lll).*metric.lengthScale + 1;
    lpg = [lpg lpgs(lll+1:end)];
  end
end

function gdist  = metric_euclidean_distg(metric, x, x2, mask)
%METRIC_EUCLIDEAN_DISTG  Evaluate the gradient of the metric function
%
%  Description
%    DISTG = METRIC_EUCLIDEAN_DISTG(METRIC, X) takes a metric
%    structure METRIC together with a matrix X of input
%    vectors and return the gradient matrices GDIST and
%    GPRIOR_DIST for each parameter.
%
%    DISTG = METRIC_EUCLIDEAN_DISTG(METRIC, X, X2) forms the
%    gradient matrices between two input vectors X and X2.
%
%    DISTG = METRIC_EUCLIDEAN_DISTG(METRIC, X, X2, MASK) forms
%    the gradients for masked covariances matrices used in sparse
%    approximations.
%
%  See also
%    METRIC_EUCLIDEAN_PAK, METRIC_EUCLIDEAN_UNPAK, METRIC_EUCLIDEAN, GP_E
%

  gdist=[];
  components = metric.components;

  n = size(x,1);
  m = length(components);
  i1=0;i2=1;

  % NOTE! Here we have already taken into account that the parameters
  % are transformed through log() and thus dK/dlog(p) = p * dK/dp

  if ~isempty(metric.p.lengthScale)
    if nargin <= 3
      if nargin == 2
        x2 = x;
      end
      ii1=0;

      dist  =  0;
      distc = cell(1,m);
      % Compute the distances for each component set
      for i=1:m
        if length(metric.lengthScale)==1
          s=1./metric.lengthScale.^2;
        else
          s=1./metric.lengthScale(i).^2;
        end
        distc{i} = 0;
        for j = 1:length(components{i})
          if metric.deltadist(i)
            distc{i} = distc{i} + double(bsxfun(@ne,x(:,components{i}(j)),x2(:,components{i}(j))'));
          else
            distc{i} = distc{i} + bsxfun(@minus,x(:,components{i}(j)),x2(:,components{i}(j))').^2;
          end
        end
        distc{i} = distc{i}.*s;
        % Accumulate to the total distance
        dist = dist + distc{i};
      end
      dist = sqrt(dist);
      % Loop through component sets
      if length(metric.lengthScale)==1
        D = -distc{1};
        D(dist~=0) = D(dist~=0)./dist(dist~=0);
        ii1 = ii1+1;
        gdist{ii1} = D;
      else
        for i=1:m
          D = -distc{i};
          ind = dist~=0;
          D(ind) = D(ind)./dist(ind);
          ii1 = ii1+1;
          gdist{ii1} = D;
        end
      end
% $$$         elseif nargin == 3
% $$$             if size(x,2) ~= size(x2,2)
% $$$                 error('metric_euclidean -> _ghyper: The number of columns in x and x2 has to be the same. ')
% $$$             end
    elseif nargin == 4
      gdist = cell(1,length(metric.lengthScale));
    end

    % Evaluate the prior contribution of gradient with respect to lengthScale
    if ~isempty(metric.p.lengthScale)
      i1=1;
      lll = length(metric.lengthScale);
      gg = -metric.p.lengthScale.fh.lpg(metric.lengthScale, metric.p.lengthScale);
      gprior(i1:i1-1+lll) = gg(1:lll).*metric.lengthScale - 1;
      gprior = [gprior gg(lll+1:end)];
    end
  end
end

function dist = metric_euclidean_dist(metric, x1, x2)
%METRIC_EUCLIDEAN_DIST  Compute the euclidean distence between
%                       one or two matrices.
%
%  Description
%    DIST = METRIC_EUCLIDEAN_DIST(METRIC, X) takes a metric
%    structure METRIC together with a matrix X of input
%    vectors and calculates the euclidean distance matrix DIST.
%
%    DIST = METRIC_EUCLIDEAN_DIST(METRIC, X1, X2) takes a
%    metric structure METRIC together with a matrices X1 and
%    X2 of input vectors and calculates the euclidean distance
%    matrix DIST.
%
%  See also
%    METRIC_EUCLIDEAN_PAK, METRIC_EUCLIDEAN_UNPAK, METRIC_EUCLIDEAN, GP_E
%
  if (nargin == 2 || isempty(x2))
    % use fast c-code for self-distance
    x2=x1;
    % force deltadist to be logical for simplified c-code
    metric.deltadist=logical(metric.deltadist);
    dist = dist_euclidean(metric,x1);
    if ~any(isnan(dist))
      % if c-code was available, result is not NaN
      return
    end
  end

  [n1,m1]=size(x1);
  [n2,m2]=size(x2);

  if m1~=m2
    error('the number of columns of X1 and X2 has to be same')
  end

  components = metric.components;
  m = length(components);
  dist  =  0;

  s=1./metric.lengthScale.^2;
  if m>numel(s)
    s=repmat(s,1,m);
  end
  for i=1:m
    for j = 1:length(components{i})
      if metric.deltadist(i)
        dist = dist + s(i).*double(bsxfun(@ne,x1(:,components{i}(j)),x2(:,components{i}(j))'));
      else
        dist = dist + s(i).*bsxfun(@minus,x1(:,components{i}(j)),x2(:,components{i}(j))').^2;
      end
    end
  end
  dist=sqrt(dist); % euclidean distance

end

function [ginput, gprior_input]  = metric_euclidean_ginput(metric, x1, x2)
%METRIC_EUCLIDEAN_GINPUT  Compute the gradient of the
%                         euclidean distance function with
%                         respect to input. [n, m]=size(x);
  ii1 = 0;
  components = metric.components;

  if nargin == 2 || isempty(x2)
    x2=x1;
  end

  [n1,m1]=size(x1);
  [n2,m2]=size(x2);

  if m1~=m2
    error('the number of columns of X1 and X2 has to be same')
  end

  s = 1./metric.lengthScale.^2;
  dist = 0;
  for i=1:length(components)
    for j = 1:length(components{i})
      if metric.deltadist(i)
        dist = dist + s(i).*double(bsxfun(@ne,x1(:,components{i}(j)),x2(:,components{i}(j))'));
      else
        dist = dist + s(i).*bsxfun(@minus,x1(:,components{i}(j)),x2(:,components{i}(j))').^2;
      end
    end
  end
  dist = sqrt(dist);

  for i=1:m1
    for j = 1:n1
      DK = zeros(n1,n2);
      for k = 1:length(components)
        if ismember(i,components{k})
          if metric.deltadist(i)
            DK(j,:) = DK(j,:)+s(k).*double(bsxfun(@ne,x1(j,i),x2(:,i)'));
          else
            DK(j,:) = DK(j,:)+s(k).*bsxfun(@minus,x1(j,i),x2(:,i)');
          end
        end
      end
      if nargin == 2
        DK = DK + DK';
      end
      DK(dist~=0) = DK(dist~=0)./dist(dist~=0);

      ii1 = ii1 + 1;
      ginput{ii1} = DK;
      gprior_input(ii1) = 0;
    end
  end
  %size(ginput)
  %ginput

end

function recmetric = metric_euclidean_recappend(recmetric, ri, metric)
%RECAPPEND  Record append
%
%  Description
%    RECMETRIC = METRIC_EUCLIDEAN_RECAPPEND(RECMETRIC, RI,
%    METRIC) takes old metric function record RECMETRIC, record
%    index RI and metric function structure. Appends the
%    parameters of METRIC to the RECMETRIC in the ri'th place.
%
%  See also
%    GP_MC and GP_MC -> RECAPPEND

% Initialize record
  if nargin == 2
    recmetric.type = 'metric_euclidean';
    metric.components = recmetric.components;

    % Initialize parameters
    recmetric.lengthScale = [];

    % Set the function handles
    recmetric.fh.pak       = @metric_euclidean_pak;
    recmetric.fh.unpak     = @metric_euclidean_unpak;
    recmetric.fh.lp        = @metric_euclidean_lp;
    recmetric.fh.lpg       = @metric_euclidean_lpg;
    recmetric.fh.dist      = @metric_euclidean_dist;
    recmetric.fh.distg     = @metric_euclidean_distg;
    recmetric.fh.ginput    = @metric_euclidean_ginput;
    recmetric.fh.recappend = @metric_euclidean_recappend;
    return
  end
  mp = metric.p;

  % record parameters
  if ~isempty(metric.lengthScale)
    recmetric.lengthScale(ri,:)=metric.lengthScale;
    recmetric.p.lengthScale = metric.p.lengthScale.fh.recappend(recmetric.p.lengthScale, ri, metric.p.lengthScale);
  elseif ri==1
    recmetric.lengthScale=[];
  end

end

author: gpstuff-dev

repo: gpstuff

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


text:
function triangle_ncc_rule_test02 ( )

%*****************************************************************************80
%
%% TEST02 tests TRIANGLE_NCC_RULE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  TRIANGLE_NCC_RULE returns the points and weights\n' );
  fprintf ( 1, '  of an NCC rule for the triangle.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this test, we simply check that the weights\n' );
  fprintf ( 1, '  sum to 1.\n' );

  rule_num = triangle_ncc_rule_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of available rules = %d\n', rule_num );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      Rule    Sum of weights\n' );
  fprintf ( 1, '\n' );

  for rule = 1 : rule_num

    order_num = triangle_ncc_order_num ( rule );

    [ xy, w ] = triangle_ncc_rule ( rule, order_num );

    w_sum = sum ( w(1:order_num) );

    fprintf ( 1, '  %8d  %14f\n', rule, w_sum );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef SwerlingIISqLawD
%%SWERLINGIISQLAWD This class implements various functions and detections
%                statistics related to the Swerling II target model when
%                used with a square law detector. All methods are static
%                (the class holds no information and does not need to be
%                instantiated).
%
%Implemented methods are: mean, var, PDF, CDF, PD4Threshold, PD4PFA,
%                         avgSNR4PDThresh, avgSNR4PDPFA, thresh4AvgSNRPD,
%                         rand
%
%Swerling models are given in [1] and in Chapter 11 of [2]. The Swerling II
%model assumes that the observed power SNR of the target varies in terms of
%an exponential distribution with rate parameter (1/avgSNR), that in a
%pulse train for detection the target SNR does fluctuate from pulse to
%pulse, and that a square law detector is used to incoherently integrate
%multiple pulses. The difference from a Swerling I model is the target
%amplitude fluctuates from pulse to pulse.
%
%The square law detector for N samples is
%y=sum_{i=1}^N r_i^2
%where r_i is the ith real amplitude. The model for the amplitude given a
%target amplitude or power is developed in Chapter 9 of [2]. The squared
%real amplitude is r^2_i=y_{I,i}^2+y_{Q,i}^2 where y_{I,i} and y_{Q,i} are
%the in-phase and quadrature components of the filter output. Define
%Rp=2*sampPowSNR, as in Equation 9.2-34. The sampled power SNR of the
%target in the Swerling II model is a different random value for each
%pulse. The sample power SNR is sampPowSNR=ExponentialD.rand(1,1/avgSNR)
%as mentioned in [1]. The value y=y_{I,i}+sqrt(-1)*y_{Q,i} conditioned on
%the sampled power SNR is modeled as being distributed circularly complex
%Gaussian with mean sqrt(Rp) and variance 2. (See, Equation 9.3-35a in
%[2]). The variance being 2 simply reflects having normalized the variance
%on the I and Q components each to 1.  As in Equation 10.4-13 of [2], the
%distribution of Y=y/2 conditioned on the target SNR values is noncentral
%chi-squared with a change of variables. The value y conditioned on the
%target SNR values is noncentral chi-squared with nu=2*N degrees of freedom
%and lambda=N*Rp as the noncentrality parameter. The expressions in [1] and
%[2] are in terms of Y and not y, so a transformation has to be performed
%to make things in terms of y.
%
%Under this model, for a given average power SNR value, one can generate a
%random sample from the distribution as
% sample=0;
% for curPulse=1:N
%     sampPowSNR=ExponentialD.rand(1,1/avgSNR);
%     Rp=2*sampPowSNR;%The peak power of the signal from the SNR, Equation
%                 %9.4-8 in [2].
%     A=sqrt(Rp);%Equation 9.4-8, the amplitude of the signal
%     sample=sample+sum(abs(ComplexGaussianD.rand(1,A,2)).^2);
% end
%
%REFERENCES:
%[1] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%[2] J. V. Di Franco and W. L. Rubin, Radar Detection. SciTech Publishing
%    Inc., Rayliegh, NC: 2004.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)
function val=mean(avgSNR,N)
%%MEAN Obtain the mean of the distribution of the detection power
%      (normalized to a noise variance of 1) under a Swerling II model in a
%      square-law detector.
%
%INPUTS:avgSNR The average power signal to noise ratio of the target.
%        N The number of pulses that are to be incoherently added for
%          detection (in a square-law detector). In a Swerling II model, a
%          different realization of the target power is used across each
%          pulse and the noise corrupting the pulses varies. If this
%          parameter is omitted or an empty matrix is passed, N=1 is used.
%
%OUTPUTS: val The value of the mean for the given parameters.
%
%The expression for the mean implemented here is E{x}, where the expected
%value is taken over the PDF in Equation III.10 of [1]. However, the first
%moment from that expression has to be multiplied by 2, because as seen in
%Equation 10.4-4 of [2], the definition in [2], which is the same as in
%[1], is in terms of a scaled version of the square law detector.
%
%EXAMPLE:
%Here, we validate the mean by generating random samples and comparing the
%computed mean to the sample mean.
% avgSNR=2;
% N=4;
% numSamples=1e5;
% meanVal=SwerlingIISqLawD.mean(avgSNR,N)
% meanSampVal=mean(SwerlingIISqLawD.rand([numSamples,1],avgSNR,N))
%One will see that both mean values are about 24.
%
%REFERENCES:
%[1] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%[2] J. V. Di Franco and W. L. Rubin, Radar Detection. SciTech Publishing
%    Inc., Rayliegh, NC: 2004.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

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

    val=2*N.*(1+avgSNR);
end

function val=var(avgSNR,N)
%%VAR Obtain the variance of the distribution of the detection power
%     (normalized to a noise variance of 1) under a Swerling II model in a
%     square-law detector.
%
%INPUTS:avgSNR The average power signal to noise ratio of the target.
%        N The number of pulses that are to be incoherently added for
%          detection (in a square-law detector). In a Swerling II model, a
%          different realization of the target power is used across each
%          pulse and the noise corrupting the pulses varies. If this
%          parameter is omitted or an empty matrix is passed, N=1 is used.
%
%OUTPUTS: val The value of the variance for the given parameters.
%
%To derive the formula implemented here, the first and second noncentral
%moments were obtained by finding E{x^2} and E{x} integrating over the PDF
%in Equation III.10 in [1].  One then uses the identity that
%variance=E{x^2}-E{x}^2
%However, the variance from that expression has to be multiplied by 4,
%because as seen in Equation 10.4-4 of [2], the definition in [2], which is
%the same as the definition in [1], is in terms of a scaled version of the
%square law detector.
%
%EXAMPLE:
%Here, we validate the variance by generating random samples and comparing
%the computed variance to the sample variance.
% avgSNR=2;
% N=4;
% numSamples=1e5;
% varVal=SwerlingIISqLawD.var(avgSNR,N)
% varSampVal=var(SwerlingIISqLawD.rand([numSamples,1],avgSNR,N))
%One will see that both variance values are about 144.
%
%REFERENCES:
%[1] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%[2] J. V. Di Franco and W. L. Rubin, Radar Detection. SciTech Publishing
%    Inc., Rayliegh, NC: 2004.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

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

    val=4*N*(1+avgSNR)^2;
end

function val=PDF(v,avgSNR,N)
%%PDF Evaluate the scalar probability density function (PDF) of the
%     distribution of the detection power (normalized to a noise variance
%     of 1) under a Swerling II model in a square-law detector.
%
%INPUTS: v The point or points at which the PDF should be evaluated. v>0
%          for nonzero PDF values.
%   avgSNR The average power signal to noise ratio of the target.
%        N The number of pulses that are to be incoherently added for
%          detection (in a square-law detector). In a Swerling II model, a
%          different realization of the target power is used across each
%          pulse and the noise corrupting the pulses varies. If this
%          parameter is omitted or an empty matrix is passed, N=1 is used.
%
%OUTPUTS: val The values of the PDF evaluated at the given points.
%
%The PDF is taken from Equation III.10 of [1]. However the PDF is derived
%based on a scaled square-law detector, as in Equation 10.4-4 in [2], which
%uses the same definition as in [1]. Thus, in implementing the PDF, the
%scaling is removed.
%
%EXAMPLE:
%Here, we validate the PDF by generating random samples and comparing the
%PDF plot with a histogram of the random samples.
% avgSNR=2;
% N=4;
% numSamples=10e3;
%
% figure(1)
% clf
% histogram(SwerlingIISqLawD.rand([numSamples,1],avgSNR,N),'Normalization','pdf')
% hold on
% numPoints=1000;
% x=linspace(0,80,numPoints);
% vals=SwerlingIISqLawD.PDF(x,avgSNR,N);
% plot(x,vals,'linewidth',2)
% h1=xlabel('x');
% h2=ylabel('PDF(x)');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%One will see that the histogram matches well with the plot.
%
%REFERENCES:
%[1] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%[2] J. V. Di Franco and W. L. Rubin, Radar Detection. SciTech Publishing
%    Inc., Rayliegh, NC: 2004.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    %Initial change of variable to undo the scaling from Equation 10.4-4 in
    %[2].
    v=v/2;

    val=exp((N-1)*log(v)-N*log(1+avgSNR)-gammaln(N)-v/(1+avgSNR));

    %Adjust for the change of variable to undo the scaling from Equation
    %10.4-4 in [2].
    val=val/2;

    val(v<0)=0;
end

function val=CDF(v,avgSNR,N)
%%CDF Evaluate the scalar cumulative distribution function (CDF) of the
%     distribution of the detection power (normalized to a noise variance
%     of 1) under a Swerling II model in a square-law detector.
%
%INPUTS: v The point or points at which the CDF should be evaluated. v>0
%          for nonzero CDF values.
%   avgSNR The average power signal to noise ratio of the target.
%        N The number of pulses that are to be incoherently added for
%          detection (in a square-law detector). In a Swerling II model, a
%          different realization of the target power is used across each
%          pulse and the noise corrupting the pulses varies. If this
%          parameter is omitted or an empty matrix is passed, N=1 is used.
%
%OUTPUTS: val The values of the CDF evaluated at the given points.
%
%The derivation of the detection probability in [1] is as 1-the value of
%the CDF. Thus, this function just evaluated
%1-SwerlingIISqLawD.PD4Threshold(avgSNR,v,N).
%
%REFERENCES:
%[1] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(N))
        N=1;
    end

    val=1-SwerlingIISqLawD.PD4Threshold(avgSNR,v,N);
    val(v<=0)=0;
end

function PD=PD4Threshold(avgSNR,thresh,N)
%%PD4THRESHOLD Determine the detection probability of a Swerling II target
%           given its signal to noise ratio and the value of the
%           detection threshold, assuming that a square-law detector is
%           used. The noise in each sample is assumed to have variance 1
%           (normalized noise power).
%
%INPUTS: avgSNR A vector or matrix of average power signal to noise ratios
%               of the target at which one wishes to evaluate the detection
%               probability.
%        thresh The scalar normalized detection threshold to use. This is
%               the threshold to use if the noise variance is 1.
%             N The number of pulses that are to be incoherently added for
%               detection (in a square-law detector). In a Swerling II
%               model, a different realization of the target power is used
%               across each pulse and the noise corrupting the pulses
%               varies. If this parameter is omitted or an empty matrix is
%               passed, N=1 is used.
%
%OUTPUTS: PD The detection probability of the target.
%
%This function implements Equation II.10 in [1]. However, the expression
%used is derived based on a scaled square-law detector, as in Equation
%10.4-4 in [2], which uses the same definition as in [1]. Thus, the
%threshold in this function is scaled appropriately.
%
%EXAMPLE:
%Here, we validate the results by comparing the PD from this function to
%the PD computed using random samples.
% avgSNR=2;
% thresh=30;
% N=4;
% PD=SwerlingIISqLawD.PD4Threshold(avgSNR,thresh,N)
% numSamples=1e5;
% PDSamp=mean(SwerlingIISqLawD.rand([numSamples,1],avgSNR,N)>=thresh)
%One will see that both PD values are about 0.265.
%
%REFERENCES:
%[1] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%[2] J. V. Di Franco and W. L. Rubin, Radar Detection. SciTech Publishing
%    Inc., Rayliegh, NC: 2004.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(N))
       N=1;
    end
    thresh=thresh/2;

    PD=1-PearsonsGammaInc(thresh./((1+avgSNR)*sqrt(N)),N-1);
end

function PD=PD4PFA(avgSNR,PFA,N)
%%PD4PFA Determine the detection probability of a Swerling II target given
%        its signal to noise ratio and the probability of false alarm
%        implied by a particular unspecified detection threshold, assuming
%        that a square law detector is used.
%
%INPUTS: avgSNR A vector or matrix of average power signal to noise ratios
%               of the target at which one wishes to evaluate the detection
%               probability.
%           PFA The probability of false alarm, 0<=PFA<1.
%             N The number of pulses that are to be incoherently added for
%               detection (in a square-law detector). In a Swerling II
%               model, a different realization of the target power is used
%               across each pulse and the noise corrupting the pulses
%               varies. If this parameter is omitted or an empty matrix is
%               passed, N=1 is used.
%
%OUTPUTS: PD The detection probability of the target, 0<=PD<=1.
%
%This function just calls PFA2SquareLawThreshold to convert the false alarm
%rate into a normalized detection threshold after which it calls
%SwerlingIISqLawD.PD4Threshold.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(N))
       N=1;
    end

    thresh=PFA2SquareLawThreshold(PFA,N);
    PD=SwerlingIISqLawD.PD4Threshold(avgSNR,thresh,N,method);
end

function avgSNR=avgSNR4PDThresh(PD,thresh,N)
%%AVGSNR4PDTHRESH Given a detection probabilty and a normalized detection
%       threshold (normalized in terms of the receiver noise having a unit
%       covariance), determine the average power signal to noise ratio
%       (SNR) of the target needed under a Swerling II model for one or
%       more pulses in a square-law detector.
%
%INPUTS: PD The detection probability of the target , 0<=PD<1.
%    thresh The scalar normalized detection threshold to use. This is the
%           threshold to use if the noise variance is 1.
%         N The number of pulses that are to be incoherently added for
%           detection (in a square-law detector). In a Swerling II model, a
%           different realization of the target power is used across each
%           pulse and the noise corrupting the pulses varies. If this
%           parameter is omitted or an empty matrix is passed, N=1 is used.
%
%OUTPUTS: avgSNR The average power signal to noise ratio needed to achieve
%                the desired probability of detection at the given
%                threshold. If this is negative, then for a given PD, the
%                threshold is so low that due to the high false alarm
%                rate, the desired PD is impossiblely low.
%
%Equation II.10 in [1] gives an expression for PD in terms of thresh and
%the average SNR. Here, we have simply inverted the expression.
%
%However, the expressions used are derived based on a scaled square-law
%detector, as in Equation 10.4-4 in [2], which is the same definition used
%in [1]. Thus, the threshold in this function is scaled appropriately.

%EXAMPLE:
%Here, we show that the results are consistent.
% N=4;
% PD=0.5;
% thresh=5;
% avgSNR=SwerlingIISqLawD.avgSNR4PDThresh(PD,thresh,N);
% PDBack=SwerlingIISqLawD.PD4Threshold(avgSNR,thresh,N)
% numSamples=1e5;
% PDSamp=mean(SwerlingIISqLawD.rand([numSamples,1],avgSNR,N)>=thresh)
%One will see that PDBack is the same as PD and about the same as PDSamp.
%
%REFERENCES:
%[1] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%[2] J. V. Di Franco and W. L. Rubin, Radar Detection. SciTech Publishing
%    Inc., Rayliegh, NC: 2004.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(N))
       N=1;
    end

    %The division by 2 deals with the scaling in Equation 10.4-4.
    avgSNR=(thresh/2)/gammaincinv(1-PD,N)-1;
end

function avgSNR=avgSNR4PDPFA(PD,PFA,N)
%%AVGSNR4PDTHRESH Given a detection probabilty and the probability of false
%       alarm, determine the average power signal to noise ratio (SNR) of
%       the target needed under a  Swerling II model for one or more pulses
%       in a square-law detector.
%
%INPUTS: PD The detection probability of the target , 0<=PD<1.
%       PFA The probability of false alarm, 0<=PFA<1.
%         N The number of pulses that are to be incoherently added for
%           detection (in a square-law detector). In a Swerling II model, a
%           different realization of the target power is used across each
%           pulse and the noise corrupting the pulses varies. If this
%           parameter is omitted or an empty matrix is passed, N=1 is used.
%
%OUTPUTS: avgSNR The average power signal to noise ratio needed to achieve
%                the desired probability of detection at the given
%                threshold. If this is negative, then for a given PD, the
%                threshold is so low that due to the high false alarm
%                rate, the desired PD is impossiblely low.
%
%Equation II.10 in [1] gives an expression for PD in terms of a normalized
%threshold and the average SNR. Here, we call the function
%PFA2SquareLawThreshold and insert the result into
%SwerlingIISqLawD.avgSNR4PDThresh, which inverts the expression in [1].
%
%REFERENCES:
%[1] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(N))
        N=1;
    end

    thresh=PFA2SquareLawThreshold(PFA,N);
    avgSNR=SwerlingIISqLawD.avgSNR4PDThresh(PD,thresh,N);
end

function thresh=thresh4AvgSNRPD(PD,avgSNR,N)
%%THRESH4AVGSNRPD Given a detection probability and the average power
%           signal to noise ratio (SNR), determine the threshold needed
%           under a Swerling II model for one or more pulses in a
%           square-law detector.
%
%INPUTS: PD The detection probability of the target , 0<=PD<1.
%    avgSNR A vector or matrix of average power signal to noise ratios
%           of the target at which one wishes to evaluate the threshold.
%         N The number of pulses that are to be incoherently added for
%           detection (in a square-law detector). In a Swerling II model, a
%           different realization of the target power is used across each
%           pulse and the noise corrupting the pulses varies. If this
%           parameter is omitted or an empty matrix is passed, N=1 is used.
%
%OUTPUTS: thresh The normnalized detection threshold(s) (assuming the noise
%                variance is 1).
%
%Equation II.10 in [1] gives an expression for PD in terms of a normalized
%threshold and the average SNR. Here, we simply invert the equation.
%
%However, the expressions used are derived based on a scaled square-law
%detector, as in Equation 10.4-4 in [2], which is the same defintiion used
%in [1]. Thus, the threshold in this function is scaled appropriately.
%
%EXAMPLE:
%Here, we show that the results are consistent.
% N=4;
% PD=0.5;
% avgSNR=2;
% thresh=SwerlingIISqLawD.thresh4AvgSNRPD(PD,avgSNR,N);
% PDBack=SwerlingIISqLawD.PD4Threshold(avgSNR,thresh,N)
% numSamples=1e5;
% PDSamp=mean(SwerlingIISqLawD.rand([numSamples,1],avgSNR,N)>=thresh)
%One will see that PDBack=0.5 and PDSamp is about the same.
%
%REFERENCES:
%[1] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%[2] J. V. Di Franco and W. L. Rubin, Radar Detection. SciTech Publishing
%    Inc., Rayliegh, NC: 2004.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(N))
        N=1;
    end

    thresh=2*gammaincinv(1-PD,N).*(avgSNR+1);
end

function vals=rand(NDims,avgSNR,N)
%%RAND Generate random variables representing a target whose sampled signal
%      to average noise power ratio is generated according to a Swerling II
%      model detected by a square-law detector.
%
%INPUTS: NDims If NDims is a scalar, then rand returns an NDimsXNDims
%          matrix of random variables. If NDims=[M,N1] is a two-element row
%          vector, then rand returns an MXN1 matrix of random variables.
%   avgSNR The average power signal to noise ratio of the target.
%        N The number of pulses that are to be incoherently added for
%          detection (in a square-law detector). In a Swerling II model, a
%          different realization of the target power is used across each
%          pulse and the noise corrupting the pulses varies. If this
%          parameter is omitted or an empty matrix is passed, N=1 is used.
%
%OUTPUTS: vals A matrix whose dimensions are determined by N of the
%              generated random variables.
%
%Following the model of [1], a sample from the exponential distribution is
%performed to obtain the complex signal to noise ratio of every single one
%of the N pulses. Then, the square law detector output conditioned on the
%signal to noise ratio is generated directly using NonFlucSqLawD.rand for
%each pulse and the values are added.
%
%REFERENCES:
%[1] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(N))
        N=1;
    end

    if(isscalar(NDims))
        dims=[NDims, NDims];
    else
        dims=NDims;
    end

    vals=zeros(dims);

    numEls=numel(vals);
    for curEl=1:numEls
        for curPulse=1:N
            %Equation I.1 for  the input signal-to-noise power ratio. It is
            %different for each sample in the pulse train in the Swerling
            %II model.
            curSNR=ExponentialD.rand(1,1/avgSNR);
            vals(curEl)=vals(curEl)+NonFlucSqLawD.rand(1,curSNR,1);
        end
    end
end
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function DEM_demo_PEB
% DEM demo for a hierarchical linear model (MFX).  This inversion is
% cross-validated with restricted maximum likelihood using and parametric
% empirical Bayes (spm_PEB). It uses a simple two level model that embodies
% empirical shrinkage priors on the first-level parameters (c.f.,
% DEM_demo_GLM, with no empirical priors)

% MFX design
%==========================================================================
X1    = kron(eye(4),kron(ones(2,1),eye(2)));
X2    = kron(ones(4,1),eye(2));

% specify parameters
%--------------------------------------------------------------------------
N     = 1;                                      % length of data sequence
h     = {log(16) log(8)};                       % precisions

% generate data
%--------------------------------------------------------------------------
DEM   = spm_DEM_generate(spm_DEM_M('HLM',X1,X2),N,{},h);

% DEM estimation
%==========================================================================
DEM.M(1).E.nE = 16;

DEM   = spm_DEM(DEM);
qU    = DEM.qU;
qP    = DEM.qP;
qH    = DEM.qH;

% compare real and estimated factor and causes (using ReML and PEB)
%==========================================================================
P            = cell(1,2);
Q1           = DEM.M(1).Q{1};
Q2           = DEM.M(2).Q{1};
[Ce,hr,W,Fr] = spm_reml_sc(DEM.Y*DEM.Y',X1*X2,{Q1; X1*Q2*X1'},N);
P{1}.X       = X1;
P{1}.C       = {Q1};
P{2}.X       = X2;
P{2}.C       = {Q2};
[C,P,Fp]     = spm_PEB(DEM.Y,P);

% graphics for variance component estimators
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 1');

subplot(3,1,1)
qh   = spm_vec(qH.h);
bar([exp(-[h{1}; h{2}]) [exp(-qh)] [C{1}.h; C{2}.h]])
legend({'true' 'DEM' 'ReML'})
title({'ReML estimators';'Parametric Empirical Bayes'},'FontSize',16)
c     = sqrt(diag(qH.C))*1.64;
for i = 1:length(h)
    line([i i],exp(-([-c(i) c(i)] + qh(i))),'color',[1 0 0],'LineWidth',2)
end
axis square
grid on

% check estimators
%--------------------------------------------------------------------------
subplot(3,2,3)
bar(qU.v{2})
axis square; title('DEM level-1','FontSize',16)
a = axis;
subplot(3,2,4)
bar(C{2}.E)
axis square; title('VB level-1','FontSize',16)
axis(a)

subplot(3,2,5)
bar(qU.v{3})
axis square; title('DEM level-2','FontSize',16)
a = axis;
subplot(3,2,6)
bar(C{3}.E)
axis square; title('VB level-2','FontSize',16)
axis(a)

author: spm


text:
function [P,model] = chrrParseModel(model)
% Parse a COBRA model into the right format for the CHRR sampler
%
% USAGE:
%
%      [P,model] = chrrParseModel(model);
%
% We are trying to sample uniformly at random from the points v that satisfy:
%
% .. math::
%                     Sv = b\\
%             ~~ l_b \leq v \leq u_b
%
% INPUTS:
%    model:    COBRA model structure with fields:
%
%               * .S - The `m x n` stoichiometric matrix
%               * .lb - `n x 1` lower bounds on fluxes
%               * .ub - `n x 1` upper bounds on fluxes
%
% OPTIONAL INPUTS:
%               * .C - 'k x n' matrix of additional inequality constraints
%               * .d - 'k x 1' rhs of the above constraints
%               * .dsense - 'k x 1' the sense of the above constraints ('L' or 'G')
%
% OUTPUTS:
%    P:        A structure with fields:
%
%               * .A_eq - Equality constraint matrix (`model.S`)
%               * .b_eq - Right hand side of equality constraints (`model.b`)
%               * .A - Inequality constraint matrix (`[I_n 0; 0 -I_n]`)
%               * .b - Right hand side of inequality constraints (`[lb; -ub]`)
%
% .. Authors:
%       - Ben Cousins and Hulda Haraldsdóttir, 10/2016
%       - Ben Cousins, 12/2017, Moved objective function handling to preprocess function
%       - Ben Cousins, 05/2019, Added support for C,d inequalities.

dim = length(model.lb);

P.A = [eye(dim); -eye(dim)];
P.b = [model.ub; -model.lb];

if isfield(model,'C') && isfield(model,'d')
   for i=1:size(model.C,1)
      if model.dsense(i)=='G'
          % convert constraint to <=
          model.C(i,:) = model.C(i,:)*-1;
          model.d(i) = model.d(i)*-1;
      elseif model.dsense{i}=='E'
          error('Equality constraints not supported in C,d fields.');
      end
   end
   P.A = [P.A; model.C];
   P.b = [P.b; model.d];
end

P.A_eq = model.S;
P.b_eq = model.b;

end

author: opencobra

repo: cobratoolbox


text:
function pass = test_isreal()
% Test chebfun3/isreal.

f = chebfun3(@(x,y,z) x+y-z);
pass(1) = isreal(f);

f = chebfun3(@(x,y,z) 1i*x+y-z);
pass(2) = ~isreal(f);

f = real(chebfun3(@(x,y,z) 1i*x+y-z));
pass(3) = isreal(f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function i4mat_transpose_print ( m, n, a, title )

%*****************************************************************************80
%
%% I4MAT_TRANSPOSE_PRINT prints an I4MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 September 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer A(M,N), an M by N matrix to be printed.
%
%    Input, string TITLE, a title.
%
  i4mat_transpose_print_some ( m, n, a, 1, 1, m, n, title );

  return
end


text:
function [pck,nseg] = eval_plot_pck(isMatchAll, pidxs, verbose, pslidx)

if (nargin < 2)
    pidxs = 1:14;
end

if (nargin < 3)
    verbose = true;
end

if (nargin < 4)
    pslidx = -1;
end

nMatchAll = 0;
nSegAll = 0;
pck = zeros(1,15);
nseg = zeros(1,15);
idxsAll = [];
for pidx = pidxs
    idxs = find(~isnan(isMatchAll(:,pidx)));
    idxsAll = union(idxsAll,idxs);
    nSeg = sum(~isnan(isMatchAll(:,pidx)));
    nMatch = sum(isMatchAll(:,pidx) == 1);
    nMatchAll = nMatchAll + nMatch;
    nSegAll = nSegAll + nSeg;
    pck(pidx) = nMatch/nSeg*100;
    nseg(pidx) = nSeg;
    if (verbose)
        %fprintf('pidx %d: %f; seg: %d %d\n',pidx,nMatch/nSeg,nMatch,nSeg);
    end

end

pck(end) = nMatchAll/nSegAll*100;
nseg(end) = nSegAll;

if pslidx == -1
    if (verbose)
        fprintf('total: %f; seg: %d %d\n',nMatchAll/nSegAll, nMatchAll, nSegAll);
        fprintf(' Head & Shoulder & Elbow & Wrist & Hip & Knee  & Ankle & UBody & Total & # images \\\n');
        if (length(pidxs) == 14)
            fprintf('&%1.1f  & %1.1f  & %1.1f  & %1.1f  & %1.1f  & %1.1f & %1.1f & %1.1f & %1.1f & %d\\\n',(pck(7)+pck(8))/2,(pck(11)+pck(12))/2,(pck(10)+pck(13))/2,(pck(9)+pck(14))/2,(pck(3)+pck(4))/2,(pck(2)+pck(5))/2,(pck(1)+pck(6))/2,mean(pck(9:14)),pck(15),length(idxsAll));
        elseif (length(pidxs) == 6)
            fprintf('& - & %1.1f & %1.1f & %1.1f  & -  & -  & - & %1.1f & %1.1f & %d\\\n',(pck(11)+pck(12))/2,(pck(10)+pck(13))/2,(pck(9)+pck(14))/2, mean(pck(9:14)), pck(15),length(idxsAll));
        end
    end
else
    if (verbose)
%         fprintf('total: %f; seg: %d %d\n',nMatchAll/nSegAll, nMatchAll, nSegAll);
%         fprintf('Poselet & Torso & Upper & Lower & Upper & Fore- & Head  & Total & Total \\\n');
%         fprintf('        &       & Leg   & Leg&  & Arm   & arm   &       & PCP   & images\\\n');
        if (length(pidxs) == 10)
            fprintf('%s & &%1.1f  & %1.1f  & %1.1f  & %1.1f  & %1.1f  & %1.1f & %1.1f & %d\\\n',padZeros(num2str(pslidx),4),pck(5),(pck(2)+pck(3))/2,(pck(1)+pck(4))/2,(pck(8)+pck(9))/2,(pck(7)+pck(10))/2,pck(6),pck(11),length(idxsAll));
        elseif (length(pidxs) == 5)
            fprintf('&%1.1f  & - & - & %1.1f  & %1.1f  & -  & %1.1f & %d\\\n',pck(5),(pck(8)+pck(9))/2,(pck(7)+pck(10))/2, pck(11),length(idxsAll));
        elseif (length(pidxs) == 4)
            fprintf('& - & - & - & %1.1f  & %1.1f  & -  & %1.1f & %d\\\n',(pck(8)+pck(9))/2,(pck(7)+pck(10))/2, pck(11),length(idxsAll));
        end
    end
end

pck = [(pck(7)+pck(8))/2,(pck(11)+pck(12))/2,(pck(10)+pck(13))/2,(pck(9)+pck(14))/2,(pck(3)+pck(4))/2,(pck(2)+pck(5))/2,(pck(1)+pck(6))/2,mean(pck(9:14)),pck(15)];

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
clc;
clear all;
close all;

file = 'data.csv'; % Dataset

% Reading training file
data = dlmread(file);
label = data(:,end);

% Extracting positive data points
idx = (label==1);
pos_data = data(idx,:);
row_pos = size(pos_data,1);

% Extracting negative data points
neg_data = data(~idx,:);
row_neg = size(neg_data,1);

% Random permuation of positive and negative data points
p = randperm(row_pos);
n = randperm(row_neg);

% 80-20 split for training and test
tstpf = p(1:round(row_pos/5));
tstnf = n(1:round(row_neg/5));
trpf = setdiff(p, tstpf);
trnf = setdiff(n, tstnf);

train_data = [pos_data(trpf,:);neg_data(trnf,:)];
test_data = [pos_data(tstpf,:);neg_data(tstnf,:)];

% Decision Tree
prediction = SMOTEBoost(train_data,test_data,'tree',false);
disp ('    Label   Probability');
disp ('-----------------------------');
disp (prediction);

author: Sable

repo: mcbench-benchmarks


text:
function res = imdirfilter(img, varargin)
%IMDIRFILTER Apply and combine several directional filters
%
%   Deprecated, use "imDirectionlFilter" instead.
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 16/02/2004.
%

%   HISTORY
%   17/02/2004 debug, added some doc, manage different types of images
%   26/02/2007 cleanup code
%   03/10/2007 update doc
%   19/08/2009 fix bug in parsing arguments
%   18/06/2010 fix initialization for 'min' operator, update doc
%   02/12/2011 deprecated and replace by imDirectionalFilter

warning('imael:deprecated', ...
    'function "imdirfilter" has been replaced by "imDirectionalFilter"');

res = imDirectionalFilter(img, varargin{:});

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function theta = vectorAngle3d(v1, v2)
%VECTORANGLE3D Angle between two 3D vectors.
%
%   THETA = vectorAngle3d(V1, V2)
%   Computes the angle between the 2 3D vectors V1 and V2. The result THETA
%   is given in radians, between 0 and PI.
%
%
%   Example
%   % angle between 2 orthogonal vectors
%   vectorAngle3d([1 0 0], [0 1 0])
%   ans =
%       1.5708
%
%   % angle between 2 parallel vectors
%   v0 = [3 4 5];
%   vectorAngle3d(3*v0, 5*v0)
%   ans =
%       0
%
%   See also
%   vectors3d, vectorNorm3d, crossProduct3d
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-10-04,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% 2011-03-10 improve computation precision

% compute angle using arc-tangent to get better precision for angles near
% zero, see the discussion in:
% http://www.mathworks.com/matlabcentral/newsreader/view_thread/151925#381952
theta = atan2(vectorNorm3d(crossProduct3d(v1, v2)), sum(bsxfun(@times, v1, v2),2));

% equivalent to:
% v1 = normalizeVector3d(v1);
% v2 = normalizeVector3d(v2);
% theta = acos(dot(v1, v2, 2));

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function imdb = vid(prefix, opts)
savefn = fullfile(opts.dataDir, [prefix '.mat']);
if exist(savefn, 'file')
    imdb = load(savefn);
    return;
end
viddir = fullfile(opts.dataDir, 'VehicleID_V1.0');

% imgs, vid
img2vid = fullfile(viddir, 'attribute', 'img2vid.txt');
[imgs, vid] = textread(img2vid, '%d %d');
n = numel(imgs);

% model info
mid = zeros(n, 1);
vid2model = fullfile(viddir, 'attribute', 'model_attr.txt');
[vid_m, mid_m] = textread(vid2model, '%d %d');
for i = 1:length(vid_m)
    mid(vid == vid_m(i)) = mid_m(i) + 1;
end

% color info
cid = zeros(n, 1);
vid2color = fullfile(viddir, 'attribute', 'color_attr.txt');
[vid_c, cid_c] = textread(vid2color, '%d %d');
for i = 1:length(vid_c)
    mid(vid == vid_c(i)) = cid_c(i) + 1;
end

% train/test
sets = zeros(n, 4);  % [train test_small test_medium test_large]

trainFile = fullfile(viddir, 'train_test_split', 'train_list.txt');
[~, vid_train] = textread(trainFile, '%d %d');
sets(ismember(vid, vid_train), 1) = 1;
logInfo('train: %d imgs', sum(sets(:, 1)));

testSmallFile = fullfile(viddir, 'train_test_split', 'test_list_800.txt');
[~, vid_s] = textread(testSmallFile, '%d %d');
sets(ismember(vid, vid_s), 2) = 1;
logInfo('test-small: %d imgs', sum(sets(:, 2)));

testMediumFile = fullfile(viddir, 'train_test_split', 'test_list_1600.txt');
[~, vid_m] = textread(testMediumFile, '%d %d');
sets(ismember(vid, vid_m), 3) = 1;
logInfo('test-medium: %d imgs', sum(sets(:, 3)));

testLargeFile = fullfile(viddir, 'train_test_split', 'test_list_2400.txt');
[~, vid_l] = textread(testLargeFile, '%d %d');
sets(ismember(vid, vid_l), 4) = 1;
logInfo('test-large: %d imgs', sum(sets(:, 4)));

% resize to 256x256
a = tic; tic;
vid256 = fullfile(viddir, 'image_256x256');
if ~exist(vid256, 'dir'), mkdir(vid256); end
fn256 = cell(n, 1);
for i = 1:n
    fn256{i} = sprintf('%s/%07d.jpg', vid256, imgs(i));
    if ~exist(fn256{i}, 'file')
        im = imread(sprintf('%s/image/%07d.jpg', viddir, imgs(i)));
        im = imresize(im, [256 256]);
        imwrite(im, fn256{i});
    end
    if toc > 20
        logInfo(fn256{i}); tic;
    end
end
toc(a)

% imdb
imdb.images.data = fn256;
imdb.images.labels = single([vid mid cid]) ;
imdb.images.set = uint8(sets) ;
imdb.meta.sets = {'train', 'test_small', 'test_medium', 'test_large'} ;
save(savefn, '-struct', 'imdb');
end

author: kunhe

repo: FastAP-metric-learning


text:
function [out] = prep_addFeatures(dat1,dat2)
% prep_addFeatures (Pre-processing procedure):
%
% Description:
%     This function add the latter feature vector(dat2)
%     to the former one(dat1)
%
% Example:
% [out] = prep_addFeatures(dat1,dat2)
% [out] = prep_addFeatures({dat1,dat2,dat3,dat4,...})
%
% Input:
%     dat1 - Feature vector in a struct form
%     dat2 - Feature vector in a struct form, with same class as dat1
%
% Returns:
%     out - Updated feature vector added aditional features
%
%
% Seon Min Kim, 04-2016
% seonmin5055@gmail.com

if iscell(dat1)
    temp = dat1{1};
    for i=2:size(dat1,2)
        temp = prep_addFeatures(temp,dat1{i});
    end
    out = temp;
    return
end

if ~isfield(dat1,'x') || ~isfield(dat2,'x')
    warning('OpenBMI: Data structure must have a field named ''x''')
    return
end
if size(dat1.x,2) ~= size(dat2.x,2)
    warning('OpenBMI: Unmatched the number of trials')
    return
end
if ~isfield(dat1,'class') || ~isfield(dat1,'y_dec')
    warning('OpenBMI: Data structure should have a field named ''class'',''y_dec''')
end

r1 = rmfield(dat1,'x');
r2 = rmfield(dat2,'x');
if ~isequal(r1,r2)
    warning('OpenBMI: Two features must be extracted from the data with same class and trials')
    return
end

out = r1;
out.x = cat(1,dat1.x,dat2.x);

author: PatternRecognition

repo: OpenBMI


text:
function [D,u,X,div_X,phi,pre,B,t] = heat_geodesic(varargin)
  % HEAT_GEODESIC  geodesic distance approximation following the method of
  % "Geodesics in Heat" [Krane et al. 2013] D from all points V in the domain
  % (V,F) to a source point/set of points, gamma. The method is motivated by
  % heat using a time parameter t to guide the heat diffusion.
  %
  % D = heat_geodesic(V,F,gamma,t)
  % [D,u,X,div_X,phi,pre] = heat_geodesic(V,F,gamma,t,'ParamName',ParamValue)
  %
  % Inputs:
  %   V  #V by 3 set of vertex positions
  %   F  #F by 3 set of face indices
  %     or
  %      #T by 4 set of tetrahedra indices
  %   gamma  #gamma list of vertex indices of source points
  %   t  time parameter
  %   Optional:
  %     'BoundaryConditions'  followed by one of the following strings
  %       {'average'}:
  %       'dirichlet'  domain boundaries set to 0 when computing heat
  %         diffusion
  %       'neumann'  heat diffusion solved with implicit neumann conditions
  %       'average'  Uniform average of neumann and Dirichlet solutions
  %       'robin'  The original ArXiv paper called the method of averaging the
  %         Dirichlet and Neumann solutions 'robin'. This was later corrected
  %         (this solution cannot be created by boundary conditions to the heat
  %         equation alone: average solutions is a non-linear operation).  In
  %         any case, 'robin' is treated as 'average'.
  %     'Precomputation' Followed by pre struct as returned by this function
  %     'IntrinsicDelaunay'  followed by whether to use intrinsic Delaunay
  %       Laplacian {false}
  %     'Legacy' followed by bool telling whether to use Alec's legacy
  %       implmentation. In particular this is useful because the original
  %       paper is unclear how the boundary of the domain should be handle with
  %       respect to the seed locations (gamma). What if they overlap? What
  %       should be the boundary conditions for the final poisson solve?
  % Outputs:
  %   D  #V list of geodesic distances from vertices to gamma
  %   u  #V list of results of heat diffusion step
  %   X  #F by 3 list of reversed normalized gradients of u
  %   div_X  #V list od divergence of X
  %   phi  #V list of solution to final poisson equation solve
  %
  D=[];
  % Note: This is Alec's previous implmentation. He is still convinced that
  % this is more exact/correct despite not what's written in the article (and
  % implying that refactoring is necessary).
  legacy = false;
  use_intrinsic = false;

  % mandatory input
  V = varargin{1};
  F = varargin{2};

  % number of domain vertices
  n = size(V,1);
  % simplex size
  ss = size(F,2);
  gamma = varargin{3};
  if nargin>=4 && ~isempty(varargin{4})
    t = varargin{4};
  else
    switch ss
    case 3
      % Section 3.1.1
      AM = sum(doublearea(V,F))/2;
      sF = size(F,1);
      c = 5;
      t = c * AM / sF;
      %t = 20*mean(doublearea(V,F));
    case 4
      AM = sum(volume(V,F));
      sF = size(F,1);
      c = 1e3;
      t = c * AM / sF;
      %t = 20*mean(volume(V,F));
    end
  end

  % option parameter default values
  bc_type = 'average';
  pre = [];
  pre.L = [];
  % precomputation for Dirichlet solve
  pre.D = [];
  % precomputation for Neumann solve
  pre.N = [];
  % precomputation for Poisson solve
  pre.poisson = [];
  u = [];

  ii = 5;
  while(ii <= nargin)
    switch varargin{ii}
    case 'BoundaryConditions'
      ii = ii + 1;
      assert(ii<=nargin);
      bc_type = varargin{ii};
    case 'Precomputation'
      ii = ii + 1;
      assert(ii<=nargin);
      % skip empty input
      if ~isempty(varargin{ii})
        pre = varargin{ii};
      end
    case 'Legacy'
      ii = ii + 1;
      assert(ii<=nargin);
      legacy = varargin{ii};
    case 'u'
      ii = ii + 1;
      assert(ii<=nargin);
      u = varargin{ii};
    case 'IntrinsicDelaunay'
      ii = ii + 1;
      assert(ii<=nargin);
      use_intrinsic = varargin{ii};
    otherwise
      error(['Unsupported parameter: ' varargin{ii}]);
    end
    ii = ii + 1;
  end

  % Algorithm 1 in Section 3:
  %   Integrate the heat flow ?? = ??u for time t
  %   Evaluate the vector field X = -???u/|???u|
  %   Solve the Poisson equation ???? = ??????X

  % Integrate the heat flow: "We discretize the heat equation from step I of
  % Algorithm 1 using a single backward Euler step"

  % "... with initial conditions u0 = ??(x)"
  % "Note that a Dirac delta appears as a literal one in this system since we
  % are effectively working with integrated quantities"
  u0 = zeros(n,1);
  u0(gamma) = 1;

  % "where ???????? is one third the area of all triangles incident on vertex ...
  % where ???? ??? R|????|??|????| is a diagonal matrix containing the vertex areas"
  if isempty(pre.L)
    if use_intrinsic
    pre.L = intrinsic_delaunay_cotmatrix(V,F);
    else
    pre.L = cotmatrix(V,F);
    end
    pre.M = massmatrix(V,F,'barycentric');
    pre.Q = pre.M - t*pre.L;
    pre.G = grad(V,F);
    pre.Div = div(V,F);
    pre.Tik = 1e-8*speye(size(V,1));
    pre.P = -pre.L*0.5 + pre.Tik;
    % get outline ("boundary") of mesh
    switch ss
    case 3
      pre.out = unique(reshape(outline(F),[],1));
    case 4
      pre.out = unique(boundary_faces(F));
    end
  end

  % 11/7/2017: Keenan writes that the mass matrix should not be here.
  % B = M*u0;
  B = u0;

  if isempty(u)
    if strcmp(bc_type,'dirichlet') || strcmp(bc_type,'average') || strcmp(bc_type,'robin')

      if legacy
        % remove any gamma from outline
        pre.out = setdiff(pre.out(:),gamma(:));
        % find all boundary vertices and append to gamma
        b = [gamma(:); pre.out];
        % "zero Dirichlet conditions"
        bc = [ones(numel(gamma),1); zeros(numel(b)-numel(gamma),1)];
        [uD,pre.D] = min_quad_with_fixed(pre.Q,B,b,bc,[],[],pre.D);
      else
        b = pre.out;
        % Q: How should we deal with gamma ??? out ?
        % This gives *reasonable* results, but doesn't look right for fixing
        % entire boundary. Probably more because of the lack of correct boundary
        % conditions for the final poisson solve.
        bc = -1*ismember(pre.out,gamma);
        %assert(~any(ismember(gamma,out)));
        [uD,pre.D] = min_quad_with_fixed(pre.Q,B,b,bc,[],[],pre.D);
      end
    end
    if strcmp(bc_type,'neumann') || strcmp(bc_type,'average') || strcmp(bc_type,'robin')
      % See Figure 9, pretty sure these are implicit ???x/???n = 0 neumann
      % conditions, though it is not stated in the text. At best, "Neumann
      % conditions prevent heat from flowing out of the domain..."
      if legacy
        b = gamma(:);
        bc = ones(numel(gamma),1);
        [uN,pre.N] = min_quad_with_fixed(pre.Q,B,b,bc,[],[],pre.N);
      else
        [uN,pre.N] = min_quad_with_fixed(pre.Q,B,[],[],[],[],pre.N);
      end
    end

    if strcmp(bc_type,'natural')
      K = keenan(V,F);
      Q = M - t*K;
      B = M*u0;
      G = grad(V,F);
      switch size(F,2)
      case 4
        vol = volume(V,F);
      case 3
        vol = doublearea(V,F);
      end
      vol = vol/sum(vol);
      % kill off affine functions
      A = [
        sum(M)/sum(M(:)); ...
        kron(speye(size(V,2)),vol)'*G];
      uT = min_quad_with_fixed(Q,B,[],[],A,[0;zeros(size(V,2),1)]);
    end

    switch bc_type
    case 'dirichlet'
      u = uD;
    case 'neumann'
      u = uN;
    case {'average','robin'}
      % "We advocate the use of the Robin boundary conditions obtained by taking
      % the mean of the Neumann solution uN and the Dirichlet solution uD, i.e.,
      % u = 0.5*(uN + uD)"
      u = 0.5*(uN+uD);
    case 'natural'
      u = uT;
    otherwise
      error(['Unsupported BoundaryCondtions value: ' bc_type]);
    end
  end

  robust_X = false;
  if robust_X
    X = robust_heat_gradients(V,F,uN,uD,t);
    warning('robust');
  else
    % Evaluate the vector field X
    grad_u = reshape(pre.G*u,size(F,1),size(V,2));
    %grad_u_norm = sqrt(sum(grad_u.^2,2));
    %% normalize grad_u
    %normalized_grad_u = bsxfun(@rdivide,grad_u,grad_u_norm);
    %% correct any zero-norm gradients
    %normalized_grad_u(grad_u_norm == 0,:) = 0;
    % normalize reverse direction; normalizerow will use robust norm that
    % avoids underflow.
    X = -normalizerow(grad_u);
    X(isnan(X)) = 0;
  end

  % Solve the Poisson equation
  % divergence of X
  div_X = pre.Div*X(:);
  if legacy
    [phi,pre.poisson] = min_quad_with_fixed( ...
      -pre.L*0.5,2*div_X,gamma,zeros(numel(gamma),1),[],[],pre.poisson);
  else
    [phi,pre.poisson] = ...
      min_quad_with_fixed(pre.P,2*div_X,[],[],[],[],pre.poisson);
    %Aeq = sparse(ones(1,size(V,1)));
    %Beq = 0;
    %[phi,pre.poisson] = ...
    %  min_quad_with_fixed(-pre.L*0.5,2*div_X,[],[],Aeq,Beq,pre.poisson);
  end
  D = phi;
  % "Note that ???? is unique only up to an additive constant and should be
  % shifted such that the smallest distance value is zero."
  %
  % D should be zero at gamma
  D = D - mean(D(gamma));
  % Flip sign so farthest point is positive distance
  [~,mi] = max(abs(D));
  D = D*sign(D(mi));

end

author: alecjacobson

repo: gptoolbox


text:
function x = reshape( x, varargin )

%   Disciplined convex/geometric programming information for RESHAPE:
%       RESHAPE imposes no convexity restrictions on its arguments.

%
% Check size arguments
%

switch nargin,
    case {0,1},
        error( 'Not enough input arguments.' );
    case 2,
        [ temp, sz ] = cvx_check_dimlist( varargin{1}, true );
        if ~temp,
            error( 'Second argument must be a valid dimension list.' );
        end
    otherwise,
        [ temp, sz ] = cvx_check_dimlist( varargin, true );
        if ~temp,
            error( 'Second and subsequent arguments must be nonnegative integers.' );
        end
end

%
% Quick exit if the size remains the same
%

sx = x.size_;
if isequal( sx, sz ),
    return;
end

%
% Confirm compatible reshape
%

px = prod( sx );
tt = isnan( sz );
if any( tt ),
    if nnz( tt ) > 1,
        error( 'Size can only have one unknown dimension.' );
    end
    tmp = px / prod( sz( ~tt ) );
    if isnan( tmp ) || tmp ~= floor( tmp ),
        error( 'Product of known dimensions, %d, not divisible into total number of elements, %d.', prod( sz( ~tt ) ), px );
    end
    sz( tt ) = tmp;
elseif px ~= prod( sz ),
    error( 'To RESHAPE the number of elements must not change.' );
end

%
% Perform the resize
%

x = cvx( sz, x.basis_ );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
classdef neural_network < handle
    properties
    train_file;
    test_file;
    rounds;
    train_data;
    test_data;
    target;
    test_target;
    first_wghts;
    op_wghts;
    op_acivation;
    phi;
    test_phi;
    attr_no;
    no_of_samples;
    no_of_units;
    unique_class;
    op_units;
    no_of_layers;
    act;
    net;
    hidden_wghts;
    net_1;
    act_1;
    op_act;
    delta;
    step;
    del;
    delta_op;
    tot_acc;
    act_res;
    acc_res;
    pred_res;
    end
    methods (Static)
        function [obj] = neural_network(training_file, test_file, layers, units_per_layer, rounds)
            obj.train_file = training_file;
            obj.test_file = test_file;
            obj.rounds = rounds;
            obj.train_data = load(obj.train_file);
            obj.test_data = load(obj.test_file);
            obj.attr_no = size(obj.train_data, 2);
            obj.no_of_samples = size(obj.train_data, 1);
            obj.no_of_units = units_per_layer;
            obj.target = obj.train_data(1: end, end);
            obj.test_target = obj.test_data(1: end, end);
            obj.unique_class = unique(obj.target);
            obj.unique_class = sort(obj.unique_class);
            obj.op_units = size(obj.unique_class, 1);
            obj.no_of_layers = layers;
            obj.step = 0.98;
            obj.tot_acc = 0;
            obj.del = zeros(obj.no_of_units, obj.no_of_layers-2);
        end

        function [obj] = initialise(obj)
            ones = zeros(obj.no_of_samples, 1);
            ones(ones == 0) = 1;

            obj.phi = obj.train_data(1: end, 1: end-1);
            maximum = max(obj.phi(:));
            obj.phi = (obj.phi)/(maximum);
            obj.phi = [obj.phi ones];

            ones = zeros(size(obj.test_data, 1), 1);
            ones(ones == 0) = 1;
            obj.test_phi = obj.test_data(1: end, 1: end-1);
            maximum = max(obj.test_phi(:));
            obj.test_phi = (obj.test_phi)/(maximum);
            obj.test_phi = [obj.test_phi ones];

            obj.act_res = zeros(size(obj.test_data, 1), 1);
            obj.acc_res = zeros(size(obj.test_data, 1), 1);
            obj.pred_res = zeros(size(obj.test_data, 1), 1);

            obj.first_wghts = -0.05 + (0.05- (-0.05)).*rand(size(obj.phi, 2), obj.no_of_units);

            obj.act = zeros(obj.no_of_units+1, obj.no_of_layers-2);
            obj.net = zeros(obj.no_of_units, obj.no_of_layers-2);

            if obj.no_of_layers > 2
                obj.op_wghts = -0.05 + (0.05- (-0.05)).*rand(obj.no_of_units+1, obj.op_units);

            else
                obj.op_wghts = -0.05 + (0.05- (-0.05)).*rand(size(obj.phi, 2), obj.op_units);

            end

            obj.hidden_wghts = -0.05 + (0.05- (-0.05)).*rand(obj.no_of_units+1, obj.no_of_units, obj.no_of_layers-3);
            obj.delta_op = zeros(obj.op_units, 1);
        end

        function [obj] = feed_forward(obj, round)
            for row = 1:size(obj.phi, 1)
                for layer = 1:obj.no_of_layers-1

                    if layer == obj.no_of_layers-1 && obj.no_of_layers > 2

                        obj.net_1 =  transpose(obj.op_wghts(1:end, 1:end)) * (obj.act(1:end, layer-1));
                        obj.op_act = logsig(obj.net_1);

                    elseif layer == obj.no_of_layers-1 && obj.no_of_layers == 2

                        obj.net_1  = transpose(obj.op_wghts(1:end, 1:end)) * transpose(obj.phi(row, 1:end));
                        obj.op_act = logsig(obj.net_1);

                    elseif layer == 1

                        obj.net(1:end, layer) = transpose(obj.first_wghts) * transpose(obj.phi(row, 1:end));
                        obj.act(1:end-1, layer) = logsig(obj.net(1:end, layer));
                        obj.act(end, layer) = 1;

                    elseif layer > 1 && layer < obj.no_of_layers-1

                        obj.net(1:end, layer) = transpose(obj.hidden_wghts(1:end, 1:end, layer-1)) * (obj.act(1:end, layer-1));
                        obj.act(1:end-1, layer) = logsig(obj.net(1:end, layer));
                        obj.act(end, layer) = 1;

                    end
                end
                obj = obj.back_prop(obj, round, row);
            end
        end

        function [obj] = back_prop(obj, round, row)
            tg = obj.target(row, 1);

            pos = find(obj.unique_class == tg);

            rate = power(obj.step,round);

            for layer = obj.no_of_layers-1: -1: 1
                if layer == obj.no_of_layers-1 && obj.no_of_layers > 2
                    for unit = 1:obj.op_units
                        if unit == pos
                             t = 1;
                        else
                             t = 0;
                        end
                        o = obj.op_act(unit, 1);
                        delta = (o-t)*o*(1-o);
                        obj.delta_op(unit, 1) = delta;
                        for wght = 1:obj.no_of_units+1
                            obj.op_wghts(wght, unit) = obj.op_wghts(wght, unit) - (rate*delta*obj.act(wght, end));

                        end
                    end

                elseif layer == obj.no_of_layers-1 && obj.no_of_layers == 2
                    for unit = 1:obj.op_units
                        if unit == pos
                             t = 1;
                        else
                             t = 0;
                        end
                        o = obj.op_act(unit, 1);
                        delta = (o-t)*o*(1-o);
                        obj.delta_op(unit, 1) = delta;
                        for wght = 1:size(obj.phi, 2)
                            obj.op_wghts(wght, unit) = obj.op_wghts(wght, unit) - (rate*delta*obj.phi(row, wght));
                            %fprintf('unit=%d,wght=%d, delta=%6.4f, obj.act(wght, end)=%6.4f \n', unit, wght, delta, obj.act(wght, end))
                        end
                    end

                elseif layer == obj.no_of_layers-2 && layer > 1 && obj.no_of_layers >3

                    for unit = 1:obj.no_of_units
                        obj.del(unit, layer) = transpose(obj.delta_op)*transpose(obj.op_wghts(unit, 1:end));
                        o = obj.act(unit, end);
                        obj.del(unit, layer) = obj.del(unit, layer)*o*(1-o);
                        obj.delta = obj.del(unit, layer);

                        for wght = 1:obj.no_of_units+1
                            obj.hidden_wghts(wght, unit, end) = obj.hidden_wghts(wght, unit, end)-rate*obj.delta*obj.act(wght, end-1);
                        end
                    end

                elseif layer < obj.no_of_layers-2 && layer > 1 && obj.no_of_layers >3

                    for unit = 1:obj.no_of_units

                        obj.del(unit, layer) = transpose(obj.del(1:end, layer+1))*transpose(obj.hidden_wghts(unit, 1:end, layer));
                        o = obj.act(unit, layer);
                        obj.del(unit, layer) = obj.del(unit, layer)*o*(1-o);
                        obj.delta = obj.del(unit, layer);

                        for wght = 1:obj.no_of_units+1
                            obj.hidden_wghts(wght, unit, layer-1) = obj.hidden_wghts(wght, unit, layer-1)-rate*obj.delta*obj.act(wght, layer-1);
                        end
                    end
                    %disp(obj.hidden_wghts(:, :, layer-1))

                elseif layer == 1 && obj.no_of_layers == 3

                    for unit = 1:obj.no_of_units
                        obj.del(unit, layer) = transpose(obj.delta_op)*transpose(obj.op_wghts(unit, 1:end));
                        o = obj.act(unit, 1);
                        obj.del(unit, layer) = obj.del(unit, layer)*o*(1-o);
                        obj.delta = obj.del(unit, layer);
                        for wght = 1:size(obj.phi, 2)
                            obj.first_wghts(wght, unit) = obj.first_wghts(wght, unit)-(rate*obj.delta*obj.phi(row, wght));

                        end
                    end

                elseif layer == 1 && obj.no_of_layers > 3

                    for unit = 1:obj.no_of_units

                        obj.del(unit, layer) = transpose(obj.del(1:end, layer+1))*transpose(obj.hidden_wghts(unit, 1:end, 1));
                        o = obj.act(unit, 1);
                        obj.del(unit, layer) = obj.del(unit, layer)*o*(1-o);
                        obj.delta = obj.del(unit, layer);
                        for wght = 1:size(obj.phi, 2)
                            obj.first_wghts(wght, unit) = obj.first_wghts(wght, unit)-rate*obj.delta*obj.phi(row, wght);
                        end
                    end

                end

            end
        end
        function [obj] = testing(obj)
            for row = 1:size(obj.test_phi, 1)
                for layer = 1:obj.no_of_layers-1

                    if layer == obj.no_of_layers-1 && obj.no_of_layers > 2

                        obj.net_1 =  transpose(obj.op_wghts(1:end, 1:end)) * (obj.act(1:end, layer-1));
                        obj.op_act = logsig(obj.net_1);

                    elseif layer == obj.no_of_layers-1 && obj.no_of_layers == 2

                        obj.net_1  = transpose(obj.op_wghts(1:end, 1:end)) * transpose(obj.test_phi(row, 1:end));
                        obj.op_act = logsig(obj.net_1);

                    elseif layer == 1

                        obj.net(1:end, layer) = transpose(obj.first_wghts) * transpose(obj.test_phi(row, 1:end));
                        obj.act(1:end-1, layer) = logsig(obj.net(1:end, layer));
                        obj.act(end, layer) = 1;

                    elseif layer > 1 && layer < obj.no_of_layers-1

                        obj.net(1:end, layer) = transpose(obj.hidden_wghts(1:end, 1:end, layer-1)) * (obj.act(1:end, layer-1));
                        obj.act(1:end-1, layer) = logsig(obj.net(1:end, layer));
                        obj.act(end, layer) = 1;

                    end
                end
                max_act = max(obj.op_act(:));
                pos = find(obj.op_act ==  max_act);
                predicted = obj.unique_class(pos);
                test_t = obj.test_target(row, 1);

                acc= 0;
                if predicted == test_t
                    acc = 1;
                end

                obj.act_res(row, 1) = max_act;
                obj.acc_res(row, 1) = acc;
                obj.pred_res(row, 1) = predicted;

                obj.tot_acc = obj.tot_acc+acc;
                fprintf('ID=%5d, predicted=%3d, true=%3d, accuracy=%4.2f \n', row-1, predicted, test_t, acc);

            end
            fprintf('classification accuracy=%6.4f  ', (sum(obj.acc_res)/size(obj.test_phi, 1)));
        end

    end
end

repo: Machine-Learning-Models


text:
%% Copyright (C) 2014-2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @deffn  Command syms
%% @deffnx Command syms @var{x}
%% @deffnx Command syms {@var{x} @var{y} @dots{}}
%% @deffnx Command syms @var{f(x)}
%% @deffnx Command syms {@var{x} @var{asm}}
%% @deffnx Command syms {@var{x} @var{asm} @var{asm2} @dots{}}
%% Create symbolic variables and symbolic functions.
%%
%% This is a convenience function.  For example:
%% @example
%% syms x y z
%% @end example
%% instead of:
%% @example
%% @group
%% x = sym('x');
%% y = sym('y');
%% z = sym('z');
%% @end group
%% @end example
%%
%% The last arguments can provide one or more assumptions (type or
%% restriction) on the variable (@pxref{sym}).
%% @example
%% @group
%% syms x y z positive
%% syms n positive even
%% @end group
%% @end example
%%
%% Symfuns represent abstract or concrete functions.  Abstract
%% symfuns can be created with @code{syms}:
%% @example
%% syms f(x)
%% @end example
%% Here @code{x} is created in the callers workspace,
%% as a @strong{side effect}.
%%
%% Called without arguments, @code{syms} displays a list of
%% all symbolic functions defined in the current workspace.
%% @example
%% @group
%% syms x y z
%% syms
%%   @print{} Symbolic variables in current scope:
%%   @print{}   x
%%   @print{}   y
%%   @print{}   z
%% @end group
%% @end example
%%
%% Using @code{syms} to create a new symbol with different assumptions
%% does not automatically replace instances of the old symbol in the
%% caller's workspace.  For example suppose we make an absolute value
%% expression:
%% @example
%% @group
%% syms x
%% f = abs (x)
%%   @result{} f = (sym) │x│
%% @end group
%% @end example
%% If we make a new positive @code{x}, @code{f} still contains the old
%% symbol:
%% @example
%% @group
%% syms x positive
%% simplify (f)
%%   @result{} (sym) │x│
%% @end group
%% @end example
%% Note: this behaviour is slightly different from Matlab's Symbolic
%% Math Toolbox, which does change existing symbols.
%% If that behaviour is desired, @pxref{@@sym/assume} which modifies
%% existing expressions in the caller's workspace:
%% @example
%% @group
%% syms x
%% assume x positive
%% f
%%   @result{} f = (sym) x
%% @end group
%% @end example
%%
%% Caution: On Matlab, you may not want to use @code{syms} within
%% functions.
%%   In particular, if you shadow a function name, you may get
%%   hard-to-track-down bugs.  For example, instead of writing
%%   @code{syms alpha} use @code{alpha = sym('alpha')} within functions.
%%   [https://www.mathworks.com/matlabcentral/newsreader/view_thread/237730]
%%
%% @seealso{sym, assume}
%% @end deffn

%% Author: Colin B. Macdonald
%% Keywords: symbolic, symbols, CAS

function syms(varargin)

  %% No inputs
  %output names of symbolic vars
  if (nargin == 0)
    S = evalin('caller', 'whos');
    disp('Symbolic variables in current scope:')
    for i=1:numel(S)
      %S(i)
      if strcmp(S(i).class, 'sym')
        disp(['  ' S(i).name])
      elseif strcmp(S(i).class, 'symfun')
        % FIXME improve display of symfun
        disp(['  ' S(i).name ' (symfun)'])
      end
    end
    return
  end

  %% Find assumptions
  valid_asm = assumptions('possible');
  last = -1;
  for n=1:nargin
    assert(ischar(varargin{n}), 'syms: expected string inputs')
    if (ismember(varargin{n}, valid_asm))
      if (last < 0)
        last = n - 1;
      end
    elseif (strcmp(varargin{n}, 'clear'))
      error ('OctSymPy:oldsyntax', ...
             'Old "syms x clear" not supported; use "assume x clear" or "assume(x, ''clear'')"')
    elseif (last > 0)
      error('syms: cannot have symbols after assumptions')
    end
  end

  if (last < 0)
    asm = {};
    exprs = varargin;
  elseif (last == 0)
    error('syms: cannot have only assumptions w/o symbols')
  else
    asm = varargin((last+1):end);
    exprs = varargin(1:last);
  end

  % loop over each input
  for i = 1:length(exprs)
    expr = exprs{i};

    % look for parenthesis: check if we're making a symfun
    if (isempty (strfind (expr, '(') ))  % no
      assert(isvarname(expr)); % help prevent malicious strings
      assignin('caller', expr, sym(expr, asm{:}))

    else  % yes, this is a symfun
      assert(isempty(asm), 'mixing symfuns and assumptions not supported')
      % regex matches: abc(x,y), f(var), f(x, y, z), f(r2d2), f( x, y )
      % should not match: Rational(2, 3), f(2br02b)
      assert(~isempty(regexp(expr, '^\w+\(\s*[A-z]\w*(,\s*[A-z]\w*)*\s*\)$')), ...
             'invalid symfun expression')

      T = regexp (expr, '^(\w+)\((.*)\)$', 'tokens');
      assert (length (T) == 1 && length (T{1}) == 2)
      name = T{1}{1};
      varnames = strtrim (strsplit (T{1}{2}, ','));

      vars = {};
      for j = 1:length (varnames)
        % is var in the workspace?
        try
          v = evalin ('caller', varnames{j});
          create = false;
        catch
          create = true;
        end
        % if var was defined, is it a symbol with the right name?
        if (~ create)
          if (~ isa (v, 'sym'))
            create = true;
          elseif (~ strcmp (v.flat, varnames{j}))
            create = true;
          end
        end
        if (create)
          v = sym (varnames{j});
          assignin ('caller', varnames{j}, v);
        end
        vars{j} = v;
      end

      cmd = { 'f, vars = _ins'
              'return Function(f)(*vars)' };
      s = pycall_sympy__ (cmd, name, vars);

      sf = symfun(s, vars);
      assignin('caller', name, sf);
    end
  end
end

%!test
%! %% assumptions
%! syms x real
%! x2 = sym('x', 'real');
%! assert (isequal (x, x2))

%!test
%! % assumptions and clearing them on a symbol
%! syms x real
%! assert (~isempty (assumptions (x)))
%! syms x
%! assert (isempty (assumptions (x)))

%!test
%! % Note SMT would clear syms in existing expressions
%! syms x real
%! f = {x {2*x} cos(x/2)};
%! assert (~isempty (assumptions (f)))
%! syms x
%! % but we do not: this would have to toggle for pure SMT compat
%! assert (~isempty (assumptions (f)))
%! % assert (isempty (assumptions (f)))

%!error <not supported>
%! syms x clear

%!error <symbols after assumptions>
%! syms x positive y

%!error <symbols after assumptions>
%! % this sometimes catches typos or errors in assumption names
%! % (if you need careful checking, use sym not syms)
%! syms x positive evne

%!error <cannot have only assumptions>
%! syms positive integer

%!test
%! % does not create a variable called positive
%! syms x positive integer
%! assert (logical(exist('x', 'var')))
%! assert (~logical(exist('positive', 'var')))

%!test
%! % Issue #885
%! syms S(x) I(x) O(x)

%!test
%! % Issue #290
%! syms FF(x)
%! syms ff(x)
%! syms Eq(x)

%!test
%! % Issue #290
%! syms beta(x)

%!test
%! syms x real
%! syms f(x)
%! assert (~ isempty (assumptions (x)))

%!test
%! syms x real
%! f(x) = symfun(sym('f(x)'), x);
%! assert (~ isempty (assumptions (x)))
%! assert (~ isempty (assumptions (argnames (f))))

repo: octsympy


text:
classdef TestStylization
    %TestStylization

    methods (Static)
        function test_1
            img = imread(fullfile(mexopencv.root(),'test','lena.jpg'));
            out = cv.stylization(img);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_2
            img = imread(fullfile(mexopencv.root(),'test','lena.jpg'));
            out = cv.stylization(img, 'SigmaS',60 ,'SigmaR',0.45);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_error_argnum
            try
                cv.stylization();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
classdef MetaDataTagsClass  < matlab.mixin.Copyable

    % SNIRF-spec class properties
    properties
        tags
    end

    % Non-SNIRF class properties
    properties (Access = private)
        filename
        fileformat
    end

    methods

        % -------------------------------------------------------
        function obj = MetaDataTagsClass(varargin)
            % Set class properties not part of the SNIRF format
            obj.fileformat = 'hdf5';

            obj.tags.SubjectID = 'default';
            obj.tags.MeasurementDate = datestr(now,29);
            obj.tags.MeasurementTime = datestr(now,'hh:mm:ss');
            obj.tags.LengthUnit = 'mm';
            obj.tags.TimeUnit = 'unknown';
            obj.tags.FrequencyUnit = 'unknown';
            obj.tags.AppName  = 'snirf-homer3';
            obj.tags.SnirfDraft = '3';

            if nargin==1
                obj.filename = varargin{1};
                obj.Load();
            end
        end

        % -------------------------------------------------------
        function err = LoadHdf5(obj, fileobj, location)
            err = 0;

            % Arg 1
            if ~exist('fileobj','var') || (ischar(fileobj) && ~exist(fileobj,'file'))
                fileobj = '';
            end

            % Arg 2
            if ~exist('location', 'var') || isempty(location)
                location = '/nirs/metaDataTags';
            elseif location(1)~='/'
                location = ['/',location];
            end

            % Error checking
            if ~isempty(fileobj) && ischar(fileobj)
                obj.filename = fileobj;
            elseif isempty(fileobj)
                fileobj = obj.filename;
            end
            if isempty(fileobj)
               err = -1;
               return;
            end

            %%%%%%%%%%%% Ready to load from file

            try
                % Reset tags
                obj.tags = struct();

                % Open group
                [gid, fid] = HDF5_GroupOpen(fileobj, location);

                metaDataStruct = h5loadgroup(gid);
                tags = fieldnames(metaDataStruct); %#ok<*PROPLC>
                for ii=1:length(tags)
                    eval(sprintf('obj.tags.%s = metaDataStruct.%s;', tags{ii}, tags{ii}));
                end

                HDF5_GroupClose(fileobj, gid, fid);

                % Detect old or invalid metadatatag data
                assert(obj.IsValid());

            catch

                err = -1;

            end

        end

        % -------------------------------------------------------
        function SaveHdf5(obj, fileobj, location) %#ok<*INUSD>
            % Arg 1
            if ~exist('fileobj', 'var') || isempty(fileobj)
                error('Unable to save file. No file name given.')
            end

            % Arg 2
            if ~exist('location', 'var') || isempty(location)
                location = '/nirs/metaDataTags';
            elseif location(1)~='/'
                location = ['/',location];
            end

            if ~exist(fileobj, 'file')
                fid = H5F.create(fileobj, 'H5F_ACC_TRUNC', 'H5P_DEFAULT', 'H5P_DEFAULT');
                H5F.close(fid);
            end
            props = propnames(obj.tags);
            for ii=1:length(props)
                eval(sprintf('hdf5write_safe(fileobj, [location, ''/%s''], obj.tags.%s);', props{ii}, props{ii}));
            end
        end

        % -------------------------------------------------------
        function b = IsValid(obj)
            b = false;

            % Use latest required fields to determine if we're loading old
            % metaDataTag format version of SNIRF spec
            if ~isproperty(obj.tags, 'FrequencyUnit')
                return;
            end

            b = true;
        end

        % -------------------------------------------------------
        function B = eq(obj, obj2)
            B = false;
            props1 = propnames(obj.tags);
            props2 = propnames(obj2.tags);
            for ii=1:length(props1)
                if ~isproperty(obj2.tags, props1{ii})
                    return;
                end
                if eval(sprintf('~strcmp(obj.tags.%s, obj2.tags.%s)', props1{ii}, props1{ii}))
                    return;
                end
            end
            for ii=1:length(props2)
                if ~isproperty(obj.tags, props2{ii})
                    return;
                end
                if eval(sprintf('~strcmp(obj.tags.%s, obj2.tags.%s)', props2{ii}, props2{ii}))
                    return;
                end
            end
            B = true;
        end

        % -------------------------------------------------------
        function Add(obj, key, value)
            key(key==' ') = '';
            eval(sprintf('obj.tags.%s = value', key));
        end

        % ----------------------------------------------------------------------------------
        function tags = Get(obj)
            fields = propnames(obj.tags);
            tags = repmat(struct('key','','value',[]),length(fields),1);
            for ii=1:length(fields)
                eval(sprintf('tags(ii).key = ''%s'';', fields{ii}));
                eval(sprintf('tags(ii).value = obj.tags.%s;', fields{ii}));
            end
        end

        % ----------------------------------------------------------------------------------
        function nbytes = MemoryRequired(obj)
            nbytes = 0;
            fields = propnames(obj.tags);
            for ii=1:length(fields)
                nbytes = nbytes + eval(sprintf('sizeof(obj.tags.%s)', fields{ii}));
            end
        end

    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function linpack_s_test04 ( )

%*****************************************************************************80
%
%% TEST04 tests SGBCO.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  ml = 1;
  mu = 1;
  lda = 2 * ml + mu + 1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04\n' );
  fprintf ( 1, '  For a general banded matrix,\n' );
  fprintf ( 1, '  SGBCO estimates the reciprocal condition number.\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Set the matrix A.
%
  m = ml + mu + 1;
  fprintf ( 1, '  The bandwidth of the matrix is %d\n', m );

  for j = 1 : n
    a(m-1,j) = -1.0;
    a(m,  j) =  2.0;
    a(m+1,j) = -1.0;
  end
%
%  Estimate the condition.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimate the condition.\n' );

  [ a, ipivot, rcond, z ] = sgbco ( a, lda, n, ml, mu );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimated reciprocal condition = %f\n', rcond );

  return
end


text:
function [suppressionIndex, pos_Vol, neg_Vol, pos_Vol_pRF, neg_Vol_pRF] = rmGetDoGSuppressionIndex(model,varargin)
%[suppressionIndex, pos_Vol, neg_Vol, pos_Vol_pRF, neg_Vol_pRF] = rmGetDoGSuppressionIndex(model,ROIcoords)
%
% rmGetDoGSuppressionIndex - calculates the suppression Index from the DoG model
% Calculates both the volume of the positive gaussian and the negative gaussian that make the
% pRF, it only takes into account the part of the gaussians that fall inside the stimulus window.
% The suppression index is the ratio of the negative and the positive
% volumes. Higher values indicate more suppression. (JOV Zuiderbaan et al.
% 2012)
%
% !! Needs at least the size of the stimulus (sts) as an input in varargin !!
% either by its direct value or by the params of the model:
% [suppressionIndex, pos_Vol, neg_Vol, pos_Vol_pRF, neg_Vol_pRF] = rmGetDoGSuppressionIndex(model,'sts',6.25)
% [suppressionIndex, pos_Vol, neg_Vol, pos_Vol_pRF, neg_Vol_pRF] = rmGetDoGSuppressionIndex(model,'rm',vw.rm)
%
% suppressionIndex          : ratio of the volumes of the negative and the positive gaussian that fall inside the stimuluswindow
% pos_Vol                   : volume of the positive gaussian that falls inside the stimuluswindow
% neg_Vol                   : volume of the negative gaussian that falls inside the stimuluswindow
% pos_Vol_pRF               : volume of the positive part of the pRF that falls inside the stimuluswindow
% neg_Vol_pRF               : volume of the negative part of the pRF that falls inside the stimuluswindow
%
% WZ  02/12: Wrote it

% set default parameters
params.sampleRate = 0.125;      % samplerate of the grid to make the gaussians
params.SI_Centered = false;     % by default the pRF has its original position on the grid, you can make it centered on (0,0)
                                % by giving this parameter the value true, using varargin

if nargin > 1,
    addArg = varargin;
    if numel(addArg) == 1,
        addArg=addArg{1};
    end;
else
    addArg = [];
end;

% parse command line inputs:
params = rmProcessVarargin(params,addArg);

% the stimulus size needs to be given
if ~isfield(params,'stimSize') || isempty(params.stimSize),
    disp('ERROR: Need size of the stimulus');
    error('Need size of the stimulus');
end

% make the stimuluswindow
mygridx = (-params.stimSize:params.sampleRate:params.stimSize);
[X2, Y2] = meshgrid(mygridx,mygridx);
Y2 = flipud(Y2);
X2 = X2(:);
Y2 = Y2(:);

ecc = sqrt(X2.^2 + Y2.^2);
keep = ecc > -params.stimSize & ecc < params.stimSize;
X = X2(keep);
Y = Y2(keep);

gridsize = params.sampleRate.^2;

sigma = model.sigma.major;
sigma2 = model.sigma2.major;
beta1 = model.beta(1,:,1);
beta2 = model.beta(1,:,2);
x = model.x0;
y = model.y0;

% if ROI indices are given, take only the data for that ROI
if isfield(params,'ROIindex') && ~isempty(params.ROIindex)
    sigma = sigma(params.ROIindex);
    sigma2 = sigma2(params.ROIindex);
    beta1 = beta1(params.ROIindex);
    beta2 = beta2(params.ROIindex);
    x = x(params.ROIindex);
    y = y(params.ROIindex);
end

suppressionIndex = zeros(1,numel(sigma2));
pos_Vol = zeros(1,numel(sigma2));
neg_Vol = zeros(1,numel(sigma2));
pos_Vol_pRF = zeros(1,numel(sigma2));
neg_Vol_pRF = zeros(1,numel(sigma2));

for k =1:numel(sigma2)
    if ~params.SI_Centered
        rfpositive = rfGaussian2d(X,Y,sigma(k),sigma(k),0,x(k),y(k));   % original position of the pRF in the stimulus window
        rfnegative = rfGaussian2d(X,Y,sigma2(k),sigma2(k),0,x(k),y(k));
    else
        rfpositive = rfGaussian2d(X,Y,sigma(k),sigma(k),0,0,0);         % center the pRF at (0,0)
        rfnegative = rfGaussian2d(X,Y,sigma2(k),sigma2(k),0,0,0);
    end
    rfposGaus = beta1(k).*rfpositive;
    rfnegGaus = beta2(k).*rfnegative;
    rftotal = rfposGaus + rfnegGaus;

    % calculate the volume of both gaussians that make the pRF
    pos_Vol(k) = (sum(rfposGaus)).*gridsize;
    neg_Vol(k) = -(sum(rfnegGaus)).*gridsize;
    ipos = rftotal > 0;
    ineg = rftotal < 0;
    % calculate the volume of the positive and negative part of the pRF
    pos_Vol_pRF(k) = (sum(rftotal(ipos))).*gridsize;
    neg_Vol_pRF(k) = -(sum(rftotal(ineg))).*gridsize;

    if beta2(k)>= 0;
        suppressionIndex(k) = 0;                % if beta2 is bigger or equal to zero, there is no suppression and the suppression index will be 0
    else
        suppressionIndex(k) = neg_Vol(k)./pos_Vol(k);
    end

end

function params = rmProcessVarargin(params,vararg)
if ~exist('vararg','var') || isempty(vararg), return; end
for n=1:2:numel(vararg),
    data = vararg{n+1};
    fprintf(1,' %s,',vararg{n});
    switch lower(vararg{n}),
        case {'ri','roi','roiindex', 'roi_index','roi index'}
            params.ROIindex = data;

        case {'sts','stimsize','size of the stimulus'}
            params.stimSize = data;

        case {'rm','rm params','params struct of the retinotopic model'}    %vw.rm
            params.stimSize = data.retinotopyParams.stim.stimSize;

        case {'sr','samplerate','samplerate of the grid'}
            params.sampleRate = data;

        case {'sic','si_centered','si centered','suppression index centered'}
            params.SI_Centered = logical(data);

        otherwise,
            fprintf(1,'[%s]:IGNORING unknown parameter: %s\n',...
                mfilename,vararg{n});
    end;
end;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft