text:
function [pvec, pstruct] = tapas_hgf_categorical_transp(r, ptrans)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Number of states whose contingencies have to be learned
no = r.c_prc.n_outcomes;

pvec    = NaN(1,length(ptrans));
pstruct = struct;

pvec(1:no)       = ptrans(1:no);                           % mu2_0
pstruct.mu2_0    = pvec(1:no);
pvec(no+1:2*no)  = exp(ptrans(no+1:2*no));                 % sa2_0
pstruct.sa2_0    = pvec(no+1:2*no);
pvec(2*no+1)     = ptrans(2*no+1);                         % mu3_0
pstruct.mu3_0    = pvec(2*no+1);
pvec(2*no+2)     = exp(ptrans(2*no+2));                    % sa3_0
pstruct.sa3_0    = pvec(2*no+2);
pvec(2*no+3)     = tapas_sgm(ptrans(2*no+3),r.c_prc.kaub); % ka
pstruct.ka       = pvec(2*no+3);
pvec(2*no+4)     = ptrans(2*no+4);                         % om
pstruct.om       = pvec(2*no+4);
pvec(2*no+5)     = tapas_sgm(ptrans(2*no+5),r.c_prc.thub); % th
pstruct.th       = pvec(2*no+5);

return;

author: translationalneuromodeling

repo: tapas


text:
function inside = polygon_contains_point_2d ( n, v, p )

%*****************************************************************************80
%
%% POLYGON_CONTAINS_POINT_2D finds if a point is inside a simple polygon in 2D.
%
%  Discussion:
%
%    A simple polygon is one whose boundary never crosses itself.
%    The polygon does not need to be convex.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    M Shimrat,
%    Position of Point Relative to Polygon,
%    ACM Algorithm 112,
%    Communications of the ACM,
%    Volume 5, Number 8, page 434, August 1962.
%
%  Parameters:
%
%    Input, integer N, the number of nodes or vertices in the polygon.
%    N must be at least 3.
%
%    Input, real V(2,N), the coordinates of the vertices of the polygon.
%
%    Input, real P(2), the coordinates of the point to be tested.
%
%    Output, logical INSIDE, is TRUE if the point is inside the polygon.
%
  inside = 1;

  for i = 1 : n

    x1 = v(1,i);
    y1 = v(2,i);

    if ( i < n )
      x2 = v(1,i+1);
      y2 = v(2,i+1);
    else
      x2 = v(1,1);
      y2 = v(2,1);
    end

    if ( ( y1 < p(2) & p(2) <= y2 ) | ( p(2) <= y1 & y2 < p(2) ) )
      if ( ( p(1) - x1 ) - ( p(2) - y1 ) * ( x2 - x1 ) / ( y2 - y1 ) < 0.0 )
        inside = ~inside;
      end
    end

  end

  inside = ~inside;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Wrapper for vl_nnvisibilitymask block
% inputs{1} :   X     : 1 x 3 x n x b (Rotated vertices)
% obj       :   faces : nfaces x 3
% outputs{1}:   y     : 112 x 112 x 3 x b

classdef visibilitymask < dagnn.Layer

    properties
        faces = [];
    end

    methods
        function outputs = forward(obj, inputs, ~)

            useGPU = isa(inputs{1}, 'gpuArray');
            if useGPU
                outputs{1} = gpuArray( vl_nnvisibilitymask(gather(inputs{1}), obj.faces) );
            else
                outputs{1} = vl_nnvisibilitymask(inputs{1}, obj.faces);
            end

        end

        function [derInputs, derParams] = backward(obj, inputs, ~, derOutputs)

            useGPU = isa(inputs{1}, 'gpuArray');
            if useGPU
                derInputs{1} = gpuArray( vl_nnvisibilitymask(gather(inputs{1}), obj.faces, gather(derOutputs{1})) );
            else
                derInputs{1} = vl_nnvisibilitymask(inputs{1}, obj.faces, derOutputs{1});
            end

            derParams = {};
        end

        function outputSizes = getOutputSizes(obj, inputSizes)
            outputSizes = {112 112 3 inputSizes{1}(4)} ; % Note: We may want to set this dynamically as well
        end

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

author: anilbas


text:
function [mrk, ev]= mrk_removeVoidClasses(mrk)
%mrk= mrk_selectEvents(mrk)
%
% requires a field 'y' containing the class labels in mrk.

misc_checkType(mrk, 'STRUCT(time y)');

nonvoidClasses= find(any(mrk.y,2));
if length(nonvoidClasses)<size(mrk.y,1),
  msg= sprintf('void classes removed, %d classes remaining', ...
                  length(nonvoidClasses));
%   util_warning(msg, 'mrk', mfilename);
  warning(msg, 'mrk');
  mrk.y= mrk.y(nonvoidClasses,:);
  if isfield(mrk, 'className'),
    mrk.className= {mrk.className{nonvoidClasses}};
  end
end

author: bbci


text:
function [CAN_POOL] = generateObjHyps( aid, numhyp, rect_score )
%GENERATESEEDS Popup object hypothesis from 2D to 3D
%   Detailed explanation goes here
global config

%% load data
rectangle = [];
rectCuboid2 = [];

load([config.bufname config.roomHypoFile]);
load([config.bufname config.objHypoFileUpdate]);
% load([config.folderName config.IMGLIST(aid).name '/' config.roomHypoFile]);
% load([config.folderName config.IMGLIST(aid).name '/' config.objHypoFileUpdate]);
% load([config.folderName config.IMGLIST(aid).name '/' config.rectangleScore]);

load(config.classifierName);
if ~exist('ClassNames','var')
    ClassNames = B.ClassNames;
end
% load('./object_classifier/objectclassifier.mat');
% load('./living_room/objectClassifier_livingroom.mat');

%% update
% ALL_HYPS = cell(numhyp, numsample);

all_opScr = [hypothesis.opScr];
[S,I] = sort(all_opScr, 'descend');
hypothesis = hypothesis(I(1:config.sampleroomnum));
hypothesis_ori = hypothesis_ori(I(1:config.sampleroomnum));

CAN_POOL = cell(min(numhyp, length(hypothesis)),1);

for hid = 1:min(numhyp, length(hypothesis))
% while hid<=numhyp
    room_hyp = hypothesis(hid);
    camera_height = 160;
    [ room_obj ] = room3DFromHyp( room_hyp.hCorner, -camera_height );

    room_obj.omapScr = hypothesis(hid).omapScr;
    room_obj.gcScr = hypothesis(hid).gcScr;
    room_obj.mgScr = hypothesis(hid).mgScr;
    room_obj.opScr = hypothesis(hid).opScr;

    out_sign = sign(room_obj.out_points_w);
    if ~all(sum(out_sign,1)==0)
        continue;
    end

    %% get all hypotheses and project to 3D
    rectHyps = repmat(struct('objects',[],'scores',[],'valid',[],'bottomscore',[]), 14, 1);

    fprintf('Generaing cuboid: %d/%d, part 1\n', aid, hid);
    parfor rid = 1:5 %parfor
        objects = repmat(struct('out_points_w',[],'name',[],'type',[],'points',[],'x_w',[]),0,1);
        bottomscore = zeros(0,1);
        num = size(rectangle(rid).xyzBox,1);
        validIDs = find(rect_score{rid}>0.6);
        for bid = validIDs % 1:num
    %         fprintf('%d/%d\n', bid, num);
            if rid<=4
                newObj = rectangle2hypothesis( rectangle(rid).xyzBox(bid,:), rid, room_obj.out_points_w );
            elseif rid==5
                newObj = floorRectangle2hypothesis( rectangle(rid).xyzBox(bid,:), rid, room_obj.out_points_w );
            end
            objects(end+1:end+length(newObj)) = newObj;
            bottomscore(end+1:end+length(newObj),1) = rectangle(rid).score(bid);
    %         figure(1); clf; viewSingleRectangle(rotImg, rectangle, rid, bid);
    %         figure(2); clf; visualObject3D( [newObj;room_obj], 4, [1 1 1] );
        end
        rectHyps(rid).objects = objects;
        rectHyps(rid).bottomscore = bottomscore;
    end

    fprintf('Generaing cuboid: %d/%d, part 2\n', aid, hid);
    objects = repmat(struct('out_points_w',[],'name',[],'type',[],'points',[],'x_w',[]),0,1);
    bottomscore = zeros(0,1);
    for cid = 1:regionCuboid.count
    %     fprintf('%d/%d\n', cid, 100);
        views = regionCuboid.views(cid,:);
        if any(views==6)
            continue;
        end
        rects = regionCuboid.xyzBox(cid,:);
        rects = [rects(1:12); rects(13:24); rects(25:36)];
    %     newObj = obj3DFromHyp( rects, views, room_obj.out_points_w );
        newObj = rectangleCombo3Hypothesis( rects, views, room_obj.out_points_w );
        objects(end+1:end+length(newObj)) = newObj;
        bottomscore(end+1:end+length(newObj),1) = regionCuboid.score(cid);
    %     figure(1); clf; viewRegionHypothesisSubset( rotImg, regionCuboid, cid );
    %     figure(2); clf; visualObject3D( [newObj;room_obj], 4 );
    end
    rectHyps(6).objects = objects;
    rectHyps(6).bottomscore = bottomscore;

    fprintf('Generaing cuboid: %d/%d, part 3\n', aid, hid);
    parfor rid = 1:8 %parfor
    %     fprintf('%d/%d\n', rid, length(rectCuboid2));
        views = rectCuboid2(rid).views;
        views(views==0) = [];
        score1 = rectangle(views(1)).score;
        score2 = rectangle(views(2)).score;
        if any(views==5)
            num = length(rectCuboid2(rid).score);
            objects = repmat(struct('out_points_w',[],'name',[],'type',[],'points',[],'x_w',[]),0,1);
            bottomscore = zeros(0,1);
            for cid = 1:num
                combo = rectCuboid2(rid).combo(cid,:);
                rects = [rectangle(views(1)).xyzBox(combo(1),1:12); rectangle(views(2)).xyzBox(combo(2),1:12)];
                newObj = floorRectangleCombo2Hypothesis( rects, views, room_obj.out_points_w );

    %             figure(1); clf; viewRectHypothesisSubset(rotImg, rectCuboid2, rectangle, rid, cid);
    %             figure(2); clf; visualObject3D( [newObj;room_obj], 4 );
                objects(end+1:end+length(newObj)) = newObj;
                bottomscore(end+1:end+length(newObj),1) = max(score1(cid), score2(cid)); %rectCuboid2(rid).score(cid);
            end
        else
            num = length(rectCuboid2(rid).score);
            objects = repmat(struct('out_points_w',[],'name',[],'type',[],'points',[],'x_w',[]),0,1);
            bottomscore = zeros(0,1);
            for cid = 1:num
                combo = rectCuboid2(rid).combo(cid,:);
                rects = [rectangle(views(1)).xyzBox(combo(1),1:12); rectangle(views(2)).xyzBox(combo(2),1:12)];
                newObj = rectangleCombo2Hypothesis( rects, views, room_obj.out_points_w );

    %             figure(1); clf; viewRectHypothesisSubset(rotImg, rectCuboid2, rectangle, rid, cid);
    %             figure(2); clf; visualObject3D( [newObj;room_obj], 4 );
                objects(end+1:end+length(newObj)) = newObj;
                bottomscore(end+1:end+length(newObj),1) = max(score1(cid), score2(cid)); %rectCuboid2(rid).score(cid);
            end
        end

        rectHyps(rid+6).objects = objects;
        rectHyps(rid+6).bottomscore = bottomscore;
    end

    %% classify: semantic classifier, should move to DDsampling
    for rid = 1:14
        [ scores, valid ] = objectClassifier(room_obj, rectHyps(rid).objects, B, ClassNames);
        rectHyps(rid).scores = scores;
        rectHyps(rid).valid = valid;
%         num = length(rectHyps(rid).objects);
%         rectHyps(rid).scores = zeros(num,28);
%         rectHyps(rid).valid = false(num,1);
    end
    all_objects = repmat(struct('out_points_w',[],'name',[],'type',[],'points',[],'x_w',[]),0,1);
    scores = zeros(0,28);
    valid = false(0,1);
    obj_xyz = zeros(0,6);
    bottomscore = zeros(0,1);
    for rid = 1:14
        num = length(rectHyps(rid).objects);
        all_objects(end+1:end+num,:) = rectHyps(rid).objects;
        scores(end+1:end+num,:) = rectHyps(rid).scores ;
        valid(end+1:end+num,1) = rectHyps(rid).valid;
        bottomscore(end+1:end+num,1) = rectHyps(rid).bottomscore;
        xyz = zeros(num,6);
        for i = 1:num
            xyz(i,:) = [min(rectHyps(rid).objects(i).out_points_w, [], 1) ...
                max(rectHyps(rid).objects(i).out_points_w, [], 1)];
        end
        obj_xyz(end+1:end+num,:) = xyz;
    end

    scores(~valid,:) = 0;
    all_hyps.objects = addBoundingBox(all_objects(valid),'both'); %all_objects(valid);
    all_hyps.scores = scores(valid,:);
    all_hyps.valid = valid(valid);
    all_hyps.obj_xyz = obj_xyz(valid,:);
    all_hyps.bottomscore = bottomscore(valid);
    fprintf('%d:%d, %d-->%d\n', aid, hid, length(valid), sum(valid));

    %% prune weak cuboid hypothesis
%     t = 0;
%     sel_hyps = [];
%     for oid = 1:28
%         valid = bottomscore>0 & scores(:,oid)>t;
%         sel_hyps(oid).h_objects = all_hyps.objects(valid);
%         sel_hyps(oid).h_scores = all_hyps.scores(valid,:);
%         sel_hyps(oid).h_obj_xyz = all_hyps.obj_xyz(valid,:);
%         sel_hyps(oid).h_bottomscore = all_hyps.bottomscore(valid,:);
%         sel_hyps(oid).h_objID = find(valid);
%     end
% %     t = 0.60;
% %     sel_hyps = [];
% %     for oid = 1:28
% %         [S,I] = sort(all_hyps.scores(:,oid),'descend');
% %         obj_xyz = all_hyps.obj_xyz(I(S>t),:);
% %         [ outxyz, valid ] = nms3D( obj_xyz, 0.5 );
% %         sel_hyps(oid).objects = all_hyps.objects(I(valid));
% %         sel_hyps(oid).scores = all_hyps.scores(I(valid),:);
% %         sel_hyps(oid).obj_xyz = all_hyps.obj_xyz(I(valid),:);
% %         sel_hyps(oid).bottomscore = all_hyps.bottomscore(I(valid),:);
% %         sel_hyps(oid).objID = I(valid);
% %
% %         v = sel_hyps(oid).bottomscore>0;
% %
% %         sel_hyps(oid).h_objects = sel_hyps(oid).objects(v);
% %         sel_hyps(oid).h_objID = sel_hyps(oid).objID(v);
% %         sel_hyps(oid).h_scores = sel_hyps(oid).scores(v,:);
% %         sel_hyps(oid).h_obj_xyz = sel_hyps(oid).obj_xyz(v,:);
% %         sel_hyps(oid).h_bottomscore = sel_hyps(oid).bottomscore(v,:);
% %     end
%
%     % pack gnd
%     [ALL_HYPS(hid,:), CAN_POOL{hid}.sel_hyps, CAN_POOL{hid}.room] = getInitialWholeLayoutB( sel_hyps, room_obj, numsample );
    CAN_POOL{hid}.sel_hyps = all_hyps;
    CAN_POOL{hid}.room = room_obj;
end

% ALL_HYPS = ALL_HYPS(:);

% %% output candidate pool
% for tid = 1:12
%
% end

end

author: zouchuhang

repo: LayoutNet


text:
function [node,elem,phiValue,oldNodeIdx,oldNode,oldElemIdx,oldElem] = interfaceadaptivemesh(node,elem,phi)
% INTERFACEMESH bisect the initial mesh to get a mesh fiting the interface
%
%   [node,elem] = INTERFACEADAPTIVEMESH(node,elem,phi,a), where
%    - node and elem represent an initial mesh enclosing the interface;
%    - phi is a level set function defining the interface as follows
%       {p | phi(p)=0} represents the interface
%       {p | phi(p) < 0} represents the interior domain enclosed by the interface
%       {p | phi(p) > 0} rpepresents the exterior domain of the interface.
%
%  Definition
%   Interface element: the element whose vertices are not in the same side
%   of the interface.
%
%   Interface node: the vertices of all the interface elements. Split the
%   set of interface nodes into two parts: interior interface nodes (phi(p)<=0) and
%   exterior interface nodes (phi(p)>=0);
%
%  Example
%   clear all; close all;
%   node = [-1 -1; 1 -1; 1 1; -1 1];
%   elem = [2 3 1; 4 1 3];
%   [node,elem] = interfacemesh(node,elem,@ellipse);
%
%  See also bisect, circle, ellipse.
%
%  Reference: H. Wei, L. Chen, Y. Huang and B. Zheng. Adaptive Mesh
%  Refinement and Superconvergence for Two Dimensional Interface Problems.
%  SIAM Journal on Scientific Computing, 2014.
%
%  Add by Huayi Wei. Based on discussion with Long Chen.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

oldNodeIdx = [];
oldNode = [];
oldElemIdx = [];
oldElem = [];

%% Step 1: Refine mesh such that the set of interface elements is non-empty
N0 = size(node,1);
phiValue = phi(node);

if all(phiValue<0)
    error('The initial mesh is in the interior of the domain the interface surrounds, please use a mesh which can cover the interface!')
end

% uniform bisect the mesh, to make sure there exist at least one mesh node
% which is inside of the interface, namely, phi(p)<0.
while all(phiValue>0)
    [node,elem] = uniformbisect(node,elem);
    N1 = size(node,1);
    phiValue(N0+1:N1) = phi(node(N0+1:N1,:)); % append the phi value
    N0 = N1;
end

display('Finish Step 1: uniform bisect all elements');

%% Step 2: estimate the discrete curvature, then bisect elements with big curvature
% find interface elements, which are elements crossed by the interface and
% restrict the computation to the interface elements
[interfaceElemIdx, ~, isInterfaceElem] = markinterface(elem,phiValue);
interfaceElem = elem(interfaceElemIdx,:);
isInterfaceNode = false(size(node,1),1);
isInterfaceNode(interfaceElem(:))=true;
interfaceNT = size(interfaceElem,1);
T = auxstructure(interfaceElem);
edge = double(T.edge);
edge2elem = T.edge2elem;
clear T;

% case: Fig 2(a) in p6
tmpAngle = ones(interfaceNT,1)*[90,45,45];
angle = accumarray(interfaceElem(:),tmpAngle(:),[size(node,1) 1]);
is360 = (angle == 360);

% case: Fig 2(b) in p6
isInterfaceBEdge = (edge2elem(:,1) == edge2elem(:,2));
interfaceBEdge = edge(isInterfaceBEdge,:);
v2vI = sparse(interfaceBEdge,interfaceBEdge(:,[2,1]),1,size(node,1),size(node,1));
isLinkNode = (v2vI*isInterfaceNode >= 3);

% case: Fig 2(c) in p6
isInteriorNode = phiValue < 0;
turnAngle = zeros(size(node,1),1);
isIn = isInterfaceNode & isInteriorNode & ~is360 & ~isLinkNode;
turnAngle(isIn) = 180 - (360 - angle(isIn));
isOut = isInterfaceNode & ~isInteriorNode & ~is360 & ~isLinkNode;
turnAngle(isOut ) = (360 - angle(isOut)) -180;
v2v = sparse(edge,edge(:,[2,1]),1,size(node,1),size(node,1));
turnAngleTotal = v2v*turnAngle+turnAngle;

% case: size of interface elements should be small enough
ve = node(elem(:,3),:) - node(elem(:,2),:);
L = sqrt(sum(ve.^2,2));
mL = max(L);
isBigSizeElem = ((L > mL*mL) & isInterfaceElem);
% mark nodes for refinement
isBigCurvNode = (is360 | isLinkNode |(abs(turnAngleTotal)> 170));

if any(isBigCurvNode) || any(isBigSizeElem)

    bigCurvElem = interfaceElemIdx(max(isBigCurvNode(interfaceElem),[],2));
    bigSizeElem = find(isBigSizeElem);

    count = 0;
    while  count < 100
        % bisect marked elements
        [node,elem] = bisect(node,elem,[bigCurvElem; bigSizeElem]);

        N1 = size(node,1);
        phiValue(N0+1:N1) = phi(node(N0+1:N1,:));
        N0 = N1;

        % find interface elements and nodes
        interfaceElemIdx = markinterface(elem,phiValue);
        interfaceElem = elem(interfaceElemIdx,:);
        isInterfaceNode = false(size(node,1),1);
        isInterfaceNode(interfaceElem(:))=true;
        T = auxstructure(interfaceElem);
        edge = double(T.edge);
        edge2elem = T.edge2elem;
        clear T;

        % find triangles containing curves with high curvature
        % case: Fig 2(a) in p6
        interfaceNT = size(interfaceElem,1);
        tmpAngle = ones(interfaceNT,1)*[90,45,45];
        angle = accumarray(interfaceElem(:),tmpAngle(:),[size(node,1) 1]);
        is360 = angle == 360;

        % case: Fig 2(b) in p6
        isInterfaceBEdge = edge2elem(:,1) == edge2elem(:,2);
        interfaceBEdge = edge(isInterfaceBEdge,:);
        v2vI = sparse(interfaceBEdge,interfaceBEdge(:,[2,1]),1,size(node,1),size(node,1));
        isLinkNode = (v2vI*isInterfaceNode >= 3);

        % case: Fig 2(c) in p6
        isInteriorNode = (phiValue < 0);
        turnAngle = zeros(size(node,1),1);
        isIn = isInterfaceNode & isInteriorNode & ~is360 & ~isLinkNode;
        turnAngle(isIn) = 180 - (360 - angle(isIn));
        isOut = isInterfaceNode & ~isInteriorNode & ~is360 & ~isLinkNode;
        turnAngle(isOut ) = (360 - angle(isOut)) -180;
        v2v = sparse(edge,edge(:,[2,1]),1,size(node,1),size(node,1));
        turnAngleTotal = v2v*turnAngle+turnAngle;

        % control size of interface elements
        ve = node(interfaceElem(:,3),:) - node(interfaceElem(:,2),:);
        L = sqrt(sum(ve.^2,2));
        bigSizeElem = interfaceElemIdx(L > mL^2);

        isBigCurvNode = is360 | isLinkNode |(abs(turnAngleTotal)> 170);
        if ~any(isBigCurvNode) && ~any(isBigSizeElem)
            break;
        end

        bigCurvElem = interfaceElemIdx(max(isBigCurvNode(interfaceElem),[],2));
    end
end

display('Finish Step 2: bisect interface with big curvature!');

%% Step 3: Update the type of interface triangles
% Two types of triangles in the bisection grids
%  Type A: two edges parallel to coordinates; see Fig 2.6
%  Type B: the longest edge parallel to x or y coordinate; see Fig 2.7

ve = node(elem(:,3),:) - node(elem(:,2),:);
elemType = (abs(ve(:,1)) > eps & abs(ve(:,2)) > eps);
[interfaceElemIdx, ~,isInterfaceElem] = markinterface(elem,phiValue);
T = auxstructure(elem);
neighbor = T.neighbor;
clear T;
isMark = false(size(elem,1),1);
isMark(interfaceElemIdx) = true;
isMark(neighbor(interfaceElemIdx,1)) = true;

t = cumsum(ones(size(elem,1),1));
tt = neighbor(:,1);
isTransitionTypeBElem = isInterfaceElem & isInterfaceElem(tt) & neighbor(tt,1) ~= t & ~elemType;

%% Deal with the transition Type B elements, see Fig 2.7 (c) and (d) and conditions (C1) and (C2)
%
%  case 1: a Type B interface element whose second neighbor is a interface
%          and  Type A element.
idx2 = find(isTransitionTypeBElem & isInterfaceElem(neighbor(:,2)) ...
            & elemType(neighbor(:,2)));
if ~isempty(idx2)
    p1 = node(elem(idx2,1),:);
    p2 = node(elem(idx2,2),:);
    p3 = node(elem(idx2,3),:);
    p4 = node(elem(neighbor(idx2,2),1),:);
    m14 = (p1 + p4)/2;
    m34 = (p3 + p4)/2;
    m23 = (p2 + p3)/2;
    m1 = (m23 + p3)/2;
    m2 = (m23 + p2)/2;
    isNotTB = ((phiValue(elem(idx2,3)) .* phi(m1)) > 0 & ...
               (phi(m14) .* phi(m34) < 0)) | phi(m1).*phi(m2) < 0;
    isTransitionTypeBElem(idx2(isNotTB)) = false;
end

%  case 2: a Type B interface element whose third neighbor is a interface
%          and  Type A element.
idx3 = find(isTransitionTypeBElem & isInterfaceElem(neighbor(:,3)) ...
            & elemType(neighbor(:,3)));
if ~isempty(idx3)
    p1 = node(elem(idx3,1),:);
    p2 = node(elem(idx3,2),:);
    p3 = node(elem(idx3,3),:);
    p4 = node(elem(neighbor(idx3,3),1),:);
    m14 = (p1 + p4)/2;
    m24 = (p2 + p4)/2;
    m23 = (p2 + p3)/2;
    m1 = (m23 + p3)/2;
    m2 = (m23 + p2)/2;
    isNotTB = ((phiValue(elem(idx3,2)).*phi(m2)) > 0 & ...
               (phi(m14).*phi(m24) < 0)) | phi(m1).*phi(m2) < 0;
    isTransitionTypeBElem(idx3(isNotTB)) = false;
end

%  case 3: a Type B interface element whose third neighbor is a interface
%          and  Type B element.
idx4 = find(isTransitionTypeBElem & isInterfaceElem(neighbor(:,2)) ...
          & ~elemType(neighbor(:,2)));
if ~isempty(idx4)
    p2 = node(elem(idx4,2),:);
    p3 = node(elem(idx4,3),:);
    m23 = (p2 + p3)/2;
    m1 = (m23 + p3)/2;
    m2 = (m23 + p2)/2;
    isNotTB =  phi(m1).*phi(m2) < 0;
    isTransitionTypeBElem(idx4(isNotTB)) = false;
end

%  case 4: a Type B interface element whose third neighbor is a interface
%          and  Type B element.
idx5 = find(isTransitionTypeBElem & isInterfaceElem(neighbor(:,3)) ...
          & ~elemType(neighbor(:,3)));
if ~isempty(idx5)
    p2 = node(elem(idx5,2),:);
    p3 = node(elem(idx5,3),:);
    m23 = (p2 + p3)/2;
    m1 = (m23 + p3)/2;
    m2 = (m23 + p2)/2;
    isNotTB =  phi(m1).*phi(m2) < 0;
    isTransitionTypeBElem(idx5(isNotTB)) = false;
end
typeB = (isMark & ~elemType & ~isTransitionTypeBElem);

% bisect some transition type B elements to type A
while any(typeB)

    [node,elem,~, ~,tree] = bisect(node,elem,typeB);
    N1 = size(node,1);
    phiValue(N0+1:N1) = phi(node(N0+1:N1,:));
    N0 = N1;

    elemType = updatetype(elemType,tree(:,[2,3]));

    [interfaceElemIdx, ~,isInterfaceElem] = markinterface(elem,phiValue);

    T = auxstructure(elem);
    neighbor = T.neighbor;
    clear T;
    isMark = false(size(elem,1),1);
    isMark(interfaceElemIdx) = true;
    isMark(neighbor(interfaceElemIdx,1)) = true;

    t = cumsum(ones(size(elem,1),1));
    tt = neighbor(:,1);
    isTransitionTypeBElem = isInterfaceElem & isInterfaceElem(tt) & ...
                            neighbor(tt,1) ~= t & ~elemType;

    %  case 1: a Type B interface element whose second neighbor is a interface
    %          and  Type A element.
    idx2 = find(isTransitionTypeBElem & isInterfaceElem(neighbor(:,2)) & ...
                elemType(neighbor(:,2)));
    if ~isempty(idx2)
        p1 = node(elem(idx2,1),:);
        p2 = node(elem(idx2,2),:);
        p3 = node(elem(idx2,3),:);
        p4 = node(elem(neighbor(idx2,2),1),:);
        m14 = (p1 + p4)/2;
        m34 = (p3 + p4)/2;
        m23 = (p2 + p3)/2;
        m1 = (m23 + p3)/2;
        m2 = (m23 + p2)/2;
        isNotTB = ((phiValue(elem(idx2,3)) .* phi(m1)) > 0 & ...
                   (phi(m14) .* phi(m34) < 0)) | phi(m1).*phi(m2) < 0;
        isTransitionTypeBElem(idx2(isNotTB)) = false;
    end
    %  case 2: a Type B interface element whose third neighbor is a interface
    %          and  Type A element.
    idx3 = find(isTransitionTypeBElem & isInterfaceElem(neighbor(:,3)) & ...
                elemType(neighbor(:,3)));
    if ~isempty(idx3)
        p1 = node(elem(idx3,1),:);
        p2 = node(elem(idx3,2),:);
        p3 = node(elem(idx3,3),:);
        p4 = node(elem(neighbor(idx3,3),1),:);
        m14 = (p1 + p4)/2;
        m24 = (p2 + p4)/2;
        m23 = (p2 + p3)/2;
        m1 = (m23 + p3)/2;
        m2 = (m23 + p2)/2;
        isNotTB = ((phiValue(elem(idx3,2)).*phi(m2)) > 0 & ...
                   (phi(m14).*phi(m24) < 0)) | phi(m1).*phi(m2) < 0;
        isTransitionTypeBElem(idx3(isNotTB)) = false;
    end

    %  case 3: a Type B interface element whose second neighbor is a interface
    %          and  Type B element.
    idx4 = find(isTransitionTypeBElem & isInterfaceElem(neighbor(:,2)) & ...
               ~elemType(neighbor(:,2)));
    if ~isempty(idx4)
        p2 = node(elem(idx4,2),:);
        p3 = node(elem(idx4,3),:);
        m23 = (p2 + p3)/2;
        m1 = (m23 + p3)/2;
        m2 = (m23 + p2)/2;
        isNotTB =  phi(m1).*phi(m2) < 0;
        isTransitionTypeBElem(idx4(isNotTB)) = false;
    end
    %  case 4: a Type B interface element whose third neighbor is a interface
    %          and  Type B element.
    idx5 = find(isTransitionTypeBElem & isInterfaceElem(neighbor(:,3)) & ...
                ~elemType(neighbor(:,3)));
    if ~isempty(idx5)
        p2 = node(elem(idx5,2),:);
        p3 = node(elem(idx5,3),:);
        m23 = (p2 + p3)/2;
        m1 = (m23 + p3)/2;
        m2 = (m23 + p2)/2;
        isNotTB =  phi(m1).*phi(m2) < 0;
        isTransitionTypeBElem(idx5(isNotTB)) = false;
    end
    typeB = (isMark & ~elemType & ~isTransitionTypeBElem);
end

% find interface elements
interfaceElemIdx= markinterface(elem,phiValue);
interfaceElem = elem(interfaceElemIdx,:);
% find interface nodes, which are the vertices of interface elements.
isInterfaceNode = false(size(node,1),1);
isInterfaceNode(interfaceElem(:))=true;
T = auxstructure(elem);
edge = T.edge;
edge2elem = T.edge2elem;
clear T;
% case: Fig 2.7 (b)
isEdge = edge2elem(:,1) ~= edge2elem(:,2);
isEdge = isEdge & (isInterfaceNode(edge(:,1)) | isInterfaceNode(edge(:,2)) |...
    isInterfaceNode(elem(edge2elem(:,1),1)) | isInterfaceNode(elem(edge2elem(:,2),1)));
isEdge = isEdge &  edge2elem(:,3) == 1 & edge2elem(:,4)==1;
isEdge = isEdge & ~elemType(edge2elem(:,1));
isMark = false(size(elem,1),1);
isMark(edge2elem(isEdge,1)) = true;
isMark(edge2elem(isEdge,2)) = true;
[node,elem,~, ~,tree] = bisect(node,elem,isMark);
N1 = size(node,1);
phiValue(N0+1:N1) = phi(node(N0+1:N1,:));
% N0 = N1;
elemType = updatetype(elemType,tree(:,[2,3]));

display('Finish Step 3: transform the type B of  element to type A!');

%% Step 4: Move some interface nodes onto the interface
%
%  Here we just move interface nodes along the shortest edges of all the
%  type A interface elements. For a such kind of edge, we get the
%  intersection points and compute the distances between the two endpoints
%  and the intersection points. An interface node maybe have two moving
%  direction, just choose a direction which move shorter distance, then the
%  qualities of the elements connecting to this node will decrease less.

interfaceElemIdx= markinterface(elem,phiValue);
interfaceElem = elem(interfaceElemIdx,:);
T = auxstructure(interfaceElem);
edge = double(T.edge);
edge2elem = T.edge2elem;
isSpecial = elemType(interfaceElemIdx(edge2elem(:,1))) & ...
    ~elemType(interfaceElemIdx(edge2elem(:,2))) & edge2elem(:,3) == 1;
isSpecial = isSpecial | (~elemType(interfaceElemIdx(edge2elem(:,1))) & ...
    elemType(interfaceElemIdx(edge2elem(:,2))) & edge2elem(:,4) == 1);
isShortEdge = (edge2elem(:,1) ~= edge2elem(:,2)) & ~(edge2elem(:,3) == 1 & ...
               edge2elem(:,4) == 1) & ~isSpecial;
A = node(edge(isShortEdge,1),:);
B = node(edge(isShortEdge,2),:);
h = sqrt(sum((A-B).^2,2));
M = findintersectbisect(phi,A,B);
h1 = sqrt(sum((A-M).^2,2));
h2 = h - h1;
t = sparse(edge(isShortEdge,[1,2]),edge(isShortEdge,[2,1]),[h1./h,h2./h], size(node,1),size(node,1));
t1 = sparse(edge(isShortEdge,[1,2]),edge(isShortEdge,[2,1]),repmat(cumsum(ones(size(M,1),1)),1,2), size(node,1),size(node,1));
[maxt,I] = max(t,[],1);
isIrregularNode = maxt >= 0.5;
idx = find(isIrregularNode);
idx1 = I(idx)+(idx - 1)*size(node,1);

oldNodeIdx = [oldNodeIdx;idx'];
oldNode = [oldNode;node(idx,:)];

node(idx,:) = M(t1(idx1),:);

phiValue(idx) = 0;

display('Finish Step 4: move the nearest nodes to interface!');
% figure
% showmesh(node,elem);
% findelem(node,elem,interfaceElemIdx,'noindex');
% plotcurve(a);

%% Step 5: Edge swap
%
%  We just consider the case that two elements share a common longest edge,
%  whose one or two of its four vertices are moved on the interface. After
%  moving, the qualities of this two elements will change(and most case,
%  the quality will decrease). If edge swapping can improve the qualities,
%  we do edge swapping for this two elements.
%
%  But if the common longest edge of the two element is an interface edge,
%  namely its two endpoints are on the interface, we can't use edge
%  swapping. In this case, we deal with it in Step 7. Just find the element
%  patchs of  the first vertices of these two element, and replace this two
%  vertices by the centroids of the element patchs, respectively.

T = auxstructure(elem);
edge = double(T.edge);
edge2elem = T.edge2elem;
clear T;

leftElemIdx = edge2elem(:,1); % the index of the first neighbor of edge
rightElemIdx = edge2elem(:,2); % the index of the second neighbor of  edge

isIrregularNode = full(isIrregularNode)';
isLongestEdge = (edge2elem(:,1) ~= edge2elem(:,2)) & ...
    edge2elem(:,3) == 1 & edge2elem(:,4)== 1 & ...
    (isIrregularNode(edge(:,1)) | isIrregularNode(edge(:,2)) | ...
    isIrregularNode(elem(leftElemIdx,1)) | isIrregularNode(elem(rightElemIdx,1)));

isCrossLongestEdge = (edge2elem(:,1) ~= edge2elem(:,2)) & ...
    edge2elem(:,3) == 1 & edge2elem(:,4)== 1 & ...
    isIrregularNode(elem(leftElemIdx,1)) & isIrregularNode(elem(rightElemIdx,1)) &...
    ~isIrregularNode(edge(:,1)) & ~isIrregularNode(edge(:,2));

isInterfaceEdge = isIrregularNode(edge(:,1)) & isIrregularNode(edge(:,2));
isInterfaceEdge = isLongestEdge & isInterfaceEdge ;
isInterfaceEdge = isInterfaceEdge & ~isIrregularNode(elem(leftElemIdx,1)) ;
isInterfaceEdge = isInterfaceEdge & ~isIrregularNode(elem(rightElemIdx,1));

isLongestEdge(isInterfaceEdge) = false;

leftElemIdx = edge2elem(isLongestEdge,1);
rightElemIdx = edge2elem(isLongestEdge,2);

leftElem = elem(leftElemIdx,:);
rightElem = elem(rightElemIdx,:);

leftElemNew = [leftElem(:,2), rightElem(:,1),leftElem(:,1)];
rightElemNew = [rightElem(:,2), leftElem(:,1),rightElem(:,1)];

l1 = triangleratio(node,leftElem);
r1 = triangleratio(node,rightElem);

l2 = triangleratio(node,leftElemNew);
r2 = triangleratio(node,rightElemNew);

notNeedSwap = min([l1,r1],[],2) >= min([l2,r2],[],2) & ~isCrossLongestEdge(isLongestEdge);

leftElemNew(notNeedSwap,:) = leftElem(notNeedSwap,:);
rightElemNew(notNeedSwap,:) = rightElem(notNeedSwap,:);

oldElemIdx =[oldElemIdx;leftElemIdx;rightElemIdx];
oldElem = [oldElem;elem([leftElemIdx;rightElemIdx],:)];

elem(leftElemIdx,:) = leftElemNew;
elem(rightElemIdx,:) = rightElemNew;

display('Finish Step 5: edge swap!')

%% Step 6: Mesh smoothing for near interface nodes
% q = minangle(node,elem);
% [minq,I] = min(q);

T = auxstructure(elem);
% edge = double(T.edge);
neighbor = double(T.neighbor);
clear T;

isInterfaceNode = (msign(phiValue) == 0);

badElem = find(~isInterfaceNode(elem(:,1)) & isInterfaceNode(elem(:,2)) ...
              & isInterfaceNode(elem(:,3)));
v1 = node(elem(badElem,3),:) - node(elem(badElem,2),:);
v2 = node(elem(badElem,1),:) - node(elem(badElem,3),:);
v3 = node(elem(badElem,2),:) - node(elem(badElem,1),:);
badElem = badElem(sum(v2.^2,2)+sum(v3.^2,2) < sum(v1.^2,2)-eps);
moveNodeIdx = elem(badElem,1);
if ~isempty(moveNodeIdx)
   N = size(node,1);
   NT = size(elem,1);
   v2t = sparse(elem,[(1:NT)',(1:NT)',(1:NT)'],1,N,NT);
   area = simplexvolume(node,elem);
   bc = (node(elem(:,1),:) + node(elem(:,2),:) + node(elem(:,3),:))/3;
   oldNodeIdx = [oldNodeIdx;moveNodeIdx];
   oldNode = [oldNode;node(moveNodeIdx,:)];
   ba = v2t(moveNodeIdx,:)*area;
   node(moveNodeIdx,:) = v2t(moveNodeIdx,:)*([area,area].*bc)./[ba,ba];
end
% case:
badElem2 = find(~isInterfaceNode(elem(:,2)) & isInterfaceNode(elem(:,1)) ...
              & isInterfaceNode(elem(:,3)) & ~elemType(neighbor(:,1)));
v1 = node(elem(badElem2,3),:) - node(elem(badElem2,2),:);
v2 = node(elem(badElem2,1),:) - node(elem(badElem2,3),:);
v3 = node(elem(badElem2,2),:) - node(elem(badElem2,1),:);
badElem2 = badElem2(sum(v2.^2,2)+sum(v3.^2,2)<sum(v1.^2,2)-eps);
% case:
badElem3 = find(~isInterfaceNode(elem(:,3)) & isInterfaceNode(elem(:,1)) ...
              & isInterfaceNode(elem(:,2)) & ~elemType(neighbor(:,1)));
v1 = node(elem(badElem3,3),:) - node(elem(badElem3,2),:);
v2 = node(elem(badElem3,1),:) - node(elem(badElem3,3),:);
v3 = node(elem(badElem3,2),:) - node(elem(badElem3,1),:);
badElem3 = badElem3(sum(v2.^2,2)+sum(v3.^2,2)<sum(v1.^2,2)-eps);
moveNodeIdx = [elem(badElem2,2);elem(badElem3,3)];
if ~isempty(moveNodeIdx)
   if ~exist('v2t','var')
       N = size(node,1);
       NT = size(elem,1);
       v2t = sparse(elem,[(1:NT)',(1:NT)',(1:NT)'],1,N,NT);
   end
   area = simplexvolume(node,elem);
   bc = (node(elem(:,1),:) + node(elem(:,2),:) + node(elem(:,3),:))/3;
   isMany = sum(v2t(moveNodeIdx,:),2)<=6;
   moveNodeIdx = moveNodeIdx(isMany);
   oldNodeIdx = [oldNodeIdx;moveNodeIdx];
   oldNode = [oldNode;node(moveNodeIdx,:)];
   ba = v2t(moveNodeIdx,:)*area;
   node(moveNodeIdx,:) = v2t(moveNodeIdx,:)*([area,area].*bc)./[ba,ba];
end
display('Finish Step 6: We got it finally!');
%
% q = minangle(node,elem);
% [minq,I] = min(q);

repo: ifem


text:
function x = rsc_encode(G,m,termination)
% Copyright 1998, Yufei Wu, MPRG lab, Virginia Tech. for academic use
% encodes a binary data block m (0/1) with a RSC (recursive systematic
% convolutional) code defined by generator matrix G, returns the output
% in x (0/1), terminates the trellis with all-0 state if termination>0
if nargin<3, termination = 0; end
[N,L] = size(G); % Number of output bits, Constraint length
M = L-1; % Dimension of the state
lu = length(m)+(termination>0)*M; % Length of the input
lm = lu-M; % Length of the message
state = zeros(1,M); % initialize the state vector
% To generate the codeword
x = [];
for i = 1:lu
   if termination<=0 | (termination>0 & i<=length(m))
     d_k = m(i);
    elseif termination>0 & i>lm
     d_k = rem(G(1,2:L)*state.',2);
   end
   a_k = rem(G(1,:)*[d_k state].',2);
   xp = rem(G(2,:)*[a_k state].',2); % 2nd output (parity) bits
   state = [a_k state(1:M-1)]; % Next sttate
   x = [x [d_k; xp]]; % since systematic, first output is input bit
end

author: Sable

repo: mcbench-benchmarks


text:
function x = triu( x, k )

%   Disciplined convex/geometric programming information for TRIU:
%       TRIU imposes no convexity restrictions on its arguments.

%
% Check inputs
%

if nargin < 2, k = 0; end
s = x.size_;
if length( s ) > 2,
    error( 'The first argument must be 2-D.' );
elseif ~isnumeric( k ) || length( k ) ~= 1,
    error( 'The second argument must be an integer scalar.' );
end

%
% Zero out the elements outside of the desired triangle
%

b = x.basis_;
b( :, ~triu(ones(s),k) ) = 0;
x = cvx( s, b );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function boxesAll = extendBBoxStickmen(stickmen)
boxesAll = zeros(length(stickmen),4);
for ridx = 1:length(stickmen)
    box = stickmen(ridx).det;
    h = box(4) - box(2);
    % 2x extension in Y
    box(4) = box(4) + h;
    % pad in from all sides
    delta1 = 0.1*h;
    delta2 = 0.1*h;
    box = box + [-delta1 -delta2 delta1 delta2];
    boxesAll(ridx,:) = box;
end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
%FONT_SIZE Set large graphic font
%
%       font_size(size)
%
% Set font size for current figure

function font_size(size)

  V = axis;
  H = get(gcf,'Children');
	c1 = [];
	for h = H'
	if strcmp(get(h,'type'),'axes')
  	set(get(h,'XLabel'), 'FontSize', size);
  	set(get(h,'YLabel'), 'FontSize', size);
  	set(get(h,'Title'),  'FontSize', size);
  	set(h, 'FontSize', size);
  	c1 = [c1; get(gca, 'Children')];
	end
	end
 	axis(V);
  for h1 = c1'
    v1 = get (h1);
    if (isfield (v1, 'FontSize'))
      set (h1, 'FontSize', size);
    end;
    c2 = get (h1, 'Children');
    for h2 = c2'
      v2 = get (h2);
      if (isfield (v2, 'FontSize'))
        set (h2, 'FontSize', size);
      end;
    end;
  end;

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
classdef testNumericalHomogenizerPrinter < testNotShowingError...
                                          & testPrintingComparator

    properties (Access = protected)
        fileOutputName = 'testNumericalHomogenizerPrinter';
        iter = 0;
    end

    properties (Access = private)
        dataBase
    end

    methods (Access = public)

        function obj = testNumericalHomogenizerPrinter()
            d.iter     = 0;
            d.file     = obj.fileOutputName;
            d.postCase = 'NumericalHomogenizer';
            tp = testPrintingPostProcessComputer(d);
            tp.print();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% RES = spyrHigh(PYR, INDICES)
%
% Access the highpass residual band from a steerable pyramid.

% Eero Simoncelli, 6/96.

function res =  spyrHigh(pyr,pind)

res  = pyrBand(pyr, pind, 1);

repo: SelfExSR


text:
function [theta,rho,itheta,prho,rhomin] = getdata(S2G)
% return index of all points in a epsilon neighborhood of a vector
%
% Input
%  S2G     - @S2Grid
%
% Output
%  theta  - double
%  rho    - double
%  itheta - double
%  prho   - double

theta = double(S2G.thetaGrid);
rho = double(S2G.rhoGrid);
itheta = cumsum([0,GridLength(S2G.rhoGrid)]);
prho = S2G.rhoGrid(1).max;
rhomin = S2G.rhoGrid(1).min;
rhomax = S2G.rhoGrid(1).max;

author: mtex-toolbox

repo: mtex

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


text:
function z = min( x, y )

% MAX   Minimum, z = min(x,y)

if isa(x,'tfocs_tuple')
    z = x;
    if isa(y,'tfocs_tuple')
        z.value_ = cellfun( @min, x.value_, y.value_, 'UniformOutput', false );
    elseif isscalar(y)
        z.value_ = cellfun( @min, x.value_, {y}, 'UniformOutput', false );
    else
        z.value_ = cellfun( @min, x.value_, {y}, 'UniformOutput', false );
    end
else
    z = y;
    if isscalar(x)
        z.value_ = cellfun( @min, {x}, y.value_, 'UniformOutput', false );
    else
        z.value_ = cellfun( @min, {x}, y.value_, 'UniformOutput', false );
    end
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
% This program is an implementation of the on-line
% Co-Active Neuro-Fuzzy Inference System (CANFIS) algorithm.
% The structure of the network is determoned by the user.
% The premise (nonlinear) parameters are estimated by Gradient Descent (GD)
% through error backpropagation.
% The consequent (linear) parameters are estimated by Recursive LSE.

clc;
clear;
close all;
colordef black;

NumOfEpochs   = 10;
NumOfSamples = 400;
NumInVars         = 6;
NumOutVars      = 3;
NumInTerms      = 2;
lamda                 = 1;
rmse                   = zeros(NumOutVars,NumOfEpochs);
NumRules          = NumInTerms^NumInVars;
Ita1                      = .0015;
Ita2                      = .0025;
alpha1                = 0.00;
alpha2                = 0.00;
P                         = 1e8*eye((NumInVars+1)*NumRules);

% INPUT-OUTPUT training pair construction.
% We will try to teach the ANFIS network the chaotic

load lorenz_data.mat;
Process  = lorenz_dat; % a (3 X 8056) Lorenz Attractor Solution [X Y Z]' Points Data Matrix.

% Target Output is formulated as a matrix whose lines represent different
% output variables and its columns different time instants.

e         = zeros(NumOutVars,NumOfSamples);
TrgtOut   = zeros(NumOutVars,NumOfSamples);
TrgtOutCh = zeros(NumOutVars,NumOfSamples);

ThetaL4  = zeros((NumInVars+1)*NumRules,NumOutVars);  % Consequent Parameters in Layer 4.
dThetaL4 = zeros((NumInVars+1)*NumRules,NumOutVars);

In1           = zeros(NumInVars,NumOfSamples);
In11         = zeros(NumInVars,NumOfSamples);
Out5        = zeros(NumOutVars,NumOfSamples);
Theta21  = zeros(NumInVars,NumInTerms);
Theta32  = zeros(NumRules,NumRules);
mean1    = zeros(NumInVars,NumInTerms);
sigma1   = zeros(NumInVars,NumInTerms);
dmean1  = zeros(NumInVars,NumInTerms);
dsigma1 = zeros(NumInVars,NumInTerms);
db1          = zeros(NumInVars,NumInTerms);

% Each row corresponds to a different Input Linguistic Variable Xi.
% Each i-j entry corresponds to a different mean (center)
% of the bell-shaped function of the j-th Term of the i-th Input Linguistic
% Variable Xi.

Process0 = Process(:,1:NumOfSamples);
Process1 = Process(:,2:NumOfSamples+1);
Process2 = Process(:,3:NumOfSamples+2);

In1 = [Process1; Process0;];
%In1 = Process1;
TrgtOut = Process2;

Process0ch = Process(:,5*NumOfSamples+1:6*NumOfSamples);
Process1ch = Process(:,5*NumOfSamples+2:6*NumOfSamples+1);
Process2ch = Process(:,5*NumOfSamples+3:6*NumOfSamples+2);

%In11 = [Process1ch; Process0ch;];
In11 = Process1ch;
TrgtOutCh = Process2ch;

for i=1:NumInVars
	 mean1(i,:) = linspace(min(In1(i,:)), max(In1(i,:)),NumInTerms);
	 sigma1(i,:) = ((mean1(i,2)-mean1(i,1))/2)*ones(1,NumInTerms);
end
b1     = 2*ones(NumInVars,NumInTerms);

% plotgenbell(mean1,sigma1,b1,[-60,80]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                			Beginning of the Main "for" loop of the ANFIS  program.        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ne=1:NumOfEpochs

for n=1:NumOfSamples
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        								 							     				                                              %
%   			                  	NETWORK FUNCTIONALITY SECTION						      %
%      																					                                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% LAYER 1 - INPUT TERM NODES
In2 = In1(:,n)*ones(1,NumInTerms);
Out1 = 1./(1 + (abs((In2-mean1)./sigma1)).^(2*b1));

% LAYER 2 - PRODUCT NODES
precond = comb(Out1);
Out2 = prod(precond,2);
S_2 = sum(Out2);

% LAYER 3 - NORMALIZATION NODES
 Out3 = Out2/S_2;

% LAYERS 4 - 5: CONSEQUENT NODES - SUMMING NODE
Aux1 = [In1(:,n); 1]*Out3';

a = reshape(Aux1,(NumInVars+1)*NumRules,1);  % New Input Training Data shaped as a column vector.

Out5(:,n) = ThetaL4'*a;

e(:,n) = TrgtOut(:,n)-Out5(:,n);

%%%%%%%%%%% END OF NETWORK FUNCTIONALITY SECTION %%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%			 					PARAMETER LEARNING SECTION	                    			 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BACKWARD PASS. Error Backpropagation

% We may adjust the Consequent Parameters by LMS.
 for m=1:NumOutVars
     deltaThetaL4 = a*e(m,n);
     dThetaL4(:,m) = Ita1*deltaThetaL4 + alpha1*dThetaL4(:,m);           % Classic LMS.
 end

% LAYER 4
for m=1:NumOutVars
        ThetaL4_mat = reshape(ThetaL4(:,m),NumInVars+1,NumRules);
        f(m,:) = [In1(:,n)' 1]*ThetaL4_mat*e(m,n);
end

% LAYER 3
e3 = sum(f,1);
denom = S_2*S_2;

% LAYER 2
ThetaE32 = zeros(NumRules,NumRules);
if denom~=0
    for k1=1:NumRules
         for k2=1:NumRules
              if k1==k2
                 ThetaE32(k1,k2) = ((S_2-Out2(k2))/denom)*e3(k2);
              else
                 ThetaE32(k1,k2) = -(Out2(k2)/denom)*e3(k2);
              end
         end
    end

else
    continue;
end

e2 = sum(ThetaE32,2);

% LAYER 1
Q = zeros(NumInVars,NumInTerms,NumRules);
for i=1:NumInVars
     for j=1:NumInTerms
          for k=1:NumRules
      	       if Out1(i,j)== precond(k,i) && Out1(i,j)~=0
      			  Q(i,j,k) = (Out2(k)/Out1(i,j))*e2(k);
               end
          end
      end
 end

ThetaE21 = sum(Q,3);

% LAYER 1 PARAMETER ADJUSTMENT BY GRADIENT DESCENT.
deltamean1   = -ThetaE21.*(2*b1./(In2-mean1)).*Out1.*(1-Out1);
deltab1          = -ThetaE21.*(-2).*log(abs((In2-mean1)./sigma1)).*Out1.*(1-Out1);
deltasigma1  = -ThetaE21.*(2*b1./sigma1).*Out1.*(1-Out1);

% dmean1 = -Ita2*deltamean1 + alpha2*dmean1;
% mean1 = mean1 + dmean1;
%
% dsigma1 = -Ita2*deltasigma1 + alpha2*dsigma1;
% sigma1 = sigma1 + dsigma1;
%
% db1 = -Ita2*deltab1 + alpha2*db1;
% b1 = b1 + db1;

% Now update the Layer 4 linear parameters.
ThetaL4_old = ThetaL4;  % Keep here the old values of ThetaL4.
% Fixing of Consequent Parameters by LMS.
ThetaL4  = ThetaL4 + dThetaL4;

% Fixing of Consequent Parameters by Recursive LSE.
% P = (1./lamda).*(P - P*a*a'*P./(lamda+a'*P*a));
% ThetaL4 = ThetaL4 + P*a*e(:,n)';

%%%%%%%%%%%% END OF PARAMETER LEARNING PROCESS %%%%%%%%%%%%

end %Of for i=1:NumOfSamples loop.

for m=1:NumOutVars
    rmse(m,ne) = norm(e(m,:))/sqrt(NumOfSamples);
end

rmse(:,ne)

end   %Of for i=1:NumOfEpochs loop.

figure
subplot(3,1,1),plot(rmse(1,:),'*r');
title('RMSE for X-coordinate');
grid on;

subplot(3,1,2),plot(rmse(2,:),'*g');
title('RMSE for Y-coordinate');
grid on;

subplot(3,1,3),plot(rmse(3,:),'*b');
title('RMSE for Z-coordinate');
grid on;

X = 1:NumOfSamples;
figure;
subplot(3,2,1)
plot(X,TrgtOut(1,:),'r',X,Out5(1,:),'m');
title('Training Target X-coordinate (Red) VS ANFIS X Output (Magenta)')
grid on;

subplot(3,2,2)
plot(X,e(1,:),'r')
title('Training Error(n) for X-coordinate')
grid on;

subplot(3,2,3)
plot(X,TrgtOut(2,:),'b',X,Out5(2,:),'c');
title('Training Target Y-coordinate (Blue) VS ANFIS Y Output (Cyan)')
grid on;

subplot(3,2,4)
plot(X,e(2,:),'b')
title('Training Error(n) for Y-coordinate')
grid on;

subplot(3,2,5)
plot(X,TrgtOut(3,:),'g',X,Out5(3,:),'y');
title('Training Target Z-coordinate (Green) VS ANFIS Z Output (Yellow)')
grid on;

subplot(3,2,6)
plot(X,e(3,:),'g')
title('Training Error(n) for Z-coordinate')
grid on

% plotgenbell(mean1,sigma1,b1,[-60,80]);

%% %%% %%%%%%% Gradient Checking Section %%%%%%%%%%%%%
ThetaJ_mean = zeros(NumInVars,NumInTerms);
ThetaJ_sigma = zeros(NumInVars,NumInTerms);
ThetaJ_b         = zeros(NumInVars,NumInTerms);
epsilon             = 1e-4;

for i=1:NumInVars
    for j=1:NumInTerms

        mean1_up = mean1;
        mean1_do = mean1;

        mean1_up(i,j) = mean1(i,j) + epsilon;
        mean1_do(i,j) = mean1(i,j) - epsilon;

        J_up_mean = sum(1/2*(TrgtOut(:,n) - canfis_grid_forward(In1(:,n),mean1_up,sigma1,b1,ThetaL4_old)).^2);
        J_do_mean = sum(1/2*(TrgtOut(:,n) - canfis_grid_forward(In1(:,n),mean1_do,sigma1,b1,ThetaL4_old)).^2);

        ThetaJ_mean(i,j) = (J_up_mean - J_do_mean)/2/epsilon;

        sigma1_up = sigma1;
        sigma1_do = sigma1;

        sigma1_up(i,j) = sigma1(i,j) + epsilon;
        sigma1_do(i,j) = sigma1(i,j) - epsilon;

        J_up_sigma = sum(1/2*(TrgtOut(:,n) - canfis_grid_forward(In1(:,n),mean1,sigma1_up,b1,ThetaL4_old)).^2);
        J_do_sigma = sum(1/2*(TrgtOut(:,n) - canfis_grid_forward(In1(:,n),mean1,sigma1_do,b1,ThetaL4_old)).^2);

        ThetaJ_sigma(i,j) = (J_up_sigma - J_do_sigma)/2/epsilon;

        b1_up = b1;
        b1_do = b1;

        b1_up(i,j) = b1(i,j) + epsilon;
        b1_do(i,j) = b1(i,j) - epsilon;

        J_up_b = sum(1/2*(TrgtOut(:,n) - canfis_grid_forward(In1(:,n),mean1,sigma1,b1_up,ThetaL4_old)).^2);
        J_do_b = sum(1/2*(TrgtOut(:,n) - canfis_grid_forward(In1(:,n),mean1,sigma1,b1_do,ThetaL4_old)).^2);

        ThetaJ_b(i,j) = (J_up_b - J_do_b)/2/epsilon;

    end
end

rel_diff_mean  = deltamean1(:) - ThetaJ_mean(:);
rel_diff_sigma = deltasigma1(:) - ThetaJ_sigma(:);
rel_diff_b         = deltab1(:) - ThetaJ_b(:);

[rel_diff_mean rel_diff_sigma rel_diff_b]

[deltamean1(:)   ThetaJ_mean(:) ]

[deltasigma1(:)  ThetaJ_sigma(:)]

[deltab1(:)           ThetaJ_b(:)        ]


text:
function [V, policy] = mdp_bellman_operator(P, PR, discount, Vprev)

% mdp_bellman_operator Applies the Bellman operator on the value function Vprev
%                      Returns a new value function and a Vprev-improving policy
% Arguments ---------------------------------------------------------------
% Let S = number of states, A = number of actions
%   P(SxSxA) = transition matrix
%              P could be an array with 3 dimensions or
%              a cell array (1xA), each cell containing a matrix (SxS) possibly sparse
%   PR(SxA) = reward matrix
%              PR could be an array with 2 dimensions or
%              a sparse matrix
%   discount = discount rate, in ]0, 1]
%   Vprev(S) = value function
% Evaluation --------------------------------------------------------------
%   V(S)   = new value function
%   policy(S) = Vprev-improving policy

% MDPtoolbox: Markov Decision Processes Toolbox
% Copyright (C) 2009  INRA
% Redistribution and use in source and binary forms, with or without modification,
% are permitted provided that the following conditions are met:
%    * Redistributions of source code must retain the above copyright notice,
%      this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright notice,
%      this list of conditions and the following disclaimer in the documentation
%      and/or other materials provided with the distribution.
%    * Neither the name of the <ORGANIZATION> nor the names of its contributors
%      may be used to endorse or promote products derived from this software
%      without specific prior written permission.
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
% INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
% OF THE POSSIBILITY OF SUCH DAMAGE.

if iscell(P); S = size(P{1},1); else S = size(P,1); end;
if discount <= 0 || discount > 1
     disp('--------------------------------------------------------')
     disp('MDP Toolbox ERROR: Discount rate must be in ]0; 1]')
     disp('--------------------------------------------------------')
elseif size(Vprev,1) ~= S
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: Vprev must have the same dimension as P')
    disp('--------------------------------------------------------')
else

    if iscell(P)
        A = length(P);
        for a=1:A
            Q(:,a) = PR(:,a) + discount*P{a}*Vprev;
        end
    else
        A = size(P,3);
        for a=1:A
            Q(:,a) = PR(:,a) + discount*P(:,:,a)*Vprev;
        end
    end
    [V, policy] = max(Q,[],2);

end;


text:
function disp(t,name)
%DISP Command window display of a ttensor.
%
%   DISP(T) displays a ttensor with no name.
%
%   DISP(T,NAME) display a ttensor with the given name.
%
%   See also TTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) 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

if ~exist('name','var')
    name = 'ans';
end

fprintf(1,'%s is a ttensor of size %s\n', name, tt_size2str(size(t)));
disp(t.core, sprintf('\t%s.core',name));

for j = 1 : ndims(t)
    fprintf('\t%s.U{%d} = \n', name, j);
    output = tt_matrix2cellstr(t.u{j});
    fprintf('\t\t%s\n',output{:});
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function fem1d_bvp_quadratic_test04 ( )

%*****************************************************************************80
%
%% FEM1D_BVP_QUADRATIC_TEST04 carries out test case #4.
%
%  Discussion:
%
%    Use A4, C4, F4, EXACT4, EXACT_UX4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dianne O'Leary,
%    Scientific Computing with Case Studies,
%    SIAM, 2008,
%    ISBN13: 978-0-898716-66-5,
%    LC: QA401.O44.
%
  n = 11;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_BVP_QUADRATIC_TEST04\n' );
  fprintf ( 1, '  Solve -( A(x) U''(x) )'' + C(x) U(x) = F(x)\n' );
  fprintf ( 1, '  for 0 < x < 1, with U(0) = U(1) = 0.\n' );
  fprintf ( 1, '  A4(X)  = 1.0 + X * X\n' );
  fprintf ( 1, '  C4(X)  = 0.0\n' );
  fprintf ( 1, '  F4(X)  = ( X + 3 X^2 + 5 X^3 + X^4 ) * exp ( X )\n' );
  fprintf ( 1, '  U4(X)  = X * ( 1 - X ) * exp ( X )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes = %d\n', n );
%
%  Geometry definitions.
%
  x_first = 0.0;
  x_last = 1.0;
  x = linspace ( x_first, x_last, n );
  x = x(:);

  u = fem1d_bvp_quadratic ( n, @a4, @c4, @f4, x );

  uexact = exact4 ( x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     I    X         U         Uexact    Error\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %4d  %8f  %8f  %8f  %8e\n', ...
      i, x(i), u(i), uexact(i), abs ( u(i) - uexact(i) ) );
  end

  e1 = l1_error ( n, x, u, @exact4 );
  e2 = l2_error_quadratic ( n, x, u, @exact4 );
  h1s = h1s_error_quadratic ( n, x, u, @exact_ux4 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  l1 norm of error  = %g\n', e1 );
  fprintf ( 1, '  L2 norm of error  = %g\n', e2 );
  fprintf ( 1, '  Seminorm of error = %g\n', h1s );

  return
end
function value = a4 ( x )

%*****************************************************************************80
%
%% A4 evaluates A function #4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of A(X).
%
  value = 1.0 + x .* x;

  return
end
function value = c4 ( x )

%*****************************************************************************80
%
%% C4 evaluates C function #4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of C(X).
%
  value = 0.0;

  return
end
function value = exact4 ( x )

%*****************************************************************************80
%
%% EXACT4 evaluates exact solution #4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of U(X).
%
  value = x .* ( 1.0 - x ) .* exp ( x );

  return
end
function value = exact_ux4 ( x )

%*****************************************************************************80
%
%% EXACT_UX4 evaluates the derivative of exact solution #4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of dUdX(X).
%
  value = ( 1.0 - x - x .* x ) .* exp ( x );

  return
end
function value = f4 ( x )

%*****************************************************************************80
%
%% F4 evaluates right hand side function #4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of F(X).
%
  value = ( x + 3.0 * x.^2 + 5.0 * x.^3 + x.^4 ) .* exp ( x );

  return
end


text:
function dres = build_graph(dres)
ov_thresh = 0.5;
dnum = length(dres.x);
time1 = tic;
len1 = max(dres.fr);
for fr = 2:max(dres.fr)
  if toc(time1) > 2
    fprintf('%0.1f%%\n', 100*fr/len1);
    time1 = tic;
  end
  f1 = find(dres.fr == fr);     %% indices for detections on this frame
  f2 = find(dres.fr == fr-1);   %% indices for detections on the previous frame
  for i = 1:length(f1)
    ovs1  = calc_overlap(dres, f1(i), dres, f2);
    inds1 = find(ovs1 > ov_thresh);                       %% find overlapping bounding boxes.

    ratio1 = dres.h(f1(i))./dres.h(f2(inds1));
    inds2  = (min(ratio1, 1./ratio1) > 0.8);          %% we ignore transitions with large change in the size of bounding boxes.

    dres.nei(f1(i),1).inds  = f2(inds1(inds2))';      %% each detction window will have a list of indices pointing to its neighbors in the previous frame.
%     dres.nei(f1(i),1).ovs   = ovs1(inds1(inds2));
  end
end

author: yuxng


text:
function [ output_maps ] = max_pooling2( input_maps, kernel_size, stride)
%POOLING Summary of this function goes here
%   Detailed explanation goes here

    orig_rows = size(input_maps,1);
    orig_cols = size(input_maps,2);

    pooled_rows = round((orig_rows - kernel_size)/stride) + 1;
    pooled_cols = round((orig_cols - kernel_size)/stride) + 1;

    if(exist('vl_nnpool', 'file') == 3)
        % Caffe and MatConvNet do pooling slightly differently, so need to
        % counter for that

        pooled_cols_vl = floor((orig_cols - kernel_size)/stride) + 1;
        pooled_rows_vl = floor((orig_rows - kernel_size)/stride) + 1;

        if(pooled_rows_vl == pooled_rows && pooled_cols_vl == pooled_cols)
            output_maps = vl_nnpool(input_maps, [kernel_size, kernel_size], 'stride', stride);
        else
            % Else need to pad right and bottom with infinities
            for x=1:kernel_size
                pooled_cols_vl = floor((orig_cols + x - kernel_size)/stride) + 1;
                if(pooled_cols_vl == pooled_cols)
                    break;
                end
            end
            for y=1:kernel_size
                pooled_rows_vl = floor((orig_rows +y - kernel_size)/stride) + 1;
                if(pooled_rows_vl == pooled_rows)
                    break;
                end
            end

            input_maps_new = -inf * ones(size(input_maps,1)+y, size(input_maps,2)+x, size(input_maps,3), size(input_maps,4));
            input_maps_new(1:size(input_maps,1),1:size(input_maps,2),:,:) = input_maps;
            output_maps = vl_nnpool(input_maps_new, [kernel_size, kernel_size], 'stride', stride);
        end
    else

        up_to_rows_out = floor((orig_rows - kernel_size)/stride) + 1;
        up_to_cols_out = floor((orig_cols - kernel_size)/stride) + 1;

        % How many full max-pooling steps are there
        up_to_cols = kernel_size + (up_to_cols_out-1) * stride;
        up_to_rows = kernel_size + (up_to_rows_out-1) * stride;

        output_maps = zeros(pooled_rows, pooled_cols, size(input_maps,3), size(input_maps,4));

        % Pick only the striding elements
        [y, x] = meshgrid(1:up_to_cols-kernel_size+1, 1:up_to_rows-kernel_size+1);
        to_keep_map = mod(y, stride) == 1 & mod(x, stride) == 1;
        to_keep = find(to_keep_map);

        inds_pooling = im2col_inds(input_maps(1:up_to_rows,1:up_to_cols,1,1), [kernel_size, kernel_size]);
        inds_pooling = inds_pooling(:, to_keep);
        for m=1:size(input_maps,4)
            for i=1:size(input_maps,3)
    %             temp = im2col(input_maps(1:up_to_rows,1:up_to_cols,i,m), [kernel_size, kernel_size], 'sliding');
    %             temp = im2col_mine(input_maps(1:up_to_rows,1:up_to_cols,i,m), [kernel_size, kernel_size]);
    %             temp = temp(:,to_keep);

                temp = input_maps(1:up_to_rows,1:up_to_cols,i,m);
                temp = temp(inds_pooling);

                max_val = max(temp);
                output_maps(1:up_to_rows_out,1:up_to_cols_out,i,m) = reshape(max_val, up_to_rows_out, up_to_cols_out);
            end
        end
        % A bit of a hack for non-even number of rows or columns
        if(orig_cols ~= up_to_cols)
            span = orig_cols - (up_to_cols - kernel_size + stride);
            inds_pooling = im2col_inds(input_maps(1:up_to_rows,end-span+1:end,i,m), [kernel_size, span]);
            inds_pooling = inds_pooling(:, 1:stride:end);
            for m=1:size(input_maps,4)
                for i=1:size(input_maps,3)
    %                 temp = im2col(input_maps(1:up_to_rows,end-span+1:end,i,m), [kernel_size, span], 'sliding');
    %                 temp = im2col_mine(input_maps(1:up_to_rows,end-span+1:end,i,m), [kernel_size, span]);
    %                 max_val = max(temp(:,1:stride:end));

                    temp = input_maps(1:up_to_rows,end-span+1:end,i,m);
                    max_val = max(temp(inds_pooling));
                    output_maps(1:up_to_rows_out,end,i,m) = max_val;
                end
            end
        end

        if(orig_rows ~= up_to_rows)
            span = orig_rows - (up_to_rows - kernel_size + stride);
            inds_pooling = im2col_inds(input_maps(end-span+1:end, 1:up_to_cols,i,m), [span, kernel_size]);
            inds_pooling = inds_pooling(:, 1:stride:end);

            for m=1:size(input_maps,4)
                for i=1:size(input_maps,3)
    %                 temp = im2col(input_maps(end-span+1:end, 1:up_to_cols,i,m), [span, kernel_size], 'sliding');
    %                 temp = im2col_mine(input_maps(end-span+1:end, 1:up_to_cols,i,m), [span, kernel_size]);
    %                 max_val = max(temp(:,1:stride:end));
                    temp = input_maps(end-span+1:end, 1:up_to_cols,i,m);
                    max_val = max(temp(inds_pooling));

                    output_maps(end, 1:up_to_cols_out,i,m) = max_val;
                end
            end
        end

        if(orig_cols ~= up_to_cols && orig_rows ~= up_to_rows)
            for m=1:size(input_maps,4)
                for i=1:size(input_maps,3)
                    tmp = input_maps(up_to_rows- kernel_size + stride + 1:end,up_to_cols - kernel_size + stride+1:end,i,m);
                    output_maps(end,end,i,m) = max(tmp(:));
                end
            end
        end

    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
%% RD_FEM solves a 1D reaction/diffusion problem using finite elements.
%
%  Discussion:
%
%    This script sets up and runs a finite element simulation code
%    for a 1D reaction/diffusion equation.
%
%    The dynamics are given by the following reaction/diffusion equation
%    for the function W(T,X):
%
%      W_t = W_xx + NL(W,c),
%
%    where NL(W,c) is a polynomial in W:
%
%       NL(W,c) = c(1) + c(2) * W + c(3) * W^2 + c(4) * W^3
%
%    with Neumann boundary conditions at X = 0.0 and X = 1.0:
%
%      W_x(T,0.0) = 0.0
%      W_x(T,1.0) = 0.0
%
%    and initial condition at T = 0.0:
%
%      W(0,X) = sin ( pi * X ).
%
%    The problem is to be solved for 0.0 <= T <= 4.0, 0.0 <= X <= 1.0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 April 2011
%
%  Author:
%
%   Eugene Cliff
%
%  Reference:
%
%    Jeffrey Borggaard, John Burkardt, John Burns, Eugene Cliff,
%    Working Notes on a Reaction Diffusion Model: a Finite Element Formulation.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RD_FEM:\n' );
  fprintf ( 1, '  Solve a 1D time-dependent reaction/diffusion problem\n' );
  fprintf ( 1, '  with a nonlinear term and Neumann boundary conditions,\n' );
  fprintf ( 1, '  using Neumann boundary conditions.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The equation is discretized using the finite element method.\n' );
%
%  Set initial condition.
%  The first option uses a centered hat.
%
  if ( 0 )
    w_0 = @(x) 1.0 + basic_hat ( 8.0 * x - 4.0 );
  else
    w_0 = @(x) sin(pi*x);
  end
%
%  Set the polynomial coefficients of the reaction term.
%  Here, c = -x + x^3 = x * ( x^2 - 1 )
%
  c = [ 0.0; -1.0; 0.0; 1.0 ];
%
%  Select 101 times for output.
%
  t = ( 0.0 : 0.04 : 4.0 );
%
%  The number of spatial grid points is N+1;
%
  n = 32;
  x = linspace ( 0.0, 1.0, n + 1 );
%
%  T is a vector of 101 times.
%  W is an array of 101x33 function values.
%
  [ T, W ] = rd_lin_spline ( w_0, t, n, c );
%
%  Plot the results.
%
  plot_rd ( t, x, W, c )
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RD_FEM:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );