X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fdata-smoothing-1.3.0%2Frgdtsmcore.m;fp=octave_packages%2Fdata-smoothing-1.3.0%2Frgdtsmcore.m;h=bb63511d561681ebcb68dbccd665095d6899bdef;hp=0000000000000000000000000000000000000000;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hpb=1705066eceaaea976f010f669ce8e972f3734b05

diff --git a/octave_packages/data-smoothing-1.3.0/rgdtsmcore.m b/octave_packages/data-smoothing-1.3.0/rgdtsmcore.m
new file mode 100644
index 0000000..bb63511
--- /dev/null
+++ b/octave_packages/data-smoothing-1.3.0/rgdtsmcore.m
@@ -0,0 +1,166 @@
+## Copyright (C) 2008 Jonathan Stickel <jonathan.stickel@nrel.gov>
+##
+## 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/>.
+
+## -*- texinfo -*-
+##@deftypefn {Function File} {[@var{yhat}, @var{v}] =} rgdtsmcore (@var{x}, @var{y}, @var{d}, @var{lambda}, [@var{options}])
+##
+## Smooths @var{y} vs. @var{x} values by Tikhonov regularization.
+## Although this function can be used directly, the more feature rich
+## function "regdatasmooth" should be used instead.  In addition to
+## @var{x} and @var{y}, required input includes the smoothing derivative
+## @var{d} and the regularization parameter @var{lambda}.  The smooth
+## y-values are returned as @var{yhat}.  The generalized cross
+## validation variance @var{v} may also be returned.
+##
+## Note:  the options have changed!
+## Currently supported input options are (multiple options are allowed):
+##
+## @table @code
+## @item "xhat", @var{vector}
+##  A vector of x-values to use for the smooth curve; must be
+##  monotonically increasing and must at least span the data
+## @item "weights", @var{vector}
+##  A vector of weighting values for fitting each point in the data.
+## @item "relative"
+##  use relative differences for the goodnes of fit term.  Conflicts
+##  with the "weights" option.
+## @item "midpointrule"
+##  use the midpoint rule for the integration terms rather than a direct
+##  sum; this option conflicts with the option "xhat"
+## @end table
+##
+## References:  Anal. Chem. (2003) 75, 3631; AIChE J. (2006) 52, 325
+## @seealso{regdatasmooth}
+## @end deftypefn
+
+
+function [yhat, v] = rgdtsmcore (x, y, d, lambda, varargin)
+
+  if (nargin < 4)
+    print_usage;
+  endif
+
+  ## Defaults if not provided
+  xhatprov = 0;
+  xhat = x;
+  weights = 0;
+  relative = 0;
+  midpr = 0;
+  
+  ## parse the provided options
+  if ( length(varargin) )
+    for i = 1:length(varargin)
+      arg = varargin{i};
+      if ischar(arg)
+        switch arg
+          case "xhat"
+            xhatprov = 1;
+            xhat = varargin{i+1};
+          case "weights"
+            weights = 1;
+            weightv = varargin{i+1};
+          case "relative"
+            relative = 1;
+          case "midpointrule"
+            midpr = 1;
+          otherwise
+            printf("Option '%s' is not implemented;\n", arg)
+        endswitch
+      endif
+    endfor
+  endif
+  if (xhatprov && midpr)
+    warning("midpointrule is currently not used if xhat is provided (since x,y may be scattered)")
+    midpr = 0;
+  endif
+  if (weights && relative)
+    warning("relative differences is not used if a weighting vector is provided")
+  endif
+  
+  N = length(x);
+  Nhat = length(xhat);
+  
+  ## test that xhat is increasing
+  if !all(diff(xhat)>0)
+    if xhatprov
+      error("xhat must be monotonically increasing")
+    else
+      error("x must be monotonically increasing if xhat is not provided")
+    endif
+  endif
+  ## test that xhat spans x
+  if ( min(x) < min(xhat) || max(xhat) < max(x) )
+    error("xhat must at least span the data")
+  endif
+
+  ## construct M, D
+  M = speye(Nhat);
+  idx = interp1(xhat,1:Nhat,x,"nearest"); # works for unequally spaced xhat
+  M = M(idx,:);
+  D = ddmat(xhat,d);
+
+  ## construct "weighting" matrices W and U
+  if (weights)
+    ## use arbitrary weighting as provided
+    W = diag(weightv);
+  elseif (relative)
+    ## use relative differences
+    Yinv = sparse(diag(1./y));
+    W = Yinv^2;
+  else
+    W = speye(N);
+  endif
+  ## use midpoint rule integration (rather than simple sums)
+  if (midpr)
+    Bhat = sparse(diag(-ones(N-1,1),-1)) + sparse(diag(ones(N-1,1),1));
+    Bhat(1,1) = -1;
+    Bhat(N,N) = 1;
+    B = 1/2*sparse(diag(Bhat*x));
+    if ( floor(d/2) == d/2 ) # test if d is even
+      dh = d/2;
+      Btilda = B(dh+1:N-dh,dh+1:N-dh);
+    else # d is odd
+      dh = ceil(d/2);
+      Btilda = B(dh:N-dh,dh:N-dh);
+    endif
+    W = W*B;
+    U = Btilda;
+  else
+    ## W = W*speye(Nhat);
+    U = speye(Nhat-d);
+  endif
+  
+  ## Smoothing
+  delta = trace(D'*D)/Nhat^(2+d);  # using "relative" or other weighting affects this!
+  yhat = (M'*W*M + lambda*delta^(-1)*D'*U*D) \ M'*W*y;
+  #[R,P,S] = splchol(M'*W*M + lambda*delta^(-1)*D'*U*D);
+  #yhat = S*(R'\(R\(S'*M'*W*y)));
+  
+  ## Computation of hat diagonal and cross-validation
+  if (nargout > 1)
+    ## from AIChE J. (2006) 52, 325
+    ## note: chol factorization does not help speed up the computation of H;
+    ## should implement Eiler's partial H computation if many point smoothing by GCV is needed
+    ##H = M*(S*(R'\(R\(S'*M'*W))));
+    H = M*((M'*W*M + lambda*delta^(-1)*D'*U*D)\M'*W);
+    ## note: this is variance, squared of the standard error that Eilers uses
+    v = (M*yhat - y)'*(M*yhat - y)/N / (1 - trace(H)/N)^2;
+  endif
+  
+  ## test mapping
+  ##figure(5)
+  ##plot(x,y,'o',x,M*yhat,'x')
+
+endfunction