text:
classdef Generic3DTrajectoryViewType < AbstractTrajectoryViewType
    %Inertial3DTrajectoryViewType Summary of this class goes here
    %   Detailed explanation goes here

    properties

    end

    methods
        function obj = Generic3DTrajectoryViewType()

        end

        function [hCBodySurf, childrenHGs] = plotStateLog(obj, orbitNumToPlot, lvdData, viewProfile, handles, app)
%             dAxes = app.dispAxes;
            dAxes = handles.dispAxes;
            hFig = app.ma_LvdMainGUI;
            celBodyData = lvdData.celBodyData;
            stateLog = lvdData.stateLog;

%             axes(dAxes);
%             cla(dAxes);
%             cla(dAxes,'reset');
            delete(dAxes.Children);
            dAxes.Color = viewProfile.backgroundColor.color;

            hFig.Renderer = viewProfile.renderer.renderer;
            if(viewProfile.renderer == FigureRendererEnum.OpenGL && ~isunix())
                d = opengl('data');
                if(strcmpi(d.HardwareSupportLevel,'full'))
                    opengl hardware;
                elseif(strcmpi(d.HardwareSupportLevel,'basic'))
                    opengl hardwarebasic;
                end
            end

            hFig.GraphicsSmoothing = 'on';

            if(stateLog.getNumberOfEntries() == 0)
                return;
            end

            viewInFrame = viewProfile.frame;
            viewCentralBody = viewInFrame.getOriginBody();
            lvdStateLogEntries = LaunchVehicleStateLogEntry.empty(1,0);
            switch viewProfile.trajEvtsViewType
                case ViewEventsTypeEnum.SoIChunk
                    entries = stateLog.getAllEntries();
%                     maStateLog = stateLog.getMAFormattedStateLogMatrix(false);

%                     chunkedStateLog = breakStateLogIntoSoIChunks(maStateLog);
                    chunkedStateLog = stateLog.breakUpStateLogBySoIChunk();
                    if(orbitNumToPlot > size(chunkedStateLog,1))
                        orbitNumToPlot = size(chunkedStateLog,1);
                        set(dAxes,'UserData',orbitNumToPlot);
                    end
                    subStateLogs = chunkedStateLog(orbitNumToPlot,:);

                    bodyId = subStateLogs{1}(1,8);
                    bodyInfo = celBodyData.getBodyInfoById(bodyId);
                    evtIds = [];
                    for(i=1:length(subStateLogs))
                        subStateLog = subStateLogs{i};

                        if(isempty(subStateLog))
                            continue;
                        end

                        bodyId = subStateLog(1,8);
                        bodyInfo = celBodyData.getBodyInfoById(bodyId);
                        inertialFrame = bodyInfo.getBodyCenteredInertialFrame();

                        elemSet = CartesianElementSet(subStateLog(:,1), subStateLog(:,2:4)', subStateLog(:,5:7)', inertialFrame);
%                         elemSet = repmat(CartesianElementSet.getDefaultElements(), [1, size(subStateLog,1)]);
%                         for(j=1:size(subStateLog,1))
%                             elemSet(j) = CartesianElementSet(subStateLog(j,1), subStateLog(j,2:4)', subStateLog(j,5:7)', inertialFrame);
%                         end
                        elemSet = convertToFrame(elemSet, viewInFrame);

                        subStateLog(:,2:4) = [elemSet.rVect]';
                        subStateLog(:,5:7) = [elemSet.vVect]';

                        evtIds = [evtIds, unique(subStateLog(:,13))']; %#ok<AGROW>

                        subStateLogs{i} = subStateLog;
                    end

                    evtIds = unique(evtIds);
                    for(i=1:length(entries))
                        entry = entries(i);

                        if(not(isempty(entry.event.getEventNum())) && ...
                           ismember(entry.event.getEventNum(), evtIds) && ...
                           entry.centralBody == bodyInfo)
                            lvdStateLogEntries(end+1) = entry; %#ok<AGROW>
                        end
                    end

                    numTotMissionSegs = size(chunkedStateLog,1);

%                     curMissionSegStr = num2str(orbitNumToPlot);
%                     totalMissionSegStr = num2str(numTotMissionSegs);

                    if(numTotMissionSegs <= 1)
                        app.decrOrbitToPlotNum.Enable = 'off';
                        app.incrOrbitToPlotNum.Enable = 'off';
                    else
                        app.decrOrbitToPlotNum.Enable = 'on';
                        app.incrOrbitToPlotNum.Enable = 'on';
                    end
                case ViewEventsTypeEnum.All
                    entries = stateLog.getAllEntries();
                    maStateLogMatrix = stateLog.getMAFormattedStateLogMatrix(false);
                    numRows = size(maStateLogMatrix,1);
%                     subStateLogsMat = NaN(numRows, 13);

                    cartesianEntry = convertToFrame(getCartesianElementSetRepresentation(entries, false),viewInFrame);
                    times = [cartesianEntry.time]';
                    rVect = [cartesianEntry.rVect]';
                    vVect = [cartesianEntry.vVect]';
                    bodyId = viewCentralBody.id + zeros(numRows,1);
                    subStateLogsMat = [times, rVect, vVect, bodyId, maStateLogMatrix(:,9:13)];
%                     for(i=1:numRows)
%                         tempMaMatrix = entries(i).getMAFormattedStateLogMatrix(false);
%
%                         cartesianEntry =
%                         cartesianEntry = cartesianEntry.convertToFrame(viewInFrame);
%
%                         subStateLogsMat(i,:) = [entries(i).time, cartesianEntry.rVect', cartesianEntry.vVect', viewCentralBody.id, maStateLogMatrix(i,9:13)];
%                     end

                    subStateLogs = {};
                    for(evtNum=1:max(subStateLogsMat(:,13)))
                        subStateLogs{evtNum} = subStateLogsMat(subStateLogsMat(:,13) == evtNum,:); %#ok<AGROW>
                    end

%                     curMissionSegStr = num2str(1);
%                     totalMissionSegStr = num2str(1);

                    lvdStateLogEntries = entries;
                otherwise
                    error('Unknown trajectory view type when plotting trajectory: %s', viewProfile.frame.name);
            end

            eventsList = [];
            minTime = Inf;
            maxTime = -Inf;
            for(i=1:length(subStateLogs))
                if(~isempty(subStateLogs{i}))
                    eventsList = [eventsList;unique(subStateLogs{i}(:,13))]; %#ok<AGROW>
                end
                if(i>1)
                    prevSubStateLog = subStateLogs{i-1};
                else
                    prevSubStateLog = NaN(1,size(subStateLogs{i},2));
                end

                if(size(subStateLogs{i},1)>1)
                    [childrenHGs] = plotSubStateLog(subStateLogs{i}, prevSubStateLog, lvdData, dAxes);

                    minTime = min([minTime, min(subStateLogs{i}(:,1))]);
                    maxTime = max([maxTime, max(subStateLogs{i}(:,1))]);
                end
            end

            if(viewInFrame.typeEnum == ReferenceFrameEnum.BodyFixedRotating && ...
               viewProfile.showLongLatAnnotations)
                plotBodyFixedGrid(dAxes, viewCentralBody);
            end

            showSoI = viewProfile.showSoIRadius;

            if(viewProfile.showThrustVectors)
                entryInc = viewProfile.thrustVectEntryIncr;
                scale = viewProfile.thrustVectScale;
                color = viewProfile.thrustVectColor.color;
                lineStyle = viewProfile.thrustVectLineType.linespec;

                subsetLvdStateLogEntries = lvdStateLogEntries(1:entryInc:length(lvdStateLogEntries));
                subsetLvdStateLogEntries = subsetLvdStateLogEntries(:)';
                cartesianEntries = convertToFrame(getCartesianElementSetRepresentation(subsetLvdStateLogEntries), viewInFrame);

                rVects = [cartesianEntries.rVect];
                tVects = [];

                [~, ~, ~, rotMatToInertial12] = viewInFrame.getOffsetsWrtInertialOrigin([cartesianEntries.time], cartesianEntries);
                for(i=1:length(subsetLvdStateLogEntries)) %#ok<*NO4LP>
                    entry = subsetLvdStateLogEntries(i);
                    cartesianEntry = cartesianEntries(i);

                    tVect = lvd_ThrottleTask(entry, 'thrust_vector', cartesianEntry.frame);

                    if(norm(tVect) > 0)
                        [~, ~, ~, rotMatToInertial32] = entry.centralBody.getBodyCenteredInertialFrame().getOffsetsWrtInertialOrigin(entry.time, cartesianEntry);

                        tVectNew = rotMatToInertial32 * rotMatToInertial12(:,:,i)' * tVect;
                    else
                        tVectNew = [0;0;0];
                    end

                    tVects = [tVects, tVectNew]; %#ok<AGROW>
                end

                tVects = scale .* tVects;

                hold(dAxes,'on');
                quiver3(dAxes, rVects(1,:),rVects(2,:),rVects(3,:), tVects(1,:),tVects(2,:),tVects(3,:), 0, 'Color',color, 'LineStyle',lineStyle);
                hold(dAxes,'off');
            end

%             lfm = LiftForceModel();
%             rVects = NaN([3, length(lvdStateLogEntries)]);
%             forceVect = NaN([3, length(lvdStateLogEntries)]);
%             for(i=1:length(lvdStateLogEntries))
%                 stateLogEntry = lvdStateLogEntries(i);
%
%                 bodyInfo = stateLogEntry.centralBody;
%                 ut = stateLogEntry.time;
%                 rVect = stateLogEntry.position;
%                 vVect = stateLogEntry.velocity;
%                 aero = stateLogEntry.aero;
%                 mass = stateLogEntry.getTotalVehicleMass();
%                 attState = stateLogEntry.attitude;
%
%                 rVects(:,i) = rVect;
%                 force = lfm.getForce(ut, rVect, vVect, mass, bodyInfo, aero, [], [], [], [], [], [], [], [], [], [], attState);
%                 forceVect(:,i) = 1000*force;
%             end
%
%             hold(dAxes,'on');
%             quiver3(dAxes, rVects(1,:),rVects(2,:),rVects(3,:), forceVect(1,:),forceVect(2,:),forceVect(3,:), 0, 'Color','c', 'LineStyle','-');
%             hold(dAxes,'off');

            if(showSoI && ~isempty(viewCentralBody.getParBodyInfo(celBodyData)))
                hold(dAxes,'on');

                r = getSOIRadius(viewCentralBody, viewCentralBody.getParBodyInfo(celBodyData));

                x = r*sin(0:0.01:2*pi);
                y = r*cos(0:0.01:2*pi);
                z = zeros(size(x));
                plot3(dAxes, x, y, z, 'k--','LineWidth',0.5);
                plot3(dAxes, y, z, x, 'k--','LineWidth',0.5);
                plot3(dAxes, z, x, y, 'k--','LineWidth',0.5);

                hold(dAxes,'off');
            end

            if(viewProfile.dispXAxis || viewProfile.dispYAxis || viewProfile.dispZAxis)
                axisLength = 2*viewCentralBody.radius;

                hold(dAxes,'on');
                if(viewProfile.dispXAxis)
                    quiver3(dAxes, 0,0,0, axisLength,0,0, 'r', 'LineWidth',2);
                end

                if(viewProfile.dispYAxis)
                    quiver3(dAxes, 0,0,0, 0,axisLength,0, 'g', 'LineWidth',2);
                end

                if(viewProfile.dispZAxis)
                    quiver3(dAxes, 0,0,0, 0,0,axisLength, 'b', 'LineWidth',2);
                end
                hold(dAxes,'off');
            end

            %plot central body
            [hCBodySurf, hCBodySurfXForm] = ma_initOrbPlot(hFig, dAxes, viewCentralBody);
            hCBodySurf.EdgeAlpha = viewProfile.meshEdgeAlpha;

            if(viewProfile.showAtmosphere && viewCentralBody.atmohgt > 0)
                hold(dAxes,'on');
                atmoRadius = viewCentralBody.radius + viewCentralBody.atmohgt;
                [X,Y,Z] = sphere(50);
                hCBodySurf = surf(dAxes, atmoRadius*X,atmoRadius*Y,atmoRadius*Z, 'BackFaceLighting','lit', 'FaceLighting','gouraud', 'FaceColor',[223 223 223]/255, 'FaceAlpha',0.2, 'EdgeLighting','gouraud', 'LineWidth',0.1, 'EdgeColor','none');
                hold(dAxes,'off');
            end

            eventsList = unique(eventsList);
            minEventNum = min(eventsList);
            maxEventNum = max(eventsList);

            if(minEventNum < maxEventNum)
                eventStr = ['Events ', num2str(minEventNum), ' - ', num2str(maxEventNum)];
            else
                eventStr = ['Event ', num2str(minEventNum)];
            end

            hDispAxisTitleLabel = handles.dispAxisTitleLabel;
            titleStr = sprintf('%s Orbit -- %s\n%s', viewCentralBody.name, eventStr, viewInFrame.getNameStr());
            hDispAxisTitleLabel.String = titleStr;
            hDispAxisTitleLabel.TooltipString = sprintf('Frame: %s', viewInFrame.getNameStr());

            set(dAxes,'LineWidth',1);
            set(dAxes,'Box','on');
            grid(dAxes,viewProfile.gridType.gridStr);
            dAxes.GridColor = viewProfile.majorGridColor.color;
            dAxes.MinorGridColor = viewProfile.minorGridColor.color;
            dAxes.GridAlpha = viewProfile.gridTransparency;
%             axis(dAxes,'equal');
%             axis(dAxes,'tight');

%             xlabel(dAxes, '');
%             ylabel(dAxes, '');

            set(dAxes,'XTickLabel',[]);
            set(dAxes,'YTickLabel',[]);
            set(dAxes,'ZTickLabel',[]);

            hold(dAxes,'off');

            if(not(viewProfile.updateViewAxesLimits))
%                 setappdata(handles.ma_LvdMainGUI,'dispOrbitXLim',xlim(dAxes));
%                 setappdata(handles.ma_LvdMainGUI,'dispOrbitYLim',ylim(dAxes));
%                 setappdata(handles.ma_LvdMainGUI,'dispOrbitZLim',zlim(dAxes));

%                 zoom reset;

%                 view(dAxes,viewProfile.viewAzEl);

                camPos = viewProfile.viewCameraPosition;
                camTgt = viewProfile.viewCameraTarget;
                camUpVec = viewProfile.viewCameraUpVector;
                camVA = viewProfile.viewCameraViewAngle;
                if(not(any(isnan(camPos))))
                    dAxes.CameraPosition = camPos;
                end

                if(not(any(isnan(camTgt))))
                    dAxes.CameraTarget = camTgt;
                end

                if(not(any(isnan(camUpVec))))
                    dAxes.CameraUpVector = camUpVec;
                end

                if(not(any(isnan(camVA))))
                    dAxes.CameraViewAngle = camVA;
                end

%                 if(any(isnan(viewProfile.viewZoomAxLims)))
%                     viewProfile.viewZoomAxLims = [xlim(dAxes);
%                                                   ylim(dAxes);
%                                                   zlim(dAxes)];
%                 else
%                     dAxes.XLim = viewProfile.viewZoomAxLims(1,:);
%                     dAxes.YLim = viewProfile.viewZoomAxLims(2,:);
%                     dAxes.ZLim = viewProfile.viewZoomAxLims(3,:);
%                 end
%
            else
                cameratoolbar('ResetCamera');
                view(dAxes, 3);
            end

            vehPosVelData = LaunchVehicleViewProfile.createVehPosVelData(subStateLogs, lvdData.script.evts, viewInFrame);
            vehAttData = LaunchVehicleViewProfile.createVehAttitudeData(vehPosVelData, lvdStateLogEntries, lvdData.script.evts, viewInFrame);

            hold(dAxes,'on');
            viewProfile.createBodyMarkerData(dAxes, subStateLogs, viewInFrame, showSoI, viewProfile.meshEdgeAlpha, lvdData.script.evts);
            viewProfile.createTrajectoryMarkerData(subStateLogs, lvdData.script.evts);
            viewProfile.createBodyAxesData(vehPosVelData, vehAttData); %lvdStateLogEntries, lvdData.script.evts, viewInFrame
            viewProfile.createSunLightSrc(dAxes, viewInFrame);
            viewProfile.createGroundObjMarkerData(dAxes, lvdStateLogEntries, vehPosVelData, lvdData.script.evts, viewInFrame, celBodyData);
            viewProfile.createCentralBodyData(viewCentralBody, hCBodySurfXForm, viewInFrame);
            viewProfile.createPointData(viewInFrame, subStateLogs, lvdData.script.evts);
            viewProfile.createVectorData(viewInFrame, subStateLogs, lvdData.script.evts);
            viewProfile.createRefFrameData(viewInFrame, subStateLogs, lvdData.script.evts);
            viewProfile.createAngleData(viewInFrame, subStateLogs, lvdData.script.evts);
            viewProfile.createPlaneData(viewInFrame, subStateLogs, lvdData.script.evts);
            viewProfile.createSensorData(lvdStateLogEntries, vehPosVelData, vehAttData, viewInFrame);
            viewProfile.createSensorTargetData(viewInFrame);

            viewProfile.configureTimeSlider(minTime, maxTime, subStateLogs, handles, app);
            hold(dAxes,'off');

            switch viewProfile.projType
                case ViewProjectionTypeEnum.Orthographic
                    camproj(dAxes, 'orthographic');

                case ViewProjectionTypeEnum.Perspective
                    camproj(dAxes, 'perspective');

                otherwise
                    error('Unknown projection type: %s', viewProfile.projType.name);

            end
        end
    end
end

function [childrenHGs] = plotSubStateLog(subStateLog, prevSubStateLog, lvdData, dAxes)
    if(isempty(subStateLog))
        childrenHGs = [];
        return;
    end

%     bodyID = subStateLog(1,8);
%     bodyInfo = getBodyInfoByNumber(bodyID, celBodyData);

    eventNum = subStateLog(1,13);
    event = lvdData.script.getEventForInd(eventNum);
    if(isempty(event))
        childrenHGs = [];
        return;
    end

    plotLineColor = event.colorLineSpec.color.color;
    plotLineStyle = event.colorLineSpec.lineSpec.linespec;
    plotLineWidth = event.colorLineSpec.lineWidth;
    plotMarkerType = event.colorLineSpec.markerSpec.shape;
    plotMarkerSize = event.colorLineSpec.markerSize;
    plotMethodEnum = event.plotMethod;

    hold(dAxes,'on');

    switch plotMethodEnum
        case EventPlottingMethodEnum.PlotContinuous
            t = [prevSubStateLog(end,1);subStateLog(1:end,1)];
            x = [prevSubStateLog(end,2);subStateLog(1:end,2)];
            y = [prevSubStateLog(end,3);subStateLog(1:end,3)];
            z = [prevSubStateLog(end,4);subStateLog(1:end,4)];

            [~,I] = sort(t);
            x = x(I);
            y = y(I);
            z = z(I);

        case EventPlottingMethodEnum.SkipFirstState
            x = subStateLog(2:end,2);
            y = subStateLog(2:end,3);
            z = subStateLog(2:end,4);
        case EventPlottingMethodEnum.DoNotPlot
            x = [];
            y = [];
            z = [];
        otherwise
            error('Unknown event plotting method enum: %s', plotMethodEnum.name);
    end

    plot3(dAxes, x, y, z, 'Color', plotLineColor, 'LineStyle', plotLineStyle, 'LineWidth',plotLineWidth, 'Marker',plotMarkerType, 'MarkerSize',plotMarkerSize, 'MarkerEdgeColor','none', 'MarkerFaceColor',plotLineColor);
    childrenHGs = cell(0,4);
end

function plotBodyFixedGrid(dAxes, bodyInfo)
    r = 1.2*bodyInfo.radius;
    rTxt = 1.3*bodyInfo.radius;

    %draw longitude circle and text
    th = linspace(0, 2*pi, 100);
    xunit = r * cos(th);
    yunit = r * sin(th);
    patch(dAxes, 'XData',xunit,'YData',yunit, 'FaceColor', 'k', 'FaceAlpha',0.15);

    th = linspace(0,2*pi - (1/12)*2*pi, 12);
    xunit = r * cos(th);
    yunit = r * sin(th);

    xToPlot = [];
    yToPlot = [];
    for(i=1:length(th))
        xToPlot = [xToPlot, 0, xunit(i), NaN]; %#ok<AGROW>
        yToPlot = [yToPlot, 0, yunit(i), NaN]; %#ok<AGROW>
    end

    plot(dAxes, xToPlot, yToPlot, 'k');

    xunitTxt = rTxt * cos(th);
    yunitTxt = rTxt * sin(th);

    for(i=1:length(xunitTxt))
        text(dAxes, xunitTxt(i),yunitTxt(i),sprintf('%.0f%s', rad2deg(th(i)), char(176)));
    end

    %draw latitude circle and text
    th = linspace(-pi/2, pi/2, 100);
    xunit = r * cos(th);
    yunit = zeros(size(th));
    zunit = r * sin(th);
    patch(dAxes, 'XData',xunit,'YData',yunit,'ZData',zunit, 'FaceColor', 'k', 'FaceAlpha',0.15);

    th = rad2deg(linspace(-pi/2, pi/2, 7));
    xunit = r * cosd(th);
    yunit = zeros(size(th));
    zunit = r * sind(th);

    xToPlot = [];
    yToPlot = [];
    zToPlot = [];
    for(i=1:length(th))
        xToPlot = [xToPlot, 0, xunit(i), NaN]; %#ok<AGROW>
        yToPlot = [yToPlot, 0, yunit(i), NaN]; %#ok<AGROW>
        zToPlot = [zToPlot, 0, zunit(i), NaN]; %#ok<AGROW>
    end

    plot3(dAxes, xToPlot, yToPlot, zToPlot, 'k');

    xunitTxt = rTxt * cosd(th);
    yunitTxt = zeros(size(th));
    zunitTxt = rTxt * sind(th);

    for(i=1:length(xunitTxt))
        text(dAxes, xunitTxt(i),yunitTxt(i),zunitTxt(i), sprintf('%.0f%s', th(i), char(176)));
    end
end

author: Arrowstar

repo: ksptot


text:
function gpmf = gpmf_constant(varargin)
%GPMF_CONSTANT  Create a constant mean function
%
%  Description
%    GPMF = GPMF_CONSTANT('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates constant mean function structure in which the named
%    parameters have the specified values. Any unspecified
%    parameters are set to default values.
%
%    GPMF = GPMF_CONSTANT(GPMF,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a mean function structure with the named parameters
%    altered with the specified values.
%
%    Parameters for constant mean function
%      constant          - constant value for the constant
%                          base function (default 1)
%      prior_mean        - prior mean (scalar or vector) for base
%                          functions' weight prior (default 0)
%      prior_cov         - prior covariances (scalar or vector)
%                          for base functions' prior corresponding
%                          each selected input dimension. In
%                          multiple dimension case prior_cov is a
%                          struct containing scalars or vectors.
%                          The covariances must all be either
%                          scalars (diagonal cov.matrix) or
%                          vectors (for non-diagonal cov.matrix)
%                          (default 100)
%
%  See also
%    GP_SET, GPMF_LINEAR, GPMF_SQUARED
%

% Copyright (c) 2010 Tuomas Nikoskinen
% Copyright (c) 2011 Jarno Vanhatalo

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'GPMF_CONSTANT';
  ip.addOptional('gpmf', [], @isstruct);
  ip.addParamValue('constant',1, @(x) isvector(x) && all(x>0));
  ip.addParamValue('prior_mean',0, @(x) isvector(x));
  ip.addParamValue('prior_cov',100, @(x) isvector(x));
  ip.addParamValue('mean_prior', [], @isstruct);
  ip.addParamValue('cov_prior', [], @isstruct);
  ip.parse(varargin{:});
  gpmf=ip.Results.gpmf;

  if isempty(gpmf)
    % Initialize a mean function
    init=true;
    gpmf.type = 'gpmf_constant';
  else
    % Modify a mean function
    if ~isfield(gpmf,'type') && isequal(gpmf.type,'gpmf_constant')
      error('First argument does not seem to be a constant mean function')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('type',ip.UsingDefaults)
    gpmf.constant = ip.Results.constant;
  end
  if init || ~ismember('prior_mean',ip.UsingDefaults)
    gpmf.b=ip.Results.prior_mean(:)';
  end
  if init || ~ismember('prior_mean',ip.UsingDefaults)
    gpmf.B=ip.Results.prior_cov(:)';
  end
  if init || ~ismember('mean_prior',ip.UsingDefaults)
    gpmf.p.b=ip.Results.mean_prior;
  end
  if init || ~ismember('cov_prior',ip.UsingDefaults)
    gpmf.p.B=ip.Results.cov_prior;
  end
  if init
    % Set the function handles to the nested functions
    gpmf.fh.geth = @gpmf_geth;
    gpmf.fh.pak = @gpmf_pak;
    gpmf.fh.unpak = @gpmf_unpak;
    gpmf.fh.lp = @gpmf_lp;
    gpmf.fh.lpg = @gpmf_lpg;
    gpmf.fh.recappend = @gpmf_recappend;
  end

end

function h = gpmf_geth(gpmf, x)
%GPMF_GETH  Calculate the base function values for a given input.
%
%  Description
%    H = GPMF_GETH(GPMF,X) takes in a mean function structure
%    GPMF and inputs X. The function returns a row vector of
%    length(X) containing the constant value which is by default
%    1.

  constant=gpmf.constant;
  h = repmat(constant,1,length(x(:,1)));

end

function [w, s] = gpmf_pak(gpmf, w)
%GPMF_PAK  Combine GP mean function parameters into one vector
%
%  Description
%    W = GPCF_LINEAR_PAK(GPCF) takes a covariance function
%    structure GPCF and combines the covariance function
%    parameters and their hyperparameters into a single row
%    vector W.
%
%       w = [ log(gpcf.coeffSigma2)
%             (hyperparameters of gpcf.coeffSigma2)]'
%
%  See also
%    GPCF_LINEAR_UNPAK

  w = []; s = {};
  if ~isempty(gpmf.p.b)
    w = gpmf.b;
    if numel(gpmf.b)>1
      s = [s; sprintf('gpmf_constant.b x %d',numel(gpmf.b))];
    else
      s = [s; 'gpmf_constant.b'];
    end
    % Hyperparameters of coeffSigma2
    [wh sh] = gpmf.p.b.fh.pak(gpmf.p.b);
    w = [w wh];
    s = [s; sh];
  end

  if ~isempty(gpmf.p.B)
    w = [w log(gpmf.B)];
    if numel(gpmf.B)>1
      s = [s; sprintf('log(gpmf_constant.B x %d)',numel(gpmf.B))];
    else
      s = [s; 'log(gpmf_constant.B)'];
    end
    % Hyperparameters of coeffSigma2
    [wh sh] = gpmf.p.B.fh.pak(gpmf.p.B);
    w = [w wh];
    s = [s; sh];
  end

end

function [gpmf, w] = gpmf_unpak(gpmf, w)
%GPCF_LINEAR_UNPAK  Sets the mean function parameters
%                   into the structure
%
%  Description
%    [GPCF, W] = GPMF_UNPAK(GPCF, W) takes a covariance
%    function structure GPCF and a hyper-parameter vector W, and
%    returns a covariance function structure identical to the
%    input, except that the covariance hyper-parameters have been
%    set to the values in W. Deletes the values set to GPCF from
%    W and returns the modified W.
%
%    Assignment is inverse of
%       w = [ log(gpcf.coeffSigma2)
%             (hyperparameters of gpcf.coeffSigma2)]'
%
%  See also
%   GPCF_LINEAR_PAK

  gpp=gpmf.p;

  if ~isempty(gpp.b)
    i2=length(gpmf.b);
    i1=1;
    gpmf.b = w(i1:i2);
    w = w(i2+1:end);

    % Hyperparameters of coeffSigma2
    [p, w] = gpmf.p.b.fh.unpak(gpmf.p.b, w);
    gpmf.p.b = p;
  end

  if ~isempty(gpp.B)
    i2=length(gpmf.B);
    i1=1;
    gpmf.B = exp(w(i1:i2));
    w = w(i2+1:end);

    % Hyperparameters of coeffSigma2
    [p, w] = gpmf.p.B.fh.unpak(gpmf.p.B, w);
    gpmf.p.B = p;
  end

end

function lp = gpmf_lp(gpmf)
%GPCF_SEXP_LP  Evaluate the log prior of covariance function parameters
%
%  Description
%
%  See also

% Evaluate the prior contribution to the error. The parameters that
% are sampled are transformed, e.g., W = log(w) where w is all
% the "real" samples. On the other hand errors are evaluated in
% the W-space so we need take into account also the Jacobian of
% transformation, e.g., W -> w = exp(W). See Gelman et.al., 2004,
% Bayesian data Analysis, second edition, p24.
  lp = 0;
  gpp=gpmf.p;

  if ~isempty(gpmf.p.b)
    lp = lp + gpp.b.fh.lp(gpmf.b, ...
                   gpp.b);
  end

  if ~isempty(gpp.B)
    lp = lp + gpp.B.fh.lp(gpmf.B, ...
                   gpp.B) +sum(log(gpmf.B));
  end
end

function [lpg_b, lpg_B] = gpmf_lpg(gpmf)
%GPCF_SEXP_LPG  Evaluate gradient of the log prior with respect
%               to the parameters.
%
%  Description
%    LPG = GPCF_SEXP_LPG(GPCF) takes a covariance function
%    structure GPCF and returns LPG = d log (p(th))/dth, where th
%    is the vector of parameters.
%
%  See also
%    GPCF_SEXP_PAK, GPCF_SEXP_UNPAK, GPCF_SEXP_LP, GP_G

  lpg_b=[]; lpg_B=[];
  gpp=gpmf.p;

  if ~isempty(gpmf.p.b)
    lll = length(gpmf.b);
    lpgs = gpp.b.fh.lpg(gpmf.b, gpp.b);
    lpg_b = [lpgs(1:lll) lpgs(lll+1:end)];  %.*gpmf.b+1
  end

  if ~isempty(gpmf.p.B)
    lll = length(gpmf.B);
    lpgs = gpp.B.fh.lpg(gpmf.B, gpp.B);
    lpg_B = [lpgs(1:lll).*gpmf.B+1 lpgs(lll+1:end)];
  end
end

function recmf = gpmf_recappend(recmf, ri, gpmf)
%RECAPPEND  Record append
%
%  Description
%
%  See also
%    GP_MC and GP_MC -> RECAPPEND

% Initialize record
  if nargin == 2
    recmf.type = 'gpmf_constant';

    % Initialize parameters
    recmf.b= [];
    recmf.B = [];

    % Set the function handles
    recmf.fh.geth = @gpmf_geth;
    recmf.fh.pak = @gpmf_pak;
    recmf.fh.unpak = @gpmf_unpak;
    recmf.fh.lp = @gpmf_lp;
    recmf.fh.lpg = @gpmf_lpg;
    recmf.fh.recappend = @gpmf_recappend;

    recmf.p=[];
    recmf.p.b=[];
    recmf.p.B=[];
    if isfield(ri.p,'b') && ~isempty(ri.p.b)
      recmf.p.b = ri.p.b;
    end
    if ~isempty(ri.p.B)
      recmf.p.B = ri.p.B;
    end
    return
  end

  gpp = gpmf.p;

  % record magnSigma2
  if ~isempty(gpmf.b)
    recmf.b(ri,:)=gpmf.b;
    if ~isempty(recmf.p.b)
      recmf.p.b = gpp.b.fh.recappend(recmf.p.b, ri, gpmf.p.b);
    end
  elseif ri==1
    recmf.b=[];
  end

  if ~isempty(gpmf.B)
    recmf.B(ri,:)=gpmf.B;
    if ~isempty(recmf.p.B)
      recmf.p.B = gpp.B.fh.recappend(recmf.p.B, ri, gpmf.p.B);
    end
  elseif ri==1
    recmf.B=[];
  end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function triangle_to_fem ( prefix )

%*****************************************************************************80
%
%% MAIN is the main program for TRIANGLE_TO_FEM.
%
%  Discussion:
%
%    TRIANGLE_TO_FEM converts mesh data from TRIANGLE format to FEM format.
%
%  Usage:
%
%    triangle_to_fem prefix
%
%    where 'prefix' is the common filename prefix:
%
%    * 'prefix'.node contains the triangle node coordinates,
%    * 'prefix'.ele contains the triangle element node connectivity.
%    * 'prefix'_nodes.txt will contain the node coordinates.
%    * 'prefix'_elements.txt will contain the element node connectivity.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 October 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGLE_TO_FEM\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '  Read a mesh description created by TRIANGLE:\n' );
  fprintf ( 1, '  * "prefix".node, node coordinates.\n' );
  fprintf ( 1, '  * "prefix".ele, element connectivity.\n' );
  fprintf ( 1, '  Write two simple FEM files listing nodes and elements.\n' );
  fprintf ( 1, '  * "prefix"_nodes.txt, node coordinates.\n' );
  fprintf ( 1, '  * "prefix"_elements.txt, element connectivity.\n' );
%
%  Get the filename prefix.
%
  if ( nargin < 1 )

    prefix = input ( 'Enter the filename prefix:  ' );

  end
%
%  Create the filenames.
%
  triangle_node_filename = strcat ( prefix, '.node' );
  triangle_element_filename = strcat ( prefix, '.ele' );
  fem_node_filename = strcat ( prefix, '_nodes.txt' );
  fem_element_filename = strcat ( prefix, '_elements.txt' );
%
%  Read the triangle node size information.
%
  [ node_num, m, node_att_num, node_marker_num ] = node_size_read ( ...
    triangle_node_filename );
%
%  Read the triangle element size information.
%
  [ element_num, element_order, element_att_num ] = element_size_read ( ...
    triangle_element_filename );
%
%  Report sizes.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Size information from TRIANGLE files:\n' );
  fprintf ( 1, '  Spatial dimension M = %d\n', m );
  fprintf ( 1, '  Number of nodes NODE_NUM = %d\n', node_num );
  fprintf ( 1, '  NODE_ATT_NUM = %d\n', node_att_num );
  fprintf ( 1, '  NODE_MARKER_NUM = %d\n', node_marker_num );
  fprintf ( 1, '  Number of elements ELEMENT_NUM = %d\n', element_num );
  fprintf ( 1, '  Element order ELEMENT_ORDER = %d\n', element_order );
  fprintf ( 1, '  ELEMENT_ATT_NUM = %d\n', element_att_num );
%
%  Read data.
%
  [ node_x, node_att, node_marker ] = node_data_read ( ...
    triangle_node_filename, node_num, m, node_att_num, node_marker_num );

  [ element_node, element_att ] = element_data_read ( ...
    triangle_element_filename, element_num, element_order, element_att_num );
%
%  Write data.
%
  r8mat_write ( fem_node_filename, m, node_num, node_x );

  i4mat_write ( fem_element_filename, element_order, element_num, ...
    element_node );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGLE_TO_FEM:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function [ element_node, element_att ] = element_data_read ( ...
  element_filename, element_num, element_order, element_att_num )

%*****************************************************************************80
%
%% ELEMENT_DATA_READ reads the data from an element file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 October 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ELEMENT_FILENAME, the name of the file.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, integer ELEMENT_ATT_NUM, number of element attributes
%    listed on each node record.
%
%    Output, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM),
%    the indices of the nodes that make up each element.
%
%    Output, real ELEMENT_ATT(ELEMENT_ATT_NUM,ELEMENT_NUM), the
%    attributes of each element.
%
  element = 0;

  input = fopen ( element_filename, 'rt' );

  if ( input < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'ELEMENT_DATA_READ - Fatal error!\n' );
    fprintf ( 1, '  Unable to open file.\n' );
    error ( 'ELEMENT_DATA_READ - Fatal error!' );
  end
%
%  Create format.
%
  format = '%d';
  for i = 1 : element_order
    format = strcat ( format, '%d' );
  end
  for i = 1 : element_att_num
    format = strcat ( format, '%g' );
  end
%
%  Create arrays.
%
  element_node = zeros ( element_order, element_num );
  element_att = zeros ( element_att_num, element_num );
%
%  Read values and place in the arrays.
%
  while ( 1 )

    text = fgets ( input );

    if ( text == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'ELEMENT_DATA_READ - Fatal error!\n' );
      fprintf ( 1, '  Unexpected end of file while reading.\n' );
      error ( 'ELEMENT_DATA_READ - Fatal error!' );
    end

    if ( s_len_trim ( text ) == 0 )
      continue
    end

    if ( text(1) == '#' )
      continue
    end
%
%  Ignore the dimension line.
%
    if ( element == 0 )

    else

      [ value, count ] = sscanf ( text, format );

      if ( count < 1 + element_order + element_att_num )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'ELEMENT_DATA_READ - Fatal error!\n' );
        fprintf ( 1, '  Line did not have enough data.\n' );
        error ( 'ELEMENT_DATA_READ - Fatal error!' );
      end

      j = 1;

      for i = 1 : element_order
        j = j + 1;
        element_node(i,element) = value(j);
      end

      for i = 1 : element_att_num
        element_att(i,element) = value(j);
      end

    end

    element = element + 1;

    if ( element_num < element )
      break
    end

  end

  fclose ( input );

  return
end
function [ element_num, element_order, element_att_num ] = ...
  element_size_read ( element_filename )

%*****************************************************************************80
%
%% ELEMENT_SIZE_READ reads the header information from an element file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 October 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ELEMENT_FILENAME, the name of the
%    element file.
%
%    Output, integer ELEMENT_NUM, the number of elements.
%
%    Output, integer ELEMENT_ORDER, the order of the elements.
%
%    Output, integer ELEMENT_ATT_NUM, the number of
%    element attributes.
%
  input = fopen ( element_filename, 'rt' );

  if ( input < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'ELEMENT_SIZE_READ - Fatal error!\n' );
    fprintf ( 1, '  Unable to open file.\n' );
    error ( 'ELEMENT_SIZE_READ - Fatal error!' );
  end

  while ( 1 )

    text = fgets ( input );

    if ( text == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'ELEMENT_SIZE_READ - Fatal error!\n' );
      fprintf ( 1, '  Unexpected end of file while reading.\n' );
      error ( 'ELEMENT_SIZE_READ - Fatal error!' );
    end

    if ( s_len_trim ( text ) == 0 )
      continue
    end

    if ( text(1) == '#' )
      continue
    end

    [ value, count ] = sscanf ( text, '%d%d%d%d' );

    if ( count < 3 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'ELEMENT_SIZE_READ - Fatal error!\n' );
      fprintf ( 1, '  Header line did not have 3 value integers.\n' );
      error ( 'ELEMENT_SIZE_READ - Fatal error!' );
    end

    element_num = value(1);
    element_order = value(2);
    element_att_num = value(3);

    break

  end

  fclose ( input );

  return
end
function i4mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Discussion:
%
%    An I4MAT is an array of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function [ node_coord, node_att, node_marker ] = node_data_read ( ...
  node_filename, node_num, node_dim, node_att_num, node_marker_num )

%*****************************************************************************80
%
%% NODE_DATA_READ reads the data from a node file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 October 2014
%
%  Author:
%
%    John Burkardt
%
% Parameters:
%
%   Input, string NODE_FILENAME, the name of the node file.
%
%   Input, integer NODE_NUM, the number of nodes.
%
%   Input, integer NODE_DIM, the spatial dimension.
%
%   Input, integer NODE_ATT_NUM, number of node attributes
%   listed on each node record.
%
%   Input, integer NODE_MARKER_NUM, 1 if every node record
%   includes a final boundary marker value.
%
%   Output, real NODE_COORD(NODE_DIM,NODE_NUM), the nodal
%   coordinates.
%
%   Output, real NODE_ATT(NODE_ATT_NUM,NODE_NUM), the nodal
%   attributes.
%
%   Output, integer NODE_MARKER(NODE_MARKER_NUM,NODE_NUM), the
%   node markers.
%
  node = 0;

  input = fopen ( node_filename, 'rt' );

  if ( input < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'NODE_DATA_READ - Fatal error!\n' );
    fprintf ( 1, '  Unable to open file.\n' );
    error ( 'NODE_DATA_READ - Fatal error!' );
  end
%
%  Create format.
%
  format = '%d';
  for i = 1 : node_dim
    format = strcat ( format, '%g' );
  end
  for i = 1 : node_att_num
    format = strcat ( format, '%g' );
  end
  for i = 1 : node_marker_num
    format = strcat ( format, '%d' );
  end
%
%  Create arrays.
%
  node_coord = zeros ( node_dim, node_num );
  node_att = zeros ( node_att_num, node_num );
  node_marker = zeros ( node_marker_num, node_num );
%
%  Read values and place in the arrays.
%
  while ( 1 )

    text = fgets ( input );

    if ( text == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'NODE_DATA_READ - Fatal error!\n' );
      fprintf ( 1, '  Unexpected end of file while reading.\n' );
      error ( 'NODE_DATA_READ - Fatal error!' );
    end

    if ( s_len_trim ( text ) == 0 )
      continue
    end

    if ( text(1) == '#' )
      continue
    end
%
%  Ignore the dimension line.
%
    if ( node == 0 )

    else

      [ value, count ] = sscanf ( text, format );

      if ( count < 1 + node_dim + node_att_num + node_marker_num )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'NODE_DATA_READ - Fatal error!\n' );
        fprintf ( 1, '  Line did not have enough data.\n' );
        error ( 'NODE_DATA_READ - Fatal error!' );
      end

      j = 1;

      for i = 1 : node_dim
        j = j + 1;
        node_coord(i,node) = value(j);
      end

      for i = 1 : node_att_num
        j = j + 1;
        node_att(i,node) = value(j);
      end

      for i = 1 : node_marker_num
        j = j + 1;
        node_marker(i,node) = value(j);
      end

    end

    node = node + 1;

    if ( node_num < node )
      break
    end

  end

  fclose ( input );

  return
end
function [ node_num, node_dim, node_att_num, node_marker_num ] = ...
  node_size_read ( node_filename )

%*****************************************************************************80
%
%% NODE_SIZE_READ reads the header information from a node file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 October 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string NODE_FILENAME, the name of the node file.
%
%    Output, integer NODE_NUM, the number of nodes.
%
%    Output, integer NODE_DIM, the spatial dimension.
%
%    Output, integer NODE_ATT_NUM, number of node attributes
%    listed on each node record.
%
%    Output, integer NODE_MARKER_NUM, 1 if every node record
%    includes a final boundary marker value.
%
  input = fopen ( node_filename, 'rt' );

  if ( input < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'NODE_SIZE_READ - Fatal error!\n' );
    fprintf ( 1, '  Unable to open file.\n' );
    error ( 'NODE_SIZE_READ - Fatal error!' );
  end

  while ( 1 )

    text = fgets ( input );

    if ( text == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'NODE_SIZE_READ - Fatal error!\n' );
      fprintf ( 1, '  Unexpected end of file while reading.\n' );
      error ( 'NODE_SIZE_READ - Fatal error!' );
    end

    if ( s_len_trim ( text ) == 0 )
      continue
    end

    if ( text(1) == '#' )
      continue
    end

    [ value, count ] = sscanf ( text, '%d%d%d%d' );

    if ( count < 4 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'NODE_SIZE_READ - Fatal error!\n' );
      fprintf ( 1, '  Header line did not have 4 value integers.\n' );
      error ( 'NODE_SIZE_READ - Fatal error!' );
    end

    node_num = value(1);
    node_dim = value(2);
    node_att_num = value(3);
    node_marker_num = value(4);

    break

  end

  fclose ( input );

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Discussion:
%
%    An R8MAT is an array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
%  Alternative print statements include:
%
%     fprintf ( output_unit, '  %24.16e', table(i,j) );
%     fprintf ( output_unit, '  %14.6e', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %g', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LENGTH, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function plotData(X, y)
%PLOTDATA Plots the data points X and y into a new figure
%   PLOTDATA(x,y) plots the data points with + for the positive examples
%   and o for the negative examples. X is assumed to be a Mx2 matrix.

% Create New Figure
figure; hold on;

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the positive and negative examples on a
%               2D plot, using the option 'k+' for the positive
%               examples and 'ko' for the negative examples.
%

% =========================================================================

hold off;

end

author: Ayatans

repo: Machine-Learning-homework


text:
function out= util_unixCmd(cmd, action)
%UTIL_UNIXCMD - Run a unix command and check for errors
%
%Synopsis:
% util_unixCcmd(CMD, <ACTION>)
%
%Arguments:
% CMD:    String, the unix command to be executed
% ACTION: String, printed in case of an error

if ~isunix,
  warning('attempt to execute unix command in non-unix system');
  return;
end

if nargin<2,
  action= '';
end

[stat,out]= unix(cmd);
if stat~=0,
  error(sprintf('error %s (%s -> %s)', action, cmd, out));
end

if nargout==0,
  clear out
end

author: bbci


text:
% Test Dag2Cpdag and Cpdag2Dag functions ...

close all;
clear all;

bnet=mk_asia_bnet ;

n=8;
names={ 'A' , 'S' , 'T' , 'L' , 'B' , 'O' , 'X' , 'D' };
carre=ones(1,n);

[xx yy] = make_layout(bnet.dag);
yy=(yy-0.2)*.8/.6+.1;
xx=(xx-0.2833)*.8/.517+.1;

cpdag=dag_to_cpdag(bnet.dag) ;
dag1 = cpdag_to_dag(cpdag) ;

figure;
subplot(1,3,1), draw_graph(bnet.dag,names,carre,xx,yy); title('original ASIA');
subplot(1,3,2), draw_graph(cpdag,names,carre,xx,yy); title('cpasia=DAGtoCPDAG(asia)');
subplot(1,3,3), draw_graph(dag1,names,carre,xx,yy); title('CPDAGtoDAG(cpasia)');

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function imgOut = GammaDrago(img, drago_gamma, drago_slope, drago_start)
%
%        imgOut = GammaDrago(img, drago_gamma, drago_slope, drago_start)
%
%        This function applies gamma correction for the Drago et al.'s TMO,
%        please see DragoTMO.m
%
%        Input:
%           -img: an LDR image tone mapped with DragoTMO.m
%           -drago_gamma: is the elevation ratio of the line passing by the
%            origin and tangent to the curve
%           -drago_slope: f-stop value
%           -drago_start: is the abscissa at the point of tangency
%
%        Output:
%           -imgOut: gamma corrected img version
%
%     Copyright (C) 2013  Francesco Banterle
%
%     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 3 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, see <http://www.gnu.org/licenses/>.
%

if(~exist('drago_gamma','var'))
    drago_gamma = 2.2;
end

if(~exist('drago_slope','var'))
    drago_slope = 4.5;
end

if(~exist('drago_start','var'))
    drago_start = 0.018;
end

%applying gamma correction as in Drago et al. 2003
indx1 = find(img<=drago_start);
indx2 = find(img> drago_start);
img(indx1) =  img(indx1)*drago_slope;
img(indx2) = (img(indx2).^(0.9/drago_gamma))*1.099-0.099;

%clamping values out of the range [0.0,1.0]
imgOut = ClampImg(img, 0.0, 1.0);

end

author: banterle


text:
% main file

close all
clear all

% read Template image
im1=imread('K.bmp');
%im1=imread('S.bmp');
%im1=imread('image1.jpg');

% read Traget Image
im2=imread('letters.bmp');
%im2=imread('image2.jpg');

% apply templete matching using power of the image
result1=tmp(im1,im2);

figure,
subplot(2,2,1),imshow(im1);title('Template');
subplot(2,2,2),imshow(im2);title('Target');
subplot(2,2,3),imshow(result1);title('Matching Result using tmp');

% apply templete matching using DC components of the image
result2=tmc(im1,im2);

figure,
subplot(2,2,1),imshow(im1);title('Template');
subplot(2,2,2),imshow(im2);title('Target');
subplot(2,2,3),imshow(result2);title('Matching Result using tmc');

author: Sable

repo: mcbench-benchmarks


text:
function subset_sum_serial_test ( )

%*****************************************************************************80
%
%% SUBSET_SUM_SERIAL_TEST tests the SUBSET_SUM_SERIAL library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 March 2011
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SUBSET_SUM_SERIAL_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the SUBSET_SUM_SERIAL library.\n' );

  subset_sum_serial_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SUBSET_SUM_SERIAL_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
% Generate various time and frequency-varying test signals
%
% Usage 1:
%
%   output = gsig( data_type, f1, f2, num_samples, sig_type);
%
%
%   Inputs
%
%	data_type
%
%          'lin' linear FM
%          'quad' quadratic FM
%          'cubic' cubic FM
%          'hyp' hyperbolic FM
%
%	f1, f2
%
%          normalised frequency bounds for the signal (i.e. taking
%   	   sampling frequency as 1 Hz) and sig_type=1 for real data
%          otherwise  the result is complex.
%
%	num_samples
%
%          is the number of generated data samples.
%
%
%      sig_type
%
%          Must be 0 (analytic) or 1 (real).
%
% Usage 2:
%
%   output = gsig( 'sin', cf, mf, num_samples, sig_type, fdev );
%
%
%   Inputs
%
%	cf and mf
%
%          central and modulation frequencies.
%
%       fdev
%
%          frequency deviation.
%
% Usage 3:
%
%   output1 = gsig( 'step', f1, f2, num_samples, sig_type, ns );
%
%   Inputs
%
%	f1, f2
%
%          The normalised frequency bounds for the signal (i.e. taking
%	   sampling frequency as 1 Hz).
%
%       ns
%
%          The number of steps.
%
%
%
%  See Also: analyt
%
% TFSAP 7.0
% Copyright Prof. B. Boashash
% Qatar University, Doha
% email: tfsap.research@gmail.com

author: HeLiangHIT


text:
function B = horzcat(varargin)
%HORZCAT   Horizontally concatenate LINOP objects.

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

% Cast the result as a linop if there are any linop arguments.
B = linop( horzcat@chebmatrix(varargin{:}) );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%  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 St, Fifth Floor, Boston, MA  02110-1301  USA
function dividert
%DIVIDERT Unit test for the function DIVIDER.

%	O. Lemoine - March 1996.

N=500;

for k=1:N,
 [Nk,Mk]=divider(k);
 err(k)=abs(Nk*Mk-k);
end

if any(err),
 error('divider test 1 failed');
end;


text:
function legendlinestyles(h,markers,linestyles,linecolors)
% legendlinestyles(h,markers,linestyles,linecolors)
%==========================================================================
% Change the linestyles in an existing plot legend.
%
% Input: Handle to an existing legend object, h
%        (optional) Cell array with the new markers, markers
%        (optional) Cell array with the new linestyles, linestyles
%        (optional) Cell array with the new line and marker colors
%
% If you do not want to change a markerstyle or linestyle, supply an empty
% array {}. See examples below.
%
% If the 'linestyles' arument is ommitted, the original lines are left
% unchanged. Note, that no lines in the actual plot are affected. If emtpy,
% {}, the styles are left unchanged but the colors are changed according to
% the fourth argument.
%
% Output: None
%
%
% Examples on how to use (see also accompanying file
% legendlinestyles_test.m):
%
%  legendlinestyles(h, markers);
%  legendlinestyles(h, markers,linestyles);
%  legendlinestyles(h, markers,linestyles,linecolors);
%  legendlinestyles(h, markers,{},linecolors);
%  legendlinestyles(h, {},linestyles,linecolors);
%  legendlinestyles(h, {},{},linecolors);
%
%==========================================================================
% Version: 1.0
% Created: October 3, 2008, by Johan E. Carlson
% Last modified: May 12, 2009, by Johan E. Carlson
%==========================================================================

lines = findobj(get(h,'children'),'type','line');
m = 1;
for k = length(lines):-2:1,
    if ~isempty(markers),
        set(lines(k-1),'marker',char(markers{m}));
    end
    if (nargin >= 3) && (~isempty(linestyles)),
        set(lines(k),'linestyle',char(linestyles{m}));
    end
    if nargin == 4,
        if ischar(linecolors{m})==1,
            set(lines(k),'color',char(linecolors{m}));
            if ~isempty(markers),
                set(lines(k-1),'MarkerEdgeColor',char(linecolors{m}));
            end
        else
            set(lines(k),'color',linecolors{m});
            if ~isempty(markers),
                set(lines(k-1),'MarkerEdgeColor',linecolors{m});
            end
        end
    end
    m = m+1;
end


text:
function dftoggleaxlimctrl(dffig)
%DFTOGGLEAXLIMCTRL Toggle x and y axis limit controls on or off

%   $Revision: 1.1.6.2 $  $Date: 2004/01/24 09:35:53 $
%   Copyright 2003-2004 The MathWorks, Inc.

% Get handle to menu item, may be current object or may not
h = gcbo;
if ~isequal(get(h,'Tag'),'showaxlimctrl')
   h = findall(dffig,'Tag','showaxlimctrl');
end

% Get new state
onoff = on2off(get(h,'Checked'));
dfgetset('showaxlimctrl',onoff);

% Add or remove controls
dfaxlimctrl(dffig,onoff)

% Remove effects of controls on layout
if isequal(onoff,'off')
   dfadjustlayout(dffig);
end

% Change menu state
set(h,'Checked',onoff);

author: zouchuhang

repo: LayoutNet


text:
%%*****************************************************************
%% primLB: compute a lower bound for the exact primal
%%         optimal value.
%%
%%*****************************************************************

  function LB = primLB(blk,At,C,b,eigX,y,mu);

  if (nargin < 7); mu = 1.1; end
  Aty = sdpnalAtyfun(blk,At,y);
  Znew = ops(C,'-',Aty);

  LB0 = b'*y;
  pert = 0;
  for p = 1:size(blk,1)
     pblk = blk(p,:);
     if strcmp(pblk{1},'s')
        eigtmp = eig(Znew{p});
        idx = find(eigtmp < 0);
        Xbar = mu*max(eigX{p});
     elseif strcmp(pblk{1},'l')
        eigtmp = Znew{p};
        idx = find(eigtmp < 0);
        Xbar = mu*max(eigX{p});
     end
     numneg = length(idx);
     if (numneg)
        mineig = min(eigtmp(idx));
        fprintf('\n numneg = %3.0d,  mineigZnew = %- 3.2e',numneg,mineig);
        pert = pert + Xbar*sum(eigtmp(idx));
     end
  end
  LB = LB0 + pert;
  fprintf('\n dobj = %-10.9e  \n LB   = %-10.9e\n',LB0,LB);
%%*****************************************************************

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function [logp, yhat, res] = tapas_logrt_linear_whatworld(r, infStates, ptrans)
% Calculates the log-probability of log-reaction times y (in units of log-ms) according to the
% linear log-RT model developed with Louise Marshall and Sven Bestmann
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2014 Christoph Mathys, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Transform zetas to their native space
be0  = ptrans(1);
be1  = ptrans(2);
be2  = ptrans(3);
be3  = ptrans(4);
ze   = exp(ptrans(5));

% Initialize returned log-probabilities, predictions,
% and residuals as NaNs so that NaN is returned for all
% irregualar trials
n = size(infStates,1);
logp = NaN(n,1);
yhat = NaN(n,1);
res  = NaN(n,1);

% Weed irregular trials out from responses and inputs
y = r.y(:,1);
y(r.irr) = [];

u = r.u(:,1);
u(r.irr) = [];

% Extract trajectories of interest from infStates
mu1hat = squeeze(infStates(:,1,:,:,1));
mu1    = squeeze(infStates(:,1,:,:,3));
mu2    = squeeze(infStates(:,2,:,:,3));
sa2    = squeeze(infStates(:,2,:,:,4));
mu3    = squeeze(infStates(:,3,1,1,3));

% Surprise
% ~~~~~~~~

% mu1 contains the actually occurring transition -> multiply with
% mu1hat to get probability of that transition (other elements are
% zero)
otp    = mu1.*mu1hat; % observed transition probabilities (3-dim)
otps3  = sum(otp, 3, 'omitnan');      % sum over 3rd dim
otps23 = sum(otps3, 2, 'omitnan');    % sum over 2nd dim

surp = -log(otps23);
surp(r.irr) = [];

% Expected uncertainty
% ~~~~~~~~~~~~~~~~~~~~
euo    = mu1.*sa2;    % expected uncertainty of observed transition (3-dim)
euos3  = sum(euo, 3, 'omitnan');      % sum over 3rd dim
euos23 = sum(euos3, 2, 'omitnan');    % sum over 2nd dim

to     = mu1.*mu2;    % tendency of observed transition (3-dim)
tos3   = sum(to, 3, 'omitnan');       % sum over 3rd dim
tos23  = sum(tos3, 2, 'omitnan');     % sum over 2nd dim

eu = tapas_sgm(tos23,1).*(1-tapas_sgm(tos23,1)).*euos23; % transform down to 1st level
eu(r.irr) = [];

% Unexpected uncertainty
% ~~~~~~~~~~~~~~~~~~~~~~
ueu = tapas_sgm(tos23,1).*(1-tapas_sgm(tos23,1)).*exp(mu3); % transform down to 1st level
ueu(r.irr) = [];

% Calculate predicted log-reaction time
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
logrt = be0 +be1.*surp +be2.*eu +be3.*ueu;

% Calculate log-probabilities for non-irregular trials
% Note: 8*atan(1) == 2*pi (this is used to guard against
% errors resulting from having used pi as a variable).
reg = ~ismember(1:n,r.irr);
logp(reg) = -1/2.*log(8*atan(1).*ze) -(y-logrt).^2./(2.*ze);
yhat(reg) = logrt;
res(reg) = y-logrt;

return;

author: translationalneuromodeling

repo: tapas


text:
function lpp_value_test ( )

%*****************************************************************************80
%
%% LPP_VALUE_TEST tests LPP_VALUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LPP_VALUE_TEST:\n' );
  fprintf ( 1, '  LPP_VALUE evaluates a Legendre product polynomial.\n' );

  m = 3;
  n = 1;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( m, -1.0, +1.0, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Evaluate at X = ( %g, %g, %g )\n', x(1:3) )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Rank  I1  I2  I3:  L(I1,X1)*L(I2,X2)*L(I3,X3)  P(X1,X2,X3)\n' );
  fprintf ( 1, '\n' );

  for rank = 1 : 20

    l = comp_unrank_grlex ( m, rank );
%
%  Evaluate the LPP directly.
%
    v1 = lpp_value ( m, n, l, x );
%
%  Convert the LPP to a polynomial, and reevaluate.
%
    o_max = prod ( floor ( ( l(1:m) + 2 ) / 2 ) );

    [ o, c, e ] = lpp_to_polynomial ( m, l, o_max );

    v2 = polynomial_value ( m, o, c, e, n, x );
%
%  Compare results.
%
    fprintf ( 1, '  %4d  %2d  %2d  %2d  %12g  %12g\n', rank, l(1:m), v1, v2 );

  end

  return
end


text:
function retVal=setProcRoundingMode(roundMode)
%%SETPROCROUNDINGMODE Set the rounding mode used in the processor. This
%                     specifies how results should be truncated to a
%                     finite number of bits when performing floating point
%                     arithmetic. Control of the rounding mode is useful
%                     when implementing routines utilizing interval
%                     algebra. Note that this does not appear to change
%                     the rounding mode used in mex files called from
%                     Matlab. When running things in parallel, the rounding
%                     mode will have to be set for each parallel process.
%                     For example, set it at the start of a parfor loop,
%                     not outside the loop.
%
%INPUTS: roundMode An integer specifying the rounding mode to use.
%                  Possible values are
%                  0 Rounding is done towards negative infinity.
%                  1 Rounding is done towards zero.
%                  2 Rounding is done to the nearest value.
%                  3 Rounding is done towards positive infinity.
%
%OUTPUTS: retVal This is zero if the rounding direction was successfully
%                set and a nonzero value otherwise.
%
%The rounding modes are standardized in [1]. The native Matlab version of
%this function calls the undocumented system_dependent('setround',dir)
%function in Matlab. The version that can be compiled in C just calls the
%standard C function fesetround.
%
%EXAMPLE:
% setProcRoundingMode(0);
% 1+eps(0)>1
% %The result should be false (0).
% setProcRoundingMode(3);
% 1+eps(0)>1
% %The result should be true (1).
%
%REFERENCES:
%[1] IEEE Standard 754-1985 for Binary Floating-point Arithmetic, IEEE,
%    (1985).
%
%January 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

switch(roundMode)
    case 0%Rounding is done towards negative infinity.
        system_dependent('setround',-Inf);
    case 1%Rounding is done towards zero.
        system_dependent('setround',0);
    case 2%Rounding is done to the nearest value.
        system_dependent('setround',0.5);
    case 3%Rounding is done towards positive infinity.
        system_dependent('setround',Inf);
    otherwise
        error('Unknown rounding mode specified.')
end

if(nargout>0)
    actualMode=getProcRoundingMode();
    retVal=(actualMode~=roundMode);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [item, inputs] = fillvals(item, inputs, infcn)

% function [item, inputs] = fillvals(item, inputs, infcn)
% If ~all_set_item, try to set item.val to the items listed in inputs{1}.
% inputs{1} should be a cell array of indices into item.values. For
% cfg_choice items, this list should only contain one item.
% Validity checks are performed through setval. If inputs{1} is not
% suitable for this item, it is discarded. If infcn is a function handle,
% [val sts] = infcn(item)
% will be called to obtain a value for this item. This call will be
% repeated until either val can be assigned to item or sts is true.
%
% This function is identical for all cfg_intree classes.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: fillvals.m 5678 2013-10-11 14:58:04Z volkmar $

rev = '$Rev: 5678 $'; %#ok

% Set item itself
if ~all_set_item(item)
    if ~isempty(inputs)
        for k = 1:numel(inputs{1})
            item = setval(item, [inputs{1}{k} Inf], false);
        end;
        inputs = inputs(2:end);
    end;
    if ~all_set_item(item) && ~isempty(infcn) && subsasgn_check_funhandle(infcn)
        sts = false;
        while ~sts && ~all_set_item(item)
            [val, sts] = feval(infcn, item);
            if sts
                for k = 1:numel(val)
                    item = setval(item, [val{k} Inf], false);
                end;
            end;
        end;
    end;
end;
% Set children
citems = subsref(item, substruct('.','val'));
for k = 1:numel(citems)
    if ~all_set(citems{k})
        [citems{k}, inputs] = fillvals(citems{k}, inputs, infcn);
    end;
end;
item = subsasgn(item, substruct('.','val'), citems);

author: spm


text:
% mycorr2 modified version of the 2D correlation
%         for the use with im2col and col2im
%         see GETPOINT
%
%
% $Id: mycorr2.m,v 2.0 2003/06/19 12:07:11 svoboda Exp $

% Note: It written in order to gain speed. The clarity of the code suffers accordingly

function R = mycorr2(X,G,Gn,Gn2)

% Gn  = G-mean(G);
% Gn2 = sqrt(sum(Gn.^2));

mX	= repmat(mean(X),size(X,1),1);
mXn = X - mX;
smX	= sum(mXn.^2);

numerator = (mXn'*Gn)';
denominator = smX*Gn2;

R = numerator./denominator;

return

author: strawlab

repo: MultiCamSelfCal


text:
function h = cutfunc(x,a,b)
c1 = (x <= a);
c2 = (x>a & x<b);
c3 = (x >= b);
h = 0*c1+ (x-a)/(b-a).*c2 + 1*c3;
% c1 = (x <= a);
% c2 = (x>a & x<b);
% c3 = (x >= b);
% h = 0*c1+ 0.5*c2 + 1*c3;

%  for i=2:nx
% %
%      h(1,i)= (x(1,i)-1)/(-1);
% %
% %     h(1,1)=1;
%  end

repo: Aero-matlab


text:
%--- help for abstvar/setOptions ---
%
%  INTERNAL FUNCTION: Sets options for the VAR
%
%  ::
%
%     var_obj = set_inputs(var_obj,varargin);
%
%  Args:
%     var_obj (var object): VAR object
%     varargin: Options to set. Must come in pairs
%
%        - **'data'**: Data for the VAR. Can be either a struct or ts object
%        - **'prior'**: Prior to use if BVAR. Available priors are
%
%           - minnesota
%           - jeffrey
%           - nw
%           - normal-wishart
%           - inw
%           - sz
%
%        - **'linear_restrictions'**: Linear restrictions. For the format, see XXXXXXXXX
%
%  Returns:
%     :
%
%        - var_obj (var object): New VAR with the given inputs updated.
%
%
%    Other functions named setOptions
%
%       AbstractPortfolio/setOptions    PortfolioCVaR/setOptions
%       Portfolio/setOptions            PortfolioMAD/setOptions
%

author: jmaih


text:
function vol = pr_blobs2vol(xyz,vals,mat)
% Take XYZ matrix and values and return SPM matrix vol struct
% FORMAT vol = pr_blobs2vol(xyz,vals,mat)
%
% Inputs
% xyz      - 3xN X Y Z coordinate matrix (in voxels)
% vals     - 1xN values, one per coordinate
% mat      - 4x4 voxel->world space transformation
%
% Outputs
% vol      - vol struct, with matrix data 'imgdata' field
%__________________________________________________________________________

% Matthew Brett
% $Id: pr_blobs2vol.m 6623 2015-12-03 18:38:08Z guillaume $

vol = [];
if ~isempty(xyz),
    rcp      = round(xyz);
    vol.dim  = max(rcp,[],2)';
    off      = rcp(1,:) + vol.dim(1)*(rcp(2,:)-1+vol.dim(2)*(rcp(3,:)-1));
    vol.imgdata = zeros(vol.dim)+NaN;
    vol.imgdata(off) = vals;
    vol.imgdata      = reshape(vol.imgdata,vol.dim);
    vol.mat = mat;
end


text:
function Z = jacobiansparsityfromnonlinear(model,includeLinear)

if nargin == 1
    % IPOPT appends linear constraints, and that sparsity structure should
    % thus be included. KNITRO does not append linears to the structure
    includeLinear = 1;
end

m = length(model.lb);
allA=[model.Anonlinineq];
if any(model.K.q)
    top = startofSOCPCone(model.K);
    allQ = [];
    for i = 1:length(model.K.q)
        allQ = [allQ;any(model.F_struc(top:top+model.K.q(i)-1,2:end))];
        top = top + model.K.q(i);
    end
    allA = [allA;allQ];
end
if any(model.K.e)
    top = startofEXPCone(model.K);
    allQ = [];
    for i = 1:(model.K.e)
        allQ = [allQ;any(model.F_struc(top:top+3-1,2:end))];
        top = top + 3;
    end
    allA = [allA;allQ];
end
if any(model.K.p)
    top = startofPOWCone(model.K);
    allQ = [];
    for i = 1:length(model.K.p)
        allQ = [allQ;any(model.F_struc(top:top+model.K.p(i)-1,2:end))];
        top = top + model.K.p(i);
    end
    allA = [allA;allQ];
end
if any(model.K.s)
    top = startofSDPCone(model.K);
    allQ = [];
    for i = 1:length(model.K.s)
        s = any(model.F_struc(top:top+model.K.s(i)^2-1,2:end));
        allQ = [allQ;repmat(s,model.K.s(i),1)];
        top = top + model.K.s(i)^2;
    end
    allA = [allA;allQ];
end

allA = [allA;model.Anonlineq];
depends = allA | allA;
if length(model.linearindicies) == size(allA,2)
    % It is linear + possibly linear cone
    jacobianstructure = double(depends);
else
    jacobianstructure = spalloc(size(allA,1),m,0);
    for i = 1:size(depends,1)
        vars = find(depends(i,:));
        [ii,vars] = find(model.deppattern(vars,:));
        vars = unique(vars);
        s = size(jacobianstructure,1);
        if 0
        [~,loc] = ismember(vars,model.linearindicies);
        jacobianstructure(i,find(loc)) = 1;
        else
        for j = 1:length(vars)
            jacobianstructure(i,find(vars(j) == model.linearindicies)) = 1;
        end
        end
    end
end
if includeLinear
    allA=[model.A; model.Aeq];
    depends = allA | allA;
    jacobianstructure = [jacobianstructure;depends];
end
Z = sparse(jacobianstructure);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function test01 ( input_filename )

%*****************************************************************************80
%
%% TEST01 tests LATINIZE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2006
%
%  Author:
%
%    John Burkardt
%

%
%  Get the row and column dimensions of the dataset.
%
  [ m, n ] = r8mat_header_read (  input_filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".\n', input_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension M = %d\n', m );
  fprintf ( 1, '  Number of points N  = %d\n', n );
%
%  Read the array.
%
  table = r8mat_data_read ( input_filename, m, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in "%s".\n', input_filename );

  r8mat_transpose_print_some ( m, n, table, 1, 1, 5, 5, ...
    '  Small portion of data read from file:' );
%
%  Latinize the array.
%
  table_latin = r8mat_latinize ( m, n, table );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Latinized the data.\n' );
%
%  Print out a small sample of the array.
%
  r8mat_transpose_print_some ( m, n, table_latin, 1, 1, 5, 5, ...
    '  Small portion of Latinized data:' );
%
%  Make up a name for the output file that is likely to be related
%  to the input file name, and unique.
%
  output_filename = file_name_ext_swap ( input_filename, 'latin.txt' );
%
%  Write the latinized array to a file.
%
  r8mat_write ( output_filename, m, n, table_latin );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Wrote the latinized data to "%s".\n', output_filename );

  return
end


text:
function E = material_boundary_facets(F,FI)
  % E = material_boundary_facets(F,FI)
  %
  % Inputs:
  %   F  #F by 3 list of triangle indices into rows of some V
  %   FI #F list of material ids
  % Outputs:
  %   E  #E by 2 list of unique edge boundaries between elements with different
  %     ids
  [uE,~,EMAP] = unique(sort([F(:,[2 3]);F(:,[3 1]);F(:,[1 2])],2),'rows');
  uE2FI = sparse(EMAP,repmat(FI,3,1),1);
  E = uE(sum(uE2FI,2)==1 | sum(uE2FI>0,2)==2,:);
end

author: alecjacobson

repo: gptoolbox


text:
function res = immedian(img, filtre, varargin)
%IMMEDIAN Compute median value in the neighboorhood of each pixel
%
%   RES = immedian(IMG, SE)
%   Compute the median filter of image IMG, using structuring element SE.
%   The goal of this function is to provide the same interface as for
%   other image filters (imopen, imerode ...), and to allow the use of
%   median filter with user-defined structuring element.
%   This function can be used for directional filtering.
%
%
%   RES = immedian(IMG, SE, PADOPT) also specify padding option. PADOPT can
%   be either 'zeros' or 'SYMMETRIC', see medfilt2 or ordfilt2 for details.
%
%   See Also: IMMEAN, IMDIRFILTER, MEDFILT2, ORDFILT2
%
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 16/02/2004.
%

%   HISTORY
%   17/02/2004: add support for 'strel' objects.
%   20/02/2004: add 'padopt' option, and documentation
%   2011-11-05 deprecate

warning('imael:deprecatedFunction', ...
    'function "immedian" has been deprecated and replaced by "imMedianFilter"');

% transform STREL object into single array
if strcmp(class(filtre), 'strel') %#ok<STISA>
    filtre = getnhood(filtre);
end

% get padopt option.
padopt = 'zeros'; % default for standard median filtering in matlab
if ~isempty(varargin)
    padopt = varargin{1};
end

% perform filtering
order = ceil(sum(filtre(:))/2);
res = ordfilt2(img, order, filtre, padopt);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function varargout = surface(varargin)
%SURFACE  Plot surface of a SPHEREFUN.
%   See SURF for a complete description of the various forms of the input.
%
% See also SPHEREFUN/SURF.

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

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

end


text:
function f=fst(x,a,b)

   f=((abs(x-a)-abs(x-b))/(b-a)+1)/2;

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
function att = attenuationWater(f, T)
%ATTENUATIONWATER Calculate ultrasound attenuation in distilled water.
%
% DESCRIPTION:
%       attenuationWater calculates the ultrasonic absorption in distilled
%       water at a given temperature and frequency using a 7th order
%       polynomial fitted to the data given by Pinkerton (1949).
%
% USAGE:
%       att = attenuationWater(f, T)
%
% INPUTS:
%       f   - array of frequency values [MHz]
%       T   - water temperature [degC]
%
% OUTPUTS:
%       att - attenuation [dB/cm]
%
% ABOUT:
%       author      - Bradley E. Treeby
%       date        - 10th November 2008
%       last udpate - 10th November 2008
%
% REFERENCES:
%   [1] Pinkerton (1949) "The Absorption of Ultrasonic Waves in Liquids and
%       its Relation to Molecular Constitution," Proceedings of the
%       Physical Society. Section B, 2, 129-141
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also speedSoundWater

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

if T < 0 || T > 60
    disp('WARNING: Temperature outside range of experimental data');
end

% conversion factor between Nepers and dB
NEPER2DB = 8.686;

% coefficients for 7th order polynomial fit
a_0 = 56.723531840522710;
a_1 = -2.899633796917384;
a_2 = 0.099253401567561;
a_3 = -0.002067402501557;
a_4 = 2.189417428917596e-005;
a_5 = -6.210860973978427e-008;
a_6 = -6.402634551821596e-010;
a_7 = 3.869387679459408e-012;

% compute attenuation
a_on_fsqr = (a_0 + a_1*T + a_2*T.^2 + a_3*T.^3 + a_4*T.^4 + a_5*T.^5 + a_6*T.^6 + a_7*T.^7)*1e-17;
att = NEPER2DB*1e12*f.^2*a_on_fsqr;


text:
function  dC=dervar3(C,n)

%DERVAR3 for core rotations
%
%function dC=dervar3(C,n)
%
%This function determines the derivative of the nth mode of
%the core rotation expression for the 3-way case w.r.t.
%maximization of the variance of the core
%
%This version has been optimized for speed, see below for a
%more easy to read scheme.

% Copyright (C) 1995-2006  Rasmus Bro & Claus Andersson
% Copenhagen University, DK-1958 Frederiksberg, Denmark, rb@life.ku.dk
%
% 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.

% $ Version 2.00 $ May 2001 $ Changed to array notation $ RB $ Not compiled $

W = size(C);
C = reshape(C,W(1),prod(W(2:end)));

mc = mean(mean(C.^2));
E = (C.^2 - mc).*C;
dC = zeros(W(n),W(n));
W1 = W(1);
W2 = W(2);
W3 = W(3);

if n==1,
   for a=1:W1,
      for b=1:W1,

         dc=0;
         for j=1:W3,
            tmp_0 = W2*(j-1);
            idxja = [1:W2] + tmp_0;
            dc = dc + sum(C(b,idxja).*E(a,idxja));
         end;

         dC(a,b) = dc;
      end;
   end;
end;

if n==2,
   for a=1:W2,
      for b=1:W2,

         dc=0;
         for j=1:W3,
            tmp_1 = W2*(j-1);
            idxja = a + tmp_1;
            idxjb = b + tmp_1;
            dc = dc + sum(C(:,idxjb).*E(:,idxja));
         end;

         dC(a,b) = dc;
      end;
   end;
end;

if n==3,
   for a=1:W3,
      tmp_2 = W2*(a-1);
      for b=1:W3,
         tmp_3 = W2*(b-1);

         dc=0;
         for j=1:W2,
            idxja = j + tmp_2;
            idxjb = j + tmp_3;
            dc = dc + sum(C(:,idxjb).*E(:,idxja));
         end;

         dC(a,b) = dc;
      end;
   end;
end;

%----------------------------------------------------------------------------------------
%function  dC=dervar3(C,W,n)
%
%%function dC=dervar3(C,W,n)
%
%%This function determines the derivative of the nth mode of
%%the core rotation expression for the 3-way case w.r.t.
%%maximization of the variance of the core
%
%mc=mean(mean(C.^2));
%dC=zeros(W(n),W(n));
%
%if n==1,
%   for a=1:W(1),
%      for b=1:W(1),
%
%         dc=0;
%         for i=1:W(2),
%            for j=1:W(3),
%               [idxia idxja]=getindxn(W,[a i j]);
%               [idxib idxjb]=getindxn(W,[b i j]);
%               dc = dc + (C(idxia,idxja)^2 - mc)*C(idxib,idxjb)*C(idxia,idxja);
%            end;
%         end;
%
%         dC(a,b) = dc;
%      end;
%   end;
%end;
%
%if n==2,
%   for a=1:W(2),
%      for b=1:W(2),
%
%         dc=0;
%         for i=1:W(1),
%            for j=1:W(3),
%               [idxia idxja]=getindxn(W,[i a j]);
%               [idxib idxjb]=getindxn(W,[i b j]);
%               dc = dc + (C(idxia,idxja)^2 - mc)*C(idxib,idxjb)*C(idxia,idxja);
%            end;
%         end;
%
%         dC(a,b) = dc;
%      end;
%   end;
%end;
%
%if n==3,
%   for a=1:W(3),
%      for b=1:W(3),
%
%         dc=0;
%         for i=1:W(1),
%            for j=1:W(2),
%               [idxia idxja]=getindxn(W,[i j a]);
%               [idxib idxjb]=getindxn(W,[i j b]);
%               dc = dc + (C(idxia,idxja)^2 - mc)*C(idxib,idxjb)*C(idxia,idxja);
%            end;
%         end;
%
%         dC(a,b) = dc;
%      end;
%   end;
%end;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [Au, idx ,idx2] = uniquecell(A)
    %function [Au, idx, idx2] = uniquecell(A)
    %For A a cell array of matrices (or vectors), returns
    %Au, which contains the unique matrices in A, idx, which contains
    %the indices of the last appearance of each such unique matrix, and
    %idx2, which contains th indices such that Au(idx2) == A
    %
    %Example usage:
    %
    %A = {[1,2,3],[0],[2,3,4],[2,3,1],[1,2,3],[0]};
    %[Au,idx,idx2] = uniquecell(A);
    %
    %Results in:
    %idx = [6,5,4,3]
    %Au  = {[0],[1,2,3],[2,3,1],[2,3,4]}
    %idx2 = [2,1,4,3,2,1]
    %
    %Algorithm: uses cellfun to translate numeric matrices into strings
    %           then calls unique on cell array of strings and reconstructs
    %           the initial matrices
    %
    %See also: unique
    B = cellfun(@(x) num2str(x(:)'),A,'UniformOutput',false);
    if nargout > 2
        [~,idx,idx2] = unique(B);
        Au = A(idx);
    else
        [~,idx] = unique(B);
        Au = A(idx);
    end
end


text:
function region = p24_region ( )

%*****************************************************************************80
%
%% P24_REGION returns the name of the integration region for problem 24.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string REGION, the name of the integration region.
%
  region = 'box';

  return
end

author: johannesgerer


text:
function [col4row,row4col,gain]=assign2DHungarian(C,maximize)
%%ASSIGN2DHUNGARIAN Solve the two-dimensional assignment problem with a
%          rectangular cost matrix C using the O(n^4) complexity Hungarian
%          algorithm. The problem being solved can be formulated as
%          minimize (or maximize)
%          \sum_{i=1}^{numRow}\sum_{j=1}^{numCol}C_{i,j}*x_{i,j}
%          subject to
%          \sum_{j=1}^{numCol}x_{i,j}<=1 for all i
%          \sum_{i=1}^{numRow}x_{i,j}=1 for all j
%          x_{i,j}=0 or 1.
%          Assuming that numCol<=numRow. If numCol>numRow, then the
%          inequality and inequality conditions are switched. The function
%          assign2D uses an O(n^3) shortest augmenting path modification of
%          the Hungarian algorithm and thus should be preferred over this
%          algorithm.
%
%INPUTS: C A numRowXnumCol cost matrix that does not contain any NaNs and
%          where the largest finite element minus the smallest element is a
%          finite quantity (does not overflow) when performing minimization
%          and where the smallest finite element minus the largest
%          element is finite when performing maximization. Forbidden
%          assignments can be given costs of +Inf for minimization and -Inf
%          for maximization.
% maximize If true, the minimization problem is transformed into a
%          maximization problem. The default if this parameter is omitted
%          or an empty matrix is passed is false.
%
%OUTPUTS: col4row A numRowX1 vector where the entry in each element is an
%                 assignment of the element in that row to a column. 0
%                 entries signify unassigned rows. If the problem is
%                 infeasible, this is an empty matrix.
%         row4col A numColX1 vector where the entry in each element is an
%                 assignment of the element in that column to a row. 0
%                 entries signify unassigned columns. If the problem is
%                 infeasible, this is an empty matrix.
%            gain The sum of the values of the assigned elements in C. If
%                 the problem is infeasible, this is -4
%
%This function implements the O(n^4) complexity Hungarian algorithm, which
%is the combination of algorithms 4.2 and 4.3 of Chapter 4.2 of [1]. The
%algorithm has been slightly modified to handle rectangular matrices and to
%transform the input matrix to have all positice elements before running
%the algorithm, so that the algorithm can work on general matrices.
%Additionally, the algorithm contains a floating-point comparison to zero,
%which has been replaced by a threshold to make the algorithm more robust
%when handling non-integer costs.
%
%See the comments to assign2D for examples. The function assign2D
%
%REFERENCES:
%[1] R. Burkard, M. Dell'Amico, and S. Martello, Assignment Problems.
%    Philadelphia: Society for Industrial and Applied Mathematics, 2009.
%
%October 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(maximize))
    maximize=false;
end

numU=size(C,1);
numV=size(C,2);

didTranspose=false;
if(numU>numV)
    C=C';

    temp=numU;
    numU=numV;
    numV=temp;

    didTranspose=true;
end

%The cost matrix must have all non-negative elements for the assignment
%algorithm to work. This forces all of the elements to be positive. The
%delta is added back in when computing the gain in the end.
if(maximize==true)
    CDelta=max(C(:));

    %If C is all negative, do not shift.
    if(CDelta<0)
        CDelta=0;
    end

    C=-C+CDelta;
else
    CDelta=min(C(:));

    %If C is all positive, do not shift.
    if(CDelta>0)
        CDelta=0;
    end

    C=C-CDelta;
end

UScanned=zeros(numU,1);%Allocate space
VLabeledUnscanned=zeros(numV,1);%Allocate space
pred=zeros(numV,1);%Allocate space

VLabeled=zeros(numV,1);%Allocate space.

unassignedU=1:numU;
numUnassignedU=numU;
%Quick indicator of which vertex is assigned.
assignedBool=false(numU,1);
numAssignedU=0;

%Allocate and initialize.
u=zeros(numU,1);
v=zeros(numV,1);
row4col=zeros(numV,1);
col4row=zeros(numU,1);

while(numAssignedU<numU)
    kIdx=1;
    k=unassignedU(kIdx);
    while(assignedBool(k)==false)
        %%%GET THE ALTERNATING PATH

        %Unlabeled vertices.
        VUnlabeled=1:numV;
        numUnlabeled=numV;
        numLabeledUnscanned=0;
        numScannedU=0;
        numLabeled=0;

        fail=false;
        sink=0;
        i=k;
        while(fail==false&&sink==0)
            numScannedU=numScannedU+1;
            UScanned(numScannedU)=i;
            jIdx=1;
            while(jIdx<=numUnlabeled)
                j=VUnlabeled(jIdx);
                cost=C(i,j)-u(i)-v(j);
                %Ad-hoc threshold for testing if zero.
                if(cost<=4*eps(C(i,j)))
                    pred(j)=i;
                    %Remove from the list of unlabeled vertices.
                    VUnlabeled(jIdx)=VUnlabeled(numUnlabeled);
                    numUnlabeled=numUnlabeled-1;

                    %Add to the list of labeled but unscanned vertices.
                    numLabeledUnscanned=numLabeledUnscanned+1;
                    VLabeledUnscanned(numLabeledUnscanned)=j;

                    %Also just add to the list of labeled vertices.
                    numLabeled=numLabeled+1;
                    VLabeled(numLabeled)=j;
                else
                    jIdx=jIdx+1;
                end
            end
            if(numLabeledUnscanned==0)
                fail=true;
            else
                %Scan the last labeled, unscanned vertex.
                j=VLabeledUnscanned(numLabeledUnscanned);
                numLabeledUnscanned=numLabeledUnscanned-1;

                if(row4col(j)==0)
                    sink=j;
                else
                    i=row4col(j);
                end
            end
        end

        %%UPDATE USING THE ALTERNATING PATH
        if(sink>0)
            %Increase the primal solution.
            assignedBool(k)=true;
            numAssignedU=numAssignedU+1;
            unassignedU(kIdx)=unassignedU(numUnassignedU);
            numUnassignedU=numUnassignedU-1;

            j=sink;
            while(1)
                i=pred(j);
                row4col(j)=i;
                h=col4row(i);
                col4row(i)=j;
                j=h;
                if(i==k)
                    break;
                end
            end
        else
            %Update the dual solution.
            deltaVal=Inf;
            for curU=1:numScannedU
                i=UScanned(curU);
                for curV=1:numUnlabeled
                    j=VUnlabeled(curV);
                    cost=C(i,j)-u(i)-v(j);

                    if(cost<deltaVal)
                        deltaVal=cost;
                    end
                end
            end

            %If the problem is infeasible.
            if(~isfinite(deltaVal))
                col4row=[];
                row4col=[];
                gain=-1;
                return;
            end

            for curU=1:numScannedU
                i=UScanned(curU);
                u(i)=u(i)+deltaVal;
            end

            for curV=1:numLabeled
                j=VLabeled(curV);
                v(j)=v(j)-deltaVal;
            end
        end
    end
end

if(nargout>2)
    gain=0;
    for curRow=1:numU
        gain=gain+C(curRow,col4row(curRow));
    end

    %Adjust the gain for the initial offset of the cost matrix.
    if(maximize)
        gain=-gain+CDelta*numU;
    else
        gain=gain+CDelta*numU;
    end
end

if(didTranspose)
    temp=col4row;
    col4row=row4col;
    row4col=temp;
end
end

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


text:
function [bodyX, bodyY, bodyZ, R_body_2_inertial] = computeInertialBodyAxesFromFrameAeroAngles(ut, rVect, vVect, bodyInfo, bankAng, angOfAttack, angOfSideslip, baseFrame)
    arguments
        ut(1,1) double
        rVect(3,1) double
        vVect(3,1) double
        bodyInfo(1,1) KSPTOT_BodyInfo
        bankAng(1,1) double
        angOfAttack(1,1) double
        angOfSideslip(1,1) double
        baseFrame(1,1) AbstractReferenceFrame
    end

    frame = bodyInfo.getBodyCenteredInertialFrame();
    ce = CartesianElementSet(ut, rVect, vVect, frame);
    ce = ce.convertToFrame(baseFrame, true);
    rVectFrame = ce.rVect;
    vVectFrame = ce.vVect;

    [R_wind_2_frame, ~, ~, ~] = computeWindFrame(rVectFrame, vVectFrame);
    R_body_2_wind = eul2rotmARH([angOfSideslip,angOfAttack,bankAng],'zyx');
    [~,~,~, R_frame_2_inertial] = baseFrame.getOffsetsWrtInertialOrigin(ut, ce);
    R_body_2_inertial = real(R_frame_2_inertial * R_wind_2_frame * R_body_2_wind);

    bodyX = R_body_2_inertial(:,1);
    bodyY = R_body_2_inertial(:,2);
    bodyZ = R_body_2_inertial(:,3);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function plot_z_coord(V,F)
%PLOT_Z_COORD  plot the z-coordinate on the input mesh
%
% plot_z_coord(V,F);
%
% Input:
%  V,F  mesh that is to be plotted

f = V(:,3);
tsurf(F,V, 'CData',f);
axis equal;
shading interp;

end

author: odedstein

repo: sgi-introduction-course


text:
function yv = dif_vals ( nd, xd, yd, nv, xv )

%*****************************************************************************80
%
%% DIF_VALS evaluates a divided difference polynomial at a set of points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carl deBoor,
%    A Practical Guide to Splines,
%    Springer, 2001,
%    ISBN: 0387953663,
%    LC: QA1.A647.v27.
%
%  Parameters:
%
%    Input, integer ND, the order of the difference table.
%
%    Input, real XD(ND), the X values of the difference table.
%
%    Input, real YD(ND), the divided differences.
%
%    Input, integer NV, the number of evaluation points.
%
%    Input, real XV(NV), the evaluation points.
%
%    Output, real YV(NV), the value of the divided difference
%    polynomial at the evaluation points.
%

%
%  Deal with the blasted curse of row vectors.
%
  yd = yd ( : );
  xd = xd ( : );
  xv = xv ( : );

  yv = zeros ( nv, 1 );

  yv(1:nv) = yd(nd);

  for i = 1 : nd - 1
    yv(1:nv) = yd(nd-i) + ( xv(1:nv) - xd(nd-i) ) .* yv(1:nv);
  end

  return
end

repo: jburkardt-m


text:
function [xVals, yVals, zVals] = getDeformXYZVals(deformStruct)
%"getDeformXYZVals"
%   Returns the x, y, and z values of the cols, rows, and slices of the
%   passed deformStruct.
%
%   REMINDER: These x,y,z values are the coordinates of the MIDDLE of the
%   voxels of the scan.  They are not the coordinates of the dividers
%   between the voxels.
%
%
% xVals yVals zVals : x,y,z Values for scan deformation grid.
%
% Usage:
%   function [xVals, yVals, zVals] = getDeformXYZVals(deformStruct)
%
% 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/>.

sizeDim1 = double(deformStruct.gridDimensions(2))-1;
sizeDim2 = double(deformStruct.gridDimensions(1))-1;
sizeDim3 = double(deformStruct.gridDimensions(3));
zStart = deformStruct.imagePositionPatient(3);
zRes = deformStruct.gridResolution(3);
deformStruct.xOffset = double(deformStruct.xOffset);
deformStruct.yOffset = double(deformStruct.yOffset);

xVals = deformStruct.xOffset - (sizeDim2*deformStruct.gridResolution(1))/2 : deformStruct.gridResolution(1) : deformStruct.xOffset + (sizeDim2*deformStruct.gridResolution(1))/2;
yVals = fliplr(deformStruct.yOffset - (sizeDim1*deformStruct.gridResolution(2))/2 : deformStruct.gridResolution(2) : deformStruct.yOffset + (sizeDim1*deformStruct.gridResolution(2))/2);
zVals = linspace(zStart, zStart + zRes * sizeDim3, sizeDim3);
zVals = -fliplr(zVals);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function u = u_hat ( x )

%*****************************************************************************80
%
%% U_HAT evaluates the target function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2011
%
%  Author:
%
%    Jeff Borggaard, John Burkardt, Catalin Trenchea, Clayton Webster
%
%  Parameters:
%
%    Input, real X(*), the evaluation points.
%
%    Output, real U(*), the function values.
%
  u = x .* ( 1 - x.^2 );

  return
end


text:
function [h,g,a,info] = wfilt_oddevena(N)
%WFILT_ODDEVENA  Kingsbury's symmetric even filters
%
%   Usage: [h,g,a] = wfilt_oddevena(N);
%
%   `[h,g,a]=wfilt_oddevena(N)` with $N \in {1}$ returns Kingsbury's
%   even filters.
%
%   Examples:
%   ---------
%   :::
%     figure(1);
%     wfiltinfo('ana:oddevena1');
%
%     figure(2);
%     wfiltinfo('syn:oddevena1');
%
%   References: king02

% AUTHOR: Zdenek Prusa

info.istight = 0;

a = [2;2];

switch(N)
 case 1
    % Example 1. from the reference. Symmetric near-orthogonal
    garr = [
             0           0
             0           0
             0          -0.0004645
             0           0.0013349
            -0.0058109   0.0022006
             0.0166977  -0.0130127
            -0.0000641   0.0015360
            -0.0834914   0.0869008
             0.0919537   0.0833552
             0.4807151  -0.4885957
             0.4807151   0.4885957
             0.0919537  -0.0833552
            -0.0834914  -0.0869008
            -0.0000641  -0.0015360
             0.0166977   0.0130127
            -0.0058109  -0.0022006
             0          -0.0013349
             0           0.0004645
             0           0
             0           0
    ];

    % This scaling is not in the reference paper, but it is here to be
    % consistent
    garr = garr*sqrt(2);
    %garr = setnorm(garr,'energy');

    offset = -10;

  otherwise
        error('%s: No such filters.',upper(mfilename));
end

    %garr = [garr(:,3:4),garr(:,1:2)];
    modrange = (-1).^((0:size(garr,1)-1) + offset+1).';
    modrange2 = (-1).^((0:size(garr,1)-1) + offset).';

    harr =       [garr(:,2).*modrange2,...
                  garr(:,1).*modrange,...
                  ];

% In the biorthogonal case, the filters do not get time reversed
garr = flipud(garr);

htmp=mat2cell(harr,size(harr,1),ones(1,size(harr,2)));
h = cellfun(@(hEl)struct('h',hEl,'offset',offset),htmp(1:2),...
                   'UniformOutput',0);

gtmp=mat2cell(garr,size(garr,1),ones(1,size(garr,2)));

g = cellfun(@(gEl)struct('h',gEl,'offset',offset),gtmp(1:2),...
                   'UniformOutput',0);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat