text:
classdef prtRvUniformImproper < prtRv
    % prtRvUniformImproper  Improper uniform random variable
    %
    %   RV = prtRvUniformImproper creates a prtRvUniformImproper object
    %   with unknown dimensionality nDimensions. nDimensions can be set
    %   manually or using the MLE method. A prtRvUniformImproper
    %   models an improper pdf that always yields a value of 1 no matter
    %   the input. prtRvUniformImproper is sometimes useful for creating
    %   one class classifiers. See the examples below for more information
    %
    %   The draw method of prtRvUniformImproper draws values uniformly
    %   distributed from realmin to realmax in each dimension.
    %
    %   RV = prtRvUniformImproper(PROPERTY1, VALUE1,...) creates a
    %   prtRvMultinomial object RV with properties as specified by
    %   PROPERTY/VALUE pairs.
    %
    %   A prtRvUniformImproper object inherits all properties from the
    %   prtRv class. In addition, it has the following properties:
    %
    %   nDimensions - dimensionality of the data modeled by this RV.
    %
    %  A prtRvUniformImproper object inherits all methods from the prtRv
    %  class. The MLE  method can be used to set the parameters from data.
    %
    %  Example:
    %
    %  % In this example we show that the PDF of a prtRvUniformImproper is
    %  % always 1
    %  dataSet = prtDataGenUnimodal;        % Load a dataset consisting of
    %                                       % 2 features
    %  dataSet = retainFeatures(dataSet,1); % Retain only the first feature
    %                                       % only for the example.
    %
    %  RV = prtRvUniformImproper;           % Create a prtRvUniform object
    %  RV = RV.mle(dataSet);                % Compute the bounds
    %
    %  RV.plotPdf([-10 10]);                % We must manually specify
    %                                       % plot limits since
    %                                       % prtRvUniformImproper does not
    %                                       % have actual plot limits
    %
    %
    %  % In this example we show how to build a one class MAP classifier
    %  dataSet = prtDataGenUnimodal;        % Load a dataset consisting of
    %                                       % 2 features
    %
    %  % Create and train a GLRT classifier that uses a
    %  % prtRvUniformImproper to model class 0 and a prtRvMvn to model
    %  % class 1
    %  glrtClass = train(prtClassGlrt('rvH0',prtRvUniformImproper,'rvH1',prtRvMvn),dataSet);
    %
    %  plot(glrtClass) % Contours only show the log-likelihood of class 1
    %
    %   See also: prtRv, prtRvMvn, prtRvGmm, prtRvMultinomial,
    %   prtRvVq, prtRvKde

    properties (SetAccess = private)
        name = 'Improper Uniform Random Variable'
        nameAbbreviation = 'RVImUnif';
    end

    properties (SetAccess = protected)
        isSupervised = false;
        isCrossValidateValid = true;
    end

    properties (Hidden = true, SetAccess = 'private', GetAccess = 'private')
        nDimensionsPrivate
    end

    properties (Hidden = true, Dependent = true)
        nDimensions
    end

    methods
        % The Constructor
        function R = prtRvUniformImproper(varargin)
            R = constructorInputParse(R,varargin{:});
        end

        function R = mle(R,X)
            X = R.dataInputParse(X); % Basic error checking etc
           R.nDimensionsPrivate = size(X,2);
        end

        function vals = pdf(R,X)
            X = R.dataInputParse(X); % Basic error checking etc

            [isValid, reasonStr] = R.isValid;
            assert(isValid,'PDF cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            vals = ones(size(X,1),1);
        end

        function vals = logPdf(R,X)
            X = R.dataInputParse(X); % Basic error checking etc
            [isValid, reasonStr] = R.isValid;
            assert(isValid,'LOGPDF cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);
            vals = log(pdf(R,X));
        end

        function vals = cdf(R,X)
            X = R.dataInputParse(X); % Basic error checking etc
            [isValid, reasonStr] = R.isValid;
            assert(isValid,'CDF cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);
            assert(size(X,2) == R.nDimensions,'Data, RV dimensionality missmatch. Input data, X, has dimensionality %d and this RV has dimensionality %d.', size(X,2), R.nDimensions)
            vals = nan(size(X,1),1);
        end

        function vals = draw(R,N)
            if nargin < 2
                N = 1;
            end

            assert(numel(N)==1 && N==floor(N) && N > 0,'N must be a positive integer scalar.')

            [isValid, reasonStr] = R.isValid;
            assert(isValid,'DRAW cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            vals = bsxfun(@times,bsxfun(@times,rand(N,R.nDimensions),realmax),sign(randn(N,1)));
        end

        function R = set.nDimensions(R,N)
            assert(numel(N)==1 && N==floor(N) && N > 0,'nDimensions must be a positive integer scalar.')
            R.nDimensionsPrivate = N;
        end

        function val = get.nDimensions(R)
            val = R.nDimensionsPrivate;
        end
    end

    methods (Hidden = true)
        function [val, reasonStr] = isValid(R)
            if numel(R) > 1
                val = false(size(R));
                for iR = 1:numel(R)
                    [val(iR), reasonStr] = isValid(R(iR));
                end
                return
            end

            val = true;
            reasonStr = '';
        end
        function val = plotLimits(R)
            [isValid, reasonStr] = R.isValid;
            if isValid
                val = zeros(2*R.nDimensions,1);
                val(1:2:end) = Inf;
                val(2:2:end) = -Inf;
                % We send backwards limits so that we don't effect plot
                % limits if you were to check the limits of a set of RVs.
                % In this case this RV would not change the resulting
                % max(upperBounds), min(lowerBounds)
            else
                error('prtRvUniformImproper:plotLimits','Plotting limits can not be determined for this RV. It is not yet valid %s.',reasonStr)
            end
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function dxdt=ideal(t,x)
% Here you define your differential equation
% dxdt=your outputs, as initialized with zero array,
dxdt=zeros(2,1);

k=0.5;
u=-k*x(1);
% the differential eq is
%   x1dot=x2;
%   x2dot=u
% which can formulized as follows for the solver
dxdt(1)=x(2);
dxdt(2)=u;

author: Sable

repo: mcbench-benchmarks


text:
function [Ep,Cp] = spm_dcm_sparse(DCM,field)
% Bayesian model reduction of all permutations of model parameters
% FORMAT [RCM,BMR] = spm_dcm_sparse(DCM,field
%
% DCM      - A single estimated DCM (or PEB) structure:
%
%  DCM.M.pE  - prior expectation
%  DCM.M.pC  - prior covariance
%  DCM.Ep    - posterior expectation
%  DCM.Cp    - posterior covariances
%  DCM.gamma - prior variance    of reduced parameters (default: 0)
%
% field      - parameter fields in DCM{i}.Ep to optimise [default: {'A','B'}]
%             'All' will invoke all fields (i.e. random effects)
%             If Ep is not a structure, all parameters will be considered
%
% Returns:
%  Ep    - (BMA) posterior expectation
%  Cp    - (BMA) posterior covariance
%
%--------------------------------------------------------------------------
% This routine searches over reduced (nested) models of a full model (DCM)
% using Bayesian model reduction and performs Bayesian Model Averaging.
% 'Reduced' means some free parameters (parameters with a non-
% zero prior covariance) are switched off by fixing their prior variance
% to zero.This version incorporates a sparsity  prior over models (with a
% Gaussian hyperprior). In other words, the free energy is taken to be the
% likelihood of some data under a given model. The prior on that model
% corresponds to a softmax function of the prior entropy. Finally, the
% softmax (Gibbs) parameter is equipped with a Gaussian prior. Using
% Bayesian model reduction, this routine evaluates the joint probability
% over model and softmax sparsity parameter. The marginals over model space
% are then used to form Bayesian model averaging.
%
% The greedy search in this version simply evaluates the log evidence of
% models with and without each parameter and then successively removes the
% parameters with the least evidence.
%
% See also: spm_dcm_bmr and spm_dcm_bmr_all
%__________________________________________________________________________
% Copyright (C) 2010-2014 Wellcome Trust Centre for Neuroimaging

% Karl Friston, Peter Zeidman
% $Id: spm_dcm_sparse.m 7082 2017-05-27 19:36:36Z karl $

%-Number of parameters to consider before invoking greedy search
%--------------------------------------------------------------------------
nmax  = 8;

%-specification of null prior covariance
%--------------------------------------------------------------------------
if isfield(DCM,'beta'),  beta  = DCM.beta;  else, beta  = 0; end
if isfield(DCM,'gamma'), gamma = DCM.gamma; else, gamma = 0; end

%-Check fields of parameter stucture
%--------------------------------------------------------------------------
if nargin < 2 || isempty(field)
    field = {'A','B'};
end
if ischar(field)
    field = {field};
end

%-dela with filenames stucture
%--------------------------------------------------------------------------
if ischar(DCM)
    DCM = load(DCM,'DCM');
    DCM = DCM.DCM;
end

% Get prior covariances
%--------------------------------------------------------------------------
if isstruct(DCM.M.pC), DCM.M.pC = diag(spm_vec(DCM.M.pC)); end
if spm_length(DCM.M.pE) ~= size(DCM.M.pC,1)
    DCM.M.pC = diag(spm_vec(DCM.M.pC));
end

% Get priors and posteriors
%--------------------------------------------------------------------------
qE  = DCM.Ep;
qC  = DCM.Cp;
pE  = DCM.M.pE;
pC  = DCM.M.pC;

% Remove (a priori) null space
%--------------------------------------------------------------------------
U   = spm_svd(pC);
qE  = U'*spm_vec(qE);
pE  = U'*spm_vec(pE);
qC  = U'*qC*U;
pC  = U'*pC*U;

%-Greedy search (GS) - eliminating parameters in a top down fashion
%==========================================================================

% Accumulated reduction vector (C)
%--------------------------------------------------------------------------
q   = diag(DCM.M.pC);
if sum(q < 1024)
    C   = double(q > mean(q(q < 1024))/1024);
else
    C   = double(q > 0);
end

%-Find free coupling parameters
%----------------------------------------------------------------------
if isstruct(DCM.Ep)
    k = spm_fieldindices(DCM.Ep,field{:});
else
    k = 1:spm_length(DCM.Ep);
end
k     = k(find(C(k))); %#ok<FNDSB>

% Model search over new prior without the i-th parameter
%------------------------------------------------------------------
nparam = length(k);
for i  = 1:nparam

    % Identify parameters to retain r and to remove s
    %--------------------------------------------------------------
    r    = C; r(k(i)) = 0; s = 1 - r;

    % Create reduced priors
    %--------------------------------------------------------------
    R    = U'*diag(r + s*gamma)*U;
    rC   = R*pC*R;
    F(i) = spm_log_evidence(qE,qC,pE,pC,pE,rC);

end

% Find parameters with the least evidence
%--------------------------------------------------------------------------
[F,i] = sort(-F);
k     = k(i);
M     = cell(0);
for i = 1:nparam

    % parameters to retain (r) and to remove (s)
    %----------------------------------------------------------------------
    r    = C; r(k(1:i)) = 0; s = 1 - r;

    % Create reduced prior covariance matrix
    %----------------------------------------------------------------------
    R    = U'*diag(r + s*gamma)*U;
    rC   = R*pC*R;

    % record
    %----------------------------------------------------------------------
    M(i).F  = spm_log_evidence(qE,qC,pE,pC,pE,rC)
    M(i).H  = spm_logdet(rC);
    M(i).rC = rC;

end

% Sparsity hyperpriors
%--------------------------------------------------------------------------
s     = (1:64)/64;
Ps    = exp(-((1:64) - 32).^2/(2*16));
Ps    = Ps/sum(Ps);

%  model likelihood, model prior and  sparsity hyperprior
%--------------------------------------------------------------------------
Lm    = spm_softmax(spm_vec(M.F));
for i = 1:numel(s)
    Pm      = spm_softmax(-s(i)*spm_vec(M.H));
    Qm(:,i) = Lm.*Pm*Ps(i);
end

% evidence and log evidence
%--------------------------------------------------------------------------
Qm    = Qm/sum(sum(Qm));
G     = log(sum(Qm,2));

%-Bayesian model average
%==========================================================================
qE    = DCM.Ep;
qC    = DCM.Cp;
pE    = DCM.M.pE;
pC    = DCM.M.pC;
pE    = spm_vec(pE);
Gmax  = max(G);
BMA   = {};
for i = 1:length(G)
    if G(i) > (Gmax - 4)
        [F,Ep,Cp]    = spm_log_evidence_reduce(qE,qC,pE,pC,pE,M(i).rC);
        BMA{end + 1} = struct('Ep',Ep,'Cp',Cp,'F',F);
    end
end

BMA   = spm_dcm_bma(BMA);
Ep    = BMA.Ep;
Cp    = BMA.Cp;
if isstruct(Cp) || (spm_length(Cp) == spm_length(Ep))
    Cp = diag(spm_vec(Cp));
end

% Show results
% -------------------------------------------------------------------------
GRAPHICS = 1;

if ~GRAPHICS, return, end

spm_figure('Getwin','BMR - all'); clf
subplot(3,2,1)
imagesc(log(Qm)')
title('Joint density','FontSize',16)
xlabel('model','FontSize',12)
ylabel('sparsity','FontSize',12)
axis square

subplot(3,2,2)
plot(DCM.Ep,Ep,'.',DCM.Ep,DCM.Ep,':')
title('Full and reduced expectations','FontSize',16)
xlabel('Full expectations','FontSize',12)
ylabel('Reduced expectations','FontSize',12)
axis square

subplot(3,2,3)
plot(s,sum(Qm,1))
title('Marginal over sparsity','FontSize',16)
xlabel('sparsity parameter','FontSize',12)
ylabel('probability','FontSize',12)
axis square

subplot(3,2,4)
plot(sum(Qm,2))
title('Marginal over  model','FontSize',16)
xlabel(' model','FontSize',12)
ylabel('probability','FontSize',12)
axis square
drawnow

author: spm


text:
% --------- DeepMIMO: A Generic Dataset for mmWave and massive MIMO ------%
% Authors: Umut Demirhan, Abdelrahman Taha, Ahmed Alkhateeb
% Date: March 17, 2022
% Goal: Encouraging research on ML/DL for MIMO applications and
% providing a benchmarking tool for the developed algorithms
% ---------------------------------------------------------------------- %

function [channel_params,channel_params_BS,BS_loc]=read_raytracing(BS_ID, params, scenario_files)
%% Loading channel parameters between current active basesation transmitter and user receivers
filename_DoD=strcat(scenario_files, '.', int2str(BS_ID),'.DoD.mat');
filename_DoA=strcat(scenario_files, '.', int2str(BS_ID),'.DoA.mat');
filename_CIR=strcat(scenario_files, '.', int2str(BS_ID),'.CIR.mat');
filename_LOS=strcat(scenario_files, '.', int2str(BS_ID),'.LoS.mat');
filename_PL=strcat(scenario_files, '.', int2str(BS_ID),'.PL.mat');
filename_Loc=strcat(scenario_files, '.Loc.mat');

DoD_array=importdata(filename_DoD);
DoA_array=importdata(filename_DoA);
CIR_array=importdata(filename_CIR);
LOS_array=importdata(filename_LOS);
PL_array=importdata(filename_PL);
Loc_array=importdata(filename_Loc);

user_ids = double(params.user_ids);
user_last = double(max(user_ids));
num_paths = double(params.num_paths);
tx_power_raytracing = double(params.transmit_power_raytracing); %The transmit power in dBm used to generate the channel data
transmit_power = 30;

total_num_users=double(DoD_array(1));
pointer=0;

DoD_array(1)=[];
DoA_array(1)=[];
CIR_array(1)=[];
LOS_array(1)=[];

channel_params_all=struct('DoD_phi',[],'DoD_theta',[],'DoA_phi',[],'DoA_theta',[],'phase',[],'ToA',[],'power',[],'num_paths',[],'loc',[],'LoS_status',[]);
channel_params_all_BS =struct('DoD_phi',[],'DoD_theta',[],'DoA_phi',[],'DoA_theta',[],'phase',[],'ToA',[],'power',[],'num_paths',[],'loc',[],'LoS_status',[]);

dc = duration_check(params.symbol_duration);

user_count = 1;
[~,user_order] = sort(user_ids);
for Receiver_Number=1:user_last
    max_paths=double(DoD_array(pointer+2));
    num_path_limited=double(min(num_paths,max_paths));
    if ismember(Receiver_Number, user_ids)
        if (max_paths>0)
            Relevant_data_length=max_paths*4;
            Relevant_limited_data_length=num_path_limited*4;

            Relevant_DoD_array=DoD_array(pointer+3:pointer+2+Relevant_data_length);
            Relevant_DoA_array=DoA_array(pointer+3:pointer+2+Relevant_data_length);
            Relevant_CIR_array=CIR_array(pointer+3:pointer+2+Relevant_data_length);

            channel_params_all(user_order(user_count)).DoD_phi=Relevant_DoD_array(2:4:Relevant_limited_data_length);
            channel_params_all(user_order(user_count)).DoD_theta=Relevant_DoD_array(3:4:Relevant_limited_data_length);
            channel_params_all(user_order(user_count)).DoA_phi=Relevant_DoA_array(2:4:Relevant_limited_data_length);
            channel_params_all(user_order(user_count)).DoA_theta=Relevant_DoA_array(3:4:Relevant_limited_data_length);
            channel_params_all(user_order(user_count)).phase=Relevant_CIR_array(2:4:Relevant_limited_data_length);
            channel_params_all(user_order(user_count)).ToA=Relevant_CIR_array(3:4:Relevant_limited_data_length);
            channel_params_all(user_order(user_count)).power=1e-3*(10.^(0.1*( Relevant_CIR_array(4:4:Relevant_limited_data_length)+(transmit_power-tx_power_raytracing) )));
            channel_params_all(user_order(user_count)).num_paths=num_path_limited;
            channel_params_all(user_order(user_count)).loc=Loc_array(Receiver_Number,2:4);
            channel_params_all(user_order(user_count)).distance=PL_array(Receiver_Number,1);
            channel_params_all(user_order(user_count)).pathloss=PL_array(Receiver_Number,2);
            channel_params_all(user_order(user_count)).LoS_status=LOS_array(Receiver_Number);

            dc.add_ToA(channel_params_all(user_order(user_count)).power, channel_params_all(user_order(user_count)).ToA);

        else
            channel_params_all(user_order(user_count)).DoD_phi=[];
            channel_params_all(user_order(user_count)).DoD_theta=[];
            channel_params_all(user_order(user_count)).DoA_phi=[];
            channel_params_all(user_order(user_count)).DoA_theta=[];
            channel_params_all(user_order(user_count)).phase=[];
            channel_params_all(user_order(user_count)).ToA=[];
            channel_params_all(user_order(user_count)).power=[];
            channel_params_all(user_order(user_count)).num_paths=0;
            channel_params_all(user_order(user_count)).loc=Loc_array(Receiver_Number,2:4);
            channel_params_all(user_order(user_count)).distance=PL_array(Receiver_Number,1);
            channel_params_all(user_order(user_count)).pathloss=[];
            channel_params_all(user_order(user_count)).LoS_status=LOS_array(Receiver_Number);
        end
        user_count = double(user_count + 1);
    end
    pointer=double(pointer+max_paths*4+2);
end

dc.print_warnings('user', BS_ID);
dc.reset();

channel_params=channel_params_all(1,:);

%% Loading channel parameters between current active basesation transmitter and all the basestation receivers
if params.enable_BS2BSchannels
    filename_BSBS_DoD=strcat(scenario_files, '.', int2str(BS_ID),'.DoD.BSBS.mat');
    filename_BSBS_DoA=strcat(scenario_files, '.', int2str(BS_ID),'.DoA.BSBS.mat');
    filename_BSBS_CIR=strcat(scenario_files, '.', int2str(BS_ID),'.CIR.BSBS.mat');
    filename_BSBS_LOS=strcat(scenario_files, '.', int2str(BS_ID),'.LoS.BSBS.mat');
    filename_BSBS_PL=strcat(scenario_files, '.', int2str(BS_ID),'.PL.BSBS.mat');
    filename_BSBS_Loc=strcat(scenario_files, '.BSBS.RX_Loc.mat');

    DoD_BSBS_array=importdata(filename_BSBS_DoD);
    DoA_BSBS_array=importdata(filename_BSBS_DoA);
    CIR_BSBS_array=importdata(filename_BSBS_CIR);
    LOS_BSBS_array=importdata(filename_BSBS_LOS);
    PL_BSBS_array=importdata(filename_BSBS_PL);
    Loc_BSBS_array=importdata(filename_BSBS_Loc);

    num_paths_BSBS = double(params.num_paths);
    total_num_BSs=double(DoD_BSBS_array(1));
    BS_pointer=0;

    DoD_BSBS_array(1)=[];
    DoA_BSBS_array(1)=[];
    CIR_BSBS_array(1)=[];
    LOS_BSBS_array(1)=[];

    BS_count = 1;
    for Receiver_BS_Number=1:total_num_BSs
        max_paths_BSBS=double(DoD_BSBS_array(BS_pointer+2));
        num_path_BS_limited=double(min(num_paths_BSBS,max_paths_BSBS));
        if sum(Receiver_BS_Number == double(params.active_BS)) == 1 %Only read the channels related to the active basestation receivers

            if (max_paths_BSBS>0)
                Relevant_data_length_BSBS=max_paths_BSBS*4;
                Relevant_limited_data_length_BSBS=num_path_BS_limited*4;

                Relevant_DoD_BSBS_array=DoD_BSBS_array(BS_pointer+3:BS_pointer+2+Relevant_data_length_BSBS);
                Relevant_DoA_BSBS_array=DoA_BSBS_array(BS_pointer+3:BS_pointer+2+Relevant_data_length_BSBS);
                Relevant_CIR_BSBS_array=CIR_BSBS_array(BS_pointer+3:BS_pointer+2+Relevant_data_length_BSBS);

                channel_params_all_BS(BS_count).DoD_phi=Relevant_DoD_BSBS_array(2:4:Relevant_limited_data_length_BSBS);
                channel_params_all_BS(BS_count).DoD_theta=Relevant_DoD_BSBS_array(3:4:Relevant_limited_data_length_BSBS);
                channel_params_all_BS(BS_count).DoA_phi=Relevant_DoA_BSBS_array(2:4:Relevant_limited_data_length_BSBS);
                channel_params_all_BS(BS_count).DoA_theta=Relevant_DoA_BSBS_array(3:4:Relevant_limited_data_length_BSBS);
                channel_params_all_BS(BS_count).phase=Relevant_CIR_BSBS_array(2:4:Relevant_limited_data_length_BSBS);
                channel_params_all_BS(BS_count).ToA=Relevant_CIR_BSBS_array(3:4:Relevant_limited_data_length_BSBS);
                channel_params_all_BS(BS_count).power=1e-3*(10.^(0.1*( Relevant_CIR_BSBS_array(4:4:Relevant_limited_data_length_BSBS)+(transmit_power-tx_power_raytracing) )));
                channel_params_all_BS(BS_count).num_paths=num_path_BS_limited;
                channel_params_all_BS(BS_count).loc=Loc_BSBS_array(Receiver_BS_Number,2:4);
                channel_params_all_BS(BS_count).distance=PL_BSBS_array(Receiver_BS_Number,1);
                channel_params_all_BS(BS_count).pathloss=PL_BSBS_array(Receiver_BS_Number,2);
                channel_params_all_BS(BS_count).LoS_status=LOS_BSBS_array(Receiver_BS_Number);

                dc.add_ToA(channel_params_all_BS(BS_count).power, channel_params_all_BS(BS_count).ToA);

            else
                channel_params_all_BS(BS_count).DoD_phi=[];
                channel_params_all_BS(BS_count).DoD_theta=[];
                channel_params_all_BS(BS_count).DoA_phi=[];
                channel_params_all_BS(BS_count).DoA_theta=[];
                channel_params_all_BS(BS_count).phase=[];
                channel_params_all_BS(BS_count).ToA=[];
                channel_params_all_BS(BS_count).power=[];
                channel_params_all_BS(BS_count).num_paths=0;
                channel_params_all_BS(BS_count).loc=Loc_BSBS_array(Receiver_BS_Number,2:4);
                channel_params_all_BS(BS_count).distance=PL_BSBS_array(Receiver_BS_Number,1);
                channel_params_all_BS(BS_count).pathloss=[];
                channel_params_all_BS(BS_count).LoS_status=LOS_BSBS_array(Receiver_BS_Number);
            end
            BS_count = double(BS_count + 1);
        end
        BS_pointer=double(BS_pointer+max_paths_BSBS*4+2);
    end

    dc.print_warnings('BS', BS_ID);
    dc.reset()

end
channel_params_BS=channel_params_all_BS(1,:);

%% Loading current active basestation location
TX_Loc_array=importdata(strcat(scenario_files, '.TX_Loc.mat')); %Reading transmitter locations
BS_loc = TX_Loc_array(BS_ID,2:4); %Select current active basestation location

end

author: DeepMIMO

repo: DeepMIMO-matlab


text:
function a = r8gd_indicator ( n, ndiag, offset )

%*****************************************************************************80
%
%% R8GD_INDICATOR sets up a R8GD indicator matrix.
%
%  Discussion:
%
%    The R8GD storage format is suitable for matrices whose only nonzero entries
%    occur along a few diagonals, but for which these diagonals are not all
%    close enough to the main diagonal for band storage to be efficient.
%
%    In that case, we assign the main diagonal the offset value 0.
%    Each successive superdiagonal gets an offset value 1 higher, until
%    the highest superdiagonal (the A(1,N) entry) is assigned the offset N-1.
%    Similarly, the subdiagonals are assigned offsets of -1 through -(N-1).
%
%    Now, assuming that only a few of these diagonals contain nonzeros,
%    then for the I-th diagonal to be saved, we stored its offset in
%    OFFSET(I), and its entries in column I of the matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer NDIAG, the number of diagonals of the matrix
%    that are stored in the array.
%    NDIAG must be at least 1, and no more than 2 * N - 1.
%
%    Input, integer OFFSET(NDIAG), the offsets for the diagonal storage.
%
%    Output, real A(N,NDIAG), the R8GD matrix.
%
  fac = 10^( i4_log_10 ( n ) + 1 );

  for i = 1 : n
    for diag = 1 : ndiag
      j = i + offset(diag);
      if ( 1 <= j & j <= n )
        a(i,diag) = fac * i + j;
      else
        a(i,diag) = 0.0;
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef GDE3 < ALGORITHM
% <multi> <real/integer> <constrained/none>
% Generalized differential evolution 3

%------------------------------- Reference --------------------------------
% S. Kukkonen and J. Lampinen, GDE3: The third evolution step of
% generalized differential evolution, Proceedings of the IEEE Congress on
% Evolutionary Computation, 2005, 443-450.
%------------------------------- 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
        function main(Algorithm,Problem)
            %% Generate random population
            Population = Problem.Initialization();

            %% Optimization
            while Algorithm.NotTerminated(Population)
                Offspring  = OperatorDE(Problem,Population,Population(randi(Problem.N,1,Problem.N)),Population(randi(Problem.N,1,Problem.N)));
                Population = EnvironmentalSelection(Population,Offspring,Problem.N);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% For all clusters, correlation with nearest neighbor
C = FindCentroid(hfig);
coeffs = corr(C');

A = zeros(1,length(coeffs));
for i = 1:length(coeffs),
   coeffs(i,i) = nan;
   A(i) = max(coeffs(i,:));
end

figure; hist(A)

%% Distribution of within-cluster correlations
U = unique(gIX);
numU = length(U);
B = zeros(numU,3);
for i=1:numU,
    i
    IX = find(gIX == U(i));
    coeffs = corr(M(IX,:)');
    m = coeffs(:);
    B(i,1) = min(m);
    B(i,2) = mean(m);
    B(i,3) = median(m);
end

figure;hist(B)
%%
figure;
hist(B(:,2))
h = findobj(gca,'Type','patch');
h.FaceColor = [0.5 0.5 0.5];
h.EdgeColor = 'w';
xlim([0,1])
xlabel('average corr. within cluster')
ylabel('count')

%% Distribution of cluster sizes
U = unique(gIX);
numU = length(U);
C = zeros(numU,1);
for i=1:numU,
    i
    IX = find(gIX == U(i));
    C(i) = length(IX);
end
[N,edges] = histcounts(C,10:10:2100);
%%
figure;
h = bar(edges,[0,N])%+1)
set(gca,'XScale','log')
set(gca,'YScale','log')
% h = findobj(gca,'Type','patch');
h.FaceColor = [0.5 0.5 0.5];
h.EdgeColor = 'w';
xlim([5,1000])
ylim([-10,100])
xlabel('cluster size')
ylabel('number of clusters (+1 to differentiate 0 from 1)')
%%
figure;hold on;
for i = 1:length(N),
    plot([edges(i),edges(i)],[0,N(i)],'color',[1,0.5,0.5],'linewidth',6)
end
set(gca,'XScale','log')
% set(gca,'YScale','log')
xlim([5,1000])
ylim([0,85])
xlabel('cluster size')
ylabel('count')

%%
set(gca,'XScale','log')
set(gca,'YScale','log')
% h = findobj(gca,'Type','patch');
h.FaceColor = [0.5 0.5 0.5];
h.EdgeColor = 'w';
xlim([5,1000])
ylim([-10,100])
xlabel('cluster size')
ylabel('number of clusters (+1 to differentiate 0 from 1)')

%% Plot Foxel Characterization

figure;
% after second kmeans
counts = hist(gIX,1:max(gIX));
hist(counts,1:1:200);
title('Fish8: cluster sizes after 2nd kmeans');
text(20,300,'mean=10;median=4;mode=1;31 clusters>200')

author: xiuyechen

repo: FishExplorer


text:
% SARImageWriter is an abstract superclass describing the family of all
% writers of SAR raster data.  It was desgined particularly with complex
% SAR data in mind, but subclasses could be built to write detected SAR
% data as well.
%
% Written by: Wade Schwartzkopf, NGA/ID
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

classdef SARImageWriter < handle
    % Public class data.  Anyone has access to these in a read-only fashion.
    properties(SetAccess=protected, GetAccess=public)
        sicdmeta; % Collection and image formation metadata
        filename;
    end

    methods
        % Constructor.  Since this superclass constructor requires
        % arguments, constructors for derived classes will typically call
        % this explicitly (MATLAB's default is to call the superclass
        % constructor without argument.) to assure that the internal
        % metatdata is correctly set on creation.
        function obj = SARImageWriter(filename, sicdmeta)
            % Fundamental properties
            obj.filename = filename;
            obj.sicdmeta = sicdmeta;
        end
    end

    % All subclasses need a way to write pixel data
    methods (Abstract)
        status = write_chip(obj, data, start_indices);
    end

    methods (Static)
        % Utility function that mimics the behavior of MATLAB's FSEEK
        % function, but allows for seeking past the end of a file (just
        % appends zeros to the end of the file).
        status = fseek(fid, offset, origin);
    end

end

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

author: ngageoint


text:
function organizeFusedTextures_HN(pathWORK,nPatient,CTinv_cell,CTweight_mat,R_mat,scale_cell,algo_cell,Ng_mat)
% -------------------------------------------------------------------------
% function organizeFusedTextures_HN(pathWORK,nPatient,CTinv_cell,CTweight_mat,R_mat,scale_cell,algo_cell,Ng_mat)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes FUSED texture features for all patients, for all
% different combinations of the following texture extraction parameters:
% - CT Inv.: Inversion of CT intensities in the PET/CT fusion process.
% - CT weight: Weight given to MRI wavelet low-pass sub-bands in the
%               PET/CT fusion process.
% - WPBF ratio 'R': Ratio specifying the ratio of weight given to band-pass
%                   coefficients to the weigth given to the rest of
%                   coefficients (HHH and LLL) in the wavelet band-pass
%                   filtering process.
% - Scale: Resolution at which the ROI is isotropically resampled.
% - Quant. Algo.: Quantization algorithm used on analyzed ROI.
% - Ng: Number of gray-levels in the quantization process.
%
% Different extraction parameters are passed as arrays or cells in the
% function in order to test all possible combinations. This function is
% used for FUSED scans specifically. See Ref. [1,2] and 'prepareVolume.m'
% for more details.
%
% Textures are organized in one cell: 'textures_PET_CT', where each cell
% entry organizes the data for all patients for a specific extraction
% parameter combination.
%
% -- > For example, size(textures_PET_CT) is equal to:
% [numel(CTinv_cell),numel(CTweight_mat),numel(R_mat),numel(scale_cell),numel(algo_cell),numel(Ng_mat)]
%
% In addition, the Spearman's rank correlation between each texture
% features and lung metastases is computed and saved in the cells for each
% texture feature of each entry. Results are saved in the STS WORKSPACE
% with the names given above.
% -------------------------------------------------------------------------
% REFERENCES:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% [2] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% - pathWORK: Full path to the HN WORKSPACE directory.
% - nPatient: Number of patients to read.
% - CTinv_cell: Cell vector specifying the different 'CT Inv.' values to test.
% - CTweight_mat: Array vector specifying the different 'CT weight' values to test.
% - R_mat: Array vector specifying the different 'R' ratio values to test.
% - scale_cell: Cell vector specifying the different 'Scale' values to test.
% - algo_cell: Cell vector specifying the different 'Quant. Algo.' values to test.
% - Ng_mat: Array vector specifying the different 'Ng' values to test.
% -------------------------------------------------------------------------
% EXAMPLE:
% CTinv_cell = {'NoInv','Inv'};
% CTweight_mat = [1/4,1/3,1/2,2/3,3/4];
% R_mat = [1/2,2/3,1,3/2,2];
% scale_cell = {'pixelW',1,2,3,4,5};
% algo_cell = {'Uniform','Equal','Lloyd'};
% Ng_mat = [8,16,32,64];
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: July 2015
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
cd(pathWORK), load('outcomes')

% INITIALIZATION
scans = {'PET_CT'};
nScans = numel(scans);
suffix = ['_TEXT'];
nameOutcomes = fieldnames(outcomes);
nOutcomes = length(nameOutcomes);

for scan = 1:numel(scans)
    cd([pathWORK,'/TEXTURES'])
    temp = load('Patient1_PET_CT_TEXT'); temp = struct2cell(temp); temp = temp{1}; % In order to get the necessary 'nameType' and 'nameFeature' fields
    nameType = fieldnames(temp.Experiment1); nType = numel(nameType); % All texture types are the same
    text = cell(numel(CTinv_cell),numel(CTweight_mat),numel(R_mat),numel(scale_cell),numel(algo_cell),numel(Ng_mat)); % Initialization of master cell
    tempText = cell(1,nPatient); % Cell used to load patient textures only once
    for p = 1:nPatient
        load(['Patient',num2str(p),'_',scans{scan},suffix]) % Variable 'textures' is now in MATLAB workspace
        tempText{p} = textures;
    end
    experiment = 0;
    for i = 1:numel(CTinv_cell)
        for w = 1:numel(CTweight_mat)
            for r = 1:numel(R_mat)
                for s = 1:numel(scale_cell)
                    for a = 1:numel(algo_cell)
                        for n = 1:numel(Ng_mat)
                            text{i,w,r,s,a,n} = struct;
                            experiment = experiment + 1;
                            strExperiment = ['Experiment',num2str(experiment)];
                            for t = 1:nType-1 % Last field is 'parameters'
                                nameFeature = fieldnames(temp.(strExperiment).(nameType{t})); nFeature = numel(nameFeature);
                                for f = 1:nFeature
                                    data = zeros(nPatient,1);
                                    for p = 1:nPatient
                                        data(p,1) = tempText{p}.(strExperiment).(nameType{t}).(nameFeature{f});
                                    end
                                    text{i,w,r,s,a,n}.(nameType{t}).(nameFeature{f}).Data = data;
                                    for o = 1:nOutcomes
                                        text{i,w,r,s,a,n}.(nameType{t}).(nameFeature{f}).Spearman.(nameOutcomes{o}) = corr(data,outcomes.(nameOutcomes{o}),'type','Spearman');
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    cd(pathWORK), save(['textures_',scans{scan}],'text')
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function makelegend (offset, namedisplay, xlabel_my, ylabel_my)
% function makelegend (offset, namedisplay, xlabel_my, ylabel_my)
% sep = .1:.1:2;
% offset = [10 100 1000];
kk = 1;
% t{1} = namedisplay;

% t{2} = 'ICA';

for ii=1:1
    for jj= 1:length(offset)
%         l{kk} = [t{ii} 'offset ' num2str(offset(jj))];
%         l{kk} = [t{ii} num2str(offset(jj,:))];
%         l{kk}=namedisplay{jj};
        l{kk}=[namedisplay  num2str(offset(jj))];
        kk = kk+1;
    end
end

legend(l)
xlabel (xlabel_my)
ylabel (ylabel_my)
% xlabel('separation [pixels]')
% ylabel('loc err [pixels]')
grid on
% xlim([0 10])
% ylabel('loc err / separation')

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text: %SUBSASGN Dataset overload

author: marianux

repo: ecg-kit

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


text:
function surf = specificSurface(img, varargin)
%SPECIFICSURFACE implementation of Ohser's algo for surface comput.
%
%   compute surface area in discrete images.
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 16/02/2005.
%

%   HISTORY

img = img~=0;

delta = [1 1 1];

% edges of the unit cell
kl = [1 2;1 3;1 5;1 4;2 3;1 6;2 5;1 7;3 5;1 8;2 7;3 6;4 5];

% unit surface for each cell
c1 = 0.045778;
c2 = 0.036981;
c3 = 0.035196;
c = [c1 c1 c1   c2 c2 c2   c2 c2 c2   c3 c3 c3 c3];

% distances between pixels
d12 = sqrt(delta(1)*delta(1) + delta(2)*delta(2));
d13 = sqrt(delta(1)*delta(1) + delta(3)*delta(3));
d23 = sqrt(delta(2)*delta(2) + delta(3)*delta(3));
d123 = sqrt(sum(delta.*delta));
r = [delta(1) delta(2) delta(3) d12 d12 d13 d13 d23 d23 d123 d123 d123 d123];

% compute gray-tone histogram of the image
h = grayHist(img);

sv = 0;

% for each type of configuration
for l=1:256

    v = l-1;
    b(1) = bitand(v,1)~=0;
    b(2) = bitand(v,2)~=0;
    b(3) = bitand(v,4)~=0;
    b(4) = bitand(v,8)~=0;
    b(5) = bitand(v,16)~=0;
    b(6) = bitand(v,32)~=0;
    b(7) = bitand(v,64)~=0;
    b(8) = bitand(v,128)~=0;

    % for each edge of configuration
    for nu=1:13
        sv = sv + h(l)*c(nu)/r(nu)*xor( b(kl(nu, 1)), b(kl(nu,2)) );
    end
end

surf = 4*sv/sum(h(:));

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [X,Y,Z] = unplot(filename)
%%To unplot data from a matlab figure (.fig) files generated
% using version 7 or later. It can be used for both 2D and 3D plots
% Example usage:
% To unplot 2D graphs
% [x,y] = unplot('example2D.fig')
% To unplot 3D graphs
% [x,y,z] = unplot('example3D.fig')

% Pradyumna
% January 2012

if nargin==1
    fig1 = load (filename,'-mat');
    a = fig1.hgS_070000.children.children(1,1).properties;
    if isfield(a,'ZData')
        % unploting 3D plot
        Y = fig1.hgS_070000.children.children(1,1).properties.YData;
        X = fig1.hgS_070000.children.children(1,1).properties.XData;
        Z = fig1.hgS_070000.children.children(1,1).properties.ZData;
    else
        % unploting 2D plot
        Y = fig1.hgS_070000.children.children(1,1).properties.YData;
        X = fig1.hgS_070000.children.children(1,1).properties.XData;
    end
else
    disp('Usage unplot(''filename.fig''). See ''help unplot'' for more details');
end


text:
function determ = ris_determinant ( n )

%*****************************************************************************80
%
%% RIS_DETERMINANT returns the determinant of the Ris matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real DETERM the determinant.
%
  top = 1.0;
  for i = 1 : n
    for j = i + 1 : n
      top = top * ( 4 * ( i - j ) * ( i - j ) );
    end
  end

  bottom = 1.0;
  for i = 1 : n
    for j = 1 : n
      bottom = bottom * ( 3 + 2 * n - 2 * i - 2 * j );
    end
  end

  determ = top / bottom;

  return
end


text:
function  MENSRB  = readMENSRB(fid)
%
% <<<<<<<<<< CLASSIFICATION: UNCLASSIFIED >>>>>>>>>>
%
% READMENSRB Create structure for metadata fields.
%   MENSRB = READMENSRB(FID) returns a structure of fields and their
%   associated values.

% Written by Matt Donath, NGA, matthew.b.donath@nga.ic.gov
% References:
%  > MIL-STD-2500A, NATIONAL IMAGERY TRANSMISSION FORMAT VERSION 2.0
%  > MIL-STD-2500C, NATIONAL IMAGERY TRANSMISSION FORMAT VERSION 2.1
%  > STDI-0001,     NATIONAL SUPPORT DATA EXTENSIONS (SDE)(VERSION 1.3/CN2)
%                   FOR THE NATIONAL IMAGERY TRANSMISSION FORMAT (NITF)
%  > STDI-0002,     THE COMPENDIUM OF CONTROLLED EXTENSIONS (CE) FOR THE
%                   NATIONAL IMAGERY TRANSMISSION FORMAT (NITF) VERSION 2.1
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Create the metadata structure for MENSRB.

MENSRB.CETAG = fread(fid,6,'uint8=>char')';
MENSRB.CEL = str2double(fread(fid,5,'uint8=>char')');

MENSRB.ACFT_LOC = fread(fid,25,'uint8=>char')';
MENSRB.ACFT_LOC_ACCY = str2double(fread(fid,6,'uint8=>char')');
MENSRB.ACFT_ALT = str2double(fread(fid,6,'uint8=>char')');
MENSRB.RP_LOC = fread(fid,25,'uint8=>char')';
MENSRB.RP_LOC_ACCY = str2double(fread(fid,6,'uint8=>char')');
MENSRB.RP_ELV = str2double(fread(fid,6,'uint8=>char')');
MENSRB.OF_PC_R = str2double(fread(fid,7,'uint8=>char')');
MENSRB.OF_PC_A = str2double(fread(fid,7,'uint8=>char')');
MENSRB.COSGRZ = str2double(fread(fid,7,'uint8=>char')');
MENSRB.RGCCRP = str2double(fread(fid,7,'uint8=>char')');
MENSRB.RLMAP = fread(fid,1,'uint8=>char')';
MENSRB.RP_ROW = str2double(fread(fid,5,'uint8=>char')');
MENSRB.RP_COL = str2double(fread(fid,5,'uint8=>char')');
MENSRB.C_R_NC = str2double(fread(fid,10,'uint8=>char')');
MENSRB.C_R_EC = str2double(fread(fid,10,'uint8=>char')');
MENSRB.C_R_DC = str2double(fread(fid,10,'uint8=>char')');
MENSRB.C_AZ_NC = str2double(fread(fid,9,'uint8=>char')');
MENSRB.C_AZ_EC = str2double(fread(fid,9,'uint8=>char')');
MENSRB.C_AZ_DC = str2double(fread(fid,9,'uint8=>char')');
MENSRB.C_AL_NC = str2double(fread(fid,9,'uint8=>char')');
MENSRB.C_AL_EC = str2double(fread(fid,9,'uint8=>char')');
MENSRB.C_AL_DC = str2double(fread(fid,9,'uint8=>char')');
MENSRB.TOTAL_TILES_COLS = str2double(fread(fid,3,'uint8=>char')');
MENSRB.TOTAL_TILES_ROWS = str2double(fread(fid,5,'uint8=>char')');

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


text:
function g = mdngrad(net, x, t)
%MDNGRAD Evaluate gradient of error function for Mixture Density Network.
%
%	Description
%	 G = MDNGRAD(NET, X, T) takes a mixture density network data
%	structure NET, a matrix X of input vectors and a matrix T of target
%	vectors, and evaluates the gradient G of the error function with
%	respect to the network weights. The error function is negative log
%	likelihood of the target data.  Each row of X corresponds to one
%	input vector and each row of T corresponds to one target vector.
%
%	See also
%	MDN, MDNFWD, MDNERR, MDNPROB, MLPBKP
%

%	Copyright (c) Ian T Nabney (1996-2001)
%	David J Evans (1998)

% Check arguments for consistency
errstring = consist(net, 'mdn', x, t);
if ~isempty(errstring)
  error(errstring);
end

[mixparams, y, z] = mdnfwd(net, x);

% Compute gradients at MLP outputs: put the answer in deltas
ncentres = net.mdnmixes.ncentres;
dim_target = net.mdnmixes.dim_target;
nmixparams = net.mdnmixes.nparams;
ntarget = size(t, 1);
deltas = zeros(ntarget, net.mlp.nout);
e = ones(ncentres, 1);
f = ones(1, dim_target);

post = mdnpost(mixparams, t);

% Calculate prior derivatives
deltas(:,1:ncentres)  = mixparams.mixcoeffs - post;

% Calculate centre derivatives
long_t = kron(ones(1, ncentres), t);
centre_err = mixparams.centres - long_t;

% Get the post to match each u_jk:
% this array will be (ntarget, (ncentres*dim_target))
long_post = kron(ones(dim_target, 1), post);
long_post = reshape(long_post, ntarget, (ncentres*dim_target));

% Get the variance to match each u_jk:
var = mixparams.covars;
var = kron(ones(dim_target, 1), var);
var = reshape(var, ntarget, (ncentres*dim_target));

% Compute centre deltas
deltas(:, (ncentres+1):(ncentres*(1+dim_target))) = ...
                       (centre_err.*long_post)./var;

% Compute variance deltas
dist2             = mdndist2(mixparams, t);
c                 = dim_target*ones(ntarget, ncentres);
deltas(:, (ncentres*(1+dim_target)+1):nmixparams) = ...
                      post.*((dist2./mixparams.covars)-c)./(-2);

% Now back-propagate deltas through MLP
g = mlpbkp(net.mlp, x, z, deltas);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function varargout = createDodecahedron()
%CREATEDODECAHEDRON Create a 3D mesh representing a dodecahedron.
%
%   [V, E, F] = createDodecahedron;
%   Create a 3D mesh representing a dodecahedron
%   V is the 20-by-3 array of vertex coordinates
%   E is the 30-by-2 array of edge vertex indices
%   F is the 12-by-5 array of face vertex indices
%
%   [V, F] = createDodecahedron;
%   Returns only the vertices and the face vertex indices.
%
%   MESH = createDodecahedron;
%   Returns the data as a mesh structure, with fields 'vertices', 'edges'
%   and 'faces'.
%
%   Example
%   [v, e, f] = createDodecahedron;
%   drawMesh(v, f);
%
%   Use values given by P. Bourke, see:
%   http://local.wasp.uwa.edu.au/~pbourke/geometry/platonic/
%   faces are re-oriented to have normals pointing outwards.
%
%   See also
%   meshes3d, drawMesh
%   createCube, createOctahedron, createIcosahedron, createTetrahedron
%

%   ---------
%   author : David Legland
%   e-mail: david.legland@inra.fr
%   INRA - TPV URPOI - BIA IMASTE
%   created the 29/07/2010.
%

%   HISTORY

% golden ratio
phi = (1+sqrt(5))/2;

% coordinates pre-computations
b = 1 / phi ;
c = 2 - phi ;

% use values given by P. Bourke, see:
% http://local.wasp.uwa.edu.au/~pbourke/geometry/platonic/
tmp = [ ...
 c  0  1 ;   b  b  b ;   0  1  c  ; -b  b  b  ; -c  0  1 ;  ...
-c  0  1 ;  -b -b  b ;   0 -1  c  ;  b -b  b  ;  c  0  1 ;   ...
 c  0 -1 ;   b -b -b ;   0 -1 -c  ; -b -b -b  ; -c  0 -1 ;  ...
-c  0 -1 ;  -b  b -b ;   0  1 -c  ;  b  b -b  ;  c  0 -1 ; ...
 0  1 -c ;   0  1  c ;   b  b  b  ;  1  c  0  ;  b  b -b ; ...
 0  1  c ;   0  1 -c ;  -b  b -b  ; -1  c  0  ; -b  b  b ; ...
 0 -1 -c ;   0 -1  c ;  -b -b  b  ; -1 -c  0  ; -b -b -b ; ...
 0 -1  c ;   0 -1 -c ;   b -b -b  ;  1 -c  0  ;  b -b  b ; ...
 1  c  0 ;   b  b  b ;   c  0  1  ;  b -b  b  ;  1 -c  0 ;  ...
 1 -c  0 ;   b -b -b ;   c  0 -1  ;  b  b -b  ;  1  c  0 ; ...
-1  c  0 ;  -b  b -b ;  -c  0 -1  ; -b -b -b  ; -1 -c  0 ; ...
-1 -c  0 ;  -b -b  b ;  -c  0  1  ; -b  b  b  ; -1  c  0 ;  ...
];

% extract coordinates of unique vertices
[verts, M, N] = unique(tmp, 'rows', 'first'); %#ok<ASGLU>

% compute indices of face vertices, put result in a 12-by-5 index array
ind0 = reshape((1:60), [5 12])';
faces = N(ind0);

% extract edges from faces
edges = [reshape(faces(:, 1:5), [60 1]) reshape(faces(:, [2:5 1]), [60 1])];
edges = unique(sort(edges, 2), 'rows');

% format output
varargout = formatMeshOutput(nargout, verts, edges, faces);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%BASE64DECODE decodes a base64 string as a vector of uint8s
% BASE64DECODE(BASE64)
%    Decodes every four BASE64 characters to three bytes.
%    If BASE64 is padded with '=', these are translated to zeros, and
%    stripped from the resulting bytes.

% (c) 2014 Bastian Bechtold
% This code is licensed under the BSD 3-clause license

function bytes = base64decode(base64)
    % strip line breaks
    base64 = strrep(base64, sprintf('\n'), '');
    % add padding if missing
    if mod(length(base64), 4) ~= 0
        base64 = [base64 repmat('=', [1 mod(length(base64), 4)])];
    end
    % remember padding
    padding = sum(base64 == '=');
    % convert from string representation to base64 bytes
    table = zeros(1, 128, 'uint8');
    table(['A':'Z' 'a':'z' '0':'9' '+' '/']+0) = 0:63;
    base64 = table(base64');

    % convert every four base64 bytes to three uint8 bytes
    bytes = zeros(length(base64)/4*3, 1, 'uint8');
    bytes(1:3:end) = bitor(bitshift(base64(1:4:end), 2), ...
                           bitshift(base64(2:4:end), -4));
    bytes(2:3:end) = bitor(bitshift(bitand(base64(2:4:end), 15), 4), ... % four LSB
                           bitshift(base64(3:4:end), -2));
    bytes(3:3:end) = bitor(bitshift(bitand(base64(3:4:end), 3), 6), ... % two LSB
                           base64(4:4:end));
    % strip padding
    bytes = bytes(1:end-padding);
end

author: bastibe

repo: transplant

save_path: github-repos/MATLAB/bastibe-transplant


text:
function fx = p20_fun ( option, nvar, x )

%*****************************************************************************80
%
%% P20_FUN evaluates the function for problem 20.
%
%  Title:
%
%    The Buckling Spring
%
%  Description:
%
%    We are given three points A, B, and C.
%    A is at the origin (0,0).
%    B has coordinates (X,Y) with Y nonnegative, and the ray from A to B
%    makes an angle of THETA with the horizontal axis.
%    C is at the point (2*X,0).
%
%    A spring extends from A to B, and is normally of length 1,
%    and is currently of length L.
%    A spring extends from B to C, and is normally of length 1,
%    and is currently of length L.
%    A spring force is also exerted, which tends to draw the two
%    springs together, proportional to the angle between the two springs.
%
%    A vertical load MU is applied at point B (downward is positive).
%    A horizontal load is applied at point C (leftware is positive).
%    The spring force is applied perpendicularly to the axes of the two springs.
%
%    If we compute F(1), the force along the axis of one spring, and
%    F(2), the force perpendicular to the axis of one spring, we have that
%    F(L,THETA,LAMBDA,MU) is given by:
%
%    F(1) =  - 2 ( 1 - L ) + 2 * LAMBDA * cos ( THETA ) + MU * sin ( THETA )
%    F(2) =  0.5 * THETA - 2 * LAMBDA * L * sin ( THETA ) + MU * L * cos ( THETA )
%
%    The user must specify a third, augmenting equation, of the form
%
%    F(3) = X(HOLD_INDEX) - HOLD_VALUE.
%
%    Typically, HOLD_INDEX is 2, for the varlable THETA, and HOLD_VALUE is
%    an angle measured in radians and strictly between 0 and PI/2.
%
%    Another choice for HOLD_INDEX would be 1, for the variable L, with
%    HOLD_VALUE greater than 0.  Values of L less than 1 represent compressed
%    springs; values greater than 1 indicate extended springs.
%
%    L represents the length of a spring.
%    THETA represents the angle that the spring makes with the horizontal axis.
%    MU is a vertical load applied at the midpoint B.
%    LAMBDA is a horizontal load applied at the right endpoint C.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Tim Poston, Ian Stewart,
%    Catastrophe Theory and its Applications,
%    Dover, 1996,
%    ISBN13: 978-0486692715,
%    LC: QA614.58.P66.
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Input, real X(NVAR), the argument of the function.
%
%    Output, real FX(NVAR-1), the value of the function at X.
%
  l      = x(1);
  theta  = x(2);
  lambda = x(3);
  mu     = x(4);

  fx = zeros ( nvar - 1, 1 );

  fx(1) = - 2.0 * ( 1.0 - l ) ...
          + 2.0 * lambda * cos ( theta ) ...
          + mu * sin ( theta );

  fx(2) = 0.5 * theta ...
          - 2.0 * lambda * l * sin ( theta ) ...
          + mu * l * cos ( theta );

  hold_index = p20_i4_get ( 'hold_index' );
  hold_value = p20_r8_get ( 'hold_value' );

  fx(3) = x(hold_index) - hold_value;

  return
end


text:
  function data = load_ascii_skip_header(file, varargin)
%|function data = load_ascii_skip_header(file)
%| read ascii file, skipping lines that start with '#'
%| output is double array of numeric values
%| option
%|	'char'	0|1	if 1 then output cell array of each line
%| Copyright 2007, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
arg.char = false;
arg = vararg_pair(arg, varargin);

fid = fopen(file, 'r');
if fid == -1
	fail('problem opening "%s"', file)
end

if arg.char
	data = {};
else
	data = [];
end

while 1
	line = fgetl(fid);
	if ~ischar(line)
		if isempty(data)
			fail('file %s ended too soon', file)
		else
			break % go to fclose
		end
	end
	if ~length(line) || line(1) == '#'
		continue
	end
	if arg.char
		data{end+1} = line;
	else
		data(end+1,:) = sscanf(line, '%f');
	end
end

tmp = fclose(fid);
if tmp ~= 0
	fail('problem with fclose for "%s"', file)
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% Test file for @chebfun/feval.m

function pass = test_feval(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebfunpref();
end

% Generate a few random points in [-1 1] to use as test values.
seedRNG(7681);
xr = 2 * rand(1000, 1) - 1;

% Test empty input.
f = chebfun(@(x) x, pref);
fx = feval(f, zeros(0, 4));
pass(1) = isequal(size(fx), [0 4]);

% Test endpoint evaluation.
pref.splitting = 0;
f = chebfun(@(x) erf(x), [-1 1], pref);

lval1 = feval(f, 'left');
lval2 = feval(f, 'start');
lval3 = feval(f, '-');
pass(2) = (abs(lval1 - erf(-1)) < 10*eps*vscale(f) ) ...
    && isequal(lval1, lval2, lval3);

rval1 = feval(f, 'right');
rval2 = feval(f, 'end');
rval3 = feval(f, '+');
pass(3) = (abs(rval1 - erf(1)) < 10*eps*vscale(f) ) ...
    && isequal(rval1, rval2, rval3);

% Test bogus input for endpoint specification.
try
    badval = feval(f, '???');
    pass(4) = false;
catch ME
    pass(4) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN:feval:strInput');
end

% Test evaluation with and without left/right limit values.
pref.splitting = 1;
f = chebfun(@(x) sign(x), [-1 1], pref);
tol = 10*eps*vscale(f);

x = [-0.5 ; 0 ; 0.5];
err = norm(feval(f,x) - [-1; 0; 1], inf);
pass(5) = err < tol;

leftlim1 = feval(f, x, 'left');
leftlim2 = feval(f, x, '-');
err1 = norm(leftlim1 - [-1; -1; 1], inf);
err2 = norm(leftlim1 - leftlim2, inf);
pass(6) = err1 < tol && err2 < tol;

rightlim1 = feval(f, x, 'right');
rightlim2 = feval(f, x, '+');
err1 = norm(rightlim1 - [-1; 1; 1], inf);
err2 = norm(rightlim1 - rightlim2, inf);
pass(7) = err1 < tol && err2 < tol;

% Test bogus input for limit specification string.
try
    badval = feval(f, x, '???');
    pass(8) = false;
catch ME
    pass(8) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN:feval:leftRightChar');
end

try
    badval = feval(f, x, []);
    pass(9) = false;
catch ME
    pass(9) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN:feval:leftRight');
end

% Spot-check a few functions.
pref.splitting = 0;
f_exact = @(x) exp(x) - 1;
f = chebfun(f_exact, [-1 1], pref);
x = xr;
pass(10) = (norm(feval(f, x) - f_exact(x), inf) < 10*eps*vscale(f));

pref.splitting = 1;
f_exact = @(x) abs(x)./(1 + x.^2);
f = chebfun(f_exact, [-2 7], pref);
x = 4.5*xr + 2.5;
pass(11) = (norm(feval(f, x) - f_exact(x), inf) < 1e2*eps*vscale(f));

pref.splitting = 0;
f_exact = @(x) cos(1e4*x);
f = chebfun(f_exact, [1 5], pref);
x = 2*xr + 3;
pass(12) = (norm(feval(f, x) - f_exact(x), inf) < 1e5*eps*vscale(f));

pref.splitting = 0;
z = exp(2*pi*1i/6);
f_exact = @(t) sinh(t*z);
f = chebfun(f_exact, [-1 1], pref);
x = xr;
pass(13) = (norm(feval(f, x) - f_exact(x), inf) < 10*eps*vscale(f));

% Check row vector and matrix input.
pref.splitting = 0;
f_exact = @(x) cos(x - 0.2);
f = chebfun(f_exact, [-1 1], pref);

x = xr;
err = feval(f, x.') - f_exact(x.');
pass(14) = isequal(size(err), [1 1000]) && (norm(err(:), inf) < ...
    10*eps*vscale(f));

x_mtx = reshape(x, [100 10]);
err = feval(f, x_mtx) - f_exact(x_mtx);
pass(15) = isequal(size(err), [100 10]) && (norm(err(:), inf) < ...
    10*eps*vscale(f));

x_3mtx = reshape(x, [5 20 10]);
err = feval(f, x_3mtx) - f_exact(x_3mtx);
pass(16) = isequal(size(err), [5 20 10]) && (norm(err(:), inf) < ...
    10*eps*vscale(f));

x_4mtx = reshape(x, [5 4 5 10]);
err = feval(f, x_4mtx) - f_exact(x_4mtx);
pass(17) = isequal(size(err), [5 4 5 10]) && (norm(err(:), inf) < ...
    10*eps*vscale(f));

% Check behavior for transposed chebfuns.
f = transpose(f);

err = feval(f, x_mtx) - f_exact(x_mtx).';
pass(18) = isequal(size(err), [10 100]) && (norm(err(:), inf) < ...
    10*eps*vscale(f));

err = feval(f, x_3mtx) - permute(f_exact(x_3mtx), [2 1 3]);
pass(19) = isequal(size(err), [20 5 10]) && (norm(err(:), inf) < ...
    10*eps*vscale(f));

f = transpose(f);

% Check evaluation at points just outside the domain.
x = [(-1 - 1e-6) ; 1e-6 + 1e-6i ; (1 + 1e-6)];
err = feval(f, x.') - f_exact(x.');
pass(20) = norm(err(:), inf) < 10*eps*vscale(f);

% Check operation for array-valued chebfuns.
pref.splitting = 1;
f_exact = @(x) [sin(x) abs(x) exp(1i*x)];
f = chebfun(f_exact, [], pref);

x = xr;
err = feval(f, x) - f_exact(x);
pass(21) = all(max(abs(err)) < 1e2*eps*vscale(f));

x_mtx = reshape(x, [100 10]);
err = feval(f, x_mtx) - f_exact(x_mtx);
pass(22) = isequal(size(err), [100 30]) && (norm(err(:), inf) < ...
    1e2*eps*vscale(f));

x_3mtx = reshape(x, [5 20 10]);
err = feval(f, x_3mtx) - f_exact(x_3mtx);
pass(23) = isequal(size(err), [5 60 10]) && (norm(err(:), inf) < ...
    1e2*eps*vscale(f));

x_4mtx = reshape(x, [5 4 5 10]);
err = feval(f, x_4mtx) - f_exact(x_4mtx);
pass(24) = isequal(size(err), [5 12 5 10]) && (norm(err(:), inf) < ...
    1e2*eps*vscale(f));

f.isTransposed = 1;
err = feval(f, x_4mtx) - permute(f_exact(x_4mtx), [2 1 3 4]);
pass(25) = isequal(size(err), [12 5 5 10]) && (norm(err(:), inf) < ...
    1e2*eps*vscale(f));

err = feval(f, 'left') - [sin(-1) 1 exp(-1i)].';
pass(26) = norm(err(:), inf) < 10*eps*vscale(f);

err = feval(f, 'right') - [sin(1) 1 exp(1i)].';
pass(27) = norm(err(:), inf) < 10*eps*vscale(f);

% Test feval(chebfun, location, 'left'/'right'):
f = chebfun({@(x) [x, 2*x], @(x) [3*x, 4*x]}, [0 1 2]);
pass(28) = norm(feval(f, 1, 'left') - [1, 2], inf) < 10*eps*vscale(f);
pass(29) = norm(feval(f, 1, 'right') - [3, 4], inf) < 10*eps*vscale(f);

% Check for dimension mismatch bug when evaluating an array-valued chebfun on a
% vector which contains breakpoint values, repeated several times.
f_exact = @(x) [sin(x) cos(x)];
f = chebfun(f_exact, [-1 1]);
x = [-1 ; 0.5 ; 1];
err = feval(f, x) - f_exact(x);
pass(30) = all(max(abs(err)) < 10*eps*vscale(f));

%% Test on singular function:

dom = [-2 7];

% Generate a few random points to use as test values.
seedRNG(6178);
x = diff(dom) * rand(100, 1) + dom(1);

pow = -0.5;
op = @(x) (x - dom(1)).^pow.*sin(200*x);
f = chebfun(op, dom, 'exps', [pow 0], 'splitting', 'on');
fval = feval(f, x);
vals_exact = feval(op, x);
err = fval - vals_exact;
pass(31) = ( norm(err, inf) < 1e5*eps*norm(vals_exact, inf) );

%% Test for function defined on unbounded domain:

% Functions on [-inf inf]:

% Set the domain:
dom = [-Inf Inf];
domCheck = [-1e2 1e2];

% Generate a few random points to use as test values:
x = diff(domCheck) * rand(100, 1) + domCheck(1);

op = @(x) (1-exp(-x.^2))./x;
f = chebfun(op, dom);
fVals = feval(f, x);
fExact = op(x);
err = fVals - fExact;
pass(32) = ( norm(err, inf) < 1e1*eps*vscale(f) ) && ...
    ( feval(f, Inf) < eps*vscale(f) ) && ...
    ( feval(f, -Inf) < eps*vscale(f) );

% Blow-up function:
op = @(x) x.^2.*(1-exp(-x.^2));
f = chebfun(op, dom, 'exps', [2 2]);
fVals = feval(f, x);
fExact = op(x);
err = fVals - fExact;
pass(33) = ( norm(err, inf) < 1e5*eps*vscale(f) ) && ...
    ( feval(f, Inf) == Inf ) && ( feval(f, -Inf) == Inf );

%% Functions on [a inf]:

% Set the domain:
dom = [1 Inf];
domCheck = [1 1e2];

% Generate a few random points to use as test values:
x = diff(domCheck) * rand(100, 1) + domCheck(1);

op = @(x) x.*exp(-x);
f = chebfun(op, dom);
fVals = feval(f, x);
fExact = op(x);
err = fVals - fExact;
pass(34) = ( norm(err, inf) < 1e1*eps*vscale(f) ) && ...
     ( feval(f, Inf) < 1e2*eps*vscale(f) );

%% Functions on [-inf b]:

% Set the domain:
dom = [-Inf -3*pi];
domCheck = [-1e6 -3*pi];

% Generate a few random points to use as test values:
x = diff(domCheck) * rand(100, 1) + domCheck(1);

% Array-valued function:
op = @(x) [exp(x) x.*exp(x) (1-exp(x))./x];
f = chebfun(op, dom);
fVals = feval(f, x);
fExact = op(x);
err = fVals - fExact;
pass(35) = ( norm(err, inf) < 1e2*max(eps*vscale(f)) ) && ...
    all( feval(f, -Inf) < 1e2*eps*vscale(f) );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun