Script 'mail_helper' called by obssrc
Hello community,
here is the log from the commit of package octave-forge-splines for
openSUSE:Factory checked in at 2021-03-17 20:16:47
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Comparing /work/SRC/openSUSE:Factory/octave-forge-splines (Old)
and /work/SRC/openSUSE:Factory/.octave-forge-splines.new.2401 (New)
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Package is "octave-forge-splines"
Wed Mar 17 20:16:47 2021 rev:6 rq:879590 version:1.3.4
Changes:
--------
---
/work/SRC/openSUSE:Factory/octave-forge-splines/octave-forge-splines.changes
2019-11-28 10:15:22.131637926 +0100
+++
/work/SRC/openSUSE:Factory/.octave-forge-splines.new.2401/octave-forge-splines.changes
2021-03-17 20:20:06.463330051 +0100
@@ -1,0 +2,6 @@
+Sat Mar 13 13:48:18 UTC 2021 - Atri Bhattacharya <[email protected]>
+
+- Update to version 1.3.4:
+ * New functions regularization, regularization2D
+
+-------------------------------------------------------------------
Old:
----
splines-1.3.3.tar.gz
New:
----
splines-1.3.4.tar.gz
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Other differences:
------------------
++++++ octave-forge-splines.spec ++++++
--- /var/tmp/diff_new_pack.eOPCTN/_old 2021-03-17 20:20:06.951330718 +0100
+++ /var/tmp/diff_new_pack.eOPCTN/_new 2021-03-17 20:20:06.955330723 +0100
@@ -1,7 +1,7 @@
#
# spec file for package octave-forge-splines
#
-# Copyright (c) 2019 SUSE LLC
+# Copyright (c) 2021 SUSE LLC
#
# All modifications and additions to the file contributed by third parties
# remain the property of their copyright owners, unless otherwise agreed
@@ -18,7 +18,7 @@
%define octpkg splines
Name: octave-forge-%{octpkg}
-Version: 1.3.3
+Version: 1.3.4
Release: 0
Summary: Additional spline functions for Octave
License: GPL-2.0-or-later AND GPL-3.0-or-later AND SUSE-Public-Domain
++++++ splines-1.3.3.tar.gz -> splines-1.3.4.tar.gz ++++++
diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn'
'--exclude=.svnignore' old/splines-1.3.3/CITATION new/splines-1.3.4/CITATION
--- old/splines-1.3.3/CITATION 2019-10-23 20:45:53.000000000 +0200
+++ new/splines-1.3.4/CITATION 2021-02-24 02:14:16.000000000 +0100
@@ -2,6 +2,6 @@
author = {Krakauer, Nir Y. and others},
title = {Splines Package for GNU Octave},
url = {http://octave.sourceforge.net/splines},
- version = {1.3.2},
- date = {2016-11-23},
+ version = {1.3.4},
+ date = {2021-02-23},
}
diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn'
'--exclude=.svnignore' old/splines-1.3.3/DESCRIPTION
new/splines-1.3.4/DESCRIPTION
--- old/splines-1.3.3/DESCRIPTION 2019-10-23 20:45:53.000000000 +0200
+++ new/splines-1.3.4/DESCRIPTION 2021-02-24 02:14:16.000000000 +0100
@@ -1,6 +1,6 @@
Name: splines
-Version: 1.3.3
-Date: 2019-10-17
+Version: 1.3.4
+Date: 2021-02-23
Author: various authors
Maintainer: Nir Krakauer <[email protected]>
Title: Splines.
diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn'
'--exclude=.svnignore' old/splines-1.3.3/INDEX new/splines-1.3.4/INDEX
--- old/splines-1.3.3/INDEX 2019-10-23 20:45:53.000000000 +0200
+++ new/splines-1.3.4/INDEX 2021-02-24 02:14:16.000000000 +0100
@@ -10,6 +10,8 @@
fnder
fnplt
fnval
+ regularization
+ regularization2D
tpaps
tps_val
tps_val_der
diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn'
'--exclude=.svnignore' old/splines-1.3.3/Makefile new/splines-1.3.4/Makefile
--- old/splines-1.3.3/Makefile 1970-01-01 01:00:00.000000000 +0100
+++ new/splines-1.3.4/Makefile 2021-02-24 02:14:16.000000000 +0100
@@ -0,0 +1,235 @@
+## Copyright 2015-2016 Carn?? Draug
+## Copyright 2015-2016 Oliver Heimlich
+## Copyright 2017 Julien Bect <[email protected]>
+## Copyright 2017 Olaf Till <[email protected]>
+##
+## Copying and distribution of this file, with or without modification,
+## are permitted in any medium without royalty provided the copyright
+## notice and this notice are preserved. This file is offered as-is,
+## without any warranty.
+
+## Some basic tools (can be overriden using environment variables)
+SED ?= sed
+TAR ?= tar
+GREP ?= grep
+CUT ?= cut
+TR ?= tr
+
+## Note the use of ':=' (immediate set) and not just '=' (lazy set).
+## http://stackoverflow.com/a/448939/1609556
+package := $(shell $(GREP) "^Name: " DESCRIPTION | $(CUT) -f2 -d" " | \
+$(TR) '[:upper:]' '[:lower:]')
+version := $(shell $(GREP) "^Version: " DESCRIPTION | $(CUT) -f2 -d" ")
+
+## These are the paths that will be created for the releases.
+target_dir := target
+release_dir := $(target_dir)/$(package)-$(version)
+release_tarball := $(target_dir)/$(package)-$(version).tar.gz
+html_dir := $(target_dir)/$(package)-html
+html_tarball := $(target_dir)/$(package)-html.tar.gz
+## Using $(realpath ...) avoids problems with symlinks due to bug
+## #50994 in Octaves scripts/pkg/private/install.m. But at least the
+## release directory above is needed in the relative form, for 'git
+## archive --format=tar --prefix=$(release_dir).
+real_target_dir := $(realpath .)/$(target_dir)
+installation_dir := $(real_target_dir)/.installation
+package_list := $(installation_dir)/.octave_packages
+install_stamp := $(installation_dir)/.install_stamp
+
+## These can be set by environment variables which allow to easily
+## test with different Octave versions.
+ifndef OCTAVE
+OCTAVE := octave
+endif
+OCTAVE := $(OCTAVE) --no-gui --silent --norc
+MKOCTFILE ?= mkoctfile
+
+## Command used to set permissions before creating tarballs
+FIX_PERMISSIONS ?= chmod -R a+rX,u+w,go-w,ug-s
+
+## Detect which VCS is used
+vcs := $(if $(wildcard .hg),hg,$(if $(wildcard .git),git,unknown))
+ifeq ($(vcs),hg)
+release_dir_dep := .hg/dirstate
+endif
+ifeq ($(vcs),git)
+release_dir_dep := .git/index
+endif
+
+
+## .PHONY indicates targets that are not filenames
+## (https://www.gnu.org/software/make/manual/html_node/Phony-Targets.html)
+.PHONY: help
+
+## make will display the command before runnning them. Use @command
+## to not display it (makes specially sense for echo).
+help:
+ @echo "Targets:"
+ @echo " dist - Create $(release_tarball) for release."
+ @echo " html - Create $(html_tarball) for release."
+ @echo " release - Create both of the above and show md5sums."
+ @echo " install - Install the package in $(installation_dir), where
it is not visible in a normal Octave session."
+ @echo " check - Execute package tests."
+ @echo " doctest - Test the help texts with the doctest package."
+ @echo " run - Run Octave with the package installed in
$(installation_dir) in the path."
+ @echo " clean - Remove everything made with this Makefile."
+
+
+##
+## Recipes for release tarballs (package + html)
+##
+
+.PHONY: release dist html clean-tarballs clean-unpacked-release
+
+## To make a release, build the distribution and html tarballs.
+release: dist html
+ md5sum $(release_tarball) $(html_tarball)
+ @echo "Upload @ https://sourceforge.net/p/octave/package-releases/new/";
+ @echo " and note the changeset the release corresponds to"
+
+## dist and html targets are only PHONY/alias targets to the release
+## and html tarballs.
+dist: $(release_tarball)
+html: $(html_tarball)
+
+## An implicit rule with a recipe to build the tarballs correctly.
+%.tar.gz: %
+ $(TAR) -c -f - --posix -C "$(target_dir)/" "$(notdir $<)" | gzip -9n >
"$@"
+
+clean-tarballs:
+ @echo "## Cleaning release tarballs (package + html)..."
+ -$(RM) $(release_tarball) $(html_tarball)
+ @echo
+
+## Create the unpacked package.
+##
+## Notes:
+## * having ".hg/dirstate" (or ".git/index") as a prerequesite means it is
+## only rebuilt if we are at a different commit.
+## * the variable RM usually defaults to "rm -f"
+## * having this recipe separate from the one that makes the tarball
+## makes it easy to have packages in alternative formats (such as zip)
+## * note that if a commands needs to be run in a specific directory,
+## the command to "cd" needs to be on the same line. Each line restores
+## the original working directory.
+$(release_dir): $(release_dir_dep)
+ -$(RM) -r "$@"
+ifeq (${vcs},hg)
+ hg archive --exclude ".hg*" --type files "$@"
+endif
+ifeq (${vcs},git)
+ git archive --format=tar --prefix="$@/" HEAD | $(TAR) -x
+ $(RM) "$@/.gitignore"
+endif
+## Don't fall back to run the supposed necessary contents of
+## 'bootstrap' here. Users are better off if they provide
+## 'bootstrap'. Administrators, checking build reproducibility, can
+## put in the missing 'bootstrap' file if they feel they know its
+## necessary contents.
+ifneq (,$(wildcard src/bootstrap))
+ cd "$@/src" && ./bootstrap && $(RM) -r "autom4te.cache"
+endif
+## Uncomment this if your src/Makefile.in has these targets for
+## pre-building something for the release (e.g. documentation).
+# cd "$@/src" && ./configure && $(MAKE) prebuild && \
+# $(MAKE) distclean && $(RM) Makefile
+##
+ ${FIX_PERMISSIONS} "$@"
+
+run_in_place = $(OCTAVE) --eval ' pkg ("local_list", "$(package_list)"); ' \
+ --eval ' pkg ("load", "$(package)"); '
+
+html_options = --eval 'options = get_html_options ("octave-forge");'
+## Uncomment this for package documentation.
+# html_options = --eval 'options = get_html_options ("octave-forge");' \
+# --eval 'options.package_doc = "$(package).texi";'
+$(html_dir): $(install_stamp)
+ $(RM) -r "$@";
+ $(run_in_place) \
+ --eval ' pkg load generate_html; ' \
+ $(html_options) \
+ --eval ' generate_package_html ("$(package)", "$@", options); ';
+ $(FIX_PERMISSIONS) "$@";
+
+clean-unpacked-release:
+ @echo "## Cleaning unpacked release tarballs (package + html)..."
+ -$(RM) -r $(release_dir) $(html_dir)
+ @echo
+
+##
+## Recipes for installing the package.
+##
+
+.PHONY: install clean-install
+
+octave_install_commands = \
+' llist_path = pkg ("local_list"); \
+ mkdir ("$(installation_dir)"); \
+ load (llist_path); \
+ local_packages(cellfun (@ (x) strcmp ("$(package)", x.name),
local_packages)) = []; \
+ save ("$(package_list)", "local_packages"); \
+ pkg ("local_list", "$(package_list)"); \
+ pkg ("prefix", "$(installation_dir)", "$(installation_dir)"); \
+ pkg ("install", "-local", "-verbose", "$(release_tarball)"); '
+
+## Install unconditionally. Maybe useful for testing installation with
+## different versions of Octave.
+install: $(release_tarball)
+ @echo "Installing package under $(installation_dir) ..."
+ $(OCTAVE) --eval $(octave_install_commands)
+ touch $(install_stamp)
+
+## Install only if installation (under target/...) is not current.
+$(install_stamp): $(release_tarball)
+ @echo "Installing package under $(installation_dir) ..."
+ $(OCTAVE) --eval $(octave_install_commands)
+ touch $(install_stamp)
+
+clean-install:
+ @echo "## Cleaning installation under $(installation_dir) ..."
+ -$(RM) -r $(installation_dir)
+ @echo
+
+
+##
+## Recipes for testing purposes
+##
+
+.PHONY: run doctest check
+
+## Start an Octave session with the package directories on the path for
+## interactice test of development sources.
+run: $(install_stamp)
+ $(run_in_place) --persist
+
+## Test example blocks in the documentation. Needs doctest package
+## https://octave.sourceforge.io/doctest/index.html
+doctest: $(install_stamp)
+ $(run_in_place) --eval 'pkg load doctest;'
\
+ --eval "targets = '$(shell (ls inst; ls src | $(GREP) .oct) | $(CUT)
-f2 -d@ | $(CUT) -f1 -d.)';" \
+ --eval "targets = strsplit (targets, ' '); doctest (targets);"
+
+
+## Test package.
+octave_test_commands = \
+' args = {"inst", "src"}; \
+ args(cellfun (@ (x) isempty (a = stat (x)) || ! S_ISDIR (a.mode), args)) =
[]; \
+ if (isempty (args)) error ("no \"inst\" or \"src\" directory"); exit (1); \
+ else cellfun(@runtests, args); endif '
+check: $(install_stamp)
+ $(run_in_place) --eval $(octave_test_commands)
+
+
+##
+## CLEAN
+##
+
+.PHONY: clean
+
+clean: clean-tarballs clean-unpacked-release clean-install
+ @echo "## Removing target directory (if empty)..."
+ -rmdir $(target_dir)
+ @echo
+ @echo "## Cleaning done"
+ @echo
+
diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn'
'--exclude=.svnignore' old/splines-1.3.3/NEWS new/splines-1.3.4/NEWS
--- old/splines-1.3.3/NEWS 2019-10-23 20:45:53.000000000 +0200
+++ new/splines-1.3.4/NEWS 2021-02-24 02:14:16.000000000 +0100
@@ -1,3 +1,9 @@
+Summary of important user-visible changes for splines 1.3.4:
+-------------------------------------------------------------------
+
+ ** New functions regularization, regularization2D
+
+
Summary of important user-visible changes for splines 1.3.3:
-------------------------------------------------------------------
diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn'
'--exclude=.svnignore' old/splines-1.3.3/inst/regularization.m
new/splines-1.3.4/inst/regularization.m
--- old/splines-1.3.3/inst/regularization.m 1970-01-01 01:00:00.000000000
+0100
+++ new/splines-1.3.4/inst/regularization.m 2021-02-24 02:14:16.000000000
+0100
@@ -0,0 +1,175 @@
+## Copyright (C) 2019 Andreas Stahel
+##
+## 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
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{grid},@var{u}] =} regularization
(@var{data}, @var{interval}, @var{N}, @var{F1})
+## @deftypefnx {Function File} {[@var{grid},@var{u}] =} regularization
(@var{data}, @var{interval}, @var{N}, @var{F1}, @var{F2})
+##
+## Apply a Tikhonov regularization, the functional to be minimized is@*
+## @var{F} = @var{FD} + @var{lambda1}*@var{F1} + @var{lambda2}*@var{F2} @*
+## = sum_(i=1)^M (y_i-u(x_i))^2
+## + @var{lambda1}*int_a^b (u'(x) - @var{g1}(x))^2 dx
+## + @var{lambda2}*int_a^b (u''(x) - @var{g2}(x))^2 dx
+##
+## With @var{lambda1} = 0 and @var{G2}(x) = 0 this leads to a smoothing spline.
+##
+## Parameters:
+## @itemize
+## @item @var{data} is a M*2 matrix with the x values in the first column and
the y values in the second column.
+## @item @var{interval} = [a,b] is the interval on which the regularization is
applied.
+## @item @var{N} is the number of subintervals of equal length. @var{grid}
will consist of @var{N+1} grid points.
+## @item @var{F1} is a structure containing the information on the first
regularization term, integrating the square of the first derivative.
+## @itemize
+## @item @var{F1.lambda} is the value of the regularization parameter
@var{lambda1}>=0.
+## @item @var{F1.g} is the function handle for the function @var{g1(x)}. If
not provided @var{G1=0} is used.
+## @end itemize
+## @item @var{F2} is a structure containing the information on the second
regularization term, integrating the square of the second derivative. If
@var{F2} is not provided @var{lambda2}=0 is assumed.
+## @itemize
+## @item @var{F2.lambda} is the value of the regularization parameter
@var{lambda2}>=0.
+## @item @var{F2.g} is the function handle for the function @var{g2(x)}. If
not provided @var{G2=0} is used.
+## @end itemize
+## @end itemize
+##
+## Return values:
+## @itemize
+## @item @var{grid} is the grid on which @var{u} is evaluated. It consists of
@var{N+1} equidistant points on the @var{interval}.
+## @item @var{u} are the values of the regularized approximation to the
@var{data} evaluated at @var{grid}.
+## @end itemize
+## @seealso{csaps, regularization2D, demo regularization}
+## @end deftypefn
+
+## Author: Andreas Stahel <[email protected]>
+## Created: 2019-03-26
+
+function [grid,u] = regularization (data, interval, N, F1, F2)
+ a = interval (1); b = interval (2);
+ grid = linspace (a, b, N + 1)';
+ dx = grid (2) - grid (1);
+ x = data (:, 1);
+ ## select points in interval only
+ ind = find ((x >= a) .* (x <= b));
+ x = x (ind);
+ y = data (:, 2); y = y (ind);
+ M = length (x);
+ Interp = sparse (M, N + 1); ## interpolation matrix
+ pos = floor ((x - a) / dx) + 1;
+ theta = mod ((x - a) / dx, 1);
+ for ii = 1:M
+ if theta (ii) > 10*eps
+ Interp (ii, pos (ii)) = 1 - theta (ii);
+ Interp (ii, pos (ii) + 1) = theta (ii);
+ else
+ Interp (ii, pos (ii)) = 1;
+ endif
+ endfor
+ mat = Interp' * Interp;
+ rhs = (Interp' * y);
+ if isfield (F1, 'lambda')
+ A1 = spdiags ([-ones(N, 1), +ones(N, 1)], [0, 1], N, N + 1) / dx;
+ mat = mat + F1.lambda * dx * A1' * A1;
+ if isfield (F1, 'g')
+ g1 = F1.g (grid (1:end - 1) + dx / 2);
+ rhs = rhs + F1.lambda * dx * A1' * g1;
+ endif
+ endif
+ if exist ('F2')
+ A2 = spdiags (ones (N, 1) * [1, -2, 1], [0, 1, 2], N - 1, N + 1) / dx^2;
+ mat = mat + F2.lambda * dx * A2' * A2;
+ if isfield (F2, 'g')
+ g2 = F2.g (grid (2:end - 1));
+ rhs = rhs + F2.lambda * dx * A2' * g2;
+ endif
+ endif
+ u = mat \ rhs;
+endfunction
+
+%!demo
+%! N = 100; interval = [0,10];
+%! x = [3.2,4,5,5.2,5.6]'; y = x;
+%!
+%! clear F1 F2
+%! %% regularize towards slope 0.1, no smoothing
+%! F1.lambda = 1e-2; F1.g = @(x)0.1*ones(size(x));
+%! [grid,u1] = regularization([x,y],interval,N,F1);
+%! %% regularize towards slope 0.1, with some smoothing
+%! F2.lambda = 1*1e-3;
+%! [grid,u2] = regularization([x,y],interval,N,F1,F2);
+%!
+%! figure(1)
+%! plot(grid,u1,'b',grid,u2,'g',x,y,'*r')
+%! xlabel('x'); ylabel('solution');
+%! legend('regular1','regular2','data','location','northwest')
+
+%!demo
+%! N = 1000; interval = [0,pi];
+%! x = linspace( pi/4,3*pi/4,15)'; y = sin(x)+ 0.03*randn(size(x));
+%!
+%! clear F1 F2
+%! %% regularize by smoothing only
+%! F1.lambda = 0; F2.lambda = 1e-3;
+%! [grid,u1] = regularization([x,y],interval,N,F1,F2);
+%! %% regularize by smoothing and aim for slope 0
+%! F1.lambda = 1*1e-2;
+%! [grid,u2] = regularization([x,y],interval,N,F1,F2);
+%!
+%! figure(1)
+%! plot(grid,u1,'b',grid,u2,'g',x,y,'*r')
+%! xlabel('x'); ylabel('solution');
+%! legend('regular1','regular2','data','location','northwest')
+
+%!demo
+%! interval = [0,1];
+%! N = 400;
+%! x = rand(200,1);
+%! %% generate the data on four line segments, add some noise
+%! y = 2 - 2*x + (x>0.25) - 2*(x>0.5).*(x<0.65)+ 0.1*randn(length(x),1);
+%! clear F1
+%! %% apply regularization with three different parameters
+%! F1.lambda = 1e-3; [grid,u1] = regularization([x,y],interval,N,F1);
+%! F1.lambda = 1e-1; [grid,u2] = regularization([x,y],interval,N,F1);
+%! F1.lambda = 3e+0; [grid,u3] = regularization([x,y],interval,N,F1);
+%!
+%! figure(1); plot(grid,u1,'b',grid,u2,'g',grid,u3,'m',x,y,'+r')
+%! xlabel('x'); ylabel('solution')
+%! legend('\lambda_1=0.001','\lambda_1=0.1','\lambda_1=3','data')
+
+%!demo
+%! %% generate a smoothing spline, see also csaps() in the package splines
+%! N = 1000; interval = [0,10.3];
+%! x = [0 3 4 6 10]'; y = [0 1 0 1 0]';
+%! clear F2
+%! F2.lambda = 1e-2;
+%! %% apply regularization, the result is a smoothing spline
+%! [grid,u] = regularization([x,y],interval,N,0,F2);
+%!
+%! figure(1);
+%! plot(grid,u,'b',x,y,'*r')
+%! legend('spline','data')
+%! xlabel('x'); ylabel('solution')
+
+
+%!test
+%! data = [0,0;1,1;2,2];
+%! F1.lambda = 0.0; F2.lambda = 0.1;
+%! [grid,u] = regularization(data,[0,1],10,F1,F2);
+%! assert(norm(grid-u),0,1e-12)
+
+%!test
+%! x = linspace(-1,1,11);
+%! F1.lambda = 0.01;F2.g = @(x) 2*ones(size(x)); F2.lambda = 0.01;
+%! [grid,u] = regularization([x;x.^2]',[-2,2],20,F1,F2);
+%! assert(u(11),7.330959483903200e-03,1e-8)
+
diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn'
'--exclude=.svnignore' old/splines-1.3.3/inst/regularization2D.m
new/splines-1.3.4/inst/regularization2D.m
--- old/splines-1.3.3/inst/regularization2D.m 1970-01-01 01:00:00.000000000
+0100
+++ new/splines-1.3.4/inst/regularization2D.m 2021-02-24 02:14:16.000000000
+0100
@@ -0,0 +1,161 @@
+## Copyright (C) 2021 Andreas Stahel
+##
+## 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
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{grid},@var{u},@var{data_valid}] =}
regularization2D (@var{data}, @var{box}, @var{N}, @var{lambda1},@var{lambda2})
+##
+## Apply a Tikhonov regularization, the functional to be minimized is@*
+## @var{F} = @var{FD} + @var{lambda1} * @var{F1} + @var{lambda2} * @var{F2} @*
+## = sum_(i=1)^M (y_i-u(x_i))^2 +@*
+## + @var{lambda1} * dintegral (du/dx)^2+(du/dy)^2 dA +@*
+## + @var{lambda2} * dintegral (d^2u/dx^2)^2+(d^2u/dy^2)^2+2*(d^2u/dxdy) dA
+##
+## With @var{lambda1} = 0 and @var{lambda2}>0 this leads to a thin plate
smoothing spline.
+##
+## Parameters:
+## @itemize
+## @item @var{data} is a M*3 matrix with the (x,y) values in the first two
columns and the y values in the third column.@*
+## Only data points strictly inside the @var{box} are used
+## @item @var{box} = [x0,x1;y0,y1] is the rectangle x0<x<x1 and y0<y<y1 on
which the regularization is applied.
+## @item @var{N} = [N1,N2] determines the number of subintervals of equal
length. @var{grid} will consist of (@var{N1+1})x(@var{N2+1}) grid points.
+## @item @var{lambda1} >= 0 is the value of the first regularization parameter
+## @item @var{lambda2} > 0 is the value of the secondregularization parameter
+## @end itemize
+##
+## Return values:
+## @itemize
+## @item @var{grid} is the grid on which @var{u} is evaluated. It consists of
+## (@var{N1+1})x(@var{N2}+1) equidistant points on the the rectangle @var{box}.
+## @item @var{u} are the values of the regularized approximation to the
@var{data} evaluated on the @var{grid}.
+## @item @var{data_valid} returns the values data points used and the values
of the regularized function at these points
+## @end itemize
+## @seealso {tpaps, regularization, demo regularization2D}
+## @end deftypefn
+
+## Author: Andreas Stahel <[email protected]>
+## Created: 2021-01-13
+
+function [grid,u,data_valid] = regularization2D (data, box, N, lambda1,
lambda2)
+%% generate the grid
+ N = N+1; %%% now N is the number of grid points in either direction
+ x = linspace(box(1,1),box(1,2),N(1));
+ y = linspace(box(2,1),box(2,2),N(2));
+ [xx,yy] = meshgrid(x,y);
+ dx = diff(box(1,:))/(N(1)-1); dy = diff(box(2,:))/(N(2)-1);
+ grid.x = xx;
+ grid.y = yy;
+ x = data (:,1); y = data(:,2); z = data(:,3);
+ ## select points in box only
+ ind = find ((x > box(1,1)) .* (x < box(1,2)) .* (y > box(2,1)) .* (y <
box(2,2)));
+ x = x (ind); y = y (ind); z = z (ind);
+ %% generate the sparse interpolation matrix
+ M = length (x);
+ x_ind = floor((x-box(1,1))/dx);
+ %% xi = mod(x,dx)/dx;
+ xi = (x-box(1,1)-x_ind*dx)/dx;
+ y_ind = floor((y-box(2,1))/dy);
+ %%nu = mod(y,dy)/dy;
+ nu = (y-box(2,1)-y_ind*dy)/dy;
+ row = ones(4,1)*[1:M];
+ index_base = N(2)*x_ind+y_ind+1;
+ index = index_base + [0,N(2),1,N(2)+1]; index = index';
+ coeff = [(1-xi).*(1-nu),xi.*(1-nu),(1-xi).*nu,xi.*nu]; coeff = coeff';
+ Interp = sparse(row(:),index(:),coeff(:),M,N(1)*N(2));
+ mat = Interp' * Interp;
+ rhs = (Interp' * z);
+
+%%% derivative with respect to x
+ Dx = kron(spdiags(ones(N(1),1)*[-1 1],[0 1],N(1)-1,N(1))/dx,speye(N(2)));
+ Wx = ones(N(2),1); Wx(1) = 1/2; Wx(N(2)) = 1/2;
+ Wx = kron(speye(N(1)-1),diag(Wx))*dx*dy;
+%%% derivative with respect to y
+
+ Wy = ones(N(1),1); Wy(1) = 1/2; Wy(N(1)) = 1/2;
+ Wy = kron(diag(Wy),speye(N(2)-1))*dx*dy;
+ if lambda1 > 0
+ Dy = kron(speye(N(1)),spdiags(ones(N(2),1)*[-1 1],[0 1],N(2)-1,N(2))/dy);
+ mat += lambda1*(Dx'*Wx*Dx + Dy'*Wy*Dy);
+ endif
+
+%%% second derivative with respect to x
+ Dxx = spdiags(ones(N(1),1)*[1 -1 -1 1],[-1 0 1 2],N(1)-1,N(1));
+ Dxx(1,1:4) = [3 -7 5 -1]; Dxx(N(1)-1,N(1)-3:N(1)) = [-1 5 -7 3];
+ Dxx = Dxx/(2*dx^2);
+ Dxx = kron(Dxx,speye(N(2)));
+%%% second derivative with respect to y
+ Dyy = spdiags(ones(N(2),1)*[1 -1 -1 1],[-1 0 1 2],N(2)-1,N(2));
+ Dyy(1,1:4) = [3 -7 5 -1]; Dyy(N(2)-1,N(2)-3:N(2)) = [-1 5 -7 3];
+ Dyy = Dyy/(2*dy^2);
+ Dyy = kron(speye(N(1)),Dyy);
+%%% mixed second derivative
+ Dy2 = kron(speye(N(1)-1),spdiags(ones(N(2),1)*[-1 1],[0 1],N(2)-1,N(2))/dy);
+ Dxy = Dy2*Dx*sqrt(dx*dy);
+ mat += lambda2*(Dxx'*Wx*Dxx + Dyy'*Wy*Dyy + 2*Dxy'*Dxy);
+
+%%% solve
+ u = reshape(mat\rhs,N(2),N(1));
+ if nargout>2
+ data_valid =[x,y,Interp*u(:)];
+ endif
+
+endfunction
+
+%!demo
+%! M = 100;
+%! lambda1 = 0; lambda2 = 0.05;
+%! x = 2*rand(M,1)-1; y = 2*rand(M,1)-1;
+%! z = x.*y + 0.1*randn(M,1);
+%! data = [x,y,z];
+%! [grid,u] = regularization2D(data,[-1 1;-1 1],[50 50],lambda1,lambda2);
+%! figure()
+%! mesh(grid.x, grid.y,u)
+%! xlabel('x'); ylabel('y');
+%! hold on
+%! plot3(data(:,1),data(:,2),data(:,3),'*b','Markersize',2)
+%! hold off
+%! view([30,30]);
+
+%!demo
+%! lambda1 = 0; lambda2 = 0.01;
+%! M = 4; angles = [1:M]/M*2*pi;
+%! data = zeros(M+1,3); data(M+1,3) = 1;
+%! data(1:M,1) = cos(angles); data(1:M,2) = sin(angles);
+%! [grid,u] = regularization2D(data,[-1.5 1.5;-1.5 1.5],[50
50],lambda1,lambda2);
+%! figure()
+%! mesh(grid.x, grid.y,u)
+%! xlabel('x'); ylabel('y');
+%! hold on
+%! plot3(data(:,1),data(:,2),data(:,3),'*b','Markersize',2)
+%! hold off
+
+%!test
+%! data = [0,0,0;1,2,3;2,0,2;0,2,2]; % data on z = x+y
+%! lambda1 = 0.0 ; lambda2 = 0.1;
+%! [grid,u,u_valid] = regularization2D(data,[-0.2 2.2;
-0.2,2.2],[11,12],lambda1,lambda2);
+%! assert(norm(data-u_valid),0,1e-12)
+%! assert(norm(grid.x(:) + grid.y(:) - u(:)),0,1e-10)
+
+%!test
+%! data = [0,0,0;1,1,3;2,1,2;0,2,-2;-1.5,-1,0;-2,2,0];
+%! lambda1 = 0.001 ; lambda2 = 0.01;
+%! [grid,u,u_valid] = regularization2D(data,[-2.1
2.1;-2.1,2.1],[40,50],lambda1,lambda2);
+%! value_at_data = [ 1.233351091378741e-01
+%! 2.694454049778803e+00
+%! 2.102786670571836e+00
+%! -1.870655633783656e+00
+%! -1.635243922246946e-02
+%! -3.356775648311422e-02];
+%! assert(norm(value_at_data-u_valid(:,3)),0,1e-12)
