text:
function pass = test_removeDeltas(pref)

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

% Test that the deltas are removed.
x = chebfun(@(x) x, [0 1], pref);
f = dirac(x - 0.2) + diff(dirac(x-1), 2);
thenorm = norm(removeDeltas(f), Inf); % Is infinite if deltas are present.
pass(1) = ( thenorm < Inf );

% Test that pointValues are changed.
f = dirac(x);
val = feval(removeDeltas(f), 0);
pass(2) = ( val == 0 );

% Test that removeDeltas is the identity for classicfuns.
f = sin(x);
g = removeDeltas(f);
pass(3) = isequal(f, g);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef SetEngineActiveStateAction < AbstractEventAction
    %SetEngineActiveStateAction Summary of this class goes here
    %   Detailed explanation goes here

    properties
        engine LaunchVehicleEngine
        activeStateToSet(1,1) logical = false;
    end

    properties(Constant)
        emptyVarArr = AbstractOptimizationVariable.empty(0,1);
    end

    methods
        function obj = SetEngineActiveStateAction(engine, activeStateToSet)
            if(nargin > 0)
                obj.engine = engine;
                obj.activeStateToSet = activeStateToSet;
            end

            obj.id = rand();
        end

        function newStateLogEntry = executeAction(obj, stateLogEntry)
            newStateLogEntry = stateLogEntry;
            stage = obj.engine.stage;

            stgState = newStateLogEntry.stageStates([newStateLogEntry.stageStates.stage] == stage);

            engineState = stgState.engineStates([stgState.engineStates.engine] == obj.engine);
            engineState.active = obj.activeStateToSet;
        end

        function initAction(obj, initialStateLogEntry)
            %nothing
        end

        function name = getName(obj)
            if(obj.activeStateToSet)
                tf = 'Active';
            else
                tf = 'Inactive';
            end

            name = sprintf('Set Engine State (%s = %s)',obj.engine.name,tf);
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = ([obj.engine] == engine);
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end

        function tf = usesExtremum(obj, extremum)
            tf = false;
        end

        function tf = usesTankToTankConn(obj, tankToTank)
            tf = false;
        end

        function [tf, vars] = hasActiveOptimVar(obj)
            tf = false;
            vars = obj.emptyVarArr;
        end
    end

    methods(Static)
        function addActionTf = openEditActionUI(action, lv)
%             addActionTf = lvd_EditActionSetEngineStateGUI(action, lv);

            output = AppDesignerGUIOutput({false});
            lvd_EditActionSetEngineStateGUI_App(action, lv, output);
            addActionTf = output.output{1};
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% converts all double variables to mp objects
% put variable names as cell strings into mp_dontmakeall if you don't want them converted
%  e.g. mp_dontmakeall={'C1','C2'};
%  if you don't want C1 and C2 converted to mp objects

tempVarList=whos;

if exist('default_precision')~=1
 default_precision=mp_Defaults(mp(0));
end
for tempLoopVar=1:length(tempVarList)
 if strcmp(tempVarList(tempLoopVar).class,'double')
  tempGoon=1;
  if exist('mp_dontmakeall')
   if any(strcmp(tempVarList(tempLoopVar).name,mp_dontmakeall))
    tempGoon=0;
   end
  end
  if tempGoon
   if ~strcmp(tempVarList(tempLoopVar).name,'default_precision')
    eval([tempVarList(tempLoopVar).name,'=mp(',tempVarList(tempLoopVar).name,...
          ',',num2str(default_precision),');']);
   end
  end
 end
end

clear tempVarList tempLoopVar tempGoon

author: opencobra

repo: cobratoolbox


text:
% Data File MAH3
% Triple Pendulum
  s = 3; % degree of freedom
%  L a g r a n g i a n   of the system
  L = [' l1^2/6*(m1 + 3*m2 + 3*m3)*qt1^2 + '        , ...
       ' l2^2/6*(m2 + 3*m3)*qt2^2 + '               , ...
       ' l3^2/6*m3*qt3^2 + '                        , ...
       ' 1/2*(m2+2*m3)*l1*l2*cos(q1-q2)*qt1*qt2 + ' , ...
       ' 1/2*m3*l1*l3*cos(q1-q3)*qt1*qt3 + '        , ...
       ' 1/2*m3*l2*l3*cos(q2-q3)*qt2*qt3 + '        , ...
       ' 9.81*l1*(m1/2 + m2 + m3)*cos(q1) + '          , ...
       ' 9.81*l2*(m2/2 + m3)*cos(q2) + '               , ...
       ' 9.81*l3*m3/2*cos(q3) '              ];
  QN{1} = '-k1*qt1'; % Generalized
  QN{2} = '-k2*qt2'; % non potential
  QN{3} = '-k3*qt3'; % forces
  qj0   = [0, 0, 0]; % Initial coordinates
  qtj0  = [5, 0, 0]; % Initial velocities
  Tend  = 20;        % Upper bound of integration
  eps   = 1e-10;     % Desirable accuracy
  np    = 9;         % Number of parameters
  P{1}  = 'm1';      % Masses of
  P{2}  = 'm2';      % the rods 1, 2, 3
  P{3}  = 'm3';
  P{4}  = 'l1';      % Lenghts of
  P{5}  = 'l2';      % the rods 1, 2, 3
  P{6}  = 'l3';
  P{7}  = 'k1';      % Coefficients of
  P{8}  = 'k2';      % damping
  P{9}  = 'k3';

author: Sable

repo: mcbench-benchmarks


text:
classdef NormDescriptor < handle

    properties (Access = protected)
        dist
    end

    methods (Access = protected, Abstract)
        computeDistance(obj)
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef xIMUdataClass < handle

    %% Public properties
    properties (SetAccess = private)
        FileNamePrefix = '';
        ErrorData = [];
        CommandData = [];
        RegisterData = [];
        DateTimeData = [];
        RawBatteryAndThermometerData = [];
        CalBatteryAndThermometerData = [];
        RawInertialAndMagneticData = [];
        CalInertialAndMagneticData = [];
        QuaternionData = [];
        RotationMatrixData = [];
        EulerAnglesData = [];
        DigitalIOdata = [];
        RawAnalogueInputData = [];
        CalAnalogueInputData = [];
        PWMoutputData = [];
        RawADXL345busData = [];
        CalADXL345busData = [];
    end

    %% Public methods
    methods (Access = public)
        function obj = xIMUdataClass(varargin)
            % Create data objects from files
            obj.FileNamePrefix = varargin{1};
            dataImported = false;
            try obj.ErrorData = ErrorDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.CommandData = CommandDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.RegisterData = RegisterDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.DateTimeData = DateTimeDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.RawBatteryAndThermometerData = RawBatteryAndThermometerDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.CalBatteryAndThermometerData = CalBatteryAndThermometerDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.RawInertialAndMagneticData = RawInertialAndMagneticDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.CalInertialAndMagneticData = CalInertialAndMagneticDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.QuaternionData = QuaternionDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.EulerAnglesData = EulerAnglesDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.RotationMatrixData = RotationMatrixDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.DigitalIOdata = DigitalIOdataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.RawAnalogueInputData = RawAnalogueInputDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.CalAnalogueInputData = CalAnalogueInputDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.PWMoutputData = PWMoutputDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.RawADXL345busData = RawADXL345busDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            try obj.CalADXL345busData = CalADXL345busDataClass(obj.FileNamePrefix); dataImported = true; catch e, end
            if(~dataImported)
                error('No data was imported.');
            end

            % Apply SampleRate from register data
            try h = obj.DateTimeData; h.SampleRate = obj.SampleRateFromRegValue(obj.RegisterData.GetValueAtAddress(67)); catch e, end
            try h = obj.RawBatteryAndThermometerData; h.SampleRate = obj.SampleRateFromRegValue(obj.RegisterData.GetValueAtAddress(68)); catch e, end
            try h = obj.CalBatteryAndThermometerData; h.SampleRate = obj.SampleRateFromRegValue(obj.RegisterData.GetValueAtAddress(68)); catch e, end
            try h = obj.RawInertialAndMagneticData; h.SampleRate = obj.SampleRateFromRegValue(obj.RegisterData.GetValueAtAddress(69)); catch e, end
            try h = obj.CalInertialAndMagneticData; h.SampleRate = obj.SampleRateFromRegValue(obj.RegisterData.GetValueAtAddress(69)); catch e, end
            try h = obj.QuaternionData; h.SampleRate = obj.SampleRateFromRegValue(obj.RegisterData.GetValueAtAddress(70)); catch e, end
            try h = obj.RotationMatrixData; h.SampleRate = obj.SampleRateFromRegValue(obj.RegisterData.GetValueAtAddress(70)); catch e, end
            try h = obj.EulerAnglesData; h.SampleRate = obj.SampleRateFromRegValue(obj.RegisterData.GetValueAtAddress(70)); catch e, end
            try h = obj.DigitalIOdata; h.SampleRate = obj.SampleRateFromRegValue(obj.RegisterData.GetValueAtAddress(78)); catch e, end
            try h = obj.RawAnalogueInputData; h.SampleRate = obj.SampleRateFromRegValue(obj.RegisterData.GetValueAtAddress(80)); catch e, end
            try h = obj.CalAnalogueInputData; h.SampleRate = obj.SampleRateFromRegValue(obj.RegisterData.GetValueAtAddress(80)); catch e, end
            try h = obj.RawADXL345busData; h.SampleRate = obj.SampleRateFromRegValue(obj.RegisterData.GetValueAtAddress(85)); catch e, end
            try h = obj.CalADXL345busData; h.SampleRate = obj.SampleRateFromRegValue(obj.RegisterData.GetValueAtAddress(85)); catch e, end

            % Apply SampleRate if specified as argument
            for i = 2:2:(nargin)
                if strcmp(varargin{i}, 'DateTimeSampleRate')
                    try h = obj.DateTimeData; h.SampleRate = varargin{i+1}; catch e, end
                elseif strcmp(varargin{i}, 'BattThermSampleRate')
                    try h = obj.RawBatteryAndThermometerData; h.SampleRate = varargin{i+1}; catch e, end
                    try h = obj.CalBatteryAndThermometerData; h.SampleRate = varargin{i+1}; catch e, end
                elseif strcmp(varargin{i}, 'InertialMagneticSampleRate')
                    try h = obj.RawInertialAndMagneticData; h.SampleRate = varargin{i+1}; catch e, end
                    try h = obj.CalInertialAndMagneticData; h.SampleRate = varargin{i+1}; catch e, end
                elseif strcmp(varargin{i}, 'QuaternionSampleRate')
                    try h = obj.QuaternionData; h.SampleRate = varargin{i+1}; catch e, end
                    try h = obj.RotationMatrixData.SampleRate; h.SampleRate = varargin{i+1}; catch e, end
                    try h = obj.EulerAnglesData; h.SampleRate = varargin{i+1}; catch e, end
                elseif strcmp(varargin{i}, 'DigitalIOSampleRate')
                    try h = obj.DigitalIOdata; h.SampleRate = varargin{i+1}; catch e, end
                elseif strcmp(varargin{i}, 'AnalogueInputSampleRate')
                    try h = obj.RawAnalogueInputData; h.SampleRate = varargin{i+1}; catch e, end
                    try h = obj.CalAnalogueInputData; h.SampleRate = varargin{i+1}; catch e, end
                elseif strcmp(varargin{i}, 'ADXL345SampleRate')
                    try h = obj.RawADXL345busData; h.SampleRate = varargin{i+1}; catch e, end
                    try h = obj.CalADXL345busData; h.SampleRate = varargin{i+1}; catch e, end
                else
                    error('Invalid argument.');
                end
            end
        end
        function obj = Plot(obj)
            try obj.RawBatteryAndThermometerData.Plot(); catch e, end
            try obj.CalBatteryAndThermometerData.Plot(); catch e, end
            try obj.RawInertialAndMagneticData.Plot(); catch e, end
            try obj.CalInertialAndMagneticData.Plot(); catch e, end
            try obj.QuaternionData.Plot(); catch e, end
            try obj.EulerAnglesData.Plot(); catch e, end
            try obj.RotationMatrixDataClass.Plot(); catch e, end
            try obj.DigitalIOdata.Plot(); catch e, end
            try obj.RawAnalogueInputData.Plot(); catch e, end
            try obj.CalAnalogueInputData.Plot(); catch e, end
            try obj.RawADXL345busData.Plot(); catch e, end
            try obj.CalADXL345busData.Plot(); catch e, end
        end
    end

    %% Private methods
    methods (Access = private)
        function sampleRate = SampleRateFromRegValue(obj, value)
            sampleRate = floor(2^(value-1));
        end
    end
end

author: xioTechnologies

repo: Gait-Tracking-With-x-IMU


text:
%Author: Moeti Ncube

%Comments: This code calibrates the heston model to any dataset of the form
%of the marketdata.txt file. It also simulates the heston model given the
%optimized parameters.

%Marketdata=[strike, maturity,impvol,time till expiry]

%Strike: Strike price of options
%Maturity: (1=Prompt month, 2=Prompt month+1,...)
%Impvol: Market Implied vol
%Time till expiry: Days until option expires

clear all

global impvol; global strike; global T; global F0; global r;

%Initial Parameter Guess for Hestion Model
%V(1), Kappa, Theta, Vol of Volatility (sig), Correlation (rho)
x0=[.5,.5,.5,.05,.5];
%Constraints (Lower and Upper bounds on Parameters)
lb = [0, 0, 0, 0, -.9];
ub = [1, 100, 1, .5, .9];
%Number of MCMC simulations of Heston Model
M=50000;
%Current Asset Price (Prompt Month Price) and Interest rate
F0=1250; r=0;

load marketdata.txt

T=marketdata(:,4)/365;
strike=marketdata(:,1);
impvol=marketdata(:,3);

%Optimization
x = lsqnonlin(@costf2,x0,lb,ub);

for k=1:length(T);
%Initial asset price
shes(1)=F0;
%Number of Time Steps,time step size
N=round(T(k)/(1/360));dt=T(k)/N;

%Heston Parameters
vhes(1)=x(1);  kappa=x(2); theta=x(3); vsigma=x(4);rho=x(5); simPath=0;

%Simulation of Heston Model
for i = 1:M
for j=1:N
%heston model
r1 = randn;
r2 = rho*r1+sqrt(1-rho^2)*randn;
shes(j+1)=shes(j)*exp((-0.5*vhes(j))*dt+sqrt(vhes(j))*sqrt(dt)*r1);
vhes(j+1)=vhes(j)*exp(((kappa*(theta - vhes(j))-0.5*vsigma^2)*dt)/vhes(j) + vsigma*(1/sqrt(vhes(j)))*sqrt(dt)*r2);
end
simPath = simPath + exp(-r*T(k)) * max(shes(j+1) - strike(k), 0);
end
simhes(k)=simPath/M;
simimpvol(k)=blkimpv(shes(1), strike(k), r, T(k), simhes(k));
modhes(k)=HestonCall(shes(1),strike(k),r,T(k),vhes(1),kappa,theta,vsigma,rho,0);
hesimpvol(k)=blkimpv(shes(1), strike(k), r, T(k), modhes(k));
end

for i=1:length(T)
bsprice(i)=blsprice(F0,strike(i),r,T(i),impvol(i));
end

%Output optimized Parameters
x

%Compare blackscholes option price, analytical heston price, and simulated
%heston prices for a given maturity and strike
pricedata=[T,strike,bsprice',modhes',simhes']

%Compare blackscholes option IV, analytical heston IV, and simulated
%heston IV for a given maturity and strike

voldata=[T,strike,impvol,hesimpvol',simimpvol']


text:
classdef GeometricPointSet < matlab.mixin.SetGet
    %GeometricPointSet Summary of this class goes here
    %   Detailed explanation goes here

    properties
        points(1,:) AbstractGeometricPoint

        lvdData LvdData
    end

    methods
        function obj = GeometricPointSet(lvdData)
            obj.lvdData = lvdData;
        end

        function addPoint(obj, point)
            obj.points(end+1) = point;
        end

        function removePoint(obj, point)
            obj.points([obj.points] == point) = [];
        end

        function [listBoxStr, points] = getListboxStr(obj)
            listBoxStr = {};

            for(i=1:length(obj.points)) %#ok<*NO4LP>
                listBoxStr{end+1} = obj.points(i).getListboxStr(); %#ok<AGROW>
            end

            points = obj.points;
        end

        function indPoints = getPointsForInds(obj, inds)
            indPoints = obj.points(inds);
        end

        function inds = getIndsForPoints(obj, indPoints)
            inds = [];
            for(i=1:length(indPoints))
                for(j=1:length(obj.points))
                    if(indPoints(i) == obj.points(j))
                        inds(end+1) = j; %#ok<AGROW>
                        break;
                    end
                end
            end
        end

        function indPoint = getPointAtInd(obj, ind)
            if(ind > 0 && ind <= length(obj.points))
                indPoint = obj.points(ind);
            else
                indPoint = AbstractGeometricPoint.empty(1,0);
            end
        end

        function [plotPoints, listBoxStr] = getPlottablePoints(obj)
            plotPoints = AbstractGeometricPoint.empty(1,0);
            listBoxStr = {};
            for(i=1:length(obj.points))
                if(obj.points(i).canBePlotted())
                    plotPoints(end+1) = obj.points(i); %#ok<AGROW>
                    listBoxStr{end+1} = obj.points(i).getListboxStr(); %#ok<AGROW>
                end
            end
        end

        function numPoints = getNumPoints(obj)
            numPoints = length(obj.points);
        end

        function tf = usesGroundObj(obj, groundObj)
            tf = false;
            for(i=1:length(obj.points))
                tf = tf || obj.points(i).usesGroundObj(groundObj);
            end
        end

        function tf = usesGeometricPoint(obj, point)
            tf = false;
            for(i=1:length(obj.points))
                tf = tf || obj.points(i).usesGeometricPoint(point);
            end
        end

        function tf = usesGeometricVector(obj, vector)
            tf = false;
            for(i=1:length(obj.points))
                tf = tf || obj.points(i).usesGeometricVector(vector);
            end
        end

        function tf = usesGeometricCoordSys(obj, coordSys)
            tf = false;
            for(i=1:length(obj.points))
                tf = tf || obj.points(i).usesGeometricCoordSys(coordSys);
            end
        end

        function tf = usesGeometricRefFrame(obj, refFrame)
            tf = false;
            for(i=1:length(obj.points))
                tf = tf || obj.points(i).usesGeometricRefFrame(refFrame);
            end
        end

        function tf = usesGeometricAngle(obj, angle)
            tf = false;
            for(i=1:length(obj.points))
                tf = tf || obj.points(i).usesGeometricAngle(angle);
            end
        end

        function tf = usesGeometricPlane(obj, plane)
            tf = false;
            for(i=1:length(obj.points))
                tf = tf || obj.points(i).usesGeometricPlane(plane);
            end
        end

        %graphical analysis tasks
        function ptGAStr = getAllPointGraphAnalysisTaskStrs(obj)
            ptGAStr = horzcat(obj.getPointPositionXGraphAnalysisTaskStrs(), ...
                              obj.getPointPositionYGraphAnalysisTaskStrs(), ...
                              obj.getPointPositionZGraphAnalysisTaskStrs(), ...
                              obj.getPointVelocityXGraphAnalysisTaskStrs(), ...
                              obj.getPointVelocityYGraphAnalysisTaskStrs(), ...
                              obj.getPointVelocityZGraphAnalysisTaskStrs());
        end

        function [ptPosXGAStr, points] = getPointPositionXGraphAnalysisTaskStrs(obj)
            points = obj.points;

            ptPosXGAStr = cell(1,length(points));
            A = length(points);
            formSpec = sprintf('%%0%uu',floor(log10(abs(A)+1)) + 1);
            for(i=1:length(points))
                ptPosXGAStr{i} = sprintf(sprintf('Point %s Position (X) - "%s"',formSpec, points(i).getName()), i);
            end
        end

        function [ptPosXGAStr, points] = getPointPositionYGraphAnalysisTaskStrs(obj)
            points = obj.points;

            ptPosXGAStr = cell(1,length(points));
            A = length(points);
            formSpec = sprintf('%%0%uu',floor(log10(abs(A)+1)) + 1);
            for(i=1:length(points))
                ptPosXGAStr{i} = sprintf(sprintf('Point %s Position (Y) - "%s"',formSpec, points(i).getName()), i);
            end
        end

        function [ptPosXGAStr, points] = getPointPositionZGraphAnalysisTaskStrs(obj)
            points = obj.points;

            ptPosXGAStr = cell(1,length(points));
            A = length(points);
            formSpec = sprintf('%%0%uu',floor(log10(abs(A)+1)) + 1);
            for(i=1:length(points))
                ptPosXGAStr{i} = sprintf(sprintf('Point %s Position (Z) - "%s"',formSpec, points(i).getName()), i);
            end
        end

        function [ptVelXGAStr, points] = getPointVelocityXGraphAnalysisTaskStrs(obj)
            points = obj.points;

            ptVelXGAStr = cell(1,length(points));
            A = length(points);
            formSpec = sprintf('%%0%uu',floor(log10(abs(A)+1)) + 1);
            for(i=1:length(points))
                ptVelXGAStr{i} = sprintf(sprintf('Point %s Velocity (X) - "%s"',formSpec, points(i).getName()), i);
            end
        end

        function [ptVelYGAStr, points] = getPointVelocityYGraphAnalysisTaskStrs(obj)
            points = obj.points;

            ptVelYGAStr = cell(1,length(points));
            A = length(points);
            formSpec = sprintf('%%0%uu',floor(log10(abs(A)+1)) + 1);
            for(i=1:length(points))
                ptVelYGAStr{i} = sprintf(sprintf('Point %s Velocity (Y) - "%s"',formSpec, points(i).getName()), i);
            end
        end

        function [ptVelZGAStr, points] = getPointVelocityZGraphAnalysisTaskStrs(obj)
            points = obj.points;

            ptVelZGAStr = cell(1,length(points));
            A = length(points);
            formSpec = sprintf('%%0%uu',floor(log10(abs(A)+1)) + 1);
            for(i=1:length(points))
                ptVelZGAStr{i} = sprintf(sprintf('Point %s Velocity (Z) - "%s"',formSpec, points(i).getName()), i);
            end
        end
    end
end


text:
%computes the spectral crest from the magnitude spectrum
%> called by ::ComputeFeature
%>
%> @param X: spectrogram (dimension FFTLength X Observations)
%> @param f_s: sample rate of audio data (unused)
%>
%> @retval vtsc spectral crest factor
% ======================================================================
function [vtsc] = FeatureSpectralCrestFactor (X, f_s)

    % get maximum
    vtsc = max(X, [], 1) ./ sum(X, 1);

    % avoid NaN for silence frames
    vtsc (sum(X, 1) == 0) = 0;
end

author: alexanderlerch

repo: ACA-Code


text:
function plot_subdivided(V,F,u,k)
%PLOT_SUBDIVIDED Plot the function u on the mesh V,F which has subdivided,
%with Loop subdivision, k times
%
% plot_subdivided(V,F,u,k);
%
% Inputs:
%  V,F  the input mesh to be subdivided
%  u  the input function to be subdivided along with the mesh
%  k  how many times to subdivide
%

Vu = ...
Fu = ...
uu = ...

t = tsurf(Fu,Vu, 'CData',uu);
shading interp;
axis equal;
axis off;
colormap(cbrewer('Blues', 500));
light('Position',[-1.5 1 1],'Style','local');
lights = camlight;
set(t, 'FaceLighting','gouraud', 'FaceColor','interp');
set(t, 'DiffuseStrength',0.5, 'SpecularStrength',0.2, 'AmbientStrength',0.3);
camproj('perspective');
add_shadow([t],lights);

end

author: odedstein

repo: sgi-introduction-course


text:
function SDP = relax_category_registration(problem,varargin)
%% Apply a sparse third-order relaxation to category registration
%% Depending on multivariate polynomial package in SPOT
%% Bounded translation is modelled as an inequality constraint, which leads
%% to an extra PSD block in the semidefinite relaxation
%% Heng Yang
%% July 05, 2021

params = inputParser;
params.CaseSensitive = false;

params.addParameter('checkMonomials',true, @(x) islogical(x));
params.addParameter('lambda',0.1, @(x) isscalar(x));

params.parse(varargin{:});

checkMonomials = params.Results.checkMonomials;
lambda         = params.Results.lambda;

fprintf('\n===================================================================')
fprintf('\nApplying SDP relaxation to category registration problem')
fprintf('\n===================================================================\n')
t0              = tic;

N               = problem.N;
K               = problem.K;
scene           = problem.scene;
shapes          = problem.shapes;
noiseBoundSq    = problem.noiseBoundSq;
tBound          = problem.translationBound;
tBoundSq        = tBound^2; % t'*t <= tBoundSq
cBoundSq        = problem.cBound^2; % should just be 1
barc2           = 1.0;

%% define POP variables
nrPrimalVars    = 9+3+K+N; % rotation: 9, translation: 3, binary: N, shape: K
p               = msspoly('p',nrPrimalVars);
r               = p(1:9);
R               = reshape(r,3,3);
col1 = R(:,1); col2 = R(:,2); col3 = R(:,3);
t               = p(10:12);
c               = p(12+1:12+K);
theta           = p(12+K+1:nrPrimalVars);

%% define cost function
shape           = combine_shapes(shapes,c);
residuals       = {};
for i = 1:N
    distance            = scene(:,i) - R * shape(:,i) - t;
    residuals{end+1}    = (distance' * distance) / noiseBoundSq;
end
f_cost = 0;
for i = 1:N
    f_cost = f_cost + (1+theta(i))/2 * residuals{i} + (1-theta(i))/2 * barc2;
end
% add regularization on c
f_cost = f_cost + lambda * (c'*c);

%% define constraints
h_r  = [1.0-col1'*col1;...
        1.0-col2'*col2;...
        1.0-col3'*col3;... % column unit length
        col1'*col2;...
        col2'*col3;...
        col3'*col1;... % colums orthogonal
        cross(col1,col2) - col3;...
        cross(col2,col3) - col1;...
        cross(col3,col1) - col2]; % columns righthandedness
% h_c  = [sum(c) - 1.0];

h_theta = [];
for i = 1:N
    h_theta =[h_theta; 1-theta(i)^2];
end

g_t = tBoundSq - t'*t; % Translation bounded
g_c = [cBoundSq - c'*c;c]; % nonnegative and bounded shape parameters

%% Formulate the sparse third-order relaxation
cr          = mykron(c,r);
x           = [r;t];
basis_p     = [1;x;c;cr;theta;mykron(theta,x);mykron(theta,cr)];
n           = length(basis_p);
basis_r     = get_multiplier_basis(p,basis_p,h_r(1));
basis_theta = get_multiplier_basis(p,basis_p,h_theta(1));
% basis_c     = get_multiplier_basis(p,basis_p,h_c(1));
basis_g_t   = [1;theta];
basis_g_t_f = mykron(basis_g_t,basis_g_t);
basis_g_c   = [1;r];
basis_g_c_f = mykron(basis_g_c,basis_g_c);
n1          = length(basis_g_t);
n2          = length(basis_g_c);
n1delta     = triangle_number(n1);
n2delta     = triangle_number(n2);

fprintf('Computing localizing and moment polynomials ...')
time_start  = tic;
pop = [mykron(basis_r,h_r);...
       mykron(basis_theta,h_theta);...
       mykron(basis_p,basis_p);...
       f_cost;...
       mykron(g_t,basis_g_t_f);...
       mykron(g_c,basis_g_c_f)];
[~,degmat,coef_all] = decomp(pop);
coef_all            = coef_all';
time_prep   = toc(time_start);
fprintf(' Done in %g seconds.\n',time_prep);

if checkMonomials
    fprintf('Checking consistency of monomials ...')
    time_check0   = tic;
    monomials_mom = mono(mykron(basis_p,basis_p));
    fprintf(' %d ... %d ...',size(degmat,1),length(monomials_mom));
    assert(size(degmat,1) == length(monomials_mom),'monomials not consistent');
    time_check    = toc(time_check0);
    fprintf('Done in %g seconds.\n',time_check);
end

dim_loc_eq      = length(basis_r)*length(h_r)+length(basis_theta)*length(h_theta);
dim_loc_ineq    = length(g_t)*n1delta+length(g_c)*n2delta;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% generate standard SDP data from degmat and coefficients
ndelta      = triangle_number(n);
nterms      = size(degmat,1);
m_mom       = ndelta - nterms;
m_loc       = dim_loc_eq;
m_loc_ineq  = dim_loc_ineq;
m           = m_mom + m_loc + m_loc_ineq + 1;

fprintf('SDP: n = %d, n1 = %d, m = %d, m_mom = %d, m_loc = %d, m_loc_ineq = %d, ndelta = %d.\n',...
        n,n1,m,m_mom,m_loc,m_loc_ineq,ndelta);

coef_mom    = coef_all(:,dim_loc_eq+1:dim_loc_eq+n^2);
coef_mom    = coef_mom';

B           = {};
B_normalize = {};
A           = {};

fprintf('Building B and A... Progress ')
for i = 1:nterms
    if rem(i,10000) == 1
        fprintf('%d/%d ',i,nterms);
    end

    [row,~,~]   = find(coef_mom(:,i));
    SDP_coli    = floor((row-1)./n) + 1;
    SDP_rowi    = mod(row-1,n) + 1;
    nnz         = length(SDP_rowi);

    Bi          = sparse(SDP_rowi,SDP_coli,ones(nnz,1),n,n);
    B{end+1}    = Bi;
    B_normalize{end+1} = Bi/nnz;

    mask_triu   = (SDP_rowi >= SDP_coli);
    si          = SDP_rowi(mask_triu);
    sj          = SDP_coli(mask_triu);

    nnz_triu    = length(si);

    if nnz_triu > 1
        [~,base_idx]        = max(sj);
        si_base             = si(base_idx);
        sj_base             = sj(base_idx);

        si_nonbase          = si;
        si_nonbase(base_idx)= [];
        sj_nonbase          = sj;
        sj_nonbase(base_idx)= [];

        is_base_diag        = (si_base == sj_base);

        if is_base_diag
            A_si            = [si_base];
            A_sj            = [sj_base];
            A_v             = [1];
        else
            A_si            = [si_base,sj_base];
            A_sj            = [sj_base,si_base];
            A_v             = [0.5,0.5];
        end

        for nonbase_idx = 1:length(si_nonbase)
            is_nonbase_diag = (si_nonbase(nonbase_idx) == sj_nonbase(nonbase_idx));
            if is_nonbase_diag
                A_sii       = [A_si,si_nonbase(nonbase_idx)];
                A_sjj       = [A_sj,sj_nonbase(nonbase_idx)];
                A_vv        = [A_v,-1];
            else
                A_sii       = [A_si,si_nonbase(nonbase_idx),sj_nonbase(nonbase_idx)];
                A_sjj       = [A_sj,sj_nonbase(nonbase_idx),si_nonbase(nonbase_idx)];
                A_vv        = [A_v,-0.5,-0.5];
            end
            A_temp          = sparse(A_sii,A_sjj,A_vv,n,n);

            A{end+1}        = A_temp;
        end
    end
end
fprintf('Done.\n')
assert(length(A) == m_mom,'length(A)+length(B) == ndelta!');

%% Now build A's associated with localizing constraints
if dim_loc_eq == 0
    % Do nothing
    A_local = {};
else
    coef_loc    = coef_all(:,1:dim_loc_eq);
    A_local     = {};
    fprintf('Building localizing constraints A_local... Progress ')
    for i = 1:dim_loc_eq
        if rem(i,10000) == 1
            fprintf('%d/%d ',i,m_loc);
        end

        [rowi,~,vi] = find(coef_loc(:,i));

        Ai      = sparse(n,n);
        for j   = 1:length(rowi)
            Ai  = Ai + vi(j) * B_normalize{rowi(j)};
        end
        A_local{end+1} = Ai;
    end
end
fprintf('Done.\n')

%% Leading A
A0 = sparse([1],[1],[1],n,n);

%% Combine all A for the main block
A = [{A0},A_local,A];

%% Now build the cost matrix
coef_cost   = coef_all(:,dim_loc_eq+n^2+1);
[row,~,v]   = find(coef_cost);
C           = sparse(n,n);
fprintf('Building cost matrix C... Progress ')
for i = 1:length(row)
    if rem(i,1000) == 1
        fprintf('%d/%d ',i,length(row));
    end
    C       = C + v(i) * B_normalize{row(i)};
end
fprintf('Done.\n')

%% Now build the sub PSD constraint g_t
coef_ineq   = coef_all(:,dim_loc_eq+n^2+1+1:dim_loc_eq+n^2+1+n1^2); % nterms by n1^2
A1          = {};
for ii = 1:n1^2
    row     = mod(ii-1,n1) + 1;
    col     = floor((ii-1)./n1) + 1;
    if row < col
        % Do nothing for upper triangular parts
    else
        if row == col
            A1tmp       = sparse([row],[col],[-1],n1,n1);
        else
            A1tmp       = sparse([row,col],[col,row],[-0.5,-0.5],n1,n1);
        end
        [termIds,~,v]   = find(coef_ineq(:,ii));
        Atmp            = sparse(n,n);
        for iii = 1:length(termIds)
            Atmp        = Atmp + v(iii) * B_normalize{termIds(iii)};
        end
        A{end+1}        = Atmp;
        A1{end+1}       = A1tmp;
    end
end

%% Now build the sub PSD constraint g_c (multiple of them)
Ac = {};
for k = 1:length(g_c)
    idx         = dim_loc_eq+n^2+1+n1^2+blkIndices(k,n2^2);
    coef_ineq   = coef_all(:,idx); % nterms by n2^2
    A2          = {};
    for ii = 1:n2^2
        row     = mod(ii-1,n2) + 1;
        col     = floor((ii-1)./n2) + 1;
        if row < col
            % Do nothing for upper triangular parts
        else
            if row == col
                A2tmp       = sparse([row],[col],[-1],n2,n2);
            else
                A2tmp       = sparse([row,col],[col,row],[-0.5,-0.5],n2,n2);
            end
            [termIds,~,v]   = find(coef_ineq(:,ii));
            Atmp            = sparse(n,n);
            for iii = 1:length(termIds)
                Atmp        = Atmp + v(iii) * B_normalize{termIds(iii)};
            end
            A{end+1}        = Atmp;
            A2{end+1}       = A2tmp;
        end
    end
    Ac{end+1} = A2;
end

assert(length(A) == m,'Total number of equality constraints wrong.')
count = length(A1);
for k = 1:length(g_c)
    count = count + length(Ac{k});
end
assert(count == m_loc_ineq,'Equality constraints from second PSD blk wrong')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
blk{1,1}        = 's';
blk{1,2}        = n;
blk{2,1}        = 's';
blk{2,2}        = n1;
for k = 1:length(g_c)
    blk{2+k,1}  = 's';
    blk{2+k,2}  = n2;
end
b               = sparse([1],[1],[1],m,1);
%% svec in sdpt3 format and output standard data
At0             = sparsesvec(blk(1,:),A);
At1             = [sparse(n1delta,m-m_loc_ineq),... % moment and localizing
                   sparsesvec(blk(2,:),A1),... % sub PSD t
                   sparse(n1delta,length(g_c)*n2delta)]; % sub PSD c
At              = {At0;At1};
for k = 1:length(g_c)
    Atc         = [sparse(n2delta,m-m_loc_ineq),... % moment and localizing
                   sparse(n2delta,n1delta),... % sub PSD t
                   sparse(n2delta,(k-1)*n2delta),... % sub PSD c
                   sparsesvec(blk(2+k,:),Ac{k}),... % sub PSD k-th c
                   sparse(n2delta,(length(g_c)-k)*n2delta)]; % sub PSD k-th c
    At          = [At;{Atc}];
end
C       = {C;sparse(n1,n1)};
for k = 1:length(g_c)
    C   = [C;{sparse(n2,n2)}];
end

SDP.blk = blk;
SDP.At  = At;
SDP.m   = m;
SDP.C   = C;
SDP.b   = b;
SDP.lam = lambda;
SDP.M   = 4+tBoundSq+cBoundSq+3*cBoundSq+N+(3+tBoundSq)*N+(3*cBoundSq)*N;

tf    = toc(t0);
fprintf('\nDone in %g seconds.\n',tf);
fprintf('===================================================================\n')
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
%  INTERNAL FUNCTION: computes chebyshev distances
%
%  ::
%
%    c=chebyshev_distance(y)
%
%  Args:
%
%     - **y** [numeric] : N x T x G array, with
%
%       - **N** [numeric] : number of simulations/replications
%       - **T** [numeric] : sample length (time series dimension)
%       - **G** [numeric] : number of variables
%
%  Returns:
%     :
%
%     - **c** [N x 1 vector] : chebyshev distances
%
%  See also:
%      - standardized_distance
%      - multivariate_chebyshev_box
%
%

author: jmaih


text:
classdef TestSuperEllipseExponent < handle

    properties (Access = private)
        dataFile = 'test_superEllipseExponent'
        m1
        m2
        qD
        q
    end

    methods (Access = public)

        function obj = TestSuperEllipseExponent()
            obj.loadData();
            obj.computeExponent();
        end

        function error = computeError(obj)
            error = norm(obj.qD - obj.q(:))/norm(obj.qD);
        end

    end

    methods (Access = private)

        function loadData(obj)
            d = load(obj.dataFile);
            obj.m1 = d.m1;
            obj.m2 = d.m2;
            obj.qD = d.q;
        end

        function computeExponent(obj)
            s.m1 = obj.m1;
            s.m2 = obj.m2;
            o = SmoothingExponentComputerOptimal(s);
            obj.q = o.compute();
        end

    end

end


text:
function [gist, param] = LMgist(D, HOMEIMAGES, param, HOMEGIST)
%
% param.imageSize = [256 256]; % it works also with non-square images (use the most common aspect ratio in your set)
% param.orientationsPerScale = [8 8 8 8]; % number of orientations per scale
% param.numberBlocks = 4;
% param.fc_prefilt = 4;
%
% [gist, param] = LMgist(D, HOMEIMAGES, param);
% [gist, param] = LMgist(filename, HOMEIMAGES, param);
% [gist, param] = LMgist(filename, HOMEIMAGES, param, HOMEGIST);
%
% For a one image or several images of the same size in a 4D array, use:
% gist = LMgist(img, [], param);
%
% When calling LMgist with a fourth argument it will store the gists in a
% new folder structure mirroring the folder structure of the images. Then,
% when called again, if the gist files already exist, it will just read
% them without recomputing them:
%
%   [gist, param] = LMgist(filename, HOMEIMAGES, param, HOMEGIST);
%   [gist, param] = LMgist(D, HOMEIMAGES, param, HOMEGIST);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Modeling the shape of the scene: a holistic representation of the spatial envelope
% Aude Oliva, Antonio Torralba
% International Journal of Computer Vision, Vol. 42(3): 145-175, 2001.

if nargin==4
    precomputed = 1;
    % get list of folders and create non-existing ones
    %listoffolders = {D(:).annotation.folder};

    %for i = 1:length(D);
    %    f{i} = D(i).annotation.folder;
    %end
    %[categories,b,class] = unique(f);
else
    precomputed = 0;
    HOMEGIST = '';
end

param.boundaryExtension = 32; % number of pixels to pad

if nargin<3
    % Default parameters
    param.imageSize = 128;
    param.orientationsPerScale = [8 8 8 8];
    param.numberBlocks = 4;
    param.fc_prefilt = 4;
    param.G = createGabor(param.orientationsPerScale, param.imageSize(1)+2*param.boundaryExtension);
else
    if ~isfield(param, 'G')
        param.G = createGabor(param.orientationsPerScale, param.imageSize(1)+2*param.boundaryExtension);
    end
end

% Precompute filter transfert functions (only need to do this once, unless
% image size is changes):
Nfeatures = size(param.G,3)*param.numberBlocks^2;

if isstruct(D)
    % [gist, param] = LMgist(D, HOMEIMAGES, param);
    Nscenes = length(D);
    typeD = 1;
end
if iscell(D)
    % [gist, param] = LMgist(filename, HOMEIMAGES, param);
    Nscenes = length(D);
    typeD = 2;
end
if isnumeric(D)
    % [gist, param] = LMgist(img, HOMEIMAGES, param);
    Nscenes = size(D,4);
    typeD = 3;
end

% Loop: Compute gist features for all scenes
gist = zeros([Nscenes Nfeatures], 'single');
for n = 1:Nscenes
    g = [];
    todo = 1;

    % if gist has already been computed, just read the file
    if precomputed==1
        filegist = fullfile(HOMEGIST, D(n).annotation.folder, [D(n).annotation.filename(1:end-4) '.mat']);
        if exist(filegist, 'file')
            load(filegist, 'g');
            todo = 0;
        end
    end

    % otherwise compute gist
    if todo==1
        if Nscenes>1
            disp([n Nscenes])
        end

        % load image
        try
            switch typeD
                case 1
                    img = LMimread(D, n, HOMEIMAGES);
                case 2
                    img = imread(fullfile(HOMEIMAGES, D{n}));
                case 3
                    img = D(:,:,:,n);
            end
        catch
            disp(D(n).annotation.folder)
            disp(D(n).annotation.filename)
            rethrow(lasterror)
        end

        % convert to gray scale
        img = single(mean(img,3));

        % resize and crop image to make it square
        img = imresizecrop(img, param.imageSize, 'bilinear');
        %img = imresize(img, param.imageSize, 'bilinear'); %jhhays

        % scale intensities to be in the range [0 255]
        img = img-min(img(:));
        img = 255*img/max(img(:));

        if Nscenes>1
            imshow(uint8(img))
            title(n)
        end

        % prefiltering: local contrast scaling
        output    = prefilt(img, param.fc_prefilt);

        % get gist:
        g = gistGabor(output, param);

        % save gist if a HOMEGIST file is provided
        if precomputed
            mkdir(fullfile(HOMEGIST, D(n).annotation.folder))
            save (filegist, 'g')
        end
    end

    gist(n,:) = g;
    drawnow
end

function output = prefilt(img, fc)
% ima = prefilt(img, fc);
% fc  = 4 (default)
%
% Input images are double in the range [0, 255];
% You can also input a block of images [ncols nrows 3 Nimages]
%
% For color images, normalization is done by dividing by the local
% luminance variance.

if nargin == 1
    fc = 4; % 4 cycles/image
end

w = 5;
s1 = fc/sqrt(log(2));

% Pad images to reduce boundary artifacts
img = log(img+1);
img = padarray(img, [w w], 'symmetric');
[sn, sm, c, N] = size(img);
n = max([sn sm]);
n = n + mod(n,2);
img = padarray(img, [n-sn n-sm], 'symmetric','post');

% Filter
[fx, fy] = meshgrid(-n/2:n/2-1);
gf = fftshift(exp(-(fx.^2+fy.^2)/(s1^2)));
gf = repmat(gf, [1 1 c N]);

% Whitening
output = img - real(ifft2(fft2(img).*gf));
clear img

% Local contrast normalization
localstd = repmat(sqrt(abs(ifft2(fft2(mean(output,3).^2).*gf(:,:,1,:)))), [1 1 c 1]);
output = output./(.2+localstd);

% Crop output to have same size than the input
output = output(w+1:sn-w, w+1:sm-w,:,:);

function g = gistGabor(img, param)
%
% Input:
%   img = input image (it can be a block: [nrows, ncols, c, Nimages])
%   param.w = number of windows (w*w)
%   param.G = precomputed transfer functions
%
% Output:
%   g: are the global features = [Nfeatures Nimages],
%                    Nfeatures = w*w*Nfilters*c

img = single(img);

w = param.numberBlocks;
G = param.G;
be = param.boundaryExtension;

if ndims(img)==2
    c = 1;
    N = 1;
    [nrows ncols c] = size(img);
end
if ndims(img)==3
    [nrows ncols c] = size(img);
    N = c;
end
if ndims(img)==4
    [nrows ncols c N] = size(img);
    img = reshape(img, [nrows ncols c*N]);
    N = c*N;
end

[ny nx Nfilters] = size(G);
W = w*w;
g = zeros([W*Nfilters N]);

% pad image
img = padarray(img, [be be], 'symmetric');

img = single(fft2(img));
k=0;
for n = 1:Nfilters
    ig = abs(ifft2(img.*repmat(G(:,:,n), [1 1 N])));
    ig = ig(be+1:ny-be, be+1:nx-be, :);

    v = downN(ig, w);
    g(k+1:k+W,:) = reshape(v, [W N]);
    k = k + W;
    drawnow
end

if c == 3
    % If the input was a color image, then reshape 'g' so that one column
    % is one images output:
    g = reshape(g, [size(g,1)*3 size(g,2)/3]);
end

function y=downN(x, N)
%
% averaging over non-overlapping square image blocks
%
% Input
%   x = [nrows ncols nchanels]
% Output
%   y = [N N nchanels]

nx = fix(linspace(0,size(x,1),N+1));
ny = fix(linspace(0,size(x,2),N+1));
y  = zeros(N, N, size(x,3));
for xx=1:N
  for yy=1:N
    v=mean(mean(x(nx(xx)+1:nx(xx+1), ny(yy)+1:ny(yy+1),:),1),2);
    y(xx,yy,:)=v(:);
  end
end

author: CSAILVision

repo: LabelMeToolbox


text:
function []=stvoldisp(beta_median,beta_std,beta_lbound,beta_ubound,sigma_median,sigma_t_median,sigma_t_lbound,sigma_t_ubound,gamma_median,X,Y,n,m,p,k,T,stvol,bex,ar,lambda1,lambda2,lambda3,lambda4,lambda5,gamma,alpha0,delta0,gamma0,zeta0,IRFt,const,endo,exo,startdate,enddate,stringdates1,decimaldates1,pref,PriorExcel)

% before displaying and saving the results, start estimating the evaluation measures for the model

% obtain first a point estimate betatilde of the VAR coefficients
% this is simply the median
betatilde=beta_median;
Btilde=reshape(betatilde,k,n);
% use this estimate to produce predicted values for the model, following (a.8.2)
Ytilde=X*Btilde;
% then produce the corresponding residuals, using (a.8.3)
EPStilde=Y-Ytilde;

% check first whether the model is stationary, using (a.7.2)
[stationary eigmodulus]=bear.checkstable(betatilde,n,p,k);

% Compute then the sum of squared residuals
% compute first the RSS matrix, defined in (a.8.4)
RSS=EPStilde'*EPStilde;
% retain only the diagonal elements to get the vector of RSSi values
rss=diag(RSS);

% Go on calculating R2
% generate Mbar
Mbar=eye(T)-ones(T,T)/T;
% then compute the TSS matrix, defined in (a.8.7)
TSS=Y'*Mbar*Y;
% generate the R2 matrix in (a.8.8)
R2=eye(n)-RSS./TSS;
% retain only the diagonal elements to get the vector of R2 values
r2=diag(R2);

% then calculate the adjusted R2, using (a.8.9)
R2bar=eye(n)-((T-1)/(T-k))*(eye(n)-R2);
% retain only the diagonal elements to get the vector of R2bar values
r2bar=diag(R2bar);

% now start displaying and saving the results

% preliminary task: create and open the txt file used to save the results

filelocation=fullfile(pref.results_path, [pref.results_sub '.txt']);
fid=fopen(filelocation,'wt');

% print toolbox header

fprintf('%s\n','');
fprintf(fid,'%s\n','');

% print the list of contributors
bear.printcontributors(fid);

% print then estimation results

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

toolboxinfo='BEAR toolbox estimates';
fprintf('%s\n',toolboxinfo);
fprintf(fid,'%s\n',toolboxinfo);

time=clock;
datestring=datestr(time);
dateinfo=['Date: ' datestring(1,1:11) '   Time: ' datestring(1,13:17)];
fprintf('%s\n',dateinfo);
fprintf(fid,'%s\n',dateinfo);

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

VARtypeinfo='Stochastic volatility BVAR';
fprintf('%s\n',VARtypeinfo);
fprintf(fid,'%s\n',VARtypeinfo);

if stvol==1
    modelinfo='Stochastic volatility model: standard';
elseif stvol==2
    modelinfo='Stochastic volatility model: random inertia';
elseif stvol==3
    modelinfo='Stochastic volatility model: large BVAR';
end
fprintf('%s\n',modelinfo);
fprintf(fid,'%s\n',modelinfo);

if IRFt==1
    SVARinfo='structural decomposition: none';
elseif IRFt==2
    SVARinfo='structural decomposition: choleski factorisation';
elseif IRFt==3
    SVARinfo='structural decomposition: triangular factorisation';
elseif IRFt==4
    SVARinfo='structural decomposition: sign restrictions';
end
fprintf('%s\n',SVARinfo);
fprintf(fid,'%s\n',SVARinfo);

temp='endogenous variables: ';
for ii=1:n
    temp=[temp ' ' endo{ii,1} ' '];
end
endoinfo=temp;
fprintf('%s\n',endoinfo);
fprintf(fid,'%s\n',endoinfo);

temp='exogenous variables: ';
if const==0 && m==0
    temp=[temp ' none'];
elseif const==1 && m==1
    temp=[temp ' constant '];
elseif const==0 && m>0
    for ii=1:m-1
        temp=[temp ' ' exo{ii,1} ' '];
    end
elseif const==1 && m>1
    temp=[temp ' constant '];
    for ii=1:m-1
        temp=[temp ' ' exo{ii,1} ' '];
    end
end
exoinfo=temp;
fprintf('%s\n',exoinfo);
fprintf(fid,'%s\n',exoinfo);

sampledateinfo=['estimation sample: ' startdate '-' enddate];
fprintf('%s\n',sampledateinfo);
fprintf(fid,'%s\n',sampledateinfo);

samplelengthinfo=['sample size (omitting initial conditions): ' num2str(T)];
fprintf('%s\n',samplelengthinfo);
fprintf(fid,'%s\n',samplelengthinfo);

laginfo=['number of lags included in regression: ' num2str(p)];
fprintf('%s\n',laginfo);
fprintf(fid,'%s\n',laginfo);

if PriorExcel==1
    arprint=[];
    for ii=1:n
        arprint=[arprint num2str(ar(ii,1)) '  '];
    end
    hyperparam2=['autoregressive coefficients (ar):                ' arprint];
else
    hyperparam2=['autoregressive coefficients (ar):                ' num2str(ar(1,1))];
end
fprintf('%s\n',hyperparam2);
fprintf(fid,'%s\n',hyperparam2);

hyperparam3=['overall tightness (lambda1):                    ' num2str(lambda1)];
fprintf('%s\n',hyperparam3);
fprintf(fid,'%s\n',hyperparam3);

if stvol==1||stvol==2
    hyperparam4=['cross-variable weighting (lambda2):             ' num2str(lambda2)];
    fprintf('%s\n',hyperparam4);
    fprintf(fid,'%s\n',hyperparam4);
end

hyperparam5=['lag decay (lambda3):                            ' num2str(lambda3)];
fprintf('%s\n',hyperparam5);
fprintf(fid,'%s\n',hyperparam5);

%hyperparam6=['exogenous variable tightness (lambda4):         ' num2str(lambda4)];
%fprintf('%s\n',hyperparam6);
%fprintf(fid,'%s\n',hyperparam6);

if bex==1
    hyperparam7=['block exogeneity shrinkage (lambda5):           ' num2str(lambda5)];
    fprintf('%s\n',hyperparam7);
    fprintf(fid,'%s\n',hyperparam7);
end

if stvol==1||stvol==3
    hyperparam8=['AR coefficient on residual variance (gamma):    ' num2str(gamma)];
    fprintf('%s\n',hyperparam8);
    fprintf(fid,'%s\n',hyperparam8);
end

hyperparam9=['IG shape on residual variance (alpha0):         ' num2str(alpha0)];
fprintf('%s\n',hyperparam9);
fprintf(fid,'%s\n',hyperparam9);

hyperparam10=['IG scale on residual variance (delta0):         ' num2str(alpha0)];
fprintf('%s\n',hyperparam10);
fprintf(fid,'%s\n',hyperparam10);

if stvol==2
    hyperparam11=['Prior mean on inertia (gamma0):                 ' num2str(gamma)];
    fprintf('%s\n',hyperparam11);
    fprintf(fid,'%s\n',hyperparam11);
end

if stvol==2
    hyperparam12=['Prior variance on inertia (zeta0):              ' num2str(gamma)];
    fprintf('%s\n',hyperparam12);
    fprintf(fid,'%s\n',hyperparam12);
end

% display coefficient estimates

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

coeffinfo=['VAR coefficients (beta): posterior estimates'];
fprintf('%s\n',coeffinfo);
fprintf(fid,'%s\n',coeffinfo);

for ii=1:n

    fprintf('%s\n','');
    fprintf(fid,'%s\n','');
    if ii~=1
        fprintf('%s\n','');
        fprintf(fid,'%s\n','');
    end

    endoinfo=['Endogenous: ' endo{ii,1}];
    fprintf('%s\n',endoinfo);
    fprintf(fid,'%s\n',endoinfo);

    coeffheader=fprintf('%25s %15s %15s %15s %15s\n','','Median','St.dev','Low.bound','Upp.bound');
    coeffheader=fprintf(fid,'%25s %15s %15s %15s %15s\n','','Median','St.dev','Low.bound','Upp.bound');

    % handle the endogenous
    for jj=1:n
        for kk=1:p
            values=[beta_median((ii-1)*k+n*(kk-1)+jj,1) beta_std((ii-1)*k+n*(kk-1)+jj,1) beta_lbound((ii-1)*k+n*(kk-1)+jj,1) beta_ubound((ii-1)*k+n*(kk-1)+jj,1)];
            fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n',strcat(endo{jj,1},'(-',int2str(kk),')'),values);
            fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n',strcat(endo{jj,1},'(-',int2str(kk),')'),values);
        end
    end

    % handle the exogenous
    % if there is no constant:
    if const==0
        % if there is no exogenous at all, obvioulsy, don't display anything
        if m==0
            % if there is no constant but some other exogenous, display them
        else
            for jj=1:m
                values=[beta_median(ii*k-m+jj,1) beta_std(ii*k-m+jj,1) beta_lbound(ii*k-m+jj,1) beta_ubound(ii*k-m+jj,1)];
                fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
                fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
            end
        end
        % if there is a constant
    else
        % display the results related to the constant
        values=[beta_median(ii*k-m+1,1) beta_std(ii*k-m+1,1) beta_lbound(ii*k-m+1,1) beta_ubound(ii*k-m+1,1)];
        fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n','Constant',values);
        fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n','Constant',values);
        % if there is no other exogenous, stop here
        if m==1
            % if there are other exogenous, display their results
        else
            for jj=1:m-1
                values=[beta_median(ii*k-m+jj+1,1) beta_std(ii*k-m+jj+1,1) beta_lbound(ii*k-m+jj+1,1) beta_ubound(ii*k-m+jj+1,1)];
                fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
                fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
            end
        end
    end

    fprintf('%s\n','');
    fprintf(fid,'%s\n','');

    % display evaluation measures
    rssinfo=['Sum of squared residuals: ' num2str(rss(ii,1),'%.2f')];
    fprintf('%s\n',rssinfo);
    fprintf(fid,'%s\n',rssinfo);

    r2info=['R-squared: ' num2str(r2(ii,1),'%.3f')];
    fprintf('%s\n',r2info);
    fprintf(fid,'%s\n',r2info);

    adjr2info=['adj. R-squared: ' num2str(r2bar(ii,1),'%.3f')];
    fprintf('%s\n',adjr2info);
    fprintf(fid,'%s\n',adjr2info);

end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display VAR stability results
eigmodulus=reshape(eigmodulus,p,n);
stabilityinfo1=['Roots of the characteristic polynomial (modulus):'];
fprintf('%s\n',stabilityinfo1);
fprintf(fid,'%s\n',stabilityinfo1);
for ii=1:p
    temp=num2str(eigmodulus(ii,1),'%.3f');
    for jj=2:n
        temp=[temp,'  ',num2str(eigmodulus(ii,jj),'%.3f')];
    end
    fprintf('%s\n',temp);
    fprintf(fid,'%s\n',temp);
end
if stationary==1;
    stabilityinfo2=['No root lies outside the unit circle.'];
    stabilityinfo3=['The estimated VAR model satisfies the stability condition'];
    fprintf('%s\n',stabilityinfo2);
    fprintf(fid,'%s\n',stabilityinfo2);
    fprintf('%s\n',stabilityinfo3);
    fprintf(fid,'%s\n',stabilityinfo3);
else
    stabilityinfo2=['Warning: at leat one root lies on or outside the unit circle.'];
    stabilityinfo3=['The estimated VAR model will not be stable'];
    fprintf('%s\n',stabilityinfo2);
    fprintf(fid,'%s\n',stabilityinfo2);
    fprintf('%s\n',stabilityinfo3);
    fprintf(fid,'%s\n',stabilityinfo3);
end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display posterior for sigma
sigmainfo1=['sigma (residual covariance matrix): posterior estimates'];
fprintf('%s\n',sigmainfo1);
fprintf(fid,'%s\n',sigmainfo1);
% calculate the (integer) length of the largest number in sigma, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(sigma_median))))));
% add a separator, a potential minus sign, and three digits (total=5) to obtain the total space for each entry in the matrix
width=width+5;
for ii=1:n
    temp=[];
    for jj=1:n
        % convert matrix entry into string
        number=num2str(sigma_median(ii,jj),'% .3f');
        % pad potential missing blanks
        while numel(number)<width
            number=[' ' number];
        end
        number=[number '  '];
        temp=[temp number];
    end
    fprintf('%s\n',temp);
    fprintf(fid,'%s\n',temp);
end
% notice about use of long-run values
sigmainfo2=['Note: estimates rely on long-run (homoskedastic) values'];
fprintf('%s\n',sigmainfo2);
fprintf(fid,'%s\n',sigmainfo2);

% display posterior for gamma (if random intertia model)
if stvol==2
    fprintf('%s\n','');
    fprintf(fid,'%s\n','');
    fprintf('%s\n','');
    fprintf(fid,'%s\n','');
    % display posterior for gamma
    gammainfo1=['gamma (autoregressive coefficient on stochastic volatility): posterior estimates'];
    fprintf('%s\n',gammainfo1);
    %fprintf(fid,'%s\n',gammainfo1);
    gammainfo2=[];
    for ii=1:n
        gammainfo2=[gammainfo2 repmat(' ',1,5-numel(ii)) 'gamma' num2str(ii)];
    end
    fprintf('%s\n',gammainfo2);
    %fprintf(fid,'%s\n',gammainfo2);
    gammainfo3=[];
    for ii=1:n
        % convert matrix entry into string
        temp=num2str(gamma_median(1,ii),'% .3f');
        % pad potential missing blanks and complete
        temp=[repmat(' ',1,10-numel(temp)) temp];
        gammainfo3=[gammainfo3 temp];
    end
    fprintf('%s\n',gammainfo3);
    %fprintf(fid,'%s\n',gammainfo3);
end
fclose(fid);

% Finally, display the results in terms of graph
if pref.plot
    % plot actual vs. fitted
    actualfitted=figure('Tag','BEARresults');
    set(actualfitted,'Color',[0.9 0.9 0.9]);
    set(actualfitted,'name','model estimation: actual vs fitted')
    ncolumns=ceil(n^0.5);
    nrows=ceil(n/ncolumns);
    for ii=1:n
        subplot(nrows,ncolumns,ii)
        hold on
        plot(decimaldates1,Y(:,ii),'Color',[0 0 0],'LineWidth',2);
        plot(decimaldates1,Ytilde(:,ii),'Color',[1 0 0],'LineWidth',2);
        hold off
        set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
        title(endo{ii,1},'FontName','Times New Roman','FontSize',10,'FontWeight','normal');
        if ii==1
            plotlegend=legend('actual','fitted');
            set(plotlegend,'FontName','Times New Roman');
        end
    end

    % plot the residuals
    residuals=figure('Tag','BEARresults');
    set(residuals,'Color',[0.9 0.9 0.9]);
    set(residuals,'name','model estimation: residuals')
    for ii=1:n
        subplot(nrows,ncolumns,ii)
        plot(decimaldates1,EPStilde(:,ii),'Color',[0 0 0],'LineWidth',2)
        set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
        title(endo{ii,1},'FontName','Times New Roman','FontSize',10,'FontWeight','normal');
    end

    % plot first the time-varying variance and covariance estimates
    varcov=figure('Tag','BEARresults');
    set(varcov,'Color',[0.9 0.9 0.9]);
    set(varcov,'name','model estimation: residual variance and covariance')
    for ii=1:n
        for jj=1:ii
            subplot(n,n,n*(ii-1)+jj)
            hold on
            Xpatch=[decimaldates1' fliplr(decimaldates1')];
            Ypatch=[sigma_t_lbound{ii,jj}' fliplr(sigma_t_ubound{ii,jj}')];
            HDpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
            set(HDpatch,'facealpha',0.6);
            set(HDpatch,'edgecolor','none');
            plot(decimaldates1,sigma_t_median{ii,jj},'Color',[0.4 0.4 1],'LineWidth',2);
            plot([decimaldates1(1,1),decimaldates1(end,1)],[0 0],'k--');
            hold off
            set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
            % top labels
            if jj==ii
                title(['var(' endo{ii,1} ')'],'FontWeight','normal');
            else
                title(['cov(' endo{jj,1} ',' endo{ii,1} ')'],'FontWeight','normal');
            end
        end
    end

end % pref.plot

% finally, save the results on excel

% compute the cell for actual and fitted
% create the cell that will be saved on excel
afcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},T+3,1);
% loop over variables (horizontal dimension)
for ii=1:n
    % create cell of actual/fitted for variable ii
    temp=['actual and fitted: ' endo{ii,1}];
    af_i=[temp {''} {''};{''} {''} {''};{''} {'sample'} {'fitted'};stringdates1 num2cell(Y(:,ii)) num2cell(Ytilde(:,ii))];
    afcell=[afcell af_i vertspace];
end
% trim
afcell=afcell(:,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),afcell,'actual fitted','B2');
end

% then compute the cell for the residuals
% create the cell that will be saved on excel
horzspace=repmat({''},1,n);
rescell=[{'residuals'} horzspace;{''} horzspace;{''} endo';stringdates1 num2cell(EPStilde)];
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),rescell,'resids','B2');
end

% finally compute the cell for the time varying variance and covariance
% create the cell that will be saved on excel
varcovcell={};
% build preliminary elements: space between the tables
horzspace=repmat({''},2,5*n);
vertspace=repmat({''},T+3,1);
% loop over variables (vertical dimension)
for ii=1:n
    tempcell={};
    % loop over shocks (horizontal dimension)
    for jj=1:ii
        % create cell of hd record for the contribution of shock jj in variable ii fluctuation
        % if a sign restriction identification scheme has been used, use the structural shock labels
        if jj==ii
            temp=['variance of ' endo{ii,1} ' residuals'];
            % otherwise, the shocks are just orthogonalised shocks from the variables: use variable names
        else
            temp=['covariance between ' endo{jj,1} ' and ' endo{ii,1} 'residuals'];
        end
        vc_ij=[temp {''} {''} {''};{''} {''} {''} {''};{''} {'lw. bound'} {'median'} {'up. bound'};stringdates1 num2cell(sigma_t_lbound{ii,jj}) num2cell(sigma_t_median{ii,jj}) num2cell(sigma_t_ubound{ii,jj})];
        tempcell=[tempcell vc_ij vertspace];
    end
    % complete with blanks (for repeated covariance values)
    for jj=1:n-ii
        tempcell=[tempcell cell(T+3,5)];
    end
    varcovcell=[varcovcell;horzspace;tempcell];
end
% trim
varcovcell=varcovcell(3:end,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),varcovcell,'time variation','B2');
end

author: european-central-bank

repo: BEAR-toolbox


text:
function rtk=udpos_rtkins(rtk,tt) %#ok

ins=rtk.ins;
rtk.x(1:15)=0;
rtk.P(1:15,1:15)=0;

% lever arm correction
pos=ins.pos+ins.Mpv*ins.Cnb*ins.lever;
vel=ins.vel+ins.Cnb*askew(ins.web)*ins.lever;

rtk.x(1:15)=[ins.att;vel;pos;ins.bg;ins.ba];
rtk.P(1:15,1:15)=ins.P;

return

repo: GINav


text:
function [V]=quad_smooth(V,Ev,IND_V_not_Ev,Lc,n,Ls,w1)

L1=Ev>0;
L2=IND_V_not_Ev>0;
w2=1-w1;

p=V;

for i=1:n
    Xp=NaN(size(Ev));  Yp=NaN(size(Ev));  Zp=NaN(size(Ev));
    Xp(L1)=p(Ev(L1),1); Yp(L1)=p(Ev(L1),2); Zp(L1)=p(Ev(L1),3);
    Vp1=[gnanmean(Xp,2) gnanmean(Yp,2) gnanmean(Zp,2)];

    Xp=NaN(size(Ev));  Yp=NaN(size(Ev));  Zp=NaN(size(Ev));
    Xp(L2)=p(IND_V_not_Ev(L2),1); Yp(L2)=p(IND_V_not_Ev(L2),2); Zp(L2)=p(IND_V_not_Ev(L2),3);
    Vp2=[gnanmean(Xp,2) gnanmean(Yp,2) gnanmean(Zp,2)];

    Vp=(w1.*Vp1+w2.*Vp2)./(w1+w2);
    p=p+Ls.*(Vp-p);

    p(Lc,:)=V(Lc,:);
    if i>1
        D=gnansum((p-p_old).^2,2);
        SSQD_new=gnansum(D(:));
        if i>2
            SSQD_ratio=SSQD_new./SSQD_old;
            disp(num2str(SSQD_ratio));
        end
        SSQD_old=SSQD_new;
    end
    p_old=p;

end
V=p;

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function slice = imMiddleSlice(img)
%IMMIDDLESLICE Extract the middle slice of a 3D stack
%
%   SLICE = imMiddleSlice(STACK)
%
%   Example
%   imMiddleSlice
%
%   See also
%   imStacks, stackSlice

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2018-06-05,    using Matlab 9.4.0.813654 (R2018a)
% Copyright 2018 INRA - Cepia Software Platform.

switch ndims(img)
    case 3
        slice = img(:,:,round(size(img, 3) / 2));
    case 4
        slice = img(:,:,:,round(size(img, 4) / 2));
    otherwise
        error('Requires an input image with 3 or 4 dimensions');
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function imdata = mcmcComputeImageData(im, imsegs)

[imh, imw] = size(imsegs.segimage);
imdata.yim = 1-repmat([(0:imh-1)/(imh-1)]', 1, imw);
imdata.xim = repmat([(0:imw-1)/(imw-1)], imh, 1);

[gx, gy] = gradient(im);
imdata.gradim = sqrt(sum(gx.^2 + gy.^2, 3));

minEdgeLen = sqrt(imh^2+imw^2)*0.02;
[vpdata.lines, vpdata.spinfo] = ...
    APPgetLargeConnectedEdges(rgb2gray(im), minEdgeLen, imsegs);
   [vpdata.v, vpdata.vars, vpdata.p, vpdata.hpos] = ...
       APPestimateVp(vpdata.lines, [imh imw], 0);
%  [vpdata.v vpdata.vars vpdata.p vpdata.hpos]=AppgetVP(vpdata.lines,[imh imw],0);

imdata.vpdata = vpdata;

% get pixels in each superpixel
stats = regionprops(imsegs.segimage, 'PixelIdxList');
imdata.pixlist = {stats.PixelIdxList};

%%by varsha%%%

for segno=1:imsegs.nseg
   BW=(imsegs.segimage==segno);
   imdata.tracedbndy{segno}= bwboundaries(BW);
end

%%%%%%%%%%%
% get boundary pixels for each superpixel
[tx, ty] = gradient(double(imsegs.segimage));
segedgeim = (tx~=0) | (ty ~=0);
stats = regionprops(double(imsegs.segimage).*segedgeim, 'PixelIdxList');
imdata.bndpixlist = {stats.PixelIdxList};
imdata.bndnpix = zeros(size(imdata.bndpixlist));
for s = 1:imsegs.nseg
    try
        imdata.bndnpix(s) = numel(imdata.bndpixlist{s});
    catch
        imdata.bndnpix(s) = 0;
    end
end

% get a random subset of pixels and boundary pixels to speed computation of
% some features
for s = 1:imsegs.nseg
    npix = imsegs.npixels(s);
    nbndpix = imdata.bndnpix(s);
    imdata.nsmpix(s) = min(npix, 1000);
    imdata.nsmbndpix(s) = floor(nbndpix/4);
    if npix>1000
        rind = randperm(npix);
        imdata.smpixlist{s} = imdata.pixlist{s}(rind(1:1000));
    else
        imdata.smpixlist(s) = imdata.pixlist(s);
    end
    rind = randperm(nbndpix);
    try
        imdata.smbndpixlist{s} = imdata.bndpixlist{s}(rind(1:imdata.nsmbndpix(s)));
    catch
        imdata.smbndpixlist{s} = [];
    end
end

author: zouchuhang

repo: LayoutNet


text:
classdef RodAlign < dagnn.Layer

    properties
        subdivisions = [6 6]
        transformation = 1/8
    end

    methods
        function outputs = forward(obj, inputs, params)
          % calculation on cpu is faster for matlab
          inputs{1} = gather(inputs{1});
          inputs{2} = gather(inputs{2});

          feat_h = size(inputs{1},1);
          feat_w = size(inputs{1},2);
          num_roi = size(inputs{2},2);
          roi_h_start = max(round(inputs{2}(2,:) * obj.transformation),1);
          roi_w_start = max(round(inputs{2}(3,:) * obj.transformation),1);
          roi_h_end = min(round(inputs{2}(4,:) * obj.transformation),feat_h);
          roi_w_end = min(round(inputs{2}(5,:) * obj.transformation),feat_w);

          scale_h = feat_h / obj.subdivisions(1);
          scale_w = feat_w / obj.subdivisions(2);
          bin_h = ((1:obj.subdivisions(1))-0.5) * scale_h;
          bin_W = ((1:obj.subdivisions(2))-0.5) * scale_w;

          floor_bin_h = max(floor(bin_h),1);
          ceil_bin_h = floor_bin_h+1;
          floor_bin_w = max(floor(bin_W),1);
          ceil_bin_w = floor_bin_w+1;
          c11 = min(ceil_bin_h-bin_h,1)'*min(ceil_bin_w-bin_W,1);
          c12 = min(ceil_bin_h-bin_h,1)'*max(bin_W-floor_bin_w,0);
          c21 = max((bin_h-floor_bin_h),0)'*min(ceil_bin_w-bin_W,1);
          c22 = max(bin_h-floor_bin_h,0)'*max(bin_W-floor_bin_w,0);
          feats = repmat(inputs{1},1,1,1,num_roi);
          for i = 1:num_roi
              feats(roi_h_start(i):roi_h_end(i),roi_w_start(i):roi_w_end(i),:,i) = 0;
          end
          outputs = c11.*feats(floor_bin_h,floor_bin_w,:,:) ...
                  + c12.*feats(floor_bin_h,ceil_bin_w,:,:) ...
                  + c21.*feats(ceil_bin_h,floor_bin_w,:,:) ...
                  + c22.*feats(ceil_bin_h,ceil_bin_w,:,:);
          outputs = {gpuArray(outputs)};

        end

        function [derInputs, derParams] = backward(obj, inputs, param, derOutputs)
          %When assigning into a GPUArray, the subscripts must contain unique values.
          inputs{1} = gather(inputs{1});
          inputs{2} = gather(inputs{2});
          derOutputs{1} = gather(derOutputs{1});

          feat_h = size(inputs{1},1);
          feat_w = size(inputs{1},2);
          ch = size(inputs{1},3);
          num_roi = size(inputs{2},2);
          roi_h_start = max(round(inputs{2}(2,:) * obj.transformation),1);
          roi_w_start = max(round(inputs{2}(3,:) * obj.transformation),1);
          roi_h_end = min(round(inputs{2}(4,:) * obj.transformation),feat_h);
          roi_w_end = min(round(inputs{2}(5,:) * obj.transformation),feat_w);

          scale_h = feat_h / obj.subdivisions(1);
          scale_w = feat_w / obj.subdivisions(2);
          bin_h = ((1:obj.subdivisions(1))-0.5) * scale_h;
          bin_W = ((1:obj.subdivisions(2))-0.5) * scale_w;

          floor_bin_h = max(floor(bin_h),1);
          ceil_bin_h = floor_bin_h+1;
          floor_bin_w = max(floor(bin_W),1);
          ceil_bin_w = floor_bin_w+1;
          c11 = min(ceil_bin_h-bin_h,1)'*min(ceil_bin_w-bin_W,1);
          c12 = min(ceil_bin_h-bin_h,1)'*max(bin_W-floor_bin_w,0);
          c21 = max((bin_h-floor_bin_h),0)'*min(ceil_bin_w-bin_W,1);
          c22 = max(bin_h-floor_bin_h,0)'*max(bin_W-floor_bin_w,0);

          derbilinear = zeros([feat_h,feat_w,ch,num_roi],'single');
          for i = 1:num_roi
              derbilinear(floor_bin_h,floor_bin_w,:,i) = derbilinear(floor_bin_h,floor_bin_w,:,i) + c11 .* derOutputs{1}(:,:,:,i);
              derbilinear(floor_bin_h,ceil_bin_w,:,i) = derbilinear(floor_bin_h,ceil_bin_w,:,i) + c12 .* derOutputs{1}(:,:,:,i);
              derbilinear(ceil_bin_h,floor_bin_w,:,i) = derbilinear(ceil_bin_h,floor_bin_w,:,i) + c21 .* derOutputs{1}(:,:,:,i);
              derbilinear(ceil_bin_h,ceil_bin_w,:,i) = derbilinear(ceil_bin_h,ceil_bin_w,:,i) + c22 .* derOutputs{1}(:,:,:,i);
              derbilinear(roi_h_start(i):roi_h_end(i),roi_w_start(i):roi_w_end(i),:,i) = 0;
          end
          derInputs{1} = gpuArray(sum(derbilinear,4));

          derInputs{2} = [];
          derParams = {} ;
        end

        function obj = RodAlign(varargin)
          obj.load(varargin);
        end
    end
end

author: HuiZeng

repo: Grid-Anchor-based-Image-Cropping


text:
function [uOutput] = yourfilterimp(nFunction, sFilename)
% Import filter template -- copy and modify to import your own delimted ascii files
%
% To use this:
%   - Copy it into an appropriate new name (still in the importfilters directory)
%   - Change it

% Assume this function is to convert a comma-separated file formatted thusly:
%
%  ID, La, Lo, Dp, Yr, Mo, Dy, Tm, Ma, MagT, cm
%  01, 46.0, 10.3, 4.0, 2015, 12, 3, 3.1, 12:00:03, Mb, comment
%  02, 45.2, 6.4, 6.0, 2015, 12, 4, 3.2, 03:00:03, Mb, comment
%  03, 45.8, 9.1, 8.0, 2015, 12, 3, 2.0, 12:00:03, Mb, comment
%   ... etc ...
%
%  The goal is to turn this into a ZmapCatalog.
%

% if the file did not have a header, you would specify the column names here:
%
% ColumnNamesInFile = {'ID','Latitude','Longitude','Depth','Year','Month','Day','Time','Magnitude','MagnitudeType', 'Comment'}

% if the names need to be remapped, create the mapping here:
%   - note, if a column is not mentioned, it will be unchanged
%   - note, if a column is mapped to [], it will be removed
fileCols2catalogCols = {'La', 'Latitude';...
                        'Lo', 'Longitude';...
                        'Dp', 'Depth';...
                        'Ma', 'Magnitude';...
                        'MagT', 'MagnitudeType';...
                        'cm', []...
                        };

% Filter function switchyard
switch nFunction
    case FilterOp.getDescription
        % return a short description which will appear in the drop-down list
        uOutput = 'Your data format - adjust the file yourfilterimp.m';
    case FilterOp.getWebpage
        % return a web-page detailing your format
        uOutput = 'yourfilterimp.html';
    case FilterOp.importCatalog

        opts = detectImportOptions(sFilename);
        tb = readtable(sFilename, opts);

        % now add or change column names, thesemay not map directly to ZmapCatalog yet
        if exist('ColumnNamesInFile','var')
            tb.Properties.VariableNames = ColumnNamesInFile;
        end

        % now calculate any columns as necessary.
        tb.Date = datetime(tb.Yr, tb.Mo, tb.Dy) + tb.Tm;  % in this case, tm was automatically converted to a duration

        % now rename and delete columns as desired
        if exist('fileCols2catalogCols','var')
            remapper = containers.Map(fileCols2catalogCols(:,1),fileCols2catalogCols(:,2));
            for idx = numel(tb.Properties.VariableNames) : -1 : 1 % going backwards avoids trouble deleting
                fcolname = tb.Properties.VariableNames{idx};
                if ismember(fcolname, remapper.keys())
                    if isempty(remapper(fcolname))
                        tb.(fcolname) = [];
                    else
                        tb.Properties.VariableNames{idx} = remapper(fcolname);
                    end
                end
            end
        end
        uOutput = ZmapCatalog.from(tb);
        return

end

author: CelsoReyes


text:
%Copyright (c) October,15 2008 by Varsha Hedau, UIUC.  All rights reserved.
function Xpnts = ComputeIntersectionPoints(lines)

%%%%%%Computing intersections of all the lines%%%%%%
p1 = [lines(:, [1 3]) ones(size(lines, 1), 1)];
p2 = [lines(:, [2 4]) ones(size(lines, 1), 1)];
% get plane normals for line segments
l = cross(p1, p2);
l = l ./ repmat(sqrt(sum(l.^2,2)), 1, 3);

[XX YY]=meshgrid(1:size(l,1));
ll1=l(XX(:),:);ll2=l(YY(:),:);
Xpnts=cross(ll1,ll2);

%[x1 y1 x2 y2 x3 y3] are colinear if x1(y2-y3)+x2(y3-y1)+x3(y1-y2)=0;
colchck=[lines(XX(:),1) lines(XX(:),3) lines(YY(:),1) lines(YY(:),3) lines(YY(:),2) lines(YY(:),4)];
colchck=colchck(:,1).*(colchck(:,4)-colchck(:,6))+colchck(:,3).*(colchck(:,6)-colchck(:,2))+...
    colchck(:,5).*(colchck(:,2)-colchck(:,4));

keepind=find(abs(colchck)>50);
Xpnts=Xpnts(keepind,:);
Xpnts=[Xpnts(:,1)./Xpnts(:,3) Xpnts(:,2)./Xpnts(:,3)];


text:
function G = filterDog2d( r, var, order, show )
% Difference of Gaussian (Dog) Filter.
%
% Adapted from code by Serge Belongie.  Takes a "Difference of Gaussian" -
% all centered on the same point but with different values for sigma. Also
% serves as an approximation to an Laplacian of Gaussian (LoG) filter (if
% order==1).
%
% USAGE
%  G = filterDog2d( r, var, order, [show] )
%
% INPUTS
%  r       - Final filter will be 2*r+1 on each side
%  var     - variance of central Gaussian
%  order   - should be either 1-LoG or 2-difference of 3 Gaussians
%  show    - [0] figure to use for optional display
%
% OUTPUTS
%  G       - filter
%
% EXAMPLE
%  G = filterDog2d( 40, 40, 1, 1 ); %order=1 (LoG)
%  G = filterDog2d( 40, 40, 2, 3 ); %order=2
%
% See also FILTERDOOG, FILTERGAUSS
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<4 || isempty(show) ); show=0; end

% create filter
N = 2*r+1;
if (order==1)
  Ga = filterGauss( [N N], [], .71*var );
  Gb = filterGauss( [N N], [], 1.14*var );
  a=1; b=-1; G = a*Ga + b*Gb;

elseif (order==2)
  Ga = filterGauss( [N N], [], 0.62*var );
  Gb = filterGauss( [N N], [], var );
  Gc = filterGauss( [N N], [], 1.6*var );
  a=-1; b=2; c=-1; G = a*Ga + b*Gb + c*Gc;

else
  error('order must be 1 or 2');
end

% normalize
G=G-mean(G(:));
G=G/norm(G(:),1);

% display
if(show); filterVisualize( G, show, 'row' ); end


text:
function ge_test1
% This is an example script on how to get data (loacations with altitude, colored image) from GoogleEarth into
% a matlab surface.
% Use this function as a script so comment the first line!
%
% The code was tested with Matlab 2010a an GoogleEarth 6 (actual version) so download it or try your actual
% GE-version:       http://www.google.com/intl/en/earth/index.html
%
% Further help for GoogleEarthAPI you find under:
% http://earth.google.com/comapi/interfaceIApplicationGE.html#9e1c8da5b36e8687fe718f2267224103
%
% You need additional MCFE-files before script will run without errors:
%       http://www.mathworks.com/matlabcentral/fileexchange/5256-pos2dist
%       http://www.mathworks.com/matlabcentral/fileexchange/28357-jmouseemu-mouse-emulator-v2-2
%       http://www.mathworks.com/matlabcentral/fileexchange/28603
%       http://www.mathworks.com/matlabcentral/fileexchange/14584-clipboardimage
%
%  ATTENTION: There is a heavy workaround to get an colored image from GE in line 70 (jmouseemu). Here it is
%  necessery to set the mousepointer in the actual GE-picture, so the x and y pixels must be in it! The
%  problem is, that it will not work in debug mode.
%  Depending of internetconnection speed you must play a little with the pause-times to get all needed data.
%
% Programmed by Sven Koerner: koerner(underline)sven(add)gmx.de
% Date: 2011/07/27

% open GR
ge_app = actxserver('GoogleEarth.ApplicationGE');   % need GoogleEarth tested with version 6

% wait while initialising
pause(7)

% get camerapostion
camPos = ge_app.GetCamera(0);

% teleport to location
camPos.FocusPointLatitude           = 51.10719595658241;
camPos.FocusPointLongitude          = 13.638689585783112;
camPos.FocusPointAltitude           = 0.0;
camPos.FocusPointAltitudeMode       = 'AbsoluteAltitudeGE';
camPos.Range                        = 11780.53549488341;
camPos.Tilt                         = 0;
camPos.Azimuth                      = 7;  % 0 < Speed <5 --> Speed> 5: Teleport
ge_app.SetCameraParams(camPos.FocusPointLatitude, camPos.FocusPointLongitude, camPos.FocusPointAltitude, camPos.FocusPointAltitudeMode, camPos.Range, camPos.Tilt, camPos.Azimuth, 6);

% measure x and y - distances
u_li   = ge_app.GetPointOnTerrainFromScreenCoords(-1,-1);    % bottom left
u_re   = ge_app.GetPointOnTerrainFromScreenCoords(1,-1);     % bottom right
o_li   = ge_app.GetPointOnTerrainFromScreenCoords(-1, 1);    % upper_left
dist_u = pos2dist(u_li.Latitude ,u_li.Longitude, u_re.Latitude , u_re.Longitude  ,1);   % distance x  need:  http://www.mathworks.com/matlabcentral/fileexchange/5256-pos2dist
dist_o = pos2dist(u_li.Latitude ,u_li.Longitude, o_li.Latitude , o_li.Longitude  ,1);   % distance y

%  generate mesh for screencoords from -1 to 1 in each direction
[X,Y] = meshgrid(-1:(1/30):1, -1:(1/30):1);
Z     = NaN(size(X));    % preallocate for altitude-data
% get altitude from GE
for i= 1: 1 : size (X,1)
    for j =1:1: size (X,2)
        pot = ge_app.GetPointOnTerrainFromScreenCoords(X(i,j), Y(i,j) );
        Z(i,j) = pot.Altitude;
    end;
end;

% calculate distance an heigt
X_dist = (X + 1).*0.5.*dist_u;   % distance in km
Y_dist = (Y + 1).*0.5.*dist_o;   % distance in km
Z_dist = Z;                      % heigt in meter

% get a colored image from GE --> heavy workaround!!
jmouseemu([600,500],'normal');      % set actual view active            need: http://www.mathworks.com/matlabcentral/fileexchange/28357-jmouseemu-mouse-emulator-v2-2
inputemu('key_alt','c');            % copy image to clipboard           need: http://www.mathworks.com/matlabcentral/fileexchange/28603
pause(2);                           % wait for image in clipboard
[imdata] = clipboardimage('jpg');   % get image from clipboard;         need: http://www.mathworks.com/matlabcentral/fileexchange/14584-clipboardimage

% Convert to indexed image
[X2, map2] = rgb2ind(imdata, 256);  % not more than 256 because of uint8 is necessary

% Create surface plot an texture it with image
figure('Color', [1 1 1 ]);
surface(X_dist,Y_dist,Z_dist,flipud(X2),'FaceColor','texturemap','EdgeColor','none', 'CDataMapping','direct')
colormap(map2);
xlim([0 dist_u]);
ylim([0 dist_o]);
view([-12,72]);
box on
grid on


text:
function [ panoout ] = combineViews( Imgs, width, height )
%COMBINEVIEWS Combine separate views to panorama
%   Imgs: same format as separatePano

panoout = zeros(height, width, size(Imgs(1).img,3));
panowei = zeros(height, width, size(Imgs(1).img,3));
imgNum = length(Imgs);
for i = 1:imgNum
    [sphereImg validMap] = ...
        im2Sphere( Imgs(i).img, Imgs(i).fov, width, height, Imgs(i).vx, Imgs(i).vy);
    sphereImg(~validMap) = 0;
    panoout = panoout + sphereImg;
    panowei = panowei + validMap;
end
panoout(panowei==0) = 0;
panowei(panowei==0) = 1;
panoout = panoout./double(panowei);

end


text:
% DEMBRENDANFGPLVM3 Use the GP-LVM to model the Frey face data with DTCVAR.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'brendan';
experimentNo = 3;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('dtcvar');

latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
modelWriteResult(model, dataSetName, experimentNo);

if exist('printDiagram') & printDiagram
  lvmPrintPlot(model, lbls, dataSetName, experimentNo);
end

% Load the results and display dynamically.
lvmResultsDynamic(model.type, dataSetName, experimentNo, 'image', [20 28], 1, 0, 1)

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function f = NLP_objFunction(x, Prob)
% f = NLP_objFunction(x, Prob)
% From tomlab quickguide

if ~isfield(Prob, 'uP')
    alpha = 100;
else
    if isempty(Prob.uP)
        alpha = 100;
    else
        alpha = Prob.uP(1);
    end
end
f = alpha * (x(2) - x(1) ^ 2) ^ 2 + (1 - x(1)) ^ 2;


text:
function [x, c, funVal, ValueL]=sgLogisticR(A, y, z, opts)
%
%%
% Function sgLogisticR
%      Logistic Loss with the
%           sparse group Lasso Regularization
%
%% Problem
%
%  min  f(x,c) = - sum_i weight_i * log (p_i)
%                + z_1 \|x\|_1 + z_2 * sum_j w_j ||x_{G_j}||
%
%  a_i denotes a training sample,
%      and a_i' corresponds to the i-th row of the data matrix A
%
%  y_i (either 1 or -1) is the response
%
%  p_i= 1/ (1+ exp(-y_i (x' * a_i + c) ) ) denotes the probability
%
%  G_j's are nodes with tree structure
%
%  For this special case,
%    we have L1 for each element
%    and the L2 for the non-overlapping group
%
%  The tree overlapping group information is contained in
%  opts.ind, which is a 3 x nodes matrix, where nodes denotes the number of
%  nodes of the tree.
%  opts.ind(1,:) contains the starting index
%  opts.ind(2,:) contains the ending index
%  opts.ind(3,:) contains the corresponding weight (w_j)
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size mx1)
%  z -        Regularization parameter (z >=0)
%  opts-      Optional inputs (default value: opts=[])
%             !!We require that opts.ind is specified.!!
%
%% Output parameters:
%  x-         The obtained weight of size n x 1
%  c-         The obtained intercept (scalar)
%  funVal-    Function value during iterations
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified on February 19, 2010.
%
%% Related papers
%
% [1] Jun Liu and Jieping Ye, Moreau-Yosida Regularization for
%     Grouped Tree Structure Learning, NIPS 2010
%
%% Related functions:
%
%  sll_opts
%
%%

%% Verify and initialize the parameters
%%
if (nargin <4)
    error('\n Inputs: A, y, z, and opts (.ind) should be specified!\n');
end

[m,n]=size(A);

if (length(y) ~=m)
    error('\n Check the length of y!\n');
end

if (z(1)<0 || z(2)<0)
    error('\n z should be nonnegative!\n');
end

lambda1=z(1);
lambda2=z(2);

opts=sll_opts(opts); % run sll_opts to set default values (flags)

% restart the program for better efficiency
%  this is a newly added function
if (~isfield(opts,'rStartNum'))
    opts.rStartNum=opts.maxIter;
else
    if (opts.rStartNum<=0)
        opts.rStartNum=opts.maxIter;
    end
end

%% Detailed initialization
%% Normalization

% Please refer to the function 'sll_opts'
%                 for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly.
%     This implicit normalization is suggested for the sparse matrix,
%                                    but not for the dense matrix.
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Group & Others

% Initialize ind
if (~isfield(opts,'ind'))
    error('\n In sgLogisticR, the field .ind should be specified');
else
    ind=opts.ind;

    if (size(ind,1)~=3)
        error('\n Check opts.ind');
    end
end

% The parameter 'weight' contains the weight for each training sample.
% See the definition of the problem above.
% The summation of the weights for all the samples equals to 1.
if (isfield(opts,'sWeight'))
    sWeight=opts.sWeight;

    if ( length(sWeight)~=2 || sWeight(1) <=0 || sWeight(2) <= 0)
        error('\n Check opts.sWeight, which contains two positive values');
    end

    % we process the weight, so that the summation of the weights for all
    % the samples is 1.

    p_flag=(y==1);                  % the indices of the postive samples
    m1=sum(p_flag) * sWeight(1);    % the total weight for the positive samples
    m2=sum(~p_flag) * sWeight(2);   % the total weight for the positive samples

    weight(p_flag,1)=sWeight(1)/(m1+m2);
    weight(~p_flag,1)=sWeight(2)/(m1+m2);
else
    weight=ones(m,1)/m;             % if not specified, we apply equal weight
end

%% Starting point initialization

p_flag=(y==1);                  % the indices of the postive samples
m1=sum(weight(p_flag));         % the total weight for the positive samples
m2=1-m1;                        % the total weight for the positive samples

% process the regularization parameter
if (opts.rFlag==0)
    lambda=z;
else % z here is the scaling factor lying in [0,1]
    if (lambda1<0 || lambda1>1 || lambda2<0 || lambda2>1)
        error('\n opts.rFlag=1, and z should be in [0,1]');
    end

    % we compute ATb for computing lambda_max, when the input z is a ratio

    b(p_flag,1)=m2;  b(~p_flag,1)=-m1;
    b=b.*weight;

    % compute AT b
    if (opts.nFlag==0)
        ATb =A'*b;
    elseif (opts.nFlag==1)
        ATb= A'*b - sum(b) * mu';  ATb=ATb./nu;
    else
        invNu=b./nu;               ATb=A'*invNu-sum(invNu)*mu';
    end

    % compute lambda1_max
    temp=abs(ATb);
    lambda1_max=max(temp);

    lambda1=lambda1*lambda1_max;

    % compute lambda2_max(lambda_1)
    temp=max(temp-lambda1,0);

    if ( min(ind(3,:))<=0 )
        error('\n In this case, lambda2_max = inf!');
    end
    lambda2_max=computeLambda2Max(temp,n,ind,size(ind,2));

    lambda2=lambda2*lambda2_max;
end

% initialize a starting point
if opts.init==2
    x=zeros(n,1); c=log(m1/m2);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if (length(x)~=n)
            error('\n Check the input .x0');
        end
    else
        x=zeros(n,1);
    end

    if isfield(opts,'c0')
        c=opts.c0;
    else
        c=log(m1/m2);
    end
end

% compute A x
if (opts.nFlag==0)
    Ax=A* x;
elseif (opts.nFlag==1)
    invNu=x./nu; mu_invNu=mu * invNu;
    Ax=A*invNu -repmat(mu_invNu, m, 1);
else
    Ax=A*x-repmat(mu*x, m, 1);    Ax=Ax./nu;
end

%% The main program

if (opts.mFlag==0 && opts.lFlag==0)

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1/m; % the intial guess of the Lipschitz continuous gradient

    weighty=weight.*y;
    % the product between weight and y

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,1);
    cp=c;         ccp=0;

    %% The Armijo Goldstein line search schemes + accelearted gradient descent

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;  sc=c + beta* ccp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute As=A*s
        As=Ax + beta* (Ax-Axp);

        % aa= - diag(y) * (A * s + sc)
        aa=- y.*(As+ sc);

        % fun_s is the logistic loss at the search point
        bb=max(aa,0);
        fun_s= weight' * ( log( exp(-bb) +  exp(aa-bb) ) + bb );

        % compute prob=[p_1;p_2;...;p_m]
        prob=1./( 1+ exp(aa) );

        % b= - diag(y.* weight) * (1 - prob)
        b= -weighty.*(1-prob);

        gc=sum(b); % the gradient of c

        % compute g= AT b, the gradient of x
        if (opts.nFlag==0)
            g=A'*b;
        elseif (opts.nFlag==1)
            g=A'*b - sum(b) * mu';  g=g./nu;
        else
            invNu=b./nu;              g=A'*invNu-sum(invNu)*mu';
        end

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;
        cp=c;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the L1/Lq-norm regularized projection
            v=s-g/L;  c= sc- gc/L;

            % tree overlapping group Lasso projection
            ind_work(1:2,:)=[ [-1, -1]', ind(1:2,:) ];
            ind_work(3,:)=[ lambda1/L, ind(3,:) * (lambda2 / L) ];

            x=altra(v, n, ind_work, size(ind_work,2));

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute A x
            if (opts.nFlag==0)
                Ax=A* x;
            elseif (opts.nFlag==1)
                invNu=x./nu; mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
            end

            % aa= - diag(y) * (A * x + c)
            aa=- y.*(Ax+ c);

            % fun_x is the logistic loss at the new approximate solution
            bb=max(aa,0);
            fun_x= weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb );

            r_sum= (v'*v + (c-sc)^2) / 2;
            l_sum=fun_x - fun_s - v'* g - (c-sc)* gc;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is fun_x <= fun_s + v'* g + c * gc
            %                           + L/2 * (v'*v + (c-sc)^2 )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                %fprintf('\n L=%e, r_sum=%e',L, r_sum);
            end
        end

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        ValueL(iterStep)=L;
        % store values for L

        xxp=x-xp;    ccp=c-cp;
        funVal(iterStep)=fun_x;

        % compute the regularization part
        ind_work(1:2,:)=[ [-1, -1]', ind(1:2,:) ];
        ind_work(3,:)=[ lambda1, ind(3,:) * lambda2 ];
        tree_norm=treeNorm(x, n, ind_work, size(ind_work,2));

        % function value = loss + regularizatioin
        funVal(iterStep)=fun_x + tree_norm;

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp);    norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end

        % restart the program every opts.rStartNum
        if (~mod(iterStep, opts.rStartNum))
            alphap=0; alpha=1;
            xp=x; Axp=Ax; xxp=zeros(n,1); L =L/2;
        end
    end
else
    error('\n The function does not support opts.mFlag neq 0 & opts.lFlag neq 0!');
end

author: thomaskuestner