text:
function linplus_test19 ( )

%*****************************************************************************80
%
%% TEST19 tests R8CBB_FA, R8CBB_PRINT, R8CBB_RANDOM, R8CBB_SL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n1 = 8;
  n2 = 2;
  n = n1 + n2;
  ml = 1;
  mu = 1;
  na = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST19\n' );
  fprintf ( 1, '  For a compressed border banded matrix:\n' );
  fprintf ( 1, '  R8CBB_RANDOM randomly generates;\n' );
  fprintf ( 1, '  R8CBB_PRINT prints;\n' );
  fprintf ( 1, '  R8CBB_FA factors (no pivoting);\n' );
  fprintf ( 1, '  R8CBB_SL solves.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N     = %d\n', n );
  fprintf ( 1, '  Matrix suborder N1 = %d\n', n1 );
  fprintf ( 1, '  Matrix suborder N2 = %d\n', n2 );
  fprintf ( 1, '  Lower bandwidth ML = %d\n', ml );
  fprintf ( 1, '  Upper bandwidth MU = %d\n', mu );
%
%  Set the matrix.
%
  [ a, seed ] = r8cbb_random ( n1, n2, ml, mu, seed );

  r8cbb_print ( n1, n2, ml, mu, a, '  The R8CBB matrix:' );
%
%  Set the desired solution.
%
  x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
  b = r8cbb_mxv ( n1, n2, ml, mu, a, x );
%
%  Factor the matrix
%
  [ a_lu, info ] = r8cbb_fa ( n1, n2, ml, mu, a );

  r8cbb_print ( n1, n2, ml, mu, a_lu, '  The factored matrix:' );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST19 - Fatal error!\n' );
    fprintf ( 1, '  R8CBB_FA claims the matrix is singular.\n' );
    fprintf ( 1, '  The value of INFO is %d\n', info );
    return
  end
%
%  Solve the system.
%
  r8vec_print ( n, b, '  The right hand side vector:' );

  x = r8cbb_sl ( n1, n2, ml, mu, a_lu, b );

  r8vec_print ( n, x, '  Solution:' );

  return
end


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
% This function computes the MCG candidates given an image.
%  INPUT:
%  - image : Input image
%  - mode  : It can be: + 'fast'     (SCG in the paper)
%                       + 'accurate' (MCG in the paper)
%  - compute_masks : Compute the candidates maks [0] or not [1]. Note that
%                    it is very time consuming. Otherwise use the labels
%                    as shown in the demo.
%
%  OUTPUT:
%  - candidates : Struct containing the following fields
%          + superpixels : Label matrix of the superpixel partition
%          + labels : Cell containing the superpixel labels that form
%                     each of the candidates
%          + scores : Score of each of the ranked candidates
%          + masks  : 3D boolean matrix containing the masks of the candidates
%                     (only if compute_masks==1)
%  - ucm2       : Ultrametric Contour Map from which the candidates are
%                 extracted
%  - bboxes     : Bounding boxes of the candidates (up,left,down,right)
%                 See 'bboxes' folder for functions to work with them
%
%  DEMO:
%  - See demos/demo_im2mcg.m
% ------------------------------------------------------------------------
function [candidates, ucm2] = im2mcg(image,mode,compute_masks)
if nargin<2
    mode = 'fast';
end
if nargin<3
    compute_masks = 0;
end

% Load pre-trained Structured Forest model
sf_model = loadvar(fullfile(root_dir, 'datasets', 'models', 'sf_modelFinal.mat'),'model');

% Level of overlap to erase duplicates
J_th = 0.95;

% Max margin parameter
theta = 0.7;

if strcmp(mode,'fast')
    % Which scales to work on (MCG is [2, 1, 0.5], SCG is just [1])
    scales = 1;

    % Get the hierarchies at each scale and the global hierarchy
    [ucm2] = img2ucms(image, sf_model, scales);
    all_ucms = ucm2;

    % Load pre-trained pareto point
    pareto_n_cands = loadvar(fullfile(root_dir, 'datasets', 'models', 'scg_pareto_point_train2012.mat'),'n_cands');

    % Load pre-trained random forest regresssor for the ranking of candidates
    rf_regressor = loadvar(fullfile(root_dir, 'datasets', 'models', 'scg_rand_forest_train2012.mat'),'rf');

elseif strcmp(mode,'accurate')
    % Which scales to work on (MCG is [2, 1, 0.5], SCG is just [1])
    scales = [2, 1, 0.5];

    % Get the hierarchies at each scale and the global hierarchy
    [ucm2,ucms] = img2ucms(image, sf_model, scales);
    all_ucms = cat(3,ucm2,ucms(:,:,3),ucms(:,:,2),ucms(:,:,1)); % Multi, 0.5, 1, 2

    % Load pre-trained pareto point
    pareto_n_cands = loadvar(fullfile(root_dir, 'datasets', 'models', 'mcg_pareto_point_train2012.mat'),'n_cands');

    % Load pre-trained random forest regresssor for the ranking of candidates
    rf_regressor = loadvar(fullfile(root_dir, 'datasets', 'models', 'mcg_rand_forest_train2012.mat'),'rf');
else
    error('Unknown mode for MCG: Possibilities are ''fast'' or ''accurate''')
end
% ------------------------------------

% Transform ucms to hierarchies (dendogram) and put them all together
n_hiers = size(all_ucms,3);
lps = [];
ms  = cell(n_hiers,1);
ths = cell(n_hiers,1);
for ii=1:n_hiers
    % Transform the UCM to a hierarchy
    curr_hier = ucm2hier(all_ucms(:,:,ii));
    ths{ii}.start_ths = curr_hier.start_ths';
    ths{ii}.end_ths   = curr_hier.end_ths';
    ms{ii}            = curr_hier.ms_matrix;
    lps = cat(3, lps, curr_hier.leaves_part);
end

% Get full cands, represented on a fused hierarchy
[f_lp,f_ms,cands,start_ths,end_ths] = full_cands_from_hiers(lps,ms,ths,pareto_n_cands);

% Hole filling and complementary candidates
if ~isempty(f_ms)
    [cands_hf, cands_comp] = hole_filling(double(f_lp), double(f_ms), cands); %#ok<NASGU>
else
    cands_hf = cands;
    cands_comp = cands; %#ok<NASGU>
end

% Select which candidates to keep (Uncomment just one line)
cands = cands_hf;                       % Just the candidates with holes filled
% cands = [cands_hf; cands_comp];         % Holes filled and the complementary
% cands = [cands; cands_hf; cands_comp];  % All of them

% Compute base features
b_feats = compute_base_features(f_lp, f_ms, all_ucms);
b_feats.start_ths = start_ths;
b_feats.end_ths   = end_ths;
b_feats.im_size   = size(f_lp);

% Filter by overlap
red_cands = mex_fast_reduction(cands-1,b_feats.areas,b_feats.intersections,J_th);

% Compute full features on reduced cands
[feats, bboxes] = compute_full_features(red_cands,b_feats);

% Rank candidates
class_scores = regRF_predict(feats,rf_regressor);
[scores, ids] = sort(class_scores,'descend');
red_cands = red_cands(ids,:);
bboxes = bboxes(ids,:);
if isrow(scores)
    scores = scores';
end

% Max margin
[new_ids, mm_scores] = mex_max_margin(red_cands-1,scores,b_feats.intersections,theta); %#ok<NASGU>
cand_labels = red_cands(new_ids,:);
candidates.scores = scores(new_ids);
bboxes = bboxes(new_ids,:);

% Change the coordinates of bboxes to be coherent with
% other results from other sources (sel_search, etc.)
candidates.bboxes = [bboxes(:,2) bboxes(:,1) bboxes(:,4) bboxes(:,3)];

% Get the labels of leave regions that form each candidates
candidates.superpixels = f_lp;
if ~isempty(f_ms)
    candidates.labels = cands2labels(cand_labels,f_ms);
else
    candidates.labels = {1};
end

% Transform the results to masks
if compute_masks
    if ~isempty(f_ms)
        candidates.masks = cands2masks(cand_labels, f_lp, f_ms);
    else
        candidates.masks = true(size(f_lp));
    end
end

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
classdef MaOPP_binary < PROBLEM
% <many> <binary> <large/none> <expensive/none>
% Many-objective pathfinding problem based on binary encoding
% xmax                    ---  10 --- xmax
% ymax                    ---  10 --- ymax
% obstacleValue           ---   0 --- obstacleValue
% nh                      ---   1 --- nh
% neighbourhood           ---   2 --- neighbourhood
% backtracking            ---   0 --- backtracking
% allowObstaclesOnPath    ---   1 --- allowObstaclesOnPath
% overheadVariablesFactor --- 1.5 --- overheadVariablesFactor

%------------------------------- Reference --------------------------------
% J. Weise and S. Mostaghim, A scalable many-objective pathfinding
% benchmark suite, IEEE Transactions on Evolutionary Computation, 2022,
% 26(1): 188-194.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% The most recent reference fronts and sets can be found at
% https://ci.ovgu.de/Publications/Scalable+Many_Objective+Pathfinding+Benchmark+Suite-p-910.html
% In case of any question, please contact me at
% jens.weise@ovgu.de

    properties(Access = protected)
        x_max;          % size [1,x_max] in x direction
        y_max;          % size [1,y_max] in y direction
        delay;          % expected delay
        nh;             % elevation function {1,2,3,M} -> 1 = 1, 2 = 2, 3 = 3, 4 = M
        neighborhood;   % 2^k, k {2,3}
        backtracking;   % True, False
        obstacle;       %{0,1,2} -> 0 = No, 1 = CH, 2 = LA
        vmax_high;
        vmax_medium;
        vmax_low;
        v_max;
        elevation;
        allowObstaclesOnPath;
        upperBoundsForObjectives;
        overheadVariablesFactor;
        numerOfBits;
    end

    methods
        %% Default settings of the problem
        function Setting(obj)
            [obj.x_max, obj.y_max, obstacleValue, obj.nh, obj.neighborhood, obj.backtracking, obj.allowObstaclesOnPath, obj.overheadVariablesFactor] = obj.ParameterSet(10,10,0,1,2,0,1,1);
            if isempty(obj.vmax_high); obj.vmax_high = 130; end
            if isempty(obj.vmax_medium); obj.vmax_medium = 100; end
            if isempty(obj.vmax_low); obj.vmax_low = 50; end
            obj.M = 5;
            obj.setBinaryEncoding();
            obj.v_max     = zeros(obj.x_max, obj.y_max);  % matrix of shape (x_max,y_max), matrix with different v values for each cell
            obj.elevation = zeros(obj.x_max, obj.y_max);  % matrix of shape (x_max,y_max)
            for x = 1 : obj.x_max
                for y = 1 : obj.y_max
                    obj.obstacle(x,y) = obstacleValue;
                    w_x_y = max(sin(x-1), cos(y-1));
                    if w_x_y > 0.9
                        obj.v_max(x, y) = obj.vmax_high;
                    elseif w_x_y < -0.4
                        obj.v_max(x, y) = obj.vmax_low;
                    else
                        obj.v_max(x, y) = obj.vmax_medium;
                    end
                    % velocity matrix obstacles for CH
                    if obj.obstacle(x,y) == 1 && (sign(sin(pi/2 + pi*x)) + sign(sin(pi/2 + pi*y)) - 2* (heaviside(x-obj.x_max+0.5)- heaviside(x-obj.x_max-0.5))* (heaviside(y-obj.y_max+0.5)- heaviside(y-obj.y_max-0.5)) == 2)
                        obj.v_max(x, y) = 0;
                    end
                    % velocity matrix obstacles for LA
                    if obj.obstacle(x,y) == 2 && (x-1-obj.x_max/2)^2 + (y-1-obj.y_max/2)^2 - (0.25*obj.x_max)^2 <0 %Q: using radius of x_max/4 from paper page 4 line 22 A: Yes.
                        obj.v_max(x, y) = 0;
                    end
                    xs = map(x,1,obj.x_max+1,-3,3);
                    ys = map(y,1,obj.y_max+1,-3,3);
                    if obj.nh == 1
                        obj.elevation(x,y) = 5*exp(-(xs-1.5)^2-(ys+1.5)^2);
                    elseif obj.nh == 2
                        obj.elevation(x,y) = 5*exp(-(xs+1.5)^2-(ys+1.5)^2) + 5*exp(-(xs-1.5)^2-(ys-1.5)^2);
                    elseif obj.nh == 3
                        obj.elevation(x,y) = 5*exp(-(xs+1.5)^2-(ys+1.5)^2) + 5*exp(-(xs-1.5)^2-(ys-1.5)^2) + 5*exp(-(xs-1.5)^2-(ys+1.5)^2);
                    elseif obj.nh == 4
                        obj.elevation(x,y) = 3*(1-xs)^2*exp(-(xs^2)-(ys+1)^2)-10*exp(-xs^2-ys^2)*(-xs^3+xs/5-ys^5)-1/3*exp(-(xs+1)^2-ys^2);
                    end
                end
            end
            manhattanSteps = (obj.x_max + obj.y_max - 2);
            obj.upperBoundsForObjectives = [manhattanSteps*1.5, manhattanSteps*1.5, 1.5*5*obj.nh, 1.5*manhattanSteps/50, 0.5*manhattanSteps*pi/2];
        end
        function setBinaryEncoding(obj)
            if obj.neighborhood == 2 && obj.backtracking == 0
                obj.D = 1*(obj.x_max + obj.y_max - 2)*obj.overheadVariablesFactor;
                obj.numerOfBits = 1;
            elseif obj.neighborhood == 2 && obj.backtracking == 1
                obj.D = 2*(obj.x_max + obj.y_max - 2)*obj.overheadVariablesFactor;
                obj.numerOfBits = 2;
            elseif obj.neighborhood == 3 && obj.backtracking == 0
                obj.D = 2*(obj.x_max + obj.y_max - 2)*obj.overheadVariablesFactor;
                obj.numerOfBits = 2;
            elseif obj.neighborhood == 3 && obj.backtracking == 1
                obj.D = 3*(obj.x_max + obj.y_max - 2)*obj.overheadVariablesFactor;
                obj.numerOfBits = 3;
            end
            obj.encoding = 4 + zeros(1,obj.D);
        end
        function [x_coords, y_coords,D] = decodePath(obj,PopDec,x_max,y_max)
            D = size(PopDec,2);
            n = obj.numerOfBits;
            for i = 0 : (D/n)-1
               temp = PopDec(:,[n*i+1:n*i+n]);
               realPopDec(:,i+1) = b2d(temp);
            end
            realPopDec = 1.0*realPopDec/(2^n);
            [x_coords, y_coords, D] = decodePathWithRealValues(obj,realPopDec,x_max,y_max);
            D = size(realPopDec,2);
        end
        function [x_coords, y_coords, D] = decodePathWithRealValues(obj,PopDec,x_max,y_max)
            [N,D] = size(PopDec);
            % convert the array into a sequence of x,y coordinates
            x_coords = zeros(N,D+1);
            y_coords = zeros(N,D+1);
            x_coords(:,1) = 1;
            y_coords(:,1) = 1;

            for i = 1 : N
                x_current = 1;
                y_current = 1;
                for d = 1 : D
                    possibleNeighours      = obj.getPossNeighbours(x_current,y_current,obj.backtracking,obj.neighborhood,obj.allowObstaclesOnPath);
                    noOfPossibleNeighbours = size(possibleNeighours,1);
                    neighbourToTake        = noOfPossibleNeighbours;
                    for a = 1 : noOfPossibleNeighbours
                        if PopDec(i,d) >= (a-1)/noOfPossibleNeighbours && PopDec(i,d) < (a)/noOfPossibleNeighbours
                            neighbourToTake = a;
                            break;
                        end
                    end
                    if neighbourToTake == 0 % Path ends if there are no neighbours available
                        break;
                    end
                    addvals   = possibleNeighours(neighbourToTake,:);
                    xAdd      = addvals(1);
                    yAdd      = addvals(2);
                    x_current = x_current + xAdd;
                    y_current = y_current + yAdd;
                    x_coords(i,d+1) = x_current;
                    y_coords(i,d+1) = y_current;
                    if x_current == obj.x_max && y_current == obj.y_max
                        break;
                    end
                end
            end
        end
        function neighbours = getPossNeighbours(obj,xCurr,yCurr,backtracking,neighborhood,allowObstaclesOnPath)
            e  = [1 0];
            se = [1 1];
            s  = [0 1];
            sw = [-1 1];
            w  = [-1 0];
            nw = [-1 -1];
            n  = [0 -1];
            ne = [1 -1];
            currentCoordinates = [xCurr yCurr];
            if backtracking == 0
                if neighborhood == 2
                    neighbours = [e;s];
                else
                    neighbours=[e;se;s];
                end
            else
                if neighborhood == 2
                    neighbours = [e;s;w;n];
                else
                    neighbours = [e;se;s;sw;w;nw;n;ne];
                end
            end
            % should we remove the outliers?
            % it is closer to the original implementation
            nextCoordinates = neighbours+currentCoordinates;
            out = nextCoordinates(:,1)<1 | nextCoordinates(:,1)>obj.x_max | nextCoordinates(:,2)<1 | nextCoordinates(:,2)>obj.y_max;
            neighbours(out,:) = [];
            if allowObstaclesOnPath == 0
                nextCoordinates = neighbours+currentCoordinates;
                obstancles = getVmaxValues(nextCoordinates(:,2),nextCoordinates(:,1),obj.v_max) == 0;
                neighbours(obstancles,:) = [];
            end
        end
        %% separate the calculations here to be able to compute the feasibility rate as well (with a lot of extra runtime, but whatever..)
        function [objectives,distToTarget,fObstacle] = objectiveValuesPath(obj,PopDec)
            % read out parameters
            N = size(PopDec,1);
            [x_coords, y_coords, D] = obj.decodePath(PopDec,obj.x_max,obj.y_max);
            % for entire population (one column) iterate over nodes (columns)
            f1 = zeros(N,1); % euclidean distance
            f2 = zeros(N,1); % expected delays
            f3 = zeros(N,1); % elevations
            f4 = zeros(N,1); % traveling time
            f5 = zeros(N,1); % smoothness
            fObstacle = zeros(N,1);
            distances = zeros(N,D); % stores distances that were calculated for f1 for later use in f5
            c = [];
            for d = 1 : D
                x1        = x_coords(:,d);
                x2        = x_coords(:,d+1);
                y1        = y_coords(:,d);
                y2        = y_coords(:,d+1);
                hasEnded  = (x2~=0&x1~=0);
                x2(x2==0) = 1;
                y2(y2==0) = 1;
                x1(x1==0) = 1;
                y1(y1==0) = 1;

                vMaxValues1 = getVmaxValues(x1,y1,obj.v_max);
                vMaxValues2 = getVmaxValues(x2,y2,obj.v_max);

                if obj.allowObstaclesOnPath==1
                    obstancles = vMaxValues1==0|vMaxValues2==0;
                    fObstacle  = fObstacle+obstancles;
                else
                    obstancles = zeros(N,1);
                    fObstacle  = fObstacle+obstancles;
                end

                %---Objective 1: euclidean distance
                distances(:,d) = sqrt((x2-x1).^2+(y2-y1).^2);
                res = (distances(:,d)).*hasEnded;
                res(res>sqrt(2)) = 0;
                f1 = f1 + res;

                %---Objective 2: delays
                inds      = vMaxValues1 ~= vMaxValues2;
                del       = inds*2;
                inds      = vMaxValues1 == vMaxValues2;
                loc       = find(inds==1);
                inds      = zeros(N,1);
                inds(loc) = vMaxValues1(loc);
                inds50    = inds == 50;
                del       = del + inds50*3;
                inds100   = inds == 100;
                del       = del + inds100*1;
                inds0     = del == 0;
                del       = del + inds0*0.2;
                res       = del.*hasEnded;
                f2        = f2 + res;

                %---Objective 3: elevations
                h1   = getVmaxValues(x1,y1,obj.elevation);
                h2   = getVmaxValues(x2,y2,obj.elevation);
                inds = h1 < h2;
                res  = (inds.*(h2-h1)).*hasEnded;
                f3   = f3 + res;

                %---Objective 4: traveling time
                res = (2*distances(:,d)./(vMaxValues1+vMaxValues2));
                res(isinf(res)|isnan(res)) = 0;
                if any(isnan(res))
                    errr;
                end
                res = res.*hasEnded;
                f4  = f4 + res;

                %---Objective 5: smoothness
                if d>1
                    x0       = x_coords(:,d-1);
                    y0       = y_coords(:,d-1);
                    x0(x0<1) = 1;
                    y0(y0<1) = 1;

                    test = transpose(dot([(x0-x1).'; (y0-y1).'],[(x1-x2).'; (y1-y2).'],1))./(distances(:,d-1).*distances(:,d));
                    test(isnan(test)) = 0;
                    try
                    	assert(all(test <= 1) && all(test >= -1))
                    catch
                    end
                    res = (acos(transpose(dot([(x0-x1).'; (y0-y1).'],[(x1-x2).'; (y1-y2).'],1))./(distances(:,d-1).*distances(:,d)))).*hasEnded;
                    res(isnan(res) )= 0;
                    f5 = f5 + res;
                end
            end
            objectives = [f1,f2,f3,f4,f5];
            try
                last = [];
                for i = 1 : N
                    last = [last find(x_coords(i,:)~=0,1,'last')];
                end
                last = last';
                lxc  = [];
                for k = 1 : N
                    lxc(k,1) = x_coords(k,last(k));
                    lxc(k,2) = y_coords(k,last(k));
                end
            catch
            end
            target       = [repmat(obj.x_max,N,1) repmat(obj.y_max,N,1)];
            distToTarget = abs(target(:,1)-lxc(:,1))+abs(target(:,2)-lxc(:,2));
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            [objectives,distToTarget,fObstacle] = objectiveValuesPath(obj,varargin{1});
            objectives(:,1)                     = objectives(:,1) + 2*distToTarget + fObstacle;
            Population = SOLUTION(varargin{1},objectives,distToTarget+fObstacle,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        % Here, the pre-computed fronts are read. If it is not available,
        % an educated-guess Nadir point is computed. Please check the
        % website above in the header to see, if there are new fronts
        % available.
        function R = GetOptimum(obj,N)
            obs = '';
            if obj.obstacle == 0
                obs = 'NO';
            elseif obj.obstacle == 1
                obs = 'CH';
            elseif obj.obstacle == 2
                obs = 'LA';
            end
            sizeX = obj.x_max;
            sizeY = obj.y_max;
            ele   = '';
            if obj.nh == 1
                ele = '1';
            elseif obj.nh == 2
                ele = '2';
            elseif obj.nh == 3
                ele = '3';
            elseif obj.nh == 4
                ele = 'M';
            end
            k = obj.neighborhood;
            if obj.backtracking
                bt = 'T';
            else
                bt = 'F';
            end
            problem = sprintf('ASLETISMAC_%s_X%i_Y%i_P%s_K%i_B%s',obs,sizeX,sizeY,ele,k,bt);
            try
                CallStack = dbstack('-completenames');
                load(fullfile(fileparts(CallStack(1).file),'PF.mat'),problem);
                order = [2 5 1 3 4];
                eval(['R=',problem,';']);
                R      = R(:,order);
                R(:,5) = deg2rad(R(:,5));
            catch
                R = repmat(obj.upperBoundsForObjectives,N,1);
            end
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

function out= map(x,in_min,in_max,out_min,out_max)
    out=(x - in_min) * (out_max - out_min) / (in_max - in_min) + out_min;
end

function vList = getVmaxValues(x1,y1,aMatrix)
    len = size(x1,1);
    vList = zeros(len,1);
    for i = 1:len
        vList(i) = aMatrix(x1(i),y1(i));
    end
end

% Converts binary to decimal (MSB 0)!
function dec = b2d(vector)
    l   = size(vector);
    dec = ones(l(1),1);
    for r = 1 : l(1)
        dTemp = 0;
        for i = 1 : l(2)
            dTemp = dTemp + vector(r,i)*2^(i-1);
        end
        dec(r,1) = dTemp;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function f = projectOntoBMCII(f)
%PROJECTONTOBMCII   Projection onto BMC-II symmetry.
%   F = projectOntoBMCII(F) is the orthogonal projection of f onto BMC-II
%   symmetry, i.e., a function that is
%   1. even in r for every even wave number in theta.
%   2. odd in r for every odd wave number in theta.
%   Additionally, for all but the k=0 mode in theta, the resulting
%   projection enforces the diskfun is zero at the origin.
%
%   The projection is orthogonal, i.e., the correction matrix to fix up the
%   structure has the smallest possible Frobenius norm.

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

% Even part
feven = f;
feven.cols = feven.cols(:, feven.idxPlus);
feven.rows = feven.rows(:, feven.idxPlus);
feven = projectOntoEvenBMCII(feven);

% Odd part
fodd = f;
fodd.cols = fodd.cols(:, fodd.idxMinus);
fodd.rows = fodd.rows(:, fodd.idxMinus);
fodd = projectOntoOddBMCII(fodd);

% Put pieces back together.
f.cols(:, f.idxPlus) = feven.cols;
f.rows(:, f.idxPlus) = feven.rows;

f.cols(:, f.idxMinus) = fodd.cols;
f.rows(:, f.idxMinus) = fodd.rows;

end

function f = projectOntoEvenBMCII(f)
% Project a diskfun to have even BMC-II symmetry, i.e., a diskfun that
% is pi-periodic in theta and even in r. The projection is orthogonal,
% i.e., the correction matrix to fix up the structure has the smallest
% possible Frobenius norm.

% Nothing to project
if isempty(f)
    return;
end

% Operate on the column coefficients first.
X = f.cols.funs{1}.onefun.coeffs;

% Get size:
[m, n] = size(X);

 % First we enforce that the every function in r is even:
 X(2:2:end, :) = 0;

% The rest of the code now needs to operate on the remaining even,
% non-zero modes in theta.
evenModes = 2:n;

if ( ~isempty(evenModes) )
    % to enforce that f(r) = 0 for each column function,
    % we want C with smallest 2-norm such that  A*(X(1:2:end, 2:end)+C) = 0,
    % where A =[1 -1 1 -1..] .
    % The solution is
    % C = A'*((A*A')\(A*X)).
    %odd modes in r are zero, they won't contribute.
    even = 1:2:m;
    Xe = X(even, evenModes);
    factor = 1/length(even)*(sum(Xe(1:2:end, :),1)-sum(Xe(2:2:end, :),1));
    C = ((-1*ones(length(even), 1)).^((2:length(even)+1)'))*factor;
    %now add C to X
    X(even, evenModes) = Xe+C;

end

ctechs = chebtech2({'',X});
f.cols.funs{1}.onefun = ctechs;

% Now operate on the rows. The coefficients for the rows of an even BMC-II
% function should only contain even wave numbers. The projection is to
% simply zero out the odd wave numbers.
X = f.rows.funs{1}.onefun.coeffs;
n = size(X, 1);
zeroMode = floor(n/2) + 1;
oddModes = [fliplr(zeroMode-1:-2:1) zeroMode+1:2:n];
X(oddModes, :) = 0;
rtechs = real(trigtech({'', X}));
f.rows.funs{1}.onefun = rtechs;

% Weird feval behavior in chebfun requires this
f.cols.pointValues = feval(ctechs, [-1; 1]);
f.rows.pointValues = feval(rtechs, [-1; 1]);

end

function f = projectOntoOddBMCII(f)
% Project a diskfun to have odd BMC-II symmetry, i.e., a diskfun that is
% pi-anti-periodic in theta and even in r. The projection is orthogonal,
% i.e., the correction matrix to fix up the structure has the smallest
% possible Frobenius norm.

% Nothing to project
if isempty(f)
    return;
end

% to enforce odd, simply zero out even coeffs:

% Operate on the column coefficients first to project them onto odd
% functions.
X = f.cols.funs{1}.onefun.coeffs;
X(1:2:end, :) = 0;

ctechs = chebtech2({'',X});
f.cols.funs{1}.onefun = ctechs;

% Now operate on the rows. The coefficients for the rows of an odd BMC-II
% function should only contain odd wave numbers. The projection is to
% simply zero out the even wave numbers.
X = f.rows.funs{1}.onefun.coeffs;
n = size(X, 1);
zeroMode = floor(n/2) + 1;
evenModes = [fliplr(zeroMode-2:-2:1) zeroMode:2:n];
X(evenModes, :) = 0;

rtechs = real(trigtech({'', X}));
f.rows.funs{1}.onefun = rtechs;

% Weird feval behavior in chebfun requires this
f.cols.pointValues = feval(ctechs, [-1; 1]);
f.rows.pointValues = feval(rtechs, [-1; 1]);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function title = p05_title ( option )

%*****************************************************************************80
%
%% P05_TITLE sets the title for problem 5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, string TITLE, the title of the problem.
%
  if ( option == 1 )
    title = 'Wacker function, A = 0.1.';
  elseif ( option == 2 )
    title = 'Wacker function, A = 0.5.';
  elseif ( option == 3 )
    title = 'Wacker function, A = 1.0.';
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [all_theta] = oneVsAll(X, y, num_labels, lambda)
%ONEVSALL trains multiple logistic regression classifiers and returns all
%the classifiers in a matrix all_theta, where the i-th row of all_theta
%corresponds to the classifier for label i
%   [all_theta] = ONEVSALL(X, y, num_labels, lambda) trains num_labels
%   logisitc regression classifiers and returns each of these classifiers
%   in a matrix all_theta, where the i-th row of all_theta corresponds
%   to the classifier for label i

% Some useful variables
m = size(X, 1);
n = size(X, 2);

% You need to return the following variables correctly
all_theta = zeros(n + 1 , num_labels);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: You should complete the following code to train num_labels
%               logistic regression classifiers with regularization
%               parameter lambda.
%
% Hint: theta(:) will return a column vector.
%
% Hint: You can use y == c to obtain a vector of 1's and 0's that tell use
%       whether the ground truth is true/false for this class.
%
% Note: For this assignment, we recommend using fmincg to optimize the cost
%       function. It is okay to use a for-loop (for c = 1:num_labels) to
%       loop over the different classes.
%
%       fmincg works similarly to fminunc, but is more efficient when we
%       are dealing with large number of parameters.
%
% Example Code for fmincg:
%
%     % Set Initial theta
%     initial_theta = zeros(n + 1, 1);
%
%     % Set options for fminunc
%     options = optimset('GradObj', 'on', 'MaxIter', 50);
%
%     % Run fmincg to obtain the optimal theta
%     % This function will return theta and the cost
%     [theta] = ...
%         fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), ...
%                 initial_theta, options);
%

classes = zeros(m, num_labels);
for i=1:num_labels
    classes(:,i) = y==i;
end

initial_theta = zeros(n + 1, 1);

% Set options for fminunc
options = optimset('GradObj', 'on', 'MaxIter', 50);

% Run fmincg to obtain the optimal theta
% This function will return theta and the cost
for i = 1:num_labels
    [all_theta(:,i)] = fmincg (@(t)(lrCostFunction(t, X, classes(:,i), lambda)), ...
                  initial_theta, options);
end

all_theta = all_theta';

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

end

repo: MachineLearningEx


text:
%vgg_X_from_xP_nonlin  Estimation of 3D point from image matches and camera matrices, nonlinear.
%   X = vgg_X_from_xP_lin(x,P,imsize) computes max. likelihood estimate of projective
%   3D point X (column 4-vector) from its projections in K images x (2-by-K matrix)
%   and camera matrices P (K-cell of 3-by-4 matrices). Image sizes imsize (2-by-K matrix)
%   are needed for preconditioning.
%   By minimizing reprojection error, Newton iterations.
%
%   X = vgg_X_from_xP_lin(x,P,imsize,X0) takes initial estimate of X. If X0 is omitted,
%   it is computed by linear algorithm.
%
%   See also vgg_X_from_xP_lin.

% werner@robots.ox.ac.uk, 2003

function X = vgg_X_from_xP_nonlin(u,P,imsize,X)

if iscell(P)
  P = cat(3,P{:});
end
K = size(P,3);
if K < 2
  error('Cannot reconstruct 3D from 1 image');
end

if nargin==3
  X = vgg_X_from_xP_lin(u,P,imsize);
end

if nargin==2
  X = vgg_X_from_xP_lin(u,P);
end

% precondition
if nargin>2
  for k = 1:K
    H = [2/imsize(1,k) 0 -1
         0 2/imsize(2,k) -1
         0 0              1];
    P(:,:,k) = H*P(:,:,k);
    u(:,k) = H(1:2,1:2)*u(:,k) + H(1:2,3);
  end
end

% Parametrize X such that X = T*[Y;1]; thus x = P*T*[Y;1] = Q*[Y;1]
[dummy,dummy,T] = svd(X',0);
T = T(:,[2:end 1]);
for k = 1:K
  Q(:,:,k) = P(:,:,k)*T;
end

% Newton
Y = [0;0;0];
eprev = inf;
for n = 1:10
  [e,J] = resid(Y,u,Q);
  if 1-norm(e)/norm(eprev) < 1000*eps
    break
  end
  eprev = e;
  Y = Y - (J'*J)\(J'*e);
end

X = T*[Y;1];

return

%%%%%%%%%%%%%%%%%%%%%%%%%%

function [e,J] = resid(Y,u,Q)
K = size(Q,3);
e = [];
J = [];
for k = 1:K
  q = Q(:,1:3,k);
  x0 = Q(:,4,k);
  x = q*Y + x0;
  e = [e; x(1:2)/x(3)-u(:,k)];
  J = [J; [x(3)*q(1,:)-x(1)*q(3,:)
           x(3)*q(2,:)-x(2)*q(3,:)]/x(3)^2];
end
return

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function f = abs( f )
%ABS Absolute value of a BALLFUN.
%   ABS(F) is the absolute value of the BALLFUN F. This function does not work
%   if the function passes through or becomes numerically close to zero.

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

if ( isempty(f) )
    return
end

if ( f.isReal )
    % Positive/negative test.
    % Returns TRUE if there is no sign change.

    tol = chebfun2eps;
    bool = false;                  % Assume false

    % Evaluate on a grid:
    X = sample(f);

    X = X(:);

    if ( all( X >= -tol * f.vscale ))   % If all values are nonnegative
        bool = true;
    elseif ( all( X <= tol * f.vscale))  % If all values are not positive
        bool = true;
    end
    if ( ~bool )
        error('CHEBFUN:BALLFUN:abs:notSmooth', ...
            'Sign change detected. Unable to represent the result.');
    end
end

% Still call the constructor in case we missed a change of sign.
f = compose( f, @abs );
end


text:
function [Y,index] = Shuffle(X)
% [Y,index] = Shuffle(X)
%
% Randomly sorts X.
% If X is a vector, sorts all of X, so Y = X(index).
% If X is an m-by-n matrix, sorts each column of X, so
%	for j=1:n, Y(:,j)=X(index(:,j),j).
%
% Also see SORT, Sample, Randi, and RandSample.
%
% xx/xx/92  dhb  Wrote it.
% 10/25/93  dhb  Return index.
% 5/25/96   dgp  Made consistent with sort and "for i=Shuffle(1:10)"
% 6/29/96	dgp	 Edited comments above.
%
% Probably from PTB.  I think this is Pelli and Brainard.

% I guess rand() never returns two identical numbers.
[null,index] = sort(rand(size(X)));
Y = X(index);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function counter=population(i,j,k,order)
%-fanDTasia ToolBox------------------------------------------------------------------
% This Matlab script is part of the fanDTasia ToolBox: a Matlab library for Diffusion
% Weighted MRI (DW-MRI) Processing, Diffusion Tensor (DTI) Estimation, High-order
% Diffusion Tensor Analysis, Tensor ODF estimation, Visualization and more.
%
% A Matlab Tutorial on DW-MRI can be found in:
% http://www.cise.ufl.edu/~abarmpou/lab/fanDTasia/tutorial.php
%
%-CITATION---------------------------------------------------------------------------
% If you use this software please cite the following work:
% A. Barmpoutis and B.C. Vemuri, "A Unified Framework for Estimating Diffusion Tensors
% of any order with Symmetric Positive-Definite Constraints",
% In the Proceedings of ISBI, 2010
%
%-DESCRIPTION------------------------------------------------------------------------
% This function computes how many identical copies of a tensor coefficient appear in
% a fully symmetric high order tensor in 3 variables of spesific order.
%
%-USE--------------------------------------------------------------------------------
% counter=population(i,j,k,order)
%
% order: is the order of the symmetric tensor in 3 variables (x,y,z)
% i,j,k: is a triplet that represent the coefficient which is weighted by the monomial
% x^i*y^j*z^k, where i+j+k=order
%
%-DISCLAIMER-------------------------------------------------------------------------
% You can use this source code for non commercial research and educational purposes
% only without licensing fees and is provided without guarantee or warrantee expressed
% or implied. You cannot repost this file without prior written permission from the
% authors. If you use this software please cite the following work:
% A. Barmpoutis and B.C. Vemuri, "A Unified Framework for Estimating Diffusion Tensors
% of any order with Symmetric Positive-Definite Constraints", In Proc. of ISBI, 2010.
%
%-AUTHOR-----------------------------------------------------------------------------
% Angelos Barmpoutis, PhD
% Computer and Information Science and Engineering Department
% University of Florida, Gainesville, FL 32611, USA
% abarmpou at cise dot ufl dot edu
%------------------------------------------------------------------------------------

    size=3^order;
    counter=0;
    for z=0:size-1
        c=populationBasis(z,order,[0 0 0]);
        if (c(1)==i)&(c(2)==j)&(c(3)==k)
            counter=counter+1;
        end
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ret=populationBasis(i,order,c)
	if order==0
        ret=c;
    else
	    j=mod(i,3);
	    c(j+1)=c(j+1)+1;
	    ret=populationBasis((i-j)/3,order-1,c);
    end

author: Sable

repo: mcbench-benchmarks


text:
%% Advanced Morphology Transformations Demo
%
% In this sample you will learn how to use the OpenCV function
% <matlab:doc('cv.morphologyEx') cv.morphologyEx> to apply
% Morphological Transformation such as:
%
% * Opening
% * Closing
% * Morphological Gradient
% * Top Hat
% * Black Hat
%
% Sources:
%
% * <https://docs.opencv.org/3.2.0/d3/dbe/tutorial_opening_closing_hats.html>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/ImgProc/Morphology_2.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/morphology2.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/python/morphology.py>
%

%% Theory
%
% We previously covered two basic Morphology operations:
%
% * Erosion
% * Dilation
%
% Based on these two we can effectuate more sophisticated transformations to
% our images. Here we discuss briefly 5 operations offered by OpenCV.
%
%% 1) Opening
%
% It is obtained by the erosion of an image followed by a dilation.
%
% $$dst = open( src, element) = dilate( erode( src, element ) )$$
%
% Useful for removing small objects (it is assumed that the objects are bright
% on a dark foreground). For instance, check out the example below. The image
% at the left is the original and the image at the right is the result after
% applying the opening transformation. We can observe that the small spaces in
% the corners of the letter tend to disappear.
%
% <<https://docs.opencv.org/3.2.0/Morphology_2_Tutorial_Theory_Opening.png>>
%
% For the sake of clarity, we have performed the opening operation (|7x7|
% rectangular structuring element) on the same original image but inverted
% such as the object in white is now the letter.
%
% <<https://docs.opencv.org/3.2.0/Morphology_2_Tutorial_Theory_Opening_2.png>>
%
% (Left image: original image inverted, right image: resulting opening).
%
%% 2) Closing
%
% It is obtained by the dilation of an image followed by an erosion.
%
% $$dst = close( src, element ) = erode( dilate( src, element ) )$$
%
% Useful to remove small holes (dark regions).
%
% <<https://docs.opencv.org/3.2.0/Morphology_2_Tutorial_Theory_Closing.png>>
%
% On the inverted image, we have performed the closing operation (|7x7|
% rectangular structuring element):
%
% <<https://docs.opencv.org/3.2.0/Morphology_2_Tutorial_Theory_Closing_2.png>>
%
% (Left image: original image inverted, right image: resulting closing).
%
%% 3) Morphological Gradient
%
% It is the difference between the dilation and the erosion of an image.
%
% $$dst = morph_{grad}( src, element )
%       = dilate( src, element ) - erode( src, element )$$
%
% It is useful for finding the outline of an object as can be seen below:
%
% <<https://docs.opencv.org/3.2.0/Morphology_2_Tutorial_Theory_Gradient.png>>
%
%% 4) Top Hat
%
% It is the difference between an input image and its opening.
%
% $$dst = tophat( src, element ) = src - open( src, element )$$
%
% <<https://docs.opencv.org/3.2.0/Morphology_2_Tutorial_Theory_TopHat.png>>
%
%% 5) Black Hat
%
% It is the difference between the closing and its input image
%
% $$dst = blackhat( src, element ) = close( src, element ) - src$$
%
% <<https://docs.opencv.org/3.2.0/Morphology_2_Tutorial_Theory_BlackHat.png>>
%
%% References
%
% * "Learning OpenCV" by Bradski and Kaehler.
%

function varargout = morphology_demo_gui2(im)
    % load source image
    if nargin < 1
        src = imread(fullfile(mexopencv.root(),'test','img001.jpg'));
    elseif ischar(im)
        src = imread(im);
    else
        src = im;
    end

    % create the UI
    h = buildGUI(src);
    if nargout > 0, varargout{1} = h; end
end

function onChange(~,~,h)
    %ONCHANGE  Event handler for UI controls

    % retrieve current values from UI controls
    ksize = round(get(h.slid, 'Value'));
    opIdx = get(h.pop(1), 'Value');
    shapeIdx = get(h.pop(2), 'Value');
    set(h.txt, 'String',sprintf('%2d',ksize));

    % options (structuring element shape and morphological operation)
    shapes = {'Rect', 'Cross', 'Ellipse'};
    types = {'Erode', 'Dilate', 'Open', 'Close', 'Gradient', 'Tophat', ...
        'Blackhat', ''};

    % apply the specified morphology operation
    elem = cv.getStructuringElement('Shape',shapes{shapeIdx}, ...
        'KSize',[ksize ksize]*2+1, 'Anchor',[ksize ksize]);
    if ~isempty(types{opIdx})
        dst = cv.morphologyEx(h.src, types{opIdx}, 'Element',elem);
    else
        dst = h.src;
    end

    % show result
    set(h.img, 'CData',dst);
    drawnow;
end

function h = buildGUI(img)
    %BUILDGUI  Creates the UI

    % parameters
    ksize = 3;
    max_ksize = 21;
    sz = size(img);
    sz(2) = max(sz(2), 250);  % minimum figure width

    % build the user interface (no resizing to keep it simple)
    h = struct();
    h.src = img;
    h.fig = figure('Name','Morphology Transformations Demo', ...
        'NumberTitle','off', 'Menubar','none', 'Resize','off', ...
        'Position',[200 200 sz(2) sz(1)+29]);
    if ~mexopencv.isOctave()
        %HACK: not implemented in Octave
        movegui(h.fig, 'center');
    end
    h.ax = axes('Parent',h.fig, 'Units','pixels', 'Position',[1 30 sz(2) sz(1)]);
    if ~mexopencv.isOctave()
        h.img = imshow(img, 'Parent',h.ax);
    else
        %HACK: https://savannah.gnu.org/bugs/index.php?45473
        axes(h.ax);
        h.img = imshow(img);
    end
    h.pop(1) = uicontrol('Parent',h.fig, 'Style','popupmenu', ...
       'Position',[5 5 70 20], 'Value',3, ...
       'String', {'Erosion','Dilation','Opening','Closing','Gradient',...
            'Top Hat','Black Hat','-None-'});
    h.pop(2) = uicontrol('Parent',h.fig, 'Style','popupmenu', ...
       'Position',[75 5 70 20], 'String',{'Rect','Cross','Ellipse'});
    h.slid = uicontrol('Parent',h.fig, 'Style','slider', 'Value',ksize, ...
        'Min',1, 'Max',max_ksize, 'SliderStep',[1 5]./(max_ksize-1), ...
        'Position',[150 5 sz(2)-150-20-5 20]);
    h.txt = uicontrol('Parent',h.fig, 'Style','text', 'FontSize',11, ...
        'Position',[sz(2)-20-5 5 20 20], 'String',num2str(ksize));

    % hook event handlers, and trigger default start
    set([h.slid, h.pop], 'Callback',{@onChange,h}, ...
        'Interruptible','off', 'BusyAction','cancel');
    onChange([],[],h);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [Aeq, beq] = getAbeq(n_seg, n_order, ts, start_cond, end_cond)
    n_coeff = n_order+1
    n_all_poly = n_seg*(n_coeff);
    %#####################################################
    % STEP 2.1 p,v,a constraint in start
    Aeq_start = zeros(3, n_all_poly);
    beq_start = zeros(3, 1);

    d1 = n_order;
    d2 = n_order * (n_order - 1);
    Aeq_start(1:1:3, 1:1:n_coeff) = [  1,    0,  0, 0, 0, 0, 0, 0;
                                       -d1,   d1,  0, 0, 0, 0, 0, 0;
                                       d2, -2*d2, d2, 0, 0, 0, 0, 0];
    beq_start = start_cond';% p,v,a

    %#####################################################
    % STEP 2.2 p,v,a constraint in end
    Aeq_end = zeros(3, n_all_poly);
    Aeq_end = zeros(3, 1);
    Aeq_end(1:1:3, n_all_poly-n_order:1:n_all_poly) = [0, 0, 0, 0, 0,  0,     0,   1;
                                                       0, 0, 0, 0, 0,  0,   -d1,  d1;
                                                       0, 0, 0, 0, 0, d2, -2*d2,  d2];
    beq_end = end_cond';% p,v,a

    %#####################################################
    % Init array for continuity constrains
    Aeq_con_p = zeros(n_seg-1, n_all_poly);
    beq_con_p = zeros(n_seg-1, 1);
    Aeq_con_v = zeros(n_seg-1, n_all_poly);
    beq_con_v = zeros(n_seg-1, 1);
    Aeq_con_a = zeros(n_seg-1, n_all_poly);
    beq_con_a = zeros(n_seg-1, 1);

    for con_index = 1:n_seg-1
        index = n_coeff + n_coeff * (con_index - 1);

        % STEP 2.3 position continuity constrain between 2 segments
        Aeq_con_p(con_index,index) = 1;% the end of previous segment
        Aeq_con_p(con_index,index+1) = -1;% the begin of next segment

        % STEP 2.4 velocity continuity constrain between 2 segments
        Aeq_con_v(con_index,index-1:index) = [-d1, d1];% the end of previous segment
        Aeq_con_v(con_index,index+1:index+2) = -[-d1, d1];% the begin of next segment

        % STEP 2.5 acceleration continuity constrain between 2 segments
        Aeq_con_a(con_index,index-2:index) = [d2, -2*d2, d2];% the end of previous segment
        Aeq_con_a(con_index,index+1:index+3) = -[d2, -2*d2, d2];% the begin of next segment
    end

    %#####################################################
    % combine all components to form Aeq and beq
    Aeq_con = [Aeq_con_p; Aeq_con_v; Aeq_con_a];
    beq_con = [beq_con_p; beq_con_v; beq_con_a];
    Aeq = [Aeq_start; Aeq_end; Aeq_con];
    beq = [beq_start; beq_end; beq_con];
end

author: Mesywang

repo: Motion-Planning-Algorithms


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%
% time invariant or time variant

%time variant
t=-5:.001:10;
p=heaviside(t)-heaviside(t-5);
y=t.*exp(-t).*p;
plot(t,y)
ylim([-.05 .4]);
legend('y(t)')

figure
plot(t+3,y)
ylim([-.05 .4]);
legend('y(t-3)')

figure
t=-5:.001:10;
p=heaviside(t-3)-heaviside(t-8);
y2=t.*exp(-t).*p;
plot(t,y2)
ylim([-.01 .2]);
legend('S[x(t-3)]')

%time invariant
figure
t=-5:.01:20;
p=heaviside(t-1)-heaviside(t-11);
y=1-2*cos(t-1).*p;
plot(t,y)
legend('y(t)')

figure
p=heaviside(t)-heaviside(t-10);
x=cos(t).*p;
plot(t+1,1-2*x)
legend('y(t)')

figure
plot(t+4,y)
legend('y(t-4)')

figure
p=heaviside(t-4)-heaviside(t-14);
x=cos(t-4).*p;
plot(t,x)

figure
p=heaviside(t-4)-heaviside(t-14);
x=cos(t-4).*p;
plot(t+1,1-2*x)

figure
t=-5:.01:20;
p=heaviside(t-5)-heaviside(t-15);
y2=1-2*cos(t-5).*p ;
plot(t,y2) ;
legend('S[x(t-4)]')

%time variant
figure
t=-5:.1:10;
x=heaviside(t+2)-heaviside(t-2);
plot(t,x);
ylim([-.1 1.1]);

figure
plot((1/2)*t,x);
ylim([-.1 1.1]);

figure
plot((1/2)*t+2,x);
ylim([-.1 1.1]);
legend('y_1(t)')

figure
plot(t+2,x);
ylim([-.1 1.1]);
legend('x(t-2)')

figure
t=-5:.1:10;
x2=heaviside(t)-heaviside(t-4);
plot(t,x2);
ylim([-.1 1.1]);

figure
plot((1/2)*t,x2);
ylim([-.1 1.1]);
legend('y_2(t)')

%shift invariant
figure
n=0:5;
x=0.8.^n;
y=x.^2;
stem(n,y);
xlim([-.1 5.1])
legend('y[n]')

figure
stem(n+2,y)
legend('y[n-2]')
xlim([1.9 7.1])

figure
n=2:7;
y2=(0.8.^(n-2)).^2;
stem(n,y2);
xlim([1.9 7.1])
legend('S[x[n-2]]')


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
function stats = eval_labels(labels_folder,database,gt_set,n_cands)

if nargin<2
    database = 'pascal2012';
end
if nargin<3
    gt_set = 'val2012';
end
if nargin<4
    % Number of sampled candidates
    n_cands = [10:5:100,125:25:1000,1500:500:6000,10000];
end

% Get the name of the folder to refer to the method
if strcmp(labels_folder(end),filesep)
    labels_folder(end) = [];
end
tmp = strfind(labels_folder,filesep);
if isempty(tmp)
    method_name = labels_folder;
else
    method_name = labels_folder(tmp(end)+1:end);
end

res_dir    = fullfile(root_dir,'results',database, method_name);
stats_file = fullfile(root_dir,'results',database, [method_name '_' gt_set '.mat']);

% Is the result already gathered?
if exist(stats_file, 'file')
    load(stats_file)
    recompute = 0;
    disp(['Loaded: ' stats_file '.'])
else
    disp(['RECOMPUTING: ' stats_file '.'])
    recompute = 1;
end

if recompute
    %% Evaluate and save each image independently to be able to parallelize
    % You can adapt the matlabpool to your system
    eval_and_save_labels(labels_folder,database,gt_set);

    %% Gather and save results
    % Load which images to consider
    im_ids = database_ids(database,gt_set);

    % Store and initialize
    stats.n_cands = n_cands;
    stats.gt_set  = gt_set;
    stats.obj_classes = [];

    % Compute statistics
    stats.num_objects = 0;
    for ii=1:numel(im_ids)
        res_file = fullfile(res_dir, [im_ids{ii} '.mat']);
        if exist(res_file,'file')
            load(res_file)
            for jj=1:length(stats.n_cands)
                curr_n_cands = min(stats.n_cands(jj), size(jaccards,2));
                to_consider = 1:curr_n_cands;
                stats.all_n_masks(ii,jj) = length(to_consider);
                if (stats.all_n_masks(ii,jj)>0)
                    for kk=1:size(jaccards,1)
                        [stats.max_J(stats.num_objects+kk,jj), which_one] = max(jaccards(kk,to_consider));
                        stats.max_indicator(stats.num_objects+kk,jj) = to_consider(which_one);
                        stats.max_fp(stats.num_objects+kk,jj) = false_pos(kk,to_consider(which_one));
                        stats.max_fn(stats.num_objects+kk,jj) = false_neg(kk,to_consider(which_one));
                        stats.max_inters(stats.num_objects+kk,jj) = inters(kk,to_consider(which_one));
                    end
                end
            end
            stats.obj_classes = [stats.obj_classes; obj_classes(1:size(inters,1))];
            stats.num_objects = stats.num_objects + size(jaccards,1);
        else
            error([res_file ' not found']);
        end
    end

    % Check
    if size(stats.obj_classes,1)~=stats.num_objects
        error('Something went wrong')
    end

    stats.mean_n_masks   = mean(stats.all_n_masks);

    % ----- Jaccard at instance level (J_i) ----
    % It is the mean best jaccard for all objects
    stats.jaccard_object = mean(stats.max_J);

    % ----- Compute jaccard at class level (J_c) ----
    class_ids = unique(stats.obj_classes);
    for ii=1:length(class_ids)
        curr_class = class_ids(ii);
        stats.per_class_results{ii}.num_objects = sum(stats.obj_classes==curr_class);

        stats.per_class_results{ii}.max_fp     = stats.max_fp(logical(stats.obj_classes==curr_class),:);
        stats.per_class_results{ii}.max_fn     = stats.max_fn(logical(stats.obj_classes==curr_class),:);
        stats.per_class_results{ii}.max_inters = stats.max_inters(logical(stats.obj_classes==curr_class),:);
        stats.per_class_results{ii}.max_J      = stats.max_J(logical(stats.obj_classes==curr_class),:);
        stats.per_class_results{ii}.meanmax    = mean(stats.per_class_results{ii}.max_J);

        stats.per_class_results{ii}.global_fp      = sum(stats.per_class_results{ii}.max_fp,1);
        stats.per_class_results{ii}.global_fn      = sum(stats.per_class_results{ii}.max_fn,1);
        stats.per_class_results{ii}.global_inters  = sum(stats.per_class_results{ii}.max_inters,1);

        % Compute per-class total inters, fp, fn
        stats.per_class_results{ii}.global_J = ...
             stats.per_class_results{ii}.global_inters ./...
            (stats.per_class_results{ii}.global_inters+...
             stats.per_class_results{ii}.global_fp    +...
             stats.per_class_results{ii}.global_fn);
    end

    % Compute global mean on all classes
    tmp = [];
    for ii=1:length(class_ids)
        tmp = [tmp; stats.per_class_results{ii}.global_J]; %#ok<AGROW>
    end
    stats.jaccard_class = mean(tmp,1);

    save(stats_file,'stats');

    % Remove temporal results folder
    rmdir(res_dir,'s');
end
end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function [nv,ev,uv] = AltVelCalc(msg)

% Copyright 2010, The MathWorks, Inc.

%AltVelCalc
% msg = '8DABDEEC99153E09802C013F4BDF';
% msg = '8DA66A1399104CAA80A406129D1D';
% msg = '8DAA8A83991502A7000411EE1091';
% msg = '8D40067899050CA680050D26CFB1';
% msg = '8D4CA74E9914C0AC80040F8D6965';
% msg = '8D80043C9904E4A7A0070D597008';
% msg = '8DABC6519904AFAC000513E3A869';
% msg = '8D4B187E9944CDAAA8040F67994E';
% msg = '8DA6CC41990488AE00070049848E';
% msg = '8DA9BF4B99948BB0A0040E5828F0';
% msg = '8D40067899050CA680050D26CFB1';
% msg = '8DA3CC3790C380976B152295D11F';  % Bad results
% msg = '8DA3CC37994483ACB004003CB62A';

a = hex2dec(msg');
b = dec2bin(a);
bin = reshape(b',1,length(msg)*4);

aircraftID = msg(3:8);

ewDir = bin(46);
if ewDir == '0'
    EW = 'East';
    ed = 1;
else
    EW = 'West';
    ed = -1;
end
ewVel = bin2dec(bin(47:56))-1;

nsDir = bin(57);
if nsDir == '0'
    NS = 'North';
    nd = 1;
else
    NS = 'South';
    nd = -1;
end
nsVel = bin2dec(bin(58:67))-1;

udDir = bin(69);
if udDir == '0'
    UD = 'Up';
    ud = 1;
else
    UD = 'Down';
    ud = -1;
end
udVel = (bin2dec(bin(70:78))-1)*64;

speed = sqrt(ewVel^2+nsVel^2);   % Speed in knots

disp(sprintf('Aircraft ID %s is traveling at %f knots\nDirection %s at %f knots, direction %s at %f knots ', aircraftID, speed, EW, ewVel, NS, nsVel));
disp(sprintf('Aircraft ID %s is going %s at %f feet/min\n', aircraftID, UD, udVel));

nv = nd*nsVel;
ev = ed*ewVel;
uv = ud*udVel;

author: analogdevicesinc


text:
function K = linear(kern,dat1,dat2,ind1,ind2,kerParam),

K = get_x(dat2,ind2)*get_x(dat1,ind1)';

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function varargout = sameDeltaLocs(f, g, tol)
%SAMEDELTALOCS   Find common locations of delta functions in F and G.
%   IDX1 = SAMEDELTALOCS(F, G) returns a logical vector of of same length as
%   that of F.DELTALOCS indicating which locations are shared with G.
%   [IDX1, IDX2] = SAMEDELTALOCS(F, G) does the same but also returns the
%   vector IDX2 locations of G.
%
%   SAMEDELTALOCS(F, G, TOL) uses the user-specified tolerance TOL for testing.
%
% See also SIMPLIFY, CLEANROWS, CLEANCOLUMNS

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

% Deal with the empty cases:
if ( isempty(f) || isempty(g) )
    if ( nargout == 2 )
        varargout{1} = [];
        varargout{2} = [];
    else
        varargout{1} = [];
    end
    return
end

% Get preferences:
if ( nargin < 3 )
    pref = chebfunpref();
    tol = pref.deltaPrefs.proximityTol;
end

% Initialize variables:
loc1 = f.deltaLoc;
idx1 = zeros(size(loc1));
loc2 = g.deltaLoc;
idx2 = zeros(size(loc2));

% Loop and compare locations:
for i = 1:length(loc1)
    for j = 1:length(loc2)
        % If both entries are equal to zero:
        p = ( (loc1(i) == 0) && (loc2(j) == 0) );
        % Or if they are very close to each other:
        maxLoc = max(abs([loc1(i), loc2(j)]));
        p = p | ( abs((loc1(i) - loc2(j)))/maxLoc < tol );
        if ( p )
            % Mark the indices:
            idx1(i) = 1;
            idx2(j) = 1;
        end
    end
end

% Convert to logical indices:
idx1 = logical(idx1);
idx2 = logical(idx2);

% Output based on the number of outputs:
if ( nargout == 2 )
    varargout{1} = idx1;
    varargout{2} = idx2;
else
    varargout{1} = idx1;
end

end


text:
 function st = nufft_init(om, Nd, Jd, Kd, varargin)
%function st = nufft_init(om, Nd, Jd, Kd, [n_shift,] ...)
%
% Initialize structure for d-dimension NUFFT using KB interpolator,
% particularly the interpolation matrix in sparse format.
% caution: this routine can require a lot of memory!
% in
%	om [M,d]	"digital" frequencies in radians
%	Nd [d]		image dimensions (N1,N2,...,Nd)
%	Jd [d]		# of neighbors used (in each direction)
%	Kd [d]		FFT sizes (should be >= N1,N2,...)
% optional arguments
%	n_shift [d]	n = 0-n_shift to N-1-n_shift (must be first)
%	'minmax:kb'	minmax interpolator with excellent KB scaling!
%				(minmax:kb is recommended, and used by default)
%	'minmax:tuned'	minmax interpolator, somewhat numerically tuned
%	'minmax:user'	minmax interpolator with user ({alpha}, {beta})
%	'uniform'	uniform scaling factors (not recommended)
%	'kaiser'	kaiser-bessel (KB) interpolator (minmax best alpha, m)
%			or 'kaiser', alpha, m to specify parameter (vectors)
%	'linear'	linear interpolator (a terrible straw man)
%	kernel		user-provided inline interpolation kernel(k,J)
%			(or a cell array of kernels, one for each dimension)
%	'table'		use table-based interpolation rather than sparse matrix.
%			this can save a lot of memory for large problems.
%			example ..., 'table', 2^11, 'minmax:kb'
%			where 2^11 is the table over-sampling factor.
% out
%	st.p		[M, *Nd]	sparse interpolation matrix
%					(or empty if table-based)
%	st.sn		[(Nd)]		scaling factors
%	st.Nd,Jd,Kd,om	copies of inputs
%
% *Nd is shorthand for prod(Nd).
% (Nd) is shorthand for (N1,N2,...,Nd)
%
% Like fft(), the NUFFT expects the signals to be x(0,0), ...
% Use n_shift = [N1/2, N2/2, ...] for x(-N1/2,-N2/2,...), ...
%
% Copyright 2002-5-30	Jeff Fessler	The University of Michigan

if nargin < 4, help(mfilename), error args, end

% dimensionality of input space (usually 2 or 3)
dd = length(Nd);
if dd ~= length(Jd) | dd ~= length(Kd)
	error 'inconsistent dim'
end
if dd ~= size(om,2), error(sprintf('omega needs %d columns', dd)), end

%
% process optional arguments
%

% n_shift argument? (must be first)
if length(varargin) > 0 & isnumeric(varargin{1})
	n_shift = varargin{1};
	if dd ~= length(n_shift)
		error(sprintf('n_shift needs %d columns', dd))
	end
	varargin = {varargin{2:end}};
else
	n_shift = zeros(size(Nd));
end
st.n_shift = n_shift;

% default/recommended interpolator is minmax with KB scaling factors
if length(varargin) == 0
	varargin = {'minmax:kb'};
end

st.alpha = {};
st.beta = {};
is_kaiser_scale = logical(0);

% table based?
if ischar(varargin{1}) & streq(varargin{1}, 'table')
	st = nufft_table_init(om, Nd, Jd, Kd, n_shift, varargin{2:end});
	return
end

ktype = varargin{1};

% cell array of kernel functions: {kernel1, kernel2, ..., kernelD}
if isa(ktype, 'cell')
	if isa(ktype{1}, 'inline') | isa(ktype{1}, 'function_handle')
		ktype = 'inline';
		if length(varargin) > 1, error 'excess arguments?', end
		if length(varargin{1}) ~= dd, error 'wrong # of kernels', end
		st.kernel = varargin{1};
	else
		error 'cell array should be inline kernels!?'
	end

% or a single inline kernel for all dimension
elseif isa(ktype, 'inline') | isa(ktype, 'function_handle')
	ktype = 'inline';
	if length(varargin) > 1, error 'excess arguments?', end
	for id = 1:dd
		st.kernel{id} = varargin{1};	% all same
	end

% or a string that describes the type of interpolator
elseif ~ischar(ktype)
	error 'non-string kernel type?'

end

st.ktype = ktype;

%
% set up whatever is needed for each interpolator
%
if streq(ktype, 'inline')
	% already did it above

% linear interpolator straw man
elseif streq(ktype, 'linear')
	ktype = 'inline';
	kernel = inline('(1 - abs(k/(J/2))) .* (abs(k) < J/2)', 'k', 'J');
	for id = 1:dd
		st.kernel{id} = kernel;
	end

% KB interpolator
elseif streq(ktype, 'kaiser')
	is_kaiser_scale = logical(1);

	% with minmax-optimized parameters
	if length(varargin) == 1
		for id = 1:dd
			[st.kernel{id} st.kb_alf(id) st.kb_m(id)] = ...
				kaiser_bessel('inline', Jd(id));
		end

	% with user-defined parameters
	elseif length(varargin) == 3
		alpha_list = varargin{2};
		m_list = varargin{3};
		if (length(alpha_list) ~= dd) | (length(m_list) ~= dd)
			error 'need #dim alpha and m'
		end
		for id = 1:dd
			[st.kernel{id} st.kb_alf(id) st.kb_m(id)] = ...
				kaiser_bessel('inline', Jd(id), ...
					alpha_list(id), m_list(id));
		end
	else
		error 'kaiser should have no arguments, or both alpha and m'
	end

% minmax interpolator with KB scaling factors (recommended default)
elseif streq(ktype, 'minmax:kb')
	for id = 1:dd
		[st.alpha{id}, st.beta{id}] = ...
			nufft_alpha_kb_fit(Nd(id), Jd(id), Kd(id));
	end

% minmax interpolator with numerically "tuned" scaling factors
elseif streq(ktype, 'minmax:tuned')
	for id = 1:dd
		[st.alpha{id}, st.beta{id}, ok] = ...
			nufft_best_alpha(Jd(id), 0, Kd(id)/Nd(id));
		if ~ok, error 'unknown J,K/N', end
	end

% minmax interpolator with user-provided scaling factors
elseif streq(ktype, 'minmax:user')
	if length(varargin) ~= 3, error 'user must provide alpha/beta', end
	st.alpha = varargin{2};
	st.beta = varargin{3};
	if length(st.alpha) ~= dd | length(st.beta) ~= dd
		error 'alpha/beta size mismatch'
	end

elseif streq(ktype, 'uniform')
	for id = 1:dd
		st.alpha{id} = 1;
		st.beta{id} = 0;
	end

else
	error 'unknown kernel type'
end

st.tol	= 0;

st.Jd	= Jd;
st.Nd	= Nd;
st.Kd	= Kd;

M = size(om,1);
st.M	= M;
st.om	= om;

%
% scaling factors: "outer product" of 1D vectors
%
st.sn = 1;
for id=1:dd
	if is_kaiser_scale
		nc = [0:Nd(id)-1]'-(Nd(id)-1)/2;
		tmp = 1 ./ kaiser_bessel_ft(...
			nc/Kd(id), Jd(id), st.kb_alf(id), st.kb_m(id), 1);
	elseif streq(ktype, 'inline')
		tmp = 1 ./ nufft_interp_zn(0, Nd(id), Jd(id), Kd(id), st.kernel{id});
	else
		tmp = nufft_scale(Nd(id), Kd(id), st.alpha{id}, st.beta{id});
	end
	st.sn = st.sn(:) * tmp';
end
if length(Nd) > 1
	st.sn = reshape(st.sn, Nd);	% [(Nd)]
else
	st.sn = st.sn(:);	% [(Nd)]
end

%
% [J?,M] interpolation coefficient vectors.  will need kron of these later
%
for id=1:dd
	N = Nd(id);
	J = Jd(id);
	K = Kd(id);
	if isvar('st.kernel')
		[c, arg] = ...
		nufft_coef(om(:,id), J, K, st.kernel{id});	% [J?,M]
	else
		alpha = st.alpha{id};
		beta = st.beta{id};
		T = nufft_T(N, J, K, st.tol, alpha, beta);	% [J?,J?]
		[r, arg] = ...
		nufft_r(om(:,id), N, J, K, alpha, beta);	% [J?,M]
		c = T * r;	clear T r
	end

	gam = 2*pi/K;
	phase_scale = 1i * gam * (N-1)/2;

	phase = exp(phase_scale * arg);	% [J?,M] linear phase
	ud{id} = phase .* c;		% [J?,M]

	%
	% indices into oversampled FFT components
	%
	koff = nufft_offset(om(:,id), J, K);	% [M,1] to leftmost near nbr
	kd{id} = mod(outer_sum([1:J]', koff'), K) + 1;	% [J?,M] {1,...,K?}
	if id > 1	% trick: pre-convert these indices into offsets!
		kd{id} = (kd{id}-1) * prod(Kd(1:(id-1)));
	end

end, clear c arg gam phase phase_scale koff N J K

%
% build sparse matrix that is [M,*Kd]
% with *Jd nonzero entries per frequency point
%
if dd >= 3, printf('Needs at least %g Gbyte RAM', prod(Jd)*M*8/10^9*2), end

kk = kd{1};	% [J1,M]
uu = ud{1};	% [J1,M]
for id = 2:dd
	Jprod = prod(Jd(1:id));
	kk = block_outer_sum(kk, kd{id});	% outer sum of indices
	kk = reshape(kk, Jprod, M);
	uu = block_outer_prod(uu, ud{id});	% outer product of coefficients
	uu = reshape(uu, Jprod, M);
end	% now kk and uu are [*Jd, M]

%
% apply phase shift
% pre-do Hermitian transpose of interpolation coefficients
%
phase = exp(1i * (om * n_shift(:))).';			% [1,M]
uu = conj(uu) .* phase(ones(1,prod(Jd)),:);		% [*Jd,M]

mm = [1:M]; mm = mm(ones(prod(Jd),1),:);		% [*Jd,M]
st.p = sparse(mm(:), kk(:), uu(:), M, prod(Kd));	% sparse matrix

%
% in
%	x1	[J1,M]
%	x2	[J2,M]
% out
%	y	[J1,J2,M]	y(i1,i2,m) = x1(i1,m) + x2(i2,m)
%
function y = block_outer_sum(x1, x2)
[J1 M] = size(x1);
[J2 M] = size(x2);
xx1 = reshape(x1, [J1 1 M]);	% [J1,1,M] from [J1,M]
xx1 = xx1(:,ones(J2,1),:);	% [J1,J2,M], emulating ndgrid
xx2 = reshape(x2, [1 J2 M]);	% [1,J2,M] from [J2,M]
xx2 = xx2(ones(J1,1),:,:);	% [J1,J2,M], emulating ndgrid
y = xx1 + xx2;			% [J1,J2,M]

function y = block_outer_prod(x1, x2)
[J1 M] = size(x1);
[J2 M] = size(x2);
xx1 = reshape(x1, [J1 1 M]);	% [J1,1,M] from [J1,M]
xx1 = xx1(:,ones(J2,1),:);	% [J1,J2,M], emulating ndgrid
xx2 = reshape(x2, [1 J2 M]);	% [1,J2,M] from [J2,M]
xx2 = xx2(ones(J1,1),:,:);	% [J1,J2,M], emulating ndgrid
y = xx1 .* xx2;			% [J1,J2,M]

author: thomaskuestner


text:
function D = meeg(varargin)
% Function for creating meeg objects.
% FORMAT m = meeg(varargin)
%
% SPM8 format consists of a binary memory - mapped data file and a header
% mat file.
%
% The header file will contain an object of class meeg called D. All
% information other than data is contained in this object and access to the
% data is via methods of the object. Also, arbitrary data can be stored
% inside the object if their field names do not conflict with existing
% methods' names.
%
% The following is a description of the internal implementation of meeg.
%
% Fields of meeg:
% .type - type of data in the file: 'continuous', 'single', 'evoked'
% .Fsample - sampling rate
% .data - substruct containing the information related to the bimary data
% file and memory mapping.
%
%   Subfields of .data
%       .y - reference to the data file_array
%       .fnamedat - name of the .dat file
%       .datatype - type of the data for file_array
%
%
% .Nsamples - length of the trial (whole data if the file is continuous).
% .timeOnset - the peri-stimulus time of the first sample in the trial (in sec)
%
% .fname, .path - strings added by spm_eeg_load to keep track of where a
%                 header struct was loaded from.
%
% .trials - this describes the segments of the epoched file and is also a
%           structure array.
%
%   Subfields of .trials
%
%       .label - user-specified string for the condition
%       .onset - time of the first sample in seconds in terms of the
%                original file
%       .bad - 0 or 1 flag to allow rejection of trials.
%       .repl - for epochs that are averages - number of replications used
%               for the average.
%       .events - this is a structure array describing events related to
%                 each trial.
%
%           Subfields of .events
%
%           .type - string (e.g. 'front panel trigger')
%           .value - number or string, can be empty (e.g. 'Trig 1').
%           .time - in seconds in terms of the original file
%           .duration - in seconds
%
% .channels - This is a structure array which is a field of meeg.
%             length(channels) should equal size(.data.y, 1) and the order
%             must correspond to the order of channels in the data.
%
%   Subfields of .channels
%
%       .label - channel label which is always a string
%       .type - a string, possible values - 'MEG', 'EEG', 'VEOG', 'HEOG',
%               'EMG' ,'LFP' etc.
%       .units - units of the data in the channel.
%       .bad - 0 or 1 flag to mark channels as bad.
%       .X_plot2D, .Y_plot2D - positions on 2D plane (formerly in ctf). NaN
%                              for channels that should not be plotted.
%
% .sensors
%
%
%   Subfields of .sensors (optional)
%       .meg - struct with sensor positions for MEG (subfields: .pnt .ori .tra .label)
%       .eeg - struct with sensor positions for MEG (subfields: .pnt .tra .label)
%
% .fiducials - headshape and fiducials for coregistration with sMRI
%
%   Subfiels of .fiducials (optional)
%       .pnt - headshape points
%       .fid.pnt - fiducial points
%       .fid.label - fiducial labels
%
% .artifacts - structure array with fields .start and .stop expressed in
%              seconds in original file time.
%
% .history - structure array describing commands that modified the file.
%
%   Subfields of .history:
%
%       .function - string, the function name
%       .arguments - cell array, the function arguments
%       .time - when function call was made
%
% .other - structure used to store other information bits, not fitting the
%          object structure at the moment,
%       for example:
%       .inv - structure array corresponding to the forw/inv problem in MEEG.
%       .ival - index of the 'inv' solution currently used.
%
% ______________________________________________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: meeg.m 3940 2010-06-21 12:51:19Z christophe $

if nargin == 1
    if isstruct(varargin{1})
        [OK D] = checkmeeg(varargin{1}, 'basic');
        if OK
            D  = class(D, 'meeg');
        else
            error('Struct not fit for conversion to meeg');
        end
        return;
    elseif isa(varargin{1},'meeg'),
        D = varargin{1};
        return;
    end;
end;

D          = [];
D.Nsamples = 0;
[res D]    = checkmeeg(D);
D          = class(D, 'meeg');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [params, result] = ortho_default(params, W, w)
% Default DSS orthogonalization function
%   W = orthof(W)     for symmetric dss
%   w = orthof(W, w)  for deflation dss
%     W Matrix with projection vectors as rows. For deflation
%       algorithm only previously calculated projections are given.
%     w Currently iterated projection. Only for deflation algorithm.

% Copyright (C) 2004, 2005 DSS MATLAB package team (dss@cis.hut.fi).
% Distributed by Laboratory of Computer and Information Science,
% Helsinki University of Technology. http://www.cis.hut.fi/projects/dss/.
% $Id$

if nargin<2
    params.name = 'Default orthogonalization';
    params.description = 'Description of this function.';
    return;
end

if nargin>2
  % per component orthogonalization
  w = w - W' * W * w;
  result = w / norm(w);
else
  % symmetric orthogonalization
  W = real(inv(W * W')^(1/2))' * W;
  result = W;
end


text:
function Zvals = cp_wfg_sparse_setup(Z,W)
%CP_WFG_SPARSE_SETUP Creates a special array.
%
%   ZVALS = CP_WFG_SPARSE_SETUP(Z,W) creates an array ZVALS that
%   contains the values of Z corresponding to the indices specified
%   by W.subs.
%
%   See also CP_WFG_SPARSE.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

Zsubs = Z.subs;
Wsubs = W.subs;
Zvals = zeros(size(W.vals));
[junk,loc] = ismember(Zsubs,Wsubs,'rows');
Zvals(loc) = Z.vals;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function outHtml = colText(inText, inColor)
% Returns a HTML string with colored font. Used for GUIs.
% Called by `metCompareGUI`, `reactionCompareGUI`.
%
% USAGE:
%
%    outHtml = colText(inText, inColor)
%
% INPUTS:
%    inText:      String
%    inColor      Color to make string.
%
% OUTPUTS:
%    outHtml:     Text with formatting
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

outHtml = ['', ...
    inText, ...
    ''];

author: opencobra

repo: cobratoolbox


text:
function [camPoints, prjPoints, Nodes, Edges] = getMatchedNodes(whiteLight, colorGrid, camCorners, prjW, prjH, verbose)
%% Get matched Nodes on both camera and projector image.
% Tthe outputs are matched camera and projector 2D point pairs.

%% License
% ACADEMIC OR NON-PROFIT ORGANIZATION NONCOMMERCIAL RESEARCH USE ONLY
% Copyright (c) 2018 Bingyao Huang
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:

% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.

% If you publish results obtained using this software, please cite our paper.

% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
% SOFTWARE.

%% Step 1. Extract grid mask on the white board
[imBWGrid, imColorGridMasked] = ImgProc.segColorGrid(whiteLight, colorGrid, camCorners, verbose);
imSize = size(imBWGrid);

%% Step 2. Skeletonize
imSkele = ImgProc.bw2skele(imBWGrid);

if (verbose)
    figure;
    imshowpair(imBWGrid, imSkele, 'Montage');
    title('Grid mask and skeletonized grid mask');
end

%% Step 3. Extract edge, node endpoint and # neighbors image
[imEdge, imNode, ~, ~] = extractStructs(imSkele);

%% Step 4. Pre-extract Nodes and Edges structures for cleaning
[Edges, Nodes] = ImgProc.skeleToStruct(imNode, imEdge);

% get rid unnormally long edges
    hEdge = Edges(logical([Edges.isH]));
    hEdge = hEdge([hEdge.Area] < mean([hEdge.Area]) + 2*std([hEdge.Area]));

    vEdge = Edges(~logical([Edges.isH]));
    vEdge = vEdge([vEdge.Area] < mean([vEdge.Area]) + 2*std([vEdge.Area]));

    Edges = [hEdge; vEdge];

% find isolated edges (area == 1) and remove those edges from imEdge
% and add to imNodes
labels = [Edges.Area] == 1;
isoEdges = Edges(labels);
vecNumNodes = cellfun(@numel, {isoEdges.nodes});
leafEdgeIdx = vecNumNodes == 1;
nodeEdgeIdx = vecNumNodes == 2;

imEdge([isoEdges(nodeEdgeIdx).PixelIdxList]) = 0;
imNode([isoEdges(nodeEdgeIdx).PixelIdxList]) = 1;

% Extract Nodes and Edges from cleaned imEdge and imNode
[Edges, Nodes] = ImgProc.skeleToStruct(imNode, imEdge);
numNodes = numel(Nodes);
numEdges = numel(Edges);

% imMix = imEdge + imNode*3;
% figure;
% title('imMix - Edges and Nodes');
% imagesc(imMix);

%% step 5. create labeled edge and node images for debugging

if (verbose)
%     imEdgeLabel = zeros(size(imNode), 'uint16');
%     imNodeLabel = zeros(size(imNode), 'uint16');
%
%     for i = 1:numEdges
%         imEdgeLabel(Edges(i).PixelIdxList) = i;
%     end
%
%     for i = 1:numNodes
%         imNodeLabel(Nodes(i).PixelIdxList) = i;
%     end
%
%     subplot(2,2,3);drawnow
%     title('imEdgeLabel - Labeled Edges');
%     imagesc(imEdgeLabel);
%     caxis([0 0.1]);
%     myMap = [0, 0, 0; 1, 1, 1];
%     colormap(myMap);
%
%     subplot(2,2,4);drawnow
%     title('imNodeLabel - Labeled Nodes');
%     imagesc(imNodeLabel);
%     caxis([0 0.1]);
%     myMap = [0, 0, 0; 1, 1, 1];
%     colormap(myMap);
end

%% step 6. create the ajacency matrix and graph from Edges and Nodes
A = zeros(numNodes, numNodes);

    for i = 1:numEdges
    curEdge = Edges(i);

    if (numel(curEdge.nodes) == 2)% ignore end point and end edge
        edgeLength = numel(curEdge.PixelIdxList);
        A(curEdge.nodes(1), curEdge.nodes(2)) = edgeLength;
        A(curEdge.nodes(2), curEdge.nodes(1)) = edgeLength;
    end

end

%% step 6-7. create the ajacency matrix and graph from Edges and Nodes and traverse the grid to assign four neighbors to Nodes
Nodes = ImgProc.traverseGrid(Nodes, Edges);

%% step 8. get horizontal and vertical edges
horiEdges = Edges([Edges.isH] == 1);
vertEdges = Edges([Edges.isH] == 0);

imHoriEdge = false(imSize);
imHoriEdge(vertcat(horiEdges.PixelIdxList)) = 1;

imVertEdge = false(imSize);
imVertEdge(vertcat(vertEdges.PixelIdxList)) = 1;

%% step 9. read color grid and convert RGB to labels
[imAllLabel, imHoriLabel, imVertLabel] = ImgProc.colorToLabel(imColorGridMasked, imNode, imHoriEdge, imVertEdge, verbose);

%% step 10. assign each node horizontal and vertical color labels
[Nodes(:).horiColor] = deal(-1);
[Nodes(:).vertColor] = deal(-1);

% imNodeWide = imdilate(imNode, strel('disk', 1));
imNodeWide = imdilate(imNode, strel('square', 3));

imHoriWide = imdilate(imHoriEdge, strel('disk', 1));
imHoriWide = logical(imHoriWide .* (~imNodeWide));

imVertWide = imdilate(imVertEdge, strel('disk', 1));
imVertWide = logical(imVertWide .* (~imNodeWide));

% labeled images (faster than imreconstruct)
imHoriWideL = bwlabel(imHoriWide);
imVertWideL = bwlabel(imVertWide);

for i = 1:numel(Nodes)
    hEdges = [Nodes(i).hEdges];
    hEdges = hEdges(hEdges > 0);
    curHoriEdgePixIdx = vertcat(Edges(hEdges).PixelIdxList);

%     imMarker = false(imSize);
%     imMarker(curHoriEdgePixIdx) = 1;
%     imWideRecon = imreconstruct(imMarker, imHoriWide, 4);
%     curHoriEdgeLabels = imHoriLabel(imWideRecon > 0);

    hLabels = imHoriWideL(curHoriEdgePixIdx);
    idx = ismember(imHoriWideL, hLabels(hLabels>0));
    curHoriEdgeLabels = imHoriLabel(idx);

    % vertical
    vEdges = [Nodes(i).vEdges];
    vEdges = vEdges(vEdges > 0);
    curVertEdgePixIdx = vertcat(Edges(vEdges).PixelIdxList);

%     imMarker = false(imSize);
%     imMarker(curVertEdgePixIdx) = 1;
%     imWideRecon = imreconstruct(imMarker, imVertWide, 4);
%     curVertEdgeLabels = imVertLabel(imWideRecon > 0);

    vLabels = imVertWideL(curVertEdgePixIdx);
    idx = ismember(imVertWideL, vLabels(vLabels>0));
    curVertEdgeLabels = imVertLabel(idx);

    % vote majority
    Nodes(i).horiColor = mode(curHoriEdgeLabels);
    Nodes(i).vertColor = mode(curVertEdgeLabels);
end

%% step 11. correct each Node's label by majority voting on each stripe

% NodesBackup = Nodes; % keep a copy of old Nodes for debug
Nodes = ImgProc.correctColors(Nodes);

%% step 12. decode the De Bruijn sequence

Nodes = ImgProc.decodeDebruijn(Nodes, prjW, prjH);
[camPoints, prjPoints] = extractNodesPosition(Nodes);

% fill missing coordinates by majority vote on the same stripe
for i = 1:numel(Nodes)
    % find all nodes on the same horizontal line
    Nnb = ImgProc.findNbInDir(Nodes, i, 'N', 0);
    Snb = ImgProc.findNbInDir(Nodes, i, 'S', 0);
    Enb = ImgProc.findNbInDir(Nodes, i, 'E', 0);
    Wnb = ImgProc.findNbInDir(Nodes, i, 'W', 0);

    % activeRow
    hNodeIdx = [flip(Wnb), i, Enb];
    activeRow = [Nodes(hNodeIdx).activeRow];
    if(nnz(activeRow > 0))
        [Nodes(hNodeIdx).activeRow] = deal(mode(activeRow(activeRow > 0)));
    end

    % activeCol
    vNodeIdx = [flip(Nnb), i, Snb];
    activeCol = [Nodes(vNodeIdx).activeCol];
    if(nnz(activeCol > 0))
        [Nodes(vNodeIdx).activeCol] = deal(mode(activeCol(activeCol > 0)));
    end
end

% fill missing node prj coordinates using linear interpolation
for i = 1:numel(Nodes)
    % find all nodes on the same horizontal line
    Nnb = ImgProc.findNbInDir(Nodes, i, 'N',0);
    Snb = ImgProc.findNbInDir(Nodes, i, 'S',0);
    Enb = ImgProc.findNbInDir(Nodes, i, 'E',0);
    Wnb = ImgProc.findNbInDir(Nodes, i, 'W',0);

    hNodeIdx = [flip(Wnb), i, Enb];
%     hStripes = [Nodes(hNodeIdx).hEdges];

    vNodeIdx = [flip(Nnb), i, Snb];
%     vStripes = [Nodes(vNodeIdx).vEdges];

    % fill -1 with linear interpolation

    % horizontal stripes
    vecCol = [Nodes(hNodeIdx).activeCol];
    vecCol = num2cell(fillmissing(vecCol, 'linear', 'MissingLocations', vecCol<0));
    [Nodes(hNodeIdx).activeCol] = vecCol{:};

    vecRow = [Nodes(hNodeIdx).activeRow];
    vecRow = num2cell(fillmissing(vecRow, 'linear', 'MissingLocations', vecRow<0));
    [Nodes(hNodeIdx).activeRow] = vecRow{:};

    % vertical stripes
    vecRow = [Nodes(vNodeIdx).activeRow];
    vecRow = num2cell(fillmissing(vecRow, 'linear', 'MissingLocations', vecRow<0));
    [Nodes(vNodeIdx).activeRow] = vecRow{:};

    vecCol = [Nodes(vNodeIdx).activeCol];
    vecCol = num2cell(fillmissing(vecCol, 'linear', 'MissingLocations', vecCol<0));
    [Nodes(vNodeIdx).activeCol] = vecCol{:};
end

%% Local functions
function [camPoints, prjPoints, validNodeIdx] = extractNodesPosition(Nodes)

validNodes = [Nodes.activeCol] > 0 & [Nodes.activeRow] > 0;
validNodeIdx = find(validNodes);

camPoints = vertcat(Nodes(validNodes).Centroid);
prjPoints = [[Nodes(validNodes).activeCol]', [Nodes(validNodes).activeRow]'];
end

%% convert skeleton image to node, edge endpoints images
function [imEdges, imNodes, imEndPoints, imNeighbor] = extractStructs(imSkele)
imNeighbor = imfilter(double(imSkele), [1, 1, 1; 1, 0, 1; 1, 1, 1]);
imNeighbor = imNeighbor .* imSkele;

imEdges = imNeighbor == 2;
imEndPoints = imNeighbor == 1;
imNodes = logical(imSkele - imEdges - imEndPoints);
end
end

author: BingyaoHuang

repo: single-shot-pro-cam-calib


text:
function [C,A,G,B] = bekk_parameter_transform(parameters,p,o,q,k,type)
% Parameter transformation for BEKK(p,o,q) multivariate volatility model simulation and estimation
%
% USAGE:
%  [C,A,G,B] = bekk_parameter_transform(PARAMETERS,P,O,Q,K,TYPE)
%
% INPUTS:
%   PARAMETERS - Vector of parameters governing the dynamics.  See BEKK or BEKK_SIMULATE
%   P          - Positive, scalar integer representing the number of symmetric innovations
%   O          - Non-negative, scalar integer representing the number of asymmetric innovations
%   Q          - Non-negative, scalar integer representing the number of conditional covariance lags
%   K          - Number of assets
%   TYPE       - Integer indicating type
%                  1: Scalar
%                  2: Diagonal
%                  3: Full
%
% OUTPUTS:
%   C - K by K covariance model intercept
%   A - K by K by P matrix of symmetric innovation parameters
%   G - K by K by O matrix of asymmetric innovation parameters
%   B - K by K by Q matrix of smoothing parameters
%
% COMMENTS:
%
% See also BEKK, BEKK_SIMULATE, BEKK_LIKELIHOOD

if type==1
    numParams = 1;
elseif type==2
    numParams = k;
else
    numParams = k*k;
end

k2 = k*(k+1)/2;
C = parameters(1:k2);
C = vec2chol(C);
C = C*C';
offset = k2;
[V,D] = eig(C);
D = diag(D);
if (min(D))<(2*eps*max(D))
    D((D/max(D))<eps) = 2*max(D)*eps;
    C = V*diag(D)*V';
    C=(C+C)/2;
end

m = p+o+q;
temp = zeros(k,k,m);
for j=1:m
    tempP = parameters(offset+(1:numParams));
    offset = offset+numParams;
    if type==1
        temp(:,:,j) = tempP*eye(k);
    elseif type==2
        temp(:,:,j) = diag(tempP);
    else
        temp(:,:,j) = reshape(tempP,k,k);
    end
end
A = temp(:,:,1:p);
G = temp(:,:,p+1:p+o);
B = temp(:,:,p+o+1:p+o+q);

author: bashtage

repo: mfe-toolbox


text:
function dX = Motor2D(X,K,isForceField)

% Parameters
m   = 2;%1;  % mass
b   = 10; % viscosity
c1  = 0.15/2;
c2  = 0.05/2;
dt  = 0.005;
tau = 0.05;
A = [ 0 0    1    0;
      0 0    0    1;
      0 0 -b/m    0;
      0 0    0 -b/m];
B = [0   0;
     0   0;
     1/m 0;
     0 1/m];
B1 = [1/tau 0;
     0  1/tau];

x = X(1:6);
z = X(7:8);

dx = x;
x = x(:);
u = -K*x;

w = randn(2,1)*sqrt(dt);
M = [c1*u(1) c2*u(2); -c2*u(1) c1*u(2)];

VF=[-10.1 -11.2; -11.2 11.1];

v = M*w; % control dependent noise
f = z;

dx(1:4) = A * x(1:4)+B*(x(5:6)+f);
dx(5:6) = B1 * (-x(5:6)+u+v./dt);

if isForceField
	dz = -1/0.01*(z-VF*[x(3);x(4)]);
else
	dz=[0;0];
end

dX=[dx;dz];
end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


%==============================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% function [vc,yc,wc,MLhis] = MLLDDMM(ML,varargin)
%
% Multi-Level for Large Deformation Diffeomorphic Metric Mapping (LDDMM)
%
% minimizes J^h(y) = D^h(T(y),R) + S^h(y) for h=coarse:fine
%
% Input:
%   ML         coarse to fine representations of the data, see getMultilevel.m
%   varagin    optional parameters, see default parameter
% Output:
%   vc          numerical optimizer (velocity field)
%   yc          optimal transformation
%   wc          optimal parameters for PIR
%   MLhis       iteration history
%   para        additional info from objective function
%
% for level=minLevel:maxLevel,
%   get data(level)
%   if level==minLevel, get wOpt running PIR; end;
%   use pre-registered Yref=trafo(wOpt,X) for regularization
%   if level==minLevel
%     v0 = trafo('get','w0'); % use this as starting guess
%   else
%     get v0 by prologating vopt from coarse to finer level
%   end;
%   get vopt using GaussNewton using v0 as starting guess
% end%For
%
%==============================================================================
%
% optimizer are controlled via [default]
%    PIRobj      [= @PIRobjFctn;         ]
%    PIRpara     [= optPara('PIR-GN');   ]
%    NPIRobj     [= @NPIRobjFctn;         ]
%    NPIRpara    [= optPara('NPIR-GN')   ]
%
% see also Ex_MLLDDMM.m for an example
%==============================================================================

function [vc,yc,wc,MLhis] = MLLDDMM(ML,varargin)

if nargin == 0
  % run minimal example
  help(mfilename);
  Ex_MLLDDMM;
  yc = 'endOfMinimalExample';
  return;
end;

%------------------------------------------------------------------------------
% initialize some default parameter:
% 1.1  initialize output
% 1.2  determine minLevel and maxLevel from the data ML
% 1.3  flag for parametric registration (PIR)
% 2.1  regularization of of PIR objective
% 2.2  regularization of Hessian approximation, H -> H + beta*speye(size(H))
% 2.3  setup stopping value wStop and initial value w0
% 2.4  setup objective in PIR  and optimization parameters
% 3    setup objective in NPIR and optimization parameters
% 4       setup defaults for plots
% 5    replace defaults by user input
% 6    initializes dimension, discretization, grids etc
% 7    convert optimization parameter sets for PIR and NPIR to lists
%------------------------------------------------------------------------------
% 1.1 initialize output
vc          = [];
yc          = [];
wc          = [];
MLhis       = [];

% 1.2  get minLevel and maxLevel from ML
[ML,minLevel,maxLevel] = getMultilevel(ML);

% 1.3  flag for parametric registration (PIR)
parametric  = 1;     % flag: if true, run parametric pre-registration

% 2.1  regularization of PIR objective
%   J_{PIR} = D(T(Y(wc)),Rc) + S(wc) with S_{PIR} = 0.5*(wc-wRef)'*M*(wc-wRef);
wRef        = [];
M           = [];

% 2.2  regularization of Hessian approximation, H -> H + beta*speye(size(H))
beta        = 0;

% 2.3  setup stopping value wStop and initial value w0
wStop       = [];           % global stopping for PIR
w0          = [];           % starting guess for PIR

% 2.4  PIR: default objective function and optimization parameter via optPara('PIR-GN')
PIRobj      = @PIRobjFctn;
PIRpara     = optPara('PIR-GN');

% 3.   NPIR: default objective function and optimization parameter via optPara('NPIR-GN')
NPIRobj     = @LDDMMobjFctn;
NPIRpara    = optPara('NPIR-GN');

% 4       setup defaults for plots
pause       = 0;        % flag for pauses
plots       = 1;        % flag for plots
plotIter    = 0;        % flag for output of iteration history each level
plotMLiter  = 0;        % flag for output of summarized iteration history

omegaV      = [];       % domain for velocity field
N           = 5;        % number of forward euler steps
mV          = @(m) m;   % resolution of velocity field
nt          = regularizer('get','nt');        % number of time discretizations for velocities
vStop       = [];           % global stopping for NPIR
vRef        = [];           % regularization: S(vc-vRef)

% 5    replace defaults by user input
for k=1:2:length(varargin) % overwrites default parameter
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

% 6    initializes dimension, discretization, grids etc
dimstr      = @(m) sprintf('[%s]',sprintf(' %d',m));
omega       = ML{end}.omega;  % spatial domain
grid        = regularizer('get','grid');
switch grid
  case 'cell-centered', getGrid = @(m) getCellCenteredGrid(omega,m);
  case 'staggered',     getGrid = @(m) getStaggeredGrid(omega,m);
  case 'nodal',         getGrid = @(m) getNodalGrid(omega,m);
  otherwise, error('nyi');
end;

% 7    convert optimization parameter sets for PIR and NPIR to lists
PO = FAIRcell2struct(PIRpara);
NO = FAIRcell2struct(NPIRpara);
%------------------------------------------------------------------------------

fprintf('%s: MultiLevel Large Deformation Diffeomorphic Metric Mapping\n',mfilename)
fprintf('>> %-20s : %s\n','omega',dimstr(omega));
fprintf('>> %-20s : %s\n','IMAGE MODEL',imgModel);
fprintf('>> %-20s : %s\n','DISTANCE',distance);
fprintf('>> %-20s : %s\n','TRAFO',trafo);
fprintf('>> %-20s : %s, alpha=%s\n','REGULARIZER',...
  regularizer,num2str(regularizer('get','alpha')));
fprintf('>> %-20s : %d\n','PRE-REGISTRATION',parametric);
fprintf('>> %-20s : objective=<%s>, method = <%s>\n',...
  'PIR',func2str(PIRobj),func2str(PIRpara.scheme));
fprintf('>> %-20s : objective=<%s>, method = <%s>\n',...
  'NPIR',func2str(NPIRobj),func2str(NPIRpara.scheme));
fprintf('>> %-20s : minLevel=%d, maxLevel=%d\n','LEVEL',minLevel,maxLevel);
FAIRmessage('#')

tic;
%--------------------------------------------------------------------------
for level=minLevel:maxLevel

  FAIRmessage(sprintf('%s: level %d from %d to %d, %s',...
    mfilename,level,minLevel,maxLevel,dimstr(ML{level}.m)));

  % save(old grid), update(m,grid,data coefficients)
  m     = ML{level}.m;
  xc    = getGrid(m);
  [T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega);
  Rc    = imgModel(R,omega,center(xc,m));
  if strcmp(func2str(NPIRobj),'MPLDDMMobjFctn')
      T = imgModel(T,omega,center(xc,m));
  end

  if level == minLevel    % on coarse level

    if parametric          % perform parametric pre-registration

      % initialize plots for PIR
      FAIRplots('set','mode','PIR','fig',minLevel-1,'plots',plots);
      FAIRplots('init',struct('Tc',T,'Rc',R,'omega',omega,'m',m));

      % ----- call Parametric Image Registration ----------------------------
      xC = getGrid(m);
      PIRfctn = @(wc) PIRobj(T,Rc,omega,m,beta,M,wRef,xC,wc);
      PIRfctn([]); % report status of objective function
      if isempty(w0),    w0    = trafo('w0'); end;
      if isempty(wStop), wStop = w0;          end;

      [wc,hisPIR] = PIRpara.scheme(PIRfctn,w0,'yStop',wStop,PO{:});
      fprintf('wc = \n');
      disp(wc');
      doPause(pause);
      % ----- end PIR --------------------------------------

      if plotIter
        hisPIR.str{1} = 'iteration history for PIR';
        plotIterationHistory(hisPIR,'J',1:4,'fh',100+minLevel-1);
      end;
    else                    % no pre-registration
      wc = [];
    end;

  end;

  % compute xc = trafo(wc,xc)
  if ~isempty(wc) % use yRef(x) = x
    xc = trafo(wc,xc);
  end

  % compute starting guess y0 for this level and stopping value yStop
  if level == minLevel
      if isempty(vRef)
          vRef = getVelocityStartingGuess(omegaV,mV(m),nt);
      end
      v0 = vRef;  % the best known so far
  else
    % prolongate vc (coarse) vRef (current)
    vcOld   = reshape(vc,[],nt+1);
    vRefOld = reshape(vRef,[],nt+1);
    v0 = []; vRef = [];
    for k=1:nt+1
        v0   = [v0 mfPu(vcOld(:,k),omega,mV(m/2))];
        vRef = [vRef mfPu(vRefOld(:,k),omega,mV(m/2))];
    end
    v0 = v0(:); vRef = vRef(:);
 end;
    vStop = 0*getVelocityStartingGuess(omegaV,mV(m),nt);

  % ----- call NPIR -------------------------------------
  % initialize plots for NPIR
  FAIRplots('reset','mode','NPIR','fig',level,'plots',plots);
  FAIRplots('init',struct('Tc',T,'Rc',R,'omega',omega,'m',m));

  NPIRfctn = @(vc) NPIRobj(T,Rc,omega,m,vRef,xc,omegaV,mV(m),N,vc);
  if level == minLevel   % report status of objective function
    NPIRfctn([]);
  end; %

  [vc,his] = NPIRpara.scheme(NPIRfctn,v0,'yStop',vStop,NO{:});
  % ----- end NPIR --------------------------------------

  if plotIter
    his.str{1} = sprintf('iteration history for LDDMM, level=%d',level);
    plotIterationHistory(his,'J',1:4,'fh',100+level);
  end;

  % update iteration history
  if level == minLevel
    MLhis.str = his.str;
    MLhis.his = his.his;
  else
    MLhis.his = [MLhis.his;his.his];
  end;
  doPause(pause);

%--------------------------------------------------------------------------
end;%For level
%--------------------------------------------------------------------------

MLhis.time = toc;

if plotMLiter
  plotMLIterationHistory(MLhis,'fh',[]);
end;

id0 = find(MLhis.his(:,1)==-1,1,'last');
MLhis.reduction = MLhis.his(end,2)/MLhis.his(id0,2);
J = find(MLhis.his(:,1)==-1);
MLhis.iter(minLevel:maxLevel) = MLhis.his([J(2:end)-1;size(MLhis.his,1)],1)';

FAIRmessage([mfilename,' : done !'],'#');

%--------------------------------------------------------------------------

function doPause(p)
if strcmp(p,'on')
  FAIRpause;
elseif p>0
  FAIRpause(p);
end;

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


text:
%   AUTHORSHIP
%   Primary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Math Lead & Secondary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%
classdef CellBasics
    methods(Static)
        function tabLines=ToTabLines(data, labels, removeHtml)
            if nargin<3
                removeHtml=true;
            end
            if removeHtml
                J=edu.stanford.facs.swing.Basics;
            end
            tabLines=java.util.ArrayList;
            [rows,cols]=size(data);
            if nargin>1
                data=[labels; data];
            end
            for row=1:rows
                line='';
                for col=1:cols
                    value=data{row, col};
                    if removeHtml && ischar(value)
                        value=J.RemoveXml(value);
                    end
                    if isnumeric(value)
                        value=num2str(value);
                    end
                    line=sprintf('%s%s\t', line, value);
                end
                tabLines.add(line);
            end
        end

        function groups=GetRowGroups(sortedRows)
            groups={};
            [R,C]=size(sortedRows);
            if R>0
                lastStartingR=1;
                for r=2:R
                    if sortedRows(r,1)>sortedRows(r-1,1)+1
                        groups{end+1}=sortedRows(lastStartingR:r-1,:);
                        lastStartingR=r;
                    end
                end
            end
            groups{end+1}=sortedRows(lastStartingR:r,:);
        end

        function rowIdxs=Find(c, colIdxs, fnc, first, rowStartIdx)
            if nargin<5
                rowStartIdx=1;
                if nargin<4
                    first=true;
                end
            end
            N=size(c,1);
            C=size(c,2);
            colIdxsN=length(colIdxs);
            rowIdxs=[];
            for i=rowStartIdx:N
                for j=1:colIdxsN
                    colIdx=colIdxs(j);
                    if colIdx>=1 && colIdx<=C
                        o=c{i,colIdx};
                        if fnc(o, i, colIdx)
                            if colIdxsN==1
                                rowIdxs(end+1)=i;
                            else
                                rowIdxs(end+1,:)=[i colIdx];
                            end
                            if first
                                return;
                            end
                        end
                    end
                end
            end
        end

        function yes=equals(o, arg)
            yes=o==arg;
        end

        function [found, idx]=Contains(cell, object)
            found=false;
            N_=length(cell);
            for idx=1:N_
                if isequal(object, cell{idx})
                    found=true;
                    return;
                end
            end
            idx=0;
        end

        function c=Java(collectionOrArray)
            try
                N=collectionOrArray.size;
                c=cell(1,N);
                it=collectionOrArray.iterator;
                for i=1:N
                    c{i}=it.next;
                end
            catch
                N=length(collectionOrArray);
                c=cell(1,N);
                for i=1:N
                    c{i}=collectionOrArray(i);
                end

            end
        end

        function u=UniqueNumbers(cellNums)
            if isnumeric(cellNums)
                u=unique(cellNums);
            else
                u=[];
                N=length(cellNums);
                for i=1:N
                    u=unique([u cellNums{i}']);
                end
            end
        end
    end
end