text:
function [SR, report] = mapsr(LRImages, model, varargin)

    if nargin > 2
        solverParams = varargin{end};
    else
        % Use default solver parameters
        solverParams = SRSolverParams;
    end

    % Assemble equation system to be solved consisting of low-resolution
    % observations and the system matrix.
    lrDim = size(LRImages(:,:,1));
    if solverParams.verbose
        disp('Assemble equation system (system matrix)...');
    end
    [W, LRImages] = composeSREquationSystem(LRImages, model);

    if isempty(model.SR)
        % Compute initial guess for the super-resolved image. We use the
        % "average image" estimated from the system matrix and photometric
        % parameters (if available).
        if solverParams.verbose
            disp('Compute average image used as initial guess...');
        end
        SR = getInitialSRImage(W, LRImages, model.photometricParams);
    else
        % Use initial guess provided by the user and reshape 2-D image into
        % parameter vector.
        SR = imageToVector(model.SR);
    end

    scgOptions = setupSCGOptions(solverParams);
    if solverParams.verbose
        disp('Minimize objective function...');
    end
    [SR, ~, flog, ~, ~]  = scg(@mapfunc, SR', scgOptions, @mapfunc_grad, model, LRImages, W, W');
    if nargout > 1
        report.numFunEvals = length(flog);
    end

    % Reshape parameter vector to a 2D image.
    SR = vectorToImage(SR, model.magFactor * lrDim);
    if solverParams.verbose
        disp('DONE!');
    end

function scgOptions = setupSCGOptions(solverParams)

    scgOptions = zeros(1,18);
    scgOptions(1) = 0;
    scgOptions(2) = solverParams.tolX;
    scgOptions(3) = solverParams.tolF;
    scgOptions(9) = solverParams.gradCheck;
    scgOptions(10) = solverParams.maxFunEvals;
    scgOptions(14) = solverParams.maxIter;

function f = mapfunc(SR, model, LR, W, ~)

    if ~iscolumn(SR)
        % Reshape to column vector.
        SR = SR';
    end

    % Evaluate the data fidelity term.
    dataTerm = mapDataTerm(SR, model, LR, W);

    % Evaluate image prior for regularization the super-resolved estimate.
    prior = model.imagePrior.function(SR, model.imagePrior.parameters{1:end});

    % Calculate objective function.
    f = dataTerm + model.imagePrior.weight * prior;

function grad = mapfunc_grad(SR, model, LR, W, Wt)

    if ~iscolumn(SR)
        % Reshape to column vector.
        SR = SR';
    end

    % Calculate gradient of the data fidelity term w.r.t. the
    % super-resolved image.
    dataTerm_grad = mapDataTerm_gradImage(SR, model, LR, W, Wt);

    % Calculate gradient of the regularization term w.r.t. the
    % super-resolved image.
    prior_grad = model.imagePrior.gradient(SR, model.imagePrior.parameters{1:end});

    % Sum up to total gradient
    grad = dataTerm_grad + model.imagePrior.weight * prior_grad;
    grad = grad';

author: thomas-koehler

repo: SupER


text:
function [x,state] = struct_diag(z,task)
%STRUCT_DIAG Diagonal matrix.
%   [x,state] = struct_diag(z) generates x as a diagonal matrix with
%   the elements of the vector z on the diagonal. The structure state
%   stores information which is reused in computing the right and left
%   Jacobian-vector products.
%
%   struct_diag(z,task) computes the right or left Jacobian-vector product
%   of this transformation, depending on the structure task. Use the
%   structure state and add the field 'r' of the same shape as z or the
%   field 'l' of the same shape as x to obtain the structure task for
%   computing the right and left Jacobian-vector products
%
%      (dF(:)/dz(:).')*task.r(:) and
%      (dF(:)/dz(:).')'*task.l(:) + conj((dF(:)/dconj(z(:)).')'*task.l(:)),
%
%   respectively. Here, F(z) represents this transormation, (:) signifies
%   vectorization and the derivative w.r.t. z (conj(z)) is a partial
%   derivative which treats conj(z) (z) as constant. The output has the
%   same shape as x or z for the right and left Jacobian-vector products,
%   respectively.
%
%   See also struct_band, struct_tridiag, struct_tril, struct_triu.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Structured data fusion,"
%       ESAT-SISTA Internal Report 13-177, KU Leuven, 2013.

if nargin < 2, task = []; end
state = [];

if isempty(task) || (isempty(task.l) && isempty(task.r))
    x = diag(z);
elseif ~isempty(task.r)
    x = diag(task.r);
elseif ~isempty(task.l)
    x = diag(task.l);
end

end

author: andrewssobral

repo: lrslibrary


text:
function polygon_solid_angle_3d_test ( )

%*****************************************************************************80
%
%% POLYGON_SOLID_ANGLE_3D_TEST tests POLYGON_SOLID_ANGLE_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 May 2015
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;
  test_num = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_SOLID_ANGLE_3D_TEST\n' );
  fprintf ( 1, '  POLYGON_SOLID_ANGLE_3D computes the solid angle\n' );
  fprintf ( 1, '  subtended by a planar polygon in 3D as viewed from\n' );
  fprintf ( 1, '  a point P.\n' );

  for test = 1 : test_num
%
%  One eighth of sphere surface, on the unit sphere surface.
%
    if ( test == 1 )

      n = 3;

      v = [ ...
        1.0, 0.0, 0.0; ...
        0.0, 1.0, 0.0; ...
        0.0, 0.0, 1.0 ]';

      p(1:3,1) = [ 0.0; 0.0; 0.0 ];
%
%  Reverse order of vertices.
%
    elseif ( test == 2 )

      n = 3;

      v = [ ...
        1.0, 0.0, 0.0; ...
        0.0, 0.0, 1.0; ...
        0.0, 1.0, 0.0 ]';

      p(1:3,1) = [ 0.0; 0.0; 0.0 ];
%
%  One eighth of sphere surface, on the unit sphere surface,
%  translated by (1,2,3).
%
    elseif ( test == 3 )

      n = 3;

      v = [ ...
        2.0, 2.0, 3.0; ...
        1.0, 3.0, 3.0; ...
        1.0, 2.0, 4.0 ]';

      p(1:3,1) = [ 1.0; 2.0; 3.0 ];
%
%  One eighth of sphere surface, but on sphere of radius 2.
%
    elseif ( test == 4 )

      n = 3;

      v = [ ...
        2.0, 0.0, 0.0; ...
        0.0, 2.0, 0.0; ...
        0.0, 0.0, 2.0 ]';

      p(1:3,1) = [ 0.0; 0.0; 0.0 ];

    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  TEST # %d\n', test );
    fprintf ( 1, '\n' );

    r8vec_print ( dim_num, p, '  The viewing point P:' );

    r8mat_transpose_print ( dim_num, n, v, '  The polygon vertices V:' );

    solid_angle = polygon_solid_angle_3d ( n, v, p );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Solid angle subtended: %f\n', solid_angle );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function b = dppsl ( ap, n, b )

%*****************************************************************************80
%
%% DPPSL solves a real symmetric positive definite system factored by DPPCO or DPPFA.
%
%  Discussion:
%
%    To compute inverse(A) * C where C is a matrix with P columns
%
%      call dppco ( ap, n, rcond, z, info )
%
%      if ( rcond is too small .or. info /= 0 ) then
%        exit
%      end if
%
%      do j = 1, p
%        call dppsl ( ap, n, c(1,j) )
%      end do
%
%    A division by zero will occur if the input factor contains
%    a zero on the diagonal.  Technically this indicates
%    singularity but it is usually caused by improper subroutine
%    arguments.  It will not occur if the subroutines are called
%    correctly and INFO == 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real AP(N*(N+1)/2), the output from DPPCO or DPPFA.
%
%    Input, integer N, the order of the matrix.
%
%    Input, real B(N), the right hand side.
%
%    Output, real B(N), the solution.
%
  kk = 0;

  for k = 1 : n
    t = ddot ( k-1, ap(kk+1:kk+k-1), 1, b(1:k-1), 1 );
    kk = kk + k;
    b(k) = ( b(k) - t ) / ap(kk);
  end

  for k = n : -1 : 1
    b(k) = b(k) / ap(kk);
    kk = kk - k;
    t = -b(k);
    b(1:k-1) = daxpy ( k-1, t, ap(kk+1:kk+k-1), 1, b(1:k-1), 1 );
  end

  return
end


text:
function B = vertcat(varargin)
%VERTCAT   Horizontally concatenate LINOP objects.

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

% Cast the result as a linop if there are any linop arguments.
B = linop( vertcat@chebmatrix(varargin{:}) );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef (Abstract) BirthModelX < BaseX
% BirthModelX Abstract class
%
% Summary of BirthModelX
%  This is the base class for all TrackingX birth models. Must be used as
%  the superclass of any custom defined TrackingX birth model.
%
% BirthModelX Properties:
%   + NumMeasDims - The number of observation dimensions.
%
% BirthModelX Methods:
%   ~ random - Process noise sample generator function
%   ~ pdf    - Function to evaluate the probability p(x_k|x_{k-1}) of
%              a set of new states, given a set of (particle) state vectors
%              e.g. eval = @(xk,xkm1) mvnpdf(xk,xkm1,Q);
%
% (+) denotes puplic properties/methods
% (~) denotes abstract properties/methods
%
% The above parameters and methods are accessed by the majority of
% existing TrackingX library components and are COMPULSORY to guarantee
% compatibility with the TrackingX framework.
%
% For examples see also the source of GenericBirthModelX

    properties
    end

    methods (Abstract)
        pdf(this);      % Function to evaluate the intensity of the birth model
        random(this);   % Generate iid samples from the birth model
    end

    methods
        function this = BirthModelX()
        % BirthModelX Constructor method
        %
        % See also GenericBirthModelX.

            this@BaseX();
        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function styleLabelAsValid(hLabel,validStr)
%styleLabelAsValid Summary of this function goes here
%   Detailed explanation goes here
    lblStr = strWrapForLabel(hLabel, validStr);

    set(hLabel,'Visible','on');
    set(hLabel,'BackgroundColor',[34,139,34]/255);
    set(hLabel,'ForegroundColor',[1,1,1]);
    set(hLabel,'String',lblStr);
    set(hLabel,'TooltipString','');
    set(hLabel,'FontWeight','bold');

    setappdata(hLabel,'tooltipstring',validStr);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function cvx_optval = sum_smallest( x, varargin )

%SUM_SMALLEST Sum of the smallest k elements of a vector.
%   For a real vector X and an integer k between 1 and length(X) inclusive,
%   y = SUM_SMALLEST(X,k) is the sum of the k smallest elements of X; e.g.,
%       temp = sort( x )
%       y = sum( temp( 1 : k ) )
%   If k=1, then SUM_SMALLEST(X,k) is equivalent to MIN(X); if k=length(X),
%   then SUM_SMALLEST(X,k) is equivalent to SUM(X).
%
%   Both X and k must be real, and k must be a scalar. But k is not, in
%   fact, constrained to be an integer between 1 and length(X); the
%   function is extended continuously and logically to all real k. For
%   example, if k <= 0, then SUM_SMALLEST(X,k)=0. If k > length(X), then
%   SUM_SMALLEST(X,k)=SUM(X). Non-integer values of k interpolate linearly
%   between their integral neighbors.
%
%   For matrices, SUM_SMALLEST(X,k) is a row vector containing the
%   application of SUM_SMALLEST to each column. For N-D arrays, the
%   SUM_SMALLEST operation is applied to the first non-singleton dimension
%   of X.
%
%   SUM_SMALLEST(X,k,DIM) performs the operation along dimension DIM of X.
%
%   Disciplined convex programming information:
%       SUM_SMALLEST(X,...) is concave and nondecreasing in X. Thus, when
%       used in CVX expressions, X must be concave (or affine). k and DIM
%       must both be constant.

narginchk(2,3);
cvx_optval = -sum_largest( -x, varargin{:} );

% 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:
% correct_mc() - compute an upper limit for the number of independant
%                time-frequency estimate in a given time-frequency image.
%                This number can be used to correct for multiple comparisons.
%
% Usage:
%   [ncorrect array] = correct_mc( EEG, cycles, maxfreq, timesout);
%
% Inputs:
%    EEG       - EEGLAB structure
%    cycles    - [float] same as the cycle input to timef(). Default is [3 0.5].
%    freqrange - [float] minimum and maximum frequency. Default is [2 50] Hz.
%    timesout  - [integer] array of number of time points to test.
%
% Output:
%    ncorrect - number of independant tf estimate in the time-freq image
%    array    - array of size (freqs x timesout) containing pvalues.
%
% Method details:
%
% Dividing by the total number of time-frequency estimate in the 2-D
% time-frequency image decomposition would be too conservative since
% spectral estimates of neighboring time-frequency points are highly
% correlated. One must thus estimate the number of independent
% time-frequency points in the TF image. Here, I used geometrical wavelets
% which are optimal in terms of image compression, so neighboring
% frequencies can be assume to carry independent spectral estimates.
% We thus had time-frequency decompositions at only X frequencies (e.g. 120,
% 60, 30, 15, 7.5, 3.25, 1.625 Hz). For each frequency, I then found
% the minimum number of time points for which there was a significant
% correlation of the spectral estimates between neighboring time points
% (for each frequency and number of time point, I computed the correlation
% from 0 to 1 for all data channel to obtain an a probability distribution
% of correlation; we then fitted this distribution using a 4th order curve
% (Ramberg, J. S., E. J. Dudewicz, et al. (1979). "A probability
% distribution and its uses in fitting data." Technometrics 21(2)) and
% assessed the probability of significance for the value 0 (no correlation)
% to be within the distribution of estimated correlation). For instance,
% using 28 time points at 120 Hz, there was no significant (p>0.05 taking
% into account Bonferoni correction for multiple comparisons) correlation
% between neighboring time-frequency power estimate, but there was a
% significant correlation using 32 time points instead of 28 (p<0.05).
% Applying the same approach for the X geometrical frequencies and summing
% the minimum number of time points for observing a significant correlation
% at all frequencies, ones obtain in general a number below 200 (with the
% defaults above and 3-second data epochs) independent estimates. In all
% the time-frequency plots, one has to used a significance mask at p<0.00025
% (0.05/200). An alternative method for correcting for multiple comparisons
% is presented in Nichols & Holmes, Human Brain Mapping, 2001.
%
% Author: Arnaud Delorme, SCCN, Jan 17, 2004

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

function [ncorrect, pval] = correct_mc( EEG, cycles, freqrange, timesout);

    if nargin < 1
        help correct_mc;
        return;
    end;
    if nargin < 2
        cycles  = [3 0.5];
    end;
    if nargin < 3
        freqrange = [2 50];
    end;
    if nargin < 4
        % possible number of time outputs
        % -------------------------------
        timesout = [5 6 7 8 9 10 12 14 16 18 20 24 28 32 36 40];
    end;
    nfreqs = ceil(log2(freqrange(2)));

    % scan times
    % ----------
    for ti = 1:length(timesout)
        clear tmpf

        % scan data channels
        % ------------------
        for index = 1:EEG.nbchan

            clf; [ersp,itc,powbase,times,freqs,erspboot,itcboot] = newtimef(EEG.data(index,:),EEG.pnts, ...
                             [EEG.xmin EEG.xmax]*1000,EEG.srate, cycles, 'timesout', timesout(ti), ...
                             'freqscale', 'log', 'nfreqs', nfreqs, 'freqrange', freqrange, 'plotitc', 'off', 'plotersp', 'off');

            % compute correlation
            % -------------------
            for fi = 1:length(freqs)
                tmp      = corrcoef(ersp(fi,1:end-1), ersp(fi,2:end));
                tmpf(index,fi) = tmp(2,1);
            end;

        end;

        % fit curve and determine if the result is significant
        % ----------------------------------------------------
        for fi = 1:length(freqs)
            pval(fi, ti) = rsfit(tmpf(:,fi)', 0);
            if pval(fi,ti) > 0.9999, pval(fi,ti) = NaN; end;
        end;
    end;

    % find minimum number of points for each frequency
    % ------------------------------------------------
    ncorrect = 0;
    threshold = 0.05 / prod(size(pval));
    for fi = 1:size(pval,1)
        ti = 1;
        while ti <= size(pval,2)
            if pval(fi,ti) < threshold
                ncorrect = ncorrect +  timesout(ti);
                ti = size(pval,2)+1;
            end;
            ti = ti+1;
        end;
    end;

author: PatternRecognition

repo: OpenBMI


text:
function bz_PreprocessExtracellEphysSession(basepath)

if ~exist('basepath','var')
    basepath = cd;
elseif isempty(basepath)
    basepath = cd;
end

%% Assuming one already did bz_SetAnimalMetadata
bz_RunSessionMetadata(basepath);

%% making an xml if needed
basename = bz_BasenameFromBasepath(basepath);
txmlname = fullfile(basepath,[basename '.xml']);
d = dir(txmlname);
load(fullfile(basepath,[basename '.SessionMetadata.mat']))
if isempty(d);%if no basename.xml already
    if isempty(SessionMetadata.ExtracellEphys.ParametersDivergentFromAnimalMetadata)%if not different from animal xml, copy the animal xml here
        ax = fullfile(SessionMetadata.AnimalMetadata.AnimalBasepath,[SessionMetadata.AnimalMetadata.AnimalName, '.xml']);
        if exist(ax,'file')
            copyfile(ax,txmlname);
        end
    end
end
d = dir(txmlname);
if isempty(d);%if no basename.xml already
    bz_MakeXML(basepath)
end

%%
disp('Concatenating .dat files')
bz_ConcatenateDats(basepath);
disp('Converting .dat to .lfp')
bz_LFPFromDat(basepath);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function value = r8vec_norm_lp ( n, a, p )

%*****************************************************************************80
%
%% R8VEC_NORM_LP returns the LP norm of an R8VEC.
%
%  Discussion:
%
%    The vector LP norm is defined as:
%
%      value = ( sum ( 1 <= I <= N ) ( abs ( A(I) ) )^P )^(1/P).
%
%    This routine allows P to have the special Matlab value Inf,
%    in which case the L-infinity norm is returned.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in A.
%
%    Input, real A(N), the vector whose LP norm is desired.
%
%    Input, real P, the index of the norm.
%
%    Output, real VALUE, the LP norm of A.
%
  if ( p <= 0.0 )
    value = -1.0;
  elseif ( p == Inf )
    value = maxval ( abs ( a(1:n) ) );
  elseif ( p == 1.0 )
    value = sum ( abs ( a(1:n) ) );
  elseif ( p == 2.0 )
    value = sqrt ( sum ( a(1:n).^2 ) );
  else
    value = ( sum ( ( abs ( a(1:n) ) ).^p ) ).^( 1.0 / p );
  end

  return
end


text:
classdef bndfun < classicfun
%BNDFUN   Represent global functions on a bounded interval [a, b].
%
%   Class for representing global functions on a bounded interval [a, b].
%   Functions are approximated via a ONEFUN object, which lives on the interval
%   [-1, 1], stored in the BNDFUN. Forward and inverse maps stored in the
%   BNDFUN object map the interval [-1, 1] to [a, b], and vice versa.
%
% Constructor inputs:
%   BNDFUN(OP) constructs a BNDFUN object from the function handle OP on the
%   domain determined by the default in CHEBFUNPREF by mapping the domain to
%   [-1, 1] and constructing a ONEFUN object to represent the function
%   prescribed by OP.  OP should be vectorised (i.e., accept a vector input)
%   and output a vector of the same length as its input.
%
%   BNDFUN(OP, DATA) constructs a BNDFUN object using the data supplied in the
%   DATA structure.  DATA fields used by BNDFUN are:
%     DATA.DOMAIN    (Default:  Determined by CHEBFUNPREF)
%         A row vector with two elements in increasing order defining the
%         construction domain.  Both elements must be finite.
%   In addition, BNDFUN may modify the following DATA fields before passing
%   them on to the ONEFUN constructor:
%     DATA.HSCALE    (Default:  Determined by DATA.DOMAIN)
%         Horizontal construction scale.
%   If any fields in DATA are empty or not supplied, or if DATA itself is empty
%   or not supplied, appropriate default values are set.  Any fields in DATA
%   which are not recognized will be passed as-is to the ONEFUN constructor.
%
%   BNDFUN(OP, DATA, PREF) overrides the default behavior with that given by
%   the preferences in the structure or CHEBFUNPREF object PREF. See
%   CHEBFUNPREF for details.
%
% See also CLASSICFUN, CHEBFUNPREF, ONEFUN.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BNDFUN Class Description:
%
% The BNDFUN class is a class for representations of globally smooth functions
% on a finite interval [a, b]. It achieves this by taking a ONEFUN on [-1, 1]
% and applying a linear mapping to re-scale its domain.
%
% Note that all binary BNDFUN operators (methods which can take two BNDFUN
% arguments) assume that the domains of the BNDFUN objects agree. The methods
% will not issue warnings if this condition is violated, but the results will
% not be meaningful.
%
% Class diagram: [<<CLASSICFUN>>] <>-- [<<onefun>>]
%                   ^
%                   |
%                [bndfun]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function obj = bndfun(op, data, pref)
            % Parse inputs.
            if ( (nargin < 1) || isempty(op) )
                obj.domain = [];
                obj.mapping = [];
                obj.onefun = [];
                return
            end

            if ( (nargin < 2) || isempty(data) )
                data = struct();
            end

            if ( (nargin < 3) || isempty(pref) )
                pref = chebfunpref();
            else
                pref = chebfunpref(pref);
            end

            data = parseDataInputs(data, pref);

            % Check the domain input.
            if ( ~all(size(data.domain) == [1, 2]) || (diff(data.domain) <= 0) )
                error('CHEBFUN:BNDFUN:bndfun:badDomain', ...
                    ['Domain argument should be a row vector with two ', ...
                    'entries in increasing order.']);
            elseif ( any(isinf(data.domain)) )
                error('CHEBFUN:BNDFUN:bndfun:unboundedDomain', ...
                    'Should not encounter unbounded domain in bndfun class.');
            end

            % TODO:  Why do we rescale the hscale like this?
            data.hscale = data.hscale / diff(data.domain);

            % Remap the OP to be a function on [-1, 1].
            linmap = bndfun.createMap(data.domain);
            if ( isa(op, 'function_handle') && ~all(data.domain == [-1, 1]) ...
                    && ~isnumeric(op) )
                op = @(x) op(linmap.For(x));
            end

            % Call the ONEFUN constructor:
            obj.onefun = onefun.constructor(op, data, pref);

            % Add the domain and mapping:
            obj.domain = data.domain;
            obj.mapping = linmap;
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% STATIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = true )

        % Linear map from [-1, 1] to the domain of the BNDFUN.
        m = createMap(domain);

        % Make a BNDFUN (constructor shortcut):
        f = make(varargin);

    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% METHODS IMPLEMENTED IN THIS FILE:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function data = parseDataInputs(data, pref)
%PARSEDATAINPUTS   Parse inputs from the DATA structure and assign defaults.

if ( ~isfield(data, 'domain') || isempty(data.domain) )
    data.domain = pref.domain;
end

if ( ~isfield(data, 'hscale') || isempty(data.hscale) )
    % TODO:  Or should this be 1?  What does chebfun pass down?
    data.hscale = norm(data.domain, inf);
end

end


text:
function [constants,general,singles,pairs] = classifyfactors(L,M,R)

nfac = length(L);
used = zeros(1,nfac);

constants = 0;
general = 0;
singles = {};
pairs = {};
for i = 1:length(L)
    if isa(M{i},'double')
        constants = constants + L{i}*M{i}*R{i};
        used(i)=1;
    end
end

for i = 1:length(L)
    if ~used(i)
        temp = struct(M{i});
        if ~(is(M{i},'sdpcone') || is(M{i},'lpcone') || isequal(temp.originalbasis,'skew') || isequal(temp.originalbasis,'diagonal'))
            general = general + L{i}*M{i}*R{i};
            used(i) = 1;
        elseif size(M{i},2)~=size(R{i},1) | size(M{i},1)~=size(L{i},2)
            % Something like scalar*matrix
            % simplifies code if we treat that as a general
            general = general + L{i}*M{i}*R{i};
            used(i) = 1;
        end
    end
end

% sameL = [];
% sameR = [];
% for i = 1:length(L)-1
%     for j = i+1:length(L)
%         if isequal(M{i},M{j})
%             if isequal(L{i},L{j})
%                 mmm(i,j) = 1;
%                 mmm(j,i) = 1;
%                 sameL = [sameL;i j];
%             elseif isequal(R{i},R{j})
%                 sameR = [sameR;i j];
%                 mmm(i,j) = 1;
%                 mmm(j,i) = 1;
%             end
%         end
%     end
% end
% if length(sameL)>1
%     for i = 1:size(sameL,1)-1
%         for j = 2:size(sameL,2)
%             AA = sameL([i j],:);
%             BB = sameR([i j],:);
%             if AA == BB'
%                 candidate = unique([AA(:);BB(:)]);
%             end
%         end
%     end
% end

for i = 1:length(L)

     if issymmetric(M{i}) && isequal(find(L{i}),find(R{i}'))
        Rttmp = R{i}';
        tmp = L{i}(L{i}~=0)./(Rttmp(Rttmp~=0));
     else
         tmp = [];
     end

    if ~used(i)
        if issymmetric(M{i}) && (isequal(L{i},R{i}'))
            singles{end+1}.L = L{i};
            singles{end}.M = M{i};
            singles{end}.R = R{i};
            singles{end}.negated = 0;
            used(i)=1;
        elseif issymmetric(M{i}) && (isequal(L{i},-R{i}'))
            singles{end+1}.L = L{i};
            singles{end}.M = M{i};
            singles{end}.R = R{i};
            singles{end}.negated = 1;
            used(i)=1;
        elseif issymmetric(M{i}) && isequal(find(L{i}),find(R{i}')) &  length(unique(tmp)) == 1
        %    Rttmp = R{i}';
        %    tmp = L{i}(L{i}~=0)./(Rttmp(Rttmp~=0));

         %   if length(unique(tmp)) == 1
                  singles{end+1}.L = L{i}/sqrt(abs(tmp(1)));
                  singles{end}.M = M{i};
                  singles{end}.R = R{i}*sqrt(abs(tmp(1)));
                  if tmp(1) > 0
                      singles{end}.negated = 0;
                  else
                      singles{end}.negated = 1;

                  end
                  used(i)=1;
        %  end

        else
            for j = i+1:length(L)
                if isa(L{i}*M{i}*R{i}-R{j}'*M{j}'*L{j}','double')
                    pairs{end+1}.L = L{i};
                    pairs{end}.M = M{i};
                    pairs{end}.R = R{i};
                    used(j)=1;
                    used(i)=1;
                    break
                end
            end
        end
    end
end

for k = find(~used)
    if isequal(size(M{k}),[1 1])
        general = general + L{k}*M{k}*R{k};
    else
        error('Couldn''t classify all factors');
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function beep2(w,t)
%plays a short tone as an audible cue
%
%USAGE:
%    beep2
%    beep2(w)    specify frequency (200-1,000 Hz)
%    beep2(w,t)     "       "       and duration in seconds

fs=8192;  %sample freq in Hz

if (nargin == 0)
   w=1000;            %default
   t = [0:1/fs:.2];   %default
elseif (nargin == 1)
   t = [0:1/fs:.2];   %default
elseif (nargin == 2)
   t = [0:1/fs:t];
end

%one possible wave form
wave=sin(2*pi*w*t);

%play sound
sound(wave,fs);

author: Sable

repo: mcbench-benchmarks


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%	SurfBox-MATLAB (c)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%	Yue M. Lu and Minh N. Do
%%
%%	Department of Electrical and Computer Engineering
%%	Coordinated Science Laboratory
%%	University of Illinois at Urbana-Champaign
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%	surf_vec2coeff.m
%%
%%	First created: 04-06-07
%%	Last modified: 04-06-07
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Y = surf_coeff2vec(c, size_info)

%% Convert the surfacelet coefficient stored in a linear vector to a nested
%% cell array.
%%
%% Input:
%%
%% c: a linear vector storing all the coefficients.
%%
%% size_info: a book keeping cell array storing the dimension information
%% of Y. This is useful when we want to convert c back to Y. See
%% surf_vec2coeff.m for details.
%%
%% Output:
%%
%% Y: an L+1 by 1 cell array containing the surfacelet coefficients. See
%% Surfdec.m for details
%%

%% we can inherit the nested structure
Y = size_info;

%% number of multiscale levels
L = length(Y) - 1;

Ncoeff = 1;

for n = 1 : L
    for m = 1 : length(Y{n})
        for k = 1 : length(Y{n}{m})
            sz = size_info{n}{m}{k};
            d = prod(sz);
            Y{n}{m}{k} = reshape(c(Ncoeff : Ncoeff + d - 1), sz);
            Ncoeff = Ncoeff + d;
        end
    end
end

%% the lowpass subband
sz = size_info{end};
Y{end} = reshape(c(Ncoeff : end), sz);

author: thomaskuestner


text:
function [data,data2,Isymbols,Qsymbols]=serial_to_parallel(N)
%Bob Gess, June 2008 for EE 473
%Extracted and modified from original code (QPSK_TX_IQ_RX) written by JC and posted on
%Mathworks 12-23-2005.

%This module generates the data and splits it into the I channel and Q
%channel
%N=1e3;		        % Number of data bits(bit rate)
fs=40*2e3;		    % Sampling frequency
Fn=fs/2;            % Nyquist frequency
Ts=1/fs;	        % Sampling time = 1/fs
T=1/N;		        % Bit time
randn('state',0);   % Keeps PRBS from changing on reruns
td=[0:Ts:(N*T)-Ts]';% Time vector(data)(transpose)
%===================================================================
% Input data
%===================================================================
data=sign(randn(N,1))';%transpose
data1=ones(T/Ts,1)*data;
data2=data1(:);
%length(data2)

%display input data bits in command window
data_2=data2';
data_2=data_2 >0;
x=0;
transmitted_data_bits=data_2(1:(fs+x)/N:end);

%Serial to parallel (alternating)
tiq = [0:Ts*2:(N*T)-Ts]';% Time vector for I and Q symbols(transpose)

bs1=data(1:2:length(data));%odd
symbols=ones(T/Ts,1)*bs1;
Isymbols=symbols(:);%I_waveform
%LI=length(Isymbols)

bs2=data(2:2:length(data));%even
symbols1=ones(T/Ts,1)*bs2;
Qsymbols=symbols1(:);%Q_waveform
%LQ=length(Qsymbols)

%=====================================================================
%Plots
%======================================================================
figure(1)
subplot(3,1,1)
plot(td,data2)
axis([0 20/N -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('Input Data')

subplot(3,1,2)
plot(tiq,Isymbols)
axis([0 20/N -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('I Channel(one bit/symbol(phase)) Data')

subplot(3,1,3)
plot(tiq,Qsymbols)
axis([0 20/N -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('Q Channel(one bit/symbol(phase)) Data')


text:
% -----------------------------------------------------------------  %
% Matlab Programs included the Appendix B in the book:               %
%  Xin-She Yang, Engineering Optimization: An Introduction           %
%                with Metaheuristic Applications                     %
%  Published by John Wiley & Sons, USA, July 2010                    %
%  ISBN: 978-0-470-58246-6,   Hardcover, 347 pages                   %
% -----------------------------------------------------------------  %
% Citation detail:                                                   %
% X.-S. Yang, Engineering Optimization: An Introduction with         %
% Metaheuristic Application, Wiley, USA, (2010).                     %
%                                                                    %
% http://www.wiley.com/WileyCDA/WileyTitle/productCd-0470582464.html %
% http://eu.wiley.com/WileyCDA/WileyTitle/productCd-0470582464.html  %
% -----------------------------------------------------------------  %
% ===== ftp://  ===== ftp://   ===== ftp:// =======================  %
% Matlab files ftp site at Wiley                                     %
% ftp://ftp.wiley.com/public/sci_tech_med/engineering_optimization   %
% ----------------------------------------------------------------   %

% Simulated Annealing for constrained optimization
% by Xin-She Yang @ Cambridge University @2008
% Usage: sa_mincon(alpha)

function [bestsol,fmin,N]=sa_mincon(alpha)

% Default cooling factor
if nargin<1,
    alpha=0.95;
end

% Display usage
disp('sa_mincon or [Best,fmin,N]=sa_mincon(0.95)');

% d dimensions

% Welded beam design optimization
Lb=[0.1 0.1  0.1 0.1];
Ub=[2.0 10.0 10.0 2.0];
u0=(Lb+Ub)/2;

if length(Lb) ~=length(Ub),
    disp('Simple bounds/limits are improper!');
    return
end

%% Start of the main program -----------------------------------------
d=length(Lb);       % Dimension of the problem

% Initializing parameters and settings
T_init = 1.0;       % Initial temperature
T_min =  1e-10;     % Finial stopping temperature
F_min = -1e+100;    % Min value of the function
max_rej=500;        % Maximum number of rejections
max_run=150;        % Maximum number of runs
max_accept = 50;    % Maximum number of accept
initial_search=500; % Initial search period
k = 1;              % Boltzmann constant
Enorm=1e-5;         % Energy norm (eg, Enorm=1e-8)

% Initializing the counters i,j etc
i= 0; j = 0; accept = 0; totaleval = 0;
% Initializing various values
T = T_init;
E_init = Fun(u0);
E_old = E_init; E_new=E_old;
best=u0;  % initially guessed values
% Starting the simulated annealling
while ((T > T_min) & (j <= max_rej) & E_new>F_min)
    i = i+1;
    % Check if max numbers of run/accept are met
    if (i >= max_run) | (accept >= max_accept)
        % reset the counters
        i = 1;  accept = 1;
      % Cooling according to a cooling schedule
        T = cooling(alpha,T);
    end

    % Function evaluations at new locations
    if totaleval<initial_search,
        init_flag=1;
        ns=newsolution(u0,Lb,Ub,init_flag);
    else
        init_flag=0;
        ns=newsolution(best,Lb,Ub,init_flag);
    end

      totaleval=totaleval+1;
      E_new = Fun(ns);
    % Decide to accept the new solution
    DeltaE=E_new-E_old;
    % Accept if improved
    if (DeltaE <0)
        best = ns; E_old = E_new;
        accept=accept+1;   j = 0;
    end
    % Accept with a small probability if not improved
    if (DeltaE>=0 & exp(-DeltaE/(k*T))>rand );
        best = ns; E_old = E_new;
        accept=accept+1;
    else
        j=j+1;
    end
    % Update the estimated optimal solution
    f_opt=E_old;
end

bestsol=best
fmin=f_opt
N=totaleval

%% New solutions
function s=newsolution(u0,Lb,Ub,init_flag)
  % Either search around
if length(Lb)>0 & init_flag==1,
  s=Lb+(Ub-Lb).*rand(size(u0));
else
  % Or local search by random walk
  s=u0+0.01*(Ub-Lb).*randn(size(u0));
end

s=bounds(s,Lb,Ub);

%% Cooling
function T=cooling(alpha,T)
T=alpha*T;

function ns=bounds(ns,Lb,Ub)
if length(Lb)>0,
% Apply the lower bound
  ns_tmp=ns;
  I=ns_tmp<Lb;
  ns_tmp(I)=Lb(I);
  % Apply the upper bounds
  J=ns_tmp>Ub;
  ns_tmp(J)=Ub(J);
% Update this new move
  ns=ns_tmp;
else
  ns=ns;
end

% d-dimensional objective function
function z=Fun(u)

% Objective
z=fobj(u);

% Apply nonlinear constraints by penalty method
% Z=f+sum_k=1^N lam_k g_k^2 *H(g_k)
z=z+getnonlinear(u);

function Z=getnonlinear(u)
Z=0;
% Penalty constatn
lam=10^15; lameq=10^15;
[g,geq]=constraints(u);

% Inequality constraints
for k=1:length(g),
    Z=Z+ lam*g(k)^2*getH(g(k));
end

% Equality constraints (when geq=[], length->0)
for k=1:length(geq),
   Z=Z+lameq*geq(k)^2*geteqH(geq(k));
end

% Test if inequalities hold
function H=getH(g)
if g<=0,
    H=0;
else
    H=1;
end

% Test if equalities hold
function H=geteqH(g)
if g==0,
    H=0;
else
    H=1;
end

% Objective functions
function z=fobj(u)
% Welded Beam Design Optimization
% K. Ragsdell and D. Phillips, Optimal design of a class of welded
% strucures using geometric programming,
% J. Eng. Ind., 98 (3):1021-1025, (1976).
% Best solution found in literature
% [0.205730, 3.470489, 9.036624, 0.205729) with the objective 1.724852
% by Cagnina et al., Solving engineering optimization problems with the
% simple constrained particle swarm optimizer, Informatica, 32 (2008)319-326

z=1.10471*u(1)^2*u(2)+0.04811*u(3)*u(4)*(14.0+u(2));

% For Rosenbrock's function
% z=(1-u(1))^2+(u(2)-u(1)^2)^2+(1-u(3))^2+(1-u(4))^2

% All constraints
function [g,geq]=constraints(x)
% Inequality constraints
Q=6000*(14+x(2)/2);
D=sqrt(x(2)^2/4+(x(1)+x(3))^2/4);
J=2*(x(1)*x(2)*sqrt(2)*(x(2)^2/12+(x(1)+x(3))^2/4));
alpha=6000/(sqrt(2)*x(1)*x(2));
beta=Q*D/J;
tau=sqrt(alpha^2+2*alpha*beta*x(2)/(2*D)+beta^2);
sigma=504000/(x(4)*x(3)^2);
delta=65856000/(30*10^6*x(4)*x(3)^3);
F=4.013*(30*10^6)/196*sqrt(x(3)^2*x(4)^6/36)*(1-x(3)*sqrt(30/48)/28);

g(1)=tau-13600;
g(2)=sigma-30000;
g(3)=x(1)-x(4);
g(4)=0.10471*x(1)^2+0.04811*x(3)*x(4)*(14+x(2))-5.0;
g(5)=0.125-x(1);
g(6)=delta-0.25;
g(7)=6000-F;

% Equality constraints
geq=[];
%% End of the program ------------------------------------------------


text:
function [V,F] = simple_cube()
%SIMPLE_CUBE  construct a triangle mesh for a simple cube
%
% [V,F] = simple_cube();
%
% Outputs:
%  V,F  the triangle mesh of a simple cube

V = [?, ?, ?; ... %dots in MATLAB allow us to insert a linebreak, but continue with the same statement
    ?, ?, ?; ...
    ?, ?, ?; ...
    ?, ?, ?; ...
    ?, ?, ?; ...
    ?, ?, ?; ...
    ?, ?, ?; ...
    ?, ?, ?];

F = [  ];

end

author: odedstein

repo: sgi-introduction-course


text:
% Demo for Edge Boxes (please see readme.txt first).

addpath(genpath(pwd)); savepath;
addpath(genpath('/home/priyanka/Documents/autonomous_systems/sem-4/edgeBoxes/toolbox/')); savepath;

%% load pre-trained edge detection model and set opts (see edgesDemo.m)
% model=load('models/forest/modelBsds'); model=model.model;
model=load('models/forest/modelNyud2Rgbd'); model=model.model;
model.opts.multiscale=0; model.opts.sharpen=2; model.opts.nThreads=4;

%% set up opts for edgeBoxes (see edgeBoxes.m)
opts = edgeBoxes;
opts.alpha = .65;     % step size of sliding window search
opts.beta  = .75;     % nms threshold for object proposals
opts.minScore = .01;  % min score of boxes to detect
opts.maxBoxes = 1e4;  % max number of boxes to detect

%dataDir='/home/priyanka/Documents/autonomous_systems/sem-4/edgeBoxes/pdollar-edges-94260b5/BSR/nyud2_dataset/data/';
dataDir='/home/priyanka/Documents/autonomous_systems/sem-4/edgeBoxes/pdollar-edges-94260b5/BSR/vocb3do/data/';
flist=dir([dataDir '/images/train/*.png']); flist={flist.name};

cnt=0;
trainfileID = fopen('train_annot.txt','a');
testfileID = fopen('test_annot.txt','a');
trainFileName=importdata('train.txt');
testFileName=importdata('val.txt');
recall=0;
t_imgs=0;
num_prop=2000;
numel(flist)
for fidx=1:numel(flist)
	tp=0;
	tpfn=0;
	I=single(imread([dataDir 'images/train/' flist{fidx}]))/255;
	%D=single(imread([dataDir 'depth/train/' flist{fidx}]))/1e4;
	D=single(imread([dataDir 'depth/train/' flist{fidx}(1:end-4) '_abs_smooth.png']))/1e4;
	tic, bbs=edgeBoxes(I,D,model,opts); toc
	display(sprintf('img %d of %d, number of proposals: %d',fidx,numel(flist),size(bbs,1)));
	%if(size(bbs,1) > num_prop)
	%	bbs=bbs(1:num_prop,:);
	%end
	%num_prop=size(bbs,1);
	%t_imgs=numel(flist);
	%%% detect Edge Box bounding box proposals (see edgeBoxes.m)
	%I = imread('peppers.png');
	%tic, bbs=edgeBoxes(I,model,opts); toc

	%%% show evaluation results (using pre-defined or interactive boxes)
	%gt=[122 248 92 65; 193 82 71 53; 410 237 101 81; 204 160 114 95; ...
	%  9 185 86 90; 389 93 120 117; 253 103 107 57; 81 140 91 63];
	try
		display('in try')
		fName=['/home/priyanka/Desktop/VOCB3DO/VOCB3DO_Annotations/all_objects/' flist{fidx}(1:end-3) 'txt'];
		%fName=['/home/priyanka/Desktop/VOCB3DO/VOCB3DO_Annotations/all_objects/img_0013.txt'];
		fileID = fopen(fName);
		C = textscan(fileID,'%d %d %d %d %s');
		gt=cat(2,C{1},C{2},C{3},C{4});
		class_name=C{5};
		fclose(fileID);
		t_imgs=t_imgs+1;
	catch
		display('in catch')
		continue;
	end
	if(0), gt='Please select an object box.'; disp(gt); figure(1); imshow(I);
	  title(gt); [~,gt]=imRectRot('rotate',0); gt=gt.getPos(); end
	gt(:,5)=0; [gtRes,dtRes]=bbGt('evalRes',gt,double(bbs),.7);
	figure(1); bbGt('showRes',I,gtRes,dtRes(dtRes(:,6)==1,:)); pause(0.1);
	title('green=matched gt  red=missed gt  dashed-green=matched detect');
	A=bbGt('showRes',I,gtRes,dtRes(dtRes(:,6)==1,:));

	iname=['imgs/' flist{fidx}(1:end-4) '.png'];
	saveas(gcf,iname);

	% getting object bounding boxes
	%fName=replace(fName,'/home/priyanka/Desktop/VOCB3DO/VOCB3DO_Annotations/all_objects/','');
	fName=strrep(fName,'/home/priyanka/Desktop/VOCB3DO/VOCB3DO_Annotations/all_objects/','');

	isTrain=any(ismember(trainFileName,flist{fidx}(1:end-4)));
	isTest=any(ismember(testFileName,flist{fidx}(1:end-4)));
	if(isTrain == 1)
		depthDir='/home/priyanka/Documents/autonomous_systems/sem-4/edgeBoxes/pdollar-edges-94260b5/b3do/train/depth/';
		rgbDir='/home/priyanka/Documents/autonomous_systems/sem-4/edgeBoxes/pdollar-edges-94260b5/b3do/train/rgb/';
		writefileID=trainfileID;
	end
	if(isTest == 1)
		depthDir='/home/priyanka/Documents/autonomous_systems/sem-4/edgeBoxes/pdollar-edges-94260b5/b3do/test/depth/';
		rgbDir='/home/priyanka/Documents/autonomous_systems/sem-4/edgeBoxes/pdollar-edges-94260b5/b3do/test/rgb/';
		writefileID=testfileID;
	end

	origD=imread([dataDir 'depth/train/' flist{fidx}(1:end-4) '_abs_smooth.png']);
	for gidx=1:size(gt,1)
		cnt=cnt+1;
		tempgt = gt(gidx,:);
		[gtResT,dtResT]=bbGt('evalRes',tempgt,double(bbs),.7);
		loc = find(dtResT(:,6)==1);
		if(~isempty(loc))
			oname = class_name{gidx};
			I2 = imcrop(I,gtResT(1:4));
			D2 = imcrop(origD,gtResT(1:4));
			r_name=[oname '_1_1_' num2str(cnt) '_crop.png'];
			d_name=[oname '_1_1_' num2str(cnt) '_depthcrop.png'];
			imwrite(I2,[rgbDir '' r_name]);
			imwrite(D2,[depthDir '' d_name]);
			fprintf(writefileID,'%s %s\n',r_name,num2str(lookup(oname)));

		end
	end

	loc = find(dtRes(:,6)==1);
	tp = tp + numel(loc);
	tpfn = tpfn+size(gtRes,1);
	%display(tp);
	%display('...');
	%display(dtRes(loc,:));
	%display(tpfn);
	%display('---');
	%disp(tp);
	%disp(tpfn);
	recall=recall+(tp/tpfn);
end

display(recall);
display(recall/t_imgs);
if 0
%% run and evaluate on entire dataset (see boxesData.m and boxesEval.m)
if(~exist('boxes/VOCdevkit/','dir')), return; end
split='val'; data=boxesData('split',split);
nm='EdgeBoxes70'; opts.name=['boxes/' nm '-' split '.mat'];
edgeBoxes(data.imgs,model,opts); opts.name=[];
boxesEval('data',data,'names',nm,'thrs',.7,'show',2);
boxesEval('data',data,'names',nm,'thrs',.5:.05:1,'cnts',1000,'show',3);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMat(data,win,params)
%
%
% Multi-taper cross-spectral matrix - another routine, this one allows for multiple trials and channels
% but does not do confidence intervals. Also this routine always averages over trials - continuous process
%
% Usage:
%
% [Sc,Cmat,Ctot,Cvec,Cent,f]=CrossSpecMat(data,win,params)
% Input:
% Note units have to be consistent. See chronux.m for more information.
%       data (in form samples x channels x trials)
%       win  (duration of non-overlapping window)
%       params: structure with fields tapers, pad, Fs, fpass
%       - optional
%           tapers (precalculated tapers from dpss, or in the form [NW K] e.g [3 5]) -- optional.
%                                                 If not specified, use [NW K]=[3 5]
%	        pad		    (padding factor for the FFT) - optional. Defaults to 0.
%			      	 e.g. For N = 500, if PAD = 0, we pad the FFT
%			      	 to 512 points; if PAD = 2, we pad the FFT
%			      	 to 2048 points, etc.
%           Fs   (sampling frequency) - optional. Default 1.
%           fpass    (frequency band to be used in the calculation in the form
%                                   [fmin fmax])- optional.
%                                   Default all frequencies between 0 and Fs/2
% Output:
%       Sc (cross spectral matrix frequency x channels x channels)
%       Cmat Coherence matrix frequency x channels x channels
%       Ctot Total coherence: SV(1)^2/sum(SV^2) (frequency)
%       Cvec leading Eigenvector (frequency x channels)
%       Cent A different measure of total coherence: GM/AM of SV^2s
%       f (frequencies)
d=ndims(data);
if d<2, error('Need multidimensional array'); end
if d==2, [N,C]=size(data); end;
if d==3, [N,C,Ntr]=size(data); end;
if nargin < 3; params=[]; end;
[tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
clear err trialave params
nwin=round(win*Fs); nfft=2^(nextpow2(nwin)+pad);
[f,findx]=getfgrid(Fs,nfft,fpass);
tapers=dpsschk(tapers,nwin,Fs); % check tapers
Sc=zeros(length(findx),C,C);

Nwins=floor(N/nwin);

if d==3, % If there are multiple trials
for iwin=1:Nwins,
    for i=1:Ntr,
        data1=squeeze(data(1+(iwin-1)*nwin:iwin*nwin,:,i));
        J1=mtfftc(detrend(data1),tapers,nfft,Fs);
        J1=J1(findx,:,:);
        for k=1:C,
            for l=1:C,
                spec=squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2));
            Sc(:,k,l)=Sc(:,k,l)+spec;
            end
        end
    end
end
Sc=Sc/(Nwins*Ntr);
end

if d==2, % only one trial
for iwin=1:Nwins,
        data1=squeeze(data(1+(iwin-1)*nwin:iwin*nwin,:));
        J1=mtfftc(data1,tapers,nfft,Fs);
        J1=J1(findx,:,:);
        for k=1:C,
            for l=1:C,
            Sc(:,k,l)=Sc(:,k,l)+squeeze(mean(conj(J1(:,:,k)).*J1(:,:,l),2));
            end
        end
end
Sc=Sc/Nwins;
end

Cmat=Sc;
Sdiag=zeros(length(findx),C);
for k=1:C,
    Sdiag(:,k)=squeeze(Sc(:,k,k));
end

for k=1:C,
    for l=1:C,
        Cmat(:,k,l)=Sc(:,k,l)./sqrt(abs(Sdiag(:,k).*Sdiag(:,l)));
    end
end

Ctot=zeros(length(findx),1); Cent=Ctot;
Cvec=zeros(length(findx),C);
for i=1:length(findx),
    [u s]=svd(squeeze(Sc(i,:,:)));s=diag(s);
    Ctot(i)=s(1).^2/sum(s.^2); Cent(i)=exp(mean(log(s.^2)))/mean(s.^2);
    Cvec(i,:)=transpose(u(:,1));

end


text:
%%*******************************************************************
%% schurmat_sblk: compute Schur complement matrix corresponding to
%%                SDP blocks.
%%
%% symm = 0, HKM
%%      = 1, NT
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*******************************************************************

   function schur = schurmat_sblk(blk,At,par,schur,p,X,Y);

   global iter smallblkdim  nnzschur  nzlistschur

   if isempty(smallblkdim); smallblkdim = 15; end
   if (nargin == 7); symm = 0; else; symm = 1; Y = X; end;
   m = length(schur);
   pblk = blk(p,:);
   if (iter == 1)
      nnzschur(size(blk,1),1) = m*m;
      nzlistschur = cell(size(blk,1),1);
   end
%%
   if (max(pblk{2}) > smallblkdim)
      %%
      %% compute schur for matrices that are very sparse.
      %%
      m1 = size(At{p,1},2);
      if issparse(schur); schur = full(schur); end;
      J = min(m1, max(find(par.nzlistA{p,1} < inf))-1);
      if (J > 0)
         if issparse(X{p}) & ~issparse(Y{p}); X{p} = full(X{p}); end
         if ~issparse(X{p}) & issparse(Y{p}); Y{p} = full(Y{p}); end
         if (iter <= 3)
            [nnzschur(p),nzlisttmp] = mexschur(pblk,At{p,1},par.nzlistA{p,1},...
            par.nzlistA{p,2},par.permA(p,:),Y{p},X{p},J,symm,schur);
            if (nnzschur(p) == mexnnz(nzlisttmp))
               nzlistschur{p} = nzlisttmp;
            else
               nzlistschur{p} = [];
            end
         else
            if isempty(nzlistschur{p})
               mexschur(pblk,At{p,1},par.nzlistA{p,1},...
               par.nzlistA{p,2},par.permA(p,:),Y{p},X{p},J,symm,schur);
            else
               mexschur(pblk,At{p,1},par.nzlistA{p,1},...
               par.nzlistA{p,2},par.permA(p,:),Y{p},X{p},J,symm,schur,nzlistschur{p});
            end
         end
      end
      %%
      %% compute schur for matrices that are not so sparse or dense.
      %%
      if (m1 < m)
         ss = [0, cumsum(pblk{3})];
         if (length(At(p,:)) > 2)
            dd = At{p,3};
         else
            dd = ones(sum(pblk{3}),1);
         end
         XVD = X{p}*At{p,2}*spdiags(dd,0,length(dd),length(dd));
         YVD = Y{p}*At{p,2}*spdiags(dd,0,length(dd),length(dd));
      end
      L = max(find(par.nzlistAsum{p,1} < inf)) -1;
      if (J < L)
         len = par.nzlistAsum{p,1}(J+1); list = par.nzlistAsum{p,2}(1:len,:);
      end
      if (m1 > 0)
         for k = J+1:m
            if (k<=m1)
               isspAk = par.isspA(p,k);
               Ak = mexsmat(blk,At,isspAk,p,k);
               if (k <= L)
                  idx1 = par.nzlistAsum{p,1}(k)+1; idx2 = par.nzlistAsum{p,1}(k+1);
                  list = [list; par.nzlistAsum{p,2}(idx1:idx2,:)];
                  list = sortrows(list,[2 1]);
                  tmp = Prod3(pblk,X{p},Ak,Y{p},symm,list);
               else
                  tmp = Prod3(pblk,X{p},Ak,Y{p},symm);
               end
            else
               idx = [ss(k-m1)+1 :ss(k-m1+1)];
               tmp = XVD(:,idx)* (Y{p}*At{p,2}(:,idx))';
            end
            if (~symm)
               tmp = 0.5*(mexsvec(pblk,tmp) + mexsvec(pblk,tmp,[],1));
            else
               tmp = mexsvec(pblk,tmp);
            end
            permk = par.permA(p,k);
            idx  = par.permA(p,1:min(k,m1));
            tmp2 = schur(idx,permk) + mexinprod(blk,At,tmp,min(k,m1),p);
            schur(idx,permk) = tmp2;
            schur(permk,idx) = tmp2';
         end
      end
      if (m1 < m)
         m2 = m - m1;
         YVtmp = At{p,2}'*YVD;
         XVtmp = XVD'*At{p,2};
         for k = 1:m2
            idx0 = [ss(k)+1 : ss(k+1)];
            tmp = XVtmp(:,idx0) .* YVtmp(:,idx0);
            tmp = tmp*ones(length(idx0),1);
            tmp3 = schur(m1+[1:m2],m1+k) + mexqops(pblk{3},tmp,ones(length(tmp),1),1);
            schur(m1+[1:m2],m1+k) = tmp3;
         end
      end
   else
      if issparse(X{p}) & ~issparse(Y{p}); Y{p} = sparse(Y{p}); end
      if ~issparse(X{p}) & issparse(Y{p}); X{p} = sparse(X{p}); end
      tmp  = mexskron(pblk,X{p},Y{p});
      Perm = spconvert([(1:m)' par.permA(p,:)' ones(m,1)]);
      schurtmp = At{p,1}'*tmp*At{p,1};
      schurtmp = 0.5*(schurtmp + schurtmp');
      schur = schur + Perm'*schurtmp*Perm;
   end
%%*******************************************************************

author: yu-jiang

repo: radpbook

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


text:
function pass = test_uminus( pref )
% This tests the basic arithmetic operations on chebfun2 objects.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e5 * pref.cheb2Prefs.chebfun2eps;
j = 1;

D = [-1 1 -1 1; -2 2 -2 2; -1 pi 0 2*pi];

for r = 1 : size(D,1)
    f = chebfun2(@(x,y) cos(x.*y), D(r,:));

    uminusF = chebfun2(@(x,y) -cos(x.*y), D(r,:));

    tolr = norm(D(r,:),inf)*tol;

    pass(j) = ( norm( (-f) - uminusF ) < tolr ); j = j + 1;

end

end


text:
function [delay,maxBitRate,energy] = dtiNeuronBitRate(radius, length, meanFiringRate, myelinated)
%
% [delay,maxBitRate,energy] = dtiNeuronBitRate(radius, length, meanFiringRate, [myelinated=true])
%
% Returns neural conduction delay in msec, the maximum bit-rate, and the energy consumption
% given the outer radius (in micrometers), the length (in millimeters), and the mean firing
% rate (in spikes/sec).
%
% Energy requirements are specified in amol glucose / action potential.
%
% See:
%   Brenner et. al. (2000). Synergy in a neural code. Neural Comput.
%   Wang et. al. (2008). Functional Trade-Offs in White Matter Axonal Scaling. J. Neurosci.
%   (The main reference that works out the bit-rate.)
%
% Also see:
%   Timing jitter vs. radius: Swadlow (2000). Time and the brain (R. Miller, ed.)
%   (We assume that timing jitter is +/-10% or the total conduction delay.)
%
% 2009.02.06 RFD wrote it.

if(~exist('myelinated','var') || isempty(myelinated))
    myelinated = true;
end

[speed,energyPerMm] = dtiNeuralConductionSpeed(radius, myelinated);

delay = 1 ./ (speed ./ length);
% delay is in ms- we need seconds, thus the 1000
energy = energyPerMm.*length;
maxBitRate = log2(1./(meanFiringRate.*0.2.*(delay./1000)));

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [outStr outNum] = gunits(before,after)
%Function to convert units using Google
%SCd 09/01/2010
%
%gunits searches using the Google convention e.g: '3ft to in'
%
%Input Arguments:
%   -before: string with number and unit such as:
%         '3ft', '17c','17lbf/in^2', '3 ft*lbf'
%   -after: unit you wish to convert such as: (respectively)
%       'mm', 'f', 'N/ft^2','Newton*m'
%
%   NOTE1: Google can be picky with the input.  It's recommended
%       to spell everything out and fill in multiplication/division/powers
%       if you're having issues.  Example:
%
%       >>gunits('3psi','psf')
%           ans = '3 psi' %if outNum was called; an error would've been produced
%
%       >>gunits('3lbf/in^2','lbf/ft^2')
%           ans = 3 (lbf / (in^2)) = 432 lbf / (ft^2)
%
%   NOTE2: The units should be distinguishable with spaces removed e.g:
%      'lbf ft' becomes 'lbfft' which is not a unit.  To fix use 'lbf*ft'
%
%Output Arguments:
%   -outStr: string with the conversion
%   -outNum: numeric value of answer (only if asked for)
%
%NOTE/DISCLAIMER:
%   The algorithm used for retrieving the string is based solely on Google's
%   web format.  As (unfortunately!), I do not own or control Google, their
%   format is subject to change, which could nullify the results of this program.
%   Please let me know if this happens.
%

%Error checking:
assert(nargin==2,'gunits() expects two and only two input arguments!');
assert(ischar(before),'gunits expected its first input argument to be a string!');
assert(ischar(after),'gunits expected its second input argument to be a string!');

%Create search string, remove spaces, read it
search_string = ['http://www.google.com/search?hl=en&source=hp&q=' before '+to+' after];
S = urlread(search_string(~isspace(search_string)));

%The first HTML bold ('') in the returned string will be the initiallization of the
%answer.  Find it and corresponding close bold ('').
strinit = strfind(lower(S),'');
strend = strfind(lower(S),'');
if isempty(strinit) || isempty(strend)
    error('Input units were not accepted: please try again!');
end
strinit = strinit(1);
strend = strend(strend>strinit);
strend = strend(1);

%Output the answer (remove the , < found earlier)
outStr = S(strinit+3:strend-1);

%If the numeric value is wanted:
if nargout == 2
    %Split into parts; find equal sign; keep value after equal sign, remove
    %spaces that Google uses to show 1000's
    parts = regexp(outStr,' ','split');
    eqsign = cellfun(@(x)strcmp(x,'='),parts);
    assert(any(eqsign),'The conversion did not work: please try again');
    outNum = parts{find(eqsign)+1};
    outNum = str2double(outNum(~isspace(outNum)));
end


text:
function [L,P,Lc,Pc]=corr_meas(rtk,obs,nav,dantr,dants,phw)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SNR test not surpport
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global glc
lam=nav.lam(obs.sat,:);
L=zeros(glc.NFREQ,1);P=zeros(glc.NFREQ,1);Lc=0;Pc=0;

for i=1:glc.NFREQ
    L(i)=0; P(i)=0;
    if lam(i)==0||obs.L(i)==0||obs.P(i)==0,continue;end

    %antenna phase center and phase windup correction
    L(i)=obs.L(i)*lam(i)-dants(i)-dantr(i)-phw*lam(i);
    P(i)=obs.P(i)       -dants(i)-dantr(i);

end

% DCB correction
[cbias,~]=getdcb(nav,obs,rtk.opt);
for i=1:glc.NFREQ
    if P(i)~=0,P(i)=P(i)-cbias(i);end
end
C1= lam(2)^2/(lam(2)^2-lam(1)^2);
C2=-lam(1)^2/(lam(2)^2-lam(1)^2);

%IFLC measurements
if L(1)~=0&&L(2)~=0,Lc=C1*L(1)+C2*L(2);end
if P(1)~=0&&P(2)~=0,Pc=C1*P(1)+C2*P(2);end

return

repo: GINav


text:
function pass = test_complex(pref)
% Test complex

if ( nargin == 0)
    pref = chebfunpref;
end
tol = 1000*pref.cheb3Prefs.chebfun3eps;

f = chebfun3(@(x,y,z) sin(x.*y.*z));
g = chebfun3(@(x,y,z) cos(x.*y.*z));
h = chebfun3(@(x,y,z) sin(x.*y.*z) + 1i*cos(x.*y.*z));
pass = norm(h - complex(f, g)) < tol;

end


text:
function [co, amp, ph] = computeCorAnalTSeries(vw, scan, tSeries)

% Variable check
%
if notDefined('vw'),        vw = getCurView;                    end
if notDefined('scan'),      scan = viewGet(vw, 'curscan');      end

% Compute Fourier transform
%
ft = fft(tSeries);
ft = ft(1:1+fix(size(ft, 1)/2), :);

% This quantity is proportional to the amplitude
%
scaledAmp = abs(ft);

% This is in fact, the correct amplitude
%
nCycles = viewGet(vw, 'nCycles', scan);
amp = 2*(scaledAmp(nCycles+1,:))/size(tSeries,1);

% We use the scaled amp here which is OK for computing the
% correlation. Note that the noiseBand defines the portion
% of the spectrum to use for the noise metric. As such, this
% calculation now corresponds to the correlation of the fundamental
% stimulus sinusoid with a FILTERED version of the data, where
% noiseBand determines the passband of a square-edged bandpass
% filter.
%
noiseBand    = GetNoiseBand(vw, scan);
noiseIndices = CreateNoiseIndices(scaledAmp, nCycles, noiseBand);
sqrtsummagsq = sqrt(sum(scaledAmp(noiseIndices, :).^2));

% (ras 06/07: sometimes sqrtsummagsq can be zero for some voxels;
% don't throw a warning for this line only.)
warning off MATLAB:divideByZero
co = scaledAmp(nCycles+1,:)./sqrtsummagsq;
warning on MATLAB:divideByZero
clear scaledAmp

% Calculate phase:
% 1) add pi/2 so that it is in sine phase.
% 2) minus sign because sin(x-phi) is shifted to the right by phi.
% 3) Add 2pi to any negative values so phases increase from 0 to 2pi.
%
ph = -(pi/2) - angle(ft(nCycles+1,:));
ph(ph<0) = ph(ph<0)+pi*2;


text:
function errors = test_filter()
%TEST_FILTERS This function test all the filters

errors = 0;

errors = errors + test_mexican_hat1();
errors = errors + test_mexican_hat2();
errors = errors + test_meyer();
errors = errors + test_abspline();
errors = errors + test_simple_tf();
errors = errors + test_itersine();
errors = errors + test_uniform_half_cosine();
errors = errors + test_design_warped_translates();
errors = errors + test_held();
errors = errors + test_simoncelli();
errors = errors + test_papadakis();
errors = errors + test_regular();

errors = errors + test_filter_lanczos();
errors = errors + test_filter_lanczos_multiple();

errors = errors + test_filter_lanczos_syn();
errors = errors + test_filter_lanczos_syn_multiple();

errors = errors + test_cheb_coeff();
errors = errors + test_cheb_op();
errors = errors + test_filter_analysis();
errors = errors + test_filter_synthesis();
errors = errors + test_filter_inverse();
errors = errors + test_localization();

errors = errors + test_multiple_dimensions();

errors = errors + test_reshape();
errors = errors + test_matrix_op();

try  %#ok<TRYNC>
    close(100)
end

end

function errors = test_regular()

errors = 0;
try

   figure(100);
   G = gsp_sensor(100);
   G = gsp_estimate_lmax(G);
   g = gsp_design_regular(G);
   gsp_plot_filter(G,g);
   close(100);

   fprintf('FILTER: regular kernel 1 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error kernel regular 1 test')
end

try

   figure(100);
   G = gsp_sensor(100);
   G = gsp_estimate_lmax(G);
   param.d = 5;
   g = gsp_design_regular(G,param);
   gsp_plot_filter(G,g);
   close(100);

   fprintf('FILTER: regular kernel 2 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error kernel regular 2 test')
end

try

   figure(100);
   G = gsp_sensor(100);
   G = gsp_estimate_lmax(G);
   param.d = 0;
   g = gsp_design_regular(G,param);
   gsp_plot_filter(G,g);
   close(100);

   fprintf('FILTER: regular kernel 3 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error kernel regular 3 test')
end

end

function errors = test_papadakis()

errors = 0;
try

   figure(100);
   G = gsp_sensor(100);
   G = gsp_estimate_lmax(G);
   g = gsp_design_papadakis(G);
   gsp_plot_filter(G,g);
   close(100);

   fprintf('FILTER: papadakis kernel 1 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error kernel papadakis 1 test')
end

try

   figure(100);
   G = gsp_sensor(100);
   G = gsp_estimate_lmax(G);
   param.a = 0.25;
   g = gsp_design_papadakis(G,param);
   gsp_plot_filter(G,g);
   close(100);

   fprintf('FILTER: papadakis kernel 2 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error kernel papadakis 2 test')
end

end

function errors = test_simoncelli()

errors = 0;
try

   figure(100);
   G = gsp_sensor(100);
   G = gsp_estimate_lmax(G);
   g = gsp_design_simoncelli(G);
   gsp_plot_filter(G,g);
   close(100);

   fprintf('FILTER: simoncelli kernel 1 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error kernel simoncelli 1 test')
end

try

   figure(100);
   G = gsp_sensor(100);
   G = gsp_estimate_lmax(G);
   param.a = 0.25;
   g = gsp_design_simoncelli(G,param);
   gsp_plot_filter(G,g);
   close(100);

   fprintf('FILTER: simoncelli kernel 2 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error kernel simoncelli 2 test')
end

end

function errors = test_held()

errors = 0;
try

   figure(100);
   G = gsp_sensor(100);
   G = gsp_estimate_lmax(G);
   g = gsp_design_held(G);
   gsp_plot_filter(G,g);
   close(100);

   fprintf('FILTER: held kernel 1 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error kernel held 1 test')
end

try

   figure(100);
   G = gsp_sensor(100);
   G = gsp_estimate_lmax(G);
   param.a = 0.25;
   g = gsp_design_held(G,param);
   gsp_plot_filter(G,g);
   close(100);

   fprintf('FILTER: held kernel 2 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error kernel held 2 test')
end

end

function errors = test_design_warped_translates()

errors = 0;
try

    figure(100);
    Nf = 10;
    G = gsp_sensor(100);
    G = gsp_estimate_lmax(G);
    G = gsp_spectrum_cdf_approx(G);
    g = gsp_design_warped_translates(G, Nf);
    gsp_plot_filter(G,g);
    close(100);

   fprintf('FILTER: test_design_warped_translates 1 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error in test_design_warped_translates 1 test')
end

try

    figure(100);
    Nf = 10;
    G = gsp_sensor(100);
    G = gsp_estimate_lmax(G);
    param.log =1 ;
    param.warping_type = 'spectrum_approximation';
    G = gsp_spectrum_cdf_approx(G);
    g = gsp_design_warped_translates(G, Nf,param);
    gsp_plot_filter(G,g);
     close(100);

   fprintf('FILTER: test_design_warped_translates 2 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error in test_design_warped_translates 2 test')
end

try

    figure(100);
    Nf = 10;
    G = gsp_sensor(100);
    G = gsp_estimate_lmax(G);
    param.log =1 ;
    param.warping_type = 'custom';
    G = gsp_spectrum_cdf_approx(G);
    g = gsp_design_warped_translates(G, Nf,param);
    gsp_plot_filter(G,g);
    close(100);

   fprintf('FILTER: test_design_warped_translates 3 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error in test_design_warped_translates 3 test')
end

try

    figure(100);
    Nf = 10;
    G = gsp_sensor(100);
    G = gsp_estimate_lmax(G);
    param.log =0 ;
    param.warping_type = 'custom';
    G = gsp_spectrum_cdf_approx(G);
    g = gsp_design_warped_translates(G, Nf,param);
    gsp_plot_filter(G,g);
    close(100);

   fprintf('FILTER: test_design_warped_translates 4 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error in test_design_warped_translates 4 test')
end

try

    figure(100);
    Nf = 10;
    G = gsp_sensor(100);
    G = gsp_compute_fourier_basis(G);
    tol=1e-8;
    [unique_E,unique_E_inds]=unique(round(G.e*1/tol)*tol);
    param.approx_spectrum.x=unique_E;
    param.approx_spectrum.y=(unique_E_inds-1)/(G.N-1);
    param.interpolation_type='monocubic';
    param.log =0 ;
    param.warping_type = 'spectrum_interpolation';
    g = gsp_design_warped_translates(G, Nf,param);
    gsp_plot_filter(G,g);
    close(100);

   fprintf('FILTER: test_design_warped_translates 5 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error in test_design_warped_translates 5 test')
end

try

    figure(100);
    Nf = 10;
    G = gsp_sensor(100);
    G = gsp_compute_fourier_basis(G);
    tol=1e-8;
    approx_spectrum_inds=1:11:100;
    N_bar=length(approx_spectrum_inds);
    [unique_subset_E, unique_subset_E_inds]=unique(round(G.e(approx_spectrum_inds)*1/tol)*tol);
    param.approx_spectrum.x=unique_subset_E;
    param.approx_spectrum.y=(unique_subset_E_inds-1)/(N_bar-1);
    param.interpolation_type='pwl';
    param.log =1 ;
    param.warping_type = 'spectrum_interpolation';
    G = gsp_spectrum_cdf_approx(G);
    g = gsp_design_warped_translates(G, Nf,param);
    gsp_plot_filter(G,g);
    close(100);

   fprintf('FILTER: test_design_warped_translates 6 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error in test_design_warped_translates 6 test')
end

end

function errors = test_mexican_hat1()

errors = 0;
try

   G = gsp_sensor(100);
   gsp_design_mexican_hat(G, 5);

   fprintf('FILTER: mexican hat 1 ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error in the mexican hat 1 test')
end

end

function errors = test_mexican_hat2()

errors = 0;
try

    figure(100);
    Nf = 4;
    G = gsp_sensor(100);
    G = gsp_estimate_lmax(G);
    g = gsp_design_mexican_hat(G, Nf);
    gsp_plot_filter(G,g);
    close(100);

   fprintf('FILTER: mexican hat 2 ok\n');
catch
    errors = 1;
    warning('FILTER: Error in the mexican hat 2 test')
end

end

function errors = test_uniform_half_cosine()

errors = 0;
try

    figure(100);
    Nf = 4;
    G = gsp_sensor(100);
    G = gsp_estimate_lmax(G);
    g = gsp_design_half_cosine(G, Nf);
    gsp_plot_filter(G,g);
    close(100);

   fprintf('FILTER: uniform half cosine ok\n');
catch
    errors = 1;
    warning('FILTER: Error in uniform half cosine test')
end

end

function errors = test_abspline()

errors = 0;
try

    figure(100);
    Nf = 4;
    G = gsp_sensor(100);
    G = gsp_estimate_lmax(G);
    g = gsp_design_abspline(G, Nf);
    gsp_plot_filter(G,g);
    close(100);

   fprintf('FILTER: abspline ok\n');
catch
    errors = 1;
    warning('FILTER: Error abspline test')
end

end

function errors = test_meyer()

errors = 0;
try

    figure(100);
    Nf = 4;
    G = gsp_sensor(100);
    G = gsp_estimate_lmax(G);
    g = gsp_design_meyer(G, Nf);
    gsp_plot_filter(G,g);
    close(100);

   fprintf('FILTER: meyer ok\n');
catch
    errors = 1;
    warning('FILTER: Error in meyer test')
end

end

function errors = test_simple_tf()

errors = 0;
try

    figure(100);
    Nf = 4;
    G = gsp_sensor(100);
    G = gsp_estimate_lmax(G);
    g = gsp_design_simple_tf(G, Nf);
    gsp_plot_filter(G,g);
    close(100);

   fprintf('FILTER: simple_tf ok\n');
catch
    errors = 1;
    warning('FILTER: Error in simple_tf test')
end

end

function errors = test_itersine()

errors = 0;
try

        Nf = 20;
        G = gsp_sensor(100);
        G = gsp_estimate_lmax(G);
        g = gsp_design_itersine(G, Nf);
        gsp_plot_filter(G,g);
        [A,B] = gsp_filterbank_bounds(G,g);

   fprintf('FILTER: itersine ok\n');
catch
    errors = errors + 1;
    warning('FILTER: Error in itersine test')
end

if norm(A-B)<eps(1000)
    fprintf('FILTER: itersine 2 ok\n');
else
    errors = errors + 1;
    warning('FILTER: Error in itersine 2 test')
end

end

function errors = test_cheb_coeff()

 errors = 0;
 try

   Nf = 4;
    G = gsp_sensor(100);
    G = gsp_estimate_lmax(G);
    g = gsp_design_meyer(G, Nf);
    c = gsp_cheby_coeff(G, g);

  fprintf('FILTER: cheb coeff ok\n');

 catch
	 errors = 1;
	 warning('FILTER: Error in cheb coeff test')
end

end

function errors = test_cheb_op()

errors = 0;
try

        Nf = 4;
        G = gsp_sensor(100);
        G = gsp_estimate_lmax(G);
        g = gsp_design_meyer(G, Nf);
        c = gsp_cheby_coeff(G, g);
        f = rand(G.N,1);
        r = gsp_cheby_op(G, c, f);

   fprintf('FILTER: cheb op ok\n');
catch
    errors = 1;
    warning('FILTER: Error in cheb op test')
end

end

function errors = test_filter_analysis()

    Nf = 4;
    G = gsp_sensor(10);
    G = gsp_estimate_lmax(G);
    G = gsp_compute_fourier_basis(G);
    g = gsp_design_meyer(G, Nf);
    f = rand(G.N,1);

    param.method = 'exact';
    param.order = 100;
    f1 = gsp_filter_analysis(G,g,f,param);
    param.method = 'cheby';
    f2 = gsp_filter_analysis(G,g,f,param);

    if norm(f1(:)-f2(:))/norm(f1) < 1e-4
        errors = 0;
       fprintf('FILTER: analysis ok\n');
    else
        errors = 1;
        warning('FILTER: Error in analysis test')
    end

end

function errors = test_filter_synthesis()

    Nf = 4;
    G = gsp_sensor(10);
    G = gsp_estimate_lmax(G);
    G = gsp_compute_fourier_basis(G);
    g = gsp_design_simple_tf(G, Nf);
    f = rand(G.N,1);
    f = f/norm(f);

    param.method = 'cheby';
    param.order = 1000;
    f1 = gsp_filter_analysis(G,g,f,param);
    f2 = gsp_filter_synthesis(G,g,f1,param);
    param.method = 'exact';
    f3 = gsp_filter_synthesis(G,g,f1,param);
    f4 = f2/norm(f2);

    if norm(f-f4) < 1e-5
        errors = 0;
       fprintf('FILTER: synthesis ok\n');
    else
        errors = 1;
        norm(f-f4)
        warning('FILTER: Error in synthesis test')
    end

    if norm(f3-f2)/norm(f2) < 1e-5
       fprintf('FILTER: synthesis 2 ok\n');
    else
        errors = errors+1;
        warning('FILTER: Error in synthesis 2 test')
    end
end

function errors = test_filter_inverse()

    Nf = 4;
     G = gsp_sensor(10);
%     G = gsp_estimate_lmax(G);
    G = gsp_compute_fourier_basis(G);

    g = gsp_design_mexican_hat(G, Nf);
    f = rand(G.N,1);

    param.method = 'exact';
    f1 = gsp_filter_analysis(G,g,f,param);
    f2 = gsp_filter_inverse(G,g,f1,param);

    if norm(f-f2) < 1e-10
        errors = 0;
       fprintf('FILTER: inverse ok\n');
    else
        errors = 1;
        warning('FILTER: Error in inverse test')
    end

end

function errors = test_localization()

   % Number of nodes
    N = 100;

    % Kernel
    tau = 0.1;

    g = @(x) exp(-tau*x);

    G1 = gsp_sensor(N);

    G1 = gsp_estimate_lmax(G1);
    s1 = sqrt(G1.N)*gsp_filter_analysis(G1,g,eye(N));

    % Do the test
    G1 = gsp_compute_fourier_basis(G1);
    shat = g(G1.e);
    s = gsp_igft(G1,shat);

    s2 = zeros(G1.N);
    for ii = 1:G1.N
        s2(:,ii) = gsp_translate_old(G1, s, ii);
    end

    if norm(s1-s2) < 1e-10
        errors = 0;
       fprintf('FILTER: localisation ok\n');
    else
        errors = 1;
        warning('FILTER: Error in localisation test')
    end

end

function errors = test_multiple_dimensions()

   % Number of nodes
    N = 100;

    G1 = gsp_sensor(N);

    Nf = 4;
    k = 6;

    G1 = gsp_estimate_lmax(G1);
    g = gsp_design_meyer(G1, Nf);

    G1 = gsp_estimate_lmax(G1);

    f = rand(N,k);

    s1 = gsp_filter_analysis(G1,g,f);

    s2 = zeros(G1.N*Nf,k);
    for ii = 1:k
        s2(:,ii) =  gsp_filter_analysis(G1,g,f(:,ii));
    end

    if norm(s1(:)-s2(:)) < 1e-10
        errors = 0;
       fprintf('FILTER: multiple_dimensions ok\n');
    else
        errors = 1;
        warning('FILTER: Error in multiple_dimensions test')
    end

    f1 = gsp_filter_synthesis(G1,g,s1);

    f2 = zeros(G1.N,k);
    for ii = 1:k
        f2(:,ii) =  gsp_filter_synthesis(G1,g,s1(:,ii));
    end

     if norm(f1(:)-f2(:)) < 1e-10
       fprintf('FILTER: multiple_dimensions 2 ok\n');
    else
        errors = errors + 1;
        warning('FILTER: Error in multiple_dimensions 2 test')
     end

        f1 = gsp_filter_inverse(G1,g,s1);

    f2 = zeros(G1.N,k);
    for ii = 1:k
        f2(:,ii) =  gsp_filter_inverse(G1,g,s1(:,ii));
    end

    if norm(f1(:)-f2(:)) < 1e-10
        fprintf('FILTER: multiple_dimensions 3 ok\n');
    else
        errors = errors + 1;
        warning('FILTER: Error in multiple_dimensions 3 test')
    end

%     if norm(f1(:)-f(:))/norm(f(:)) < 1e-6
%         fprintf('FILTER: multiple_dimensions 4 ok\n');
%     else
%         errors = errors + 1;
%         norm(f1(:)-f(:))/norm(f(:))
%         warning('FILTER: Error in multiple_dimensions 4 test')
%     end

end

function errors = test_reshape()

   % Number of nodes
    N = 10;
    M = 4;

    f = rand(M*N,N);

    ft = gsp_vec2mat(f,M);
    [f1, M2] = gsp_mat2vec(ft);

    errors = 0;

    if norm(f1(:)-f(:))/norm(f(:)) < 1e-12
        fprintf('FILTER: reshape ok\n');
    else
        errors = errors + 1;
        warning('FILTER: erros reshape test')
    end

    if norm(M2-M) < 1e-12
        fprintf('FILTER: reshape 2 ok\n');
    else
        errors = errors + 1;
        warning('FILTER: erros reshape test 2')
    end
end

function errors = test_matrix_op()

    errors = 0;
    % Number of nodes
    N = 100;

    G1 = gsp_sensor(N);

    Nf = 4;

    G1 = gsp_estimate_lmax(G1);
    g = gsp_design_meyer(G1, Nf);

    G1 = gsp_estimate_lmax(G1);

    x = rand(N,1);

    f = gsp_filter_analysis(G1,g,x);
    F = gsp_filterbank_matrix(G1,g);
    f1 = F'*x;

    if norm(f1(:)-f(:))/norm(f(:)) < 1e-12
        fprintf('FILTER: operator matrix ok\n');
    else
        errors = errors + 1;
        warning('FILTER: error in operator matrix test')
    end

end

function errors = test_filter_lanczos()

    Nf = 4;
    G = gsp_sensor(10);
    G = gsp_estimate_lmax(G);
    G = gsp_compute_fourier_basis(G);
    g = gsp_design_meyer(G, Nf);
    f = rand(G.N,1);

    param.method = 'exact';
    param.order = 100;
    f1 = gsp_filter_analysis(G,g,f,param);
    param.method = 'lanczos';
    f2 = gsp_filter_analysis(G,g,f,param);

    if norm(f1(:)-f2(:))/norm(f1(:)) < 1e-10
        errors = 0;
       fprintf('FILTER: lanczos ok\n');
    else
        errors = 1;
        norm(f1(:)-f2(:))/norm(f1(:))
        warning('FILTER: Error in lanczos test')
    end

end

function errors = test_filter_lanczos_multiple()

    Nf = 4;
    G = gsp_sensor(10);
    G = gsp_estimate_lmax(G);
    G = gsp_compute_fourier_basis(G);
    g = gsp_design_meyer(G, Nf);
    f = rand(G.N,5);

    param.method = 'exact';
    param.order = 100;
    f1 = gsp_filter_analysis(G,g,f,param);
    param.method = 'lanczos';
    f2 = gsp_filter_analysis(G,g,f,param);

    if norm(f1(:)-f2(:))/norm(f1(:)) < 1e-10
        errors = 0;
       fprintf('FILTER: lanczos multiple ok\n');
    else
        errors = 1;
        norm(f1(:)-f2(:))/norm(f1(:))
        warning('FILTER: Error in lanczos multiple test')
    end

end

function errors = test_filter_lanczos_syn()

    Nf = 4;
    G = gsp_sensor(10);
    G = gsp_estimate_lmax(G);
    G = gsp_compute_fourier_basis(G);
    g = gsp_design_meyer(G, Nf);
    f = rand(G.N*Nf,1);

    param.method = 'exact';
    param.order = 100;
    f1 = gsp_filter_synthesis(G,g,f,param);
    param.method = 'lanczos';
    f2 = gsp_filter_synthesis(G,g,f,param);

    if norm(f1(:)-f2(:))/norm(f1(:)) < 1e-10
        errors = 0;
       fprintf('FILTER: lanczos synthesis ok\n');
    else
        errors = 1;
        norm(f1(:)-f2(:))/norm(f1(:))
        warning('FILTER: Error in lanczos synthesis test')
    end

end

function errors = test_filter_lanczos_syn_multiple()

    Nf = 4;
    G = gsp_sensor(10);
    G = gsp_estimate_lmax(G);
    G = gsp_compute_fourier_basis(G);
    g = gsp_design_meyer(G, Nf);
    f = rand(G.N*Nf,5);

    param.method = 'exact';
    param.order = 100;
    f1 = gsp_filter_synthesis(G,g,f,param);
    param.method = 'lanczos';
    f2 = gsp_filter_synthesis(G,g,f,param);

    if norm(f1(:)-f2(:))/norm(f1(:)) < 1e-10
        errors = 0;
       fprintf('FILTER: lanczos synthesis multiple test ok\n');
    else
        errors = 1;
        norm(f1(:)-f2(:))/norm(f1(:))
        warning('FILTER: Error in lanczos synthesis multiple test')
    end

end

repo: gspbox


text:
classdef (Abstract) nme_load < mp.nm_element

%   MATPOWER
%   Copyright (c) 2019-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end     %% properties

    methods
        function name = name(obj)
            name = 'load';
        end

        function np = np(obj)
            np = 1;     %% this is a 1 port element
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% pop_crossf() - Return estimates and plots of event-related spectral coherence
%
% Usage:
%       >> pop_crossf(EEG, typeproc, num1, num2, tlimits,cycles,
%                                 'key1',value1,'key2',value2, ... );
% Inputs:
%   INEEG    - Input EEG dataset
%   typeproc - Type of processing:
%                  1 = process two raw-data channels,
%                  0 = process two ICA components
%   num1     - First component or channel number
%   num2     - Second component or channel number
%   tlimits  - [mintime maxtime] Sub-epoch time limits in ms
%   cycles   -   >0 -> Number of cycles in each analysis wavelet
%                 0 -> Use FFTs (with constant window length)
%
% Optional inputs: As for crossf().  See >> help crossf
%
% Outputs: Same as crossf(). No outputs are returned when a
%          window pops-up to ask for additional arguments
%
% Author: Arnaud Delorme, CNL / Salk Institute, 11 March 2002
%
% See also: timef(), eeglab()

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

% 03-18-02 added title -ad & sm
% 04-04-02 added outputs -ad & sm

function varargout = pop_crossf(EEG, typeproc, num1, num2, tlimits, cycles, varargin );

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

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

	geometry = { [1 0.5 0.5] [1 0.5 0.5] [1 0.5 0.5] [1 0.5 0.5] [0.92 0.1 0.78] [1 0.5 0.5] [1 0.8 0.2] [1] [1 1]};
    uilist = { { 'Style', 'text', 'string', fastif(typeproc, 'First channel number', 'First component number'), 'fontweight', 'bold'  } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,3,[],'1') } {} ...
			   { 'Style', 'text', 'string', fastif(typeproc, 'Second channel number', 'Second component number'), 'fontweight', 'bold'  } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,4,[],'2') } {} ...
			   { 'Style', 'text', 'string', 'Epoch time range [min max] (msec)', 'fontweight', 'bold', ...
				 'tooltipstring', 'Sub epoch time limits' } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,5,[],[ int2str(EEG.xmin*1000) ' ' int2str(EEG.xmax*1000) ]) } {} ...
			   { 'Style', 'text', 'string', 'Wavelet cycles (0->FFT, see >> help timef)', 'fontweight', 'bold', ...
				 'tooltipstring', context('cycles',vars,txt) } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,6,[],'3 0.5') } {} ...
			   { 'Style', 'text', 'string',  '[set]->Linear coher / [unset]->Phase coher', 'fontweight', 'bold', ...
				 'tooltipstring', ['Compute linear inter-trial coherence (coher)' 10 ...
					'OR Amplitude-normalized inter-trial phase coherence (phasecoher)'] } ...
			   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'phasecoher','present',1) } { } ...
			   { 'Style', 'text', 'string', 'Bootstrap significance level (Ex: 0.01 -> 1%)', 'fontweight', 'bold', ...
				 'tooltipstring', context('alpha',vars,txt) } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,'alpha') } {} ...
			   { 'Style', 'text', 'string', 'Optional timef() arguments (see Help)', 'fontweight', 'bold', ...
				 'tooltipstring', 'See crossf() help via the Help button on the right...' } ...
			   { 'Style', 'edit', 'string', '''padratio'', 4' } ...
			   { 'Style', 'pushbutton', 'string', 'Help', 'callback',  'pophelp(''crossf'');' } ...
			   {} ...
			   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotersp','present',0), 'string', ...
				 'Plot coherence amplitude', 'tooltipstring', ...
				 'Plot coherence ampltitude image in the upper panel' } ...
			   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotphase','present',0), 'string', ...
				 'Plot coherence phase', 'tooltipstring', ...
				 'Plot coherence phase image in the lower panel' } ...
			 };

	result = inputgui( geometry, uilist, 'pophelp(''pop_crossf'');', ...
					   fastif(typeproc, 'Plot channel cross-coherence -- pop_crossf()', ...
							  'Plot component cross-coherence -- pop_crossf()'));
	if length( result ) == 0 return; end;

	num1     = eval( [ '[' result{1} ']' ] );
	num2     = eval( [ '[' result{2} ']' ] );
	tlimits	 = eval( [ '[' result{3} ']' ] );
	cycles	 = eval( [ '[' result{4} ']' ] );
    if result{5}
    	options = [',''type'', ''coher''' ];
    else
		options = [',''type'', ''phasecoher''' ];
    end;

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

    % add title
    % ---------
	if isempty( findstr(  'title', result{7}))
        if ~isempty(EEG.chanlocs) & typeproc
            chanlabel1 = EEG.chanlocs(num1).labels;
            chanlabel2 = EEG.chanlocs(num2).labels;
        else
            chanlabel1 = int2str(num1);
            chanlabel2 = int2str(num2);
        end;
		if result{5}
            options = [options ', ''title'',' fastif(typeproc, '''Channel ', '''Component ') chanlabel1 '-' chanlabel2 ...
					' Coherence'''];
        else
            options = [options ', ''title'',' fastif(typeproc, '''Channel ', '''Component ') chanlabel1 '-' chanlabel2 ...
					' Phase Coherence''' ];
		end;
	end;
	if ~isempty( result{6} )
		options      = [ options ', ''alpha'',' result{6} ];
	end;
	if ~isempty( result{7} )
		  options = [ options ',' result{7} ];
	end;
	if ~result{8}
		options = [ options ', ''plotersp'', ''off''' ];
	end;
	if ~result{9}
		options = [ options ', ''plotphase'', ''off''' ];
	end;
    figure; try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end;
else
	options = [ ',' vararg2str(varargin) ];
end;

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

% call function sample either on raw data or ICA data
% ---------------------------------------------------
if typeproc == 1
	tmpsig1 = EEG.data(num1,pointrange,:);
	tmpsig2 = EEG.data(num2,pointrange,:);
else
	if ~isempty( EEG.icasphere )
        tmpsig1 = eeg_getdatact(EEG, 'component', num1, 'samples', pointrange);
        tmpsig2 = eeg_getdatact(EEG, 'component', num2, 'samples', pointrange);
	else
		error('You must run ICA first');
	end;
end;
tmpsig1 = reshape( tmpsig1, 1, size(tmpsig1,2)*size(tmpsig1,3));
tmpsig2 = reshape( tmpsig2, 1, size(tmpsig2,2)*size(tmpsig2,3));

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

% plot the datas and generate output command
% --------------------------------------------
if length( options ) < 2
    options = '';
end;
varargout{1} = sprintf('figure; pop_crossf( %s, %d, %d, %d, [%s], [%s] %s);', ...
          inputname(1), typeproc, num1, num2, int2str(tlimits), num2str(cycles), options);

%options = [ options ', ''ydir'', ''norm''' ];
com = sprintf( '%s crossf( tmpsig1, tmpsig2, length(pointrange), [tlimits(1) tlimits(2)], EEG.srate, cycles %s);', outstr, options);
eval(com)

return;

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


text:
% dtft_test.m
% test mex versions of DTFT by comparing to m-file versions

for1 = @(om, x) jf_mex('dtft,forward', om', x, int32(1));
adj1 = @(om, Xm, N) jf_mex('dtft,adjoint', om', Xm, int32(N), int32(1));

%
% 1d tests
%
if 1, printm '1d dtft for'
	N1 = 32;
	M = 10;
	om = 2*pi*[0:(M-1)]'/M; % test with uniform frequency locations
	x = [1:N1]';
	Xm = dtft1(x, om);
	Xc = for1(om, x);
	printm('1d forward real mpd = %g%%', max_percent_diff(Xm,Xc))

	X2 = for1([om 0*om], x); % trick: fake as 2d!
	if any(X2 ~= Xc), error bug, end

	x = x - 2i*flipud(x);
	Xm = dtft1(x, om);
	Xc = for1(om, x);
	printm('1d forward cplx mpd = %g%%', max_percent_diff(Xm,Xc))

	X2 = for1(om, [x x]); % L=2 reps
	if any(X2 ~= [Xc Xc]), error bug, end
end

if 1, printm '1d dtft adj'
	xm = dtft_adj(Xm, om, N1, [], 0);
	xc = adj1(om, Xm, N1);
	printm('1d adjoint cplx mpd = %g%%', max_percent_diff(xm,xc))

	x2 = dtft_adj([Xm Xm], om, N1, [], 0);
%	if any(x2 ~= [xm xm]), error bugxm2, end % strangely nonzero!
	x2 = adj1(om, [Xm Xm], N1);
	if any(x2 ~= [xc xc]), error bug, end
end

%
% 2d tests
%
if 1, printm '2d dtft for'
	N1 = 4; N2 = 6;
	x = [[1:N1]'*ones(1,3), ones(N1,N2-3)]; % test signal

	o1 = 2*pi*[0:(N1-1)]'/N1; % test with uniform frequency locations
	o2 = 2*pi*[0:(N2-1)]'/N2; % test with uniform frequency locations
	[o1 o2] = ndgrid(o1, o2);
	omega = [o1(:) o2(:)];
	omega = [omega; [2 7]];

	Xm = dtft2(x, omega, [], 1);
	Xc = for1(omega, x);
	printm('2d forward real mpd = %g%%', max_percent_diff(Xm,Xc))

	x = x - 2i * flipud(x); % complex
	Xm = dtft(x, omega, 'how', 'loop');
	Xc = for1(omega, x);
	printm('2d forward cplx mpd = %g%%', max_percent_diff(Xm,Xc))

	X2 = for1(omega, stackup(x,x));
	if any(X2 ~= [Xc Xc]), error bug, end
end

if 1, printm '2d dtft adj'
	xm = dtft2_adj(Xm, omega, N1, N2, [], 0);
	xc = adj1(omega, Xm, [N1 N2]);
	printm('2d adjoint cplx mpd = %g%%', max_percent_diff(xm,xc))

	if 1 % stack 2d adjoint
		X = [Xm Xm]; % L=2
		x2 = dtft_adj(X, omega, [N1 N2], [], 0);
	%	if any(x2 ~= stackup(xm, xm)), error bug, end % strangely nonzero!
		x2 = adj1(omega, X, [N1 N2]);
		if any(x2 ~= stackup(xc, xc)), error bug, end
	end

	if 0 % real 2d spectrum (not implemented, and probably not needed)
		X = real(Xm);
		xm = dtft2_adj(X, omega, N1, N2, [], 0);
		xc = adj1(omega, X, [N1 N2]);
		printm('2d adjoint real mpd = %g%%', max_percent_diff(xm,xc))
	end
end

%
% 3d tests
%
if 1, printm '3d dtft for'
	N1 = 4; N2 = 6; N3 = 5;
	rng(3)
	x = rand(N1,N2,N3);

	o1 = 2*pi*[0:(N1-1)]'/N1; % test with uniform frequency locations
	o2 = 2*pi*[0:(N2-1)]'/N2;
	o3 = 2*pi*[0:(N3-1)]'/N3;
	[o1 o2 o3] = ndgrid(o1, o2, o3);
	omega = [o1(:) o2(:) o3(:)];
	omega = [omega; [2 7 1]];

	Xm = dtft(x, omega, 'how', 'loop');
	Xc = for1(omega, x);
	printm('3d forward real mpd = %g%%', max_percent_diff(Xm,Xc))

	x = x - 2i * flipdim(x,3); % complex
	Xm = dtft(x, omega, 'how', 'loop');
	Xc = for1(omega, x);
	printm('3d forward cplx mpd = %g%%', max_percent_diff(Xm,Xc))

	X3 = for1(omega, stackup(x,x));
	if any(X3 ~= [Xc Xc]), error bug, end
end

if 1, printm '3d dtft adj'
	xm = dtft_adj(Xm, omega, [N1 N2 N3], [], 0);
	xc = adj1(omega, Xm, [N1 N2 N3]);
	printm('3d adjoint cplx mpd = %g%%', max_percent_diff(xm,xc))

	if 1 % stack 3d adjoint
		X = [Xm Xm]; % L=2
		x2 = dtft_adj(X, omega, [N1 N2 N3], [], 0);
	%	if any(x2 ~= stackup(xm, xm)), error bug, end % strangely nonzero!
		x2 = adj1(omega, X, [N1 N2 N3]);
		if any(x2 ~= stackup(xc, xc)), error bug, end
	end

	if 0 % real 3d spectrum (not implemented, and probably not needed)
		X = real(Xm);
		xm = dtft_adj(X, omega, [N1 N2 N3], [], 0);
		xc = adj1(omega, X, [N1 N2 N3]);
		printm('3d adjoint real mpd = %g%%', max_percent_diff(xm,xc))
	end
end

%
% now test thread speedup
%
if 1 % thread setup
	clear
	N1 = 2^6;
	N2 = 2^5;
	M = 2^12;
	rng(0)
	omega = rand(M, 2);
	s0 = rand(M,1) + 1i * rand(M,1);
	x0 = rand(N1,N2) + 1i * rand(N1,N2);
	nthread = int32(8);

	arg1 = {'dtft,forward', omega', x0};
	arg2 = {'dtft,adjoint', omega', s0, int32([N1 N2]')};
end

if 1, printm 'threads forw'
	tic
	s1 = jf_mex(arg1{:}, int32(1));
	t1 = toc;
	tic
	s2 = jf_mex(arg1{:}, int32(2));
	t2 = toc;

	if any(s1 ~= s2), error bug, end
	printm('thread for time %g %g, speedup = %g', t1, t2, t1 / t2)
	if 1 % vs exact
		se = dtft(x0, omega, 'how', 'loop');
		printm('2d forward exct mpd = %g%%', max_percent_diff(se,s1))
	end
end

if 1, printm 'threads adj'
	tic
	x1 = jf_mex(arg2{:}, int32(1));
	t1 = toc;
	tic
	x2 = jf_mex(arg2{:}, int32(2));
	t2 = toc;
%	printm('2d adjoint mpd2 = %g%%', max_percent_diff(xm,x2))

	if nrms(x2,x1) > 1e-11, warning 'bug-adj-thread'
		max_percent_diff(x2,x1), end
	printm('thread adj time %g %g, speedup = %g', t1, t2, t1 / t2)
	if 1 % vs exact
		xm = dtft2_adj(s0, omega, N1, N2, [], 0);
		printm('2d adjoint exct mpd = %g%%', max_percent_diff(xm,x1))
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [RF] = applyEnsembleRF_table(table,Y,categories,cost)

% INITIALIZATION
nNeg = sum(~Y); % Number of negative instances
nPos = sum(Y); % Number of positive instances
indPos = find(Y); % Indexes of positive instances
indNeg = find(~Y); % Indexes of negative instances
nSub = round(nNeg/nPos); % Number of subsets

% RANDOM NUMBER GENERATOR SEED
if ~RandStream.getGlobalStream.Seed
    rng('shuffle')
end

if nargin < 4
    cost = 1;
end

% FIND THE NUMBER OF NEGATIVE INSTANCES IN EACH SUBSET
nNegS = nNeg / nSub; nSup = ceil(nNegS); nInf = floor(nNegS);
if nSup ~= nInf
    nSubInf = nSub - 1; nSubSup = 1; total = nSubInf*nInf + nSubSup*nSup;
    while total ~= nNeg
        nSubInf = nSubInf - 1; nSubSup = nSubSup + 1;
        total = nSubInf*nInf + nSubSup*nSup;
    end
    nNegS = [repmat(nInf,[1,nSubInf]),repmat(nSup,[1,nSubSup])];
else % The number of negative instances in all partitions will be the same
    nNegS = repmat(nSup,[1,nSub]);
end

% FINDING THE INDEXES OF NEGATIVE INSTANCES IN EACH SUBSET
indNegSub = cell(1,nSub);
for i = 1:nSub-1
    indNegSub{i} = zeros(nNegS(i),1);
    indTemp = ceil(numel(indNeg)*rand(nNegS(i),1));
    indTemp = unique(indTemp);
    total = numel(indTemp);
    while total ~= nNegS(i)
        indMore = ceil(numel(indNeg)*rand(nNegS(i)-total,1));
        indTemp = [indTemp;unique(indMore)];
        indTemp = unique(indTemp);
        total = numel(indTemp);
    end
    indNegSub{i}(:) = indNeg(indTemp);
    indNeg(indTemp) = [];
end
indNegSub{end} = indNeg;

% COMPUTING A RANDOM FOREST FOR EACH SUBSET AND APPENDING
ind = [indNegSub{1};indPos];
Xtemp = table(ind,:);
Ytemp = Y(ind);
[Xtemp,Ytemp] = shufflePartition(Xtemp,Ytemp);
RF = TreeBagger(1,Xtemp,Ytemp,'CategoricalPredictors',categories,'SampleWithReplacement','on','Cost',[0,1/cost;1,0]);
for i = 2:nSub
    ind = [indNegSub{i};indPos];
    Xtemp = table(ind,:);
    Ytemp = Y(ind);
    [Xtemp,Ytemp] = shufflePartition(Xtemp,Ytemp);
    RFtemp = TreeBagger(1,Xtemp,Ytemp,'CategoricalPredictors',categories,'SampleWithReplacement','on','Cost',[0,1/cost;1,0]);
    RF = append(RF,RFtemp);
end

end

author: mvallieres

repo: radiomics


text:
function params = parseVarArgs(params,varargin)
% Parse variable input arguments supplied in name/value format.
%
%    params = parseVarArgs(params,'property1',value1,'property2',value2) sets
%    the fields propertyX in p to valueX.
%
%    params = parseVarArgs(params,varargin{:},'strict') only sets the field
%    names already present in params. All others are ignored.
%
% AE 2007-06-01

if isempty(varargin)
    return
end

% check if correct number of inputs
if mod(length(varargin),2)
    if ~strcmp(varargin{end},'strict')
        err.message = 'Name and value input arguments must come in pairs.';
        err.identifier = 'parseVarArgs:wrongInputFormat';
        error(err)
    else
        % in 'strict' case, remove all fields that are not already in params
        fields = fieldnames(params);
        ndx = find(~ismember(varargin(1:2:end-1),fields));
        varargin([2*ndx-1 2*ndx end]) = [];
    end
end

% parse arguments
for i = 1:2:length(varargin)
    if ischar(varargin{i})
        params.(varargin{i}) = varargin{i+1};
    else
        err.message = 'Name and value input arguments must come in pairs.';
        err.identifier = 'parseVarArgs:wrongInputFormat';
        error(err)
    end
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% QPSK demonstration packet-based transceiver for Chilipepper
% Toyon Research Corp.
% http://www.toyon.com/chilipepper.php
% Created 10/17/2012
% embedded@toyon.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%#codegen
function [d_i_out, d_q_out] = qpsk_rx_srrc(d_i_in, d_q_in)

persistent buf_i buf_q

OS_RATE = 8;
f = SRRC;

if isempty(buf_i)
    buf_i = zeros(1,OS_RATE*2+1);
    buf_q = zeros(1,OS_RATE*2+1);
end

buf_i = [buf_i(2:end) d_i_in];
buf_q = [buf_q(2:end) d_q_in];

d_i_out = buf_i*f;
d_q_out = buf_q*f;


text:
function[cmap] = cmapThingsStuff()
% [cmap] = cmapThingsStuff()
%
% Returns the color map for stuff and thing labels in CocoStuff.
%
% Copyright by Holger Caesar, 2017

% Get stuff and thing color
stuffColors = cmapStuff();
thingColors = cmapThings();

% Remove duplicate backgrund
stuffColors(1, :) = [];
thingColors(1, :) = [];

% Combine unlabeled (black), things and stuff
cmap = [0, 0, 0; thingColors; stuffColors];

author: nightrome


text:
clc;
close all;
clearvars;

rng default;

A = mat_selector();

k  = 30;
spx_bdpro(A, k);

author: indigits

repo: sparse-plex


text:
%% Import/Export
%
%%

author: mtex-toolbox

repo: mtex

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


text:
function ref = randref(low, high, N)

% Create roughly N random integers between between
% low and high. Each value is unique, the values
% are monotonically increasing.
% These integer values can be used as reference indices
% for statistical purposes
% -1 for N means : return low:high
%
% length(randref(low, high, N)) will not be exactly equal to N !!!

narginchk(3,3);

if (high<=low)
	error('Cannot create random indices when upper limit is smaller than lower limit');
end

if N < 1 | N > high-low
    ref = low:high;
    return
end

m_target = (high-low) / N;

if m_target < 10
    m = (0.99+log10(m_target)) *  m_target;
    N  = ceil(1.5 * N);     % calculate more indices than necesary
else
    m = 1.99 * m_target;
    N = ceil(1.2*N);        % calculate more indices than necesary
end

ref = cumsum([low ; ceil(rand(N,1)*m)]);    % compute indices
ref = ref(find(ref <= high));               % remove out of range indices

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa