This file is indexed.

/usr/share/octave/packages/optim-1.3.0/expfit.m is in octave-optim 1.3.0-1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
## Copyright (C) 2000 Gert Van den Eynde <na.gvandeneynde@na-net.ornl.gov>
## Copyright (C) 2002 Rolf Fabian <fabian@tu-cottbus.de>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.

## USAGE  [alpha,c,rms] = expfit( deg, x1, h, y )
##
## Prony's method for non-linear exponential fitting
##
## Fit function:   \sum_1^{deg} c(i)*exp(alpha(i)*x)
##
## Elements of data vector y must correspond to
## equidistant x-values starting at x1 with stepsize h
##
## The method is fully compatible with complex linear
## coefficients c, complex nonlinear coefficients alpha
## and complex input arguments y, x1, non-zero h .
## Fit-order deg  must be a real positive integer.
##
## Returns linear coefficients c, nonlinear coefficients
## alpha and root mean square error rms. This method is
## known to be more stable than 'brute-force' non-linear
## least squares fitting.
##
## Example
##    x0 = 0; step = 0.05; xend = 5; x = x0:step:xend;
##    y = 2*exp(1.3*x)-0.5*exp(2*x);
##    error = (rand(1,length(y))-0.5)*1e-4;
##    [alpha,c,rms] = expfit(2,x0,step,y+error)
##
##  alpha =
##    2.0000
##    1.3000
##  c =
##    -0.50000
##     2.00000
##  rms = 0.00028461
##
## The fit is very sensitive to the number of data points.
## It doesn't perform very well for small data sets.
## Theoretically, you need at least 2*deg data points, but
## if there are errors on the data, you certainly need more.
##
## Be aware that this is a very (very,very) ill-posed problem.
## By the way, this algorithm relies heavily on computing the
## roots of a polynomial. I used 'roots.m', if there is
## something better please use that code.
##
## Demo for a complex fit-function:
## deg= 2; N= 20; x1= -(1+i), x= linspace(x1,1+i/2,N).';
## h = x(2) - x(1)
## y= (2+i)*exp( (-1-2i)*x ) + (-1+3i)*exp( (2+3i)*x );
## A= 5e-2; y+= A*(randn(N,1)+randn(N,1)*i); % add complex noise
## [alpha,c,rms]= expfit( deg, x1, h, y )

function [a,c,rms] = expfit(deg,x1,h,y)

% Check input
if deg<1, error('expfit: deg must be >= 1');     end
if  ~h,   error('expfit: vanishing stepsize h'); end

if ( N=length( y=y(:) ) ) < 2*deg
   error('expfit: less than %d samples',2*deg);
end

% Solve for polynomial coefficients
A = hankel( y(1:N-deg),y(N-deg:N) );
s = - A(:,1:deg) \ A(:,deg+1);

% Compose polynomial
p = flipud([s;1]);

% Compute roots
u = roots(p);

% nonlinear coefficients
a = log(u)/h;

% Compose second matrix A(i,j) = u(j)^(i-1)
A = ( ones(N,1) * u(:).' ) .^ ( [0:N-1]' * ones(1,deg) );

% Solve linear system
f = A\y;

% linear coefficients
c = f./exp( a*x1 );

% Compute rms of y - approx
% where approx(i) = sum_k ( c(k) * exp(x(i)*a(k)) )
if nargout > 2
   x = x1+h*[0:N-1];
   approx = exp( x(:) * a(:).' ) * c(:);
   rms = sqrt( sumsq(approx - y) );
end

endfunction

%!demo	 % same as in help - part
%! deg= 2; N= 20; x1= -(1+i), x= linspace(x1,1+i/2,N).';
%! h = x(2) - x(1)
%! y= (2+i)*exp( (-1-2i)*x ) + (-1+3i)*exp( (2+3i)*x );
%! A= 5e-2; y+= A*(randn(N,1)+randn(N,1)*i); % add complex noise
%! [alpha,c,rms]= expfit( deg, x1, h, y )

%!demo	 % demo for stepsize with negative real part
%! deg= 2; N= 20; x1= +3*(1+i), x= linspace(x1,1+i/3,N).';
%! h = x(2) - x(1)
%! y= (2+i)*exp( (-1-2i)*x ) + (-1+3i)*exp( (2+3i)*x );
%! A= 5e-2; y+= A*(randn(N,1)+randn(N,1)*i); % add complex noise
%! [alpha,c,rms]= expfit( deg, x1, h, y )

%!demo
%! x0 = 1.5; step = 0.05; xend = 5;
%! a = [1.3, 2]';
%! c = [2, -0.5]';
%! v = 1e-4;
%! 
%! x = x0:step:xend;
%! y = exp (x(:) * a(:).') * c(:);
%! err = randn (size (y)) * v;
%! plot (x, y + err);
%! 
%! [a_out, c_out, rms] = expfit (2, x0, step, y+err)