1## Copyright (C) 2013 Nir Krakauer
2##
3## This program is free software; you can redistribute it and/or modify
4## it under the terms of the GNU General Public License as published by
5## the Free Software Foundation; either version 2 of the License, or
6## (at your option) any later version.
7##
8## This program is distributed in the hope that it will be useful,
9## but WITHOUT ANY WARRANTY; without even the implied warranty of
10## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11## GNU General Public License for more details.
12##
13## You should have received a copy of the GNU General Public License
14## along with this program; If not, see <http://www.gnu.org/licenses/>.
15
16## -*- texinfo -*-
17## @deftypefn{Function File}{[@var{yi}] =} tps_val(@var{x}, @var{coefs}, @var{xi}, @var{vectorize}=true)
18##
19## Evaluates a thin plate spline at given points @*
20## @var{xi}
21##
22## @var{coefs} should be the vector of fitted coefficients returned from @code{tpaps(x, y, [p])}
23##
24## @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{coefs} should be @var{n} + @var{d} + 1 by 1; @var{xi} should be @var{k} by @var{d}
25##
26## The logical argument @var{vectorize} controls whether an @var{k} by @var{n} by @var{d} intermediate array is formed to speed up computation (the default) or whether looping is used to economize on memory
27##
28## The returned @var{yi} will be @var{k} by 1
29##
30## See the documentation to @code{tpaps} for more information
31##
32## @end deftypefn
33## @seealso{tpaps, tps_val_der}
34
35## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu>
36
37function [yi]=tps_val(x,coefs,xi,vectorize=true)
38
39
40  [n d] = size(x); #d: number of dimensions; n: number of points
41  k = size(xi, 1); #number of points for which to find the spline function value
42
43  #form of the Green's function for solutions
44  G = @(r) merge(r == 0, 0, r .^ 2 .* log(r));
45
46 a = coefs(1:n);
47 b = coefs((n+1):end);
48
49
50 yi = [ones(k, 1) xi] * b;
51 if vectorize
52   if d == 1
53    yi = yi + G(abs(x' - xi)) * a;
54   else
55    yi = yi + G(sqrt(sumsq((reshape(x, 1, n, d) - reshape(xi, k, 1, d)), 3))) * a;
56   endif
57 else
58   dist = @(x1, x2) norm(x2 - x1, 2, "rows"); #Euclidian distance between points in d-dimensional space
59   warn_state = warning ("query", "Octave:broadcast").state;
60   warning ("off", "Octave:broadcast"); #turn off warning message for automatic broadcasting when dist is called
61   unwind_protect
62  if k > n ##choose from either of two ways of computing the values of the thin plate spline at xi
63    for i = 1:n
64      yi = yi + a(i)*G(dist(x(i, :), xi));
65    endfor
66  else
67    for i = 1:k
68      yi(i) = yi(i) + dot(a, G(dist(x, xi(i, :))));
69    endfor
70  endif
71   unwind_protect_cleanup
72   warning (warn_state, "Octave:broadcast");
73   end_unwind_protect
74 endif
75
76
77endfunction
78
79%!shared x,y,c,xi
80%! x = ([1:10 10.5 11.3])'; y = sin(x);
81%! c = tpaps(x,y,1);
82%!assert (tpaps(x,y,1,x), tps_val(x,c,x), 100*eps);
83%! x = ((1 ./ (1:100))' - 0.5) * ([0.2 0.6]);
84%! y = x(:, 1) .^ 2 + x(:, 2) .^ 2;
85%! c = tpaps(x,y,1);
86%!assert (tpaps(x,y,1,x), tps_val(x,c,x), 100*eps);
87%!assert (tps_val(x,c,x,true), tps_val(x,c,x,false), 100*eps);
88