text:
function soITrans = findSoITransitions(initialState, maxSearchUT, soiSkipIds, massLoss, orbitDecay, celBodyData)
%findSoITransitions Summary of this function goes here
%   Detailed explanation goes here

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %findSoITransitions() returns a matrix of information that contains
    % data on upcoming SoI transitions for a given state.  SoI transition
    % information is given in the following format:
    %
    % [transType, transUT, fromSoI, toSoI, fromRx, fromRy, fromRz, fromVx, fromVy, fromVz, toRx, toRy, toRz, toVx, toVy, toVz]
    %
    % transType - type of SoI transition.  1 for upwards, -1 for downwards
    % transUT   - UT of the SoI transition
    % fromSoI   - ID number of the body transitioning out of
    % toSoI     - ID number of the body transitioning into
    %
    % fromRx    - The following three are the position components of the s/c prior to leaving, relative to the fromSoI body
    % fromRy
    % fromRz
    % fromVx    - The following three are the velocity components of the s/c prior to leaving, relative to the fromSoI body
    % fromVy
    % fromVz
    %
    % toRx    - The following three are the position components of the s/c after moving to new SoI, relative to the toSoI body
    % toRy
    % toRz
    % toVx    - The following three are the velocity components of the s/c after moving to new SoI, relative to the toSoI body
    % toVy
    % toVz
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    global num_SoI_search_revs use_selective_soi_search soi_search_tol num_soi_search_attempts_per_rev;

    if(isempty(use_selective_soi_search))
        use_selective_soi_search = 0;
    end

    if(isempty(num_SoI_search_revs))
        num_SoI_search_revs = 3;
    end

    if(isempty(soi_search_tol))
        soi_search_tol = 1E-12;
    end

    if(isempty(num_soi_search_attempts_per_rev))
        num_soi_search_attempts_per_rev = 1000;
    end

    ut = initialState(1);

    if(isempty(maxSearchUT))
        maxSearchUT = Inf;
    else
        if(maxSearchUT <= ut)
            maxSearchUT = Inf;
        end
    end

    bodyID = initialState(8);
    bodyInfo = getBodyInfoByNumber(bodyID, celBodyData);

    gmu = bodyInfo.gm;
    rVect = initialState(2:4)';
    vVect = initialState(5:7)';

    [sma, ecc, inc, raan, arg, truINI] = getKeplerFromState(rVect,vVect,gmu);
%     meanINI = computeMeanFromTrueAnom(truINI, ecc);
%     scBodyInfo = getBodyInfoStructFromOrbit([sma, ecc, inc, raan, arg, meanINI, ut]);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %There are two cases to consider: moving up the heirarchy and moving
    % down the heirarchy.  For example:
    % Up Heirarchy: Kerbin -> Sun
    % Down Heirarchy Kerbin -> Mun
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %This segment looks for SoI transitions upwards.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    upSoITrans = [];
    parentBodyInfo = bodyInfo.getParBodyInfo(celBodyData);
    if(~isempty(parentBodyInfo))
        soiRadius = getSOIRadius(bodyInfo, parentBodyInfo);
        if(ecc < 1)
            [rAp, ~] = computeApogeePerigee(sma, ecc);
        else
            rAp = Inf;
        end

        if(ecc >= 1 || rAp > soiRadius) %definite SoI transition upwards exists!
            if(abs(norm(rVect)-soiRadius)<0.001 && dot(rVect, vVect) > 0)
                truSoI = truINI;
                tempEventLog = initialState;
            else
                truSoI = computeTrueAFromRadiusEcc(soiRadius, sma, ecc);
                tempEventLog = ma_executeCoast_goto_tru(truSoI, initialState, -1, false, soiSkipIds, [], massLoss, orbitDecay, celBodyData);
            end

            upSoITransUT = tempEventLog(end,1);

            if(upSoITransUT <= maxSearchUT)
                [rVectUp, vVectUp] = convertRVVectOnUpwardsSoITransition(bodyInfo, celBodyData, upSoITransUT, tempEventLog(end,2:4), tempEventLog(end,5:7));

                upSoITrans = [1, upSoITransUT, bodyInfo.id, parentBodyInfo.id, ...
                              tempEventLog(end,2:7), ...
                              rVectUp', vVectUp'];
            else
                upSoITrans = [];
            end
        else
            upSoITrans = [];
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %This segment looks for SoI transitions downwards.  This is the
    % trickier of the two algorithms.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    parentBodyInfo = bodyInfo.getParBodyInfo(celBodyData);
    [childBodies] = getChildrenOfParentInfo(celBodyData, lower(bodyInfo.name));
    downSoITrans = [];

    searchableChildBodies = KSPTOT_BodyInfo.empty(1,0);
    childBodySoIRadii = [];
    for(child = childBodies)
        childBodyInfo = child{1};

        if(use_selective_soi_search && ismember(childBodyInfo.id,soiSkipIds))
            continue;
        else
            searchableChildBodies(end+1) = childBodyInfo; %#ok<AGROW>
            childBodySoIRadii(end+1) = getSOIRadius(childBodyInfo, bodyInfo);
        end
    end

    if(not(isempty(searchableChildBodies)))
        meanMotion = computeMeanMotion(sma, gmu);
        odefun = @(ut, mean) soiSearchOdeFun(ut, mean, meanMotion);

        [maxSearchUTComputed, maxStepSize] = getMaxSoISearchTime(ut, sma, ecc, truINI, bodyInfo, gmu, parentBodyInfo, num_SoI_search_revs, num_soi_search_attempts_per_rev);
        if(isempty(maxSearchUT) || not(isfinite(maxSearchUT)) || isnan(maxSearchUT))
            maxSearchUT = maxSearchUTComputed;
        end

        tspan = [ut, maxSearchUT];
        if(isempty(maxSearchUT))
            fprintf('%0.9f - %0.9f - %0.9f - %0.9f - %0.9f - %0.9f - %0.9f\n\n', ut, sma, ecc, truINI, gmu, num_SoI_search_revs, num_soi_search_attempts_per_rev);
        end

        meanIni = computeMeanFromTrueAnom(truINI, ecc);
        if(ecc < 1)
            meanIni = AngleZero2Pi(meanIni);
        end
        y0 = meanIni;

        odeEvtFcn = @(ut,mean) getSoITransitionOdeEvents(ut, mean, sma, ecc, inc, raan, arg, bodyInfo, gmu, searchableChildBodies, childBodySoIRadii, celBodyData);
        options = odeset('RelTol',soi_search_tol, 'AbsTol',soi_search_tol, 'Events',odeEvtFcn, 'MaxStep',maxStepSize, 'Refine',1, 'InitialStep',maxStepSize);

        [~,~,te,~,ie] = ode113(odefun,tspan,y0,options);

        if(not(isempty(ie)))
            [crossingUT, I] = min(te);
            childBodyInfo = searchableChildBodies(ie(I));

            tempEventLog = ma_executeCoast_goto_ut(crossingUT, initialState, -1, false, soiSkipIds, massLoss,  orbitDecay, celBodyData);

            [rVectDown, vVectDown] = convertRVVectOnDownwardsSoITransition(childBodyInfo, celBodyData, crossingUT, tempEventLog(end,2:4), tempEventLog(end,5:7));

            if(dot(rVectDown, vVectDown) < 0)
                downSoITrans(end+1,:) = [-1, crossingUT, bodyInfo.id, childBodyInfo.id, ...
                                         tempEventLog(end,2:7), ...
                                         rVectDown', vVectDown'];
            end
        end
    end

    soITrans = [upSoITrans;downSoITrans];
end

function [maxSearchUT, maxStep] = getMaxSoISearchTime(utINI, sma, ecc, truINI, bodyInfo, gmu, parentBodyInfo, num_SoI_search_revs, num_soi_search_attempts_per_rev)
    numStepsPer = num_soi_search_attempts_per_rev;

    if(ecc < 1)
        oPeriod = computePeriod(sma,gmu);

        maxSearchUT = utINI + num_SoI_search_revs .* oPeriod;

        maxStep = abs(maxSearchUT - utINI)/num_SoI_search_revs/numStepsPer;
    else
        soiRadius = getSOIRadius(bodyInfo, parentBodyInfo);
        maxHypTru = AngleZero2Pi(real(computeTrueAFromRadiusEcc(soiRadius, sma, ecc)));

        meanMotion = computeMeanMotion(sma, gmu);
        meanIni = computeMeanFromTrueAnom(truINI, ecc);
        maxMean = computeMeanFromTrueAnom(maxHypTru, ecc);

        maxSearchUT = utINI + (maxMean - meanIni)./meanMotion;

        maxStep = abs(maxSearchUT - utINI)./numStepsPer;
    end
end

function meandot = soiSearchOdeFun(~, ~, meanMotion)
    meandot = meanMotion;
end

function [value, isterminal, direction] = getSoITransitionOdeEvents(ut, mean, sma, ecc, inc, raan, arg, bodyInfo, gmu, childBodies, searchableChildBodies, celBodyData)
    value = [];
    isterminal = [];
    direction = [];

    tru = computeTrueAnomFromMean(mean, ecc);
    [rVect, vVect] = getStatefromKepler(sma, ecc, inc, raan, arg, tru, gmu, false);

    for(i=1:length(childBodies)) %#ok<*NO4LP>
        childBodyInfo = childBodies(i);

        dVect = getAbsPositBetweenSpacecraftAndBody(ut, rVect, bodyInfo, childBodyInfo, celBodyData);
        distToChild = norm(dVect);

        rSOI = searchableChildBodies(i);

        val = distToChild - rSOI;

        value(end+1) = val; %#ok<AGROW>
        direction(end+1) = -1; %#ok<AGROW>
        isterminal(end+1) = 1; %#ok<AGROW>
    end
end

%     for(child = childBodies)
%         childBodyInfo = child{1};
%
%         if(use_selective_soi_search && ismember(childBodyInfo.id,soiSkipIds))
%             continue;
%         end
%
%         soiRadius = getSOIRadius(childBodyInfo, bodyInfo);
%         if(~isempty(parentBodyInfo))
%             soiRadiusParent = getSOIRadius(bodyInfo, parentBodyInfo);
%         else
%             soiRadiusParent = Inf;
%         end
%
%         [rAp, rPe] = computeApogeePerigee(sma, ecc);
%         if(ecc > 1)
%             rAp = Inf;
%         end
%
%         [rApBody, rPeBody] = computeApogeePerigee(childBodyInfo.sma, childBodyInfo.ecc);
%
%         %%%%%%%%%%%
%         %These are the two cases in which an SoI transition downwards is
%         %impossible:
%         % 1) S/C apogee is less than the child body perigee minus the SoI
%         %    radius (we don't get up the child body).
%         % 2) S/C perigee is greater than the child body apogee plus the SoI
%         %    radius (we don't get down to the child body).
%         % 3) Position vectors at the asc/desc nodes must be within some
%         %    small number of SoI radii.
%         %%%%%%%%%%%
%         if(rAp < rPeBody-soiRadius)
%             continue;
%         end
%         if(rPe > rApBody+soiRadius)
%             continue;
%         end
%
%         if(strict_SoI_search)
%             sma1=sma;
%             ecc1=ecc;
%             inc1=inc;
%             raan1=raan;
%             arg1=arg;
%
%             sma2=childBodyInfo.sma;
%             ecc2=childBodyInfo.ecc;
%             inc2=deg2rad(childBodyInfo.inc);
%             raan2=deg2rad(childBodyInfo.raan);
%             arg2=deg2rad(childBodyInfo.arg);
%
%             [TF] = intersectCheck(sma1, ecc1, inc1, raan1, arg1,...
%                                   sma2, ecc2, inc2, raan2, arg2,...
%                                   gmu, soiRadius);
%
%             if(TF == false)
%                 continue;
%             end
%         end
%
%         soiSf = 0;
%
%         distScChildBody = norm(getAbsPositBetweenSpacecraftAndBody(ut, rVect, scBodyInfo, childBodyInfo, celBodyData));
%         if(distScChildBody < (soiRadius-2*soiSf)) %The spacecraft starts out inside the SoI for some reason!
%                 [rVectDown, vVectDown] = convertRVVectOnDownwardsSoITransition(childBodyInfo, celBodyData, ut, initialState(end,2:4), initialState(end,5:7));
%                 if(dot(rVectDown,vVectDown) < 0) %we're already headed down, so we can transition downwards.  But if we're headed up at the SoI boundary, don't transition!  Will cause SoI dithering.
%                     downSoITrans(end+1,:) = [-1, ut, bodyInfo.id, childBodyInfo.id, ...
%                                              initialState(end,2:7), ...
%                                              rVectDown', vVectDown']; %#ok<AGROW>
%                 end
%         end
%
%         n1 = computeMeanMotion(sma, gmu);
%         n2 = computeMeanMotion(childBodyInfo.sma, gmu);
%         maxSearchRadius = min(rApBody+soiRadius,rAp);
%         minSearchRadius = max(rPeBody-soiRadius,rPe);
%         if(ecc < 1)
%             sPeriod = computeSynodicPeriod(n1, n2);
%             oPeriod = computePeriod(sma,gmu);
%             if(sPeriod / oPeriod > 3)
%                 sPeriod = 3*computePeriod(sma,gmu);
%             end
%             numPeriods = ceil(sPeriod/oPeriod);
%             numPeriods = num_SoI_search_revs;
%
%             truMax1=real(computeTrueAFromRadiusEcc(maxSearchRadius, sma, ecc));
%             truMin1=real(computeTrueAFromRadiusEcc(minSearchRadius, sma, ecc));
%
%             truMax2 = 2*pi-truMin1;
%             truMin2 = 2*pi-truMax1;
%
%             %%%%%%%%%Set 1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             MA1 = computeMeanFromTrueAnom(truMin1, ecc);
%             MA2 = computeMeanFromTrueAnom(truMax1, ecc);
%
%             dt = (MA1 - meanINI)/n1;
%             if(dt < 0)
% %                 if(meanINI>=MA1 && meanINI<=MA2)
% %                     dt = (2*pi-meanINI+MA1)/n1;
% %                 else
% %                     dt = (2*pi-meanINI+MA1)/n1;
% %                 end
%                 dt = (2*pi-meanINI+MA1)/n1;
%             end
%
%             minT = ut+dt-100;
%             if(minT < ut)
%                 minT = ut;
%             end
%
%             dMA = MA2-MA1;
%             sPeriod = dMA/n1;
%             tArr = [];
%             for(i=1:numPeriods)
%                 t1 = minT + (i-1)*oPeriod;
%                 t2 = minT + (i-1)*oPeriod + sPeriod + 100;
%
%                 if(t1 <= maxSearchUT && t2 <= maxSearchUT)
%                     tArr(end+1) = t1;
%                     tArr(end+1) = t2;
%                 elseif(t1 <= maxSearchUT && t2 > maxSearchUT)
%                     tArr(end+1) = t1;
%                     tArr(end+1) = maxSearchUT; %we don't want to exceed maxSearchUT
%                 else
%                     continue; %both t1 and t2 are greater than maxSearchUT, don't add them
%                 end
%             end
%
%             if(truINI >= truMin1 && truINI <= truMax1)
%                 tArr(end+1) = ut;
%                 tArr(end+1) = minT;
%             end
%
%             %%%%%%%%%Set 2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             MA1 = computeMeanFromTrueAnom(truMin2, ecc);
%             MA2 = computeMeanFromTrueAnom(truMax2, ecc);
%
%             dt = (MA1 - meanINI)/n1;
%             if(dt < 0)
% %                 dt = 0;
%                 dt = (2*pi-meanINI+MA1)/n1;
%             end
%
%             minT = ut+dt-100;
%             if(minT < ut)
%                 minT = ut;
%             end
%
%             dMA = MA2-MA1;
%             sPeriod = dMA/n1;
%             for(i=1:numPeriods)
%                 t1 = minT + (i-1)*oPeriod;
%                 t2 = minT + (i-1)*oPeriod + sPeriod + 100;
%
%                 if(t1 <= maxSearchUT && t2 <= maxSearchUT)
%                     tArr(end+1) = t1;
%                     tArr(end+1) = t2;
%                 elseif(t1 <= maxSearchUT && t2 > maxSearchUT)
%                     tArr(end+1) = t1;
%                     tArr(end+1) = maxSearchUT; %we don't want to exceed maxSearchUT
%                 else
%                     continue; %both t1 and t2 are greater than maxSearchUT, don't add them
%                 end
%             end
%
%             if(truINI >= truMin2 && truINI <= truMax2)
%                 tArr(end+1) = ut;
%                 tArr(end+1) = minT;
%             end
%
%         else
%             iniHyTruMax1 = AngleZero2Pi(real(computeTrueAFromRadiusEcc(maxSearchRadius, sma, ecc)));
%             iniHyTruMin1 = AngleZero2Pi(real(computeTrueAFromRadiusEcc(minSearchRadius, sma, ecc)));
%
%             iniHyTruMax2 = -iniHyTruMin1;
%             iniHyTruMin2 = -iniHyTruMax1;
%
%             %%%%%%%%%Set 1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             MA1 = computeMeanFromTrueAnom(iniHyTruMin1, ecc);
%             MA2 = computeMeanFromTrueAnom(iniHyTruMax1, ecc);
%
%             dt = (MA1 - meanINI)/n1;
%             if(dt < 0)
%                 dt = 0;
%             end
%
%             dMA = MA2-MA1;
%             sPeriod = dMA/n1;
%
%             minT = ut+dt-100;
%             if(minT < ut)
%                 minT = ut;
%             end
%             tArr = linspace(minT, ut+dt+sPeriod+100, 2);
%
%             %%%%%%%%%Set 2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             MA1 = computeMeanFromTrueAnom(iniHyTruMin2, ecc);
%             MA2 = computeMeanFromTrueAnom(iniHyTruMax2, ecc);
%
%             dt = (MA1 - meanINI)/n1;
%             if(dt < 0)
%                 dt = 0;
%             end
%
%             dMA = MA2-MA1;
%             sPeriod = dMA/n1;
%
%             minT = ut+dt-100;
%             if(minT < ut)
%                 minT = ut;
%             end
%
%             tArrTmp = linspace(minT, ut+dt+sPeriod+100, 2);
%             tArr(end+1) = min(tArrTmp);
%             tArr(end+1) = max(tArrTmp);
%         end
%
%         %These two lines sort the SoI calculations so that they're in
%         %time-order, meaning that earlier SoI transitions will not be
%         %missed if there are multiple possible SoI transitions.
%         tArr = sortrows(reshape(tArr, [2,length(tArr)/2])', 1);
%         tArr = reshape(tArr',1,[]);
%
%         findChildSoITransFunc = @(ut) isSoICrossing(ut, scBodyInfo, childBodyInfo, bodyInfo, soiRadius-soiSf, true, celBodyData);
%         findChildSoITransFuncNoAbs = @(ut) isSoICrossing(ut, scBodyInfo, childBodyInfo, bodyInfo, soiRadius-soiSf, false, celBodyData);
%         findChildSoITransFuncRealSoI = @(ut) isSoICrossing(ut, scBodyInfo, childBodyInfo, bodyInfo, soiRadius-10/1000, false, celBodyData);
%
%         options = optimset('TolX',soi_search_tol);
%         tolX2 = soi_search_tol;
%         for(i=1:length(tArr)/2) %#ok<*NO4LP>
%             tInds = [2*(i-1)+1, 2*i];
%
%             if(tArr(tInds(1)) >= tArr(tInds(2)))
%                 continue;
%             end
%
%             %%%%%%%%%
%             %This block of code is designed to weed out extraneous calls to
%             %fminbnd by avoiding looking for SoI transitions where the
%             %child body in question is clearly off on the other side of the
%             %orbit.
%             %%%%%%%%%
%             [rVectSc1] = getStateAtTime(scBodyInfo, tArr(tInds(1)), gmu);
%             [rVectSc2] = getStateAtTime(scBodyInfo, tArr(tInds(2)), gmu);
%
%             [rVectC1] = getStateAtTime(childBodyInfo, tArr(tInds(1)), gmu);
%             [rVectC2] = getStateAtTime(childBodyInfo, tArr(tInds(2)), gmu);
%
%             distT1 = norm(rVectSc1-rVectC1)/childBodyInfo.sma;
%             distT2 = norm(rVectSc2-rVectC2)/childBodyInfo.sma;
%
%             if(strict_SoI_search)
%                 smaFrac = 0.5;
%             else
%                 smaFrac = 60.0; %used to be 0.9
%             end
%
%             if(min(distT1, distT2) > smaFrac)
%                 if(strict_SoI_search)
%                     fprintf('While propagating Mission Architect orbit, strict SoI search setting skipped over a possible SoI transition between UT = %f sec and UT = %f sec.\n', tArr(tInds(1)), tArr(tInds(2)))
%                 end
%                 continue;
%             end
%
%             try
%                 minUT = max(tArr(tInds(1)), initialState(1));
%                 maxUT = tArr(tInds(2));
%
%                 if(minUT > maxUT)
%                     continue;
%                 end
%
%                 if(num_soi_search_attempts_per_rev > 1)
%                     startPts = linspace(minUT,maxUT,num_soi_search_attempts_per_rev+1);
%                 else
%                     startPts = [minUT,maxUT];
%                 end
%                 crossingUTs = NaN(length(startPts)-1,1);
%                 minDistFromSOI2s = crossingUTs;
%                 for(j=1:length(startPts)-1)
%                     minUT_Opt = startPts(j);
%                     maxUT_Opt = startPts(j+1);
%
%                     [crossingUTs(j), ~,~,~] = fminbnd(findChildSoITransFunc,minUT_Opt,maxUT_Opt, options);
%                     minDistFromSOI2s(j) = findChildSoITransFuncRealSoI(crossingUTs(j));
%
%                     if(minDistFromSOI2s(j) < 1)
%                         break; %we found SoI transition, no need to keep going
%                     end
%                 end
%
%                 [minDistFromSOI2, minI] = min(minDistFromSOI2s);
%                 crossingUT = crossingUTs(minI);
%
%                 if(minDistFromSOI2 > 100)
%                     continue;
%                 end
%             catch ME %#ok<*NASGU>
%                 continue;
%             end
%
%             if(abs(minDistFromSOI2) > 0)
%                 minUT = crossingUT-1200;
%                 maxUT = crossingUT+1200;
%                 unrefinedCrossingUt = crossingUT;
%
%                 [crossingUT,minDistFromSOI2, exitFlag] = bisection(findChildSoITransFuncNoAbs, minUT, maxUT, 0.0, tolX2);
%                 if(minDistFromSOI2 > 0.01)
%                     continue;
%                 end
%
%                 if(minDistFromSOI2 < 0.01 && isnan(crossingUT))
%                     options = optimset('Display','off', 'TolX',tolX2);
%                     [crossingUT,minDistFromSOI2, exitFlag] = fzero(findChildSoITransFuncNoAbs, unrefinedCrossingUt, options);
%
%                     if(minDistFromSOI2 > 0.01)
%                         continue;
%                     end
%                 end
%             end
%
%             if(isempty(crossingUT) || isnan(crossingUT) || ~isreal(crossingUT) || ~isfinite(crossingUT))
%                 continue;
%             end
%
%             if(abs(crossingUT-initialState(1)) <= 0.1)
%                 crossingUT = initialState(1);
%             elseif(abs(crossingUT-initialState(1)) > 0.1)
%                 if(crossingUT < initialState(1))
%                     continue;
%                 end
%             end
%
%             if(abs(maxSearchUT - crossingUT) <= 0.1)
%                 crossingUT = maxSearchUT;
%             elseif(abs(maxSearchUT - crossingUT) > 0.1)
%                 if(crossingUT > maxSearchUT)
%                     continue;
%                 end
%             end
%
%             tempEventLog = ma_executeCoast_goto_ut(crossingUT, initialState, -1, false, soiSkipIds, massLoss,  orbitDecay, celBodyData);
%
%             [rVectDown, vVectDown] = convertRVVectOnDownwardsSoITransition(childBodyInfo, celBodyData, crossingUT, tempEventLog(end,2:4), tempEventLog(end,5:7));
%             if(dot(rVectDown, vVectDown) > 0) %we found the outgoing part of the SoI transition
%                 crossingUTTemp = crossingUT;
%
%                 [smaTd, eccTd, ~, ~, ~, truTd] = getKeplerFromState(rVectDown, vVectDown, childBodyInfo.gm);
%                 [meanTd] = computeMeanFromTrueAnom(truTd, eccTd);
%                 [meanMotionTd] = computeMeanMotion(smaTd, childBodyInfo.gm);
% %                 [~,vVectTd] = getStatefromKepler(smaTd, eccTd, incTd, raanTd, argTd, 0.0, childBodyInfo.gm);
%
%                 maxUT = crossingUT - meanTd/meanMotionTd+100;
%                 minUT = crossingUT - 6*meanTd/meanMotionTd;
%                 minUT = max(minUT,initialState(1));
%
%                 if(maxUT < initialState(1))
%                     continue;
%                 end
%
%                 try
%                     [crossingUT, ~,~,~] = fminbnd(findChildSoITransFunc,minUT,maxUT, options);
%                     minDistFromSOI2 = findChildSoITransFuncRealSoI(crossingUT);
%                     if(minDistFromSOI2 > 100)
%                         continue;
%                     end
%                 catch ME %#ok<*NASGU>
%                     continue;
%                 end
%
%                 if(minDistFromSOI2 > 0)
%                     minUT = crossingUT-1200;
%                     minUT = max(minUT,initialState(1));
%                     maxUT = crossingUT+1200;
%                     [crossingUT,minDistFromSOI2] = bisection(findChildSoITransFuncRealSoI, minUT, maxUT, 0.0, tolX2);
%                     if(minDistFromSOI2 > 0.01)
%                         continue;
%                     end
%                 end
%
%                 if(isempty(crossingUT) || isnan(crossingUT) || ~isreal(crossingUT) || ~isfinite(crossingUT))
%                     continue;
%                 end
%
%                 if(abs(crossingUT-initialState(1)) <= 1)
%                     crossingUT = initialState(1);
%                 elseif(abs(crossingUT-initialState(1)) > 1)
%                     if(crossingUT < initialState(1))
%                         continue;
%                     end
%                 end
%
%                 if(abs(maxSearchUT - crossingUT) <= 1)
%                     crossingUT = maxSearchUT;
%                 elseif(abs(maxSearchUT - crossingUT) > 1)
%                     if(crossingUT > maxSearchUT)
%                         continue;
%                     end
%                 end
%
%                 tempEventLog = ma_executeCoast_goto_ut(crossingUT, initialState, -1, false, soiSkipIds, massLoss, orbitDecay, celBodyData);
%                 [rVectDown, vVectDown] = convertRVVectOnDownwardsSoITransition(childBodyInfo, celBodyData, crossingUT, tempEventLog(end,2:4), tempEventLog(end,5:7));
%             end
% %             if(abs(norm(rVectDown) - soiRadius) >= 0.01)
% %                 continue; %something wierd happened, go around
% %             end
%             if(crossingUT <= maxSearchUT)
%                 downSoITrans(end+1,:) = [-1, crossingUT, bodyInfo.id, childBodyInfo.id, ...
%                                          tempEventLog(end,2:7), ...
%                                          rVectDown', vVectDown']; %#ok<AGROW>
%             end
%             break;
%         end
%     end
%
%     soITrans = [upSoITrans;downSoITrans];
% end
%

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function lambda = g_lambda(L, mu)
%G_LAMBDA Computes lambda
%   Input: L, mu
%   Output: lambda

q = mu/L;

if (1/q > 10^6)
    lambda = 10 * mu;
elseif (1/q > 10^4)
    lambda = mu;
else
    lambda = mu/10;
end

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function header = egi_writeRawHeader( fid, header )
% Write EGI RAW header structure to an open file.
%
% USAGE:  header = egi_writeRawHeader( fid, header )
%
% INPUT:
%    - fid    : Matlab file handle to the EGI RAW file (file must be already open)
%    - header : Header structure of the same file
% OUTPUT:
%    - epochLength : length (in number of time samples) of each epoch in the RAW file

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2008
% ----------------------------- Script History ---------------------------------
% FT  24-Jun-2008  Creation
% ------------------------------------------------------------------------------

% WRITE RAW HEADER
fwrite(fid, header.versionNumber,           'integer*4');
fwrite(fid, header.recordingTime.Year     , 'integer*2');
fwrite(fid, header.recordingTime.Month    , 'integer*2');
fwrite(fid, header.recordingTime.Day      , 'integer*2');
fwrite(fid, header.recordingTime.Hour     , 'integer*2');
fwrite(fid, header.recordingTime.Minute   , 'integer*2');
fwrite(fid, header.recordingTime.Second   , 'integer*2');
fwrite(fid, header.recordingTime.Millisec , 'integer*4');
fwrite(fid, header.samplingRate           , 'integer*2');
fwrite(fid, header.numChans               , 'integer*2');
fwrite(fid, header.boardGain              , 'integer*2');
fwrite(fid, header.numConvBits            , 'integer*2');
fwrite(fid, header.ampRange               , 'integer*2');
fwrite(fid, header.numSamples             , 'integer*4');
fwrite(fid, header.numEvents              , 'integer*2');

% WRITE EVENTS
if (header.numEvents ~= 0)  % File contains event info.
    for i = 1:header.numEvents
        fwrite(fid, header.eventCodes{i}, 'uchar');
    end
end

author: brainstorm-tools


text:
function f1 = p02_f1 ( x )

%*****************************************************************************80
%
%% P02_F1 evaluates the first derivative for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the value of the variable.
%
%    Output, real F1, the first derivative of the
%    objective function.
%
  f1 = 2.0 * x - exp ( -x );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ output_im ] = patch_sew(input_patches, output_sz, patch_stsz)
% function [ output_im ] = patch_sew(input_patches, output_sz, patch_stsz)
%   Combine evenly spaced patches into a single image, averaging at voxels
%   in more than one patch. Allows for the input of patches from multiple
%   output images (assuming each image has been divided into similarly
%   sized and spaced patches.
%
%   Note: patch_create and patch_sew differ from im2col and col2im by
%   allowing a custom stride, and were designed to output a stack of 2d
%   patches for direct input to a 2d neural network.
%
% Inputs:
%   input_patches   [nlin_patch, ncol_patch, npatch]  stack of input patches
%   output size     2x1   dimensions of output 2D image (or images)
%   patch_stsz      1x1   step size btwn patches (stride)
%
% Outputs:
%   output_im       [nlin, ncol, nims]   sewed together images
%
%  Melissa Haskell, University of Michigan, 2021-09-08

if nargin < 3, ir_usage, end

%% Find patch size and determine the indices of each patch using createpatches

patch_sz = size(input_patches,1);

% input 2D image to get number of patches for 2d image of that size
[~, patch_indices] = patch_create(zeros(output_sz), patch_sz, patch_stsz);

%% Calculate total number of images and initialize output

npatch_per_image = size(patch_indices,3);
npatch = size(input_patches,3);
nrow = output_sz(1);
ncol = output_sz(2);
nims = npatch / npatch_per_image;
output_im = zeros(nrow,ncol,nims);

%% Sew all the patches together
for ii = 1:nims
    psf_img = zeros(nrow, ncol);  % for tracking how many patches are at each voxel
    sum_img = zeros(nrow, ncol);  % for adding all the patch values together
    for jj = 1:npatch_per_image

        patch_indx = (ii-1)*npatch_per_image + jj;
        patch = reshape(input_patches(:,:,patch_indx),[patch_sz, patch_sz]);

        ind_2dim = patch_indices(:,:,jj);
        r1 = ind_2dim(1); c1 = ind_2dim(2);
        r2 = ind_2dim(3); c2 = ind_2dim(4);

        psf_img(r1:r2,c1:c2) = psf_img(r1:r2,c1:c2) + 1;
        sum_img(r1:r2,c1:c2) = sum_img(r1:r2,c1:c2) + patch;

    end
    % scale summation image by the psf image
    output_im(:,:,ii) = sum_img ./ (psf_img + 1e-12);
end

end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
classdef MimDatasetStackItem < handle
    % MimDatasetStackItem. Part of the internal framework of the TD MIM Toolkit.
    %
    %     You should not use this class within your own code. It is intended to
    %     be used internally within the framework of the TD MIM Toolkit.
    %
    %     Used to store the dependency information of a plugin as it is being
    %     build up during execution. MimDatasetStackItem are created
    %     temporarily by the class MimDatasetStack and used to build up the
    %     dependency list. See MimDatasetStack for more information.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (SetAccess = private)
        InstanceIdentifier     % A unique identifier for this plugin result
        IgnoreDependencyChecks % Certain types of plugin do not cached their results and so are exempt from dependency checking
        DependencyList         % Current list of plugin results this plugin depends on
        IsEdited
        IsManual
        IsMarkerSet
        Schema                 % The disk cache version
    end

    properties (SetAccess = private, Transient = true)
        ExecutionTimer         % Times the execution time for this plugin
    end

    properties (Access = private, Transient)
        DiskCacheSchema = '0.1'
    end

    methods
        function obj = MimDatasetStackItem(instance_id, dependency_list, ignore_dependency_checks, start_timer, reporting)
            if nargin > 0
                if ~isa(instance_id, 'PTKDependency')
                    reporting.Error('MimDatasetStackItem:BadInstanceID', 'instance_id must be a PTKDependency object');
                end

                if ~isa(dependency_list, 'PTKDependencyList')
                    reporting.Error('MimDatasetStackItem:BadDependencyList', 'dependency_list must be a PTKDependencyList object');
                end

                obj.InstanceIdentifier = instance_id;
                obj.DependencyList = dependency_list;
                obj.IgnoreDependencyChecks = ignore_dependency_checks;
                obj.Schema = obj.DiskCacheSchema;
                obj.IsEdited = false;
                obj.IsManual = false;
                obj.IsMarkerSet = false;

                if start_timer
                    obj.ExecutionTimer = MimTimer(reporting);
                    obj.ExecutionTimer.Start;
                end
            end
        end

        function AddDependencies(obj, dependencies, reporting)
            % Adds more plugin result which this particular plugin result depends on.

            obj.DependencyList.AddDependenciesList(dependencies, reporting);
        end

        function StopAndDeleteTimer(obj)
            if ~isempty(obj.ExecutionTimer)
                obj.ExecutionTimer.Stop;
                obj.InstanceIdentifier.Attributes.ExecutionTime = obj.ExecutionTimer.TotalElapsedTime;
                obj.InstanceIdentifier.Attributes.SelfTime = obj.ExecutionTimer.SelfTime;
                obj.ExecutionTimer = [];
            end
        end

        function MarkEdited(obj)
            obj.IsEdited = true;
        end

        function MarkManual(obj)
            obj.IsManual = true;
        end

        function MarkMarkerSet(obj)
            obj.IsMarkerSet = true;
        end
    end

    methods (Static)
        function obj = loadobj(a)
            % This method is called when the object is loaded from disk.

            if isa(a, 'MimDatasetStackItem')
                obj = a;
            else
                % In the case of a load error, loadobj() gives a struct
                obj = MimDatasetStackItem;
                for field = fieldnames(a)'
                    if isprop(obj, field{1})
                        mp = findprop(obj, (field{1}));
                        if (~mp.Constant) && (~mp.Dependent) && (~mp.Abstract)
                            obj.(field{1}) = a.(field{1});
                        end
                    end
                end
            end

        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [A,B,D,E,I]=bpf1(X)
% Multiple Feedback Bandpass Filter
% Note that complex frequency variable s is not used here.
% These arrays are real, not complex.
%
R1=X(1);R2=X(2);R3=X(3);C1=X(4);C2=X(5);
%
N=2;M=1;U=1;Y=1;
%
A1=[0 1 1;
   1/R3 1 0;
   1 0 0];
%
B2=[-(1/R1+1/R2) 0 1/R1;
   0 0 0;
   1 -1 0];
%
P=diag([C1 C2]);
%
% The following statements never change from circuit to circuit.
%
V=A1\B2;H=V(U+1:U+N,1:N+M);I=eye(N);
AB=P\H;A=AB(1:N,1:N);B=AB(1:N,N+1:N+M);
D=V(Y:Y,1:N);E=V(Y:Y,N+1:N+M);

author: Sable

repo: mcbench-benchmarks


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_func10N(Temp,Event,handles)
handles = guidata(handles.MU_matrix_display);

MergeM=get(handles.Matrix_name_edit,'String');
if isfloat(handles.TMatrix(1))
    set(handles.Matrix_name_edit,'String',['10.^[' MergeM ']']);
else
    set(handles.Matrix_name_edit,'String',['10.^double([' MergeM '])']);
end
MU_calc_matrix(handles);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
classdef IMMOEA_F7 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing IM-MOEA

%------------------------------- Reference --------------------------------
% R. Cheng, Y. Jin, K. Narukawa, and B. Sendhoff, A multiobjective
% evolutionary algorithm using Gaussian process-based inverse modeling,
% IEEE Transactions on Evolutionary Computation, 2015, 19(6): 838-856.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            t = X(:,2:obj.D).^(1./(1+3*repmat(2:obj.D,size(X,1),1)/obj.D)) - repmat(X(:,1),1,obj.D-1);
            g = 1 + 9*mean(t.^2,2);
            PopObj(:,1) = 1 - exp(-4*X(:,1)).*sin(6*pi*X(:,1)).^6;
            PopObj(:,2) = g.*(1-(PopObj(:,1)./g).^2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            minf1  = min(1-exp(-4*(0:1e-6:1)).*(sin(6*pi*(0:1e-6:1))).^6);
            R(:,1) = linspace(minf1,1,N)';
            R(:,2) = 1 - R(:,1).^2;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [model] = SpectralHashing(X, maxbits)
%
% Input
%   X = features matrix [Nsamples, Nfeatures]
%   maxbits = number of bits (nbits do not need to be a multiple of 8)
%
%
% Spectral Hashing
% Y. Weiss, A. Torralba, R. Fergus.
% Advances in Neural Information Processing Systems, 2008.

% 1) PCA

C = cov(X);
sizeC = size(C,1);
npca = min(maxbits, sizeC);
if maxbits > sizeC
    maxbits = sizeC;
end

if npca > sizeC/2
    [pc, eigvalue] = eig(C);
    eigvalue = diag(eigvalue);
    [~, index] = sort(-eigvalue);
    index=index(1:npca);
    pc = pc(:,index);
else
    [pc, ~] = eigs(C, npca);
end

% The above line can be replaced as it is very slow in comparison to
% the PCA code by Deng Cai
% [pc, l] = PCA(X./sqrt(Nsamples-1),struct('ReducedDim',npca));

X = X * pc; % no need to remove the mean

% 2) fit uniform distribution

mn = min(X)-eps;
mx = max(X)+eps;

% 3) enumerate eigenfunctions

R=(mx-mn);
maxMode = ceil((maxbits+1)*R/max(R));

nModes = sum(maxMode)-length(maxMode)+1;
modes = ones([nModes npca]);
m = 1;
for i=1:npca
    modes(m+1:m+maxMode(i)-1,i) = 2:maxMode(i);
    m = m+maxMode(i)-1;
end
modes = modes - 1;
omega0 = pi./R;
omegas = modes.*repmat(omega0, [nModes 1]);
eigVal = -sum(omegas.^2,2);
[yy,ii]= sort(-eigVal);
modes=modes(ii(2:maxbits+1),:);

% 4) store model paramaters

model.pc = pc;
model.mn = mn;
model.mx = mx;
% model.mx = mx;
model.modes = modes;

end

author: ZJULearning

repo: MatlabFunc


text:
% QUATERNION MULTIPLICATION
function [qv] = qMultiply(q,v)
% Calculate the product of two quaternions
% Associated block:
% "Quaternion Multiplication"
% Multiply the quaternion elements

assert(size(q,1) == 4 && size(v,1) == 4,...
    'Both quaternion must be provided as 4x1 column vectors')
% Quaternion projection matrix
qv = [v(1), -v(2), -v(3), -v(4);
      v(2),  v(1), -v(4),  v(3);
      v(3),  v(4),  v(1), -v(2);
      v(4), -v(3),  v(2),  v(1)]*q; % Confirmed with matlab website
% Re-normalise the output
qv = unit(qv);
end

repo: OpenMAS


text:
%TROTX Rotation about X axis
%
% T = TROTX(THETA) is a homogeneous transformation (4x4) representing a rotation
% of THETA radians about the x-axis.
%
% T = TROTX(THETA, 'deg') as above but THETA is in degrees.
%
% Notes::
% - Translational component is zero.
%
% See also ROTX, TROTY, TROTZ, TROT2.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function T = trotx(t, varargin)
	T = [rotx(t, varargin{:}) [0 0 0]'; 0 0 0 1];

author: Allopart

repo: rbpf-gmapping


text:
function byteStream = mrMeshTcpSerializeParams(params)
%
%    byteStream = mrMeshTcpSerializeParams(params)
%
% Turns the Matlab structure params into a byte-stream suitable for
% TCP transmission using the mrMesh protocol.  That protocol uses pnet.
%
% See also:  mrMesh, pnet_WriteMsg and pnet_ReadMsg
%
% HISTORY
% 2007.04.12 RFD wrote it.

if(~isstruct(params))
    byteStream=params;
    return;
end

byteStream = [];
if(~isempty(params))
    % the termination char is a new-line
    termSeq = uint8(sprintf('\n'));
    eq = uint8(' = ');
    % arrays (non-scalars) are sent in raw for and are enclodes in
    % single-quotes. Scalars are sent in text (sprintf) form.
    sq = uint8('''');
    fnStr = fieldnames(params);

    for(ii=1:length(fnStr))
        data = getfield(params,fnStr{ii});
        sz = size(data);
        data = data(:)';
        fn = uint8(fnStr{ii});
        if(ischar(data))
            fn = uint8(fnStr{ii});
            byteStream = [byteStream fn eq sq data sq termSeq];

        elseif(length(data)>1)
            % *** FIX ME ***
            % The current mrMesh server expects all numeric arrays to be
            % sent as doubles. This sends 8x more data than necessary when
            % sending uint8's! We should fix the server so that it can
            % accept various data types.
            dimStr = sprintf('%d,',sz);
            fn = uint8(sprintf('%s[%s]',fnStr{ii},dimStr(1:end-1)));
            byteStream = [byteStream fn eq sq typecast(double(data),'uint8') sq termSeq];
        else
            fn = uint8(fnStr{ii});
            if(isinteger(data))
                byteStream = [byteStream fn eq uint8(sprintf('%d',data)) termSeq];
            else
                byteStream = [byteStream fn eq uint8(sprintf('%f',data)) termSeq];
            end
        end
    end
end

byteStream(end+1) = 0;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% dataset      -> (1=>user data) or (2=>toy example)
% type         -> (1=> Regression model) or (2=>Classification model)
% num_glevel   -> number of hidden nodes in the net (gating levels)
% num_exp      -> number of experts in the net
% branch_fact  -> dimension of the hidden nodes in the net
% cov_dim      -> root node dimension
% res_dim      -> output node dimension
% nodes_info   -> 4 x num_glevel+2 matrix that contain all the info about the nodes:
%                 nodes_info(1,:) = nodes type: (0=>gaussian)or(1=>softmax)or(2=>mlp)
%                 nodes_info(2,:) = nodes size: [cov_dim   num_glevel x branch_fact   res_dim]
%                 nodes_info(3,:) = hidden units number (for mlp nodes)
%                                  |- optimizer iteration number (for softmax & mlp CPD)
%                 nodes_info(4,:) =|- covariance type (for gaussian CPD)->
%                                  | (1=>Full)or(2=>Diagonal)or(3=>Full&Tied)or(4=>Diagonal&Tied)
% fh1 -> Figure: data & decizion boundaries; fh2 -> confusion matrix; fh3 -> LL trace
% test_data    -> test data matrix
% train_data   -> training data matrix
% ntrain       -> size(train_data,2)
% ntest        -> size(test_data,2)
% cases        -> (cell array) training data formatted for the learning engine
% bnet         -> bayesian net before learning
% bnet2        -> bayesian net after learning
% ll           -> log-likelihood before learning
% LL2          -> log-likelihood trace
% onodes       -> obs nodes in bnet & bnet2
% max_em_iter  -> maximum number of interations of the EM algorithm
% train_result -> prediction on the training set (as test_result)
%
% IMPORTANT: CHECK the loading path (lines 64 & 364)
% ----------------------------------------------------------------------------------------------------
% -> pierpaolo_b@hotmail.com   or   -> pampo@interfree.it
% ----------------------------------------------------------------------------------------------------

error('this no longer works with the latest version of BNT')

clear all;
clc;
disp('---------------------------------------------------');
disp('  Hierarchical Mixtures of Experts models builder ');
disp('---------------------------------------------------');
disp(' ')
disp('   Using this script you can build both an HME model')
disp('as in [Wat94] and [Jor94] i.e. with ''softmax'' gating')
disp('nodes and ''gaussian'' ( for regression ) or ''softmax''')
disp('( for classification ) expert node, and its variants')
disp('called ''gated nets'' where we use ''mlp'' models in')
disp('place of a number of ''softmax'' ones [Mor98], [Wei95].')
disp('  You can decide to train and test the model on your')
disp('datasets  or  to evaluate its  performance on  a toy')
disp('example.')
disp(' ')
disp('Reference')
disp('[Mor98] P. Moerland (1998):')
disp('        Localized mixtures of experts. (http://www.idiap.ch/~perry/)')
disp('[Jor94] M.I. Jordan, R.A. Jacobs (1994):')
disp('        HME and the EM algorithm. (http://www.cs.berkeley.edu/~jordan/)')
disp('[Wat94] S.R. Waterhouse, A.J. Robinson (1994):')
disp('        Classification using HME. (http://www.oigeeza.com/steve/)')
disp('[Wei95] A.S. Weigend, M. Mangeas (1995):')
disp('        Nonlinear gated experts for time series.')
disp(' ')

if 0
disp('(See the figure)')
pause(5);
%%%%%WARNING!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
im_path=which('HMEforMatlab.jpg');
fig=imread(im_path, 'jpg');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure('Units','pixels','MenuBar','none','NumberTitle','off', 'Name', 'HME model');
image(fig);
axis image;
axis off;
clear fig;
set(gca,'Position',[0 0 1 1])
disp('(Press any key to continue)')
pause
end

clc
disp('---------------------------------------------------');
disp('              Specify the Architecture             ');
disp('---------------------------------------------------');
disp(' ');
disp('What kind of model do you need?')
disp(' ')
disp('1) Regression ')
disp('2) Classification')
disp(' ')
type=input('1 or 2?: ');
if (isempty(type)|(~ismember(type,[1 2]))), error('Invalid value'); end
clc
disp('----------------------------------------------------');
disp('               Specify the Architecture             ');
disp('----------------------------------------------------');
disp(' ')
disp('Now you have to set the number of experts and gating')
disp('levels in the net.  This script builds only balanced')
disp('hierarchy with the same branching factor (>1)at each')
disp('(gating) level. So remember that: ')
disp(' ')
disp('         num_exp = branch_fact^num_glevel           ')
disp(' ')
disp('with branch_fact >=2.')
disp('You can also set to zeros the number of gating level')
disp('in order to obtain a classical GLM model.           ')
disp(' ')
disp('----------------------------------------------------');
disp(' ')
num_glevel=input('Insert the number of gating levels {0,...,20}: ');
if (isempty(num_glevel)|(~ismember(num_glevel,[0:20]))), error('Invalid value'); end
nodes_info=zeros(4,num_glevel+2);
if num_glevel>0, %------------------------------------------------------------------------------------
    for i=2:num_glevel+1,
        clc
        disp('----------------------------------------------------');
        disp('               Specify the Architecture             ');
        disp('----------------------------------------------------');
        disp(' ')
        disp(['-> Gating network ', num2str(i-1), ' is a: '])
        disp(' ')
        disp('   1) Softmax model');
        disp('   2) Two layer perceptron model')
        disp(' ')
        nodes_info(1,i)=input('1 or 2?: ');
        if (isempty(nodes_info(1,i))|(~ismember(nodes_info(1,i),[1 2]))), error('Invalid value'); end
        disp(' ')
        if nodes_info(1,i)==2,
           nodes_info(3,i)=input('Insert the number of units in the hidden layer: ');
           if (isempty(nodes_info(3,i))|(floor(nodes_info(3,i))~=nodes_info(3,i))|(nodes_info(3,i)<=0)),
              error(['Invalid value: ', num2str(nodes_info(3,i)), ' is not a positive integer!']);
           end
           disp(' ')
        end
        nodes_info(4,i)=input('Insert the optimizer iteration number: ');
        if (isempty(nodes_info(4,i))|(floor(nodes_info(4,i))~=nodes_info(4,i))|(nodes_info(4,i)<=0)),
           error(['Invalid value: ', num2str(nodes_info(4,i)), ' is not a positive integer!']);
        end
    end
    clc
    disp('---------------------------------------------------------');
    disp('                 Specify the Architecture                ');
    disp('---------------------------------------------------------');
    disp(' ')
    disp('Now you have to set the number  of experts in the network');
    disp('The value will be adjusted in order to obtain a hierarchy');
    disp('as said above.')
    disp(' ');
    num_exp=input(['Insert the approximative number of experts (>=', num2str(2^num_glevel), '): ']);
    if (isempty(num_exp)|(num_exp<=0)|(num_exp<2^num_glevel)),
        error('Invalid value');
    end
    app1=0; base=2;
    while app1<num_exp,
        app1=base^num_glevel;
        base=base+1;
    end
    app2=(base-2)^num_glevel;
    branch_fact=base-1;
    if app2>=(2^num_glevel)&(abs(app2-num_exp)<abs(app1-num_exp)),
        branch_fact=base-2;
    end
    clear app1 app2 base;
    disp(' ')
    disp(['The effective number of experts in the net is: ', num2str(branch_fact^num_glevel), '.'])
    disp(' ');
else
    clc
    disp('---------------------------------------------------------');
    disp('        Specify the Architecture (GLM model)             ');
    disp('---------------------------------------------------------');
    disp(' ')
end % END of: if num_glevel>0-------------------------------------------------------------------------

if type==2,
    disp(['-> Expert node is a: '])
    disp(' ')
    disp('   1) Softmax model');
    disp('   2) Two layer perceptron model')
    disp(' ')
    nodes_info(1,end)=input('1 or 2?: ');
    if (isempty(nodes_info(1,end))|(~ismember(nodes_info(1,end),[1 2]))),
        error('Invalid value');
    end
    disp(' ')
    if nodes_info(1,end)==2,
       nodes_info(3,end)=input('Insert the number of units in the hidden layer: ');
       if (isempty(nodes_info(3,end))|(floor(nodes_info(3,end))~=nodes_info(3,end))|(nodes_info(3,end)<=0)),
           error(['Invalid value: ', num2str(nodes_info(3,end)), ' is not a positive integer!']);
       end
       disp(' ')
    end
    nodes_info(4,end)=input('Insert the optimizer iteration number: ');
    if (isempty(nodes_info(4,end))|(floor(nodes_info(4,end))~=nodes_info(4,end))|(nodes_info(4,end)<=0)),
        error(['Invalid value: ', num2str(nodes_info(4,end)), ' is not a positive integer!']);
    end
elseif type==1,
    disp('What kind of covariance matrix structure do you want?')
    disp(' ')
    disp('   1) Full');
    disp('   2) Diagonal')
    disp('   3) Full & Tied');
    disp('   4) Diagonal & Tied')

    disp(' ')
    nodes_info(4,end)=input('1, 2, 3 or 4?: ');
    if (isempty(nodes_info(4,end))|(~ismember(nodes_info(4,end),[1 2 3 4]))),
        error('Invalid value');
    end
end
clc
disp('----------------------------------------------------');
disp('                    Specify the Input               ');
disp('----------------------------------------------------');
disp(' ')
disp('Do you want to...')
disp(' ')
disp('1) ...use your own dataset?')
disp('2) ...apply the model on a toy example?')
disp(' ')
dataset=input('1 or 2?: ');
if (isempty(dataset)|(~ismember(dataset,[1 2]))), error('Invalid value'); end
if dataset==1,
    if type==1,
        clc
        disp('-------------------------------------------------------');
        disp('        Specify the Input - Regression problem         ');
        disp('-------------------------------------------------------');
        disp(' ')
        disp('Be sure that each row of your data matrix is an example');
        disp('with the covariate values that precede the respond ones')
        disp(' ')
        disp('-------------------------------------------------------');
        disp(' ')
        cov_dim=input('Insert the covariate space dimension: ');
        if (isempty(cov_dim)|(floor(cov_dim)~=cov_dim)|(cov_dim<=0)),
          error(['Invalid value: ', num2str(cov_dim), ' is not a positive integer!']);
        end
        disp(' ')
        res_dim=input('Insert the dimension of the respond variable: ');
        if (isempty(res_dim)|(floor(res_dim)~=res_dim)|(res_dim<=0)),
            error(['Invalid value: ', num2str(res_dim), ' is not a positive integer!']);
        end
        disp(' ');
    elseif type==2
        clc
        disp('-------------------------------------------------------');
        disp('      Specify the Input - Classification problem       ');
        disp('-------------------------------------------------------');
        disp(' ')
        disp('Be sure that each row of your data matrix is an example');
        disp('with the covariate values that precede the class labels');
        disp('(integer value >=1).                                   ');
        disp(' ')
        disp('-------------------------------------------------------');
        disp(' ')
        cov_dim=input('Insert the covariate space dimension: ');
        if (isempty(cov_dim)|(floor(cov_dim)~=cov_dim)|(cov_dim<=0)),
          error(['Invalid value: ', num2str(cov_dim), ' is not a positive integer!']);
        end
        disp(' ')
        res_dim=input('Insert the number of classes: ');
        if (isempty(res_dim)|(floor(res_dim)~=res_dim)|(res_dim<=0)),
          error(['Invalid value: ', num2str(res_dim), ' is not a positive integer!']);
        end
        disp(' ')
    end
    % ------------------------------------------------------------------------------------------------
    % Loading training data --------------------------------------------------------------------------
    % ------------------------------------------------------------------------------------------------
    train_path=input('Insert the complete (with extension) path of the training data file:\n >> ','s');
    if isempty(train_path), error('You must specify a data set for training!'); end
    if ~isempty(findstr('.mat',train_path)),
        ap=load(train_path); app=fieldnames(ap); train_data=eval(['ap.', app{1,1}]);
        clear ap app;
    elseif ~isempty(findstr('.txt',train_path)),
        train_data=load(train_path, '-ascii');
    else
        error('Invalid data format: not a .mat or a .txt file')
    end
    if (size(train_data,2)~=cov_dim+res_dim)&(type==1),
        error(['Invalid data matrix size: ', num2str(size(train_data,2)), ' columns rather than ',...
            num2str(cov_dim+res_dim),'!']);
    elseif (size(train_data,2)~=cov_dim+1)&(type==2),
        error(['Invalid data matrix size: ', num2str(size(train_data,2)), ' columns rather than ',...
            num2str(cov_dim+1),'!']);
    elseif (~isempty(find(ismember(intersect([train_data(:,end)' 1:res_dim],...
            train_data(:,end)'),[1:res_dim])==0)))&(type==2),
        error('Invalid class label');
    end
    ntrain=size(train_data,1);
    train_d=train_data(:,1:cov_dim);
    if type==2,
        train_t=zeros(ntrain, res_dim);
        for m=1:res_dim,
            train_t((find(train_data(:,end)==m))',m)=1;
        end
    else
        train_t=train_data(:,cov_dim+1:end);
    end
    disp(' ')
    % ------------------------------------------------------------------------------------------------
    % Loading test data ------------------------------------------------------------------------------
    % ------------------------------------------------------------------------------------------------
    disp('(If you don''t want to specify a test-set press ''return'' only)');
    test_path=input('Insert the complete (with extension) path of the test data file:\n >> ','s');
    if ~isempty(test_path),
        if ~isempty(findstr('.mat',test_path)),
            ap=load(test_path); app=fieldnames(ap); test_data=eval(['ap.', app{1,1}]);
            clear ap app;
        elseif ~isempty(findstr('.txt',test_path)),
            test_data=load(test_path, '-ascii');
        else
            error('Invalid data format: not a .mat or a .txt file')
        end
        if (size(test_data,2)~=cov_dim)&(size(test_data,2)~=cov_dim+res_dim)&(type==1),
            error(['Invalid data matrix size: ', num2str(size(test_data,2)), ' columns rather than ',...
                num2str(cov_dim+res_dim), ' or ', num2str(cov_dim), '!']);
        elseif (size(test_data,2)~=cov_dim)&(size(test_data,2)~=cov_dim+1)&(type==2),
            error(['Invalid data matrix size: ', num2str(size(test_data,2)), ' columns rather than ',...
                num2str(cov_dim+1), ' or ', num2str(cov_dim), '!']);
        elseif (~isempty(find(ismember(intersect([test_data(:,end)' 1:res_dim],...
                test_data(:,end)'),[1:res_dim])==0)))&(type==2)&(size(test_data,2)==cov_dim+1),
            error('Invalid class label');
        end
        ntest=size(test_data,1);
        test_d=test_data(:,1:cov_dim);
        if (type==2)&(size(test_data,2)>cov_dim),
            test_t=zeros(ntest, res_dim);
            for m=1:res_dim,
                test_t((find(test_data(:,end)==m))',m)=1;
            end
        elseif (type==1)&(size(test_data,2)>cov_dim),
            test_t=test_data(:,cov_dim+1:end);
        end
        disp(' ');
    end
else
    clc
    disp('----------------------------------------------------');
    disp('                  Specify the Input                 ');
    disp('----------------------------------------------------');
    disp(' ')
    ntrain = input('Insert the number of examples in training (<500): ');
    if (isempty(ntrain)|(floor(ntrain)~=ntrain)|(ntrain<=0)|(ntrain>500)),
          error(['Invalid value: ', num2str(ntrain), ' is not a positive integer <500!']);
    end
    disp(' ')
    test_path='toy';
    ntest = input('Insert the number of examples in test (<500): ');
    if (isempty(ntest)|(floor(ntest)~=ntest)|(ntest<=0)|(ntest>500)),
          error(['Invalid value: ', num2str(ntest), ' is not a positive integer <500!']);
    end

    if type==2,
        cov_dim=2;
        res_dim=3;
        seed = 42;
        [train_d, ntrain1, ntrain2, train_t]=gen_data(ntrain, seed);
        for m=1:ntrain
            q=[]; q = find(train_t(m,:)==1);
            train_data(m,:)=[train_d(m,:) q];
        end
        [test_d, ntest1, ntest2, test_t]=gen_data(ntest);
        for m=1:ntest
            q=[]; q = find(test_t(m,:)==1);
            test_data(m,:)=[test_d(m,:) q];
        end
    else
        cov_dim=1;
        res_dim=1;
        global HOME
        %%%%%WARNING!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        load([HOME '/examples/static/Misc/mixexp_data.txt'], '-ascii');
        %%%%%WARNING!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        train_data = mixexp_data(1:ntrain, :);
        train_d=train_data(:,1:cov_dim); train_t=train_data(:,cov_dim+1:end);
        test_data = mixexp_data(ntrain+1:ntrain+ntest, :);
        test_d=test_data(:,1:cov_dim);
        if size(test_data,2)>cov_dim,
            test_t=test_data(:,cov_dim+1:end);
        end
    end
end
% Set the nodes dimension-----------------------------------
if num_glevel>0,
    nodes_info(2,2:num_glevel+1)=branch_fact;
end
nodes_info(2,1)=cov_dim; nodes_info(2,end)=res_dim;
%-----------------------------------------------------------
% Prepare the training data for the learning engine---------
%-----------------------------------------------------------
cases = cell(size(nodes_info,2), ntrain);
for m=1:ntrain,
    cases{1,m}=train_data(m,1:cov_dim)';
    cases{end,m}=train_data(m,cov_dim+1:end)';
end
%-----------------------------------------------------------------------------------------------------
[bnet onodes]=hme_topobuilder(nodes_info);
engine = jtree_inf_engine(bnet, onodes);
clc
disp('---------------------------------------------------------------------');
disp('                         L  E  A  R  N  I  N  G                      ');
disp('---------------------------------------------------------------------');
disp(' ')
ll = 0;
for l=1:ntrain
  scritta=['example number: ', int2str(l),'---------------------------------------------'];
  disp(scritta);
  ev = cases(:,l);
  [engine, loglik] = enter_evidence(engine, ev);
  ll = ll + loglik;
end
disp(' ')
disp(['Log-likelihood before learning: ', num2str(ll)]);
disp(' ')
disp('(Press any key to continue)');
pause
%-----------------------------------------------------------
clc
disp('---------------------------------------------------------------------');
disp('                         L  E  A  R  N  I  N  G                      ');
disp('---------------------------------------------------------------------');
disp(' ')
max_em_iter=input('Insert the maximum number of the EM algorithm iterations: ');
if (isempty(max_em_iter)|(floor(max_em_iter)~=max_em_iter)|(max_em_iter<=1)),
          error(['Invalid value: ', num2str(ntest), ' is not a positive integer >1!']);
end
disp(' ')
disp(['Log-likelihood before learning: ', num2str(ll)]);
disp(' ')

[bnet2, LL2] = learn_params_em(engine, cases, max_em_iter);
disp(' ')
fprintf('HME: loglik before learning %f, after %d iters %f\n', ll, length(LL2),  LL2(end));
disp(' ')
disp('(Press any key to continue)');
pause
%-----------------------------------------------------------------------------------
% Classification problem: plot data & decision boundaries if the input data size = 2
% Regression problem: plot data & prediction if the input data size = 1
%-----------------------------------------------------------------------------------
if (type==2)&(nodes_info(2,1)==2)&(~isempty(test_path)),
    fh1=hme_class_plot(bnet2, nodes_info, train_data, test_data);
    disp(' ');
    disp('(See the figure)');
elseif (type==2)&(nodes_info(2,1)==2)&(isempty(test_path)),
    fh1=hme_class_plot(bnet2, nodes_info, train_data);
    disp(' ');
    disp('(See the figure)');
elseif (type==1)&(nodes_info(2,1)==1)&(~isempty(test_path)),
    fh1=hme_reg_plot(bnet2, nodes_info, train_data, test_data);
    disp(' ');
    disp('(See the figure)');
elseif (type==1)&(nodes_info(2,1)==1)&(isempty(test_path)),
    fh1=hme_reg_plot(bnet2, nodes_info, train_data);
    disp(' ')
    disp('(See the figure)');
end
%-----------------------------------------------------------------------------------
% Classification problem: plot confusion matrix
%-----------------------------------------------------------------------------------
if (type==2)
    ztrain=fhme(bnet2, nodes_info, train_d, size(train_d,1));
    [Htrain, trainRate]=confmat(ztrain, train_t); % CM on the training set
    fh2=figure('Name','Confusion matrix', 'MenuBar', 'none', 'NumberTitle', 'off');
    if (~isempty(test_path))&(size(test_data,2)>cov_dim),
        ztest=fhme(bnet2, nodes_info, test_d, size(test_d,1));
        [Htest, testRate]=confmat(ztest, test_t);   % CM on the test set
        subplot(1,2,1);
    end
    plotmat(Htrain,'b','k',12)
    tick=[0.5:1:(0.5+nodes_info(2,end)-1)];
    set(gca,'XTick',tick)
    set(gca,'YTick',tick)
    grid('off')
    ylabel('True')
    xlabel('Prediction')
    title(['Confusion Matrix: training set (' num2str(trainRate(1)) '%)'])
    if (~isempty(test_path))&(size(test_data,2)>cov_dim),
        subplot(1,2,2)
        plotmat(Htest,'b','k',12)
        set(gca,'XTick',tick)
        set(gca,'YTick',tick)
        grid('off')
        ylabel('True')
        xlabel('Prediction')
        title(['Confusion Matrix: test set (' num2str(testRate(1)) '%)'])
    end
    disp(' ')
    disp('(Press any key to continue)');
    pause
end
%-----------------------------------------------------------------------------------
% Regression & Classification problem: calculate the predictions & plot the LL trace
%-----------------------------------------------------------------------------------
train_result=fhme(bnet2,nodes_info,train_d,size(train_d,1));
if ~isempty(test_path),
    test_result=fhme(bnet2,nodes_info,test_d,size(test_d,1));
end
fh3=figure('Name','Log-likelihood trace', 'MenuBar', 'none', 'NumberTitle', 'off')
plot(LL2,'-ro',...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor',[1 1 0],...
                'MarkerSize',4)
title('Log-likelihood trace')
%-----------------------------------------------------------------------------------
% Regression & Classification problem: save the predictions
%-----------------------------------------------------------------------------------
clc
disp('------------------------------------------------------------------');
disp('                           Save the results                       ');
disp('------------------------------------------------------------------');
disp(' ')
%-----------------------------------------------------------------------------------
save_quest_m=input('Do you want to save the HME model (Y/N)? [Y default]: ', 's');
if isempty(save_quest_m),
    save_quest_m='Y';
end
if ~findstr(save_quest_m, ['Y', 'N']), error('Invalid input'); end
if save_quest_m=='Y',
    disp(' ');
    m_save=input('Insert the complete path for save the HME model (.mat):\n >> ', 's');
    if isempty(m_save), error('You must specify a path!'); end
    save(m_save, 'bnet2');
end
%-----------------------------------------------------------------------------------
disp(' ')
save_quest=input('Do you want to save the HME predictions (Y/N)? [Y default]: ', 's');
disp(' ')
if isempty(save_quest),
    save_quest='Y';
end
if ~findstr(save_quest, ['Y', 'N']), error('Invalid input'); end
if save_quest=='Y',
    tr_save=input('Insert the complete path for save the training data prediction (.mat):\n >> ', 's');
    if isempty(tr_save), error('You must specify a path!'); end
    save(tr_save, 'train_result');
    if ~isempty(test_path),
        disp(' ')
        te_save=input('Insert the complete path for save the test data prediction (.mat):\n >> ', 's');
        if isempty(te_save), error('You must specify a path!'); end
        save(te_save, 'test_result');
    end
end
clc
disp('----------------------------------------------------');
disp('                      B  Y  E !                     ');
disp('----------------------------------------------------');
pause(2)
%clear
clc

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function perr = p_err(t, pc, pk, qp, tt, ts)
    %p_err calculate p error
    sumln=sum(log(t+pc));
    qsumln=pk/qp^2;
    qsumln=qsumln*(((tt+pc)^qp)*(1-qp*log(tt+pc))-((ts+pc)^qp)*(1-qp*log(ts+pc)));
    esumln=qsumln+sumln;
    perr=esumln;
end

author: CelsoReyes


text:
function [Gi, ng] = Multi_Process(I, g, n)

if nargin < 3
    n = 6;
end
ng = g;
for i = 1:n
    ng = imdilate(ng, g);
end

Gi1 = imopen(I, ng);
Gi1 = imdilate(Gi1, ng);
Gi2 = imclose(I, ng);
Gi2 = imerode(Gi2, ng);
Gi = imsubtract(Gi1, Gi2);

author: decouples


text:
function result = MimReduceResultToContext(full_result, child_context, image_templates, dataset_stack, reporting)
    % MimReduceResultToContext. Part of the internal framework of the TD MIM Toolkit.
    %
    %     This function is used to convert a plugin result from one context
    %     to a lower context
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    % If the result is a composite result, then get the result for this
    % context
    if isa(full_result, 'MimCompositeResult') && full_result.IsField(char(child_context))
        result = full_result.(char(child_context));
        return
    end

    % Otherwise, we can only perform the reduction if the result it a
    % PTKImage
    if ~isa(full_result, 'PTKImage')
        result = full_result;
        return
    end

    % Resize an image to match the template for a particular context,
    % and extract out only the part of the image specified by the
    % context map

    template_image = image_templates.GetTemplateImage(child_context, dataset_stack, reporting);

    % Make a copy before we resize
    result = full_result.Copy;
    result.ResizeToMatch(template_image);

    if template_image.ImageExists
        result.Clear;
        result.ChangeSubImageWithMask(full_result, template_image, true);
    end
end


text:
function varargout=lpcifilt(varargin)
%
% For calling details please see v_lpcifilt.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpcifilt.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpcifilt.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpcifilt(varargin{:});
else
    v_lpcifilt(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [nScoutProj, destSurfMat, sAtlasProj] = bst_project_scouts( srcSurfFile, destSurfFile, sAtlas, isSingleHemi, isSave )
% BST_PROJECT_SCOUTS: Project scouts on a different surface (need the FreeSurfer registered spheres).
%
% USAGE:  [nScoutProj, destSurfMat, sAtlasProj] = bst_project_scouts( srcSurfFile, destSurfFile, sAtlas=[all], isSingleHemi=0, isSave=1 )

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2015-2023

% ===== PARSE INPUTS ======
if (nargin < 5) || isempty(isSave)
    isSave = 1;
end
if (nargin < 4) || isempty(isSingleHemi)
    isSingleHemi = 0;
end
if (nargin < 3) || isempty(sAtlas)
    sAtlas = [];
end
nScoutProj = 0;
destSurfMat = [];
sAtlasProj = sAtlas;

% ===== GET INTERPOLATION =====
% Make sure files are different
if file_compare(srcSurfFile, destSurfFile)
    disp('BST> Error: Source and destination surfaces are the same.');
    return;
end
% Compute interpolation
[Wmat, sSrcSubj, sDestSubj, srcSurfMat, destSurfMat] = tess_interp_tess2tess(srcSurfFile, destSurfFile, 1, [], isSingleHemi);
% Source subject and destination subject are the same
isSameSubject = file_compare(sSrcSubj.FileName, sDestSubj.FileName);
% If no scouts in input, project everything
if isempty(sAtlas)
    sAtlas = srcSurfMat.Atlas;
end
% Check if there are things to project
if isempty(sAtlas)
    disp('BST> Error: No scouts to project.');
    return;
end
% Ratio of vertex number
ratio = size(destSurfMat.Vertices,1) ./ size(srcSurfMat.Vertices,1);

% ===== PROCESS ATLAS/SCOUTS =====
for iAtlas = 1:length(sAtlas)
    % Initialize probability maps
    scoutIndex = zeros(size(destSurfMat.Vertices,1),1);
    scoutProba = zeros(size(destSurfMat.Vertices,1),1);
    % Project scouts one by one and keep for each vertex only the maximum probability
    for iScout = 1:length(sAtlas(iAtlas).Scouts)
        % Vertex map on the original surface
        vMap = zeros(size(srcSurfMat.Vertices,1),1);
        vMap(sAtlas(iAtlas).Scouts(iScout).Vertices) = 1;
        % Project to destination surface
        vMapProj = full(Wmat * vMap);
        % Keep only the projections that have a higher probability than the previous scouts
        isHigherProba = vMapProj > scoutProba;
        scoutProba(isHigherProba) = vMapProj(isHigherProba);
        scoutIndex(isHigherProba) = iScout;
    end

% DISABLED 2018
%         % If the number of vertices does not decrease: force the selection of the closest vertex to each input vertex
%         if (ratio > 0.9)
%             for iVert = 1:length(sScout.Vertices)
%                 % Get the closest projected vertex for iVert
%                 [tmp, iVertClosest] = max(Wmat(:, sScout.Vertices(iVert)));
%                 % Force the selection by setting the interpolated value higher than 1
%                 vMapProj(iVertClosest) = 2;
%             end
%         end

    % Create all the scouts in the destination surface
    for iScout = 1:length(sAtlas(iAtlas).Scouts)
        % Current scout
        sScout = sAtlas(iAtlas).Scouts(iScout);

        % Get vertices identified in this scout
        iVertices = find(scoutIndex == iScout);
        if isempty(iVertices)
            sAtlasProj(iAtlas).Scouts(iScout).Vertices = [];
            sAtlasProj(iAtlas).Scouts(iScout).Seed     = [];
            continue;
        end
        % Limit the growth to extra vertices when not projecting an entire atlas
        if (length(sAtlas(iAtlas).Scouts) < 10)
            % Sort the projected values and keep the highest ones, up to desired number of vertices
            [tmp,I] = sort(scoutProba(iVertices), 1, 'descend');
            % Keep the highest values
            nVertices = round(ratio * length(sScout.Vertices));
            iVertices = iVertices(I(1:min(nVertices,length(I))));
        end

        % Identify seed (closest point to the center of mass of the scout)
        c = mean(destSurfMat.Vertices(iVertices,:),1);
        distC = sqrt((destSurfMat.Vertices(iVertices,1)-c(1)).^2 + (destSurfMat.Vertices(iVertices,2)-c(2)).^2 + (destSurfMat.Vertices(iVertices,3)-c(3)).^2);
        [distMin,iMin] = min(distC);
        iSeed = iVertices(iMin);

        % Get destination atlas
        iAtlasDest = find(strcmpi({destSurfMat.Atlas.Name}, sAtlas(iAtlas).Name));
        if isempty(iAtlasDest)
            iAtlasDest = length(destSurfMat.Atlas) + 1;
            destSurfMat.Atlas(iAtlasDest).Name = sAtlas(iAtlas).Name;
        end
        % Destination scout name
        if ~isSave
            ScoutLabel = sScout.Label;
        elseif isSameSubject
            ScoutLabel = sScout.Label;
        else
            ScoutLabel = [sSrcSubj.Name '_' sScout.Label];
            ScoutLabel = strrep(ScoutLabel, '@default_subject', 'Default');
        end
        if ~isempty(destSurfMat.Atlas(iAtlasDest).Scouts)
            ScoutLabel = file_unique(ScoutLabel, {destSurfMat.Atlas(iAtlasDest).Scouts.Label});
        end
        % Create new scout
        iScoutDest = length(destSurfMat.Atlas(iAtlasDest).Scouts) + 1;
        destSurfMat.Atlas(iAtlasDest).Scouts(iScoutDest).Vertices = iVertices(:)';
        destSurfMat.Atlas(iAtlasDest).Scouts(iScoutDest).Seed     = iSeed;
        destSurfMat.Atlas(iAtlasDest).Scouts(iScoutDest).Color    = sScout.Color;
        destSurfMat.Atlas(iAtlasDest).Scouts(iScoutDest).Label    = ScoutLabel;
        destSurfMat.Atlas(iAtlasDest).Scouts(iScoutDest).Function = sScout.Function;
        destSurfMat.Atlas(iAtlasDest).Scouts(iScoutDest).Region   = sScout.Region;
        % Report projected scouts
        nScoutProj = nScoutProj + 1;
        sAtlasProj(iAtlas).Scouts(iScout).Vertices = iVertices(:)';
        sAtlasProj(iAtlas).Scouts(iScout).Seed     = iSeed;
    end
end

% Save destination surface (append the atlas to existing file)
if isSave && (nScoutProj > 0)
    s.Atlas = destSurfMat.Atlas;
    bst_save(file_fullpath(destSurfFile), s, 'v7', 1);
end


text:
% Spherical Harmonics basis of order lmax.
% NB: the definition now matches the one in DIPY
function Ylm = AMICO_CreateYlm( lmax, colatitude, longitude )

Ylm = zeros( size(longitude,1), (lmax+2)*(lmax+1)/2 );
for l = 0:2:lmax
    Pm = legendre(l,cos(colatitude'))';
    lconstant = sqrt((2*l + 1)/(4*pi));

    center = (l+1)*(l+2)/2 - l;
    Ylm(:,center) = lconstant*Pm(:,1);
    for m=1:l
        precoeff = lconstant * sqrt(factorial(l - m)/factorial(l + m));

		Ylm(:, center - m) = sqrt(2)*precoeff*Pm(:,m+1).*cos(m*longitude);
		Ylm(:, center + m) = sqrt(2)*precoeff*Pm(:,m+1).*sin(m*longitude);
	end
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function title = p06_title ( )

%*****************************************************************************80
%
%% P06_TITLE returns the title of problem P06.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Problem 6, Enright and Pryce #B1';

  return
end


text:
function approxhess = getApproxHessian(problem, x, d, storedb, key)
% Computes an approximation of the Hessian of the cost fun. at x along d.
%
% function approxhess = getApproxHessian(problem, x, d)
% function approxhess = getApproxHessian(problem, x, d, storedb)
% function approxhess = getApproxHessian(problem, x, d, storedb, key)
%
% Returns an approximation of the Hessian at x along d of the cost function
% described in the problem structure.
%
% storedb is a StoreDB object, key is the StoreDB key to point x.
%
% If no approximate Hessian was provided, this call is redirected to
% getHessianFD.
%
% See also: getHessianFD canGetApproxHessian

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.

    % Allow omission of the key, and even of storedb.
    if ~exist('key', 'var')
        if ~exist('storedb', 'var')
            storedb = StoreDB();
        end
        key = storedb.getNewKey();
    end

    if isfield(problem, 'approxhess')
    %% Compute the approximate Hessian using approxhess.

        % Check whether this function wants to deal with storedb or not.
        switch nargin(problem.approxhess)
            case 2
                approxhess = problem.approxhess(x, d);
            case 3
                % Obtain, pass along, and save the store for x.
                store = storedb.getWithShared(key);
                [approxhess, store] = problem.approxhess(x, d, store);
                storedb.setWithShared(store, key);
            case 4
                % Pass along the whole storedb (by reference), with key.
                approxhess = problem.approxhess(x, d, storedb, key);
            otherwise
                up = MException('manopt:getApproxHessian:badapproxhess', ...
                    'approxhess should accept 2, 3 or 4 inputs.');
                throw(up);
        end

    else
    %% Try to fall back to a standard FD approximation.

        approxhess = getHessianFD(problem, x, d, storedb, key);

    end

end

author: NicolasBoumal

repo: manopt


text:
if frame == 1
    bias = min(w(:));
end
if (frame==1||mod(frame,update_interval)==0)&&params.show_regularization
    figure(2);
    I = im;
    init_rect = rect_position(loop_frame,:);
    I = insertShape(I, 'Rectangle', init_rect, 'LineWidth', 2, 'Color', 'red');
    padding = 2.0;
    rect_paddding = [init_rect(1)-init_rect(3)*padding/2,...
    init_rect(2)-init_rect(4)*padding/2,...
    init_rect(3)*(1+padding),init_rect(4)*(1+padding)];
    crop = imresize(imcrop(I, rect_paddding), [125, 125]);
    [h, width, c] = size(crop);

    startw=round(size(w,1)*0.4/2)+1;
    starth=round(size(w,2)*0.4/2)+1;
    w_=w(startw:startw+round(size(w,1)*0.6)-1,starth:starth+round(size(w,2)*0.6)-1);
    w_=imresize(w_,[h,width]);

    w_=(w_-bias+0.1);
    surf(w_);
    %surf(X,Y,(Z-min(Z(:))),'FaceAlpha',0.4);colormap(hsv);hold on

    g = hgtransform('Matrix',makehgtform('translate',[0 0 0]));
    image(g, crop)

    axis off
    view(205, 30)
    set(gcf, 'position', [100 100 900 900], 'Color',[0,0,0]);
end

author: Daikenan

repo: ASRCF

save_path: github-repos/MATLAB/Daikenan-ASRCF


text:
% Function which propagates and adds noise to the 3D points and popoulates
% other variables of the structs passed to the function
% motion = [trans, rot_yaw_theta]
function [detectionsQ, detectionsT] = propagateDetections(detectionsQ, detectionsT, params_2D3D, params_carCuboid, motion, first_time)

%% for all detections in Query propagate the box to 3D box in F2 frame

    for i = 1:length(detectionsQ)

        b1Q = [];
        B1Q = [];
        bvolumeQ1 = [];

        % check if the detection is first time i.e. we have to start with
        % the canonical cuboid. Else we already have the cuboid.

        if(length(detectionsQ{i}.bvolume) == 0 )
            %disp(sprintf('first_time - id :%d\n',detectionsQ{i}.dno));
            canonicalCuboid = getCanonicalCuboid( params_carCuboid.avg_Sz);

            % find center of the bounding box bottom line
            b1Q = [detectionsQ{i}.bbox(1) + (detectionsQ{i}.bbox(3) - detectionsQ{i}.bbox(1))/2 ;
                  detectionsQ{i}.bbox(4);
                  1.0];

            % project it to 3D
            B1Q = (params_2D3D.h * inv(params_2D3D.K) * b1Q) / (params_2D3D.n' * inv(params_2D3D.K) * b1Q);

            % apply offset which is a function of yaw and get car's origin (remember approx.)
            offset_Z = getOffsetBasedOnYaw([params_carCuboid.avg_Sz(1); params_carCuboid.avg_Sz(3)], detectionsQ{i}.yaw);

            % car's origin in Query Frame
            B1Q = B1Q + [0; 0; offset_Z];

            % translate canonical cuboid
            canonicalCuboid = canonicalCuboid + repmat(B1Q', 8,1);

            % BOUNDING VOLUME IN QUERY FRAME
            [bvolumeQ1, k1] = getBoundingVolume(B1Q, canonicalCuboid, detectionsQ{i}.yaw, ...
                                                detectionsQ{i}.sigma_3D, ...
                                                params_carCuboid.sz_ub, params_carCuboid.sz_lb ...
                                               );

            bvolumeQ1 = bvolumeQ1 - repmat([0 params_carCuboid.avg_Sz(2)/2 0], size(bvolumeQ1,1), 1);
        else

            B1Q = detectionsQ{i}.origin;
            bvolumeQ1 = detectionsQ{i}.bvolume;

        end

        % TODO
             % propagate noise and find new noise .... for now same noise
        %

        % car's origin in Train frame
        B2Q = motion(1:3) + B1Q ;
        bvolumeQ1 = bvolumeQ1 + repmat(motion(1:3)', size(bvolumeQ1,1),1);
        % CUBOID IN TRAIN FRAME

        %dummy_sigma = eye(4,4);
        %dumm_sigma(4,4) = 0;

        [bvolumeQ2, k2] = getBoundingVolume(B2Q, bvolumeQ1, motion(4), ...
                                             detectionsQ{i}.sigma_3D, ...
                                             params_carCuboid.sz_ub, params_carCuboid.sz_lb ...
                                            );

       % bvolumeQ2 = bvolumeQ2 - repmat([0 params_carCuboid.avg_Sz(2)/2 0], size(bvolumeQ2,1), 1);    % offset the h/2 as car ar to be on road

        % compute projection of bvolume and store in bvolume_proj variable
        bvolume_proj = (params_2D3D.K*bvolumeQ2')';
        bvolume_proj(:,1:3) = bvolume_proj(:,1:3)./repmat(bvolume_proj(:,3), 1,3);
        kbvolume_proj = convhull(bvolume_proj(:,1:2));
        bvolume_proj = bvolume_proj(kbvolume_proj,:);

        % update the fields of the structure

        detectionsQ{i}.bvolume = bvolumeQ2;
        detectionsQ{i}.bvolume_proj = bvolume_proj;
        detectionsQ{i}.yaw  = detectionsQ{i}.yaw+motion(4);
        detectionsQ{i}.origin  = B2Q;
        detectionsQ{i}.k = k2;

    end

    %% for all detections in Train propagate the box to 3D box in its frame

    for i = 1:length(detectionsT)

        canonicalCuboid = getCanonicalCuboid( params_carCuboid.avg_Sz);

        % find center of the bounding box bottom line
        b1T = [detectionsT{i}.bbox(1) + (detectionsT{i}.bbox(3) - detectionsT{i}.bbox(1))/2 ;
              detectionsT{i}.bbox(4);
              1.0];

        % project it to 3D
        B1T = (params_2D3D.h * inv(params_2D3D.K) * b1T) / (params_2D3D.n' * inv(params_2D3D.K) * b1T);

        % apply offset which is a function of yaw and get car's origin (remember approx.)
        offset_Z = getOffsetBasedOnYaw([params_carCuboid.avg_Sz(1); params_carCuboid.avg_Sz(3)], detectionsT{i}.yaw);

        % car's origin in Query Frame
        B1T = B1T + [0; 0; offset_Z];

        % translate canonical cuboid
        canonicalCuboid = canonicalCuboid + repmat(B1T', 8,1);

        % BOUNDING VOLUME IN QUERY FRAME
        [bvolumeT1, k1] = getBoundingVolume(B1T, canonicalCuboid, detectionsT{i}.yaw, ...
                                            detectionsT{i}.sigma_3D, ...
                                            params_carCuboid.sz_ub, params_carCuboid.sz_lb ...
                                           );
        % update the fields of the structure
        detectionsT{i}.origin  = B1T;
        detectionsT{i}.k = k1;

        detectionsT{i}.bvolume = bvolumeT1 - repmat([0 params_carCuboid.avg_Sz(2)/2 0], size(bvolumeT1,1), 1);    % offset the h/2 as car ar to be on road             ;

        bvolume_proj = detectionsT{i}.bvolume;
        bvolume_proj = (params_2D3D.K*bvolume_proj')';
        bvolume_proj(:,1:3) = bvolume_proj(:,1:3)./repmat(bvolume_proj(:,3), 1,3);
        kbvolume_proj = convhull(bvolume_proj(:,1:2));
        bvolume_proj = bvolume_proj(kbvolume_proj,:);

        detectionsT{i}.bvolume_proj = bvolume_proj;
    end

end

repo: MOTBeyondPixels


text:
function [ X ] = special_add_left( X, S )
% add on Lie group-Lie Algebra space
    s_theta=S(1:3);
    s_p=S(4:6);
    sizeS=size(S,1);
    NumberOfLandmarks=(sizeS-6)/3;
    X.position=X.position+X.orientation*jaco_r(-s_theta)*s_p;
    if NumberOfLandmarks>=1
        s_landmarksMatrix=reshape(S(7:end),3,NumberOfLandmarks);
        X.landmarks(1:3,:)=X.landmarks(1:3,:)+X.orientation*jaco_r(-s_theta)*s_landmarksMatrix;
    end
    X.orientation=X.orientation*so3_exp(s_theta);
end

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
function dist = codedist_bay(C1,C2,Py)
% compute the distance using a Bayesian metric between 'C1' and  'C2'
%
%    dist = codedist_bay(C1,C2,{Py})
%
% 'C1' contains the result code,
% 'C2' contains the codebooks prototype. '0' or  an infine small number 'eps'  represents the don't care
% 'Py' is a matrix containing the moderated outputs of the binary classifiers
%
%
% An example:
%  >> Ye = [1 1; 1 1];
%  >> codebook = [1 2 3];
%  >> old_codebook = [1 1 -1; 1 -1 -1];
%  >> Py = [.9 -.1; -.1 .13];
%  >> code(Ye, codebook, [],old_codebook,'codedist_bay',{Py})
% in this call, 'Ye' is not used explicitly.
%
% To use this distance measure in LS-SVMlab, the following
% procedure is to be followed, assume input data 'X' and multiclass
% output 'Y':
%
%  >> [Ycode,codebook,old_codebook] = code(Y,'code_MOC');
%  >> [alpha,b] = trainlssvm({X,Yc,'c',gam,sig2});
%  >> Yhc = simlssvm({X,Yc,'c',gam,sig2},{alpha,b},Xt);
%
%  The moderated output for the LS-SVM can be computed using the bayesian inference
%  framework for the LS-SVM:
%  >> Ymod = bay_modoutClass(model,Xt);
%  >> Yh = code(Yhc,old_codebook,[],codebook,'codedist_bay',{Ymod});
%
% see also:
%   bay_modoutClass, codedist_hamming, code_ECOC

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

% % encode for training
% >> model = initlssvm(X,Y,'classification',gam,sig2,'preprocess','RBF_kernel');
% >> model = changelssvm(model,'codetype','code_MOC');
% >> model = changelssvm(model,'codedist_fct','codedist_hamming');
% >> model = trainlssvm(model);
%
% % decode for simulating
% >> model = changelssvm(model,'codedist_fct','codedist_bay');
% >> model = changelssvm(model,'codedist_args',{bay_modoutClass(model,Xt)});
% >>   Yt  = simlssvm(model,Xt);

if nargin<3,
  error(['moderated output needed as function arguments' ...
	 '(for LS-SVM, model.codedist_args =' ...
	 ' bay_modoutClass(model,X)).']);
end

[nb,nbin] = size(Py);
[~,dim] = size(C2);
dist = zeros(nb,dim);

for d = 1:dim,
  for n= 1:nb,
    dist(n,d) = sum((1-Py(n,:).*C2(:,d)'))-sum(C2(:,d)==eps);
  end
end

author: peterhcharlton

repo: RRest


text:
function pass = test_vdpIVP(~)

%% Setup
vdpFun = @(mu) @(u) diff(u, 2) - mu.*(1-u.^2).*diff(u) + u;
dom = [0 10];
bcFun =  @(u) [u-2; diff(u)];
opts = odeset('abstol', 1e5*eps, 'reltol', 100*eps);
tol = 1e-10;

%% Solve with CHEBOPs, compare with standard MATLAB, mu = 1
N = chebop(vdpFun(1), dom);
N.lbc = bcFun;
u = N\0;
uend = u(dom(end));
sol =ode45(@vdp1, dom, [2 0], opts);
yend = deval(sol, dom(end)); yend = yend(1);
err(1) = abs(uend - yend);

%% Solve with CHEBOPs, compare with standard MATLAB, mu = 1000
N = chebop(vdpFun(1000), dom);
N.lbc = bcFun;
u = N\0;
uend = u(dom(end));
sol =ode45(@vdp1000, dom, [2 0], opts);
yend = deval(sol, dom(end)); yend = yend(1);
err(2) = abs(uend - yend);

%% Happy?
pass = err < tol;

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function m = uplus(k)
%  kronMatrix/uplus
%
%     +k = k
%

m = k;

repo: IRtools