text:
classdef ShFunc_StressNorm3 < ShFunWithElasticPdes

    properties (Access = private)
        adjointProb
        fAdjoint
        stressNorm
        stressNormSquare
        stressNorml2LpNorm
        isElementToOptimize
        integralSigmaP

        ngaus
        nelem
        sigma
        sigmaA
        sigma2
        sigmaP

        dsigHPdsigH2v
        PpdSigHP_DsigH2M
    end

    methods (Access = public)

        function obj = ShFunc_StressNorm3(cParams)
            cParams.filterParams.quadratureOrder = 'LINEAR';
            obj.init(cParams);
            fileName = cParams.femSettings.fileName;
            obj.createEquilibriumProblem(fileName);
            obj.createAdjointProblem(fileName);
            obj.createOrientationUpdater();
            obj.createElementsToOptimize();
            obj.updateHomogenizedMaterialProperties();
            obj.nelem = obj.physicalProblem.mesh.nelem;
            obj.ngaus = obj.physicalProblem.element.quadrature.ngaus;
        end

        function f = getPdesVariablesToPrint(obj)
            f{1} = obj.getPdeVariableToPrint(obj.physicalProblem);
            f{2} = obj.getPdeVariableToPrint(obj.adjointProb);
        end

       function fP = addPrintableVariables(obj)
            fP{1}.value = obj.getPdeVariableToPrint(obj.physicalProblem);
            fP{2}.value = obj.getPdeVariableToPrint(obj.adjointProb);
            fP{3}.value = sqrt(obj.sigma2');
            fP{4}.value = obj.designVariable.alpha;
            fP{5}.value = abs(obj.designVariable.alpha);
            fP = obj.addHomogVariables(fP);
       end

        function fP = createPrintVariables(obj)
            types = {'Elasticity','Elasticity','ScalarGauss'...
                        'VectorGauss','VectorGauss'};
            names = {'Primal','Adjoint','StressNormGauss',...
                        'AlphaGauss','AlphaAbsGauss'};
            fP = obj.obtainPrintVariables(types,names);
            fP = obj.addHomogPrintVariablesNames(fP);
        end

        function v = getVariablesToPlot(obj)
            v{1} = max(sqrt(obj.sigma2));
            sNorm2 = obj.computeClosedSigmaProduct(obj.sigma,obj.sigma);
            v{2} = max(sqrt(sNorm2));
        end

        function t = getTitlesToPlot(obj)
            t{1} = 'Max amplified stress';
            t{2} = 'Max stress';
        end

        function setPnorm(obj,pNorm)
            obj.target_parameters.stressNormExponent = pNorm;
        end

    end

    methods (Access = protected)

        function updateHomogenizedMaterialProperties(obj)
            obj.filterDesignVariable();
            obj.homogenizedVariablesComputer.computeCtensor(obj.regDesignVariable);
            pNorm = obj.obtainPnorm();
            obj.homogenizedVariablesComputer.computePtensor(obj.regDesignVariable,pNorm);
        end

        function computeFunctionValue(obj)
            p  = obj.obtainPnorm();
            Pp = obj.homogenizedVariablesComputer.Pp;
            obj.sigma  = obj.computeSigmaTensor();
            obj.sigmaA = obj.amplifySigma(Pp,obj.sigma);
            obj.sigma2 = obj.computeClosedSigmaProduct(obj.sigmaA,obj.sigmaA);
            obj.sigmaP = obj.computeSigmaP(obj.sigma2);
            obj.integralSigmaP = obj.integrateSigmaP();
            obj.value = obj.integralSigmaP.^(1/p);
        end

        function solveState(obj)
            obj.physicalProblem.setC(obj.homogenizedVariablesComputer.C);
            obj.physicalProblem.computeVariables();
        end

        function computeGradientValue(obj)
            dJdx = obj.computeDJdx();
            dHdx = obj.computeDHdx();
            obj.gradient = dJdx + dHdx;
        end

        function solveAdjoint(obj)
            obj.computeFadjoint();
            obj.adjointProb.setC(obj.homogenizedVariablesComputer.C);
            obj.adjointProb.computeVariablesWithBodyForces(obj.fAdjoint);
        end

    end

    methods (Access = private)

        function g = initG(obj)
            g = zeros(obj.nelem,obj.ngaus,obj.nVariables);
        end

        function euR = rotateStrain(obj,eu)
            rot = obj.homogenizedVariablesComputer.rotator;
            euR  = rot.rotateStrain(eu);
        end

        function sR = rotateStress(obj,s)
            rot = obj.homogenizedVariablesComputer.rotator;
            sR  = rot.rotateStress(s);
        end

        function ds2ds = computedSigma2dSigma(obj,ds)
            Pp    = obj.homogenizedVariablesComputer.Pp;
            s     = obj.sigma;
            sA    = obj.amplifySigma(Pp,s);
            dsA   = obj.amplifySigma(Pp,ds);
            sXds  = obj.computeClosedSigmaProduct(sA,dsA);
            dsXs  = obj.computeClosedSigmaProduct(dsA,sA);
            ds2ds = sXds + dsXs;
        end

        function ds2dP = computedSigma2dP(obj,dPp)
            s     = obj.sigma;
            Pp    = obj.homogenizedVariablesComputer.Pp;
            sA    = obj.amplifySigma(Pp,s);
            dsA   = obj.amplifySigma(dPp,s);
            sXds  = obj.computeClosedSigmaProduct(sA,dsA);
            dsXs  = obj.computeClosedSigmaProduct(dsA,sA);
            ds2dP = sXds + dsXs;
        end

        function s = computeSigmaTensor(obj)
            eu = obj.physicalProblem.variables.strain;
            C  = obj.homogenizedVariablesComputer.Cref;
            s  =  obj.computeSigma(C,eu);
        end

        function sA = amplifySigma(obj,Pp,s)
            sA = obj.computeSigma(Pp,s);
        end

        function computeSigma2(obj)
            s = obj.sigmaA;
            obj.sigma2  = obj.computeClosedSigmaProduct(s,s);
        end

        function sNorm = computeClosedSigmaProduct(obj,sa,sb)
            nStre = size(sa,2);
            factor = [1 1 2];
            sNorm = zeros(obj.ngaus,obj.nelem);
            for iStre = 1:nStre
                sAi = obj.squeezeParticular(sa(:,iStre,:),2);
                sBi = obj.squeezeParticular(sb(:,iStre,:),2);
                sNorm = sNorm + factor(iStre)*(sAi.*sBi);
            end
        end

        function intSigmaP = integrateSigmaP(obj)
            dvolum = obj.physicalProblem.geometry.dvolu';
            int    = obj.sigmaP.*dvolum;
            intSigmaP = int(:,obj.isElementToOptimize);
            intSigmaP = sum(intSigmaP(:));
        end

        function sP = computeSigmaP(obj,s2)
            p  = obj.obtainPnorm;
            sP = s2.^(p/2);
        end

        function computeFadjoint(obj)
            dSigmaPdu    = obj.computedSigmaPdu();
            dJdu         = obj.computedJdSigmaP(dSigmaPdu);
            obj.fAdjoint = obj.assambleVector(-dJdu);
        end

        function dJdx = computeDJdx(obj)
            dSigmaPdx = obj.computedSigmaPdx();
            dJdx      = obj.computedJdSigmaP(dSigmaPdx);
            dJdx(~obj.isElementToOptimize,:) = 0;
        end

        function s = computeSigma(obj,C,eu)
            %eu = obj.rotateStrain(eu);
            s  = obj.computeStress(C,eu);
            %s = obj.rotateStress(s);
        end

        function dSigmaPdx = computedSigmaPdx(obj)
            eu = obj.physicalProblem.variables.strain;
            dC = obj.homogenizedVariablesComputer.dCref;
            dP = obj.homogenizedVariablesComputer.dPp;
            nX = size(dC,3);
            dSigmaPdx = zeros(nX,obj.ngaus,obj.nelem);
            for ix = 1:nX
                dCdx       = squeeze(dC(:,:,ix,:));
                dPdx       = squeeze(dP(:,:,ix,:));
                dSigmadx   = obj.computeSigma(dCdx,eu);
                dSigma2dxA = obj.computedSigma2dSigma(dSigmadx);
                dSigma2dxB = obj.computedSigma2dP(dPdx);
                dSigma2dx  = dSigma2dxA + dSigma2dxB;
                dSigmaPdx(ix,:,:) = obj.computedSigmaPdSigma2(dSigma2dx);
            end
            dSigmaPdx = permute(dSigmaPdx,[3 2 1]);
        end

        function dSigmaPdu = computedSigmaPdu(obj)
            C    = obj.homogenizedVariablesComputer.Cref;
            dEps = obj.computedEps();
            nV   = size(dEps,1);
            dSigmaPdu = zeros(nV,obj.ngaus,obj.nelem);
            for iv = 1:nV
                dEpsDu    = obj.squeezeParticular(dEps(iv,:,:,:),1);
                dSigmadu  = obj.computeSigma(C,dEpsDu);
                dSigma2du = obj.computedSigma2dSigma(dSigmadu);
                dSigmaPdu(iv,:,:) = obj.computedSigmaPdSigma2(dSigma2du);
            end
        end

        function dSigmaPdx = computedSigmaPdSigma2(obj,dSigma2dx)
            p = obj.obtainPnorm();
            dSigmaPdSigma2 = p/2*(obj.sigma2).^(p/2-1);
            dSigmaPdx = dSigmaPdSigma2.*dSigma2dx;
        end

        function s = computeStress(obj,C,e)
            s = zeros(size(e));
            for igaus = 1:obj.ngaus
                eG(1,:,:) = squeeze(e(igaus,:,:));
                sGt = bsxfun(@times,C,eG);
                s(igaus,:,:) = obj.squeezeParticular(sum(sGt,2),2);
            end
        end

        function g = computeDHdx(obj)
            eu = obj.physicalProblem.variables.strain;
            %eu = obj.rotateStrain(eu);
            dC = obj.homogenizedVariablesComputer.dCref;
            ep = obj.adjointProb.variables.strain;
            g = obj.initG();
            for ivar = 1:obj.nVariables
                dCiv = squeeze(dC(:,:,ivar,:));
                ds = obj.computeSigma(dCiv,ep);
                g(:,:,ivar) = obj.squeezeParticular(sum(eu.*ds,2),2)';
            end
        end

        function p = obtainPnorm(obj)
            p = obj.target_parameters.stressNormExponent;
        end

        function createAdjointProblem(obj,fileName)
            obj.adjointProb = FEM.create(fileName);
        end

        function dEps = computedEps(obj)
            dvolum = obj.physicalProblem.geometry.dvolu';
            nstre = obj.physicalProblem.element.getNstre();
            nnode = obj.physicalProblem.element.nnode;
            nunkn = obj.physicalProblem.element.dof.nunkn;
            nv    = nnode*nunkn;
            dEps = zeros(nv,obj.ngaus,nstre,obj.nelem);
            for igaus = 1:obj.ngaus
               B  = obj.physicalProblem.element.computeB(igaus);
               Bm = permute(B,[2 1 3]);
               dvG(1,1,:) = squeeze(dvolum(igaus,:));
               dvGm = repmat(dvG,nv,nstre,1);
               dEps(:,igaus,:,:) = Bm.*dvGm;
            end
        end

        function fAdjoint = assambleVector(obj,dsigmaPdu)
            eforce = dsigmaPdu;
            eforce(:,:,~obj.isElementToOptimize) = 0;
            Fvol = obj.physicalProblem.element.AssembleVector({eforce});
            fAdjoint = Fvol;
        end

        function dJdx = computedJdSigmaP(obj,dSigmaPdx)
            p = obj.obtainPnorm();
            intSigmaP = obj.integralSigmaP;
            dJ_dSigmaP = 1/p*intSigmaP^(1/p-1);
            dJdx = dJ_dSigmaP*dSigmaPdx;
        end

        function createElementsToOptimize(obj)
             phy = obj.getPhysicalProblems();
             m = phy{1}.mesh;
             xG = transpose(m.computeBaricenter);
             x = xG(:,1);
             y = xG(:,2);
%             isDirichletPartX = x > -1e-12 & x < 0.05;
%             isDirichletPartY = y > 0.20 & y < 0.80;
%             isDirichletPart = isDirichletPartX & isDirichletPartY;
%             isNeumannPartX = x > (1-0.05) & x < (1+1e-12);
%             isNeumannPartY = y > 0.40 & y < 0.60;
%             isNeumannPart = isNeumannPartX & isNeumannPartY;
%             isForOptimizing = ~isDirichletPart & ~isNeumannPart;
%             obj.isElementToOptimize = isForOptimizing;

%             s.connec = m.connec(isForOptimizing,:);
%             s.coord  = m.coord;
%             m2 = Mesh(s);
%             m2.plot();
           obj.isElementToOptimize = true(size(x));
        end

    end

    methods (Access = private, Static)

        function squeezed = squeezeParticular(t,dim)
            z = size(t);
            index = setdiff(1:length(z),dim);
            squeezed = reshape(t,[z(index) z(dim)]);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function nearest_interp_1d_test ( )

%*****************************************************************************80
%
%% NEAREST_INTERP_1D_TEST tests the NEAREST_INTERP_1D library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2012
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../r8lib' )
  addpath ( '../test_interp' )

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NEAREST_INTERP_1D_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the NEAREST_INTERP_1D library.\n' );
  fprintf ( 1, '  The R8LIB library is needed.\n' );
  fprintf ( 1, '  The test needs the TEST_INTERP library.\n' );

  prob_num = p00_prob_num ( );

  ni = 11;
  for prob = 1 : prob_num
    nearest_interp_1d_test01 ( prob, ni );
  end

  for prob = 1 : prob_num
    nearest_interp_1d_test02 ( prob );
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NEAREST_INTERP_1D_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../r8lib' )
  rmpath ( '../test_interp' )

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Psi, keepIDs] = relabelFeatIDs( Psi, keepIDs )
% Obtain Left-Ordered Form ordering on the features
%   and remove all features not allocated to any time series obj.
% As a result, F is reindexed by swapping columns
%    so it is in Left-Ordered Form and any columns of all zeros are dropped
% Relabelling must also apply to HMM transition and emission params
%    and also to the HMM hidden state sequence and some cached suff. stats
% These are just *relabelled*... their value (in log prob.) does not change

origF = Psi.F;
Psi.F = origF(:, keepIDs);
if isfield( Psi,'cache') && isfield( Psi.cache, 'logSoftEv')
    for ii = 1:size( Psi.F, 1 )
        [Kmax,Tii] = size( Psi.cache.logSoftEv{ii} );
        if max( keepIDs ) <= Kmax
            Psi.cache.logSoftEv{ii} = Psi.cache.logSoftEv{ii}( keepIDs, : );
        else
            kIDs = keepIDs( keepIDs <= Kmax );
            SoftEvTmp = -inf( length(keepIDs), Tii );
            SoftEvTmp( keepIDs <= Kmax, : ) = Psi.cache.logSoftEv{ii}( kIDs, : );
            Psi.cache.logSoftEv{ii} = SoftEvTmp;
        end
    end
end
if isfield( Psi, 'TransM' )
Psi.TransM = Psi.TransM.reallocateFeatIDs( origF, keepIDs );
end
if isfield( Psi, 'ThetaM' )
Psi.ThetaM = Psi.ThetaM.reallocateFeatIDs( keepIDs );
end

if isfield( Psi, 'stateSeq' )
    stateSeq = Psi.stateSeq;
    for ii = 1:size(Psi.F,1)
        for jj = 1:length( keepIDs )
            stateSeq(ii).z( Psi.stateSeq(ii).z == keepIDs(jj) ) = jj;
        end
    end
    Psi.stateSeq = stateSeq;
end

if isfield( Psi, 'activeFeatIDs' )
   for kk = 1:length( Psi.activeFeatIDs )
        Psi.activeFeatIDs(kk) = find( keepIDs == Psi.activeFeatIDs(kk)  );
   end
end

author: michaelchughes

repo: NPBayesHMM


text:
function fgsegment = dtiFiberExcludeDiscontinuous(fg, threshold)

% Remove the fascicles which has discontinuity by VOI-based fiber clipping
%
% Sometimes we create fascicles/streamlines which is discontinuous when we
% clip the fascicles using VOI/ROI (e.g. feConnectomeClip.m in LiFE). This function
% excludes the discontinuous fascicles from the connectome.
%
% INPUT:
% fg: fg structure for input
% threshold: If the nodal spacing exceeds this threshold, the fascicle will
%           be removed
%
% (C) Hiromasa Takemura, 2015 CiNet HHS/Stanford Vista Team

% The matrix to define fascicles to exclude
fascicletoremove = ones(length(fg.fibers),1);

for i=1:length(fg.fibers)

    x = fg.fibers{i}(1,:);
    y = fg.fibers{i}(2,:);
    z = fg.fibers{i}(3,:);

    try
        % Compute the 3d distance on each node
        for k=1:(length(x)-1)
            nodalspace(k) = sqrt((x(k+1)-x(k))^2 + (y(k+1)-y(k))^2 + (z(k+1)-z(k))^2);

            % Remove fascicles when the nodal spacing exceeds the threshold
            if nodalspace(k) > threshold
                fascicletoremove(i) = 0;
            end

        end

    catch
        % Remove fascicles if the code could not execute the process above
        % This happens when fascicle only contains one node as a consequence of the
        % clipping
        fascicletoremove(i) = 0;

    end
    clear nodalspace x y z k

end

% Remove fascicles having longer nodal spacing
fgsegment = fgExtract(fg, logical(fascicletoremove), 'keep');

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%Maps any angle to the equivalent between -pi and pi
function a=modAngle(a)
a=mod(a+pi,2*pi)-pi;

author: NicolasBoumal

repo: manopt


text:
function S_scalar2d = poynting(varargin)

chkarg(nargin == 5 || nargin == 7, 'five or seven arguments are required.')

iarg = 0;
iarg = iarg + 1; polarization = varargin{iarg};
chkarg(istypesizeof(polarization, 'Axis'), '"argument %d should be "polarization" (instance of Axis).', iarg);

if istypesizeof(varargin{iarg+1}, 'Scalar3d')
	iarg = iarg + 1; Ep3d = varargin{iarg};
	chkarg(istypesizeof(Ep3d, 'Scalar3d'), '"argument %d should be "Ep3d" (instance of Scalar3d)."', iarg);

	iarg = iarg + 1; Eq3d = varargin{iarg};
	chkarg(istypesizeof(Eq3d, 'Scalar3d'), '"argument %d should be "Eq3d" (instance of Scalar3d)."', iarg);

	iarg = iarg + 1; Hp3d = varargin{iarg};
	chkarg(istypesizeof(Hp3d, 'Scalar3d'), '"argument %d should be "Hp3d" (instance of Scalar3d)."', iarg);

	iarg = iarg + 1; Hq3d = varargin{iarg};
	chkarg(istypesizeof(Hq3d, 'Scalar3d'), '"argument %d should be "Hq3d" (instance of Scalar3d)."', iarg);

	iarg = iarg + 1; normal_axis = varargin{iarg};
	chkarg(istypesizeof(normal_axis, 'Axis'), 'argument %d should be "normal_axis" (instance of Axis).', iarg);

	iarg = iarg + 1; intercept = varargin{iarg};
	chkarg(istypesizeof(intercept, 'real'), 'argument %d should be "intercept" (real).', iarg);

	chkarg(isequal(Ep3d.grid3d, Eq3d.grid3d, Hp3d.grid3d, Hq3d.grid3d), ...
		'instances of Scalar3d do not have same grid3d.');

	% This makes the Poynting vector calculation independnet of whether the
	% primary grid is the E-field grid or H-field grid, when normal_axis is
	% normal to p and q.
	Ep2d = slice_scalar3d(Ep3d, normal_axis, intercept);
	Eq2d = slice_scalar3d(Eq3d, normal_axis, intercept);
	Hp2d = slice_scalar3d(Hp3d, normal_axis, intercept);
	Hq2d = slice_scalar3d(Hq3d, normal_axis, intercept);

	grid2d = Ep2d.grid2d;
else
	iarg = iarg + 1; Ep2d = varargin{iarg};
	chkarg(istypesizeof(Ep2d, 'Scalar2d'), '"argument %d should be "Ep2d" (instance of Scalar2d)."', iarg);

	iarg = iarg + 1; Eq2d = varargin{iarg};
	chkarg(istypesizeof(Eq2d, 'Scalar2d'), '"argument %d should be "Eq2d" (instance of Scalar2d)."', iarg);

	iarg = iarg + 1; Hp2d = varargin{iarg};
	chkarg(istypesizeof(Hp2d, 'Scalar2d'), '"argument %d should be "Hp2d" (instance of Scalar2d)."', iarg);

	iarg = iarg + 1; Hq2d = varargin{iarg};
	chkarg(istypesizeof(Hq2d, 'Scalar2d'), '"argument %d should be "Hq2d" (instance of Scalar2d)."', iarg);

	chkarg(isequal(Ep2d.grid2d, Eq2d.grid2d, Hp2d.grid2d, Hq2d.grid2d), ...
		'instances of Scalar2d do not have same grid2d.');  % Grid2d has normal_axis, so passing this test implies shared normal axis
	chkarg(isequal(Ep2d.intercept, Eq2d.intercept, Hp2d.intercept, Hq2d.intercept), ...
		'instances of Scalar2d do not have same intercept.');

	grid2d = Ep2d.grid2d;
	normal_axis = grid2d.normal_axis;
	intercept = Ep2d.intercept;
end

if polarization == normal_axis
	assert(all(Ep2d.gt_array==Hq2d.gt_array) && all(Eq2d.gt_array==Hp2d.gt_array));

	pi = grid2d.l{Dir.h,GT.dual};
	qi = grid2d.l{Dir.v,GT.dual};
	[PI, QI] = ndgrid(pi, qi);  % centers of grid cell faces

	[ep, l1] = Ep2d.data_expanded();
	hq = Hq2d.data_expanded();  % hq and ep have same location
	[eq, l2] = Eq2d.data_expanded();
	hp = Hp2d.data_expanded();  % hp and eq have same location

	sr1 = ep .* conj(hq);
	sr2 = eq .* conj(hp);

	[P1, Q1] = ndgrid(l1{:});  % locations of sr1
	[P2, Q2] = ndgrid(l2{:});  % locations of sr2

	% Interpolate sr1 and sr2 at the centers of grid cell faces.
	sr1 = interpn(P1, Q1, sr1, PI, QI);
	sr2 = interpn(P2, Q2, sr2, PI, QI);

	array = real(sr1 - sr2) / 2;

	% Attach extra points.
	for d = Dir.elems
		array = attach_extra_S(array, d, grid2d);
	end

	osc = Ep2d.osc;
	physQ = PhysQ.S;
	gt_array = [GT.dual, GT.dual];  % face centers
else  % polarization ~= normal_axis
	% To be implemented.  See an attempted usage in example/2d/pc_2d_basic.m.
	% Note that in this case, Ep and Hq are not co-located.  Neither are Eq and
	% Hp.  (Draw their locations, and compare it with Ew2d.gt_array and
	% Hw2d.gt_array.)
end

% The attached values to array should be the same as the ones inside the array
% if BC is not periodic.
S_scalar2d = Scalar2d(array, grid2d, gt_array, osc, physQ, ['<', physQ.symbol, '_', char(polarization), '>'], intercept);

function array = attach_extra_S(array, d, grid2d)
ind_n = {':', ':'};
ind_p = {':', ':'};
bc_d = grid2d.bc(d);
if bc_d == BC.p
	ind_n{d} = grid2d.N(d);
	ind_p{d} = 1;
else  % bc_d == BC.e or BC.m
	ind_n{d} = 1;
	ind_p{d} = grid2d.N(d);
end
array = cat(int(d), array(ind_n{:}), array, array(ind_p{:}));  % Bloch phases in S are ignored due to conj(H)

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function sr=comp_gabreassign(s,tgrad,fgrad,a);
%COMP_GABREASSIGN  Reassign time-frequency distribution.
%   Usage:  sr = comp_gabreassign(s,tgrad,fgrad,a);
%
%   `comp_gabreassign(s,tgrad,fgrad,a)` will reassign the values of the positive
%   time-frequency distribution *s* using the instantaneous time and frequency
%   *fgrad* and *ifdummy*. The lattice is determined by the time shift *a* and
%   the number of channels deduced from the size of *s*.
%
%   See also: gabreassign
%
%   References: aufl95

%   AUTHOR : Peter L. Søndergaard.
%   TESTING: OK
%   REFERENCE: OK

[M,N,W]=size(s);
L=N*a;
b=L/M;

freqpos=fftindex(M);
tgrad=bsxfun(@plus,tgrad/b,freqpos);

timepos=fftindex(N);
fgrad=bsxfun(@plus,fgrad/a,timepos.');

tgrad=round(tgrad);
fgrad=round(fgrad);

tgrad=mod(tgrad,M);
fgrad=mod(fgrad,N);

sr=zeros(M,N,W,assert_classname(s,tgrad,fgrad));

fgrad=fgrad+1;
tgrad=tgrad+1;

for w=1:W
    for ii=1:M
        for jj=1:N
            sr(tgrad(ii,jj),fgrad(ii,jj),w) = sr(tgrad(ii,jj),fgrad(ii,jj),w)+s(ii,jj,w);
        end;
    end;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%%
%
function show_2dresults(m,n,varargin)
  k = nargin;
  vars = varargin;
  if(k < 3)
    error('Not supported operation!');
  end
  for i = 1:size(vars{1},2)
    for j = 1:k-2
      T = vars{j};
      T = reshape(T(:,i),m,n);
      switch(k)
        case 1
          imshow(T,[],'InitialMagnification','fit');
        case {2,3}
          subplot(1,k,j), imshow(T,[],'InitialMagnification','fit');
        case 4
          subplot(2,2,j), imshow(T,[],'InitialMagnification','fit');
        case {5,6}
          subplot(2,3,j), imshow(T,[],'InitialMagnification','fit');
        case {8}
          subplot(2,4,j), imshow(T,[],'InitialMagnification','fit');
        case {7,9}
          subplot(3,3,j), imshow(T,[],'InitialMagnification','fit');
        otherwise
          error('Not supported operation!');
      end
    end
    pause(0.01);
  end
end

author: andrewssobral

repo: lrslibrary


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function p = pnig(xnew, xold, t, r, alpha, beta, delta)
% probability distribution of the NIG model
d = delta * t;                          % time dependent parameter

gamma = sqrt(alpha^2 - beta^2);

% martingale correction removes the drift!
omega = delta * gamma - delta*sqrt(alpha^2-(1+beta)^2); % martingale correction

x = xnew -xold - (r - omega) * t;           %
arg = sqrt(1+(x./delta).^2);

p = real(alpha / pi * besselk(1,alpha .* d .* arg)./arg.* exp(d.*(gamma+beta.*(x./d))));
end

author: Sable

repo: mcbench-benchmarks


text:
function p = preprocess_bilinear_bounds(p)

if ~isempty(p.integer_variables)
    for i = 1:size(p.bilinears,1)
        if ismember(p.bilinears(i,2),p.integer_variables)
            if ismember(p.bilinears(i,3),p.integer_variables)
                p.integer_variables = [p.integer_variables p.bilinears(i,1)];
            end
        end
    end
    if p.K.f > 0
        for i = 1:p.K.f
            if all(ismember(p.F_struc(i,:),[0 1 -1]))
                involved = find(p.F_struc(i,2:end));
                % One variable is linear combination of  integer variables
                if (nnz(ismember(involved,p.integer_variables)) == length(involved)-1) & length(involved)>1
                    p.integer_variables = [p.integer_variables involved];
                end
            end
        end
        p.integer_variables  = unique(p.integer_variables );
    end
end

if ~isempty(p.binary_variables)
    for i = 1:size(p.bilinears,1)
        if ismember(p.bilinears(i,2),p.binary_variables)
            if ismember(p.bilinears(i,3),p.binary_variables)
                p.binary_variables = [p.binary_variables p.bilinears(i,1)];
            end
        end
    end
    for i = 1:p.K.f
        if all(p.F_struc(i,:) == fix(p.F_struc(i,:)))
            involved = find(p.F_struc(i,2:end));
            % One variable is linear combination of binary variables
            if (nnz(ismember(involved,p.binary_variables)) == length(involved)-1) & length(involved)>1
                p.integer_variables = [p.integer_variables involved];
            end
        end
    end
    p.binary_variables  = unique(p.binary_variables );

end

if isempty(p.ub)
    p.ub = repmat(inf,length(p.c),1);
end
if isempty(p.lb)
    p.lb = repmat(-inf,length(p.c),1);
end
if ~isempty(p.F_struc)
    [lb,ub,used_rows_eq,used_rows_lp] = findulb(p.F_struc,p.K);
    if ~isempty([used_rows_eq;used_rows_lp])
        lower_defined = find(~isinf(lb));
        if ~isempty(lower_defined)
            p.lb(lower_defined) = max(p.lb(lower_defined),lb(lower_defined));
        end
        upper_defined = find(~isinf(ub));
        if ~isempty(upper_defined)
            p.ub(upper_defined) = min(p.ub(upper_defined),ub(upper_defined));
        end
        % Remove linear bound inequalities
        if ~isempty(used_rows_lp)
            used_rows_lp = used_rows_lp(find(~any(p.F_struc(p.K.f+used_rows_lp,1+p.nonlinears),2)));
            not_used_rows = setdiff(1:p.K.l,used_rows_lp);
             newKCutl = [];
            for i = 1:length(p.KCut.l)
                 newKCutl  = [newKCutl  find(not_used_rows==p.KCut.l(i))];
               % p.KCut.l(i) = find(not_used_rows == p.KCut.l(i));
               % p.originalModel.KCut.l(i) = find(not_used_rows == p.originalModel.KCut.l(i) );
            end
             p.KCut.l = newKCutl;
            if ~isempty(used_rows_lp)
                p.F_struc(p.K.f+used_rows_lp,:)=[];
              %  p.originalModel.F_struc(p.originalModel.K.f+used_rows_lp,:)=[];
                p.K.l = p.K.l - length(used_rows_lp);
               % p.originalModel.K.l = p.originalModel.K.l - length(used_rows_lp);
            end
        end
        % Remove linear bound inequalities
        if ~isempty(used_rows_eq)
            used_rows_eq = used_rows_eq(find(~any(p.F_struc(used_rows_eq,1+p.nonlinears),2)));
            not_used_rows = setdiff(1:p.K.f,used_rows_eq);
            newKCutf = [];
            for i = 1:length(p.KCut.f)
                newKCutf  = [newKCutf  find(not_used_rows==p.KCut.f(i))];
               % p.KCut.f(i) = find(not_used_rows==p.KCut.f(i));
              %  p.originalModel.KCut.f(i) = find(not_used_rows==p.originalModel.KCut.f(i));
            end
            p.KCut.f = newKCutf;
            if ~isempty(used_rows_eq)
                p.F_struc(used_rows_eq,:)=[];
              %  p.originalModel.F_struc(used_rows_eq,:)=[];
                p.K.f = p.K.f - length(used_rows_eq);
              %  p.originalModel.K.f = p.originalModel.K.f - length(used_rows_eq);
            end
        end
    end
end
p.lb(p.binary_variables) = max(0,p.lb(p.binary_variables));
p.ub(p.binary_variables) = min(1,p.ub(p.binary_variables));
p.lb(p.integer_variables) = ceil(p.lb(p.integer_variables));
p.ub(p.integer_variables) = floor(p.ub(p.integer_variables));
p = clean_bounds(p);
if ~isempty(p.bilinears)
    p = updatemonomialbounds(p);
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
 function y = tomo_filter(omega, ob)
%function y = tomo_filter(omega, ob)
% build the sinogram-spectrum-sized matrix
% that is .* multiplied after 2D FT, before iFFT or NUiFFT
% Copyright 2001-1, Jeff Fessler, The University of Michigan
% Extend to fan-beam geometry, 2003-11, Yingying Zhang

%
% effect of image-domain shift
%
y = exp(1i * (omega * ob.nxy_shift));

%
% effect of image basis function (extension to non-rect by S. Matej)
% basis: b(x/Dx,y/Dy) <-FT-> Dx*Dy * B(Dx*u,Dy*v)
%
y = y .* (ob.pixel_size).^2;

if isempty(ob.basis.type) | streq(ob.basis.type, 'pixel')
	% sinc_2 due to square pixels
	if ob.chat, printf('pixel basis, pixel size=%g',ob.pixel_size), end
	y = y .* sinc(omega(:,1)/(2*pi)) .* sinc(omega(:,2)/(2*pi));

elseif streq(ob.basis.type, 'no')
	if ob.chat, printf('no image basis modeled'), end

elseif streq(ob.basis.type, 'KB')
	if ob.chat, printf('KB basis: J=%g, alpha=%g, m=%g, n=%g', ...
		ob.basis.diam, ob.basis.shape, ob.basis.m, ob.basis.dim); end
	norm = kaiser_bessel_ft(0, ...
		ob.basis.diam, ob.basis.shape, ob.basis.m, ob.basis.dim);
	nd = sqrt(omega(:,1)/(2*pi) .* omega(:,1)/(2*pi) + ...
		omega(:,2)/(2*pi) .* omega(:,2)/(2*pi));
	y = y .* kaiser_bessel_ft(nd, ...
		ob.basis.diam, ob.basis.shape, ob.basis.m, ob.basis.dim) / norm;

elseif streq(ob.basis.type, 'Gauss')
	error('Gaussian basis not yet implemented')

else
	error(sprintf('basis function %s not implemented', ob.basis.type))
end

% form into sinogram shape
K = ob.Krho;
y = reshape(y, K, ob.na);

%
% detector radial response resolution-loss effect in frequency domain
% extension to parallel-beam non-rect case by S. Matej
%

kk = [-K/2 : K/2-1]';
if isempty(ob.beam.type) | streq(ob.beam.type, 'rect')

	% Trick: in fan-beam, the detector blur is not shift-invariant.
	% We approximate by the blur at the center of object.
	% The effective strip_width / ray_spacing at the center is the same.

	strip_ray = ob.strip_width / ob.ray_spacing;
	blur = sinc(strip_ray * kk / K);
	if ob.chat
		printf('strip integrals, relative beam width=%g', strip_ray)
	end

elseif streq(ob.beam.type, 'line')
	blur = ones(size(kk));
	if ob.chat, printf('line integrals modeled'), end

elseif streq(ob.beam.type, 'KB')
	if ob.chat, printf('KB beam shape: J=%g, alpha=%g, m=%g', ...
		ob.beam.diam, ob.beam.shape, ob.beam.m), end
	rel_bdiam = ob.beam.diam;	% KB_diameter relative to ob.ray_spacing
	norm = kaiser_bessel_ft(0, rel_bdiam, ob.beam.shape, ob.beam.m, 1);
	blur = kaiser_bessel_ft(...
		kk/K, rel_bdiam, ob.beam.shape, ob.beam.m, 1) / norm;

elseif streq(ob.beam.type, 'Gauss')
	g_sigma = ob.beam.shape / sqrt(8*log(2));
	blur = exp(-(2*pi*kk/K).^2 .* g_sigma^2/2);

else
	error(sprintf('beam shape %s not implemented', ob.beam.type))

end

y = y .* repmat(blur, [1 ob.na]);	% include blur effect

%
% phase "shift" to effect a half-pixel shift in each row
% corresponding to "offset=0" in tomographic projection.
% build in the post-fft shift too while at it.
%
if streq(ob.geometry, 'par') % do we need this in fan-beam ???
	phase = exp(1i*2*pi*(K+ob.is.shift0)/2 * kk / K);
	y = y .* repmat(phase, [1 ob.na]);
	y = y ./ ob.ray_spacing;	% see JF tech. report
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function totalRCS = RCS_Compute(complex_data, oversample_ratio, cal_sf, mask)
%RCS_COMPUTE Demonstrates area-based computation of calibrated RCS for an ROI
%
% USAGE:
%   totalRCS = RCS_Compute(complex_data, oversample_ratio, cal_sf, mask)
%
% INPUTS:
%   complex_data     - required : Complex valued SAR dataset in the image
%                                 domain.  First dimension is azimuth,
%                                 second range.  Third dimension could be
%                                 for multi-channel (i.e. polarimetric)
%                                 data.  This code assumes the frequency
%                                 support is centered and constant across
%                                 the image.
%   oversample_ratio - optional : Overample or zeropad factor. (Required
%                                 for calibrated RCS.)
%   cal_sf           - optional : Calibration scale factor (linear).
%                                 Either a constant or an array the same
%                                 size as complex_data with per-pixel
%                                 values. (Required for calibrated RCS.)
%   mask             - optional : Binary image which is ones over the
%                                 region of interest. (Default is an image
%                                 of all ones.)
%
% OUTPUTS:
%   totalRCS         - required : Total calibrated RCS of ROI (calibrated
%                                    if oversample_ratio and cal_sf were
%                                    given.)
%
% Author: Tim Cox, NRL; Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Default parameter values
if ~exist('oversample_ratio','var')
    oversample_ratio = [1 1]; % RCS values will be uncalibrated if this is unknown
end
if ~exist('cal_sf','var')
    cal_sf = 1; % RCS values will be uncalibrated if this is unknown
end
if ~exist('mask','var')
    mask = ones(size(complex_data));
end

% Apply shape mask to rectangular data
filtimg = complex_data.*repmat(mask,[1 1 size(complex_data,3)]);
if ~isscalar(cal_sf)
    % Only works if the radiometric scale factors are same for both channels
    cal_sf = repmat(cal_sf,[1 1 size(complex_data,3)]);
end

% Image domain computation of total RCS (Reference Adam Bryant, NGA/PL)
totalRCS = (1/prod(oversample_ratio)) * sum(sum(cal_sf.*abs(filtimg).^2));
totalRCS = squeeze(totalRCS); % Multi-channel (polarimetric) data
% The oversample ratio (or zeropad) factor is the ratio between the sum of
% squared (power detected) samples of an ideal sinc function (which is what
% we are measuring, at least for an ROI which we assume contains nearly all
% the energy of that sinc) and the peak of that ideal sinc^2 function
% (which is how the calibration constant is defined).
% We put the scale factor inside the sum so that it can be applied per
% pixel, rather than only as a constant.

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
%% Code Division Multiple Access Transmitter
% CDMAt - Function
% 1. (s) Input the data
% Input data must be (+/- 1's) (this is a PSK modulation (BPSK))
% 2. (hl) Hadamard matrix length
% 3. (cn) code number to be used for this user (row - number of H - matrix)
% 4. Spread the data by multiplying s by cn.
% 6. Outpot of the function is spread symbol of user cn.
% Montadar Abas Taher
% 11/03/2011
function [outcdmat]=cdmat(s,hl,cn)
if cn>hl
    errordlg('The input code number must be equal or less than the Hadamard length','File Error');
end
%% Generate Hadamard Matrix of length (hl)
h=hadamard(hl);
%% Spread the input sequence
outcdmat=kron(s,h(cn,:));


text:
function [img, mmPerVoxNew, dimOrder, dimFlip] = applyCannonicalXform(img, img2std, mmPerVox, insertMarkerFlag)
% Applies the cannonical xform specified in img2std to img
%
%  [img, mmPerVoxNew, dimOrder, dimFlip] = ...
%   applyCannonicalXform(img, img2std, mmPerVox, insertMarkerFlag)
%
% Applies the cannonical xform specified in img2std. NOTE! This function
% does not do a full affine transform- it assumes that the transform
% specified is a set of simple rotations and flips.
%
% When used with computeCannonicalXformFromIfile, img is reoriented to a
% standard axial orientation. That is, the image will be reoriented so that
% left-right is along the x-axis, anterior-posterior is along the y-axis,
% superior-inferior is along the z-axis, and the leftmost, anterior-most,
% superior-most point is at 0,0,0 (which, for Analyze, is the lower
% left-hand corner of the last slice).
%
% To help you find the 0,0,0 point, a 4-pixel rectange is drawn there with a pixel
% value equal to the maximum image intensity. Set insertMarkerFlag to
% false to disbale this. Also, if ndims(img)>3, this is automatically
% diabled.
%
% To adjust an img-to-standard space xofrm, use:
%
%    newXform = inv(canXform*inv(imToScanXform));
%
% SEE ALSO: computeCannonicalXformFromIfile
%
% HISTORY:
%   2003.06.19 RFD (bob@white.stanford.edu): wrote it.
%%

if(nargin<2)
    help(mfilename);
    return;
end
if(~exist('insertMarkerFlag','var')||isempty(insertMarkerFlag))
    insertMarkerFlag = true;
end

% Extract rotation & scale matrix. We have set things up so that the scales
% should be 1.
img2stdRot = round(img2std(1:3,1:3));

% Note that we have constructed this transform matrix so that it will only
% involve cannonical rotations. We did this by specifying corresponding
% points from cannonical locations (the corners of the volume- see
% stdCoords and volCoords).

% We use shortcuts to apply the transform. Since all rotations are
% cannonical, we can achieve them efficiently by swapping dimensions with
% 'permute'. The dimension permutation logic- we want to know which of the
% current dimensions should be x, which should be y, and which should be z:
xdim = find(abs(img2stdRot(1,:))==1);
ydim = find(abs(img2stdRot(2,:))==1);
zdim = find(abs(img2stdRot(3,:))==1);
dimOrder = [xdim, ydim, zdim];
dimFlip = [0 0 0];
if exist('mmPerVox','var')
    mmPerVoxNew = [mmPerVox(xdim), mmPerVox(ydim), mmPerVox(zdim)];
else
    mmPerVoxNew = [];
end

% Allow >3 dims- just leave the extra dims alone.
img = permute(img, [dimOrder, 4, 5]);

% Now do any necessary mirror flips (indicated by negative rotation matrix values).
if img2stdRot(1,xdim)<0
    dimFlip(xdim) = 1;
    % flip each slice ud (ie. along matlab's first dimension- our x-axis)
    %for(jj=1:size(img,3)) img(:,:,jj) = flipud(squeeze(img(:,:,jj))); end
    img = flipdim(img, 1);
end
if img2stdRot(2,ydim)<0
    dimFlip(ydim) = 1;
    % flip each slice lr (ie. along matlab's second dimension- our y-axis)
    %for(jj=1:size(img,3)) img(:,:,jj) = fliplr(squeeze(img(:,:,jj))); end
    img = flipdim(img, 2);
end
if img2stdRot(3,zdim)<0
    dimFlip(zdim) = 1;
    % reorder slices (ie. flip along the 3rd dim)
    %for(jj=1:size(img,1)) img(jj,:,:) = fliplr(squeeze(img(jj,:,:))); end
    img = flipdim(img, 3);
end

if insertMarkerFlag&&ndims(img)==3
    % insert a marker at 1,end,end (should be left, anterior, superior
    % given stdCoords of [0,0,0; 0,-1,0; 0,0,-1; -1,0,0])
    img(1,end,end) = max(img(:));
    img(2,end,end) = img(end,1,end);
    img(1,end-1,end) = img(end,1,end);
    img(2,end-1,end) = img(end,1,end);
end

return;

%%


text:
% HLINE - Plot 2D lines defined in homogeneous coordinates.
%
% Function for ploting 2D homogeneous lines defined by 2 points
% or a line defined by a single homogeneous vector
%
% Usage:   hline(p1,p2)   where p1 and p2 are 2D homogeneous points.
%          hline(p1,p2,'colour_name')  'black' 'red' 'white' etc
%          hline(l)       where l is a line in homogeneous coordinates
%          hline(l,'colour_name')
%

%  Peter Kovesi
%  School of Computer Science & Software Engineering
%  The University of Western Australia
%  pk @ csse uwa edu au
%  http://www.csse.uwa.edu.au/~pk
%
%  April 2000

function hline(a,b,c)

col = 'blue';  % default colour

if nargin >= 2 & isa(a,'double')  & isa(b,'double')   % Two points specified

  p1 = a./a(3);        % make sure homogeneous points lie in z=1 plane
  p2 = b./b(3);

  if nargin == 3 & isa(c,'char')  % 2 points and a colour specified
    col = c;
  end

elseif nargin >= 1 & isa(a,'double')       % A single line specified

  a = a./a(3);   % ensure line in z = 1 plane (not needed??)

  if abs(a(1)) > abs(a(2))   % line is more vertical
    ylim = get(get(gcf,'CurrentAxes'),'Ylim');
    p1 = hcross(a, [0 1 0]');
    p2 = hcross(a, [0 -1/ylim(2) 1]');
  else                       % line more horizontal
    xlim = get(get(gcf,'CurrentAxes'),'Xlim');
    p1 = hcross(a, [1 0 0]');
    p2 = hcross(a, [-1/xlim(2) 0 1]');
  end

  if nargin == 2 & isa(b,'char') % 1 line vector and a colour specified
    col = b;
  end

else
  error('Bad arguments passed to hline');
end

line([p1(1) p2(1)], [p1(2) p2(2)], 'color', col);

author: jianxiongxiao

repo: ProfXkit


text:
function data_num = p07_data_num ( )

%*****************************************************************************80
%
%% P07_DATA_NUM returns the number of data points for problem p07.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer DATA_NUM, the number of data points.
%
  data_num = 4;

  return
end


text:
classdef m_22_vic_10p_3s < MARRMoT_model
% Class for hydrologic conceptual model:  VIC

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Model references
% Clark, M. P., Slater, A. G., Rupp, D. E., Woods, R. A., Vrugt, J. A.,
% Gupta, H. V., � Hay, L. E. (2008). Framework for Understanding Structural
% Errors (FUSE): A modular framework to diagnose differences between
% hydrological models. Water Resources Research, 44(12), W00B02.
% http://doi.org/10.1029/2007WR006735
%
% Liang, X., Lettenmaier, D. P., Wood, E. F., & Burges, S. J. (1994). A
% simple hydrologically based model of land surface water and energy fluxes
% for general circulation models. Journal of Geophysical Research, 99,
% 14415�14428.

    properties
        % model-specific attributes
        aux_theta     % auxiliary parameter set
    end
    methods

        % creator method
        function obj = m_22_vic_10p_3s()
            obj.numStores = 3;                                             % number of model stores
            obj.numFluxes = 11;                                            % number of model fluxes
            obj.numParams = 10;

            obj.JacobPattern  = [1,0,0;
                                 1,1,0;
                                 1,1,1];                                   % Jacobian matrix of model store ODEs

            obj.parRanges = [   0.1 , 5;        % ibar,     Mean interception capacity [mm]
                                0   , 1;        % idelta,   Seasonal interception change as fraction of mean [-]
                                1   , 365;      % ishift,   Maximum interception peak timing [-]
                                1   , 2000;     % stot,    Maximum soil moisture capacity [mm]
                                0.01, 0.99;     % fsm, Fraction of stot that constitutes maximum soil moisture smmax [-]
                                0   ,10;        % b,        Infiltration excess shape parameter [-]
                                0   , 1;        % k1,       Percolation time parameter [d-1]
                                0   ,10;        % c1,       Percolation non-linearity parameter [-]
                                0   ,1;         % k2,       Baseflow time parameter [d-1]
                                1   , 5];       % c2,       Baseflow non-linearity parameter

            obj.StoreNames = {"S1", "S2" "S3"};                             % Names for the stores
            obj.FluxNames  = {"ei",   "peff", "iex", "qie",  "inf", "et1",...
                              "qex1", "pc",   "et2", "qex2", "qb"};        % Names for the fluxes

            obj.FluxGroups.Ea = [1 6 9];                                   % Index or indices of fluxes to add to Actual ET
            obj.FluxGroups.Q  = [4 7 10 11];                               % Index or indices of fluxes to add to Streamflow

        end

        % INITialisation function
        function obj = init(obj)
            %Parameters
            theta = obj.theta;
            stot    = theta(4);     % Total available storage [mm]
            fsm     = theta(5);     % Fraction of stot that constitutes maximum soil mositure storage [-]

            % Auxiliary parameter
            smmax   = fsm*stot;     % Maximum soil moisture capacity [mm]
            gwmax   = (1-fsm)*stot; % Maximum groundwater storage [mm]
            tmax    = 365.25;       % Length of one growing cycle [d]
            obj.aux_theta = [smmax; gwmax; tmax];

        end

        % MODEL_FUN are the model governing equations in state-space formulation
        function [dS, fluxes] = model_fun(obj, S)
            % parameters
            theta   = obj.theta;
            ibar    = theta(1);     % Mean interception capacity [mm]
            idelta  = theta(2);     % Seasonal interception change as fraction of mean [-]
            ishift  = theta(3);     % Maximum interception peak timing [-]
            b       = theta(6);     % Infiltration excess shape parameter [-]
            k1      = theta(7);     % Percolation time parameter [d-1]
            c1      = theta(8);     % Percolation non-linearity parameter [-]
            k2      = theta(9);     % Baseflow time parameter [d-1]
            c2      = theta(10);    % Baseflow non-linearity parameter

            aux_theta = obj.aux_theta;
            smmax = aux_theta(1);
            gwmax = aux_theta(2);
            tmax = aux_theta(3);

            % delta_t
            delta_t = obj.delta_t;

            % stores
            S1 = S(1);
            S2 = S(2);
            S3 = S(3);

            % climate input
            t = obj.t;                             % this time step
            climate_in = obj.input_climate(t,:);   % climate at this step
            P  = climate_in(1);
            Ep = climate_in(2);
            T  = climate_in(3);

            % fluxes functions
            aux_imax  = phenology_2(ibar,idelta,ishift,obj.t,tmax,delta_t);
            flux_ei   = evap_7(S1,aux_imax,Ep,delta_t);
            flux_peff = interception_1(P,S1,aux_imax);
            flux_iex  = excess_1(S1,aux_imax,delta_t);
            flux_qie  = saturation_2(S2,smmax,b,flux_peff+flux_iex);
            flux_inf  = effective_1(flux_peff+flux_iex,flux_qie);
            flux_et1  = evap_7(S2,smmax,max(0,Ep-flux_ei),delta_t);
            flux_qex1 = saturation_1(flux_inf,S2,smmax);
            flux_pc   = percolation_5(k1,c1,S2,smmax,delta_t);
            flux_et2  = evap_7(S3,gwmax,max(0,Ep-flux_ei-flux_et1),delta_t);
            flux_qex2 = saturation_1(flux_pc,S3,gwmax);
            flux_qb   = baseflow_5(k2,c2,S3,gwmax,delta_t);

            % stores ODEs
            dS1 = P - flux_ei - flux_peff - flux_iex;
            dS2 = flux_inf - flux_et1 - flux_qex1 - flux_pc;
            dS3 = flux_pc  - flux_et2 - flux_qex2 - flux_qb;

            % outputs
            dS = [dS1 dS2 dS3];
            fluxes = [flux_ei  flux_peff flux_iex  flux_qie ...
                      flux_inf flux_et1  flux_qex1 flux_pc ...
                      flux_et2 flux_qex2 flux_qb];
        end

        % STEP runs at the end of every timestep
        function obj = step(obj)
        end
    end
end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function X = calculateHuberMean(Y, rho, iters)
% Perform a robust mean under the Huber loss function.
% x = calculateRobust(Y, rho, iters)
%
% Input:
%   Y : MxN matrix over which to average (columnwise)
%   rho : augmented Lagrangian variable (default: 1)
%   iters : number of iterations to perform (default: 1000)
%
% Output:
%   x : 1xN vector that is the roust mean of Y
%
% Based on the ADMM Matlab codes also found at:
%   http://www.stanford.edu/~boyd/papers/distr_opt_stat_learning_admm.html
%
% Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
% 2013-09-26

if ~exist('rho', 'var') || isempty(rho)
    rho = 1;
end
if ~exist('iters', 'var') || isempty(iters)
    iters = 1000;
end

m = size(Y,1);
if m==1
    X = Y;
else
    mu = sum(Y)/m;
    Z = zeros(size(Y)); U = Z;
    for k = 1:iters
        X = mu + sum(Z - U)/m;
        D = bsxfun(@minus, X, Y - U);
        Z = (rho/(1+rho) + (1/(1+rho))*max(0,(1-(1+1/rho)./abs(D)))).*D;
        U = D - Z;
    end
end

author: VisLab

repo: EEG-Clean-Tools


text:
function [normF, normLoc] = norm(f, n)
%NORM   Norm of a CHEBFUN object.
%   For scalar-valued column CHEBFUN objects:
%       NORM(F) = sqrt(integral of abs(F)^2).
%       NORM(F, 2) is the same as NORM(F).
%       NORM(F, 'fro') is also the same as NORM(F).
%       NORM(F, 1) = integral of abs(F).
%       NORM(F, P) = (integral of abs(F)^P)^(1/P).
%       NORM(F, inf) = max(abs(F)).
%       NORM(F, -inf) = min(abs(F)).
%
%   For array-valued column CHEBFUN objects:
%       NORM(F) is the Frobenius norm, sqrt(sum(svd(F).^2)).
%       NORM(F, 'fro') is the same as NORM(F).
%       NORM(F, 1) is the maximum of the 1-norms of the columns of F.
%       NORM(F, 2) is the largest singular value of F.
%       NORM(F, inf) is the maximum of the 1-norms of the rows of F.
%       NORM(F, -inf) is the minimum of the 1-norms of the rows of F.
%       NORM(F, P) is the P-th root of the maximum of the sum of the P-th
%                  powers of the magnitudes of the columns of F.
%
% Furthermore, the +\-inf norms for scalar-valued CHEBFUN objects may also
% return a second output, giving the position where the max/min occurs. For
% array-valued CHEBFUN objects, the 1 norm can return as its 2nd output the
% index of the column with the largest norm, while the inf and -inf norms
% can return as their 2nd output the point in the domain of the CHEBFUN at
% which the norm is attained.
%
% If F is a row CHEBFUN, NORM(F, TYPE) is equal to NORM(F.', TYPE).

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

% Empty CHEBFUN has norm 0:
if ( isempty(f) )
    normF = 0;
    return
end

if ( nargin == 1 )
    n = 'fro'; 	% Frobenius norm is the default (2 norm would be much slower).
end

% Initialise:
normLoc = [];
if ( numel(f) > 1 )
    numCols = numel(f);
else
    numCols = size(f.funs{1}, 2);
end

%% %%%%%%%%%%%%%%%%%%%%%%%%% SCALAR-VALUED CHEBFUNS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ( numCols == 1 )
    switch n
        case 1
            if ( nargout == 2 )
                error('CHEBFUN:CHEBFUN:norm:argout', ...
                        'Cannot return two outputs for 1-norms');
            end
            normF = sum(abs(f));

        case {2, 'fro'}
            if ( nargout == 2 )
                error('CHEBFUN:CHEBFUN:norm:argout', ...
                        'Cannot return two outputs for ''fro''-norms');
            end
            f.isTransposed = 0;
            normF = sqrt(abs(innerProduct(f, f)));

        case {inf, 'inf'}
            if ( isreal(f) )
                [normF, normLoc] = minandmax(f);
                [normF, idx] = max(abs(normF));
                normLoc = normLoc(idx);
            else
                [normF, normLoc] = max(conj(f).*f);
                normF = sqrt(normF);
            end

        case {-inf, '-inf'}
            [normF, normLoc] = min(conj(f).*f);
            normF = sqrt(normF);

        otherwise
            if ( isnumeric(n) && isreal(n) )
                if ( nargout == 2 )
                    error('CHEBFUN:CHEBFUN:norm:argout', ...
                            'Cannot return two outputs for p-norms.');
                end
                if ( mod(n, 2) == 0 )
                    normF = sum((conj(f).*f).^(n/2))^(1/n);
                else
                    normF = sum(abs(f).^n)^(1/n);
                end
            else
                error('CHEBFUN:CHEBFUN:norm:unknownNorm', ...
                 'The only matrix norms available are 1, 2, inf, and ''fro''.');
            end
    end

%% %%%%%%%%%%%%%%%%%%%%%%%%%% ARRAY-VALUED CHEBFUNS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
    if ( f(1).isTransposed )
        f = f.';
    end

    switch n
        case 1
            f = mat2cell(f);
            normF = zeros(1, numel(f));
            for k = 1:numel(f)
                normF(k) = norm(f{k}, 1);
            end
            [normF, normLoc] = max(normF);

        case 2
            if (nargout == 2 )
                error('CHEBFUN:CHEBFUN:norm:argout', ...
                    ['Cannot return two outputs for 2-norms of ' ...
                     'array-valued CHEBFUNs.']);
            end
            s = svd(f, 0);
            normF = s(1);

        case 'fro'
            if ( nargout == 2 )
                error('CHEBFUN:CHEBFUN:norm:argout', ...
                    ['Cannot return two outputs for ''fro''-norms of ' ...
                     'array-valued CHEBFUNs.']);
            end
            normF = sqrt(sum(sum(f.*conj(f))));

        case {'inf', inf}
            [normF, normLoc] = max(sum(abs(f), 2));

        case {'-inf', -inf}
            [normF, normLoc] = min(sum(abs(f), 2));

        otherwise
            if ( isnumeric(n) && isreal(n) )
                [normF, normLoc] = max(sum(abs(f).^n, 2));
                normF = normF^(1/n);
            else
                error('CHEBFUN:CHEBFUN:norm:unknownNorm', ...
                 'The only matrix norms available are 1, 2, inf, and ''fro''.');
            end

    end

end

% Discard possible imaginary rounding errors:
normF = real(normF);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ l, u ] = l1pp_lu ( n, h )

%*****************************************************************************80
%
%% L1PP_LU computes the LU factors of the 1D PP Laplacian.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%    N must be at least 3.
%
%    Input, real H, the spacing between points.
%
%    Output, real L(N,N), U(N,n), the LU factors.
%
  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'L1PP_LU - Fatal error!\n' );
    fprintf ( 1, '  N < 3.\n' );
    error ( 'L1PP_LU - Fatal error!' );
  end

  l = zeros ( n, n );

  for i = 1 : n
    l(i,i) = 1.0;
  end

  for i = 2 : n - 1
    l(i,i-1) = - ( i - 1 ) / i;
    l(n,i-1) =       - 1   / i;
  end
  l(n,n-1) = - 1.0;

  u = zeros ( n, n );

  for i = 1 : n - 2
    u(i,i)   = ( i + 1 ) / i;
    u(i,i+1) = - 1.0;
    u(i,n) =     - 1   / i;
  end

  i = n - 1;
  u(i,i)   =   ( i + 1 ) / i;
  u(i,i+1) = - ( i + 1 ) / i;

  i = n;
  u(i,i) = 0.0;

  u(1:n,1:n) = u(1:n,1:n) / h / h;

  return
end


text:
% 自由增长模型演示程序

ts = 0 : 20;   % 时间天数
lambda = 0.3;  % 每个病人每天感染人数
x0 = 1;   % 初始病人数
infective = x0 * exp(lambda * ts);   % 病人数指数增长
plot(ts, infective);
title('自由增长模型');
grid;

repo: NumericalAnalysis


text:
function [f, g] = hsvargplvmObjectiveGradient(params, model)

% SVARGPLVMOBJECTIVEGRADIENT Wrapper function for shared VARGPLVM objective and gradient.
% FORMAT
% DESC returns the negative log likelihood of a Gaussian process
% model given the model structure and a vector of parameters. This
% allows the use of NETLAB minimisation functions to find the model
% parameters.
% ARG params : the parameters of the model for which the objective
% will be evaluated.
% ARG model : the model structure for which the objective will be
% evaluated.
% RETURN f : the negative log likelihood of the SVARGPLVM model.
% RETURN g : the gradient of the negative log likelihood of the SVARGPLVM
% model with respect to the parameters.
%
% SEEALSO : minimize, svargplvmModelCreate, svargplvmGradient, svargplvmLogLikelihood, svargplvmOptimise
%
% COPYRIGHT :  Andreas C. Damianou, 2011

% SVARGPLVM

% Check how the optimiser has given the parameters
if size(params, 1) > size(params, 2)
  % As a column vector ... transpose everything.
  transpose = true;
  model = hsvargplvmExpandParam(model, params');
else
  transpose = false;
  model = hsvargplvmExpandParam(model, params);
end

f = - hsvargplvmLogLikelihood(model);
% fprintf(1,'# F: %.13f\n',f); %%% DEBUG
if nargout > 1
  g = - hsvargplvmLogLikeGradients(model);
%  fprintf(1,'# G: %.13f .\n',sum(abs(g))); %%% DEBUG
end
if transpose
  g = g';
end

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
function [f, mergedPts] = merge(f, index, pref)
%MERGE   Remove unnecessary breakpoints in from a CHEBFUN.
%   F = MERGE(F, PREF) removes unnecessary breakpoints from a CHEBFUN F. In
%   particular the kth breakpoint is removed if the resulting FUN on the
%   interval [x_{k-1}, x_{k+1}] can be represented to the same accuracy as F
%   with a fewer than PREF.MAXLENGTH points when PREF.SPLITTING = 0 and
%   PREF.SPLITPREFS.SPLITLENGTH points when PREF.SPLITTING = 1. If a PREF is
%   not passed, then the default CHEBFUN.PREF() is used.
%
%   [F, MERGEDPTS] = MERGE(F) returns the index of the merged endpoints in the
%   vector MERGEDPTS.
%
%   MERGE(F, INDEX) or MERGE(F, INDEX, PREF) attempts to eliminate the endpoints
%   specified in INDEX. MERGE(F, 'all') is equivalent to MERGE(F,
%   [2:length(F.domain)-1]). (Note that it doesn't make sense to consider
%   merging the first and final breakpoints.)
%
%   In all cases, elimination is attempted from left to right, and non-trivial
%   pointValues will prevent merging at the corresponding breakpoints.
%
%   Example:
%       f = chebfun(@(x) abs(x), 'splitting','on');
%       [g, mergedPts] = merge(f.^2);
%
% See also CHEBFUNPREF.

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

if ( numel(f) > 1 )
    % TODO:  Implement this.
    error('CHEBFUN:CHEBFUN:merge:quasi', ...
        'MERGE does not support quasimatrices.');
end

% Convert to a column CHEBFUN so that feval(f, x) returns a column instead of a
% row.  (Needed to construct a CHEBFUN out of @(x) feval(f, x) later on.)
isTrans = f.isTransposed;
if ( isTrans )
    f = f.';
end

% Parse the inputs:
if ( nargin == 1 )
    % Choose all indices by default:
    index = 2:numel(f.funs);
    % Obtain preferences:
    pref = chebfunpref();
elseif ( nargin == 2 )
    if ( isstruct(index) || isa(index, 'chebfunpref') ) % MERGE(F, PREF)
        % index actually is a struct of preferences or a CHEBFUNPREF
        pref = chebfunpref(index);
        % Choose all indices by default:
        index = 2:numel(f.funs);
    else                   % MERGE(F, INDEX)
        % indices passed, obtain preferences:
        pref = chebfunpref();
    end
end

[f, mergedPts] = mergeColumn(f, index, pref);

% Convert back to a row CHEBFUN if we started with one.
if ( isTrans )
    f.isTransposed = true;
end

end

function [f, mergedPts] = mergeColumn(f, index, pref)

% Deal with input arguments:
if ( isempty(index) )
    % No indices requested.
    mergedPts = [];
    return

elseif ( ischar(index) )
    % 'all' indices requested:
    index = 2:numel(f.funs);

else
    % Index of endpoints was provided:
    index = unique(index);
    % Break points must be in range 2:numel(f.funs)
    index((index <= 1) | (index > numel(f.funs))) = [];

    if ( isempty(index) )
        % All the requested indices were trivial.
        mergedPts = [];
        return
    end

end

% Determine the maximum length of the merged pieces:
if ( ~pref.splitting )
    maxn = pref.techPrefs.maxLength;
else
    maxn = pref.splitPrefs.splitLength;
end
pref.techPrefs.maxLength = maxn;

% Splitting forces extrapolate:
if ( pref.splitting )
    pref.techPrefs.extrapolate = true;
end

% Obtain scales of the CHEBFUN:
vs = vscale(f);
hs = hscale(f);
pref.chebfuneps = max(eps, pref.chebfuneps);
mergedPts = [];

% Store data from input CHEBFUN:
oldPointVals = f.pointValues;
oldDom = f.domain;
oldFuns = f.funs;
newPointVals = oldPointVals;
newDom = oldDom;
newFuns = oldFuns;

% Loop through the index:
for k = index
    % Find corresponding break:
    j = find(oldDom(k) == newDom, 1, 'first');
    % And lengths of funs on either side:
    lengthPrevFun = length(newFuns{j-1});
    lengthThisFun = length(newFuns{j});

    % Prevent merge if existing FUN lengths add to more than 1.2*maxn:
    if ( lengthPrevFun + lengthThisFun >= 1.2*maxn )
        % Skip to next breakpoint:
        continue
    end

    % Prevent merging if there are jumps:
    v = [oldPointVals(k,:); get(oldFuns{k-1}, 'rval'); get(oldFuns{k}, 'lval')];
    jumps = norm(v([1, 1],:) - v(2:3,:), inf)./vs;
    if ( all( jumps >= 1e3*pref.chebfuneps ) || any(isinf(v(:))) )
        % Skip to next breakpoint:
        continue
    end

    % Call merge at the FUN level:
    [mergedFun, ishappy] = merge(newFuns{j-1}, newFuns{j}, vs, hs, pref);

    % Prevent merge if the result is not happy:
    if ( ~ishappy )
        % Skip to next breakpoint:
        continue
    end

    % Merging successful!
    mergedPts = [mergedPts, k]; %#ok<AGROW> % Store index of the removed point.
    newFuns{j-1} = mergedFun;   % Insert new fun.
    newFuns(j) = [];            % Remove unneeded fun.
    newDom = [newDom(1:j-1), newDom(j+1:end)]; % Update domain.
    % Update pointValues.
    newPointVals = [newPointVals(1:j-1,:) ; newPointVals(j+1:end,:)];

end

% Assign data to CHEBFUN:
f.domain = newDom;
f.funs = newFuns;
f.pointValues = newPointVals;

end


text:
 function [yik, scalefactor] = ir_mri_field_map_reg_scale(yik, etime, varargin)
%function [yik, scalefactor] = ir_mri_field_map_reg_scale(yik, etime, varargin)
%|
%| Scale images to account for R2* effects and differences in absolute value
%| using median(ri), where ri = sum_j sum_k |yik_j yik_k|^2 (t_k - t_j)^2
%|
%| in
%|	yik	[N nset]	scan images
%|	etime	[1 nset]	echo times (units of sec if fieldmap is in Hz)
%|
%| option
%|	fmax			threshold for absolute yik value (default 0.1)
%|	dmax			threshold for absolute rj value (default 0.1)
%|	show	0|1		1 to show result (default 0)
%|
%| out
%|	yik	[N nset]	scaled scan images
%|	scalefactor		sqrt(median(rj)) * effect_of_fmax
%|
%| MJ Allison
%| 2015-06-27 JF cosmetic changes

if nargin < 4, ir_usage, end

arg.fmax = 0.1;
arg.dmax = 0.1;
arg.show = 0;
arg = vararg_pair(arg, varargin);

dim_yik = size(yik);
yik = reshapee(yik, [], dim_yik(end)); % [*N nset]

[nn, nset] = size(yik);

% Scale by median of first set of data to get rid of large mag_j
% effects (not actually needed, but left in for consistency)
if arg.fmax > 0
	y1 = abs(yik(:,1));
	scalefactor = median(y1(y1(:) > arg.fmax * max(y1(:))));
	if scalefactor == 0
		fail 'median is zero?'
	end
	yik = yik / scalefactor;
else
	scalefactor = 1;
end

% Try to compensate for R2 effects on effective regularization.

d = zeros(nn,1);
for j = 1:nset
	for k = 1:nset
		tmp = abs(yik(:,j) .* yik(:,k)).^2 * (etime(k)-etime(j)).^2;
		d = d + tmp;
	end
end

% divide by numerator of wj^mn -> sum(abs(y)^2)
d = div0(d, sum(abs(yik).^2,2));

if arg.show
	im(reshape(d, [dim_yik(1:2)]))
end

% compute dtyp
dtyp = median(d(d(:) > arg.dmax * max(d(:))));
% pr dtyp

% now uniformly scale by the square root of dtyp
yik = yik / sqrt(dtyp);

scalefactor = scalefactor * sqrt(dtyp); % 2022-11-27 fix
yik = reshape(yik, dim_yik); % [(N) nset]


text:
function bnet = mk_alarm_bnet()

% Written by Qian Diao <qian.diao@intel.com> on 11 Dec 01

N = 37;
dag = zeros(N,N);
dag(21,23) = 1 ;
dag(21,24) = 1 ;
dag(1,24) = 1 ;
dag(1,23) = 1 ;
dag(2,26) = 1 ;
dag(2,25) = 1 ;
dag(2,24) = 1 ;
dag(2,13) = 1 ;
dag(2,23) = 1 ;
dag(13,30) = 1 ;
dag(30,31) = 1 ;
dag(3,14) = 1 ;
dag(3,19) = 1 ;
dag(4,36) = 1 ;
dag(14,35) = 1 ;
dag(32,33) = 1 ;
dag(32,35) = 1 ;
dag(32,34) = 1 ;
dag(32,36) = 1 ;
dag(15,21) = 1 ;
dag(5,31) = 1 ;
dag(27,30) = 1 ;
dag(28,31) = 1 ;
dag(28,29) = 1 ;
dag(26,28) = 1 ;
dag(26,27) = 1 ;
dag(16,31) = 1 ;
dag(16,37) = 1 ;
dag(23,26) = 1 ;
dag(23,29) = 1 ;
dag(23,25) = 1 ;
dag(6,15) = 1 ;
dag(7,27) = 1 ;
dag(8,21) = 1 ;
dag(19,20) = 1 ;
dag(19,22) = 1 ;
dag(31,32) = 1 ;
dag(9,14) = 1 ;
dag(9,17) = 1 ;
dag(9,19) = 1 ;
dag(10,33) = 1 ;
dag(10,34) = 1 ;
dag(11,16) = 1 ;
dag(12,13) = 1 ;
dag(12,18) = 1 ;
dag(35,37) = 1 ;

node_sizes = 2*ones(1,N);
node_sizes(2) = 3;
node_sizes(6) = 3;
node_sizes(14) = 3;
node_sizes(15) = 4;
node_sizes(16) = 3;
node_sizes(18) = 3;
node_sizes(19) = 3;
node_sizes(20) = 3;
node_sizes(21) = 4;
node_sizes(22) = 3;
node_sizes(23) = 4;
node_sizes(24) = 4;
node_sizes(25) = 4;
node_sizes(26) = 4;
node_sizes(27) = 3;
node_sizes(28) = 3;
node_sizes(29) = 4;
node_sizes(30) = 3;
node_sizes(32) = 3;
node_sizes(33) = 3;
node_sizes(34) = 3;
node_sizes(35) = 3;
node_sizes(36) = 3;
node_sizes(37) = 3;

bnet = mk_bnet(dag, node_sizes);

bnet.CPD{1} = tabular_CPD(bnet, 1,[0.96 0.04 ]);
bnet.CPD{2} = tabular_CPD(bnet, 2,[0.92 0.03 0.05 ]);
bnet.CPD{3} = tabular_CPD(bnet, 3,[0.8 0.2 ]);
bnet.CPD{4} = tabular_CPD(bnet, 4,[0.95 0.05 ]);
bnet.CPD{5} = tabular_CPD(bnet, 5,[0.8 0.2 ]);
bnet.CPD{6} = tabular_CPD(bnet, 6,[0.01 0.98 0.01 ]);
bnet.CPD{7} = tabular_CPD(bnet, 7,[0.01 0.99 ]);
bnet.CPD{8} = tabular_CPD(bnet, 8,[0.95 0.05 ]);
bnet.CPD{9} = tabular_CPD(bnet, 9,[0.95 0.05 ]);
bnet.CPD{10} = tabular_CPD(bnet, 10,[0.9 0.1 ]);
bnet.CPD{11} = tabular_CPD(bnet, 11,[0.99 0.01 ]);
bnet.CPD{12} = tabular_CPD(bnet, 12,[0.99 0.01 ]);
bnet.CPD{13} = tabular_CPD(bnet, 13,[0.95 0.95 0.05 0.1 0.1 0.01 0.05 0.05 0.95 0.9 0.9 0.99 ]);
bnet.CPD{14} = tabular_CPD(bnet, 14,[0.05 0.95 0.5 0.98 0.9 0.04 0.49 0.01 0.05 0.01 0.01 0.01 ]);
bnet.CPD{15} = tabular_CPD(bnet, 15,[0.01 0.01 0.01 0.97 0.01 0.01 0.01 0.97 0.01 0.01 0.01 0.97 ]);
bnet.CPD{16} = tabular_CPD(bnet, 16,[0.3 0.98 0.4 0.01 0.3 0.01 ]);
bnet.CPD{17} = tabular_CPD(bnet, 17,[0.99 0.1 0.01 0.9 ]);
bnet.CPD{18} = tabular_CPD(bnet, 18,[0.05 0.01 0.9 0.19 0.05 0.8 ]);
bnet.CPD{19} = tabular_CPD(bnet, 19,[0.05 0.98 0.01 0.95 0.9 0.01 0.09 0.04 0.05 0.01 0.9 0.01 ]);
bnet.CPD{20} = tabular_CPD(bnet, 20,[0.95 0.04 0.01 0.04 0.95 0.29 0.01 0.01 0.7 ]);
bnet.CPD{21} = tabular_CPD(bnet, 21,[0.97 0.97 0.01 0.97 0.01 0.97 0.01 0.97 0.01 0.01 0.97 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.97 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.97 0.01 ]);
bnet.CPD{22} = tabular_CPD(bnet, 22,[0.95 0.04 0.01 0.04 0.95 0.04 0.01 0.01 0.95 ]);
bnet.CPD{23} = tabular_CPD(bnet, 23,[0.97 0.97 0.97 0.97 0.97 0.97 0.01 0.95 0.97 0.97 0.01 0.95 0.01 0.4 0.97 0.97 0.01 0.5 0.01 0.3 0.97 0.97 0.01 0.3 0.01 0.01 0.01 0.01 0.01 0.01 0.97 0.03 0.01 0.01 0.97 0.03 0.01 0.58 0.01 0.01 0.01 0.48 0.01 0.68 0.01 0.01 0.01 0.68 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.97 0.01 0.01 0.01 0.97 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.97 0.01 0.01 0.01 0.97 0.01 ]);
bnet.CPD{24} = tabular_CPD(bnet, 24,[0.97 0.97 0.97 0.97 0.97 0.97 0.01 0.01 0.4 0.1 0.01 0.01 0.01 0.01 0.2 0.05 0.01 0.01 0.01 0.01 0.2 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.97 0.49 0.58 0.84 0.9 0.29 0.01 0.01 0.75 0.25 0.01 0.01 0.01 0.01 0.7 0.15 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.3 0.01 0.05 0.08 0.3 0.97 0.08 0.04 0.25 0.38 0.08 0.01 0.01 0.09 0.25 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.2 0.01 0.01 0.01 0.4 0.01 0.9 0.01 0.45 0.6 0.9 0.97 0.97 0.01 0.59 0.97 0.97 ]);
bnet.CPD{25} = tabular_CPD(bnet, 25,[0.97 0.97 0.97 0.01 0.6 0.01 0.01 0.5 0.01 0.01 0.5 0.01 0.01 0.01 0.01 0.97 0.38 0.97 0.01 0.48 0.01 0.01 0.48 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.97 0.01 0.97 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.97 0.01 0.97 ]);
bnet.CPD{26} = tabular_CPD(bnet, 26,[0.97 0.97 0.97 0.01 0.01 0.03 0.01 0.01 0.01 0.97 0.01 0.01 0.01 0.01 0.01 0.97 0.97 0.95 0.01 0.01 0.94 0.01 0.01 0.88 0.01 0.01 0.01 0.01 0.01 0.01 0.97 0.97 0.04 0.01 0.01 0.1 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.97 0.97 0.01 ]);
bnet.CPD{27} = tabular_CPD(bnet, 27,[0.98 0.98 0.98 0.98 0.95 0.01 0.95 0.01 0.01 0.01 0.01 0.01 0.04 0.95 0.04 0.01 0.01 0.01 0.01 0.01 0.01 0.04 0.01 0.98 ]);
bnet.CPD{28} = tabular_CPD(bnet, 28,[0.01 0.01 0.04 0.9 0.01 0.01 0.92 0.09 0.98 0.98 0.04 0.01 ]);
bnet.CPD{29} = tabular_CPD(bnet, 29,[0.97 0.01 0.01 0.01 0.97 0.01 0.01 0.01 0.97 0.01 0.01 0.01 0.01 0.97 0.97 0.97 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.97 0.97 0.97 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.97 0.97 0.43 ]);
bnet.CPD{30} = tabular_CPD(bnet, 30,[0.98 0.98 0.01 0.98 0.01 0.69 0.01 0.01 0.98 0.01 0.01 0.3 0.01 0.01 0.01 0.01 0.98 0.01 ]);
bnet.CPD{31} = tabular_CPD(bnet, 31,[0.05 0.01 0.05 0.01 0.05 0.01 0.05 0.01 0.05 0.01 0.05 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.1 0.01 0.95 0.01 0.95 0.05 0.1 0.01 0.95 0.01 0.95 0.05 0.1 0.01 0.3 0.01 0.3 0.01 0.95 0.01 0.99 0.05 0.95 0.05 0.95 0.01 0.99 0.05 0.99 0.05 0.3 0.01 0.99 0.01 0.3 0.01 0.95 0.99 0.95 0.99 0.95 0.99 0.95 0.99 0.95 0.99 0.95 0.99 0.99 0.99 0.99 0.99 0.99 0.99 0.9 0.99 0.05 0.99 0.05 0.95 0.9 0.99 0.05 0.99 0.05 0.95 0.9 0.99 0.7 0.99 0.7 0.99 0.05 0.99 0.00999999 0.95 0.05 0.95 0.05 0.99 0.01 0.95 0.01 0.95 0.7 0.99 0.01 0.99 0.7 0.99 ]);
bnet.CPD{32} = tabular_CPD(bnet, 32,[0.1 0.01 0.89 0.09 0.01 0.9 ]);
bnet.CPD{33} = tabular_CPD(bnet, 33,[0.98 0.33333334 0.01 0.33333334 0.01 0.33333334 0.01 0.33333334 0.98 0.33333334 0.01 0.33333334 0.01 0.33333334 0.01 0.33333334 0.98 0.33333334 ]);
bnet.CPD{34} = tabular_CPD(bnet, 34,[0.98 0.33333334 0.01 0.33333334 0.01 0.33333334 0.01 0.33333334 0.98 0.33333334 0.01 0.33333334 0.01 0.33333334 0.01 0.33333334 0.98 0.33333334 ]);
bnet.CPD{35} = tabular_CPD(bnet, 35,[0.98 0.95 0.3 0.95 0.04 0.01 0.8 0.01 0.01 0.01 0.04 0.69 0.04 0.95 0.3 0.19 0.04 0.01 0.01 0.01 0.01 0.01 0.01 0.69 0.01 0.95 0.98 ]);
bnet.CPD{36} = tabular_CPD(bnet, 36,[0.98 0.98 0.01 0.4 0.01 0.3 0.01 0.01 0.98 0.59 0.01 0.4 0.01 0.01 0.01 0.01 0.98 0.3 ]);
bnet.CPD{37} = tabular_CPD(bnet, 37,[0.98 0.98 0.3 0.98 0.1 0.05 0.9 0.05 0.01 0.01 0.01 0.6 0.01 0.85 0.4 0.09 0.2 0.09 0.01 0.01 0.1 0.01 0.05 0.55 0.01 0.75 0.9 ]);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
 function [opt, extra] = vararg_pair(opt, varargs, varargin)
%function [opt, extra] = vararg_pair(opt, varargs, [options])
%|
%| Process name / value pairs, replacing the "default" field values
%| of the opt structure with the user-specified values.
%| This allows flexible argument order and "named arguments" somewhat like IDL.
%| This newer version allows option names to differ from the field names.
%| If two output arguments are given, then any "extra" name value pairs
%| that are not associated with the fields in structure opt will be
%| returned as name/value pairs in the "extra" output.
%|
%| in
%|	opt	struct	opt.name1 opt.name2 etc., containing default values
%|	varargs	{'name1', 'value1', ... } as a cell array
%|
%| options
%|	'subs'	Nx2 cell	{'arg_name1', 'struct_name1'; ...}
%|		short argument name followed by longer structure name
%|	'allow_new'	0|1	allow new structure elements? (default: 0)
%|
%| out
%|	opt		opt.(name) = value
%|	extra		{'name_i', 'value_i', ...} all "unused" pairs
%|
%| Copyright 2005-06-18, Jeff Fessler, University of Michigan
%| Modified 2016-04-06, Gopal Nataraj, University of Michigan

if nargin < 1, ir_usage, end
if nargin == 1 && streq(opt, 'test'), vararg_pair_test, clear, return, end

if isempty(varargs) && numel(varargin) == 0
	extra = {};
	return
end

base = [caller_name ' arg: ']; % for printing

local.subs = {};
local.allow_new = 0;
local = vararg_pair_do(local, varargin, {}, 0, 0);

allow_extra = nargout > 1;
[opt, extra] = vararg_pair_do(opt, varargs, local.subs, local.allow_new, ...
	allow_extra, base);

% Structure elements of opt with value "NaN" are checked as being required.
% NO, this just adds overhead. Removed 2007-2-9.
% vararg_pair_check(opt);

% vararg_pair_do()
function [opt, extra] = ...
	vararg_pair_do(opt, varargs, subs, allow_new, allow_extra, base)

if allow_extra && allow_new
	error 'only one of "new" or "extra" is sensible for vararg_pair'
end

npair = floor(length(varargs) / 2);
if 2*npair ~= length(varargs), fail('need names and values in pairs'), end
args = varargs(1:2:end);
vals = varargs(2:2:end);

if ~isempty(subs)
	if size(subs,2) ~= 2, fail('subs must be Nx2'), end
%	subs1 = strvcat(subs(:,1));
	subs1 = char(subs(:,1)); % strmatch doc says fastest for char array!
	subs2 = subs(:,2);
	% todo: apply 'subs' all at once here?
end

extra = {};
for ii=1:npair
	arg = args{ii};
	val = vals{ii};

	if ~ischar(arg)
		fail('unknown option of class %s', class(arg))
	end

	if ~isempty(subs)
		arg = vararg_pair_sub(arg, subs1, subs2);
	end

	[,opt ex] = vararg_pair_put(opt, arg, val, allow_new, allow_extra);
	if allow_extra && ~isempty(ex)
		extra = [extra(:); ex(:)]'; % faster than {extra{:}, ex{:}};
	end

	if isfield(opt, 'chat') && opt.chat && isempty(ex)
		show_pair(arg, val, base)
	end
end

% vararg_pair_sub()
% substitute option name for field name, if a Nx2 cell array is provided.
%
function arg = vararg_pair_sub(arg, subs1, subs2)
ii = strmatch(arg, subs1, 'exact');
if isempty(ii)
	return
elseif length(ii) == 1
	arg = subs2{ii};
else
	printm([subs1(ii) ' ' subs2{ii}])
	error 'double match? only one substitution allowed'
end

% vararg_pair_put()
% opt.(arg) = value
% with extra effort to handle opt.(arg1.arg2) = value
% because matlab does not handle that case by itself, except via subsagn
%
function [opt, extra] = vararg_pair_put(opt, arg, value, allow_new, allow_extra)

extra = {};
if isfield(opt, arg) % simple case of name / value pair
	opt.(arg) = value;
	return
end

idot = strfind(arg, '.');

% old: jf, 2005
%{
if isempty(idot)
	if allow_new
		opt.(arg) = value;
		return
	elseif allow_extra
		extra = {arg, value};
		return
	else
		fail('unknown option name "%s"', arg)
	end
end

% tricky case of a.b.c
if length(idot) > 1, error 'a.b.c.d not done', end

arg1 = arg([1:(idot-1)]);
arg2 = arg([(idot+1):end]);
if ~isfield(opt, arg1)
	fail('unknown option1 name "%s"', arg1)
end
if ~isfield(opt.(arg1), arg2)
	fail('unknown option2 name %s', arg2)
end
s = struct('type', {'.', '.'}, 'subs', ...
	{ arg([1:(idot-1)]), arg([(idot+1):end]) });
try
	opt = subsasgn(opt, s, value);
catch
	printf(lasterr)
	fail('subsasgn? unknown option name %s', arg)
end
%}

% new: gn, 2016
switch length(idot)
	case 0
		if allow_new
			opt.(arg) = value;
			return
		elseif allow_extra
			extra = {arg, value};
			return
		else
			fail('unknown option name "%s"', arg)
		end

	case 1
		arg1 = arg([1:(idot-1)]);
		arg2 = arg([(idot+1):end]);
		if ~isfield(opt, arg1)
			fail('unknown option1 name "%s"', arg1)
		end
		if ~isfield(opt.(arg1), arg2)
			fail('unknown option2 name "%s"', arg2)
		end
		s = struct('type', {'.', '.'}, 'subs', {arg1, arg2});
		try
			opt = subsasgn(opt, s, value);
		catch
			 printm(lasterr)
			 fail('subsasgn? unknown option name %s', arg)
		end

	case 2
		arg1 = arg([1:(idot(1)-1)]);
		arg2 = arg([(idot(1)+1):(idot(2)-1)]);
		arg3 = arg([(idot(2)+1):end]);
		if ~isfield(opt, arg1)
			fail('unknown option1 name "%s"', arg1)
		end
		if ~isfield(opt.(arg1), arg2)
			fail('unknown option2 name "%s"', arg2)
		end
		if ~isfield(opt.(arg1).(arg2), arg3)
			fail('unknown option3 name "%s"', arg3)
		end
		s = struct('type', {'.', '.', '.'},...
			'subs', {arg1, arg2, arg3});
		try
			opt = subsasgn(opt, s, value);
		catch
			printf(lasterr)
			fail('subsasgn? unknown option name %s', arg)
		end

	otherwise
		error('a.b.c.d.e not done');
end

% vararg_pair_check()
% check that the user supplied the required arguments,
% which are signified by 'nan' values.
%
function vararg_pair_check(opt)
names = fieldnames(opt);
ok = 1;
for ii=1:length(names)
	t = opt.(names{ii});
	if isnumeric(t) && any(isnan(t(:)))
		printf('Required option "%s" missing', names{ii});
		ok = 0;
	end
end
if ~ok, fail 'missing required "option(s)"', end

% show_pair()
function show_pair(name, value, base)

pri = @(varargin) printf([base varargin{1}], varargin{2:end});

if isnumeric(value) || islogical(value)
	if max(size(value)) == 1	% only print scalars
		if ~streq(name, 'chat')
			pri('%s = %g', name, value)
		end
	else	% otherwise print sum (e.g., for mask)
		pri('sum(%s) = %g', name, sum(value(:)))
	end

elseif ischar(value)
	pri('%s = %s', name, value)

elseif isa(value, 'cell')
	pri('%s {cell}', name)

elseif isa(value, 'function_handle')
	pri('%s = @%s', name, func2str(value))

elseif isstruct(value)	% show first structure element, if string
	fields = fieldnames(value);
%	name = getfield(arg, fields{1});
	name1 = value.(fields{1});
	if ischar(name1)
		pri('%s.%s = %s', name, fields{1}, name1)
	end

else
	warn('display type "%s" not implemented', class(value))
end

% test routine
function vararg_pair_test
args = {'a', 1, 'b', 2, 'c.d', 3, 'g', 9, 'l', 5, 'll', 6};%, 'c.f', 4};
opt.a = 0;
opt.b = 0;
opt.f = 0;
opt.long1 = 0;
opt.long2 = 0;
%opt.g = nan;
opt.c.d = 0;
opt.c.e = 0;
sub = {'g', 'f'; 'l', 'long1'; 'll', 'long2'};
opt = vararg_pair(opt, args, 'subs', sub);
jf_equal(opt.a, 1)
jf_equal(opt.b, 2)
jf_equal(opt.c.d, 3)
jf_equal(opt.c.e, 0)
jf_equal(opt.f, 9)
jf_equal(opt.long1, 5)
jf_equal(opt.long2, 6)

opt = vararg_pair(opt, {'z', 'new'}, 'allow_new', 1);
if ~isequal(opt.z, 'new'), fail('bug'), end
try % typo in name should yield error
	opt = vararg_pair(opt, {'verify_bad_test!', 'new'}, 'allow_new', 0)
catch
	printm 'test passed'
	return
end
fail 'should not get here'


text:
classdef BT1 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP with bias feature

%------------------------------- Reference --------------------------------
% H. Li, Q. Zhang, and J. Deng, Biased multiobjective optimization and
% decomposition algorithm, IEEE Transactions on Cybernetics, 2017, 47(1):
% 52-66.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            [N,D] = size(X);
            I1    = 2 : 2 : D;
            I2    = 3 : 2 : D;
            Y     = X - sin(repmat(1:D,N,1)*pi/2/D);
            PopObj(:,1) = X(:,1)         + sum(Y(:,I1).^2+(1-exp(-Y(:,I1).^2/1e-10))/5,2);
            PopObj(:,2) = 1-sqrt(X(:,1)) + sum(Y(:,I2).^2+(1-exp(-Y(:,I2).^2/1e-10))/5,2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - R(:,1).^0.5;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO