text:
function [p,w] = nrbeval(nurbs,tt)
%
% Function Name:
%
%   nrbeval - Evaluate a NURBS at parameteric points
%
% Calling Sequence:
%
%   [p,w] = nrbeval(crv,ut)
%   [p,w] = nrbeval(srf,{ut,vt})
%
% Parameters:
%
%   crv		: NURBS curve, see nrbmak.
%
%   srf		: NURBS surface, see nrbmak.
%
%   ut		: Parametric evaluation points along U direction.
%
%   vt		: Parametric evaluation points along V direction.
%
%   p		: Evaluated points on the NURBS curve or surface as cartesian
% 		coordinates (x,y,z). If w is included on the lhs argument list
% 		the points are returned as homogeneous coordinates (wx,wy,wz).
%
%   w		: Weights of the homogeneous coordinates of the evaluated
% 		points. Note inclusion of this argument changes the type
% 		of coordinates returned in p (see above).
%
% Description:
%
%   Evaluation of NURBS curves or surfaces at parametric points along the
%   U and V directions. Either homogeneous coordinates are returned if the
%   weights are requested in the lhs arguments, or as cartesian coordinates.
%   This function utilises the 'C' interface bspeval.
%
% Examples:
%
%   Evaluate the NURBS circle at twenty points from 0.0 to 1.0
%
%   nrb = nrbcirc;
%   ut = linspace(0.0,1.0,20);
%   p = nrbeval(nrb,ut);
%
% See:
%
%     bspeval
%

%  D.M. Spink
%  Copyright (c) 2000.

if nargin < 2
  error('Not enough input arguments');
end

foption = 1;    % output format 3D cartesian coordinates
if nargout == 2
  foption = 0;  % output format 4D homogenous coordinates
end

if ~isstruct(nurbs)
  error('NURBS representation is not structure!');
end

if ~strcmp(nurbs.form,'B-NURBS')
  error('Not a recognised NURBS representation');
end

if iscell(nurbs.knots)
  % NURBS structure represents a surface

  num1 = nurbs.number(1);
  num2 = nurbs.number(2);
  degree = nurbs.order-1;

  if iscell(tt)
    % Evaluate over a [u,v] grid
    % tt{1} represents the u direction
    % tt{2} represents the v direction

    nt1 = length(tt{1});
    nt2 = length(tt{2});

    % Evaluate along the v direction
    val = reshape(nurbs.coefs,4*num1,num2);
    val = bspeval(degree(2),val,nurbs.knots{2},tt{2});
    val = reshape(val,[4 num1 nt2]);

    % Evaluate along the u direction
    val = permute(val,[1 3 2]);
    val = reshape(val,4*nt2,num1);
    val = bspeval(degree(1),val,nurbs.knots{1},tt{1});
    val = reshape(val,[4 nt2 nt1]);
    val = permute(val,[1 3 2]);

    w = val(4,:,:);
    p = val(1:3,:,:);
    if foption
      p = p./repmat(w,[3 1 1]);
    end

  else

    % Evaluate at scattered points
    % tt(1,:) represents the u direction
    % tt(2,:) represents the v direction

    nt = size(tt,2);

    val = reshape(nurbs.coefs,4*num1,num2);
    val = bspeval(degree(2),val,nurbs.knots{2},tt(2,:));
    val = reshape(val,[4 num1 nt]);

    % evaluate along the u direction
    pnts = zeros(4,nt);
    for v = 1:nt
      coefs = squeeze(val(:,:,v));
      pnts(:,v) = bspeval(degree(1),coefs,nurbs.knots{1},tt(1,v));
    end

    w = pnts(4,:);
    p = pnts(1:3,:);
    if foption
       p = p./w;
    end

  end

else

  % NURBS structure represents a curve
  %  tt represent a vector of parametric points in the u direction

  val = bspeval(nurbs.order-1,nurbs.coefs,nurbs.knots,tt);

  w = val(4,:);
  p = val(1:3,:);
  if foption
    p = p./repmat(w,3,1);
  end

end


text:
function hh = rotateXLabels( ax, angle, varargin )
%rotateXLabels: rotate any xticklabels
%
%   hh = rotateXLabels(ax,angle) rotates all XLabels on axes AX by an angle
%   ANGLE (in degrees). Handles to the resulting text objects are returned
%   in HH.
%
%   hh = rotateXLabels(ax,angle,param,value,...) also allows one or more
%   optional parameters to be specified. Possible parameters are:
%     'MaxStringLength'   The maximum length of label to show (default inf)
%
%   Examples:
%   >> bar( hsv(5)+0.05 )
%   >> days = {'Monday','Tuesday','Wednesday','Thursday','Friday'};
%   >> set( gca(), 'XTickLabel', days )
%   >> rotateXLabels( gca(), 45 )
%
%   See also: GCA, BAR

%   Copyright 2006-2013 The MathWorks Ltd.

error( nargchk( 2, inf, nargin ) );
if ~isnumeric( angle ) || ~isscalar( angle )
    error( 'RotateXLabels:BadAngle', 'Parameter ANGLE must be a scalar angle in degrees' )
end
angle = mod( angle, 360 );

% From R2014b, rotating labels is built-in using 'XTickLabelRotation'
if ~verLessThan('matlab','8.4.0')
    set(ax, 'XTickLabelRotation', angle)
    if nargout > 1
        hh = [];
    end
    return;
end

[maxStringLength] = parseInputs( varargin{:} );

% Get the existing label texts and clear them
[vals, labels] = findAndClearExistingLabels( ax, maxStringLength );

% Create the new label texts
h = createNewLabels( ax, vals, labels, angle );

% Reposition the axes itself to leave space for the new labels
repositionAxes( ax );

% If an X-label is present, move it too
repositionXLabel( ax );

% Store angle
setappdata( ax, 'RotateXLabelsAngle', angle );

% Only send outputs if requested
if nargout
    hh = h;
end

%-------------------------------------------------------------------------%
    function [maxStringLength] = parseInputs( varargin )
        % Parse optional inputs
        maxStringLength = inf;
        if nargin > 0
            params = varargin(1:2:end);
            values = varargin(2:2:end);
            if numel( params ) ~= numel( values )
                error( 'RotateXLabels:BadSyntax', 'Optional arguments must be specified as parameter-value pairs.' );
            end
            if any( ~cellfun( 'isclass', params, 'char' ) )
                error( 'RotateXLabels:BadSyntax', 'Optional argument names must be specified as strings.' );
            end
            for pp=1:numel( params )
                switch upper( params{pp} )
                    case 'MAXSTRINGLENGTH'
                        maxStringLength = values{pp};

                    otherwise
                        error( 'RotateXLabels:BadParam', 'Optional parameter ''%s'' not recognised.', params{pp} );
                end
            end
        end
    end % parseInputs
%-------------------------------------------------------------------------%
    function [vals,labels] = findAndClearExistingLabels( ax, maxStringLength )
        % Get the current tick positions so that we can place our new labels
        vals = get( ax, 'XTick' );

        % Now determine the labels. We look first at for previously rotated labels
        % since if there are some the actual labels will be empty.
        ex = findall( ax, 'Tag', 'RotatedXTickLabel' );
        if isempty( ex )
            % Store the positions and labels
            labels = get( ax, 'XTickLabel' );
            if isempty( labels )
                % No labels!
                return
            else
                if ~iscell(labels)
                    labels = cellstr(labels);
                end
            end
            % Clear existing labels so that xlabel is in the right position
            set( ax, 'XTickLabel', {}, 'XTickMode', 'Manual' );
            setappdata( ax, 'OriginalXTickLabels', labels );
        else
            % Labels have already been rotated, so capture them
            labels = getappdata( ax, 'OriginalXTickLabels' );
            set(ex, 'DeleteFcn', []);
            delete(ex);
        end
        % Limit the length, if requested
        if isfinite( maxStringLength )
            for ll=1:numel( labels )
                if length( labels{ll} ) > maxStringLength
                    labels{ll} = labels{ll}(1:maxStringLength);
                end
            end
        end

    end % findAndClearExistingLabels
%-------------------------------------------------------------------------%
    function restoreDefaultLabels( ax )
        % Restore the default axis behavior
        removeListeners( ax );

        % Try to restore the tick marks and labels
        set( ax, 'XTickMode', 'auto', 'XTickLabelMode', 'auto' );
        rmappdata( ax, 'OriginalXTickLabels' );

        % Try to restore the axes position
        if isappdata( ax, 'OriginalAxesPosition' )
            set( ax, 'Position', getappdata( ax, 'OriginalAxesPosition' ) );
            rmappdata( ax, 'OriginalAxesPosition' );
        end
    end
%-------------------------------------------------------------------------%
    function textLabels = createNewLabels( ax, vals, labels, angle )
        % Work out the ticklabel positions
        zlim = get(ax,'ZLim');
        z = zlim(1);

        % We want to work in normalised coords, but this doesn't print
        % correctly. Instead we have to work in data units even though it
        % makes positioning hard.
        y = getYPositionToUse( ax );

        % Now create new text objects in similar positions.
        textLabels = -1*ones( numel( vals ), 1 );
        for ll=1:numel(vals)
            textLabels(ll) = text( ...
                'Units', 'Data', ...
                'Position', [vals(ll), y, z], ...
                'String', labels{ll}, ...
                'Parent', ax, ...
                'Clipping', 'off', ...
                'Rotation', angle, ...
                'Tag', 'RotatedXTickLabel', ...
                'UserData', vals(ll));
        end
        % So that we can respond to CLA and CLOSE, attach a delete
        % callback. We only attach it to one label to save massive numbers
        % of callbacks during axes shut-down.
        set(textLabels(end), 'DeleteFcn', @onTextLabelDeleted);

        % Now copy font properties into the texts
        updateFont();
        % Update the alignment of the text
        updateAlignment();

    end % createNewLabels

%-------------------------------------------------------------------------%
    function repositionAxes( ax )
        % Reposition the axes so that there's room for the labels
        % Note that we only do this if the OuterPosition is the thing being
        % controlled
        if ~strcmpi( get( ax, 'ActivePositionProperty' ), 'OuterPosition' )
            return;
        end

        % Work out the maximum height required for the labels
        labelHeight = getLabelHeight(ax);

        % Remove listeners while we mess around with things, otherwise we'll
        % trigger redraws recursively
        removeListeners( ax );

        % Change to normalized units for the position calculation
        oldUnits = get( ax, 'Units' );
        set( ax, 'Units', 'Normalized' );

        % Not sure why, but the extent seems to be proportional to the height of the axes.
        % Correct that now.
        set( ax, 'ActivePositionProperty', 'Position' );
        pos = get( ax, 'Position' );
        axesHeight = pos(4);
        % Make sure we don't adjust away the axes entirely!
        heightAdjust = min( (axesHeight*0.9), labelHeight*axesHeight );

        % Move the axes
        if isappdata( ax, 'OriginalAxesPosition' )
            pos = getappdata( ax, 'OriginalAxesPosition' );
        else
            pos = get(ax,'Position');
            setappdata( ax, 'OriginalAxesPosition', pos );
        end
        if strcmpi( get( ax, 'XAxisLocation' ), 'Bottom' )
            % Move it up and reduce the height
            set( ax, 'Position', pos+[0 heightAdjust 0 -heightAdjust] )
        else
            % Just reduce the height
            set( ax, 'Position', pos+[0 0 0 -heightAdjust] )
        end
        set( ax, 'Units', oldUnits );
        set( ax, 'ActivePositionProperty', 'OuterPosition' );

        % Make sure we find out if axes properties are changed
        addListeners( ax );

    end % repositionAxes

%-------------------------------------------------------------------------%
    function repositionXLabel( ax )
        % Try to work out where to put the xlabel
        removeListeners( ax );
        labelHeight = getLabelHeight(ax);

        % Use the new max extent to move the xlabel. We may also need to
        % move the title
        xlab = get(ax,'XLabel');
        titleh = get( ax, 'Title' );
        set( [xlab,titleh], 'Units', 'Normalized' );
        if strcmpi( get( ax, 'XAxisLocation' ), 'Top' )
            titleExtent = get( xlab, 'Extent' );
            set( xlab, 'Position', [0.5 1+labelHeight-titleExtent(4) 0] );
            set( titleh, 'Position', [0.5 1+labelHeight 0] );
        else
            set( xlab, 'Position', [0.5 -labelHeight 0] );
            set( titleh, 'Position', [0.5 1 0] );
        end
        addListeners( ax );
    end % repositionXLabel

%-------------------------------------------------------------------------%
    function height = getLabelHeight(ax)
        height = 0;
        textLabels = findall( ax, 'Tag', 'RotatedXTickLabel' );
        if isempty(textLabels)
            return;
        end
        oldUnits = get( textLabels(1), 'Units' );
        set( textLabels, 'Units', 'Normalized' );
        for ll=1:numel(vals)
            ext = get( textLabels(ll), 'Extent' );
            if ext(4) > height
                height = ext(4);
            end
        end
        set( textLabels, 'Units', oldUnits );
    end % getLabelExtent

%-------------------------------------------------------------------------%
    function updateFont()
        % Update the rotated text fonts when the axes font changes
        properties = {
            'FontName'
            'FontSize'
            'FontAngle'
            'FontWeight'
            'FontUnits'
            };
        propertyValues = get( ax, properties );
        textLabels = findall( ax, 'Tag', 'RotatedXTickLabel' );
        set( textLabels, properties, propertyValues );
    end % updateFont

    function updateAlignment()
        textLabels = findall( ax, 'Tag', 'RotatedXTickLabel' );
        angle = get( textLabels(1), 'Rotation' );
        % Depending on the angle, we may need to change the alignment. We change
        % alignments within 5 degrees of each 90 degree orientation.
        if strcmpi( get( ax, 'XAxisLocation' ), 'Top' )
            if 0 <= angle && angle < 5
                set( textLabels, 'HorizontalAlignment', 'Center', 'VerticalAlignment', 'Bottom' );
            elseif 5 <= angle && angle < 85
                set( textLabels, 'HorizontalAlignment', 'Left', 'VerticalAlignment', 'Bottom' );
            elseif 85 <= angle && angle < 95
                set( textLabels, 'HorizontalAlignment', 'Left', 'VerticalAlignment', 'Middle' );
            elseif 95 <= angle && angle < 175
                set( textLabels, 'HorizontalAlignment', 'Left', 'VerticalAlignment', 'Top' );
            elseif 175 <= angle && angle < 185
                set( textLabels, 'HorizontalAlignment', 'Center', 'VerticalAlignment', 'Top' );
            elseif 185 <= angle && angle < 265
                set( textLabels, 'HorizontalAlignment', 'Right', 'VerticalAlignment', 'Top' );
            elseif 265 <= angle && angle < 275
                set( textLabels, 'HorizontalAlignment', 'Right', 'VerticalAlignment', 'Middle' );
            elseif 275 <= angle && angle < 355
                set( textLabels, 'HorizontalAlignment', 'Right', 'VerticalAlignment', 'Bottom' );
            else % 355-360
                set( textLabels, 'HorizontalAlignment', 'Center', 'VerticalAlignment', 'Bottom' );
            end
        else
            if 0 <= angle && angle < 5
                set( textLabels, 'HorizontalAlignment', 'Center', 'VerticalAlignment', 'Top' );
            elseif 5 <= angle && angle < 85
                set( textLabels, 'HorizontalAlignment', 'Right', 'VerticalAlignment', 'Top' );
            elseif 85 <= angle && angle < 95
                set( textLabels, 'HorizontalAlignment', 'Right', 'VerticalAlignment', 'Middle' );
            elseif 95 <= angle && angle < 175
                set( textLabels, 'HorizontalAlignment', 'Right', 'VerticalAlignment', 'Bottom' );
            elseif 175 <= angle && angle < 185
                set( textLabels, 'HorizontalAlignment', 'Center', 'VerticalAlignment', 'Bottom' );
            elseif 185 <= angle && angle < 265
                set( textLabels, 'HorizontalAlignment', 'Left', 'VerticalAlignment', 'Bottom' );
            elseif 265 <= angle && angle < 275
                set( textLabels, 'HorizontalAlignment', 'Left', 'VerticalAlignment', 'Middle' );
            elseif 275 <= angle && angle < 355
                set( textLabels, 'HorizontalAlignment', 'Left', 'VerticalAlignment', 'Top' );
            else % 355-360
                set( textLabels, 'HorizontalAlignment', 'Center', 'VerticalAlignment', 'Top' );
            end
        end
    end

%-------------------------------------------------------------------------%
    function onAxesFontChanged( ~, ~ )
        updateFont();
        repositionAxes( ax );
        repositionXLabel( ax );
    end % onAxesFontChanged

%-------------------------------------------------------------------------%
    function onAxesPositionChanged( ~, ~ )
        % We need to accept the new position, so remove the appdata before
        % redrawing
        if isappdata( ax, 'OriginalAxesPosition' )
            rmappdata( ax, 'OriginalAxesPosition' );
        end
        if isappdata( ax, 'OriginalXLabelPosition' )
            rmappdata( ax, 'OriginalXLabelPosition' );
        end
        repositionAxes( ax );
        repositionXLabel( ax );
    end % onAxesPositionChanged

%-------------------------------------------------------------------------%
    function onXAxisLocationChanged( ~, ~ )
        updateAlignment();
        repositionAxes( ax );
        repositionXLabel( ax );
    end % onXAxisLocationChanged

%-------------------------------------------------------------------------%
    function onAxesLimitsChanged( ~, ~ )
        % The limits have moved, so make sure the labels are still ok
        textLabels = findall( ax, 'Tag', 'RotatedXTickLabel' );
        xlim = get( ax, 'XLim' );
        pos = [0 getYPositionToUse( ax )];
        for tt=1:numel( textLabels )
            xval = get( textLabels(tt), 'UserData' );
            pos(1) = xval;
            % If the tick is off the edge, make it invisible
            if xval<xlim(1) || xval>xlim(2)
                set( textLabels(tt), 'Visible', 'off', 'Position', pos )
            elseif ~strcmpi( get( textLabels(tt), 'Visible' ), 'on' )
                set( textLabels(tt), 'Visible', 'on', 'Position', pos )
            else
                % Just set the position
                set( textLabels(tt), 'Position', pos );
            end
        end

        repositionXLabel( ax );
    end % onAxesPositionChanged

%-------------------------------------------------------------------------%
    function onTextLabelDeleted( ~, ~ )
        % The final text label has been deleted. This is likely from a
        % "cla" or "close" call, so we should remove all of our dirty
        % hacks.
        restoreDefaultLabels(ax);
    end

%-------------------------------------------------------------------------%
    function addListeners( ax )
        % Create listeners. We store the array of listeners in the axes to make
        % sure that they have the same life-span as the axes they are listening to.
        axh = handle( ax );
        listeners = [
            handle.listener( axh, findprop( axh, 'FontName' ), 'PropertyPostSet', @onAxesFontChanged )
            handle.listener( axh, findprop( axh, 'FontSize' ), 'PropertyPostSet', @onAxesFontChanged )
            handle.listener( axh, findprop( axh, 'FontWeight' ), 'PropertyPostSet', @onAxesFontChanged )
            handle.listener( axh, findprop( axh, 'FontAngle' ), 'PropertyPostSet', @onAxesFontChanged )
            handle.listener( axh, findprop( axh, 'FontUnits' ), 'PropertyPostSet', @onAxesFontChanged )
            handle.listener( axh, findprop( axh, 'OuterPosition' ), 'PropertyPostSet', @onAxesPositionChanged )
            handle.listener( axh, findprop( axh, 'XLim' ), 'PropertyPostSet', @onAxesLimitsChanged )
            handle.listener( axh, findprop( axh, 'YLim' ), 'PropertyPostSet', @onAxesLimitsChanged )
            handle.listener( axh, findprop( axh, 'XAxisLocation' ), 'PropertyPostSet', @onXAxisLocationChanged )
            ];
        setappdata( ax, 'RotateXLabelsListeners', listeners );
    end % addListeners

%-------------------------------------------------------------------------%
    function removeListeners( ax )
        % Rempove any property listeners whilst we are fiddling with the axes
        if isappdata( ax, 'RotateXLabelsListeners' )
            delete( getappdata( ax, 'RotateXLabelsListeners' ) );
            rmappdata( ax, 'RotateXLabelsListeners' );
        end
    end % removeListeners

%-------------------------------------------------------------------------%
    function y = getYPositionToUse( ax )
        % Use the direction and XAxisLocation properties to work out which
        % Y-value to draw the labels at.
        whichYLim = 1;
        % If YDir is reversed, switch where we position
        if strcmpi( get( ax, 'YDir' ), 'Reverse' )
            whichYLim = 3-whichYLim;
        end
        % If set to "top", then switch (again?)
        if strcmpi( get( ax, 'XAxisLocation' ), 'Top' )
            whichYLim = 3-whichYLim;
        end
        % Now get the value
        ylim = get( ax, 'YLim' );
        y = ylim(whichYLim);
    end % getYPositionToUse

end % EOF

author: yuxng


text:
function varargout = cgal_meshseg(varargin)
% CGAL_MESHSEG  Surface meshing of an isosurface from a segmentation or
% grayscale image.
%
% This function is a Matlab wrapper of the CGAL 3D Surface Mesh Generation
% for a grayscale input image.
%
% http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Surface_mesher/Chapter_main.html
%
% Note that even though CGAL internally associates x <-> rows, y <-> cols,
% this functions presents to the user the usual Matlab convention of
% x <-> cols, y <-> rows.
%
% [TRI, X] = cgal_meshseg(IM, ISOVAL)
%
%   IM is a 3D array or a SCIMAT struct with a segmentation or a grayscale
%   image.
%
%   ISOVAL is a scalar value that defines the isosurface to be meshed.
%
%   Note that if you have a binary segmentation (background=0, segmented
%   voxels=1), ISOVAL=1 will give a very tight surface connecting the
%   centres of the segmented voxels on the boundary. With small values,
%   ISOVAL=0.01, the surface will be close to the centres of the adjacent
%   background voxels. With ISOVAL=0.5, the surface will be halfway between
%   the centres of the background and segmented voxels. The latter is
%   usually the desired result.
%
%   TRI is a 3-column matrix. Each row represents the indices of the three
%   vertices that form a triangle. TRI as a whole represents the closed
%   surface.
%
%   X is a 3-column matrix. Each row represents the Cartesian coordinates
%   of a vertex on the surface, indexed by TRI values.
%
% ... = cgal_meshseg(IM, ISOVAL, MINALPHA, MAXRAD, MAXD, C, MANIFOLD)
%
%   MINALPHA, MAXRAD, MAXD are scalars that implement the three meshing
%   criteria in CGAL::Surface_mesh_default_criteria_3
%
%   http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Surface_mesher_ref/Class_Surface_mesh_default_criteria_3.html#Cross_link_anchor_1519
%
%   MINALPHA is "a lower bound on the minimum angle in degrees of the
%   surface mesh facets". By default, MINALPHA = 30.0.
%
%   MAXRAD is "an upper bound on the radius of surface Delaunay balls. A
%   surface Delaunay ball is a ball circumscribing a facet, centered on the
%   surface and empty of vertices. Such a ball exists for each facet of the
%   current surface mesh. Indeed the current surface mesh is the Delaunay
%   triangulation of the current sampling restricted to the surface which
%   is just the set of facets in the three dimensional Delaunay
%   triangulation of the sampling that have a Delaunay surface ball". By
%   default, MAXRAD is 1/2 of the minimum voxel size dimension. For
%   example, if voxels have size [0.1 0.2 0.5], then by default
%   MAXRAD=0.05.
%
%   MAXD is "an upper bound on the center-center distances of the surface
%   mesh facets. The center-center distance of a surface mesh facet is the
%   distance between the facet circumcenter and the center of its surface
%   Delaunay ball". By default, MAXD is computed the same as MAXRAD.
%
%   C is a 3-vector with the coordinates of the centre of the bounding
%   sphere used by the meshing algorithm. This is an important parameter.
%   If C is close to the surface, it can produce lots of little triangles
%   in that area. If C is outside the segmentation, the computed mesh may
%   be incomplete.
%
%   MANIFOLD is a boolean flag.
%
%     MANIFOLD=false, the mesher uses CGAL::Non_manifold_tag: "When
%     instantiated with the tag Non_manifold_tag the function template
%     make_surface_mesh does not ensure that the output mesh is a manifold
%     surface. The manifold property of output mesh may nevertheless result
%     from the choice of appropriate meshing criteria".
%
%     MANIFOLD=true, the mesher uses CGAL::Manifold_tag: "When instantiated
%     with the tag Manifold_tag the function template make_surface_mesh
%     ensures that the output mesh is a manifold surface without boundary".
%
% Important!
%
% Note that this function can produce meshes with (1) stray vertices that
% belong to no triangle, (2) triangles not oriented with respect to the
% manifold, and (3) holes in the surface. These problems can be solved with
% the following functions available in Gerardus:
%
%   [tri, x] = tri_squeeze(tri, x); % remove stray vertices
%   [x, tri] = meshcheckrepair(x, tri, 'deep'); % correct orientation
%   tri = cgal_tri_fillholes(tri, x); % fill holes in the surface
%
% See also: bwmesh, tri_squeeze, meshcheckrepair, cgal_tri_fillholes.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2013 University of Oxford
% Version: 0.1.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

error('MEX function not found')

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function lengths = meshEdgeLength(vertices, edges, faces) %#ok<INUSD>
%MESHEDGELENGTH Lengths of edges of a polygonal or polyhedral mesh
%
%   output = meshEdgeLength(V, E, F)
%
%   Example
%   meshEdgeLength
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-10-04,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% extract vertices
p1 = vertices(edges(:, 1), :);
p2 = vertices(edges(:, 2), :);

% compute euclidean distance betwenn the two vertices
lengths = sqrt(sum((p2-p1).^2, 2));

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function airy_hi_values_test ( )

%*****************************************************************************80
%
%% AIRY_HI_VALUES_TEST demonstrates the use of AIRY_HI_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'AIRY_HI_VALUES_TEST:\n' );
  fprintf ( 1, '  AIRY_HI_VALUES stores values of \n' );
  fprintf ( 1, '  the Airy Hi function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = airy_hi_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_suite = test_regression_meanf

%   Run specific demo, save values and compare the results to the expected.
%   Works for both xUnit Test Framework package by Steve Eddins and for
%   the built-in Unit Testing Framework (as of Matlab version 2013b).
%
%   See also
%     TEST_ALL, DEMO_REGRESSION_MEANF
%
% Copyright (c) 2011-2012 Ville Tolvanen
% Copyright (c) 2014 Tuomas Sivula

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  % Check which package to use
  if exist('initTestSuite', 'file')
    initTestSuite_path = which('initTestSuite');
  else
    initTestSuite_path = '';
  end
  if ~isempty(initTestSuite_path) && ...
     exist([initTestSuite_path(1:end-15) 'runtests'], 'file')
    % xUnit package
    initTestSuite;
  else
    % Built-in package
    % Use all functions except the @setup
    tests = localfunctions;
    tests = tests(~cellfun(@(x)strcmp(func2str(x),'setup'),tests));
    test_suite = functiontests(tests);
  end
end

% -------------
%     Tests
% -------------

function testRunDemo(testCase)
  % Run the correspondin demo and save the values. Note this test has to
  % be run at lest once before the other test may succeed.
  rundemo(getName())
end

function testPredictions(testCase)
  verifyVarsEqual(testCase, getName(), {'Eft'}, {@mean}, ...
    'RelTolElement', 0.1, 'RelTolRange', 0.02)
end

% ------------------------
%     Helper functions
% ------------------------

function testCase = setup
  % Helper function to suply empty array into variable testCase as an
  % argument for each test function, if using xUnit package. Not to be
  % used with built-in test framework.
  testCase = [];
end

function name = getName
  % Helperfunction that returns the name of the demo, e.g. 'binomial1'.
  name = mfilename;
  name = name(6:end);
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function genBox( obb, filename )
%GENBOX
cornerpoints = OBBrep2cornerpoints(obb);

fid = fopen(filename,'w');
for i = 1:size(cornerpoints,1)
	fprintf(fid, 'v  %f  %f  %f \n', cornerpoints(i,3), cornerpoints(i,1), cornerpoints(i,2));
end
fprintf(fid, 'f  %d  %d  %d %d \n', 4, 8, 7, 3);
fprintf(fid, 'f  %d  %d  %d %d \n', 8, 6, 5, 7);
fprintf(fid, 'f  %d  %d  %d %d \n', 6, 2, 1, 5);
fprintf(fid, 'f  %d  %d  %d %d \n', 2, 4, 3, 1);
fprintf(fid, 'f  %d  %d  %d %d \n', 4, 2, 6, 8);
fprintf(fid, 'f  %d  %d  %d %d \n', 7, 5, 1, 3);
fclose(fid);

end

repo: GRAINS


text:
function [llh] = tapas_sem_prosa_llh(y, u, theta, ptheta)
%% Computes the likelihood of the data.
%
% Input
%
%   y -- Observed behavioral data. A structure with fields 't' times and 'a'
%       action
%   u -- Experimental input. A structure with fields: 'tt' trial type, either
%       prosaccade or antisaccade
%   theta -- Model parameters
%   ptheta -- Priors
%
% Output
%
%   llh -- Log likelilhood
%
% aponteeduardo@gmail.com
% copyright (C) 2015
%

% Compute the likelihood of antisaccades and prosaccades

llh = zeros(1, numel(theta));

it = y.i;

method = ptheta.method;
ptrans = ptheta.ptrans;

for i = 1:numel(theta)
    llh(i) = sum(tapas_sem_prosa_cllh(y.t(~it), y.a(~it), u.tt(~it), ...
        ptrans(theta{i}), method, 1), 1);
end

end

author: translationalneuromodeling

repo: tapas


text:
function value = p38_autonomous ( )

%*****************************************************************************80
%
%% P38_AUTONOMOUS reports whether problem P38 is autonomous.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer VALUE, is 1 if the equation is autonomous.
%
  value = 1;

  return
end


text:
function [bodyInfo] = getBodyInfoStructFromOrbit(inputOrbit)
%getBodyInfoStructFromOrbit Summary of this function goes here
%   Detailed explanation goes here

%     bodyInfo = struct();
    bodyInfo = KSPTOT_BodyInfo();
    bodyInfo.epoch = inputOrbit(7);
    bodyInfo.sma = inputOrbit(1);
    bodyInfo.ecc = inputOrbit(2);
    bodyInfo.inc = rad2deg(inputOrbit(3));
    bodyInfo.raan = rad2deg(inputOrbit(4));
    bodyInfo.arg = rad2deg(inputOrbit(5));
    bodyInfo.mean = rad2deg(inputOrbit(6));
    bodyInfo.id = rand();
    bodyInfo.propTypeEnum = BodyPropagationTypeEnum.TwoBody;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function xi = calcProjection( X, Y )
%CALCPROJECTION Projection into the tangent space.
%   XI = CALCPROJECTION( X, Y ) computes the orth. projection of the Tucker
%   tensor (ttensor) Y into the tangent space at X (ttensor).
%   Of the resulting tangent tensor XI, only the variations are stored.
%   Hence, XI is a struct with the fields
%
%       XI.Y_tilde  (var. in the core)
%       XI.U1_tilde (var. in the factors...
%       XI.U2_tilde
%       XI.U3_tilde  ...)
%
%   For efficiency, parts of the calculation are performed by the mex
%   routine calcProjection_mex
%
%   See also calcGradient.
%

%   GeomCG Tensor Completion. Copyright 2013 by
%   Michael Steinlechner
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt

    [temp,temp1,temp2,temp3] = calcProjection_mex( Y.subs', Y.vals, X.U{1}', X.U{2}', X.U{3}' );

    [n1,k1] = size(X.U{1});
    [n2,k2] = size(X.U{2});
    [n3,k3] = size(X.U{3});

    xi.Y_tilde = tensor(reshape( temp, [k1 k2 k3]), [k1, k2, k3]);

    Y23 = tensor(reshape(temp1, [n1 k2 k3]));
    Y13 = tensor(reshape(temp2, [k1 n2 k3]));
    Y12 = tensor(reshape(temp3, [k1 k2 n3]));

    S1_inv = pinv( double( tenmat( X.core, 1 ) ));
    S2_inv = pinv( double( tenmat( X.core, 2 ) ));
    S3_inv = pinv( double( tenmat( X.core, 3 ) ));

    U1_tilde = double( tenmat( Y23, 1 )) * S1_inv;
    U2_tilde = double( tenmat( Y13, 2 )) * S2_inv;
    U3_tilde = double( tenmat( Y12, 3 )) * S3_inv;

    xi.U1_tilde = U1_tilde - X.U{1} * ( X.U{1}' * U1_tilde );
    xi.U2_tilde = U2_tilde - X.U{2} * ( X.U{2}' * U2_tilde );
    xi.U3_tilde = U3_tilde - X.U{3} * ( X.U{3}' * U3_tilde );

end

author: andrewssobral


text:
function [xY,word,NI] = spm_voice_get_xY(PATH)
% Create word arrays from sound file exemplars
% FORMAT [xY,word,NI] = spm_voice_get_xY(PATH)
%
% PATH      -  directory containing sound files of exemplar words
%
% xY(nw,ns) -  structure array for ns samples of nw words
% word(nw)  -  cell array of word names
% NI(nw,ns) -  numeric array of number of minima
%
%  This routine uses a library of sound files, each containing 32 words
%  spoken with varying prosody. The name of the sound file labels the word
%  in question. These exemplars are then transformed (using a series of
%  discrete cosine and Hilbert transforms) into a set of parameters, which
%  summarise the lexical content and prosody. The inverse transform
%  generates  timeseries that can be played to articulate a word. The
%  transform operates on a word structure xY to create lexical and prosody
%  parameters (Q and P respectively).
%__________________________________________________________________________
% Copyright (C) 2019 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_voice_get_xY.m 7766 2020-01-05 21:37:39Z karl $

% get corpus
%==========================================================================
rng('default')

% get the list of words and sampling frequency (FS)
%--------------------------------------------------------------------------
spm_figure('GetWin','voice'); clf
cd(PATH)
wfile    = dir('*.wav');
try
    audioread(wfile(1).name,[1,1]);
    read = @audioread;
catch
    read = @wavread;
end
[Y,FS]   = read(wfile(1).name,[1,1]);

% assemble cell array of word structures for subsequent characterisation
%==========================================================================
nw    = numel(wfile);                               % number of words
ns    = 32;                                         % number of samples
nj    = 16;                                         % number of jitters
sj    = FS/64;                                      % s.d. of jitter
for w = 1:nw

    % get lexicon name and create structure
    %----------------------------------------------------------------------
    wname   = wfile(w).name;                        % name of word
    [d,str] = fileparts(wname);
    word{w} = str; disp(str)

    % get the midpoint of words from the (maxima) of successive exemplars
    %----------------------------------------------------------------------
    G     = spm_voice_check(read(wname),FS,1/4);
    I     = find((diff(G(1:end - 1)) > 0) & (diff(G(2:end)) < 0));
    [i,j] = sort(G(I),'descend');
    I     = sort(I(j(1:ns)));

    % and plot
    %----------------------------------------------------------------------
    subplot(2,1,1), plot(G), hold on, plot(I,G(I),'ro'), hold off
    title(sprintf('Power peaks - %s',word{w}),'FontSize',16)
    xlabel('frequency (hertz)'), ylabel('power'), box off
    drawnow

    for s = 1:ns

        % retrieve (one second) epoch around midpoint and transform
        %------------------------------------------------------------------
        i     = round([-1/2 1/2]*FS + I(s));
        if i(1) < 1
            i = i - i(1) + 1;
        end
        Y     = read(wname,i);
        i     = spm_voice_onsets(Y,FS);
        i     = i{end};
        ni    = numel(Y);
        for j = 1:nj

            % jitter interval
            %--------------------------------------------------------------
            k       = (s - 1)*nj + j;
            ik      = fix((i(1) + sj*randn):(i(end) + sj*randn));
            ik      = ik(ik < ni & ik > 1);
            xY(w,k) = spm_voice_ff(Y(ik),FS);

            % record interval
            %--------------------------------------------------------------
            xY(w,k).i(1) = ik(1)/FS   - 1/2;
            xY(w,k).i(2) = ik(end)/FS - 1/2;
        end

        % apply inverse transform and play, if requested
        %------------------------------------------------------------------
        spm_voice_iff(xY(w,k));

        % record number of spectral minima intervals
        %------------------------------------------------------------------
        [Y,IS]  = spm_voice_get_next(Y);
        j       = fix((0:FS + FS/4) + IS);
        j       = j(j < ni & j > 1);
        NI(w,s) = numel(spm_voice_onsets(Y(j),FS));

    end

end

author: spm


text:
function clearDoseMeshes(varargin)
%function clearDoseMeshes()
%This function clears all the surface meshes for passed dose uidC.
%Note that the contour levels which are cleared depend on the current
%selections.
%
%APA, 06/21/07
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC
indexS = planC{end};

if isempty(varargin)
    doseNumToClear = 1:length(planC{indexS.dose});
else
    uidC = varargin{1};
    doseNumToClear = getAssociatedDose(uidC);
end

%Set Matlab path to directory containing the library
currDir = cd;

if ispc
    meshDir = fileparts(which('libMeshContour.dll'));
    cd(meshDir);
    loadlibrary('libMeshContour','MeshContour.h');
elseif isunix
    meshDir = fileparts(which('libMeshContour.so'));
    cd(meshDir);
    loadlibrary('libMeshContour.so','MeshContour.h');
end
%Clear existing dose meshes
waitbarH = waitbar(0,'Clearing surface meshes for dose...');
try
    contourLevels = getIsoDoseLevels;
    for doseNum = doseNumToClear
        doseUID = planC{indexS.dose}(doseNum).doseUID;
        for level = 1:length(contourLevels)
            calllib('libMeshContour','clear',[doseUID,'_',num2str(contourLevels(level))])
        end
        waitbar(doseNum/length(planC{indexS.dose}),waitbarH)
    end
end
close(waitbarH)
%switch back the current irectory
cd(currDir)

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function grad = fif2grad(filename);

% FIF2GRAD constructs a gradiometer definition from a Neuromag *.fif file
% The resulting gradiometer definition can be used by Fieldtrip for forward
% and inverse computations.
%
% Use as
%   grad = fif2grad(filename);
%
% See also READ_HEADER, CTF2GRAD, BTI2GRAD

% Copyright (C) 2004, Joachim Gross
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: fif2grad.m 2885 2011-02-16 09:41:58Z roboos $

% this try-catch construct ensures that missing gradiometer information is
% handeled in a "graceful" way
grad = [];
try
   megmodel('head',[0 0 0],filename);
   [n,s,t]=chaninfo;
   [TY,NA]=chaninfo('type');
   nCoils=sum(TY+1);      % number of coils
   nSensors=length(TY);   % number of sensors
   grad.pnt=zeros(nCoils,3);
   grad.ori=zeros(nCoils,3);
   grad.tra=zeros(nSensors,nCoils);
   grad.unit='cm';
   % define coils
   kCoil=1;
   for k=1:nSensors,
      if (TY(k)==0),  %magnetometer
         grad.pnt(kCoil,:)=100*(t{k}(1:3,4));
         grad.ori(kCoil,:)=t{k}(1:3,3);
         grad.tra(k,kCoil)=1;
         kCoil=kCoil+1;
         grad.label{k}=deblank(s(k,:));
      elseif (TY(k) == 1), %planar gradiometer
         grad.pnt(kCoil,:)=100*(t{k}(1:3,4)-0.008*t{k}(1:3,1)); % multiply with 100 to get cm
         grad.ori(kCoil,:)=t{k}(1:3,3);
         grad.tra(k,kCoil)= -1;
         kCoil=kCoil+1;
         grad.pnt(kCoil,:)=100*(t{k}(1:3,4)+0.008*t{k}(1:3,1));
         grad.ori(kCoil,:)=t{k}(1:3,3);
         grad.tra(k,kCoil)= 1;
         kCoil=kCoil+1;
         grad.label{k}=deblank(s(k,:));
      else
         error('unknown sensor type');
      end
   end
catch
   warning(['gradiometer information could not be extracted from file']);
   warning(['returning an empty grad structure']);
   return;
end

author: PatternRecognition

repo: OpenBMI


text:
function cell2file(file,str,flag)
% write cellstring to file

if nargin == 2, flag = 'w'; end
fid = efopen(file,flag);

write_cell(fid,str);
fclose(fid);

author: mtex-toolbox

repo: mtex

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


text:
function monom = recovermonoms(newton_m,x);

if isempty(newton_m)
    monom = 1;
else
    error('Report this bug (call to recovermonoms with double)');
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [Population,Dec,Mask,FrontNo,CrowdDis] = SPEA2_EnvironmentalSelection(Population,Dec,Mask,N)
% The environmental selection of MSKEA

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Lei Chen

    %% Delete duplicated solutions
    [~,uni] = unique(Population.objs,'rows');
    Population = Population(uni);
    Dec        = Dec(uni,:);
    Mask       = Mask(uni,:);
    N          = min(N,length(Population));
    %% Calculate the fitness of each solution

    [FrontNo,MaxFNo] = NDSort(Population.objs,N);
    Next = false(1,length(FrontNo));
    Next(FrontNo<MaxFNo) = true;

    PopObj = Population.objs;
    fmax   = max(PopObj(FrontNo==1,:),[],1);
    fmin   = min(PopObj(FrontNo==1,:),[],1);
    PopObj = (PopObj-repmat(fmin,size(PopObj,1),1))./repmat(fmax-fmin,size(PopObj,1),1);

    %% Environmental selection
    Last = find(FrontNo==MaxFNo);
    del  = Truncation(PopObj(Last,:),length(Last)-N+sum(Next));
    Next(Last(~del)) = true;
    % Population for next generation
    Population = Population(Next);
    Dec        = Dec(Next,:);
    Mask       = Mask(Next,:);
    FrontNo    = FrontNo(Next);
    CrowdDis = CrowdingDistance(Population.objs,FrontNo);
end

function Del = Truncation(PopObj,K)
% Select part of the solutions by truncation

    %% Truncation
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Del = false(1,size(PopObj,1));
    while sum(Del) < K
        Remain   = find(~Del);
        Temp     = sort(Distance(Remain,Remain),2);
        [~,Rank] = sortrows(Temp);
        Del(Remain(Rank(1))) = true;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [] = bihist(x1,x2,n)
% BIHIST(X1,X2,N) open a new figure an plot a back to back histogram
% X1: variable 1 (right side, colour blue)
% X2: variable 2 (right side, colour green)
% n : number of classes for each histogramm

% make two dummy-histogram
[temp, d1] = hist(x1,n);  %#ok<ASGLU>
[temp, d2] = hist(x2,n);  %#ok<ASGLU>
clear temp
% calculate classes
d = min([d1 d2]):((max([d1 d2])-min([d1 d2]))/(n-1)) : max([d1 d2]);
% create the two histograms
templ = hist(x1,d);
tempr = hist(x2,d);
% scaling the histograms
templ=templ*-1;
% create bar plot
left=barh(templ,'style','hist');
hold on
right = barh(tempr,'style','hist');
plot([0 0], [1 n],'k')
% change color of left histogram
set(left,'FaceColor',[0 0 1]);
% change color of right histogram
set(right,'FaceColor',[0 1 0]);
% distance factor to frame
dt=1.2;
xlim([-max(abs([templ tempr]))*dt max(abs([templ tempr]))*dt])
% mark the y-axis
yl = get(gca,'YTick');
b = (max(d)-min(d))/(n-1);
a = max(d)-b*n;
yl2 = yl*b+a;
ds = cell(size(yl2));
for k=1:length(yl2)
    ds{k} = num2str(yl2(k),3);
end
set(gca,'YTickLabel',ds)

% rewrite the x-axis with positiv values
xl = get(gca,'XTickLabel');
xl = regexprep(cellstr(xl),'-','');
set(gca,'XTickLabel',regexprep(cellstr(xl),'-',''))

author: Sable

repo: mcbench-benchmarks


text:
function lambda = rodman_eigenvalues ( n, alpha )

%*****************************************************************************80
%
%% RODMAN_EIGENVALUES returns the eigenvalues of the RODMAN matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real ALPHA, the parameter.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  lambda(1:n-1,1) = 1.0 - alpha;

  lambda(n,1) = 1.0 + alpha * ( n - 1 );

  return
end


text:
function [varargout]=cunique(varargin)

% function [A_uni,ind1,ind2,Ac]=cunique(A)
%-------------------------------------------------------------------------
%This function is similar to MATLAB's unique function. There are three
%differences: 1) An additional 4th optional output is available providing
%the count, or number of occurances, for each element in the input array.
%2) The 2nd output mathes the size of the first input, 3) The 3rd output
%is reshaped to be the size of the input variable.
%
% See also: unique
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% Change log:
% 2018/03/21 Created
% 2019/07/02 Adding varargin handling (pass unique options) including
% 'rows' support.
% 2020/01/10 Fixed bug in accumarray array size specification
%-------------------------------------------------------------------------

%%

A=varargin{1};

[A_uni,ind1,ind2]=unique(varargin{:});

if any(strcmp(varargin,'rows'))
    typeOpt=2;
else
    typeOpt=1;
end

switch typeOpt
    case 1
        varargout{1}=A_uni;
        varargout{2}=reshape(ind1,size(A_uni));
        varargout{3}=reshape(ind2,size(A));
    case 2
        varargout{1}=A_uni;
        varargout{2}=reshape(ind1,[size(A_uni,1) 1]);
        varargout{3}=reshape(ind2,[size(A,1) 1]);
end

if nargout==4
    [subInd] = ind2subn(size(A_uni),ind2);
%     Ac=accumarray(subInd,ones(numel(ind2),1),[size(A_uni,1),1]);
    Ac=accumarray(subInd,ones(numel(ind2),1),size(A_uni));
    switch typeOpt
        case 1
            Ac=reshape(Ac(ind2),size(A));
        case 2
            Ac=reshape(Ac(ind2),[size(A,1) 1]);
    end
    varargout{4}=Ac;
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [pf_, opf_] = t_node_test(quiet)
%T_NODE_TEST  Tests for network model with multipe node-creating elements.

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if nargin < 1
    quiet = 0;
end

cases = {'t_case9_gizmo', 'case4gs', 'case4_dist', 'case9', 'case14', 'case57', 'case300'};

t_begin(12*length(cases), quiet);

define_constants;
if quiet
    verbose = 0;
else
    verbose = 1;
end

mpopt = mpoption('out.all', 0, 'verbose', 0, 'pf.tol', 1e-10);
mpopt = mpoption(mpopt, 'opf.ignore_angle_lim', 1);
mpopt0 = mpopt;
mpopt0.exp.use_legacy_core = 1;
mpopt.exp.mpx = mp.xt_node_test();

for k = 1:length(cases)
    t = sprintf('PF - %s - ', cases{k});
    mpc = loadcase(cases{k});
    if strcmp(cases{k}, 't_case9_gizmo')
        mpc.bus(2, BS) = 1;
    end
    if ~isfield(mpc, 'gencost')
        ng = size(mpc.gen, 1);
        mpc.gencost = ones(ng, 1) * [2 0 0 2 1 0];
        if strcmp(cases{k}, 'case4gs')
            mpc.gen(:, PMAX) = 320;
            mpc.gen(:, QMAX) = 200;
        end
    end

    r = runpf(mpc, mpopt0);
    evm = r.bus(:, VM);
    eva = r.bus(:, VA);
    epg = r.gen(:, PG);
    eqg = r.gen(:, QG);
    t_ok(r.success, [t 'success 1']);

    pf = run_pf(mpc, mpopt);
    r2 = pf.dmc.export(pf.dm, pf.dm.source);
    va = r2.bus(:, VA);
    vm = r2.bus(:, VM);
    pg = r2.gen(:, PG);
    qg = r2.gen(:, QG);
    t_ok(pf.success, [t 'success 2']);
    t_is(va, eva, 9, [t 'va']);
    t_is(vm, evm, 9, [t 'vm']);
    t_is(pg, epg, 9, [t 'pg']);
    t_is(qg, eqg, 9, [t 'qg']);

    t = sprintf('OPF - %s - ', cases{k});
    r = runopf(mpc, mpopt0);
    t_ok(r.success, [t 'success 1']);
    evm = r.bus(:, VM);
    eva = r.bus(:, VA);
    epg = r.gen(:, PG);
    eqg = r.gen(:, QG);

    opf = run_opf(mpc, mpopt);
    r2 = opf.dmc.export(opf.dm, opf.dm.source);
    va = r2.bus(:, VA);
    vm = r2.bus(:, VM);
    pg = r2.gen(:, PG);
    qg = r2.gen(:, QG);
    t_ok(opf.success, [t 'success 2']);
    t_is(va, eva, 9, [t 'va']);
    t_is(vm, evm, 9, [t 'vm']);
    t_is(pg, epg, 9, [t 'pg']);
    t_is(qg, eqg, 9, [t 'qg']);
end

t_end;

if nargout
    pf_ = pf;
    if nargout > 1
        opf_ = opf;
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function ps_parms_default()
%PS_PARMS_DEFAULT set parms to default value if not already set
%
%   Andy Hooper, June 2006
%
%   ======================================================================
%   04/2008 AH: SB default processing corrected
%   06/2009 AH: scla_deramp added
%   08/2009 AH: unwrap_gold_alpha added
%   09/2009 AH: merge_resample_size added
%   09/2009 AH: select_method added
%   09/2009 AH: density_rand added
%   11/2009 AH: pixel_aspect_ratio removed
%   01/2010 AH: add check for write permission before saving
%   02/2010 AH: weed_alpha replaced by weed_time_win
%   02/2010 AH: weed_max_noise added
%   02/2010 AH: plot_pixel_size replaced by plot_scatterer_size (in m)
%   02/2010 AH: plot_dem_posting added
%   02/2010 AH: ref_centre_lonlat and ref_radius added
%   02/2010 AH: unwrap_ifg_index replaced by drop_ifg_index
%   02/2010 AH: weed_time_win and unwrap_time_win default changed to 730
%   03/2010 AH: add logging
%   09/2010 MA: oversampling factor
%   09/2010 AH: plot_pixels_scatterer added
%   11/2010 AH: (sb_)recalc_index replaced by (sb_)scla_drop_index
%   06/2011 AH: add weed_neighbours
%   12/2012 AH: add insar_processor
%   06/2013 DB: add parameters to include tropopshere in scla estimation
%   01/2014 DB: Add wavelength and heading to parameters list. This in case
%               not all stamps steps are ran.
%   01/2014 DB: Add a multi-core option for step 1-5
%   02/2014 AH: add unwrap_hold_good_values
%   08/2014 DB: explicit search if there is not a processor.txt file
%   05/2015 AH: add gamma_max_iterations
%   09/2015 AH: change default weed_neighbours = 'n'
%   09/2015 AH: change default max_topo_err = 20
%   10/2016 DB: Change other processor type to a warning not an error.
%   ======================================================================

parmfile='parms.mat';
parent_flag=0;

if exist('./parms.mat','file')
    parms=load(parmfile);
elseif exist('../parms.mat','file')
    parmfile='../parms.mat';
    parms=load(parmfile);
    parent_flag=1;
else
    parms=struct('Created',date);
    parms.small_baseline_flag='n';
end

parmfields_before=fieldnames(parms);
num_fields=size(parmfields_before,1);

if ~isfield(parms,'max_topo_err')
    parms.max_topo_err=20;
end

if ~isfield(parms,'quick_est_gamma_flag')
    parms.quick_est_gamma_flag='y';
end

if ~isfield(parms,'select_reest_gamma_flag')
    parms.select_reest_gamma_flag='y';
end

if ~isfield(parms,'filter_grid_size')
    parms.filter_grid_size=50;
end

if ~isfield(parms,'filter_weighting')
    parms.filter_weighting='P-square'; % filter weighting strategy
end

if ~isfield(parms,'gamma_change_convergence')
    parms.gamma_change_convergence=0.005; % change in change in gamma that signals convergence
end

if ~isfield(parms,'gamma_max_iterations')
    parms.gamma_max_iterations=3; % maximum number of iterations for gamma estimation
end

if ~isfield(parms,'slc_osf')
    parms.slc_osf=1;  % [MA] SLC oversampling factor
end

if ~isfield(parms,'clap_win')
    parms.clap_win=32;
end

if ~isfield(parms,'clap_low_pass_wavelength')
    parms.clap_low_pass_wavelength=800;
end

if ~isfield(parms,'clap_alpha')
    parms.clap_alpha=1;
end

if ~isfield(parms,'clap_beta')
    parms.clap_beta=0.3;
end

if ~isfield(parms,'select_method')
    parms.select_method='DENSITY'; %'DENSITY' or 'PERCENT'
end

if ~isfield(parms,'density_rand')
    if strcmpi(parms.small_baseline_flag,'y')
        parms.density_rand=2;  % Random-phase pixels per km2 (before weeding)
    else
        parms.density_rand=20; % Random-phase pixels per km2 (before weeding)
    end
end

if ~isfield(parms,'percent_rand')
    if strcmpi(parms.small_baseline_flag,'y')
        parms.percent_rand=1; % Percent random-phase pixels (before weeding)
    else
        parms.percent_rand=20;% Percent random-phase pixels (before weeding)
    end
end

if ~isfield(parms,'gamma_stdev_reject')
    parms.gamma_stdev_reject=0;
end

if isfield(parms,'weed_alpha')
    parms=rmfield(parms,'weed_alpha');
end

if ~isfield(parms,'weed_time_win')
    parms.weed_time_win=730;    % weeding smoothing window alpha (days)
end

if ~isfield(parms,'weed_max_noise')
    parms.weed_max_noise=inf;   % maximum arc noise in any ifg (rad)
end

if ~isfield(parms,'weed_standard_dev')
    if strcmpi(parms.small_baseline_flag,'y')
        parms.weed_standard_dev=inf;
    else
        parms.weed_standard_dev=1.0;
    end
end

if ~isfield(parms,'weed_zero_elevation')
    parms.weed_zero_elevation='n';
end

if ~isfield(parms,'weed_neighbours')
    parms.weed_neighbours='n';
end

if ~isfield(parms,'unwrap_method')
    if strcmpi(parms.small_baseline_flag,'y')
        parms.unwrap_method='3D_QUICK';
    else
        parms.unwrap_method='3D';
    end
end

if ~isfield(parms,'unwrap_patch_phase')
    parms.unwrap_patch_phase='n';
end

if isfield(parms,'unwrap_ifg_index')
    try
       ps=load('ps2.mat');
    catch
       try
          ps=load('ps1.mat');
       catch
       end
    end
    if exist('ps','var') & ~strcmpi(parms.unwrap_ifg_index,'all')
       parms.drop_ifg_index=setdiff([1:ps.n_ifg],parms.unwrap_ifg_index);
    end
    parms=rmfield(parms,'unwrap_ifg_index');
    num_fields=0;
end

if ~isfield(parms,'drop_ifg_index')
    parms.drop_ifg_index=[];
end

if ~isfield(parms,'unwrap_la_error_flag')
    parms.unwrap_la_error_flag='y';
end

if ~isfield(parms,'unwrap_spatial_cost_func_flag')
    parms.unwrap_spatial_cost_func_flag='n';
end

if ~isfield(parms,'unwrap_prefilter_flag')
    parms.unwrap_prefilter_flag='y';
end

if ~isfield(parms,'unwrap_grid_size')
    parms.unwrap_grid_size=200; % prefilter grid size
end

if ~isfield(parms,'unwrap_gold_n_win')
    parms.unwrap_gold_n_win=32; % prefilter goldstein filtering window size
end

if ~isfield(parms,'unwrap_alpha')
    parms.unwrap_alpha=8;    % unwrapping smoothing window alpha
end

if ~isfield(parms,'unwrap_time_win')
    parms.unwrap_time_win=730;    % unwrapping smoothing window alpha
end

if ~isfield(parms,'unwrap_gold_alpha')
    parms.unwrap_gold_alpha=0.8;  % unwrapping goldstein filter alpha
end

if ~isfield(parms,'unwrap_hold_good_values')
    parms.unwrap_hold_good_values='n';       % fix unwrapped values judged to be good
end

if isfield(parms,'recalc_index')
    try
       ps=load('ps2.mat');
    catch
       try
          ps=load('ps1.mat');
       catch
       end
    end
    if exist('ps','var') & ~strcmpi(parms.recalc_index,'all')
        if strcmpi(parms.small_baseline_flag,'y')
           parms.scla_drop_index=setdiff([1:ps.n_image],parms.recalc_index);
        else
           parms.scla_drop_index=setdiff([1:ps.n_ifg],parms.recalc_index);
        end
    end
    parms=rmfield(parms,'recalc_index');
    if isfield(parms,'sb_recalc_index')
        if exist('ps','var') & ~strcmpi(parms.sb_recalc_index,'all')
            parms.sb_scla_drop_index=setdiff([1:ps.n_ifg],parms.sb_recalc_index);
        end
        parms=rmfield(parms,'sb_recalc_index');
    end
    num_fields=0;
end

if ~isfield(parms,'scla_drop_index')
    parms.scla_drop_index=[];
end

if ~isfield(parms,'scn_wavelength')
    parms.scn_wavelength=100; % spatially correlated noise wavelength
end

if ~isfield(parms,'scn_time_win')
    parms.scn_time_win=365; % 1 year time window
end

if ~isfield(parms,'scn_deramp_ifg')
    parms.scn_deramp_ifg=[];  % deramp these ifgs and add to estimate of scn
end

if ~isfield(parms,'scn_kriging_flag')
    parms.scn_kriging_flag='n';  % use kriging, 'y' or 'n'
end

if ~isfield(parms,'ref_lon')
    parms.ref_lon=[-inf,inf];  % low and high longitude for ref ps
end

if ~isfield(parms,'ref_lat')
    parms.ref_lat=[-inf,inf];  % low and high latitude for ref ps
end

if ~isfield(parms,'ref_centre_lonlat')
    parms.ref_centre_lonlat=[0,0];  % centre lon/lat for ref ps
end

if ~isfield(parms,'ref_radius')
    parms.ref_radius=inf;  % radius from centre for ref ps
end

if ~isfield(parms,'ref_velocity')
    parms.ref_velocity=0;  % velocity of reference point if known
end

if ~isfield(parms,'n_cores')
    parms.n_cores=1;  % n_cores for the muti-core option (step 1-5)
end

if ~isfield(parms,'plot_dem_posting')
    parms.plot_dem_posting=90;
end

if isfield(parms,'plot_pixel_size')
    parms.plot_scatterer_size=parms.plot_pixel_size*25;
    num_fields=0;
    parms=rmfield(parms,'plot_pixel_size');
end

if ~isfield(parms,'plot_scatterer_size')
    parms.plot_scatterer_size=120;
end

if ~isfield(parms,'plot_pixels_scatterer')
    parms.plot_pixels_scatterer=3;
end

if ~isfield(parms,'plot_color_scheme')
    parms.plot_color_scheme='inflation';
end

if isfield(parms,'pixel_aspect_ratio')
    parms=rmfield(parms,'pixel_aspect_ratio');
end

if ~isfield(parms,'shade_rel_angle')
    parms.shade_rel_angle=[90,45]; % look angle for dem shaded relief
end

if ~isfield(parms,'lonlat_offset')
    parms.lonlat_offset=[0,0]; % offset of PS in degrees from dem
end

if ~isfield(parms,'merge_resample_size')
    if strcmpi(parms.small_baseline_flag,'y')
        parms.merge_resample_size=100; % grid size (in m) to resample to during merge of patches
    else
        parms.merge_resample_size=0;   % (0=no resampling)
    end
end

if ~isfield(parms,'merge_standard_dev')
    parms.merge_standard_dev=inf;
end

if ~isfield(parms,'scla_method')
    parms.scla_method='L2'; % method for estmating SCLA, L1- or L2-norm
end

if ~isfield(parms,'scla_deramp')
    parms.scla_deramp='n'; % estimate an orbital ramp before SCLA
end

lambdaname=['lambda.1.in'];         % wavelength
if ~isfield(parms,'lambda')
    if ~exist(lambdaname,'file')
        lambdaname= ['../',lambdaname];
    end
    if ~exist(lambdaname,'file')
        lambdaname= ['../',lambdaname];
    end
    if ~exist(lambdaname,'file')
        parms.lambda=NaN; % Add wavelength
    else
        lambda=load(lambdaname);
        parms.lambda=lambda; % Add wavelength
    end
end

headingname=['heading.1.in'];       % satellite heading
if ~isfield(parms,'heading')
    if ~exist(headingname,'file')
        headingname= ['../',headingname];
    end
    if ~exist(headingname,'file')
        headingname= ['../',headingname];
    end
    if ~exist(headingname,'file')
        parms.heading=NaN; % Add heading
    else
        heading=load(headingname);
        parms.heading=heading; % Add heading
    end
end

if ~isfield(parms,'scla_deramp')
    parms.scla_deramp='n'; % estimate an orbital ramp before SCLA
end

if ~isfield(parms,'sb_scla_drop_index')
    if strcmpi(parms.small_baseline_flag,'y')
        parms.sb_scla_drop_index=[];
    end
end

if ~isfield(parms,'insar_processor')
    processor_file = 'processor.txt';
    if exist(processor_file,'file')~=2
       if exist(['..' filesep processor_file],'file')==2
           processor_file = ['..' filesep processor_file];
           if exist(['..' filesep processor_file],'file')==2
               processor_file = ['..' filesep processor_file];
           end
       end
    end

   if exist(processor_file,'file')~=2
       parms.insar_processor='doris'; %
   else
       processor = fileread(processor_file);
       processor = strtrim(processor);
       parms.insar_processor=processor; %

       if ~strcmpi(processor,'gamma') & ~strcmpi(processor,'doris')
           fprintf('WARNING: This processor is not supported (doris and gamma)')
       end
   end
end

if ~isfield(parms,'subtr_tropo')
    parms.subtr_tropo='n';         % remove tropospheric estimate
end

if ~isfield(parms,'tropo_method')
    parms.tropo_method='a_l';       % method for tropopsheric estimate
end

parmfields=fieldnames(parms);
if size(parmfields,1)~=num_fields
    try
        save(parmfile,'-struct','parms')
        for i=1:size(parmfields,1)
            if isempty(strmatch(parmfields{i},parmfields_before))
               parmname=parmfields{i};
               value=getfield(parms,parmname);
               if isempty(value)
                   value='[]';
               end
               if isnumeric(value)
                   logit([parmname,' = ',num2str(value)],0,parent_flag);
               else
                   logit([parmname,' = ',value],0,parent_flag);
               end
            end
        end

    catch
        fprintf('Warning: missing parameters could not be updated (no write access)\n')
    end
end

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
lw=2;
fz=20;

data_fb_real =load('longer_fb_real.log');
data_fac_real=load('longer_fac_real.log');
data_ola_real=load('longer_ola_real.log');

% Columns in data, fb : a M L W gl time
% Columns in data, fac: a M L W time
Ls=data_fac_real(:,3);
t_fb_real =data_fb_real(:,6);
t_fac_real=data_fac_real(:,5);
t_ola_real=data_ola_real(:,7);

if 0
  % Color legend
  l1='b';
  l2='b--';
  l3='r';
  l4='r--';
else
  % bw legend
  l1='b';
  l2='b--';
  l3='b-.';
  l4='b:';
end;

figure(1);

a=data_fac_real(1,1);
M=data_fac_real(1,2);
L=data_fac_real(:,3);

plot(Ls,t_fb_real,l1,...
     Ls,t_fac_real,l2,...
     Ls,t_ola_real,l3,'LineWidth',lw);
set(gca,'Fontsize',fz);

h=legend('Portnoff, real','Fac, real','Fac-OLA, real',...
       'Location','NorthWest');
% Grow the box a little, otherwise the export to .eps is messed up.
q=get(h,'Position');
set(h,'Position',[q(1)*0.90 q(2)*.95 q(3)*1.8 q(4)]);

xlabel('Signal length / samples','Fontsize',fz);
ylabel('Running time / seconds','Fontsize',fz);

print -deps plot_longer_1.eps

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the following exercises:
%   A) Exercise A: Use the Runge-Kutta 4 algorithm to integrate the
%   differential equation:
%            dy/dt = 2*t
%      Integrate from t0=0, to tfinal = 10 s.
%      Compare the solution with the algebraic integration.
%   B) Exercise B: Use the Runge-Kutta 4 algorithm to integrate the
%   second order differential equation:
%            d2y/dt^2+5*dy/dt + 4*y(t) = 0
%
%   C) Exercise C: Use the Runge-Kutta 4 algorithm to simulate the movement
%   of a 1 dof robot arm with friction under the effect of gravity and with
%   zero torque applied.
%
%   D) Exercise D: Use the Runge-Kutta 4 algorithm to simulate the movement
%   of a 2 dof robot arm with friction under the effect of gravity and with
%   zero torques applied.
%
% Help: Function prototype
% [y, t] = runge_kutta(f, y0, [t0 tfinal], timestep)
% where f is the function being integrated as dy/dt = f(t, y).
% y0 are the initial conditions
% t0: initial time
% tfinal: final time.
% h: time step for the calculations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function runge_kutta_exercises()
close all;

%uncomment to execute each of the exercises
% exerciseA()%
exerciseB()
%exerciseC()
%exerciseD()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Integrate a simple time function. dy/dt = 2*ts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function exerciseA()
t0 = 0;
tfinal = 10;
[t, y] = runge_kutta(@line, 0, [t0 tfinal], 1);

%Compare with the integral of 2*t
error = y(:)-t(:).^2;
mean(error)
plot(t, y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use Runge-Kutta to integrate a second order equation of the form.
% d2y/dt^2+5*dy/dt + 4*y(t) = 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function exerciseB()
t0 = 0;
tfinal = 10;
%                                          initial values
[t, y] = runge_kutta(@second_order_system, [10 1]', [t0 tfinal], 0.01);

y = y(:, 1:length(t));
%plot results
plot(t, y(1,:)), hold
plot(t, y(2,:)), hold

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now use Runge-kutta to integrate the movement of a 1 dof robot arm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function exerciseC()
%these variables are shared by the forward_dynamic_robot1 function defined
%below
global robot tau g
t0 = 0;
tfinal = 10;
robot = load_robot('example','1dofplanar')
robot.dynamics.friction=0
tau = [0];
g = [0 -9.81 0]';
[t, y] = runge_kutta(@forward_dynamic_robot1, [0 0]', [t0 tfinal], 0.01);

figure, plot(y')
% Animate the movement. Change speed from 1-30-100-200
speed = 5
animate(robot,[y(1,1:speed:length(y)); y(2,1:speed:length(y))])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now use Runge-Kutta to simulate the movement of a 2 DOF robot arm.
function exerciseD()
global robot tau g
t0 = 0;
tfinal = 10;
tau = [0 0];
g = [0 -9.81 0]';
robot = load_robot('example','2dofplanar')
robot.dynamics.friction=1

[t, y] = runge_kutta(@forward_dynamic_robot2, [0 0 0 0]', [t0 tfinal], 0.01);

%speed, change to 10, 30, 50, 100
speed = 5
animate(robot,[y(1,1:speed:length(y)); y(2,1:speed:length(y))])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function.
% The function returns dy/dt = 2*t. Function called from exerciseA()
%
% Integrate a line in time. dy/dt = 2*t. Obviously, the integration should
% yield. y(t) = t^2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy = line(t, y)
dy = 2*t;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function to solve a second order differential equation.
% Called from function exerciseB()
%
% In order to solve for d2y/dt^2 + 5*dy/dt + 4*y(t) = 0. We can use:
% x1 = y
% x2 = dy/dt
% thus,
% dx1/dt = x2
% dx2/dt = d2y/dt^2=-5*dy/dt-4*y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xd = second_order_system(t, y)
% We must return the solution of
% [dx1/dt; dx2/dt]
xd(1) = y(2);
xd(2) = 0 - 5*y(2)-4*y(1);
xd = xd(:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function to simulate the movement of a 1 DOF robot
% Called from function exerciseC()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xd = forward_dynamic_robot1(t, y)
global tau g robot

qdd = accel(robot, y(1), y(2), tau, g);
%return qd, qdd
xd = [y(2); qdd];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function to simulate the movement of a 2 DOF robot
% Called from function exerciseD()
% t is only used to plot time on screen during simulation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xd = forward_dynamic_robot2(t, y)
global tau g robot
%we must return the solution of
% [dx1/dt; dx2/dt]
t
% caution, g in the forwarddynamics_2dofplanar takes the opposite sign
qdd=forwarddynamics_2dofplanar(robot, y(1:2,1), y(3:4,1), tau', -sum(g), [0 0 0 0 0 0]);
xd = [y(3:4,1); qdd];

repo: ARTE


text:
function prob = mosekProb(H,f,A,rl,ru,Q,l,qrl,qru,lb,ub,xint,sdcone,x0,opts)
%MOSEKPROB  From general inputs check then build MOSEK problem
%
%   prob = mosekProb(H,f,A,rl,ru,Q,l,qrl,qru,lb,ub,xint,sdcone,x0,opts)

%   This function is based in parts on examples from the MOSEK Toolbox,
%   Copyright (c) 1998-2011 MOSEK ApS, Denmark.

%   Copyright (C) 2012 Jonathan Currie (IPL)

%Get warning status
warn = strcmpi(opts.warnings,'all');

%Get sizes
ndec = length(f);
[nlin,Ad] = size(A);

%Error checks
if(~isempty(H))
    [r1,c1] = size(H);
    if(r1 ~= ndec || c1 ~= ndec), error('H is the wrong size, expected %d x %d',ndec,ndec); end
    %Symmetry check is done in OPTI, won't repeat here
end
if(isempty(f)), error('You must supply the f vector to a MOSEK problem!'); end
if(xor(isempty(A),isempty(rl)) || xor(isempty(A),isempty(ru))), error('You must supply A + rl + ru for linear constraints'); end
if(~isempty(A))
    if(Ad ~= ndec), error('A is the wrong size, expected %d x %d',nineq,ndec); end
    if(nlin ~= length(rl)), error('The sizes of A and rl do not correspond'); end
    if(nlin ~= length(ru)), error('The sizes of A and ru do not correspond'); end
end
if(~isempty(lb) && length(lb) ~= ndec), error('lb is the wrong length, expected %d x 1',ndec); end
if(~isempty(ub) && length(ub) ~= ndec), error('ub is the wrong length, expected %d x 1',ndec); end
if(~isempty(xint))
    if(length(xint) ~= ndec), error('xint is the wrong length, expected %d x 1',ndec); end
    if(~ischar(xint)), error('xint should be a char array with characters ''C'', ''I'' or ''B'''); end
end
if(~isempty(x0) && length(x0) ~= ndec), error('x0 is the wrong length, expected %d x 1',ndec); end

%Quadratic Objective
if(~isempty(H))
    if(~issparse(H))
        if(warn)
            optiwarn('mosek:sparse','The H matrix should be sparse, correcting: [sparse(H)]');
        end
        [prob.qosubi,prob.qosubj,prob.qoval] = find(sparse(tril(H)));
    else
        [prob.qosubi,prob.qosubj,prob.qoval] = find(tril(H));
    end
end

%Linear Objective (if not sd)
if(isempty(sdcone))
    prob.c = f;
else
    prob.c = zeros(size(f));
end

%Linear constraints
prob.a = A;
if(isempty(prob.a))
    prob.a = sparse(0,length(f));
elseif(~issparse(prob.a))
    if(warn)
        optiwarn('mosek:sparse','The A matrix should be sparse, correcting: [sparse(A)]');
    end
    prob.a = sparse(prob.a);
end
if(size(f,2) > 1), f = f'; end
if(size(rl,2) > 1), rl = rl'; end
if(size(ru,2) > 1), ru = ru'; end
prob.blc = rl;
prob.buc = ru;

%Quadratic constraints
if(~isempty(Q))
    prob = genQC(prob,Q,l,qrl,qru,warn);
end

%Semidefinite Constraints
if(~isempty(sdcone))
    prob = genSDCONE(prob,f,sdcone);
end

%Integer variables
if(~isempty(xint))
    xint = lower(xint); [r,c] = size(xint);
    if(r > c), xint = xint'; end
    cind = strfind(xint,'c');
    bind = strfind(xint,'b');
    iind = strfind(xint,'i');
    if(length([cind bind iind]) ~= ndec) %double check no other chars
        error('xint should be a char array with characters ''C'', ''I'' or ''B''');
    end
    %Only need to setup if integer vars exist
    if(~isempty(bind) || ~isempty(iind))
        prob.ints.sub = sort([bind iind]); %sorted indices of all binary and integer vars
        %If binary, need to setup binary bounds
        if(~isempty(bind))
            if(isempty(lb)), lb = -Inf(size(f)); end %create false bounds if required
            lb(bind) = 0;
            if(isempty(ub)), ub = Inf(size(f)); end %create false bounds if required
            ub(bind) = 1;
        end
    end
end

%Bounds (note must be done after binary stuff above)
prob.blx = lb;
prob.bux = ub;

%Initial guess
if(~isempty(x0))
    if(isfield(prob,'ints'))
        prob.sol.int.xx = x0;
    else
        prob.sol.itr.xx = x0;
    end
end

%Determine problem type
if(isempty(H)) %linear or sdp
    if(isempty(xint))
        if(isempty(sdcone))
            prob.names.name = 'OPTI MOSEK LP Problem';
        else
            prob.names.name = 'OPTI MOSEK SDP Problem';
        end
    elseif(all(lower(xint) == 'b'))
        prob.names.name = 'OPTI MOSEK BILP Problem';
    else
        prob.names.name = 'OPTI MOSEK MILP Problem';
    end
else %quadratic
    if(isempty(Q))
        if(isempty(xint))
            prob.names.name = 'OPTI MOSEK QP Problem';
        else
            prob.names.name = 'OPTI MOSEK MIQP Problem';
        end
    else %quadratically constrained
        if(isempty(xint))
            prob.names.name = 'OPTI MOSEK QCQP Problem';
        else
            prob.names.name = 'OPTI MOSEK MIQCQP Problem';
        end
    end
end

function prob = genQC(prob,Q,l,qrl,qru,warn)
%Generate MOSEK quadratic constraints

%Create default fields
prob.qcsubi = [];
prob.qcsubj = [];
prob.qcsubk = [];
prob.qcval = [];

ndec = length(prob.c);

%Error checks
if(~isempty(Q))
    if(isempty(l) || isempty(qrl) || isempty(qru))
        error('You must supply Q, l, qrl and qru for quadratic constraints');
    end
end

%Process cell by cell
if(iscell(Q))
   if(iscell(l) || iscell(qrl) || iscell(qru))
       error('Only Quadratic Constraints Q may be a cell!');
   end
   %Check overall sizes
   [r0,c0] = size(Q);
   [r1,c1] = size(l);
   [r2,c2] = size(qrl);
   [r3,c3] = size(qru);
   if(r0 > c0)
       Q = Q';
       [~,c0] = size(Q);
   end
   if(r1 ~= c1) %not square
       if(r1 == c0)
           l = l';
           c1 = size(l,2);
       end
   end
   if(r2 > c2)
       qrl = qrl';
       [~,c2] = size(qrl);
   end
   if(r3 > c3)
       qru = qru';
       [~,c3] = size(qru);
   end
   if(c0 ~= c1 || c0 ~= c2 || c0 ~= c3)
       error('Quadratic Constraints Q + l + qrl + qru are not the same length!');
   end
   %For each cell
   for i = 1:c2
        Qt = Q{i};
        lt = l(:,i);
        rlt = qrl(i);
        rut = qru(i);
        %Error Check
        checkQC(Qt,lt,rlt,rut,ndec,sprintf('in cell %d',i),sprintf('in column %d',i));
        %Check sparsity
        if(~issparse(Qt))
            if(warn)
                optiwarn('mosek:sparse','The Q matrix in cell %d should be sparse, correcting: [sparse(Q)]',i);
            end
            Qt = sparse(Qt);
        end
        %Build structure fields
        [qcsubi,qcsubj,qcval] = find(tril(Qt));
        qcval = 2*qcval; %MOSEK includes 1/2x'Qx
        %Concatenate to original vectors
        prob.qcsubi = [prob.qcsubi; qcsubi];
        prob.qcsubj = [prob.qcsubj; qcsubj];
        prob.qcval = [prob.qcval; qcval];
        %Insert linear part into a matrix
        prob = insertQClin(prob,lt,rlt,rut,length(qcval));
   end
else  %Single Quadratic Constraint
    if(size(l,2) > size(l,1))
        l = l';
    end
    %Error Check
    checkQC(Q,l,qrl,qru,ndec,'','');
    %Check sparsity
    if(~issparse(Q))
        if(warn)
            optiwarn('mosek:sparse','The Q matrix should be sparse, correcting: [sparse(Q)]');
        end
        Q = sparse(Q);
    end
    %Build structure fields
    [prob.qcsubi,prob.qcsubj,prob.qcval] = find(tril(Q));
    prob.qcval = 2*prob.qcval; %MOSEK includes 1/2x'Qx
    %Insert linear part into a matrix
    prob = insertQClin(prob,l,qrl,qru,length(prob.qcval));
end

function prob = genSDCONE(prob,f,sdcone)
%Generate MOSEK semidefinite constraints
if(~iscell(sdcone)), sdcone = {sdcone}; end
ndec = length(f);
%Create default fields
prob.bardim    = [];
prob.barc.subj = []; %which cone
prob.barc.subk = []; %row
prob.barc.subl = []; %col
prob.barc.val  = [];
prob.bara.subi = []; %which A matrix
prob.bara.subj = []; %which cone
prob.bara.subk = []; %row
prob.bara.subl = []; %col
prob.bara.val  = [];

%For each sdcone
for i = 1:length(sdcone)
    dim = chkSDDim(sdcone{i},ndec,i);
    %Extract C, enter
    prob.bardim = [prob.bardim dim];
    [ci,cj,cv] = find(tril(reshape(sdcone{i}(:,1),dim,dim)));
    prob.barc.subj = [prob.barc.subj repmat(i,1,length(ci))];
    prob.barc.subk = [prob.barc.subk ci'];
    prob.barc.subl = [prob.barc.subl cj'];
    prob.barc.val = [prob.barc.val -cv'];
    %For each A, enter
    for j = 1:ndec
        [ai,aj,av] = find(tril(reshape(sdcone{i}(:,j+1),dim,dim)));
        prob.bara.subi = [prob.bara.subi repmat(j,1,length(ai))];
        prob.bara.subj = [prob.bara.subj repmat(i,1,length(ai))];
        prob.bara.subk = [prob.bara.subk ai'];
        prob.bara.subl = [prob.bara.subl aj'];
        prob.bara.val = [prob.bara.val av'];
    end
end
%Modify a
prob.a = [prob.a; sparse(ndec,ndec)];
%Set bounds
prob.blc = [prob.blc;f];
prob.buc = [prob.buc;f];

function checkQC(Q,l,qrl,qru,ndec,msgC,msgV)
%Complete QC checks
if((ndec ~= size(Q,1)) || (ndec ~= size(Q,2)))
    error('Constraint Q matrix %s is the wrong size! Expected %d x %d',msgC,ndec,ndec);
end
if(ndec ~= size(l,1))
   error('Constraint l vector %s is the wrong size! Expected %d x 1',msgV,ndec);
end
if(~isscalar(qrl))
   error('Constraint qrl %s is the wrong size! Expected 1 x 1 ',msgV);
end
if(~isscalar(qru))
   error('Constraint qru %s is the wrong size! Expected 1 x 1 ',msgV);
end
try
    e = eig(Q);
catch
    e = eigs(Q);
end
if(any(e < 0))
   error('Constraint Q matrix %s is not positive semidefinite!',msgC);
end
if(all(e == 0))
   error('Constraint Q matrix %s is indefinite!',msgC);
end
sym = abs(tril(Q,-1)-triu(Q,1)') > 1e-9; %tol may need to be adjusted
if(any(any(sym)))
    error('Constraint Q matrix %s is not symmetric',msgC);
end

function dim = chkSDDim(cone,ndec,i)
if(size(cone,2) ~= ndec+1)
    error('Semidefinite Cone %d does not have the correct number of columns',i);
end
dim = sqrt(size(cone,1));
if(floor(dim) ~= dim)
    error('Semidefinite Cone %d does not have the correct number of rows to form a square matrix');
end

function prob = insertQClin(prob,l,qrl,qru,n)
%Insert linear components of qc into prob structure

%Generate k vector
k = size(prob.a,1) + 1;
prob.qcsubk = [prob.qcsubk; k*ones(n,1)];
%Augment to linear a and constraint bounds
prob.a = [prob.a; l'];
prob.blc = [prob.blc; qrl];
prob.buc = [prob.buc; qru];


text:
function v = subSet(v,ind)
% subindex vector3d

v.x = v.x(ind);
v.y = v.y(ind);
v.z = v.z(ind);

end


text:
function [y, b] = remove_baseline(y, sn)
% estiamte baseline of the calcium traces and subtract it
if ~exist('sn', 'var') || isempty(sn)
    sn = get_noise_fft(reshape(y, 1, []));
end
sz = size(y);
y = reshape(y, 1, []);
y_diff = [-1, diff(y)];
b = median(y(and(y_diff>=0, y_diff<sn)));
y = reshape(y-b, sz);

author: flatironinstitute

repo: CaImAn-MATLAB


text:
% function calvinNNDetection()
%
% Copyright by Holger Caesar, 2016

% Global variables
global glDatasetFolder glFeaturesFolder;
assert(~isempty(glDatasetFolder) && ~isempty(glFeaturesFolder));

%%% Settings
% Dataset
vocYear = 2010;
trainName = 'train';
testName  = 'val';

% Specify paths
vocName = sprintf('VOC%d', vocYear);
datasetDir = [fullfile(glDatasetFolder, vocName), '/'];
outputFolder = fullfile(glFeaturesFolder, 'CNN-Models', 'FRCN', vocName, sprintf('%s-testRelease', vocName));
netPath = fullfile(glFeaturesFolder, 'CNN-Models', 'matconvnet', 'imagenet-vgg-verydeep-16.mat');
logFilePath = fullfile(outputFolder, 'log.txt');

% Fix randomness
randSeed = 42;
rng(randSeed);

% Setup dataset specific options and check validity
setupDataOpts(vocYear, testName, datasetDir);
global DATAopts;
assert(~isempty(DATAopts), 'Error: Dataset not initialized properly!');

% Task-specific
nnOpts.testFn = @testDetection;
nnOpts.misc.overlapNms = 0.3;
nnOpts.derOutputs = {'objective', 1, 'regressObjective', 1};

% General
nnOpts.batchSize = 2;
nnOpts.numSubBatches = nnOpts.batchSize; % 1 image per sub-batch
nnOpts.weightDecay = 5e-4;
nnOpts.momentum = 0.9;
nnOpts.numEpochs = 16;
nnOpts.learningRate = [repmat(1e-3, 12, 1); repmat(1e-4, 4, 1)];
nnOpts.misc.netPath = netPath;
nnOpts.expDir = outputFolder;
nnOpts.gpus = 1; % for automatic selection use: SelectIdleGpu();

% Create outputFolder
if ~exist(outputFolder, 'dir')
    mkdir(outputFolder);
end

% Start logging
diary(logFilePath);

%%% Setup
% Start from pretrained network
net = load(nnOpts.misc.netPath);

% Setup imdb
imdb = setupImdbDetection(trainName, testName, net);

% Create calvinNN CNN class
% By default, network is transformed into fast-rcnn with bbox regression
calvinn = CalvinNN(net, imdb, nnOpts);

%%% Train
calvinn.train();

%%% Test
stats = calvinn.test();

%%% Eval
evalDetection(testName, imdb, stats, nnOpts);

author: nightrome

repo: matconvnet-calvin


text:
function out=proj_Lorentz(x)
%PROJ_LORENTZ computes the orthogonal projection onto the Lorentz cone {x:||x(1,..,n)||<=x(n+1)}
%
%  Usage:
%  out = PROJ_LORENTZ(x)
%  ===========================================
%  Input:
%  x - (n+1)-length vector to be projected
%  ===========================================
%  Output:
%  out - projection vector

% This file is part of the FOM package - a collection of first order methods for solving convex optimization problems
% Copyright (C) 2017 Amir and Nili Beck
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

%reading the user
if (nargin < 1)
    error ('usage: proj_Lorentz(x)') ;
end
n = length(x)-1;
s = x(n+1) ;
x = x(1:n) ;

if (norm(x,s) >= abs(s))
    outx = (norm(x,2) + s)/(2*norm(x,2)) * x;
    outs = (norm(x,2) + s)/2 ;
else
    if (norm(x,2) <= s)
        outx = x ;
        outs = s ;
    else
        % s < norm(x,2) < -s
        outx = zeros(n,1) ;
        outs = 0 ;
    end
end

out = [outx; outs] ;

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function a = jacobi ( m, n )

%*****************************************************************************80
%
%% JACOBI returns the JACOBI matrix.
%
%  Formula:
%
%    if ( J = I - 1 )
%      A(I,J) = 0.5 * sqrt ( ( 4 * J^2 ) / ( 4 * J^2 - 1 ) )
%    else if ( J = I + 1 )
%      A(I,J) = 0.5 * sqrt ( ( 4 * (J-1)^2 ) / ( 4 * (J-1)^2 - 1 ) )
%    else
%      A(I,J) = 0
%
%  Example:
%
%    M = 4, N = 4
%
%    0            0.577350269  0            0
%    0.577350269  0            0.516397779  0
%    0            0.516397779  0            0.507092553
%    0            0            0.507092553  0
%
%  Properties:
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A has a zero diagonal.
%
%    The eigenvalues of A are the zeros of the Legendre polynomial
%    of degree N.  They lie symmetrically in [-1,1], and are also
%    the nodes of Gauss-Legendre quadrature.  For the case of N = 4,
%    these eigenvalues are
%
%      [ -0.861136312, -0.339981044, +0.339981044, +0.861136312 ].
%
%    It follows that A is singular when N is odd.
%
%    The J-th Gauss-Legendre weight is twice the square of the first
%    component of the J-th eigenvector of A.  For the case of N = 4,
%    the eigenvector matrix is:
%
%      -0.417046     -0.571028     -0.571028    0.417046
%       0.622037      0.336258     -0.336258    0.622038
%      -0.571028      0.417046      0.417046    0.571028
%       0.336258     -0.622037      0.622038    0.336258
%
%    and the corresponding weights are
%
%      [ 0.347854845, 0.652145155, 0.652145155, 0.347854845 ]
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Lloyd Trefethen, David Bau,
%    Numerical Linear Algebra,
%    SIAM, 1997, pages 287-292.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of A.
%
%    Output, real A(M,N), the matrix.
%
  a = zeros ( m, n );

  for i = 1 : m
    for j = 1 : n

      if ( j == i - 1 )
        a(i,j) = 0.5 * sqrt ( ( 4 * j * j ) / ( 4 * j * j - 1 ) );
      elseif ( j == i + 1 )
        a(i,j) = 0.5 * sqrt ( ( 4 * ( j - 1 ) * ( j - 1 ) ) ...
          / ( 4 * ( j - 1 ) * ( j - 1 ) - 1 ) );
      else
        a(i,j) = 0.0;
      end

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text: bgest

author: petercorke

repo: machinevision-toolbox-matlab


text:
classdef matrix < handle

    methods(Static)
        function result = is_square(A)
            [m, n]  = size(A);
            result = (m == n);
        end

        function result = is_symmetric(A)
            result = spx.matrix.is_square(A);
            if ~result
                return;
            end
            result =  all(all(A==A.'));
        end

        function result = is_hermitian(A)
            result = spx.matrix.is_square(A);
            if ~result
                return;
            end
            result =  all(all(A==A'));
        end

        function result = is_positive_definite(A)
            result = spx.matrix.is_hermitian(A);
            if ~result
                return;
            end
            % Perform Cholesky decomposition of A
            [R, p] = chol(A);
            result = (p == 0);
        end

        function result = is_orthogonal(A)
            G = A' * A;
            [m, n] = size(G);
            idx = eye(m, n);
            G = abs(G(~idx));
            result = all(G <= 1e-12);
        end

        function result = is_identity(A)
            [m, n]  = size(A);
            if m ~= n
                result = false;
                return;
            end
            I = eye(m);
            diff = abs(A - I);
            result = all(all(diff <= 1e-12));
        end

        function result = is_orthonormal(A)
            if ~isreal(A)
                result = false;
                return;
            end
            G = A' * A;
            m = size(G, 1);
            I = eye(m);
            diff = abs(G - I);
            result = all(all(diff <= 1e-12));
        end

        function X = off_diagonal_elements(X)
            % Returns a column vector of off diagonal elements
            [m, n] = size(X);
            idx = eye(m, n);
            X = X(~idx);
        end

        function X = off_diagonal_matrix(X)
            % Returns the matrix of off diagonal elements
            [m, n] = size(X);
            idx = eye(m, n);
            X = (1 - idx).* X;
        end

        function X = off_diag_upper_tri_elements(X)
            % Returns the upper triangular off diagonal elements in a vector
            [m, n] = size(X);
            idx = tril(ones(m, n));
            X = X(~idx);
        end

        function X = off_diag_upper_tri_matrix(X)
            % Returns the upper triangular off diagonal elements in a matrix
            [m, n] = size(X);
            idx = tril(ones(m, n));
            X = (1 - idx).* X;
        end

        function result = nonzero_density(X)
            % Density of nonzero entries in the matrix
            result = nnz(X) / numel(X);
        end

        function [ result ] = is_diagonally_dominant( A, strict)
        %ISDIAGONALLYDOMINANT Returns if A is a diagonally dominant
        %matrix
        if ~exist('strict', 'var')
            strict = true;
        end
        A = abs(A);
        % Extract the diagonal
        d = diag(A);
        % Set the diagonal elements to 0
        A(logical(eye(size(A)))) = 0;
        % Now sum over the rows
        s = sum(A, 2);
        % We now check whether
        if strict
            result = all(d > s);
        else
            result = all(d >= s);
        end
        end

        function [ A ] = make_diagonally_dominant( A, strict )
        %MAKEDIAGONALLYDOMINANT Makes a matrix diagonally dominant
        if ~exist('strict', 'var')
            strict = true;
        end
        B = abs(A);
        % Set the diagonal elements to 0
        B(logical(eye(size(B)))) = 0;
        % Now sum over the rows
        s = sum(B, 2);
        if strict
            d = s + 1;
        else
            d = s;
        end
        % Extract the diagonal elements of A
        dd = diag(A);
        % Identify the ones which are not dominant
        requiredChanges = abs(dd) < d;
        % Assign the updated values
        dd(requiredChanges) = d(requiredChanges);
        % Update the matrix
        A(logical(eye(size(A)))) = dd;
        end

        function compare(A, B)
            % Compares two matrices
            C = abs(A - B);
            max_diff = max(max(abs(C)));
            fprintf('Maximum difference: %.4e\n', max_diff);
        end

        function df = dof_lowrank(m, n, r)
            % Returns the degrees of freedom of a low rank matrix
            df = r * (m + n - r);
        end

    end

end

author: indigits

repo: sparse-plex


text:
% Authors: A. Iscen, G. Tolias, Y. Avrithis, T. Furon, O. Chum. 2017.
% function to perform diffusion
% solving system A*f = y
function f = dfs(A, y, tol, it)
	if nargin < 4, it = 20;	end
	if nargin < 3, tol = 1e-10; end

   [f,~,~,~] = pcg(A,y,tol,it);

author: ahmetius

repo: diffusion-retrieval


text:
% Directly plots a simple figure using plotPub function.

clear all;
addpath('../lib');

%% lets plot 3 cycles of 50Hz AC voltage
f = 50;
Vm = 10;
phi = 0;

% generate the signal
t = [0:0.0001:3/f];
th = 2*pi*f*t;
v = Vm*sin(th+phi);

%% plot now
plotx{1} = t*1E3; %convert time in ms and create a cell array
ploty{1} = v; % assign v to a cell array
opt.XLabel = 'Time, t (ms)'; % xlabel
opt.YLabel = 'Voltage, V (V)'; %ylabel

% Save? comment the following line if you do not want to save
opt.FileName = 'plotDirectSimple1.eps';

% create the plot
plotPub(plotx, ploty, 1, opt);

author: masumhabib

repo: PlotPub

save_path: github-repos/MATLAB/masumhabib-PlotPub


text:
function [ newPositions, zerosInARow ] = CalcSimpleOpticalFlowHists( centroids, img1, train, neighborSize, windowSize, w, thresh, zerosInARow)
%  calcSimpleOpticalFlowHists - Used to track object of interest
%--------------------------------------------------------------------------
%   Params: centroids: n x 3 matrix where n is the number of centroids showing
%           img1: image to be looked at for new centroid position
%           train: training histograms
%           neighborSize: size of neighborhood to look at.  Should be odd
%           windowSize: size of windows in centroid neighborhood
%           w - the width of the bins for the RGB color
%               histograms
%           thresh: distance threshold
%           zerosInARow: zeros in a row
%
%   Returns: newPositions - new positions of centroid.  0 for x and y of
%               centroid if its out of view bounds or doesn't
%               have a good enough match in the view
%            zerosInARow - number of frames straight where zeros
%               detected
%
%   Assumes: brightness constant, small motion changes betw frames
%--------------------------------------------------------------------------

newPositions = zeros(size(centroids,1), size(centroids,2));
imgWidth = size(img1,2);
imgHeight = size(img1,1);
for i = 1:size(centroids,1)
    currCentr = centroids(i,:);
    row = currCentr(1);
    col = currCentr(2);
%     if useOrigColor == 1
%         pixRef = [currCentr(4) currCentr(5) currCentr(6)]';
%     else
%         pixRef = double(refImg(row,col,:));
%         pixRef = pixRef(:);
%     end
    indent = floor(neighborSize / 2);
    bestDist = 0;
    for kRow = -indent:indent
        for kCol = -indent:indent
            testRow = max(1, row + kRow);
            testRow = min(imgHeight, testRow);
            testCol = max(1, col + kCol);
            testCol = min(imgWidth, testCol);
            hist = SimpleHist1D(img1(max(1,testRow - windowSize):min(imgHeight, testRow + windowSize),...
                max(1,testCol- windowSize):min(imgWidth, testCol + windowSize)...
                ,:), w);
            dist = Score1D(hist, train);
            %prefer points closer to old centroid so put penalty
            manhDist2center = (abs(testRow - centroids(i,1)) + ...
              abs(testCol - centroids(i,2))) / 2;
            %good vals is - 15, - 20,
            if (dist > exp(thresh - 25 + manhDist2center))
                if (dist > bestDist)
                    newPositions(i,1) = testRow;
                    newPositions(i,2) = testCol;
                    bestDist = dist;
                end
            end
        end
    end
    if (newPositions(i,1) == 0)
        zerosInARow = zerosInARow + 1;
        if (zerosInARow <= 15)
            newPositions(i,1) = centroids(i,1);
            newPositions(i,2) = centroids(i,2);
        end
    else
        zerosInARow = 0;
    end
    %newPositions
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function rot = fit(r,l,varargin)
% find rotation rot such that l = rot * r
%
% Syntax
%
%   rot = rotation.fit(r,l)
%   rot = rotation.fit(r,l,'weights',w)
%
% Input
%  r, l - @vector3d
%
% Output
%  rot - @rotation
%
% Description
% Find the rotation that best maps all the vectors |r| onto the vectors |l|.
%
% See also
% rotation/rotation rotation/byMatrix rotation/byAxisAngle
% rotation/byEuler
%
% References
%
% * W. Kabsch, A solution for the best rotation to relate two sets of
% vectors, Acta Cryst. (1976). A32, 922.
%
% * B. K. P. Horn, Closed-form solution of absolute orientation using unit
% quaternions, Journal of the Optical Society of America A (1987), Vol, 4,
% 629.
%

if check_option(varargin,'antipodal') || r.antipodal || l.antipodal
  warning('antipodal symmetry not yet supported in rotation.fit');
end

w = get_option(varargin,'weights');
if ~isempty(w), r = r .* w; end

M = r.normalize * l.normalize;

switch lower(get_option(varargin,'method','horn'))

  case 'horn' % quaternion based method

    MS = M + M.';
    MA = M - M.';

    N = [trace(M) MA(2,3)             MA(3,1)              MA(1,2);
      MA(2,3)    M(1,1)-M(2,2)-M(3,3) MS(1,2)              MS(1,3);
      MA(3,1)    MS(1,2)              M(2,2)-M(1,1)-M(3,3) MS(2,3);
      MA(1,2)    MS(1,3)              MS(2,3)              M(3,3)-M(1,1)-M(2,2)];

    [V,d] = eig(N,'vector');
    [~, ind] = sort(d);

    rot = rotation(V(:,ind(4)).');

  case 'kabsch' % Kabsch algorithm

    [U,~,V] = svd(M);

    if length(r)<3
      S = ones(3);
      S(3,3) = det(V*U');
      rot = rotation.byMatrix(V*S*U');
    else
      rot = rotation.byMatrix(V*U');
      if r.antipodal || l.antipodal, rot.i = false; end
    end
end
end

function check

r = rotation.rand;
% u = [xvector;yvector;zvector];
u = vector3d.rand(randi(100,1));

v = (r * u).';

% add some noise
vn = v + 0.1 .* vector3d.rand(size(v));

% fit a rotation
%rec = rotation.fit(u,vn,'method','kabsch');
rec = rotation.fit(u,vn,'method','horn');

% the distance to the initial rotation
angle(r,rec)./degree

% the fit
f = sum(dot(rec * u,vn.'));

% check for local minimum
S3G = localOrientationGrid(specimenSymmetry,specimenSymmetry,1*degree,'resolution',0.5*degree);
all(f > sum(dot(rec * S3G * u,vn.'),2))

end

author: mtex-toolbox

repo: mtex