function [lgraph,hyperParams,numsNetParams,FLOPs] = importDarknetWeights(cfgfile,weightsfile,cutoffModule)
% IMPORTDARKNETWEIGHTS 功能：指定导入部分module的darknet模型
% 输入：cfgfile, 指定的cfg后缀的模型描述文件
%       weighfile, 指定的.weights后缀的二进制文件
%       cutoffModule,(可选项)1*1的正整数，(可选项)1*1的正整数，指定导出darknet前cutoffModule个module。以cfg文件中第一个非[net]开始的module为0开始的计数，没有该项则导出整个网络
% 输出：lgraph， matlab深度学习模型计算图
%      hyperParams,结构体，超参配置文件
%      numsReadParams,权重参数个数
%      FLOPs， 模型计算力
% 注意：1、适合2019b版本及以上
%       2、leaky阈值darknet为0.1
%       3、如果某个module中有bn层，则conv的bias为0，因为darknet是这种存储形式
%      4、darknet weights保存顺序依次为BN层offset,scale,mean,variance,Conv层的bias,weights
%      特征图输出output Size = (Input Size – ((Filter Size – 1)*Dilation Factor + 1) + 2*Padding)/Stride + 1
% 参考：1、官方文档，Specify Layers of Convolutional Neural Network
%      2、https://www.zhihu.com/question/65305385
%       3、https://github.com/ultralytics/yolov3/blob/master/models.py
% cuixingxing150@gmail.com
% 2019.8.19
% 2019.9.4修改，由原来的darknet中[net]为0开始的索引改为以cfg文件中第一个非[net]开始的module为0开始的计数的索引
% 2020.4.25修改函数默认输入参数
%
arguments
    cfgfile (1,:) char
    weightsfile (1,:) char
    cutoffModule {mustBeNonnegative} = 0 % 默认导入所有的层
end

[lgraph,hyperParams,numsNetParams,FLOPs,...
    moduleTypeList,moduleInfoList,layerToModuleIndex] = importDarkNetLayers(cfgfile);% only weights需要cutoff,layer都导入
assert(length(moduleTypeList)==length(moduleInfoList));

%% 读取权重参数文件
fid_w = fopen(weightsfile,'rb');
if fid_w == -1
  error('Author:Function:OpenFile', 'Cannot open file: %s', weightsfile);
end
header = fread(fid_w, 3, '*int32');
if header(2) > 1
    header2 = fread(fid_w, 1, '*int64'); % int64占用8个字节
else
    header2 = fread(fid_w, 1, '*int32'); % int32占用4个字节
end
fprintf('Major :%d, Minor :%d,Revision :%d,number of images during training:%d,reading params...\n',...
    header(1),header(2),header(3),header2);
weights = fread(fid_w,'*single');
fclose(fid_w);

% numsWeightsParams = numel(weights);
readSize = 1;
numsModule = length(moduleTypeList);
if cutoffModule>0
    numsModule = cutoffModule+1;% [net] plus 1
end

for i = 1:numsModule
    currentModuleType = moduleTypeList{i};
    currentModuleInfo = moduleInfoList{i};
    if strcmp(currentModuleType,'[convolutional]')
        currentModule = lgraph.Layers(i==layerToModuleIndex);
        filterSize = str2double(currentModuleInfo.size);
        numFilters = str2double(currentModuleInfo.filters);
        channels_in = moduleInfoList{i-1}.channels;

        if isfield(currentModuleInfo,'batch_normalize')
            % bn bias
            bn_bias = weights(readSize:readSize+numFilters-1);
            bn_bias = reshape(bn_bias,[1,1,numFilters]);
            currentModule(2).Offset = bn_bias;
            readSize = readSize+numFilters;
            % bn weight
            bn_weight = weights(readSize:readSize+numFilters-1);
            bn_weight = reshape(bn_weight,[1,1,numFilters]);
            currentModule(2).Scale = bn_weight;
            readSize = readSize+numFilters;
             % bn trainedMean
            bn_mean = weights(readSize:readSize+numFilters-1);
            bn_mean = reshape(bn_mean,[1,1,numFilters]);
            currentModule(2).TrainedMean = bn_mean;
            readSize = readSize+numFilters;
             % bn trainedVariance
            bn_var = weights(readSize:readSize+numFilters-1);
            bn_var = reshape(bn_var,[1,1,numFilters]);
            if any(bn_var<-0.01)
                error("方差应该大于0！");
            end
            currentModule(2).TrainedVariance = abs(bn_var); % 防止接近于0的数是负数
            readSize = readSize+numFilters;
            % conv bias 为0
            if isfield(currentModuleInfo,'groups')
                numGroups = str2double(currentModuleInfo.groups);
                numFiltersPerGroup_out = numFilters/numGroups;
                currentModule(1).Bias = zeros(1,1,numFiltersPerGroup_out,numGroups,'single');
            else
                currentModule(1).Bias = zeros(1,1,numFilters,'single');
            end
        else
            % load conv bias
            conv_bias = weights(readSize:readSize+numFilters-1);
            if isfield(currentModuleInfo,'groups')
                numGroups = str2double(currentModuleInfo.groups);
                numFiltersPerGroup_out = numFilters/numGroups;
                conv_bias = reshape(conv_bias,1,1,numFiltersPerGroup_out,numGroups);
            else
                conv_bias = reshape(conv_bias,1,1,numFilters);
            end
            currentModule(1).Bias = conv_bias;
            readSize = readSize+numFilters;
        end % end of is bn
        % load conv weights
        if isfield(currentModuleInfo,'groups')
            numGroups = str2double(currentModuleInfo.groups);
            numFiltersPerGroup_out = numFilters/numGroups;
            nums_conv_w = filterSize*filterSize*channels_in/numGroups*numFiltersPerGroup_out*numGroups;
            conv_weights = weights(readSize:readSize+nums_conv_w-1);
            conv_weights = reshape(conv_weights,filterSize,filterSize,channels_in/numGroups,numFiltersPerGroup_out,numGroups);
            conv_weights = permute(conv_weights,[2,1,3,4,5]);
            currentModule(1).Weights = conv_weights;
            readSize = readSize+nums_conv_w;
        else
            nums_conv_w = filterSize*filterSize*channels_in*numFilters;% weights
            conv_weights = weights(readSize:readSize+nums_conv_w-1);
            conv_weights = reshape(conv_weights,filterSize,filterSize,channels_in,numFilters);
            conv_weights = permute(conv_weights,[2,1,3,4]);
            currentModule(1).Weights = conv_weights;
            readSize = readSize+nums_conv_w;
        end % end of load conv weights
        % 更新参数
        % lgraph.Layers(i==layerToModuleIndex) = currentModule;
        for replaceInd = 1:length(currentModule)
            layerName = currentModule(replaceInd).Name;
            lgraph = replaceLayer(lgraph,layerName,currentModule(replaceInd));
        end
    end % end of module '[convolutional]'

    % fullyConnectedLayer
    if  strcmp(currentModuleType,'[connected]')
        currentModule = lgraph.Layers(i==layerToModuleIndex);
        numFilters = str2double(currentModuleInfo.output);
        % load fc bias
        numBias = numFilters;
        fl_bias = weights(readSize:readSize+numBias-1);
        fl_bias = reshape(fl_bias,numBias,1);
        currentModule(1).Bias = fl_bias;
        readSize = readSize+numBias;
         % load fc weights
        input_all_neurons = prod(moduleInfoList{i-1}.mapSize)*moduleInfoList{i-1}.channels;
        numWeights = numFilters*input_all_neurons; % help fullyConnectedLayer weights
        fl_weights = weights(readSize:readSize+numWeights-1);
        fl_weights = reshape(fl_weights,input_all_neurons,numFilters);
        fl_weights = permute(fl_weights,[2,1]);% fc不需要permute?
        currentModule(1).Weights = fl_weights;
        readSize = readSize+numWeights;
        % 更新参数
        for replaceInd = 1:length(currentModule)
            layerName = currentModule(replaceInd).Name;
            lgraph = replaceLayer(lgraph,layerName,currentModule(replaceInd));
        end
    end % end of module '[connected]'
end % end of nums of module

if isa(lgraph.Layers(end),'nnet.cnn.layer.SoftmaxLayer')
    lastLayerName = lgraph.Layers(end).Name;
    classifyLayer = classificationLayer('Name','classify');
    lgraph = addLayers(lgraph,classifyLayer);
    lgraph = connectLayers(lgraph,lastLayerName,'classify');
end

fprintf('Load parameters succfully!\n')


text:
function [sigma2,Sigma,didOverflow]=spherHarmonicCov(CStdDev,SStdDev,point,a,c,scalFactor)
%%SPHERHARMONICCOV Evaluate the variance of a potential or the covariance
%                  matrix of a gradient that one might compute using a
%                  spherical harmonic coefficient model with
%                  spherHarmonicEval. Here, one provides the standard
%                  deviations associated with the coefficients in the
%                  model, such as those obtained using the function
%                  getEGMGravCoeffs for the EGM2008 gravity model.
%
%INPUTS: CStdDev A length (M+2)*(M+1)/2 array holding the standard
%               deviations of the coefficient terms that are multiplied by
%               cosines in the harmonic expansion. The coefficients must be
%               fully normalized using the type of full normalization that
%               is used in the EGM2008 model. Their normalization type can
%               be changed using the changeSpherHarmonicNorm function. If
%               given to a CountingClusterSet class, C(n+1,m+1) is the
%               coefficient of degree n and order m. When a maximum degree
%               of M is used, all C must have values for all n from 0 to M
%               and for all m from 0 to n for each n. If coefficients are
%               not present for certain degrees, then insert a 0. It is
%               assumed that M>=3. For magnetic field models, the units of
%               C should generally be Tesla. For gravitational models, C is
%               generally unitless with the units being determined entirely
%               by the parameters a and c.
%       SStdDev A length (M+2)*(M+1)/2 array holding the standard
%               deviations of the coefficient terms that are multiplied by
%               sines in the harmonic expansion. The requirements on
%               SStdDev are the same as those on CStdDev.
%         point The 3XN set of N points at which the potential and/or
%               gradient should be evaluated given in SPHERICAL, ECEF
%               coordinates consisting of [r;azimuth;elevation]; When
%               evaluating points on a grid, the algorithm will be fastest
%               if the points are provided presorted by range and then by
%               azimuth. This reduces the amount of recomputation of
%               certain values. Alternatively, if C and S are for
%               evaluating terrain heights, then points are 2XN having the
%               format [azimuth;elevation] and it is best if the points are
%               sorted by azimuth.
%             a The numerator in the (a/r)^n term in the spherical harmonic
%               sum. Normally, this is some type of a reference radius. For
%               example, when using most gravitational models, a is the
%               semi-major axis of the reference ellipsoid. If this
%               parameter is omitted, it is assumed that one is using the
%               spherical harmonics with something like the National
%               Geospatial Intelligence Agency's (NGA's) EGM96 or EGM2008
%               models, in which case a=Constants.EGM2008SemiMajorAxis is
%               used unless point is 2D, in which case c=1 is used.
%             c The constant value by which the spherical harmonic series
%               is multiplied. For example, for gravitational potentials,
%               the value is usually GM where G is the universal
%               gravitational constant and M is the mass of the Earth. When
%               using the International Geomagnetic Reference Field (IGRF),
%               the constant is a^2, where a is the same as the numerator
%               in the a/r term. If this parameter is omitted, it is
%               assumed that one is using the spherical harmonics with
%               something like the NGA's EGM96 or EGM2008 models, in which
%               case c=Constants.EGM2008GM is used unless point is 2D, in
%               which case c=1 is used.
%    scalFactor An optional scale factor used in computing the normalized
%               associated Legendre polynomials if fullyNormalized=true.
%               Generally, the default value (is if the scalFactor
%               parameter is omitted) of 2^(-470) is sufficient. When very
%               high-order models are used, this scale factor prevents
%               overflows. However overflows (and a loss of precision) are
%               unavoidable when using the full EGM2008 model. These
%               effects are worse near the poles.
%
%OUTPUTS: sigma2 The NX1 vector of variances (squared standard deviations)
%                of the potential estimate at the given points.
%          Sigma The covariance matrix of the gradient of the potential at
%                the given points.
%    didOverFlow This indicates whether during the computation of sigma2 or
%                Sigma there were any overflow errors leading to term
%                being dropped. Note that this does not indicate potential
%                losses of precision due to underflow errors making terms
%                zero, which becomes more common the smaller scalFactor is.
%                Also, if scaling is extremely bad, it is possible for
%                NaNs or Inf terms to still be returned in sigma2 and
%                Sigma.
%
%The algorithm used here is described in [1].
%
%Since Matlab uses double precision arithmetic, when using high degree and
%order models, such as the full 2190 degree EGM2008 model, the precision of
%Sigma will be reduced as higher-order terms can experience overflow
%problems and are thus discarded to avoid NaNs from occurring. Lower-order
%models will not suffer from the same problem.
%
%EXAMPLE:
%Here, we evaluate the function at a point on the reference ellipsoid and
%then at a point at the pole using the full degree 2190 EGM2008 model. Note
%that the Matlab implementation will bee too slow for this high a degree
%and the C++ implementation should be compiled.
% latLongAlt=[-30*(pi/180);45*(pi/180);0];
% spherLoc=Cart2Sphere(ellips2Cart(latLongAlt));
%
% [C,S,a,c,CStdDev,SStdDev]=getEGMGravCoeffs(2190,false);
% [sigma21,Sigma1,didOverflow1]=spherHarmonicCov(CStdDev,SStdDev,spherLoc,a,c)
% spherLoc(end)=pi/2;%90 degree elevation --the pole.
% [sigma22,Sigma2,didOverflow2]=spherHarmonicCov(CStdDev,SStdDev,spherLoc,a,c)
%One will see that in both instances, finite values are returned, but for
%the point at the pole, overflows are flagged indicating the loss of
%prevision due to some terms being dropped due to overflows. A model with
%fewer spherical harmonic coefficients would not be as susceptible to such
%overflow problems.
%
%REFERENCES:
%[1] D. F. Crouse, "An Overview of Major Terrestrial, Celestial, and
%    Temporal Coordinate Systems for Target Tracking," Formal Report, Naval
%    Research Laboratory, no. NRL/FR/5344--16-10,279, 10 Aug. 2016, 173
%    pages.
%
%April 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<6||isempty(scalFactor))
    scalFactor=2^(-470);
end

if(nargin<5||isempty(c))
    if(size(point,1)==2)
        c=1;
    else
        c=Constants.EGM2008GM;
    end
end

if(nargin<4||isempty(a))
    if(size(point,1)==2)
        a=1;
    else
        a=Constants.EGM2008SemiMajorAxis;
    end
end

M=(1/2)*(sqrt(1+8*length(CStdDev))-1)-1;

if(M<3)
    error('The coefficients must be provided to at least degree 3. To use a lower degree, one can insert zero coefficients.');
end

numPoints=size(point,2);
%If we are evaluating terrain heights.
switch(size(point,1))
    case 2
        point=[ones(1,numPoints);point(1,:);point(2,:)];
    case 3
    otherwise
        error('Invalid point length');
end

CStdDev=CountingClusterSet(CStdDev);
SStdDev=CountingClusterSet(SStdDev);

sigma2=zeros(numPoints,1);
Sigma=zeros(3,3,numPoints);

didOverflow=false;
rPrev=Inf;
thetaPrev=Inf;
for curPoint=1:numPoints
    r=point(1,curPoint);
    thetaCur=point(3,curPoint);

    rChanged=r~=rPrev;
    rPrev=r;

    if(rChanged)
        nCoeff=zeros(M+1,1);
        nCoeff(1)=1;
        for n=1:M
            nCoeff(n+1)=nCoeff(n)*(a/r);
        end
    end

    %The non-singular algorithm of Pines using the fully normalized
    %Helmholtz equations from Fantino and Casotto is used. The algorithm
    %has been slightly modified so that the c/r term is out front and the
    %fully normalized Helmholtz polynomials can be scaled. Also, lumped
    %coefficients are not used. The Pines algorithm can suffer a loss of
    %precision near the equator. However, it is simple to just square the
    %terms in the sum.
    CartPoint=spher2Cart(point(:,curPoint));

    x=CartPoint(1);
    y=CartPoint(2);
    z=CartPoint(3);

    %Get the direction cosines used by Pines' algorithm.
    s=x/r;
    t=y/r;
    u=z/r;

    %Compute the fully normalized Helmholtz polynomials.
    if(thetaPrev~=thetaCur)
        [HBar,dHBardu]=normHelmholtz(u,M,scalFactor);
        thetaPrev=thetaCur;
    end

    %Recursively compute the rm and im terms for the sums.
    rm=zeros(M+1,1);
    im=zeros(M+1,1);
    rm(0+1)=1;
    im(0+1)=0;
    for m=1:M
        %These are equation 49 in the Fantino and Casotto paper.
        rm(m+1)=s*rm(m-1+1)-t*im(m-1+1);
        im(m+1)=s*im(m-1+1)+t*rm(m-1+1);
    end

    %Perform the sum for the potential from Equation 44 in the Fantino and
    %Casotto paper, but square all of the terms to represent a sigma.
    %All cross terms are expected to be zero, so the sum is very similar to
    %the sum for the potential in spherHarmonicEval.
    sigma2(curPoint)=0;
    for n=0:M
        innerTerm=0;
        for m=0:n
            innerTerm=innerTerm+(CStdDev(n+1,m+1)*rm(m+1)*HBar(n+1,m+1))^2+(SStdDev(n+1,m+1)*im(m+1)*HBar(n+1,m+1))^2;
        end
        if(isfinite(innerTerm))
            sigma2(curPoint)=sigma2(curPoint)+nCoeff(n+1)^2*innerTerm;
        else
            didOverflow=true;
        end
    end

    sigma2(curPoint)=(c/r)^2*sigma2(curPoint)/scalFactor^2;

    %Now, compute the cosigma matrix of the gradient, if requested.
    if(nargout>1)
        a11=0;
        a22=0;
        a33=0;
        a44=0;
        a12=0;
        a13=0;
        a14=0;
        a23=0;
        a24=0;
        a34=0;

        %The equations in these loops are from Table 10.
        for n=0:M
            a11Loop=0;
            a22Loop=0;
            a12Loop=0;
            a13Loop=0;
            a14Loop=0;
            a23Loop=0;
            a24Loop=0;

            %The m=0 case only applies to a3 and a4, so that means only to
            %a33, a34, and a44.
            m=0;
            HVal=HBar(n+1,m+1);
            dHVal=dHBardu(n+1,m+1);
            CProdMN=CStdDev(n+1,m+1)^2*rm(m+1)^2+SStdDev(n+1,m+1)^2*im(m+1)^2;
            Lmn=(n+m+1)*HVal+u*dHVal;%Defined in Table 14.

            a33Loop=CProdMN*dHVal^2;
            a44Loop=CProdMN*Lmn^2;
            a34Loop=-CProdMN*Lmn*dHVal;

            for m=1:n
                HVal=HBar(n+1,m+1);
                dHVal=dHBardu(n+1,m+1);

                CProdMN=CStdDev(n+1,m+1)^2*rm(m+1)^2+SStdDev(n+1,m+1)^2*im(m+1)^2;
                Lmn=(n+m+1)*HVal+u*dHVal;

                %These if-statements are to deal with numerical precision
                %problems near the poles. We want to avoid 0*Inf terms due
                %to limitations in the valid range of double precision
                %numbers. Of course, the loss of the terms where overflow
                %occurs means that the covariance matrix will be
                %underestimated.
                if(isfinite(HVal))
                    a11Loop=a11Loop+m^2*(CStdDev(n+1,m+1)^2*rm(m-1+1)^2+SStdDev(n+1,m+1)^2*im(m-1+1)^2)*HVal^2;
                    a12Loop=a12Loop+m^2*rm(m-1+1)*im(m-1+1)*(SStdDev(n+1,m+1)^2-CStdDev(n+1,m+1)^2)*HVal^2;
                    a22Loop=a22Loop+m^2*(SStdDev(n+1,m+1)^2*rm(m-1+1)^2+CStdDev(n+1,m+1)^2*im(m-1+1)^2)*HVal^2;
                else
                    didOverflow=true;
                end

                if(isfinite(Lmn))
                    a44Loop=a44Loop+CProdMN*Lmn^2;
                    if(isfinite(HVal))
                        a14Loop=a14Loop-m*(CStdDev(n+1,m+1)^2*rm(m-1+1)*rm(m+1)+SStdDev(n+1,m+1)^2*im(m-1+1)*im(m+1))*HVal*Lmn;
                        a24Loop=a24Loop-m*(-CStdDev(n+1,m+1)^2*im(m-1+1)*rm(m+1)+SStdDev(n+1,m+1)^2*rm(m-1+1)*im(m+1))*HVal*Lmn;
                    end
                    if(isfinite(dHVal))
                        a34Loop=a34Loop-CProdMN*Lmn*dHVal;
                    end
                else
                    didOverflow=true;
                end

                if(isfinite(dHVal))
                    a33Loop=a33Loop+CProdMN*dHVal^2;
                    if(isfinite(HVal))
                        a13Loop=a13Loop+m*(CStdDev(n+1,m+1)^2*rm(m-1+1)*rm(m+1)+SStdDev(n+1,m+1)^2*im(m-1+1)*im(m+1))*HVal*dHVal;
                        a23Loop=a23Loop+m*(-CStdDev(n+1,m+1)^2*im(m-1+1)*rm(m+1)+SStdDev(n+1,m+1)^2*rm(m-1+1)*im(m+1))*HVal*dHVal;
                    end
                else
                    didOverflow=true;
                end
            end

            a11=a11+nCoeff(n+1)^2*a11Loop;
            a22=a22+nCoeff(n+1)^2*a22Loop;
            a33=a33+nCoeff(n+1)^2*a33Loop;
            a44=a44+nCoeff(n+1)^2*a44Loop;
            a12=a12+nCoeff(n+1)^2*a12Loop;
            a13=a13+nCoeff(n+1)^2*a13Loop;
            a14=a14+nCoeff(n+1)^2*a14Loop;
            a23=a23+nCoeff(n+1)^2*a23Loop;
            a24=a24+nCoeff(n+1)^2*a24Loop;
            a34=a34+nCoeff(n+1)^2*a34Loop;
        end

%These are based on squaring the terms in equation 70, removing cross
%terms. However, an additional 1/r (squared) term has been added,
%which the original paper omitted when going from Equation 68 to 70.

        s11=a11+2*s*a14+s^2*a44;
        s12=a12+s*a24+t*a14+s*t*a44;
        s13=a13+s*a34+u*a14+s*u*a44;
        s22=a22+2*t*a24+t^2*a44;
        s23=a23+t*a34+u*a24+t*u*a44;
        s33=a33+2*u*a34+u^2*a44;

        temp=c/(r^2*scalFactor);

        Sigma(:,:,curPoint)=temp*(temp*[s11,s12,s13;
                                        s12,s22,s23;
                                        s13,s23,s33]);
    end

end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function triangle_symq_rule_test02 ( degree, numnodes, vert1, vert2, vert3 )

%*****************************************************************************80
%
%% TRIANGLE_SYMQ_RULE_TEST02 calls TRIASYMQ for a quadrature rule of given order and region.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    27 June 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer DEGREE, the desired total polynomial degree exactness
%    of the quadrature rule.  0 <= DEGREE <= 50.
%
%    Input, integer NUMNODES, the number of nodes to be used by the rule.
%
%    Input, real VERT1(2), VERT2(2), VERT3(2), the
%    vertices of the triangle.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGLE_SYMQ_RULE_TEST02\n' );
  fprintf ( 1, '  Symmetric quadrature rule for a triangle.\n' );
  fprintf ( 1, '  Polynomial exactness degree DEGREE = %d\n', degree );

  area = triangle_area ( vert1, vert2, vert3 );
%
%  Retrieve and print a symmetric quadrature rule.
%
  [ rnodes, weights ] = triasymq ( degree, vert1, vert2, vert3, numnodes );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  NUMNODES = %d', numnodes );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     J      W               X               Y\n' );
  fprintf ( 1, '\n' );
  for j = 1 : numnodes
    fprintf ( 1, '  %4d  %14.6g  %14.6g  %14.6g\n', ...
      j, weights(j), rnodes(1,j), rnodes(2,j) );
  end

  d = sum ( weights(1:numnodes) );

  fprintf ( 1, '   Sum  %14.6g\n', d );
  fprintf ( 1, '  Area  %14.6g\n', area );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function p=rmReconParams(vw)

% this parameters is needed to load all raw-data
p.wData = 'all';

% convert data to percent BOLD
p.analysis.calcPC  = true;

% number of trends
p.analysis.nDCT = 4;

% stimulus parameters
ns = viewGet(vw,'nscans');
for n=1:ns
    p.stim(n).nFrames    = viewGet(vw,'numframes',n);
    p.stim(n).nUniqueRep = 1;
    p.stim(n).nDCT = p.analysis.nDCT;
end

% thresholds above which to use voxels
p.thresh.rmVarexp = 0.15; % goodness of fit of pRF model
p.thresh.stimTmap = 1.96; % goodness of fit of stimuli

% Number of events in a par file
p.numOfEvents = 37;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function drawAzBeam(azCenter,beamHalfwidth,range,systemType,lRx,M,closeEnds,varargin)
%%DRAWANGLEBEAM Draw the outline in 2D Cartesian coordinates of a polar
%           beam in a particular direction. This does not shape the beam
%           based upon gain. It just outlines a region, drawing a curve on
%           the end at a constant range, if desired. This can be useful for
%           determining whether plotted targets are near various beams in
%           2D.
%
%INPUTS: azCenter The scalar center angle of the beam in radians.
%   beamHalfwidth The distance from the center of the beam to one side of
%                 the beam in radians. This must be >0.
%           range This can be a scalar or a length-2 vector. If a scalar,
%                 it is the maximum extent of the range to which the beam
%                 is drawn. If this is a length-2 vector, then range(1) is
%                 the minimum range to draw and range(2) is the maximum
%                 range to be covered by the beam. This is a one-way range.
%      systemType An optional parameter specifying the axis from which the
%                 angles are measured. Possible values are:
%                 0 (The default if omitted or an empty matrix is passed)
%                   The azimuth angle is counterclockwise from the x axis.
%                 1 The azimuth angle is measured clockwise from the y axis.
%             lRx The 2X1 [x;y] location vector of the receiver in
%                 Cartesian coordinates.  If this parameter is omitted or
%                 an empty matrix is passed, then the receiver is assumed
%                 to be at the origin.
%               M A 2X2 rotation matrices to go from the alignment of the
%                 global coordinate system to that at the receiver. If
%                 omitted or an empty matrix is passed, then it is assumed
%                 that the local coordinate system is aligned with the
%                 global and M=eye(2) --the identity matrix is used.
%
%OUTPUTS: None.
%
%EXAMPLE:
%This function draws 60-degree beams in three directions making an image
%reminiscent of a radioactive symbol.
% figure(1)
% clf
% hold on
% %60 degree beam.
% beamHalfwidth=30*(pi/180);
% range=[2;10];
% systemType=0;
% azCenter=30*(pi/180);
% drawAzBeam(azCenter,beamHalfwidth,range,[],[],[],[],'-k','linewidth',2)
% azCenter=pi-30*(pi/180);
% drawAzBeam(azCenter,beamHalfwidth,range,[],[],[],[],'-k','linewidth',2)
% azCenter=-pi/2;
% drawAzBeam(azCenter,beamHalfwidth,range,[],[],[],[],'-k','linewidth',2)
%
%January 2023 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(systemType))
    systemType=0;
end

if(nargin<5||isempty(lRx))
    lRx=[0;0];
end

if(nargin<6||isempty(M))
    M=eye(2,2);
end

if(nargin<7||isempty(closeEnds))
    closeEnds=true;
end

numPts=1000;

if(isscalar(range))
   %If only a maximum range is given, set the minimum range to 0.
   range=[0;range];
end

if(systemType==1)
    azCenter=pi/2-azCenter;
end

if(beamHalfwidth<=0)
    error('beamHalfwidth must be positive.')
end

%Save the value of hold on the plot so that it can be reset to its previous
%state after plotting multiple lines.
holdVal=ishold();

systemType=0;
useHalfRange=true;

if(closeEnds==false)
    %Draw the sides
    az1=azCenter-beamHalfwidth;
    startPoint=pol2Cart([range(1);az1],systemType,useHalfRange,lRx,lRx,M);
    endPoint=pol2Cart([range(2);az1],systemType,useHalfRange,lRx,lRx,M);
    plot([startPoint(1);endPoint(1)],[startPoint(2);endPoint(2)],varargin{:})

    hold on
    az2=azCenter+beamHalfwidth;
    startPoint=pol2Cart([range(1);az2],systemType,useHalfRange,lRx,lRx,M);
    endPoint=pol2Cart([range(2);az2],systemType,useHalfRange,lRx,lRx,M);
    plot([startPoint(1);endPoint(1)],[startPoint(2);endPoint(2)],varargin{:})

    %Restore the hold value.
    if(~holdVal)
        hold off
    end
    return
end

if(range(1)~=0)
    totalNumPoints=2*numPts+1;
else
    totalNumPoints=numPts+2;
end
zPts=zeros(2,totalNumPoints);

az1=azCenter-beamHalfwidth;
az2=azCenter+beamHalfwidth;
ang=[az1+linspace(0,2*beamHalfwidth,numPts-1),az2];
zPts(:,1:numPts)=pol2Cart([range(2)*ones(1,numPts);ang],systemType,useHalfRange,lRx,lRx,M);

if(range(1)~=0)
    ang=flip(ang);
    zPts(:,(numPts+1):(2*numPts))=pol2Cart([range(1)*ones(1,numPts);ang],systemType,useHalfRange,lRx,lRx,M);
else
    zPts(:,numPts+1)=lRx;
end
zPts(:,totalNumPoints)=zPts(:,1);
plot(zPts(1,:),zPts(2,:),varargin{:})

%Restore the hold value to its original setting.
if(~holdVal)
    hold off
end
end

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


text:
function strip_svn_dirs()
% Removes the .svn dirs from the path
    strip_path_dirs('\.svn');
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function plotbernstein(B,a,b)
%PLOTBERNSTEIN  Plots Berstein points with connecting lines
%
%For univariate Bernstein points,
%
%   plotbernstein(B)
%
%plots Bernstein points in default interval [0,1]. Correspondingly,
%
%   plotbernstein(B,a,b)
%
%plots Bernstein points in the interval [a,b].
%
%For given polynomial P and an interval [a,b], a typical call for plotting P together with
%its Bernstein points is
%
%   B = bernsteincoeff(ptrans(P,a,b,0,1));
%   plotpoly(P,a,b), hold on, plotbernstein(B,a,b), hold off
%

% written  12/25/02     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if ~isreal(B.c)
    error('polynomial must be real (point or interval)')
  end

  if size(B.e,2)==1                   % univariate polynomial

    if nargin==1
      a = 0;
      b = 1;
    end
    n = length(B.c);
    X = linspace(a,b,n);

    if isa(B.c,'intval')

      Bc = fliplr([B.c.inf B.c.sup]);
      XX = [X X];
      if n>1
        index = convhull(XX,Bc);
      else
        index = 1:2;
      end
      plot(XX,Bc,'o',XX(index),Bc(index),'-o')

    else

      Bc = fliplr(B.c);
      if n>2
        index = convhull(X,Bc);
      else
        index = 1:length(Bc);
      end
      plot(X,Bc,'o',X(index),Bc(index),'-o')

    end

  elseif size(B.e,2)==2                  % multivariate polynomial in two unknowns
    error('not yet implemented')
  else
    error('Bernstein plot only for polynomials in one or two unknowns')
  end

  setround(rndold)

repo: OpenMAS


text:
function display(ker)

% DISPLAY
%
% Display a textual representation of a linear kernel object.
%
%    display(ker);

%
% File        : @linear/display.m
%
% Date        : Tuesday 12th September
%
% Author      : Dr Gavin C. Cawley
%
% Description : Display method for a class implementing a linear kernel.
%               Part of an object-oriented implementation of Vapnik's Support
%               Vector Machine, as described in [1].
%
% References  : [1] V.N. Vapnik,
%                   "The Nature of Statistical Learning Theory",
%                   Springer-Verlag, New York, ISBN 0-387-94559-8,
%                   1995.
%
% History     : 12/09/2000 - v1.00
%
% Copyright   : (c) Dr Gavin C. Cawley, September 2000.
%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
%

fprintf(1,'\nlinear kernel\n\n');

% bye bye...


text:
function flag = isProblem(problem)
% Check the input problem has all the mandatory fields
mandatory_fields = {'type', 'N', 'outliers', 'priors', 'dof'};
flag = all(isfield(problem, mandatory_fields));
end

author: MIT-SPARK

repo: GNC-and-ADAPT


text:
%PARZENML Optimum smoothing parameter in Parzen density estimation.
%
%   H = PARZENML(A)
%
% INPUT
%   A    Input dataset
%
% OUTPUT
%   H    Scalar smoothing parameter (in case of crisp labels)
%        Vector with smoothing parameters (in case of soft labels)
%
% DESCRIPTION
% Maximum likelihood estimation for the smoothing parameter H in the
% Parzen denstity estimation of the data in A. A leave-one out
% maximum likelihood estimation is used.
%
% The dataset A can either be crisp or soft labeled. In case of crisp
% labeling the class information is not used and a single smoothing
% parameter is estimated. In case of soft labels a smoothing parameter
% for every class is estimated and objects are weighted in relation to
% their class weigthts (soft label value).
% It may be profitable to scale the data before calling it. eg.
% WS = SCALEM(A,'variance'); A = A*WS.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, SCALEM, SELDAT, PARZENM, PARZENDC, PRPROGRESS

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: parzenml.m,v 1.11 2010/03/25 15:39:46 duin Exp $

function h = parzenml(A,fid)

	if nargin < 2, fid = []; end

	if isdouble(A), A = prdataset(A); end

	A = testdatasize(A);
	A = testdatasize(A,'objects');

	if islabtype(A,'crisp')
		h = parzenmlc(A,fid);
	elseif islabtype(A,'soft')
		h = parzenmls(A,fid);
	else
		error('Label type should be either ''crisp'' or ''soft''')
	end

	return

function h = parzenmlc(A,fid) %crisp version

	[m,k] = size(A);
	DD= distm(+A) + diag(1e70*ones(1,m));
	E = min(DD);

	h1 = sqrt(max(E));    % initial estimate of h
	F1 = derlc(DD,E,h1,k); % derivative

	prprogress(fid,'parzenml:\n');
	prprogress(fid,' %6.4f   %6.3e\n',h1,F1);
	if abs(F1) < 1e-70
		h = h1;
		prwarning(4,'jump out\n');
		return;
	end

	a1 = (F1+m*k)*h1*h1;
	h2 = sqrt(a1/(m*k));  % second guess
	F2 = derlc(DD,E,h2,k); % derivative

	prprogress(fid,' %6.4f   %6.3e\n',h2,F2);
	if (abs(F2) < 1e-70) | (abs(1e0-h1/h2) < 1e-6)
		h = h2;
		prwarning(4,'jump out\n');
		return
	end

	% find zero-point of derivative to optimize h^2
	% stop if improvement is small, or h does not change significantly

	alf = 1;
	prwaitbar(100,'parzenml: Optimizing smoothing parameter',m > 100)
	iter = 0;
	while abs(1e0-F2/F1) > 1e-4 & abs(1e0-h2/h1) > 1e-3 & abs(F2) > 1e-70
		iter = iter+1;
		h3 = (h1*h1*h2*h2)*(F2-F1)/(F2*h2*h2-F1*h1*h1);
		if h3 < 0 % this should not happen
			h3 = sqrt((F2+m*k)*h2*h2/(m*k));
		else
			h3 = sqrt(h3);
		end
		prwaitbar(100,100-100*exp(-iter/10));
		h3 = h2 +alf*(h3-h2);
		F3 = derlc(DD,E,h3,k);
		prprogress(fid,' %6.4f   %6.3e\n',h3,F3);
		F1 = F2; F2 = F3;
		h1 = h2; h2 = h3;
		alf = alf*0.99; % decrease step size
	end
	h = h2;
	prwaitbar(0);

return

function F = derlc(DD,E,h,k) % crisp version
	% computation of the likelihood derivative for Parzen density
	% given distances D and their object minima E (for increased accuracy)
	m = size(DD,1);
	warning off MATLAB:divideByZero;
		Y = (DD-repmat(E,m,1))/(2*h*h); % correct for minimum distance to save accuracy
	warning on MATLAB:divideByZero;
	IY = find(Y<20);                % take small distance only, others don't contribute
	P = zeros(m,m);
	P(IY) = exp(-Y(IY));
	PP = sum(P,2)';
	FU = repmat(realmax,1,m);
	J = find(PP~=0);
	FU(J) = 1./PP(J);
	FF = sum(DD.*P,2);
	warning off MATLAB:divideByZero;
		F = (FU*FF)./(h*h) - m*k;
	warning on MATLAB:divideByZero;
return

function h = parzenmls(A,fid) %soft version

	SS = gettargets(setlabtype(A,'soft'));
	[m,k,c] = getsize(A);
	DD= distm(+A) + diag(1e70*ones(1,m));
	E = min(DD);
	h = zeros(c,1);
	h0 = sqrt(max(E));    % initial estimate of h

	s = sprintf('parzenml: runover classes');
	prwaitbar(c,s,m > 100);
	iter = 0;

	for j=1:c
		prwaitbar(c,j)
		S = SS(:,j);
		h1 = h0;
		F1 = derls(DD,E,h1,k,S); % derivative

	  prprogress(fid,'parzenml: class %i : \n',j);
		prprogress(fid,' %6.4f   %6.3e\n',h1,F1);
		if abs(F1) < 1e-70
			h(j) = h1;
			prwarning(4,'jump out\n');
			break;
		end

		a1 = (F1+m*k)*h1*h1;
		h2 = sqrt(a1/(m*k));  % second guess
		F2 = derls(DD,E,h2,k,S); % derivative

		prprogress(fid,' %6.4f   %6.3e\n',h2,F2);
		if (abs(F2) < 1e-70) | (abs(1e0-h1/h2) < 1e-6)
			h(j) = h2;
			prwarning(4,'jump out\n');
			break;
		end

		% find zero-point of derivative to optimize h^2
		% stop if improvement is small, or h does not change significantly

		prwaitbar(100,'parzenml: Optimizing smoothing parameter',m > 100)
		iter = 0;
		alf = 1;
		while abs(1e0-F2/F1) > 1e-4 & abs(1e0-h2/h1) > 1e-3 & abs(F2) > 1e-70
			iter = iter+1;
			prwaitbar(100,100-100*exp(-iter/10));
			h3 = (h1*h1*h2*h2)*(F2-F1)/(F2*h2*h2-F1*h1*h1);
			if h3 < 0 % this should not happen
				h3 = sqrt((F2+m*k)*h2*h2/(m*k));
			else
				h3 = sqrt(h3);
			end
			h3 = h2 +alf*(h3-h2);
			F3 = derls(DD,E,h3,k,S);
			prprogress(fid,' %6.4f   %6.3e\n',h3,F3);
			F1 = F2; F2 = F3;
			h1 = h2; h2 = h3;
			alf = alf*0.99; % decrease step size
		end
		prwaitbar(0)
		h(j) = h2;
	end
  prwaitbar(0)
return

function F = derls(DD,E,h,k,S) %soft version
	% computation of the likelihood derivative for Parzen density
	% given distances D and their object minima E (for increased accuracy)
	% S are the object weigths
	c = size(S,2);                  % number of classes
	m = size(DD,1);
	Y = (DD-repmat(E,m,1))/(2*h*h); % correct for minimum distance to save accuracy
	IY = find(Y<20);                % take small distance only, others don't contribute
	F = 0;
	for j=1:c
		P = zeros(m,m);
		P(IY) = exp(-Y(IY));
		PP = S(:,j)'*P';
		FU = repmat(realmax,1,m);
		J = find(PP~=0);
		FU(J) = S(J,j)'./PP(J);
		K = find(S(:,j)==0);
		FU(K) = zeros(1,length(K));
		FF = (DD.*P)*S(:,j);
		F = F + (FU*FF)./(h*h);
	end
	F = F - sum(S(:))*k;
return

author: marianux

repo: ecg-kit

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


text:
function cin_lta()
    %

    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    report_this_filefun();

    var1 = zeros(1,ncu);
    var2 = zeros(1,ncu);
    lta = zeros(1,ncu);
    maxlta = zeros(1,ncu);
    maxlta = maxlta -5;

    cu = [cumuall(1:ti-1,:) ; cumuall(ti+winlen_days+1:len,:)];
    mean1 = mean(cu(:,:));
    mean2 = mean(cumuall(it:it+winlen_days,:));
    it

    for i = 1:ncu
        var1(i) = cov(cu(:,i));
    end     % for i

    for i = 1:ncu
        var2(i) = cov(cumuall(it:it+winlen_days,i));
    end     % for i

    lta = (mean1 - mean2)./(sqrt(var1/(len-winlen_days)+var2/(winlen_days)));
    valueMap = reshape(lta,length(gy),length(gx));

    % define size of the plot etc.
    %
    % set values gretaer ZG.tresh_km = nan
    %
    %[len, ncu] = size(cumuall);
    s = cumuall(len,:);
    r = reshape(s,length(gy),length(gx));
    re4 = valueMap;
    re4(r > ZG.tresh_km) = nan;

    figure(tmp);
    clf reset
    rect = [0.10 0.30 0.55 0.50 ];
    rect1 = rect;

    % find max and min of data for automatic scaling
    %

    % plot image
    %
    orient landscape
    axes('position',rect)
    pco1 = pcolor(gx,gy,re4);
    shading interp
    caxis([ZG.minc ZG.maxc]);
    colormap(jet)
    set(gca,'NextPlot','add')
    % plot overlay
    %
    overlay

    tx2 = text(0.07,0.85 ,['ti=' char(it*ZG.bin_dur+t0b)  ] ,...
        'Units','Norm','FontSize',14,'Color','k','FontWeight','bold');

    tx = text(0.07,0.95 ,['LTA;' char(ZG.compare_window_dur_v3)] ,...
        'Units','Norm','FontSize',14,'Color','k','FontWeight','bold');

    has = gca;

end

author: CelsoReyes


text:
function f = uminus(f)
%UMINUS   Unary minus of a SINGFUN.
%   UMINUS(F) is the negative of F.
%
% See also UPLUS, MINUS.

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

% Negate the smooth part:
f.smoothPart = -f.smoothPart;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function res = bf_output_image_powcorr(BF, S)
% Computes phase-amplitude coupling
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak using code bits from OSL library
% $Id: bf_output_image_powcorr.m 7703 2019-11-22 12:06:29Z guillaume $

%--------------------------------------------------------------------------
if nargin == 0
    all = cfg_const;
    all.tag = 'all';
    all.name = 'All';
    all.val  = {1};

    condlabel = cfg_entry;
    condlabel.tag = 'condlabel';
    condlabel.name = 'Condition label';
    condlabel.strtype = 's';
    condlabel.val = {''};

    conditions = cfg_repeat;
    conditions.tag = 'conditions';
    conditions.name = 'Conditions';
    conditions.help = {'Specify the labels of the conditions to be included in the inversion'};
    conditions.num  = [1 Inf];
    conditions.values  = {condlabel};
    conditions.val = {condlabel};

    whatconditions = cfg_choice;
    whatconditions.tag = 'whatconditions';
    whatconditions.name = 'What conditions to include?';
    whatconditions.values = {all, conditions};
    whatconditions.val = {all};

    sametrials = cfg_menu;
    sametrials.tag = 'sametrials';
    sametrials.name = 'Trials same as for filters';
    sametrials.labels = {'yes', 'no'};
    sametrials.values = {true, false};
    sametrials.val = {false};
    sametrials.help = {'Take the same trials as used for filter computation',...
        'This is useful for bootstrap.'};

    woi = cfg_entry;
    woi.tag = 'woi';
    woi.name = 'Time window of interest';
    woi.strtype = 'r';
    woi.num = [1 2];
    woi.val = {[-Inf Inf]};
    woi.help = {'Time windows (in ms)'};

    freqref         = cfg_entry;
    freqref.tag     = 'freqref';
    freqref.name    = 'Reference frequencies';
    freqref.strtype = 'r';
    freqref.num     = [1 Inf];
    freqref.val     = {20};
    freqref.help    = {'Frequencies in the reference channel'};

    resref         = cfg_entry;
    resref.tag     = 'resref';
    resref.name    = 'Reference resolutions';
    resref.strtype = 'r';
    resref.num     = [1 Inf];
    resref.val     = {5};
    resref.help    = {'Frequency resolution for reference frequencies. Single value or vector'};

    freq         = cfg_entry;
    freq.tag     = 'freq';
    freq.name    = 'Data frequencies';
    freq.strtype = 'r';
    freq.num     = [1 Inf];
    freq.val     = {20};
    freq.help    = {'First set of frequencies'};

    res         = cfg_entry;
    res.tag     = 'res';
    res.name    = 'Data resolutions';
    res.strtype = 'r';
    res.num     = [1 Inf];
    res.val     = {5};
    res.help    = {'Frequency resolution for data frequencies. Single value or vector'};

    refchan = cfg_entry;
    refchan.tag = 'refchan';
    refchan.name = 'Reference channel';
    refchan.strtype = 's';
    refchan.num = [1 Inf];
    refchan.help = {'Reference channel name.'};

    movavg         = cfg_entry;
    movavg.tag     = 'movavg';
    movavg.name    = 'Moving average window';
    movavg.strtype = 'r';
    movavg.num     = [1 1];
    movavg.val     = {100};
    movavg.help    = {'Time window for moving average of power envelope (ms).',...
        'Specify 0 to not average'};

    movavg         = cfg_entry;
    movavg.tag     = 'movavg';
    movavg.name    = 'Moving average window';
    movavg.strtype = 'r';
    movavg.num     = [1 1];
    movavg.val     = {100};
    movavg.help    = {'Time window for moving average of power envelope (ms).',...
        'Specify 0 to not average'};

    subsample = cfg_entry;
    subsample.tag = 'subsample';
    subsample.name = 'Subsample';
    subsample.strtype = 'n';
    subsample.num = [1 1];
    subsample.val = {1};
    subsample.help = {'Set to N to subsample the power to every Nth sample'};

    shuffle         = cfg_entry;
    shuffle.tag     = 'shuffle';
    shuffle.name    = 'Shuffle';
    shuffle.strtype = 'w';
    shuffle.num     = [1 1];
    shuffle.help    = {'Shuffle the reference channel to produce the null case.',...
        'Specify the number of shufflings'};
    shuffle.val = {0};

    modality         = cfg_menu;
    modality.tag     = 'modality';
    modality.name    = 'Modality';
    modality.help    = {'Specify modality'};
    modality.labels  = {
        'MEG'
        'MEGPLANAR'
        'EEG'
        }';
    modality.values  = {
        'MEG'
        'MEGPLANAR'
        'EEG'
        }';
    modality.val = {'MEG'};

    image_powcorr      = cfg_branch;
    image_powcorr.tag  = 'image_powcorr';
    image_powcorr.name = 'Power correlations image';
    image_powcorr.val  = {whatconditions, sametrials, shuffle, woi, refchan, freqref, ....
        resref, freq, res, movavg, subsample, modality};

    res = image_powcorr;

    return
elseif nargin < 2
    error('Two input arguments are required');
end

D = BF.data.D;

S.woi = 1e-3*S.woi; % ms -> s

samples =  D.indsample(S.woi(1)):D.indsample(S.woi(2));
nsamples = length(samples);
times = D.time(samples);

if isfield(S.whatconditions, 'all')
    S.whatconditions.condlabel = D.condlist;
end

for i = 1:numel(S.whatconditions.condlabel)
    if S.sametrials
        trials{i} = BF.features.trials(strmatch(S.whatconditions.condlabel{i},...
            D.conditions(BF.features.trials)));
    else
        trials{i} = D.indtrial(S.whatconditions.condlabel{i}, 'GOOD');
    end

    if isempty(trials{i})
        error('No trials matched the selection.');
    end

end

if isempty(trials)
    error('No trials matched the selection, check the specified condition labels');
end

channels = BF.features.(S.modality).chanind;
U        = BF.features.(S.modality).U;
nchan    = size(U, 2);

alltrials = spm_vec(trials);
ntrials   = length(alltrials);

nref = length(S.freqref);
nfreq    = length(S.freq);

W = BF.inverse.(S.modality).W;
nvert = numel(W);

Y = U'*reshape(D(channels, samples, alltrials), nchan, []);
Y = reshape(Y, size(Y, 1), nsamples, ntrials);

Yr = squeeze(D(D.indchannel(S.refchan), samples, alltrials));
if size(Yr, 1) == 1
    Yr = Yr';
end

spectrum = ft_specest_hilbert(Yr', times,...
    'freqoi', S.freqref, 'width', S.resref, 'filttype', 'but', 'filtorder', 2,...
    'filtdir', 'twopass', 'verbose', 0);

spectrum = reshape(permute(spectrum, [3 1 2]), nsamples, ntrials*nref);
spectrum = abs(spectrum);

if S.movavg
    avwin =  spm_hanning(1e-3*(S.movavg)*D.fsample);
    spectrum = conv2(avwin, 1, spectrum, 'same');
end

spectrum = spectrum(1:S.subsample:end, :);

spectrum = detrend(spectrum);
spectrum = spectrum./repmat(std(spectrum), size(spectrum, 1), 1);

refsig   = reshape(spectrum, [], ntrials, nref);

powc   = nan(length(S.freqref), length(S.freq), nvert);

if S.shuffle
    spowc = repmat(powc, [1 1 1 S.shuffle]);
    sind = zeros(S.shuffle, ntrials);
    for s = 1:S.shuffle
        sind(s, :) = randperm(ntrials);
    end
end

for f = 1:length(S.freq)

    spm_progress_bar('Init', ntrials, ...
        sprintf('Computing data spectra')); drawnow;
    if ntrials > 100, Ibar = floor(linspace(1, ntrials,100));
    else Ibar = 1:ntrials; end

    Yh = 0*Y;
    for i = 1:ntrials
        Yh(: , : ,i) = spm_squeeze(ft_specest_hilbert(squeeze(Y(:,:, i)), times,...
            'freqoi', S.freq(f), 'width', S.res(f), 'filttype', 'but', ...
            'filtorder', 2,  'filtdir', 'twopass', 'verbose', 0), 2);

        if ismember(i, Ibar)
            spm_progress_bar('Set', i); drawnow;
        end
    end

    Yh = reshape(Yh, nchan, []);

    spm_progress_bar('Clear');

    spm_progress_bar('Init', nvert, ...
        sprintf('Scanning grid points image')); drawnow;
    if nvert > 100, Ibar = floor(linspace(1, nvert,100));
    else Ibar = 1:nvert; end

    for i = 1:nvert
        if ~isnan(W{i})
            w    = W{i};

            sYh  = w*Yh;

            sYh  = reshape(abs(sYh), nsamples, ntrials);

            if S.movavg
                sYh = conv2(avwin, 1, sYh, 'same');
            end

            sYh = sYh(1:S.subsample:end, :);

            sYh = detrend(sYh);

            sYh = sYh./repmat(std(sYh), size(sYh, 1), 1);

            x = sYh(:);

            for j = 1:nref

                rYh = spm_squeeze(refsig(:,:, j), 3);

                for shuffle = 0:S.shuffle
                    if shuffle
                        rYh = rYh(:, sind(shuffle, :));
                    end

                    y = rYh(:);

                    pinvx = pinv(x);
                    pe  = pinvx*y;
                    r   = y-x*pe;
                    vr  = diag(r'*r/(size(y,1)-size(x,2)));
                    vrp = pinv(x'*x)*vr;
                    cs  = pe/sqrt(vrp);

                    if shuffle
                        spowc(j, f, i, shuffle) = cs;
                    else
                        powc(j, f, i) = cs;
                    end
                end
            end
        end

        if ismember(i, Ibar)
            spm_progress_bar('Set', i); drawnow;
        end
    end
end

spm_progress_bar('Clear');

if max(nfreq, nref)>1
    image(1).val     = squeeze(sum(sum(powc, 2), 1));
    image(1).label   = ['sumpowc_'  spm_file(D.fname, 'basename')];
    c = 2;
else
    c = 1;
end

for f = 1:nref
    for g = 1:nfreq
        image(c).val     = squeeze(powc(f, g, :));
        image(c).label   = ['powc_ref_' num2str(S.freqref(f)) 'Hz_meg_'...
            num2str(S.freq(g)) 'Hz_' spm_file(D.fname, 'basename')];
        c = c+1;
    end
end

for shuffle = 1:S.shuffle
    if max(nfreq, nref)>1
        image(c).val     = squeeze(sum(sum(spowc(:,:,:, shuffle), 2), 1));
        image(c).label   = ['shuffled' num2str(shuffle) '_sumpowc_'  spm_file(D.fname, 'basename')];
        c = c+1;
    end

    for f = 1:nref
        for g = 1:nfreq
            image(c).val     = squeeze(spowc(f, g, :, shuffle));
            image(c).label   = ['shuffled' num2str(shuffle) '_powc_ref_' num2str(S.freqref(f)) 'Hz_meg_'...
                num2str(S.freq(g)) 'Hz_' spm_file(D.fname, 'basename')];
            c = c+1;
        end
    end
end

res = image;

author: spm


text:
function result = clara(x,kclus,vtype,stdize,metric,nsamp,sampsize)

%CLARA is the 'Clustering Large Applications' clustering algorithm.
% It returns a list representing a clustering of the data
% into kclus clusters following the clara algorithm which is
% designed for large data sets.
%
%The algorithm is fully described in:
%   Kaufman, L. and Rousseeuw, P.J. (1990),
%   "Finding groups in data: An introduction to cluster analysis",
%   Wiley-Interscience: New York (Series in Applied Probability and
%   Statistics), ISBN 0-471-87876-6.
%
% Required input arguments:
%       x : Data matrix (rows = observations, columns = variables)
%   kclus : The number of desired clusters
%   vtype : Variable type vector (length equals number of variables)
%           Possible values are 1  Asymmetric binary variable (0/1)
%                               2  Nominal variable (includes symmetric binary)
%                               3  Ordinal variable
%                               4  Interval variable
%
% Optional input arguments:
%     stdize : standardise the variables given by the x-matrix
%              Possible values are 0 : no standardisation (default)
%                                  1 : standardisation by the mean
%                                  2 : standardisation by the median
%     metric : Metric to be used
%              Possible values are 'eucli' Euclidian (all interval variables, default)
%                                  'manha' Manhattan
%                                  'mixed' Mixed (not all interval variables, default)
%      nsamp : Number of samples to be drawn from the data set
%   sampsize : Number of observations in each sample (should be higher
%              than the number of clusters and lower than the number of
%              observations)
%
% I/O:
%   result=clara(x,kclus,vtype,'eucli',5,40+2*kclus)
%
% Example (subtracted from the referenced book)
%   load obj200.mat
%   result=clara(obj200,3,[4 4]);
%
% The output of CLARA is a structure containing:
%   result.dysobs     : dissimilarities for each observation with the medoids
%   result.metric     : metric used
%   result.number     : number of observations
%   result.idmed      : Id of medoid observations
%   result.ncluv      : A vector with length equal to the number of observations,
%                       giving for each observation the number of the cluster to
%                       which it belongs
%   result.obj        : Objective function for the best subsample
%   result.clusinf    : Matrix, each row gives numerical information for
%                       one cluster. These are the cardinality of the cluster
%                       (number of observations), the maximal and average
%                       dissimilarity between the observations in the cluster
%                       and the cluster's medoid, the diameter of the cluster
%                       (maximal dissimilarity between two observations of the
%                       cluster), and the separation of the cluster (minimal
%                       dissimilarity between an observation of the cluster
%                       and an observation of another cluster).
%   result.sylinf     : Matrix based on the best subsample, with for each
%                       observation i of this subsample the cluster to
%                       which i belongs, as well as the neighbor cluster of i
%                       (the cluster, not containing i, for which the average
%                       dissimilarity between its observations and i is minimal),
%                       and the silhouette width of i.
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by Guy Brys (May 2006)

%Checking and filling out the inputs
if (nargin<3)
    error('Three input arguments required')
elseif (nargin<4)
    stdize = 0;
    if (sum(vtype)~=4*size(x,2))
        metric = 'mixed';
    else
        metric = 'eucli';
    end
    nsamp=5;
    sampsize=40+2*kclus;
elseif (nargin<5)
    if (sum(vtype)~=4*size(x,2))
        metric = 'mixed';
    else
        metric = 'eucli';
    end
    nsamp=5;
    sampsize=40+2*kclus;
elseif (nargin<6)
    nsamp=5;
    sampsize=40+2*kclus;
elseif (nargin<7)
    sampsize=40+2*kclus;
end

%Standardization
if (stdize==1)
    x = ((x - repmat(mean(x),size(x,1),1))./(repmat(std(x),size(x,1),1)));
elseif (stdize==2)
    x = ((x - repmat(median(x),size(x,1),1))./(repmat(mad(x),size(x,1),1)));
end

%Actual calculations
obj = Inf;
for i=1:nsamp
    sampindex = randperm(size(x,1));
    restemp = pam(x(sampindex(1:sampsize),:),kclus,vtype,metric);
    if (restemp.obj(1)<obj)
        obj = restemp.obj(1);
        idmed = sampindex(restemp.idmed);
    end
end

%Calculating some extra dissimilarities for output
for i=1:size(x,1)
    for j=1:kclus
        distemp = daisy(x([i idmed(j)],:),vtype,metric);
        disv(i,j) = distemp.disv(1);
    end
    [zz,clu(i)] = min(disv(i,:));
end
mindisv=[];
for j=1:kclus
    clusinf(j,1) = length(find(clu==j));
    clusinf(j,2) = max(disv(find(clu==j),j));
    clusinf(j,3) = mean(disv(find(clu==j),j));
end
for i=1:kclus
    for j=1:kclus
        distemp = daisy(x([idmed(i) idmed(j)],:),vtype,metric);
        mindisv = [mindisv distemp.disv(1)];
    end
end
clusinf(:,4) = clusinf(:,2)/min(mindisv(mindisv~=0));

%Putting things together
result = struct('dysobs',disv,'metric',metric,'number',size(x,1),...
    'idmed',idmed,'ncluv',clu,'obj',obj,'clusinf',clusinf,...
    'sylinf',restemp.sylinf,'x',x);


text:
%% Outlier demo
%Copyright (c) 2011, The MathWorks, Inc.

clear all
clc

% Generate a vector of X and Y variables
[CleanX, CleanY] = meshgrid(1:10);
CleanX = reshape(CleanX,100,1);
CleanY = reshape(CleanY,100,1);

%% Specify Z as a function of X and Y
A = 3; B = 4;
CleanZ = A*CleanX + B*CleanY;

%% Add in a noise vector
noise = 2*randn(100,1);
Noisy_Z = CleanZ + noise;

% Plot
scatter3(CleanX, CleanY, Noisy_Z, 'filled')
hold on

%% Predict Z as a function of X and Y
[foo gof] = fit([CleanX, CleanY], Noisy_Z, 'poly11')
plot(foo)

%% Use Robust Regression

bar = fit([CleanX, CleanY], Noisy_Z, 'poly11', 'Robust', 'LAR')

%% Use robustfit

b = robustfit([CleanX, CleanY], Noisy_Z)

author: Sable

repo: mcbench-benchmarks


text:
filename = 'CantileverSquare';%'CantileverSquareSYmmetricMesh';%'CantileverSquareNewFine';%'LshapeTriFine';%'CantileverSquareNew';'ArchTriFine';%'CantileverSquare';%'ArchTriFine';%'CantileverSquare';%'Lshape';%'LshapeTriFine';%'CantileverSquareSmall';%'';%'Lshape';'CantileverSquare';'Lshape';%'CantileverSquare';%'LshapeFine';%'Bridge_quad_coarse';%'Arch_quad_coarse';%'BridgeCool_Quadrilateral_Bilinear_Structured_Coarse';%'Bridge';%'CantileverSquareSmall';
ptype = 'MACRO';
initial_case = 'given';
m1 = 0.0101;
m2 = 0.0101;
%cost = {'compliance'};
cost = {'stressNorm'};
%cost = {'stressNorm','compliance'};
%weights = [0.55,0.45];
weights = 1;
constraint = {'volumeConstraint'};
filterType = 'P1';
constraint_case = 'EQUALITY';

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

Vfrac_final = 0.3;
optimality_final = 1e-5;
constr_final = 1e-5;

stressNormExponent_initial = 2;
stressNormExponent_final = 64;

optimizer = 'DualNestedInPrimal';
optimizerUnconstrained = 'PROJECTED GRADIENT';

designVariable = 'MicroParams';
ub = 0.989;
lb = 0.011;
homegenizedVariablesComputer = 'ByVademecum';
%
%vademecumFileName = 'SuperEllipseQMax';
%vademecumFileName = 'SuperEllipseQ2';
vademecumFileName = 'SuperEllipseQOptAnalytic';
%
% designVariable = 'Density';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';

kfrac = 2;
nsteps = 20;

plotting = false;
printing = true;
monitoring = true;
monitoring_interval = 1;
maxiter = 500;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%This Matlab script can be used to reproduce Figures 4.5-6 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%Number of BSs
L = 16;

%Number of UEs per BS
K = 10;

%Define the range of BS antennas
Mrange = 10:10:100;

%Extract maximum number of BS antennas
Mmax = max(Mrange);

%Define the range of pilot reuse factors
fRange = [1 2 4];

%Select the number of setups with random UE locations
nbrOfSetups = 100;

%Select the number of channel realizations per setup
nbrOfRealizations = 100;

%% Propagation parameters

%Communication bandwidth
B = 20e6;

%Total uplink transmit power per UE (mW)
p = 100;

%Noise figure at the BS (in dB)
noiseFigure = 7;

%Compute noise power
noiseVariancedBm = -174 + 10*log10(B) + noiseFigure;

%Select length of coherence block
tau_c = 200;

%Use the approximation of the Gaussian local scattering model
accuracy = 2;

%Angular standard deviation in the local scattering model (in degrees)
ASDdeg = 10;

%Prepare to save simulation results
sumSE_MR = zeros(length(Mrange),length(fRange),nbrOfSetups);
sumSE_ZF = zeros(length(Mrange),length(fRange),nbrOfSetups);
sumSE_RZF = zeros(length(Mrange),length(fRange),nbrOfSetups);
sumSE_SMMSE = zeros(length(Mrange),length(fRange),nbrOfSetups);
sumSE_MMMSE = zeros(length(Mrange),length(fRange),nbrOfSetups);

%% Go through all setups
for n = 1:nbrOfSetups

    %Output simulation progress
    disp([num2str(n) ' setups out of ' num2str(nbrOfSetups)]);

    %Compute channel statistics for one setup
    [R,channelGaindB] = functionExampleSetup(L,K,Mmax,accuracy,ASDdeg);

    %Compute the normalized average channel gain, where the normalization
    %is based on the noise power
    channelGainOverNoise = channelGaindB - noiseVariancedBm;

    %Go through all number of antennas
    for m = 1:length(Mrange)

        %Output simulation progress
        disp([num2str(m) ' antennas out of ' num2str(length(Mrange))]);

        %Go through all pilot reuse factors
        for s = 1:length(fRange)

            %Extract pilot reuse factor
            f = fRange(s);

            %Generate channel realizations with estimates and estimation
            %error correlation matrices
            [Hhat,C,tau_p,Rscaled] = functionChannelEstimates(R(1:Mrange(m),1:Mrange(m),:,:,:),channelGainOverNoise,nbrOfRealizations,Mrange(m),K,L,p,f);

            %Compute SEs using Theorem 4.1
            [SE_MR,SE_RZF,SE_MMMSE,SE_ZF,SE_SMMSE] = functionComputeSE_UL(Hhat,C,Rscaled,tau_c,tau_p,nbrOfRealizations,Mrange(m),K,L,p);

            %Save average sum SE per cell
            sumSE_MR(m,s,n) = mean(sum(SE_MR,1));
            sumSE_ZF(m,s,n) = mean(sum(SE_ZF,1));
            sumSE_SMMSE(m,s,n) = mean(sum(SE_SMMSE,1));
            sumSE_RZF(m,s,n) = mean(sum(SE_RZF,1));
            sumSE_MMMSE(m,s,n) = mean(sum(SE_MMMSE,1));

            %Delete large matrices
            clear Hhat C Rscaled;

        end

    end

    %Delete large matrices
    clear R;

end

%% Plot the simulation results
for s = 1:length(fRange)

    figure(s);
    hold on; box on;

    plot(Mrange,mean(sumSE_MMMSE(:,s,:),3),'rd-','LineWidth',1);
    plot(Mrange,mean(sumSE_SMMSE(:,s,:),3),'b:','LineWidth',1);
    plot(Mrange,mean(sumSE_RZF(:,s,:),3),'k-.','LineWidth',1);
    plot(Mrange,mean(sumSE_ZF(:,s,:),3),'r--','LineWidth',1);
    plot(Mrange,mean(sumSE_MR(:,s,:),3),'bs-','LineWidth',1);

    xlabel('Number of antennas (M)');
    ylabel('Average sum SE [bit/s/Hz/cell]');
    legend('M-MMSE','S-MMSE','RZF','ZF','MR','Location','NorthWest');
    ylim([0 60]);

end

author: emilbjornson

repo: massivemimobook


text:
classdef CoreCachedFileReader < handle
    % CoreCachedFileReader. A helper class to assist with reading a file
    %
    %
    %     Licence
    %     -------
    %     Part of CoreMat. https://github.com/tomdoel/coremat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        % The size of the text buffer to be read in. It is safe to change this
        % value while file reading is in progress
        MaxFileCacheSize = uint32(600000)
    end

    properties (Access = private)
        FileText
        FileId

        FilePointer
        FileCache
        SizeOfCurrentCache
        FlipEndian

        EofReached
    end

    methods
        function obj = CoreCachedFileReader(file_path, file_name, transfer_syntax, reporting)

            full_file_name = fullfile(file_path, file_name);

            computer_endian = CoreSystemUtilities.GetComputerEndian;
            obj.FlipEndian = (computer_endian ~= transfer_syntax.Endian);

            switch transfer_syntax.Endian
                case CoreEndian.BigEndian
                    file_encoding = 'b';
                case CoreEndian.LittleEndian
                    file_encoding = 'l';
                otherwise
                    error('Unknown file encoding');
            end

            switch transfer_syntax.CharacterEncoding
                case CoreCharacterEncoding.UTF8
                    character_encoding = 'UTF-8';
                otherwise
                    error('Unknown character encoding');
            end

            file_id = fopen(full_file_name, 'r', file_encoding, character_encoding);

            if (file_id == -1)
                reporting.Error('CoreCachedFileReader:CannotOpenFile', ['Unable to open file ' full_file_name]);
                obj.EofReached = true;
                obj.FileId = [];
            else
                obj.FileId = file_id;
                obj.EofReached = false;
                obj.ReadNextFileBlock;
            end
        end

        function data = ReadData(obj, data_type, number_of_bytes)
            data = ReadBytesFromCache(obj, number_of_bytes);
            data = typecast(data, data_type);
            if obj.FlipEndian
                data = swapbytes(data);
            end
        end

        function data = ReadString(obj, number_of_bytes)
            data = char(ReadBytesFromCache(obj, number_of_bytes));
        end

        function wrapper = ReadDataIntoWrapper(obj, data_type, number_of_bytes)
            wrapper = CoreWrapper;
            wrapper.RawImage = ReadBytesFromCache(obj, number_of_bytes);

            if strcmp(data_type, 'char')
                wrapper.RawImage = char(wrapper.RawImage);
            else
                wrapper.RawImage = typecast(wrapper.RawImage, data_type);
                if obj.FlipEndian
                    wrapper.RawImage = swapbytes(wrapper.RawImage);
                end
            end
        end

        function Close(obj)
            if ~isempty(obj.FileId)
                fclose(obj.FileId);
                obj.FileId = [];
            end
        end

        function delete(obj)
            obj.Close;
        end

        function more_data = MoreDataToRead(obj)
            more_data_in_cache = ~isempty(obj.FileCache) && (obj.FilePointer <= obj.SizeOfCurrentCache);
            more_data_to_read_from_disk = ~obj.EofReached;
            more_data = more_data_in_cache || more_data_to_read_from_disk;
        end

        function Skip(obj, number_of_bytes)
            bytes_left_to_read = uint32(number_of_bytes);
            bytes_left_in_cache = 1 + obj.SizeOfCurrentCache - obj.FilePointer;

            if bytes_left_in_cache >= bytes_left_to_read
                obj.FilePointer = obj.FilePointer + bytes_left_to_read;
            else

                while (bytes_left_to_read > 0) && (obj.MoreDataToRead)
                    bytes_left_in_cache = 1 + obj.SizeOfCurrentCache - obj.FilePointer;
                    next_bytes_read = min(bytes_left_in_cache, bytes_left_to_read);
                    bytes_left_to_read = bytes_left_to_read - next_bytes_read;
                    obj.FilePointer = obj.FilePointer + next_bytes_read;
                    if obj.FilePointer > obj.SizeOfCurrentCache
                        obj.ReadNextFileBlock;
                    end
                end

            end
        end
    end

    methods (Access = private)
        function data_bytes = ReadBytesFromCache(obj, number_of_bytes)
            bytes_left_to_read = uint32(number_of_bytes);
            bytes_left_in_cache = 1 + obj.SizeOfCurrentCache - obj.FilePointer;

            if bytes_left_in_cache >= bytes_left_to_read
                data_bytes = obj.FileCache(obj.FilePointer : obj.FilePointer + bytes_left_to_read - 1)';
                obj.FilePointer = obj.FilePointer + bytes_left_to_read;
            else
                data_bytes = zeros(1, number_of_bytes, 'uint8');
                data_bytes_pointer = uint32(1);

                while (bytes_left_to_read > 0) && (obj.MoreDataToRead)
                    bytes_left_in_cache = 1 + obj.SizeOfCurrentCache - obj.FilePointer;
                    next_bytes_read = min(bytes_left_in_cache, bytes_left_to_read);
                    data_bytes(data_bytes_pointer : data_bytes_pointer + next_bytes_read - 1) = obj.FileCache(obj.FilePointer : obj.FilePointer + next_bytes_read - 1);
                    data_bytes_pointer = data_bytes_pointer + next_bytes_read;
                    bytes_left_to_read = bytes_left_to_read - next_bytes_read;
                    obj.FilePointer = obj.FilePointer + next_bytes_read;
                    if obj.FilePointer > obj.SizeOfCurrentCache
                        obj.ReadNextFileBlock;
                    end
                end

                % If the end of the file has been reached, we may not have the full
                % number of bytes
                if bytes_left_to_read > 0
                    number_of_bytes_actually_read = number_of_bytes - bytes_left_to_read;
                    data_bytes = data_bytes(1 : number_of_bytes_actually_read);
                end
            end
        end

        function ReadNextFileBlock(obj)
            if obj.EofReached
                obj.FileCache = [];
            else
                obj.FileCache = fread(obj.FileId, obj.MaxFileCacheSize, 'uint8=>uint8');
                obj.SizeOfCurrentCache = uint32(numel(obj.FileCache));
                obj.FilePointer = uint32(1);
                if obj.SizeOfCurrentCache < obj.MaxFileCacheSize
                    obj.EofReached = true;
                    obj.Close;
                end
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
classdef MimScript < handle
    % MimScript. Base class for an API script used by the Pulmonary Toolkit.
    %
    %     MimScript are classes you create to run routines using the
    %     Pulmonary Toolkit API. Whereas Plugins execute for a single
    %     dataset and produce a deterministic result, Scripts can load and
    %     operate over multiple datasets. Unike Gui Plugins, Scripts do not
    %     access the user interface, so they can be run from the command
    %     window using the API's RunScript method.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Abstract = true)

        % This specifies the version of Script interface implemented by this script
        InterfaceVersion

        % A version number for this script implementation
        Version

        % Set to an apropriate value for grouping related scripts
        Category
    end

    methods (Abstract, Static)

        % Called when the user has clicked the button for this gui plugin.
        % Implement this method with the code you with to run.
        output = RunScript(ptk_obj, parameters)
    end
end


text:
% This make.m is for MATLAB and OCTAVE under Windows, Mac, and Unix

try
	Type = ver;
	% This part is for OCTAVE
	if(strcmp(Type(1).Name, 'Octave') == 1)
		mex libsvmread.c
		mex libsvmwrite.c
		mex train.c linear_model_matlab.c ../linear.cpp ../tron.cpp ../blas/*.c
		mex predict.c linear_model_matlab.c ../linear.cpp ../tron.cpp ../blas/*.c
	% This part is for MATLAB
	% Add -largeArrayDims on 64-bit machines of MATLAB
	else
		mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims libsvmread.c
		mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims libsvmwrite.c
		mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims train.c linear_model_matlab.c ../linear.cpp ../tron.cpp "../blas/*.c"
		mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims predict.c linear_model_matlab.c ../linear.cpp ../tron.cpp "../blas/*.c"
	end
catch
	fprintf('If make.m fails, please check README about detailed instructions.\n');
end

author: rbgirshick

repo: rcnn

save_path: github-repos/MATLAB/rbgirshick-rcnn


text:
function [ c, ic1, ic2, nc, an1, an2, ncp, d, itran, live, indx ] = ...
  optra ( a, m, n, c, k, ic1, ic2, nc, an1, an2, ncp, d, itran, live, indx )

%*****************************************************************************80
%
%% OPTRA carries out the optimal transfer stage.
%
%  Discussion:
%
%    This is the optimal transfer stage.
%
%    Each point is re-allocated, if necessary, to the cluster that
%    will induce a maximum reduction in the within-cluster sum of
%    squares.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2008
%
%  Author:
%
%    Original FORTRAN77 version by John Hartigan, Manchek Wong.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    John Hartigan, Manchek Wong,
%    Algorithm AS 136:
%    A K-Means Clustering Algorithm,
%    Applied Statistics,
%    Volume 28, Number 1, 1979, pages 100-108.
%
%  Parameters:
%
%    Input, real A(M,N), the points.
%
%    Input, integer M, the number of points.
%
%    Input, integer N, the number of spatial dimensions.
%
%    Input/output, real C(K,N), the cluster centers.
%
%    Input, integer K, the number of clusters.
%
%    Input/output, integer IC1(M), the cluster to which each
%    point is assigned.
%
%    Input/output, integer IC2(M), used to store the cluster
%    which each point is most likely to be transferred to at each step.
%
%    Input/output, integer NC(K), the number of points in
%    each cluster.
%
%    Input/output, real AN1(K).
%
%    Input/output, real AN2(K).
%
%    Input/output, integer NCP(K).
%
%    Input/output, real D(M).
%
%    Input/output, integer ITRAN(K).
%
%    Input/output, integer LIVE(K).
%
%    Input/output, integer INDX, the number of steps since a
%    transfer took place.
%

%
%  If cluster L is updated in the last quick-transfer stage, it
%  belongs to the live set throughout this stage.   Otherwise, at
%  each step, it is not in the live set if it has not been updated
%  in the last M optimal transfer steps.
%
  for l = 1 : k
    if ( itran(l) == 1)
      live(l) = m + 1;
    end
  end

  for i = 1 : m

    indx = indx + 1;
    l1 = ic1(i);
    l2 = ic2(i);
    ll = l2;
%
%  If point I is the only member of cluster L1, no transfer.
%
    if ( 1 < nc(l1)  )
%
%  If L1 has not yet been updated in this stage, no need to
%  re-compute D(I).
%
      if ( ncp(l1) ~= 0 )
        de = 0.0;
        for j = 1 : n
          df = a(i,j) - c(l1,j);
          de = de + df * df;
        end
        d(i) = de * an1(l1);
      end
%
%  Find the cluster with minimum R2.
%
     da = 0.0;
      for j = 1 : n
        db = a(i,j) - c(l2,j);
        da = da + db * db;
      end
      r2 = da * an2(l2);

      for l = 1 : k
%
%  If LIVE(L1) <= I, then L1 is not in the live set.   If this is
%  true, we only need to consider clusters that are in the live set
%  for possible transfer of point I.   Otherwise, we need to consider
%  all possible clusters.
%
        if ( ( i < live(l1) | i < live(l2) ) & l ~= l1 & l ~= ll )

          rr = r2 / an2(l);

          dc = 0.0;
          for j = 1 : n
            dd = a(i,j) - c(l,j);
            dc = dc + dd * dd;
          end

          if ( dc < rr )
            r2 = dc * an2(l);
            l2 = l;
          end

        end

      end
%
%  If no transfer is necessary, L2 is the new IC2(I).
%
      if ( d(i) <= r2 )

        ic2(i) = l2;
%
%  Update cluster centers, LIVE, NCP, AN1 and AN2 for clusters L1 and
%  L2, and update IC1(I) and IC2(I).
%
      else

        indx = 0;
        live(l1) = m + i;
        live(l2) = m + i;
        ncp(l1) = i;
        ncp(l2) = i;
        al1 = nc(l1);
        alw = al1 - 1.0;
        al2 = nc(l2);
        alt = al2 + 1.0;
        for j = 1 : n
          c(l1,j) = ( c(l1,j) * al1 - a(i,j) ) / alw;
          c(l2,j) = ( c(l2,j) * al2 + a(i,j) ) / alt;
        end
        nc(l1) = nc(l1) - 1;
        nc(l2) = nc(l2) + 1;
        an2(l1) = alw / al1;
        if ( 1.0 < alw )
          an1(l1) = alw / ( alw - 1.0 );
        else
          an1(l1) = Inf;
        end
        an1(l2) = alt / al2;
        an2(l2) = alt / ( alt + 1.0 );
        ic1(i) = l2;
        ic2(i) = l1;

      end

    end

    if ( indx == m )
      return
    end

  end
%
%  ITRAN(L) = 0 before entering QTRAN.   Also, LIVE(L) has to be
%  decreased by M before re-entering OPTRA.
%
  for l = 1 : k
    itran(l) = 0;
    live(l) = live(l) - m;
  end

  return
end


text:
function pass = test_wronskian(pref)

% Choose a tolerance:
tol = 1e-11;

% Check different input styles for a simple problem:
op = @(x,u) diff(u, 2) + (u);
L = chebop(op, [-pi, pi]);
x = chebfun('x', [-pi, pi]);
w = wronskian(L, cos(x), sin(x));
pass(1) = norm(w - 1, inf) < tol;
w = wronskian(L, [cos(x), sin(x)]);
pass(2) = norm(w - 1, inf) < tol;
w = wronskian(L, chebmatrix({cos(x), sin(x)}));
pass(3) = norm(w - 1, inf) < tol;

% A slightly harder problem:
op = @(u) diff(u,2) + 4*diff(u) - 5*u;
dom = [0, 5];
L = chebop(op, dom);
x = chebfun(@(x) x, dom);
g = exp(x);
f = exp(-5*x);
w = wronskian(L, f, g);
% Compare with the exact wronskian:
pass(4) = norm(w - 6*exp(-4*x), inf) < tol;

end


text:
function [datapt, unitStr] = lvd_GrdObjTasks(stateLogEntry, subTask, grdObj, inFrame)
%lvd_GrdObjTasks Summary of this function goes here
%   Detailed explanation goes here

    stateLogEntry = stateLogEntry.deepCopy();
    cartElem = stateLogEntry.getCartesianElementSetRepresentation().convertToFrame(inFrame);
    stateLogEntry.setCartesianElementSet(cartElem);

    switch subTask
        case 'azimuth'
            [az, ~, ~] = getAzElRngOfScFromGrdObj(stateLogEntry, grdObj);

            datapt = rad2deg(az);
            unitStr = 'deg';

        case 'elevation'
            [~, elev, ~] = getAzElRngOfScFromGrdObj(stateLogEntry, grdObj);

            datapt = rad2deg(elev);
            unitStr = 'deg';

        case 'range'
            [~, ~, r] = getAzElRngOfScFromGrdObj(stateLogEntry, grdObj);

            datapt = r;
            unitStr = 'km';

        case 'LoS'
            maStateLogEntry = stateLogEntry.getMAFormattedStateLogMatrix(false);
            bodyInfo = stateLogEntry.centralBody;
            celBodyData = bodyInfo.celBodyData;

            targetBodyInfo = grdObj.centralBodyInfo;

            allBodyInfo = celBodyData.getAllBodyInfo();
            hasLoSAll = true;
            for(j=1:length(allBodyInfo))
                eclipseBodyInfo = allBodyInfo(j);

                hasLoS = LoS2Target(maStateLogEntry, bodyInfo, eclipseBodyInfo, targetBodyInfo, celBodyData, grdObj);

                if(hasLoS == 0)
                    hasLoSAll = false;
                    break;
                end
            end

            datapt = double(hasLoSAll);
            unitStr = '';

        otherwise
            error('Unknown sub task string: %s', subTask);
    end
end

function [az, elev, r] = getAzElRngOfScFromGrdObj(stateLogEntry, grdObj)
    time = stateLogEntry.time;
    scCartElem = stateLogEntry.getCartesianElementSetRepresentation();
    grdObjElemSet = grdObj.getStateAtTime(time);

    grdObjBodyInfo = grdObj.centralBodyInfo;
    grdObjParentBodyInertialFrame = grdObj.centralBodyInfo.getBodyCenteredInertialFrame();
    scCartElemStnFrame = scCartElem.convertToFrame(grdObjParentBodyInertialFrame).convertToCartesianElementSet();

    grdObjElemSetInertialFrame = grdObjElemSet.convertToFrame(grdObjParentBodyInertialFrame).convertToCartesianElementSet();

    rVectScToTarget = grdObjElemSetInertialFrame.rVect - scCartElemStnFrame.rVect;
    stnRVectECIRelToParent = grdObjElemSetInertialFrame.rVect;

    R_ned_2_inert = computeNedFrame(time, stnRVectECIRelToParent, grdObjBodyInfo);
    rVectTargetToSc = -rVectScToTarget;
    rVectTargetToScNed = R_ned_2_inert' * rVectTargetToSc;
    [az, elev, r] = getAzElRngFromNedPosition(rVectTargetToScNed);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function Fmem = PQinitFMem (Nc, PCinit)
% Initialize the filter memories

% P. Kabal $Revision: 1.1 $  $Date: 2003/12/07 13:34:10 $

Fmem.TDS.Ef(1:2,1:Nc) = 0;
Fmem.Adap.P(1:2,1:Nc) = 0;
Fmem.Adap.Rn(1:Nc) = 0;
Fmem.Adap.Rd(1:Nc) = 0;
Fmem.Adap.PC(1:2,1:Nc) = PCinit;
Fmem.Env.Ese(1:2,1:Nc) = 0;
Fmem.Env.DE(1:2,1:Nc) = 0;
Fmem.Env.Eavg(1:2,1:Nc) = 0;

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
function f = powell ( m, x )

%*****************************************************************************80
%
%% POWELL computes the Powell singular quartic function.
%
%  Discussion:
%
%    This function has a global minimizer:
%
%      X* = ( 0.0, 0.0, 0.0, 0.0 ), F(X*) = 0.
%
%    Start the search at
%
%      X = ( 3.0, -1.0, 0.0, 1.0 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Michael Powell,
%    An Iterative Method for Finding Stationary Values of a Function
%    of Several Variables,
%    Computer Journal,
%    Volume 5, 1962, pages 147-151.
%
%  Parameters:
%
%    Input, integer M, the number of variables.
%
%    Input, real X(M), the argument of the function.
%
%    Output, real F, the value of the function at X.
%
  f1 = x(1) + 10.0 * x(2);
  f2 =                            x(3) - x(4);
  f3 =               x(2) - 2.0 * x(3);
  f4 = x(1)                            - x(4);

  f = f1 * f1 + f2 * f2 + f3 * f3 + f4 * f4;

  return
end


text:
%% Default template
%% Plot some Pole Figures

plotPDF(odf,[Miller(1,0,0),Miller(1,1,0)],'antipodal')

author: mtex-toolbox

repo: mtex

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


text:
function y = medfilt1m( x, r, z )
% One-dimensional adaptive median filtering with missing values.
%
% Applies a width s=2*r+1 one-dimensional median filter to vector x, which
% may contain missing values (elements equal to z). If x contains no
% missing values, y(j) is set to the median of x(j-r:j+r). If x contains
% missing values, y(j) is set to the median of x(j-R:j+R), where R is the
% smallest radius such that sum(valid(x(j-R:j+R)))>=s, i.e. the number of
% valid values in the window is at least s (a value x is valid x~=z). Note
% that the radius R is adaptive and can vary as a function of j.
%
% This function uses a modified version of medfilt1.m from Matlab's 'Signal
% Processing Toolbox'. Note that if x contains no missing values,
% medfilt1m(x) and medfilt1(x) are identical execpt at boundary regions.
%
% USAGE
%  y = medfilt1m( x, r, [z] )
%
% INPUTS
%  x      - [nx1] length n vector with possible missing entries
%  r      - filter radius
%  z      - [NaN] element that represents missing entries
%
% OUTPUTS
%  y      - [nx1] filtered vector x
%
% EXAMPLE
%  x=repmat((1:4)',1,5)'; x=x(:)'; x0=x;
%  n=length(x); x(rand(n,1)>.8)=NaN;
%  y = medfilt1m(x,2); [x0; x; y; x0-y]
%
% See also MODEFILT1, MEDFILT1
%
% Piotr's Image&Video Toolbox      Version 2.35
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

% apply medfilt1 (standard median filter) to valid locations in x
if(nargin<3 || isempty(z)), z=NaN; end; x=x(:)'; n=length(x);
if(isnan(z)), valid=~isnan(x); else valid=x~=z; end; v=sum(valid);
if(v==0), y=repmat(z,1,n); return; end
if(v<2*r+1), y=repmat(median(x(valid)),1,n); return; end
y=medfilt1(x(valid),2*r+1);

% get radius R needed at each location j to span s=2r+1 valid values
% get start (a) and end (b) locations and map back to location in y
C=[0 cumsum(valid)]; s=2*r+1; R=find(C==s); R=R(1)-2; pos=zeros(1,n);
for j=1:n, R0=R;
  R=R0-1; a=max(1,j-R); b=min(n,j+R);
  if(C(b+1)-C(a)<s), R=R0; a=max(1,j-R); b=min(n,j+R);
    if(C(b+1)-C(a)<s), R=R0+1; a=max(1,j-R); b=min(n,j+R); end
  end
  pos(j)=(C(b+1)+C(a+1))/2;
end
y=y(floor(pos));

end

function y = medfilt1( x, s )
% standard median filter (copied from medfilt1.m)
n=length(x); r=floor(s/2); indr=(0:s-1)'; indc=1:n;
ind=indc(ones(1,s),1:n)+indr(:,ones(1,n));
x0=x(ones(r,1))*0; X=[x0'; x'; x0'];
X=reshape(X(ind),s,n); y=median(X,1);
end

% function y = medfilt1( x, s )
% % standard median filter (slow)
% % get unique values in x
% [vals,disc,inds]=unique(x); m=length(vals); n=length(x);
% if(m>256), warning('x takes on large number of diff vals'); end %#ok<WNTAG>
% % create quantized representation [H(i,j)==1 iff x(j)==vals(i)]
% H=zeros(m,n); H(sub2ind2([m,n],[inds; 1:n]'))=1;
% % create histogram [H(i,j) is count of x(j-r:j+r)==vals(i)]
% H=localSum(H,[0 s],'same');
% % compute median for each j and map inds back to original vals
% [disc,inds]=max(cumsum(H,1)>s/2,[],1); y=vals(inds);
% end

author: zouchuhang

repo: LayoutNet