This file is indexed.

/usr/share/octave/packages/splines-1.2.6/tpaps.m is in octave-splines 1.2.6-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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
## Copyright (C) 2013 Nir Krakauer
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn{Function File}{[@var{yi} @var{p}] =} tpaps(@var{x}, @var{y}, @var{p}, @var{xi})
## @deftypefnx{Function File}{[@var{coefs} @var{p}] =} tpaps(@var{x}, @var{y}, @var{p}, [])
##
## Thin plate smoothing of scattered values in multi-D @*
## approximately interpolate [@var{x},@var{y}] at @var{xi}
##
## The chosen thin plate spline minimizes the sum of squared deviations from the given points plus a penalty term proportional to the curvature of the spline function
##
## @var{x} should be @var{n} by @var{d} in size, where @var{n} is the number of points and @var{d} the number of dimensions; @var{y} and @var{w} should be @var{n} by 1; @var{xi} should be @var{k} by @var{d}; the points in @var{x} should be distinct
##
## @table @asis
## @item @var{p}=0
##       maximum smoothing: flat surface
## @item @var{p}=1
##       no smoothing: interpolation
## @item @var{p}<0 or not given
##       an intermediate amount of smoothing is chosen (such that the smoothing term and the interpolation term are of the same magnitude)
## @end table
## 
## If @var{xi} is not specified, returns a vector @var{coefs} of the @var{n} + @var{d} + 1 fitted thin plate spline coefficients.
## Given @var{coefs}, the value of the thin-plate spline at any @var{xi} can be determined with @code{tps_val}
##  
## Note: Computes the pseudoinverse of an @var{n} by @var{n} matrix, so not recommended for very large @var{n}
##
## Example usages:
## @example
## x = ([1:10 10.5 11.3])'; y = sin(x); xi = (0:0.1:12)';
## yi = tpaps(x, y, 0.5, xi); 
## plot(x, y, xi, yi)
## @end example
## 
## @example
## x = rand(100, 2)*2 - 1; 
## y = x(:, 1) .^ 2 + x(:, 2) .^ 2;
## scatter(x(:, 1), x(:, 2), 10, y, "filled")
## [x1 y1] = meshgrid((-1:0.2:1)', (-1:0.2:1)');
## xi = [x1(:) y1(:)];
## yi = tpaps(x, y, 1, xi);
## contourf(x1, y1, reshape(yi, 11, 11))
## @end example
## 
## Reference: 
## David Eberly (2011), Thin-Plate Splines, www.geometrictools.com/Documentation/ThinPlateSplines.pdf
## Bouhamidi, A. (2005) Weighted thin plate splines, Analysis and Applications, 3: 297-324
##
## @end deftypefn
## @seealso{csaps, tps_val}

## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu>

function [ret, p]=tpaps(x,y,p,xi)


  if(nargin < 4)
    xi = [];
    if(nargin < 3) || p < 0
      p = [];
    endif
  endif


  [n d] = size(x); #d: number of dimensions; n: number of points [y should be n*1]

  dist = @(x1, x2) norm(x2 - x1, 2, "rows"); #Euclidian distance between points in d-dimensional space
  
  #form of the Green's function for solutions
  G = @(r) merge(r == 0, 0, r .^ 2 .* log(r));
 
  N = [ones(n, 1) x]; 
  
  if p == 0 #infinite regularization (no curvature allowed), solution is a regression plane
  
    b = N \ y;
    
    a = zeros(n, 1);
    
  else  
  
    #coefficient matrices
    #need pairwise distances between points
    M = zeros(n);
     warn_state = warning ("query", "Octave:broadcast").state;
     warning ("off", "Octave:broadcast"); #turn off warning message for automatic broadcasting when dist is called
     unwind_protect
    for i = 1:(n-1) #M is symmetric, so only need to compute half
      M(i, (i+1):n) = G(dist(x(i, :), x((i+1):n, :)));
    endfor
     unwind_protect_cleanup
     warning (warn_state, "Octave:broadcast");
     end_unwind_protect    
    M = M + M';
    if isempty(p) #choose an intermediate value for the regularization parameter
      lambda = mean(spdiags(M, 1:min(n, 3))(:));
      p = 1 / (lambda + 1);
    else #use the given value
      lambda = (1 - p) / p;
    endif
    
    M = M + lambda*eye(n); #add regularization term
   
    M_inv = pinv(M);
   
    b = pinv(N' * M_inv * N) * N' * M_inv * y;
    
    a = M_inv * (y - N*b);
  
  endif
  
    
  if isempty(xi) #return the coefficients
      ret = [a' b']';
  else #return the thin plate spline values at xi
      k = size(xi, 1);
      ret = [ones(k, 1) xi] * b;
      if p ~= 0     
         warn_state = warning ("query", "Octave:broadcast").state;
         warning ("off", "Octave:broadcast"); #turn off warning message for automatic broadcasting when dist is called
         unwind_protect
        if k > n ##choose from either of two ways of computing the values of the thin plate spline at xi
          for i = 1:n
            ret = ret + a(i)*G(dist(x(i, :), xi));
          endfor
        else
          for i = 1:k
            ret(i) = ret(i) + dot(a, G(dist(x, xi(i, :))));
          endfor
        endif
         unwind_protect_cleanup
         warning (warn_state, "Octave:broadcast");
         end_unwind_protect  
      endif
  endif


endfunction

%!shared x,y
%! x = ([1:10 10.5 11.3])'; y = sin(x);
%!assert (tpaps(x,y,1,x), y, 1E-13);
%! x = rand(100, 2)*2 - 1; 
%! y = x(:, 1) .^ 2 + x(:, 2) .^ 2;
%!assert (tpaps(x,y,1,x), y, 1E-11);