text:
%
% Test to make sure the acceptance ratios are correct for the MH moves.
% This compares the ratios from the relevant conditional distributions
% to the overall generative score
%
% If this does not throw an assert, then it should be working.
%
ntry = 1;
verbose = true;
debug = true;
ps = defaultps;
load(ps.libname,'lib');
for t=1:ntry
    fprintf(1,'TRY %d\n',t);
    motor_program = generate_all_rel_character(lib);
    M = motor_program();
    mcmc_all(M,lib,1,'verbose','debug','type','token','relation');
end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
function fb = h_goNogo(yt,t,in)
% generates the feedback (go/nogo task) given the agent's last choice
if yt==1 % 'go'
    fb = in.u(t);
else % 'nogo'
    fb = 0;
end

author: MBB-team

repo: VBA-toolbox


text:
classdef HPOLD < Algorithm
    %HPOLD Hierarchical Partial Order Label Decomposition performs partial
    %order classification with a hierarchical model [1]. Class 1 versus the
    %rest is first classified with a binary model and the remainder labels
    %(2,3,...,Q) with an ordinal model. The available methods are logistic
    %regression and suport vector machine. For additional details see [1].
    %
    %   SVR methods:
    %      fitpredict               - runs the corresponding algorithm,
    %                                   fitting the model and testing it in a dataset.
    %      fit                        - Fits a model from training data
    %      predict                    - Performs label prediction
    %
    %   References:
    %     [1] J. Sánchez-Monedero, M. Pérez-Ortiz, A. Sáez, P.A. Gutiérrez,
    %         and C. Hervás-Martínez. "Partial order label decomposition
    %         approaches for melanoma diagnosis". Applied Soft Computing.
    %         Volume 64, March 2018, Pages 341-355.
    %         https://doi.org/10.1016/j.asoc.2017.11.042
    %
    %   This file is part of ORCA: https://github.com/ayrna/orca
    %   Original authors: Pedro Antonio Gutiérrez, María Pérez Ortiz, Javier Sánchez Monedero
    %   Citation: If you use this code, please cite the associated paper http://www.uco.es/grupos/ayrna/orreview
    %   Copyright:
    %       This software is released under the The GNU General Public License v3.0 licence
    %       available at http://www.gnu.org/licenses/gpl-3.0.html
    properties
        description = 'Hierarchical Partial Order Label Decomposition';
        %C penalty coefficient and the kernel parameters (both for the binary
        %and ordinal methods).
        parameters = struct('C', 0.1, 'k', 0.1);
        binaryMethod = 'SVC1V1';
        ordinalMethod = 'SVORIM';
    end
    properties(Access = private)
        objBI;
        objOR;
    end

    methods

        % TODO: update doc
        function obj = HPOLD(varargin)
            %HPOLD constructs an object of the class HPOLD and sets its default
            %   characteristics
            %   OBJ = HPOLD() builds HPOLD
            obj.parseArgs(varargin);
        end

        function [projectedTrain, predictedTrain] = privfit( obj, train, param)
            %PRIVFIT trains the model for the HPOLD method with TRAIN data and
            %vector of parameters PARAM.

            projectedTrain = -1*ones(length(train.targets),1);% dummy value
            predictedTrain = -1*ones(length(train.targets),1);% dummy value

            % Create binary dataset
            trainTargetsBi = ones(size(train.targets));
            trainTargetsBi(train.targets~=1)=2;

            trainBi.patterns = train.patterns;
            trainBi.targets = trainTargetsBi;

            % Create ordinal dataset by removing class obj.binClass patterns
            % and relabelling the rest of the labels
            trainOr.patterns = train.patterns(train.targets~=1,:);
            trainOr.targets = train.targets(train.targets~=1,:) - 1;

            % Create and train binary classifier
            switch(lower(obj.binaryMethod))
                case 'svc1v1'
                    parambi.C = param.C;
                    parambi.k = param.k;
                    obj.objBI = SVC1V1();
                    obj.objBI.fit(trainBi, parambi);
                case 'csvc1v1'
                    error('TODO')
                    parambi.C = param.C;
                    parambi.k = param.k;
                    obj.objBI = CSVC();
                    obj.objBI.fit(trainBi, parambi);
                case 'liblinear'
                    parambi.C = param.C;
                    obj.objBI = LIBLINEAR();
                    obj.objBI.fit(trainBi, parambi);
                otherwise
                    error(['Unknown binary classifier method:', obj.binaryMethod])
            end

            % Create and train ordinal classifier
            switch(lower(obj.ordinalMethod))
                case 'svorim'
                    obj.objOR = SVORIM();
                    paramor.C = param.C;
                    paramor.k = param.k;
                    obj.objOR.fit(trainOr,paramor);
                    %obj.objOR.model = obj.objOR.fit(trainOr,paramor);
                case 'pom'
                    obj.objOR = POM();
                    obj.objOR.fit(trainOr);
                    %obj.objOR.model = obj.objOR.fit(trainOr);
                otherwise
                    error(['Unknown ordinal classifier method:', obj.ordinalMethod])
            end

            % Save model and parameters
            model.parameters = param;
            model.modelBI = obj.objBI.getModel();
            model.modelOR = obj.objOR.getModel();
            obj.model = model;
            [projectedTrain, predictedTrain] = obj.predict(train.patterns);
        end

        function [projected, predTargets]= privpredict(obj, testPatterns)
            %PREDICT predict labels of TEST patterns labels using MODEL.
            projected = -1*ones(size(testPatterns,1),1);% dummy value
            % Binary prediction: classes 1/2
            [projectedTest_bi,predTargets] = obj.objBI.predict(testPatterns);
            % Ordinal prediction for patterns of class ~= class 1
            [projectedTest_or, predictedTest_or] = obj.objOR.predict(testPatterns(predTargets~=1,:));
            % +1 to correct label numbering
            predictedTest_or = predictedTest_or + 1;
            predTargets(predTargets~=1,:) = predictedTest_or;
        end

    end

end

author: ayrna

repo: orca

save_path: github-repos/MATLAB/ayrna-orca


text:
function v= visutil_goodContourValues(mi, ma, ct)
% to be optimized

if ct<0,
  spacing= (ma-mi)/(-ct+1);
  gro= floor(log10(spacing));
  res= spacing*10^-gro;
  if res<1.5,
    resi= 1;
  elseif res<3.5,
    resi= 2;
  elseif res<7,
    resi= 5;
  else
    resi= 1;
    gro= gro+1;
  end
  ct= resi*10^gro;
end

v= ct*ceil(mi/ct):ct:ct*floor(ma/ct);

author: bbci


text:
%% gviridis
% Below is a demonstration of the features of the |gviridis| function

%%
clear; close all; clc;

%% Syntax
% |[cMap]=gviridis(varargin);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function value = r8_exponential_01_pdf ( rval )

%*****************************************************************************80
%
%% R8_EXPONENTIAL_01_PDF: PDF of a standard exponential distribution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2013
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Input, real RVAL, the point where the PDF is evaluated.
%
%    Output, real VALUE, the value of the PDF.
%
  if ( rval < 0.0 )
    value = 0.0;
  else
    value = exp ( - rval );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [h,s,v] = hsl2hsv(hh,ss,ll)

h = hh;
ll = ll*2;

ind = ll <= 1;
ss(ind) = ss(ind) .* ll(ind);
ss(~ind) = ss(~ind) .* (2 - ll(~ind));

v = (ll + ss) ./ 2;
s = (2 * ss) ./ (ll + ss);
s(isnan(s)) = 0;

end

% r = linspace(1,2,30);
% omega = linspace(0,2*pi);
% [r,omega] = meshgrid(r,omega);
% x = r .* cos(omega);
% y = r .* sin(omega);
% z = zeros(size(r));
% rgb = hsv2rgb(omega./2./pi,ones(size(omega)),ones(size(omega)));
% surf(x,y,z,rgb)
% axis equal

author: mtex-toolbox

repo: mtex

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


text:
function [An Bn] = mieLayeredTerms(mu, epsilon, radius, isPEC, frequency, nmax)

% Compute the sphere coefficients An and Bn for a sphere having a number
% of homogeneous dielectric layers, with or without a perfectly
% electrically conducting (PEC) core. Follows the treatment in
% Chapter 3 of
%
% Ruck, et. al. "Radar Cross Section Handbook", Plenum Press, 1970.
%
% The incident electric field is in the -z direction (theta = 0) and is
% theta-polarized. The time-harmonic convention exp(jwt) is assumed, and
% the Green's function is of the form exp(-jkr)/r.
%
% Inputs:
%   mu: Relative complex permeability for each dielectric region
%   epsilon: Relative complex permittivity for each dielectric region
%       (the unbounded (free space) region should be the first entry)
%   radius: Radius of each spherical dielectric interface, from outermost
%           to innermost
%   isPEC: A flag for when the innermost sphere is conducting
%   frequency: Operating frequency (Hz)
%   nmax: Maximum mode for computing Bessel functions
% Outputs:
%   An: Array of Mie solution constants (used in mieScatteredField)
%   Bn: Array of Mie solution constants (used in mieScatteredField)
%
%   Author: Walton C. Gibson, email: kalla@tripoint.org

% speed of light
c = 299792458.0;

% radian frequency
w = 2.0*pi*frequency;

% wavenumber
k0 = w/c;

% total number of dielectric interfaces
numInterfaces = length(mu);

% free space impedance
eta0 = 376.7303134617707;

% impedance of each dielectric region
eta = eta0 * sqrt(mu./epsilon);

% scale factor
m = sqrt(mu.*epsilon);

% mode numbers
mode = 1:nmax;

Z = zeros(numInterfaces, length(mode));
Y = zeros(numInterfaces, length(mode));

numIFC = numInterfaces - 1;

for L = numIFC:-1:1

    if L == numIFC
        if isPEC == 0
            x = k0 * m(L + 1) * radius(L);
            [J JP] = SphericalBessel_JJP(mode, x);
            % Ruck, et. al. (3.4-23)
            P = (J ./ JP);
            % Ruck, et. al. (3.4-26), at innermost interface
            Z(L,:) = eta(L + 1) * P;
            % Ruck, et. al. (3.4-26), at innermost interface
            Y(L,:) = P / eta(L + 1);
        end
    else
        x1 = k0 * m(L + 1) * radius(L);
        x2 = k0 * m(L + 1) * radius(L + 1);
        [J1 JP1] = SphericalBessel_JJP(mode, x1);
        [H1 HP1] = SphericalHankel_HHP(mode, 2, x1);
        [J2 JP2] = SphericalBessel_JJP(mode, x2);
        [H2 HP2] = SphericalHankel_HHP(mode, 2, x2);
        % Ruck, et. al. (3.4-24)
        U = JP2 .* HP1 ./ (JP1 .* HP2);
        % Ruck, et. al. (3.4-25)
        V = (J2 .* H1) ./ (J1 .* H2);
        % Ruck, et. al. (3.4-23)
        P1 = J1 ./ JP1;
        % Ruck, et. al. (3.4-23)
        P2 = J2 ./ JP2;
        % Ruck, et. al. (3.4-23)
        Q2 = H2 ./ HP2;
        if (L == numIFC - 1) && (isPEC == 1)
            % Ruck, et. al. (3.4-27) for PEC boundary condition
            Z(L,:) = eta(L + 1) * P1 .* (1.0 - V) ./ (1.0 - U .*  P2 ./ Q2);
            % Ruck, et. al. (3.4-28) for PEC boundary condition
            Y(L,:) = (P1 / eta(L + 1)) .* (1.0 - V .*  Q2 ./ P2) ./ (1.0 - U);
        else
            % Ruck, et. al. (3.4-27)
            Z(L,:) = eta(L + 1) * P1 .* (1.0 - V .*((1.0 - Z(L+1,:)./(eta(L+1)*P2))./(1.0 - Z(L+1,:)./(eta(L+1)*Q2)))) ./ ...
                (1.0 - U .* ((1.0 - eta(L+1)*P2./Z(L+1,:))./(1.0 - eta(L+1)*Q2./Z(L+1,:))));
            % Ruck, et. al. (3.4-28)
            Y(L,:) = (P1 / eta(L + 1)) .* (1.0 - V .*((1.0 - eta(L+1)* Y(L+1,:)./ P2)./(1.0 - eta(L+1)*Y(L+1,:)./Q2))) ./ ...
                (1.0 - U .* ((1.0 - P2./(eta(L+1)*Y(L+1,:)))./(1.0 - Q2./(eta(L + 1)*Y(L+1,:)))));
        end
    end

end

% Ruck, et. al. (3.4-29)
Zn = i*Z(1,:)/eta0;
% Ruck, et. al. (3.4-29)
Yn = i*eta0*Y(1,:);

x = k0*radius(1);
[J JP] = SphericalBessel_JJP(mode, x);
[H HP] = SphericalHankel_HHP(mode, 2, x);

% Ruck, et. al. (3.4-1)
An = -((i).^(mode)) .* (2*mode + 1) ./ (mode.*(mode + 1)) .* (J + i*Zn .* JP) ./ (H + i*Zn .* HP);

% Ruck, et. al. (3.4-2) - there is an error in Ruck which is fixed here
Bn = ((i).^(mode+1)) .* (2*mode + 1) ./ (mode.*(mode + 1)) .* (J + i*Yn .* JP) ./ (H + i*Yn .* HP);

    function [J JP] = SphericalBessel_JJP(mode, x)

        % compute spherical bessel functions and their derivatives

        s = sqrt(0.5*pi/x);

        [J] = besselj(mode + 1/2, x); J = J*s;
        [J2] = besselj(mode - 1/2, x); J2 = J2*s;

        JP = (x * J2 - mode .* J);
        J = x*J;

    end

    function [H HP] = SphericalHankel_HHP(mode, arg, x)

        % compute spherical hankel functions and their derivatives

        s = sqrt(0.5*pi/x);

        [H] = besselh(mode + 1/2, arg, x); H = H*s;
        [H2] = besselh(mode - 1/2, arg, x); H2 = H2*s;

        HP = (x * H2 - mode .* H);
        H = x*H;
    end

end

author: Sable

repo: mcbench-benchmarks


text:
% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
function qq = mrdivide(q1, q2)
%Quaternion.mrdivide Compute quaternion quotient.
%
% Q1/Q2   is a quaternion formed by Hamilton product of Q1 and inv(Q2)
% Q/S     is the element-wise division of quaternion elements by by the scalar S

    if isa(q2, 'Quaternion')
        % qq = q1 / q2
        %    = q1 * qinv(q2)

        qq = q1 * inv(q2);
    elseif isa(q2, 'double')
        qq = Quaternion( double(q1) / q2 );
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [ut, itn] = project_lin_sys( work, data, n, m, u, v, rho_x, i, use_indirect, cg_rate, extra_verbose, h, g, gTh)
    ut = u+v;
    ut(1:n) = rho_x*ut(1:n);
    ut(1:n+m) = ut(1:n+m) - ut(end)*h;
    ut(1:n+m) = ut(1:n+m) - h*((g'*ut(1:n+m))/(gTh+1));
    warm_start = u(1:n+m);
    ut(n+1:end-1) = -ut(n+1:end-1);
    [ut(1:n+m), itn] = solve_lin_sys(work, data, ut(1:n+m), n, m, warm_start, rho_x, i, use_indirect, cg_rate, extra_verbose);
    ut(end) = (ut(end) + h'*ut(1:n+m));
end

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [H,info] = freqwavelet(name,L, varargin)
%FREQWAVELET  Wavelet in the freq. domain
%   Usage: H=freqwavelet(name,L)
%          H=freqwavelet(name,L,scale)
%          [H,info]=freqwavelet(...)
%
%   Input parameters:
%         name  : Name of the wavelet
%         L     : System length
%   Output parameters:
%         H     : Frequency domain wavelet
%         info  : Struct with additional outputs
%
%   `freqwavelet(name,L)` returns peak-normalized "mother" frequency-domain
%   wavelet `name` for system length *L*. The basic scale is selected such
%   that the peak is positioned at the frequency 0.1 relative to the Nyquist
%   rate (fs/2).
%
%   The supported wavelets that can be used in place of `name` (the actual
%   equations can be found at the end of this help):
%
%     'cauchy'    Cauchy wavelet with alpha=100. Custom order=(alpha-1)/2
%                 (with alpha>1) can be set by `{'cauchy',alpha}`. The
%                 higher the order, the narrower and more symmetric the
%                 wavelet. A numericaly stable formula is used in order
%                 to allow support even for very high `alpha` e.g.
%                 10^6.
%
%     'morse'     Morse wavelet with alpha=100 and gamma=3. Both parameters
%                 can be set directly by `{'morse',alpha,gamma}`. `alpha`
%                 has the same role as for 'cauchy', `gamma` is the
%                 'skewness' parameter. 'morse' becomes 'cauchy' for
%                 gamma=1.
%
%     'morlet'    Morlet wavelet with sigma=4. The parameter `sigma`
%                 is the center frequency to standard deviation ratio of
%                 the main generating Gaussian distribution. Note that the
%                 true peak and standard deviation of the Morlet wavelet
%                 differ, in particular for low values of `sigma` (<5).
%                 This is a consequence of the correction factor. The
%                 parameter can be set directly by `{'morlet',sigma}`.
%
%    'fbsp'       Frequency B-spline wavelet of order m=3, center frequency
%                 to bandwidth ratio fb = 2. The parameters can be set
%                 by `{'fbsp',m,fb}`. Note that `m` must be integer and
%                 greater than or equal to 1, and `fb` must be greater than
%                 or equal to 2.
%
%    'analyticsp' Positive frequency part of cosine-modulated B-spline
%                 wavelet of  order order=3, with center frequency to main
%                 lobe width ratio fb = 1. The parameters can be set by
%                 `{'analyticsp',order,fb}`. Note that `order` and `fb`
%                 must be integer and greater than or equal to 1.
%
%    'cplxsp'     Complex-modulated B-spline of order order=3, with center
%                 frequency to main lobe width ratio fb = 1. The parameters
%                 can be set by `{'cplxsp',order,fb}`. Note that `order`
%                 and `fb` must be integer and greater than or equal to 1.
%
%   `freqwavelet(name,L,scale)` returns a "dilated" version of the wavelet.
%   The nonzero scalar `scale` controls both the bandwidth and the center
%   frequency. Values greater than 1 make the wavelet wider (and narrower
%   in the frequency domain), values lower than 1 make the wavelet narrower.
%   The center frequency is moved to`0.1/scale`. The center frequency is
%   limited to the range of "positive" frequencies ]0,1] (up to Nyquist rate).
%   If `scale` is a vector, the output is a `L x numel(scale)` matrix with
%   the individual wavelets as columns.
%
%   The following additional flags and key-value parameters are available:
%
%     'scale'           Wavelet scale (relative to basic scale)
%
%     'waveletParams'   a vector containing the respective wavelet parameters
%                       [alpha, beta, gamma] for cauchy and morse wavelets
%                       [sigma] for morlet wavelets
%                       [order, fb] for splines
%
%     'basefc',fc      Normalized center frequency of the mother wavelet
%                      (scale=1). The default is 0.1.
%
%     'bwthr',bwthr    The height at which the bandwidth is computed.
%                      The default value is 10^(-3/10) (~0.5).
%
%     'efsuppthr',thr  The threshold determining the effective support of
%                      the wavelet. The default value is 10^(-5).
%
%     'scal',s         Scale the filter by the constant *s*. This can be
%                      useful to equalize channels in a filter bank.
%
%     'delay',d        Set the delay of the filter. Can be either a scalar or
%                      a vector of the same length as *scale*. Default value is zero.
%
%   The admissible range of scales can be adjusted to handle different
%   scenarios:
%
%     'positive'       Enables the construction of wavelets at postive
%                      center frequencies ]0,1]. If `basefc=0.1`, this
%                      corresponds to scales larger than or equal to 0.1.
%                      This is the default.
%
%     'negative'       Enables the construction of wavelets at negative
%                      center frequencies [-1,0[. If `basefc=0.1`, this
%                      corresponds to scales smaller than or equal to -0.1.
%
%     'analytic'       Enables the construction of wavelets with center
%                      frequencies in ]0,2] for analysis of analytic signals.
%                      [! This feature is currently experimental and may not
%                      always work as intended. !]
%
%   The format of the output is controlled by the following flags:
%   'full' (default),'econ','asfreqfilter':
%
%     'full'           The output is a `L x numel(scale)` matrix.
%
%     'econ'           The output is a `numel(scale)` cell array with
%                      individual freq. domain wavelets truncated to the
%                      length of the effective support given by parameter 'efsuppthr'.
%                      Does not run stably for system lengths > 2000
%
%     'asfreqfilter'   As 'econ', but the elements of the cell-array are
%                      filter structs with fields .H and .foff as in
%                      |blfilter| to be used in |filterbank| and related.
%
%   `[H,info]=freqwavelet(...)` additionally returns a struct with the
%   following fields:
%
%     .fc             Normalized center frequency.
%
%     .foff           Index of the first sample above the effective
%                     support threshold (minus one).
%                     It can be directly used to convert the 'econ'
%                     output to 'full' by `circshift(postpad(H,L),foff)`.
%
%     .fsupp          Length of the effective support (with values above
%                     efsuppthr).
%
%     .basefc         Center frequency of the implied mother wavelet.
%
%     .scale          The scale used.
%
%     .dilation       The actual dilation used in the formula.
%
%     .bw             Relative bandwidth at -3 dB (half of the height).
%
%     .tfr            Time-frequency ratio of a Gaussian with the same
%                     bandwidth as the wavelet.
%
%     .aprecise       Exact natural subsampling factors (not rounded).
%
%     .a_natural      Fractional natural subsampling factors in the
%                     format acceptable by |filterbank| and related.
%
%     .cauchyAlpha    Alpha value of closest Cauchy wavelet [NOTE: Not
%                     implemented for non-Morse wavelets.]
%
%   Additionally, the function accepts flags to normalize the output.
%   Please see the help of |setnorm|. By default, no normaliazation is
%   applied.
%
%   Wavelet definitions
%   -------------------
%
%   `C` is a normalization constant.
%
%   Cauchy wavelet
%
%       H = C \xi^{\frac{\alpha-1}{2}} exp( -2\pi\xi )
%
%       .. math:: H = C \xi^{\frac{\alpha-1}{2}} exp( -2\pi\xi )
%
%   Morse wavelet
%
%       H = C \xi^{\frac{\alpha-1}{2\gamma}} exp( -2\pi\xi^{\gamma} )
%
%       .. math:: H = C \xi^{\frac{\alpha-1}{2\gamma}} exp( -2\pi\xi^{\gamma} )
%
%   Morlet wavelet
%
%       H = C \xi^{\frac{\alpha-1}{2\gamma}} exp( -2\pi\xi^{\gamma} )
%
%       .. math:: H = C \xi^{\frac{\alpha-1}{2\gamma}} exp( -2\pi\xi^{\gamma} )
%
%   Frequency bandlimited spline wavelet
%
%       H = C B (\xi - m \frac{\xi}{4})
%
%       .. math:: H = C B (\xi - m \frac{\xi}{4})
%
%   Analytic spline wavelet
%
%       H = C exp(-j \omega x) A(-exp(j \omega)) H(exp(-j \omega)
%
%       .. math:: H = C exp(-j \omega x) A(-exp(j \omega)) H(exp(-j \omega))
%
%   Complex spline wavelet
%
%       H = C exp(-j \omega x + \xi )
%
%       .. math:: H = C exp(-j \omega x + \xi )
%
%
%   References: rioul92 unalsc94
%
%   See also: setnorm, filterbank, blfilter

complainif_notenoughargs(nargin,2,upper(mfilename));

%set default input parameters
if ~isscalar(L)
    error('%s: L must be a scalar',upper(mfilename));
end

if ~iscell(name), name = {name}; end

freqwavelettypes = getfield(arg_freqwavelet(),'flags','wavelettype');

if ~ischar(name{1}) || ~any(strcmpi(name{1},freqwavelettypes))
  error('%s: First input argument must be the name of a supported window.',...
        upper(mfilename));
end

definput.import={'setnorm', 'freqwavelet'};
definput.importdefaults={'null'};
definput.keyvals.scale = 1;
definput.keyvals.scal = [];
definput.keyvals.basefc = 0.1;
definput.keyvals.delay=0;
definput.keyvals.fc = [];
definput.keyvals.bwthr = 10^(-3/10);
definput.keyvals.efsuppthr = 10^(-5);
definput.flags.freqrange = {'positive','negative','analytic'};
definput.flags.outformat = {'full','econ','asfreqfilter'};
definput.keyvals.fs = 2;
definput.keyvals.alphaStep = definput.keyvals.fs/L;

if ~iscell(name)
    name = {name};
end

[flags,kv,scale]=ltfatarghelper({'scale'},definput,varargin,'freqwavelet');

%check L
if ~isscalar(L)
    error('%s: L must be a scalar',upper(mfilename));
end

%check basefc
if ~isscalar(kv.basefc)
    error('%s: basefc must be a positive scalar',upper(mfilename));
end

scale = scale(:).'; % Scale should always be a row vector
if ~isnumeric(scale), error('%s: scale must be numeric',upper(mfilename)); end

if isempty(kv.scal)
    kv.scal = scale;
elseif ~isnumeric(kv.scal)
    error('%s: scal must be numeric',upper(mfilename));
elseif numel(kv.scal) ~= numel(scale)
    error('%s: scal must have exactly as many entries as scale',upper(mfilename));
end

if kv.alphaStep > 0.02
    error('%s: wavelet sampling too small. increase system length.',upper(mfilename));
end

%if L/scale > 10
%   error('%s: scale too large.',upper(mfilename));
%end

% Check range of scales
if flags.do_positive && (any(scale <= 0) || any(kv.basefc./scale > 1))
    error('%s: Frequency range flag is set to positive. scale must be positive and not smaller than basefc.', upper(mfilename));
end
if flags.do_negative && (any(scale >= 0) || any(kv.basefc./scale < -1))
    error('%s: Frequency range flag is set to negative. scale must be negative and not larger than -basefc.', upper(mfilename));
end
if flags.do_analytic && (any(scale <= 0) || any(kv.basefc./scale > 2))
    error('%s: Frequency range flag is set to analytic. scale must be positive and not smaller than 2*basefc.', upper(mfilename));
end

% Check delay vector
if numel(kv.delay) == 1
    kv.delay = repmat(kv.delay, 1, numel(scale));
end
if numel(kv.delay) ~= numel(scale)
    error('%s: delay must have exactly as many entries as scale',upper(mfilename));
end
% Check other parameters
if kv.efsuppthr < 0, error('%s: efsuppthr must be nonnegative',upper(mfilename)); end
if kv.bwthr < 0, error('%s: bwthr must be nonnegative',upper(mfilename)); end
if kv.bwthr < kv.efsuppthr, error('%s: efsuppthr must be lower than bwthr.',upper(mfilename)); end

M = numel(scale);

if flags.do_full
    H = zeros(L,M);
else
    H = cell(1,M);
end

if flags.do_negative
    wltype = 'negative';
else
    wltype = 'positive';
end

%% generate the wavelet prototype
[fun, fsupp_, peakpos, cauchyAlpha] = helper_waveletgeneratorfunc(name, wltype);

%% calculate the support as f(scale)
fsupp = zeros(5,M);
fsupp(5,:) = kv.fs;
basedil = peakpos/(kv.basefc);

if kv.efsuppthr > 0
    fsupp(1,:) = max(0,(fsupp_(1)/basedil)./scale);
    fsupp(5,:) = min(kv.fs,(fsupp_(5)/basedil)./scale);
end
fsupp(2,:) = max(0,(fsupp_(2)/basedil)./scale);
fsupp(3,:) = (fsupp_(3)/basedil)./scale;
fsupp(4,:) = min(kv.fs,(fsupp_(4)/basedil)./scale);

if flags.do_negative
    fsupp = flip(fsupp);
end

fsuppL = fsuppL_inner(fsupp,kv.fs,L,1:5);

%% calculate H
if flags.do_full
    if ~flags.do_negative
        y = ((0:L-1)').*basedil*kv.alphaStep*scale;
    else
        y = ([0;(L-1:-1:1)']).*basedil*kv.alphaStep*scale;
    end
    H = abs(kv.scal).*setnorm(fun(y), flags.norm);
elseif flags.do_econ
    for ii = 1:numel(scale)
        y = ((fsuppL(1,ii):fsuppL(end,ii)-1)').*basedil*kv.alphaStep*abs(scale(ii));
        H{ii} = abs(kv.scal(ii)).*setnorm(fun(y), flags.norm);%TODO: check output format: should be cell
    end
elseif flags.do_asfreqfilter
    for m = 1:numel(scale)
        y = @(L) ((fsuppL_inner(fsupp(:,m),kv.fs,L,1):fsuppL_inner(fsupp(:,m),kv.fs,L,5)-1)').*basedil*abs(scale(m))*kv.fs/L;
        H{m} = struct('H',@(L) abs(kv.scal(m))*setnorm(fun(y(L)),flags.norm),'foff',@(L)fsuppL_inner(fsupp(:,m),kv.fs,L,1),'realonly',0, 'delay', kv.delay(m));
   end
end

%% write info struct
info.basefc = kv.basefc;
info.fc = fsupp(3,:);
info.scale = scale;%';
info.dilation = basedil.*scale;%';
info.fsupp = fsuppL(end,:) - fsuppL(1,:) + ones(1,numel(scale));
if info.fsupp <= 0, info.fsupp = 0; end
info.bw  = (fsupp(4,:) - fsupp(2,:));
bwinsamples = info.bw./kv.alphaStep;
info.aprecise = L./bwinsamples;
info.a_natural(:,2) = ceil(bwinsamples);
info.a_natural = info.a_natural';
info.tfr = (cauchyAlpha - 1)./(pi*info.fc.^2*L);
info.cauchyAlpha = cauchyAlpha;
info.foff = fsuppL(1,:);

%if M==1 && iscell(H)
%    H = H{1};
%end

end

function fsuppL = fsuppL_inner(fsupp,fs,L,idx)
    fsuppL_all = [ ceil(fsupp(1:2,:)/fs*L); round(fsupp(3,:)/fs*L); floor(fsupp(4:5,:)/fs*L) ];
    fsuppL = fsuppL_all(idx,:);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function result = cvMcd(data,kmax,resMCD,h)

%CVMCD calculates the robust cross-validated PRESS (predicted residual error sum of squares)
% curve for the MCD method in a fast way.
%
% Input arguments:
%   data   : the full data set
%   kmax   : the maximal number of components to be considered (mostly kmax = p).
%   resMCD : the result of mcdcov(data,'plots',0,'factor',1)
%   h      : the quantile used in MCD.
%
% output:
%   result.press     : vector of length kmax with the press values
%   result.weights   : the weights for all observations
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by Sanne Engelen
% Last Update: 01/07/2004

% Some initialisations:
n = size(data,1);
p = size(data,2);
r = rank(data);
Pk = [];
Lk = [];
teller_if_lus = 0;

if nargin < 4
    alfa = 0.75;
    h=floor(2*floor((n+p+1)/2)-n+2*(n-floor((n+p+1)/2))*alfa);
end

outWeights = weightscvMcd(data,r,kmax,resMCD,h);

w_min = outWeights.w_min;

Hopt = resMCD.Hsubsets.Hopt;
inputH0.H0 = Hopt;
Tfull = mean(data(Hopt,:));
Sfull = cov(data(Hopt,:));

for i = 1:n
    % deciding which index should be removed from H0.
    inputH0.same = 0;
    if isempty(find(inputH0.H0 == i))
        inputH0.j = h;
        if teller_if_lus >= 1
            inputH0.same = 1;
        end
        teller_if_lus = teller_if_lus + 1;
    else
        inputH0.j = find(inputH0.H0 == i);
    end

    % assigning the input variables:
    inputFull.T = Tfull;
    inputFull.S = Sfull;

    if ~inputH0.same
        res = removeObsMcd(data,i,inputH0,inputFull);
    end
    if (isempty(find(inputH0.H0 == i))) & (teller_if_lus == 1)
        resfixed = res;
    end
    if isempty(find(inputH0.H0 == i)) & (teller_if_lus ~= 1)
        res = resfixed;
    end

    P_min_i = res.P_min_i;
    L_min_i = res.L_min_i;
    mu_min_i = res.mu_min_i;

     for k = 1:kmax
        clear Pk Lk;
        Pk = P_min_i(:,1:k);
        Lk = L_min_i(1:k,1:k);
        Xhoedk_min_i(i,(k-1)*p + 1:k*p) = (data(i,:) - mu_min_i)*Pk*Pk' + mu_min_i;

        if k~=r
            odk(i,k) = norm(data(i,:) - Xhoedk_min_i(i,(k-1)*p + 1:k*p));
        else
            odk(i,k) = 0;
        end
    end
end

for k = 1:kmax
    press_min(k) = 1/sum(w_min)*w_min*odk(:,k).^2;
end

result.press = press_min;
result.weights = outWeights;

%----------------------------------------------------------------------------------
function out = weightscvMcd(data,r,kmax,resMCD,h)

% computes the weights used to calculate the robust PRESS values.
%
% input:
%   data   : the whole data
%       r  : the rank of the data
%   kmax   : the maximal number of components to be considered
%   resMCD : the result of mcdcov(data,'plots',0,'factor',1)
%   h      : the number of observations on which the computations are based.
%
% output:
%   out.w_min    : the weights computed by taken the minimum over all k

% Some initialisations:
n = size(data,1);
p = size(data,2);
Pk = [];
Lk = [];
Tik = [];

[P,L] = eig(resMCD.cov);
[L,I] = greatsort(diag(L));
P = P(:,I);

for k = 1:kmax
    Pk = P(:,1:k);
    if h==n
        Lk=L(1:k);
    else
    Lk = chi2inv(h/n,k)/chi2inv(h/n,kmax/2)*L(1:k);% with correction for the factor
    end
    muk = resMCD.center;
    Xhoedk(:,(k-1)*p + 1:k*p) = (data - repmat(muk,n,1))*Pk*Pk' + repmat(muk,n,1);
    Tk = (data - repmat(muk,n,1))*Pk;

    for i =1:n
        % defining the sd for the observation that is left:
        sdk(i,k) = sqrt(libra_mahalanobis(Tk(i,:),zeros(1,k),'cov',diag(Lk)));

        % defining the od for the observation that is left:
        if k~=r
            odk(i,k) = norm(data(i,:) - Xhoedk(i,(k-1)*p + 1:k*p));
        else
            odk(i,k) = 1;
        end
    end

    % defining weights for odk and sdk:
    if k~=r
        [m,s]=unimcd(odk(:,k).^(2/3),h);
        cutoff(k)=sqrt(norminv(0.975,m,s).^3);
        wod(:,k) = (odk(:,k) <= cutoff(k));
    else
        cutoff(k)= 0;
        wod(:,k) = 1;
    end
    wsd(:,k) = (sdk(:,k) <= sqrt(chi2inv(0.975,k)));
end

% determine the weights for every observation:
wk = wsd & wod;

if size(wk,1) == 1 | size(wk,2) == 1
    w_min = wk';
else
    w_min = min(wk');
end

out.w_min = w_min;

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function out = slidingavg(in, N)
%   OUTPUT_ARRAY = SLIDINGAVG(INPUT_ARRAY, N)
%
%  The function 'slidingavg' implements a one-dimensional filtering, applying a sliding window to a sequence. Such filtering replaces the center value in
%  the window with the average value of all the points within the window. When the sliding window is exceeding the lower or upper boundaries of the input
%  vector INPUT_ARRAY, the average is computed among the available points. Indicating with nx the length of the the input sequence, we note that for values
%  of N larger or equal to 2*(nx - 1), each value of the output data array are identical and equal to mean(in).
%
%  *  The input argument INPUT_ARRAY is the numerical data array to be processed.
%  *  The input argument N  is the number of neighboring data points to average over for each point of IN.
%
%  *  The output argument OUTPUT_ARRAY is the output data array.
%
%     Š 2002 - Michele Giugliano, PhD and Maura Arsiero
%     (Bern, Friday July 5th, 2002 - 21:10)
%    (http://www.giugliano.info) (bug-reports to michele@giugliano.info)
%
% Two simple examples with second- and third-order filters are
%  slidingavg([4 3 5 2 8 9 1],2)
%  ans =
%   3.5000  4.0000  3.3333  5.0000  6.3333  6.0000  5.0000
%
%  slidingavg([4 3 5 2 8 9 1],3)
%  ans =
%   3.5000  4.0000  3.3333  5.0000  6.3333  6.0000  5.0000
%

if (isempty(in)) | (N<=0)                                              % If the input array is empty or N is non-positive,
 disp(sprintf('SlidingAvg: (Error) empty input data or N null.'));     % an error is reported to the standard output and the
 return;                                                               % execution of the routine is stopped.
end % if

if (N==1)                                                              % If the number of neighbouring points over which the sliding
 out = in;                                                             % average will be performed is '1', then no average actually occur and
 return;                                                               % OUTPUT_ARRAY will be the copy of INPUT_ARRAY and the execution of the routine
end % if                                                               % is stopped.

nx   = length(in);             % The length of the input data structure is acquired to later evaluate the 'mean' over the appropriate boundaries.

if (N>=(2*(nx-1)))                                                     % If the number of neighbouring points over which the sliding
 out = mean(in)*ones(size(in));                                        % average will be performed is large enough, then the average actually covers all the points
 return;                                                               % of INPUT_ARRAY, for each index of OUTPUT_ARRAY and some CPU time can be gained by such an approach.
end % if                                                               % The execution of the routine is stopped.

out = zeros(size(in));         % In all the other situations, the initialization of the output data structure is performed.

if rem(N,2)~=1                 % When N is even, then we proceed in taking the half of it:
 m = N/2;                      % m = N     /  2.
else                           % Otherwise (N >= 3, N odd), N-1 is even ( N-1 >= 2) and we proceed taking the half of it:
 m = (N-1)/2;                  % m = (N-1) /  2.
end % if

for i=1:nx,                                                 % For each element (i-th) contained in the input numerical array, a check must be performed:
 if ((i-m) < 1) & ((i+m) <= nx)                             % If not enough points are available on the left of the i-th element..
  out(i) = mean(in(1:i+m));                                 % then we proceed to evaluate the mean from the first element to the (i + m)-th.
 elseif ((i-m) >= 1) & ((i+m) <= nx)                        % If enough points are available on the left and on the right of the i-th element..
  out(i) = mean(in(i-m:i+m));                               % then we proceed to evaluate the mean on 2*m elements centered on the i-th position.
elseif ((i-m) >= 1) & ((i+m) > nx)                          % If not enough points are available on the rigth of the i-th element..
  out(i) = mean(in(i-m:nx));                                % then we proceed to evaluate the mean from the element (i - m)-th to the last one.
 elseif ((i-m) < 1) & ((i+m) > nx)                          % If not enough points are available on the left and on the rigth of the i-th element..
  out(i) = mean(in(1:nx));                                  % then we proceed to evaluate the mean from the first element to the last.
 end % if
end % for i


text:
function bz_RunAnimalMetadata(basepath)

%% Initial variable parsing
if ~exist('basepath','var')
    basepath = cd;
elseif isempty(basepath)
    basepath = cd;
end
basename = bz_BasenameFromBasepath(basepath);

%%
bz_EditAnimalMetadata(basepath)

%%
notesname = [basename,'_AnimalMetadataText'];
notesfullpath = fullfile(basepath,notesname);
prompt = 'Push any key in this window when done editing the NoteText file ';
str = input(prompt,'s');
run(notesfullpath);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [mCat, mCatSmStdr] = calc_fpfitqual(mCatalog, fSm, fStdr, fUncertMin, fUncertMax)
% [mCat, mCatSmStdr] = calc_fpfitqual(mCatalog,fSm, fStdr, fUncertMin, fUncertMax)
% --------------------------------------------------------------------------------
% Create quality restricted catalog using quality parameters from FPFIT
% Reference: Oppenheimer & Reasenberg, Open-File Report 85-739, 1985
% Based on catalogs imported using ncedc_fomech_imp.m
%
% Input variables
% mCatalog   : Focal meachnism catalog
% fSm        : Maximum solution misfit
% fStdr      : Minimum station distribution ratio
% fUncertMin : Min uncertainty for strike (dip direction), rake, dip
% fUncertMax : Max uncertainty for strike (dip direction), rake, dip
%
% Output:
% mCat       : Quality catalog selected according to the input
% mCatSmStdr : Quality only selected for Stdr and SM
%
% J. Woessner, 04.2004

% Selection of solution misfit and station distribution ratio
vSel = (mCatalog(:,15) >= fStdr & mCatalog(:,13) <= fSm);
mCatSmStdr = mCatalog(vSel,:);

% Selection using the uncertainties for strike, dip and slip
vSel = (mCatSmStdr(:,17) >= fUncertMin & mCatSmStdr(:,17) < fUncertMax & mCatSmStdr(:,18) >= fUncertMin & mCatSmStdr(:,18) < fUncertMax & mCatSmStdr(:,19) >= fUncertMin & mCatSmStdr(:,19) < fUncertMax);
mCat = mCatSmStdr(vSel,:);

author: CelsoReyes


text:
clear;
clc;

ori_root = './train/';
ori_images = dir(fullfile([ori_root 'ShadowImages/'],'*.jpg'));
ori_gts = dir(fullfile([ori_root 'ShadowMasks/'],'*.png'));

save_path = './train_agru/';

for i=1: length(ori_images)
    i
    ro_degree = 0;%rand(1,1)*360-180;

    w_e = rand(1,1)*-0.3;
    h_e = rand(1,1)*-0.3;

    %  b_h = rand(1,1)*90+10;
    % b_w = rand(1,1)*90+10;

    image = imread([ori_root 'ShadowImages/' ori_images(i).name]);
    gt = imread([ori_root 'ShadowMasks/' ori_gts(i).name]);

    [h,w,~]=size(image);

    new_im = imrotate(image,ro_degree);
    new_gt = imrotate(gt,ro_degree);
    [new_h,new_w,~]=size(new_im);

    if (ro_degree<=45 && ro_degree>=-45) || (ro_degree>=135 && ro_degree<=-135)
        f_new_im = new_im(max(new_h/2-h/2,1):min(h/2+new_h/2,size(new_im,1)),max(new_w/2-w/2,1):min(w/2+new_w/2,size(new_im,2)),:);
        f_new_gt = new_gt(max(new_h/2-h/2,1):min(h/2+new_h/2,size(new_im,1)),max(new_w/2-w/2,1):min(w/2+new_w/2,size(new_im,2)),:);
    else
        f_new_im = new_im(max(new_w/2-w/2,1):min(w/2+new_w/2,size(new_im,1)),max(new_h/2-h/2,1):min(h/2+new_h/2,size(new_im,2)),:);
        f_new_gt = new_gt(max(new_w/2-w/2,1):min(w/2+new_w/2,size(new_im,1)),max(new_h/2-h/2,1):min(h/2+new_h/2,size(new_im,2)),:);
    end

%     f_new_im = image;
%     f_new_gt = gt;

    f_new_im = imresize(f_new_im, [h*(1+h_e), w*(1+w_e)]);
    f_new_gt = imresize(f_new_gt, [h*(1+h_e), w*(1+w_e)]);

  %   figure(2),imshow(f_new_im);
% %
  %   figure(1),imshow(f_new_gt);

%%
    imwrite(f_new_im,[save_path 'ShadowImages/a_' ori_images(i).name]);
    imwrite(f_new_gt,[save_path 'ShadowMasks/a_' ori_gts(i).name]);

end

author: xw-hu

repo: DSC

save_path: github-repos/MATLAB/xw-hu-DSC


text:
function images = rmfilter_energy(images, display)

% simlpe energy filter:
%   (1) subtract background intenstity
%   (2) rectify
%   (3) divide by maximum
%       (max = greatest absolute distance from background)

try
    % look up min, max, and bk from screen calibration
    Cmap.min  = min(display.stimRgbRange);
    Cmap.max  = max(display.stimRgbRange);
    Cmap.bk   = display.backColorIndex;
catch
    % if not found assume them from image matrix
    Cmap.min  = min(images(:));
    Cmap.max  = max(images(:));
    Cmap.bk   = round(mean(images(:)));
end

images = double(images);
images = images - Cmap.bk;
images = sqrt(images.^2);
maxEnergy = double(max(Cmap.max - Cmap.bk, Cmap.bk - Cmap.min));
maxEnergy = double(maxEnergy);
images = images  / maxEnergy;

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function s = logspace_len(a, b, n, l)

if nargin < 4
    l = n
end

r = exp((log(10^b)-log(10^a)) / (n-1));
s = 10^a * r.^(0:l-1);

end

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet


text:
%wrap - Set radian angles in range [0,2pi] or [-pi,pi].
%
%  USAGE
%
%    y = wrap(x,offset)
%
%    x              angles in radians
%    offset         phase angle [0-2*pi] that circularly shifts output data
%
%
%  SEE ALSO
%
%    See also isradians, clinspace.
%

% Copyright (C) 2010-2011 by Michaël Zugaro
% edited by David Tingley, 2017 to add phase offset as input
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

function y = bz_wrap(x,offset)

% Check number of parameters
if nargin < 1,
  error('Incorrect number of parameters (type ''help wrap'' for details).');
end

% if nargin < 2,
% 	range = 1;
% end

if ~isa(x,'double'), y = []; return; end

% Determine angle in [0,2*pi]
y = mod(mod(x,2*pi)+offset,2*pi);

% Change range if necessary
% if range == 1,
% 	change = y > pi;
% 	y(change) = y(change)-2*pi + offset;
% end


text:
%%
% Auther : Gaul Swapnil Narhari,
% Developed @ IIT Kharagpur, India.
% In collaborationwith TaraNG, India.
%%
% for optimization see 'pso.m' file for fitness function defination see 'simple_fitness.m'
%%
clc;
clear all;
ObjectiveFunction = @simple_fitness;
nvars = 3;    % Number of variables
LB = [0 0 0];   % Lower bound
UB = [3 5 10];  % Upper bound
Popln=20;
Genrtn=50;
WByn=1;
[x,fvalue] = pso(ObjectiveFunction,'NumVar',nvars,'LowerBound',LB,'UpperBound',UB,'Population',Popln,'Genrations',Genrtn,'IncludeWB',WByn);
disp(x);
disp(fvalue);


text:
function [surf, labels] = imSurfaceEstimate(img, varargin)
% Estimate surface area of a binary 3D structure
%
%   Deprecated, use 'imSurfaceAreaEstimate' instead.
%
%   Usage
%   S = imSurfaceEstimate(IMG)
%   Estimate the surface area of the structure within the image, without
%   measuring surface area of borders.
%   The aim of this function is to be called by the "imSurfaceDensity"
%   function, for providing an estimate of surface area density within a
%   representative volume of interest.
%
%
%   Example
%   imSurfaceEstimate
%
%   See also
%   imSurfaceAreaEstimate, imSurfaceDensity
%

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

%% Process input arguments

warning('MatImage:deprecated', ...
    'function imSurfaceEstimate is obsolete, use imSurfaceAreaEstimate instead');

% check image dimension
if ndims(img) ~= 3
    error('first argument should be a 3D image');
end

% in case of a label image, return a vector with a set of results
if ~islogical(img)
    labels = unique(img);
    labels(labels==0) = [];
    surf = zeros(length(labels), 1);
    for i=1:length(labels)
        surf(i) = imSurfaceEstimate(img==labels(i), varargin{:});
    end
    return;
end

%% Process input arguments

% in case of binary image, compute only one label...
labels = 1;

% default number of directions
nDirs = 13;

% default image resolution
delta = [1 1 1];

% Process user input arguments
while ~isempty(varargin)
    var = varargin{1};
    if ~isnumeric(var)
        error('option should be numeric');
    end

    % option is either connectivity or resolution
    if isscalar(var)
        nDirs = var;
    else
        delta = var;
    end
    varargin(1) = [];
end

%% Use the LUT to estimate surface area

lut = imSurfaceLut(delta, nDirs);
bch = imBinaryConfigHisto(img);
surf = sum(bch .* lut);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function ADE = estimate_average_displacement_error(displacementField1,displacementField2,mask)
% ESTIMATE_AVERAGE_DISPLACEMENT_ERROR Estimates the average displacement error
%
% ADE = estimate_average_displacement_error(displacementField1,displacementField2,mask)
%
% INPUT ARGUMENTS
% displacementField1    - Displacement field 1
% displacementField2    - Displacement field 2
% mask                  - Mask, to remove certain points in the data
%
% OPTIONAL INPUT ARGUMENTS
% N/A
%
% OUTPUT ARGUMENTS
% ADE                   - Average displacement error between the two
%                         provided displacement fields

% Copyright (c) 2012 Daniel Forsberg
% danne.forsberg@outlook.com
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

dims = length(displacementField1);
numberOfElements = sum(mask(:) ~= 0);
switch dims
    case 1
        x_error = mask.*(displacementField1{1}-displacementField2{1});
        ADE = sum(sqrt(x_error(:).^2))/numberOfElements;
    case 2
        x_error = mask.*(displacementField1{1}-displacementField2{1});
        y_error = mask.*(displacementField1{2}-displacementField2{2});
        ADE = sum(sqrt(x_error(:).^2 + y_error(:).^2))/numberOfElements;
    case 3
        x_error = mask.*(displacementField1{1}-displacementField2{1});
        y_error = mask.*(displacementField1{2}-displacementField2{2});
        z_error = mask.*(displacementField1{3}-displacementField2{3});
        ADE = sum(sqrt(x_error(:).^2 + y_error(:).^2 + z_error(:).^2))/numberOfElements;
    otherwise
        error('Not implemented for dimensions higher than three.')
end

author: fordanic

repo: image-registration


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym polylog (@var{s}, @var{z})
%% Symbolic polylogarithm function.
%%
%% Returns the polylogarithm of order @var{s} and argument @var{z}.
%%
%% Example:
%% @example
%% @group
%% syms z s
%% polylog(s, z)
%%   @result{} ans = (sym) polylog(s, z)
%% diff(ans, z)
%%   @result{} (sym)
%%       polylog(s - 1, z)
%%       ─────────────────
%%               z
%% @end group
%% @end example
%%
%% The polylogarithm satisfies many identities, for example:
%% @example
%% @group
%% syms s positive
%% polylog (s+1, 1)
%%   @result{} (sym) ζ(s + 1)
%% zeta (s+1)
%%   @result{} (sym) ζ(s + 1)
%% @end group
%% @end example
%%
%% @seealso{@@sym/dilog, @@sym/zeta}
%% @end defmethod

function L = polylog(s, z)
  if (nargin ~= 2)
    print_usage ();
  end

  L = elementwise_op ('polylog', sym(s), sym(z));
end

%!assert (isequal (polylog (sym('s'), 0), sym(0)))

%!assert (isequal (double (polylog (1, sym(-1))), -log(2)))

%!assert (isequal (double (polylog (0, sym(2))), -2))
%!assert (isequal (double (polylog (-1, sym(2))), 2))
%!assert (isequal (double (polylog (-2, sym(3))), -1.5))
%!assert (isequal (double (polylog (-3, sym(2))), 26))
%!assert (isequal (double (polylog (-4, sym(3))), -15))

%!assert (isequal (double (polylog (1, sym(1)/2)), log(2)))

%!test
%! % round trip
%! syms s z
%! f = polylog (s, z);
%! h = function_handle (f, 'vars', [s z]);
%! A = h (1.1, 2.2);
%! B = polylog (1.1, 2.2);
%! assert (A, B)

repo: octsympy


text:
function om = add_quad_cost(om, name, idx, Q, c, k, varsets)
%ADD_QUAD_COST  Adds a set of user costs to the model.
%   OM.ADD_QUAD_COST(NAME, Q, C);
%   OM.ADD_QUAD_COST(NAME, Q, C, K);
%   OM.ADD_QUAD_COST(NAME, Q, C, K, VARSETS);
%   OM.ADD_QUAD_COST(NAME, IDX_LIST, Q, C);
%   OM.ADD_QUAD_COST(NAME, IDX_LIST, Q, C, K);
%   OM.ADD_QUAD_COST(NAME, IDX_LIST, Q, C, K, VARSETS);
%
%   Adds a named block of quadratic costs to the model. Costs are of the
%   form
%       F(X) = 1/2 * X'*Q*X + C'*X + K
%   where Q is an NX x NX matrix (possibly sparse), C is an NX x 1 vector,
%   K is a scalar and NX is the number of elements in X. Here X is the vector
%   formed by combining the specified VARSETS (the full optimization vector
%   by default). Alternatively, if Q is an NX x 1 vector or empty, then F(X)
%   is also NX x 1, and K can be either NX + 1 or scalar.
%       F(X) = 1/2 * Q .* X.^2 + C .* X + K
%
%   Indexed Named Sets
%       A cost set can be identified by a single NAME, as described
%       above, such as 'PgCost', or by a name that is indexed by one
%       or more indices, such as 'PgCost(3,4)'. For an indexed named
%       set, before adding the cost sets themselves, the dimensions
%       of the indexed set must be set by calling INIT_INDEXED_NAME.
%
%       The constraints are then added using the following, where
%       all arguments are as described above, except IDX_LIST is a cell
%       array of the indices for the particular cost set being added.
%
%       OM.ADD_QUAD_COST(NAME, IDX_LIST, Q, C, K);
%       OM.ADD_QUAD_COST(NAME, IDX_LIST, Q, C, K, VARSETS);
%
%   Examples:
%       om.add_quad_cost('quad_cost1', Q1, c1, 0);
%       om.add_quad_cost('lin_cost2',  [], c2, k2, {'Vm', 'Pg', 'z'});
%
%       om.init_indexed_name('c', {2, 3});
%       for i = 1:2
%         for j = 1:3
%           om.add_quad_cost('c', {i, j}, Q{i,j}, ...);
%         end
%       end
%
%   See also OPT_MODEL, PARAMS_QUAD_COST, EVAL_QUAD_COST.

%   MP-Opt-Model
%   Copyright (c) 2008-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

%% initialize input arguments
if iscell(idx)          %% indexed named set
    if nargin < 7
        varsets = {};
    end
else                    %% simple named set
    if nargin < 6
        varsets = {};
    else
        varsets = k;
    end
    if nargin < 5
        k = 0;
    else
        k = c;
    end
    c = Q;
    Q = idx;
    idx = {};
end

%% convert varsets from cell to struct array if necessary
varsets = om.varsets_cell2struct(varsets);
nv = om.varsets_len(varsets);   %% number of variables

%% check sizes
[MQ, NQ] = size(Q);
[Mc, Nc] = size(c);
if MQ
    if NQ ~= MQ && NQ ~= 1
        error('@opt_model/add_quad_cost: Q (%d x %d) must be square or a column vector (or empty)', MQ, NQ);
    end
end
if Mc && Nc ~= 1
    error('@opt_model/add_quad_cost: c (%d x %d) must be a column vector (or empty)', Mc, Nc);
end
if MQ
    if Mc && Mc ~= MQ
        error('@opt_model/add_quad_cost: dimensions of Q (%d x %d) and c (%d x %d) are not compatible', MQ, NQ, Mc, Nc);
    end
    nx = MQ;
else
    if ~Mc
        error('@opt_model/add_quad_cost: Q and c cannot both be empty');
    end
    nx = Mc;
end
if nx ~= nv
    error('@opt_model/add_quad_cost: dimensions of Q (%d x %d) and c (%d x %d) do not match\nnumber of variables (%d)\n', MQ, NQ, Mc, Nc, nv);
end

%% size of named cost set
if NQ == 1 || isempty(Q)    %% if Q is a column vector or empty
    N = nx;                 %%   cost is element-wise, i.e. a vector
else                        %% otherwise Q is a square matrix
    N = 1;                  %%   cost is scalar
end

%% add the named quadratic cost set
om.add_named_set('qdc', name, idx, N, Q, c, k, varsets);

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
classdef  VideoMaker_TopOpt_densityBridge < VideoMaker_TopOpt
  properties

  end

  methods (Access = public)

        function Make_video_design_variable(obj,output_video_name)
            post = Postprocess_TopOpt_density();
            field2print = post.density_name;
            componentfield = post.density_name_component;
            obj.Make_video_density_field(field2print,componentfield,output_video_name)
        end

        function Make_video_design_variable_reg(obj,output_video_name)
            post = Postprocess_TopOpt_density();
            field2print = post.density_name_reg;
            componentfield = post.density_name_component_reg;
            obj.Make_video_density_field(field2print,componentfield,output_video_name)
        end

  end

  methods (Access = private)
        function Make_video_density_field(obj,field2print,componentfield,output_video_name)
            file_tcl_name = 'tcl_gid.tcl';
            file_list = obj.create_file_list(obj.iterations_to_print,obj.file_name,obj.files_folder);
            [output_video_name] = obj.replace_special_character(output_video_name);
            [output_photo]=obj.replace_special_character(strcat(obj.files_folder,'\',obj.file_name,'.png'));
            [file_list] = obj.replace_special_character(file_list);

            file_tcl_name_with_path = fullfile(obj.files_folder,file_tcl_name);
            file_path_in = fullfile(pwd,'FEM','PostProcess','Make_Video_densityBridge.tcl');
            filepath = obj.replace_special_character(file_path_in);

            fid = fopen(file_tcl_name_with_path,'w+');

            fprintf(fid,'GiD_Process PostProcess \n');
            fprintf(fid,'%s\n',['set arg1 "',file_list,'"']);
            fprintf(fid,'%s\n',['set arg2 "',output_video_name,'"']);
            fprintf(fid,'%s\n',['set arg3 "',field2print,'"']);
            fprintf(fid,'%s\n',['set arg4 "',componentfield,'"']);
            fprintf(fid,'%s\n',['set arg5 "',output_photo,'"']);

            fprintf(fid,'%s\n',['source "',filepath,'"']);
            fprintf(fid,'%s\n',['Make_Video_density $arg1 $arg2 $arg3 $arg4 $arg5']);
            fprintf(fid,'%s\n',['GiD_Process Mescape Quit']);
            fclose(fid);
            obj.execute_tcl_files(obj.gidPath,file_tcl_name_with_path)
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef prtUiTableEditStructure < prtUiManagerPanel
    % prtUiTableEditStructure
    %

    properties

        enableRowDeletion = 'on';
        enableColumnCreation = 'on';
        dataStructure
        dataCell = {};

        tableFieldNames = {};
        handleStruct

        sortingInds

        fontSizeTable = 8;
        fontSizeEdit = 12;

        editableArray = [];
    end

    methods
        function self = prtUiTableEditStructure(varargin)

            self = prtUtilAssignStringValuePairs(self,varargin{:});

            if nargin == 0
                self.create();
            elseif nargin~=0 && ~self.hgIsValid
               self.create()
            end

            init(self);
        end

        function create(self)
            % Overload the default prtUiManagerPanel
            % This will "create" a bigger figure and turn off toolbars and
            % stuff.

            self.handleStruct.figure = gcf;

            set(self.handleStruct.figure,...
                'Toolbar','none',...
                'MenuBar','none',...
                'Name','prtDataSetSelector',...
                'NumberTitle','off',...
                'WindowStyle','Normal',...
                'DockControls','off');

            if isempty(self.parent)
                self.managedHandle = uipanel(self.handleStruct.figure,...
                    'BackgroundColor',get(0,'DefaultFigureColor'),...
                    'BorderType','none');
            else
                self.managedHandle = uipanel('Parent',self.parent,...
                    'BackgroundColor',get(0,'DefaultFigureColor'),...
                    'BorderType','none');
            end
        end
        function init(self)

            self.handleStruct.table = uitable('parent',self.managedHandle,...
                'units','normalized',...
                'position',[0.05 0.05 0.9 0.8],...
                'fontUnits','points',...
                'fontSize',self.fontSizeTable,...
                'visible','on');
            drawnow;
            self.handleStruct.jScrollPane = findjobj(self.handleStruct.table);
            if ~isempty(self.handleStruct.jScrollPane)
                self.handleStruct.jTable = self.handleStruct.jScrollPane(1).getViewport.getView;
            end
            set(self.handleStruct.table,'visible','on');

            self.handleStruct.tableContextMenu = uicontextmenu;
            self.handleStruct.tableContextMenuItemNewSelection = uimenu(self.handleStruct.tableContextMenu, 'Label', 'Sort By', 'Callback', @(myHandle,eventData)self.sortBy(myHandle,eventData));
            self.handleStruct.tableContextMenuItemAndSelection = uimenu(self.handleStruct.tableContextMenu, 'Label', 'Set Values', 'Callback', @(myHandle,eventData)self.setValues(myHandle,eventData));
            self.handleStruct.tableContextMenuItemOrSelection = uimenu(self.handleStruct.tableContextMenu, 'Label', 'String Search Sort', 'Callback', @(myHandle,eventData)self.sortByStringSearch(myHandle,eventData));
            self.handleStruct.tableContextMenuCreateColumn = uimenu(self.handleStruct.tableContextMenu, 'Label', 'Add Column', 'Callback', @(myHandle,eventData)self.createColumn(myHandle,eventData),'Enable',self.enableColumnCreation,'Separator','on');
            self.handleStruct.tableContextMenuDeleteSelection = uimenu(self.handleStruct.tableContextMenu, 'Label', 'Delete Current Rows', 'Callback', @(myHandle,eventData)self.deleteSelectedRows(myHandle,eventData),'Enable',self.enableRowDeletion);

            %             self.handleStruct.tableContextMenuItemNoSelection = uimenu(self.handleStruct.tableContextMenu, 'Label', 'Remove From Table', 'Callback', @(myHandle,eventData)self.uiMenuNoSelection(myHandle,eventData));
            self.updateDataStructure;
        end

        function updateDataStructureFromGui(self,h,e)

            cellData = get(self.handleStruct.table,'data');

            selected = get(self.handleStruct.table,'UserData');
            if isempty(selected)
                warndlg('No rows selected');
                return;
            end
            selectedRow = selected(1);
            selectedCol = selected(2);

            obsInfo = self.dataStructure;
            obsInfo(self.sortingInds(selectedRow)).(self.tableFieldNames{selectedCol}) = cellData{selectedRow,selectedCol};
            self.dataStructure = obsInfo;

            self.updateDataStructure(obsInfo);
            self.updateTableDisplay;

            jUIScrollPane = findjobj(self.handleStruct.table);
            jUITable = jUIScrollPane.getViewport.getView;
            jUITable.changeSelection(selectedRow,selectedCol-1, false, false);
        end

        function appendDataStructure(self,obsInfo)
            if ~isempty(self.dataStructure)
                self.updateDataStructure(cat(1,obsInfo(:),self.dataStructure(:)));
            else
                self.updateDataStructure(obsInfo(:));
            end
        end

        function updateDataStructure(self,obsInfo)
            if nargin > 1
                obsInfo = obsInfo(:);
                self.dataStructure = obsInfo;
            end
            if ~isempty(self.dataStructure)
                self.tableFieldNames = fieldnames(self.dataStructure(:));
                self.dataCell = struct2cell(self.dataStructure(:))';

                if isempty(self.editableArray)
                    self.editableArray = false(1,size(self.dataCell,2));
                end

                set(self.handleStruct.table,'data',self.dataCell,...
                    'ColumnName',self.tableFieldNames,...
                    'ColumnEditable',self.editableArray,...
                    'ColumnWidth','auto','RearrangeableColumns','off',...
                    'RowStriping','off','RowName',cellstr(num2str((1:size(self.dataCell,1))')),...
                    'uicontextmenu',self.handleStruct.tableContextMenu,...
                    'CellSelectionCallback',@(src,evnt)set(src,'UserData',evnt.Indices),...
                    'CellEditCallback',@(src,event)self.updateDataStructureFromGui(src,event));

                if length(self.sortingInds) ~= numel(self.dataStructure)
                    self.sortingInds = [];
                end
                % I do not think we have to do this
                if isempty(self.sortingInds)
                    self.sortingInds = (1:size(self.dataCell,1))';
                end

                set(self.managedHandle,'ResizeFcn',@self.resizeFunction);
                resizeFunction(self);
                set(self.handleStruct.table,'visible','on');
            else
                set(self.handleStruct.table,'data',[],...
                    'ColumnName',self.tableFieldNames,...
                    'ColumnEditable',self.editableArray,...
                    'ColumnWidth','auto','RearrangeableColumns','off',...
                    'RowStriping','off','RowName',cellstr(num2str((1:size(self.dataCell,1))')),...
                    'uicontextmenu',self.handleStruct.tableContextMenu,...
                    'CellSelectionCallback',@(src,evnt)set(src,'UserData',evnt.Indices),...
                    'CellEditCallback',@(src,event)error('Not implemented'));

                self.sortingInds = [];

                set(self.managedHandle,'ResizeFcn',@self.resizeFunction);
                resizeFunction(self);
                set(self.handleStruct.table,'visible','on');
            end
        end

        function setValues(self,e,s)

            selected = get(self.handleStruct.table,'UserData');
            if isempty(selected)
                warndlg('No rows selected');
                return;
            end
            selectedRows = selected(:,1);

            strs = fieldnames(self.dataStructure);
            strOut = prtUiListDlg(strs);
            selectedColInd = strOut.enteredSelection;
            if isempty(selectedColInd) || selectedColInd == 0;
                warndlg('No field selected');
                return;
            end
            selectedColData = self.dataCell(:,selectedColInd);
            answer = inputdlg(sprintf('Please input a new value for field "%s" in the selected cells',strs{selectedColInd}));
            answer = answer{1};
            if isempty(answer)
                warndlg('No value entered');
                return;
            end
            if isnumeric(selectedColData{1})
                answer = str2double(answer);
                if isempty(answer)
                    warndlg('Improper numeric value provided');
                    return;
                end
            elseif islogical(selectedColData{1})
                answer = logical(str2double(answer));
                if isempty(answer)
                    warndlg('Improper numeric value provided');
                    return;
                end
            elseif isa(selectedColData{1},'char')
                if ~isa(answer,'char');
                    answer = num2str(answer);
                end
            end
            relativeInds = self.sortingInds(selectedRows);
            if ~isempty(answer)
                self.dataCell(relativeInds,selectedColInd) = {answer};
                for i = 1:length(relativeInds)
                    self.dataStructure(relativeInds(i)).(strs{selectedColInd}) = answer;
                end
            end

            self.updateTableDisplay;
        end

        function sortByStringSearch(self,e,s)
            strs = fieldnames(self.dataStructure);
            strOut = prtUiListDlg(strs);
            selectedColInd = strOut.enteredSelection;
            if isempty(selectedColInd) || selectedColInd == 0;
                return;
            end
            sortVecData = self.dataCell(:,selectedColInd);
            if isnumeric(sortVecData{1}) || islogical(sortVecData{1});
                warndlg('Numeric or logical field specified');
                return;
            end
            str = inputdlg('Input a search string, use * for wildcards');
            regexpStr = strrep(lower(str),'*','.*');

            found = regexp(lower(sortVecData),regexpStr);
            sortVecBool = cellfun(@(x)~isempty(x),found);
            [~,inds] = sort(~sortVecBool);
            self.sort(inds);
        end

        function sortBy(self,e,s)
            strs = fieldnames(self.dataStructure);
            strOut = prtUiListDlg(strs);
            selectedColInd = strOut.enteredSelection;
            if isempty(selectedColInd) || selectedColInd == 0;
                return;
            end
            sortVec = self.dataCell(:,selectedColInd);
            if isnumeric(sortVec{1}) || islogical(sortVec{1});
                sortVec = cat(1,sortVec{:});
            end
            [~,inds] = sort(sortVec);
            self.sort(inds);
        end

        function createColumn(self,h,e)

            answer=inputdlg('Enter a name for the new column','New Column Name');
            if isempty(answer)
                return;
            end
            varName = genvarname(answer{1});
            obsInfo = self.dataStructure;
            if isfield(obsInfo,varName);
                h = errordlg('%s is already a field of the database');
                uiwait(h);
                return
            end

            ButtonName = questdlg('Add a numeric or character?', ...
                'Numeric or Character', ...
                'Numeric', 'Character','Numeric');

            switch ButtonName
                case 'Numeric'
                    vals = repmat({nan},size(obsInfo));
                otherwise
                    vals = repmat({''},size(obsInfo));
            end
            for i = 1:length(obsInfo)
                obsInfo(i).(varName) = vals{i};
            end

            self.editableArray = [self.editableArray,true];
            self.updateDataStructure(obsInfo);
            self.updateTableDisplay;
        end

        function deleteSelectedRows(self,h,e,force)

            if nargin < 4
                force = false;
            end
            if ~force
                ButtonName = questdlg('Remove selected rows from database?', ...
                    'Clear Current?', ...
                    'Yes, erase', 'No, go back','No, go back');
            else
                ButtonName = 'Yes, erase';
            end
            switch ButtonName
                case '';
                    return;
                case 'No, go back'
                    return;
                otherwise

                    selected = get(self.handleStruct.table,'UserData');
                    if isempty(selected)
                        warndlg('No rows selected');
                        return;
                    end
                    selectedRowsUnsorted = selected(:,1);
                    selectedRows = self.sortingInds(selectedRowsUnsorted);
                    selectedBool = false(size(self.dataStructure));
                    selectedBool(selectedRows) = true;
                    selectedBoolUnsorted = false(size(self.dataStructure));
                    selectedBoolUnsorted(selectedRowsUnsorted) = true;

                    obsInfo = self.dataStructure(~selectedBool);
                    [~,~,ic] = unique(self.sortingInds(~selectedBoolUnsorted));
                    self.sortingInds = ic;

                    self.updateDataStructure(obsInfo);
                    self.updateTableDisplay;
            end
        end

        function sort(self,inds)
            self.sortingInds = inds;
            self.updateTableDisplay;
        end

        function resizeFunction(self,varargin)

            rowHeightInPixels = self.fontSizeEdit*2; % A decent rule of thumb for edit boxes?
            %border = max(rowHeightInPixels/8,1); % A decent rule of thumb ?
            border = 0;

            parentPosPixels = getpixelposition(self.managedHandle);
            parentHeight = parentPosPixels(4);
            parentWidth = parentPosPixels(3);

            editPos = [1 parentHeight-rowHeightInPixels, parentWidth, rowHeightInPixels];

            % Re-adjust height of the uitable
            tableTopPixels = editPos(2)-border;
            tablePos = [1 1 parentWidth tableTopPixels-1];

            % Set everything
            set(self.handleStruct.table,'units','pixels'); % Reset below
            set(self.handleStruct.table,'position',tablePos);
            set(self.handleStruct.table,'units','normalized');
        end

        function updateTableDisplay(self)
            set(self.handleStruct.table,'Data',self.dataCell(self.sortingInds,:));
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [fcomp] = qsubcompile(fname, varargin)

% QSUBCOMPILE compiles your function into an standalone executable that can easily
% be distributed on a cluster by QSUBCELLFUN. Running a compiled version of your
% function does not take any additional MATLAB licenses. Note that it does require
% that the corresponding MATLAB run-time environment (MCR) is installed on your
% cluster.
%
% Use as
%   compiledfun = qsubcompile(fname)
%   argout      = qsubcellfun(compiledfun, argin, ...)
% or
%   compiledfun = qsubcompile(fname)
%   jobid       = qsubfeval(compiledfun, argin, ...)
%   argout      = qsubget(jobid)
%
% Optional input arguments should be specified in key-value pairs
% and can include
%   batchid     = string that is used for the compiled application filename
%                 and to identify the jobs in the queue, the default is
%                 automatically determined and looks like user_host_pid_batch.
%   toolbox     = string or cell-array with strings, additional Mathworks
%                 toolboxes to include (see below).
%   executable  = string with the name of a previous compiled executable
%                 to start, which usually takes the form "run_xxx.sh". This
%                 implies that compilation does not have to be done.
%   numthreads  = number of threads, can be 1 or inf (default = 1)
%
% When executing a single batch of jobs using QSUBCELLFUN, you can also
% compile your function on the fly with the compile flag like this
%   argout      = qsubcellfun(fname, argin, ..., 'compile', 'yes')
% Using this syntax, the compiled function will be automatically cleaned
% up immediately after execution.
%
% If you need to include additional functions that are not automatically
% detected as dependencies by the MATLAB compiler, e.g. because using
% constructs like feval(sprintf(...)), you can specify fname as a
% cell-array. For example
%   compiledfun = qsubcompile({@ft_definetrial, @trialfun_custom})
%
% If you need to include Mathworks toolboxes that are not automatically
% detected as dependencies by the MATLAB compiler, you can specify them
% like this
%   compiledfun = qsubcompile(fname, 'toolbox', {'signal', 'image', 'stats'})
%
% A common problem for compilation is caused by the use of addpath in
% your startup.m file. Please change your startup.m file into
%   if ~isdeployed
%    % here goes the original content of your startup file
%    % ...
%   end
%
% If you want to execute the same function multiple times with different input
% arguments, you only have to compile it once. The name of the executable can be
% specified as input parameter, and the specified function within the executable
% can be re-execured. An example is specyfying the executable as run_fieldtrip.sh,
% which is a compiled version of the complete FieldTrip toolbox.
%
% See also QSUBCELLFUN, QSUBFEVAL, MCC, ISDEPLOYED

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

% Note that the user-function will be wrapped into qsubexec, which takes care of
% the loading and saving of the input and output arguments.

% It is possible for the user to have multiple MATLAB sessions running at the same
% time (on the same or different computers) and to start multiple instances of
% qsubcellfun. To ensure that those do not interfere with each others, each batch
% of jobs (i.e. instance of qsubcellfun) should get a unique identifier that is
% used in the filename of the temporary mat files.

% get the optional input arguments
batch      = ft_getopt(varargin, 'batch',   getbatch());               % this is a number that is automatically incremented
batchid    = ft_getopt(varargin, 'batchid', generatebatchid(batch));   % this is a string like user_host_pid_batch
toolbox    = ft_getopt(varargin, 'toolbox', {});
executable = ft_getopt(varargin, 'executable');
numthreads = ft_getopt(varargin, 'numthreads', 1); % our experience is that in many cases it is better to have multiple jobs on a single machine instead of a single multi-threaded job

if ischar(toolbox)
  % this should be a cell-array with strings
  toolbox = {toolbox};
end

% some temporary filse are made during compilation, these flags determine
% whether they can be cleaned up afterwards
hasreadme = exist('./readme.txt', 'file');
hasmcclog = exist('./mccExcludedFiles.log', 'file');

if iscell(fname)
  fdeps = fname(2:end);  % this remains a cell-array
  fname = fname{1};      % this is a handle or string
else
  fdeps = {};
end

if isa(fname, 'function_handle')
  % convert the function handle back into a string (e.g. @plus should be 'plus')
  fname = func2str(fname);
end

for i=1:length(fdeps)
  if isa(fdeps{i}, 'function_handle')
    fdeps{i} = func2str(fdeps{i});
  end
end

if numthreads==1
  % pass the singleCompThread option to the MCR
  runtimeopt = '-nodisplay,-singleCompThread';
elseif numthreads>0 && isinf(numthreads)
  % the default behavior of the MCR is to start as many computational threads as cores available
  runtimeopt = '-nodisplay';
else
  error('unsupported number of threads (%d)', numthreads);
end

if ~isempty(toolbox)
  % each toolbox should be added to the mcc command line as -p <full_path>
  toolboxopt = cell(1,length(toolbox));
  for i=1:length(toolbox)
    % find the directory where the toolbox is to be found
    toolboxpath = fileparts(which(fullfile(toolbox{i}, 'Contents.m')));
    if isempty(toolboxpath)
      error('the Mathworks toolbox "%s" could not be found', toolbox{i});
    else
      fprintf('including %s\n', toolboxpath);
    end
    toolboxopt{2*(i-1)+1} = '-p';
    toolboxopt{2*(i-1)+2} = toolboxpath;
  end
else
  toolboxopt = {};
end

if ~isempty(executable)
  % compilation is not needed, just wrap the user-specified executable and the function name
  % into the output structure for use by qsubcellfun or qsubfexec
  [p, f, x] = fileparts(executable);
  if isempty(p) || p(1)~='/'
    p = fullfile(pwd, p);
  end
  % precisely locate the executable, including the full path
  executable = fullfile(p, [f x]);
  if ~exist(executable, 'file')
    error('executable "%s" does not exist', executable);
  end
  fprintf('using "%s"\n', executable);
else
  % try to compile into a stand-allone application
  fprintf('compiling %s into %s\n', fname, batchid);
  % ensure that cellfun is included, it might be needed for stacked jobs
  mcc('-N', '-R', runtimeopt, '-o', batchid, toolboxopt{:}, '-m', 'qsubexec', 'cellfun', fname, fdeps{:});
  fprintf('finished compiling\n');
  executable = fullfile(pwd, sprintf('run_%s.sh', batchid));
end

if ~hasreadme
  delete('./readme.txt');
end

if ~hasmcclog
  delete('./mccExcludedFiles.log');
end

% remember all details
fcomp.fname       = fname;
fcomp.fdeps       = fdeps;
fcomp.batch       = batch;
fcomp.batchid     = batchid;
fcomp.executable  = executable;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%%
%% Run this script in Matlab command window
%%

   function Installmex

   curdir = pwd;
   fprintf(' current directory is:  %s\n',curdir);
%%
%% generate mex files in Mexfun
%%
   if strcmp(computer,'PCWIN64') | strcmp(computer,'GLNXA64')
      computer_model = 64;
   else
      computer_model = 32;
   end
   matlabversion = sscanf(version,'%f');
   matlabversion = matlabversion(1);
   fsp = filesep;
%%
%%
%%
   src = [curdir,fsp,'Solver',fsp,'Mexfun'];
   eval(['cd ','Solver',fsp,'Mexfun']);
   fprintf ('\n Now compiling the mexFunctions in:\n');
   fprintf (' %s\n',src);

   fname{1} = 'mexProd2';
   fname{2} = 'mexProd2nz';
   fname{3} = 'mexinprod';
   fname{4} = 'mexmat';
   fname{5} = 'mexsmat';
   fname{6} = 'mexsvec';
   fname{7} = 'mexschur';
   fname{8} = 'mexqops';
   fname{9} = 'mexexpand';
   fname{10} = 'mexskron';
   fname{11} = 'mexnnz';
   fname{12} = 'mexschurfun';
   fname{13} = 'mexMatvec';
   fname{14} = 'mextriang';
   fname{15} = 'mextriangsp';

   if (matlabversion < 7.3)
      mexcmd = 'mex  -O  ';
   else
      mexcmd = 'mex  -O -largeArrayDims  ';
   end
   for k = 1:length(fname)
      cmd([mexcmd,fname{k},'.c']);
   end
   cd ..
   cd ..
%%
%% generate mex files in spchol
%%
   if (matlabversion < 7.3)
      clear fname
      src = [curdir,fsp,'Linsysolver',fsp,'spchol'];
      eval(['cd ','Linsysolver',fsp,'spchol']);
      fprintf ('\n Now compiling the mexFunctions in:\n');
      fprintf (' %s\n',src);

      fname{1} = 'mexordmmd.c ordmmd.c';
      fname{2} = 'mexsymbfct.c  symbfct.c';
      fname{3} = 'choltmpsiz.c';
      fname{4} = 'cholsplit.c';
      fname{5} = 'mexsparchol.c sparchol2.c sdmauxFill.c sdmauxScalarmul.c';
      fname{6} = 'mexfwblkslv.c sdmauxScalarmul.c';
      fname{7} = 'mexbwblkslv.c sdmauxFill.c sdmauxRdot.c';
      mexcmd = 'mex  -O  ';
      for k = 1:length(fname)
         cmd([mexcmd,fname{k}]);
      end
      cd ..
      cd ..
   end
%%***********************************************
   function cmd(s)

   fprintf(' %s\n',s);
   eval(s);
%%***********************************************

author: yu-jiang

repo: radpbook

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


text:
function [y, C] = noisy_joint_density(s, w, mu, sigma, noise_dev)

% Evaluate the joint density of hidden and observed variables
% Note that the joint density is *not* normalized, i.e., it is only
% evaluated up to proportionality. This is because it can assume Inf values
% if mu << 0 and sigma >> 0 and the proportionalty factor is not needed in
% the EM algorithm.
%
% Syntax:
%   [y, C] = joint_density(s, w, mu, sd, sn)
%
% Input:
%   s  : Value of hidden variable.
%
%   w  : Value of wavelet coefficient. Must have the same size as s
%
%   mu : Mean of hidden variable.
%
%   sd : Standard deviation of hidden variable.
%
%   sn : Standard deviation of noise.
%
%
% Output:
%   y  : Value of density
%
%   C  : Normalizing factor for density.
%
% See also: JOINT_DENSITY

if isscalar(w)
    w2 = w^2 * ones(size(s));
else
    w2 = w.^2;
end

% The 'upper' exponential
logy = -0.5*((s - mu)/sigma).^2;

non_zero_idx = w2 ~= 0;
observed_var = exp(s) + noise_dev^2;
logy(non_zero_idx) = logy(non_zero_idx) ...
    - 0.5*w2(non_zero_idx)./observed_var(non_zero_idx);

y = exp( logy ) ./ sqrt(observed_var);

% The normalization factor
if nargout == 2
	C = 1/(2*pi*sigma);
end

end

author: Sable

repo: mcbench-benchmarks


text:
function Sp = querytp(S, U, domain, p)
%QUERYTP Query TP model at a given parameter point
%	Sp = QUERYTP(S, U, domain, p)
%
%	S        - core tensor of the TP model
%	U        - (discrete) weight function data of the TP model
%	domain   - intervals for each dimension
%	p        - parameter point where we query the model
%
%	Sp       - system at p (linear interpollated weights)

% TODO: method param to interp1 (spline,..)

W = queryw1(U, domain, p);
Sp = squeeze(tprod(S, W));


text:
% Is In Range function written by NKN(C)-2006
% is c E [a,b]
% True=1;
% False=0;
%
% Example: isinrange(1,2,3)
% ans=0;
function g=isinrange(a,b,c)
bb=max(a,b);
aa=min(a,b);
if (c<=bb) & (c>=aa)
%% disp('true');
g=1;
else
%% disp('false');
g=0;
end


text:
function [C] = K2C(K)
% Convert temperature from Kelvin to degrees Celsius.
C = K-273.15;


text:
function [Sf,AI,AIrows,AIcols] = imageStats1(f)
%imageStats1 Sample function used in Chapter 2.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt
Sf = size(f);
AI = mean2(f);
AIrows = mean(f,2);
AIcols = mean(f,1);

author: dipum

repo: dipum-toolbox

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


text:
function valid_dims = computeValidDimensions(data,options)
% Evaluate a matrix of connections S and return which elements are to be
% analysed for variance and standardisation.
S = options.S;
if ~isfield(options,'distribution') || strcmp(options.distribution,'Gaussian')
    if all(S(:)==1)
        valid_dims = [1:length(S)]; % note this omits first dimension to allow TUDA style use of intercepts
    else
        %implies a TUDA/TUCA style setup:
        lastdatadim = find(diff(any(S==1)),1);
        % check if intercept term present:
        if isstruct(data)
             X = data.X;
        else
             X = data;
        end
        if var(X(:,lastdatadim))==0
            lastdatadim = lastdatadim-1;
        end
        valid_dims = [1:lastdatadim];
    end
else
    valid_dims = [];
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function [nodes,edges,classData] = mrgGrowGray(classData,numLayers,layer0,hemisphere)
%
%     [nodes,edges,classData] = mrgGrowGray(classInfo,[numLayers],[layer0],[hemisphere='left'])
%
% This routine calls the mex function grow_gray with the parameters
% specified by the user and loads the class file if necessary
%
% INPUTS :
%     - classInfo : a variable loaded from a class file or
%                   a path string pointing to the class file
%     - numLayers : number of layers to grow the gray matter
%     - layer0 : optional input argument. (Default 0)
%                   - 0 : no white matter included in the gray graph
%                   - 1 : the white boundary included in the gray graph.
%                         In this case, the layer number of the boundary will be 0.
%                   - 2 : All the white matter fully connected is included
%                         in the gray graph
%
% OUTPUTS :
%     - nodes : structure containing nodes of the gray matter
%     - edges : structure containing edges of the gray matter
%     - classData : matlab class file structure containing the output
%
% Example:
%  To visualize a mesh with, say, two layers of gray superimposed do this.
%   fName=fullfile(mrvDataRootPath,'anatomy','anatomyV','left','left.Class');
%   [nodes,edges,classData] = mrgGrowGray(fName,2);
%   wm = uint8( (classData.data == classData.type.white) | (classData.data == classData.type.gray));
%   msh = meshColor(meshSmooth(meshBuildFromClass(wm,[1 1 1])));
%   meshVisualize(msh,2);
%
%  To see the gray matter and white matter in a slice, along with their
%  connectivity, use this.  Change the third argument of mrGrowGray to
%  change the type of connectivity you visualize.
%
%   fName ='X:\anatomy\nakadomari\left\20050901_fixV1\left.Class';
%   [nodes,edges,classData] = mrgGrowGray(fName,2,2);
%   mrgDisplayGrayMatter(nodes,edges,80,[120 140 120 140]);
%
% NOTE:  M. Schira says that the gray nodes should be more densely
% connected to reduce the error in the distance measurements.  We should
% look into completing the connections in mrGrowGray.
%
% HISTORY:
%  GB 01/11/06 wrote it.
%  2008.07.07 RFD: fixed bug in VOI clipping.
%  2008.12.19 DY: fixed hemisphere variable setting bug
%
% (c) Stanford VISTA Team, 2006

if notDefined('numLayers'),  numLayers = 3; end
if notDefined('layer0'), layer0 = 0; end
if notDefined('hemisphere'), hemisphere = 'left'; end
if ischar(classData), classData = readClassFile(classData,0,0,hemisphere); end

voi  = classData.header.voi;
data = classData.data;
% RFD: added the "| ... 160". I think this is the type code for 'selected
% gray matter', but it's not listed in
data(data(:) == classData.type.gray) = classData.type.unknown;

fprintf('Growing %i gray layers...\n',numLayers);
if layer0,
    fprintf('White matter included in the gray graph...\n');
else
    fprintf('White matter excluded in the gray graph...\n');
end

% Growgray should ignore extra labels, but doesn't-= so we 'clean' it's
% class data.
cleanData = data;
cleanData(cleanData>=classData.type.other) = classData.type.unknown;
[nodes,edges] = grow_gray(cleanData,numLayers,voi,layer0);

grayMatter = nodes(1:3,nodes(6,:) ~= 0)+1;
outOfVoi = grayMatter(1,:)<=voi(1) | grayMatter(1,:)>voi(2) ...
         | grayMatter(2,:)<=voi(3) | grayMatter(2,:)>voi(4) ...
         | grayMatter(3,:)<=voi(5) | grayMatter(3,:)>voi(6);
grayMatter(:,outOfVoi) = [];
data(sub2ind(size(data),...
    grayMatter(1,:) - voi(1),...
    grayMatter(2,:) - voi(3),...
    grayMatter(3,:) - voi(5))) = classData.type.gray;

classData.data = data;

return;


text:
function sp = esprit_1d(R, n, k, varargin)
%ESPRIT_1D 1D ESPRIT for ULAs.
%Syntax:
%   sp = ESPRIT_1D(R, n, design, wavelength, grid_size, ...);
%Inputs:
%   R - Sample covariance matrix.
%   n - Number of sources.
%   k - 2*pi*inter_element_spacing/wavelength.
%   ... - Options:
%           'Unit' - Can be 'radian', 'degree', or 'sin'. Default value is
%                   'radian'.
%           'Displacement' - The displacement between the two overlapping
%                            subarrays measured in number of element
%                            spacings. Default value is 1.
%                            Note: increasing this value will lead to
%                            smaller unambiguous range. Make sure your DOAs
%                            falls within the unambiguous range.
%           'Formulation' - Either 'TLS' (Total Lease Squares) or 'LS'
%                           (Least Squares). Default value is 'TLS'.
%           'RowWeights' - Specifies the row weights with a vector or
%                          string. Default value is 'Default', which
%                          generates the following weight vector:
%                           [1 sqrt(2) sqrt(3) ... sqrt(3) sqrt(2) 1]
%                          You can disable row weighting by passing in
%                          'Identity' or 'Off'.
%Output:
%   sp - Spectrum structure with the following fields:
%           x - An 1 x grid_size vector.
%           y - An 1 x grid_size vector. Calling `plot(x, y)` will plot the
%               spectrum.
%           x_est - An 1 x n vector storing the estimated DOAs.
%           x_unit - The same as the unit specified by 'Unit'.
%           resolved - Constant value true.
%           discrete - Constant value true.
%Reference:
%   [1] H. L. Van Trees, Optimum array processing. New York: Wiley, 2002.
use_tls = true;
ds = 1;
row_weights = [];
use_row_weights = false;
unit = 'radian';
for ii = 1:2:nargin-3
    option_name = varargin{ii};
    option_value = varargin{ii+1};
    switch lower(option_name)
        case 'unit'
            unit = option_value;
        case 'formulation'
            switch lower(option_value)
                case 'ls'
                    use_tls = false;
                case 'tls'
                    use_tls = true;
                otherwise
                    error('Formulation must be either ''LS'' or ''TLS''.');
            end
        case 'displacement'
            if option_value < 1 || mod(option_value, 1) ~= 0
                error('Displacement must be an integer that is greater or equal to one.');
            end
            ds = option_value;
        case 'rowweights'
            if ischar(option_value)
                switch lower(option_value)
                    case 'default'
                        use_row_weights = true;
                    case {'off', 'identity'}
                        use_row_weights = false;
                    otherwise
                        error('Either specify the row weights manually or pass in ''Default'' to use the default weights, or ''Identity'', ''Off'' to disable row weighting.');
                end
            else
                use_row_weights = true;
                row_weights = option_value(:);
            end
        otherwise
            error('Unknow option "%s".', option_name);
    end
end
m = size(R, 1);
if n > m - ds
    error('Too many sources.');
end
if ~isempty(row_weights) && length(row_weights) ~= m - ds
    error('The dimension of the row weights vector is not equal to (m - displacement).');
end
% ESPRIT
[E, ~] = eig(0.5*(R + R'), 'vector');
Es = E(:,end - n + 1:end);
% apply weights if necessary
if use_row_weights
    if isempty(row_weights)
        % default weights
        if mod(m - ds, 2) == 1
            w_max = (m - ds - 1)/2;
            row_weights = sqrt([1:w_max w_max + 1 w_max:-1:1])';
        else
            w_max = (m - ds)/2;
            row_weights = sqrt([1:w_max w_max:-1:1])';
        end
    end
    Es1 = bsxfun(@times, row_weights, Es(1:end - ds,:));
    Es2 = bsxfun(@times, row_weights, Es(ds + 1:end,:));
else
    Es1 = Es(1:end - ds,:);
    Es2 = Es(ds + 1:end,:);
end
if use_tls
    % TLS estimate
    C = [Es1 Es2];
    C = C'*C;
    C = 0.5*(C + C');
    [V, l] = eig(C, 'vector');
    [~, idx] = sort(real(l), 'descend');
    V = V(:,idx);
    V12 = V(1:n,n + 1:end);
    V22 = V(n + 1:end,n + 1:end);
    Phi = -V12/V22;
else
    % LS estimate
    Phi = (Es1'*Es1)\(Es1'*Es2);
end
% convert z to spectrum
z = eig(Phi);
sp = struct();
sp.x_est = sort(cm2doa(z, k*ds, unit));
sp.x = sp.x_est;
sp.x_unit = unit;
sp.y = ones(1, n);
sp.resolved = true;
sp.discrete = true;
end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools