text:
function inputs = getBatch_LapSRN(opts, imdb, batch, mode)
% -------------------------------------------------------------------------
%   Description:
%       get one batch for training LapSRN
%
%   Input:
%       - opts  : options generated from init_opts()
%       - imdb  : imdb file generated from make_imdb()
%       - batch : array of ID to fetch
%       - mode  : 'train' or 'val'
%
%   Output:
%       - inputs: input for dagnn (include LR and HR images)
%
%   Citation:
%       Deep Laplacian Pyramid Networks for Fast and Accurate Super-Resolution
%       Wei-Sheng Lai, Jia-Bin Huang, Narendra Ahuja, and Ming-Hsuan Yang
%       IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017
%
%   Contact:
%       Wei-Sheng Lai
%       wlai24@ucmerced.edu
%       University of California, Merced
% -------------------------------------------------------------------------

    %% get images
    image_batch = imdb.images.img(batch);

    %% crop patches
    HR = zeros(opts.patch_size, opts.patch_size, 1, length(batch), 'single');

    for i = 1:length(batch)

        img = image_batch{i};

        ratio = 1;

        if( opts.scale_augmentation && strcmp(mode, 'train') )
            % randomly resize between 0.5 ~ 1.0
            ratio = randi([5, 10]) * 0.1;

        end

        eps = 1e-3;

        % min width/height should be larger than patch size
        if size(img, 1) < size(img, 2)
            if size(img, 1) * ratio < opts.patch_size
                ratio = opts.patch_size / size(img, 1) + eps;
            end
        else
            if size(img, 2) * ratio < opts.patch_size
                ratio = opts.patch_size / size(img, 2) + eps;
            end
        end

        img = imresize(img, ratio);

        % random crop
        H = size(img, 1);
        W = size(img, 2);

        y = randi(H - opts.patch_size + 1);
        x = randi(W - opts.patch_size + 1);
        HR(:, :, :, i) = img(y : y + opts.patch_size - 1, x : x + opts.patch_size - 1, :);

    end

    %% data augmentation
    if( opts.data_augmentation && strcmp(mode, 'train') )

        % rotate
        rotate = rand;
        if( rotate < 0.25 )
            HR = rot90(HR, 1);
        elseif( rotate < 0.5 )
            HR = rot90(HR, 2);
        elseif( rotate < 0.75 )
            HR = rot90(HR, 3);
        end

        % horizontally flip
        if( rand > 0.5 )
            HR = fliplr(HR);
        end

    end % end of data augmentation

    %% dagnn input
    inputs = {};
    inputs{end+1} = 'level1_HR';
	inputs{end+1} = HR;

    for i = 2 : opts.level
        ratio = 1 / 2^(i - 1);
        inputs{end+1} = sprintf('level%d_HR', i);
        inputs{end+1} = imresize(HR, ratio);
    end

    inputs{end+1} = 'LR';
	inputs{end+1} = imresize(HR, 1 / opts.scale);

    %% convert to GPU array
    if( opts.gpu > 0 )
        for i = 2:2:length(inputs)
            inputs{i} = gpuArray(inputs{i});
        end
    end

end

repo: LapSRN


text:
function homofil(im,d,r,c,n)
%%%%%%%%%%Butterworth high pass filter %%%%%%%%%%%%%%
A=zeros(r,c);
for i=1:r
    for j=1:c
        A(i,j)=(((i-r/2).^2+(j-c/2).^2)).^(.5);
        H(i,j)=1/(1+((d/A(i,j))^(2*n)));
    end
end
%%%%%%%%%%%%%Using it for my application as homomorphic filtering is
%%%%%%%%%%%%%application specific, taking the value of alphaL and alphaH
%%%%%%%%%%%%%values accordingly.
alphaL=.0999;
aplhaH=1.01;
H=((aplhaH-alphaL).*H)+alphaL;
H=1-H;
%%%%%log of image
im_l=log2(1+im);
%%%%%DFT of logged image
im_f=fft2(im_l);
%%%%%Filter Applying DFT image
im_nf=H.*im_f;
%%%%Inverse DFT of filtered image
im_n=abs(ifft2(im_nf));
%%%%%Inverse log
im_e=exp(im_n);
% subplot(1,2,2);
subplot(122)
imshow((im_e),[])

% figure
% imshow(H)

author: Sable

repo: mcbench-benchmarks


text:
%SerialLink.plot Graphical display and animation
%
% R.plot(Q, options) displays a graphical animation of a robot based on
% the kinematic model.  A stick figure polyline joins the origins of
% the link coordinate frames. The robot is displayed at the joint angle Q (1xN), or
% if a matrix (MxN) it is animated as the robot moves along the M-point trajectory.
%
% Options::
% 'workspace', W    Size of robot 3D workspace, W = [xmn, xmx ymn ymx zmn zmx]
% 'floorlevel',L    Z-coordinate of floor (default -1)
%-
% 'delay',D         Delay betwen frames for animation (s)
% 'fps',fps         Number of frames per second for display, inverse of 'delay' option
% '[no]loop'        Loop over the trajectory forever
% '[no]raise'       Autoraise the figure
% 'movie',M         Save frames as files in the folder M
% 'trail',L         Draw a line recording the tip path, with line style L
%-
% 'scale',S         Annotation scale factor
% 'zoom',Z          Reduce size of auto-computed workspace by Z, makes
%                   robot look bigger
% 'ortho'           Orthographic view
% 'perspective'     Perspective view (default)
% 'view',V          Specify view V='x', 'y', 'top' or [az el] for side elevations,
%                   plan view, or general view by azimuth and elevation
%                   angle.
% 'top'             View from the top.
%-
% '[no]shading'     Enable Gouraud shading (default true)
% 'lightpos',L      Position of the light source (default [0 0 20])
% '[no]name'        Display the robot's name
%-
% '[no]wrist'       Enable display of wrist coordinate frame
% 'xyz'             Wrist axis label is XYZ
% 'noa'             Wrist axis label is NOA
% '[no]arrow'       Display wrist frame with 3D arrows
%-
% '[no]tiles'       Enable tiled floor (default true)
% 'tilesize',S      Side length of square tiles on the floor (default 0.2)
% 'tile1color',C   Color of even tiles [r g b] (default [0.5 1 0.5]  light green)
% 'tile2color',C   Color of odd tiles [r g b] (default [1 1 1] white)
%-
% '[no]shadow'      Enable display of shadow (default true)
% 'shadowcolor',C   Colorspec of shadow, [r g b]
% 'shadowwidth',W   Width of shadow line (default 6)
%-
% '[no]jaxes'       Enable display of joint axes (default false)
% '[no]jvec'        Enable display of joint axis vectors (default false)
% '[no]joints'      Enable display of joints
% 'jointcolor',C    Colorspec for joint cylinders (default [0.7 0 0])
% 'jointdiam',D     Diameter of joint cylinder in scale units (default 5)
%-
% 'linkcolor',C     Colorspec of links (default 'b')
%-
% '[no]base'        Enable display of base 'pedestal'
% 'basecolor',C     Color of base (default 'k')
% 'basewidth',W     Width of base (default 3)
%
% The options come from 3 sources and are processed in order:
% - Cell array of options returned by the function PLOTBOTOPT (if it exists)
% - Cell array of options given by the 'plotopt' option when creating the
%   SerialLink object.
% - List of arguments in the command line.
%
% Many boolean options can be enabled or disabled with the 'no' prefix.  The
% various option sources can toggle an option, the last value is taken.
%
% Graphical annotations and options::
%
% The robot is displayed as a basic stick figure robot with annotations
% such as:
% - shadow on the floor
% - XYZ wrist axes and labels
% - joint cylinders and axes
% which are controlled by options.
%
% The size of the annotations is determined using a simple heuristic from
% the workspace dimensions.  This dimension can be changed by setting the
% multiplicative scale factor using the 'mag' option.
%
% Figure behaviour::
%
% - If no figure exists one will be created and the robot drawn in it.
% - If no robot of this name is currently displayed then a robot will
%   be drawn in the current figure.  If hold is enabled (hold on) then the
%   robot will be added to the current figure.
% - If the robot already exists then that graphical model will be found
%   and moved.
%
% Multiple views of the same robot::
%
% If one or more plots of this robot already exist then these will all
% be moved according to the argument Q.  All robots in all windows with
% the same name will be moved.
%
%  Create a robot in figure 1
%         figure(1)
%         p560.plot(qz);
%  Create a robot in figure 2
%         figure(2)
%         p560.plot(qz);
%  Now move both robots
%         p560.plot(qn)
%
% Multiple robots in the same figure::
%
% Multiple robots can be displayed in the same plot, by using "hold on"
% before calls to robot.plot().
%
%  Create a robot in figure 1
%         figure(1)
%         p560.plot(qz);
%  Make a clone of the robot named bob
%         bob = SerialLink(p560, 'name', 'bob');
%  Draw bob in this figure
%         hold on
%         bob.plot(qn)
%
%  To animate both robots so they move together:
%         qtg = jtraj(qr, qz, 100);
%         for q=qtg'
%           p560.plot(q');
%           bob.plot(q');
%         end
%
% Making an animation movie::
% - The 'movie' options saves frames as files NNNN.png into the specified folder
% - The specified folder will be created
% - To convert frames to a movie use a command like:
%        ffmpeg -r 10 -i %04d.png out.avi
%
% Notes::
% - The options are processed when the figure is first drawn, to make different options come
%   into effect it is neccessary to clear the figure.
% - The link segments do not neccessarily represent the links of the robot, they are a pipe
%   network that joins the origins of successive link coordinate frames.
% - Delay betwen frames can be eliminated by setting option 'delay', 0 or
%   'fps', Inf.
% - By default a quite detailed plot is generated, but turning off labels,
%   axes, shadows etc. will speed things up.
% - Each graphical robot object is tagged by the robot's name and has UserData
%   that holds graphical handles and the handle of the robot object.
% - The graphical state holds the last joint configuration
% - The size of the plot volume is determined by a heuristic for an all-revolute
%   robot.  If a prismatic joint is present the 'workspace' option is
%   required.  The 'zoom' option can reduce the size of this workspace.
%
% See also SerialLink.plot3d, plotbotopt, SerialLink.animate, SerialLink.teach, SerialLink.fkine.

% HANDLES:
%
% A robot comprises a bunch of individual graphical elements and these are
% kept in a structure:
%
%   h.link     the graphical elements that comprise each joint/link
%   h.wrist     the coordinate frame marking the wrist frame
%   h.shadow    the robot's shadow
%   h.floorlevel the z-coordinate of the tiled floor
%   h.robot pointer to the robot object
%   h.opt   the final options structure
%
% The h.link graphical element is tagged with the robot's name and has this
% struct as its UserData.
%
%  h.links -> h -> robot
%
%  This enables us to find all robots with a given name, in all figures,
% and update them.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

% TODO
% deal with base transform and tool
% more consistent option names, scale, mag etc.
function plot(robot, qq, varargin)

    % check the joint angle data matches the robot
    n = robot.n;
    if numcols(qq) ~= n
        error('Insufficient columns in q')
    end

    % process options, these come from:
    %  - passed arguments
    %  - the robot object itself
    %  - the file plotopt.m
    opt = plot_options(robot, varargin);

    % logic to handle where the plot is drawn, are old figures updated or
    % replaced?
    %  calls create_floor() and create_robot() as required.

    if strcmp(get(gca,'Tag'), 'RTB.plot')
        % this axis is an RTB plot window

        rhandles = findobj('Tag', robot.name);

        if isempty(rhandles)
            % this robot doesnt exist here, create it or add it

            if ishold
                % hold is on, add the robot, don't change the floor
                handle = create_robot(robot, opt);

                % tag one of the graphical handles with the robot name and hang
                % the handle structure off it
%                 set(handle.joint(1), 'Tag', robot.name);
%                 set(handle.joint(1), 'UserData', handle);
            else
                % create the robot and floor
                newplot();

                if opt.tiles
                    create_tiled_floor(opt);
                end
                handle = create_robot(robot, opt);
                set(gca, 'Tag', 'RTB.plot');
            end

        end

    else
        % this axis never had a robot drawn in it before, let's use it
        if opt.tiles
            create_tiled_floor(opt);
        end
        handle = create_robot(robot, opt);
        set(gca, 'Tag', 'RTB.plot');
        set(gcf, 'Units', 'Normalized');
        pf = get(gcf, 'Position');
        if strcmp( get(gcf, 'WindowStyle'), 'docked') == 0
            set(gcf, 'Position', [0.1 1-pf(4) pf(3) pf(4)]);
        end
    end

    % deal with a few options that need to be stashed in the SerialLink object
    % movie mode has not already been flagged
    if opt.movie
        robot.framenum = 0;
        robot.moviepath = opt.movie;
    else
        robot.framenum = [];
    end
    robot.delay = opt.delay;
    robot.loop = opt.loop;

    if opt.raise
        % note this is a very time consuming operation
        figure(gcf);
    end

    if strcmp(opt.projection, 'perspective')
        set(gca, 'Projection', 'perspective');
    end

    if isstr(opt.view)
        switch opt.view
            case 'top'
                view(0, 90);
            case 'x'
                view(0, 0);
            case 'y'
                view(90, 0)
            otherwise
                error('rtb:plot:badarg', 'view must be: x, y, top')
        end
    elseif isnumeric(opt.view) && length(opt.view) == 2
        view(opt.view)
    end

    % enable mouse-based 3D rotation
    rotate3d on

    if ~isempty(opt.movie)
        mkdir(opt.movie);
        framenum = 1;
    end
    robot.animate(qq);

end

% Create a new graphical robot in the current figure.
% Returns a structure of handles that describe the various graphical entities in the robot model
% Make extensive use of hgtransform, all entities are defined at the origin, then moved to their
% proper pose
%
% The graphical hiearchy is:
%  hggroup: Tag = robot name
%     hgtransform: Tag = 'link#'
%
% The top-level group has user data which is the handle structure.

function h = create_robot(robot, opt)

    %disp('creating new robot');

    links = robot.links;
    s = opt.scale;

    % create an axis
    ish = ishold();
    if ~ishold
        % if hold is off, set the axis dimensions
        axis(opt.workspace);
        hold on
    end

    N = robot.n;

    % create the base
    if opt.base
        bt = transl(robot.base);
        bt = [bt'; bt'];
        bt(1,3) = opt.floorlevel;
        line(bt(:,1), bt(:,2), bt(:,3), 'LineWidth', opt.basewidth, 'Color', opt.basecolor);
    end

    % add the robot's name
    if opt.name
        b = transl(robot.base);
        bz = 0;
        if opt.base
            bz = 0.5*opt.floorlevel;
        end
        text(b(1), b(2)-s, bz, [' ' robot.name], 'FontAngle', 'italic', 'FontWeight', 'bold')
    end
    group = hggroup('Tag', robot.name);
    h.group = group;

    % create the graphical joint and link elements
    for L=1:N
        if opt.debug
            fprintf('create graphics for joint %d\n', L);
        end

        % create the transform for displaying this element (joint cylinder + link)
        h.link(L) = hgtransform('Tag', sprintf('link%d', L), 'Parent', group);

        % create a joint cylinder
        if opt.joints
            % create the body of the joint
            if links(L).isrevolute
                cyl('z', 2*s, opt.jointdiam/2*s*[-1 1], opt.jointcolor, [], 'Parent', h.link(L));
            else
                % create an additional hgtransform for positioning and scaling the prismatic
                % element.  The element is created with unit length.
                h.pjoint(L) = hgtransform('Tag', 'prismatic', 'Parent', h.link(L));
                if links(L).mdh
                    % make the box extend in negative z-dir because scaling factor in animate
                    % must be positive
                    box('z', s, [0 -1], opt.jointcolor, [], 'Parent', h.pjoint(L));
                else
                    box('z', s, [0 1], opt.jointcolor, [], 'Parent', h.pjoint(L));
                end
            end
        end

        % create the body of the link

        % create elements to represent translation between joint frames, ie. the link itself.
        % This is drawn to resemble orthogonal plumbing.
        if robot.mdh
            % modified DH convention
            if L < N
                A = links(L+1).A(0);
                t = transl(A);
                if t(1) ~= 0
                    cyl('x', s, [0 t(1)], opt.linkcolor, [], 'Parent', h.link(L));
                end
                if t(2) ~= 0
                    cyl('y', s, [0 t(2)], opt.linkcolor, [t(1) 0 0], 'Parent', h.link(L));
                end
                if t(3) ~= 0
                    cyl('z', s, [0 t(3)], opt.linkcolor, [t(1) t(2) 0], 'Parent', h.link(L));
                end
            end
        else
            % standard DH convention
            if L > 1
                Ainv = inv(links(L-1).A(0));
                t = transl(Ainv);
                if t(1) ~= 0
                    cyl('x', s, [0 t(1)], opt.linkcolor, [], 'Parent', h.link(L));
                end
                if t(2) ~= 0
                    cyl('y', s, [s t(2)], opt.linkcolor, [t(1) 0 0], 'Parent', h.link(L));
                end
                if t(3) ~= 0
                    cyl('z', s, [s t(3)], opt.linkcolor, [t(1) t(2) 0], 'Parent', h.link(L));
                end
                %line([0 t(1)]', [0 t(2)]', [0 t(3)]', 'Parent', h.link(L));
            end
        end

        if opt.jaxes && opt.jvec
            error('RTB:plot:badopt', 'Can''t specify ''jaxes'' and ''jvec''')
        end
        % create the joint axis line
        if opt.jaxes
            line('XData', [0 0], ...
                'YData', [0 0], ...
                'ZData', 12*s*[-1 1], ...
                'LineStyle', ':', 'Parent', h.link(L));

            % create the joint axis label
            text(0, 0, 12*s, sprintf('q%d', L), 'Parent', h.link(L))
        end
         % create the joint axis vector
        if opt.jvec
            daspect([1 1 1]);
            ha = arrow3([0 0 -12*s], [0 0 20*s]);
            set(ha, 'Parent', h.link(L));

            % create the joint axis label
            text(0, 0, -12*s, sprintf('q%d', L), 'Parent', h.link(L))
        end
    end
    if opt.debug
        fprintf('create graphics for tool\n');
    end
    % display the tool transform if it exists
    h.link(N+1) = hgtransform('Tag', sprintf('link%d', N+1), 'Parent', group);
    tool = eye(4,4);
    if ~robot.mdh
        tool = links(L).A(0);
    end
    if ~isempty(robot.tool)
        tool = tool * robot.tool;
    end
    t = transl(inv(tool));
    if t(1) ~= 0
        cyl('x', s, [0 t(1)], 'r', [], 'Parent', h.link(N+1));
    end
    if t(2) ~= 0
        cyl('y', s, [s t(2)], 'r', [t(1) 0 0], 'Parent', h.link(N+1));
    end
    if t(3) ~= 0
        cyl('z', s, [s t(3)], 'r', [t(1) t(2) 0], 'Parent', h.link(N+1));
    end

    % display the wrist coordinate frame
    if opt.wrist
        if opt.arrow
            h.wrist = trplot(eye(4,4), 'labels', upper(opt.wristlabel), ...
                'arrow', 'rgb', 'length', 15*s);
        else
            h.wrist = trplot(eye(4,4), 'labels', upper(opt.wristlabel), ...
                'rgb', 'length', 15*s);
        end
    else
        h.wrist = [];
    end

    % display a shadow on the floor
    if opt.shadow
        % create the polyline which is the shadow on the floor
        h.shadow = line('LineWidth', opt.shadowwidth, 'Color', opt.shadowcolor);
    end

    if opt.trail
        h.trail = plot(0, 0, opt.trail);
        robot.trail = [];
    end

    % deal with some display options
    if opt.shading
        lighting gouraud
        light('position', opt.lightpos)
    end
    xlabel('X')
    ylabel('Y')
    zlabel('Z')
    grid on

    % restore hold setting
    if ~ish
        hold off
    end
    h.floorlevel = opt.floorlevel;
    h.robot = robot;
    h.opt = opt;

    % attach the handle structure to the top graphical element
    set(group, 'UserData', h);
end

% draw a cylinder of radius r in the direction specified by ax, with an
% extent from extent(1) to extent(2)

function cyl(ax, r, extent, color, offset, varargin)
    if abs(extent(1) - extent(2)) < eps
        return
    end

    if isempty(offset)
        offset = [0 0 0];
    end

    %fprintf('   cyl: %s, r=%f, extent=[%g, %g]\n', ax, r, extent);

    n = 20;

    r = [r;r];

    theta = (0:n)/n*2*pi;
    sintheta = sin(theta); sintheta(n+1) = 0;

    switch ax
        case 'x'
            y = r * cos(theta) + offset(2);
            z = r * sintheta + offset(3);
            x = extent(:) * ones(1,n+1) + offset(1);
        case 'y'
            x = r * cos(theta) + offset(1);
            z = r * sintheta + offset(3);
            y = extent(:) * ones(1,n+1) + offset(2);
        case 'z'
            x = r * cos(theta) + offset(1);
            y = r * sintheta + offset(2);
            z = extent(:) * ones(1,n+1) + offset(3);
    end

    % walls of the cylinder
    surf(x,y,z, 'FaceColor', color, 'EdgeColor', 'none', varargin{:})

    % put the ends on
    patch(x', y', z', color, 'EdgeColor', 'none', varargin{:});
end

% draw a cylinder of radius r in the direction specified by ax, with an
% extent from extent(1) to extent(2)

function box(ax, r, extent, color, offset, varargin)
    if abs(extent(1) - extent(2)) < eps
        return
    end
    %fprintf('   box: %s, r=%f, extent=[%g, %g]\n', ax, r, extent);
    n = 4;

    r = [r;r];

    theta = (0:n)/n*2*pi;
    sintheta = sin(theta); sintheta(n+1) = 0;

    switch ax
        case 'x'
            y = r * cos(theta);
            z = r * sintheta;
            x = extent(:) * ones(1,n+1);
        case 'y'
            x = r * cos(theta);
            z = r * sintheta;
            y = extent(:) * ones(1,n+1);
        case 'z'
            y = r * cos(theta);
            x = r * sintheta;
            z = extent(:) * ones(1,n+1);
    end

    % walls of the cylinder
    surf(x,y,z, 'FaceColor', color, 'EdgeColor', 'none', varargin{:})

    % put the ends on
    patch(x', y', z', color, 'EdgeColor', 'none', varargin{:});
end

% draw a tiled floor in the current axes
function create_tiled_floor(opt)

    xmin = opt.workspace(1);
    xmax = opt.workspace(2);
    ymin = opt.workspace(3);
    ymax = opt.workspace(4);

    % create a colored tiled floor
    xt = xmin:opt.tilesize:xmax;
    yt = ymin:opt.tilesize:ymax;
    Z = opt.floorlevel*ones( numel(yt), numel(xt));
    C = zeros(size(Z));
    [r,c] = ind2sub(size(C), 1:numel(C));
    C = bitand(r+c,1);
    C = reshape(C, size(Z));
    C = cat(3, opt.tile1color(1)*C+opt.tile2color(1)*(1-C), ...
        opt.tile1color(2)*C+opt.tile2color(2)*(1-C), ...
        opt.tile1color(3)*C+opt.tile2color(3)*(1-C));
    [X,Y] = meshgrid(xt, yt);
    surface(X, Y, Z, C, ...
        'FaceColor','texturemap',...
        'EdgeColor','none',...
        'CDataMapping','direct');
end

    % process a cell array of options and return a struct
    % define all possible options and their default values

function opt = plot_options(robot, optin)

    % timing/looping
    opt.delay = 0.1;
    opt.fps = [];
    opt.loop = false;

    opt.raise = false;

    % general appearance
    opt.scale = 1;
    opt.zoom = 1;
    opt.trail = [];

    opt.workspace = [];
    opt.name = true;
    opt.projection = {'ortho', 'perspective'};
    opt.view = [];
    opt.top = false;

    % 3D rendering
    opt.shading = true;
    opt.lightpos = [0 0 20];

    % tiled floor
    opt.tiles = true;
    opt.tile1color = [0.5 1 0.5];  % light green
    opt.tile2color = [1 1 1];  % white
    opt.floorlevel = [];
    opt.tilesize = 0.2;

    % shadow on the floor
    opt.shadow = true;
    opt.shadowcolor = [0.5 0.5 0.5];
    opt.shadowwidth = 6;

    % the base or pedestal
    opt.base = true;
    opt.basewidth = 3;
    opt.basecolor = 'k';

    % wrist
    opt.wrist = true;
    opt.wristlabel = {'xyz', 'noa'};
    opt.arrow = true;

    % joint rotation axes
    opt.jaxes = false;
    opt.jvec = false;

    % joint cylinders
    opt.joints = true;
    opt.jointdiam = 5;
    opt.jointcolor = [0.7 0 0];

    % links
    opt.linkcolor = 'b';
    opt.toolcolor = 'r';

    % misc
    opt.movie = [];

    % build a list of options from all sources
    %   1. the M-file plotbotopt if it exists
    %   2. robot.plotopt
    %   3. command line arguments
    if exist('plotbotopt', 'file') == 2
        options = [plotbotopt robot.plotopt optin];
    else
        options = [robot.plotopt optin];
    end

    % parse the options
    [opt,args] = tb_optparse(opt, options);
    if ~isempty(args)
        error(['unknown option: ' args{1}]);
    end

    if opt.top
        opt.view = 'top';
    end
    if ~isempty(opt.projection)
        opt.projection = 'ortho';
    end

    % figure the size of the figure
    if isempty(opt.workspace)
        %
        % simple heuristic to figure the maximum reach of the robot
        %
        L = robot.links;
        if any(L.isprismatic)
            error('Prismatic joint(s) present: requires the ''workspace'' option');
        end
        reach = 0;
        for i=1:robot.n
            reach = reach + abs(L(i).a) + abs(L(i).d);
        end
        reach = reach + sum(abs(transl(robot.tool)));
        reach = reach/opt.zoom;

        % if we have a floor, quantize the reach to a tile size
        if opt.tiles
            reach = opt.tilesize * ceil(reach/opt.tilesize);
        end

        % now create a 3D volume based on this reach
        opt.workspace = [-reach reach -reach reach -reach reach];

        % if a floorlevel has been given, ammend the 3D volume
        if ~isempty(opt.floorlevel)
            opt.workspace(5) = opt.floorlevel;
        else
            opt.floorlevel = -reach;
        end
    else
        reach = min(abs(diff(reshape(opt.workspace, [2 3]))));
        if opt.tiles
            % set xy limits to be integer multiple of tilesize
            opt.workspace(1:4) = opt.tilesize * round(opt.workspace(1:4)/opt.tilesize);
            opt.floorlevel = opt.workspace(5);
        end
    end

    % update the fundamental scale factor (given by the user as a multiplier) by a length derived from
    % the overall workspace dimension
    %  we need that a lot when creating the robot model
    opt.scale = opt.scale * reach/40;

    if ~isempty(opt.fps)
        opt.delay = 1/opt.fps;
    end

end

author: Allopart

repo: rbpf-gmapping


text:
function f = dot( F, G )
%DOT   Vector dot product.
%   DOT(F, G) returns the dot product of the DISKFUNV objects F and G.
%   DOT(F, G) is the same as F'*G.
%
% See also CROSS.

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

% Empty check:
if ( isempty( F ) || isempty( G ) )
    f = diskfun();
    return
end

% Extract components:
Fc = F.components;
Gc = G.components;

% Calculate dot-product:
f = times(Fc{1}, Gc{1})+times(Fc{2}, Gc{2});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function test3 (nmat)
%TEST3 test for BTF
% Requires UFget
% Example:
%   test3
% See also btf, maxtrans, strongcomp, dmperm, UFget,
%   test1, test2, test3, test4, test5.

% Copyright 2007, Timothy A. Davis, Univ. of Florida

doplot = 1 ;
dopause = 0 ;
dostrong = 1 ;

index = UFget ;
f = find (index.nrows == index.ncols) ;
[ignore i] = sort (index.nnz (f)) ;
f = f (i) ;
clear i

figure (1)

% short test set: seg faults, lots of blocks, lots of work, and so on:
nasty = [
        % --- various test matrices (no seg fault, quick run time)
    -(1:8)'  % generated matrices
    904 % vanHeukelum/cage3 (5-by-5)
    819 % Simon/raefsky6 (permuted triangular matrix)
        %
        % --- older seg faults:
    264 % HB/west0156, causes older strongcomp_recursive to fail
    824 % TOKAMAK/utm300 (300-by-300), causes older code to fail
    868 % Pothen/bodyy4
        %
        % --- seg faults in old MATLAB dmperm
    290 % Averous/epb3
    983 % Sanghavi/ecl32
    885 % Pothen/tandem_dual
    879 % Pothen/onera_dual
    955 % Schenk_IBMSDS/2D_54019_highK
    957 % Schenk_IBMSDS/3D_51448_3D
    958 % Schenk_IBMSDS/ibm_matrix_2
    912 % vanHeukelum/cage11
    924 % Andrews/Andrews
    960 % Schenk_IBMSDS/matrix-new_3
    862 % Kim/kim1
    544 % Hamm/scircuit
    897 % Norris/torso2
    801 % Ronis/xenon1
     53 % HB/bcsstk31
    958 % Schenk_IBMSDS/matrix_9
    844 % Cunningham/qa8fk
    845 % Cunningham/qa8fk
    821 % Simon/venkat25
    822 % Simon/venkat50
    820 % Simon/venkat01
    812 % Simon/bbmat
    804 % Rothberg/cfd1
     54 % HB/bcsstk32
    913 % vanHeukelum/cage12
    846 % Boeing/bcsstk39
    972 % Schenk_IBMSDS/para-10
    974 % Schenk_IBMSDS/para-5
    975 % Schenk_IBMSDS/para-6
    976 % Schenk_IBMSDS/para-7
    977 % Schenk_IBMSDS/para-8
    978 % Schenk_IBMSDS/para-9
    961 % Schenk_ISEI/barrier2-10
    962 % Schenk_ISEI/barrier2-11
    963 % Schenk_ISEI/barrier2-12
    964 % Schenk_ISEI/barrier2-1
    965 % Schenk_ISEI/barrier2-2
    966 % Schenk_ISEI/barrier2-3
    967 % Schenk_ISEI/barrier2-4
    968 % Schenk_ISEI/barrier2-9
    851 % Chen/pkustk05
    979 % Kamvar/Stanford
    374 % Bova/rma10
        %
        % --- lots of time:
    395 % DRIVCAV/cavity16
    396 % DRIVCAV/cavity17
    397 % DRIVCAV/cavity18
    398 % DRIVCAV/cavity19
    399 % DRIVCAV/cavity20
    400 % DRIVCAV/cavity21
    401 % DRIVCAV/cavity22
    402 % DRIVCAV/cavity23
    403 % DRIVCAV/cavity24
    404 % DRIVCAV/cavity25
    405 % DRIVCAV/cavity26
    1109 % Sandia/mult_dcop_01
    1110 % Sandia/mult_dcop_02
    1111 % Sandia/mult_dcop_03
    376 % Brethour/coater2
    284 % ATandT/onetone2
    588 % Hollinger/mark3jac100
    589 % Hollinger/mark3jac100sc
    452 % Grund/bayer01
    920 % Hohn/sinc12
    590 % Hollinger/mark3jac120
    591 % Hollinger/mark3jac120sc
    809 % Shyy/shyy161
    448 % Graham/graham1
    283 % ATandT/onetone1
    445 % Garon/garon2
    541 % Hamm/bcircuit
    592 % Hollinger/mark3jac140
    593 % Hollinger/mark3jac140sc
    435 % FIDAP/ex40
    912 % Hohn/sinc15
    894 % Norris/lung2
    542 % Hamm/hcircuit
    752 % Mulvey/finan512
    753 % Mulvey/pfinan512
    564 % Hollinger/g7jac180
    565 % Hollinger/g7jac180sc
    566 % Hollinger/g7jac200
    567 % Hollinger/g7jac200sc
    748 % Mallya/lhr34
    749 % Mallya/lhr34c
    922 % Hohn/sinc18
    447 % Goodwin/rim
    807 % Rothberg/struct3
    286 % ATandT/twotone
    982 % Tromble/language
    953 % Schenk_IBMNA/c-73
    890 % Norris/heart1
    750 % Mallya/lhr71
    751 % Mallya/lhr71c
    925 % FEMLAB/ns3Da
    827 % Vavasis/av41092
    931 % FEMLAB/sme3Db
   1297 % GHS_index/boyd2
   1301 % GHS_indef/cont-300
        %
        % --- lots of time, and seg faults:
    285 % ATandT/pre2
        % --- huge matrix, turn off plotting
    940 % Shenk/af_shell1, memory leak in plot, after call to btf, once.
        % ----
]' ;

% maxtrans_recursive causes a seg fault on these matrices, because of
% stack overflow (this is expected)
skip_list_maxtrans_recursive = 285 ;

% p = dmperm (A) in MATLAB 7.4 causes a seg fault on these matrices:
skip_list_dmperm = [285 1301 1231 1251 1232 1241] ;

% [p,q,r] = dmperm (A) in MATLAB 7.4 causes a seg fault on these matrices:
skip_list_dmperm_btf = ...
[    285 879 885 290 955 957 958     924 960         897        959 844 845 ...
 821 822 820     804    913 846 972 974:978 961:968  979 940 ...
 1422 1513 1412 1510 1301 1231 1251 1434 1213 1232 1241 1357 1579 1431 1281] ;
% length(skip_list_dmperm_btf)

% time intensive
skip_costly = [1514 1297 1876 1301] ;

% strongcomp (recursive) causes a seg fault on these matrices because of
% stack overflow (this is expected).
skip_list_strongcomp_recursive     = ...
[983 285 879 885 290 955 957 958 912 924 960 862 544 897 801 53 959 844 845 ...
 821 822 820 812 804 54 913 846 972 974:978 961:968 851     374 940] ;
skip_list_strongcomp_recursive     = ...
[ skip_list_strongcomp_recursive 592 593 752 753 807 286 982 855 566 567 ] ;

% matrices with the largest # of nonzeros in the set (untested)
toobig = [
928   853   852   356   761   368   973   895   805   849   932 ...
803   854   936   802   850   537   856   898   857   859   971   937 ...
914   858   980   896   806   538   863   369   938   860   941   942 ...
943   944   945   946   947   948   915   939   916 ] ;

f = [ -(1:8) f ] ;
% f = nasty ;

h = waitbar (0, 'BTF test 3 of 6') ;

if (nargin < 1)
    nmat = 1000 ;
end
nmat = min (nmat, length (f)) ;
f = f (1:nmat) ;

try

    for matnum = 1:nmat

        waitbar (matnum/nmat, h) ;

        j = f (matnum) ;

        if (any (j == toobig) | any (j == skip_costly))                     %#ok
            fprintf ('\n%4d: %3d %s/%s too big\n', ...
                matnum, j, index.Group{j}, index.Name{j}) ;
            continue ;
        end

        rand ('state', 0) ;

        % clear all unused variables.
        % nothing here is left that is proportional to the matrix size
        clear A p1 p2 p3 q3 r3 match1 match2 match4 pa ra sa qa B C pb rb pc rc
        clear jumble B11 B12 B13 B21 B22 B23 B31 B32 B33 pjumble qjumble ans
        clear c kbad kgood
        % whos
        % pause

        if (j > 0)
            Problem = UFget (j, index) ;
            name = Problem.name ;
            A = Problem.A ;
            clear Problem
        else
            % construct the jth test matrix
            j = -j ;
            if (j == 1 | j == 2)                                            %#ok
                B11 = UFget ('Grund/b1_ss') ;       % 7-by-7 diagonal block
                B11 = B11.A ;
                B12 = sparse (zeros (7,2)) ;
                B12 (3,2) = 1 ;
                B13 = sparse (ones  (7,5)) ;
                B21 = sparse (zeros (2,7)) ;
                B22 = sparse (ones  (2,2)) ;        % 2-by-2 diagonal block
                B23 = sparse (ones  (2,5)) ;
                B31 = sparse (zeros (5,7)) ;
                B32 = sparse (zeros (5,2)) ;
                B33 = UFget ('vanHeukelum/cage3') ;    % 5-by-5 diagonal block
                B33 = B33.A ;
                A = [ B11 B12 B13 ; B21 B22 B23 ; B31 B32 B33 ] ;
                name = '(j=1 test matrix)' ;
            end
            if (j == 2)
                pjumble = [ 10 7 11 1 13 12 8 2 5 14 9 6 4 3 ] ;
                qjumble = [ 3 14 2 11 1 8 5 7 10 12 4 13 9 6 ] ;
                A = A (pjumble, qjumble) ;
                name = '(j=2 test matrix)' ;
            elseif (j == 3)
                A = sparse (1) ;
            elseif (j == 4)
                A = sparse (0) ;
            elseif (j == 5)
                A = sparse (ones (2)) ;
            elseif (j == 6)
                A = sparse (2,2) ;
            elseif (j == 7)
                A = speye (2) ;
            elseif (j == 8)
                A = sparse (2,2) ;
                A (2,1) = 1 ;
            end
            if (j > 2)
                full (A)
            end
        end

        [m n] = size (A) ;
        if (m ~= n)
            continue ;
        end
        fprintf ('\n%4d: ', matnum) ;
        fprintf (' =========================== Matrix: %3d %s\n', j, name) ;
        fprintf ('n: %d nz: %d\n', n, nnz (A)) ;

        if (nnz (A) > 6e6)
            doplot = 0 ;
        end

        %-----------------------------------------------------------------------
        % now try maxtrans
        tic
        match1 = maxtrans (A) ;
        t = toc ;
        s1 = sum (match1 > 0) ;
        fprintf ('n-sprank: %d\n', n-s1) ;
        fprintf ('maxtrans:                %8.2f seconds\n', t) ;
        singular = s1 < n ;

        if (doplot)
            % figure (1)
            clf
            subplot (2,4,1)
            spy (A)
            title (name) ;
        end

        p1 = match1 ;
        if (any (p1 <= 0))
            % complete the permutation
            badrow = find (p1 <= 0) ;

            badcol = ones (1,n) ;
            badcol (p1 (p1 > 0)) = 0 ;
            badcol = find (badcol) ;

            p1 (badrow) = badcol ;

            % construct the older form of match1
            match1 (badrow) = -p1 (badrow) ;
        end
        if (any (sort (p1) ~= 1:n))
            error ('!!') ;
        end

        B = A (:,p1) ;

        if (doplot)
            subplot (2,4,2)
            hold off
            spy (B)
            hold on
            badcol = find (match1 < 0) ;
            Junk = sparse (badcol, badcol, ones (length (badcol), 1), n, n) ;
            % if (~isempty (A))
                % spy (Junk, 'ro') ;
            % end
            title ('maxtrans') ;
        end

        d = nnz (diag (B)) ;
        if (d ~= s1)
            error ('bad sprank') ;
        end
        clear B

        %-----------------------------------------------------------------------
        % try p = dmperm(A)
        skip_dmperm = any (j == skip_list_dmperm) ;

        if (~skip_dmperm)
            tic
            match4 = dmperm (A) ;
            t = toc ;
            fprintf ('p=dmperm(A):             %8.2f seconds\n', t) ;
            s4 = sum (match4 > 0) ;
            singular4 = (s4 < n) ;

            if (doplot)
                if (~singular4)
                    subplot (2,4,3)
                    spy (A (match4,:))
                    title ('dmperm') ;
                end
            end
            if (singular ~= singular4)
                error ('s4?') ;
            end
            if (s1 ~= s4)
                error ('bad sprank') ;
            end
        else
            fprintf ('p=dmperm(A): skip\n') ;
        end

        %-----------------------------------------------------------------------
        nblocks = -1 ;
        skip_dmperm_btf = any (j == skip_list_dmperm_btf) ;
        if (~skip_dmperm_btf)
            % get btf form
            tic
            [pa,qa,ra,sa] = dmperm (A) ;
            t = toc ;
            fprintf ('[p,q,r,s]=dmperm(A):     %8.2f seconds\n', t) ;
            nblocks = length (ra) - 1 ;
            fprintf ('nblocks: %d\n', nblocks) ;
            if (~singular4)
                checkbtf (A, pa, qa, ra) ;
                if (doplot)
                    subplot (2,4,4)
                    drawbtf (A, pa, qa, ra)
                    title ('dmperm blocks')
                end
            end
        else
            fprintf ('[p,q,r,s]=dmperm(A): skip\n') ;
        end

        jumble = randperm (n) ;

        %-----------------------------------------------------------------------
        % try strongcomp, non-recursive version

            %-------------------------------------------------------------------
            % try strongcomp on original matrix
            B = A (:,p1) ;
            tic ;
            [pb,rb] = strongcomp (B) ;
            t = toc ;
            fprintf ('strongcomp               %8.2f seconds\n', t) ;
            if (~singular & ~skip_dmperm_btf & (length (rb) ~= nblocks+1))  %#ok
                error ('BTF:invalid (rb)') ;
            end
            checkbtf (B, pb, pb, rb) ;
            if (doplot)
                subplot (2,4,5)
                drawbtf (B, pb, pb, rb) ;
                title ('strongcomp') ;
            end

            %-------------------------------------------------------------------
            % try btf on original matrix
            tic ;
            [pw,qw,rw] = btf (A) ;
            t = toc ;
            fprintf ('btf                      %8.2f seconds nblocks %d\n', ...
                t, length (rw)-1) ;

            if (any (pw ~= pb))
                error ('pw') ;
            end
            if (any (rw ~= rb))
                error ('rw') ;
            end
            if (any (abs (qw) ~= p1 (pw)))
                error ('qw') ;
            end
            c = diag (A (pw,abs (qw))) ;
            if (~singular & ~skip_dmperm_btf & (length (rw) ~= nblocks+1))  %#ok
                error ('BTF:invalid (rw)') ;
            end
            checkbtf (A, pw, abs (qw), rw) ;

            kbad  = find (qw < 0) ;
            kgood = find (qw > 0) ;
            if (any (c (kbad) ~= 0))
                error ('kbad') ;
            end
            if (any (c (kgood) == 0))       %#ok
                error ('kgood') ;
            end

            if (doplot)
                subplot (2,4,6)
                drawbtf (A, pw, abs (qw), rw) ;
                if (n < 500)
                    for k = kbad
                        plot ([k (k+1) (k+1) k k]-.5, ...
                            [k k (k+1) (k+1) k]-.5, 'r') ;
                    end
                end
                title ('btf') ;
            end

            %-------------------------------------------------------------------
            % try [p,q,r] = strongcomp (A, qin) form
            tic
            [pz,qz,rz] = strongcomp (A, match1) ;
            t = toc ;
            fprintf ('[p,q,r]=strongcomp(A,qin)%8.2f seconds\n', t) ;
            if (any (pz ~= pb))
                error ('pz') ;
            end
            if (any (rz ~= rb))
                error ('rz') ;
            end
            if (any (abs (qz) ~= p1 (pz)))
                error ('qz') ;
            end
            c = diag (A (pz,abs (qz))) ;
            if (~singular & ~skip_dmperm_btf & (length (rz) ~= nblocks+1))  %#ok
                error ('BTF:invalid (rz)') ;
            end
            checkbtf (A, pz, abs (qz), rz) ;

            kbad  = find (qz < 0) ;
            kgood = find (qz > 0) ;
            if (any (c (kbad) ~= 0))
                error ('kbad') ;
            end
            if (any (c (kgood) == 0))                                       %#ok
                error ('kgood') ;
            end

            if (doplot)
                subplot (2,4,7)
                drawbtf (A, pz, abs (qz), rz) ;
                if (n < 500)
                    for k = kbad
                        plot ([k (k+1) (k+1) k k]-.5, ...
                            [k k (k+1) (k+1) k]-.5, 'r') ;
                    end
                end
                title ('strongcomp(A,qin)') ;
            end

            %-------------------------------------------------------------------
            % try strongcomp again, on a randomly jumbled matrix
            C = sparse (B (jumble, jumble)) ;
            tic ;
            [pc,rc] = strongcomp (C) ;
            t = toc ;
            fprintf ('strongcomp       (rand)  %8.2f seconds\n', t) ;
            if (~singular & ~skip_dmperm_btf & (length (rc) ~= nblocks+1))  %#ok
                error ('BTF:invalid (rc)') ;
            end
            checkbtf (C, pc, pc, rc) ;
            if (doplot)
                subplot (2,4,8)
                drawbtf (C, pc, pc, rc) ;
                title ('strongcomp(rand)') ;
            end

            if (length (rc) ~= length (rb))
                error ('strongcomp random mismatch') ;
            end

        %-----------------------------------------------------------------------
        if (doplot)
            drawnow
        end

        if (matnum ~= nmat & dopause)                                       %#ok
            input ('Hit enter: ') ;
        end

    end

catch
    % out-of-memory is OK, other errors are not
    disp (lasterr) ;
    if (isempty (strfind (lasterr, 'Out of memory')))
        error (lasterr) ;                                                   %#ok
    else
        fprintf ('test terminated early, but otherwise OK\n') ;
    end
end

close (h) ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function dim = getDimension(~)
%GETDIMENSION   Returns 3 for a SPINOP3.
%
% See also SPINOP3.

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

dim = 3;

end


text:
function LOCS = tapas_physio_create_LOCS_from_VOLLOCS(VOLLOCS, t, sqpar)
% Compute slice scan events as locations in time vector (LOCS) from volume
% scan locations and sequence parameters (TR, nSlices)
%
% LOCS = tapas_physio_create_LOCS_from_VOLLOCS(VOLLOCS, t, sqpar);
%
%
% IN
%   NOTE: The detailed description of all input structures can be found as
%   comments in tapas_physio_new
%
%   VOLLOCS - index locations in time vector (of physiological recordings),
%                             when volume scan events started
%   t       - time vector of phys time course
%
%   sqpar   - sequence timing parameters, used for computation
%             of scan events from 'nominal' timing
%           .Nslices        - number of slices per volume in fMRI scan
%           .TR             - repetition time in seconds
%           .time_slice_to_slice - time between the acquisition of 2 subsequent
%                             slices; typically TR/Nslices or
%                             minTR/Nslices, if minimal temporal slice
%                             spacing was chosen
%
% OUT
%   LOCS    - locations in time vector, when slice or volume scan
%             events started
%
%   See also tapas_physio_new tapas_physio_main_create_regresssors

% Author: Lars Kasper
% Created: 2013-08-23
% Copyright (C) 2016 Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO 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/>.

NallVols        = sqpar.Ndummies+sqpar.Nscans;
Nslices         = sqpar.Nslices;
TR              = sqpar.TR;
LOCS            = NaN(NallVols*Nslices,1);

%default for equidistantly spaced slices
if isempty(sqpar.time_slice_to_slice)
    sqpar.time_slice_to_slice = TR/Nslices;
end

iFirstValidVolume = find(~isnan(VOLLOCS), 1);
tRef = t(VOLLOCS(iFirstValidVolume));

for n = iFirstValidVolume:NallVols

    % first, restrict search interval to within TR
    LOCVOLSTART = VOLLOCS(n);

    if n == NallVols
        LOCVOLEND = numel(t);
    else
        LOCVOLEND = VOLLOCS(n+1);
    end

    % Then, find slice-by-slice nearest neighbour of actual and
    % reference time
    tSearchInterval = t(LOCVOLSTART:LOCVOLEND) - ...
        (tRef + sqpar.TR*(n-iFirstValidVolume));

    % Note: Why not the same as other slices? Some jitter problem...
    % if TR is not given with high enough precision, therefore reset for
    % each volume!
    LOCS((n-1)*Nslices + 1) = LOCVOLSTART;
    for s = 2:sqpar.Nslices
        [tmp, RELATIVELOC] = min(abs(tSearchInterval - ...
            sqpar.time_slice_to_slice*(s-1)) );

        LOCS((n-1)*Nslices + s) = LOCVOLSTART + RELATIVELOC - 1;
    end
end

author: translationalneuromodeling

repo: tapas


text:
function [ x, w ] = ncoh_compute ( n )

%*****************************************************************************80
%
%% NCOH_COMPUTE computes a Newton-Cotes "open half" quadrature rule.
%
%  Discussion:
%
%    The input value N is used to define N equal subintervals of [-1,+1].
%    The I-th abscissa is the center of the I-th subinterval.
%
%    The integral:
%
%      Integral ( X_MIN <= X <= X_MAX ) F(X) dx
%
%    The quadrature rule:
%
%      Sum ( 1 <= I <= N ) W(I) * F ( X(I) ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%
  x = zeros ( n, 1 );

  x_min = -1.0;
  x_max =  1.0;

  for i = 1 : n
    x(i) = ( ( 2 * n - 2 * i + 1 ) * x_min   ...
           + (         2 * i - 1 ) * x_max ) ...
           / ( 2 * n             );
  end

  w = nc_compute ( n, x_min, x_max, x );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function xfac=dasp(lat)
% DASP sets data aspect ratio to 1
% or if optional lat is supplied, crude projection
% Usage: xfac=dasp(lat);
% where: lat = latitude in degrees
% xfac : the ratio lon/lat
if(nargin==1),
   if (isnan(lat)),
      set(gca,'DataAspectRatioMode','auto');
   else
      xfac=cos(lat*pi/180);
      set (gca, 'DataAspectRatio', [1 xfac 1000] );
   end
else
   set (gca, 'DataAspectRatio', [1 1 1000] );
end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function line_cvt_lloyd_test02 ( )

%*****************************************************************************80
%
%% LINE_CVT_LLOYD_TEST02 tests the constrained computation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_CVT_LLOYD_TEST02:\n' );
  fprintf ( 1, '  Test the constrained computation.\n' );

  n = 25;
  a = 0.0;
  b = 1.0;
  it_num = 200;
  x = a + ( b - a ) * rand ( n, 1 );
  header = 'test02';

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use %d points in the interval [%g,%g]\n', n, a, b );
  fprintf ( 1, '  Take %d iterations.\n', it_num );
  fprintf ( 1, '  Call this calculation "%s"\n', header );
  fprintf ( 1, '  Expect a uniform spacing of %g\n', ( b - a ) / n );

  x = line_ccvt_lloyd ( n, a, b, it_num, header, x );

  return
end


text:
%MDL_TWOLINK_SYM Create symbolic model of a simple 2-link mechanism
%
% MDL_TWOLINK_SYM is a script that creates the workspace variable twolink
% which describes in symbolic form the kinematic and dynamic
% characteristics of a simple planar 2-link mechanism moving in the
% xz-plane, it experiences gravity loading.  The symbolic parameters are:
%  - link lengths: a1, a2
%  - link masses: m1, m2
%  - link CoMs in the link frame x-direction: c1, c2
%  - gravitational acceleration: g
%  - joint angles: q1, q2
%  - joint angle velocities: qd1, qd2
%  - joint angle accelerations: qdd1, qdd2
%
% Notes::
% - It is a planar mechanism operating in the vertical plane and is
%   therefore affected by gravity (unlike mdl_planar2 in the horizontal
%   plane).
% - Gear ratio is 1 and motor inertia is 0.
% - Link inertias Iyy1, Iyy2 are 0.
% - Viscous and Coulomb friction is 0.
%
% References::
%  - Based on Fig 3-6 (p73) of Spong and Vidyasagar (1st edition).
%
% See also mdl_puma560, mdl_stanford, SerialLink.

% MODEL: generic, planar, dynamics, 2DOF, symbolic, standard_DH

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

syms a1 a2 g real
syms c1 c2 m1 m2 real

%syms Iyy1 Iyy2 b1 b2 real
Iyy1 = 0;
Iyy2 = 0;
b1 = 0;
b2 = 0;

twolink = SerialLink([
    Revolute('d', 0, 'a', a1, 'alpha', 0, 'm', m1, 'r', [c1 0 0], 'I', [0 Iyy1 0], 'B', b1, 'G', 1, 'Jm', 0, 'standard')
    Revolute('d', 0, 'a', a2, 'alpha', 0, 'm', m2, 'r', [c2 0 0], 'I', [0 Iyy2 0], 'B', b2, 'G', 1, 'Jm', 0, 'standard')
    ], ...
    'name', 'two link', ...
    'comment', 'from Spong, Hutchinson, Vidyasagar');
twolink = twolink.sym();
twolink.gravity = [0; 0; g];
twolink.base = trotx(pi/2);

syms q1 q2 q1d q2d q1dd q2dd real

author: petercorke

repo: robotics-toolbox-matlab


text:
function forest = forestTrain( data, hs, varargin )
% Train random forest classifier.
%
% Dimensions:
%  M - number trees
%  F - number features
%  N - number input vectors
%  H - number classes
%
% USAGE
%  forest = forestTrain( data, hs, [varargin] )
%
% INPUTS
%  data     - [NxF] N length F feature vectors
%  hs       - [Nx1] or {Nx1} target output labels in [1,H]
%  varargin - additional params (struct or name/value pairs)
%   .M          - [1] number of trees to train
%   .H          - [max(hs)] number of classes
%   .N1         - [5*N/M] number of data points for training each tree
%   .F1         - [sqrt(F)] number features to sample for each node split
%   .split      - ['gini'] options include 'gini', 'entropy' and 'twoing'
%   .minCount   - [1] minimum number of data points to allow split
%   .minChild   - [1] minimum number of data points allowed at child nodes
%   .maxDepth   - [64] maximum depth of tree
%   .dWts       - [] weights used for sampling and weighing each data point
%   .fWts       - [] weights used for sampling features
%   .discretize - [] optional function mapping structured to class labels
%                    format: [hsClass,hBest] = discretize(hsStructured,H);
%
% OUTPUTS
%  forest   - learned forest model struct array w the following fields
%   .fids     - [Kx1] feature ids for each node
%   .thrs     - [Kx1] threshold corresponding to each fid
%   .child    - [Kx1] index of child for each node
%   .distr    - [KxH] prob distribution at each node
%   .hs       - [Kx1] or {Kx1} most likely label at each node
%   .count    - [Kx1] number of data points at each node
%   .depth    - [Kx1] depth of each node
%
% EXAMPLE
%  N=10000; H=5; d=2; [xs0,hs0,xs1,hs1]=demoGenData(N,N,H,d,1,1);
%  xs0=single(xs0); xs1=single(xs1);
%  pTrain={'maxDepth',50,'F1',2,'M',150,'minChild',5};
%  tic, forest=forestTrain(xs0,hs0,pTrain{:}); toc
%  hsPr0 = forestApply(xs0,forest);
%  hsPr1 = forestApply(xs1,forest);
%  e0=mean(hsPr0~=hs0); e1=mean(hsPr1~=hs1);
%  fprintf('errors trn=%f tst=%f\n',e0,e1); figure(1);
%  subplot(2,2,1); visualizeData(xs0,2,hs0);
%  subplot(2,2,2); visualizeData(xs0,2,hsPr0);
%  subplot(2,2,3); visualizeData(xs1,2,hs1);
%  subplot(2,2,4); visualizeData(xs1,2,hsPr1);
%
% See also forestApply, fernsClfTrain
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.24
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% get additional parameters and fill in remaining parameters
dfs={ 'M',1, 'H',[], 'N1',[], 'F1',[], 'split','gini', 'minCount',1, ...
  'minChild',1, 'maxDepth',64, 'dWts',[], 'fWts',[], 'discretize','' };
[M,H,N1,F1,splitStr,minCount,minChild,maxDepth,dWts,fWts,discretize] = ...
  getPrmDflt(varargin,dfs,1);
[N,F]=size(data); assert(length(hs)==N); discr=~isempty(discretize);
minChild=max(1,minChild); minCount=max([1 minCount minChild]);
if(isempty(H)), H=max(hs); end; assert(discr || all(hs>0 & hs<=H));
if(isempty(N1)), N1=round(5*N/M); end; N1=min(N,N1);
if(isempty(F1)), F1=round(sqrt(F)); end; F1=min(F,F1);
if(isempty(dWts)), dWts=ones(1,N,'single'); end; dWts=dWts/sum(dWts);
if(isempty(fWts)), fWts=ones(1,F,'single'); end; fWts=fWts/sum(fWts);
split=find(strcmpi(splitStr,{'gini','entropy','twoing'}))-1;
if(isempty(split)), error('unknown splitting criteria: %s',splitStr); end

% make sure data has correct types
if(~isa(data,'single')), data=single(data); end
if(~isa(hs,'uint32') && ~discr), hs=uint32(hs); end
if(~isa(fWts,'single')), fWts=single(fWts); end
if(~isa(dWts,'single')), dWts=single(dWts); end

% train M random trees on different subsets of data
prmTree = {H,F1,minCount,minChild,maxDepth,fWts,split,discretize};
for i=1:M
  if(N==N1), data1=data; hs1=hs; dWts1=dWts; else
    d=wswor(dWts,N1,4); data1=data(d,:); hs1=hs(d);
    dWts1=dWts(d); dWts1=dWts1/sum(dWts1);
  end
  tree = treeTrain(data1,hs1,dWts1,prmTree);
  if(i==1), forest=tree(ones(M,1)); else forest(i)=tree; end
end

end

function tree = treeTrain( data, hs, dWts, prmTree )
% Train single random tree.
[H,F1,minCount,minChild,maxDepth,fWts,split,discretize]=deal(prmTree{:});
N=size(data,1); K=2*N-1; discr=~isempty(discretize);
thrs=zeros(K,1,'single'); distr=zeros(K,H,'single');
fids=zeros(K,1,'uint32'); child=fids; count=fids; depth=fids;
hsn=cell(K,1); dids=cell(K,1); dids{1}=uint32(1:N); k=1; K=2;
while( k < K )
  % get node data and store distribution
  dids1=dids{k}; dids{k}=[]; hs1=hs(dids1); n1=length(hs1); count(k)=n1;
  if(discr), [hs1,hsn{k}]=feval(discretize,hs1,H); hs1=uint32(hs1); end
  if(discr), assert(all(hs1>0 & hs1<=H)); end; pure=all(hs1(1)==hs1);
  if(~discr), if(pure), distr(k,hs1(1))=1; hsn{k}=hs1(1); else
      distr(k,:)=histc(hs1,1:H)/n1; [~,hsn{k}]=max(distr(k,:)); end; end
  % if pure node or insufficient data don't train split
  if( pure || n1<=minCount || depth(k)>maxDepth ), k=k+1; continue; end
  % train split and continue
  fids1=wswor(fWts,F1,4); data1=data(dids1,fids1);
  [~,order1]=sort(data1); order1=uint32(order1-1);
  [fid,thr,gain]=forestFindThr(data1,hs1,dWts(dids1),order1,H,split);
  fid=fids1(fid); left=data(dids1,fid)<thr; count0=nnz(left);
  if( gain>1e-10 && count0>=minChild && (n1-count0)>=minChild )
    child(k)=K; fids(k)=fid-1; thrs(k)=thr;
    dids{K}=dids1(left); dids{K+1}=dids1(~left);
    depth(K:K+1)=depth(k)+1; K=K+2;
  end; k=k+1;
end
% create output model struct
K=1:K-1; if(discr), hsn={hsn(K)}; else hsn=[hsn{K}]'; end
tree=struct('fids',fids(K),'thrs',thrs(K),'child',child(K),...
  'distr',distr(K,:),'hs',hsn,'count',count(K),'depth',depth(K));
end

function ids = wswor( prob, N, trials )
% Fast weighted sample without replacement. Alternative to:
%  ids=datasample(1:length(prob),N,'weights',prob,'replace',false);
M=length(prob); assert(N<=M); if(N==M), ids=1:N; return; end
if(all(prob(1)==prob)), ids=randperm(M,N); return; end
cumprob=min([0 cumsum(prob)],1); assert(abs(cumprob(end)-1)<.01);
cumprob(end)=1; [~,ids]=histc(rand(N*trials,1),cumprob);
[s,ord]=sort(ids); K(ord)=[1; diff(s)]~=0; ids=ids(K);
if(length(ids)<N), ids=wswor(cumprob,N,trials*2); end
ids=ids(1:N)';
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function test_failed = test_libltfat_solvehermitiansystem(varargin)
test_failed = 0;

fprintf(' ===============  %s ================ \n',upper(mfilename));

definput.flags.complexity={'double','single'};
[flags]=ltfatarghelper({},definput,varargin);
dataPtr = [flags.complexity, 'Ptr'];

Larr =    [301, 9,11,110, 9, 8, 11, 10, 301];

    complexstring = 'complex';
    funname_init = makelibraryname('hermsystemsolver_init',flags.complexity,1);
    funname_execute = makelibraryname('hermsystemsolver_execute',flags.complexity,1);
    funname_done = makelibraryname('hermsystemsolver_done',flags.complexity,1);

    p = libpointer();
    calllib('libltfat',funname_init,max(Larr),p);

for Lidx = 1:numel(Larr)
    L = Larr(Lidx);

    D = randn(L,flags.complexity)+1i*randn(L,flags.complexity);
   % D(:,1) = D(:,end);
    A = D*D';
    Amessedup = A;
    for n=2:size(A,2)
        for m=1:n-1
              Amessedup(m,n) = randn(1);
        end
    end
    Aint = complex2interleaved(Amessedup);
    APtr = libpointer(dataPtr,Aint);

    b = randn(L,1,flags.complexity)+ 1i*randn(L,1,flags.complexity);
    bint = complex2interleaved(b);
    bPtr = libpointer(dataPtr,bint);

    tic
    trueres = A\b;
    toc

    tic
    status = calllib('libltfat',funname_execute,p,APtr,L,bPtr);
    toc

    res = norm(trueres - interleaved2complex(bPtr.Value));

    [test_failed,fail]=ltfatdiditfail(res+status,test_failed,1e-8);
    fprintf(['SOLVEHERM L:%3i, %s %s %s %s\n'],L,flags.complexity,complexstring,ltfatstatusstring(status),fail);

end

    calllib('libltfat',funname_done,p);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function sbmlModel = convertCobraToSBML(model, sbmlLevel, sbmlVersion, compSymbolList, compNameList)
% Converts a cobra structure to an sbml structure - kept for backward Compatability
%
% USAGE:
%
%    sbmlModel = convertCobraToSBML(model, sbmlLevel, sbmlVersion, compSymbolList, compNameList)
%
% INPUT:
%    model:             COBRA model structure
%
% OPTIONAL INPUTS:
%    sbmlLevel:         SBML Level (ignored! level 3 is used)
%    sbmlVersion:       SBML Version (ignored! version 1 is used)
%    compSymbolList:    List of compartment symbols
%    compNameList:      List of copmartment names correspoding to `compSymbolList`
%
% OUTPUT:
%    sbmlModel:         SBML MATLAB structure
%
% .. Authors:
%    - Longfei Mao 24/09/15 FBCv2 support added
%    - Thomas Pfau 19/09/2016 reinstantiated for backward compatability
% NOTE:
%
%    sbmlLevel and sbmlVersion are ignored by this function and it is just
%    kept for backward compatability and might be deprecated in future
%    versions.

warning('This function is likely to be deprecated in future COBRA Toolbox version.')

TemporaryFileName = [tempname '.xml'];

if ~exist('compSymbolList','var') || ~exist('compNameList','var')
    sbmlModel = writeSBML(model,TemporaryFileName);
else
    sbmlModel = writeSBML(model,TemporaryFileName,compSymbolList,compNameList);
end

delete(TemporaryFileName)

author: opencobra

repo: cobratoolbox


text:
function net = res_imagenet_init(n, varargin)

opts.batchNormalization = true;
opts.networkType = 'resnet'; % 'plain' | 'resnet'
opts.bottleneck = false; % only used when n is an array
opts.nClasses = 1000;
opts.reLUafterSum = true;
opts = vl_argparse(opts, varargin);
nClasses = opts.nClasses;

net = dagnn.DagNN();

% n -> specific configuration
if numel(n)==4,
  Ns = n;
else
  switch n,
    case 18, Ns = [2 2 2 2]; opts.bottleneck = false;
    case 34, Ns = [3 4 6 3]; opts.bottleneck = false;
    case 50, Ns = [3 4 6 3]; opts.bottleneck = true;
    case 101, Ns = [3 4 23 3]; opts.bottleneck = true;
    case 152, Ns = [3 8 36 3]; opts.bottleneck = true;
    otherwise, error('No configuration found for n=%d', n);
  end
end
if strcmpi(opts.networkType, 'plain') && opts.bottleneck,
  error('plain network cannot be built with bottleneck layers');
end

% Meta parameters
net.meta.inputSize = [224 224 3] ;
net.meta.trainOpts.weightDecay = 0.0001 ;
net.meta.trainOpts.momentum = 0.9;
net.meta.trainOpts.batchSize = 256 ;
if opts.batchNormalization;
  net.meta.trainOpts.learningRate = [0.1*ones(1,30) 0.01*ones(1,30) 0.001*ones(1,50)] ;
else
  net.meta.trainOpts.learningRate = [0.01*ones(1,45) 0.001*ones(1,45) 0.0001*ones(1,75)] ;
end
net.meta.trainOpts.numEpochs = numel(net.meta.trainOpts.learningRate) ;

% First conv layer
switch n
    case {18, 34, 50}
        block = dagnn.Conv('size',  [7 7 3 64], 'hasBias', true, ...
            'stride', 2, 'pad', [3 3 3 3]);
        lName = 'conv0';
net.addLayer(lName, block, 'data', lName, {[lName '_f'], [lName '_b']});
    case {101, 152}
        block = dagnn.Conv('size',  [7 7 3 64], 'hasBias', false, ...
            'stride', 2, 'pad', [3 3 3 3]);
        lName = 'conv0';
    net.addLayer(lName, block, 'data', lName, {[lName '_f']});
end

add_layer_bn(net, 64, lName, 'bn0', 0.1);
block = dagnn.ReLU('leak',0);
net.addLayer('relu0',  block, 'bn0', 'relu0');

%add_block_conv(net, '0', 'image', [7 7 3 64], 2, opts.batchNormalization, true);
block = dagnn.Pooling('poolSize', [3 3], 'method', 'max', 'pad', [0 1 0 1], 'stride', 2);
net.addLayer('pool0', block, 'relu0', 'pool0');

info.lastNumChannel = 64;
info.lastIdx = 0;
info.lastName = 'pool0';

% Four groups of layers
info = add_group(opts.networkType, net, Ns(1), info, 3, 64,  1, opts.bottleneck, opts.batchNormalization, opts.reLUafterSum);
info = add_group(opts.networkType, net, Ns(2), info, 3, 128, 2, opts.bottleneck, opts.batchNormalization, opts.reLUafterSum);
info = add_group(opts.networkType, net, Ns(3), info, 3, 256, 2, opts.bottleneck, opts.batchNormalization, opts.reLUafterSum);
info = add_group(opts.networkType, net, Ns(4), info, 3, 512, 2, opts.bottleneck, opts.batchNormalization, opts.reLUafterSum);

% Prediction & loss layers
block = dagnn.Pooling('poolSize', [7 7], 'method', 'avg', 'pad', 0, 'stride', 1);
if opts.reLUafterSum
    net.addLayer('pool_final', block, sprintf('relu%d',info.lastIdx), 'pool_final');
else
    net.addLayer('pool_final', block, sprintf('sum%d',info.lastIdx), 'pool_final');
end
block = dagnn.Conv('size', [1 1 info.lastNumChannel nClasses], 'hasBias', true, ...
                   'stride', 1, 'pad', 0);
lName = sprintf('fc%d', info.lastIdx+1);
net.addLayer(lName, block, 'pool_final', lName, {[lName '_f'], [lName '_b']});

net.addLayer('softmax', dagnn.SoftMax(), lName, 'softmax');
net.addLayer('loss', dagnn.Loss('loss', 'log'), {'softmax', 'label'}, 'loss');
net.addLayer('error', dagnn.Loss('loss', 'classerror'), {'softmax','label'}, 'error') ;
net.addLayer('error5', dagnn.Loss('loss', 'topkerror', 'opts', {'topK', 5}), ...
  {'softmax','label'}, 'error5') ;

net.initParams();

net.meta.normalization.imageSize = net.meta.inputSize;
net.meta.normalization.border = 256 - net.meta.inputSize(1:2) ;
net.meta.normalization.interpolation = 'bicubic' ;
net.meta.normalization.averageImage = [] ;
net.meta.normalization.keepAspect = true ;
net.meta.augmentation.rgbVariance = zeros(0,3) ;
net.meta.augmentation.transformation = 'stretch' ;

% Add a group of layers containing 2n/3n conv layers
function info = add_group(netType, net, n, info, w, ch, stride, bottleneck, bn, reLUafterSum)
if strcmpi(netType, 'plain'),
  if isfield(info, 'lastName'),
    lName = info.lastName;
    info = rmfield(info, 'lastName');
  else
    lName = sprintf('relu%d', info.lastIdx);
  end
  add_block_conv(net, sprintf('%d', info.lastIdx+1), lName, ...
    [w w info.lastNumChannel ch], stride, bn, true);
  info.lastIdx = info.lastIdx + 1;
  info.lastNumChannel = ch;
  for i=2:2*n,
    add_block_conv(net, sprintf('%d', info.lastIdx+1), sprintf('relu%d', info.lastIdx), ...
      [w w ch ch], 1, bn, true);
    info.lastIdx = info.lastIdx + 1;
  end
elseif strcmpi(netType, 'resnet'),
    info = add_block_res(net, info, [w w info.lastNumChannel ch], stride, bottleneck, bn, 1, reLUafterSum);
    for i=2:n,
        if bottleneck,
            info = add_block_res(net, info, [w w 4*ch ch], 1, bottleneck, bn, 0, reLUafterSum);
        else
            info = add_block_res(net, info, [w w ch ch], 1, bottleneck, bn, 0, reLUafterSum);
        end
    end
end

% Add a smallest residual unit (2/3 conv layers)
function info = add_block_res(net, info, f_size, stride, bottleneck, bn, isFirst, reLUafterSum)

if isfield(info, 'lastName'),
  lName0 = info.lastName;
  info = rmfield(info, 'lastName');
elseif reLUafterSum || info.lastIdx == 0
  lName0 = sprintf('relu%d',info.lastIdx);
else
    lName0 = sprintf('sum%d',info.lastIdx);
end

lName01 = lName0;
if stride > 1 || isFirst,
  if bottleneck,
      ch = 4*f_size(4);
  else
      ch = f_size(4);
  end
  block = dagnn.Conv('size',[1 1 f_size(3) ch], 'hasBias',false,'stride',stride, ...
    'pad', 0);
  lName_tmp = lName0;
  lName0 = [lName_tmp '_down2'];
  net.addLayer(lName0, block, lName_tmp, lName0, [lName0 '_f']);

  pidx = net.getParamIndex([lName0 '_f']);
  net.params(pidx).learningRate = 0;

  add_layer_bn(net, ch, lName0, [lName01 '_d2bn'], 0.1);
  lName0 = [lName01 '_d2bn'];
end

if bottleneck,
  add_block_conv(net, sprintf('%d',info.lastIdx+1), lName01, [1 1 f_size(3) f_size(4)], stride, bn, true);
  info.lastIdx = info.lastIdx + 1;
  info.lastNumChannel = f_size(4);
  add_block_conv(net, sprintf('%d',info.lastIdx+1), sprintf('relu%d',info.lastIdx), ...
    [f_size(1) f_size(2) info.lastNumChannel info.lastNumChannel], 1, bn, true);
  info.lastIdx = info.lastIdx + 1;
  add_block_conv(net, sprintf('%d',info.lastIdx+1), sprintf('relu%d',info.lastIdx), ...
    [1 1 info.lastNumChannel info.lastNumChannel*4], 1, bn, true);
  info.lastIdx = info.lastIdx + 1;
  info.lastNumChannel = info.lastNumChannel*4;
else
  add_block_conv(net, sprintf('%d',info.lastIdx+1), lName01, f_size, stride, bn, true);
  info.lastIdx = info.lastIdx + 1;
  info.lastNumChannel = f_size(4);
  add_block_conv(net, sprintf('%d',info.lastIdx+1), sprintf('relu%d',info.lastIdx), ...
    [f_size(1) f_size(2) info.lastNumChannel info.lastNumChannel], 1, bn, true);
  info.lastIdx = info.lastIdx + 1;
end
if bn,
  lName1 = sprintf('bn%d', info.lastIdx);
else
  lName1 = sprintf('conv%d', info.lastIdx);
end

info.lastIdx = info.lastIdx + 1;
net.addLayer(sprintf('sum%d',info.lastIdx), dagnn.Sum(), {lName0,lName1}, ...
sprintf('sum%d',info.lastIdx));

% relu
if reLUafterSum
    block = dagnn.ReLU('leak', 0);
    net.addLayer(sprintf('relu%d', info.lastIdx), block, sprintf('sum%d', info.lastIdx), ...
    sprintf('relu%d', info.lastIdx));
end

% Add a conv layer (followed by optional batch normalization & relu)
function net = add_block_conv(net, out_suffix, in_name, f_size, stride, bn, relu)
block = dagnn.Conv('size',f_size, 'hasBias',false, 'stride', stride, ...
                   'pad',[ceil(f_size(1)/2-0.5) floor(f_size(1)/2-0.5) ...
                   ]);
lName = ['conv' out_suffix];
net.addLayer(lName, block, in_name, lName, {[lName '_f']});

if bn,
  add_layer_bn(net, f_size(4), lName, strrep(lName,'conv','bn'), 0.1);
  lName = strrep(lName, 'conv', 'bn');
end
if relu,
  block = dagnn.ReLU('leak',0);
  net.addLayer(['relu' out_suffix], block, lName, ['relu' out_suffix]);
end

% Add a batch normalization layer
function net = add_layer_bn(net, n_ch, in_name, out_name, lr)
block = dagnn.BatchNorm('numChannels', n_ch);
net.addLayer(out_name, block, in_name, out_name, ...
  {[out_name '_g'], [out_name '_b'], [out_name '_m']});
pidx = net.getParamIndex({[out_name '_g'], [out_name '_b'], [out_name '_m']});
net.params(pidx(1)).weightDecay = 0;
net.params(pidx(2)).weightDecay = 0;
net.params(pidx(3)).learningRate = lr;
net.params(pidx(3)).trainMethod = 'average';

repo: ResNet-Matconvnet


text:
function X=rectsample(N)

X=rand(N,2);
X=X+repmat([-0.5,-0.5],N,1);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


function varargout = plot_line(varargin)
% PLOT_LINE M-file for plot_line.fig
%      PLOT_LINE, by itself, creates a new PLOT_LINE or raises the existing
%      singleton*.
%
%      H = PLOT_LINE returns the handle to a new PLOT_LINE or the handle to
%      the existing singleton*.
%
%      PLOT_LINE('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in PLOT_LINE.M with the given input arguments.
%
%      PLOT_LINE('Property','Value',...) creates a new PLOT_LINE or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before plot_line_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to plot_line_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help plot_line

% Last Modified by GUIDE v2.5 27-Mar-2004 18:44:35

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @plot_line_OpeningFcn, ...
                   'gui_OutputFcn',  @plot_line_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before plot_line is made visible.
function plot_line_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to plot_line (see VARARGIN)

% Choose default command line output for plot_line
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes plot_line wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = plot_line_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --------------------------------------------------------------------
function FileOpen_Callback(hObject, eventdata, handles)
% hObject    handle to FileOpen (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get the file to open
[filename, pathname] = uigetfile('*.dat','Load Data');
if filename ~= 0

   % Open the input file
   filename = [pathname filename];
   [fid,msg] = fopen(filename,'rt');

   % Check to see of the open failed.
   if fid < 0

      % There was an error--tell user.
      str = ['File ' filename ' could not be opened.'];
      title = 'File Open Failed';
      errordlg(str,title,'modal');

   else

      % File opened successfully. Read the (x,y) pairs from
      % the input file.  Get first (x,y) pair before the
      % loop starts.
      [in,count] = fscanf(fid,'%g',2);
      ii = 0;

      while ~feof(fid)
         ii = ii + 1;
         x(ii) = in(1);
         y(ii) = in(2);

         % Get next (x,y) pair
         [in,count] = fscanf(fid,'%g',2);
      end

      % Data read in.  Close file.
      fclose(fid);

      % Now plot the data.
      hline = plot(x,y,'LineWidth',3);
      xlabel('x');
      ylabel('y');
      grid on;

      % Associate the context menu with line
      set(hline,'Uicontextmenu',handles.ContextMenu1);

      % Save the line's handle as application data
      handles.hline = hline;
      guidata(gcbf, handles);

   end
end

% --------------------------------------------------------------------
function FileExit_Callback(hObject, eventdata, handles)
% hObject    handle to FileExit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

close(gcbf);

% --------------------------------------------------------------------
function LineSolid_Callback(hObject, eventdata, handles)
% hObject    handle to LineSolid (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

set(handles.hline,'LineStyle','-');

% --------------------------------------------------------------------
function LineDashed_Callback(hObject, eventdata, handles)
% hObject    handle to LineDashed (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

set(handles.hline,'LineStyle','--');

% --------------------------------------------------------------------
function LineDotted_Callback(hObject, eventdata, handles)
% hObject    handle to LineDotted (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

set(handles.hline,'LineStyle',':');

% --------------------------------------------------------------------
function LineDashDot_Callback(hObject, eventdata, handles)
% hObject    handle to LineDashDot (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

set(handles.hline,'LineStyle','-.');

% --------------------------------------------------------------------
function File_Callback(hObject, eventdata, handles)
% hObject    handle to File (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function ContextMenu1_Callback(hObject, eventdata, handles)
% hObject    handle to ContextMenu1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


text:
function params = init_default_params(params)

% Initialize default parameters

if ~isfield(params, 'interpolation_method')
    params.interpolation_method = 'none';
end
if ~isfield(params, 'interpolation_centering')
    params.interpolation_centering = false;
end
if ~isfield(params, 'interpolation_windowing')
    params.interpolation_windowing = false;
end
if ~isfield(params, 'clamp_position')
    params.clamp_position = false;
end
if ~isfield(params, 'update_projection_matrix')
    params.update_projection_matrix = true;
end
if ~isfield(params, 'proj_init_method')
    params.proj_init_method = 'pca';
end
if ~isfield(params, 'use_detection_sample')
    params.use_detection_sample = true;
end
if ~isfield(params, 'use_projection_matrix')
    params.use_projection_matrix = true;
end
if ~isfield(params, 'use_sample_merge')
    params.use_sample_merge = true;
end
if ~isfield(params, 'CG_use_FR')
    params.CG_use_FR = true;
end
if ~isfield(params, 'CG_standard_alpha')
    params.CG_standard_alpha = true;
end
if ~isfield(params, 'neglect_higher_frequency')
    params.neglect_higher_frequency = false;
end

if params.use_projection_matrix == false
    params.proj_init_method = 'none';
    params.update_projection_matrix = false;
end

repo: CFWCR


text:
function sim = gaussianKernel(x1, x2, sigma)
    %% RBFKERNEL returns a radial basis function kernel between x1 and x2
    %   sim = gaussianKernel(x1, x2) returns a gaussian kernel between x1 and x2
    %   and returns the value in sim

    % Ensure that x1 and x2 are column vectors
    x1 = x1(:); x2 = x2(:);

    % ====================== YOUR CODE HERE ======================
    % Instructions: Fill in this function to return the similarity between x1
    %               and x2 computed using a Gaussian kernel with bandwidth
    %               sigma
    value = - sum((x1 - x2) .^ 2) / (2 * sigma ^ 2);
    sim = exp(value);
end

author: worldveil

repo: coursera-ml


text:
% SPGL1: A solver for large-scale sparse reconstruction
% Version 1.7 (May 20, 2009)
%
% Files
%   spg_bp           - Solve the basis pursuit (BP) problem
%   spg_bpdn         - Solve the basis pursuit denoise (BPDN) problem
%   spg_lasso        - Solve the LASSO problem
%   spg_mmv          - Solve the multiple measurement vector (MMV) problem
%   spg_group        - Solve the group-sparse BPDN problem
%   spgdemo          - Demonstrates the use of the SPGL1 solver
%   spgl1            - Expert driver for BP, BPDN, and LASSO problems
%   spgSetParms      - Set options for SPGL1
%   spgsetup         - Compile the SPGL1 MEX interfaces
%
% $Id: Contents.m 1384 2009-05-20 19:09:37Z ewout78 $

author: mpf


text:
disp('Starting tests...')
test_digamma
%test_gammaln
test_trigamma
test_int_hist
%test_inv_posdef
test_java
test_logmulexp
test_logsumexp
test_mutable
%test_ndsum
%test_normcdf
test_normpdf
test_randgamma
test_randbeta
test_randbinom
test_randwishart
test_repmat
test_row_sum
test_sameobject
test_scale
test_solve_tri
test_sorted
test_sqdist
disp('All tests completed')

author: Cloud-CV

repo: object-proposals


text:
function c = tapas_softmax_wld_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the softmax observation model for multinomial responses with phasic
% volatility exp(mu3) as the decision temperature and parameters accounting for win- and
% loss-distortion of state values.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2019 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/>.

% Config structure
c = struct;

% Is the decision based on predictions or posteriors? Comment as appropriate.
c.predorpost = 1; % Predictions
%c.predorpost = 2; % Posteriors

% Model name
c.model = 'softmax_wld';

% Sufficient statistics of Gaussian parameter priors

% Beta
c.logbemu = log(1);
c.logbesa = 4^2;

% Win-distortion
c.la_wdmu = 0;
c.la_wdsa = 2^-2;

% Loss-distortion
c.la_ldmu = 0;
c.la_ldsa = 2^-2;

% Gather prior settings in vectors
c.priormus = [
    c.logbemu,...
    c.la_wdmu,...
    c.la_ldmu,...
         ];

c.priorsas = [
    c.logbesa,...
    c.la_wdsa,...
    c.la_ldsa,...
         ];

% Model filehandle
c.obs_fun = @tapas_softmax_wld;

% Handle to function that transforms observation parameters to their native space
% from the space they are estimated in
c.transp_obs_fun = @tapas_softmax_wld_transp;

return;


text:
%% (Internal) Convert matrix of ECG wave annotations to a struct position format, used in wavedet algorithm
%
%   this_position = matrix2positions(pos_matrix, fields)
%
% Arguments:
%
%      + pos_matrix: cell with the fields defined in fields variable
%
% Output:
%
%      + this_position: struct with fields vectors of annotations
%
%      + fields: field names in pos_matrix
%
% Example:
%
% See also positions2matrix
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/05/2014
% Birthdate  : 18/05/2016
% Copyright 2008-2016
%
function this_position = matrix2positions(pos_matrix, fields)

    this_position = [];

    if(~iscell(pos_matrix))
        pos_matrix = {pos_matrix};
    end

    fields = rowvec(cellstr(fields));

    cant_leads = length(pos_matrix);
    cant_fields = length(fields);

    for ii = 1:cant_leads

        aux_matrix = pos_matrix{ii};
        aux_struct = [];

        for count = 1:cant_fields
            aux_struct.(fields{count}) = colvec(aux_matrix(:,count));
        end

        if( isempty(this_position) )
            this_position = aux_struct;
        else
            this_position(ii) = aux_struct;
        end

    end

end

author: marianux

repo: ecg-kit

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


text:
function [epochs] = fiff_read_epochs(fname)
%
% [epochs] = fiff_read_epochs(fname,setno)
%
% Read eochs from file
%
%
%   Author : Martin Luessi, MGH Martinos Center
%   License : BSD 3-clause

global FIFF;
if isempty(FIFF)
    FIFF = fiff_define_constants();
end

me='MNE:fiff_read_epochs';

%
%   Open the file
%
fprintf(1,'Reading %s ...\n',fname);
[ fid, tree ] = fiff_open(fname);
%
%   Read the measurement info
%
[ info, meas ] = fiff_read_meas_info(fid,tree);
info.filename = fname;

%
%   Read events
%
[ events, mappings ] = fiff_read_events(fid,tree);

%
% Locate the data of interest
%
processed = fiff_dir_tree_find(meas, FIFF.FIFFB_PROCESSED_DATA);
if length(processed) == 0
        fclose(fid);
    error(me,'Could not find epochs data');
end

ep = fiff_dir_tree_find(meas, FIFF.FIFFB_MNE_EPOCHS);
if length(ep) == 0
        fclose(fid);
    error(me,'Could not find epochs data');
end

comment = '';
selection = '';
drop_log = '';
for k = 1:ep.nent
    kind = ep.dir(k).kind;
    pos  = ep.dir(k).pos;
    switch kind
        case FIFF.FIFF_FIRST_SAMPLE
            tag = fiff_read_tag(fid,pos);
            first = tag.data;
        case FIFF.FIFF_LAST_SAMPLE
            tag = fiff_read_tag(fid,pos);
            last = tag.data;
        case FIFF.FIFF_COMMENT
            tag = fiff_read_tag(fid,pos);
            comment = tag.data;
        case FIFF.FIFF_EPOCH
            tag = fiff_read_tag(fid,pos);
            epoch = tag.data;
        case FIFF.FIFF_MNE_BASELINE_MIN
            tag = fiff_read_tag(fid,pos);
            bmin = tag.data;
        case FIFF.FIFF_MNE_BASELINE_MAX
            tag = fiff_read_tag(fid,pos);
            bmax = tag.data;
        case FIFF.FIFFB_MNE_EPOCHS_SELECTION
            tag = fiff_read_tag(fid,pos);
            selection = tag.data;
        case FIFF.FIFFB_MNE_EPOCHS_DROP_LOG
            tag = fiff_read_tag(fid,pos);
            drop_log = tag.data;
    end
end

if ~exist('epoch','var')
    fclose(fid);
    error(me,'Epochs data not found');
end

if ~exist('bmin','var')
    bmin = double(first)/info.sfreq;
end

if ~exist('bmax','var')
    bmax = double(last)/info.sfreq;
end

baseline = [bmin, bmax];

nsamp = last-first+1;
fprintf(1,'\tFound the data of interest:\n');
fprintf(1,'\t\tt = %10.2f ... %10.2f ms (%s)\n',...
    1000*double(first)/info.sfreq,1000*double(last)/info.sfreq,comment);
if ~isempty(info.comps)
    fprintf(1,'\t\t%d CTF compensation matrices available\n',length(info.comps));
end

if size(epoch,1) ~= size(events,1)
    fclose(fid);
    error(me,'Incorrect number of trials (%d instead of %d)',...
        size(epoch,1),size(events,1));
end

if size(epoch,2) ~= info.nchan
    fclose(fid);
    error(me,'Incorrect number of channels (%d instead of %d)',...
        size(epoch,2),info.nchan);
end

if size(epoch,3) ~= nsamp
    fclose(fid);
    error(me,'Incorrect number of samples (%d instead of %d)',...
        size(epoch,3),nsamp);
end

%
%   Calibrate
%
for k = 1:info.nchan
    cals(k) = info.chs(k).cal;
end

nepochs = size(epoch, 1);
epoch = repmat(cals, [nepochs, 1, nsamp]) .* epoch;

times = double(first:last) / info.sfreq;
tmin = times(1);
tmax = times(end);

%
% Put it all together
%
epochs.info = info;
epochs.events = events;
epochs.name = comment;
epochs.times = times;
epochs.tmin = tmin;
epochs.tmax = tmax;
epochs.data = epoch;
epochs.baseline = baseline;
epochs.event_id = mappings;
epochs.selection = selection;
epochs.drop_log = drop_log;

fclose(fid);

return;


text:
function bsv_test ( )

%*****************************************************************************80
%
%% BSV_TEST tests the BURGERS_STEADY_VISCOUS (BSV) library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BSV_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the BURGERS_STEADY_VISCOUS (BSV) library.\n' );

  bsv_test01 ( );
  bsv_test02 ( );
  bsv_test03 ( );
  bsv_test04 ( );
  bsv_test05 ( );
  bsv_test06 ( );
  bsv_test07 ( );
  bsv_test08 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BSV_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function [fathom] = ft2fathom(ft)
% Convert length from feet to fathoms.
% Chad A. Greene 2012
fathom = ft*0.1666666666667;

repo: mcbench-benchmarks


text:
function [pairs,pOri,cOri] = getP2CPairs(job,varargin)
%

% all parent - child neighbors
pairs = neighbors(job.grains);

isParent = ismember(pairs, job.grains.id(job.isParent));
isChild = ismember(pairs, job.grains.id(job.isChild));

ind = any(isParent,2) & any(isChild,2);
pairs = pairs(ind,:);

% ensure parent is first
doFlip = isChild(ind,1);
pairs(doFlip,:) = fliplr(pairs(doFlip,:));

if  check_option(varargin,'quick') && length(pairs) > 10000
  ind = unique(randi(length(pairs),10000,1));
  pairs = pairs(ind,:);
end

% maybe there is nothing to do
if isempty(pairs)
  pOri = orientation(job.csParent);
  cOri = orientation(job.csChild);
  return
end

% compute the corresponding mean orientations
if job.useBoundaryOrientations

  % identify boundaries by grain pairs
  [gB,pairId] = job.grains.boundary.selectByGrainId(pairs);

  % extract boundary child orientations
  % TODO: this can be done faster
  pOri = job.ebsd('id',gB.ebsdId(:,1)).orientations;
  cOri = job.ebsdPrior('id',gB.ebsdId(:,2)).orientations;

  % average child orientations along the boundaries
  pOri = accumarray(pairId,pOri);
  cOri = accumarray(pairId,cOri);

else

  % simply the mean orientations of the grains
  pOri = job.grains('id',pairs(:,1)).meanOrientation;
  cOri = job.grains('id',pairs(:,2)).meanOrientation;

end

% translate to index if required
if check_option(varargin,'index'), pairs = job.grains.id2ind(pairs); end

end

author: mtex-toolbox

repo: mtex

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


text:
            % Diferite reprezentari grafice speciale
            % **************************************
			% ***** Grafica nr.1. *****
	subplot(4,2,1)
			% Desenarea unui poligon
			% - desemnarea coordonatelor varfurilor
	x=[0 2 4 3 1 0];
	y=[0 1 2 4 3 0];
			% - reprezentarea poligonului
	fill(x,y,'y')
			% - plasarea unui titlu
	title('FILL')
			% ***** Grafica nr.2. *****
	subplot(4,2,2)
	x=0:pi/10:2*pi;
	y=sin(x);
			% Reprezentarea grafica utilizand comanda bar
	bar(x,y)
	title('BAR')
			% - fixarea limitelor axelor de coordonate
	axis([0 2*pi -1.1 1.1])
			% ***** Grafica nr.3. *****
	subplot(4,2,3)
			% Reprezentarea grafica utilizand comanda stem
	stem(x,y)
	title('STEM')
	axis([0 2*pi -1.1 1.1])
			% ***** Grafica nr.4. *****
	subplot(4,2,4)
			% Reprezentarea grafica utilizand comanda stairs
	stairs(x,y)
	title('STAIRS')
	axis([0 2*pi -1.1 1.1])
			% ***** Grafica nr.5. *****
	subplot(4,2,5)
			% Generarea unui vector de eroare aleatoare
	e=rand(size(x)).*y./5;
			% Reprezentarea grafica utilizand comanda errorbar
	errorbar(x,y,e)
	title('ERRORBAR')
	axis([0 2*pi -1.1 1.1])
			% ***** Grafica nr.6. *****
	subplot(4,2,6)
			% Regenerarea vectorului x cu un pas mai mic
	x=0:pi/100:2*pi;
			% Reincarcarea vectorului y
	y=sin(x);
			% Reprezentarea grafica utilizand comanda hist
	hist(y,20)
	title('HIST')
            % ***** Grafica nr.7. *****
	subplot(4,2,7)
			% Reprezentarea grafica utilizand comanda pie
	pie([2 4 3 6],{'Nord','Sud','Est','Vest'})
	title('PIE')
            % ***** Grafica nr.8. *****
	subplot(4,2,8)
			% Reprezentarea grafica utilizand comanda area
	z = [1,5,6;3,4,7;2,8,3;2,6,10];
	area(z)
	title('AREA')

author: Sable


text:
clear all, close all, clc
load ../DATA/VORTALL;
% VORTALL contains flow fields reshaped into column vectors
X = VORTALL;
[Phi, Lambda, b] = DMD(X(:,1:end-1),X(:,2:end),21);
% Code to plot the second mode
fhandle = plotCylinder(real(reshape(Phi(:,2),199,449)));

author: dynamicslab