text:
%RangeBearingSensor Range and bearing sensor class
%
% A concrete subclass of the Sensor class that implements a range and bearing
% angle sensor that provides robot-centric measurements of landmark points in
% the world. To enable this it holds a references to a map of the world (LandmarkMap object)
% and a robot (Vehicle subclass object) that moves in SE(2).
%
% The sensor observes landmarks within its angular field of view between
% the minimum and maximum range.
%
% Methods::
%
% reading   range/bearing observation of random landmark
% h         range/bearing observation of specific landmark
% Hx        Jacobian matrix with respect to vehicle pose dh/dx
% Hp        Jacobian matrix with respect to landmark position dh/dp
% Hw        Jacobian matrix with respect to noise dh/dw
%-
% g         feature position given vehicle pose and observation
% Gx        Jacobian matrix with respect to vehicle pose dg/dx
% Gz        Jacobian matrix with respect to observation dg/dz
%
% Properties (read/write)::
% W            measurement covariance matrix (2x2)
% interval     valid measurements returned every interval'th call to reading()
% landmarklog  time history of observed landmarks
%
% Reference::
%
%   Robotics, Vision & Control, Chap 6,
%   Peter Corke,
%   Springer 2011
%
% See also Sensor, Vehicle, LandmarkMap, EKF.

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

classdef RangeBearingSensor < Sensor

    properties
        W           % measurment covariance
        r_range     % range limits
        theta_range % angle limits

        randstream  % random stream just for Sensors

        landmarklog  % time history of observed landmarks
    end

    properties (SetAccess = private)
        count       % number of reading()s
    end

    methods

        function s = RangeBearingSensor(robot, map, varargin)
            %RangeBearingSensor.RangeBearingSensor Range and bearing sensor constructor
            %
            % S = RangeBearingSensor(VEHICLE, MAP, OPTIONS) is an object
            % representing a range and bearing angle sensor mounted on the Vehicle
            % subclass object VEHICLE and observing an environment of known landmarks
            % represented by the LandmarkMap object MAP.  The sensor covariance is W
            % (2x2) representing range and bearing covariance.
            %
            % The sensor has specified angular field of view and minimum and maximum
            % range.
            %
            % Options::
            % 'covar',W               covariance matrix (2x2)
            % 'range',xmax            maximum range of sensor
            % 'range',[xmin xmax]     minimum and maximum range of sensor
            % 'angle',TH              angular field of view, from -TH to +TH
            % 'angle',[THMIN THMAX]   detection for angles betwen THMIN
            %                         and THMAX
            % 'skip',K                return a valid reading on every K'th call
            % 'fail',[TMIN TMAX]      sensor simulates failure between
            %                         timesteps TMIN and TMAX
            % 'animate'               animate sensor readings
            %
            % See also options for Sensor constructor.
            %
            % See also RangeBearingSensor.reading, Sensor.Sensor, Vehicle, LandmarkMap, EKF.

            % call the superclass constructor
            s = s@Sensor(robot, map, varargin{:});

            s.randstream = RandStream.create('mt19937ar');

            opt.range = [];
            opt.angle = [];
            opt.covar = zeros(2,2);

            [opt,args] = tb_optparse(opt, varargin);

            s.W = opt.covar;
            if ~isempty(opt.range)
                if length(opt.range) == 1
                    s.r_range = [0 opt.range];
                elseif length(opt.range) == 2
                    s.r_range = opt.range;
                end
            end
            if ~isempty(opt.angle)
                if length(opt.angle) == 1
                    s.theta_range = [-opt.angle opt.angle];
                elseif length(opt.angle) == 2
                    s.theta_range = opt.angle;
                end
            end

            s.count = 0;
            s.verbose = opt.verbose;
        end

        function init(s)
            s.landmarklog = [];
        end

        function k = selectFeature(s)
            k = s.randstream.randi(s.map.nlandmarks);
        end

        function [z,jf] = reading(s)
            %RangeBearingSensor.reading Choose landmark and return observation
            %
            % [Z,K] = S.reading() is an observation of a random visible landmark where
            % Z=[R,THETA] is the range and bearing with additive Gaussian noise of
            % covariance W (property W). K is the index of the map feature that was
            % observed.
            %
            % The landmark is chosen randomly from the set of all visible landmarks,
            % those within the angular field of view and range limits.  If no valid
            % measurement, ie. no features within range, interval subsampling enabled
            % or simulated failure the return is Z=[] and K=0.
            %
            % Notes::
            % - Noise with covariance W (property W) is added to each row of Z.
            % - If 'animate' option set then show a line from the vehicle to the
            %   landmark
            % - If 'animate' option set and the angular and distance limits are set
            %   then display that region as a shaded polygon.
            % - Implements sensor failure and subsampling if specified to constructor.
            %
            % See also RangeBearingSensor.h.

            % TODO probably should return K=0 to indicated invalid

            % model a sensor that emits readings every interval samples
            s.count = s.count + 1;

            % check conditions for NOT returning a value
            z = [];
            jf = 0;
            % sample interval
            if mod(s.count, s.interval) ~= 0
                return;
            end
            % simulated failure
            if ~isempty(s.fail) && (s.count >= s.fail(1)) && (s.count <= s.fail(2))
                return;
            end

            % create a polygon to indicate the active sensing area based on range+angle limits
            if s.animate && ~isempty(s.theta_range) && ~isempty(s.r_range)
                h = findobj(gca, 'tag', 'sensor-area');
                if isempty(h)

                    th=linspace(s.theta_range(1), s.theta_range(2), 20);
                    x = s.r_range(2) * cos(th);
                    y = s.r_range(2) * sin(th);
                    if s.r_range(1) > 0
                        th = flip(th);
                        x = [x s.r_range(1) * cos(th)];
                        y = [y s.r_range(1) * sin(th)];
                    else
                        x = [x 0];
                        y = [y 0];
                    end
                    % no sensor zone, create one
                    plot_poly([x; y], 'fillcolor', 'r', 'alpha', 0.1, 'edgecolor', 'none', 'animate', 'tag', 'sensor-area');
                else
                    %hg = get(h, 'Parent');
                    plot_poly(h, s.robot.x);

                end
            end

            if ~isempty(s.r_range)  || ~isempty(s.theta_range)
                % if range and bearing angle limits are in place look for
                % any landmarks that match criteria

                % get range/bearing to all landmarks, one per row
                z = s.h(s.robot.x');
                jf = 1:numcols(s.map.map);

                if ~isempty(s.r_range)
                    % find all within range
                    k = find( z(:,1) >= s.r_range(1) & z(:,1) <= s.r_range(2) );
                    z = z(k,:);
                    jf = jf(k);
                end
                if ~isempty(s.theta_range)
                    % find all within angular range as well
                    k = find( z(:,2) >= s.theta_range(1) & z(:,2) <= s.theta_range(2) );
                    z = z(k,:);
                    jf = jf(k);
                end

                % deal with cases for 0 or > 1 features found
                if isempty(z)
                    % no landmarks found
                    jf = 0;
                elseif length(k) >= 1
                    % more than 1 in range, pick a random one
                    i = s.randstream.randi(length(k));
                    z = z(i,:);
                    jf = jf(i);
                end

            else
                % randomly choose the feature
                jf = s.selectFeature();

                % compute the range and bearing from robot to feature
                z = s.h(s.robot.x', jf);
            end

            if s.verbose
                if isempty(z)
                    fprintf('Sensor:: no features\n');
                else
                    fprintf('Sensor:: feature %d: %.1f %.1f\n', jf, z);
                end
            end
            if s.animate
                s.plot(jf);
            end

            z = z';

            % add the reading to the landmark log
            s.landmarklog = [s.landmarklog jf];
        end

        function z = h(s, xv, jf)
            %RangeBearingSensor.h Landmark range and bearing
            %
            % Z = S.h(X, K) is a sensor observation (1x2), range and bearing, from vehicle at
            % pose X (1x3) to the K'th landmark.
            %
            % Z = S.h(X, P) as above but compute range and bearing to a landmark at coordinate P.
            %
            % Z = s.h(X) as above but computes range and bearing to all
            % map features.  Z has one row per landmark.
            %
            % Notes::
            % - Noise with covariance W (propertyW) is added to each row of Z.
            % - Supports vectorized operation where XV (Nx3) and Z (Nx2).
            % - The landmark is assumed visible, field of view and range liits are not
            %   applied.
            %
            % See also RangeBearingSensor.reading, RangeBearingSensor.Hx, RangeBearingSensor.Hw, RangeBearingSensor.Hp.

            % get the landmarks, one per row
            if nargin < 3
                % s.h(XV)
                xlm = s.map.map';
            elseif length(jf) == 1
                % s.h(XV, JF)
                xlm = s.map.map(:,jf)';
            else
                % s.h(XV, XF)
                xlm = jf(:)';
            end

            % Straightforward code:
            %
            % dx = xf(1) - xv(1); dy = xf(2) - xv(2);
            %
            % z = zeros(2,1);
            % z(1) = sqrt(dx^2 + dy^2);       % range measurement
            % z(2) = atan2(dy, dx) - xv(3);   % bearing measurement
            %
            % Vectorized code:

            % compute range and bearing
            dx = xlm(:,1) - xv(:,1); dy = xlm(:,2) - xv(:,2);
            z = [sqrt(dx.^2 + dy.^2) angdiff(atan2(dy, dx), xv(:,3)) ];   % range & bearing measurement

            % add noise with covariance W
            z = z + s.randstream.randn(size(z)) * sqrtm(s.W) ;
        end

        function J = Hx(s, xv, jf)
            %RangeBearingSensor.Hx Jacobian dh/dx
            %
            % J = S.Hx(X, K) returns the Jacobian dh/dx (2x3) at the vehicle
            % state X (3x1) for map landmark K.
            %
            % J = S.Hx(X, P) as above but for a landmark at coordinate P.
            %
            % See also RangeBearingSensor.h.
            if length(jf) == 1
                xf = s.map.map(:,jf);
            else
                xf = jf;
            end
            if isempty(xv)
                xv = s.robot.x;
            end
            Delta = xf - xv(1:2)';
            r = norm(Delta);
            J = [
                -Delta(1)/r,    -Delta(2)/r,        0
                Delta(2)/(r^2),  -Delta(1)/(r^2),   -1
            ];
        end

        function J = Hp(s, xv, jf)
            %RangeBearingSensor.Hp Jacobian dh/dp
            %
            % J = S.Hp(X, K) is the Jacobian dh/dp (2x2) at the vehicle
            % state X (3x1) for map landmark K.
            %
            % J = S.Hp(X, P) as above but for a landmark at coordinate P (1x2).
            %
            % See also RangeBearingSensor.h.
            if length(jf) == 1
                xf = s.map.map(:,jf);
            else
                xf = jf;
            end
            Delta = xf - xv(1:2)';
            r = norm(Delta);
            J = [
                Delta(1)/r,         Delta(2)/r
                -Delta(2)/(r^2),    Delta(1)/(r^2)
            ];
        end

        function J = Hw(s, xv, jf)
            %RangeBearingSensor.Hx Jacobian dh/dw
            %
            % J = S.Hw(X, K) is the Jacobian dh/dw (2x2) at the vehicle
            % state X (3x1) for map landmark K.
            %
            % See also RangeBearingSensor.h.
            J = eye(2,2);
        end

        function xf = g(s, xv, z)
            %RangeBearingSensor.g Compute landmark location
            %
            % P = S.g(X, Z) is the world coordinate (2x1) of a feature given
            % the observation Z (1x2) from a vehicle state with X (3x1).
            %
            % See also RangeBearingSensor.Gx, RangeBearingSensor.Gz.

            range = z(1);
            bearing = z(2) + xv(3); % bearing angle in vehicle frame

            xf = [xv(1)+range*cos(bearing); xv(2)+range*sin(bearing)];
        end

        function J = Gx(s, xv, z)
            %RangeBearingSensor.Gxv Jacobian dg/dx
            %
            % J = S.Gx(X, Z) is the Jacobian dg/dx (2x3) at the vehicle state X (3x1) for
            % sensor observation Z (2x1).
            %
            % See also RangeBearingSensor.g.
            theta = xv(3);
            r = z(1);
            bearing = z(2);
            J = [
                1,   0,   -r*sin(theta + bearing);
                0,   1,    r*cos(theta + bearing)
                ];
        end

        function J = Gz(s, xv, z)
            %RangeBearingSensor.Gz Jacobian dg/dz
            %
            % J = S.Gz(X, Z) is the Jacobian dg/dz (2x2) at the vehicle state X (3x1) for
            % sensor observation Z (2x1).
            %
            % See also RangeBearingSensor.g.
            theta = xv(3);
            r = z(1);
            bearing = z(2);
            J = [
                cos(theta + bearing),   -r*sin(theta + bearing);
                sin(theta + bearing),    r*cos(theta + bearing)
                ];
        end

        function str = char(s)
            str = char@Sensor(s);
            str = char(str, ['W = ', mat2str(s.W, 3)]);

            str = char(str, sprintf('interval %d samples', s.interval) );
            if ~isempty(s.r_range)
                str = char(str, sprintf('range: %g to %g', s.r_range) );
            end
            if ~isempty(s.theta_range)
                str = char(str, sprintf('angle: %g to %g', s.theta_range) );
            end
        end

    end % method
end % classdef

author: petercorke

repo: robotics-toolbox-matlab


text:
function k = triangle_to_i4 ( i, j )

%*****************************************************************************80
%
%% TRIANGLE_TO_I4 converts a triangular coordinate to an integer.
%
%  Discussion:
%
%    Triangular coordinates are handy when storing a naturally triangular
%    array (such as the lower half of a matrix) in a linear array.
%
%    Thus, for example, we might consider storing
%
%    (1,1)
%    (2,1) (2,2)
%    (3,1) (3,2) (3,3)
%    (4,1) (4,2) (4,3) (4,4)
%
%    as the linear array
%
%    (1,1) (2,1) (2,2) (3,1) (3,2) (3,3) (4,1) (4,2) (4,3) (4,4)
%
%    Here, the quantities in parenthesis represent the natural row and
%    column indices of a single number when stored in a rectangular array.
%
%    Thus, our goal is, given the row I and column J of the data,
%    to produce the value K which indicates its position in the linear
%    array.
%
%    The triangular numbers are the indices associated with the
%    diagonal elements of the original array, T(1,1), T(2,2), T(3,3)
%    and so on.
%
%  Formula:
%
%    K = J + ( (I-1) * I ) / 2
%
%  First Values:
%
%     I  J  K
%
%     0  0  0
%     1  1  1
%     2  1  2
%     2  2  3
%     3  1  4
%     3  2  5
%     3  3  6
%     4  1  7
%     4  2  8
%     4  3  9
%     4  4 10
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, J, the row and column indices.  I and J must
%    be nonnegative, and J must not be greater than I.
%
%    Output, integer K, the linear index of the (I,J) element.
%
  if ( i < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGLE_TO_I4 - Fatal error!\n' );
    fprintf ( 1, '  I < 0.\n' );
    fprintf ( 1, '  I = %d\n', i );
    error ( 'TRIANGLE_TO_I4 - Fatal error!' );
  elseif ( j < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGLE_TO_I4 - Fatal error!\n' );
    fprintf ( 1, '  J < 0.\n' );
    fprintf ( 1, '  J = %d\n', j );
    error ( 'TRIANGLE_TO_I4 - Fatal error!' );
  elseif ( i < j )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGLE_TO_I4 - Fatal error!\n' );
    fprintf ( 1, '  I < J.\n' );
    fprintf ( 1, '  I = %d\n', i );
    fprintf ( 1, '  J = %d\n', j );
    error ( 'TRIANGLE_TO_I4 - Fatal error!' );
  end

  k = j + ( ( i - 1 ) * i ) / 2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ result, node_num ] = sphere01_triangle_quad_icos2v ( a_xyz, ...
  b_xyz, c_xyz, factor, fun )

%*****************************************************************************80
%
%% SPHERE01_TRIANGLE_QUAD_ICOS2V: vertex rule, subdivide then project.
%
%  Discussion:
%
%    This function estimates an integral over a spherical triangle on the
%    unit sphere.
%
%    This function sets up an icosahedral grid, and subdivides each
%    edge of the icosahedron into FACTOR subedges.  These edges define a grid
%    within each triangular icosahedral face.   All of these calculations are
%    done, essentially, on the FLAT faces of the icosahedron.  Only then are
%    the triangle vertices projected to the sphere.
%
%    This function uses a more sophisticated projection scheme than
%    SPHERE01_TRIANGLE_QUAD_ICOS1V, but this does not seem to improve
%    the results significantly.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A_XYZ(3), B_XYZ(3), C_XYZ(3), the vertices
%    of the spherical triangle.
%
%    Input, integer FACTOR, the subdivision factor, which must
%    be at least 1.
%
%    Input, function v = FUN ( x ), evaluates the integrand at the point X.
%
%    Output, real RESULT, the estimated integral.
%
%    Output, integer NODE_NUM, the number of evaluation points.
%

%
%  Initialize the integral data.
%
  result = 0.0;
  area_total = 0.0;
  node_num = 0;
%
%  Deal with subtriangles that have same orientation as face.
%
  for f1 = 0 : factor - 1
    for f2 = 0 : factor - f1 - 1
      f3 = factor - f1 - f2;

      a2_xyz = sphere01_triangle_project2 ( a_xyz, b_xyz, c_xyz, f1 + 1, f2,     f3 - 1 );
      b2_xyz = sphere01_triangle_project2 ( a_xyz, b_xyz, c_xyz, f1,     f2 + 1, f3 - 1 );
      c2_xyz = sphere01_triangle_project2 ( a_xyz, b_xyz, c_xyz, f1,     f2,     f3 );

      area = sphere01_triangle_vertices_to_area ( a2_xyz, b2_xyz, c2_xyz );

      node_num = node_num + 3;
      va = fun ( a2_xyz );
      vb = fun ( b2_xyz );
      vc = fun ( c2_xyz );
      result = result + area * ( va + vb + vc ) / 3.0;
      area_total = area_total + area;

      end
    end
%
%  Deal with subtriangles that have opposite orientation as face.
%
  for f3 = 0 : factor - 2
    for f2 = 1 : factor - f3 - 1
      f1 = factor - f2 - f3;

      a2_xyz = sphere01_triangle_project2 ( a_xyz, b_xyz, c_xyz, f1 - 1, f2,     f3 + 1 );
      b2_xyz = sphere01_triangle_project2 ( a_xyz, b_xyz, c_xyz, f1,     f2 - 1, f3 + 1 );
      c2_xyz = sphere01_triangle_project2 ( a_xyz, b_xyz, c_xyz, f1,     f2,     f3 );

      area = sphere01_triangle_vertices_to_area ( a2_xyz, b2_xyz, c2_xyz );

      node_num = node_num + 3;
      va = fun ( a2_xyz );
      vb = fun ( b2_xyz );
      vc = fun ( c2_xyz );
      result = result + area * ( va + vb + vc ) / 3.0;
      area_total = area_total + area;

    end
  end

  return
end


text:
function rad = deg2rad(deg)
%
% Function Name:
%
%   deg2rad - Convert degrees to radians.
%
% Calling Sequence:
%
%   rad = deg2rad(deg);
%
% Parameters:
%
%   deg		: Angle in degrees.
%
%   rad		: Angle in radians.
%
% Description:
%
%   Convenient utility function for converting degrees to radians, which are
%   often the required angular units for functions in the NURBS toolbox.
%
% Examples:
%
%   // Convert 35 degrees to radians
%   rad = deg2rad(35);

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

rad = pi*deg/180.0;

author: Sable

repo: mcbench-benchmarks


text:
function best_acc = D2L2R2_top(dataset, N_train, k, lambda1, lambda2, alpha)
% function D2L2R2_top(dataset, N_train, k, lambda1, lambda2, alpha)
% * The top function of D2L2R2
% * INPUT:
%   + `dataset`: name of the dataset stored in `.mat` file in `data` folder.
%     Note that `dataset` is the file name of the `.mat`, excluding `.mat`.
%   + `N_train`: number of training samples in each class
%   + `k`: number of atoms in EACH dictionary
%   + `lambda1, lambda2, alpha`: regularization parameters.
% * To run an small example, type `D2L2r2-top` without input in MATLAB
%           command window.
% -----------------------------------------------
% Author: Tiep Vu, thv102@psu.edu, 5/13/2016 10:29:48 PM
%         (http://www.personal.psu.edu/thv102/)
% -----------------------------------------------
    addpath(genpath('utils'));
    addpath('ODL');
    addpath('LRSDL_FDDL');
    addpath('D2L2R2');
    %% test mode
    if nargin == 0
        dataset = 'myYaleB';
        N_train = 10;
        k = 8;
        lambda1 = 0.001;
        lambda2 = 0.01;
        alpha = 0.01;
    end
    %% Data preparation
    t = getTimeStr();
    [dataset, Y_train, Y_test, label_train, label_test] = train_test_split(...
        dataset, N_train);
    %% main
    [acc, rt] = D2L2R2_wrapper(Y_train, label_train, Y_test, label_test,...
                        k, lambda1, lambda2, alpha);
    %% save results
    if ~exist('results', 'dir')
        mkdir('results');
    end
    if ~exist(fullfile('results', 'D2L2R2'), 'dir')
        mkdir('results', 'D2L2R2');
    end
    fn = fullfile('results', 'D2L2R2', strcat(dataset, '_N_', ...
        num2str(N_train), '_k_', num2str(k), '_l1_', num2str(lambda1),...
        '_l2_', num2str(lambda2), '_a_', num2str(alpha), '_', t, '.mat'));
    disp(fn);

    save(fn, 'acc', 'rt');
    best_acc = max(acc);
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
classdef testStoredComputedChecker < handle

    properties (Abstract, Access = protected)
        variablesToStore
        storedVar
        error
    end

    properties (Access = protected)
        computedVar
    end

    methods (Access = protected)
        function computeError(obj)
            d = numel(obj.variablesToStore);
            err = ones(d,1);
            for ivar = 1:d
                sV = obj.storedVar{ivar};
                cV = obj.computedVar{ivar};
                err(ivar) = norm(sV - cV)/norm(sV);
            end
            obj.error = norm(err);
        end
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function out = normalize_rows(in)
    out = scale_rows(in,1./row_sum(in));
% rowSums = sum(x,2);
% [I,J] = size(x);
% p = zeros(I,J);
% for i=1:I
%     p(i,:) = x(i,:) / rowSums(i);
% end
%
end

author: jacobeisenstein

repo: SAGE


text:
function [Eft, Varft, lpyt, Eyt, Varyt] = gpsvi_predgrad(gp,x,y,varargin)
%GPSVI_PREDGRAD  Make predictions with SVI GP
%
%  Description
%    [EFT, VARFT] = GPSVI_PREDGRAD(GP, X, Y, XT, OPTIONS)
%    takes a GP structure together with matrix X of training inputs and
%    vector Y of training targets, and evaluates the predictive
%    distribution at test inputs XT. Returns a posterior mean EFT and
%    variance VARFT of latent variables. Each row of X corresponds to one
%    input vector and each row of Y corresponds to one output vector.
%
%    [EFT, VARFT, LPYT] = GPSVI_PREDGRAD(GP, X, Y, XT, 'yt', YT, OPTIONS)
%    returns also logarithm of the predictive density LPYT of the
%    observations YT at test input locations XT. This can be used
%    for example in the cross-validation. Here Y has to be a vector.
%
%    [EFT, VARFT, LPYT, EYT, VARYT] = GPSVI_PREDGRAD(GP, X, Y, XT, OPTIONS)
%    returns also the posterior predictive mean EYT and variance VARYT.
%
%    [EF, VARF, LPY, EY, VARY] = GPSVI_PREDGRAD(GP, X, Y, OPTIONS)
%    evaluates the predictive distribution at training inputs X
%    and logarithm of the predictive density LPY of the training
%    observations Y.
%
%    OPTIONS is optional parameter-value pair
%      predcf - an index vector telling which covariance functions are
%               used for prediction. Default is all (1:gpcfn).
%               See additional information below.
%      yt     - optional observed yt in test points (see below)
%      z      - optional observed quantity in triplet (x_i,y_i,z_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, expected value
%               for ith case.
%      zt     - optional observed quantity in triplet (xt_i,yt_i,zt_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, the expected
%               value for the ith case.
%
%  See also
%    GPSVI_PRED, SVIGP, DEMO_SVI*
%

% Copyright (c) 2014 Ville Tolvanen

% 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.

ip=inputParser;
ip.FunctionName = 'GP_PRED';
ip.addRequired('gp',@(x) isstruct(x) || iscell(x));
ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addOptional('xt', [], @(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))))
ip.addParamValue('yt', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('zt', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('predcf', [], @(x) isempty(x) || ...
                 isvector(x) && isreal(x) && all(isfinite(x)&x>=0))
if numel(varargin)==0 || isnumeric(varargin{1})
  % inputParser should handle this, but it doesn't
  ip.parse(gp, x, y, varargin{:});
else
  ip.parse(gp, x, y, [], varargin{:});
end
xt=ip.Results.xt;
yt=ip.Results.yt;
zt=ip.Results.zt;
z=ip.Results.z;
predcf=ip.Results.predcf;
if isempty(xt)
  xt=x;
  if isempty(yt)
    yt=y;
  end
  if isempty(zt)
    zt=z;
  end
end

tn = size(x,1);
if nargout > 2 && isempty(yt)
  lpyt=[];
end

%     [tmp,tmp,tmp,param]=gpsvi_e(gp_pak(gp),gp,x,y);

% Check if the variational parameters has been set
if ~isfield(gp, 'm') || ~isfield(gp, 'S')
  error('Variational parameters has not been set. Call SVIGP first.')
end

u = gp.X_u;
m=gp.m;
S=gp.S;
if size(u,2) ~= size(x,2)
  % Turn the inducing vector on right direction
  u=u';
end
% Calculate some help matrices
K_uu = gp_trcov(gp, u);     % u x u, noiseles covariance K_uu
K_nu = gp_dcov2(gp, u, [], xt, xt);
K_nu = (K_nu(:,size(xt,1)+1:end))';
%     K_nu = gp_cov(gp,xt,u);       % n x u
[Luu, notpositivedefinite] = chol(K_uu,'lower');
if notpositivedefinite
  Eft=NaN; Varft=NaN; lpyt=NaN; Eyt=NaN; Varyt=NaN;
  return
end

Eft = K_nu*(Luu'\(Luu\m));

if nargout > 1
  [Knn_v, Cnn_v] = gp_trvar(gp,xt,predcf);
  B2=Luu\(K_nu');
  B3=K_uu\(K_nu');

  Varft = Knn_v - sum(B2.*B2)' + sum(B3.*(S*B3))';
end
s2=gp.lik.sigma2;

if nargout > 2
  if isequal(gp.lik.type, 'Gaussian')
    Eyt = Eft;
    Varyt = Varft + Cnn_v - Knn_v;
    if ~isempty(yt)
      lpyt = norm_lpdf(yt, Eyt, sqrt(Varyt));
    end
  else
    if nargout>3
      [lpyt, Eyt, Varyt] = gp.lik.fh.predy(gp.lik, Eft, Varft+s2, yt, zt);
    else
      lpyt = gp.lik.fh.predy(gp.lik, Eft, Varft+s2, yt, zt);
    end
  end
end
end

function [C, Cinv] = gp_dcov2(gp, x1, xv1, x2, xv2, predcf)
% Replacement for GP_DCOV where virtual inputs are given as arguments.

% Split the training data for normal latent input and gradient inputs
x12=x1;
%x11=gp.xv;
x11=xv1;
x3=xv2;

% Derivative observations
[n,m]=size(x1);
[n4,m4]=size(x2);
ncf=length(gp.cf);
if isfield(gp, 'nvd')
  % Only specific dimensions
  ii1=abs(gp.nvd);
else
  % All dimensions
  ii1=1:m;
end
for i1=1:ncf
  gpcf = gp.cf{i1};
  if m==1
    Gset2 = gpcf.fh.ginput4(gpcf, x3, x12);
    Kff = gpcf.fh.cov(gpcf, x12, x2);
    if ~isempty(x11)
      Gset1 = gpcf.fh.ginput4(gpcf, x11, x2);
      Kdd = gpcf.fh.ginput2(gpcf, x11, x3);
      Kdf=Gset1{1};
    else
      Kdd={[]};
      Kdf=[];
    end

    Kfd=Gset2{1};
    C = [Kff Kfd'; Kdf Kdd{1}];

    % Input dimension is >1
  else
    [n,m]=size(x11);
    [n2,m2]=size(x3);

    Kff = gpcf.fh.cov(gpcf, x12, x2);
    Gset2 = gpcf.fh.ginput4(gpcf, x3, x12);

    %Gather matrices from Gset (d k(x1,x2) /d x1)
    Kfd22=cat(2,Gset2{ii1});
    Kdf22=cat(1,Gset2{ii1})';

    if ~isempty(x11)
      Gset1 = gpcf.fh.ginput4(gpcf, x11,x2);
      Kfd=cat(2,Gset1{ii1});
      Kdf=cat(1,Gset1{ii1});
      % both x derivatives, same dimension (to diagonal blocks)
      D = gpcf.fh.ginput2(gpcf, x11, x3);
      % both x derivatives, different dimension (non-diagonal blocks)
      Kdf2 = gpcf.fh.ginput3(gpcf, x11 ,x3);

      % Now build up Kdd m*n x m*n2 matrix, which contains all the
      % both partial derivative" -matrices

      % Add the diagonal matrices
      Kdd=blkdiag(D{1:m});
      % Add the non-diagonal matrices to Kdd
      ii3=0;
      for j=0:m-2
        for i=1+j:m-1
          ii3=ii3+1;
          Kdd(i*n+1:(i+1)*n,j*n2+1:j*n2+n2) = Kdf2{ii3};
          Kdd(j*n+1:j*n+n,i*n2+1:(i+1)*n2) = Kdf2{ii3};
        end
      end
      if isfield(gp, 'nvd')
        % Collect the correct gradient dimensions, i.e. select the blocks
        % that correspond to the input dimensions for which we want the
        % gradients to be monotonic
        Kddtmp=[];
        for ii2=1:length(ii1)
          for ii3=ii2:length(ii1)
            Kddtmp((ii2-1)*n+1:ii2*n, (ii3-1)*n2+1:ii3*n2) = ...
              Kdd((ii1(ii2)-1)*n+1:ii1(ii2)*n,(ii1(ii3)-1)*n2+1:ii1(ii3)*n2);
            if ii2~=ii3
              Kddtmp((ii3-1)*n+1:ii3*n, (ii2-1)*n2+1:ii2*n2) = ...
                Kdd((ii1(ii3)-1)*n+1:ii1(ii3)*n,(ii1(ii2)-1)*n2+1:ii1(ii2)*n2);
            end
          end
        end
        Kdd=Kddtmp;
      end
    else
      Kdf=[];
      Kdd=[];
    end

    % Gather all the matrices into one final matrix K which is the
    % training covariance matrix
    C = [Kff Kdf22; Kdf Kdd];
    %   C = [Kff; Kdf];
  end
  if i1==1
    CC=C;
  else
    CC=CC+C;
  end
end
C=CC;
end

author: gpstuff-dev

repo: gpstuff

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


text:
%% housekeeping
clear
close all
clc()

%% Create dataset
clc

do_plot=true;

scale=100;

[db,varlist0]=create_dataset(scale,do_plot);

varlist=fieldnames(varlist0);
%% Choose a model type: see cell "create the structural VAR model" below

model_type=0;

%% set up the restrictions
close()

% create restrictions on parameters as well as markov chains
%------------------------------------------------------------
switch model_type
    case 0
        % constant-parameter model
        [restrictions,markov_chains,switch_prior]=create_restrictions_and_markov_chains0();
    case 1
        % Coefficients are switching regimes across all equations
        % (synchronized case)
        [restrictions,markov_chains,switch_prior]=create_restrictions_and_markov_chains1();
    case 2
        % Coefficients and variances have different chains, different
        % regimes, and different durations
        [restrictions,markov_chains,switch_prior]=create_restrictions_and_markov_chains2();
    case 3
        % Only coefficients in monetary policy equation are changing
        [restrictions,markov_chains,switch_prior]=create_restrictions_and_markov_chains3();
    case 4
        % Only variance in monetary policy equation is changing
        [restrictions,markov_chains,switch_prior]=create_restrictions_and_markov_chains4();
    case 5
        % Both coefficients and variances in monetary policy equation
        % change with two independent Markov processes
        [restrictions,markov_chains,switch_prior]=create_restrictions_and_markov_chains5();
    case 6 % ok
        % Only variances in ALL three equations switch
        [restrictions,markov_chains,switch_prior]=create_restrictions_and_markov_chains6();
    otherwise
        error('the coded model types are 0, 1, 2, 3, 4 and 6')
end

%% Create the VAR

clc

nlags=2;

exog={};

panel=[];

constant=true;

% first we create a template structure
% ------------------------------------
sv0=svar(varlist,exog,nlags,constant,panel,markov_chains);

%% set priors % prior=[];

var_prior=svar.prior_template();

var_prior.type='sz';

prior=struct('var',var_prior,'nonvar',switch_prior);

is_prior=true;

if ~is_prior

    prior=rmfield(prior,'var');

end

%% Find posterior mode
clc

sv=sv0;

sv=estimate(sv,db,{'1960Q1','2015Q2'},prior,restrictions);

%% estimates
clc

pmode=posterior_mode(sv)

%% Printing estimates
clc

print_structural_form(sv)

%% Printing solution
clc

print_solution(sv)

%% plot smoothed state and regime probabilities
clc
close all

plot_probabilities(sv)

%% plots probabilities against data
close all

plot_data_against_probabilities(sv,'regime')

%% Impulse responses

myirfs=irf(sv);

%% Posterior sampling

%% Marginal data density

%% Out-of sample forecasts

%% Conditional forecast

author: jmaih


text:
function m = rnsubset(k,n)
%RNSUBSET choose k distinct random integers from 1:n M=(K,N)
%
%  Inputs:
%
%    K is number of disinct integers required from the range 1:N
%    N specifies the range - we must have K<=N
%
%  Outputs:
%
%    M(1,K) contains the output numbers

%      Copyright (C) Mike Brookes 2006
%      Version: $Id: rnsubset.m,v 1.2 2007/05/04 07:01:39 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 2 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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if k>n
    error('rnsubset: k must be <= n');
end
% We use two algorithms according to the values of k and n
[f,e]=log2(n);
if k>0.03*n*(e-1)
[v,m]=sort(rand(1,n)); % for large k, just do a random permutation
else
    v=ceil(rand(1,k).*(n:-1:n-k+1));
    m=1:n;
    for i=1:k
        j=v(i)+i-1;
        x=m(i);
        m(i)=m(j);
        m(j)=x;
    end
end
m=m(1:k);

author: decouples


text:
%% Example of using KST class for interfacing with KUKA iiwa robots

% Used to test switching between various operation modes.

% First start the server on the KUKA iiwa controller
% Then run this script using Matlab

% Tested under Matlab2013b, for higher versions, the realtime plot might
% not work.

% This example works with Sunrise application version KST_1.7  and higher.
% Copyright Mohammad SAFEEA, 1st-March-2019

% Operation modes:
% 1st impedance motion
% 2nd (PTP) move on lines
% 3rd (PT) move on circles
% 4th direct servo
% 5th impedance again
% 6th (PTP) move lines again
close all;clear;clc;
warning('off');
%% Create the robot object
ip='172.31.1.147'; % The IP of the controller
arg1=KST.LBR7R800; % choose the robot iiwa7R800 or iiwa14R820
arg2=KST.Medien_Flansch_elektrisch; % choose the type of flange
Tef_flange=eye(4); % transofrm matrix of EEF with respect to flange
iiwa=KST(ip,arg1,arg2,Tef_flange); % create the object

%% Start a connection with the server
flag=iiwa.net_establishConnection();
if flag==0
return;
end
pause(1);
disp('Testing KST under various operation modes');

%% Move point to point to an initial position
jPos={0,0,0,-pi/2,0,pi/2,0};
relVel=0.15;
iiwa.movePTPJointSpace(jPos, relVel); % move to initial configuration

%% Pause for 3 seocnds
pause(3);

%% 1st
disp('First test, real time impedance mode');
% Tool/Impedance parameters
massOfTool=0.5; % the mass of the tool attached to flange in Kg
cOMx=0; % X coordinate of the center of mass of the tool in (mm)
cOMy=0; % Y coordinate of the center of mass of the tool in (mm)
cOMz=40; % Z coordinate of the center of mass of the tool in (mm)
cStiness=900; % cartizian stifness
rStifness=80; % rotational stifness
nStifness=50; % null space stifness

% Start the soft realtime control with impedance
iiwa.realTime_startImpedanceJoints(massOfTool,cOMx,cOMy,cOMz,...
cStiness,rStifness,nStifness);

w=0.6; % motion constants, frequency rad/sec
A=0.2; % motion constants, amplitude of motion

a=datevec(now);
t0=a(6)+a(5)*60+a(4)*60*60; % calculate initial time

dt=0;
tstart=t0;
counter=0;
duration=0.2*60; %0.2 minutes

% Control loop
% real time plot handles
 colors={'k','b','r','g','c','m','y'};
 figureHandle=figure('Units','inches','Position',[0 0 5 3.75]);
 numOfSamples=120; % number of samples to show in the plot
 timeVec=1:numOfSamples;
 tawVec=zeros(7,numOfSamples);
 plotHandle=[];
 for i=1:7
    plotHandle=[plotHandle,plot(timeVec,tawVec(i,:),colors{i},'LineWidth',2)];
    hold on;
 end
% format the plot
ylim([0,8]);
temp=title('External torques');
set(temp,'FontSize',16);
temp=xlabel('Time (seconds)');
set(temp,'FontSize',14);
temp=ylabel('Distance (m)');
set(temp,'FontSize',14);

while(dt<duration)
% perform trajectory calculation here
  a=datevec(now);
  time=a(6)+a(5)*60+a(4)*60*60;
  dt=time-t0;

  temp=A*(1-cos(w*dt));
    for dacount=1:7
        jPosCommand{dacount}=jPos{dacount}+temp;
    end
  counter=counter+1;
  % Send joint positions to robot
  taw=iiwa.sendJointsPositionsExTorque(jPosCommand);

  tawTemp=zeros(7,1);
  for i=1:7
      tawTemp(i)=taw{i};
      torque_array(i,counter)=taw{i};
  end
  % delete first measurment
  tawVec(:,1)=[];
  % add new measurment to end of arrays
  tawVec=[tawVec,tawTemp];
  % update plots data
  for i=1:6
    set(plotHandle(i),'YData',tawVec(i,:));
  end
  set(plotHandle(7),'YData',tawVec(i,:));
  % show plot data
  set(figureHandle,'Visible','on');

end

% Stop the realtime control with impedence motion
iiwa.realTime_stopImpedanceJoints();
fprintf('\n impedance motion stopped \n');
pause(2);

%% 2nd
disp('Second test, point to point motion mode, move on lines');
% move to initial position
pinit={0,pi*20/180,0,-pi*70/180,0,pi*90/180,0}; % initial confuguration
relVel=0.15; % relative velocity
iiwa.movePTPJointSpace(pinit, relVel); % point to point motion in joint space

% Get position orientation of end-effector
disp('Cartesian position')
Pos=iiwa.getEEFPos();
disp(Pos)

% Move the end-effector in the X direction
Vel=50; % velocity of the end effector, mm/sec
displacment=50;
index=1;

Pos{index}=Pos{index}+displacment;
iiwa.movePTPLineEEF(Pos, Vel)

pause(0.1)

Pos{index}=Pos{index}-displacment;
iiwa.movePTPLineEEF(Pos, Vel)

% Move the end-effector in the z direction
index=3;

Pos{index}=Pos{index}+displacment;
iiwa.movePTPLineEEF(Pos, Vel)

pause(0.1)

Pos{index}=Pos{index}-displacment;
iiwa.movePTPLineEEF(Pos, Vel);

% Move the end-effector to a distination frame,

Pos={400,0,580,-pi,0,-pi}; % first configuration
% the coordinates are:
% x=400mm, y=0mm, z=580 mm.
% the rotation angles are:
% alpha=-180 degrees, beta=0 degrees, gama=-180 degrees
iiwa.movePTPLineEEF(Pos, Vel);

Pos={500,0,580,-pi,0,-pi+pi/4}; %% second configuration
% the coordinates are:
% x=400mm, y=0mm, z=580 mm.
% the rotation angles are:
% alpha=-180 degrees, beta=0 degrees, gama=-165 degrees,
iiwa.movePTPLineEEF(Pos, Vel);

%% 3rd
disp('Third test, point to point motion mode, move on circles');
% move on cirlces
relVel=0.25; % override relative joint velocities

pos={0, -pi / 180 * 10, 0, -pi / 180 * 100, pi / 180 * 90,pi / 180 * 90, 0};   % initial cofiguration

iiwa.movePTPJointSpace( pos, relVel); % go to initial position
% Move in an arc, the orientation of EEF changes while performing the motion,
% The function utilized (movePTPCirc1OrintationInter)
% f2 is the final frame, at which the arc motion ends
% f1 is an intermidiate frame, through wich the robot passes while
% performing the motion.

f1=iiwa.getEEFPos();
f2=f1;
r=75; % radius of arc
f1{2}=f1{2}+r;
f1{3}=f1{3}-r;
f1{6}=f1{6}+pi/8;

f2{3}=f2{3}-2*r;
f2{6}=f2{6}+pi/2;

vel=150; % linear velocity of end-effector mm/sec
iiwa.movePTPCirc1OrintationInter( f1,f2, vel)

% Move robot in joint space to some initial configuration
pinit={0,pi*20/180,0,-pi*70/180,0,pi*90/180,0}; % joint angles of initial confuguration
relVel=0.15; % the relative velocity
iiwa.movePTPJointSpace(pinit, relVel); % point to point motion in joint space

% Move EEF -100 mm in Z direction
deltaX=0.0;deltaY=0;deltaZ=-100.; % relative displacemnets of end-effector
Pos{1}=deltaX;
Pos{2}=deltaY;
Pos{3}=deltaZ;
iiwa.movePTPLineEefRelBase(Pos, vel);

% Store the current position in the memory
Cen=iiwa.getEEFPos(); % Concider the current position as the center of the arcs

% Move EEF 50mm in X direction
deltaX=100;deltaY=0;deltaZ=0.;
Pos{1}=deltaX;
Pos{2}=deltaY;
Pos{3}=deltaZ;
iiwa.movePTPLineEefRelBase(Pos, vel);
% Store the current position in the memory
circle_Starting_Point=iiwa.getEEFPos(); % Consider the current point as circle starting point

% Move in an arc, the arc is drawn on an incliend plane
% using the function ((movePTPArc_AC))
theta=pi/2; % the angle subtended by the arc at the center ((c))
k=[1;1;1]; % normal vector of the plane, on which the circle is drawn
c=[Cen{1};Cen{2};Cen{3}]; % the center of the arc
vel=100; % the motion velocity mm/sec
iiwa.movePTPArc_AC(theta,c,k,vel)

% Go back to ((circle_Starting_Point)) coordinates
vel=150;
iiwa.movePTPLineEEF(circle_Starting_Point, vel);
% Move in an arc, the arc is drawn in XY plane
% using the function ((movePTPArcXY_AC))
theta=1.98*pi; % the angle subtended by the arc at the center ((c))
c=[Cen{1};Cen{2}]; % the XY coordinate of the center of the arc
vel=150; % the motion velocity mm/sec
iiwa.movePTPArcXY_AC(theta,c,vel)

%% 4th
disp('Fourth test, realtime motion using direct servo mode');
% Move to an intial position
jPos={0,0,0,-pi/2,0,pi/2,0};
relVel=0.15;
iiwa.movePTPJointSpace(jPos, relVel); % move to initial configuration

% Get Cartesian position of EEF
fprintf('Cartesian position')
eefpos=iiwa.getEEFPos();
eefposDist=eefpos;
disp(eefpos)

% Start direct servo in Cartesian space
iiwa.realTime_startDirectServoCartesian();
disp('Starting direct servo in Cartesian space');
w=2; % motion constants, frequency rad/sec
A=75; % motion constants, amplitude of motion (mm)
tic;
deltaT=0;
counter=0;
initiationFlag=0;
disp('Entering control loop, streaming EEF positions')

% Control loop
while(deltaT<(6*pi/w))
if(initiationFlag==0)
    initiationFlag=1;
    t_0=toc;
    t0=t_0;
else
    time=toc;
    deltaT=time-t0;
    %........................................
    % Perform trajectory calculation here
    eefposDist{3}=eefpos{3}-A*(1-cos(w*deltaT));
    %........................................
    % Send EEF position to robot
    if(toc-t_0>0.003)
        counter=counter+1;
        iiwa.sendEEfPositionf(eefposDist);
        t_0=toc;
    end
end
end
tstart=t0;
tend=time;
rate=counter/(tend-tstart);
% Stop the direct servo motion
iiwa.realTime_stopDirectServoCartesian( );
fprintf('\nThe rate of update per second is: \n');
disp(rate);
fprintf('\n')
pause(2);

%% 5th
disp('Fifth test, realtime impedance motion mode, repeat procedure for 3 times');
for i=1:3
    massOfTool=0.5; % the mass of the tool attached to flange in Kg
    cOMx=0; % X coordinate of the center of mass of the tool in (mm)
    cOMy=0; % Y coordinate of the center of mass of the tool in (mm)
    cOMz=40; % Z coordinate of the center of mass of the tool in (mm)
    cStiness=900; % cartizian stifness
    rStifness=80; % rotational stifness
    nStifness=50; % null space stifness

    % Start the soft realtime control with impedance
    iiwa.realTime_startImpedanceJoints(massOfTool,cOMx,cOMy,cOMz,...
    cStiness,rStifness,nStifness);

    w=0.6; % motion constants, frequency rad/sec
    A=0.2; % motion constants, amplitude of motion

    a=datevec(now);
    t0=a(6)+a(5)*60+a(4)*60*60; % calculate initial time

    dt=0;
    tstart=t0;
    counter=0;
    duration=10; %10 seconds

    % Control loop
    while(dt<duration)
    % perform trajectory calculation here
      a=datevec(now);
      time=a(6)+a(5)*60+a(4)*60*60;
      dt=time-t0;

      temp=A*(1-cos(w*dt));
        for dacount=1:7
            jPosCommand{dacount}=jPos{dacount}+temp;
        end
      counter=counter+1;
      % Send joint positions to robot
      iiwa.sendJointsPositions(jPosCommand);
    end
    % Stop the realtime control with impedence motion
    iiwa.realTime_stopImpedanceJoints();
    fprintf('\n Impedance motion stopped \n');
    pause(2);

    daPos_temp=iiwa.getEEFPos();
end

%% 6th
disp('Sexth test, repeated from point to point motion mode, move on lines');
% move to initial position
pinit={0,pi*20/180,0,-pi*70/180,0,pi*90/180,0}; % initial confuguration
relVel=0.15; % relative velocity
iiwa.movePTPJointSpace(pinit, relVel); % point to point motion in joint space

% Get position orientation of end effector
disp('Cartesian position')
Pos=iiwa.getEEFPos();
disp(Pos)

% Move the end-effector in the X direction
Vel=50; % velocity of the end effector, mm/sec
displacment=50;
index=1;

Pos{index}=Pos{index}+displacment;
iiwa.movePTPLineEEF(Pos, Vel)

pause(0.1)

Pos{index}=Pos{index}-displacment;
iiwa.movePTPLineEEF(Pos, Vel)

% Move the end-effector in the z direction
index=3;

Pos{index}=Pos{index}+displacment;
iiwa.movePTPLineEEF(Pos, Vel)

pause(0.1)

Pos{index}=Pos{index}-displacment;
iiwa.movePTPLineEEF(Pos, Vel);

% Move the end-effector to a distination frame,

Pos={400,0,580,-pi,0,-pi}; % first configuration
% the coordinates are:
% x=400mm, y=0mm, z=580 mm.
% the rotation angles are:
% alpha=-180 degrees, beta=0 degrees, gama=-180 degrees
iiwa.movePTPLineEEF(Pos, Vel);

Pos={500,0,580,-pi,0,-pi+pi/4}; %% second configuration
% the coordinates are:
% x=400mm, y=0mm, z=580 mm.
% the rotation angles are:
% alpha=-180 degrees, beta=0 degrees, gama=-165 degrees,
iiwa.movePTPLineEEF(Pos, Vel);

%% Turn off the Sunrise.Workbench application
iiwa.net_turnOffServer

repo: KST-Kuka-Sunrise-Toolbox


text:
function Tr_Im=rotate3D(Im, T, R)
% rotate 3D volume by T & R

[row,col,layer] = size(Im);
Tr_Im = tformarray(Im, T, R, [1 2 3], [1 2 3], [row, col,layer], [], []);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function x=randiscr(p,n,a)
%RANDISCR Generate discrete random numbers with specified probabiities [X]=(P,N,A)
%
% Usage: (1) randiscr([],10)        % generate 10 uniform random binary values
%        (2) randiscr(2:6,10)       % generate 10 random numbers in the range 1:5
%                                     with probabilities [2 3 4 5 6]/20
%        (3) randiscr([],10,'abcd') % generate a string of 10 random
%                                     characters equiprobable from 'abcd'
%
% Inputs: P  vector of probabilities (not necessarily normalized) [default = uniform]
%         N  number of random values to generate [default = 1]
%         A  output alphabet [default = 1:length(p) or 0:1 if p is empty]
%
% Outputs: X  vector of not necessarily distinct values taken from alphabet A
%
% The vector P is internally normalized by dividing by its sum.
% If P is an M-dimensional matrix (and A is unspecified), then X will
% have dimensions (N,M) with the corresponding indices for each dimension.

% Somewhat similar in function to RANDSRC in the comms toolbox

%   Copyright (c) 2005-2012 Mike Brookes,  mike.brookes@ic.ac.uk
%      Version: $Id: randiscr.m 2189 2012-07-20 13:47:00Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 2 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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gota=nargin>2;
if nargin<1 || ~numel(p)
    if gota
        p=ones(1,length(a));
    else
        p=ones(1,2);
        a=(0:1)';
        gota=1;
    end
end
if nargin<2 || ~numel(n)
    n=1;
end
d=length(p(:)); % size of output alphabet
z=zeros(d+n-1,1); % array to hold random numbers
z(1:d)=cumsum(p(:)/sum(p(:))); % last value is actually overwritten in the next line
z(d:d+n-1)=rand(n,1);
[y,iy]=sort(z);
y(iy)=(1:d+n-1)';
m=zeros(d+n-1,1);
m(y(1:d-1))=1;
m(1)=m(1)+1;
mc=cumsum(m);
x=mc(y(d:d+n-1));
if gota
    x=a(x);
elseif length(p(:))>length(p) % need multiple dimensions
    v=x-1;
    s=cumprod(size(p));
    m=length(s);
    s(2:end)=s(1:end-1);
    s(1)=1;
    x=zeros(n,m);
    for i=m:-1:1
        x(:,i)=1+floor(v/s(i));
        v=rem(v,s(i));
    end
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function [mu ul ll] = circ_mean(alpha, w, dim)
%
% mu = circ_mean(alpha, w)
%   Computes the mean direction for circular data.
%
%   Input:
%     alpha	sample of angles in radians
%     [w		weightings in case of binned angle data]
%     [dim  compute along this dimension, default is 1]
%
%     If dim argument is specified, all other optional arguments can be
%     left empty: circ_mean(alpha, [], dim)
%
%   Output:
%     mu		mean direction
%     ul    upper 95% confidence limit
%     ll    lower 95% confidence limit
%
% PHB 7/6/2008
%
% References:
%   Statistical analysis of circular data, N. I. Fisher
%   Topics in circular statistics, S. R. Jammalamadaka et al.
%   Biostatistical Analysis, J. H. Zar
%
% Circular Statistics Toolbox for Matlab

% By Philipp Berens, 2009
% berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html

if nargin < 3
  dim = 1;
end

if nargin < 2 || isempty(w)
  % if no specific weighting has been specified
  % assume no binning has taken place
	w = ones(size(alpha));
else
  if size(w,2) ~= size(alpha,2) || size(w,1) ~= size(alpha,1)
    error('Input dimensions do not match');
  end
end

% compute weighted sum of cos and sin of angles
r = nansum(w.*exp(1i*alpha),dim);

% obtain mean by
mu = angle(r);

% confidence limits if desired
if nargout > 1
  t = circ_confmean(alpha,0.05,w,[],dim);
  ul = mu + t;
  ll = mu - t;
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [a]=add(a,d)

% function [a]=add(a,d)
%
% add a single algorithm or a cell array algorithms to an algs
% object

  d=make_cell(d);
  r=a.child;
  if ~isempty(r) r={r{1:length(r)} d{1:length(d)}}; else r=d; end
  a.child=r;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [varargout, peval]=updates_variational(peval, varargin)
% [varargout, peval]=variationalupdates(peval, varargin)
% Computes variational updates
% dvec = varargin{1};   % data (N x T)
% w = varargin{2};      % initialization for basis (N x peval.ncomp)
% h = varargin{3};      % initialization for a (peval.ncomp x T) and b
% a = h.a;      % initialization for a
% b = h.b;
mfprintf(peval.fid, 'Variational updates.\n')
peval=setDefaultValuesPeval(peval);
dvec = varargin{1};
w = varargin{2};        % initialization for basis (N x peval.ncomp)
h = varargin{3};        % initialization for a (peval.ncomp x T) and b
a = h.a;                % initialization for a
b = h.b;

alph=repmat(peval.alpha,1,peval.ncomp); % prior parametes for Gamma distribution (peval.ncomp x T)
beta=repmat(peval.beta,1,peval.ncomp);  % prior parametes for Gamma distribution (peval.ncomp x T)

[dovec_w, dovec_a]=setDoVec(peval); %sets which components used in updates
lb=zeros(1,peval.maxiter);
lb(1)=sum(lowerbound(dvec, w, alph, beta, a, b));
for updateindex=1:peval.maxiter;
    n=update_n(w, a, b);
    a(dovec_a,:)=update_a(dvec,n(:,dovec_a,:),alph(:,dovec_a));
    w(:,dovec_w)=update_w(dvec, n(:,dovec_w,:));
    lb(updateindex+1)=sum(lowerbound(dvec, w, alph, beta, a, b));
    if peval.showprogress
        plotprogress(w,lb(2:updateindex),peval)
    end
end

varargout = struct('w',w,'a',a,'lb',lb(end));

end % Main function

% Nested fucntions:
function peval=setDefaultValuesPeval(peval)
if ~isfield(peval,'showprogress')
    peval.showprogress = 0;
end

if ~isfield(peval, 'fix_bg_w')
    peval.fix_bg_w=0;
    if peval.addbgcomp
        peval.fix_bg_w=1;
    end
end
if ~isfield(peval, 'fix_bg_a')
    peval.fix_bg_a=0;
    if peval.addbgcomp
        peval.fix_bg_a=1;
    end
end
if peval.fix_bg_w
    mfprintf(peval.fid, '''w'' fixed for background component [%g].\n',peval.ncomp)
end
if peval.fix_bg_a
    mfprintf(peval.fid, '''a'' fixed for background component [%g].\n',peval.ncomp)
end
end

function [dovec_w, dovec_a]=setDoVec(peval)
dovec_w=1:peval.ncomp;
dovec_a=1:peval.ncomp;
if peval.fix_bg_w
    dovec_w=1:peval.ncomp-1;
end
if peval.fix_bg_a
    dovec_a=1:peval.ncomp-1;
end
if isfield(peval, 'fix_w')
    dovec_w = removerows(dovec_w',peval.fix_w)';
end
if isfield(peval, 'fix_a')
    dovec_a = removerows(dovec_a',peval.fix_a)';
end
mfprintf(peval.fid, '''w'' will be updated for components: ')
mfprintf(peval.fid, '['); mfprintf (peval.fid, '%g ', dovec_w); mfprintf (peval.fid, ']\n');
mfprintf(peval.fid, '''a'' will be updated for components: ')
mfprintf(peval.fid, '['); mfprintf (peval.fid, '%g ', dovec_a); mfprintf (peval.fid, ']\n');
end

function plotprogress(w,lb,peval)
imstiled(reshape(w, peval.nx, peval.ny, peval.ncomp),10, 'gray',[],[],1)
plot(lb)
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% y-position of nose
function [data,units] = compute_ynose_mm(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);
  data{i} = trx(fly).y_mm + 2*trx(fly).a_mm.*sin(trx(fly).theta_mm);
end
units = parseunits('mm');

author: kristinbranson

repo: JAABA


text:
function vis_people( expidx, firstidx, nImgs)

if (nargin < 3)
    nImgs = 1;
end

colors = {'r','g','b','c','m','y'};
markerSize = 6;
lineWidth = 4;
edges = [1 2; 2 3; 3 13; 4 13; 4 5; 5 6; 7 8; 8 9; 10 11; 11 12; 9 13; 10 13; 13 14];

p = exp_params(expidx);
load(p.testGT,'annolist');

lastidx = firstidx+nImgs-1;

figure;

for imgidx = firstidx:lastidx

    imgname = annolist(imgidx).image.name;
    img = imread(imgname);

    clf;
    imagesc(img); axis equal; hold on;

    load([p.multicutDir '/prediction_' padZeros(num2str(imgidx),4)], 'people');

    for j = 1:length(people)
        person = people{j};
        person_color = colors{mod(j-1,6)+1};

        for i = 1:size(edges,1)
            pos1 = person(edges(i,1), :);
            pos2 = person(edges(i,2), :);
            if (~isnan(pos1(1)) && ~isnan(pos2(1)))
                plot([pos1(1);pos2(1)],[pos1(2);pos2(2)],[person_color '-'],'linewidth',lineWidth);
            end
        end

        for i = 1:size(person, 1)
            if isnan(person(i, 1))
                continue;
            end
            pos = person(i, :);
            cp = colors{mod(i-1,6)+1};
            plot(pos(:,1),pos(:,2),[cp 'o'],'MarkerSize',markerSize,'MarkerFaceColor',cp,'MarkerEdgeColor','k');
        end

    end
end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [] = plot_skeleton(pts, A, options)

if nargin<3
    options.sizep=400;
end
sizep = getoptions(options, 'sizep', 400);
sizee = getoptions(options, 'sizee', 2);
colorp = getoptions(options, 'colorp', [1, .8, .8]);
colore = getoptions(options, 'colore', [1, .0, .0]);

scatter3(pts(:,1),pts(:,2),pts(:,3),sizep,'.','MarkerEdgeColor', colorp);  hold on;
plot_connectivity(pts, A, sizee, colore);

author: taiya

repo: cloudcontr

save_path: github-repos/MATLAB/taiya-cloudcontr


text:
function CPad = directFilter(A,D,tensorOrder,lambda)
% bsarray\directFilter: compute B-spline coefficients from data array
% usage: C = directFilter(A,D,tensorOrder,lambda);
%
% arguments:
%   A - N-dimensional array (vector, image, volume, etc.)
%   D (1xN vector) - degree of each dimension of tensor product BSpline.
%   tensorOrder (scalar) - number of nonsingleton dimensions of A
%   lambda (1xN vector) - smoothing factor for each dimension
%
%   C - array of BSpline coefficients
%

% author: Nathan D. Cahill
% email: ndcahill@gmail.com
% date: 18 April 2008

% if D = 0 or 1, and lambda = 0, return A, with reflection padding
if all(D<2) && all(lambda==0)
    padNum = floor(D/2);
    idx   = cell(1,tensorOrder);
    for k = 1:tensorOrder
        M = size(A,k);
        dimNums = [1:M (M-1):-1:2];
        p = padNum(k);
        idx{k}   = dimNums(mod(-p:M+p, 2*M-2) + 1);
    end
    CPad = A(idx{:});
    return;
end

% get coefficients for BSpline filters
F = cell(tensorOrder,1);
F0 = cell(tensorOrder,1);
for i=1:tensorOrder
    % NEED TO UPDATE THIS TO SCALE LAMBDA FOR ELEMENT SPACING
    [F{i},F0{i}] = getBSplineFiltCoeffs(D(i),true,'direct',lambda(i));
end

% initialize output
C = double(A);

% choose tolerance for K0
K0Tol = eps;

switch tensorOrder
    case 1 % X is a vector - easy case

        CLen = numel(C);

        % loop through poles of direct BSpline filter
        for i = 1:length(F{1})

            % compute K0 for current pole
            K0 = ceil(log(K0Tol)./log(abs(F{1}(i))));

            % now compute vector of indices of length K0 that runs back and
            % forth through the length of the data vector (to mimic
            % reflection of the data at each end)
            indReflect = [1:CLen (CLen-1):-1:2];
            numReps = ceil(K0/(2*CLen-2));
            KVec = repmat(indReflect,[1 numReps]);
            KVec = KVec(1:K0);

            % compute scaling factor for current pole
            C0 = -F{1}(i)./(1-F{1}(i)^2);

            % apply symmetric exponential filter
            C(:) = symExpFilt(C(:),CLen,C0,F{1}(i),K0,KVec);

        end

        % multiply by numerator of direct BSpline filter
        C = real(C)*F0{1};

    case 2 % X is a matrix - a little more difficult

        % first perform direct filtering over each column
        [CRows,CCols] = size(C);

        % loop through poles of direct BSpline filter
        for i = 1:length(F{1})

            % compute K0 for current pole
            K0 = ceil(log(K0Tol)./log(abs(F{1}(i))));

            % now compute vector of indices of length K0 that runs back and
            % forth through the length of the data vector (to mimic
            % reflection of the data at each end)
            indReflect = [1:CRows (CRows-1):-1:2];
            numReps = ceil(K0/(2*CRows-2));
            KVec = repmat(indReflect,[1 numReps]);
            KVec = KVec(1:K0);

            % compute scaling factor for current pole
            C0 = -F{1}(i)./(1-F{1}(i)^2);

            % apply symmetric exponential filter for each column
            for k = 1:CCols
                C(:,k) = symExpFilt(C(:,k),CRows,C0,F{1}(i),K0,KVec);
            end
        end

        % multiply by numerator of direct BSpline filter
        C = real(C)*F0{1};

        % now perform direct filtering across each row
        % loop through poles of direct BSpline filter
        for i = 1:length(F{2})

            % compute K0 for current pole
            K0 = ceil(log(K0Tol)./log(abs(F{2}(i))));

            % now compute vector of indices of length K0 that runs back and
            % forth through the length of the data vector (to mimic
            % reflection of the data at each end)
            indReflect = [1:CCols (CCols-1):-1:2];
            numReps = ceil(K0/(2*CCols-2));
            KVec = repmat(indReflect,[1 numReps]);
            KVec = KVec(1:K0);

            % compute scaling factor for current pole
            C0 = -F{2}(i)./(1-F{2}(i)^2);

            % apply symmetric exponential filter for each column
            for k = 1:CRows
                C(k,:) = symExpFilt(C(k,:),CCols,C0,F{2}(i),K0,KVec);
            end
        end

        % multiply by numerator of direct BSpline filter
        C = real(C)*F0{2};

    case 3 % X is a volume - a bit more difficult

        % first perform direct filtering over each column
        [CRows,CCols,CSlices] = size(C);

        % loop through poles of direct BSpline filter
        for i = 1:length(F{1})

            % compute K0 for current pole
            K0 = ceil(log(K0Tol)./log(abs(F{1}(i))));

            % now compute vector of indices of length K0 that runs back and
            % forth through the length of the data vector (to mimic
            % reflection of the data at each end)
            indReflect = [1:CRows (CRows-1):-1:2];
            numReps = ceil(K0/(2*CRows-2));
            KVec = repmat(indReflect,[1 numReps]);
            KVec = KVec(1:K0);

            % compute scaling factor for current pole
            C0 = -F{1}(i)./(1-F{1}(i)^2);

            % apply symmetric exponential filter for each column
            for k = 1:CSlices
                for j = 1:CCols
                    C(:,j,k) = symExpFilt(C(:,j,k),CRows,C0,F{1}(i),K0,KVec);
                end
            end
        end

        % multiply by numerator of direct BSpline filter
        C = real(C)*F0{1};

        % now perform direct filtering across each row
        % loop through poles of direct BSpline filter
        for i = 1:length(F{2})

            % compute K0 for current pole
            K0 = ceil(log(K0Tol)./log(abs(F{2}(i))));

            % now compute vector of indices of length K0 that runs back and
            % forth through the length of the data vector (to mimic
            % reflection of the data at each end)
            indReflect = [1:CCols (CCols-1):-1:2];
            numReps = ceil(K0/(2*CCols-2));
            KVec = repmat(indReflect,[1 numReps]);
            KVec = KVec(1:K0);

            % compute scaling factor for current pole
            C0 = -F{2}(i)./(1-F{2}(i)^2);

            % apply symmetric exponential filter for each column
            for k = 1:CSlices
                for j = 1:CRows
                    C(j,:,k) = symExpFilt(C(j,:,k),CCols,C0,F{2}(i),K0,KVec);
                end
            end
        end

        % multiply by numerator of direct BSpline filter
        C = real(C)*F0{2};

        % now perform direct filtering across each slice
        % loop through poles of direct BSpline filter
        for i = 1:length(F{3})

            % compute K0 for current pole
            K0 = ceil(log(K0Tol)./log(abs(F{3}(i))));

            % now compute vector of indices of length K0 that runs back and
            % forth through the length of the data vector (to mimic
            % reflection of the data at each end)
            indReflect = [1:CSlices (CSlices-1):-1:2];
            numReps = ceil(K0/(2*CSlices-2));
            KVec = repmat(indReflect,[1 numReps]);
            KVec = KVec(1:K0);

            % compute scaling factor for current pole
            C0 = -F{3}(i)./(1-F{3}(i)^2);

            % apply symmetric exponential filter for each column
            for k = 1:CCols
                for j = 1:CRows
                    C(j,k,:) = symExpFilt(C(j,k,:),CSlices,C0,F{3}(i),K0,KVec);
                end
            end
        end

        % multiply by numerator of direct BSpline filter
        C = real(C)*F0{3};

    otherwise % X is multidimensional - the most difficult

        % get size of coefficients array
        Csz = size(C);

        % set up vector of indices into each dimension
        idx = cell(1,tensorOrder);
        for k = 1:tensorOrder
            idx{k} = 1:size(C,k);
        end

        % loop through each dimension
        for d = 1:tensorOrder

            % copy original version of idx
            idxCurrent = idx;

            % construct vector of indices which will be used in computing
            % the initial element of the filtered data
            indReflect = [1:Csz(d) (Csz(d)-1):-1:2];

            % construct list of dimensions without d
            dimList = setdiff(1:tensorOrder,d);

            % construct index into entries of remaining dimensions
            [idxCurrent{dimList}] = ndgrid(idx{dimList});
            dimListInd = sub2ind(Csz(dimList),idxCurrent{dimList});
            dimListInd = dimListInd(:);

            % loop through poles of BSpline filter
            for i = 1:length(F{d})

                % compute K0 for current pole
                K0 = ceil(log(K0Tol)./log(abs(F{d}(i))));

                % now compute vector of indices of length K0 that runs back and
                % forth through the length of the data vector (to mimic
                % reflection of the data at each end)
                numReps = ceil(K0/(2*Csz(d)-2));
                KVec = repmat(indReflect,[1 numReps]);
                KVec = KVec(1:K0);

                % compute scaling factor for current pole
                C0 = -F{d}(i)./(1-F{d}(i)^2);

                % now loop through remaining dimensions, applying symmetric
                % exponential filter
                for k=1:numel(dimListInd)
                    [idxCurrent{dimList}] = ind2sub(Csz(dimList),dimListInd(k));
                    C(idxCurrent{:}) = symExpFilt(C(idxCurrent{:}),Csz(d),C0,F{d}(i),K0,KVec);
                end
            end

            % multiply by numerator of direct BSpline filter
            C = real(C)*F0{d};

        end

end

% now pad dimensions by reflection
padNum = floor(D/2);

% Form index vectors to subsasgn input array into output array.
% Also compute the size of the output array.
idx   = cell(1,tensorOrder);
for k = 1:tensorOrder
  M = size(C,k);
  dimNums = [1:M (M-1):-1:2];
  p = padNum(k);
  idx{k}   = dimNums(mod(-p:M+p, 2*M-2) + 1);
end
CPad = C(idx{:});


text:
function LPproblem = liftModel(model, BIG, printLevel,fileName,directory)
% Lifts a COBRA model with badly-scaled stoichiometric and
% coupling constraints of the form:
% :math:`max c*v`  subject to: :math:`Sv = 0, x, Cv <= 0`
% Converts it into a COBRA LPproblem structure, which can be used with
% solveCobraLP. Fluxes for the reactions should stay the same i.e.
% sol.full(1:nRxns) should yield an optimal flux vector.
%
% USAGE:
%
%    LPproblem = liftModel(model, BIG, printLevel,fileName,directory)
%
% INPUTS:
%    model:     COBRA LPproblem Structure containing the original LP to be solved. The format of
%                   this struct is described in the documentation for `solveCobraLP.m`
%
% OPTIONAL INPUTS:
%    BIG:           A parameter the controls the largest entries that appear in the
%                   reformulated problem (default = 1000).
%    printLevel:    printLevel = 1 enables printing of problem statistics (default);
%                   printLevel = 0 silent
%    fileName:      name of th file to load
%    directory:     file directory (if `model` is empty, you can load it using `fileName` and `directory`)
%
%
% OUTPUTS:
%    LPproblem:         COBRA Structure contain the reformulated LP to be solved.
%
% .. Authors:
%       - Michael Saunders, saunders@stanford.edu
%       - Yuekai Sun, yuekai@stanford.edu, Systems Optimization Lab (SOL), Stanford University
%       - Ronan Fleming   (updated interface to take COBRA model structure)
%       - Thomas Pfau - Updated information, that the lifted problem is
%                       converted to a COBRA LP structure

if ~exist('BIG','var')
    BIG=1000;
end
if ~exist('printLevel','var')
    printLevel=1;
end
if exist('fileName','var') && exist('directory','var') && isempty(model)
    model = loadIdentifiedModel(fileName,directory);
end

%save original model
LPproblem = buildLPproblemFromModel(model);

% Assume constraint matrix is S if no A provided.
if ~isfield(LPproblem,'A')
    if isfield(LPproblem,'S')
        LPproblem.A = LPproblem.S;
    end
end

% Assume constraint S*v = b if csense not provided
if ~isfield(LPproblem,'csense')
    % If csense is not declared in the model, assume that all
    % constraints are equalities.
    LPproblem.csense(:,1) = 'E';
end

% Assume constraint S*v = 0 if b not provided
if ~isfield(LPproblem,'b')
    warning('LP problem has no defined b in S*v=b. b should be defined, for now we assume b=0')
    LPproblem.b=zeros(size(LPproblem.A,1),1);
end

% Assume max c'v s.t. S v = b if osense not provided
if ~isfield(LPproblem,'osense')
    LPproblem.osense = -1;
end

%call the LP reformulate script by Michael and Yuekai
LPproblem = reformulate(LPproblem, BIG, printLevel);

if exist('fileName','var') && exist('directory','var')
    save([directory filesep 'L_' fileName '.mat'],'LPproblem');
end

author: opencobra

repo: cobratoolbox


text:
function cheval_segmentation
% This code was written by Muhammet Balcilar , France,
% muhammetbalcilar@gmail.com, inspired from following reference
%
%  This Matlab code implements an edge-based active contour model as an
%  application of the Distance Regularized Level Set Evolution (DRLSE) formulation in Li et al's paper:
%
%      C. Li, C. Xu, C. Gui, M. D. Fox, "Distance Regularized Level Set Evolution and Its Application to Image Segmentation",
%        IEEE Trans. Image Processing, vol. 19 (12), pp.3243-3254, 2010.

close all

%% step1, read grayscale image
Img=imread('Inputs/out_meanshift_big.jpg');
Img=rgb2gray(Img);
frm=0;
%% step2, set params
timestep=30;  % time step
mu=0.2;  % coefficient of the distance regularization term R(phi)
lambda=5; %coefficient of the weighted length term L(phi)
alfa= -3;  %  coefficient of the weighted area term A(phi)
epsilon=1.5; % papramater that specifies the width of the DiracDelta function
c0=2;
maxiter=220;
sigma=2.0;    % scale parameter in Gaussian kernel

%% step3 smooth image with gaussian filter
G=fspecial('gaussian',30,sigma); % 15 Caussian kernel
Img_smooth=conv2(Img,G,'same');  % smooth image by Gaussiin convolution
figure(1);
imagesc(Img_smooth,[0, 255]); axis off; axis equal; colormap(gray);
title('Smoothed image');

%% step4 calculate edge indicator according to Eq23
[Ix,Iy]=gradient(Img_smooth);
f=Ix.^2+Iy.^2;
%load hednms
%f=double(E);
g=1./(1+f);  % edge indicator function.
g=exp(-f);
figure(2);
imagesc(g); axis off; axis equal;
title('g,  edge indicator');

%% step5, set initial phi

im=imread('Inputs/out_meanshift_big_draw4.jpg');
bw=im(:,:,1)>220;
L=bwlabel(bw);
ul=unique(L(:));
a=hist(L(:),ul);
a=sortrows([a' ul],-1);
bw=L==a(2,2);

bww=imfill(bw, 'holes')-bw;

L=bwlabel(bww);
ul=unique(L(:));
a=hist(L(:),ul);
a=sortrows([a' ul],-1);
bw=L==a(2,2);

%bw=zeros(size(bw));
%bw(510:1060,280:950)=1;

phi = -c0*ones(size(Img));
phi(bw==1)=c0;
%phi(350:370,200:220)=-c0;
figure(3);
imagesc(phi);
axis off; axis equal;colormap(jet);
title('initial phi matrix');

[vx, vy]=gradient(g);
figure(4);
subplot(1,2,1);imagesc(vx); title('x directioned gradient of g');
subplot(1,2,2);imagesc(vy); title('y directioned gradient of g');

for k=1:maxiter

    if mod(k,1)==0
        frm=frm+1;
        %close all;
        h=figure(5);
        set(gcf,'color','w');
        %subplot(1,2,1);
        II=Img;
        II(:,:,2)=Img;II(:,:,3)=Img;
        imshow(II); axis off; axis equal; hold on;
        q=contour(phi, [0,0], 'r');
        msg=['contour result , iteration number=' num2str(k)];
        title(msg);

        frame = getframe(h);
        im = frame2im(frame);
        [imind,cm] = rgb2ind(im,256);
        %Write to the GIF File
        if frm == 1
            imwrite(imind,cm,'Outputs/cheval2.gif','gif', 'Loopcount',inf);
        else
            imwrite(imind,cm,'Outputs/cheval2.gif','gif','WriteMode','append');
        end
    end

    %% step6, check boundary conditions
    phi=NeumannBoundCond(phi);

    %% step 7 calculate differential of regularized term in Eq.30
    distRegTerm=distReg_p2(phi);

    %% step8 calculate differential of area term in Eq.30
    diracPhi=Dirac(phi,epsilon);
    areaTerm=diracPhi.*g;

    %% step9 calculate differential of length term in Eq.30
    [phi_x,phi_y]=gradient(phi);
    s=sqrt(phi_x.^2 + phi_y.^2);
    Nx=phi_x./(s+1e-10); % add a small positive number to avoid division by zero
    Ny=phi_y./(s+1e-10);
    edgeTerm=diracPhi.*(vx.*Nx+vy.*Ny) + diracPhi.*g.*div(Nx,Ny);

    %% step 10 update phi according to Eq.20
    phi=phi + timestep*(mu/timestep*distRegTerm + lambda*edgeTerm + alfa*areaTerm);

    %% show result in every 50 iteration

    %% step 11 if maxiter done then finish, else return step6
end
%% Step 12. show last iteration results
figure(6);
imagesc(Img,[0, 255]); axis off; axis equal; colormap(gray); hold on;  contour(phi, [0,0], 'r');
msg=['phi result , iteration number=' num2str(k)];
title(msg);

function f = distReg_p2(phi)
% compute the distance regularization term with the double-well potential p2 in eqaution (16)
[phi_x,phi_y]=gradient(phi);
s=sqrt(phi_x.^2 + phi_y.^2);
a=(s>=0) & (s<=1);
b=(s>1);
ps=a.*sin(2*pi*s)/(2*pi)+b.*(s-1);  % compute first order derivative of the double-well potential p2 in eqaution (16)
dps=((ps~=0).*ps+(ps==0))./((s~=0).*s+(s==0));  % compute d_p(s)=p'(s)/s in equation (10). As s-->0, we have d_p(s)-->1 according to equation (18)
f = div(dps.*phi_x - phi_x, dps.*phi_y - phi_y) + 4*del2(phi);

function f = div(nx,ny)
[nxx,junk]=gradient(nx);
[junk,nyy]=gradient(ny);
f=nxx+nyy;

function f = Dirac(x, sigma)
f=(1/2/sigma)*(1+cos(pi*x/sigma));
b = (x<=sigma) & (x>=-sigma);
f = f.*b;

function g = NeumannBoundCond(f)
% Make a function satisfy Neumann boundary condition
[nrow,ncol] = size(f);
g = f;
g([1 nrow],[1 ncol]) = g([3 nrow-2],[3 ncol-2]);
g([1 nrow],2:end-1) = g([3 nrow-2],2:end-1);
g(2:end-1,[1 ncol]) = g(2:end-1,[3 ncol-2]);

author: balcilar

repo: DRLSE-Image-Segmentation


text:
%GRAD_AND_DIV
% Generates pressure gradient matrices Pu, Pv and divergence matrix D
% Vectorized matrix generation by evaluation of stencil coefficients
% Output: Pu, Pv, Du, Dv

J= length(dx); K = length(dy);

%.................Pressure gradient matrix Pu ....................

%        	  |    0    |
% Stencil: [Pu] = |p1 p2  0 |
%        	  |    0    |

% Indexing convention in staggered grid:
%  +--k+1--+
%  |       |
%  j  jk  j+1
%  |       |
%  +---k---+
tic
p2 = 1./DXU(1,:)'; p1 = zeros(size(p2)); p1(1:J) = -p2(2:J+1);
Puu = spdiags([p1  p2],[-1;0], J+1, J); Pu = kron(speye(K),Puu);

%......................Boundary corrections....................................
for seg = 1:length(ybc(1,:))
  if (ybc(1,seg) == 1)|(ybc(1,seg) == 2)  % No-slip or inflow at left boundary
    k = 1+kseg(seg):kseg(seg+1); Pu(1+(k-1)*(J+1),:) = 0;
  end
  if (ybc(2,seg) == 1)|(ybc(2,seg) == 2)  % No-slip or inflow at right boundary
    k = 1+kseg(seg):kseg(seg+1);  Pu(k*(J+1),:) = 0;
  end
end
tijd = toc; disp(['Breakdown of grad_and_div time'])
disp(['  Pu time = ',num2str(tijd)])

%.................Pressure gradient matrix Pv ....................

%        	  |     0   |
% Stencil: [Pv] = |0   p2  0|
%        	  |    p1   |

tic
pp1 = zeros(size(XV)); pp2 = pp1;    		% Diagonals of Pv
pp1(2:K+1,:) = -1./DYV(2:K+1,:); 	pp2 = -pp1;
pp2(1,:) = 1./DYV(1,:); 		pp2(K+1,:) = 0;

%......................Boundary corrections....................................
for seg = 1:length(xbc(1,:))
  if (xbc(1,seg) == 1)|(xbc(1,seg) == 2)  % No-slip or inflow at lower boundary
    j = 1+jseg(seg):jseg(seg+1);     	pp1(1,j) = 0; pp2(1,j) = 0;
  end
  if (xbc(2,seg) == 1)|(xbc(2,seg) == 2)  % No-slip or inflow at upper boundary
    j = 1+jseg(seg):jseg(seg+1);	pp1(K+1,j) = 0; pp2(K+1,j) = 0;
  end
end
n = J*(K+1);   p1 = reshape(pp1',n,1);    p2 = reshape(pp2',n,1);
p1 = [p1(J+1:n); zeros(J,1)];		% Shift to accomodate spdiags
Pv = spdiags([p1  p2],[-J;0], n,n-J);
tijd = toc; disp(['  Pv time = ',num2str(tijd)])

%.................u divergence matrix Du ....................

%		  |   0     |
% Stencil: [Du] = |0  p1  p2|
%		  |   0     |

tic
Duu = spdiags([-ones(J,1) ones(J,1)], [0;1], J,J+1);
Du = kron(spdiags(dy',0,K,K),Duu);
tijd = toc; disp(['  Du time = ',num2str(tijd)])

%.................v divergence matrix Dv ....................

%		  |   p2   |
% Stencil: [Dv] = |0  p1  0|
%		  |   0    |

tic
Dvv = spdiags([-ones(K,1) ones(K,1)], [0;1], K,K+1);
Dv = kron(Dvv,spdiags(dx',0,J,J));
tijd = toc; disp(['  Dv time = ',num2str(tijd)])

clear  pp1 pp2 p1 p2 Puu Duu Dvv 			% Save storage

repo: Aero-matlab


text:
function [A, b, x, ProbInfo] = PRdiffusion(varargin)
% PRdiffusion Generates data for use in an inverse diffusion problems
%
% [A, b, x, ProbInfo] = PRdiffusion
% [A, b, x, ProbInfo] = PRdiffusion(n)
% [A, b, x, ProbInfo] = PRdiffusion(n, options)
%
% This function generates a 2D inverse diffusion problem on a uniform
% finite-element mesh with 2*(n-1)^2 triangular elements; think of the
% domain as an n-by-n pixel grid with two triangular elements in each pixel.
%
% The data represents a function u(t) on the grid at time t = Tfinal.
% The goal of the inverse problem is to compute the function u(t) at
% time t = 0 which, when diffused, produces the data.
%
% The forward problem is a very simple diffusion problem
%    du/dt = u_xx + u_yy,   0 < t < Tfinal
% with homogenous Neumann boundary conditions.  It is solved by means of
% the Crank-Nicolson-Galerkin finite-element method.  The inverse problem
% is to compute the initial condition, the function u(0) at time t = 0,
% from the solution u(t) at time t = Tfinal.
%
% The functions u(0) and u(Tfinal) are represented by the vecors x and b,
% resp., which hold the values at the corners of the finite elements.
% Both vectors have a total of n^2 elements.
%
% The severity of the problem increases with Tfinal and does not depend on
% Tsteps.
%
% Input:
%  n       - The grid has n-1 pixels in each direction, and n must be an
%            integer.  Default: n = 128.  There are 2*(n-1)^2 finite
%            elements, and the number of unknowns is n^2.
%  options - Parameters that define the problem
%              Tfinal : diffusion time; default Tfinal = 0.01.
%              Tsteps : number of time steps; default Tsteps = 100.
%
% Output:
%  A -  Function handle for the forward and adjoint problem.
%  b -  Vector that represents the solution u(t) at t = Tfinal.
%  x -  Vector that represents the initial condition u(t) at t = 0.
%  ProbInfo -  Structure containing some information about problem
%                problemType : type of test problem generated
%                              (in this case: 'diffusion')
%                xType       : solution type (in this case 'fem')
%                bType       : data type (in this case 'fem')
%                elmtab      : table of FEM elements
%                elmX        : X-coordinates of element grid
%                elmY        : Y-coordinates of element grid
%
% See also: PRblur, PRinvinterp2, PRnmr, PRseismic, PRspherical, PRtomo,
% PRnoise, PRshowb, PRshowx

% Silvia Gazzola, University of Bath
% Per Christian Hansen, Technical University of Denmark
% James G. Nagy, Emory University
% April, 2018.

% Based on original code by Allan P. Engsig-Karup, DTU Compute.

% This file is part of the IR Tools package and is distributed under the
% 3-Clause BSD Licence. A separate license file should be provided as part
% of the package.

% Set default values for options.
defaultopt = struct('Tfinal', 0.01, 'Tsteps', 100);

% If input is 'defaults,' return the default options in A.
if nargin == 1 && nargout <= 1 && strcmp(varargin,'defaults')
    A = defaultopt;
    return;
end

% Check for acceptable number of optional input arguments.
switch length(varargin)
    case 0
        n = []; options = [];
    case 1
        if isa(varargin{1}, 'double')
            n = varargin{1}; options = [];
        else
            n = []; options = varargin{1};
        end
    case 2
        if isa(varargin{1}, 'double')
            n = varargin{1}; options = varargin{2};
        else
            n = varargin{2}; options = varargin{1};
        end
    otherwise
        error('Too many input parameters')
end

if isempty(options)
    options = defaultopt;
end

options = PRset(defaultopt, options);
Tfinal = PRget(options, 'Tfinal', [], 'fast');
Tsteps = PRget(options, 'Tsteps', [], 'fast');

if isempty(n), n = 128; end

% Create rectangular mesh consisting of uniformly distributed triangles.
[elmX,elmY] = meshgrid(linspace(0,1,n));
elmX = elmX(:);
elmY = elmY(:);

% Element table for uniformly distributed triangles in a rectangular domain.
elmtab = zeros(2*(n-1)^2,3);
count = 0;
for j = 1:n-1
    for i = 1:n-1
        elmtab(count+(1:2),:) = [i   i+1+n i+n;
                                 i+1 i+1+n i     ] + (j-1)*n;
        count = count+2;
    end
end

% Initial condition.
x = 0.7*exp( -((elmX-0.4)/0.12).^2 - ((elmY-0.5)/0.15).^2) + ...
        exp( -((elmX-0.7)/0.1 ).^2 - ((elmY-0.4)/0.08).^2);

% Create the matrices.
dt = Tfinal/Tsteps;
[AA,CC] = assembly(elmX,elmY,elmtab);
S = CC - 0.5*dt*AA;
R = CC + 0.5*dt*AA;

% A is a function handle to the function that carries out the forward
% operation, i.e., interpolates from the regular grid to the randon points,
% or the corresponding adjoint operation (representing A').
A = @(xx,tflag) OPdiffusion(xx,R,S,Tsteps,tflag);

% The right-hand side (the data) ...
b = A(x,'notransp');

ProbInfo.problemType = 'diffusion';
ProbInfo.bType = 'fem';
ProbInfo.xType = 'fem';
ProbInfo.elmtab = elmtab;
ProbInfo.elmX = elmX;
ProbInfo.elmY = elmY;
ProbInfo.R = R;                  %PCH
ProbInfo.S = S;

% SUBFUNCTIONS -----------------------------------------------------------

function [A,C] = assembly(x,y,elmtab)
% Generate linear system coefficient matrix and right hand side vector.
% Symmetry is not exploited.

P = size(elmtab,1);
iL = zeros(3*3*P,1);
jL = iL;
AL = iL;
CL = iL;
count = 0;
for p = 1:P
    [delta,abc] = basfun(p,x,y,elmtab);
    for r = 1:3
        i = elmtab(p,r);
        for s = 1:3   % Does not exploit symmetry.
            j = elmtab(p,s);
            Ars = 1/(4*abs(delta))*(abc(r,2)*abc(s,2) + abc(r,3)*abc(s,3));
            if r == s
                Crs = abs(delta)/6;
            else
                Crs = abs(delta)/12;
            end
            count = count + 1;
            iL(count) = i;
            jL(count) = j;
            AL(count) = Ars;
            CL(count) = Crs;
        end
    end
end
iL = iL(1:count);
jL = jL(1:count);
AL = AL(1:count);
A = sparse(iL,jL,AL);
C = sparse(iL,jL,CL);

function [delta,abc] = basfun(p,x,y,elmtab)
% For a given element compute the geometric properties of the element.

xc = x(elmtab(p,:));
yc = y(elmtab(p,:));
delta = 0.5*( xc(2)*yc(3) - yc(2)*xc(3) - (xc(1)*yc(3) - ...
              yc(1)*xc(3)) + xc(1)*yc(2) - yc(1)*xc(2)      );

% Vectorized output.
j   = [2 3 1]';
k   = [3 1 2]';
at  = xc(j).*yc(k) - xc(k).*yc(j);
bt  = yc(j) - yc(k);
ct  = xc(k) - xc(j);
abc = [at bt ct];

repo: IRtools


text:
function datapath = metoffice_get_path()

os = system_dependent('getos');
if ~isempty(strfind(os,'Ubuntu'))
  % Finland
  datapath = '/share/climate/data/UK_Met_Office/RecTest/DATASETS/';
elseif ~isempty(strfind(os,'Linux'))
  % ICS GridEngine
  datapath = '/share/climate/data/UK_Met_Office/RecTest/DATASETS/';
else
  % MetOffice
  datapath = '/net/home/h03/hadia/data/RecTest/';
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
clear all;

%%%%%%%%%%%%%This example will need to be reorganized

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%provide parameters and inputs below
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

addpath('../')

addpath(genpath('../CoreModules'));

addpath('./lm_data');

n_epoch=20;%20 %%training epochs
dataset_name='char'; % dataset name
network_name='qrnn';%'gru';'rnn','lstm'
use_gpu=0; %%use gpu or not
opts.use_nntoolbox=0;

PrepareDataFunc=@PrepareData_Char_RNN; %%function handler to prepare your data

if strcmp(network_name,'lstm')
    NetInit=@net_init_char_lstm;%_bn  %% function to initialize the network
end

if strcmp(network_name,'gru')
    NetInit=@net_init_char_gru;  %% function to initialize the network
end

if strcmp(network_name,'rnn')
    NetInit=@net_init_char_rnn;  %% function to initialize the network
    opts.parameters.Id_w=1;%vanilla rnn:0, rnn with skip links: 1
end

if strcmp(network_name,'qrnn')
    NetInit=@net_init_char_qrnn;  %% function to initialize the network
end

use_selective_sgd=0; %automatically select learning rates
%%selective-sgd parameters
%ssgd_search_freq=10; %select new coarse-scale learning rates every n epochs

learning_method=@adam; %training method: @rmsprop;

%sgd parameter (unnecessary if selective-sgd is used)
sgd_lr=1e-2;

opts.parameters.batch_size=100;
opts.parameters.n_hidden_nodes=30;
opts.parameters.n_cell_nodes=30;
opts.parameters.n_input_nodes=67;
opts.parameters.n_output_nodes=67;
if strcmp(network_name,'lstm')
    opts.parameters.n_gates=3;
end
if strcmp(network_name,'gru')
    opts.parameters.n_gates=2;
end

if strcmp(network_name,'qrnn')
    opts.parameters.n_gates=1;
    use_gpu=(gpuDeviceCount>0);
    opts.use_nntoolbox=license('test','neural_network_toolbox');
end

opts.parameters.n_frames=64;%%%%max sentence length

opts.parameters.lr =sgd_lr;
opts.parameters.mom =0.9;
opts.parameters.learning_method=learning_method;
opts.parameters.selective_sgd=use_selective_sgd;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%provide parameters and inputs above
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%stupid settings below (so please ignore)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

opts.n_epoch=n_epoch; %training epochs
opts.dataset_name=dataset_name; %dataset name
opts.network_name=network_name; %network name
opts.use_gpu=use_gpu; %use gpu or not

opts.results=[];
opts.results.TrainEpochError=[];
opts.results.TestEpochError=[];
opts.results.TrainEpochLoss=[];
opts.results.TestEpochLoss=[];
opts.RecordStats=1;
opts.results.TrainLoss=[];
opts.results.TrainError=[];

opts.plot=1;

opts.dataDir=['./',opts.dataset_name,'/'];
opts=PrepareDataFunc(opts);

net=NetInit(opts);

opts=generate_output_filename(opts);

if(opts.use_gpu)
    for i=1:length(net)
        net{i}=SwitchProcessor(net{i},'gpu');
    end
else
    for i=1:length(net)
        net{i}=SwitchProcessor(net{i},'cpu');
    end
end

opts.n_batch=floor(opts.n_train/opts.parameters.batch_size);
opts.n_test_batch=floor(opts.n_test/opts.parameters.batch_size);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%stupid settings above
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%training goes below
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

opts.parameters.current_ep=1;

start_ep=opts.parameters.current_ep;
if opts.plot
    figure1=figure;
end
for ep=start_ep:opts.n_epoch

    [net,opts]=train_rnn(net,opts);
    [opts]=test_rnn(net,opts);
    opts.parameters.current_ep=opts.parameters.current_ep+1;
    disp(['Epoch ',num2str(ep),' testing error: ',num2str(opts.results.TestEpochError(end)), ' testing loss: ',num2str(opts.results.TestEpochLoss(end))])

    %
    if opts.plot
        subplot(1,2,1); plot(opts.results.TrainEpochError);hold on;plot(opts.results.TestEpochError);hold off;title('Error Rate per Epoch')
        subplot(1,2,2);plot(opts.results.TrainEpochLoss);hold on;plot(opts.results.TestEpochLoss);hold off;title('Loss per Epoch')
        drawnow;
    end
    %}
    parameters=opts.parameters;
    results=opts.results;
    save([fullfile(opts.output_dir,[opts.output_name,num2str(ep),'.mat'])],'net','parameters','results');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%training goes above
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
% RANDSAMPLE Random sample, with or without replacement.
%    Y = RANDSAMPLE(N,K) returns Y as a column vector of K values sampled
%    uniformly at random, without replacement, from the integers 1:N.
%
%    Y = RANDSAMPLE(POPULATION,K) returns K values sampled uniformly at random,
%    without replacement, from the values in the vector POPULATION.  Y is a
%    vector of the same type as POPULATION.  NOTE:  When POPULATION is a
%    numeric vector containing only non-negative integer values, and it might
%    have length 1, use
%
%       Y = POPULATION(RANDSAMPLE(LENGTH(POPULATION),K)
%
%    instead of Y = RANDSAMPLE(POPULATION,K).
%
%    Y = RANDSAMPLE(N,K,REPLACE) or RANDSAMPLE(POPULATION,K,REPLACE) returns a
%    sample taken with replacement if REPLACE is true, or without replacement
%    if REPLACE is false (the default).
%
%    Y = RANDSAMPLE(N,K,true,W) or RANDSAMPLE(POPULATION,K,true,W) returns a
%    weighted sample, using positive weights W, taken with replacement.  W is
%    often a vector of probabilities. This function does not support weighted
%    sampling without replacement.
%
%    Y = RANDSAMPLE(S,...) uses the random number stream S for random number
%    generation.  RANDSAMPLE uses the MATLAB default random number stream by
%    default.
%
%    Examples:
%
%    Draw a single value from the integers 1:10.
%       n = 10;
%       x = randsample(n,1);
%
%    Draw a single value from the population 1:n, when n > 1.
%       y = randsample(1:n,1);
%
%    Generate a random sequence of the characters ACGT, with
%    replacement, according to specified probabilities.
%       R = randsample('ACGT',48,true,[0.15 0.35 0.35 0.15])
%
%    See also RAND, RANDPERM, RANDSTREAM.
%
%    Reference page in Doc Center
%       doc randsample
%
%    Other functions named randsample
%
%       generic/randsample
%


text:
function r=rsi2(x,N)
L = length(x);
dx = diff([0;x]);
up=dx;
down=abs(dx);
% up and down moves
I=dx<=0;
up(I) = 0;
down(~I)=0;
% calculate exponential moving averages
m1 = ema(up,N); m2 = ema(down,N);
warning off
r = 100*m1./(m1+m2);
%r(isnan(r))=50;
I2=~((up+down)>0);
r(I2)=50;
warning on

author: Sable

repo: mcbench-benchmarks


text:
function qUnNormalised = unnormaliseImagePoints(qNormalised,K,kc)
%function to unnormalise a set of image points according to the perspective
%intrinsic calibration matrix K and distortion parameters kc. This is the model
%from Oulu university and described clearly in Bouguet's Matlab camera calibration toolbox
% (http://www.vision.caltech.edu/bouguetj/calib_doc/). This is done using
% code from Bouguet's toolbox.
%
%inputs:
%q: a 2*N matrix of normalised points in an image.
%
%K: a 3x3 perspective intrinsic calibration matrix
%
%kc: a 5x1 vector of radial and tangential distortion coefficients.
%
%outputs:
%qUnNormalised: the un-normalised points

%first compensate for radial and tangential distortion:
xn = qNormalised;
x = xn(1,:);
y = xn(2,:);

r = sqrt(xn(1,:).^2 + xn(2,:).^2);

dx(1,:) = 2*kc(3)*x.*y + kc(4)*(r.^2 + 2*x.^2);
dx(2,:) = kc(3)*(r.^2 + 2*y.^2) + 2*kc(4)*x.*y;

s = (1+kc(1)*r.^2 + kc(2)*r.^4 + kc(5)*r.^6);
xd = [s;s].*xn + dx;
xd(3,:) = 1;

%now apply the intrinsic matrix
qUnNormalised = K*xd;
qUnNormalised(1,:) = qUnNormalised(1,:)./qUnNormalised(3,:);
qUnNormalised(2,:) = qUnNormalised(2,:)./qUnNormalised(3,:);
qUnNormalised = qUnNormalised(1:2,:);

author: tobycollins

repo: IPPE

save_path: github-repos/MATLAB/tobycollins-IPPE


text:
function [yMeas] = stereoCamProject(p_fc_c, calibParams)
% triangulate Triangulates 3D point from stereo camera measurement

c_u = calibParams.c_u;
c_v = calibParams.c_v;
f_u = calibParams.f_u;
f_v = calibParams.f_v;
b = calibParams.b;

x = p_fc_c(1);
y = p_fc_c(2);
z = p_fc_c(3);

yMeas = (1/z)*[f_u*x; f_v*y; f_u*(x-b); f_v*y] + [c_u;c_v;c_u;c_v];
end

repo: MSCKF