text:
classdef REMO < ALGORITHM
% <multi/many> <real> <expensive>
% Expensive multiobjective optimization by relation learning and prediction
% k    ---    6 --- Number of reference solutions
% gmax --- 3000 --- Number of solutions evaluated by surrogate model

%------------------------------- Reference --------------------------------
% H. Hao, A. Zhou, H. Qian, and H. Zhang, Expensive multiobjective
% optimization by relation learning and prediction, IEEE Transactions on
% Evolutionary Computation, 2022.
%------------------------------- 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)
            %% Parameterr setting
            [k,gmax] = Algorithm.ParameterSet(6,3000);

            %% Initalize the population by Latin hypercube sampling
            if Problem.D <= 10
                N = 11*Problem.D-1;
            else
                N = 100;
            end
            PopDec     = UniformPoint(N,Problem.D,'Latin');
            Population = Problem.Evaluation(repmat(Problem.upper-Problem.lower,N,1).*PopDec+repmat(Problem.lower,N,1));
            Archive    = Population;

            %% Optimization
            while Algorithm.NotTerminated(Archive)
                % Select reference solutions and preprocess the data
                Ref       = RefSelect(Population,k);
                Input     = Population.decs;
                Catalog   = GetOutput_PBI(Population.objs,Ref.objs);
                [XXs,YYs] = GetRelationPairs(Input,Catalog);
                [TrainIn,TrainOut,TestIn,TestOut] = DataProcess(XXs,YYs);
                xDim = size(TrainIn,2);

                % Train relation model
                [TrainIn_nor,TrainIn_struct] = mapminmax(TrainIn');
                TrainIn_nor     = TrainIn_nor';
                TrainOut_onehot = onehotconv(TrainOut,1);
                net = patternnet([ceil(xDim*1.5),xDim*1,ceil(xDim/2)]);
                net.trainParam.showWindow =0;
                net        = train(net,TrainIn_nor',TrainOut_onehot');
                TestIn_nor = mapminmax('apply',TestIn',TrainIn_struct)';
                TestPre    = onehotconv(net(TestIn_nor')',2);
                p_err      = sum(TestPre ~= TestOut)/size(TestPre,1);
                Smodel.X   = Input;
                Smodel.Y   = Catalog;
                Smodel.mp_struct = TrainIn_struct;
                Smodel.net       = net;
                Smodel.p_err     = p_err;
                Next = RSurrogateAssistedSelection(Problem,Ref,Population.decs,gmax,Smodel);
                if ~isempty(Next)
                    Archive = [Archive,Problem.Evaluation(Next)];
                end
                Population = RefSelect(Archive,Problem.N);
            end
        end
	end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
classdef quadcopter < agent
    properties
        % GLOBAL.position   - The position of the object in the global axes
        % GLOBAL.velocity   - The velocity of the object in the global axes
        % GLOBAL.quaternion - The quaternion representing the earth to rotated body

        % DYNAMICS - All the models parameters are held in the DYNAMICS
        %            container field.
    end
    %% ////////////////////////// MAIN METHODS /////////////////////////////
    methods
        % Constructor
        function [this] = quadcopter(varargin)
            % Call the super class
            this@agent(varargin);

            % Create the dynamics of a quadcopter
            this.DYNAMICS = this.CreateDYNAMICS();

            % //////////////// Check for user overrides ///////////////////
            this = this.ApplyUserOverrides(varargin); % Recursive overrides
            % /////////////////////////////////////////////////////////////
        end
        % Setup
        function [this] = setup(this,localXYZVelocity,localXYZrotations)
            % This function calculates the intial state for a quadrotor
            % object.

            % For a state vector defined as x_k:
            % [x,y,z,dx,dy,dz,R11,R12,R13,R21,R22,R23,R31,R32,R33,wx,wy,wz]

            % THIS MODEL OPERATES IN THE GLOBAL FRAME
            p0 = this.GetGLOBAL('position');   % True global position
            v0 = this.GetGLOBAL('velocity');   % True global velocity
            % GET THE ROTATION MATRIX fixed local axes -> rotated global pose
            R0 = OMAS_geometry.eulersToRotationMatrix(localXYZrotations);
            % Build an initial (global) state vector
            vecR0   = reshape(R0',9,1);        % Defines local vector to global vector
            omega0  = zeros(3,1);
            x0      = [p0;v0;vecR0;omega0];
            % ASSIGN THE LOCAL FRD STATE
            this.SetGLOBAL('priorState',x0);
            this.localState = x0;
        end
        % Main
        function [this] = main(this,ENV,varargin)

            % //////////// CHECK FOR NEW INFORMATION UPDATE ///////////////
            % Update the agent with the environmental data
            [this,~,~] = this.GetAgentUpdate(ENV,varargin{1});
            % /////////////////////////////////////////////////////////////

            % The state reference
            desiredPosition = [0;0;0];

            % Call the controller loop
            this = this.Controller(ENV,desiredPosition);

            % /////////////// RECORD THE AGENT-SIDE DATA //////////////////
            this.DATA.inputNames = {'$\dot{x}$ (m/s)','$\dot{y}$ (m/s)','$\dot{z}$ (m/s)',...
                '$p$ (rad/s)','$q$ (rad/s)','$r$ (rad/s)'};
            this.DATA.inputs(1:length(this.DATA.inputNames),ENV.currentStep) = [this.localState(4:6);this.localState(16:end)];          % Record the control inputs
        end
    end
    %% //////////////////////// AUXILLARY METHODS //////////////////////////
    methods
        % The quadcopter controller
        function [this] = Controller(this,ENV,y_k)

            % Input sanity checks
            assert(IsColumn(y_k,3),"Expecting a numeric velocity vector.");
            assert(1 == size(this.localState,2),"The length of the objects state update must match the its local state.");

            % ///////////// UPDATE THE LOCAL STATE VECTOR /////////////////
            x_k_plus = this.UpdateLocalState(ENV,this.localState,y_k);
            % Update the global properties
            this = this.GlobalUpdate(x_k_plus);
        end
        % Get the state update (using ode45)
        function [x_k_plus] = UpdateLocalState(this,ENV,x_k,y_desired)
            % This function computes the state update for the current agent
            % using the ode45 function.

            % Integrate across the time delta "dt"
            opts = odeset('RelTol',1e-2,'AbsTol',ENV.dt*1E-1);

            [~,Xset] = ode45(@(t,X) this.ClosedLoopDynamics_position(X,y_desired),[0 ENV.dt],x_k,opts);
%             [~,Xset] = ode45(@(t,X) this.ClosedLoopDynamics_velocity(X,y_desired),[0 ENV.dt],x_k,opts);

            % Assign the integral state
            x_k_plus = Xset(end,:)';
        end
        % Global update from the new state
        function [this] = GlobalUpdate(this,x_k_plus)
            % This function updates the global structure from the new state
            % definition: [p_k;v_k;vecR_k;omega_k]

            % Extract the rotation matrix components
            R_k_plus = reshape(x_k_plus(7:15),3,3);
            % Define the new quaternion from R
            q_k_plus = OMAS_geometry.rotationMatrixToQuaternion(R_k_plus');

            % //////////////// UPDATE GLOBAL PROPERTIES ///////////////////
            [this] = this.GlobalUpdate_direct(...
                x_k_plus(1:3),...   % The global position is within the state
                x_k_plus(4:6),...   % The global velocity is within the state
                q_k_plus);          % Append the global quaternion

            % Ensure the local state is re-assigned
            this.localState = x_k_plus;
        end
    end
    % Dynamic functions
    methods
        % Quadcopter Dynamics (Closed-loop)
        function [dxdt] = ClosedLoopDynamics_velocity(this,x_k,v_desired)

            % Extract the current state properties
            p_k     = x_k(1:3);
            v_k     = x_k(4:6);
            R_k     = reshape(x_k(7:15),3,3); % The rotation matrix components
            eta_k   = OMAS_geometry.rotationMatrixToEulers(R_k);
            omega_k = x_k(16:18);

            % The state reference
            y_k = [p_k;v_k;eta_k;omega_k];

            % Extract the parameters from the state
            psi_desired = pi/2;
            y_desired = [zeros(3,1);v_desired;[0;0;psi_desired];zeros(3,1)];

            % Extract DYNAMIC parameters
            g  = this.DYNAMICS.g;
            m  = this.DYNAMICS.m;

            % State matrix
            A = this.DYNAMICS.A;
            A(4:6,7:9) = g*[ sin(psi_desired), cos(psi_desired), 0;
                            -cos(psi_desired), sin(psi_desired), 0;
                                            0,                0, 0];
            % Input matrix
            B = this.DYNAMICS.B;

            % LQR control gain
            [K,~,~] = lqr(A,B,eye(12),eye(4));
            K(:,1:3) = 0;
            % Calculate the error
            e_k  = y_k - y_desired;
            du_k = -K*e_k;
            % Calculate the inputs
            u_k = [m*g;0;0;0] + du_k;

            % Provide the inputs to the open-loop dynamics
            [dxdt] = this.OpenLoopDynamics(x_k,u_k);
        end
        % Quadcopter Dynamics (Closed-loop)
        function [dxdt] = ClosedLoopDynamics_position(this,x_k,p_desired)

            % Extract the current state properties
            p_k     = x_k(1:3);
            v_k     = x_k(4:6);
            R_k     = reshape(x_k(7:15),3,3); % The rotation matrix components
            eta_k   = OMAS_geometry.rotationMatrixToEulers(R_k);
            omega_k = x_k(16:18);

            % The state reference
            y_k = [p_k;v_k;eta_k;omega_k];

            % Extract the parameters from the state
            psi_desired = 0;
            y_desired = [p_desired(1:3);zeros(3,1);[0;0;psi_desired];zeros(3,1)];

            % Extract DYNAMIC parameters
            e3 = this.DYNAMICS.e3;
            g  = this.DYNAMICS.g;
            m  = this.DYNAMICS.m;
            I  = this.DYNAMICS.I;

            % State matrix
            A = this.DYNAMICS.A;
            A(1:3,4:6) = eye(3);
            A(4:6,7:9) = g*[ sin(psi_desired), cos(psi_desired), 0;
                -cos(psi_desired), sin(psi_desired), 0;
                0,                0, 0];
            A(7:9,10:12)=eye(3);

            % Input matrix
            B = this.DYNAMICS.B;
            B(4:6,1) = e3/m;
            B(10:12,2:4) = inv(I);

            % LQR control gain
            [K,~,~] = lqr(A,B,eye(12),eye(4));
            % Calculate the error
            e_k  = y_k - y_desired;
            du_k = -K*e_k;
            % Calculate the inputs
            u_k = [m*g;0;0;0] + du_k;

            % Provide the inputs to the open-loop dynamics
            [dxdt] = this.OpenLoopDynamics(x_k,u_k);
        end
        % Quadcopter Dynamics (Open-loop)
        function [dxdt] = OpenLoopDynamics(this,x_k,u_k)

            % Extract the state parameters
            v_k = x_k(4:6);
            R_k = reshape(x_k(7:15),3,3);
            w_k = x_k(16:18);
            % Extract the inputs
            f   = u_k(1);
            tau = u_k(2:4);

            % Get the constants
            m  = this.DYNAMICS.m;
            g  = this.DYNAMICS.g;
            I  = this.DYNAMICS.I;
            e3 = this.DYNAMICS.e3;

            % nonlinear dynamics based on rotation dynamics
            p_dot   = v_k;
            v_dot   = f/m*R_k*e3-g*e3;
            R_dot    = R_k*skew(w_k);
            vecR_dot = reshape(R_dot,9,1);
            w_dot    = inv(I)*(tau-skew(w_k)*I*w_k);

            % THE STATE DIFFERENCE
            dxdt = [p_dot;v_dot;vecR_dot;w_dot];
        end
    end
    % Dynamics container
    methods
        % Get the (generic) quadcopter dynamic properties
        function [DYNAMICS] = CreateDYNAMICS(this)
            % Simple quadcopter, use the default properties
            DYNAMICS = this.CreateDYNAMICS_default();
        end
        % Create (quadcopter) default dynamic structure
        function [DYNAMICS] = CreateDYNAMICS_default(this)
            % Define an infinite "throw" range by default
            DYNAMICS = struct();
            % Define kinematic constraints
            DYNAMICS.maxLinearVelocity      = ones(3,1)*this.v_max;	% Limits on the agents linear velocity
            DYNAMICS.maxLinearAcceleration  = inf(3,1);            	% Limits on the agents linear acceleration
            DYNAMICS.maxAngularVelocity     = ones(3,1)*this.w_max;	% Limits on the agents angular velocity
            DYNAMICS.maxAngularAcceleration = inf(3,1);            	% Limits on the agents angular acceleration
            % Dynamical properties
            DYNAMICS.e3 = [0;0;1];                                   % Local z-axis
            DYNAMICS.g  = 300978377846937375/30680772461461504;      % Gravitational constant
            DYNAMICS.m  = 1;                                         % Quadcopter mass
            DYNAMICS.I  = 2*diag([4.856, 4.856, 8.801])*10^(-3);     % Quadcopter inertia
            % Control properties
            % Plant matrix
            DYNAMICS.A = zeros(12);
            DYNAMICS.A(4:6,7:9)   = eye(3)*DYNAMICS.g;
            DYNAMICS.A(1:3,4:6)   = eye(3);
            DYNAMICS.A(7:9,10:12) = eye(3);
            % Input matrix
            DYNAMICS.B = zeros(12,4);
            DYNAMICS.B(4:6,1)     = DYNAMICS.e3/DYNAMICS.m;
            DYNAMICS.B(10:12,2:4) = inv(DYNAMICS.I);
            % Observation matrix
            DYNAMICS.C = eye(12);
            % Feed-forward matrix
            DYNAMICS.D = eye(4);
        end
    end
end

repo: OpenMAS


text:
classdef KuwaharaFilter < EBSDFilter

  properties
    numNeighbours % number of neigbours to consider (default 1)
  end

  methods

    function F = KuwaharaFilter(varargin)
      %

      F.numNeighbours = get_option(varargin,'neighbours',1);

      addlistener(F,'isHex','PostSet',@check);
      function check(varargin)
        if F.isHex
          warning(['Hexagonal grids are not yet fully supportet for the KuwaharaFilter. ' ...
            'It might give reasonable results anyway']);
        end
      end

    end

    function ori = smooth(F,ori,quality)

      ori(quality==0) = nan;

      % map to mean
      [qmean,q] = mean(ori);
      q = reshape(inv(qmean)*q,size(ori)); %#ok<MINV>

      % prepare the result
      tqMean = nan(length(q),3);
      stdOpt = inf(size(q));

      n = F.numNeighbours;
      % make q a bit larger
      q = [quaternion.nan(n,size(q,2)+2*n);...
        [quaternion.nan(size(q,1),n),...
        q,quaternion.nan(size(q,1),n)];...
        quaternion.nan(n,size(q,2)+2*n)];

      % map quaternions into tangential space
      tq = double(log(q));

      for d = 0:3

        % decide for one of the quadrants
        dir = (1+1i)^(2*d+1);
        xdir = sign(real(dir));
        ydir = sign(imag(dir));

        % compute the mean
        meanLocal = zeros([size(ori),3]);
        count = zeros(size(meanLocal));

        for i = 0:n
          for j = 0:n
            [meanLocal,count] = nanplus(meanLocal, ...
              tq((1+n:end-n)+i*xdir,(1+n:end-n)+j*ydir,:),count);
          end
        end
        meanLocal = meanLocal ./ count;

        % compute variance
        stdLocal = zeros(size(ori));
        for i = 0:n
          for j = 0:n
            stdLocal = nanplus(stdLocal, ...
              sum((meanLocal-tq((1+n:end-n)+i*xdir,(1+n:end-n)+j*ydir,:)).^2,3));
          end
        end
        stdLocal = (stdLocal+eps) ./ max(0,(count(:,:,1)-1));

        % keep mean with smallest variance
        ind = stdLocal < stdOpt;
        ind3 = repmat(ind,[1,1,3]);
        tqMean(ind3) = meanLocal(ind3);
        stdOpt(ind) = stdLocal(ind);

      end

      % map back to orientation space
      q = quaternion(qmean) * reshape(expquat(tqMean),size(ori));
      ori.a = q.a; ori.b = q.b; ori.c = q.c; ori.d = q.d;

    end
  end
end

author: mtex-toolbox

repo: mtex

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


text:
function label = LSC(data,k,opts)
% label = LSC(data,k,opts): Landmark-based Spectral Clustering
% Input:
%       - data: the data matrix of size nSmp x nFea, where each row is a sample
%               point
%       - k: the number of clusters
%       opts: options for this algorithm
%           - p: the number of landmarks picked (default 1000)
%           - r: the number of nearest landmarks for representation (default 5)
%           - numRep: the number of replicates for the final kmeans (default 10)
%           - maxIter: the maximum number of iterations for final kmeans (default 100)
%           - mode: landmark selection method, currently support
%               - 'kmeans': use centers of clusters generated by kmeans (default)
%               - 'random': use randomly sampled points from the original
%                           data set
%           The following parameters are effective ONLY in mode 'kmeans'
%           - kmNumRep: the number of replicates for initial kmeans (default 1)
%           - kmMaxIter: the maximum number of iterations for initial kmeans (default 5)
% Output:
%       - label: the cluster assignment for each point
% Requre:
%       litekmeans.m
% Usage:
%       data = rand([100,50]);
%       label = LSC(data,10);
%Reference:
%
%	[1] Xinlei Chen, Deng Cai, "Large Scale Spectral Clustering with
%	Landmark-Based Representation," AAAI 2011.
%
%   [2] Deng Cai, Xinlei Chen: Large Scale Spectral Clustering Via
%   Landmark-Based Sparse Representation. IEEE Trans. Cybernetics 45(8):
%   1669-1680 (2015)
%
%   version 2.0 --Dec./2011
%   version 1.0 --Oct./2010
%
%   Written by Xinlei Chen (endernewton AT gmail.com)
%              Deng Cai (dengcai AT gmail.com)

% Set and parse parameters
if (~exist('opts','var'))
   opts = [];
end

p = 1000;
if isfield(opts,'p')
    p = opts.p;
end

r = 5;
if isfield(opts,'r')
    r = opts.r;
end

maxIter = 100;
if isfield(opts,'maxIter')
    maxIter = opts.maxIter;
end

numRep = 10;
if isfield(opts,'numRep')
    numRep = opts.numRep;
end

mode = 'kmeans';
if isfield(opts,'mode')
    mode = opts.mode;
end

nSmp=size(data,1);

% Landmark selection
if strcmp(mode,'kmeans')
    kmMaxIter = 5;
    if isfield(opts,'kmMaxIter')
        kmMaxIter = opts.kmMaxIter;
    end
    kmNumRep = 1;
    if isfield(opts,'kmNumRep')
        kmNumRep = opts.kmNumRep;
    end
    [~,marks]=litekmeans(data,p,'MaxIter',kmMaxIter,'Replicates',kmNumRep);
    clear kmMaxIter kmNumRep
elseif strcmp(mode,'random')
    indSmp = randperm(nSmp);
    marks = data(indSmp(1:p),:);
    clear indSmp
else
    error('mode does not support!');
end

% Z construction
D = EuDist2(data,marks);

if isfield(opts,'sigma')
    sigma = opts.sigma;
else
    sigma = mean(mean(D));
end

dump = zeros(nSmp,r);
idx = dump;
for i = 1:r
    [dump(:,i),idx(:,i)] = min(D,[],2);
    temp = (idx(:,i)-1)*nSmp+[1:nSmp]';
    D(temp) = 1e100;
end

dump = exp(-dump/(2*sigma^2));
sumD = sum(dump,2);
Gsdx = bsxfun(@rdivide,dump,sumD);
Gidx = repmat([1:nSmp]',1,r);
Gjdx = idx;
Z=sparse(Gidx(:),Gjdx(:),Gsdx(:),nSmp,p);

% Graph decomposition
feaSum = full(sqrt(sum(Z,1)));
feaSum = max(feaSum, 1e-12);
Z = Z./feaSum(ones(size(Z,1),1),:);
U = mySVD(Z,k+1);
U(:,1) = [];

U=U./repmat(sqrt(sum(U.^2,2)),1,k);

% Final kmeans
label=litekmeans(U,k,'MaxIter',maxIter,'Replicates',numRep);

author: ZJULearning

repo: MatlabFunc


text:
% Using Chebyshev Polynomials
clear all; close all; clc;

% Define:
f = @(x) x.^4;
%f = @(x) 4*(x.^2-x.^4).*exp(-x./2);

x = (-1:0.1:1)'; % for -1 <= x <= 1

repo: Aero-matlab


text:
function [z, logPrQ_z] = sampleStateSeq_PrecompSoftEv( ii, Etaii, LL, TargetPsi)

logPrQ_z = 0;
T = size( LL,2 );

if exist( 'TargetPsi', 'var' ) && ~isempty( TargetPsi )
    % Translate from TargetPsi's feature IDs to
    %    corresponding IDs in F, theta, prevStateSeq, etc.
    z = TargetPsi.stateSeq(ii).z;
    for aa = 1:length( TargetPsi.activeFeatIDs )
        z( z == TargetPsi.activeFeatIDs(aa) ) = TargetPsi.externalFeatIDs(aa);
    end
    TargetPsi.stateSeq(ii).z = z;
end

availFeatIDs = Etaii.availFeatIDs;
if length( availFeatIDs ) == 1
    z = availFeatIDs(1) * ones(1,T);
    return;
end
pi_z = Etaii.eta;
pi_z = bsxfun( @rdivide, pi_z, sum(pi_z,2)  );

pi_init = 1/length(pi_z)*ones( 1, length(pi_z) );

% Safely convert provided *log* soft evidence matrix
%  into proper probabilities
%  (up to a prop. constant, so stays in range of CPU)
normC = max( LL, [], 1);
LL = bsxfun( @minus, LL, normC );
likelihood = exp( LL );

if ~exist( 'TargetPsi', 'var') || isempty( TargetPsi )
    % ------------------------------------------- Actually sample z(t)
    [js, logQ] = SampleHMMStateSeqWithQsC( pi_z, likelihood, pi_init, -1, randi([1 100000]) );
    z = availFeatIDs( js );
else % ----------------------------- Calc prob of moving to Target's z
    jseq = zeros(1,T);
    for jj = 1:length( availFeatIDs )
        jseq( TargetPsi.stateSeq(ii).z == availFeatIDs(jj) ) = jj;
    end
    [js,logQ] = SampleHMMStateSeqWithQsC( pi_z, likelihood, pi_init, jseq, randi([1 100000]) );
    z = availFeatIDs( js );
end

if nargout > 1
    if exist( 'logQ', 'var' )
        logPrQ_z = logPrQ_z + logQ;
    else
        logPrQ_z = logPrQ_z + sum( log( qs ) );
    end
end

end % main function

%
% % Initialize state and sub-state sequences:
% qs = zeros(1,T);

author: michaelchughes

repo: NPBayesHMM


text:
function dataOut=prtUtilSpectralOutOfSampleExtension(dataOld,dataNew,eigVectors,eigValues,sigma)

%   dataOut=prtUtilSpectralOutOfSampleExtension(dataOld,dataNew,eigVectors,eigValues,sigma)
%   Performs an out of sample spectral dimensionality reduction.
%
%   Parameters:
%
%   dataOld - Features of in-sample data
%
%   dataNew - Features of out of sample data
%
%   eigVectors - Spectral features of in-sample data
%
%   eigValues - eigen values corresponding to eigen vectors of in-sample data
%
%   sigma - sigma used for radial basis function (RBF)
%
%
%       Example usage:
%       ds=rt(prtPreProcZmuv,prtDataGenMoon);
%       nEigs=2;
%       sigma=.2;
%
%       imagesc(ds)                          %Plot Moon data in feature space
%       title('prtDataGenMoon in Feature Space')
%
%       nPicked=30;
%       picked=randi(size(ds.X,1),[nPicked,1]);
%
%       dsNew=ds.retainObservations(picked);
%       [eigValues, eigVectors] = prtUtilSpectralDimensionalityReduction(ds.X, nEigs,'sigma',sigma);      %transform Data to Spectral Space
%       dsNew.X=prtUtilSpectralOutOfSampleExtension(ds.X,dsNew.X,eigVectors,eigValues,sigma);             %transform out of sample data to spectral space
%
%
%       figure;
%       imagesc(dsNew)                      %Plot Moon data in spectral space
%       title('prtDataGenMoon out of sample in Spectral Space')

%Reference: Bengio, Y., Paiement, J. F., Vincent, P., Delalleau, O., Le Roux, N., & Ouimet, M. (2004).
%           Out-of-sample extensions for lle, isomap, mds, eigenmaps, and spectral clustering.
%           Advances in neural information processing systems, 16, 177-184.

            n=size(dataOld,1);
            k=prtUtilRbfDist(dataNew,dataOld,'sigma',sigma);
            kNorm1=sqrt(mean(prtUtilRbfDist(dataNew,dataOld,'sigma',sigma),2));
            kNorm2=sqrt(mean(prtUtilRbfDist(dataOld,dataOld,'sigma',sigma),2));

            normKern = k./(bsxfun(@times,kNorm1,kNorm2'));

            normKernTransformed=repmat(normKern,[1,1,length(eigValues)]);
            eigTransformed=permute(repmat(eigVectors,[1,1,size(dataNew,1)]),[3,1,2]);

            sumOut=squeeze(sum(normKernTransformed.*eigTransformed,2));

            dataOut=bsxfun(@rdivide,sumOut,(n*eigValues)');

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function organizeRFexperiments_OldMitaineWithAerts(pathWORK,pathExperimentTextures,pathRF)

% IMPORTANT! CETTE FONCTION DOIT ABSOLUMENT ETRE AUTOMATISEE!!!

startpath = pwd;

cd(pathWORK)
load('clinical')
load('subTypes')

% LOADING VARIABLES
cd(pathExperimentTextures)
load('permTrain'), load('permTest')
trainText = load('training'); trainText = struct2cell(trainText); trainText = trainText{1};
testText = load('testing'); testText = struct2cell(testText); testText = testText{1};
PT_loco = load('testResultsModels_PET_Locoregional'); PT_loco = struct2cell(PT_loco); PT_loco = PT_loco{1};
CT_loco = load('testResultsModels_CT_Locoregional'); CT_loco = struct2cell(CT_loco); CT_loco = CT_loco{1};
PT_distant = load('testResultsModels_PET_Distant'); PT_distant = struct2cell(PT_distant); PT_distant = PT_distant{1};
CT_distant = load('testResultsModels_CT_Distant'); CT_distant = struct2cell(CT_distant); CT_distant = CT_distant{1};
PT_death = load('testResultsModels_PET_Death'); PT_death = struct2cell(PT_death); PT_death = PT_death{1};
CT_death = load('testResultsModels_CT_Death'); CT_death = struct2cell(CT_death); CT_death = CT_death{1};

% ORGANIZING OUTCOMES
training.outcomes.Locoregional = trainText.Locoregional.outcome;
training.outcomes.Distant = trainText.Distant.outcome;
training.outcomes.Death = trainText.Death.outcome;
testing.outcomes.Locoregional = testText.Locoregional.outcome;
testing.outcomes.Distant = testText.Distant.outcome;
testing.outcomes.Death = testText.Death.outcome;
training.timeToEvents.Locoregional = trainText.Locoregional.timeToEvent;
training.timeToEvents.Distant = trainText.Distant.timeToEvent;
training.timeToEvents.Death = trainText.Death.timeToEvent;
testing.timeToEvents.Locoregional = testText.Locoregional.timeToEvent;
testing.timeToEvents.Distant = testText.Distant.timeToEvent;
testing.timeToEvents.Death = testText.Death.timeToEvent;

% ORGANIZING CLINICAL DATA
age_train = [clinical.HGJ.Age;clinical.CHUS.Age]; age_test = [clinical.HMR.Age;clinical.CHUM.Age];
subType_train = [subTypes.HGJ;subTypes.CHUS]; subType_test = [subTypes.HMR;subTypes.CHUM];
tStage_train = [clinical.HGJ.T_stage;clinical.CHUS.T_stage]; tStage_test = [clinical.HMR.T_stage;clinical.CHUM.T_stage];
nStage_train = [clinical.HGJ.N_stage;clinical.CHUS.N_stage]; nStage_test = [clinical.HMR.N_stage;clinical.CHUM.N_stage];
tnmStage_train = [clinical.HGJ.TNM_stage;clinical.CHUS.TNM_stage]; tnmStage_test = [clinical.HMR.TNM_stage;clinical.CHUM.TNM_stage];
age_train = age_train(permTrain); age_test = age_test(permTest);
subType_train = subType_train(permTrain); subType_test = subType_test(permTest);
tStage_train = tStage_train(permTrain); tStage_test = tStage_test(permTest);
nStage_train = nStage_train(permTrain); nStage_test = nStage_test(permTest);
tnmStage_train = tnmStage_train(permTrain); tnmStage_test = tnmStage_test(permTest);

Age = age_train; SubType = subType_train; T_Stage = tStage_train; N_Stage = nStage_train; TNM_Stage = tnmStage_train;
training.clinical.table = table(Age,SubType,T_Stage,N_Stage,TNM_Stage);
training.clinical.categories = [0,1,1,1,1];
Age = age_test; SubType = subType_test; T_Stage = tStage_test; N_Stage = nStage_test; TNM_Stage = tnmStage_test;
testing.clinical.table = table(Age,SubType,T_Stage,N_Stage,TNM_Stage);
testing.clinical.categories = [0,1,1,1,1];

% ORGANIZING TEXTURES (change names if necessary)

% For Locoregional
PET_GLV = PT_loco.trainData.data(:,1); PET_GLN = PT_loco.trainData.data(:,2); PET_LGRE = PT_loco.trainData.data(:,3);
training.textures.Locoregional.PET = table(PET_GLV,PET_GLN,PET_LGRE);
PET_GLV = PT_loco.testData.data(:,1); PET_GLN = PT_loco.testData.data(:,2); PET_LGRE = PT_loco.testData.data(:,3);
testing.textures.Locoregional.PET = table(PET_GLV,PET_GLN,PET_LGRE);

CT_LGZE = CT_loco.trainData.data(:,1); CT_Correlation = CT_loco.trainData.data(:,2); CT_Busyness = CT_loco.trainData.data(:,3);
training.textures.Locoregional.CT = table(CT_LGZE,CT_Correlation,CT_Busyness);
CT_LGZE = CT_loco.testData.data(:,1); CT_Correlation= CT_loco.testData.data(:,2); CT_Busyness = CT_loco.testData.data(:,3);
testing.textures.Locoregional.CT = table(CT_LGZE,CT_Correlation,CT_Busyness);

PET_GLV = PT_loco.trainData.data(:,1); PET_GLN = PT_loco.trainData.data(:,2); PET_LGRE = PT_loco.trainData.data(:,3);
CT_LGZE = CT_loco.trainData.data(:,1); CT_Correlation = CT_loco.trainData.data(:,2); CT_Busyness = CT_loco.trainData.data(:,3);
training.textures.Locoregional.PETCT = table(PET_GLV,PET_GLN,PET_LGRE,CT_LGZE,CT_Correlation,CT_Busyness);
PET_GLV = PT_loco.testData.data(:,1); PET_GLN = PT_loco.testData.data(:,2); PET_LGRE = PT_loco.testData.data(:,3);
CT_LGZE = CT_loco.testData.data(:,1); CT_Correlation= CT_loco.testData.data(:,2); CT_Busyness = CT_loco.testData.data(:,3);
testing.textures.Locoregional.PETCT = table(PET_GLV,PET_GLN,PET_LGRE,CT_LGZE,CT_Correlation,CT_Busyness);

% For Distant
PT_SRHGE = PT_distant.trainData.data(:,1); PT_SUVmax = PT_distant.trainData.data(:,2); PT_ZSV = PT_distant.trainData.data(:,3);
training.textures.Distant.PET = table(PT_SRHGE,PT_SUVmax,PT_ZSV);
PT_SRHGE = PT_distant.testData.data(:,1); PT_SUVmax = PT_distant.testData.data(:,2); PT_ZSV = PT_distant.testData.data(:,3);
testing.textures.Distant.PET = table(PT_SRHGE,PT_SUVmax,PT_ZSV);

CT_LRHGE = CT_distant.trainData.data(:,1); CT_ZSV = CT_distant.trainData.data(:,2); CT_ZSN = CT_distant.trainData.data(:,3);
training.textures.Distant.CT = table(CT_LRHGE,CT_ZSV,CT_ZSN);
CT_LRHGE = CT_distant.testData.data(:,1); CT_ZSV = CT_distant.testData.data(:,2); CT_ZSN = CT_distant.testData.data(:,3);
testing.textures.Distant.CT = table(CT_LRHGE,CT_ZSV,CT_ZSN);

PT_SRHGE = PT_distant.trainData.data(:,1); PT_SUVmax = PT_distant.trainData.data(:,2); PT_ZSV = PT_distant.trainData.data(:,3);
CT_LRHGE = CT_distant.trainData.data(:,1); CT_ZSV = CT_distant.trainData.data(:,2); CT_ZSN = CT_distant.trainData.data(:,3);
training.textures.Distant.PETCT = table(PT_SRHGE,PT_SUVmax,PT_ZSV,CT_LRHGE,CT_ZSV,CT_ZSN);
PT_SRHGE = PT_distant.testData.data(:,1); PT_SUVmax = PT_distant.testData.data(:,2); PT_ZSV = PT_distant.testData.data(:,3);
CT_LRHGE = CT_distant.testData.data(:,1); CT_ZSV = CT_distant.testData.data(:,2); CT_ZSN = CT_distant.testData.data(:,3);
testing.textures.Distant.PETCT = table(PT_SRHGE,PT_SUVmax,PT_ZSV,CT_LRHGE,CT_ZSV,CT_ZSN);

% For Death
PT_SZLGE = PT_death.trainData.data(:,1); PT_LGZE = PT_death.trainData.data(:,2); PT_Contrast = PT_death.trainData.data(:,3);
training.textures.Death.PET = table(PT_SZLGE,PT_LGZE,PT_Contrast);
PT_SZLGE = PT_death.testData.data(:,1); PT_LGZE = PT_death.testData.data(:,2); PT_Contrast = PT_death.testData.data(:,3);
testing.textures.Death.PET = table(PT_SZLGE,PT_LGZE,PT_Contrast);

CT_GLN = CT_death.trainData.data(:,1); CT_SZE = CT_death.trainData.data(:,2); CT_Energy = CT_death.trainData.data(:,3);
training.textures.Death.CT = table(CT_GLN,CT_SZE,CT_Energy);
CT_GLN = CT_death.testData.data(:,1); CT_SZE = CT_death.testData.data(:,2); CT_Energy = CT_death.testData.data(:,3);
testing.textures.Death.CT = table(CT_GLN,CT_SZE,CT_Energy);

PT_SZLGE = PT_death.trainData.data(:,1); PT_LGZE = PT_death.trainData.data(:,2); PT_Contrast = PT_death.trainData.data(:,3);
CT_GLN = CT_death.trainData.data(:,1); CT_SZE = CT_death.trainData.data(:,2); CT_Energy = CT_death.trainData.data(:,3);
training.textures.Death.PETCT = table(PT_SZLGE,PT_LGZE,PT_Contrast,CT_GLN,CT_SZE,CT_Energy);
PT_SZLGE = PT_death.testData.data(:,1); PT_LGZE = PT_death.testData.data(:,2); PT_Contrast = PT_death.testData.data(:,3);
CT_GLN = CT_death.testData.data(:,1); CT_SZE = CT_death.testData.data(:,2); CT_Energy = CT_death.testData.data(:,3);
testing.textures.Death.PETCT = table(PT_SZLGE,PT_LGZE,PT_Contrast,CT_GLN,CT_SZE,CT_Energy);

% For Death (Aerts signature)
PT_death = load('testResultsAerts_PET_Death'); PT_death = struct2cell(PT_death); PT_death = PT_death{1};
CT_death = load('testResultsAerts_CT_Death'); CT_death = struct2cell(CT_death); CT_death = CT_death{1};

PT_Energy = PT_death.trainData.data(:,1); PT_Compactness = PT_death.trainData.data(:,2); PT_GLN = PT_death.trainData.data(:,3); PT_GLN_HLH = PT_death.trainData.data(:,4);
training.textures.DeathSign.PET = table(PT_Energy,PT_Compactness,PT_GLN,PT_GLN_HLH);
PT_Energy = PT_death.testData.data(:,1); PT_Compactness = PT_death.testData.data(:,2); PT_GLN = PT_death.testData.data(:,3); PT_GLN_HLH = PT_death.testData.data(:,4);
testing.textures.DeathSign.PET = table(PT_Energy,PT_Compactness,PT_GLN,PT_GLN_HLH);

CT_Energy = CT_death.trainData.data(:,1); CT_Compactness = CT_death.trainData.data(:,2); CT_GLN = CT_death.trainData.data(:,3); CT_GLN_HLH = CT_death.trainData.data(:,4);
training.textures.DeathSign.CT = table(CT_Energy,CT_Compactness,CT_GLN,CT_GLN_HLH);
CT_Energy = CT_death.testData.data(:,1); CT_Compactness = CT_death.testData.data(:,2); CT_GLN = CT_death.testData.data(:,3); CT_GLN_HLH = CT_death.testData.data(:,4);
testing.textures.DeathSign.CT = table(CT_Energy,CT_Compactness,CT_GLN,CT_GLN_HLH);

PT_Energy = PT_death.trainData.data(:,1); PT_Compactness = PT_death.trainData.data(:,2); PT_GLN = PT_death.trainData.data(:,3); PT_GLN_HLH = PT_death.trainData.data(:,4);
CT_Energy = CT_death.trainData.data(:,1); CT_Compactness = CT_death.trainData.data(:,2); CT_GLN = CT_death.trainData.data(:,3); CT_GLN_HLH = CT_death.trainData.data(:,4);
training.textures.DeathSign.PETCT = table(PT_Energy,PT_Compactness,PT_GLN,PT_GLN_HLH,CT_Energy,CT_Compactness,CT_GLN,CT_GLN_HLH);
PT_Energy = PT_death.testData.data(:,1); PT_Compactness = PT_death.testData.data(:,2); PT_GLN = PT_death.testData.data(:,3); PT_GLN_HLH = PT_death.testData.data(:,4);
CT_Energy = CT_death.testData.data(:,1); CT_Compactness = CT_death.testData.data(:,2); CT_GLN = CT_death.testData.data(:,3); CT_GLN_HLH = CT_death.testData.data(:,4);
testing.textures.DeathSign.PETCT = table(PT_Energy,PT_Compactness,PT_GLN,PT_GLN_HLH,CT_Energy,CT_Compactness,CT_GLN,CT_GLN_HLH);

% SAVING ORGANIZED DATA
cd(pathRF)
save('training','training'), save('testing','testing')
save('permTrain','permTrain'), save('permTest','permTest')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [a, b] = VBA_guessHyperpriors (y, rangeEV)
% // VBA toolbox //////////////////////////////////////////////////////////
%
% [a, b] = VBA_guessHyperpriors (y, rangeEV)
% propose shape and scale parameters to use as hyperpriors for the
% measurement noise precision, based on the part of the variance one can
% expect the model to explain.
%
% IN:
%   - y: data to be fitted
%   - rangeEV: range of the fraction of variance to be explained.
%       default = [0.1, 0.9], .ie between 10% and 90% of the total variance
% OUT:
%   - a, b: shape and scale parameters of the Gamma hyperprior
%
% /////////////////////////////////////////////////////////////////////////

%% check parameters
% =========================================================================
% fill in defaults
if nargin == 1
    rangeEV = [0.1, 0.9];
end

% check inputs
assert (numel (y) > 3, '*** VBA_guessHyperpriors: no enough datapoints!');
assert (numel (rangeEV) == 2, '*** VBA_guessHyperpriors: rangeEV must have 2 values.');
assert (VBA_isInRange (rangeEV, [1e-8, 1 - 1e-8]), '*** VBA_guessHyperpriors: rangeEV must be between 0 and 1.');

% cleanup
y(isnan (y)) = [];

%% Compute expected precision
% =========================================================================
% reformulate statistics of interest
varianceTotal = var (VBA_vec (y));
varianceUnexplained = varianceTotal * (1 - rangeEV);
precisionExpected = 1 ./ varianceUnexplained;

% shortcuts
sumPrecision = sum (precisionExpected);
diffPrecision = max (precisionExpected) - min (precisionExpected);

% precision range as 98% confidence interval
a = 6 * (sumPrecision / diffPrecision) ^ 2;
b = 12 * sumPrecision / (diffPrecision ^ 2);

author: MBB-team

repo: VBA-toolbox


text:
function [P3, S_bar, V, RO, Tr, Z, sigma_sq, phi, Q, mu0, sigma0, c] =...
    em_sfm(P, MD, K, mask, imsize, bbox, tol, max_em_iter, rhandedCoords, rotP3d)

%  Non-Rigid Structure From Motion with Gaussian/LDS Deformation Model
%  Copyright (c) by Lorenzo Torresani, Stanford University
%
%  Based on the following paper:
%
%  Lorenzo Torresani, Aaron Hertzmann and Christoph Bregler,
%     Learning Non-Rigid 3D Shape from 2D Motion, NIPS 16, 2003
%  http://cs.stanford.edu/~ltorresa/projects/learning-nr-shape/
%
%  Please refer to this publication if you use this program for
%  research or for technical applications.
%
%
%  INPUT:
%
%  P           - (2*T) x J tracking matrix:          P([t t+T],:) contains the 2D projections of the J points at time t
%  MD          - T x J missing data binary matrix:   MD(t, j)=1 if no valid data is available for point j at time t, 0 otherwise
%  K           - number of deformation basis
%  use_lds     - set to 1 to model deformations using a linear dynamical system; set to 0 otherwise
%  tol         - termination tolerance (proportional change in likelihood)
%  max_em_iter - maximum number of EM iterations
%  mask        - mask.poly_x (Tx1 cell array of x coordinates of mask) and
%                mask.poly_y (Tx1 cell array of y coordinates of mask)
%  imsize      - Tx2 array with height and width of images
%  bbox        - Bboxes of instances in image (Nx4)
%  params.normdim - Normalization dimension for object bounding boxes
%
%  OUTPUT:
%
%  P3          - (3*T) x J 3D-motion matrix:                    ( P3([t t+T t+2*T],:) contains the 3D coordinates of the J points at time t )
%  S_bar       - shape average:            3 x J matrix
%  V           - deformation shapes:       (3*K) x J matrix     ( V((n-1)*3+[1:3],:) contains the n-th deformation basis )
%  RO          - rotation:                 cell array           ( RO{t} gives the rotation matrix at time t )
%  Tr          - translation:              T x 2 matrix
%  Z           - deformation weights:      T x K matrix
%  sigma_sq    - variance of the noise in feature position
%  phi         - LDS transition matrix
%  Q           - LDS state noise matrix
%  mu0         - initial state mean
%  sigma0      - initial state variance
%  c - scale coefficients                   T X 1 matrix

%% Parsing passed parameters
params = get_params();

%% Checking Dimensions
if mod(size(P,1), 1) ~= 0,
   fprintf('Error: size(P) must be (2*T)xJ\n');
   return;
end

if (size(P,1)/2 ~= size(MD,1)) || (size(P,2) ~= size(MD,2))
   fprintf('Error: Size incompatibility between P and MD\n');
   return;
end

if mod(K, 1) ~= 0,
   fprintf('Error: K must be an integer value\n');
   return;
end
%% Initialization for Tr, S_bar, R, c

[T, J] = size(MD);

P_hat = P; % if any of the points are missing, P_hat will be updated during the M-step

% uses rank 3 factorization to get a first initialization for rotation and S_bar
[R_init, Trvect, S_bar] = rigidfac(P_hat, MD); % P_hat = R_init*S_bar + Trvect

% initialize c also here !!
c = ones(T,1);
%if scaling
    for t = 1:T
       c(t) = mean([norm(R_init(t,:)),norm(R_init(t+T,:))]); % just an estimate
    end
%end

Tr(:,1) = Trvect(1:T);
Tr(:,2) = Trvect(T+1:2*T);

R = zeros(2*T, 3);

%% enforces rotation constraints for initilizations
for t = 1:T,
   Ru = R_init(t,:);
   Rv = R_init(T+t,:);
   Rz = cross(Ru,Rv); if det([Ru;Rv;Rz])<0, Rz = -Rz; end;
   RO_approx = apprRot([Ru;Rv;Rz]);
   RO{t} = RO_approx;
   R(t,:) = RO_approx(1,:);
   R(t+T,:) = RO_approx(2,:);
end

%% Initializations for deformation shapes and weights
% given the initial estimates of rotation, translation and shape average, it initializes
% deformation shapes and weights through LSQ minimization of the reprojection error

[V, Z] = init_SB(P_hat, Tr, R, c , S_bar, K); %handled c inside this function

% initializes sigma_sq
E_zz_init = cov(Z);
E_zz_init = repmat(E_zz_init, T, 1);
sigma_sq = mstep_update_noisevar(P_hat, S_bar, V, Z', E_zz_init, RO, Tr, c); %handled c inside this function

phi = [];
mu0 = [];
sigma0 = [];
Q = [];

%% EM iterations

loglik = 0;
annealing_const = 60;
max_anneal_iter = round(max_em_iter/4);

pBar =  TimedProgressBar( max_em_iter, 30, ...
    'Training NRSfM ', ', completed ', 'Training NRSfM concluded in ' );
for em_iter=1:max_em_iter,
    %mean(c)
    %std(c)
    %mean(mean(abs(S_bar)))
    %% compute the hidden variables distributions
    [E_z, E_zz] = estep_compute_Z_distr(P_hat, S_bar, V, R, c, Tr, sigma_sq);     % (Eq 17-18) %handled c inside this function
    Z = E_z';

    %% updates shape basis
    [S_bar, V] = mstep_update_shapebasis(P_hat, E_z, E_zz, R, c, Tr);   % (Eq 21)

    %% fill in missing points
    if sum(MD(:))>0,
        P_hat = mstep_update_missingdata(P_hat, MD, S_bar, V, E_z, RO, c, Tr, mask, imsize, bbox, params.nrsfm.norm_dim,(em_iter>max_em_iter/2));    % Occlusion reasoning after half iterations over
    end

    %% updates rotation
    [RO, R] = mstep_update_rotation(P_hat, S_bar, V, E_z, E_zz, RO, c, Tr);       % (Eq 24)

    %RO = rotP3d;
    %for t = 1:T
    %    R(t,:) = RO{t}(1,:);
    %    R(t+T,:) = RO{t}(2,:);
    %end

    %% update C here
    %if scaling
        c = mstep_update_c(P_hat, S_bar, V, E_z, E_zz, RO, Tr); % Eq 48 (PAMI paper)
        frac = mean(c);
        c = c/frac;
        S_bar = frac*S_bar;
        V = frac*V;
    %end

    %% updates translation
    Tr = mstep_update_transl(P_hat, S_bar, V, E_z, RO, c);                        % (Eq 23)

    %% updates noise variance
    sigma_sq = mstep_update_noisevar(P_hat, S_bar, V, E_z, E_zz, RO, Tr, c);      % (Eq 22)
    if em_iter < max_anneal_iter,
        sigma_sq = sigma_sq * (1 + annealing_const*(1 - em_iter/max_anneal_iter));
    end

    %visualize_3D_point(S_bar,params.part_names,wireframe_car(),h);figure(h);clf;
    %% computes log likelihood
    oldloglik = loglik;
    loglik = compute_log_lik(P_hat, S_bar, V, E_z, E_zz, RO, c, Tr, sigma_sq);

%     if(mod(em_iter,40)==0)
%         fprintf('Iteration %d/%d: Error %f\n', em_iter, max_em_iter, -loglik/(size(P_hat,1)/2));
%     end

    if (em_iter <= 2),
     loglikbase = loglik;
    elseif (loglik < oldloglik)
     %fprintf('Violation');
     %keyboard;
    elseif 0 && ((loglik-loglikbase)<(1 + tol)*(oldloglik-loglikbase)),
     fprintf('\n');
     break;
    end
    pBar.progress();
end
pBar.stop();
%% Check for right handed coordinate flips
vec1 = S_bar(:, rhandedCoords(1,2)) - S_bar(:, rhandedCoords(1,1));
vec2 = S_bar(:, rhandedCoords(2,2)) - S_bar(:, rhandedCoords(2,1));
vec3 = cross(vec1,vec2);
v = S_bar(:, rhandedCoords(3,2)) - S_bar(:, rhandedCoords(3,1));

badCt = 0;
for t = 1:T
   z_t = Z(t,:);
   S = S_bar;
   for kk = 1:K,
      S = S+z_t(kk)*V((kk-1)*3+[1:3],:);
   end
   S = c(t,1)*RO{t}*S;
   Zs = S(3,:);
   badZ = mean(Zs(MD(t,:)));
   goodZ = mean(Zs(~MD(t,:)));
   if(goodZ > badZ)
       badCt = badCt+1;
   end
end
thresh = badCt/T;

isFlipped = 0;
if(thresh > 0.5 || (thresh >0.2 && dot(vec3,v)<0))
    refZ = diag([1 1 -1]);
    rotZ = diag([-1 -1 1]);
    isFlipped = 1;
    S_bar = -S_bar;
    V = -V;
    for i=1:length(RO)
        RO{i} = rotZ*RO{i};
    end
end

%% Changing to PASCAL 3D
fprintf('Converting to PASCAL 3D co-ordinate frame\n');
if(~isempty(rotP3d))
    tformRot = averageRotationTransform(RO,rotP3d);
else
    fprintf('PASCAL 3D rotations not found. Might be a problem with evaluation\n');
    tformRot = averageRotationTransform(RO,RO);
end

%tformRot = eye(3);
S_bar = tformRot'*S_bar;
for t = 1:T
    RO{t}  = RO{t}*tformRot;
end
for kk = 1:K
  V((kk-1)*3+[1:3],:) = tformRot'*V((kk-1)*3+[1:3],:);
end

%% compute the final answer
P3 = zeros(3*T, J);
for t = 1:T,
   z_t = Z(t,:);
   S = S_bar;
   for kk = 1:K,
      S = S+z_t(kk)*V((kk-1)*3+[1:3],:);
   end;
   S = c(t,1)*RO{t}*S; %% changed code here !

   P3([t t+T t+2*T], :) = S + [Tr(t, [1 2]) -mean(S(3,:))]'*ones(1,J);
end

repo: CategoryShapes


text:
function run_experiments_seg_summary()
rcnn.name = 'rcnn' ;
rcnn.opts = {...
  'type', 'rcnn', ...
  'model', 'data/models/imagenet-vgg-m.mat', ...
  'layer', 19} ;

rcnnvd.name = 'rcnnvd' ;
rcnnvd.opts = {...
  'type', 'rcnn', ...
  'model', 'data/models/imagenet-vgg-verydeep-19.mat', ...
  'layer', 41} ;

dcnn.name = 'dcnn' ;
dcnn.opts = {...
  'type', 'dcnn', ...
  'model', 'data/models/imagenet-vgg-m.mat', ...
  'layer', 13, ...
  'numWords', 64} ;

dcnnvd.name = 'dcnnvd' ;
dcnnvd.opts = {...
  'type', 'dcnn', ...
  'model', 'data/models/imagenet-vgg-verydeep-19.mat', ...
  'layer', 35, ...
  'numWords', 64} ;

dsift.name = 'dsift' ;
dsift.opts = {...
  'type', 'dsift', ...
  'numWords', 256, ...
  'numPcaDimensions', 80} ;

datasetList = { 'msrc'} ; % add 'os'
numSplits = [1, 1] ;

setupNameList = {'rcnn', 'dcnn', 'rdcnn'} ;
encoderList = {{rcnn}, {dcnn}, {rcnn dcnn}} ;
if 1
  setupNameList =  horzcat(setupNameList, {'rcnnvd', 'dcnnvd', 'rdcnnvd'}) ;
  encoderList = horzcat(encoderList, {{rcnnvd}, {dcnnvd}, {rcnnvd dcnnvd}}) ;
end

t = LatexTable() ;
t.begin().hline().vline().pf('dataset').pf('measure (\\%%)').vline() ;
for ee = 1: numel(setupNameList)
  t.pf('%s', upper(setupNameList{ee})) ;
end
t.vline().pf('SoA') ;
t.vline().endl().hline() ;

for ii = 1 : numel(datasetList)
  switch datasetList{ii}
    case 'msrc'
      dn = 'MSRC' ; m = 'msrcAcc' ;
      art = '86.5~\cite{ladicky10what}' ;
    case 'os'
      dn = 'OS' ; m = 'acc' ;
      art = '--' ;
    otherwise, assert(false) ;
  end
  t.pf(dn).pf(m) ;
  for ee = 1: numel(encoderList)
    score = [] ;
    for jj = 1 : numSplits(ii)
      opts =  os_setup(...
        'dataset', datasetList{ii}, 'seed', jj, ...
        'encoders', encoderList{ee}, ...
        'prefix', 'exp01', ...
        'suffix', setupNameList{ee}) ;
      fprintf('loading %s\n', opts.segResultPath) ;
      res = load(opts.segResultPath) ;
      switch m
        case 'map11', score(end+1) = res.test.map11 ;
        case 'msrcAcc', score(end+1) = res.test.msrcAcc ;
        case 'acc', score(end+1) = res.test.acc ;
        otherwise, assert(false) ;
      end
    end
    if numel(score) > 1
      t.pf('%5.1f \\tiny$\\pm%3.1f$', 100*mean(score), 100*std(score)) ;
    else
      t.pf('%5.1f', 100*mean(score)) ;
    end
    if jj == 1
      figure(1) ; clf ;
      switch m
        case 'acc'
          imagesc(res.test.confusion) ;
          title(sprintf('%s accuracy: %.1f %%', upper(opts.suffix), res.test.acc*100)) ;
          vl_xmkdir(opts.segPublishDir) ;
          vl_figaspect(1) ;
          vl_printsize(.7) ;
          print(fullfile(opts.segPublishDir, ['confusion.pdf']), '-dpdf') ;
      end
    end
  end
  t.pf('%s',art).endl() ;
end
t.hline() ;
str = t.end() ;
disp(str) ;

author: mcimpoi

repo: deep-fbanks

save_path: github-repos/MATLAB/mcimpoi-deep-fbanks


text:
clear all; close all;
m=256; n=256;
a=0.04;
k=-1/a;
I=k*log(1-rand(m, n));
figure;
subplot(121);  imshow(uint8(I));
subplot(122);  imhist(uint8(I));

author: HuangCongQing


text:
%achievability bound for AWGN channel under a long-term power constraint
%

function rate_ach_lt=awgn_ach_lt(snr_db, nn, error)

snr = 10^(snr_db);

rate_ach_lt = zeros(1,length(nn)); %achievable rate (kappa beta)

%Check whether a long-term power constraint will be beneficial
C=log(1+snr);
V = 1-1/(1+snr)^2;
n_min = ( (1+snr)/snr*sqrt(2*pi*snr)*(1-error)*exp(qfuncinv(error)^2/2) + qfuncinv(error)/(1+snr)^2/sqrt(V))^2;

n_index=find(nn <= n_min);

if length(n_index)>0
    rate_ach_lt(n_index) = kappabeta_ach(2*nn(n_index),error,snr,1)./nn(n_index); %use the kappa beta bound for short-term power constraint
end

index_large_n = find(nn>n_min);

number_taus=20; %increase this number of gain tightness (at the cost of complexity)
number_errors=20; %increase this number of gain tightness (at the cost of complexity)

for index_n= index_large_n
    n=nn(index_n);
    %tic
    disp(['ach_awgn_lt():  n=', num2str(n)]);

    temp_rate_arr=[];
    for error_n = (error/number_errors): (error /number_errors): (error * (number_errors-1)/number_errors)
        %ach SISO kappa beta bound
        P_n= snr*(1-error_n)./(1-error);
        tau_array = (error_n/number_taus) : (error_n/number_taus): (error_n*(number_taus-1)/number_taus);
        for tau = tau_array
            temp_rate = (log2(kappa_inf(tau, P_n)) - betaq_up_v2(1-error_n+tau, 2*n, P_n))/n;
            temp_rate_arr=[temp_rate_arr,temp_rate];
        end
        rate_ach_lt(index_n)=max(temp_rate_arr);
    end
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
classdef FilterOp < int32
    % enumeration used with the import routines
    % Each import routine takes a first argument that describes the behavior
    enumeration
        getDescription (0)
        getWebpage     (2)
        importCatalog  (1)
    end
end

author: CelsoReyes


text:
function conMaskLengthyStartVoxels(timeMaskFileOut, startVoxelMaskFile, pathsRootFile)
% Load path files and find average start time per start point defined in a
% start voxel mask.  We assume that the pathway file contains time required
% to create a path per path.
%
% conMaskLengthyStartVoxels(timeMaskFileOut, startVoxelMaskFile, pathsRootFile)
%
%
%
%
% HISTORY:
% 2007.12.11 AJS: wrote it.
%

% Get startpoint mask
niMask = niftiRead(startVoxelMaskFile);
if(isempty(niMask.data))
    error('Could not find the start voxel mask.');
end

% Create output structure
timeMask = zeros(size(niMask.data));
% Create num pathways structure
npData = zeros(size(niMask.data));
% Get pathway directory
[pathDir] = fileparts(pathsRootFile);

pathNames = dir(pathsRootFile);
disp(['Processing ' num2str(length(pathNames)) ' pathway files...']);
ticPercent = 10;
lastPercent = 0;
for nn=1:length(pathNames)
    curName = pathNames(nn).name;
    %disp(['For ' curName '...']);
    pathFile = fullfile(pathDir,curName);
    fid = fopen(pathFile,'rb','b');
    d = fread(fid,'float');fclose(fid);
    fp=1;
    %timeV = zeros(23000,1)-1;
    %nPathsV = zeros(23000,1)-1;
    pp=1;
    while(fp <= length(d))
        numstats    = d(fp);
        npoints     = d(fp+1);
        %timeV(pp)   = d(fp+2);
        %nPathsV(pp) = d(fp+3);
        [i,j,k] = xyz2ijk(d(fp+4),d(fp+5),d(fp+6),niMask);
        timeMask(i,j,k) = timeMask(i,j,k)+d(fp+2);
        npData(i,j,k) = npData(i,j,k)+1;
        fp = fp+2+numstats+3*npoints;
        pp=1+pp;
    end
    % Squeeze vectors
    %if(find(timeV==-1,1,'first')>1)
    %    timeV = timeV(1:find(timeV==-1,1,'first')-1);
    %    nPathsV = nPathsV(1:find(nPathsV==-1,1,'first')-1);
    %end
    % Replace zero times with minimum
    %timeV(timeV==0) = min(timeV(timeV>0));
    %totalTime = (sum(timeV))/3600;
    %disp(['Total time (per path and per cpu): ' num2str(totalTime) ' hours.']);
    %totalPaths = sum(nPathsV);
    %disp(['Num paths sampled to find terminating (per connecting path): ' num2str(totalPaths) ' paths.']);
    %maxPaths = max(nPathsV);
    %disp(['Max number of paths sampled to find terminating: ' num2str(maxPaths) ' paths.']);
    [lastPercent] = progress(nn,lastPercent,ticPercent,length(pathNames));
end
disp(['Total time required by pathway computation: ' num2str(sum( timeMask(:) )) ' seconds.']);
disp(['Total GM covered: ' num2str(sum( npData(:)>0 )) ' voxels.']);
avgTimePP = zeros(size(timeMask));
avgTimePP(npData>0) = timeMask(npData>0) ./ npData(npData>0);
disp(['Max. average time required per GM point ' num2str(max( avgTimePP(:) )) ' seconds']);

% Write out time mask file
timeMask(timeMask>0) = 1;
dtiWriteNiftiWrapper (timeMask, niMask.qto_xyz, timeMaskFileOut);

return;

function [i,j,k] = xyz2ijk(x,y,z,niData)
% Assume that the input paths are in mm space so we don't need to use ACPC
% offset
index = inv(niData.qto_xyz(1:3,1:3))*[x;y;z]+1;
i = round(index(1)); j=round(index(2)); k=round(index(3));
if(outBounds(i,j,k,niData.data))
    error(['Pathway position ' num2str(i) ', ' num2str(j) ', ' num2str(k) ' is out of bounds.']);
end
return;

function [bIsOut] = outBounds(i,j,k,data)
bIsOut = i<1 || i>size(data,1) || j<1 || j>size(data,2) || k<1 || k>size(data,3);
return;

function [lastPercent] = progress(curV,lastPercent,ticPercent,maxV)
curPercent = curV/maxV*100;
if(curPercent - lastPercent > ticPercent)
    disp([ num2str(curPercent) '% complete...']);
    lastPercent = curPercent;
end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [VMH2IDmappingAll,VMH2IDmappingPresent,VMH2IDmappingMissing]=getIDfromMetStructure(metabolite_structure,idName)

% INPUT
% metabolite_structure      Structure containg metabolite related informations
% and ID's
% idName                    Name of the ID as used in the metabolite structure
%                           to be retrieved (e.g., 'pubChemId')
%
% OUTPUT
% VMH2IDmappingAll          Mapping of all VMH metabolites present in the
%                           metabolite structure (including NaNs)
% VMH2IDmappingPresent      Mapping of all VMH metabolites present in the
%                           metabolite structure (excluding NaNs)
% VMH2IDmappingMissing      Abbreviations of metabolites that are NaN's in the
%                           metabolite structure
%
% IT, Aug 2020

Mets = fieldnames(metabolite_structure);
a = 1;
b = 1;
c = 1;
VMH2IDmappingAll = '';
VMH2IDmappingPresent = '';
VMH2IDmappingMissing = '';
for i = 1 : length(Mets)
    if length(find(isnan(metabolite_structure.(Mets{i}).(idName)))) == 0
        VMH2IDmappingPresent{b,1} =Mets{i};
        if isstr(metabolite_structure.(Mets{i}).(idName))
            id = (metabolite_structure.(Mets{i}).(idName));
        else
            id = num2str(metabolite_structure.(Mets{i}).(idName));
        end
        VMH2IDmappingPresent{b,2} = id;
        VMH2IDmappingPresent{b,3} = metabolite_structure.(Mets{i}).VMHId;
        if isfield(metabolite_structure.(Mets{i}),'metNames')
            VMH2IDmappingPresent{b,4} = metabolite_structure.(Mets{i}).metNames;
        else
            VMH2IDmappingPresent{b,4} = metabolite_structure.(Mets{i}).rxnNames;
        end
        b = b + 1;
    elseif isnan(metabolite_structure.(Mets{i}).(idName))
        VMH2IDmappingMissing{c,1} =Mets{i};
        c = c + 1;
    end

    VMH2IDmappingAll{a,1} =Mets{i};
    if isstr(metabolite_structure.(Mets{i}).(idName))
        id = (metabolite_structure.(Mets{i}).(idName));
    else
        id = num2str(metabolite_structure.(Mets{i}).(idName));
    end
    VMH2IDmappingAll{a,2} = id;
    a = a + 1;
end

author: opencobra

repo: cobratoolbox


text:
function areaBB = computeArea(bb)
%computes area of the bb=[xmin ymin xmax ymax]

if ((bb(1) > bb(3)) || (bb(2) > bb(4)))
    areaBB = 0;
else
    areaBB = (bb(3) - bb(1) + 1) * (bb(4) - bb(2) + 1);
end

author: Cloud-CV

repo: object-proposals


text:
%% Dense vs Sparse SNLE
clc
fun = @(x) [10*(x(2) - x(1)^2)
            sqrt(90)*(x(4)-x(3)^2)
            sqrt(10)*(x(2) + x(4) - 2)
            (1/sqrt(10))*(x(2) - x(4))];

x0 = [-30 -10 -30 -10]';

Opt = opti('fun',fun,'x0',x0,'options',optiset('solver','auto'))

[x,f,e,i] = solve(Opt)

%% Sparse above
fun = @(x) -1;
nleq = @(x) [10*(x(2) - x(1)^2)
            sqrt(90)*(x(4)-x(3)^2)
            sqrt(10)*(x(2) + x(4) - 2)
            (1/sqrt(10))*(x(2) - x(4))];
cl = zeros(4,1);
cu = zeros(4,1);

x0 = [-30 -10 -30 -10]';

Opt = opti('fun',fun,'nl',nleq,cl,cu,'x0',x0,'options',optiset('solver','auto'))

[x,f,e,i] = solve(Opt)

%% As above but new construct [no grad]
clc
nleq = @(x) [10*(x(2) - x(1)^2)
            sqrt(90)*(x(4)-x(3)^2)
            sqrt(10)*(x(2) + x(4) - 2)
            (1/sqrt(10))*(x(2) - x(4))];

x0 = [-30 -10 -30 -10]';

Opt = opti('nleq',nleq,'x0',x0,'options',optiset('solver','auto'))

[x,f,e,i] = solve(Opt)

%% As above but new construct [w grad DENSE]
clc
nleq = @(x) [10*(x(2) - x(1)^2)
            sqrt(90)*(x(4)-x(3)^2)
            sqrt(10)*(x(2) + x(4) - 2)
            (1/sqrt(10))*(x(2) - x(4))];
if (exist('syms.m','file'))
    [nljac,nljacstr] = symJac(nleq);

    x0 = [-30 -10 -30 -10]';

    Opt = opti('nleq',nleq,'nljac',nljac,'x0',x0,'options',optiset('solver','auto'))

    [x,f,e,i] = solve(Opt)
end

%% As above but new construct [w grad SPARSE]
clc
nleq = @(x) [10*(x(2) - x(1)^2)
            sqrt(90)*(x(4)-x(3)^2)
            sqrt(10)*(x(2) + x(4) - 2)
            (1/sqrt(10))*(x(2) - x(4))];
if (exist('syms.m','file'))
    [nljac,nljacstr] = symJac(nleq); nljac = @(x) sparse(nljac(x));

    x0 = [-30 -10 -30 -10]';

    Opt = opti('nleq',nleq,'nljac',nljac,'nljacstr',nljacstr,'x0',x0,'options',optiset('solver','auto'))

    [x,f,e,i] = solve(Opt)
end

%% As above but new construct [w grad SPARSE alt nl format]
clc
nleq = @(x) [10*(x(2) - x(1)^2)
            sqrt(90)*(x(4)-x(3)^2)
            sqrt(10)*(x(2) + x(4) - 2)
            (1/sqrt(10))*(x(2) - x(4))];
if (exist('syms.m','file'))
    [nljac,nljacstr] = symJac(nleq); nljac = @(x) sparse(nljac(x));

    x0 = [-30 -10 -30 -10]';

    Opt = opti('nleq',nleq,'nljac',nljac,'nljacstr',nljacstr,'x0',x0,'options',optiset('solver','auto'))

    [x,f,e,i] = solve(Opt)
end

%% As above but new construct [w grad SPARSE alt nl format w lin]
clc
nleq = @(x) [10*(x(2) - x(1)^2)
            sqrt(90)*(x(4)-x(3)^2)
            sqrt(10)*(x(2) + x(4) - 2)
            (1/sqrt(10))*(x(2) - x(4))];
if (exist('syms.m','file'))
    [nljac,nljacstr] = symJac(nleq); nljac = @(x) sparse(nljac(x));
    A = [0 1 0 1]; b = 2;

    x0 = [-30 -10 -30 -10]';

    Opt = opti('nleq',nleq,'ineq',A,b,'nljac',nljac,'nljacstr',nljacstr,'x0',x0,'options',optiset('solver','auto'))

    [x,f,e,i] = solve(Opt)
end

%% As above but new construct [w grad SPARSE alt nl format w lin + int]
clc
nleq = @(x) [10*(x(2) - x(1)^2)
            sqrt(90)*(x(4)-x(3)^2)
            sqrt(10)*(x(2) + x(4) - 2)
            (1/sqrt(10))*(x(2) - x(4))];
if (exist('syms.m','file'))
    [nljac,nljacstr] = symJac(nleq); nljac = @(x) sparse(nljac(x));
    A = [0 1 0 1]; b = 2;

    x0 = [-30 -10 -30 -10]';

    Opt = opti('nleq',nleq,'ivars',2,'lin',A,-Inf,b,'nljac',nljac,'nljacstr',nljacstr,'x0',x0,'options',optiset('solver','auto'))

    [x,f,e,i] = solve(Opt)
end

%% As above but new construct [w grad SPARSE alt nl format w ineq]
clc
nleq = @(x) [10*(x(2) - x(1)^2)
            sqrt(90)*(x(4)-x(3)^2)
            sqrt(10)*(x(2) + x(4) - 2)
            (1/sqrt(10))*(x(2) - x(4))
            x(3)];
if (exist('syms.m','file'))
    [nljac,nljacstr] = symJac(nleq); nljac = @(x) sparse(nljac(x));
    nlrhs = zeros(5,1);
    nle = [zeros(4,1);-1];
    A = [0 1 0 1]; b = 2;

    x0 = [-30 -10 -30 -10]';

    Opt = opti('nlmix',nleq,nlrhs,nle,'ineq',A,b,'nljac',nljac,'nljacstr',nljacstr,'x0',x0,'options',optiset('solver','auto'))

    [x,f,e,i] = solve(Opt)
end

%% Wiki Ex 1
clc
% System of Nonlinear Equations
 nleq = @(x) [ 2*x(1) - x(2) - exp(-x(1));
             -x(1) + 2*x(2) - exp(-x(2))];

% Starting Guess
 x0 = [-5;5];

% Create OPTI Object
 Opt = opti('nleq',nleq,'x0',x0)

% Solve the SNLE problem
[x,fval,exitflag,info] = solve(Opt)

%% Online Example
clc
% System of Nonlinear Equations
 nleq = @(x) [10*(x(2) - x(1)^2)
             sqrt(90)*(x(4) - x(3)^2)
             sqrt(10)*(x(2) + x(4) - 2)
             (1/sqrt(10))*(x(2) - x(4))];

% Nonlinear Equations Jacobian
 nlJac = @(x) sparse([-20*x(1),10,0,0
                      0,0,-6*10^(1/2)*x(3),3*10^(1/2)
                      0,10^(1/2),0,10^(1/2)
                      0,10^(1/2)/10,0,-10^(1/2)/10]);

% Jacobian Sparsity Pattern
 nlJacstr = @() sparse([1 1 0 0
                        0 0 1 1
                        0 1 0 1
                        0 1 0 1]);

% Starting Guess
 x0 = [-30;-10;-30;-10];

% Sparse SNLE OPTI Problem
 Opt = opti('nleq',nleq,'nlJac',nlJac,'nlJacstr',nlJacstr,'x0',x0)

  %Solve
[x,f,e,i] = solve(Opt)

%% SCNLE Example
clc

% Objective (Nonlinear Equations) Function
 fun = @(x) [ 2*x(1) - x(2) - exp(-x(1));
             -x(1) + 2*x(2) - exp(-x(2))];

% Bounds
 lb = [0.6;0];
 ub = [1;1];

% Starting Guess
 x0 = [-5;5];

% Create OPTI Object
 Opt = opti('nleq',fun,'bounds',lb,ub,'x0',x0,'options',optiset('display','iter'))

% Solve the SCNLE problem
[x,fval,exitflag,info] = solve(Opt)

%% Problem SCNLE
clc

A = randn(91);
B12 = randn(91,145);
B1 = randn(145,91); B2 = randn(145,91);
Constant = randn(91,1);
A_ieq = randn(12,91); b_ieq = randn(12,1);
lb = zeros(91,1); ub = 100*ones(91,1);
z0 = lb;

nleq = @(x) A*x + B12*((B1*x).*(B2*x)) + Constant;

Opt=opti('nleq',nleq,'ineq',A_ieq,b_ieq,'bounds',lb,ub,'x0',z0,'options',optiset('display','iter'))

solve(Opt)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [estimate] = estimate_completion_time(sequences, varargin)
% estimate_completion_time Estimates completion time for a tracker
%
% Estimates time for a pass over a given set of sequences given some basic performance characteristics.
%
% Input:
% - sequences: Cell array of sequence structures.
% - varargin[FPS]: Average speed of a tracker (frames per second).
% - varargin[Failures]: Average number of failures per sequence.
% - varargin[Repeats]: Required number of repetitions for each sequence.
%
% Output:
% - estimate: Completion time estimate given in number of seconds

failures = 8;
fps = 0.5;
repeats = 1;

args = varargin;
for j=1:2:length(args)
    switch varargin{j}
        case 'fps', fps = max(0, args{j+1});
        case 'failures', failures = max(0, args{j+1});
        case 'repeats', repeats = max(1, args{j+1});
        otherwise, error(['unrecognized argument ' args{j}]);
    end
end

total_frames = 0;

for i = 1:length(sequences)
    total_frames = total_frames + sequences{i}.length;
end;

average_frames = total_frames / length(sequences);

lens = linspace(0, 1, failures + 1);
reruns = sum(lens);

estimate = (average_frames * reruns * repeats) * length(sequences) / fps;

author: votchallenge

repo: toolkit-legacy


text:
% This subroutine assigns creates a 3D grid with
% spacing dx,dy, dz (in degreees). The size will
% be selected interactiVELY. The pvalue in each
% volume around a grid point containing ni earthquakes
% will be calculated as well as the magnitude
% of completness
%   Stefan Wiemer 1/98

report_this_filefun(mfilename('fullpath'));
global no1 bo1 inb1 inb2

if sel == 'in'
    % get the grid parameter
    % initial values
    %
    dx = 0.1;
    dy = 0.1 ;
    dz = 5.00 ;
    ni = 300;
    R = 10000;

    def = {'1982.6','3','200','0.1','0.1',num2str(dz),num2str(max(a.Depth)), num2str(min(a.Depth))};

    tit ='Three dimesional z-value analysis';
    prompt={'Time of analysis?',...
        'Window length in years ?',...
        'Sample size N?',...
        'Spacing in Longitude (dx in [deg])',...
        'Sapcing in Latitude  (dy in [deg])',...
        'Spacing in Depth    (dz in [km ])',...
        'Depth Range: deep limit [km] ',...
        'Depth Range: shallow limit',...
        };

    ni2 = inputdlg(prompt,tit,1,def);

    l = ni2{1}; ti  = str2double(l);
    l = ni2{2}; ZG.compare_window_yrs= str2double(l);
    l = ni2{3}; ni= str2double(l);
    l = ni2{4}; dx= str2double(l);
    l = ni2{5}; dy= str2double(l);
    l = ni2{6}; dz= str2double(l);
    l = ni2{7}; z1= str2double(l);
    l = ni2{8}; z2= str2double(l);

    sel = 'ca'; zgrid3db

end   % if sel == 'in'

% get the grid-size interactively and
% calculate the b-value in the grid by sorting
% thge seimicity and selectiong the ni neighbors
% to each grid point

if sel == 'ca'
    selgp3dB
    zvect=[z2:dz:z1];
    gz = zvect;
    itotal = length(gx)*length(gz)*length(gy);
    zmap_message_center.set_info(' ','Running... ');think
    %  make grid, calculate start- endtime etc.  ...
    %
    zvg = ones(length(gx),length(gy),length(gz))*nan;
    ra  = ones(length(gx),length(gy),length(gz));

    t0b = min(a.Date)  ;
    n = a.Count;
    teb = a(n,3) ;
    tdiff = round((teb - t0b)*365/par1);
    loc = zeros(3, length(gx)*length(gy));

    % loop over  all points
    %
    i2 = 0.;
    i1 = 0.;
    allcount = 0.;
    wai = waitbar(0,' Please Wait ...  ');
    set(wai,'NumberTitle','off','Name',' 3D gridding - percent done');;
    drawnow
    %
    %
    tl = teb - t0b;
    iwl = floor(ZG.compare_window_yrs*365/par1);
    t = floor((ti-t0b)*365/par1);

    z0 = 0; x0 = 0; y0 = 0; dt = 1;
    % loop over all points
    for x = min(gx):dx:max(gx)
        x0 = x0+1;
        for y = min(gy):dy:max(gy)
            y0 = y0+1;
            for z = min(gz):dz:max(gz)
                z0 = z0+1;
                allcount = allcount + 1.;
                i2 = i2+1;

                % calculate distance from center point and sort wrt distance
                l = sqrt(((a.Longitude-x)*cos(pi/180*y)*111).^2 + ((a.Latitude-y)*111).^2 + ((a.Depth - z)).^2 ) ;
                [s,is] = sort(l);
                b = a(is(:,1),:) ;       % re-orders matrix to agree row-wise

                % take first ni points
                b = b(1:ni,:);      % new data per grid point (b) is sorted in distance

                [bv magco stan av me mer me2,  pr] =  bvalca3(b,inb1,inb2);
                l2 = sort(l);
                zvg(x0,y0,z0) = bv;
                ra(x0,y0,z0) = l2(ni);

                zvg(x0,y0,z0) = (mean1 - mean2)/(sqrt(var1/(ncu-iwl)+var2/iwl));
                ra(x0,y0,z0) = l(ni);
                waitbar(allcount/itotal)
            end  % for z
            z0 = 0;
        end  % for y
        y0 = 0;
    end  % for x
    x0 = 0;

    % save data
    %
    catSave3 =...
        [ 'zmap_message_center.set_info(''Save Grid'',''  '');think;',...
        '[file1,path1] = uiputfile(fullfile(hodi, ''eq_data'', ''*.mat''), ''Grid Datafile Name?'') ;',...
        ' sapa2 = [''save '' path1 file1 '' zvg ra gx gy gz dx dy dz d par1 tdiff t0b teb a main faults mainfault coastline yvect xvect tmpgri ll''];',...
        ' if length(file1) > 1, eval(sapa2),end , done']; eval(catSave3)

    close(wai)
    watchoff

    gz = -gz;
    zv2 = zvg;
    ac2 = 'new'; myslicer;

end  % if cal


text:
% Chapter 17 - Neural Networks.
% Programs_17e - Bifurcation Diagram for a Simple Bistable Neuromodule.
% Copyright Birkhauser 2013. Stephen Lynch.

% Bifurcation diagram for a two-neuron module. (See Figure 17.15).
% Vary bias b1.
clear all
format long;
halfN=1000;N=2*halfN+1;N1=1+halfN;Max=10;a=1;alpha=0.3;
x=zeros(1,N);y=zeros(1,N);
b2=3;w11=7;w12=-4;w21=5;start=-5;
x(1)=-10;y(1)=-3;

% Ramp the power up
for n=1:halfN
    b1=(start+n*Max/halfN);
    x(n+1)=b1+w11*(exp(a*x(n))-exp(-a*x(n)))/(exp(a*x(n))+exp(-a*x(n)))+w12*(exp(alpha*y(n))-exp(-alpha*y(n)))/(exp(alpha*y(n))+exp(-alpha*y(n)));
    y(n+1)=b2+w21*(exp(a*x(n))-exp(-a*x(n)))/(exp(a*x(n))+exp(-a*x(n)));
end

% Ramp the power down
for n=N1:N
    b1=(start+2*Max-n*Max/halfN);
    x(n+1)=b1+w11*(exp(a*x(n))-exp(-a*x(n)))/(exp(a*x(n))+exp(-a*x(n)))+w12*(exp(alpha*y(n))-exp(-alpha*y(n)))/(exp(alpha*y(n))+exp(-alpha*y(n)));
    y(n+1)=b2+w21*(exp(a*x(n))-exp(-a*x(n)))/(exp(a*x(n))+exp(-a*x(n)));
end

% Plot the bifurcation diagrams
fsize=14;
subplot(2,1,1)
hold on
set(gca,'XTick',0:halfN/2:N,'FontSize',fsize);
set(gca,'YTick',-Max:5:Max,'FontSize',fsize);
plot(x(1:N),'-','MarkerSize',1,'color','k')
xlabel('Number of Iterations','FontSize',fsize);
ylabel('x_n','FontSize',fsize);
hold off
x1=zeros(1,N);w=zeros(1,N);

for n=1:halfN
    x1(n)=x(N+1-n);
    w(n)=start+n*Max/halfN;
end

subplot(2,1,2)
hold on
set(gca,'XTick',start:5:start+Max,'FontSize',fsize);
set(gca,'YTick',-Max:5:Max,'FontSize',fsize);
plot(w(1:halfN),x(1:halfN),'-','MarkerSize',1,'color','k');
plot(w(1:halfN),x1(1:halfN),'-','MarkerSize',1,'color','k');
xlabel('b_1','FontSize',fsize);
ylabel('x_n','FontSize',fsize);
hold off

% End of Programs_17e.

author: Sable

repo: mcbench-benchmarks


text:
function [candidates, exact] = project_candidates(candidates, H)
  n = size(candidates, 1);

  % project all 4 coordinates of the bbox
  %   1 ---- 2
  %   |      |
  %   4 ---- 3
  points = [cat(1, ...
    candidates(:,1:2), candidates(:,[3,2]), ...
    candidates(:,3:4), candidates(:,[1,4])), ...
    ones(4*n, 1)];
  projected_points = points * H';
  projected_points = projected_points ./ repmat(projected_points(:,3), [1 3]);
%   p1 = projected_points(1:n, 1:2);
%   p2 = projected_points((n+1):(2*n), 1:2);
%   p3 = projected_points((2*n+1):(3*n), 1:2);
%   p4 = projected_points((3*n+1):end, 1:2);
  xs = reshape(projected_points(:,1), [n 4]);
  ys = reshape(projected_points(:,2), [n 4]);
  exact = [xs(:,1), ys(:,1), xs(:,2), ys(:,2), xs(:,3), ys(:,3), xs(:,4), ys(:,4)];

  xs = sort(xs, 2);
  ys = sort(ys, 2);
  candidates(:,1) = mean(xs(:,1:2), 2);
  candidates(:,3) = mean(xs(:,3:4), 2);
  candidates(:,2) = mean(ys(:,1:2), 2);
  candidates(:,4) = mean(ys(:,3:4), 2);
end

author: hosang

repo: detection-proposals


text:
function dataout = ft_nirs_referencechannelsubtraction(cfg, datain)

% FT_NIRS_REFERENCECHANNELSUBTRACTION performs NIRS reference channel subtraction
%
% Use as
%   dataout = ft_nirs_referencechannelsubtraction(cfg, datain)
% where cfg is a configuration structure and indata is raw NIRS-data that is
% represented according to the output of FT_PREPROCESSING.
%
% The configuration should contain the following options
%   cfg.shortdistance = scalar, below which distance a channel is regarded
%                       as short in cm (default = 1.5)
%   cfg.closedistance = scalar, defines the maximal distance between a
%                       long and a short channel in cm (default = 15).
%                       NOT APPLIED CURRENTLY!
%   cfg.method        = string, 'regstat2', 'QR' or 'OLS' (default = 'QR')
%   cfg.verbose       = boolean, whether text output is desired (default = false)
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_NIRS_SCALPCOUPLINGINDEX, FT_NIRS_SIGNALQUALITYINDEX, FT_NIRS_TRANSFORM_ODS

% You are using the FieldTrip NIRS toolbox developed and maintained by
% Artinis Medical Systems (http://www.artinis.com). For more information
% on FieldTrip, see http://www.fieldtriptoolbox.org
%
% This work is licensed under a Creative Commons Attribution-ShareAlike 4.0
% International License. To view a copy of this license, visit
% http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to
% Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
%
% Creative Commons Attribution-ShareAlike 4.0 International License:
% -----------------------------------
% You are free to:
%
%     Share - copy and redistribute the material in any medium or format
%     Adapt - remix, transform, and build upon the material
%     for any purpose, even commercially.
%
%     The licensor cannot revoke these freedoms as long as you follow the
%     license terms.
%
% Under the following terms:
%
%     Attribution - You must give appropriate credit, provide a link to
%                    the license, and indicate if changes were made. You
%                    may do so in any reasonable manner, but not in any way
%                    that suggests the licensor endorses you or your use.
%
%     ShareAlike - If you remix, transform, or build upon the material,
%                   you must distribute your contributions under the same
%                   license as the original.
%
%     No additional restrictions - You may not apply legal terms or
%                                   technological measures that legally
%                                   restrict others from doing anything the
%                                   license permits.
%
% -----------------------------------
%
% This toolbox is not to be used for medical or clinical purposes.
%
% Copyright (c) 2016 by Artinis Medical Systems.
% Contact: askforinfo@artinis.com
%
% Main programmer:
% Marc van Wanrooij, DCN, http://www.neural-code.com
% Jörn M. Horschig, Artinis Medical Systems BV, http://www.artinis.com
% $Id$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% the initial part deals with parsing the input options and data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function

% the ft_preamble function works by calling a number of scripts from
% fieldtrip/utility/private that are able to modify the local workspace

ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar    datain
ft_preamble provenance datain

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  % do not continue function execution in case the outputfile is present and the user indicated to keep it
  return
end

% ensure that the input data is raw NIRS-data, this will also do
% backward-compatibility conversions of old data that for example was
% read from an old *.mat file
datain = ft_checkdata(datain, 'datatype', 'raw', 'senstype', 'nirs');

% get the options
cfg.shortdistance = ft_getopt(cfg, 'shortdistance', 1.5);
cfg.closedistance = ft_getopt(cfg, 'closedistance', 15);
cfg.method        = ft_getopt(cfg, 'method', 'QR');
cfg.verbose       = ft_getopt(cfg, 'verbose', false);

% check the options
cfg = ft_checkopt(cfg, 'method', 'char', {'regstat2', 'QR', 'OLS'});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% the actual computation is done in the middle part
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute distances
distances   = ft_nirs_optodedistance(datain);

% compute correlation
selshallow	= distances<cfg.shortdistance;
shallowIdx  = find(selshallow);
seldeep     = distances>=cfg.shortdistance;
deepIdx     = find(seldeep);

dataout			 = datain;
dataout.label = datain.label(seldeep);

for tr=1:numel(datain.trial)
  dat = datain.trial{tr};
  shallow		= dat(selshallow,:);
  shallow		= bsxfun(@minus,shallow,mean(shallow,2)); % mean detrend
  shallowlabel = datain.label(selshallow);

  deep		= dat(seldeep,:);
  deep		= bsxfun(@minus,deep,mean(deep,2)); % mean detrend
  deeplabel = datain.label(seldeep);

  time		= datain.time{tr};

  %% Reference channel subtraction
  ndeep		= size(deep,1);
  signal		= NaN(size(deep));
  x			= shallow';
  for dpIdx	= 1:ndeep
    y				= deep(dpIdx,:)';
    switch (cfg.method)
      case 'regstat2'
        b				= regstats2(y,x,'linear',{'beta','r'});
        beta    = b.beta;
        res			= b.r;

      case 'QR'
        [Q,R] = qr(x,0);
        beta = R\(Q'*y);
        yhat = x*beta;
        res = y - yhat;

      case 'OLS'
        x2 = [ones(size(x, 1), 1) x];
        beta = x2\y;
        yhat = x2*beta;
        res  = y - yhat;
        beta(1) = [];

      otherwise % it should never come here as we use ft_checkopt
        error('unrecognized method');
    end

    signal(dpIdx,:) = res';

    % sanity check of results
    if cfg.verbose
      fprintf('Found the following meaningful shallow channels for deep channel %s:"', deeplabel{dpIdx});
      shIdx = find(beta>0.5);
      deepPos = datain.opto.chanpos(deepIdx(dpIdx), :);
      shallowPos = datain.opto.chanpos(shallowIdx(shIdx), :);
      posDist = sqrt(sum((shallowPos - repmat(deepPos, size(shallowPos, 1), 1)).^2, 2));
      closeChan = posDist < cfg.closedistance;
      for s=1:numel(shIdx)
        if (closeChan(s))
          fprintf('\n\t%s (d=%.2fcm)', shallowlabel{shIdx(s)}, posDist(s))
        else
          fprintf('\n\t#%s (d=%.2fcm)', shallowlabel{shIdx(s)}, posDist(s))
        end
      end
      fprintf('\n\n')
    end
  end

  % overwrite
  dataout.time{tr}	 = time;
  dataout.trial{tr} = signal;
end

% Auto-correlation
% rshallow	= corrcoef(shallow'); % correlation matrix shallow channels
% rdeep     = corrcoef(deep');    % correlation matrix deep channels
% rsignal		= corrcoef(signal');  % correlation matrix signal channels

% this might involve more active checking of whether the input options
% are consistent with the data and with each other

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% deal with the output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% do the general cleanup and bookkeeping at the end of the function

% the ft_postamble function works by calling a number of scripts from
% fieldtrip/utility/private that are able to modify the local workspace

ft_postamble debug

ft_postamble previous   datain
ft_postamble provenance dataout
ft_postamble history    dataout
ft_postamble savevar    dataout

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function tree = prtUtilRecursiveCapTreeFast(Obj,tree,x,y,index)
%tree = recursiveCapTree(Obj,tree,x,y,index)

nFeatures = size(x,2);

uniqueY = find(sum(y,1)>0);

if index > tree.maxReservedLen
    tree.W = memorySaverAppendNulls(tree.W,Obj.Memory.nAppend,Obj.nFeatures);
    tree.threshold = memorySaverAppendNulls(tree.threshold,Obj.Memory.nAppend,1);
    tree.featureIndices = memorySaverAppendNulls(tree.featureIndices,Obj.Memory.nAppend,Obj.nFeatures);
    tree.treeIndices = memorySaverAppendNulls(tree.treeIndices,Obj.Memory.nAppend,1);
    tree.terminalVote = memorySaverAppendNulls(tree.terminalVote,Obj.Memory.nAppend,1);
    tree.maxReservedLen = tree.maxReservedLen + Obj.Memory.nAppend;
end

%Base cases; if there is only one class left, we must return

if length(uniqueY) == 1
    tree.W(:,index) = inf;  %place holder for completed processing
    tree.treeIndices(index) = tree.father;
    tree.terminalVote(index) = uniqueY;
    return;
end

%Choose random subspace projection:
if Obj.featureSelectWithReplacement  %allow redundant features
    tree.featureIndices(:,index) = ceil(rand([1,Obj.nFeatures])*nFeatures)';
    tree.treeIndices(index) = tree.father;
else
    locInd = randperm(nFeatures)'; %do NOT allow redundant features
    tree.featureIndices(:,index) = locInd(1:Obj.nFeatures);
    tree.treeIndices(index) = tree.father;
end

xTrain = x(:,tree.featureIndices(:,index));

% Generate a CAP classifier (internal code, don't use prtClassCap for
% speed issues)
[w,thresholdValue,yOut] = recursiveCapTreeGenerateCap(xTrain,y);

if any(~isfinite(w)) || ~isfinite(thresholdValue)
    exitNow = true;
else

    yOut = yOut >= thresholdValue;

    % True if only one class in the output
    exitNow = ~any(yOut) || all(yOut);
end

if exitNow
    % If we get here we have the same mean under both class labels
    % or have infs in the data
    % So we exit

    % Tree params on an exit
    tree.W(:,index) = inf;  %place holder for completed processing
    tree.threshold(:,index) = nan;
    tree.treeIndices(index) = tree.father;

    % Figure out the dominant class
    classCounts = sum(y,1);
    [maxClassCounts, maxClassInd] = max(classCounts);
    isGoodClass = classCounts == maxClassCounts;

    nGoodInds = sum(isGoodClass);

    % If nGoodInds > 1 We have ties and need to decide randomly
    if nGoodInds > 1
        goodInds = find(isGoodClass);
        maxClassInd = goodInds(ceil(rand*nGoodInds));
    end

    %tree.terminalVote(index) = uniqueY(maxClassInd);
    tree.terminalVote(index) = maxClassInd;
    return
end

% Continue on...
% Store node info into tree

tree.W(:,index) = w(:);
tree.threshold(:,index) = thresholdValue;
tree.father = index;

% Split the data

% Left
xLeft = x(~yOut,:);
yLeft = y(~yOut,:);
tree = prtUtilRecursiveCapTreeFast(Obj,tree,xLeft,yLeft,index + 1);
tree.father = index;

% Right
xRight = x(yOut,:);

yRight = y(yOut,:);
maxLen = length(find(~isnan(tree.W(1,:))));
tree = prtUtilRecursiveCapTreeFast(Obj,tree,xRight,yRight,maxLen + 1);

    function M = memorySaverAppendNulls(M,nAppend,nFeats)
        M = cat(2,M,nan(nFeats,nAppend));
    end
end

function [w,thresholdValue,yOut] = recursiveCapTreeGenerateCap(x,y)
%[w,thresholdValue,yOut] = recursiveCapGenerateCap(x,y,uniqueY,nRocEvals)
% Internal function to quicly generate a CAP classifier without prt
% overhead

if size(y,2) > 2
    classCounts = sum(y,1);
    [sortedClassCounts, sortedClassInds] = sort(classCounts,'descend'); %#ok<ASGLU>
else
    sortedClassInds = [1 2];
end
% Ties go to the lower classInds

mean0 = mean(x(y(:,sortedClassInds(1)),:),1);
mean1 = mean(x(y(:,sortedClassInds(2)),:),1);
w = mean1 - mean0;

w = w./norm(w);

if isnan(w)
    w = nan(size(w));
    thresholdValue = nan;
    yOut = nan;
    return
end

selectedClasses = y(:,sortedClassInds(1)) | y(:,sortedClassInds(2));

yOut = (w*x')';
% yOut = yOut(selectedClasses,:);
yTest = y(selectedClasses,sortedClassInds(2));
tempY = yOut(selectedClasses,:);
threshold = median(tempY);

pE = 1-sum(tempY >= threshold == yTest)/length(tempY);
if pE >= .5
    %     threshold = -threshold;
    %     w = -w;
    %     yOut = -yOut;
    %
    %     pE = 1-sum(tempY >= threshold == yTest)/length(tempY);
    %     %     if pE > .5
    %         %We tried the median, then -w, and the -median, and prob(error) was
    %         %still too high. That says that the median just won't work.  So use
    %         %the ROC curve:
    %
    %         w = -w; %un-negate (!)
    %         yOut = (w*x')';
    %         yTest = y(selectedClasses,sortedClassInds(2));
    %
    %         % figure out the threshold: (1)
    %         [pf,pd,thresh] = prtScoreRoc(yOut(selectedClasses,:),yTest,'uniqueLabels',[0 1]);
    %
    %         % figure out the threshold: (2)
    %         % uY = unique(yTest);
    %         % binLabels = yTest ~= uY(1);
    %         % [~,~,pf,pd,thresh] = prtUtilMinPeThreshold(yOut(selectedClasses,:),binLabels);
    %
    %         pE = prtUtilPfPd2Pe(pf,pd);
    %         [minPe,I] = min(pE);
    %         if numel(I) > 1
    %             I = unique(I);
    %         end
    %
    %         if I < length(thresh)
    %             %for categorical variables, the threshold should live between the
    %             %samples we've seen; this is true for continuous variables too.  Makes
    %             %a big difference when multiple values take the same value
    %             thresholdValue = mean([thresh(I),thresh(I+1)]);
    %         else
    %             thresholdValue = thresh(I);
    %         end
    %
    %         if minPe >= 0.5
    %             w = -w;
    %             yOut = -yOut;
    %             % % figure out the threshold:
    %             [pf,pd,thresh] = prtScoreRoc(yOut(selectedClasses,:),yTest,'uniqueLabels',[0 1]);
    %             pE = prtUtilPfPd2Pe(pf,pd);
    %             [minPe,I] = min(pE);
    %             if numel(I) > 1
    %                 I = unique(I);
    %             end
    %             if I < length(thresh)
    %                 %for categorical variables, the threshold should live between the
    %                 %samples we've seen; this is true for continuous variables too.  Makes
    %                 %a big difference when multiple values take the same value
    %                 thresholdValue = mean([thresh(I),thresh(I+1)]);
    %             else
    %                 thresholdValue = thresh(I);
    %             end
    %
    %         end
    %         threshold = thresholdValue;
    %         pE = minPe;
    %     end
end
if pE > .5
    threshold = -inf;
end
thresholdValue = threshold;
end


text:
function colors = colorSamples(lst)
%
% lst:  A vector of integers between 1 and 20 to draw some samples
%
% RF/BW

% Check parameters

% This is static.
% This is a good list of colors that Renzo made from D3.  We use these for
% various tasks like selecting the fiber colors
colorList = [ ...
    31   119   180
    174   199   232
    255   127    14
    255   187   120
    44   160    44
    152   223   138
    214    39    40
    255   152   150
    148   103   189
    197   176   213
    140    86    75
    196   156   148
    227   119   194
    247   182   210
    127   127   127
    199   199   199
    188   189    34
    219   219   141
    23   190   207
    158   218   229];

% Here they are
colors = colorList(lst,:);
end


text:
function MTEXmenu
% show up MTEX menu

try
  g = getGitInfo;
  v = ['MTEX (' g.branch ')'];
catch
  v = getMTEXpref('version');
end

disp(' ');
if isempty(javachk('desktop'))
  if verLessThan('matlab','7.13')
    s = v;
  else
    s = ['' v ''];
  end
  disp([ s  ...
    ' (show documentation)'])
  disp('  Import pole figure data')
  disp('  Import EBSD data')
  disp('  Import ODF data')
  disp(' ');

  if isappdata(0,'MTEXInstalled') && getappdata(0,'MTEXInstalled')
    disp('  Uninstall MTEX')
    disp(' ');
  else
    disp('  Install MTEX for future sessions')
    disp(' ');
  end

end

author: mtex-toolbox

repo: mtex


text:
cur_dir = pwd;
cd(fileparts(mfilename('fullpath')));

fprintf('Downloading Selective Search IJCV code...\n');
urlwrite('http://huppelen.nl/publications/SelectiveSearchCodeIJCV.zip', 'SelectiveSearchCodeIJCV.zip');

fprintf('Unzipping...\n');
unzip('SelectiveSearchCodeIJCV.zip');

fprintf('Done.\n');
system('del SelectiveSearchCodeIJCV.zip');

cd(cur_dir);

author: ShaoqingRen


text:
function [D,DJacob,DHess,pDpt]=DPolarLin2D(x,q0Dir,param4)
%%DPOLARCV2D The diffusion matrix for a continuous-time constant heading
%            motion model where the target state is given in terms of
%            position, direction of velocity in radians counterclockwise
%            from the x-axis, and speed.
%
%INPUTS: x The 4XN set of target state vectors in 2D space in the order of
%          [2D position;direction angle;speed] for each column where the
%          heading angle is measured in radians counterclockwise from the
%          x-axis. Alternatively, a 5XN state vector of
%          [2D position;direction angle;speed; speed derivative] can be
%          used to cover a linearlly accelerating target. Only the
%          dimensionality of x matters; the components are ignored. The
%          dimensionality is used to determine whether linear acceleration
%          is modeled and what N is.
%    q0Dir The power spectral density of the process noise for the heading
%          direction. The units are rad^2/s.
%   param4 If x is only 4D, that is, does not model a constant linear
%          acceleration, then param4 is the power spectral density of the
%          process noise for the speed. The units are m^2/s^3. Otherwise,
%          param4 is the power spectral density of the noise for the
%          linear acceleration and its units are m^2/s^5.
%
%OUTPUT: D The 4X2XN set of N diffusion matrices of a continuous-time
%          linear additive noise model where the velocity is given in terms
%          of direction and speed and the noise is only added to the
%          direction and speed components or the 5X2XN set of matrices
%          where noise is added to the direction and acceleration
%          components if a linear acceleration is assumed.
%   DJacob, DHess The xDimX2XxDim, xDimX2XxDimXxDim matrices of first and
%           second partial derivatives of the elements of D with respect to
%           x. These are all zero, since D is a constant. it is the same
%           for all D and is not repeated N times.
%      pDpt The xDimX2 partial derivative of D with respect to time. This
%           is all zeros, because D is a constant.
%
%Ideas for setting q0Dir  can be obtained using processNoiseSuggest with
%'PolyKal-ROT' and order 0. Ideas for setting param4 can
%be obtained using processNoiseSuggest with 'PolyKal-ROT' and order 1 or 2
%respectively depending on whether x is constant velocity or constant
%acceleration.
%
%The idea of decomposing the velocity into direction and speed components
%is presented in [1], where the direction component is differently defined
%and additional derivatives are used.
%
%This diffusion matrix goes with the drift function aPolarLin2D.
%
%REFERENCES:
%[1] J. L. Gertz, "Multisensor surveillance for improved aircraft
%    tracking," The Lincoln Laboratory Journal, vol. 2, no. 3, pp. 381-396,
%    1989.
%
%July 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

N=size(x,2);

q0DirR=sqrt(q0Dir);
param4R=sqrt(param4);

numDim=size(x,1);

if(numDim==4)%If it is a constant velocity model
    D=repmat([0,      0;
                0,      0;
                q0DirR, 0;
                0,     param4R],[1,1,N]);
elseif(numDim==5)%If it is a constant acceleration model.
    D=repmat([0,      0;
                0,      0;
                q0DirR, 0;
                0,      0;
                0,      param4R],[1,1,N]);
else
    error('The length of x is neither 4 nor 5.');
end

if(nargout>1)
    DJacob=zeros(numDim,2,numDim);
    if(nargout>2)
        DHess=zeros(numDim,2,numDim,numDim);
        if(nargout>3)
            pDpt=zeros(numDim,2);
        end
    end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function vf_out = get_vector_field(deformS,toPlanC,toScanNum,fromPlanC,fromScanNum)
% function vf = get_vector_field(deformS,toPlanC,toScanNum,fromPlanC,fromScanNum)
%
% APA, 09/24/2012

% Obtain base and moving scan UIDs
baseScanUID = deformS.baseScanUID;
movScanUID  = deformS.movScanUID;

% Figure out whether an inverse Vector field is required
indexToS = toPlanC{end};
indexFromS = fromPlanC{end};

toScanUID = toPlanC{indexToS.scan}(toScanNum).scanUID;
fromScanUID = fromPlanC{indexFromS.scan}(fromScanNum).scanUID;
if isequal(baseScanUID, fromScanUID)
    calc_inv_vf_flag = 0;
else
    calc_inv_vf_flag = 1;
end

% Create b-spline coefficients file
bspFileName = fullfile(getCERRPath,'ImageRegistration','tmpFiles',['bsp_coeffs_',baseScanUID,'_',movScanUID,'.txt']);
% bspFileName = fullfile(tempdir,'tmpFiles',['bsp_coeffs_',baseScanUID,'_',movScanUID,'.txt']);
success     = write_bspline_coeff_file(bspFileName,deformS.algorithmParamsS);

% Obtain Vf from b-splice coefficients
vfFileName = fullfile(getCERRPath,'ImageRegistration','tmpFiles',['vf_',baseScanUID,'_',movScanUID,'.mha']);
%vfFileName = fullfile(tempdir,'tmpFiles',['vf_',baseScanUID,'_',movScanUID,'.mha']);
system(['plastimatch xf-convert --input ',escapeSlashes(bspFileName), ' --output ', escapeSlashes(vfFileName), ' --output-type vf'])
%system(['plastimatch convert --xf ',escapeSlashes(bspFileName), ' --output-vf=', escapeSlashes(vfFileName)])
delete(bspFileName)

if calc_inv_vf_flag
    % Get dims, origin, spacing for toScan
    [uniformCT, uniformScanInfoS] = getUniformizedCTScan(0,toScanNum,toPlanC);
    uniformCT = permute(uniformCT, [2 1 3]);
    uniformCT = flipdim(uniformCT,3);
    % Change data type to int16 to allow (-)ve values
    uniformCT = int16(uniformCT) - int16(toPlanC{indexToS.scan}(toScanNum).scanInfo(1).CTOffset);
    % [dx, dy, dz]
    resolution = [uniformScanInfoS.grid2Units, uniformScanInfoS.grid1Units, uniformScanInfoS.sliceThickness] * 10;
    [xVals, yVals, zVals] = getUniformScanXYZVals(toPlanC{indexToS.scan}(toScanNum));
    offset = [xVals(1) -yVals(1) -zVals(end)] * 10;
    img_size = size(uniformCT);
    system(['vf_invert --input ', escapeSlashes(vfFileName), ' --output ', escapeSlashes(vfFileName), ' --dims="',num2str(img_size),'" --origin="',num2str(offset),'" --spacing="',num2str(resolution),'"'])
end

% infoS  = mha_read_header(vfFileName);
% vf = mha_read_volume(infoS);
[vf,infoS] = readmha(vfFileName);
%vf = flipdim(permute(vf,[2,1,3]),3);
delete(vfFileName)

vf_out(:,:,:,1) = flipdim(permute(vf(:,:,:,1),[2,1,3]),3)/10;
vf_out(:,:,:,2) = flipdim(permute(-vf(:,:,:,2),[2,1,3]),3)/10;
vf_out(:,:,:,3) = flipdim(permute(-vf(:,:,:,3),[2,1,3]),3)/10;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function NewGoal = GeneGoal(PopObj,NGoal)
% Generate new goals

%------------------------------- 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".
%--------------------------------------------------------------------------

    Gmax    = max(PopObj,[],1)*1.2;
    Gmin    = min(PopObj,[],1);
    NewGoal = unifrnd(repmat(Gmin,NGoal,1),repmat(Gmax,NGoal,1));
end

author: BIMK

repo: PlatEMO