text:
function [xsol,objp,info,msg,times,x,y,z,s,w] = lipsol(A,b,c,lb,ub,opts)
% LIPSOL        - Main program for LIPSOL

% Yin Zhang, January, 1995
% Department of Mathematics and Statistics
% University of Maryland  Baltimore County
% Modified J.Currie AUT May 2013

global probData

if(nargin < 6)
    opts = lipsolset;
    opts.maxiter = 100;
    opts.maxtime = 1000;
    opts.tol = 1e-5;
end
if(nargin < 5), ub = []; end
if(nargin < 4), lb = []; end
if(nargin < 3), error('You must supply at least 3 arguments to lipsol(A,b,c)'); end

times = zeros(5,1);
t0 = tic;

%Problem Data Structure to avoid so many Global Variables
probData = struct('NNZA',0,'NNZL',0,'XLNZ',[],'PERM',[],'INVP',[],'LNZ',[],'LOOP_LEVEL',[],...
                  'XSUPER',[],'XLINDX',[],'LINDX',[],'SNODE',[],'SPLIT',[],'TMPSIZ',[],...
                  'LNZ0',[],'b_orig',b,'c_orig',c,'data_changed',0,'isFeasible',1,...
                  'Lbounds_non0',[],'Ubounds_exist',0,'nub',0,'Fixed_exist',0,...
                  'ifix',[],'infx',[],'xfix',[],'Zrcols_exist',0,'izrcol',[],...
                  'inzcol',[],'xzrcol',[],'Sgtons_exist',0,'isolved',[],...
                  'insolved',[],'xsolved',[],'nt',0,'message',[],...
                  'col_scaled',0,'colscl',[],'Dense_cols_exist',0,'Mdense',[],...
                  'idense',[],'ispars',[],'Hist',[],'Sherman_OK',[],'backed',[]);
%Initial Time
opts.t0 = tic;

%Save loop_level into problem structure
probData.LOOP_LEVEL = opts.loop_level;

%Pretty Printing
if(opts.verb), fprintf('----------------------------------------------\n'); end

%Preprocess the Input Data
[A,b,c,lb,ub] = ls_preprocess(A,b,c,lb,ub,opts);
times(1) = toc(t0);

%Check if preproccesor indicates infeasible
if(~probData.isFeasible)
    xsol = []; objp = []; info = [-1 0 0]; msg = 'Presolve Indicated Infeasible'; times = [];
    x = []; y = []; z = []; s = []; w = [];
    return;
end

%Scale the Input Data if Required
[A,b,c,ub] = scaling(A,b,c,ub,opts);
%Determine density of the problem
checkDense(A,opts);
times(2) = toc(t0) - times(1);

%Solve the Problem
[x,y,z,s,w,info] = ls_miip(A,b,c,ub,opts);
times(3) = toc(t0) - times(2);

%Display Solve Status
if(opts.verb && ~isempty(probData.message)), fprintf('\n  Status: %s\n',probData.message); end

%Post Process the Solution
[xsol, objp, times ] = ls_postprocess(x,y,lb,opts.verb,times,t0);
msg = probData.message;

%Pretty Printing
if(opts.verb), fprintf('----------------------------------------------\n'); end

function [A,b,c,ub] = scaling(A,b,c,ub,opts)
% SCALING     - Scale matrix A. Usage: [A,b,c,ubounds] = scaling(A,b,c,ubounds)
% Yin Zhang, January, 1995
% Department of Mathematics and Statistics
% University of Maryland  Baltimore County

global probData

[m, n] = size(A);
badscl = 1.e-4;

absnzs = abs(nonzeros(A));
thescl = min(absnzs)/max(absnzs);
clear absnzs

if (thescl < badscl)
    if(opts.verb), fprintf('Scaling ...\n'); end

    % ----- scaling vectors ------
    AA = abs(A);
    colscl = full(sqrt(max(AA)'));
    rowscl = full(sqrt(max(AA,[],2)));
    clear AA;

    % ----- column scaling -----
    if (probData.Ubounds_exist), ub = ub.*colscl; end
    colscl = reciprocal(colscl);
    A = A*sparse(1:n,1:n,colscl);
    c = c.*colscl;
    probData.col_scaled = 1;

    % ----- row scaling -----
    rowscl = reciprocal(rowscl);
    A = sparse(1:m,1:m,rowscl)*A;
    b = b.*rowscl;
    bnrm = norm(b);
    if (bnrm > 0)
      q = median([1 norm(c)/bnrm 1.e+8]);
      if (q > 10), A = q*A; b = q*b; end
    end
    probData.data_changed = 1;
    probData.colscl = colscl;
end

function checkDense(A,opts) %#ok<INUSD>
% CHECKDENSE - Determine and locate dense columns of matrix A.
% Yin Zhang, January, 1995
% Department of Mathematics and Statistics
% University of Maryland  Baltimore County

%PROBLEM IN PCG - SKIPPING DENSE COLUMN MANAGEMENT [J.C. 19/5/13]

% global probData
%
% [m, n] = size(A);
% probData.ispars = 1:n;
% nzratio = 1;
% if (m >  500), nzratio = 0.20; end;
% if (m > 1000), nzratio = 0.10; end;
% if (m > 5000), nzratio = 0.05; end;
%
% if (nzratio < 1)
%    checking = sum(spones(A))/m <= nzratio;
%    if any(checking == 0)
%       probData.Dense_cols_exist = 1;
%       probData.idense = find(sparse(1-checking));   % Dense  column indices
%       probData.ispars = find(checking);             % Sparse column indices
%       if(opts.verb), fprintf('Dense columns (nnz/m > %g): %i\n',nzratio,length(probData.idense)); end
%    else
%       if(opts.verb), fprintf('No Significant Dense Columns\n'); end
%    end
%    clear checking
% else
%     if(opts.verb), fprintf('Problem is Small Enough to Ignore Dense Columns\n'); end
% end

function Y = reciprocal(X)
% RECIPROCAL  - Invert the nonzero entries of a matrix elementwise.
%             Y = RECIPROCAL(X) has the same sparsity pattern as X
%	      (except possibly for underflow).

% Yin Zhang, January, 1995
% Department of Mathematics and Statistics
% University of Maryland  Baltimore County

if issparse(X)
   [m, n]  = size(X);
   [i,j,Y] = find(X);
   Y = sparse(i,j,1./Y,m,n);
else
   Y = 1./X;
end
ceiling = 1.e+16; Y = min(ceiling,Y);


text:
function mask = PTKGetLobeROI(lobe_mask, context, reporting)
    % PTKGetLobeROI. Extracts a region of interest for one of the lobes.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    mask = lobe_mask.BlankCopy;

    if context == PTKContext.LeftUpperLobe
        colormap_index = PTKColormapLabels.LeftUpperLobe;
    elseif context == PTKContext.LeftLowerLobe
        colormap_index = PTKColormapLabels.LeftLowerLobe;
    elseif  context == PTKContext.RightUpperLobe
        colormap_index = PTKColormapLabels.RightUpperLobe;
    elseif  context == PTKContext.RightMiddleLobe
        colormap_index = PTKColormapLabels.RightMiddleLobe;
    elseif  context == PTKContext.RightLowerLobe
        colormap_index = PTKColormapLabels.RightLowerLobe;
    else
        reporting.Error('PTKGetLobeROI:UnknownContext', ['The context ' char(context) ' is unknown.']);
    end
    mask.ChangeRawImage(lobe_mask.RawImage == colormap_index);
    mask.CropToFitWithBorder(5);
end

author: tomdoel

repo: pulmonarytoolkit


text:
function B = bound2im(b, M, N, x0, y0)
%BOUND2IM Converts a boundary to an image.
%   B = BOUND2IM(b) converts b, an np-by-2 or 2-by-np array
%   representing the integer coordinates of a boundary, into a binary
%   image with 1s in the locations defined by the coordinates in b
%   and 0s elsewhere.
%
%   B = BOUND2IM(b, M, N) places the boundary approximately centered
%   in an M-by-N image. If any part of the boundary is outside the
%   M-by-N rectangle, an error is issued.
%
%   B = BOUND2IM(b, M, N, X0, Y0) places the boundary in an image of
%   size M-by-N, with the topmost boundary point located at X0 and
%   the leftmost point located at Y0. If the shifted boundary is
%   outside the M-by-N rectangle, an error is issued. XO and X0 must
%   be positive integers.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.6 $  $Date: 2003/06/14 16:21:28 $

[np, nc] = size(b);
if np < nc
   b = b'; % To convert to size np-by-2.
   [np, nc] = size(b);
end

% Make sure the coordinates are integers.
x = round(b(:, 1));
y = round(b(:, 2));

% Set up the default size parameters.
x = x - min(x) + 1;
y = y - min(y) + 1;
B = false(max(x), max(y));
C = max(x) - min(x) + 1;
D = max(y) - min(y) + 1;

if nargin == 1
   % Use the preceding default values.
elseif nargin == 3
   if C > M | D > N
      error('The boundary is outside the M-by-N region.')
   end
   % The image size will be M-by-N. Set up the parameters for this.
   B = false(M, N);
   % Distribute extra rows approx. even between top and bottom.
   NR = round((M - C)/2);
   NC = round((N - D)/2); % The same for columns.
   x = x + NR; % Offset the boundary to new position.
   y = y + NC;
elseif nargin == 5
   if x0 < 0 | y0 < 0
      error('x0 and y0 must be positive integers.')
   end
   x = x + round(x0) - 1;
   y = y + round(y0) - 1;
   C = C + x0 - 1;
   D = D + y0 - 1;
   if C > M | D > N
      error('The shifted boundary is outside the M-by-N region.')
   end
   B = false(M, N);
else
   error('Incorrect number of inputs.')
end

B(sub2ind(size(B), x, y)) = true;

repo: weiyanmin


text:
function eev(obj, eventnum)
    % CATALOG.EEV - Browse an Catalog object one event at a time.
    %  catalogObject.EEV() Browse through an Catalog object one event
    %  at a time in a similar way to the Seisan program 'eev'.

    if ~exist('eventnum','var')
        eventnum = 1;
    end

    while 1,

        % don't beyond start or end of this catalogObject object
        if eventnum<1
            eventnum=1;
        end
        if eventnum>numel(obj.otime)
            eventnum=numel(obj.otime);
        end
        % display line for this event
        dstr=datestr(obj.otime(eventnum),31);
        subclass=obj.etype{eventnum};
        mag=obj.mag(eventnum);
        outstr=sprintf('%s %7.2f %7.2f %7.2 %5.1f %s %s',dstr, obj.lon(eventnum), obj.lat(eventnum), obj.depth(eventnum), mag, obj.magtype{eventnum}, subclass);
        choice=input([outstr,':  ?'],'s');

        % process choice
        if isempty(choice)
            eventnum=eventnum+1; % ENTER goes to next event

        elseif (choice(1)=='c') % CLASSIFY
            classify_event(datenum);
% ^^ does this line need obj.otime?

        elseif (choice(1)=='f') % FORWARD N EVENTS
            num=1;
            if length(choice)>1
                num=str2num(choice(2:end));
            end
            eventnum=eventnum+num;

        elseif (choice(1)=='b') % BACKWARD N EVENTS
            num=1;
            if length(choice)>1
                num=str2num(choice(2:end));
            end
            eventnum=eventnum-num;

        elseif (choice(1)=='t') % JUMP TO TIME
            month=1;dd=1;hr=0;
            if length(choice)>4
                year=str2num(choice(2:5));
            end
            if length(choice)>6
                month=str2num(choice(6:7));
            end
            if length(choice)>8
                dd=str2num(choice(8:9));
            end
            if length(choice)>10
                hr=str2num(choice(10:11));
            end
            jumptime=datenum(year,month,dd,hr,0,0);
            eventnum = min(find(obj.otime >= jumptime));

        elseif (choice(1)=='s') % SUMMARISE - SHOW S FILE or similar data
            fprintf('\nTime:\t\t%s\n',dstr);
            fprintf('Longitude:\t%7.2f degrees\n',obj.lon(eventnum));
            fprintf('Latitude:\t%7.2f degrees\n',obj.lat(eventnum));
            fprintf('Depth:\t\t%7.2f km\n',obj.depth(eventnum));
            fprintf('Magnitude:\t%7.2f\n',obj.mag(eventnum));
            fprintf('Magnitude Type:\t%s\n',obj.magtype{eventnum});
            fprintf('Event Type:\t%s\n',obj.etype{eventnum});
            fprintf('\n');

        elseif (choice(1)=='x') % CLOSE ALL
            close all;

        elseif (choice=='q') % QUIT
            break;

        elseif (choice(1)=='h') % HELP
            disp(' ');
            disp('Options:');
            disp('________');
            disp(' ');
            disp('b[num]            - go backward 1 event (or N events)');
            disp('c                 - classify');
            %disp('e                 - edit/generate S-file');
            disp('f[num]            - go forward 1 event (or N events)');
            disp('h                 - this help');
            %disp('p                 - plot');
            disp('s                 - summarise');
            disp('tYYYY[MM[DD[HH]]] - jump to date/hour specified');
            disp('x                 - close all figure windows');
            disp('q                 - quit');
            disp(' ');
        end
    end
end

author: geoscience-community-codes

repo: GISMO


text:
classdef shape2d < grain2d

  properties (Dependent=true)
    Vs      % list of vertices
    rho     % radius of polar coords
    theta   % angle of polar coords
  end

  % 1) should be constructed from calcTDF / circdensity (density function from a set of lines)
  % characteristicshape
  % surfor, paror
  % 2) purpose: take advantage of grain functions (long axis direction, aspect ratio..)
  %
  % additional functions I will try to put here: measure of asymmetry
  % nice plotting wrapper (replacing plotTDF)

  methods

    function shape = shape2d(V,CS)
      % list of vertices [x y]

      if nargin == 0, return;end

      N = size(V,1);
      shape.poly = {[1:N,1].'};
      shape.inclusionId = zeros(N,1);

      if nargin>=2
        shape.CSList = {CS};
      else
        shape.CSList = {'notIndexed'};
      end

      shape.phaseId = 1;
      shape.phaseMap = 1;
      shape.id = 1;
      shape.grainSize = 1;

      if isa(V,'grainBoundary') % grain boundary already given
        shape.boundary = V;
      else % otherwise compute grain boundary

        F = [1:N;[2:N 1]].';
        grainId = [zeros(N,1),ones(N,1)];

        shape.boundary = grainBoundary(V,F,grainId,1,...
          1,nan,shape.CSList,shape.phaseMap,1,'noTriplePoints');

      end

    end

    function Vs = get.Vs(shape)
      Vs = shape.boundary.V;
    end

    function theta = get.theta(shape)
      theta = atan2(shape.boundary.V(:,2),shape.boundary.V(:,1));
    end

    function rho = get.rho(shape)
      rho = sqrt(shape.boundary.V(:,2).^2 + shape.boundary.V(:,1).^2);
    end
  end

  methods (Static = true)
    shape = byRhoTheta(rho,theta)
    shape = byFV(F,V,varargin)
  end

end

author: mtex-toolbox

repo: mtex

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


text:
%instrBrickIO Instr interface between MATLAB and the brick (uses the
%instrumentation and control toolbox)
%
% Methods::
%
%  instrBrickIO     Constructor, initialises and opens the instr connection
%  delete           Destructor, closes the instr connection
%
%  open             Open a instr connection to the brick
%  close            Close the instr connection to the brick
%  read             Read data from the brick through instr
%  write            Write data to the brick through instr
%
% Example::
%           instrbtbrick = instrBrickIO(1,'bt','EV3',1)
%           instrwfbrick = instrBrickIO(1,'wf','192.168.1.104',5555,'0016533dbaf5')
%
% Notes::
% - Can connect to the brick through bluetooth or wifi.
% - Uses the instrumentation and control toolbox so only works in windows.

classdef instrBrickIO < BrickIO
    properties
        % connection handle
        handle
        % debug input
        debug = 0;
        % serial connection interface
        serConn;
        % bluetooth brick device name
        btDevice = 'EV3';
        % bluetooth connection channel
        btChannel = 1;
        % wifi brick ip address
        wfAddr = '192.168.1.104';
        % wifi brick tcp port
        wfPort = 5555;
        % wifi brick serial number
        wfSN = '0016533dbaf5';
    end

    methods

        function brickIO = instrBrickIO(varargin)
            %instrBrickIO.btBrickIO Create a instrBrickIO object
            %
            % instrBrickIO = btBrickIO(varargin) is an object which
            % initialises and opens a instr connection between MATLAB
            % and the brick using the instrumentation and control toolbox.
            %
            % Notes::
            % - For bluetooth; instrBrickIO(debug,btDevice,btChannel) where
            % debug is a flag specifying output printing (0 or 1), btDevice
            % is the bluetooth device name of the brick and btChannel is
            % the bluetooth channel.
            % - For wifi; instrBrickIO(debug,wfAddr,wfPort,wfSn) where
            % debug is a flag specifying output printing(0 or 1), wfAddr is
            % the IP address of the brick, wfPort is the TCP connection
            % port and wfSN is the serial number of the brick.

            if nargin > 0
                brickIO.debug = varargin{1};
                brickIO.serConn = varargin{2};
                if brickIO.debug > 0
                    fprintf('instrBrickIO init\n');
                end
            end
            % bluetooth
            % brickIO = instrBrickIO(1,'bt','EV3',1)
            if(strcmp(brickIO.serConn,'bt'))
               brickIO.btDevice = varargin{3};
               brickIO.btChannel = varargin{4};
               % set the bt handle
               brickIO.handle = Bluetooth(brickIO.btDevice,brickIO.btChannel);
               % open the bt handle
               brickIO.open;
            end
            % wifi
            % brickIO = instrBrickIO(1,'wf','192.168.1.104',5555,'0016533dbaf5')
            if(strcmp(brickIO.serConn,'wf'))
               brickIO.wfAddr = varargin{3};
               brickIO.wfPort = varargin{4};
               brickIO.wfSN = varargin{5};
               % set the wf handle
               brickIO.handle = tcpip(brickIO.wfAddr,brickIO.wfPort);
               % connection message to the brick
               wmsg = ['GET /target?sn=' brickIO.wfSN ' VMTP1.0' char(10) 'Protocol: EV3'];
               % open the wf handle
               brickIO.open;
               % send the connection message
               brickIO.write(wmsg);
               pause(0.5);
               % receive the reply
               rmsg = fread(brickIO.handle,16);
               if brickIO.debug > 0
                 fprintf('%s',char(rmsg));
               end
            end
        end

        function delete(brickIO)
            %instrBrickIO.delete Delete the instrBrickIO object
            %
            % delete(brickIO) closes the instr connection handle

            if brickIO.debug > 0
                fprintf('instrBrickIO delete\n');
            end
            % delete the instr handle
            brickIO.close;
        end

        function open(brickIO)
            %instrBrickIO.open Open the instrBrickIO object
            %
            % instrBrickIO.open() opens the instr connection to the brick
            % using fopen.

            if brickIO.debug > 0
                fprintf('instrBrickIO open\n');
            end
            % open the instr handle
            fopen(brickIO.handle);
        end

        function close(brickIO)
            %instrBrickIO.close Close the instrBrickIO object
            %
            % btBrickIO.close() closes the instr connection the brick
            % using fclose.
            if brickIO.debug > 0
                fprintf('instrBrickIO close\n');
            end
            % close the instr handle
            fclose(brickIO.handle);
        end

        function rmsg = read(brickIO)
            %instrBrickIO.read Read data from the instrBrickIO object
            %
            % rmsg = instrBrickIO.read() reads data from the brick through
            % instr (either bluetooth or wifi) via fread and returns the
            % data in uint8 format.
            if brickIO.debug > 0
                fprintf('instrBrickIO read\n');
            end
            % get the number of bytes to be read from the instrCtrl handle
            nLength = fread(brickIO.handle,2);
            % read the remaining bytes
            rmsg = fread(brickIO.handle,double(typecast(uint8(nLength),'uint16')));
            % append the reply size to the return message
            rmsg = uint8([nLength' rmsg']);
        end

        function write(brickIO,wmsg)
            %instrBrickIO.write Write data to the instrBrickIO object
            %
            % instrBrickIO.write(wmsg) writes data to the brick through
            % instr (either bluetooth or wifi).
            %
            % Notes::
            % - wmsg is the data to be written to the brick via instr
            % in uint8 format.
            if brickIO.debug > 0
                fprintf('instrBrickIO write\n');
            end
            fwrite(brickIO.handle,wmsg);
        end
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function value = p37_param ( action, name, value )

%*****************************************************************************80
%
%% P37_PARAM defaults, sets, or gets the parameters for problem 37.
%
%  Discussion:
%
%    In typical use, the user calls
%      p37_param ( 'DEFAULT' )
%    Then, once the variables have been given default values, any
%    value can be retrieved, as for example by:
%      value = p37_param ( 'GET', 'A' )
%    or set to a new value, as in:
%      p37_param ( 'SET', 'A', 2.0 )
%
%    A is the amplitude on the forcing term;
%    K is the damping coefficient.
%    W is the frequency coefficient on the forcing term;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 March 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION:
%    'DEFAULT' sets the variables to default values.
%    'SET' sets the variable chosen by NAME to have the value VALUE.
%    'GET' gets the variable chosen by NAME.
%
%    Input, string NAME, the name of the variable, which should
%    be 'A', 'K', or 'W'.  NAME is ignored if the ACTION is 'DEFAULT'.
%
%    Input, real VALUE, the value to be assigned the variable NAME,
%    if the value of ACTION was 'SET'.  Otherwise VALUE is ignored.
%
%    Output, real VALUE, depends on the input value of ACTION:
%    'DEFAULT', returns the number of variables;
%    'SET', returns the input value of VALUE;
%    'GET', returns the value of the variable NAME.
%
  persistent a;
  persistent k;
  persistent w;

  if ( s_eqi ( action, 'DEFAULT' ) )

    a = 0.3;
    k = 0.2;
    w = 1.0;

    value = 3;

  elseif ( s_eqi ( action, 'SET' ) )

    if ( s_eqi ( name, 'A' ) )
      a = value;
    elseif ( s_eqi ( name, 'K' ) )
      k = value;
    elseif ( s_eqi ( name, 'W' ) )
      w = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P37_PARAM - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized variable name: "%s".\n', name );
      error ( 'P37_PARAM - Fatal error!' );
    end

  elseif ( s_eqi ( action, 'GET' ) )

    if ( s_eqi ( name, 'A' ) )
      value = a;
    elseif ( s_eqi ( name, 'K' ) )
      value = k;
    elseif ( s_eqi ( name, 'W' ) )
      value = w;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P37_PARAM - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized variable name: "%s".\n', name );
      error ( 'P37_PARAM - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P37_PARAM - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action: "%s".\n', action );
    error ( 'P37_PARAM - Fatal error!' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
f = @(x) objfun(x);
A = [-2 1;-1 -1;1 0;0 1];
b = [-1;-2;0;0];
E = [];
e = [];

[xstar,ystar] = zoutendijk(f,A,b,E,e)

repo: Optimization-Theory-and-Algorithm


text:
function [c, s, b, g, active_set] = foopsi_oasisAR1(y, g, lam, smin, optimize_b,...
    optimize_g, decimate, maxIter, tau_range, gmax)
%% Infer the most likely discretized spike train underlying an AR(1) fluorescence trace
% Solves the sparse non-negative deconvolution problem
%  min 1/2|c-y|^2 + lam |s|_1 subject to s_t = c_t-g c_{t-1} >= 0

%% inputs:
%   y:  T*1 vector, One dimensional array containing the fluorescence intensities
%withone entry per time-bin.
%   g:  scalar, Parameter of the AR(1) process that models the fluorescence ...
%impulse response.
%   lam:  scalar, sparsity penalty parameter
%   optimize_b: bool, optimize baseline if True
%   optimize_g: integer, number of large, isolated events to consider for
%       optimizing g
%   decimate: int, decimation factor for estimating hyper-parameters faster
%       on decimated data
%   maxIter:  int, maximum number of iterations
%   active_set: npool x 4 matrix, warm stared active sets
%   gmax:  scalar, maximum value of g
%   tau_range: [tau_min, tau_max]
%% outputs
%   c: T*1 vector, the inferred denoised fluorescence signal at each time-bin.
%   s: T*1 vector, discetized deconvolved neural activity (spikes)
%   b: scalar, fluorescence baseline
%   g: scalar, parameter of the AR(1) process
%   active_set: npool x 4 matrix, active sets

%% Authors: Pengcheng Zhou, Carnegie Mellon University, 2016
% ported from the Python implementation from Johannes Friedrich

%% References
% Friedrich J et.al., NIPS 2016, Fast Active Set Method for Online Spike Inference from Calcium Imaging

%% input arguments
y = reshape(y, [], 1);
T = length(y);

if ~exist('g', 'var') || isempty(g)
    g = estimate_time_constant(y, 1);
end
if ~exist('lam', 'var') || isempty(lam);   lam = 0; end
if ~exist('smin', 'var') || isempty(smin)
    smin = 0;
elseif smin<0
    smin = abs(smin) * GetSn(y);     % use a threshold that is proportional to the noise level
end
if ~exist('optimize_b', 'var') || isempty(optimize_b)
    optimize_b = false;
end
if ~exist('optimize_g', 'var') || isempty(optimize_g)
    optimize_g = 0;
end
if ~exist('decimate', 'var') || isempty(decimate)
    decimate = 1;
else
    decimate = max(1, round(decimate));
end
if ~exist('maxIter', 'var') || isempty(maxIter)
    maxIter = 10;
end
if ~exist('tau_range', 'var') || isempty(tau_range)
   g_range = [0, 1];
else
   g_range = exp(-1./tau_range);
   g = min(max(g, g_range(1)), g_range(2));
end

% change parameters due to downsampling
if decimate>1
    decimate = 1;  %#ok<NASGU>
    disp('to be done');
    %     fluo = y;
    %     y = resample(y, 1, decimate);
    %     g = g^decimate;
    %     thresh = thresh / decimate / decimate;
    %     T = length(y);
end

%% optimize parameters
if ~optimize_b   %% don't optimize the baseline b
    %% initialization
    b = 0;
    [solution, spks, active_set] = oasisAR1(y, g, lam, smin);

    %% iteratively update parameters g
    if  optimize_g     % update g
        [solution, active_set, g, spks] = update_g(y, active_set,lam, smin, g_range);
    end
else
    %% initialization
    b = quantile(y, 0.15);
    [solution, spks, active_set] = oasisAR1(y-b, g, lam, smin);

    %% iteratively update parameters g and b
    for m=1:maxIter
        b = mean(y-solution);
        if  optimize_g     % update g
            if size(active_set, 1) == 0
                break;
            end
            g0 = g;
            if g>gmax  % spike counts are too small. stop
                g = estimate_time_constant(y, 1);
                [solution, spks, active_set] = oasisAR1(y-b, g, lam, smin);
                break;
            end
            [solution, active_set, g, spks] = update_g(y-b, active_set,lam, smin, g_range);
            if abs(g-g0)/g0 < 1e-3  % g is converged
                optimize_g = false;
            end
        else
            break;
        end
    end
end
c = solution;
s = spks;
end
%update the AR coefficient: g
function [c, active_set, g, s] = update_g(y, active_set, lam, smin, g_range)
%% inputs:
%   y:  T*1 vector, One dimensional array containing the fluorescence intensities
%withone entry per time-bin.
%   active_set: npools*4 matrix, previous active sets.
%   lam:  scalar, curret value of sparsity penalty parameter lambda.

%% outputs
%   c: T*1 vector
%   s: T*1 vector, spike train
%   active_set: npool x 4 matrix, active sets
%   g: scalar

%% Authors: Pengcheng Zhou, Carnegie Mellon University, 2016
% ported from the Python implementation from Johannes Friedrich

%% References
% Friedrich J et.al., NIPS 2016, Fast Active Set Method for Online Spike Inference from Calcium Imaging

%% initialization
if ~exist('g_range', 'var') || isempty(g_range)
   g_range = [0, 1];
end
len_active_set = size(active_set, 1);  %number of active sets
y = reshape(y,[],1);    % fluorescence data
maxl = max(active_set(:, 4));   % maximum ISI
c = zeros(size(y));     % the optimal denoised trace

%% find the optimal g and get the warm started active_set
h = [];
g = fminbnd(@rss_g, g_range(1), g_range(2));
yp = y - lam*(1-g);
for m=1:len_active_set
    tmp_h = exp(log(g)*(0:maxl)');   % response kernel
    tmp_hh = cumsum(h.*h);        % hh(k) = h(1:k)'*h(1:k)
    li = active_set(m, 4);
    ti = active_set(m, 3);
    idx = ti:(ti+li-1);
    active_set(m,1) = (yp(idx))'*tmp_h(1:li);
    active_set(m,2) = tmp_hh(li);
end
[c,s,active_set] = oasisAR1(y, g, lam, smin, active_set);

%% nested functions
    function rss = rss_g(g)
        h = exp(log(g)*(0:maxl)');   % response kernel
        hh = cumsum(h.*h);        % hh(k) = h(1:k)'*h(1:k)
        yp = y - lam*(1-g);     % include the penalty term
        for ii=1:len_active_set
            li = active_set(ii, 4);
            ti = active_set(ii, 3);
            idx = ti:(ti+li-1);
            tmp_v = max(yp(idx)' * h(1:li) / hh(li), 0);
            c(idx) = tmp_v*h(1:li);
        end
        res = y-c;
        rss = res'*res;     % residual sum of squares
    end
end

author: zhoupc


text:
function VecInd  = bst_closest(VecGuess, VecRef)
% BST_CLOSEST: Find entries of closest elements between two vectors.
%
% USAGE:  VecInd  = bst_closest(VecGuess, VecRef);
%
% DESCRIPTION:
%     VecGuess is a vector for which one wants to find the closest entries in vector VecRef
%     VecInd is the vector of indices pointing atr the entries in vector VecRef that are the closest to VecWin
%     VecInd is of the length of VecGuess
%
%     In other words, VecRef(VecInd(i)) is the element of VecRef closest to VecGuess(j)
%
%     VecRef and VecGuess do not need to be the same length

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@

if size(VecRef,1) == 1
    VecRef = VecRef';
end

tmp = repmat(VecRef,1,length(VecGuess));
[minn VecInd] = min(abs(repmat(VecGuess,length(VecRef),1) - tmp));

author: brainstorm-tools


text:
function rcov = cov_blockgeom(X,blocksize)
% like cov(), just robust (using the blockwise geometric median)
% The blocksize allows to reduce the memory requirements, at the cost of reduced robustness
% against outliers that occupy fewer samples than the blocksize.

% Copyright (C) Christian Kothe, SCCN, 2013, christian@sccn.ucsd.edu
%
% 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 should have received a copy of the GNU General Public License along with this program; if not,
% write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
% USA

if ~exist('blocksize','var')
    blocksize = 10; end

[n,m] = size(X);
X = bsxfun(@minus,X,median(X));
U = zeros(length(1:blocksize:n),m*m);
for k=1:blocksize
    range = min(n,k:blocksize:(n+k-1));
    U = U + reshape(bsxfun(@times,reshape(X(range,:),[],1,m),reshape(X(range,:),[],m,1)),size(U));
end
rcov = real(reshape(geometric_median(U/blocksize),m,m));

repo: REST


text:
function [logml]=nwmlikgrid(X,n,k,phi0,S0,Sbar,alphabar)

% compute the first determinant part
% create the square root matrix of phi0
% because phi0 is diagonal, this is simply the square root of the diagonal terms of phi0
Fphi=spdiags(diag(phi0).^0.5,0,k,k);
% compute the product
product=Fphi'*X'*X*Fphi;
% compute the eigenvalues of the product
eigenvalues=eig(product);
% now compute the full determinant term
temp1=(-n/2)*log(prod(diag(eye(k)+diag(eigenvalues))));

% compute the second determinant part
% create the square root matrix of inv(S0)
invS0=spdiags(1./diag(S0),0,n,n);
Fs=spdiags(diag(invS0).^0.5,0,n,n);
% compute the summation
summ=Fs'*(Sbar-S0)*Fs;
% compute the eigenvalues of the summation
eigenvalues=eig(summ);
% now compute the full determinant term
temp2=(-alphabar/2)*log(prod(diag(eye(n)+diag(eigenvalues))));

% compute the marginal likelihood
logml=real(temp1+temp2);

author: european-central-bank

repo: BEAR-toolbox


text: %% S2Kernel


text:
function h = thermometer(varargin)
% Create a Thermometer plot.
% THERMOMETER(DATA, MAX) - Create a thermometer
%    where DATA is a set of 1D data to plot in the thermometer.
%    Each entry in DATA adds onto the previous number, thus:
%    [ 3 4 5 ] will show the '3' between 0 and 3, and '4' between 3
%    and 7, and '5' between 7 and 12.
% THERMOMETER(DATA, MIN, MAX)
%    MIN and MAX represent the minimum and maximum thermometer
%    values.
%    Values in DATA add onto MIN. Thus, if MIN is 10, and DATA is
%    [ 3 4 ] then the '3' is shown between 10 and 13, and the '4'
%    between 13 and 17.
% THERMOMETER(..., { 'name1', 'name2' })
%    Area labels.
% THERMOMETER(..., { GOAL 'label' })
%    At a specific GOAL, add LABEL.
%    The GOAL value is absolute, and does not change based on MIN
%    or MAX.
%
% The thermometer code has been tweeked to look good printed from
% painters.  Using Zbuffer or OpenGL will make it looks strange.
%
% To adjust the tick labels, modify GCA's yticks.
%
% To make the zone labels alternate sides, uncomment the section
% that contains the word 'left'.
%
% Example:
% thermometer([ 3 4 ], 10, { 'tasks' 'bugs' }, { 8 'Almost There' })

% History:
%  Version 1.2: Update doc to better describe data.
%
%  Version 1.1:  Fixed problem using color-order for shading
%  elements so more than 7 entries can be used.
%
%
% Eric Ludlam
% Copyright (c) 2006 The MathWorks Inc.

  data = varargin{1};

  if any(data < 0)
      error('DATA Values for thermometer must be > 0.');
  end

  if nargin >= 3 && isnumeric(varargin{3})
    inputy = [ varargin{2} varargin{3} ];
    base = 4;
  else
    inputy = [ 0 varargin{2} ];
    base = 3;
  end

  names = {};
  goals = {};

  while nargin >= base

    nexta = varargin{base};

    if ~iscell(nexta)
      error('Arguments after MAX must be cell arrays');
    end

    if ischar(nexta{1})

      names = cell(size(data,1),1);
      names(1:length(varargin{base})) = varargin{base};

    elseif isnumeric(nexta{1})

      goals = nexta;

    else

      error('Unknown cell array format for thermometer.', nexta); %#ok

    end

    base = base + 1;
  end

  lw = 2;
  ax = newplot;

  set(ax,'color','white');
  total = inputy(2) - inputy(1);
  yl = [ inputy(1)-(total/15) inputy(2) ];
  set(ax,'ylim', yl)
  set(ax,'xlim', [ 0 1 ])
  set(ax,'xtick', [] )
  set(ax,'plotboxaspectratio',[ 1 20 1 ]);
  set(ax,'box','on');
  set(ax,'tickdir','out');
  set(ax,'ygrid','on');
  set(ax,'gridlinestyle','-');
  set(ax,'box','off');

  co = get(gca,'colororder');

  l = get(gca,'loose');
  l(1) = l(end);
  set(gca,'loose',l);

  % Draw a circle at bottom of thermometer

  rectangle('position',[-.5 inputy(1)-(total/15*1.5) 2 total/15*1.6 ],...
            'curvature',[ 1 1 ],...
            'facecolor',co(1,:),'edgecolor','k',...
            'clipping','off',...
            'linewidth',3);

  % Data Patchs

  % One number being plotted

  surface('xdata', [ 0 1 ] , 'ydata', [ yl(1) data(1)+inputy(1) ], ...
          'zdata', [ 1 1 ; 1 1 ],  ...
          'facecolor', co(1,:) ,'edgecolor','none');

  if ~isempty(names)
    yzonelabel_t(inputy(1), data(1)+inputy(1), names{1}, 'right');
  end

  bottom = data(1)+inputy(1);

  % A stacked data plot
  for i = 2:length(data)

    cidx = i-1;
    color = co(mod(cidx,size(co,1))+1,:);

    surface('xdata', [ 0 1 ] , 'ydata', [ bottom bottom+data(i) ], ...
            'zdata', [ 0 0 ; 0 0 ],  ...
            'facecolor', color ,'edgecolor','none');
    %if ~mod(i,2)
    %  lr = 'left';
    %else
    lr = 'right';
    %end

    if ~isempty(names) && length(names) >= i
      yzonelabel_t(bottom, bottom+data(i), names{i}, lr);
    end

    bottom = bottom+data(i);

  end

  while ~isempty(goals)

    pointlabel(goals{1}, goals{2}, false);

    goals = goals(3:end);

  end

    % Put a new box onto the axes.
  line('xdata',[0 0 1 1],'ydata', [ inputy(1) yl(2) yl(2) inputy(1) ],...
       'zdata',[1 1 1 1],...
       'color','k','linewidth',lw,'clipping','off');

  if nargout > 0
      h = ax;
  end

end

function pointlabel(yval, str, left)
% Create a fancy label at yval

  if left
    xdata = [ 0 -.7 ];
    txtpos = -.8;
    ha = 'right';
  else
    xdata = [ 1 1.7 ];
    txtpos = 1.8;
    ha = 'left';
  end
  ydata = [ yval yval ];

  line(xdata,ydata,'color','k','clipping','off');
  line(1,yval,'color','k',...
       'markersize',6,...
       'marker','diamond','markerfacecolor','k');
  line(0,yval,'color','k',...
       'markersize',6,...
       'marker','diamond','markerfacecolor','k');
  line([0 1],ydata,'color','k','linestyle',':');
  text(txtpos,yval,str,'horizontalalign',ha,...
       'verticalalign','middle');

end

function yzonelabel_t(start, fin, str, side)
% Create a small label thingy between START and END
% STR is an optional string to add to the label.  The text
% representing the total Y value encompassed is also included.
% SIDE is a string, either 'left', or 'right'.

  xl = xlim;
  xr = xl(2);
  xl = xl(1);
  factor = abs(xl-xr);

  switch side
   case 'left'
    in = xl-1.2*factor;
    out = xl-2*factor;
    xdata = [ in out out in ];
    txtpos = xl-2.1*factor;
    rot = 90;
   case 'right'
    in = xr + .5*factor;
    out = xr + .7*factor;
    xdata = [ in out out in ];
    txtpos = xr + .8*factor;
    rot = -90;
  end
  ydata = [ start start fin fin ];

  mid = (start + fin) / 2;

  numstr = num2str(fin-start);

  if isempty(str)
    newstr = numstr;
  else
    newstr = [ numstr char(10) str ];
  end

  l = line(xdata,ydata,'color','k');
  t = text(txtpos, mid, newstr,'rotation',rot,...
           'fontsize',10,...
           'horizontalalign','center',...
           'verticalalign','bottom');
  set([ l t ], 'xliminclude','off',...
               'clipping','off');

end

author: Sable

repo: mcbench-benchmarks


text:
%SeDuMi installation script
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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 should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

% Modified J.Currie April 2013
%
% My build platform:
% - Windows 7 SP1 x64
% - Visual Studio 2010
% - Intel Math Kernel Library

%List of Files to Compile (and dependencies)
sourcelist={...
    'bwblkslv.c sdmauxFill.c sdmauxRdot.c',...
    'choltmpsiz.c',...
    'cholsplit.c',...
    'dpr1fact.c auxfwdpr1.c sdmauxCone.c  sdmauxCmp.c sdmauxFill.c sdmauxScalarmul.c sdmauxRdot.c blkaux.c',...
    'symfctmex.c symfct.c',...
    'ordmmdmex.c ordmmd.c',...
    'quadadd.c',...
    'eigK.c sdmauxCone.c sdmauxRdot.c',...
    'sqrtinv.c sdmauxCone.c',...
    'givensrot.c auxgivens.c sdmauxCone.c',...
    'urotorder.c auxgivens.c sdmauxCone.c sdmauxTriu.c sdmauxRdot.c',...
    'psdframeit.c reflect.c sdmauxCone.c sdmauxRdot.c sdmauxTriu.c sdmauxScalarmul.c',...
    'psdinvjmul.c reflect.c sdmauxCone.c sdmauxRdot.c sdmauxTriu.c sdmauxScalarmul.c blkaux.c',...
    'bwdpr1.c sdmauxCone.c sdmauxRdot.c',...
    'fwdpr1.c auxfwdpr1.c sdmauxCone.c sdmauxScalarmul.c',...
    'fwblkslv.c sdmauxScalarmul.c',...
    'qblkmul.c sdmauxScalarmul.c',...
    'blkchol.c blkchol2.c sdmauxFill.c sdmauxScalarmul.c',...
    'vecsym.c sdmauxCone.c',...
    'qrK.c sdmauxCone.c sdmauxRdot.c sdmauxScalarmul.c',...
    'finsymbden.c sdmauxCmp.c',...
    'symbfwblk.c',...
    'statsK.c sdmauxCone.c',...
    'whichcpx.c sdmauxCone.c',...
    'eyeK.c sdmauxCone.c',...
    'ddot.c sdmauxCone.c sdmauxRdot.c sdmauxScalarmul.c',...
    'makereal.c sdmauxCone.c sdmauxCmp.c',...
    'partitA.c sdmauxCmp.c',...
    'getada1.c sdmauxFill.c',...
    'getada2.c sdmauxCone.c sdmauxRdot.c sdmauxFill.c',...
    'getada3.c spscale.c sdmauxCone.c sdmauxRdot.c sdmauxScalarmul.c sdmauxCmp.c',...
    'adendotd.c sdmauxCone.c',...
    'adenscale.c',...
    'extractA.c',...
    'vectril.c sdmauxCone.c sdmauxCmp.c',...
    'qreshape.c sdmauxCone.c sdmauxCmp.c',...
    'sortnnz.c sdmauxCmp.c',...
    'iswnbr.c',...
    'incorder.c',...
    'findblks.c sdmauxCone.c sdmauxCmp.c',...
    'invcholfac.c triuaux.c sdmauxCone.c sdmauxRdot.c sdmauxTriu.c sdmauxScalarmul.c blkaux.c',...
    };

sourcelist = {'psdscale.c triuaux.c sdmauxCone.c sdmauxRdot.c sdmauxTriu.c sdmauxScalarmul.c blkaux.c'...
              'psdfactor.c triuaux.c sdmauxCone.c sdmauxRdot.c sdmauxTriu.c sdmauxScalarmul.c blkaux.c'...
              'psdeig.c triuaux.c sdmauxCone.c sdmauxRdot.c sdmauxTriu.c sdmauxScalarmul.c blkaux.c'...
              'psdinvscale.c triuaux.c sdmauxCone.c sdmauxRdot.c sdmauxTriu.c sdmauxScalarmul.c blkaux.c'};

%Generate function list
funclist = cell(size(sourcelist));
for i = 1:length(sourcelist)
    files = regexp(sourcelist{i},' ','split');
    ind = strfind(files{1},'.');
    funclist{i} = files{1}(1:ind-1);
end

%Clear existing functions in memory
clear(funclist{:});

% Modify below function if it cannot find Intel MKL on your system.
mkl_link = opti_FindMKL();

fprintf('\n------------------------------------------------\n');
fprintf('SeDuMi MEX FILES INSTALL\n\n');

%Get Headers and Libraries
post = [' -IInclude -DPC ' mkl_link];

%CD to Source Directory
cdir = cd;
cd '../Downloads/sedumi';

%Setup compiler command
pre = 'mex -largeArrayDims ';

try
%Compile & Move each file
for i = 1:length(funclist)
    fprintf('Compiling %s.%s....',funclist{i},'c');
    eval([pre sourcelist{i} post])
    fprintf('Done!\n');
end
catch ME
    cd(cdir);
    error('opti:sedumi','Error Compiling SeDuMi!\n%s',ME.message);
end
cd(cdir);
fprintf('\nSeDuMi Compilation Complete\n');
fprintf('------------------------------------------------\n');

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function f = mat2cell(f, varargin)
%MAT2CELL   Wrap a DELTAFUN in a cell.
%   DELTAFUN objects are not array-valued, so G = MAT2CELL(F) simply wraps F
%   in a cell array and is equivalent to G = {F}.

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

if ( any(cellfun(@(v) v ~= 1, varargin)) )
    error('CHEBFUN:DELTAFUN:mat2cell:inputs', ...
        'DELTAFUN objects are not array-valued, so M and n must equal 1.');
end

f = {f};

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function wedge_grid_plot ( n, ng, g, header )

%*****************************************************************************80
%
%% WEDGE_GRID_PLOT sets up a GNUPLOT plot of a unit wedge grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 August 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of subintervals.
%
%    Input, integer NG, the number of nodes.
%
%    Input, real G(3,NG), the grid point coordinates.
%
%    Input, string HEADER, the header for the files.
%

%
%  Create the vertex file.
%
  [ v1, v2, v3, v4, v5, v6 ] = wedge_vertices ( );

  vertex_filename = strcat ( header, '_vertices.txt' );
  vertex_unit = fopen ( vertex_filename, 'wt' );

  fprintf ( vertex_unit, '%g  %g  %g\n', v1(1), v1(2), v1(3) );
  fprintf ( vertex_unit, '%g  %g  %g\n', v2(1), v2(2), v2(3) );
  fprintf ( vertex_unit, '%g  %g  %g\n', v3(1), v3(2), v3(3) );
  fprintf ( vertex_unit, '%g  %g  %g\n', v1(1), v1(2), v1(3) );
  fprintf ( vertex_unit, '\n' );

  fprintf ( vertex_unit, '%g  %g  %g\n', v4(1), v4(2), v4(3) );
  fprintf ( vertex_unit, '%g  %g  %g\n', v5(1), v5(2), v5(3) );
  fprintf ( vertex_unit, '%g  %g  %g\n', v6(1), v6(2), v6(3) );
  fprintf ( vertex_unit, '%g  %g  %g\n', v4(1), v4(2), v4(3) );
  fprintf ( vertex_unit, '\n' );

  fprintf ( vertex_unit, '%g  %g  %g\n', v1(1), v1(2), v1(3) );
  fprintf ( vertex_unit, '%g  %g  %g\n', v4(1), v4(2), v4(3) );
  fprintf ( vertex_unit, '\n' );

  fprintf ( vertex_unit, '%g  %g  %g\n', v2(1), v2(2), v2(3) );
  fprintf ( vertex_unit, '%g  %g  %g\n', v5(1), v5(2), v5(3) );
  fprintf ( vertex_unit, '\n' );

  fprintf ( vertex_unit, '%g  %g  %g\n', v3(1), v3(2), v3(3) );
  fprintf ( vertex_unit, '%g  %g  %g\n', v6(1), v6(2), v6(3) );
  fprintf ( vertex_unit, '\n' );

  fclose ( vertex_unit );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Created vertex file "%s"\n', vertex_filename );
%
%  Create the node file.
%
  node_filename = strcat ( header, '_nodes.txt' );
  node_unit = fopen ( node_filename, 'wt' );
  for j = 1 : ng
    fprintf ( node_unit, '%g  %g  %g\n', g(1:3,j) );
  end
  fclose ( node_unit );
  fprintf ( 1, '  Created node file "%s"\n', node_filename );
%
%  Create the command file.
%
  command_filename = strcat ( header, '_commands.txt' );
  command_unit = fopen ( command_filename, 'wt' );
  fprintf ( command_unit, '# %s\n', command_filename );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, '# Usage:\n' );
  fprintf ( command_unit, '#  gnuplot < %s\n', command_filename );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, 'set term png\n' );
  plot_filename = strcat ( header, '.png' );
  fprintf ( command_unit, 'set output "%s"\n', plot_filename );
  fprintf ( command_unit, 'set xlabel "<--- X --->"\n' );
  fprintf ( command_unit, 'set ylabel "<--- Y --->"\n' );
  fprintf ( command_unit, 'set zlabel "<--- Z --->"\n' );
  fprintf ( command_unit, 'set title "%s"\n', header );
  fprintf ( command_unit, 'set grid\n' );
  fprintf ( command_unit, 'set key off\n' );
  fprintf ( command_unit, 'set view 80, 85\n' );
  fprintf ( command_unit, 'set style data lines\n' );
  fprintf ( command_unit, 'set timestamp\n' );
  fprintf ( command_unit, 'splot "%s" with lines lw 3, \\\n', vertex_filename );
  fprintf ( command_unit, '      "%s" with points pt 7 lt 0\n', node_filename );
  fclose ( command_unit );

  fprintf ( 1, '  Created command file "%s"\n', command_filename );

  return
end


text:
function Z=SoftThreshold_S( A, thres )

tmp=A;
S=sign(tmp);
tmp=(abs(tmp)-thres);
tmp(tmp<=0)=0;
Z=(S.*tmp);

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
function  cmap = GetColormap(clrmap_name,numK)
numK = double(numK);
data_dir = GetCurrentDataDir();
load(fullfile(data_dir,'CustomColormaps.mat'));
if strcmp(clrmap_name,'hsv_new'),
    cmap64 = hsv_new;
    cmap = InterpColormap(cmap64,numK);
elseif strcmp(clrmap_name,'jet'),
    cmap = flipud(jet(numK));
elseif strcmp(clrmap_name,'greedy_hsv'),
%     cmap64 = greedy_hsv;
%     cmap = InterpColormap(cmap64,numK);
% choose random colors!
%  cmap64 = greedy_hsv;
rng(1);
H = horzcat(rand(numK,1),rand(numK,1),0.5+0.5*rand(numK,1));
cmap = hsv2rgb(H);

elseif strcmp(clrmap_name,'hsv_old'),
    cmap = hsv(round(numK*1.1));
%     cmap64 = hsv_old;
% %     same as:
%     temp = hsv(70);
%     cmap64 = temp(1:64,:);

else
    cmap64 = hsv_new;
    cmap = InterpColormap(cmap64,numK);
end

% cmap = [1,0.956763267517090,0.00606909440830350;0.393025994300842,1,0.00697167776525021;0,1,1;0.402252703905106,0.00566019536927342,1;1,0,0;1,0.477354049682617,0.00717086857184768];

cmap(cmap<0) = 0;
end

function cmap = InterpColormap(cmap64,numK) % this doesn't work well for very small numK - too close to edge
k = double(numK);
c1 = interp1(1:64,cmap64(:,1),linspace(1,64,k),'pchip');
c2 = interp1(1:64,cmap64(:,2),linspace(1,64,k),'pchip');
c3 = interp1(1:64,cmap64(:,3),linspace(1,64,k),'pchip');
% c1 = interp1(1:64,cmap64(:,1),linspace(1,64,k),'spline');
% c2 = interp1(1:64,cmap64(:,2),linspace(1,64,k),'spline');
% c3 = interp1(1:64,cmap64(:,3),linspace(1,64,k),'spline');
cmap = [c1',c2',c3'];
end

author: xiuyechen

repo: FishExplorer


text:
%% Load data
clear; close all; clc;

% 500-1000
% swf10 = load('swf_500_1000_10_dataset3.mat');
% swf50 = load('swf_500_1000_50_dataset3.mat');
% swf100 = load('swf_500_1000_100_dataset3.mat');
% msckf5 = load('msckf_500_1000_min5_maxInf');
% msckf10 = load('msckf_500_1000_min10_max50');
% msckf20 = load('msckf_500_1000_min20_max100');
% imu = load('imu_500_1000.mat');
% kStart = 500; kEnd = 1000;

% 1215-1715
swf10 = load('swf_1215_1715_10_dataset3.mat');
swf50 = load('swf_1215_1715_50_dataset3.mat');
swf100 = load('swf_1215_1715_100_dataset3.mat');
msckf5 = load('msckf_1215_1715_min5_maxInf');
msckf10 = load('msckf_1215_1715_min10_max50');
msckf20 = load('msckf_1215_1715_min20_max100');
imu = load('imu_1215_1715.mat');
kStart = 1215; kEnd = 1715;

%% Compute RMSE
swf10_trans_rmse = sqrt(mean(swf10.swf_trans_err.^2,1));
swf10_rot_rmse = sqrt(mean(swf10.swf_rot_err.^2,1));

swf50_trans_rmse = sqrt(mean(swf50.swf_trans_err.^2,1));
swf50_rot_rmse = sqrt(mean(swf50.swf_rot_err.^2,1));

swf100_trans_rmse = sqrt(mean(swf100.swf_trans_err.^2,1));
swf100_rot_rmse = sqrt(mean(swf100.swf_rot_err.^2,1));

msckf5_trans_rmse = sqrt(mean(msckf5.msckf_trans_err.^2,1));
msckf5_rot_rmse = sqrt(mean(msckf5.msckf_rot_err.^2,1));

msckf10_trans_rmse = sqrt(mean(msckf10.msckf_trans_err.^2,1));
msckf10_rot_rmse = sqrt(mean(msckf10.msckf_rot_err.^2,1));

msckf20_trans_rmse = sqrt(mean(msckf20.msckf_trans_err.^2,1));
msckf20_rot_rmse = sqrt(mean(msckf20.msckf_rot_err.^2,1));

imu_trans_rmse = sqrt(mean(imu.msckf_trans_err.^2,1));
imu_rot_rmse = sqrt(mean(imu.msckf_rot_err.^2,1));

%% Plot stuff
figure(1); clf;
fontSize = 14;
lineWidth = 2;
pos = [200,200,640,400];
xLim = [kStart,kEnd];
k = kStart:kEnd;
xticks = linspace(kStart,kEnd,11);

set(gcf, 'Position', pos);

% Translational RMSE
subplot(2,1,1);
plot(k,imu_trans_rmse, '-k', 'LineWidth', lineWidth); hold on;
plot(k,msckf5_trans_rmse, '-b', 'LineWidth', lineWidth);
plot(k,msckf10_trans_rmse, '--b', 'LineWidth', lineWidth);
plot(k,msckf20_trans_rmse, '-.b', 'LineWidth', lineWidth);
plot(k,swf10_trans_rmse, '-g', 'LineWidth', lineWidth);
plot(k,swf50_trans_rmse, '--g', 'LineWidth', lineWidth);
plot(k,swf100_trans_rmse, '-.g', 'LineWidth', lineWidth);
xlim(xLim);
set(gca,'XTick',xticks);
legend('IMU Only','MSCKF 5-Inf', 'MSCKF 10-50', 'MSCKF 20-100', 'SWF 10', 'SWF 50', 'SWF 100','Location', 'northwest');
title('Window Size Comparison');
ylabel('Trans. RMSE (m)')
set(gca,'FontSize',fontSize)
grid minor; box on;

% Rotational RMSE
subplot(2,1,2);
plot(k,imu_rot_rmse, '-k', 'LineWidth', lineWidth); hold on;
plot(k,msckf5_rot_rmse, '-b', 'LineWidth', lineWidth);
plot(k,msckf10_rot_rmse, '--b', 'LineWidth', lineWidth);
plot(k,msckf20_rot_rmse, '-.b', 'LineWidth', lineWidth);
plot(k,swf10_rot_rmse, '-g', 'LineWidth', lineWidth);
plot(k,swf50_rot_rmse, '--g', 'LineWidth', lineWidth);
plot(k,swf100_rot_rmse, '-.g', 'LineWidth', lineWidth);
xlim(xLim);
set(gca,'XTick',xticks);
legend('IMU Only','MSCKF 5-Inf', 'MSCKF 10-50', 'MSCKF 20-100', 'SWF 10', 'SWF 50', 'SWF 100','Location', 'northwest');
ylabel('Rot. RMSE (Axis-Angle)');
xlabel('Timestep');
set(gca,'FontSize',fontSize)
grid minor; box on;

%% Export figure
fileName = sprintf('RMSE-Comparison-WindowSize-%d-%d.pdf',kStart,kEnd);
export_fig(gcf, fileName, '-transparent');

repo: MSCKF


text:
%% DEMOGRADIENT  Short demonstration of gradients
%

%% Some sample applications of the gradient toolbox
% Gradients implement automatic differentiation in forward mode, which is
% conveniently to implement using the Matlab operator concept.
%

%% Initialization of gradients
% In order to use automatic differentiation, the independent variables need
% to be identified and values have to be assigned. This is performed by
% the function "gradientinit", for example

format compact short _
u = gradientinit([ -3.1 ; 4e-3 ])

%%
% The total size of the input is the number of independent variables,
% in the example 2, hence u represents a column vector of length 2 and
% defines two independent variables u(1) and u(2) with gradients [1 0]
% and [0 1], respectively.

%% Operations between gradients
% If at least one operand is of type gradient, operations are executed as
% gradient operations.
% For example,

x = gradientinit(3.5);
y = exp(3*x-sqrt(x))

%%
% For f(x):=exp(3*x-sqrt(x)), the result y contains in y.x the function value f(3.5)
% and in y.dx the derivative f'(3.5):
%

y.x, y.dx

%% Complex arguments
% When evaluating the expression for another argument, use the same
% statement as before with new values.

x = gradientinit(-3.5+.2i);
y = exp(3*x-sqrt(x))

%% Access to the gradient
% The principle works for functions in several unknowns the same way. Define, for
% example, the following function from R^3->R^3 :

f = @(x)( [ -2*x(1)*x(2)+4*x(3)^2 ; sin(x(2))/sqrt(pi-x(1)) ; atan(x(2)-x(3)) ] )
f([1.5;-1;0.7])

%%
% then the function value and gradient at [1.5;-1;0.7] is computed by

y = f(gradientinit([1.5;-1;0.7]))

%%
% where y.x contains the function value and y.dx the gradient, which is
% in this case the Jacobian. The gradient with respect the third unknown
% x(3) can be accessed by

y.dx(3,:)

%%
% However, it is recommended to use

y(3).dx

%%
% that is not to access the components of the gradient (Jacobian) but the
% gradient of the component. The advantage is visible when redefining the input function
% as a row vector:

f = @(x)( [ -2*x(1)*x(2)+4*x(3)^2  sin(x(2))/sqrt(pi-x(1))  atan(x(2)-x(3)) ] )
f([1.5;-1;0.7])

%%
% Then the "Jacobian" is a three-dimensional array
% because the gradient is always stored in the "next" dimension:

y = f(gradientinit([1.5;-1;0.7]))

%%
% It is problematic to access the components of y.dx, while accessing the
% gradient of the component works as expected:

y(3).dx

%% An example in one unknown: The Gamma function
% According to Stirling's formula it is for u -> inf,
%
%                         1      1       139        571
%   Gamma(u) ~ C * ( 1 + --- + ----- - ------- - --------- + ... )
%                        12u       2         3           4
%                              288u    51840u    2488320u
%
% with
%
%        -u  u-0.5
%   C = e   u      sqrt(2*pi) .
%
% The following function evaluates Stirling's formula. It is also
% suited for vector input.
%
%   function y = g(u)
%       C = exp(-u) .* ( u.^(u-0.5) ) * sqrt(2.0*pi) ;
%       v = (((( -571.0/2488320.0 ./ u - 139.0/51840.0 ) ./ u ...
%            + 1.0/288.0) ./ u ) + 1.0/12.0 ) ./ u + 1.0;
%       y = C .* v;
%
% A corresponding inline function is

format long e
g = @(u) ( ( exp(-u) .* ( u.^(u-0.5) ) * sqrt(2.0*pi) ) .* ...
           ( (((( -571.0/2488320.0 ./ u - 139.0/51840.0 ) ./ u ...
             + 1.0/288.0) ./ u ) + 1.0/12.0 ) ./ u + 1.0 ) )
u = [ 3.5 61 5 ]
g(u)

%% The inverse Gamma function
% Next we calculate the inverse Gamma function. For example, compute u such
% that g(u) = 100. Consider the following simple Newton procedure with starting
% value u=5.

u = gradientinit(5);
uold = u;
k = 0;
while abs(u.x-uold.x) > 1e-12*abs(u.x) | k < 1
  uold = u;
  k = k+1;
  y = g(u) - 100;
  u = u - y.x/y.dx;
end
k
u.x
g(u.x)

%%
% Due to the approximation error in Stirling''s formula, about six figures are correct.

%% The inverse Gamma function with complex arguments
% The same is possible for complex arguments. We use the same Gamma function
% and the same Newton procedure except that some u is searched with
% g(u) = 100 + 100i. We use the same starting value u=5.
%

 u = gradientinit(5);
 uold = u;
 k = 0;
 while abs(u.x-uold.x) > 1e-12*abs(u.x) | k < 1
  uold = u;
  k = k+1;
  y = g(u) - 100 - 100i;
  u = u - y.x/y.dx;
 end
 k
 u.x
 g(u.x)

%%
% Due to approximation error in Stirling''s formula, about six figures are correct.

%% Automatic differentiation with several unknowns
% Automatic differentiation with several unknowns works the same way.
% Consider the following example by Broyden:
%
%                               .5*sin(x1*x2) - x2/(4*pi) - x1/2  =  0
% (1-1/(4*pi))*(exp(2*x1)-exp(1)) + exp(1)*x2/pi - 2*exp(1)*x1 )  =  0
%
% with initial approximation [ .6 ; 3 ] and one solution [ .5 ; pi ].
% The following inline function evaluates Broyden's function.

f = @(x) ( [ .5*sin(x(1)*x(2)) - x(2)/(4*pi) - x(1)/2 ; ...
             (1-1/(4*pi))*(exp(2*x(1))-exp(1)) + exp(1)*x(2)/pi - 2*exp(1)*x(1) ] )

%% Solution of a nonlinear system
% The nonlinear system defined by Broyden's function is solved by Newton's procedure as follows:

 x = gradientinit([ .6 ; 3 ]);
 for i=1:5
  y = f(x);
  x = x - y.dx\y.x;
 end
 x

%%
% For simplicity, we omitted the stopping criterion (see above).
% Here, y.dx is the Jacobian, y.x the function value at x.x, and -y.dx\y.x
% is the correction obtained by the (approximate) solution of a linear system.
%

%% Verified solution of the nonlinear system
% For verified solution of the nonlinear system, we need a correct definition
% of the function. The main point is to make sure that a function evaluation with
% interval argument computes an inclusion of the function value. So first the
% transcendental number pi has to be replaced by an interval containing pi, for example

cPi = midrad(3.141592653589793,1e-15)

%%
% Second, Broyden's function contains exp(1), which would be computed in pure
% floating-point without extra care. This can be cured using exp(intval(1)).
%
% However, a new problem arises.
% When replacing "pi" in the function by "cPi" and 1 by intval(1), the function is
% _always_ evaluated in interval arithmetic; a pure floating point iteration
% is no longer possible.
%
% To solve this problem, we have to know the type of
% the incoming unknown "x". If "x" is double, replace "cPi" and intval(1) by its midpoint,
% if "x" is an interval, use "cPi" and intval(1) as is. This is done as follows.
%
%  function  y = f(x)
%    y = x;
%    c1 = typeadj( 1 , typeof(x) );
%    cpi = typeadj( midrad(3.14159265358979323,1e-16) , typeof(x) );
%    y(1) = .5*sin(x(1)*x(2)) - x(2)/(4*cpi) - x(1)/2;
%    y(2) = (1-1/(4*cpi))*(exp(2*x(1))-exp(c1)) + exp(c1)*x(2)/cpi - 2*exp(c1)*x(1);
%
% This code is implemented in the function test.m .

%% Real function evaluation
% Consider the following two function evaluations. First, f(x) is evaluated
% for real argument:

x = [ .6 ; 3 ];
test(x)

%% Interval function evaluation
% Second, f(x) is evaluated with interval argument:

x = [intval('.6') ; 3 ]
y = test(x)

%%
% The mathematical statement is the following. First, x is an interval vector
% such that x(1) is an inclusion of 0.6 and x(2)=3. Second, cPi is an interval
% containg the transcendental number pi. Third, y is an interval vector
% containing the exact value of Broyden's function evaluated at [ .6 ; 3 ].

%% Interval gradient function evaluation
% Finally, we may define the interval x to be of type gradient:

x = gradientinit([intval('.6') ; 3 ])
Y = test(x)

%%
% The mathematical statement is that Y is an interval vector such that Y.x
% contains the exact value of Broyden's function evaluated at [ .6 ; 3 ], and
% Y.dx is an interval matrix containing the Jacobian of Broyden's function
% evaluated at [ .6 ; 3 ].

%% Verified solution of the nonlinear system with Broyden's function
% The nonlinear system with Broyden's function and the given starting value
% [ .6 ; 3 ] can be solved with verification by

Y = verifynlss(@test,[ .6 ; 3 ])

%%
% The first parameter gives the name of the function such that test(x)
% evaluates the function at "x". The result vector Y is verified to contain
% a real vector X such that f(X)=0. This solution X of the nonlinear system
% is proved to be unique within Y. This statement is mathematically true,
% it is taken care of all procedural, approximation and rounding errors.
%
% It follows that an inclusion is not possible is roots are very close together:
% Since uniqueness of the root is proved, an inclusion is only possible if roots
% can be separated. An escape of that is described in DEMOSLOPE.

%% Verified solution of a nonlinear system with sparse gradients
% Up to now we considered only toy examples to explain how the nonlinear
% system solver works. For a larger example consider the following example
% proposed by Abbot and Brent, which is implemented in the function test.
%

%%
%  function y = test(x);
%  % Abbot/Brent     3 y" y + y'^2 = 0;    y(0)=0; y(1)=20;
%  % approximation   10*ones(n,1)
%  % solution        20*x^.75
%    y = x;
%    n = length(x); v=2:n-1;
%    y(1) = 3*x(1)*(x(2)-2*x(1)) + x(2)*x(2)/4;
%    y(v) = 3*x(v).*(x(v+1)-2*x(v)+x(v-1)) + (x(v+1)-x(v-1)).^2/4;
%    y(n) = 3*x(n).*(20-2*x(n)+x(n-1)) + (20-x(n-1)).^2/4;
%
% An inclusion of the solution for 1000 unknowns is computed by
%

format short
sparsegradient(50)
n = 1000;
tic
x = verifynlss(@test,10*ones(n,1));
toc
max(relerr(x))

%% Verified solution of a nonlinear system with full gradients
% Here we specified that gradients with 50 unknowns and more are stored
% in sparse mode. This is the defauls when calling "sparsegradient".
%
% Forcing gradients to use full storage results in a significantly increase
% of computing time.

sparsegradient(inf)
n = 1000;
tic
x = verifynlss(@test,10*ones(n,1));
toc
max(relerr(x))

%% Verified solution of a nonlinear system with 5000 uknowns
% Note that the inclusion is of high accuracy. The results for a larger
% nonlinear system with 5000 unknowns is as follows.

sparsegradient(0)
n = 5000;
tic
x = verifynlss(@test,10*ones(n,1));
toc
max(relerr(x))

%% Non-differentiable functions
% The given function need not be differentiable everywhere. Consider, for example,

f = inline('abs(x)')
f(gradientinit(infsup(-.1,2)))

%% Verified solution of non-differentiable functions
% The inclusion of a root is searched for near the given approximation. Consider

f = vectorize(inline('x*sinh(x)-3*exp(abs(x)-.5)+x*cos(x+1)+2*cosh(x)'))
close
x=linspace(-1,2.3);
plot(x,f(x),x,0*x)

%%
% It seems there are three roots. Note that the function contains abs(x). Indeed,

format long
verifynlss(f,-.5)
verifynlss(f,.5)
verifynlss(f,2)

%%
% there are three roots, and the inlcusions are of high accuracy. One can
% also calculate an inclusion of multiple roots. Since the problem is ill-posed,
% it has t be regularized. Consider

X = verifynlss2(@(x)(sin(x)-1),1.5)

%%
% It is proved that there exists some parameter e in X(2) such that the function
% g(x):=f(x)-e has a truely double root in X(1). Note the accuracy of the inclusions.
% It looks like the inclusion of e is a true zero. However, this is due to the "_"-format:
% add +/-1 to the last visible digit produces a valid inclusion:

format long e infsup
X

%%
% The function "verifynlss2" is applicable to multivariate functions as well.
%
% For more details, see "help verifynlss" or "help verifynlss2" or
%
%  S.M. Rump: Verification methods: Rigorous results using floating-point arithmetic.
%    Acta Numerica, 19:287-449, 2010.
%
% to be downloaded from "www.ti3.tuhh.de/rump" and the literature cited over there.

%% Enjoy INTLAB
% INTLAB was designed and written by S.M. Rump, head of the Institute for Reliable Computing,
% Hamburg University of Technology. Suggestions are always welcome to rump (at) tuhh.de

repo: OpenMAS


text:
%
% Smooth point-set registration method using neighboring constraints
% -------------------------------------------------------------------
%
% Authors: Gerard Sanromà, René Alquézar and Francesc Serratosa
%
% Contact: gsanorma@gmail.com
% Date: 15/02/2012
%
% Sets arcs among the Eave*n/2 closest pairs of points, where 'n' is the
% number of nodes and 'Eave' is the average number of arcs that we want at
% each node
%
% Input:
%   g: coordinates of the graph nodes
%   Eave: average number of arcs that we want at each node
%
% Output:
%   a1,a2: edges so that i-th edge joins a1(i)-th and a2(i)-th nodes
%

function [a1,a2] = adjZD_eff(g,Eave)

l = length(g);
D = zeros(l);
for it_dim = 1:2
    D = D +(g(:,it_dim)*ones(1,l) - ones(l,1)*g(:,it_dim)').^2;
end
D = sqrt(D);
[v,I] = sort(D(:),'ascend');
[N1,N2] = ind2sub(size(D),I);
indx = find(N1 > N2);
%
[a1,a2] = ind2sub([l,l],I(indx(1:floor(Eave*l/2))));


text:
function [] = lexmkpackage(fileSpec,varargin)
%
% lexmkpackage.m--Makes a software package by lexically checking the
% dependencies of the specified m-file(s). Uses the lexdepfun.m function,
% rather than Mathwork's depfun.m program, so can be used for assembling
% packages for GUI applications (Matlab's depfun.m does not detect callback
% dependencies).
%
% The package created by lexmkpackage.m is placed in a directory named
% "mkpackageDir" under the current working directory.
%
% The user specifies the m-files whose dependencies are to be searched
% using the input argument "fileSpec". fileSpec can be the fully-qualified
% filename of an m-file. Alternatively, if the m-file is on Matlab's search
% path, then just the name of the m-file (without path information) can be
% provided instead.
%
% Multiple filenames can be specified either as a cell array of strings or
% using the '*' wildcard character.
%
% lexdepfun.m detects only user-defined m-files.
%
% lexdepfun.m cannot detect dependant functions whose m-file names are
% generated dynamically at runtime with, for example, eval.m.
%
% Options:
%
%    LEXMKPACKAGE(fileSpec,...,'-quiet') suppresses progress output.
%
%    LEXMKPACKAGE(fileSpec,...,'-preserve') preserves the directory
%    structure of the dependent m-files when creating the package. This,
%    the default behaviour, is necessary to avoid namespace collisions when
%    multiple m-files with the same name are encountered or when "private"
%    directories contain special, modified versions of standard Matlab
%    routines. An additional m-file, package_path.m is created. Running
%    package_path.m will add the paths of the package sub-directories to
%    Matlab's search path.
%
%    LEXMKPACKAGE(fileSpec,...,'-consolidate') puts all dependent m-files
%    into the same directory. This method gives more convenient results
%    than the '-preserve' method, but will fail if multiple m-files with
%    the same name are encountered.
%
%    LEXMKPACKAGE(fileSpec,...,'-file') puts the code from all dependent
%    m-files into a single m-file. This method will only work if the
%    fileSpec input argument refers to a single "root" m-file. Note that
%    this method should not be used for GUI packages, as GUI callback
%    routines must be visible from the root workspace. The '-file' method
%    will not produce a useable file when any of the dependent m-files are
%    scripts, rather than functions.
%
% N.B., the package created by lexmkpackage.m will frequently contain extra
% files in addition to the ones actually needed. This is an unavoidable
% side-effect of the lexical dependency search algorithm.
%
% Syntax: lexmkpackage(fileSpec,<'-file'>,<'-preserve'>,<'-consolidate'>,<'-quiet'>,)
%
% e.g.,   lexmkpackage('/home/bartlett/matlab/units/*.m')
% e.g.,   lexmkpackage('date2yearday','-file','-quiet')
% e.g.,   lexmkpackage('acquire_vmp')
% e.g.,   lexmkpackage({'lsm' 'newfig'},'-consolidate')
% e.g.,   lexmkpackage({'/home/bartlett/matlab/FileOps/lsm.m' '/home/bartlett/matlab/plotting/newfig.m'})

% Developed in Matlab 7.3.0.298 (R2006b) on GLNX86.
% Kevin Bartlett (kpb@uvic.ca), 2007-04-27
%-------------------------------------------------------------------------

% Handle input arguments.
makeSingleFileSpecified = 0;
preserveSpecified = 0;
consolidateSpecified = 0;
isQuiet = 0;

if nargin > 1

    for arg = varargin

        arg = lower(arg);

        if strcmp(arg,'-file') == 1
            %doMakeSingleFile = 1;
            makeSingleFileSpecified = 1;
        elseif strcmp(arg,'-preserve') == 1
            preserveSpecified = 1;
        elseif strcmp(arg,'-consolidate') == 1
            consolidateSpecified = 1;
        elseif strcmp(arg,'-quiet') == 1
            isQuiet = 1;
        else
            error([mfilename '.m--Unrecognised input argument ''' varargin{1}  '''.'])
        end % if

    end % for

end % if

% Determine whether directory structure is to be preserved or whether
% dependent mfiles will be consolidated in a single directory or m-file.
% Preservation, consolidation and single-file output are mutually
% exclusive.
doPreserve = 1;
doMakeSingleFile = 0;
doConsolidate = 0;

if consolidateSpecified == 1
    doPreserve = 0;
    doMakeSingleFile = 0;
    doConsolidate = 1;
end % if

if makeSingleFileSpecified == 1
    doPreserve = 0;
    doMakeSingleFile = 1;
    doConsolidate = 0;
end % if

% Make the new directory to contain the package.
mkpackageDirName = 'mkpackageDir';
mkpackageDir = fullfile(pwd,mkpackageDirName);

if exist(mkpackageDir,'dir')==7

    %disp([mfilename '.m--Directory ' mkpackageDir ' already exists.']);
    disp([mfilename '.m--Directory already exists:']);
    disp(mkpackageDir);
    resp = input('Delete it and start over? Y/N:   ','s');

    if ~isempty(strmatch(lower(resp),'yes'))

        [rmdir_success,rmdir_mssg,rmdir_mssgID] = rmdir(mkpackageDir,'s');

        if rmdir_success == 0
            error([mfilename '.m--Failed to remove existing directory: ' mkpackageDir '.']);
        end % if

    else
        disp([mfilename '.m--Execution halted; package not created.']);
        return;
    end % if

end % if

% List the files specified by the user.

% % ...If fileSpec was entered as a cell array, assume it consists of m-file
% % names.
% if iscell(fileSpec)
%     rootFileList = fileSpec;
% else
%     % Otherwise, treat fileSpec as a string indicating the location of an
%     % m-file or files (either a filename or a wildcard-containing file
%     % specifier).
%     rootFileList = list_root_files(fileSpec);
% end % if
%

% Convert fileSpec to a cell array if necessary.
if ~iscell(fileSpec)
    fileSpec = {fileSpec};
end % if

% Assemble list of "root" files to search.
rootFileList = {};

for iFile = 1:length(fileSpec)

    thisFileSpec = fileSpec{iFile};

    % Treat this file specifier string as a sort of regular expression, and
    % look for files that match it.
    thisRootFileList = list_root_files(thisFileSpec);

    % If no file was found to match the file specifier as a regular
    % expression, assume that the user has entered the basename of an
    % m-file and look for it on the Matlab search path.
    if isempty(thisRootFileList)
        thisRootFileList = {which(thisFileSpec)};
    end % if

    % Append any mfiles found with this file specifier to the list of root
    % files.
    if ~isempty(thisRootFileList)
        %rootFileList{end+1} = thisRootFileList;
        [rootFileList{length(rootFileList)+1:length(rootFileList)+length(thisRootFileList)}] = deal(thisRootFileList{:});
    end % if

end % for

if isempty(rootFileList)
    %error([mfilename '.m--Must specify m-file(s) on Matlab''s search path.'])
    error([mfilename '.m--fileSpec argument matched no m-file or any function on Matlab''s search path. Package not created.'])
end % if

% Creating a package inside an m-file is only allowed when a single m-file
% has been specified by the user. This is because if multiple m-files are
% consolidated into a single m-file, only one of the functions is
% accessible from the root workspace, making it pointless to include
% multiple files unless they are dependents of the main function.
if doMakeSingleFile == 1 & length(rootFileList)>1
    error([mfilename '.m--Cannot create a single-file package for multiple m-files.'])
end % if

numFiles = length(rootFileList);

% Find the dependencies for the specified m-files.
depList = {};

for iFile = 1:numFiles

    if ~isQuiet
        disp([mfilename '.m--Finding dependencies of m-file ' num2str(iFile) ' of ' num2str(numFiles)]);
    end % if

    thisFile = rootFileList{iFile};

    [dummy,thisBaseName] = fileparts(thisFile);

    % Find the dependents of this m-file.
    if isQuiet
        thisDepList = lexdepfun(thisFile,'-quiet');
    else
        thisDepList = lexdepfun(thisFile);
    end % if

    % Combine with previous dependencies. It is okay to use union.m here
    % because the dependency list contains fully-qualified mfile names, so
    % the only duplicates that will be removed are TRUE duplicates (not
    % just mfiles with the same basenames).
    depList = union(depList,thisDepList);

end % for

% Test for duplicate mfile basenames. Duplicate names are not necessarily a
% problem if the directory structure is being preserved. If all the m-files
% are going to be consolidated in a single directory or in a single, large
% m-file, on the other hand, duplicate mfilenames will have unpredictable
% results.
[dummy,depBaseNames] = mfileparts(depList);
duplicatesExist = length(unique(depBaseNames)) ~= length((depBaseNames));

if duplicatesExist & (doMakeSingleFile == 1 | doPreserve == 0)
    error([mfilename '.m--Cannot consolidate all dependent m-files in a single directory or single m-file because there are duplicate m-file names.'])
end % if

% Create directory to contain package.

% Mkdir uses only 2 output arguments on linux version of matlab.
%[mkdir_success,mkdir_mssg,mkdir_mssgID] = mkdir(pwd,mkpackageDirName);
[mkdir_success,mkdir_mssg] = mkdir(pwd,mkpackageDirName);

if mkdir_success ~= 1
    error([mfilename '.m--Error creating directory: ' mkdir_mssg]);
end % if

% The list of dependencies contains the name(s) of the root file(s), but
% they may be mixed up with the dependent files. Create a new file list
% containing both the root file(s) and the dependent files, but with the
% root file(s) at the start of the list.
depList = setdiff(depList,rootFileList);
numSourceFiles = length(rootFileList) + length(depList);
[sourceFiles{1:numSourceFiles}] = deal(rootFileList{:},depList{:});
[sourceDirs,sourceBaseNames,sourceExts] = mfileparts(sourceFiles);

% Output the package.
if doMakeSingleFile == 1

    % Output package in the form of a single m-file containing (possibly)
    % multiple functions.

    % There is only one root file allowed when putting all the dependent
    % functions in a single mfile.
    rootFile = rootFileList{1};

    [dummy,rootBaseName,rootExt,dummy] = fileparts(rootFile);
    rootFileName = [rootBaseName rootExt];
    outFile = fullfile(mkpackageDir,rootFileName);
    outFid = fopen(outFile,'at');

    if outFid < 0
        error([mfilename '.m--Failed to open output file for writing.']);
    end % if

    % Output the code from the source files into the output m-file.
    for iSource = 1:length(sourceFiles)

        thisSource = sourceFiles{iSource};

        sourceFid = fopen(thisSource,'rt');

        if thisSource < 0
            error([mfilename '.m--Failed to open mfile for reading.']);
        end % if

        % Specify empty character for whitespace parameter when calling
        % textscan.m, or leading whitespace is discarded.
        %sourceCode = textscan(sourceFid,'%s','delimiter','\n');
        sourceCode = textscan(sourceFid,'%s','delimiter','\n','whitespace','');
        sourceCode = sourceCode{1};

        for iLine = 1:length(sourceCode)
            fprintf(outFid,'%s\n',sourceCode{iLine});
        end % for

        fclose(sourceFid);

        % Print out a divider line to separate functions.
        if iSource < length(sourceFiles)
            fprintf(outFid,'%s\n\n','%-------------------------------------------------------------------------');
        end % if

    end % for each dependent file

    fclose(outFid);

elseif doPreserve == 1

    % Separate the root path from the subdirectory for all the
    % directories to be copied (if all the files are in, say,
    % /home/bartlett/matlab/oceanography/mapping/, can copy them all to
    % a directory named "mapping"; there is no need to preserve the
    % /home/bartlett/matlab/oceanography/ part).

    %for i = 1:144, if strfind(sourceFiles{i},'private'),disp(sourceFiles{i}),end,end
    %[depDirs,depBaseNames,depExts] = mfileparts(sourceFiles);
    [parentPath,sourceSubDirs] = find_parent_path(sourceDirs);

    for iSource = 1:length(sourceFiles)

        thisSource = sourceFiles{iSource};
        thisSubDir = sourceSubDirs{iSource};
        thisBaseName = sourceBaseNames{iSource};
        thisExt = sourceExts{iSource};
        thisTargDir = fullfile(mkpackageDir,thisSubDir);

        if exist(thisTargDir,'dir') == 0

            mkdirSuccess = mkdir(thisTargDir);

            if mkdirSuccess ~= 1
                error([mfilename '.m--Failed to make directory ' thisTargDir '.']);
            end % if

        end % if

        thisTarget = fullfile(thisTargDir,[thisBaseName thisExt]);
        copySuccess = copyfile(thisSource,thisTarget);

        if copySuccess ~= 1
            error([mfilename '.m--Failed to copy file ' thisTarget '.']);
        end % if

    end % for each dependent file

    % Create m-file, package_path.m, to contain code for adding the package
    % contents to Matlab's path.
    pathFileFid = fopen(fullfile(mkpackageDir,'package_path.m'),'wt');

    if pathFileFid < 0
        error([mfilename '.m--Failed to open ''package_path.m'' file for writing.']);
    end % if

    fprintf(pathFileFid,'%s\n','function [] = package_path()');
    fprintf(pathFileFid,'%s\n','%');
    fprintf(pathFileFid,'%s\n','% package_path.m--Adds the package to the Matlab path.');
    fprintf(pathFileFid,'%s\n','%');
    fprintf(pathFileFid,'%s\n',['% Created automatically by ' mfilename '.m, ' datestr(now,31)]);
    fprintf(pathFileFid,'%s\n','%');
    fprintf(pathFileFid,'%s\n','%-------------------------------------------------------------------------');
    fprintf(pathFileFid,'%s\n','');
    fprintf(pathFileFid,'%s\n','% Get absolute path to package_path.m.');
    fprintf(pathFileFid,'%s\n','thisMFile = mfilename(''fullpath'');');
    fprintf(pathFileFid,'%s\n','thisDir = fileparts(thisMFile);');
    fprintf(pathFileFid,'%s\n','');
    fprintf(pathFileFid,'%s\n','% Add the paths of the mfiles in the package to the Matlab path.');

    pathsToAdd = unique(sourceSubDirs);

    for iPath = 1:length(pathsToAdd)

        thisPartialPath = pathsToAdd{iPath};

        % "Private" directories are not allowed in Matlab's search path.
        [dummy,bottomDirName] = fileparts(thisPartialPath);

        if strcmp(bottomDirName,'private') == 0
            thisFullPathStr = ['[thisDir ''' filesep thisPartialPath ''']'];
            addPathStr = ['addpath(' thisFullPathStr ');'];
            fprintf(pathFileFid,'%s\n',addPathStr);
        end % if

    end % for

    fclose(pathFileFid);

elseif doConsolidate == 1

    % Consolidate all mfiles in one directory.
    for iSource = 1:length(sourceFiles)

        thisSource = sourceFiles{iSource};
        thisBaseName = sourceBaseNames{iSource};
        thisExt = sourceExts{iSource};

        thisTarget = fullfile(mkpackageDirName,[thisBaseName thisExt]);
        copySuccess = copyfile(thisSource,thisTarget);

        if copySuccess ~= 1
            error([mfilename '.m--Failed to copy file ' thisTarget '.']);
        end % if

    end % for each dependent file

end % if making single file package

if isQuiet == 0
    disp([mfilename '.m--Package created.']);
end % if

%------------------------------------------------------------------------------
function fileList = list_root_files(fileSpec)
%
% list_root_files.m--Makes a list of the "root" files from the
% file-specifier string provided by the user and returns it as a CELL
% array.
%
% LIST_ROOT_FILES(fileSpec) returns a list of files according to the string
% found in fileSpec. fileSpec may include the wildcard character "*" and/or
% path information.
%
% Syntax: fileList = list_root_files(fileSpec);
%
% e.g., fileList = list_root_files('D:\MATLAB\toolbox\matlab\graph2d\plot*.m')

% Developed in Matlab 6.5.0.180913a (R13) on SUN OS 5.8.
% Kevin Bartlett(kpb@hawaii.edu), 2003/12/18, 11:04
%------------------------------------------------------------------------------

% If path not specified, default to current directory.
[filePath,name,extension,ver] = fileparts(fileSpec);

if isempty(filePath),
    filePath = pwd;
    fileSpec = fullfile(filePath,[name,extension,ver]);
end % if

% If fileSpec represents an existing directory with no file name or wild
% cards appended to it, then return an empty list.
if exist(fileSpec,'dir') == 7
    fileList = {};
else

    % Get list of files conforming to file specifier string.
    eval(['dirStruct = dir(''' fileSpec ''');'])
    fileList = char(dirStruct.name);

    % Exit if no files.
    if isempty(fileList)
        fileList = {};
        return;
    end % if

    % Discard directories from the list.
    isDir = {};
    [isDir{1:length(dirStruct),1}] = deal(dirStruct.isdir);
    isDir = cat(1,isDir{:});
    fileList = fileList(isDir == 0,:);

    % Convert to cell array.
    if isempty(fileList)
        fileList = {};
    else
        fileList = cellstr(fileList);
    end % if

    % Prepend path information to the list.
    numFiles = length(fileList);

    newFileList = cell(numFiles,1);

    for iFile = 1:numFiles
        newFileList{iFile} = fullfile(filePath,fileList{iFile});
    end % for

    fileList = newFileList;
    clear newFileList;

    % Sort the file list to get consistent file order.
    fileList = sort(fileList);

end % if

%-------------------------------------------------------------------------
function [pathStr,baseName,ext,versn] = mfileparts(fileList)
%
% mfileparts.m--Runs Mathworks' fileparts.m function on multiple filenames
% (fileparts.m only works on one file at a time). Result is returned in the
% form of cell arrays.
%
% If input argument is a single filename in the form of a string, rather
% than a cell array, mfileparts.m acts just like fileparts.m and
% returns strings rather than cell arrays.
%
% Syntax: [pathStr,baseName,ext,versn] = mfileparts(fileList)
%
% e.g.,   fileList = {'/home/bartlett/matlab/FileOps/where.m',...
%                     '/home/bartlett/matlab/bearings/count_rots.m',...
%                     '/home/bartlett/matlab/bearings/plot_rots.m',...
%                     '/home/bartlett/matlab/colours/getcolours.m',...
%                     '/home/bartlett/matlab/filters/decimate_me.m'};
%         [pathStr,baseName,ext,versn] = mfileparts(fileList)
%
% e.g.,   [pathStr,baseName,ext,versn] = mfileparts('/home/bartlett/matlab/FileOps/where.m')

% Developed in Matlab 7.3.0.298 (R2006b) on GLNX86.
% Kevin Bartlett (kpb@uvic.ca), 2007-04-27 16:08
%-------------------------------------------------------------------------

if iscell(fileList)
    isCellInput = 1;
elseif isstr(fileList)
    isCellInput = 0;
    fileList = cellstr(fileList);
else
    error([mfilename '.m--Input must be a cell array of strings or a single string.']);
end % if

nFiles = length(fileList);

[pathStr{1:nFiles}] = deal('');
baseName = pathStr;
ext = pathStr;
versn = pathStr;

for iFile = 1:nFiles
    thisFile = fileList{iFile};
    [thisPathStr,thisBaseName,thisExt,thisVersn] = fileparts(thisFile);
    pathStr{iFile} = thisPathStr;
    baseName{iFile} = thisBaseName;
    ext{iFile} = thisExt;
    versn{iFile} = thisVersn;
end % for

% If input was a string, make output just like fileparts.m.
if ~isCellInput
    pathStr = pathStr{1};
    baseName = baseName{1};
    ext = ext{1};
    versn = versn{1};
end % if

%-------------------------------------------------------------------------
function [parentPath,relativePaths] = find_parent_path(fileList)
%
% find_parent_path.m--Given a list of files (or directories),
% find_parent_path.m returns their shared parent path and their individual,
% relative paths.
%
% This program was developed on a Unix machine, but could presumably be
% adapted to Windows fairly easily.
%
% Syntax: [parentPath,relativePaths] = find_parent_path(fileList)
%
% e.g.,   fileList = {'/home/bartlett/matlab/FileOps/where.m',...
%                     '/home/bartlett/matlab/bearings/count_rots.m',...
%                     '/home/bartlett/matlab/bearings/plot_rots.m',...
%                     '/home/bartlett/matlab/colours/getcolours.m',...
%                     '/home/bartlett/matlab/colours/modcolours/redshift.m',...
%                     '/home/bartlett/matlab/filters/decimate_me.m'};
%         [parentPath,relativePaths] = find_parent_path(fileList)
%
% e.g., (Windows):
%        fileList = {'C:\Documents and Settings\bartlett\filesFileOps\where.m',...
%                    'C:\Documents and Settings\bartlett\bearings\count_rots.m',...
%                    'C:\Documents and Settings\bartlett\bearings\plot_rots.m',...
%                    'C:\Documents and Settings\bartlett\colours\getcolours.m',...
%                    'C:\Documents and Settings\bartlett\colours\redshift.m',...
%                    'C:\Documents and Settings\bartlett\filters\decimate_me.m'}
%         [parentPath,relativePaths] = find_parent_path(fileList)

% Developed in Matlab 7.3.0.298 (R2006b) on GLNX86.
% Kevin Bartlett (kpb@uvic.ca), 2007-04-27 16:32
%-------------------------------------------------------------------------

fileSep = filesep;

if strcmp(fileSep,'\')
   textscan_fileSep = '\\';
else
   textscan_fileSep = fileSep;
end % if

parentPathComponents = {};
nFiles = length(fileList);

% If any of the filenames is an empty string, then there is no shared
% parent directory.
for iFile = 1:nFiles

    thisFile = fileList{iFile};

    if isempty(thisFile)
        parentPath = '';
        relativePaths = fileList;
        return;
    end % if

end % for

% Make a cell array, with each cell containing the path components of one
% member of the file list.
%[pathComponents{1:nFiles}] = deal('');
pathComponents = cell(nFiles,1);

for iFile = 1:nFiles
    thisFile = fileList{iFile};
    pathComponents(iFile) = textscan(thisFile,'%s','delimiter',textscan_fileSep);
end % for

% Find the common parent directories.
commonParentExists = 1;

while commonParentExists

    [currTopLevelDir{1:nFiles}] = deal('');

    for iFile = 1:nFiles

        thisFilePathComponents = pathComponents{iFile};

        if length(thisFilePathComponents) > 0
            currTopLevelDir{iFile} = thisFilePathComponents{1};
        end % if

    end % for

    % If all the directories share a leading directory name, strip it
    % off each of the filenames.
    if length(unique(currTopLevelDir)) == 1

        for iFile = 1:nFiles

            thisFilePathComponents = pathComponents{iFile};

            if length(thisFilePathComponents) == 1
                thisFilePathComponents = {};
                pathComponents{iFile} = thisFilePathComponents;
            elseif length(thisFilePathComponents) > 1
                thisFilePathComponents = thisFilePathComponents(2:end);
                pathComponents{iFile} = thisFilePathComponents;
            end % if

            % If there are no remaining path components for this filename,
            % then it cannot share a common parent with the other
            % filenames.
            if isempty(thisFilePathComponents)
                commonParentExists = 0;
            end % if

        end % for

        % If the string giving the current top level is not empty, add it
        % to the components of the parent path (in unix file systems, the
        % first element returned by textscan.m is an empty string).
        if ~isempty(currTopLevelDir{1})
            parentPathComponents{end+1} = currTopLevelDir{1};
        end % if

    else
        commonParentExists = 0;
    end % if

end % while

% Assemble the components of the shared parent path into a single string.
if length(parentPathComponents) == 0
    parentPath = '';
elseif length(parentPathComponents) == 1
    %parentPath = [filesep parentPathComponents{1}];
    parentPath = [parentPathComponents{1}];
else
    %parentPath = [filesep fullfile(parentPathComponents{:})];
    parentPath = [fullfile(parentPathComponents{:})];
end % if

if isunix
    parentPath = [filesep parentPath];
end % if

% Assemble the components of each relative path into a single string.
relativePaths = cell(nFiles,1);

for iFile = 1:nFiles

    thisFilePathComponents = pathComponents{iFile};

    if length(thisFilePathComponents) == 0
        relativePaths{iFile} = '';
    elseif length(thisFilePathComponents) == 1
        relativePaths{iFile} = thisFilePathComponents{:};
    else
        relativePaths{iFile} = fullfile(thisFilePathComponents{:});
    end % if

end % for


text:
classdef FixedInCoordSysSensorSteeringModel < AbstractSensorSteeringModel
    %FixedInCoordSysSensorSteeringModel Summary of this class goes here
    %   Detailed explanation goes here

    properties
        rhtAsc(1,1) double
        dec(1,1) double
        roll(1,1) double
        coordSys AbstractGeometricCoordSystem

        lvdData LvdData
    end

    methods
        function obj = FixedInCoordSysSensorSteeringModel(rhtAsc, dec, roll, coordSys, lvdData)
            obj.rhtAsc = rhtAsc;
            obj.dec = dec;
            obj.roll = roll;
            obj.coordSys = coordSys;

            obj.lvdData = lvdData;
        end

        function [boreDir] = getBoresightVector(obj, time, vehElemSet, ~, inFrame)
            rotMat = obj.coordSys.getCoordSysAtTime(time, vehElemSet, inFrame);

            [x,y,z] = sph2cart(obj.rhtAsc, obj.dec, 1);
            v = [x;y;z];

            boreDir = rotMat * v;
        end

        function rollAngle = getBoresightRollAngle(obj)
            rollAngle = obj.roll;
        end

        function parentDcm = getSensorParentDcmToInertial(obj, time, vehElemSet, dcm, inFrame)
            parentDcm = obj.coordSys.getCoordSysAtTime(time, vehElemSet, inFrame);
        end

        %body to inertial
        function sensorDcm = getSensorDcmToInertial(obj, time, vehElemSet, dcm, inFrame)
%             sensorToParentDcm = eul2rotmARH([obj.rhtAsc,obj.dec,obj.roll],'zyx');
%             parentToInertialDcm = obj.getSensorParentDcmToInertial(time, vehElemSet, dcm, inFrame);
%
%             sensorDcm = parentToInertialDcm * sensorToParentDcm;

            boreDir = obj.getBoresightVector(time, vehElemSet, dcm, inFrame);
            M1 = vrrotvec2mat(vrrotvec([1;0;0], boreDir));
            M2 = rotx(rad2deg(obj.getBoresightRollAngle()));

            sensorDcm = M2 * M1;
        end

        function tf = isVehDependent(obj)
            tf = obj.coordSys.isVehDependent();
        end

        function useTf = openEditDialog(obj)
            output = AppDesignerGUIOutput({false});
            lvd_EditFixedInCoordSysSensorSteeringModelGUI_App(obj, obj.lvdData, output);
            useTf = output.output{1};
        end

        function enum = getEnum(obj)
            enum = SensorSteeringModelEnum.FixedInCoordSys;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% This class is derived from:
%
% An implementation of direct collocation
% Joel Andersson, 2016
% https://github.com/casadi/casadi/blob/master/docs/examples/matlab/direct_collocation.m
%
% CasADi -- A symbolic framework for dynamic optimization.
% Copyright (C) 2010-2014 Joel Andersson, Joris Gillis, Moritz Diehl,
%                         K.U. Leuven. All rights reserved.
% Copyright (C) 2011-2014 Greg Horn
% Under GNU Lesser General Public License

classdef Collocation < handle

  properties

    daefun
    pathcostfun

    coefficients
    coeff_eval
    coeff_der
    coeff_int

    vars
    num_x
    num_z
    num_u
    num_p
    num_t

    num_i

    tau_root
    order
  end

  methods

    function self = Collocation(states, algvars, controls, parameters, statesOrder, daefun, pathcostfun, d)

      nx = length(states);
      nz = length(algvars);
      nu = length(controls);
      np = length(parameters);
      nt = d;

      v = ocl.types.Structure();
      v.addRepeated({'states', 'algvars'}, {states, algvars}, d);

      tau = [0 ocl.collocation.collocationPoints(d)];

      coeff = ocl.collocation.coefficients(tau);

      self.daefun = daefun;
      self.pathcostfun = pathcostfun;

      self.coefficients = coeff;
      self.coeff_eval = ocl.collocation.evalCoefficients(coeff, d, 1.0);
      self.coeff_der = ocl.collocation.evalCoefficientsDerivative(coeff, tau, d);
      self.coeff_int = ocl.collocation.evalCoefficientsIntegral(coeff, d, 1.0);

      self.vars = v;
      self.num_x = nx;
      self.num_z = nz;
      self.num_u = nu;
      self.num_p = np;
      self.num_t = nt;

      self.num_i = length(v);

      self.tau_root = tau;
      self.order = d;

    end

  end
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function bbci_log_write(fid, format, varargin)
%BBCI_LOG_WRITE - Print information to screen or file
%
%Synopsis:
%  bbci_log_write(FID, FORMAT, PARAM1, ...)
%  bbci_log_write(DATA, FORMAT, PARAM1, ...)
%
%Arguments:
%  FID    - Vector of file identifier, e.g. [1 file_id] for printing log
%           log information to the screen and to the file of ID file_id.
%  DATA   - DATA struct as in bbci_calibate and bbci_apply. The FID of
%           the log file is in the subfield DATA.log.fid.
%  FORMAT - As for fprintf.
%  PARAMx - Parameters passed to fprintf.

if isstruct(fid),  % First input argument is given as DATA
  fid= fid.log.fid;
end

for k= 1:length(fid),
  fprintf(fid(k), [format '\n'], varargin{:});
end

author: bbci


text:
function X = COPAR_updateX(Y, Y_range, D, X, opts)
% function X = COPAR_updateX(Y, Y_range, D, X, opts)
% updating X in COPAR.
% -----------------------------------------------
% Author: Tiep Vu, thv102@psu.edu, 5/11/2016
%         (http://www.personal.psu.edu/thv102/)
% -----------------------------------------------
    if nargin == 0
        addpath(fullfile('..', 'utils'));
        fprintf('Test most\n');
        d = 3000;
        N = 30;
        k = 20;
        k0 = 50;
        C = 100 ;
        opts.k0 = k0;
        opts.lambda = 0.01;
        opts.eta = 0.1;
        Y = normc(rand(d, C*N));
        D = normc(rand(d, C*k + k0));
        X = zeros(size(D,2), size(Y,2));
        Y_range = N*(0:C);
        opts.D_range = k* (0:C);
        opts.D_range_ext = [opts.D_range opts.D_range(end)+k0];
        opts.verbose = false;
        opts.max_iter = 30;
    end
    %%
    C = numel(Y_range) - 1;
    DtD = D'*D;
    DtY = D'*Y;
    L = 2*max(eig(DtD)) + 10;
    optsX = opts;
    optsX.max_iter = 100;
    %%
    for c = 1: C
        Xc = get_block_col(X, c, Y_range);
        X(:, Y_range(c)+1: Y_range(c+1)) = ...
            COPAR_updateXc(DtD, DtY, Y_range, Xc, c, L, optsX);
        if opts.verbose
            costXc = COPAR_cost(Y, Y_range, D, X, opts);
            fprintf('class = %3d | costXc: %5f\n', c, costXc);
        end
    end
    %%
    if nargin == 0
        X = [];
    end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function cmap = overlapCmap(numGrays,numColors,rgyFlag)
%
% cmap = overlapCmap(numGrays,numColors,rgyFlag)
%
% Makes color map array with numGrays grays,
% and a numColors sub-map that's divided into
% 3 color bands (red, green, yellow, or red,
% blue purple, depending on the setting of rgyFlag).
% Used for visualizing overlap maps, though may be
% useful for other purposes.
%
% ras 2/04
if ieNotDefined('rgyFlag')
    rgyFlag = 1;
end

if ~exist('numGrays','var')
    numGrays=128;
end
if ~exist('numColors','var')
    numColors=128;
end

if rgyFlag==1
    % red/green/yellow colors
    colA = [.9 0 0];
    colB = [0 .9 0];
    colC = [.9 .9 0];
else
    % red/blue/purple colors
    colA = [.9 0 0];
    colB = [0 0 .9];
    colC = [.9 0 .9];
end

rngA = 1:round(numColors/3);
rngB = round(numColors/3)+1:ceil((2/3)*numColors);
rngC = ceil((2/3)*numColors)+1:numColors;

colors = zeros(numColors,3);
colors(rngA,:) = repmat(colA,[length(rngA) 1]);
colors(rngB,:) = repmat(colB,[length(rngB) 1]);
colors(rngC,:) = repmat(colC,[length(rngC) 1]);

cmap = [gray(numGrays); colors];

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [path, j1, j2] = hmmViterbiCM(logpi, logA, logB)
% Find the most-probable (Viterbi) path through the HMM state trellis.
% logpi(j) = log of initial state distribution
% logA(i,j) = log of transition matrix
% logB(k,t) = log of soft evidence
% * we use log of inputs for compatability with .mex version *
% * called hmmViterbiC since Matlab has an hmmViterbi function already *
%%

% This file is from pmtk3.googlecode.com

pi = exp(logpi);
A = exp(logA);
B = exp(logB);

[K T] = size(B);
delta = zeros(K,T);
psi = zeros(K,T);
path = zeros(1,T);
t=1;
delta(:,t) = normalize(pi(:) .* B(:,t));
psi(:,t) = 0; % arbitrary value, since there is no predecessor to t=1
for t=2:T
    for j=1:K
        [delta(j,t), psi(j,t)] = max(delta(:,t-1) .* A(:,j));
        delta(j,t) = delta(j,t) * B(j,t);
    end
    delta(:,t) = normalize(delta(:,t));
end

% Traceback
[p, path(T)] = max(delta(:,T));
for t=T-1:-1:1
    path(t) = psi(path(t+1),t+1);
end

j1 = []; % for .mex compatability
j2 = [];
end

author: marianux

repo: ecg-kit

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