text:
function submitWeb(partId)
%SUBMITWEB Generates a base64 encoded string for web-based submissions
%   SUBMITWEB() will generate a base64 encoded string so that you can submit your
%   solutions via a web form

  fprintf('==\n== [ml-class] Submitting Solutions | Programming Exercise %s\n==\n', ...
          homework_id());
  if ~exist('partId', 'var') || isempty(partId)
    partId = promptPart();
  end

  % Check valid partId
  partNames = validParts();
  if ~isValidPartId(partId)
    fprintf('!! Invalid homework part selected.\n');
    fprintf('!! Expected an integer from 1 to %d.\n', numel(partNames));
    fprintf('!! Submission Cancelled\n');
    return
  end

  [login] = loginPrompt();
  if isempty(login)
    fprintf('!! Submission Cancelled\n');
    return
  end

  [result] = submitSolution(login, partId, output(partId), ...
                            source(partId));
  result = base64encode(result);

  fprintf('\nSave as submission file [submit_ex%s_part%d.txt]: ', ...
          homework_id(), partId);
  saveAsFile = input('', 's');
  if (isempty(saveAsFile))
    saveAsFile = sprintf('submit_ex%s_part%d.txt', homework_id(), partId);
  end

  fid = fopen(saveAsFile, 'w');
  if (fid)
    fwrite(fid, result);
    fclose(fid);
    fprintf('\nSaved your solutions to %s.\n\n', saveAsFile);
    fprintf(['You can now submit your solutions through the web \n' ...
             'form in the programming exercises. Select the corresponding \n' ...
             'programming exercise to access the form.\n']);

  else
    fprintf('Unable to save to %s\n\n', saveAsFile);
    fprintf(['You can create a submission file by saving the \n' ...
             'following text in a file: (press enter to continue)\n\n']);
    pause;
    fprintf(result);
  end

end

% ================== CONFIGURABLES FOR EACH HOMEWORK ==================

function id = homework_id()
  id = '8';
end

function [partNames] = validParts()
  partNames = { 'Estimate Gaussian Parameters', ...
                'Select Threshold' ...
                'Collaborative Filtering Cost', ...
                'Collaborative Filtering Gradient', ...
                'Regularized Cost', ...
                'Regularized Gradient' ...
                };
end

function srcs = sources()
  % Separated by part
  srcs = { { 'estimateGaussian.m' }, ...
           { 'selectThreshold.m' }, ...
           { 'cofiCostFunc.m' }, ...
           { 'cofiCostFunc.m' }, ...
           { 'cofiCostFunc.m' }, ...
           { 'cofiCostFunc.m' }, ...
           };
end

function out = output(partId)
  % Random Test Cases
  n_u = 3; n_m = 4; n = 5;
  X = reshape(sin(1:n_m*n), n_m, n);
  Theta = reshape(cos(1:n_u*n), n_u, n);
  Y = reshape(sin(1:2:2*n_m*n_u), n_m, n_u);
  R = Y > 0.5;
  pval = [abs(Y(:)) ; 0.001; 1];
  yval = [R(:) ; 1; 0];
  params = [X(:); Theta(:)];
  if partId == 1
    [mu sigma2] = estimateGaussian(X);
    out = sprintf('%0.5f ', [mu(:); sigma2(:)]);
  elseif partId == 2
    [bestEpsilon bestF1] = selectThreshold(yval, pval);
    out = sprintf('%0.5f ', [bestEpsilon(:); bestF1(:)]);
  elseif partId == 3
    [J] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                       n, 0);
    out = sprintf('%0.5f ', J(:));
  elseif partId == 4
    [J, grad] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                             n, 0);
    out = sprintf('%0.5f ', grad(:));
  elseif partId == 5
    [J] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                       n, 1.5);
    out = sprintf('%0.5f ', J(:));
  elseif partId == 6
    [J, grad] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                             n, 1.5);
    out = sprintf('%0.5f ', grad(:));
  end
end

% ========================= SUBMIT HELPERS =========================

function src = source(partId)
  src = '';
  src_files = sources();
  if partId <= numel(src_files)
      flist = src_files{partId};
      for i = 1:numel(flist)
          fid = fopen(flist{i});
          while ~feof(fid)
            line = fgets(fid);
            src = [src line];
          end
          fclose(fid);
          src = [src '||||||||'];
      end
  end
end

function ret = isValidPartId(partId)
  partNames = validParts();
  ret = (~isempty(partId)) && (partId >= 1) && (partId <= numel(partNames));
end

function partId = promptPart()
  fprintf('== Select which part(s) to submit:\n', ...
          homework_id());
  partNames = validParts();
  srcFiles = sources();
  for i = 1:numel(partNames)
    fprintf('==   %d) %s [', i, partNames{i});
    fprintf(' %s ', srcFiles{i}{:});
    fprintf(']\n');
  end
  fprintf('\nEnter your choice [1-%d]: ', ...
          numel(partNames));
  selPart = input('', 's');
  partId = str2num(selPart);
  if ~isValidPartId(partId)
    partId = -1;
  end
end

function [result, str] = submitSolution(email, part, output, source)

  result = ['a:5:{' ...
            p_s('homework') p_s64(homework_id()) ...
            p_s('part') p_s64(part) ...
            p_s('email') p_s64(email) ...
            p_s('output') p_s64(output) ...
            p_s('source') p_s64(source) ...
            '}'];

end

function s = p_s(str)
   s = ['s:' num2str(numel(str)) ':"' str '";'];
end

function s = p_s64(str)
   str = base64encode(str, '');
   s = ['s:' num2str(numel(str)) ':"' str '";'];
end

% =========================== LOGIN HELPERS ===========================

function [login] = loginPrompt()
  % Prompt for password
  [login] = basicPrompt();
end

function [login] = basicPrompt()
  login = input('Login (Email address): ', 's');
end

% =========================== Base64 Encoder ============================
% Thanks to Peter John Acklam
%

function y = base64encode(x, eol)
%BASE64ENCODE Perform base64 encoding on a string.
%
%   BASE64ENCODE(STR, EOL) encode the given string STR.  EOL is the line ending
%   sequence to use; it is optional and defaults to '\n' (ASCII decimal 10).
%   The returned encoded string is broken into lines of no more than 76
%   characters each, and each line will end with EOL unless it is empty.  Let
%   EOL be empty if you do not want the encoded string broken into lines.
%
%   STR and EOL don't have to be strings (i.e., char arrays).  The only
%   requirement is that they are vectors containing values in the range 0-255.
%
%   This function may be used to encode strings into the Base64 encoding
%   specified in RFC 2045 - MIME (Multipurpose Internet Mail Extensions).  The
%   Base64 encoding is designed to represent arbitrary sequences of octets in a
%   form that need not be humanly readable.  A 65-character subset
%   ([A-Za-z0-9+/=]) of US-ASCII is used, enabling 6 bits to be represented per
%   printable character.
%
%   Examples
%   --------
%
%   If you want to encode a large file, you should encode it in chunks that are
%   a multiple of 57 bytes.  This ensures that the base64 lines line up and
%   that you do not end up with padding in the middle.  57 bytes of data fills
%   one complete base64 line (76 == 57*4/3):
%
%   If ifid and ofid are two file identifiers opened for reading and writing,
%   respectively, then you can base64 encode the data with
%
%      while ~feof(ifid)
%         fwrite(ofid, base64encode(fread(ifid, 60*57)));
%      end
%
%   or, if you have enough memory,
%
%      fwrite(ofid, base64encode(fread(ifid)));
%
%   See also BASE64DECODE.

%   Author:      Peter John Acklam
%   Time-stamp:  2004-02-03 21:36:56 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   if isnumeric(x)
      x = num2str(x);
   end

   % make sure we have the EOL value
   if nargin < 2
      eol = sprintf('\n');
   else
      if sum(size(eol) > 1) > 1
         error('EOL must be a vector.');
      end
      if any(eol(:) > 255)
         error('EOL can not contain values larger than 255.');
      end
   end

   if sum(size(x) > 1) > 1
      error('STR must be a vector.');
   end

   x   = uint8(x);
   eol = uint8(eol);

   ndbytes = length(x);                 % number of decoded bytes
   nchunks = ceil(ndbytes / 3);         % number of chunks/groups
   nebytes = 4 * nchunks;               % number of encoded bytes

   % add padding if necessary, to make the length of x a multiple of 3
   if rem(ndbytes, 3)
      x(end+1 : 3*nchunks) = 0;
   end

   x = reshape(x, [3, nchunks]);        % reshape the data
   y = repmat(uint8(0), 4, nchunks);    % for the encoded data

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Split up every 3 bytes into 4 pieces
   %
   %    aaaaaabb bbbbcccc ccdddddd
   %
   % to form
   %
   %    00aaaaaa 00bbbbbb 00cccccc 00dddddd
   %
   y(1,:) = bitshift(x(1,:), -2);                  % 6 highest bits of x(1,:)

   y(2,:) = bitshift(bitand(x(1,:), 3), 4);        % 2 lowest bits of x(1,:)
   y(2,:) = bitor(y(2,:), bitshift(x(2,:), -4));   % 4 highest bits of x(2,:)

   y(3,:) = bitshift(bitand(x(2,:), 15), 2);       % 4 lowest bits of x(2,:)
   y(3,:) = bitor(y(3,:), bitshift(x(3,:), -6));   % 2 highest bits of x(3,:)

   y(4,:) = bitand(x(3,:), 63);                    % 6 lowest bits of x(3,:)

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Now perform the following mapping
   %
   %   0  - 25  ->  A-Z
   %   26 - 51  ->  a-z
   %   52 - 61  ->  0-9
   %   62       ->  +
   %   63       ->  /
   %
   % We could use a mapping vector like
   %
   %   ['A':'Z', 'a':'z', '0':'9', '+/']
   %
   % but that would require an index vector of class double.
   %
   z = repmat(uint8(0), size(y));
   i =           y <= 25;  z(i) = 'A'      + double(y(i));
   i = 26 <= y & y <= 51;  z(i) = 'a' - 26 + double(y(i));
   i = 52 <= y & y <= 61;  z(i) = '0' - 52 + double(y(i));
   i =           y == 62;  z(i) = '+';
   i =           y == 63;  z(i) = '/';
   y = z;

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Add padding if necessary.
   %
   npbytes = 3 * nchunks - ndbytes;     % number of padding bytes
   if npbytes
      y(end-npbytes+1 : end) = '=';     % '=' is used for padding
   end

   if isempty(eol)

      % reshape to a row vector
      y = reshape(y, [1, nebytes]);

   else

      nlines = ceil(nebytes / 76);      % number of lines
      neolbytes = length(eol);          % number of bytes in eol string

      % pad data so it becomes a multiple of 76 elements
      y = [y(:) ; zeros(76 * nlines - numel(y), 1)];
      y(nebytes + 1 : 76 * nlines) = 0;
      y = reshape(y, 76, nlines);

      % insert eol strings
      eol = eol(:);
      y(end + 1 : end + neolbytes, :) = eol(:, ones(1, nlines));

      % remove padding, but keep the last eol string
      m = nebytes + neolbytes * (nlines - 1);
      n = (76+neolbytes)*nlines - neolbytes;
      y(m+1 : n) = '';

      % extract and reshape to row vector
      y = reshape(y, 1, m+neolbytes);

   end

   % output is a character array
   y = char(y);

end


text:
function [ y, j ] = j_borrow_english ( y, j )

%*****************************************************************************80
%
%% J_BORROW_ENGLISH borrows year-days from years in an English date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, J, a YJ date.
%
  while ( j <= 0 )

    y = y - 1;

    days = year_length_english ( y );

    j = j + days;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [month, day, year] = gdate (jdate)

% convert Julian date to Gregorian (calendar) date

% input

%  jdate = julian day

% output

%  month = calendar month [1 - 12]
%  day   = calendar day [1 - 31]
%  year  = calendar year [yyyy]

%  note: day may include fractional part

% Celestial Computing with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

jd = jdate;

z = fix(jd + .5);
fday = jd + .5 - z;

if (fday < 0)
   fday = fday + 1;
   z = z - 1;
end

if (z < 2299161)
   a = z;
else
   alpha = floor((z - 1867216.25) / 36524.25);
   a = z + 1 + alpha - floor(alpha / 4);
end

b = a + 1524;
c = fix((b - 122.1) / 365.25);
d = fix(365.25 * c);
e = fix((b - d) / 30.6001);
day = b - d - fix(30.6001 * e) + fday;

if (e < 14)
   month = e - 1;
else
   month = e - 13;
end

if (month > 2)
   year = c - 4716;
else
   year = c - 4715;
end

author: Sable

repo: mcbench-benchmarks


text:
function [ i, j ] = r8mat_min_index ( m, n, a )

%*****************************************************************************80
%
%% R8MAT_MIN_INDEX returns the location of the minimum entry of an R8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, real A(M,N), the M by N matrix.
%
%    Output, integer I, J, the indices of the minimum entry of A.
%
  i = -1;
  j = -1;

  for jj = 1 : n
    for ii = 1 : m
      if ( ii == 1 && jj == 1 )
        i = ii;
        j = jj;
      elseif ( a(ii,jj) < a(i,j) )
        i = ii;
        j = jj;
      end
    end
  end

  return
end


text:
function [matchIndsReal, matchIndsEst, meanDist] = computeCentroidMatches(realCellParams, estCellParams, varargin)
    % Matches real and estimates cells based on centroid location.
    % Lacey Kitsch & Biafra Ahanonu
    % started: updating 2016.04.24
    % inputs
        %
    % outputs
        %

    [matchIndsReal, matchIndsEst, meanDist] = ciapkg.image.computeCentroidMatches(realCellParams, estCellParams,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [ftlb] = eV2ftlb(eV)
% Convert energy or work from electron volts to foot-pounds.
% Chad A. Greene 2012
ftlb = eV*1.181705375e-19;


text:
% OtTableBof Check whether the first record was read.
%
% Usage:
%    bof = OtTableBof(htable)

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-04-11 Created (Johan Meijdam, VORtech)

function bof = OtTableBof(htable)
   bof = Ot('table.bof', htable);
end


text:
function f = real(f)
%REAL   Real part of a SINGFUN.
%   REAL(F) is the real part of F.
%
% See also IMAG, ISREAL, CONJ.

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

% Check for empty arguments:
if ( isempty(f) )
    return;
end

% Compute the real part of the smooth part of F.
f.smoothPart = real(f.smoothPart);

% Return a SMOOTHFUN object if F is smooth:
if ( issmooth(f) )
    f = f.smoothPart;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function v=dlyapsq(a,b)
% Solves the discrete Lyapunov equation AV'VA' - V'V +BB' =0
% V is upper triangular with real non-negative diagonal entries
% this is equivalent to v=chol(dlyap(a,b*b')) but better conditioned numerically

%      Copyright (C) Mike Brookes 2002
%      Version: $Id: dlyapsq.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[q,s]=schur(a');
[q,s]=rsf2csf(q,s);
[qd,r]=qr(b'*q,0);
% save r for testing
r0=r;
[m,n]=size(r);
u=zeros(n,n);
if m==1
   for i=1:n-1
      in=i+1:n;
      si=s(i,i);
      aa=sqrt(1-si*si');
      u(i,i)=r(1)/aa;
      u(i,in)=(u(i,i)*si'*s(i,in)+aa*r(2:end))/(eye(n-i)-si'*s(in,in));
      r=aa*(u(i,i)*s(i,in)+u(i,in)*s(in,in))-si*r(2:end);
   end
   u(n,n)=r/sqrt(1-s(n,n)*s(n,n)');

else
   w=zeros(m,1); w(m)=1;
   em=eye(m);
   for i=1:n-m
      in=i+1:n;
      si=s(i,i);
      aa=sqrt(1-si*si');
      u(i,i)=r(1,1)/aa;
      u(i,in)=(u(i,i)*si'*s(i,in)+aa*r(1,2:end))/(eye(n-i)-si'*s(in,in));
      vv=aa*(u(i,i)*s(i,in)+u(i,in)*s(in,in))-si*r(1,2:end);
      rr=zeros(m,n-i);
      rr(1:m-1,:)=r(2:end,2:end);
      [qq,r]=qrupdate(em,rr,w,vv');
   end
   for i=n-m+1:n-1
      in=i+1:n;
      si=s(i,i);
      aa=sqrt(1-si*si');
      u(i,i)=r(1,1)/aa;
      u(i,in)=(u(i,i)*si'*s(i,in)+aa*r(1,2:end))/(eye(n-i)-si'*s(in,in));
      vv=aa*(u(i,i)*s(i,in)+u(i,in)*s(in,in))-si*r(1,2:end);
      rr=zeros(n-i+1,n-i);
      rr(1:n-i,:)=r(2:end,2:end);
      [qq,rr]=qrupdate(eye(n-i+1),rr,w(m-n+i:end),vv');
      r=rr(1:n-i,:);
   end

   u(n,n)=r/sqrt(1-s(n,n)*s(n,n)');

end

v=triu(qr(u*q'));
dv=diag(v);
ix=dv~=0;
v(ix,:)=diag(abs(dv(ix))./dv(ix))*v(ix,:);
if isreal(a) & isreal(b)
   v=real(v);
end

author: decouples


text:
function [J, grad] = costFunction(theta, X, y)
%COSTFUNCTION Compute cost and gradient for logistic regression
%   J = COSTFUNCTION(theta, X, y) computes the cost of using theta as the
%   parameter for logistic regression and the gradient of the cost
%   w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Note: grad should have the same dimensions as theta
%
h_theta = sigmoid(X*theta);
J = (1 / m) * ((-y' * log(h_theta)) - (1 - y)' * log(1 - h_theta));

grad = (1 / m) * (h_theta - y)' * X;

% =============================================================

end

author: atinesh-s

repo: Coursera-Machine-Learning-Stanford


text:
classdef AnnotationMarkerHandle < handle
    properties
        lineHandle;
        textHandle;
    end

    properties (Dependent)
        visible;
    end

    methods
        function set.visible(obj,visibleStr)
            obj.textHandle.Visible = visibleStr;
            obj.lineHandle.Visible = visibleStr;
        end
    end

    methods (Access = public)
        function obj = AnnotationMarkerHandle(lineHandle,textHandle)
            if nargin ~= 0
                obj.lineHandle = lineHandle;
                obj.textHandle = textHandle;
            end
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function S_str = cell2str_v3(S,delimiter)
% S_str = cell2str_v3(S)
% Example: S_str = cell2str_v3(varargin)
if nargin<2, delimiter=', '; end
S = S(:)';
S_str = ''; for k=S, if isnumeric(k{1}), S_str =[S_str ', [' num2str(k{:}) ']']; else S_str =[S_str delimiter sprintf('%s',k{:})]; end; end
%S_str = S_str(3:end);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function geometry_test0493 ( )

%*****************************************************************************80
%
%% TEST0493 tests PARABOLA_EX, PARABOLA_EX2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0493\n' );
  fprintf ( 1, '  PARABOLA_EX finds the extreme value of a parabola\n' );
  fprintf ( 1, '    determined by three points.\n' );
  fprintf ( 1, '  PARABOLA_EX2 finds the extreme value of a parabola\n' );
  fprintf ( 1, '    determined by three points.\n' );

  a =  2.0;
  b = -4.0;
  c = 10.0;

  x1 = 1.0;
  y1 = a * x1 * x1 + b * x1 + c;
  x2 = 2.0;
  y2 = a * x2 * x2 + b * x2 + c;
  x3 = 3.0;
  y3 = a * x3 * x3 + b * x3 + c;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Parabolic coefficients (A,B,C) = %f  %f  %f\n', a, b, c );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  X, Y data\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  X1, Y1 = %f  %f\n', x1, y1 );
  fprintf ( 1, '  X2, Y2 = %f  %f\n', x2, y2 );
  fprintf ( 1, '  X3, Y3 = %f  %f\n', x3, y3 );

  a = 0.0;
  b = 0.0;
  c = 0.0;

  [ xmin, ymin, ierror ] = parabola_ex ( x1, y1, x2, y2, x3, y3 );

  if ( ierror == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  PARABOLA_EX returns (XMIN,YMIN) = %f  %f\n', xmin, ymin );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  PARABOLA_EX returns error code %d\n', ierror );
  end

  [ xmin, ymin, a, b, c, ierror ] = parabola_ex2 ( x1, y1, x2, y2, x3, y3 );

  if ( ierror == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  PARABOLA_EX2 returns (XMIN,YMIN) = %f  %f\n', xmin, ymin );
    fprintf ( 1, '  and (A,B,C) = %f  %f  %f\n', a, b, c );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  PARABOLA_EX2 returns error code %d\n', ierror );
  end

  return
end


text:
function [gamma,phase,dta,alpha,xmesh_lv,ymesh_lv,zmesh_lv]=dicomrt_GAMMAcal2D(evalm,refm,dose_xmesh,dose_ymesh,dose_zmesh,slice,resf,range,dta_criteria,dd_criteria,voi,voiselect,pbopt)
% dicomrt_GAMMAcal2D(evalm,refm,dose_xmesh,dose_ymesh,dose_zmesh,slice,resf,range,dta_criteria,dd_criteria,voi,voiselect,pbopt);
%
% Calculates 2D gamma distribution for a 2D dataset using a 2D algorithm.
%
% evalm is the 2D matrix to evaluate.
% refm is the reference 2D matrix.
%   Both eval and ref can be a TPS generated dataset or a MC generated dataset.
%   Doses are normalised to the target dose whenever is possible.
%   If no normalization dose is provided within the data it is assumed that matrices are already normalised.
% dose_xmesh, dose_ymesh, are the coordinates of the center of the pixels for eval and ref.
% slice is the number of the slice on which gamma should be calculated
% resf is the "resolution factor". Each voxel is divided resf times to allow dose to be interpolated
%   and quantities calculated. The higher resf the more accurate the calculation is, the slower the
%   function will be.
% range is the "search range". Range is the number of pixel about (j,i) that is considered for calculation.
%   If a dta or a dd match is not found within the search range gamma and all the other quantities in that point
%   will be set to infinite value.
% dta_criteria is the Distance-to-agreement criteria in cm (e.g. 0.3).
% dd_criteria:
%   1) IF matrices ARE NOT NORMALIZED before calling the function
%      dta_criteria is the percentage dose difference (e.g. =0.03);
%   2) IF matrices ARE already NORMALIZED before calling the function
%      dd_criteria must be given accordingly to the normalization applyed
%      (e.g. =0.03 for 3% over dose norm =1Gy, or =1.98 for 3% over dose norm=66Gy).
% voi and voiselect are the vois' cell array and the # of the voi to be used for the gamma calculation respectively.
%   They have to be specified together. Both matrices will be masked and reduced in size accordingly with the selected
%   voi's dimensions. This reduce calculation time especially for the 3D algorithm.
% pbopt progress bar option (OPTIONAL): =0 (default) no progress bar is displayed, ~0 progress bar is displayed
%
% Example:
%
% [gamma1_2,phase1_2,dta1_2,alpha1_2,xmesh,ymesh]=dicomrt_GAMMAcal2D(image_eval,image_ref, ...
%                                                             xmesh_red,ymesh_red,40,2,4,0.3,0.03);
%
% returns the calculated gamma function for image_eval vs image_ref at slice # 40 in gamma1_2.
% Gamma is calculated with 3%-3mm DD-DTA criteria.
% The phase is returned in phase1_2, the dta matrix in dta1_2.
% The function return also the direction cosines in alpha1_2. These are the
% direction cosines on the gamma vector in the two dimensional space xy. The direction
% cosines may provide useful information in detecting systematic spatial displacements between eval and ref.
% Coordinates of the xy location where gamma is defined are returned in xmesh and ymesh.
%
% The concept of gamma function was developed by Low et al Med. Phys. 25 656-661.
%
% See also dicomrt_gvhcal, dicomrt_gahcal
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check number of argument and set-up some parameters and variables
error(nargchk(10,13,nargin))

if nargin >10 & (exist('voi')==1 & exist('voiselect')~=1) | (exist('voi')~=1 & exist('voiselect')==1)
    error('dicomrt_GAMMAcal2D: Check VOI and voi2plot: they must be both present or both absent. Exit now!');
end

if exist('pbopt')~=1
    pbopt=0;
end

% Suppress warnings
warning off MATLAB:divideByZero

% Check case and set-up some parameters and variables
[eval_temp,type_dose_one,label1,PatientPosition]=dicomrt_checkinput(evalm);
[ref_temp,type_dose_two,label2,PatientPosition]=dicomrt_checkinput(refm);
local_eval=dicomrt_varfilter(eval_temp);
ref=dicomrt_varfilter(ref_temp);

% Get DICOM-RT toolbox dataset info
local_eval_pointer=eval_temp{1,1};
local_eval_header=local_eval_pointer{1};
ref_pointer=ref_temp{1,1};
ref_header=ref_pointer{1};

[voi_temp]=dicomrt_checkinput(voi);
voi=dicomrt_varfilter(voi_temp);

if strcmpi(type_dose_one,'rtplan')==1 & (strcmpi(type_dose_two,'mc')==1 | strcmpi(type_dose_two,'unknown')==1)
    targetdose_one=getfield(local_eval_header.DoseReferenceSequence,'Item_1','TargetPrescriptionDose'); % normalisation factor
    % Mask arrays
    temp=ref;
    % temp(find(local_eval==0))=0;
    % local_eval(find(temp==0))=0;
    % Normalise dose and calculate dose difference
    eval_norm_temp=local_eval(:,:,slice)/targetdose_one*100;
    ref_norm_temp=temp(:,:,slice)/targetdose_one*100;
elseif strcmpi(type_dose_two,'rtplan')==1 & (strcmpi(type_dose_one,'mc')==1 | strcmpi(type_dose_one,'unknown')==1)
    targetdose_two=getfield(ref{1,1}.DoseReferenceSequence,'Item_1','TargetPrescriptionDose'); % normalisation factor
    % Mask arrays
    temp=local_eval;
    % temp(find(ref==0))=0;
    % ref(find(temp==0))=0;
    % Normalise dose and calculate dose difference
    eval_norm_temp=temp(:,:,slice)/targetdose_two*100;
    ref_norm_temp=ref(:,:,slice)/targetdose_two*100;
elseif strcmpi(type_dose_one,'rtplan')==1 & strcmpi(type_dose_two,'rtplan')==1
    targetdose_one=getfield(local_eval_header.DoseReferenceSequence,'Item_1','TargetPrescriptionDose'); % normalisation factor
    targetdose_two=getfield(ref_header.DoseReferenceSequence,'Item_1','TargetPrescriptionDose'); % normalisation factor
    % no mask is performed. We assume TPS dose are calculated on the same patient using the same patient outline outline
    % Normalise dose and calculate dose difference
    eval_norm_temp=local_eval(:,:,slice)/targetdose_one*100;
    ref_norm_temp=ref(:,:,slice)/targetdose_two*100;
elseif (strcmpi(type_dose_one,'mc')==1 | strcmpi(type_dose_one,'unknown')==1) & ...
        (strcmpi(type_dose_one,'mc')==1 | strcmpi(type_dose_two,'unknown')==1)
    eval_norm_temp=local_eval(:,:,slice);
    ref_norm_temp=ref(:,:,slice);
else
    error('dicomrt_GAMMAcal2D: Cannot determine dose arrays format. Exit now!');
end

% Mask matrices as appropriate
if exist('voi')==1 & exist('voiselect')==1
    [locate_voi_min_x,locate_voi_max_x,locate_voi_min_y,locate_voi_max_y] = dicomrt_voiboundaries(...
        dose_xmesh,dose_ymesh,dose_zmesh,voi_temp,voiselect,PatientPosition);
    eval_norm=eval_norm_temp(locate_voi_min_y:locate_voi_max_y,locate_voi_min_x:locate_voi_max_x);
    ref_norm=ref_norm_temp(locate_voi_min_y:locate_voi_max_y,locate_voi_min_x:locate_voi_max_x);
else
    locate_voi_min_x=1;
    locate_voi_max_x=length(dose_xmesh);
    locate_voi_min_y=1;
    locate_voi_max_y=length(dose_ymesh);
    eval_norm=eval_norm_temp;
    ref_norm=ref_norm_temp;
end

%
% Define parameters
% radians to degrees conversion factor
r2d=360/(2*pi);
d2r=2*pi/360;

% Retrieve basic information
pixel_spacing_x=dicomrt_mmdigit(dose_xmesh(2)-dose_xmesh(1),7);
pixel_spacing_y=dicomrt_mmdigit(dose_ymesh(2)-dose_ymesh(1),7);

% 2D algorithm
%
% create grid for matrix interpolation
% this is done using the resolution factor resf
%[xmesh,ymesh]=dicomrt_build2dgrid(dose_xmesh-pixel_spacing_x/resf,dose_ymesh-pixel_spacing_y/resf);

% Start calculating elapsed time
% tic

disp('(+) Initializing ...');

[xmesh,ymesh]=dicomrt_build2dgrid(dose_xmesh(locate_voi_min_x:locate_voi_max_x),...
    dose_ymesh(locate_voi_min_y:locate_voi_max_y));
xmesh_res=imresize(xmesh,resf,'bilinear');
ymesh_res=imresize(ymesh,resf,'bilinear');

% These variables are exported to use with other functions
xmesh_lv=xmesh(1,:);
ymesh_lv=ymesh(:,1);
zmesh_lv=dose_zmesh(slice);

% Normalize quantities
eval_norm=eval_norm./dd_criteria;
ref_norm=ref_norm./dd_criteria;
xmesh=xmesh./dta_criteria;
ymesh=ymesh./dta_criteria;
xmesh_res=xmesh_res./dta_criteria;
ymesh_res=ymesh_res./dta_criteria;

% interpolate reference matrix
eval_norm_interp(:,:)=interp2(xmesh,ymesh,eval_norm(:,:),xmesh_res,ymesh_res);

% Define output size
gamma=zeros(size(ref_norm));
phase=zeros(size(ref_norm));
dta=zeros(size(ref_norm));
alpha=zeros(size(ref_norm));

% logical steps being taken:
%
% 0) So far the reference matrix was modified to comply with the resolution parameter resf input by the user.
%    Two new matrices which incorporate the dose difference criteria (dd_criteria) were also built.
% 1) The search_range parameter will be used now to define an area around each voxel of the new reference matrices.
%    This area is called Search Area (SA) and will be temporary stored onto a matrix called Transit Area (TA).
%    Distance to agreemenet (DTA) will be then searched within TA for each point of the evaluation matrix.
% 2) If one or more DTA matches are found the dose difference will be calculated for the minimum of the DTA matches.
%
% 1a) Calculate voxels that will define the search volume in 3D

disp('(+) Calculating ...');

if pbopt~=0
    h = waitbar(0,'Calculation progress');
    set(h,'Name','dicomrt_GAMMAcal2D: calculates gamma function');
end

for j=1:size(ref_norm,1)       % loop over y
    for i=1:size(ref_norm,2)   % loop over x
        % all voxels in the new matrices are used for the DTA search. A portion of the new matrices will be copied
        % in the temporary matrix called volume. In order to do this we need to calculate the index that refer to
        % the portion of the matrix to copy.
        if i==1
            iSAmax=range*resf+resf;
            iSAmin=1;
        elseif i>1 & i< range +1
            iSAmax=resf*i+range*resf;
            iSAmin=1;
        elseif i>=range +1 & i<size(ref_norm,2)-range
            iSAmax=resf*i+range*resf;
            iSAmin=i*resf-(resf-1)-range*resf;
        elseif i>=size(ref_norm,2)-range
            iSAmax=size(eval_norm_interp,2);
            iSAmin=i*resf-(resf-1)-range*resf;
        end

        if j==1
            jSAmax=range*resf+resf;
            jSAmin=1;
        elseif j>1 & j< range +1
            jSAmax=resf*j+range*resf;
            jSAmin=1;
        elseif j>=range +1 & j<size(ref_norm,1)-range
            jSAmax=resf*j+range*resf;
            jSAmin=j*resf-(resf-1)-range*resf;
        elseif j>=size(ref_norm,1)-range
            jSAmax=size(eval_norm_interp,1);
            jSAmin=j*resf-(resf-1)-range*resf;
        end

        % debug start
        %display(['i is: ',num2str(i),' - j is: ', num2str(j)]);
        %display(['iSAmin is: ',num2str(iSAmin),' - iSAmax is: ', num2str(iSAmax)]);
        %display(['jSAmin is: ',num2str(jSAmin),' - jSAmax is: ', num2str(jSAmax)]);
        %if j==64 & i==22
        %    disp('debug');
        %end

        temp_xmesh_res=xmesh_res(jSAmin:jSAmax,iSAmin:iSAmax);
        temp_ymesh_res=ymesh_res(jSAmin:jSAmax,iSAmin:iSAmax);
        temp_eval_norm=eval_norm_interp(jSAmin:jSAmax,iSAmin:iSAmax);

        temp_dta=zeros(size(temp_eval_norm));
        %temp_dta(:,:)=nan;
        temp_dta(:,:)=inf;

        temp_dd=zeros(size(temp_eval_norm));
        %temp_dd(:,:)=nan;
        temp_dd(:,:)=inf;

        temp_gamma=zeros(size(temp_eval_norm));
        %temp_gamma(:,:)=nan;
        temp_gamma(:,:)=inf;

        dd_spot=eval_norm(j,i)-ref_norm(j,i);

        [lo]=find(temp_eval_norm<=ref_norm(j,i)+1 & temp_eval_norm>ref_norm(j,i)-1);

        % initialize variable: lo will contain the location (number) of the pixel
        % where the match was found. Pixel numbering is done following the example
        % below for a 2d matrix:
        %
        % A=
        %
        % 1 1 1   1|  /|  /|
        % 2 2 2    | / | / |
        % 1 3 1    |/  |/  |9
        %
        % 6 is the number of the pixel which contains the value 3 and it is the
        % result of the following command:
        %
        % find(A(:,:)==3)

        if isempty(lo)~=1
            temp_dta(lo)=sqrt((temp_xmesh_res(lo)-xmesh(j,i)).^2+(temp_ymesh_res(lo)-ymesh(j,i)).^2);
            temp_dd(lo)=temp_eval_norm(lo)-ref_norm(j,i);
            temp_gamma(lo)=sqrt(temp_dta(lo).^2+temp_dd(lo).^2);
            [temp_gamma_min,temp_gamma_min_index]=min(temp_gamma(lo));
            [gamma(j,i),index]=min([temp_gamma_min sqrt(dd_spot.^2)]);
            % the following if is due because dicomrt_mask, which is used to calculate DVHs GVHs and GAHs,
            % mask matrices to zero. This makes impossible to disguish between a point outside the VOI from
            % a point with dose or gamma equal to 0.
            % In practice this is not a problem for DVHs but can represent a problem for gamma
            % calculation. A value 0 or 0.001 do not change the meaning of gamma.
            %
            if gamma(j,i)==0
                gamma(j,i)=0.01;
            end
            %
            if index==2
                phase(j,i)=0;
                dta(j,i)=0;
                alpha(j,i)=0;
            else
                phase(j,i)=asin(temp_dd(lo(temp_gamma_min_index))./gamma(j,i))*r2d;
                dta(j,i)=temp_dta(lo(temp_gamma_min_index));
                alpha(j,i)=acos(((temp_xmesh_res(lo(temp_gamma_min_index))-xmesh(j,i))./dta(j,i)))*r2d;
                %alpha(j,i)=acos((temp_xmesh_res(lo(temp_gamma_min_index))-xmesh(j,i))/...
                %    (gamma(j,i)*cos(d2r*phase(j,i))))*r2d;
                %alpha(j,i)=acos((temp_xmesh_res(lo(temp_gamma_min_index))-xmesh(j,i))/...
                %    temp_dta(lo(temp_gamma_min_index)))*r2d;
                %alpha(j,i)=acos(temp_dta(lo(temp_gamma_min_index))/(gamma(j,i)*cos(d2r*phase(j,i))))*r2d;
                %alpha(j,i)=acos((temp_xmesh_res(lo(temp_gamma_min_index))-xmesh(j,i))/...
                %    (gamma(j,i)*cos(phase(j,i))))*r2d;
                %alpha(j,i)=asin(temp_dd(lo(temp_gamma_min_index))*sin(phase(j,i)*d2r))*r2d;
            end
        else
            %gamma(j,i)=nan;
            gamma(j,i)=inf;
            %phase(j,i)=nan;
            phase(j,i)=inf;
            %dta(j,i)=nan;
            dta(j,i)=inf;
            %alpha(j,i)=nan;
            alpha(j,i)=inf;
        end
        if pbopt~=0
            waitbar(j/size(ref_norm,1),h);
        end
    end
end
disp('(=) Calculation completed.');
if pbopt~=0
    close(h)
end

% Returns elapsed time
% toc

% Plot
%figure;
%set(gcf,'Name',['dicomrt_GAMMAcal2D: local_eval= ',inputname(1),', ref= ',inputname(2)]);
%surf(gamma,phase,'XData',dose_xmesh(locate_voi_min_x:locate_voi_max_x),'Ydata',dose_ymesh(locate_voi_min_y:locate_voi_max_y));
%colormap jet;
%shading interp;
%title(['gamma/phase',' Z= ',num2str(dose_zmesh(slice))],'FontSize',16);
%xlabel('X axis (cm)','FontSize',12);
%ylabel('Y axis (cm)','FontSize',12);
%grid on;
%colorbar;
%set(gca,'XLim',[min(dose_xmesh(locate_voi_min_x:locate_voi_max_x)) max(dose_xmesh(locate_voi_min_x:locate_voi_max_x))]);
%set(gca,'YLim',[min(dose_ymesh(locate_voi_min_y:locate_voi_max_y)) max(dose_ymesh(locate_voi_min_y:locate_voi_max_y))]);
%set(gca,'ZLim',[min(min(gamma)) max(max(gamma))]);
%set(gca,'ZLim',[min(min(gamma)) 1]);

%figure;
%set(gcf,'Name',['dicomrt_GAMMAcal2D: local_eval= ',inputname(1),', ref= ',inputname(2)]);
%surf(gamma,alpha,'XData',dose_xmesh(locate_voi_min_x:locate_voi_max_x),'Ydata',dose_ymesh(locate_voi_min_y:locate_voi_max_y));
%colormap jet;
%shading interp;
%title(['gamma/alpha',' Z= ',num2str(dose_zmesh(slice))],'FontSize',16);
%xlabel('X axis (cm)','FontSize',12);
%ylabel('Y axis (cm)','FontSize',12);
%grid on;
%colorbar;
%set(gca,'XLim',[min(dose_xmesh(locate_voi_min_x:locate_voi_max_x)) max(dose_xmesh(locate_voi_min_x:locate_voi_max_x))]);
%set(gca,'YLim',[min(dose_ymesh(locate_voi_min_y:locate_voi_max_y)) max(dose_ymesh(locate_voi_min_y:locate_voi_max_y))]);
%set(gca,'ZLim',[min(min(gamma)) max(max(gamma))]);
%set(gca,'ZLim',[min(min(gamma)) 1]);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function d = rcnn_load_cached_pool5_features(cache_name, imdb_name, id)
% d = rcnn_load_cached_pool5_features(cache_name, imdb_name, id)
%   loads cached pool5 features from:
%   feat_cache/[cache_name]/[imdb_name]/[id].mat

% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

file = sprintf('./feat_cache/%s/%s/%s', cache_name, imdb_name, id);

if exist([file '.mat'], 'file')
  d = load(file);
else
  warning('could not load: %s', file);
  d = create_empty();
end

% standardize boxes to double (for overlap calculations, etc.)
d.boxes = double(d.boxes);

% ------------------------------------------------------------------------
function d = create_empty()
% ------------------------------------------------------------------------
d.gt = logical([]);
d.overlap = single([]);
d.boxes = single([]);
d.feat = single([]);
d.class = uint8([]);

author: rbgirshick

repo: rcnn

save_path: github-repos/MATLAB/rbgirshick-rcnn


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

%*****************************************************************************80
%
%% P20_R8_STORE sets or gets the value of a real parameter for problem 20.
%
%  Discussion:
%
%    The only legal input name is 'hold_value', which indicates the
%    value of the variable to be held fixed.  This variable is given
%    the default value of pi/8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION, either 'get' or 'set'.
%
%    Input, string NAME, the name of the variable.
%
%    Input, real VALUE, the value of the variable to be set,
%    if ACTION is 'set'.
%
%    Output, real VALUE, the value of the variable to get,
%    if ACTION is 'get'
%
  persistent hold_value;

  if ( s_eqi ( action, 'get' ) )

    if ( isempty ( hold_value ) )
      hold_value = 0.39269908169872415481;
    end

    if ( s_eqi ( name, 'hold_value' ) )
      value = hold_value;
    else
      value = [];
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P20_R8_STORE - Fatal error!\n' );
      fprintf ( 1, '  Unexpected "get" of NAME = "%s".\n', name );
      error ( 'P20_R8_STORE - Fatal error!' );
    end

  else if ( s_eqi ( action, 'set' ) )

    if ( s_eqi ( name, 'hold_value' ) )
      hold_value = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P20_R8_STORE - Fatal error!\n' );
      fprintf ( 1, '  Unexpected "set" of NAME = "%s".\n', name );
      error ( 'P20_R8_STORE - Fatal error!' );
    end

  else

    value = [];
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P20_R8_STORE - Fatal error!\n' );
    fprintf ( 1, '  Unexpected ACTION = "%s".\n', action );
    error ( 'P20_R8_STORE - Fatal error!' );

  end

  return
end


text:
function y = geo_lpdf(x,p)
%GEO_LPDF     Geometric log probability density function (lpdf).
%   Y = GEO_LPDF(X,P) returns the log of geometric pdf with probability, P,
%   at the values in X.
%
%   The size of Y is the common size of X and P. A scalar input
%   functions as a constant matrix of the same size as the other input.
%
% Copyright (c) 1999-2000 Aki Vehtari

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

y = log(p) + log(1-p) * x;

author: gpstuff-dev

repo: gpstuff

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


text:
%% Copyright (C) 2016, 2019 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

function r = uniop_bool_helper(x, scalar_fcn, opt, varargin)
% SCALAR_FCN can either be the name of a function or a lambda
% definition of a new function.  It can be multiline cell array
% defining a a new function called "sf".
%
% Many SymPy commands return True/False/None (or S.true/S.false/None),
% where None is used for "don't know" or "can't tell".
% The string OPT describes what to do about this.
%   'bool': the default, always return Octave bool (logical) array,
%           maps None to false.
%   'sym': always return symbolic array, which might contain True,
%          False, None.
%   'symIfAnyNone': if the result has any None, return a symbolic
%                   array, otherwise return a bool array.
%
% If you'd like to raise an error on None, do so in SCALAR_FCN.
% See isprime.m for example.
%
% Any additional arguments after OPT are passed to SCALAR_FCN as-is.

  if (iscell(scalar_fcn))
    cmd = scalar_fcn;
  else
    cmd = {['sf = ' scalar_fcn]};
  end

  if ((nargin < 3) || isempty(opt))
    opt = 'bool';
  end

  switch opt
    case 'bool'
      cmd = [ cmd
              'x = _ins[0]'
              'pp = _ins[1:]'
              'if x is not None and x.is_Matrix:'
              '    # bool will map None to False'
              '    return [bool(sf(a, *pp)) for a in x.T],'
              'return bool(sf(x, *pp))' ];

      r = pycall_sympy__ (cmd, x, varargin{:});

      if (~isscalar(x))
        r = reshape(cell2mat(r), size(x));
      end

    case 'sym'
      warning('FIXME: not working for scalars')

      cmd = [ cmd
              'x = _ins[0]'
              'pp = _ins[1:]'
              'if x if not None and x.is_Matrix:'
              '    return x.applyfunc(sf, *pp)'
              'return sf(x, *pp)' ];

      r = pycall_sympy__ (cmd, x, varargin{:});

    case 'symIfAnyNone'
      error('FIXME: not implemented')

    otherwise
      error('uniop_bool_helper: option "%s" not implemented', opt)

  end
end

repo: octsympy


text:
function logLoss = prtEvalLogLoss(classifier,dataSet,nFolds)
% prtEvalLogLoss    Calculate log-loss after applying classifier to dataSet
%
%   logLoss = prtEvalLogLoss(prtClassifier, prtDataSet) returns the
%   log-Loss for prtDataSet when classified by prtClassifier. prtDataSet
%   must be a labeled, prtDataSetStandard object. prtClassifier must be a
%   prtClass object.
%
%   logLoss = prtEvalLogLoss(prtClassifier, prtDataSet, nFolds)  returns
%   the log-Loss of prtDataSet when classified by prtClassifier with K-fold
%   cross-validation. prtDataSet must be a labeled, prtDataSetStandard
%   object. prtClassifier must be a prtClass object. nFolds is the number
%   of folds in the K-fold cross-validation.
%
%   logLoss = prtEvalLogLoss(prtClassifier, prtDataSet, xValInds) same as
%   above, but use crossValidation with specified indices instead of random
%   folds.
%
%       Note: since a lower log-loss is better, if log-loss is used in
%       feature selection, for example, you should optimize over the *negative*
%       log-loss.
%
%   See the help for prtScoreLogLoss for more information.
%
%   Example:
%   dataSet = prtDataGenSpiral;
%   classifier = prtClassDlrt;
%   ll = prtEvalLogLoss(classifier, dataSet)
%
%   See Also: prtScoreLogLoss, prtEvalPdAtPf, prtEvalPfAtPd, prtEvalAuc,
%   prtEvalMinCost, prtEvalPercentCorrect

if nargin < 3 || isempty(nFolds)
    nFolds = 1;
end
results = prtUtilEvalParseAndRun(classifier,dataSet,nFolds);
logLoss = prtScoreLogLoss(results);

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [best,curSel] = findWeakRuleSamples(data,labels,dist,binVals,bins,params)

numDim = size(data,2);

best.dim = 1;
best.error = 0.5;
best.dir = 1;
best.tr = 0;

% curBestErr = 0.5*ones(1,numDim);
% binNo = ones(1,numDim);
% bestDir = ones(1,numDim);

numS = params.numSample;

% KB: only sample if it is helpful
dosample = numS < numel(dist);

% always helpful to normalize
dist = dist / sum(dist);

% fot testing
% for dosample = [false,true],
%
%   niters = 100;
%
%   tic;
%   for iter = 1:niters,

if dosample,

  % KB: copied from randsample
  edges = min([0 cumsum(dist')],1); % protect against accumulated round-off
  edges(end) = 1; % get the upper edge exact
  [~,curSel] = histc(rand(numS,1),edges);

  if any(curSel==0),
    warning('Sanity check: some samples are not being chosen');
    curSel(curSel==0) = size(data,1);
  end

  curBins = bins(:,curSel);
  curLabels = labels(curSel);

else

  curBins = bins;
  curLabels = labels;
  curSel = true(size(labels));
end

% KB: precompute these
numBins = size(binVals,1)+1;
% indices with positive labels
idxpos = curLabels > 0;

% always normalize histograms by Z
if dosample,
  Z = size(curBins,2);
  Zpos = nnz(idxpos);
else
  Z = sum(dist);
  Zpos = sum(dist(idxpos));
end
Zneg = Z - Zpos;

fracpos = Zpos/Z;
fracneg = Zneg/Z;

if dosample,

  posCount = accummatrix(curBins(:,idxpos)',ones(Zpos,1),numBins)'/Z;
  negCount = accummatrix(curBins(:,~idxpos)',ones(Zneg,1),numBins)'/Z;

%   posCount = histc(curBins(:,idxpos),edges,2);
%   posCount = posCount(:,1:end-1) / Z;
%   negCount = histc(curBins(:,~idxpos),edges,2);
%   negCount = negCount(:,1:end-1) / Z;

else

  % weighted histogram, loop-free version

  posCount = accummatrix(curBins(:,idxpos)',dist(idxpos),numBins)';
  negCount = accummatrix(curBins(:,~idxpos)',dist(~idxpos),numBins)';

end

posLeft = cumsum(posCount,2);
posRight = fracpos - posLeft;
negLeft = cumsum(negCount,2);
negRight = fracneg - negLeft;

binErr = posRight+negLeft-posLeft-negRight;
negErr = binErr<0;
binErr(negErr) = -binErr(negErr);
dir = ones(numDim,numBins);
dir(negErr) = -1;
err = 0.5-binErr/2;
[curBestErr,binNo]= min(err(:,1:end-1),[],2);
curBestErr = curBestErr'; binNo = binNo';
bestDir = dir(sub2ind([numDim,numBins],1:numDim,binNo));

[minError minDim] = min(curBestErr); %#ok<UDIM>
best.error = minError;   best.dim = minDim;
best.dir = bestDir(minDim);   best.tr = binVals(binNo(minDim),minDim);

%   end
%
%   fprintf('dosample = %d, time = %f\n',dosample,toc);
%
% end

end

%{
function [predError dir] = getError(data,label,tr, dist)

dir = 1;
predLabel = 2*(data>tr)-1;
predError = sum(( (predLabel.*label) ~=1).*dist);

if(predError>0.5); predError = 1-predError; dir = -1; end
end
%}

author: kristinbranson

repo: JAABA


text:
function [output_signal,output_time] = mvstat(Inputdata_signal,Inputdata_time,windowrange,reductionrange,fn)
%------------ Introduction ----------------------
%This file is created to calculate the mathematical operation like 'max'
%on a pre-defined window range and with a pre-defined reduction range of an vector data.
%This function will need the following inputs in order to calculate the
%wished results:
%   1- Inputdata_signal: A vector containg your Y-Values that need to be
%   filtered.
%   2- Inputdata_time: A vector containg your x-Values that need to be
%   filtered (Normally is a time signal where the sampling rate can be
%   extracted).
%   3- windowrange: in unit of time (s) is the range to calculate the needed
%   mathematical function (The moving range ex. mvmean(y,x,0.02,0.02) or mvrms(y,x,0.02,0.02))
%   4- reductionrange: in unit of time (s) is the scalling reference for
%   the new data
%   5- fn : is used to define the wished mathmatical operation (ex. @max, @min, @mean, @rms, @int)
%The relation between the windowrange input and the reductionrange input
%must be an integer value of 1:1.2:1..10:1. If not the reductionrange would
%be corrected to get one of those relation. this specification give us the
%possibility to calculate the 'mvstat_final' without using a for-Loop which will save a whole amount of time
% for example:
% fn = @mean
% [output_signal,out_put_time] = mvstat_final (y,x,10,10)
% the moving average will be calculated as a result of this function. every
% 10 s the mean value for the last 10 s will be calculated.
% fn = @max
% [output_signal,out_put_time] = mvstat_final (y,x,20,10)
% the moving max will be calculated as a result of this function. every
% 10 s the maximum value for the last 20 s will be calculated.
% fn = @min
% [output_signal,out_put_time] = mvstat_final (y,x,5,10)
% the moving min will be calculated as a result of this function. every
% 10 s the minmum value for the last 5 s will be calculated.
%Notice: when the input values of windowrange and reductionrange are not
%the same then at the beginning of the calculation the fn value of the a
%available data will be calculated.
% Is m.file will also allows you to calculate the intgration of your input
% data in a spesific windowrange. in this case the reductionrange input
% should be defined as in the example:
%fn = @int;
%[Ua_sin,Ua_sin_time] = mvstat(Ua_sin,U1_time,(1/f_grid),delta_t,fn); Ua_sin = Ua_sin*2*f_grid;
% in this example the integration of your input data (in this case U1_sin)
% will be calculated ever one period (normally this mean 0.02 s or 1/50 or
% 1/f_grid in our case). and then the reductionrange is in this case equal
% to the sampling time of the inputdata (delta = U1_time(2)-U1_time(1)).
% this work was done in order to bring FAMOS (IMC) software function into
% matlab.
%Auther: Msc.Eng. Aubai Alkhatib
%Datum: 10.09.2013
%-------- End of introduction ------------------
delta = Inputdata_time(2)-Inputdata_time(1);
length_window = windowrange/delta;
length_reduction = reductionrange/delta;
time = Inputdata_time(end);
if windowrange > reductionrange
    number = windowrange/reductionrange;
    integer = floor(number);
    fract = number-integer;
    if fract >= 0.5
        A = windowrange/(integer+1);
    else
        A = windowrange/(integer);
    end
    reductionrange_new = A;
    integer_A = floor(A);
    fract_A = A-integer_A;
    [~,b] = rat(A);
    if b > 1
        length_reduction_new = round(reductionrange_new/delta);
    else
        length_reduction_new = reductionrange_new/delta;
    end
    i = length_reduction_new;
    ii = 1;
    output = zeros();
    while i < length(Inputdata_signal)
        if i >= length_window
            if fract_A == 0
                switch func2str(fn)
                    case 'mean'
                        output(ii) = mean(Inputdata_signal((i-length_window)+1:i));
                    case 'max'
                        output(ii) = max(Inputdata_signal((i-length_window)+1:i));
                    case 'sum'
                        output(ii) = sum(Inputdata_signal((i-length_window)+1:i));
                    case 'min'
                        output(ii) = min(Inputdata_signal((i-length_window)+1:i));
                    case 'int'
                        output(ii) = sum(Inputdata_signal((i-length_window)+1:i))*reductionrange_new;
                    otherwise
                        output(ii) = rms(Inputdata_signal((i-length_window)+1:i));
                end
            else
                switch func2str(fn)
                    case 'mean'
                        output(ii) = mean(Inputdata_signal((i-length_window):i));
                    case 'max'
                        output(ii) = max(Inputdata_signal((i-length_window):i));
                    case 'sum'
                        output(ii) = sum(Inputdata_signal((i-length_window):i));
                    case 'min'
                        output(ii) = min(Inputdata_signal((i-length_window):i));
                    case 'int'
                        output(ii) = sum(Inputdata_signal((i-length_window)+1:i))*reductionrange_new;
                    otherwise
                        output(ii) = rms(Inputdata_signal((i-length_window):i));
                end
            end
        else
            switch func2str(fn)
                case 'mean'
                    output(ii) = mean(Inputdata_signal(1:i));
                case 'max'
                    output(ii) = max(Inputdata_signal(1:i));
                case 'sum'
                    output(ii) = sum(Inputdata_signal(1:i));
                case 'int'
                    output(ii) = sum(Inputdata_signal(1:i))*reductionrange_new;
                case 'min'
                    output(ii) = min(Inputdata_signal(1:i));
                otherwise
                    output(ii) = rms(Inputdata_signal(1:i));
            end
        end
        i = i + length_reduction_new;
        ii = ii + 1;
    end

    if strcmp(func2str(fn),'int')
        output_signal = output(length_window:length(output));
        output_time = windowrange:reductionrange_new:(length(output)*reductionrange_new);
    elseif isinteger(reductionrange_new)
        output_signal = output;
        output_time = 0:reductionrange_new:(length(output)*reductionrange_new)-reductionrange_new;
    else
        output_signal = output;
        output_time = 0:(length_reduction_new*delta):(length(output)*(length_reduction_new*delta))-(length_reduction_new*delta);
    end
elseif reductionrange > windowrange
    Length_Output = fix(time/windowrange);
    Length_Input_New = Length_Output*windowrange/delta;
    delta_new = length_window;
    delta_time = windowrange;
    Input_reshaped = reshape(Inputdata_signal(1:Length_Input_New),delta_new,Length_Input_New/delta_new);
    switch func2str(fn)
        case 'mean'
            output_signal = mean(Input_reshaped);
        case 'max'
            output_signal = max(Input_reshaped);
        case 'sum'
            output_signal = sum(Input_reshaped);
        case 'int'
            output = sum(Input_reshaped)*delta_time;
        case 'min'
            output_signal = min(Input_reshaped);
        otherwise
            output_signal = rms(Input_reshaped);
    end
    if strcmp(func2str(fn),'int')
        output_signal = output(length_window:length(output));
        output_time = windowrange:delta_time:(Length_Output*delta_time)-delta_time;
    else
        output_time = 0:delta_time:(Length_Output*delta_time)-delta_time;
    end
elseif windowrange == reductionrange
    Length_Output = fix(time/windowrange);
    Length_Input_New = Length_Output*reductionrange/delta;
    delta_new = length_reduction;
    delta_time = windowrange;
    Input_reshaped = reshape(Inputdata_signal(1:Length_Input_New),delta_new,Length_Input_New/delta_new);
    switch func2str(fn)
        case 'mean'
            output_signal = mean(Input_reshaped);
        case 'max'
            output_signal = max(Input_reshaped);
        case 'sum'
            output_signal = sum(Input_reshaped);
        case 'int'
            output = sum(Input_reshaped)*delta_time;
        case 'min'
            output_signal = min(Input_reshaped);
        otherwise
            output_signal = rms(Input_reshaped);
    end
    if strcmp(func2str(fn),'int')
        output_signal = output(length_window:length(output));
        output_time = windowrange:delta_time:(Length_Output*delta_time)-delta_time;
    else
        output_time = 0:delta_time:(Length_Output*delta_time)-delta_time;
    end
end
end


text:
%%
% Test for c-transform in 1D

addpath('../toolbox/');

rep = 'results/c-transform/';
[~,~] = mkdir(rep);

SetAR = @(ar)set(gca, 'PlotBoxAspectRatio', [1 ar 1]);

kx = 10; % #samples
ky = 12;
n = 512/2; % for dense display.

randn('state', 66);
rand('state', 66);

x = rescale( cumsum( rand(kx,1)+.5 ), .05,.95 );
y = rescale( cumsum( rand(ky,1)+.5 ), .05,.95 );
g = linspace(0,1,n);
gx = sort(union(g,x));
gy = sort(union(g,y));

u = rand(kx,1);
u(end) = 1;
u(end/2) = 1;
u(1) = .9;

%%
% Just checking

r = 1;
c = @(x,y)5*abs(x-y).^r;
epsilon = 0;
%
v = ctransform(c, x,u,y, epsilon);
vg = ctransform(c, x,u,gy, epsilon);
%
uu = ctransform(c, y,v,x, epsilon);
uug = ctransform(c, y,v,gx, epsilon);
%
vv = ctransform(c, x,uu,y, epsilon);
vvg = ctransform(c, x,uu,gy, epsilon); % should match vg

rlist = [.5 1 1.5 2];
rlist = [1];

m = .2;

eps_list = linspace(0,.5,7);

for i=1:length(rlist)
    r = rlist(i);
    c = @(x,y)5*abs(x-y).^r;

    clf; hold on;
    for ieps=1:length(eps_list)
        epsilon = eps_list(ieps);
        mc = (ieps-1)/(length(eps_list)-1);
        col = [0 mc 1-mc];
        %
        v = ctransform(c, x,u,y, epsilon);
        vg = ctransform(c, x,u,gy, epsilon);
        %
        plot(gy,vg, '-', 'LineWidth', 2, 'color', col);
    end
    v0 = ctransform(c, x,u,y, 0);
    plot(y,v0, 'b.', 'MarkerSize', 25);
    set(gca, 'XTick', [], 'YTick', []);
    box on; axis tight; axis([0 1 min(v0)-m max(v0)+m]);
    SetAR(1/2);
    saveas(gcf, [rep 'c-transf-v-' num2str(round(10*r)) '.eps'], 'epsc');

    clf; hold on;
    for ieps=1:length(eps_list)
        epsilon = eps_list(ieps);
        mc = (ieps-1)/(length(eps_list)-1);
        col = [1-mc mc 0];
        %
        v = ctransform(c, x,u,y, epsilon);
        vg = ctransform(c, x,u,gy, epsilon);
        %
        uu = ctransform(c, y,v,x, epsilon);
        uug = ctransform(c, y,v,gx, epsilon);
        %
        vv = ctransform(c, x,uu,y, epsilon);
        vvg = ctransform(c, x,uu,gy, epsilon); % should match vg
        %
        plot(gx,uug, '-', 'LineWidth', 2, 'color', col);
    end
    plot(x,u, 'r.', 'MarkerSize', 25);
    set(gca, 'XTick', [], 'YTick', []);
    box on; axis tight; axis([0 1 0-m 1+m]);
    SetAR(1/2);
    saveas(gcf, [rep 'c-transf-u-' num2str(round(10*r)) '.eps'], 'epsc');
end

author: optimaltransport


text:
function [J, flag, new_data] = jacobianFunction(~, x, ~, ~, cj, data)
%	jacobianFunction is used to evaluate the Jacobian Matrix of the P2D model
%	according to the specifications of the IDA numerical solver. Please refer
%	to the IDS user's guide for additional information on this function

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

% Extract the function object previously obtained using CasADi
fJ          = data.fJ;

% Evaluate the Jacobian with respect to the current values of the states
% and their time derivatives.
J           = full(fJ(x,cj));

% Return the values
flag        = 0;
new_data    = [];
end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function [ H_output ] = homography_linearization( H, point, anchor_points )
vega = 5;

% Obtain kernel: Student’s t-weighting
alpha = (1 + pdist2(point,anchor_points)./vega).^(-(vega+1)/2);
alpha = alpha./sum(alpha);

% Linearization using Taylor series
H_output = zeros(3,3);

for i = 1:size(anchor_points,1)
    A = taylor_series( H, anchor_points(i,:));

    H_output = H_output + alpha(i).*A;
end

end

function [ A ] = taylor_series( H, p)
% The first two terms of Taylor series provide the best linearization of H
% H:3*3    p:1*2

h1 = H(1,1); h2 = H(1,2); h3 = H(1,3);
h4 = H(2,1); h5 = H(2,2); h6 = H(2,3);
h7 = H(3,1); h8 = H(3,2); h9 = H(3,3);
x1 = p(1,1);x2 = p(1,2);

dy1dx1 = h1/(h9 + h7*x1 + h8*x2) - (h7*(h3 + h1*x1 + h2*x2))/(h9 + h7*x1 + h8*x2)^2;
dy1dx2 = h2/(h9 + h7*x1 + h8*x2) - (h8*(h3 + h1*x1 + h2*x2))/(h9 + h7*x1 + h8*x2)^2;
dy2dx1 = h4/(h9 + h7*x1 + h8*x2) - (h7*(h6 + h4*x1 + h5*x2))/(h9 + h7*x1 + h8*x2)^2;
dy2dx2 = h5/(h9 + h7*x1 + h8*x2) - (h8*(h6 + h4*x1 + h5*x2))/(h9 + h7*x1 + h8*x2)^2;

y1 = (h3 + h1*x1 + h2*x2)/(h9 + h7*x1 + h8*x2);
y2 = (h6 + h4*x1 + h5*x2)/(h9 + h7*x1 + h8*x2);

A = [dy1dx1 dy1dx2 (y1 - x1*dy1dx1 - x2*dy1dx2);
     dy2dx1 dy2dx2 (y2 - x1*dy2dx1 - x2*dy2dx2);
          0      0                           1];
end

author: YaqiLYU

repo: AANAP

save_path: github-repos/MATLAB/YaqiLYU-AANAP


text:
classdef MimRemapGrey < MimGuiPlugin
    % MimRemapGrey. Gui Plugin for setting paint colour
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimRemapGrey is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Multiple'
        SelectedText = 'Multiple'
        ToolTip = 'Marks airway as supplying multiple lobes'
        Category = 'Airway label'
        Visibility = 'Dataset'
        Mode = 'Edit'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 5
        ButtonHeight = 1
        Icon = 'paint.png'
        IconColour = GemMarkerPoint.DefaultColours{7}
        Location = 37
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.PaintBrushColour = 7;
        end

        function enabled = IsEnabled(gui_app)
             enabled = gui_app.IsDatasetLoaded && gui_app.ImagePanel.OverlayImage.ImageExists && ...
                isequal(gui_app.ImagePanel.SelectedControl, 'Map');
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.PaintBrushColour == 7;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% Mathematics Q561696
% https://math.stackexchange.com/questions/561696
% Solving Non Negative Least Squares by Analogy with Least Squares (MATLAB)
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     05/08/2017
%   *   First release.

%% General Parameters

run('InitScript.m');

figureIdx           = 0; %<! Continue from Question 1
figureCounterSpec   = '%04d';

generateFigures = OFF;

%% Simulation Parameters

numRows = 4;
numCols = 3; %<! Number of Vectors - i (K in the question)

numIterations   = 10000;
stepSize        = 0.075;

%% Generate Data

mA = randn([numRows, numCols]);
vB = randn([numRows, 1]);

%% Solution by CVX

cvx_begin('quiet')
    cvx_precision('best');
    variable vX(numCols)
    minimize( square_pos(  norm(mA * vX - vB, 2) ) );
    subject to
        vX >= 0;
cvx_end

disp([' ']);
disp(['CVX Solution Summary']);
disp(['The CVX Solver Status - ', cvx_status]);
disp(['The Optimal Value Is Given By - ', num2str(cvx_optval)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Solution by Projected Gradient Descent

vX = zeros([numCols, 1]);

for ii = 1:numIterations
    vX = vX - ((stepSize / sqrt(ii)) * mA.' * (mA * vX - vB));
    vX = max(vX, 0);
end

objVal = sum((mA * vX - vB) .^ 2);

disp([' ']);
disp(['Projected Gradient Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(objVal)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function [zPred,PzPred,otherInfo]=sqrtKalmanMeasPred(xPred,SPred,H)
%%SQRTKALMANMEASPRED Perform the measurement prediction part of the
%           measurement update step of the square-root Kalman filter. The
%           function sqrtKalmanUpdateWithPred can be used to complete the
%           measurement update. Separating the measurement prediction step
%           from the rest of the update step can make the creation of
%           multiple measurement association hypotheses from a single
%           target prediction more efficient. The full measurement update
%           function is sqrtKalmanUpdate.
%
%INPUTS: xPred The xDimXnumComp set of numComp predicted target states.
%        SPred The xDimXxDimXnumComp lower-triangular predicted square-root
%              state covariance matrices.
%            H The zDimXxDim measurement matrix. The measurement is modeled
%              as z=H*x+noise and the model is the same for all components
%              given in xPred and SPred.
%
%OUTPUTS: zPred The zDimXnumComp measurement predictions from the filter.
%        PzPred The zDimXzDimXnumComp covariance matrix associated with
%               zPred.
%     otherInfo A structure containing members of intermediate results of
%               this function that can be passed to sqrtKalmanUpdate
%               when updating with a measurement.
%
%The mathematics behind the specific square root implementation used here
%are described in Appendix G of [1].
%
%EXAMPLE:
%With this example, we demonstrate that one gets the same result using
%KalmanUpdate versus calling sqrtKalmanMeasPred and then
%sqrtKalmanUpdateWithPred.
% xPred=[1e3;-2e3;100;200];
% SPred=chol([28,   3.5,    6,  8.5;
%            3.5,    23,  8.5,   11;
%              6,   8.5,   18, 13.5;
%            8.5,    11, 13.5,   13],'lower');
% z=1e3*[-5.498856156296510;
%        1.199241491470584];
% SR=eye(2);
% H=[0, 4, 9, 8;
%    6, 3, 0, 6];
%
% %The update in one step.
% [xUpdate,SUpdate,innov,Szz,W]=sqrtKalmanUpdate(xPred,SPred,z,SR,H);
% %The update in two steps.
% [zPred,PzPred,otherInfo]=sqrtKalmanMeasPred(xPred,SPred,H);
% [xUpdate1,SUpdate1,innov1,Szz1,W1]=sqrtKalmanUpdateWithPred(z,SR,zPred,otherInfo);
% %One will see that the one and two step updates agree.
% max(abs([xUpdate1-xUpdate;SUpdate1(:)-SUpdate(:);innov1(:)-innov;Szz1(:)-Szz(:);W1(:)-W(:)]))
%
%REFERENCES
%[1] David F. Crouse , "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%
%June 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

zDim=size(H,1);
xDim=size(H,2);
numComp=size(xPred,2);

zPred=H*xPred;

PzPred=zeros(zDim,zDim,numComp);
Pxz=zeros(xDim,zDim,numComp);
for k=1:numComp
    temp=H*SPred(:,:,k);
    PzPred(:,:,k)=temp*temp';

    Pxz(:,:,k)=SPred(:,:,k)*SPred(:,:,k)'*H';
    %Pxz is not needed for the measurement prediction, but we compute it
    %here, so that it need not be recomputed again and again if
    %sqrtKalmanUpdateWithPred is called for multiple measurements.
end
otherInfo.xPred=xPred;
otherInfo.SPred=SPred;
otherInfo.Pxz=Pxz;
otherInfo.H=H;

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% Notes for Laplace scheme and dynamic updating of parameters
%==========================================================================
n  = 128;                       % number of samples (time bins)
P  = 4;                         % true parameter
k  = n;                         % precision on fluctuations
pp = 2;                         % prior on parameter
m  = 8;                         % number of data
iS = eye(m,m)*2;                % error precision
s  = sqrtm(inv(iS));

p  = [0;0];                     % initial parameter estimates
for i = 1:n

    x      = 4 + randn(m,1);          % exogenous input
    z      = s*randn(m,1);            % error
    y      = P*x + z;                 % response
    e      = y - p(1)*x;              % prediction error
    Lp     = -e'*iS*x + pp*p(1);
    Lpp    = x'*iS*x  + pp;
    f      = [p(2); (-Lp -k*p(2))];
    dfdx   = [0 1;   -Lpp -k];
    p      = p + spm_dx(dfdx,f,1);
    X(:,i) = p;
    LP(i)  = Lp;
end

% results
%--------------------------------------------------------------------------
subplot(2,1,1)
plot(1:n,X)
title('generalised parameter estimates','FontSize',16)

subplot(2,1,2)
plot(1:n,LP)
title('energy gradient','FontSize',16)

return

% Notes
%==========================================================================
M.f = inline('[x(2); (u - K(1)*x(2))]','x','u','K','M');
M.m = 1;
M.n = 2;
M.l = 2;
M.x = [0;0];
M.u = 0;

N   = 128;
dt  = 1/64;
K   = 4;
[K0,K1,K2] = spm_kernels(M,K,N,dt);

subplot(2,1,1)
plot([1:N]*dt,K1);
xlabel('time (s)')
title('Kernels','FontSize',16)
legend('drive','trace')
axis square

author: spm


text:
function tsss = tbx_cfg_tsss
% Configuration file for toolbox 'TSSS'
%_______________________________________________________________________
% Copyright (C) 2014 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: tbx_cfg_tsss.m 7703 2019-11-22 12:06:29Z guillaume $

tbxdir = fileparts(mfilename('fullpath'));

if ~isdeployed, addpath(tbxdir); end

components = {
    'tsss_config'
    'tsss_config_momentspace'
    };

tsss = cfg_choice;
tsss.tag  = 'tsss';
tsss.name = 'TSSS';
tsss.help = {'Temporal Signal Space Separation (TSSS) toolbox'};

for i = 1:numel(components)
  tsss.values{i} = feval(components{i});
end


text:
function linpack_c_test33 ( )

%*****************************************************************************80
%
%% TEST33 tests CSPFA and CSPDI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 3;

  fprintf ( 1, '\n' )
  fprintf ( 1, 'TEST33\n' );
  fprintf ( 1, '  For a single precision complex (C)\n' );
  fprintf ( 1, '  symmetric matrix in packed storage (SP)\n' );
  fprintf ( 1, '  CSPFA factors the matrix.\n' );
  fprintf ( 1, '  CSPDI computes the determinant or inverse.\n' );
  fprintf ( 1, '\n' )
  fprintf ( 1, '  The matrix order is N = %d\n', n );
%
%  Set the values of the packed matrix A.
%
  k = 0;
  seed = 123456789;

  for j = 1 : n

    for i = 1 : j-1
      k = k + 1;
      [ a(k), seed ] = c4_uniform_01 ( seed );
    end

    k = k + 1;
    [ a(k), seed ] = c4_uniform_01 ( seed );

  end
%
%  Copy the packed matrix into a "normal" matrix.
%
  k = 0;
  for j = 1 : n
    for i = 1 : j
      k = k + 1;
      a_save(i,j) = a(k);
    end
  end

  for j = 1 : n
    a_save(j+1:n,j) = transpose ( a_save(j,j+1:n) );
  end

  fprintf ( 1, '\n' )
  fprintf ( 1, '  The matrix A is\n' );
  fprintf ( 1, '\n' )

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  (%8f  %8f)', real ( a_save(i,j) ), imag ( a_save(i,j) ) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Factor the matrix A.
%
  [ a, ipvt, info ] = cspfa ( a, n );

  if ( info ~= 0 )
    fprintf ( 1, '\n' )
    fprintf ( 1, '  CSPFA returned an error flag INFO = %d\n', info );
    return
  end
%
%  Get the determinant.
%
  job = 10;
  [ a, det ] = cspdi ( a, n, ipvt, job );

  fprintf ( 1, '\n' )
  fprintf ( 1, '  Determinant = (%8f  %8f)*10^(%8f)\n', ...
    real ( det(1) ), imag ( det(1) ), real ( det(2) ) );
%
%  Get the inverse.
%
  job = 01;
  [ a, det ] = cspdi ( a, n, ipvt, job );
%
%  Copy the packed matrix into a "normal" matrix.
%
  k = 0;
  for j = 1 : n
    for i = 1 : j
      k = k + 1;
      b_save(i,j) = a(k);
    end
  end

  for j = 1 : n
    b_save(j+1:n,j) = transpose ( b_save(j,j+1:n) );
  end

  c(1:n,1:n) = b_save(1:n,1:n) * a_save(1:n,1:n);

  fprintf ( 1, '\n' )
  fprintf ( 1, '  The product inv(A) * A is\n' );
  fprintf ( 1, '\n' )

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  (%8f  %8f)', real ( c(i,j) ), imag ( c(i,j) ) );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Local Regression and Likelihood, Figure 1.1.
%
% A linear regression and plot for Spencer's mortality
% data. Unlike the book, I use locfit for this!
%
% Setting 'ev','none' means that locfit computes only
% the parametric fit. In this case, 'deg',1; i.e. linear.

load spencer;
fit = locfit(age,mortality,'deg',1,'ev','none');
figure('Name','fig1_1: locfit linear regression');
lfplot(fit);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function c = tapas_hgf_categorical_norm_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the Hierarchical Gaussian Filter (HGF) for categorical inputs
% restricted to 3 levels, no drift, and no inputs at irregular intervals, in the absence of
% perceptual uncertainty.
%
% This model deals with the situation where an agent has to determine the probability of categorical
% outcomes. The tendencies of these outcomes are modeled as independent Gaussian random processes at
% the second level of the HGF. They are transformed into predictive probabilities at the first level
% by a softmax function (i.e., a logistic sigmoid). This amounts to the assumption that the
% probabilities are performing a Gaussian random walk in logit space. The volatility of all of
% these walks is determined by the same higher-level state x_3 in standard HGF fashion.
%
% The HGF is the model introduced in
%
% Mathys C, Daunizeau J, Friston, KJ, & Stephan KE. (2011). A Bayesian foundation for individual
% learning under uncertainty. Frontiers in Human Neuroscience, 5:39.
%
% and elaborated in
%
% Mathys, C, Lomakina, EI, Daunizeau, J, Iglesias, S, Brodersen, KH, Friston, KJ, & Stephan, KE
% (2014). Uncertainty in perception and the Hierarchical Gaussian Filter. Frontiers in Human
% Neuroscience, 8:825.
%
% This file refers to CATEGORICAL inputs (Eqs 1-3 in Mathys et al., (2011)); for continuous inputs,
% refer to tapas_hgf_config.m, for binary inputs, refer to tapas_hgf_binary.m
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The HGF configuration consists of the priors of parameters and initial values. All priors are
% Gaussian in the space where the quantity they refer to is estimated. They are specified by their
% sufficient statistics: mean and variance (NOT standard deviation).
%
% Quantities are estimated in their native space if they are unbounded (e.g., omega). They are
% estimated in log-space if they have a natural lower bound at zero (e.g., sigma2).
%
% Kappa and theta are estimated in 'logit-space' because bounding them above (in addition to
% their natural lower bound at zero) is an effective means of preventing the exploration of
% parameter regions where the assumptions underlying the variational inversion (cf. Mathys et
% al., 2011) no longer hold.
%
% 'Logit-space' is a logistic sigmoid transformation of native space with a variable upper bound
% a>0:
%
% logit(x) = ln(x/(a-x)); x = a/(1+exp(-logit(x)))
%
% Parameters can be fixed (i.e., set to a fixed value) by setting the variance of their prior to
% zero. Aside from being useful for model comparison, the need for this arises whenever the scale
% and origin of x3 are arbitrary. This is the case if the observation model does not contain the
% representations mu3 and sigma3 from the third level. A choice of scale and origin is then
% implied by fixing the initial value mu3_0 of mu3 and either kappa or omega.
%
% Kappa and theta can be fixed to an arbitrary value by setting the upper bound to twice that
% value and the mean as well as the variance of the prior to zero (this follows immediately from
% the logit transform above).
%
% Fitted trajectories can be plotted by using the command
%
% >> tapas_hgf_categorical_plotTraj(est)
%
% where est is the stucture returned by fitModel. This structure contains the estimated
% perceptual parameters in est.p_prc and the estimated trajectories of the agent's
% representations (cf. Mathys et al., 2011). Their meanings are:
%
%         est.p_prc.mu2_0      initial values of the mu2s
%         est.p_prc.sa2_0      initial values of the sigma2s
%         est.p_prc.mu3_0      initial value of mu3
%         est.p_prc.sa3_0      initial value of sigma3
%         est.p_prc.ka         kappa
%         est.p_prc.om         omega
%         est.p_prc.th         theta
%
%         est.traj.mu          mu
%         est.traj.sa          sigma
%         est.traj.muhat       prediction mean
%         est.traj.sahat       prediction variance
%         est.traj.v           inferred variances of random walks
%         est.traj.w           weighting factor of informational and environmental uncertainty at the 2nd level
%         est.traj.da          prediction errors
%         est.traj.ud          updates with respect to prediction
%         est.traj.psi         precision weights on prediction errors
%         est.traj.epsi        precision-weighted prediction errors
%         est.traj.wt          full weights on prediction errors (at the first level,
%                                  this is the learning rate)
%
% Tips:
% - When analyzing a new dataset, take your inputs u and use
%
%   >> est = tapas_fitModel([], u, 'tapas_hgf_categorical_config', 'tapas_bayes_optimal_categorical_config');
%
%   to determine the Bayes optimal perceptual parameters (given your current priors as defined in
%   this file here, so choose them wide and loose to let the inputs influence the result). You can
%   then use the optimal parameters as your new prior means for the perceptual parameters.
%
% - If you get an error saying that the prior means are in a region where model assumptions are
%   violated, lower the prior means of the omegas, starting with the highest level and proceeding
%   downwards.
%
% - Alternatives are lowering the prior mean of kappa, if they are not fixed, or adjusting
%   the values of the kappas or omegas, if any of them are fixed.
%
% - If the log-model evidence cannot be calculated because the Hessian poses problems, look at
%   est.optim.H and fix the parameters that lead to NaNs.
%
% - Your guide to all these adjustments is the log-model evidence (LME). Whenever the LME increases
%   by at least 3 across datasets, the adjustment was a good idea and can be justified by just this:
%   the LME increased, so you had a better model.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013-2014 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'hgf_categorical';

% Number of states
c.n_outcomes = 3;

% Upper bound for kappa and theta (lower bound is always zero)
c.kaub = 2;
c.thub = 0.1;

% Sufficient statistics of Gaussian parameter priors

% Initial mu2
c.mu2_0mu = repmat(tapas_logit(1/c.n_outcomes,1),1,c.n_outcomes);
c.mu2_0sa = zeros(1,c.n_outcomes);

% Initial sigma2
c.logsa2_0mu = repmat(log(1),1,c.n_outcomes);
c.logsa2_0sa = zeros(1,c.n_outcomes);

% Initial mu3
% Usually best kept fixed to 1 (determines origin on x3-scale).
c.mu3_0mu = 1;
c.mu3_0sa = 0;

% Initial sigma3
c.logsa3_0mu = log(0.1);
c.logsa3_0sa = 1;

% Kappa
% This should be fixed (preferably to 1) if the observation model
% does not use mu3 (kappa then determines the scaling of x3).
c.logitkamu = 0; % If this is 0, and
c.logitkasa = 0; % this is 0, and c.kaub = 2 above, then kappa is fixed to 1

% Omega
c.ommu =  -4;
c.omsa = 5^2;

% Theta
c.logitthmu = 0;
c.logitthsa = 2;

% Gather prior settings in vectors
c.priormus = [
    c.mu2_0mu,...
    c.logsa2_0mu,...
    c.mu3_0mu,...
    c.logsa3_0mu,...
    c.logitkamu,...
    c.ommu,...
    c.logitthmu,...
             ];

c.priorsas = [
    c.mu2_0sa,...
    c.logsa2_0sa,...
    c.mu3_0sa,...
    c.logsa3_0sa,...
    c.logitkasa,...
    c.omsa,...
    c.logitthsa,...
             ];

% Model function handle
c.prc_fun = @tapas_hgf_categorical;

% Handle to function that transforms perceptual parameters to their native space
% from the space they are estimated in
c.transp_prc_fun = @tapas_hgf_categorical_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
function [fname,qname] = usgsdems(latlim,lonlim)

    %USGSDEMS  USGS 1-Degree (3-arc-sec resolution) DEM file names
    %
    % [fname,qname] = USGSDEMS(latlim,lonlim) returns cellarrays of the file
    % names and quadrangle names covering the geographic region for 1-degree
    % USGS digital elevation maps (also referred to as "3-arc second" or
    % "1:250,000 scale" DEMs).  The region is specified by scalar latitude and
    % longitude points, or two element vectors of latitude and longitude
    % limits in units of degrees.
    %
    %
    % The digital elevation map data files are available from the U.S.
    % Geological Survey over the internet from
    % <ftp://edcftp.cr.usgs.gov/pub/data/DEM/250/>.
    %
    % See also: USGSDEM

    %  Copyright 1996-2000 Systems Planning and Analysis, Inc. and The MathWorks, Inc.
    % $Revision: 1399 $
    %  Written by:  A. Kim, W. Stumpf

    if nargin~=2
        error('Incorrect number of arguments')
    end

    if  isequal(size(latlim),[1 1])
        latlim = latlim*[1 1];
    elseif ~isequal(size(latlim),[1 2])
        error('Latitude limit input must be a scalar or 2 element vector')
    end

    if isequal(sort(size(lonlim)),[1 1])
        lonlim = lonlim*[1 1];
    elseif ~isequal(sort(size(lonlim)),[1 2])
        error('Longitude limit input must be a scalar or 2 element vector')
    end

    fid = fopen('usgsdems.dat','r');
    if fid==-1
        error('Couldn''t open usgsdems.dat')
    end

    % preallocate bounding rectangle data for speed

    YMIN = zeros(1,924); YMAX = YMIN;
    XMIN = YMIN; XMAX = YMIN;

    % read names and bounding rectangle limits

    for n=1:924
        fnames{n,1} = fscanf(fid,'%s',1);
        YMIN(n) = fscanf(fid,'%d',1);
        YMAX(n) = fscanf(fid,'%d',1);
        XMIN(n) = fscanf(fid,'%d',1);
        XMAX(n) = fscanf(fid,'%d',1);
        qnames{n,1} = fscanf(fid,'%s',1);
    end
    fclose(fid);

    do = ...
        find( ...
        (...
        (latlim(1) <= YMIN & latlim(2) >= YMAX) | ... % tile is completely within region
        (latlim(1) >= YMIN & latlim(2) <= YMAX) | ... % region is completely within tile
        (latlim(1) >  YMIN & latlim(1) <  YMAX) | ... % min of region is on tile
        (latlim(2) >  YMIN & latlim(2) <  YMAX)   ... % max of region is on tile
        ) ...
        &...
        (...
        (lonlim(1) <= XMIN & lonlim(2) >= XMAX) | ... % tile is completely within region
        (lonlim(1) >= XMIN & lonlim(2) <= XMAX) | ... % region is completely within tile
        (lonlim(1) >  XMIN & lonlim(1) <  XMAX) | ... % min of region is on tile
        (lonlim(2) >  XMIN & lonlim(2) <  XMAX)   ... % max of region is on tile
        )...
        );

    if ~isempty(do)
        fname = fnames(do);
        qname = qnames(do);
    else
        fname = [];
        qname = [];
    end

author: CelsoReyes


text:
function pos2 = nutate (tjd, pos1)

% this function nutates equatorial rectangular coordinates from
% the mean dynamical equator and equinox of epoch to the true
% equator and equinox of epoch. see explanatory supplement to the
% astronomical almanac, pp. 114-115.

% input

%  tjd  = tdb julian date of epoch

%  pos1 = position vector, geocentric equatorial rectangular
%         coordinates, referred to mean dynamical equator and
%         equinox of epoch

% output

%  pos2 = position vector, geocentric equatorial rectangular
%         coordinates, referred to true equator and equinox of epoch

% note:  if tjd is negative, inverse nutation (true to mean) is applied.

% ported from NOVAS 3.0

%%%%%%%%%%%%%%%%%%%%%%%

seccon = 180.d0 * 3600.d0 / pi;

tjd1 = abs(tjd);

[oblm, oblt, eqeq, dpsi, deps] = etilt (tjd1);

oblm = oblm * 3600.d0 / seccon;

oblt = oblt * 3600.d0 / seccon;

dpsi = dpsi / seccon;

deps = deps / seccon;

cobm = cos(oblm);

sobm = sin(oblm);

cobt = cos(oblt);

sobt = sin(oblt);

cpsi = cos(dpsi);

spsi = sin(dpsi);

% compute elements of nutation rotation matrix

xx =  cpsi;
yx = -spsi * cobm;
zx = -spsi * sobm;

xy =  spsi * cobt;
yy =  cpsi * cobm * cobt + sobm * sobt;
zy =  cpsi * sobm * cobt - cobm * sobt;

xz =  spsi * sobt;
yz =  cpsi * cobm * sobt - sobm * cobt;
zz =  cpsi * sobm * sobt + cobm * cobt;

if (tjd < 0.0d0)

    % perform rotation from true to mean

    pos2(1) = xx * pos1(1) + xy * pos1(2) + xz * pos1(3);

    pos2(2) = yx * pos1(1) + yy * pos1(2) + yz * pos1(3);

    pos2(3) = zx * pos1(1) + zy * pos1(2) + zz * pos1(3);

else

    % perform rotation from mean to true

    pos2(1) = xx * pos1(1) + yx * pos1(2) + zx * pos1(3);

    pos2(2) = xy * pos1(1) + yy * pos1(2) + zy * pos1(3);

    pos2(3) = xz * pos1(1) + yz * pos1(2) + zz * pos1(3);

end

author: goGPS-Project