text:
function [mCatDecluster, mCatAfter, vCluster, vCl, vMainCluster] = calc_decluster_ver1(mCatalog,nMethod)
% function [mCatDecluster, mCatAfter, vCluster, vCl, vMainCluster] = calc_decluster_ver1(mCatalog,nMethod)
% ----------------------------------------------------------------------------------------------------------
%
% Function to decluster earthquake catalog using the windowing technique in space and time by
% Knopoff & Gardner, GJR astr. Soc, 28, 311-313, 1972
% Gardner & Knopoff, BSSA, 64,5, 1363-1367, 1974
% using different windows
%
% Incoming variables
% mCatalog : Incoming earthquake catalog (ZMAP format)
% nMethod  : Window length for declustering (see calc_windows.m)
%            1: Gardener & Knopoff, 1974
%            2: Gruenthal pers. communication
%            3: Urhammer, 1986
%
% Outgoing variables:
% mCatDecluster : Declustered earthquake catalog
% mCatAfter     : Catalog of aftershocks (foreshocks)
% vCluster      : Vector indicating only aftershocks/foreshocls in cluster using a cluster number
% vCl           : Vector indicating all events in clusters using a cluster number
% vMainCluster  : Vector indicating mainshocks of clusters using a cluster number
%
% J. Woessner, woessner@seismo.ifg.ethz.ch
% last update: 31.07.02

%% Added:
% 31.07.02 Correction for problem of mainshocks with a cluster number as aftershocks belong to two sequences
% 31.07.02 Corrected fMaxClusterMag(nMagCount) to fMaxClusterMag, since counting variable not needed

%%% Remember: Improve zero length cluster problem, no waitbars yet, time difference mainshock-original event not implemented
%% This is the version before calc_decluster.m

%% Initialize Vectors and Matrices
mCatDecluster = [];
mCatAfter = [];
vCluster = zeros(length(mCatalog),1); % Initialize all events as mainshock
vCl = zeros(length(mCatalog),1); % Initialize all events as mainshock
vSel = zeros(length(mCatalog),1); % Initialize all events as mainshock
vMainCluster = zeros(length(mCatalog),1); % Initialize

[nXSize, nYsize] = size(mCatalog);
if nXSize == 0
    disp('Load new catalog');
    return
end

vDecDate = mCatalog(:,3);
nCount = 0;    % Variable of cluster number

fMagThreshold = min(mCatalog(:,6)); % Set Threshold to minimum magnitude of catalog
for nEvent=1:length(mCatalog(:,6))
    %nEvent
    if vCluster(nEvent) == 0
        fMagnitude(nEvent) = mCatalog(nEvent, 6);
        if fMagnitude(nEvent) >= fMagThreshold
            %% Define first aftershock zone and determine magnitude of strongest aftershock
            fMag = fMagnitude(nEvent);
            [fSpace, fTime] = calc_windows(fMagnitude(nEvent), nMethod);
            fSpaceDeg = km2deg(fSpace);
            %% This first if is for events with no location given
            if isnan(mCatalog(nEvent,1))
                %vSel = (vDecDate(:,1)-vDecDate(nEvent,1) >= 0) & (vDecDate(:,1)-vDecDate(nEvent,1) <= fTime  & vCluster(nEvent) == 0);
                vSel = (mCatalog(:,3) == mCatalog(nEvent,3));
            else
                vSel = (abs(mCatalog(:,1) - mCatalog(nEvent,1)) <= fSpaceDeg) & (abs(mCatalog(:,2) - mCatalog(nEvent,2)) <= fSpaceDeg)...
                    & (vDecDate(:,1)-vDecDate(nEvent,1) >= 0) & (vDecDate(:,1)-vDecDate(nEvent,1) <= fTime & vCluster(nEvent) == 0);
            end
            mTmp = mCatalog(vSel,:);
            %nMagCount = 1; % Variable to determine max. magnitude in a cluster
            if length(mTmp(:,1)) == 1  % Only one event thus no cluster
                fMaxClusterMag = fMag;
            else
                fMaxClusterMag = max(mTmp(:,6));
                % Search for event with bigger magnitude in cluster and resize windows if necessary
                while fMaxClusterMag-fMag > 0
                    [fSpace, fTime] = calc_windows(fMaxClusterMag, nMethod);
                    fSpaceDeg = km2deg(fSpace);
                    vSel = (abs(mCatalog(:,1) - mCatalog(nEvent,1)) <= fSpaceDeg) &...
                        (abs(mCatalog(:,2) - mCatalog(nEvent,2)) <= fSpaceDeg)...
                        & (vDecDate(:,1)-vDecDate(nEvent,1) >= 0) & (vDecDate(:,1)-vDecDate(nEvent,1) <= fTime & vSel == 1);
                    mTmp = mCatalog(vSel,:);
                    if isempty(mTmp) % no events in aftershock zone
                        break;
                    end
                    fMag = fMaxClusterMag;
                    %nMagCount = nMagCount+1;
                    fMaxClusterMag = max(mTmp(:,6));
                    if fMaxClusterMag - fMag == 0 % no bigger event in aftershock zone
                        break;
                    end
                end  % End of while
                nCount = nCount + 1; % Set cluster number
            end % End of if
            % Define final aftershock zone if different from original one
            if fMaxClusterMag > fMagnitude(nEvent)
                [fSpace, fTime] = calc_windows(fMaxClusterMag, nMethod);
                fSpaceDeg = km2deg(fSpace);
                % Select events with the Event having the number nEvent, since it is a foreshock (time condition)
                vSel = (abs(mCatalog(:,1) - mCatalog(nEvent,1)) <= fSpaceDeg) & (abs(mCatalog(:,2) - mCatalog(nEvent,2)) <= fSpaceDeg)...
                    & (vDecDate(:,1)-vDecDate(nEvent,1) >= 0) & (vDecDate(:,1)-vDecDate(nEvent,1) <= fTime & vCluster(nEvent) == 0);
            end
            [vIndice]=find(vSel); % Vector of indices with Clusters
            vTmpCluster(vIndice,:) = nCount;
            length(vTmpCluster(vIndice,:));
            nI=1; % Variable counting the length of the cluster
            % Select the right numbers for the cluster using the indice vector vIndice
            % First: Insert cluster number after check for length
            % Second: Check if it's a mainshock
            % Third: Keep the former cluster indice;
            while nI <= length(vIndice)
                if (~isempty(vTmpCluster(vIndice(nI))) & length(vTmpCluster(vIndice,:)) > 1 & vCluster(vIndice(nI)) == 0)
                    vCluster(vIndice(nI)) = vTmpCluster(vIndice(nI));
                    %vEventnr(vIndice,:) = nEvent;
                elseif  (~isempty(vTmpCluster(vIndice(nI)))  &&  length(vTmpCluster(vIndice,:)) == 1  &&  vCluster(vIndice(nI)) == 0)
                    vCluster(vIndice(nI)) = 0;
                else
                    vCluster(vIndice(nI)) = vCluster(vIndice(nI));
                end
                nI=nI+1;
            end %End of while nI
            %%% Check if the Cluster is not just one event which can happen in case of keeping the former
            %%% cluster number in preceeding while-Loop
            vSelSingle = (vCluster == nCount);
            [vIndiceSingle] = find(vSelSingle);
            %vTmpSingle(vIndiceSingle,:);
            if length(vIndiceSingle) == 1
                nCount
                vIndiceSingle
                vCluster(vIndiceSingle)=0; % Set the event as mainsock
                nCount = nCount-1; % Correct the cluster number down by one
            end
        end % End of if checking magnitude threshold fMagThreshold
    end % End of if chekcing if vCluster == 0
end % End of FOR over mCatalog
%nCount
%% vCL Cluster vector with mainshocks in it; vCluster is now modified to get rid of mainshocks
vCl = vCluster;

%% Matrix with cluster indice, magnitude and time
mTmpCat = [vCluster mCatalog(:,6) mCatalog(:,3)];
%% Delete largest event from cluster series and add to mainshock catalog
for nCevent = 1:nCount
    %nCevent
    vSel4 = (mTmpCat(:,1) == nCevent); % Select cluster
    mTmpCat2 = mCatalog(vSel4,:);
    fTmpMaxMag = max(mTmpCat2(:,6)); % Select max magnitude of cluster
    vSelMag = (mTmpCat2(:,6) == fTmpMaxMag);
    [nMag] = find(vSelMag);
    if length(nMag) == 1
        vSel5 = (mTmpCat(:,1) == nCevent & mTmpCat(:,2) == fTmpMaxMag); % Select the event
        [vIndiceMag] = find(vSel5); % Find indice
        vCluster(vIndiceMag) = 0;  % Set cluster value to zero, so it is a mainshock
        vMainCluster(vIndiceMag) = nCevent; % Set mainshock vector to cluster number
    elseif isempty(nMag)
        disp('Nothing in ')
        nCevent
    else
        vSel = (mTmpCat(:,1) == nCevent & mTmpCat(:,2) == fTmpMaxMag);
        mTmpCat3 = mCatalog(vSel,:);
        [vIndiceMag] = min(find(vSel)); % Find minimum indice of event with max magnitude in cluster
        vCluster(vIndiceMag) = 0;  % Set cluster value to zero, so it is a mainshock
        vMainCluster(vIndiceMag) = nCevent;  % Set mainshock vector to cluster number
    end
end % End of For nCevent

%% Create a catalog of aftershocks (mCatAfter) and of declustered catalog (mCatDecluster)
vSel = (vCluster(:,1) > 0);
mCatDecluster=mCatalog(~vSel,:);
mCatAfter = mCatalog(vSel,:);

author: CelsoReyes


text:
function mappedX = cfa(X, no_dims, no_analyzers, no_iterations)
%CFA Performs manifold charting on dataset X
%
%   mappedX = cfa(X, no_dims, no_analyzers, no_iterations)
%
% Performs manifold charting on dataset X to reduce its dimensionality to
% no_dims dimensions. The variable no_analyzers determines the number of
% local factor analyzers that is used in the mixture of factor analyzers
% (default = 40). The variable no_iterations sets the number of
% iterations that is employed in the EM algorithm (default = 200).
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    if ~exist('no_dims', 'var')
        no_dims = 2;
    end
    if ~exist('no_analyzers', 'var')
        no_analyzers = 40;
    end
    if ~exist('no_iterations', 'var')
        no_iterations = 200;
    end

    % Make sure data is zero-mean, unit-variance
    X = X - repmat(mean(X, 1), [size(X, 1) 1]);
    X = X ./ repmat(var(X, 1), [size(X, 1) 1]);

    % Initialize some parameters
    min_var = 1e-5;                     % minimum STD of Gaussians
    X = X';
    [D n] = size(X);
    c = no_analyzers;
    d = no_dims;

    % Randomly initialize the parameters of the CFA model
    SigmaN = repmat(eye(d), [1 1 n]);
    Z = rand(d, n) * .1;
    Pi = repmat(1 / c, [1 c]);
    Kappa = rand(d, c) * .1;
    Mu = randn(D, c) * .1;
    SigmaC = repmat(eye(d), [1 1 c]);
    Lambda = rand(D, d, c) * .1;
    Psi = zeros(D, D, c);
    for j=1:c
        tmp = zeros(D, D);
        tmp(1:size(tmp, 2) + 1:end) = rand(D, 1) * .05;
        Psi(:,:,j) = tmp + tmp';
    end

    % Perform the EM algorithm for the optimization
    iter = 0;
    while iter < no_iterations

        % E-step
        % ====================================================

        % Do some precomputations for speed
        invPsi = zeros(size(Psi));
        detSigmaC = zeros(1, c);
        detPsi = zeros(1, c);
        logPi = zeros(1, c);
        for j=1:c
            invPsi(:,:,j) = inv(Psi(:,:,j) + eye(D));
            detSigmaC(j) = det(SigmaC(:,:,j));
            detPsi(j) = det(Psi(:,:,j));
            logPi(j) = log(Pi(j));
        end

        % Compute matrices Epsilon, Vc, and m
        Eps = zeros(n, c);
        Vc = zeros(d, d, c);
        m = zeros(d, c);
        const = ((D + d) / 2) * log(2 * pi);
        for j=1:c

            % Precomputations
            Xnc = X - repmat(Mu(:,j), [1 n]);
            Znc = Z - repmat(Kappa(:,j), [1 n]);
            tmpProd = Lambda(:,:,j)' * invPsi(:,:,j) * Lambda(:,:,j);

            % Compute Epsilon
            for i=1:n
                normX = (Xnc(:,i) - Lambda(:,:,j) * Znc(:,i));
                Eps(i, j) = -logPi(j) + const ...
                            + (.5 * log(detSigmaC(j))) + (.5 * detPsi(j)) ...
                            + (.5 * trace(SigmaC(:,:,j) * (SigmaN(:,:,i) + Znc(:,i) * Znc(:,i)'))) ...
                            + (.5 * trace(SigmaN(:,:,i) * tmpProd)) ...
                            + (.5 * normX' * invPsi(:,:,j) * normX);
            end

            % Compute Vc and m
            Vc(:,:,j) = inv(SigmaC(:,:,j) + eps * eye(d)) + tmpProd;
            m(:,j) = Kappa(:,j) + (Vc(:,:,j) \ Lambda(:,:,j)') * invPsi(:,:,j) * mean(Xnc, 2);
        end

        % Update estimate of Q
        Q = (repmat(sum(exp(-Eps), 2), [1 c]) .^ -1) .* exp(-Eps);

        % Update estimate of SigmaN
        for i=1:n
            tmp = zeros(d, d);
            for j=1:c
                tmp = tmp + Q(i, j) * Vc(:,:,j);
            end

            % Compute covariance matrix
            SigmaN(:,:,i) = inv(tmp + eps * eye(d));           % code above gave us inv(SigmaN)
        end

        % Update estimate of Z (= mappedX)
        for i=1:n
            tmp = zeros(d, 1);
            for j=1:c
                tmp = tmp + (Q(i, j) * Vc(:,:,j) * m(:,j));
            end
            Z(:,i) = SigmaN(:,:,i) * tmp;
        end

        % M-step
        % ====================================================

        % Update estimate of Pi
        Pi = sum(Q, 1) ./ n;

        % Update estimate of Mu and Kappa
        for j=1:c
            tmpQ = Q(:,c) ./ sum(Q(:,c));
            Mu(:,c) = sum(repmat(tmpQ', [D 1]) .* X, 2);
            Kappa(:,c) = sum(repmat(tmpQ', [d 1]) .* Z, 2);
        end

        % Update estimate of SigmaC
        for j=1:c
            tmp = 0;
            tmpQ = Q(:,c) ./ sum(Q(:,c));
            for i=1:n
                Znc = Z(:,i) - Kappa(:,j);
                tmp = tmp + (tmpQ(i) * (SigmaN(:,:,i) + Znc * Znc'));
            end

            % Enforce some variance
            tmp(1:size(tmp, 1) + 1:end) = max(min_var, tmp(1:size(tmp, 1) + 1:end));
            SigmaC(:,:,j) = tmp;
        end

        % Update estimate of Lambda
        for j=1:c
            Sc = zeros(D, d);
            tmpQ = Q(:,j) ./ sum(Q(:,j));
            Xnc = X - repmat(Mu(:,j), [1 n]);
            Znc = Z - repmat(Kappa(:,j), [1 n]);
            for i=1:n
                Sc = Sc + tmpQ(i) * (Xnc(:,i) * Znc(:,i)');
            end
            Lambda(:,:,j) = Sc / (SigmaC(:,:,j) + eps * eye(d));
        end

        % Update estimate of Psi
        for j=1:c
            tmpPsi = zeros(D, D);
            tmpQ = Q(:,j) ./ sum(Q(:,j));
            Xnc = X - repmat(Mu(:,j), [1 n]);
            Znc = Z - repmat(Kappa(:,j), [1 n]);
            tmpProd = Lambda(:,:,j) * SigmaN(:,:,i) * Lambda(:,:,j)';
            for i=1:n
                tmpPsi(1:size(tmpPsi, 2) + 1:end) = tmpPsi(1:size(tmpPsi, 2) + 1:end) + ...
                    tmpQ(i) * (((Xnc(:,i) - Lambda(:,:,j) * Znc(:,i)) .^ 2)' + tmpProd(1:size(tmpProd, 2) + 1:end));
            end
            Psi(:,:,c) = tmpPsi;
        end

        % Update number of iterations
        iter = iter + 1;
        if rem(iter, 5) == 0
            fprintf('.');
        end
    end

    % Transpose to get lowdimensional data representation
    mappedX = Z';

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function isodoseClosedContour = getBeamCrossSection()
%isodoseClosedContour50 = getBeamCrossSection()
%
%Output: isodoseClosedContour - representation of beam contour
%dimension (3xN). 1st, 2nd and 3rd row represent x, y and z-coords of the
%contour respectively. It is assumed that the ray starts at (0,0,0) and
%moves along +(ve) x-axis. If collimatorType does not match the list, an
%empty isodoseClosedContour50 is returned.
%
%APA 01/21/09
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

% Get a 10x10 section at 100cm distance. In future, these points should
% come from Leaf-sequences or beamlet-fluence map.
isodoseClosedContour(3,:) = [-10 10 10 -10 -10]/2;
isodoseClosedContour(2,:) = [-10 -10 10 10 -10]/2;
isodoseClosedContour(1,:) = [100 100 100 100 100];
return

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef ParamsMonitor < ParamsMonitor_Interface

    properties (Access = private)
        frame
        figures
        nFigs
        nRows
        nCols

        dualVariable

        problemID
        cost
        nCostFuncs
        constraint
        nConstraints
        refreshInterval
        convergenceVars
        iteration
        hasFinished
        iRefresh
        iStep
        nStep

        namingManager
    end

    methods (Access = public)

        function obj = ParamsMonitor(cParams)
            obj.init(cParams);
            drawnow
        end

        function refresh(obj,it,hasFinished,iStep,nStep)
            obj.updateParams(it,hasFinished,iStep,nStep);
            obj.updateFigures();
            if obj.shallRefresh()
                obj.refreshFigures();
                obj.refreshFrameTitle();
            end
        end

        function f = saveFigure(obj)
            f = obj.frame;
            fNameFig = fullfile('Output',obj.problemID,'Monitoring.fig');
            fNamePng = fullfile('Output',obj.problemID,'Monitoring.png');
            saveas(f,fNameFig)
            saveas(f,fNamePng)
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.problemID = cParams.problemID;
            obj.refreshInterval = cParams.refreshInterval;
            obj.nCostFuncs      = length(cParams.costFuncNames);
            obj.nConstraints    = length(cParams.constraintFuncs);

            obj.dualVariable     = cParams.dualVariable;
            obj.cost             = cParams.cost;
            obj.constraint       = cParams.constraint;
            obj.convergenceVars  = cParams.convergenceVars;
            obj.createNamingManager(cParams);
            obj.createFigures();
        end

        function createFigures(obj)
            obj.figures = {};
            obj.initCostFigures();
            obj.initConstraintFigures();
            obj.initConvergenceVarsFigures();
            obj.initShapeFunctionsVariablesFigures();
            obj.displayFigures();
        end

        function initCostFigures(obj)
            obj.initFigure('Cost');
            for i = 1:obj.nCostFuncs
                title = obj.namingManager.getCostFuncFigureTitle(i);
                obj.initFigure(title);
            end
        end

        function initConstraintFigures(obj)
            for i = 1:obj.nConstraints
                obj.initConstraintFigure(i);
                obj.initLambdaFigure(i);
            end
        end

        function initConstraintFigure(obj,i)
            title = obj.namingManager.getConstraintFigureTitle(i);
            obj.initFigure(title);
        end

        function initLambdaFigure(obj,i)
            title = obj.namingManager.getLambdaFigureTitle(i);
            obj.initFigure(title);
        end

        function initConvergenceVarsFigures(obj)
            for i = 1:obj.convergenceVars.nVar
                title = obj.namingManager.getConvVarNames(i);
                obj.initFigure(title);
            end
        end

        function initShapeFunctionsVariablesFigures(obj)
            obj.initShapeFunctionsVariablesFromFunctional(obj.cost);
            obj.initShapeFunctionsVariablesFromFunctional(obj.constraint);
        end

        function initShapeFunctionsVariablesFromFunctional(obj,functional)
            nShapes = numel(functional.shapeFunctions);
            for i = 1:nShapes
                shape  = functional.shapeFunctions{i};
                titles = shape.getTitlesToPlot();
                for iv = 1:numel(titles)
                    title = titles{iv};
                    obj.initFigure(title);
                end
            end
        end

        function obj = initFigure(obj,title)
            chartType = obj.getChartType(title);
            newFig = DisplayFactory.create(chartType,title);
            obj.appendFigure(newFig);
        end

        function appendFigure(obj,fig)
            obj.figures{end+1} = fig;
        end

        function displayFigures(obj)
            obj.createFrame();
            obj.computeDistribution();
            for i = 1:obj.nFigs
                obj.figures{i}.show(obj.nRows,obj.nCols,i,[0.06 0.03]);
            end
        end

        function createFrame(obj)
            obj.frame = figure;
            drawnow;
            obj.frame.WindowState = "maximized";
        end

        function computeDistribution(obj)
            if obj.nFigs <= 4
                obj.nRows = 1;
            else
                obj.nRows = 2;
            end
            obj.nCols = round(obj.nFigs/obj.nRows);
        end

        function updateParams(obj,it,hasFinished,iStep,nStep)
            obj.iteration = it;
            obj.hasFinished = hasFinished;
            obj.iRefresh = 1;
            obj.iStep = iStep;
            obj.nStep = nStep;
        end

        function updateFigures(obj)
            obj.updateCost();
            obj.updateConstraints();
            obj.updateConvergenceVars();
            obj.updateShapeFunctionsVar();
        end

        function updateCost(obj)
            obj.updateFigure(obj.cost.value);
            for i = 1:obj.nCostFuncs
                shape = obj.cost.shapeFunctions{i};
                value = shape.value;
                obj.updateFigure(value);
            end
        end

        function updateConstraints(obj)
            for i = 1:obj.nConstraints
                shape = obj.constraint.shapeFunctions{i};
                value = shape.value;
                obj.updateFigure(value);
                if ~isempty(obj.dualVariable.value)
                    obj.updateFigure(obj.dualVariable.value(i))
                end
            end
        end

        function updateConvergenceVars(obj)
            for i = 1:obj.convergenceVars.nVar
                if ~isempty(obj.convergenceVars.values(i))
                    obj.updateFigure(obj.convergenceVars.values(i));
                end
            end
        end

        function updateShapeFunctionsVar(obj)
            obj.updateShapeFunctionsVarFromFunctional(obj.cost);
            obj.updateShapeFunctionsVarFromFunctional(obj.constraint);
        end

        function updateShapeFunctionsVarFromFunctional(obj,functional)
            nShapes = numel(functional.shapeFunctions);
            for i = 1:nShapes
                shape  = functional.shapeFunctions{i};
                values = shape.getVariablesToPlot();
                for iv = 1:numel(values)
                    value = values{iv};
                    obj.updateFigure(value);
                end
            end
        end

        function updateFigure(obj,value)
            obj.figures{obj.iRefresh}.updateParams(obj.iteration,value,obj);
        end

        function refreshFigures(obj)
            for iFig = 1:obj.nFigs
                obj.figures{iFig}.refresh();
            end
        end

        function refreshFrameTitle(obj)
            title = obj.namingManager.getFrameTitle(obj.iStep,obj.nStep,obj.iteration);
            set(obj.frame,'NumberTitle','off','Name',title)
        end

        function itShall = shallRefresh(obj)
            itShall = (obj.isTimeToRefresh() || obj.hasFinished) && obj.iteration > 0;
        end

        function itIs = isTimeToRefresh(obj)
            itIs = mod(obj.iteration,obj.refreshInterval) == 0;
        end

        function createNamingManager(obj,cParams)
           s.costFuncNames   = cParams.costFuncNames;
           s.costWeights     = cParams.costWeights;
           s.constraintFuncs = cParams.constraintFuncs;
           s.optimizerNames   = cParams.optimizerNames;
           obj.namingManager = NamingManager(s);
        end

    end

    methods (Access = private, Static)

        function type = getChartType(title)
            switch title
                case {'Line Search','outit','Line Search trials'}
                    type = 'bar';
                case {'kktnorm','Norm L2','inf_{pr}','inf_{du}','C - C not scaled'}
                    type = 'log';
                otherwise
                    type = 'plot';
            end
        end

    end

    methods (Access = ?Display_Abstract)

        function increaseRefreshIterator(obj)
            obj.iRefresh = obj.iRefresh + 1;
        end

    end

    methods

        function n = get.nFigs(obj)
            n = length(obj.figures);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% this function solves the eigenvalue problem omega=f(epsi,kx,ky) for
%%% in-plane propagation (i.e. z=0) in a 2D-PhC; the problem can be
%%% separated in two orthogonal polarizations, E-pol & H-pol
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [omegaE, omegaH]=eigsEH(kGx,kGy,epsi,N,bands)
eta=inv(epsi);
A=inv(eta*(kGx^2+kGy^2)); %%% matrix for E-pol
B=inv(kGx*eta*kGx+kGy*eta*kGy); %%% matrix for H-pol
%%% options used in "eigs" calculations - tolerance 1e-12, intermediate results of iterations will
%%% not be displayed
opts.tol=1e-12; opts.disp=0;
%%% calculate the largest eigs, which become the smallest after inversion
De=eigs(A,bands,'lr',opts);
Dh=eigs(B,bands,'lr',opts);
%%% eigenvalues "omega" are put in a column vector and sorted in ascending order
omegaE=sqrt(sort(1./De));
omegaH=sqrt(sort(1./Dh));

author: Sable

repo: mcbench-benchmarks


text:
function imb=regPutBorde(im,Nx,Ny,method);
% regPutBorde -  Puts a border to an image:
%	method=1 -> consider the image periodical
%	method=2 -> specular
%	method=3 -> repeat the edge pixel
%
%	imb = regPutBorde(im,Nx,Ny,method);
%
% ON - 10/96 (from putborde)
% Oscar Nestares - 5/99 (renamed to regPutBorde)
%

[sy sx]=size(im);
imb=zeros(sy+2*Ny,sx+2*Nx);
imb(1+Ny:sy+Ny,1+Nx:sx+Nx)=im;

if method == 1
	imb(1:Ny,1+Nx:sx+Nx)=im(sy-Ny+1:sy,:);
	imb(sy+Ny+1:sy+2*Ny,1+Nx:sx+Nx)=im(1:Ny,:);
	imb(1+Ny:sy+Ny,1:Nx)=im(:,sx-Nx+1:sx);
	imb(1+Ny:sy+Ny,sx+Nx+1:sx+2*Nx)=im(:,1:Nx);
	imb(1:Ny,1:Nx)=im(sy-Ny+1:sy,sx-Nx+1:sx);
	imb(1:Ny,sx+Nx+1:sx+2*Nx)=im(sy-Ny+1:sy,1:Nx);
	imb(sy+Ny+1:sy+2*Ny,1:Nx)=im(1:Ny,sx-Nx+1:sx);
	imb(sy+Ny+1:sy+2*Ny,sx+Nx+1:sx+2*Nx)=im(1:Ny,1:Nx);
elseif method == 2
	imb(Ny:-1:1,1+Nx:sx+Nx)=im(2:Ny+1,:);
	imb(sy+2*Ny:-1:sy+Ny+1,1+Nx:sx+Nx)=im(sy-Ny:sy-1,:);
	imb(1+Ny:sy+Ny,Nx:-1:1)=im(:,2:Nx+1);
	imb(1+Ny:sy+Ny,sx+2*Nx:-1:sx+Nx+1)=im(:,sx-Nx:sx-1);
	imb(Ny:-1:1,Nx:-1:1)=im(2:Ny+1,2:Nx+1);
	imb(Ny:-1:1,sx+2*Nx:-1:sx+Nx+1)=im(2:Ny+1,sx-Nx:sx-1);
	imb(sy+2*Ny:-1:sy+Ny+1,Nx:-1:1)=im(sy-Ny:sy-1,2:Nx+1);
	imb(sy+2*Ny:-1:sy+Ny+1,sx+2*Nx:-1:sx+Nx+1)=im(sy-Ny:sy-1,sx-Nx:sx-1);
elseif method==3
	for k=1:Nx
		imb(Ny+1:sy+Ny,k)=im(:,1);
		imb(Ny+1:sy+Ny,k+sx+Nx)=im(:,sx);
	end
	for k=1:Ny
		imb(k,Nx+1:sx+Nx)=im(1,:);
		imb(k+sy+Ny,Nx+1:sx+Nx)=im(sy,:);
	end
else
	error('Not a valid value for method')
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function pdf = f_pdf ( x, m, n )

%*****************************************************************************80
%
%% F_PDF evaluates the F central PDF.
%
%  Discussion:
%
%    PDF(X)(M,N) = M**(M/2) * X**((M-2)/2)
%      / ( Beta(M/2,N/2) * N**(M/2) * ( 1 + (M/N) * X )**((M+N)/2)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%    0.0 <= X
%
%    Input, integer M, N, the parameters of the PDF.
%    1 <= M,
%    1 <= N.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x < 0.0 )

    pdf = 0.0;

  else

    a = m;
    b = n;

    top = sqrt ( m^m * n^n * x^( m - 2 ) );
    bot1 = beta ( m / 2.0, n / 2.0 );
    bot2 =  sqrt ( ( n + m * x )^( m + n ) );

    pdf = top / ( bot1 * bot2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ a, b ] = p06_ab ( m )

%*****************************************************************************80
%
%% P06_AB returns bounds for problem 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Harald Niederreiter, Kevin McCurley,
%    Optimization of functions by quasi-random search methods,
%    Computing,
%    Volume 22, Number 2, 1979, pages 119-123.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Output, real A(M,1), B(M,1), lower and upper bounds.
%
  a(1:m,1) = 0.0;
  b(1:m,1) = 1.0;

  return
end


text:
function t_qps_matpower(quiet)
%T_QPS_MATPOWER  Tests of QPS_MATPOWER QP solvers.

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

if nargin < 1
    quiet = 0;
end

algs = {'DEFAULT', 'BPMPD', 'MIPS', 250, 'IPOPT', 'OT', 'CPLEX', 'MOSEK', 'GUROBI', 'CLP', 'GLPK'};
names = {'DEFAULT', 'BPMPD_MEX', 'MIPS', 'sc-MIPS', 'IPOPT', 'linprog/quadprog', 'CPLEX', 'MOSEK', 'Gurobi', 'CLP', 'glpk'};
check = {[], 'bpmpd', [], [], 'ipopt', 'quadprog', 'cplex', 'mosek', 'gurobi', 'clp', 'glpk'};
does_qp = [1 1 1 1 1 1 1 1 1 1 0];

n = 36;
nqp = 28;
t_begin(n*length(algs), quiet);

diff_alg_warn_id = 'optim:linprog:WillRunDiffAlg';
if have_feature('quadprog') && have_feature('quadprog', 'vnum') == 7.005
    s1 = warning('query', diff_alg_warn_id);
    warning('off', diff_alg_warn_id);
end

for k = 1:length(algs)
    if ~isempty(check{k}) && ~have_feature(check{k})
        t_skip(n, sprintf('%s not installed', names{k}));
    else
        opt = struct('verbose', 0, 'alg', algs{k});
        mpopt = struct( ...
            'verbose', 0, ...
            'opf', struct( ...
                'violation', 5e-6 ), ...
            'cplex', struct( ...
                'lpmethod', 0, ...
                'qpmethod', 0, ...
                'opt', 0 ), ...
            'mosek', struct( ...
                'lp_alg', 0, ...
                'max_it', 0, ...
                'gap_tol', 0, ...
                'max_time', 0, ...
                'num_threads', 0, ...
                'opt', 0 ) ...
        );
        if strcmp(names{k}, 'DEFAULT') || strcmp(names{k}, 'MIPS') || strcmp(names{k}, 'sc-MIPS')
            opt.mips_opt.comptol = 1e-8;
        end
%         if strcmp(names{k}, 'linprog/quadprog')
%             opt.verbose = 2;
%             opt.linprog_opt.Algorithm = 'interior-point';
%             opt.linprog_opt.Algorithm = 'active-set';
%             opt.linprog_opt.Algorithm = 'simplex';
%             opt.linprog_opt.Algorithm = 'dual-simplex';
%         end
        if strcmp(names{k}, 'CPLEX')
%           alg = 0;        %% default uses barrier method with NaN bug in lower lim multipliers
            alg = 2;        %% use dual simplex
            mpopt.cplex.lpmethod = alg;
            mpopt.cplex.qpmethod = min([4 alg]);
            opt.cplex_opt = cplex_options([], mpopt);
        end
        if strcmp(names{k}, 'MOSEK')
%             sc = mosek_symbcon;
%             alg = sc.MSK_OPTIMIZER_DUAL_SIMPLEX;    %% use dual simplex
%             alg = sc.MSK_OPTIMIZER_INTPNT;          %% use interior point
%             mpopt.mosek.lp_alg = alg;
            mpopt.mosek.gap_tol = 1e-10;
%             mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_PFEAS = 1e-10;
%             mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_DFEAS = 1e-10;
%             mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_INFEAS = 1e-10;
%             mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_REL_GAP = 1e-10;
            vnum = have_feature('mosek', 'vnum');
            if vnum >= 8
%                 mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_PFEAS = 1e-10;
%                 mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_DFEAS = 1e-10;
%                 mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_INFEAS = 1e-10;
%                 mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_MU_RED = 1e-10;
                mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_REL_GAP = 1e-10;
            end
%             opt.verbose = 3;
            opt.mosek_opt = mosek_options([], mpopt);
        end

        t = sprintf('%s - 3-d LP : ', names{k});
        %% based on example from 'doc linprog'
        c = [-5; -4; -6];
        A = [1 -1  1;
             -3  -2  -4;
             3  2  0];
        l = [-Inf; -42; -Inf];
        u = [20; Inf; 30];
        xmin = [0; 0; 0];
        x0 = [];
        [x, f, s, out, lam] = qps_matpower([], c, A, l, u, xmin, [], [], opt);
        t_is(s, 1, 12, [t 'success']);
        t_is(x, [0; 15; 3], 6, [t 'x']);
        t_is(f, -78, 6, [t 'f']);
        t_is(lam.mu_l, [0;1.5;0], 9, [t 'lam.mu_l']);
        t_is(lam.mu_u, [0;0;0.5], 9, [t 'lam.mu_u']);
        if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
            t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
        else
            t_is(lam.lower, [1;0;0], 9, [t 'lam.lower']);
            t_is(lam.upper, zeros(size(x)), 9, [t 'lam.upper']);
        end

        if does_qp(k)
            t = sprintf('%s - unconstrained 3-d quadratic : ', names{k});
            %% from http://www.akiti.ca/QuadProgEx0Constr.html
            H = [5 -2 -1; -2 4 3; -1 3 5];
            c = [2; -35; -47];
            x0 = [0; 0; 0];
            [x, f, s, out, lam] = qps_matpower(H, c, [], [], [], [], [], [], opt);
            t_is(s, 1, 12, [t 'success']);
            t_is(x, [3; 5; 7], 8, [t 'x']);
            t_is(f, -249, 11, [t 'f']);
            t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
            t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
            t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
            t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);

            t = sprintf('%s - constrained 2-d QP : ', names{k});
            %% example from 'doc quadprog'
            H = [   1   -1;
                    -1  2   ];
            c = [-2; -6];
            A = [   1   1;
                    -1  2;
                    2   1   ];
            l = [];
            u = [2; 2; 3];
            xmin = [0; 0];
            x0 = [];
            [x, f, s, out, lam] = qps_matpower(H, c, A, l, u, xmin, [], x0, opt);
            t_is(s, 1, 12, [t 'success']);
            t_is(x, [2; 4]/3, 7, [t 'x']);
            t_is(f, -74/9, 6, [t 'f']);
            t_is(lam.mu_l, [0;0;0], 13, [t 'lam.mu_l']);
            t_is(lam.mu_u, [28;4;0]/9, 4, [t 'lam.mu_u']);
            if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
                t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
            else
                t_is(lam.lower, zeros(size(x)), 7, [t 'lam.lower']);
                t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
            end

            t = sprintf('%s - constrained 4-d QP : ', names{k});
            %% from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm
            H = [   1003.1  4.3     6.3     5.9;
                    4.3     2.2     2.1     3.9;
                    6.3     2.1     3.5     4.8;
                    5.9     3.9     4.8     10  ];
            c = zeros(4,1);
            A = [   1       1       1       1;
                    0.17    0.11    0.10    0.18    ];
            l = [1; 0.10];
            u = [1; Inf];
            xmin = zeros(4,1);
            x0 = [1; 0; 0; 1];
            [x, f, s, out, lam] = qps_matpower(H, c, A, l, u, xmin, [], x0, opt);
            t_is(s, 1, 12, [t 'success']);
            t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
            t_is(f, 3.29/3, 6, [t 'f']);
            t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
            t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
            if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
                t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
            else
                t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
                t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
            end

            t = sprintf('%s - (struct) constrained 4-d QP : ', names{k});
            p = struct('H', H, 'A', A, 'l', l, 'u', u, 'xmin', xmin, 'x0', x0, 'opt', opt);
            [x, f, s, out, lam] = qps_matpower(p);
            t_is(s, 1, 12, [t 'success']);
            t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
            t_is(f, 3.29/3, 6, [t 'f']);
            t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
            t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
            if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
                t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
            else
                t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
                t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
            end
        else
            t_skip(nqp, sprintf('%s does not handle QP problems', names{k}));
        end

        t = sprintf('%s - infeasible LP : ', names{k});
        p = struct('A', sparse([1 1]), 'c', [1;1], 'u', -1, 'xmin', [0;0], 'opt', opt);
        [x, f, s, out, lam] = qps_matpower(p);
        t_ok(s <= 0, [t 'no success']);
    end
end

if have_feature('quadprog') && have_feature('quadprog', 'vnum') == 7.005
    warning(s1.state, diff_alg_warn_id);
end

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% plot masterthres sweeps
load('C:\Users\Xiu\Dropbox\FishExplorer2\AK Test Scripts\mydata.mat');

%% number of clusters
m = masterThresh_data;
data = zeros(size(m));
for i = 1:size(m,1),
    for j = 1:size(m,2),
        data(i,j) = m(i,j).nClus;
    end
end

figure('Position',[500,400,250,200]);
plot(0.5:0.05:0.9,data)
ylabel('# of clusters')
xlabel('clustering threshold (~correlation)')
xlim([0.5,0.9])
ylim([0,360])
title('custom tuning of stringency')

%% CV
data = zeros(size(m));
for i = 1:size(m,1),
    for j = 1:size(m,2),
        data(i,j) = m(i,j).CVscore;
    end
end

figure('Position',[500,400,300,200]);
plot(0.5:0.05:0.9,data)
ylabel('cross-val. score')
xlabel('clustering thresh. (~corr.)')
xlim([0.5,0.9])
ylim([0,1])
title('cross-val. (overlapping cell %)')

%% number of cells included
data = zeros(size(m));
for i = 1:size(m,1),
    for j = 1:size(m,2),
        data(i,j) = m(i,j).nCells;
    end
end

figure('Position',[500,400,300,200]);
plot(0.5:0.05:0.9,data)
ylabel('# of cells')
xlabel('clustering thresh. (~corr.)')
xlim([0.5,0.9])
title('total cell # in clusters')

author: xiuyechen

repo: FishExplorer


text:
function [predict_label, predict_prob] = CNN_SVM_prediction(imdb,OP,retrain_net,nfeatures,data_mean)
%
% predict_label:
%   if OP.class = 1 then: predict_label = 1:Wake, 2:REM sleep, 3:NREM light sleep, 4:NREM deep sleep
%   if OP.class = 2 then: predict_label = 1:NREM sleep (light + deep), 2:REM sleep, 3: wake
%   if OP.class = 3 then: predict_label = 1:NREM deep sleep, 2:NREM light sleep, 3:Wake + REM sleep
%   if OP.class = 4 then: predict_label = 1:NREM sleep (light + deep), 2: Wake + REM sleep

if nargin<5
    data_mean=[];
end

if nargin<4 || isempty(nfeatures)
    nfeatures=13;
end

imdb.meta.features(isinf(imdb.meta.features))=NaN;
sqi=nanmean(imdb.meta.sqi,1);
sqi(find(isnan(sqi)))=nanmean(sqi);

% load model
load(OP.model);

if nargin<3 || isempty(retrain_net)
    net=stat.net;
else
    net=retrain_net;
end

net.layers{1,end}.type='softmax';
mean_features=stat.mean_features;
std_features=stat.std_features;

mean_sqi=stat.mean_sqi;
std_sqi=stat.std_sqi;

model=stat.model;

% standardization
for i=8:nfeatures
    for j=1:size(imdb.meta.features,1)
        imdb.meta.features(j,i)=log(imdb.meta.features(j,i));
    end
end
imdb.meta.features(find(isinf(imdb.meta.features)))=NaN;

% subtract data_mean
if ~isempty(data_mean)
    stat.data_mean=data_mean;
end

for k=1:size(imdb.images.data,4)
    imdb.images.data(:,:,1,k)=imdb.images.data(:,:,1,k)-stat.data_mean;
end

for i=1:nfeatures
    imdb.meta.features(find(isnan(imdb.meta.features(:,i))),i)=mean_features(i);
    imdb.meta.features(:,i)=bsxfun(@minus, imdb.meta.features(:,i),mean_features(i));
    imdb.meta.features(:,i)=bsxfun(@rdivide, imdb.meta.features(:,i),std_features(i));
end
sqi=bsxfun(@minus,sqi,mean_sqi);
sqi=bsxfun(@rdivide,sqi,std_sqi);

% CNN prediction
if size(imdb.images.data,4)>1000
    split=round(size(imdb.images.data,4)/1000);
else
    split=1;
end

retall=[];
for kk=1:split
    each=ceil(size(imdb.images.data,4)/split);
    res=vl_simplenn(net,imdb.images.data(:,:,:,(kk-1)*each+1:min(kk*each,size(imdb.images.data,4))));

    ret=[];
    for i=1:OP.nclass
        for j=1:size(res(end).x,4)
            ret(i,j)=res(end).x(1,1,i,j);
        end
    end
    retall=[retall,ret];
end
% probability of CNN output
ret=retall;

xtest=[];
for i=1:size(ret,1)
    xtest(:,end+1)=ret(i,:);
end
for i=1:nfeatures
    if i~=6 % remove DFA feature
        xtest(:,end+1)=imdb.meta.features(:,i);
    end
end
xtest(:,end+1)=sqi;

% SVM prediction
[predict_label, accuracy_test, predict_prob] = svmpredict(ones(size(xtest,1),1), xtest, model,'-b 1');

[~,si]=sort(model.Label);
predict_prob=predict_prob(:,si);
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function F = coeffs2spherefunv( X, Y, Z )
%COEFFS2SPHEREFUNV   Convert componentwise matrices of 2D Fourier
%   coefficients to a spherefunv.
%
%   F = coeffs2spherefunv( X, Y, Z ) returns a spherefunv object F that has
%   matrices of 2D Fourier coefficients X, Y, and Z for each component.
%
% See also SPHEREFUN/COEFFS2SPHEREFUN

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

U = spherefun.coeffs2spherefun(X);
V = spherefun.coeffs2spherefun(Y);
W = spherefun.coeffs2spherefun(Z);
F = spherefunv(U, V, W);
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function varargout = departDisplayGUI(varargin)
% DEPARTDISPLAYGUI MATLAB code for departDisplayGUI.fig
%      DEPARTDISPLAYGUI, by itself, creates a new DEPARTDISPLAYGUI or raises the existing
%      singleton*.
%
%      H = DEPARTDISPLAYGUI returns the handle to a new DEPARTDISPLAYGUI or the handle to
%      the existing singleton*.
%
%      DEPARTDISPLAYGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in DEPARTDISPLAYGUI.M with the given input arguments.
%
%      DEPARTDISPLAYGUI('Property','Value',...) creates a new DEPARTDISPLAYGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before departDisplayGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to departDisplayGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help departDisplayGUI

% Last Modified by GUIDE v2.5 02-Dec-2013 19:58:32

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @departDisplayGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @departDisplayGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before departDisplayGUI is made visible.
function departDisplayGUI_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to departDisplayGUI (see VARARGIN)

% Choose default command line output for departDisplayGUI
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

%Possible fix for people with display issues.
% checkForCharUnitsInGUI(hObject);

% UIWAIT makes departDisplayGUI wait for user response (see UIRESUME)
% uiwait(handles.departDispGUI);

% --- Outputs from this function are returned to the command line.
function [dummy, varargout] = departDisplayGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

dummy{1} = handles.departAxis;
dummy{2} = handles.hyperOrbitText;
dummy{3} = handles.mainOrbitText;
dummy{4} = handles.burnInfoText;
dummy{5} = handles.departDispGUI;
dummy{6} = handles.hyperOrbitLabel;
dummy{7} = handles.transOrbitLabel;
dummy{8} = handles.departureOrbitRadio;

function hyperOrbitText_Callback(hObject, eventdata, handles)
% hObject    handle to hyperOrbitText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of hyperOrbitText as text
%        str2double(get(hObject,'String')) returns contents of hyperOrbitText as a double

% --- Executes during object creation, after setting all properties.
function hyperOrbitText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to hyperOrbitText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function burnInfoText_Callback(hObject, eventdata, handles)
% hObject    handle to burnInfoText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of burnInfoText as text
%        str2double(get(hObject,'String')) returns contents of burnInfoText as a double

% --- Executes during object creation, after setting all properties.
function burnInfoText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to burnInfoText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function mainOrbitText_Callback(hObject, eventdata, handles)
% hObject    handle to mainOrbitText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of mainOrbitText as text
%        str2double(get(hObject,'String')) returns contents of mainOrbitText as a double

% --- Executes during object creation, after setting all properties.
function mainOrbitText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to mainOrbitText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes when selected object is changed in plotDisplayPanel.
function plotDisplayPanel_SelectionChangeFcn(hObject, eventdata, handles)
% hObject    handle to the selected object in plotDisplayPanel
% eventdata  structure with the following fields (see UIBUTTONGROUP)
%	EventName: string 'SelectionChanged' (read only)
%	OldValue: handle of the previously selected object or empty if none was selected
%	NewValue: handle of the currently selected object
% handles    structure with handles and user data (see GUIDATA)
hPlotBtnGrp = findobj(handles.departDispGUI,'Tag','plotDisplayPanel');
hTypeBtn = get(hPlotBtnGrp,'SelectedObject');
plotType = get(hTypeBtn, 'Tag');

departAxis = handles.departAxis;
axisUserData = get(departAxis,'UserData');
reset(departAxis);
set(departAxis,'UserData',axisUserData);
hLegend = findobj(handles.departDispGUI,'Type','axes','Tag','legend');
delete(hLegend);

switch plotType
    case 'departureOrbitRadio'
        plotBodyDepartOrbit(departAxis)
    case 'transferOrbitRadio'
        plotTransferOrbit(departAxis)
end

% --------------------------------------------------------------------
function generateReport_Callback(hObject, eventdata, handles)
% hObject    handle to generateReport (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[hReport, handlesR] = genericReport(1, 'SINGLE DELTA-V MANEUVER', 'BALLISTIC TRANSFER APPLICATION');

    set(handlesR.rightTextBoxLabel,'String', 'Transfer Orbit Details');

    results = get(handles.departDispGUI,'UserData');
    results = results{1};

    paddLen = 29;
    form = '%9.5f';
    burnInfoTextstr = {};
    burnInfoTextstr{end+1} = [paddStr('Total Delta-V = ',paddLen), num2str(norm(results.burnVect), form), ' km/s'];
    burnInfoTextstr{end+1} = [paddStr('Prograde Delta-V = ',paddLen), num2str(1000*results.burnVect(1), form), ' m/s'];
    burnInfoTextstr{end+1} = [paddStr('Radial Delta-V = ',paddLen), num2str(1000*results.burnVect(3), form), ' m/s'];
    burnInfoTextstr{end+1} = [paddStr('Orbit Normal Delta-V = ',paddLen), num2str(1000*results.burnVect(2), form), ' m/s'];
    burnInfoTextstr{end+1} = '---------------------';
    burnInfoTextstr{end+1} = [paddStr('Departure True Anomaly = ',paddLen), num2str(rad2deg(results.burnTA), form), ' deg'];

    if(isfield(results,'burnTimeBeforePeri'))
        burnInfoTextstr{end+1} = [paddStr('Burn Time Before Peri. = ',paddLen), num2str(results.burnTimeBeforePeri, form), ' sec'];
    end

    if(isfield(results,'burnTimePastPeri'))
        burnInfoTextstr{end+1} = [paddStr('Burn Time Past Peri. = ',paddLen), num2str(results.burnTimeBeforePeri, form), ' sec'];
    end

    if(isfield(results,'iniOrbitPeriod'))
        burnInfoTextstr{end+1} = [paddStr('Ini. Orbit Period = ',paddLen), num2str(results.iniOrbitPeriod, form), ' sec'];
    end
    set(handlesR.leftTextBox, 'String', burnInfoTextstr);

    paddLen = 30;
    xferOrbitText = {};
    xferOrbitText{end+1} = ['Transfer Orbit about ', cap1stLetter(results.xferOrbitCBName)];
    xferOrbitText{end+1} = '-----------------------------------------------';
    xferOrbitText{end+1} = [paddStr('Semi-major Axis = ',paddLen), num2str(results.xferOrbit(1), form), ' km'];
    xferOrbitText{end+1} = [paddStr('Eccentricity = ', paddLen), num2str(results.xferOrbit(2))];
    xferOrbitText{end+1} = [paddStr('Inclination = ',paddLen), num2str(rad2deg(AngleZero2Pi(results.xferOrbit(3))), form), ' deg'];
    xferOrbitText{end+1} = [paddStr('Right Ascension of AN = ',paddLen), num2str(rad2deg(AngleZero2Pi(results.xferOrbit(4))), form), ' deg'];
    xferOrbitText{end+1} = [paddStr('Argument of Periapse = ',paddLen), num2str(rad2deg(AngleZero2Pi(results.xferOrbit(5))), form), ' deg'];
    xferOrbitText{end+1} = [paddStr([cap1stLetter(results.departBodyName), ' Depart True Anomaly = '],paddLen), num2str(rad2deg(AngleZero2Pi(results.xferOrbit(6))), form), ' deg'];
    xferOrbitText{end+1} = [paddStr([cap1stLetter(results.arrivalBodyName), ' Arrive True Anomaly = '],paddLen), num2str(rad2deg(AngleZero2Pi(results.xferOrbit(7))), form), ' deg'];

    set(handlesR.rightTextBox, 'String', xferOrbitText);

    set(handlesR.leftAxis,'UserData',results.axisUserData);
    set(handlesR.rightAxis,'UserData',results.axisUserData);

    plotBodyDepartOrbit(handlesR.leftAxis);
    plotTransferOrbit(handlesR.rightAxis, false);

% --------------------------------------------------------------------
function uploadManeuver_Callback(hObject, eventdata, handles)
% hObject    handle to uploadManeuver (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
guiData = get(handles.burnInfoText,'UserData');
uploadManeuverToKSP(guiData);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Fig. 7.2  Feedback Control of Dynamic Systems, 6e
%            Franklin, Powell, Emami
% script to generate Fig. 7.2

clear all;
close all;

F = [0 1; 0 -0.05];
G = [0; 0.001];
H = [0 1];
J = 0;
sys=ss(F,500*G,H,J);
step(sys);
title('Fig. 7.2: Response of car velocity to a step in u');
nicegrid;


text:
function lines = month_cal_store_common ( y, m )

%*****************************************************************************80
%
%% MONTH_CAL_STORE_COMMON stores a Common month calendar.
%
%  Discussion:
%
%    The "common" calendar is meant to be the calendar which is Julian before
%    the transition date, and Gregorian afterwards, with the transition date
%    best specified as as JED = 2299160.
%
%  Format:
%
%           1  2  3  4  5
%     6  7  8  9 10 11 12
%    13 14 15 16 17 18 19
%    20 21 22 23 24 25 26
%    27 28 29 30
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y,  M, the YM date.
%
%    Output, string LINES(6,20), the lines of the calendar.
%
  lines = [];
%
%  Make local copies of the input.
%
  m2 = m;
  y2 = y;
%
%  Check the month and year.  After this call, month is
%  guaranteed to be between 1 and 12.
%
  [ y2, m2, ierror ] = ym_check_common ( y2, m2 );

  if ( ierror ~= 0 )
    return
  end
%
%  Find the day of the week for Y M 1.
%
  d2 = 1;
  f2 = 0.0;

  w = ymdf_to_weekday_common ( y2, m2, d2, f2 );

  days = month_length_common ( y2, m2 );
%
%  Find the appropriate label for the first box in the calendar.
%
  iday = 2 - w;
%
%  Print out a line of day numbers.
%  IDAY keeps track of the numerical day,
%  JDAY keeps track of the label for the day, which differed in October 1582.
%
  d2 = iday;
  f2 = 0.0;

  n_line = 0;

  while ( n_line < 6 )
% while ( iday <= days )

    n_line = n_line + 1;
    s = '                    ';

    for w = 1 : 7

      i1 = 3 * ( w - 1 ) + 1;
      i2 = i1 + 1;

      if ( 1 <= iday && iday <= days )
        s(i1:i2) = sprintf ( '%2d', d2 );
      else
        s(i1:i2) = '  ';
      end

      iday = iday + 1;

      [ y2, m2, d2, f2 ] = ymdf_next_common ( y2, m2, d2, f2 );

    end

    lines = [ lines; s ];

  end

  return
end


text:
classdef TestRollingGuidanceFilter
    %TestRollingGuidanceFilter

    methods (Static)
        function test_1
            img = cv.imread(fullfile(mexopencv.root(),'test','lena.jpg'), ...
                'Color',true, 'ReduceScale',2);
            dst = cv.rollingGuidanceFilter(img, 'SigmaSpace',3.0);
            validateattributes(dst, {class(img)}, {'size',size(img)});
        end

        function test_error_argnum
            try
                cv.rollingGuidanceFilter();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function varargout = chebcoeffs2(varargin)
%CHEBCOEFFS2    Bivariate expansion coefficients
%   X = CHEBCOEFFS2(F) returns the matrix of bivariate coefficients such that
%       F= sum_{i=0}^{n-1} ( sum_{j=0}^{n-1} X(i+1,j+1) T_i(r) T_j(theta) ).
%   where -1 <= r <= 1 and -pi <= theta <= pi are the polar and radial
%   coordinates, respectively.
%
%   [A, D, B] = CHEBCOEFFS2( f ) returns the same coefficients keeping them in
%   low rank form, i.e., X = A * D * B'.
%
%
% See also DISKFUN/PLOTCOEFFS2, DISKFUN/COEFFS2.

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

[varargout{1:nargout}] = chebcoeffs2@separableApprox(varargin{:});

end


text:
function a = mass_matrix_t3 ( node_num, element_num, element_node, x, y )

%*****************************************************************************80
%
%% MASS_MATRIX_T3 computes the mass matrix, using 3-node triangles.
%
%  Discussion:
%
%    The mass matrix to be estimated has the form:
%
%      A(I,J) = integral ( PHI(I)(X,Y) * PHI(J)(X,Y) ) d Region
%
%    where PHI(I) and PHI(J) are the shape functions associated with
%    the I-th and J-th variables.
%
%  Referene Element T6:
%
%    |
%    1  3
%    |  |\
%    |  | \
%    S  |  \
%    |  |   \
%    |  |    \
%    0  1-----2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(3,ELEMENT_NUM), the nodes that make
%    up each element.
%
%    Input, real X(NODE_NUM), Y(NODE_NUM), the coordinates
%    of the nodes.
%
%    Output, real A(NODE_NUM,NODE_NUM), the mass matrix.
%
  element_order = 3;
%
%  Zero out the matrix.
%
  a(1:node_num,1:node_num) = 0.0;
%
%  Get the weights and abscissas for a unit triangle.
%
  rule = 4;
  quad_num = triangle_unit_size ( rule );
  [ rtab, stab, weight ] = triangle_unit_set ( rule );
%
%  For each element.
%
  for element = 1 : element_num

    p1 = element_node(1,element);
    p2 = element_node(2,element);
    p3 = element_node(3,element);

    area = 0.5 * abs ( ...
        x(p1) * ( y(p2) - y(p3) ) ...
      + x(p2) * ( y(p3) - y(p1) ) ...
      + x(p3) * ( y(p1) - y(p2) ) );

    if ( area == 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'MASS_MATRIX_T3 - Fatal error!\n' );
      fprintf ( 1, '  Zero area for element %d\n', element );
      error ( 'MASS_MATRIX_T3 - Fatal error!' );
    end
%
%  For each quadrature point in the element...
%
    for quad = 1 : quad_num

      r = rtab(quad);
      s = stab(quad);

      [ w, dwdr, dwds ] = shape_t3 ( r, s );
%
%  For each basis function PHI(I) associated with a node in the element,
%
      for iq = 1 : element_order

        ip = element_node(iq,element);
%
%  For each "neighbor" basis function PHI(J) associated with a node in
%  the element.
%
        for jq = 1 : element_order

          jp = element_node(jq,element);

          a(ip,jp) = a(ip,jp) + area * weight(quad) * w(iq) * w(jq);

        end
      end
    end
  end

  return
end


text:
%% Mixture of linear regression
close all; clear
d = 1;
k = 2;
n = 500;
[X,y] = mixLinRnd(d,k,n);
plot(X,y,'.');
[label,model,llh] = mixLinReg(X, y, k);
plotClass([X;y],label);
figure
plot(llh);
[y_,z,p] = mixLinPred(model,X,y);
figure;
plotClass([X;y],label);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
% sdae - training a stacked DAE (finetuning)
% Copyright (C) 2011 KyungHyun Cho, Tapani Raiko, Alexander Ilin
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
%
function [S] = sdae(S, patches, valid_patches, valid_portion);

if nargin < 3
    early_stop = 0;
    valid_patches = [];
    valid_portion = 0;
else
    early_stop = 1;
    valid_err = -Inf;
    valid_best_err = -Inf;
end

actual_lrate = S.learning.lrate;

n_samples = size(patches, 1);

layers = S.structure.layers;
n_layers = length(layers);

if layers(1) ~= size(patches, 2)
    error('Data is not properly aligned');
end

minibatch_sz = S.learning.minibatch_sz;
n_minibatches = ceil(n_samples / minibatch_sz);

n_epochs = S.iteration.n_epochs;

momentum = S.learning.momentum;
weight_decay = S.learning.weight_decay;

biases_grad = cell(n_layers, 1);
W_grad = cell(n_layers, 1);
biases_grad_old = cell(n_layers, 1);
W_grad_old = cell(n_layers, 1);
for l = 1:n_layers
    biases_grad{l} = zeros(size(S.biases{l}))';
    if l < n_layers
        W_grad{l} = zeros(size(S.W{l}));
    end
    biases_grad_old{l} = zeros(size(S.biases{l}))';
    if l < n_layers
        W_grad_old{l} = zeros(size(S.W{l}));
    end
end

min_recon_error = Inf;
min_recon_error_update_idx = 0;
stopping = 0;

do_normalize = S.do_normalize;
do_normalize_std = S.do_normalize_std;

if S.data.binary == 0
    if do_normalize == 1
        % make it zero-mean
        patches_mean = mean(patches, 1);
        patches = bsxfun(@minus, patches, patches_mean);
    end

    if do_normalize_std ==1
        % make it unit-variance
        patches_std = std(patches, [], 1);
        patches = bsxfun(@rdivide, patches, patches_std);
    end
end

anneal_counter = 0;
actual_lrate0 = actual_lrate;

if S.debug.do_display == 1
    figure(S.debug.display_fid);
end

try
    use_gpu = gpuDeviceCount;
catch errgpu
    use_gpu = false;
    disp(['Could not use CUDA. Error: ' errgpu.identifier])
end

for step=1:n_epochs
    if S.verbose
        fprintf(2, 'Epoch %d/%d: ', step, n_epochs)
    end
    if use_gpu
        % push
        for l = 1:n_layers
            if l < n_layers
                S.W{l} = gpuArray(single(S.W{l}));
            end
            S.biases{l} = gpuArray(single(S.biases{l}));
        end

        if S.adagrad.use
            for l = 1:n_layers
                if l < n_layers
                    S.adagrad.W{l} = gpuArray(single(S.adagrad.W{l}));
                end
                S.adagrad.biases{l} = gpuArray(single(S.adagrad.biases{l}));
            end
        elseif S.adadelta.use
            for l = 1:n_layers
                if l < n_layers
                    S.adadelta.gW{l} = gpuArray(single(S.adadelta.gW{l}));
                    S.adadelta.W{l} = gpuArray(single(S.adadelta.W{l}));
                end
                S.adadelta.gbiases{l} = gpuArray(single(S.adadelta.gbiases{l}));
                S.adadelta.biases{l} = gpuArray(single(S.adadelta.biases{l}));
            end
        end
    end

    for mb=1:n_minibatches
        S.iteration.n_updates = S.iteration.n_updates + 1;

        % p_0
        v0 = patches((mb-1) * minibatch_sz + 1:min(mb * minibatch_sz, n_samples), :);
        mb_sz = size(v0,1);

        if use_gpu > 0
            v0 = gpuArray(single(v0));
        end

        % add error
        v0_clean = v0;

        if S.data.binary == 0 && S.noise.level > 0
            if use_gpu
                v0 = v0 + S.noise.level * gpuArray(randn(size(v0)));
            else
                v0 = v0 + S.noise.level * randn(size(v0));
            end
        end

        if S.noise.drop > 0
            mask = binornd(1, 1 - S.noise.drop, size(v0));
            v0 = v0 .* mask;
            clear mask;
        end

        h0e = cell(n_layers, 1);
        h0e{1} = v0;

        for l = 2:n_layers
            h0e{l} = bsxfun(@plus, h0e{l-1} * S.W{l-1}, S.biases{l}');

            if l < n_layers || S.bottleneck.binary
                h0e{l} = sigmoid(h0e{l}, S.hidden.use_tanh);
            end
        end

        h0d = cell(n_layers, 1);
        h0d{end} = h0e{end};

        for l = n_layers-1:-1:1
            h0d{l} = bsxfun(@plus, h0d{l+1} * S.W{l}', S.biases{l}');
            if l == 1 && S.data.binary
                h0d{l} = sigmoid(h0d{l});
            end
            if l > 1
                h0d{l} = sigmoid(h0d{l}, S.hidden.use_tanh);
            end
        end

        % compute reconstruction error
        hr = sdae_get_hidden(v0_clean, S);
        vr = sdae_get_visible(hr, S);

        if S.data.binary
            rerr = -mean(sum(v0_clean .* log(max(vr, 1e-16)) + (1 - v0_clean) .* log(max(1 - vr, 1e-16)), 2));
        else
            rerr = mean(sum((v0_clean - vr).^2,2));
        end
        if use_gpu > 0
            rerr = gather(rerr);
        end
        S.signals.recon_errors = [S.signals.recon_errors rerr];

        % reset gradients
        for l = 1:n_layers
            biases_grad{l} = 0 * biases_grad{l};
            if l < n_layers
                W_grad{l} = 0 * W_grad{l};
            end
        end

        % backprop
        deltad = cell(n_layers, 1);
        deltad{1} = h0d{1} - v0_clean;
        biases_grad{1} = mean(deltad{1}, 1);

        for l = 2:n_layers
            deltad{l} = deltad{l-1} * S.W{l-1};
            if l < n_layers || S.bottleneck.binary
                deltad{l} = deltad{l} .* dsigmoid(h0d{l}, S.hidden.use_tanh);
            end
            biases_grad{l} = mean(deltad{l}, 1);
            W_grad{l-1} = (deltad{l-1}' * h0d{l}) / (size(v0, 1));
        end

        deltae = cell(n_layers, 1);
        deltae{end} = deltad{end};

        for l = n_layers-1:-1:1
            deltae{l} = deltae{l+1} * S.W{l}';
            if l == 1 && S.data.binary
                deltae{l} = deltae{l} .* dsigmoid(h0e{l});
            end
            if l > 1
                deltae{l} = deltae{l} .* dsigmoid(h0e{l}, S.hidden.use_tanh);
                biases_grad{l} = biases_grad{l} + mean(deltae{l}, 1);
            end
            W_grad{l} = W_grad{l} + (h0e{l}' * deltae{l+1}) / (size(v0, 1));
        end

        % learning rate
        if S.adagrad.use
            % update
            for l = 1:n_layers
                biases_grad_old{l} = (1 - momentum) * biases_grad{l} + momentum * biases_grad_old{l};
                if l < n_layers
                    W_grad_old{l} = (1 - momentum) * W_grad{l} + momentum * W_grad_old{l};
                end
            end

            for l = 1:n_layers
                if l < n_layers
                    S.adagrad.W{l} = S.adagrad.W{l} + W_grad_old{l}.^2;
                end

                S.adagrad.biases{l} = S.adagrad.biases{l} + biases_grad_old{l}.^2';
            end

            for l = 1:n_layers
                S.biases{l} = S.biases{l} - S.learning.lrate * (biases_grad_old{l}' + ...
                    weight_decay * S.biases{l}) ./ sqrt(S.adagrad.biases{l} + S.adagrad.epsilon);
                if l < n_layers
                    S.W{l} = S.W{l} - S.learning.lrate * (W_grad_old{l} + ...
                        weight_decay * S.W{l}) ./ sqrt(S.adagrad.W{l} + S.adagrad.epsilon);
                end
            end

        elseif S.adadelta.use
            % update
            for l = 1:n_layers
                biases_grad_old{l} = (1 - momentum) * biases_grad{l} + momentum * biases_grad_old{l};
                if l < n_layers
                    W_grad_old{l} = (1 - momentum) * W_grad{l} + momentum * W_grad_old{l};
                end
            end

            if S.iteration.n_updates == 1
                adamom = 0;
            else
                adamom = S.adadelta.momentum;
            end

            for l = 1:n_layers
                if l < n_layers
                    S.adadelta.gW{l} = adamom * S.adadelta.gW{l} + (1 - adamom) * W_grad_old{l}.^2;
                end

                S.adadelta.gbiases{l} = adamom * S.adadelta.gbiases{l} + (1 - adamom) * biases_grad_old{l}.^2';
            end

            for l = 1:n_layers
                dbias = -(biases_grad_old{l}' + ...
                    weight_decay * S.biases{l}) .* (sqrt(S.adadelta.biases{l} + S.adadelta.epsilon) ./ ...
                    sqrt(S.adadelta.gbiases{l} + S.adadelta.epsilon));
                S.biases{l} = S.biases{l} + dbias;

                S.adadelta.biases{l} = adamom * S.adadelta.biases{l} + (1 - adamom) * dbias.^2;
                clear dbias;

                if l < n_layers
                    dW = -(W_grad_old{l} + ...
                        weight_decay * S.W{l}) .* (sqrt(S.adadelta.W{l} + S.adadelta.epsilon) ./ ...
                        sqrt(S.adadelta.gW{l} + S.adadelta.epsilon));
                    S.W{l} = S.W{l} + dW;

                    S.adadelta.W{l} = adamom * S.adadelta.W{l} + (1 - adamom) * dW.^2;

                    clear dW;
                end

            end
        else
            if S.learning.lrate_anneal > 0 && (step >= S.learning.lrate_anneal * n_epochs)
                anneal_counter = anneal_counter + 1;
                actual_lrate = actual_lrate0 / anneal_counter;
            else
                if S.learning.lrate0 > 0
                    actual_lrate = S.learning.lrate / (1 + S.iteration.n_updates / S.learning.lrate0);
                else
                    actual_lrate = S.learning.lrate;
                end
                actual_lrate0 = actual_lrate;
            end

            S.signals.lrates = [S.signals.lrates actual_lrate];

            % update
            for l = 1:n_layers
                biases_grad_old{l} = (1 - momentum) * biases_grad{l} + momentum * biases_grad_old{l};
                if l < n_layers
                    W_grad_old{l} = (1 - momentum) * W_grad{l} + momentum * W_grad_old{l};
                end
            end

            for l = 1:n_layers
                S.biases{l} = S.biases{l} - actual_lrate * (biases_grad_old{l}' + weight_decay * S.biases{l});
                if l < n_layers
                    S.W{l} = S.W{l} - actual_lrate * (W_grad_old{l} + weight_decay * S.W{l});
                end
            end
        end

        if S.verbose == 1
            fprintf(2, '.');
        end

        if use_gpu > 0
            clear v0 h0d h0e v0_clean vr hr deltae deltad
        end

        if early_stop
            n_valid = size(valid_patches, 1);
            rndidx = randperm(n_valid);
            if use_gpu > 0
                v0valid = gpuArray(single(valid_patches(rndidx(1:round(n_valid * valid_portion)),:)));
            else
                v0valid = single(valid_patches(rndidx(1:round(n_valid * valid_portion)),:));
            end

            hr = sdae_get_hidden(v0valid, S);
            vr = sdae_get_visible(hr, S);

            if S.data.binary
                rerr = -mean(sum(v0valid .* log(max(vr, 1e-16)) + (1 - v0valid) .* log(max(1 - vr, 1e-16)), 2));
            else
                rerr = mean(sum((v0valid - vr).^2,2));
            end
            if use_gpu > 0
                rerr = gather(rerr);
            end

            S.signals.valid_errors = [S.signals.valid_errors rerr];

            if valid_err == -Inf
                valid_err = rerr;
                valid_best_err = rerr;
            else
                prev_err = valid_err;
                valid_err = 0.99 * valid_err + 0.01 * rerr;

                if step > S.valid_min_epochs && (1.1 * valid_best_err) < valid_err
                    fprintf(2, 'Early-stop! %f, %f\n', valid_err, prev_err);
                    stopping = 1;
                    break;
                end

                if valid_err < valid_best_err
                    valid_best_err = valid_err;
                end
            end
        else
            if S.stop.criterion > 0
                if S.stop.criterion == 1
                    if min_recon_error > S.signals.recon_errors(end)
                        min_recon_error = S.signals.recon_errors(end);
                        min_recon_error_update_idx = S.iteration.n_updates;
                    else
                        if S.iteration.n_updates > min_recon_error_update_idx + S.stop.recon_error.tolerate_count
                            fprintf(2, '\nStopping criterion reached (recon error) %f > %f\n', ...
                                S.signals.recon_errors(end), min_recon_error);
                            stopping = 1;
                            break;
                        end
                    end
                else
                    error ('Unknown stopping criterion %d', S.stop.criterion);
                end
            end
        end

        if length(S.hook.per_update) > 1
            err = S.hook.per_update{1}(S, S.hook.per_update{2});

            if err == -1
                stopping = 1;
                break;
            end
        end

        if S.debug.do_display == 1 && mod(S.iteration.n_updates, S.debug.display_interval) == 0
            S.debug.display_function (S.debug.display_fid, S, v0, v1, h0, h1, W_grad, vbias_grad, hbias_grad);
            drawnow;
        end
    end

    if use_gpu > 0
        % pull
        for l = 1:n_layers
            if l < n_layers
                S.W{l} = gather(S.W{l});
            end
            S.biases{l} = gather(S.biases{l});
        end

        if S.adagrad.use
            for l = 1:n_layers
                if l < n_layers
                    S.adagrad.W{l} = gather(S.adagrad.W{l});
                end
                S.adagrad.biases{l} = gather(S.adagrad.biases{l});
            end
        elseif S.adadelta.use
            for l = 1:n_layers
                if l < n_layers
                    S.adadelta.W{l} = gather(S.adadelta.W{l});
                    S.adadelta.gW{l} = gather(S.adadelta.gW{l});
                end
                S.adadelta.biases{l} = gather(S.adadelta.biases{l});
                S.adadelta.gbiases{l} = gather(S.adadelta.gbiases{l});
            end
        end
    end

    if length(S.hook.per_epoch) > 1
        err = S.hook.per_epoch{1}(S, S.hook.per_epoch{2});

        if err == -1
            stopping = 1;
        end
    end

    if stopping == 1
        break;
    end

    if S.verbose == 1
        fprintf(2, '\n');
    end

    fprintf(2, 'Epoch %d/%d - recon_error: %f\n', step, n_epochs, ...
        S.signals.recon_errors(end));
end

if use_gpu > 0
    % pull
    for l = 1:n_layers
        if l < n_layers
            S.W{l} = gather(S.W{l});
        end
        S.biases{l} = gather(S.biases{l});
    end

    if S.adagrad.use
        for l = 1:n_layers
            if l < n_layers
                S.adagrad.W{l} = gather(S.adagrad.W{l});
            end
            S.adagrad.biases{l} = gather(S.adagrad.biases{l});
        end
    elseif S.adadelta.use
        for l = 1:n_layers
            if l < n_layers
                S.adadelta.W{l} = gather(S.adadelta.W{l});
                S.adadelta.gW{l} = gather(S.adadelta.gW{l});
            end
            S.adadelta.biases{l} = gather(S.adadelta.biases{l});
            S.adadelta.gbiases{l} = gather(S.adadelta.gbiases{l});
        end
    end
end

author: kyunghyuncho

repo: deepmat


text:
function convnet=getConvNetPara(net)
convLayers=[];
for l=1:numel(net.layers)
    if(strcmp(net.layers{l}.type,'conv')||strcmp(net.layers{l}.type,'conv_aNet')||strcmp(net.layers{l}.type,'conv_mask'))
        convLayers(end+1)=l;
    end
end
len=length(convLayers);
convnet.targetLayers=convLayers+1;
convnet.targetScale=zeros(1,len);
convnet.targetStride=zeros(1,len);
convnet.targetCenter=zeros(1,len);
for i=1:len
    tarLay=convLayers(i);
    layer=net.layers{tarLay};
    pad=layer.pad(1);
    scale=size(layer.weights{1},1);
    stride=layer.stride(1);
    if(i==1)
        convnet.targetStride(i)=stride;
        convnet.targetScale(i)=scale;
        convnet.targetCenter(i)=(1+scale-pad*2)/2;
    else
        IsPool=false;
        poolStride=0;
        poolSize=0;
        poolPad=0;
        for j=convLayers(i-1)+1:tarLay-1
            if(strcmp(net.layers{j}.type,'pool'))
                IsPool=true;
                poolSize=net.layers{j}.pool(1);
                poolStride=net.layers{j}.stride(1);
                poolPad=net.layers{j}.pad(1);
            end
        end
        convnet.targetStride(i)=(1+IsPool*(poolStride-1))*stride*convnet.targetStride(i-1);
        convnet.targetScale(i)=convnet.targetScale(i-1)+IsPool*(poolSize-1)*convnet.targetStride(i-1)+convnet.targetStride(i)*(scale-1);
        if(IsPool)
            convnet.targetCenter(i)=(scale-pad*2-1)*poolStride*convnet.targetStride(i-1)/2+(convnet.targetCenter(i-1)+convnet.targetStride(i-1)*(poolSize-2*poolPad-1)/2);
        else
            convnet.targetCenter(i)=(scale-pad*2-1)*convnet.targetStride(i-1)/2+convnet.targetCenter(i-1);
        end
    end
end
end

repo: interpretableCNN


text:
function x = r83s_cg ( n, a, b, x )

%*****************************************************************************80
%
%% R83S_CG uses the conjugate gradient method on an R83S system.
%
%  Discussion:
%
%    The R83S storage format is used for a tridiagonal scalar matrix.
%    The vector A(3) contains the subdiagonal, diagonal, and superdiagonal
%    values that occur on every row.
%
%    The matrix A must be a positive definite symmetric band matrix.
%
%    The method is designed to reach the solution after N computational
%    steps.  However, roundoff may introduce unacceptably large errors for
%    some problems.  In such a case, calling the routine again, using
%    the computed solution as the new starting estimate, should improve
%    the results.
%
%  Example:
%
%    Here is how an R83S matrix of order 5, stored as (A1,A2,A3), would
%    be interpreted:
%
%      A2  A3   0   0   0
%      A1  A2  A3   0   0
%       0  A1  A2  A3   0
%       0   0  A1  A2  A3
%       0   0   0  A1  A2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Frank Beckman,
%    The Solution of Linear Equations by the Conjugate Gradient Method,
%    in Mathematical Methods for Digital Computers,
%    edited by John Ralston, Herbert Wilf,
%    Wiley, 1967,
%    ISBN: 0471706892,
%    LC: QA76.5.R3.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, real A(3), the matrix.
%
%    Input, real B(N), the right hand side vector.
%
%    Input/output, real X(N).
%    On input, an estimate for the solution, which may be 0.
%    On output, the approximate solution vector.
%
  b = b(:);
  x = x(:);
%
%  Initialize
%    AP = A * x,
%    R  = b - A * x,
%    P  = b - A * x.
%
  ap = r83s_mv ( n, n, a, x );

  r(1:n,1) = b(1:n,1) - ap(1:n,1);
  p(1:n,1) = b(1:n,1) - ap(1:n,1);
%
%  Do the N steps of the conjugate gradient method.
%
  for it = 1 : n
%
%  Compute the matrix*vector product AP=A*P.
%
    ap = r83s_mv ( n, n, a, p );
%
%  Compute the dot products
%    PAP = P*AP,
%    PR  = P*R
%  Set
%    ALPHA = PR / PAP.
%
    pap = p' * ap;
    pr = p' * r;

    if ( pap == 0.0 )
      return
    end

    alpha = pr / pap;
%
%  Set
%    X = X + ALPHA * P
%    R = R - ALPHA * AP.
%
    x(1:n,1) = x(1:n,1) + alpha * p(1:n,1);
    r(1:n,1) = r(1:n,1) - alpha * ap(1:n,1);
%
%  Compute the vector dot product
%    RAP = R*AP
%  Set
%    BETA = - RAP / PAP.
%
    rap = r' * ap;

    beta = - rap / pap;
%
%  Update the perturbation vector
%    P = R + BETA * P.
%
    p(1:n,1) = r(1:n,1) + beta * p(1:n,1);

  end

  return
end


text:
function [gKern, gVarmeans, gVarcovars] = kernVardistPsi0Gradient(kern, vardist, covGrad)

% KERNVARDISTPSI0GRADIENT description.
% VARGPLVM

    if ~strcmp(kern.type,'cmpnd')
       %
       fhandle = str2func([kern.type 'VardistPsi0Gradient']);
       [gKern, gVarmeans, gVarcovars] = fhandle(kern, vardist, covGrad);

       % Transformations
       gKern = paramTransformPsi0(kern, gKern);
       %
    else % the kernel is cmpnd
       %
       fhandle = str2func([kern.comp{1}.type 'VardistPsi0Gradient']);
       [gKern, gVarmeans, gVarcovars] = fhandle(kern.comp{1}, vardist, covGrad);
       % Transformations
       gKern = paramTransformPsi0(kern.comp{1}, gKern);
       %
       for i = 2:length(kern.comp)
           %
           fhandle = str2func([kern.comp{i}.type 'VardistPsi0Gradient']);
           [gKerni, gVarmeansi, gVarcovarsi] = fhandle(kern.comp{i}, vardist, covGrad);

           % Transformations
           gKerni = paramTransformPsi0(kern.comp{i}, gKerni);

           gVarmeans = gVarmeans + gVarmeansi;
           gVarcovars = gVarcovars + gVarcovarsi;
           gKern = [gKern gKerni];
           %
       end
       %
    end
end

% variational variances are positive (This should rather go to
% vargplvmLogLikeGradients)
%gVarcovars = (gVarcovars(:).*vardist.covars(:))';

%-----------------------------------------------------
% This applies transformations
% This must be done similar to kernGradient at some point
% function gK = paramTransformPsi0(kern, gK)
% %
% %
%     if strcmp(kern.type,'rbfardjit') | strcmp(kern.type,'rbfard2') | strcmp(kern.type,'bias') | strcmp(kern.type,'white')
%         gK(1) = gK(1)*kern.variance;
%     elseif strcmp(kern.type,'linard2')
%         gK(1:end) = gK(1:end).*kern.inputScales;
%     end
% end

function gKern = paramTransformPsi0(kern, gKern)

    fhandle = str2func([kern.type 'KernExtractParam']);
    params  = fhandle(kern);

    if length(kern.transforms) > 1
        for i=1:length(kern.transforms)
            fhandle = str2func([kern.transforms(i).type 'Transform']);
            if ~isfield(kern.transforms(i), 'transformsettings')
                gKern(kern.transforms(i).index) = gKern(kern.transforms(i).index).*fhandle(params(kern.transforms(i).index), 'gradfact');
            else
                gKern(kern.transforms(i).index) = gKern(kern.transforms(i).index).*fhandle(params(kern.transforms(i).index), 'gradfact', kern.transforms(i).transformsettings);
            end
        end
    else
        fhandle = str2func([kern.transforms.type 'Transform']);
        if ~isfield(kern.transforms, 'transformsettings')
            gKern(kern.transforms.index) = gKern(kern.transforms.index).*fhandle(params(kern.transforms.index), 'gradfact');
        else
            gKern(kern.transforms.index) = gKern(kern.transforms.index).*fhandle(params(kern.transforms.index), 'gradfact', kern.transforms.transformsettings);
        end
    end

end

author: SheffieldML

repo: vargplvm


text:
function [A] = au2A(au)
% Convert length from astronomical units to angstroms.
% Chad A. Greene 2012
A = au*1.49597870691e+21;


text:
function [Prefix, Ext] = Remove1DExtension (fname)
%
%   [Prefix, Ext] = Remove1DExtension (fname)
%
%Purpose:
%   removes known .1D filename extensions.
%   Look at code more more details.
%
%Input Parameters:
%   fname: filename
%
%
%Output Parameters:
%   Prefix: filename without extension
%   Ext: recognized extension of 1D filename, if any
%        if Ext is empty then no know extensions were found
%
%
%Key Terms:
%
%More Info :
%   RemoveExtension
%
%
%
%     Author : Ziad Saad
%     Date : Tue Nov 9 09:59:58 EST 2004
%     SSCC/NIMH/ National Institutes of Health, Bethesda Maryland

%Define the function name for easy referencing
FuncName = 'Remove1DExtension';

%Debug Flag
DBG = 1;

Prefix = ''; Ext = '';

if (isempty(fname)) return; end
[Prefix, Ext] = RemoveExtension(fname,...
            '.1D|.1D.dset|.1D.col|.1D.NodeList|.1D.FaceSetList|.1D.cmap');

return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ r,flag ] =  reward_grid_world_continuous(s,x,beta,s_goal,x_goal,x_bad,r_bad,r_goal)
%
%
%
if s == s_goal

    r    = r_goal;
    flag = true;

elseif ismember(s,x_bad)

    r    = r_bad;
    flag = true;

else
    r    = r_goal .* exp(-beta * pdist([x;x_goal]).^2);
    flag = false;

end

end

author: epfl-lasa


text:
%%
%An example to RMA from the book: Carrara:"Spotlight Synthetic Aperture
%Radar: Signal Processing Algortithms" Artech House 1995.
%Signal Data from: Tab10_1
%%
clear all,close all; clc
Nf=2^10; %number of points for freq. range

%array parameter
dx=0.6;L=760.8; %interelement space in synthetic array (length L)
Na=L/dx+1; %number of array elements
X=linspace(-L/2,L/2,Na).'; %position of array elements
Kx=linspace(-pi/dx,pi/dx,Na).';%azimuthal spacial wavenumber (page 410)
Kx=repmat(Kx,1,Nf);
Rs=1000; %distance from origin to the array plane
X=X+1j*Rs; % XY antenna position as complex number. Here, image processing
%in XY plane only, therefore, complex description for cartesian coordinates
%is convinient.

%frequency range
co=299792458; %free space wave velocity
fc=242.4*1E6; %last frequency in the Band of width B (in the book it should
%be center freq., but the Range Freq. in Fig. 10.11 doesn't fit then.)
B=133.5*1E6; %bandwidth
freq=linspace(fc-B,fc,Nf);
Kr=4*pi*freq/co;%Frequency wavenumber (page 410)
Kr=repmat(Kr,Na,1);

%target position:
tar=[-200+1j*100, 200+1j*100, 0+1j*100,...
    -200, 0, 200,...
    -200-1j*100, 200-1j*100, 0-1j*100]; %3x3 grid
%tar=[0, 200, -1j*100]; %(table 10.1)

%simulated data, dechirped (compressed) and deskewed in range:
data=my_sinthSISO_planresp( X,freq,tar );
%phase reference acc. to image origin
%(see also: Sanchez "3-D Radar Imaging Using Range Migration Techniques",
%IEEE; or:
%http://www.lluisvives.com/servlet/SirveObras/01604307547816095212257/003402_12.pdf):
data=data.*exp(1j*Rs*Kr);
%
%fig10.7 and 10.8
data_compr=fftshift(ifft(data,[],2),2);
imagesc(abs(data_compr));
%fig10.7

%azimuth fft
data=fftshift(data,1); %due to -Xmin..Xmax
data=fft(data,[],1);%azimuth fft
data=fftshift(data,1); %due to -Kx_min..Kx_max
%azimuth fft

%Kx accord to fig10.11
plot(Kx,Kr,'.')
Xax=-5:1:5;Yax=5:1:10;
axis([ min(Xax) max(Xax) min(Yax) max(Yax)])
set(gca,'XTick',Xax,'YTick',Yax)
grid on
%Kx accord to fig10.11

%fig10.11
imagesc(Kr(1,:),Kx(:,1),abs(data))
%fig10.11

%range compressed signal in Fig10.12
data_compr=fftshift(ifft(data,[],2),2);
%fig10.12
imagesc(abs(data_compr));
%fig10.12

%Ky wavenumber, squared, accordingly to (10.9)
Ky_sq=Kr.^2-Kx.^2;
%only positive Ky is physical
Neg=logical(Ky_sq < 0); %
Ky_sq(Neg)=NaN; data(Neg)=NaN; Kr(Neg)=NaN; %only positive Ky is acceptable
%matched filter, phase acc. to (10.8)
Fmf=-Rs*Kr+Rs*sqrt(Ky_sq);
%figure 10.13
mesh(Kr,Kx,Fmf)
axis([4 10 -2.8 2.8 -2000 0])
%figure 10.13

data=data.*exp(1j*Fmf); %data matching by matched filter
%fig. 10.14 middle
data_compr=fftshift(ifft(data,[],2),2);
imagesc(abs(data_compr));
%fig. 10.14 middle

%stolt interpolation
%Fig 10.15 (a)
plot(Kx,Kr,'.')
%Fig 10.15 (a)
Ky=sqrt(Ky_sq); %cross-range wavenumber by (10.9)
%Fig 10.15 (b)
plot(Kx,Ky,'.')
%Fig 10.15 (b)
Ky_int=Ky(Na/2-0.5,:);%data interpolation to Ky_int cross-range wavenumber
data=my_spline(Ky,data,Ky_int);%1D data interpolation by spline function
%stolt interpolation

%processing window in Ky dimension
Ny=750;
data=data(:,1:Ny);
Ky_int=Ky_int(1:Ny);

%range compression
data=fftshift(ifft(data,[],2),2);
%fig10.20 (a)
imagesc(abs(data));
%fig10.20 (a)

%azimuth compression
data=fftshift(data,1);%compare to azimuth fft above
data=ifft(data,[],1);
data=fftshift(data,1);
%%
dy=2*pi/(Ky_int(end)-Ky_int(1)); %step in cross-range dimension
Y=-Ny/2*dy:dy:(Ny/2-1)*dy;%cross-range axis
%fig10.20 (b)
imagesc(Y,X,abs(data));
ylabel('Cross-Range, m')
xlabel('Range, m')
axis([ -300 300 -300 300])
colormap('pink')
%fig10.20 (b)


text:
  function [net, info] = cnn_single_rgb(varargin)
%CNN_SINGLE_RGB Demonstrates fine-tuning a pre-trained CNN with static
% RGB frames (SI in pami journal) on UCF101 dataset

run(fullfile(fileparts(mfilename('fullpath')), ...
  '..', 'matconvnet', 'matlab', 'vl_setupnn.m')) ;

addpath Layers Datasets

opts.dataDir = fullfile('data','UCF101') ;
opts.expDir  = fullfile('exp', 'UCF101') ;
opts.modelPath = fullfile('models','resnext_50_32x4d-pt-mcn.mat');
opts.datasetFn = @cnn_ucf101_setup_data ;
opts.networkFn = @cnn_init_resnext ;
opts.pool1Type = 'none' ;
opts.pool1Layer = 'conv1' ;
opts.pool2Layer = '' ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.numFetchThreads = 8 ;

opts.lite = false ;
opts.imdbPath = fullfile(opts.dataDir, 'imdb-rgb.mat');
opts.ARPoolLayer = 'conv0'; % before conv1
opts.DropOutRate = 0.5 ;
opts.epochFactor = 5 ;

opts.split = 1; % data split
opts.reverseDyn = 0; % reverse video frames e.g.[N:-1:1]
opts.train = struct() ;
opts.train.gpus = [];
opts.train.batchSize = 128 ;
opts.train.numSubBatches = 16 ;
opts.train.solver = [] ;
opts.train.prefetch = true ;
opts.train.numEpochs = 30 ;
% resnet50
opts.train.learningRate = 1e-2 * [ones(1,2), 0.1*ones(1,1), 0.01*ones(1,1)];
% caffe-ref
opts.train.learningRate = 1e-4 * [ones(1,2), 0.1*ones(1,1), 0.01*ones(1,1)];

opts = vl_argparse(opts, varargin) ;
if ~isfield(opts.train, 'gpus'), opts.train.gpus = []; end;
% opts.train.numEpochs = numel(opts.train.learningRate);

% -------------------------------------------------------------------------
%                                                              Prepare data
% -------------------------------------------------------------------------

if exist(opts.imdbPath,'file')
  imdb = load(opts.imdbPath) ;
else
  imdb = opts.datasetFn('dataDir', opts.dataDir, 'lite', opts.lite) ;
  mkdir(opts.expDir) ;
  save(opts.imdbPath, '-struct', 'imdb') ;
end

% UCF101 has 3 data splits
if opts.split>3
  error('split should be <=3');
end
imdb.images.set = imdb.images.sets(opts.split,:);

% reverse frame order
if opts.reverseDyn
  for i=1:numel(imdb.images.names)
    imdb.images.names{i} = imdb.images.names{i}(end:-1:1);
  end
end

% -------------------------------------------------------------------------
%                                                             Prepare model
% -------------------------------------------------------------------------
net = load(opts.modelPath);
if isfield(net,'net')
  net = net.net;
end
opts.nCls = max(imdb.images.label) ;
net = opts.networkFn(net,opts);

if numel(net.meta.normalization.averageImage)>3
  sz = size(net.meta.normalization.averageImage) ;
  net.meta.normalization.averageImage = ...
    mean(reshape(net.meta.normalization.averageImage,[sz(1)*sz(2) sz(3)]),1) ;
end

% Set the class names in the network
net.meta.classes.name = imdb.classes.name ;
net.meta.classes.description = imdb.classes.name ;
% -------------------------------------------------------------------------
%                                                                     Learn
% -------------------------------------------------------------------------
if opts.epochFactor>0
  opts.train.train = repmat(find(imdb.images.set==1),[1 opts.epochFactor]) ;
else
  opts.train.train = NaN ;
end
opts.train.val = find(imdb.images.set==3) ;

[net, info] = cnn_train_dicnn_dag(net, imdb, getBatchFn(opts, net.meta), ...
                      'expDir', opts.expDir, ...
                      opts.train) ;

% -------------------------------------------------------------------------
%                                                          Report accuracy
% -------------------------------------------------------------------------
errlayer = net.getLayerIndex('errMC') ;

if ~isnan(errlayer)
  cats = imdb.classes.name ;
  accs = net.layers(errlayer).block.accuracy ;

  if numel(cats)~=numel(accs)
    error('wrong number of classes\n') ;
  end

  for i=1:numel(cats)
    fprintf('%s acc %.1f\n',cats{i},100*accs(i)) ;
  end
  fprintf('Mean accuracy %.1f\n',100*mean(accs)) ;
end

% -------------------------------------------------------------------------
function fn = getBatchFn(opts, meta)
% -------------------------------------------------------------------------
useGpu = numel(opts.train.gpus) > 0 ;

bopts.numThreads = opts.numFetchThreads ;
bopts.imageSize = meta.normalization.imageSize ;
if isfield(meta.normalization,'border')
  bopts.border = meta.normalization.border ;
else
  bopts.border = meta.normalization.imageSize(1:2) ./ ...
    meta.normalization.cropSize - meta.normalization.imageSize(1:2);

end

% bopts.averageImage = [];
bopts.averageImage = meta.normalization.averageImage ;
bopts.interpolation = meta.normalization.interpolation ;
bopts.keepAspect = meta.normalization.keepAspect ;
% bopts.rgbVariance = meta.augmentation.rgbVariance ;
% bopts.transformation = meta.augmentation.transformation ;

fn = @(x,y) getDagNNBatch(bopts,useGpu,x,y) ;

% -------------------------------------------------------------------------
function inputs = getDagNNBatch(opts, useGpu, imdb, batch)
% -------------------------------------------------------------------------

% batch refers to videos (not for frames)
if isempty(batch)
  inputs = {'input', [], 'label', [], 'VideoId1', [], 'VideoId2', []};
  return;
end

isVal = ~isempty(batch) && imdb.images.set(batch(1)) ~= 1 ;

% if ~isVal, transformation='stretch'; else transformation='none';end
if ~isVal, transformation='multiScaleRegular'; else transformation='none';end

names = imdb.images.names(batch);

% images = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;

namesM = {};
nVids = numel(batch);

VideoId1 = [];
VideoId2 = [];

% step-size
stepSize = 6;
% pool nFrames into a dynamic image
nFrames = 1;
% number of dynamic images to be max pooled later
nDynImgs = 10;

c1 = 1;
for v=1:nVids

  name = names{v};
  nFrms = numel(name);

  nSample = nFrames;
  nr = numel(1:stepSize:nFrms);

  % jitter by removing 50 % and limit a batch to nMaxs * nSamples images
  if nr > 1 && (~isVal && nr>nDynImgs)
    rat = min(nDynImgs,ceil(0.50*nr));
    ri = randperm(nr);
    ri = ri(1:rat);
    r = zeros(1,nr);
    r(ri) = 1;
  else
    r = ones(1,nr);
  end

  c3 = 1;
  c2 = 0;

  for f=1:stepSize:nFrms
    if r(c3)
      idx = f:min(f+nSample-1,nFrms) ;
      if numel(idx)<nFrames
        idx = [idx idx(end) * ones(1,nFrames-numel(idx))];
      end
      namesM{end+1} = name(idx);
      VideoId1 = [VideoId1 c1 * ones(1,numel(idx))];
      c1 = c1 + 1;
      c2 = c2 + 1;
    end
    c3 = c3 + 1;
  end
  VideoId2 = [VideoId2 v * ones(1,c2) ] ;
end

images = strcat([imdb.imageDir filesep], horzcat(namesM{:}) ) ;

im = cnn_video_rgb_get_batch(images, VideoId1, opts, ...
  'transformation', transformation, 'prefetch', nargout == 0, ...
  'subMean', false) ;

if nargout > 0
  if useGpu
    im = gpuArray(im) ;
  end
  inputs = {'input', im, 'label', imdb.images.label(batch), ...
    'VideoId2', VideoId2};
end

author: hbilen

repo: dynamic-image-nets


text:
function q = reshape(q,varargin)
% overloads reshape

q.a = reshape(q.a,varargin{:});
q.b = reshape(q.b,varargin{:});
q.c =  reshape(q.c,varargin{:});
q.d = reshape(q.d,varargin{:});
q.i = reshape(q.i,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
function padua_test02 ( )

%*****************************************************************************80
%
%% PADUA_TEST02 tests PADUA_POINTS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 June 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PADUA_TEST02\n' );
  fprintf ( 1, '  PADUA_POINTS returns the points of a Padua rule.\n' );

  for l = 0 : 10
    n = padua_order ( l );
    xy = padua_points ( l );
    label = sprintf ( '  Level %d Padua points:', l );
    r8mat_transpose_print ( 2, n, xy, label );
  end

  return
end


text:
function p = completeInitial(p)

if p.options.warmstart && any(isnan(p.x0))

    p_reduced = p;

    % Fix the given variables
    p_reduced.lb(~isnan(p.x0)) = p.x0(~isnan(p.x0));
    p_reduced.ub(~isnan(p.x0)) = p.x0(~isnan(p.x0));

    % Propagate these values a bit
    p_reduced = propagate_bounds_from_equalities(p_reduced);
    p_reduced = propagate_bounds_from_evaluations(p_reduced);
    p_reduced = propagate_bounds_from_equalities(p_reduced);

    % Keep only equalities and elementwise inerqualities (more?)
    p_reduced.F_struc = p.F_struc(1:p.K.f + p.K.l,:);
    p_reduced.K.q = 0;
    p_reduced.K.r = 0;
    p_reduced.K.p = 0;
    p_reduced.K.s = 0;

    % No objective
    p_reduced.c = p_reduced.c*0;
    p_reduced.Q = p_reduced.Q*0;

    % Kill nonlinearity information
    p_reduced.variabletype = p_reduced.variabletype*0;

    % Solve LP
    output = feval(p.solver.lpcall,p_reduced);
    if output.problem == 0
        p.x0 = output.Primal;
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function jed = yjf_to_jed_islamic_a ( y, j, f )

%*****************************************************************************80
%
%% YJF_TO_JED_ISLAMIC_A converts an Islamic-A YJF date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, J, real F, the YJF date.
%
%    Output, real JED, the Julian Ephemeris Date.
%

%
%  Copy the input.
%
  y1 = y;
  j1 = j;
  f1 = f;
%
%  Check the input.
%
  [ y1, j1, ierror ] = yj_check_islamic ( y1, j1 );

  if ( ierror ~= 0 )
    jed = -1.0;
    return
  end
%
%  Convert the input.
%
  [ y2, m2, d2, f2 ] = yjf_to_ymdf_islamic ( y1, j1, f1 );

  jed = ymdf_to_jed_islamic_a ( y2, m2, d2, f2 );

  return
end