diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..0e020b230f6d06f224446a553d8882298cdcd299
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,156 @@
+# Added Tiz
+tb_logs/
+*.npy
+
+
+# Created by .ignore support plugin (hsz.mobi)
+### Python template
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+env/
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+*.egg-info/
+.installed.cfg
+*.egg
+
+# PyInstaller
+# Usually these files are written by a python script from a template
+# before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*,cover
+.hypothesis/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# IPython Notebook
+.ipynb_checkpoints
+
+# pyenv
+.python-version
+
+# celery beat schedule file
+celerybeat-schedule
+
+# dotenv
+.env
+
+# virtualenv
+venv/
+ENV/
+
+# Spyder project settings
+.spyderproject
+
+# Rope project settings
+.ropeproject
+### VirtualEnv template
+# Virtualenv
+# http://iamzed.com/2009/05/07/a-primer-on-virtualenv/
+.Python
+[Bb]in
+[Ii]nclude
+[Ll]ib
+[Ll]ib64
+[Ll]ocal
+[Ss]cripts
+pyvenv.cfg
+.venv
+pip-selfcheck.json
+### JetBrains template
+# Covers JetBrains IDEs: IntelliJ, RubyMine, PhpStorm, AppCode, PyCharm, CLion, Android Studio and Webstorm
+# Reference: https://intellij-support.jetbrains.com/hc/en-us/articles/206544839
+
+# User-specific stuff:
+.idea/workspace.xml
+.idea/tasks.xml
+.idea/dictionaries
+.idea/vcs.xml
+.idea/jsLibraryMappings.xml
+
+# Sensitive or high-churn files:
+.idea/dataSources.ids
+.idea/dataSources.xml
+.idea/dataSources.local.xml
+.idea/sqlDataSources.xml
+.idea/dynamic.xml
+.idea/uiDesigner.xml
+
+# Gradle:
+.idea/gradle.xml
+.idea/libraries
+
+# Mongo Explorer plugin:
+.idea/mongoSettings.xml
+
+.idea/
+
+## File-based project format:
+*.iws
+
+## Plugin-specific files:
+
+# IntelliJ
+/out/
+
+# mpeltonen/sbt-idea plugin
+.idea_modules/
+
+# JIRA plugin
+atlassian-ide-plugin.xml
+
+# Crashlytics plugin (for Android Studio and IntelliJ)
+com_crashlytics_export_strings.xml
+crashlytics.properties
+crashlytics-build.properties
+fabric.properties
diff --git a/example.py b/example.py
new file mode 100644
index 0000000000000000000000000000000000000000..42062e526c9ee3f26ff8f965426f5423122b776e
--- /dev/null
+++ b/example.py
@@ -0,0 +1,73 @@
+# Example file to demonstrate the use of the VKOGA-PDE code:
+# We use \Omega = [0, 1]^2, L = -LaPlace and the Gaussian kernel.
+
+
+# Some imports
+import numpy as np
+from matplotlib import pyplot as plt
+
+from vkoga_pde.kernels_PDE import Gaussian_laplace
+from vkoga_pde.vkoga_PDE import VKOGA_PDE
+
+np.random.seed(1)
+
+# Create data set: Square domain
+dim = 2
+
+scale_domain = 1
+X1 = scale_domain * np.random.rand(int(1e4), dim)
+X2 = scale_domain * np.random.rand(400, dim)
+X2[:100, 0] = 0
+X2[100:200, 0] = scale_domain
+X2[200:300, 1] = 0
+X2[300:400, 1] = scale_domain
+
+# Define some functions
+u = lambda x: x[:, [1]] * x[:, [0]]**3
+f = lambda x: -6 * x[:, [0]] * x[:, [1]]
+
+# Define right hand side values
+y1 = f(X1)
+y2 = u(X2)
+
+# Initialize and run models
+kernel = Gaussian_laplace(dim=2)
+maxIter = 200
+
+model_f = VKOGA_PDE(kernel=kernel, greedy_type='f_greedy')
+model_p = VKOGA_PDE(kernel=kernel, greedy_type='p_greedy')
+
+_ = model_f.fit(X1, y1, X2, y2, maxIter=maxIter)
+_ = model_p.fit(X1, y1, X2, y2, maxIter=maxIter)
+
+# Plot some results
+plt.figure(1)
+plt.clf()
+plt.plot(np.array(model_f.train_hist['f'])**.5, 'r')
+plt.plot(np.array(model_p.train_hist['f'])**.5, 'b')
+plt.yscale('log')
+plt.xscale('log')
+plt.legend(['f-greedy', 'p-greedy'])
+plt.title('residuals')
+plt.draw()
+
+plt.figure(2)
+plt.clf()
+plt.plot(model_f.ctrs_[model_f.ind_ctrs1, 0], model_f.ctrs_[model_f.ind_ctrs1, 1], 'rx')
+plt.plot(model_f.ctrs_[model_f.ind_ctrs2, 0], model_f.ctrs_[model_f.ind_ctrs2, 1], 'ro')
+plt.title('Selected centers of f-greedy')
+plt.draw()
+
+plt.figure(3)
+plt.clf()
+plt.plot(model_p.ctrs_[model_p.ind_ctrs1, 0], model_p.ctrs_[model_p.ind_ctrs1, 1], 'bx')
+plt.plot(model_p.ctrs_[model_p.ind_ctrs2, 0], model_p.ctrs_[model_p.ind_ctrs2, 1], 'bo')
+plt.title('Selected centers of P-greedy')
+plt.draw()
+
+plt.figure(4)
+plt.clf()
+plt.plot(X1[:, 0], X1[:, 1], 'm.')
+plt.plot(X2[:, 0], X2[:, 1], 'r.')
+
+
diff --git a/kernels_PDE.py b/kernels_PDE.py
new file mode 100644
index 0000000000000000000000000000000000000000..5caac657bbc141b1f45ff04f796fae119c2f4b3b
--- /dev/null
+++ b/kernels_PDE.py
@@ -0,0 +1,237 @@
+#!/usr/bin/env python3
+
+from abc import ABC, abstractmethod
+from scipy.spatial import distance_matrix
+import numpy as np
+
+
+# Abstract kernel
+class Kernel(ABC):
+ @abstractmethod
+ def __init__(self):
+ super().__init__()
+
+ @abstractmethod
+ def eval(self, x, y):
+ pass
+
+ def eval_prod(self, x, y, v, batch_size=100):
+ N = x.shape[0]
+ num_batches = int(np.ceil(N / batch_size))
+ mat_vec_prod = np.zeros((N, 1))
+ for idx in range(num_batches):
+ idx_begin = idx * batch_size
+ idx_end = (idx + 1) * batch_size
+ A = self.eval(x[idx_begin:idx_end, :], y)
+ mat_vec_prod[idx_begin:idx_end] = A @ v
+ return mat_vec_prod
+
+ @abstractmethod
+ def diagonal(self, X):
+ pass
+
+ @abstractmethod
+ def __str__(self):
+ pass
+
+ @abstractmethod
+ def set_params(self, params):
+ pass
+
+
+# Abstract RBF
+class L_RBF(Kernel):
+ @abstractmethod
+ def __init__(self):
+ super(L_RBF, self).__init__()
+
+ def eval(self, x, y):
+ return self.rbf(self.ep, distance_matrix(np.atleast_2d(x), np.atleast_2d(y)))
+
+ def diagonal(self, X):
+ return np.ones(X.shape[0]) * self.rbf(self.ep, 0)
+
+ def __str__(self):
+ return self.name + ' [gamma = %2.2e]' % self.ep
+
+ def set_params(self, par):
+ self.ep = par
+
+ def dd_diagonal(self, X):
+ # Evaluation of (L^x - shift) (L^y - shift) k(x,y) for x=y: Likely some constant
+
+ return np.ones(X.shape[0]) * self.dd_eval(np.zeros((1, 1)), np.zeros((1, 1))).item()
+
+ def dd_eval(self, x, y):
+ # Evaluation of (L^x - shift) (L^y - shift) k(x,y)
+
+ array_dists = distance_matrix(np.atleast_2d(x), np.atleast_2d(y))
+
+ return self.rbf_utility(self.ep, array_dists) \
+ - 2 * self.shift * self.d2_eval(x, y) - self.shift ** 2 * self.rbf(self.ep, array_dists)
+
+ def d2_eval(self, x, y):
+ # Evaluation of (L^z - shift * Id) k(*,z)|_{x,y}
+ # Required exactly like this for "compute the nth basis", see also 14.12.21\A1
+
+ array_dists = distance_matrix(np.atleast_2d(x), np.atleast_2d(y))
+
+ array_result = - ((self.dim - 1) * self.rbf1_divided_by_r(self.ep, array_dists)
+ + self.rbf2(self.ep, array_dists)) \
+ - self.shift * self.rbf(self.ep, array_dists)
+
+ return array_result
+
+ def d1_eval(self, x, y):
+ # Evaluation of (L - shift * Id) k(*,y)|_{x}
+ # Required exactly like this for "compute the nth basis", see also 14.12.21\A1
+
+ # Since the LaPlace is radial, the formula is the same as d2_eval
+ return self.d2_eval(x, y)
+
+ def mixed_eval_L(self, x, ctrs, list_ind_interior):
+ # Kernel evaluation L^x k(x, y), whereby y1 := y[ind_interior, :] are centers in the interior,
+ # y2 := y[~ind_interior, :] are centers on the boundary.
+ # Therefore ddeval needs to be used with respect to y1, but deval with respect to y2.
+
+ list_not_ind_interior = [idx for idx in range(x.shape[0]) if idx not in list_ind_interior]
+
+ array_result = np.zeros((x.shape[0], ctrs.shape[0]))
+
+ array_result[list_ind_interior, :] = self.dd_eval(x[list_ind_interior, :], ctrs)
+ array_result[list_not_ind_interior, :] = self.d1_eval(x[list_not_ind_interior, :], ctrs)
+
+ return array_result
+
+ def mixed_eval(self, x, ctrs, list_ind_interior):
+ # Kernel evaluation k(x, y), whereby y1 := y[ind_interior, :] are centers in the interior,
+ # y2 := y[~ind_interior, :] are centers on the boundary.
+ # Therefore deval needs to be used with respect to y1, but eval with respect to y2.
+
+ list_not_ind_interior = [idx for idx in range(x.shape[0]) if idx not in list_ind_interior]
+
+ array_result = np.zeros((x.shape[0], ctrs.shape[0]))
+
+ array_result[list_ind_interior, :] = self.d2_eval(x[list_ind_interior, :], ctrs)
+ array_result[list_not_ind_interior, :] = self.eval(x[list_not_ind_interior, :], ctrs)
+
+ return array_result
+
+
+# Implementation of concrete RBFs
+class Gaussian_laplace(L_RBF):
+ # # Gaussian kernel phi(r) = exp(-r^2),
+ # used in conjunction with L = -LaPlace - shift * Id, whereby shift is some constant
+ # (e.g. an eigenvalue of -LaPlace).
+
+ def __init__(self, dim, shift=0, ep=1):
+ super().__init__()
+ self.dim = dim
+ self.shift = shift
+ self.ep = 1 # ToDo: Generalize this one!!!
+ self.name = 'gauss'
+
+ # Implement not only the rbf, but also its derivatives
+ self.rbf = lambda ep, r: np.exp(-r**2)
+
+ self.rbf1 = lambda ep, r: -2 * r * np.exp(-r**2)
+ self.rbf1_divided_by_r = lambda ep, r: -2 * np.exp(-r**2)
+
+ self.rbf2 = lambda ep, r: (4 * r ** 2 - 2) * np.exp(-r**2)
+
+ # Utility RBF for dd_eval: Double LaPlacian
+ self.rbf_utility = lambda ep, r: \
+ np.exp(-r ** 2) * (16 * r ** 4 - (16 * self.dim + 32) * r ** 2 + (4 * self.dim ** 2 + 8 * self.dim))
+
+
+class wendland32_laplace(L_RBF):
+ # Wendland kernel 3, 2: phi(r) = (1-r)_+^6 * (35r**2 + 18r + 3)
+ # used in conjunction with L = -LaPlace - shift * Id, whereby shift is some constant
+ # (e.g. an eigenvalue of -LaPlace).
+
+ # Only for dim <= 3!
+
+ def __init__(self, dim, shift=0, ep=1):
+ super().__init__()
+ self.dim = dim
+ self.shift = shift
+ self.ep = 1 # ep ToDo: Implement general case!
+ self.name = 'Wendland 3, 2'
+
+ # Implement not only the rbf, but also its derivatives
+ self.rbf = lambda ep, r: 1/3 * np.maximum(1-r, 0)**6 * (35 * r**2 + 18 * r + 3)
+
+ self.rbf1 = lambda ep, r: -1/3 * np.maximum(1 - r, 0) ** 5 * 56 * r * (5 * r + 1)
+ self.rbf1_divided_by_r = lambda ep, r: -1/3 * np.maximum(1 - r, 0) ** 5 * 56 * (5 * r + 1)
+
+ self.rbf2 = lambda ep, r: 1/3 * np.maximum(1 - r, 0) ** 4 * 56 * (35 * r ** 2 - 4 * r - 1)
+
+ # Utility RBF for dd_eval: Double LaPlacian
+ self.rbf_utility = lambda ep, r: \
+ 1/3 * (1680*self.dim**2 + 3360*self.dim
+ + r**4*(1680*self.dim**2 + 16800*self.dim + 40320)
+ + r**3*(-6720*self.dim**2 - 53760*self.dim - 100800)
+ + r**2*(10080*self.dim**2 + 60480*self.dim + 80640)
+ + r*(-6720*self.dim**2 - 26880*self.dim - 20160)) * (r < 1)
+
+
+class wendland33_laplace(L_RBF):
+ # Wendland kernel 3, 3: phi(r) = (1-r)_+^8 * (32r**3+ 25r^2 + 8r + 1)
+ # used in conjunction with L = -LaPlace - shift * Id, whereby shift is some constant
+ # (e.g. an eigenvalue of -LaPlace).
+
+ # Only for dim <= 3!
+
+ def __init__(self, dim, shift=0, ep=1):
+ super().__init__()
+ self.dim = dim
+ self.shift = shift
+ self.ep = 1 # ep ToDo: Implement general case!
+ self.name = 'wendland 3, 3'
+
+ # Implement not only the rbf, but also its derivatives
+ self.rbf = lambda ep, r: np.maximum(1 - r, 0) ** 8 * (32 * r ** 3 + 25 * r ** 2 + 8 * r + 1)
+
+ self.rbf1 = lambda ep, r: - np.maximum(1 - r, 0) ** 7 * 22 * r * (16 * r ** 2 + 7 * r + 1)
+ self.rbf1_divided_by_r = lambda ep, r: - np.maximum(1 - r, 0) ** 7 * 22 * (16 * r ** 2 + 7 * r + 1)
+
+ self.rbf2 = lambda ep, r: np.maximum(1 - r, 0) ** 6 * 22 * (160 * r ** 3 + 15 * r ** 2 - 6 * r - 1)
+
+ # Utility RBF for dd_eval: Double LaPlacian
+ self.rbf_utility = lambda ep, r: \
+ (528 * self.dim ** 2 + 1056 * self.dim + r ** 7 * (3168 * self.dim ** 2 + 50688 * self.dim + 199584)
+ + r ** 6 * (-18480 * self.dim ** 2 - 258720 * self.dim - 887040)
+ + r ** 5 * (44352 * self.dim ** 2 + 532224 * self.dim + 1552320)
+ + r ** 4 * (-55440 * self.dim ** 2 - 554400 * self.dim - 1330560)
+ + r ** 3 * (36960 * self.dim ** 2 + 295680 * self.dim + 554400)
+ + r ** 2 * (-11088 * self.dim ** 2 - 66528 * self.dim - 88704)) * (r < 1)
+
+
+
+
+class quadrMatern_laplace(L_RBF):
+ # Quadratic matern kernel phi(r) = (3 + 3r + r^2)*exp(-r),
+ # used in conjunction with L = -LaPlace
+
+ def __init__(self, dim, shift=0, ep=1):
+ super().__init__()
+ self.dim = dim
+ self.shift = shift
+ self.ep = 1 # ep ToDo: Implement general case!
+ self.name = 'quadratic Matern'
+
+ # Implement not only the rbf, but also its derivatives
+ self.rbf = lambda ep, r: np.exp(-r) * (3 + 3 *r + r ** 2)
+
+ self.rbf1 = lambda ep, r: np.exp(-r) * (-r) * (r+1)
+ self.rbf1_divided_by_r = lambda ep, r: np.exp(-r) * (-1) * (r+1)
+
+ self.rbf2 = lambda ep, r: np.exp(-r) * (r ** 2 - r - 1)
+
+ # Utility RBF for dd_eval: Double LaPlacian
+ self.rbf_utility = lambda ep, r: \
+ np.exp(-r) * (r**2 - (2 * self.dim + 3) * r + self.dim ** 2 + 2 * self.dim)
+
+
+
+
diff --git a/vkoga_PDE.py b/vkoga_PDE.py
new file mode 100644
index 0000000000000000000000000000000000000000..534434682dacf791e1e9af8ef5611810d242dac8
--- /dev/null
+++ b/vkoga_PDE.py
@@ -0,0 +1,594 @@
+# First approach to develope code for greedy PDE.
+# ToDo: Not sure whether vector valued output makes sense, but let's try to implement it ..
+
+from vkoga_pde.kernels_PDE import wendland32_laplace
+import numpy as np
+from sklearn.base import BaseEstimator
+from sklearn.utils.validation import check_X_y, check_array
+
+
+# VKOGA implementation
+class VKOGA_PDE(BaseEstimator):
+
+ def __init__(self, kernel=wendland32_laplace(dim=2), kernel_par=1,
+ verbose=True, n_report=10,
+ greedy_type='p_greedy', reg_par=0, restr_par=0,
+ tol_f=1e-10, tol_p=1e-10, beta=None):
+
+ # Set the verbosity on/off
+ self.verbose = verbose
+
+ # Set the frequency of report
+ self.n_report = n_report
+
+ # Set the params defining the method
+ self.kernel = kernel
+ self.kernel_par = kernel_par
+ self.reg_par = reg_par
+ self.restr_par = restr_par
+
+ # Define the greedy method properly
+ if beta is None:
+ self.greedy_type = greedy_type
+ if self.greedy_type == 'f_greedy':
+ self.beta = 1.0
+ elif self.greedy_type == 'p_greedy':
+ self.beta = 0.0
+ else:
+ self.beta = float(beta)
+ if self.beta == 1.0:
+ self.greedy_type = 'f_greedy'
+ elif self.beta == 0.0:
+ self.greedy_type = 'p_greedy'
+ else:
+ self.greedy_type = 'beta'
+
+ # ToDo: Did never think about the regularization parameter ..
+ assert reg_par == 0, 'reg_par > 0 might not be implemented correctly!'
+
+ self.flag_val = None
+
+ # Set the stopping values
+ self.tol_f = tol_f
+ self.tol_p = tol_p
+
+ # Some further settings
+ self.ctrs_ = None # centers for the interior
+ self.ind_ctrs1 = [] # self.ctrs1_ = self.ctrs_[self.indI_ctrs1], i.e. interior
+ self.ind_ctrs2 = [] # self.ctrs2_ = self.ctrs_[self.indI_ctrs2], i.e. boundary
+
+ self.Cut_ = None
+ self.c = None
+
+ # Initialize the convergence history
+ self.train_hist = {}
+ self.train_hist['n'] = []
+ self.train_hist['f'] = []
+ self.train_hist['f1'] = []
+ self.train_hist['f2'] = []
+ self.train_hist['p'] = []
+ self.train_hist['p1'] = []
+ self.train_hist['p2'] = []
+ self.train_hist['p selected'] = [] # list of selected power vals
+ self.train_hist['f val'] = []
+ self.train_hist['p val'] = []
+
+ # Initialize variables for the maximal kernel diagonal value in the interior and the boundary
+ # ToDo: So far only f-greedy and P-greedy is implemented in a meaningful way,
+ # for all other \beta values the weighting of the power function (via the kernel
+ # diagonal values) is not implemented.
+ self.max_kernel_diag1 = None
+ self.max_kernel_diag2 = None
+
+ def selection_rule(self, f1, f2, p1, p2):
+ # Generalization of selection rule: We have residual values and power function values
+ # for both the interior (f1, p1) as well as for the boundary (p1, p2).
+
+ # ToDo: So far only f-greedy and P-greedy is implemented in a meaningful way,
+ # for all other \beta values the weighting of the power function (via the kernel
+ # diagonal values) is not implemented.
+
+ if self.restr_par > 0:
+ # ToDo: Not sure whether this will really work, restr_idx might be empty
+
+ p_ = max([np.max(p1), np.max(p2)])
+ restr_idx1 = np.nonzero(p1 >= self.restr_par * p_)[0]
+ restr_idx2 = np.nonzero(p2 >= self.restr_par * p_)[0]
+ else:
+ restr_idx1 = np.arange(len(p1))
+ restr_idx2 = np.arange(len(p2))
+
+ f1 = np.sum(f1 ** 2, axis=1)
+ f2 = np.sum(f2 ** 2, axis=1)
+
+ # Case distinction required for fp_greedy
+ if self.greedy_type == 'fp_greedy':
+ array_weighted1 = f1[restr_idx1] / p1[restr_idx1]
+ idx1 = np.argmax(array_weighted1)
+
+ if len(f2[restr_idx2]) != 0: # required especially in 1D
+ array_weighted2 = f2[restr_idx2] / p2[restr_idx2]
+ idx2 = np.argmax(array_weighted2)
+
+ if array_weighted1[idx1] > array_weighted2[idx2]:
+ one_or_two = 1
+ idx = restr_idx1[idx1]
+ else:
+ one_or_two = 2
+ idx = restr_idx2[idx2]
+ else:
+ one_or_two = 1
+ idx = restr_idx1[idx1]
+ elif self.greedy_type == 'p_greedy':
+ # extra implementation for P-greedy to address weighting of powerfunc (interior, boundary)
+ array_weighted1 = p1[restr_idx1]
+ idx1 = np.argmax(array_weighted1)
+
+ if len(f2[restr_idx2]) != 0: # required especially in 1D
+ array_weighted2 = p2[restr_idx2]
+ idx2 = np.argmax(array_weighted2)
+
+ # Pay attention at the weighting of the power functions
+ if array_weighted1[idx1] / self.max_kernel_diag1 \
+ > array_weighted2[idx2] / self.max_kernel_diag2:
+ one_or_two = 1
+ idx = restr_idx1[idx1]
+ else:
+ one_or_two = 2
+ idx = restr_idx2[idx2]
+ else:
+ one_or_two = 1
+ idx = restr_idx1[idx1]
+
+ else:
+ array_weighted1 = f1[restr_idx1] ** self.beta * p1[restr_idx1] ** (1 - self.beta)
+ idx1 = np.argmax(array_weighted1)
+
+ if len(f2[restr_idx2]) != 0: # required especially in 1D
+ array_weighted2 = f2[restr_idx2] ** self.beta * p2[restr_idx2] ** (1 - self.beta)
+ idx2 = np.argmax(array_weighted2)
+
+ if array_weighted1[idx1] > array_weighted2[idx2]:
+ one_or_two = 1
+ idx = restr_idx1[idx1]
+ else:
+ one_or_two = 2
+ idx = restr_idx2[idx2]
+ else:
+ one_or_two = 1
+ idx = restr_idx1[idx1]
+
+ # ToDo: This implementation is quite laborious ..
+ if self.greedy_type == 'p_greedy':
+ f_max1, p_max1 = np.max(f1), np.max(p1) / self.max_kernel_diag1
+ else:
+ f_max1, p_max1 = np.max(f1), np.max(p1)
+
+ if len(f2[restr_idx2]) != 0:
+ if self.greedy_type == 'p_greedy':
+ f_max2, p_max2 = np.max(f2), np.max(p2) / self.max_kernel_diag2
+ else:
+ f_max2, p_max2 = np.max(f2), np.max(p2)
+ else:
+ f_max2, p_max2 = 0, 0
+
+ return (one_or_two, idx), f_max1, f_max2, p_max1, p_max2
+
+ def fit(self, X1, y1, X2, y2,
+ X1_val=None, y1_val=None, X2_val=None, y2_val=None, maxIter=None):
+ # Points in X1 are considered as interior points, points in X2 as boundary points.
+ # They can coincide (partly), but for X1 we use the L-evaluation, for X2 we use the point evaluation.
+
+ # Compute the maximal kernel diagonal value in the interior and the boundary
+ self.max_kernel_diag2 = np.max(self.kernel.diagonal(X2))
+ self.max_kernel_diag1 = np.max(self.kernel.dd_diagonal(X1))
+
+ list_ctrs = [] # not that beautiful, but required to store the final centers
+
+ # Check the datasets
+ X1, y1 = check_X_y(X1, y1, multi_output=True)
+ X2, y2 = check_X_y(X2, y2, multi_output=True) # ToDo: Check for consistency of X1, X2, X1_val, X2_val
+
+ if len(y1.shape) == 1:
+ y1 = y1[:, None]
+ if len(y2.shape) == 1:
+ y2 = y2[:, None]
+
+ # Check the validation datasets (if provided, otherwise set to None)
+ if X1_val is None or y1_val is None or X2_val is None or y2_val is None:
+ X1_val = None
+ y1_val = None
+ X2_val = None
+ y2_val = None
+ self.flag_val = False
+ else:
+ self.flag_val = True
+ X1_val, y1_val = check_X_y(X1_val, y2_val, multi_output=True)
+ X2_val, y2_val = check_X_y(X2_val, y2_val, multi_output=True)
+ # We will assume in the following that X_val and X are disjoint
+
+ if len(y1_val.shape) == 1:
+ y1_val = y1_val[:, None]
+ if len(y2_val.shape) == 1:
+ y2_val = y2_val[:, None]
+
+ # Check whether already fitted
+ if self.ctrs_ is None:
+ self.ctrs_ = np.zeros((0, X1.shape[1]))
+ self.Cut_ = np.zeros((0, 0))
+ self.c = np.zeros((0, y1.shape[1]))
+
+ # ToDo: If self.ctrs_ is not None, does this work?
+
+ # Check whether "new X1" and previously chosen centers overlap
+ list_truly_new_X1 = []
+ if self.ctrs_.shape[0] > 0:
+ for idx_x in range(X1.shape[0]):
+ if min(np.linalg.norm(self.ctrs_ - X1[idx_x, :], axis=1)) < 1e-10:
+ continue
+ else:
+ list_truly_new_X1.append(idx_x)
+ else:
+ list_truly_new_X1 = list(range(X1.shape[0]))
+ X1 = X1[list_truly_new_X1, :]
+ y1 = y1[list_truly_new_X1, :]
+
+ # Check whether "new X2" and previously chosen centers overlap
+ list_truly_new_X2 = []
+ if self.ctrs_.shape[0] > 0:
+ for idx_x in range(X2.shape[0]):
+ if min(np.linalg.norm(self.ctrs_ - X2[idx_x, :], axis=1)) < 1e-10:
+ continue
+ else:
+ list_truly_new_X2.append(idx_x)
+ else:
+ list_truly_new_X2 = list(range(X2.shape[0]))
+
+ X2 = X2[list_truly_new_X2, :]
+ y2 = y2[list_truly_new_X2, :]
+
+ # Initialize the residual and update the given y values by substracting the current model
+ y1 = y1 - self.predict_Ls(X1)
+ y2 = y2 - self.predict_s(X2)
+
+ if self.flag_val:
+ y1_val = y1_val - self.predict_Ls(X1_val)
+ y2_val = y2_val - self.predict_s(X2_val)
+
+ # Get the data dimension
+ N, q = y1.shape[0] + y2.shape[0], y1.shape[1]
+ N1, N2 = y1.shape[0], y2.shape[0]
+ if self.flag_val:
+ N_val = y1_val.shape[0] + y2_val.shape[0]
+ N1_val, N2_val = y1_val.shape[0], y2_val.shape[0]
+
+ # Set maxIter_continue
+ if maxIter is None or maxIter > N:
+ self.maxIter = min([N, 100])
+ else:
+ self.maxIter = maxIter
+
+ # Check compatibility of restriction
+ if self.beta == 0.0: # in case of P-greedy
+ self.restr_par = 0
+ if not self.reg_par == 0:
+ self.restr_par = 0
+
+ # Initialize list for the chosen and non-chosen indices - seperately for interior and boundary
+ indI1_, indI2_ = [], []
+ notIndI1, notIndI2 = list(range(y1.shape[0])), list(range(y2.shape[0]))
+ c = np.zeros((self.maxIter, q))
+
+ # Compute the Newton basis values (related to the old centers) on the new X
+ # ToDo: I think restart training etc does not work properly, since it is not debugged propely!!
+ # Thus the following lines of codes may be irrelevant
+ Vx_new_X1_old_ctrs = self.kernel.mixed_eval(X1, self.ctrs_, self.ind_ctrs1) @ self.Cut_.transpose()
+ Vx_new_X2_old_ctrs = self.kernel.mixed_eval(X2, self.ctrs_, self.ind_ctrs1) @ self.Cut_.transpose()
+ if self.flag_val:
+ Vx_val_new_X1_old_ctrs = self.kernel.mixed_eval(
+ X1_val, self.ctrs_, self.ind_ctrs1) @ self.Cut_.transpose()
+ Vx_val_new_X2_old_ctrs = self.kernel.mixed_eval(
+ X2_val, self.ctrs_, self.ind_ctrs1) @ self.Cut_.transpose()
+
+ # Compute the application of L to the Newton basis values (related to the old centers) on the new X
+ LVx_new_X1_old_ctrs = self.kernel.mixed_eval_L(X1, self.ctrs_, self.ind_ctrs1) @ self.Cut_.transpose()
+ LVx_new_X2_old_ctrs = self.kernel.mixed_eval_L(X2, self.ctrs_, self.ind_ctrs1) @ self.Cut_.transpose()
+ if self.flag_val:
+ LVx_val_new_X1_old_ctrs = self.kernel.mixed_eval_L(
+ X1_val, self.ctrs_, self.ind_ctrs1) @ self.Cut_.transpose()
+ LVx_val_new_X2_old_ctrs = self.kernel.mixed_eval_L(
+ X2_val, self.ctrs_, self.ind_ctrs1) @ self.Cut_.transpose()
+
+ # Initialize arrays for the (L-evaluation of the) Newton basis values ..
+ # .. (related to the new centers) on the new X1, X2
+ Vx1 = np.zeros((N1, self.maxIter))
+ Vx2 = np.zeros((N2, self.maxIter))
+ LVx1 = np.zeros((N1, self.maxIter))
+ LVx2 = np.zeros((N2, self.maxIter))
+
+ if self.flag_val:
+ Vx1_val = np.zeros((N1_val, self.maxIter))
+ Vx2_val = np.zeros((N2_val, self.maxIter))
+ LVx1_val = np.zeros((N1_val, self.maxIter))
+ LVx2_val = np.zeros((N2_val, self.maxIter))
+
+ # Compute the powervals on X1, X2 and X1_val, X2_val
+ p1 = self.kernel.dd_diagonal(X1) + self.reg_par
+ p1 = p1 - np.sum(LVx_new_X1_old_ctrs ** 2, axis=1)
+
+ p2 = self.kernel.diagonal(X2) + self.reg_par
+ p2 = p2 - np.sum(Vx_new_X2_old_ctrs ** 2, axis=1)
+
+ if self.flag_val:
+ p1_val = self.kernel.dd_diagonal(X1_val) + self.reg_par
+ p1_val = p1_val - np.sum(LVx_val_new_X1_old_ctrs ** 2, axis=1)
+
+ p2_val = self.kernel.diagonal(X2_val) + self.reg_par
+ p2_val = p2_val - np.sum(Vx_val_new_X2_old_ctrs ** 2, axis=1)
+
+ # Extend Cut_ matrix, i.e. continue to build on old self.Cut_ matrix
+ N_ctrs_so_far, N_ctrs_so_far1, N_ctrs_so_far2 = \
+ self.Cut_.shape[0], len(self.ind_ctrs1), len(self.ind_ctrs2)
+
+ Cut_ = np.zeros((N_ctrs_so_far + self.maxIter, N_ctrs_so_far + self.maxIter))
+ Cut_[:N_ctrs_so_far, :N_ctrs_so_far] = self.Cut_
+
+ # Iterative selection of new points
+ self.print_message('begin')
+
+ n1, n2 = 0, 0 # number of centers in interior or boundary
+ for n in range(self.maxIter):
+ # prepare
+ self.train_hist['n'].append(self.ctrs_.shape[0] + n + 1)
+ self.train_hist['f'].append([])
+ self.train_hist['f1'].append([])
+ self.train_hist['f2'].append([])
+ self.train_hist['p'].append([])
+ self.train_hist['p1'].append([])
+ self.train_hist['p2'].append([])
+ self.train_hist['p selected'].append([])
+ if self.flag_val:
+ self.train_hist['p val'].append([])
+ self.train_hist['f val'].append([])
+
+ # select the current index
+ (one_or_two, idx), self.train_hist['f1'][-1], self.train_hist['f2'][-1], \
+ self.train_hist['p1'][-1], self.train_hist['p2'][-1] \
+ = self.selection_rule(y1[notIndI1], y2[notIndI2], p1[notIndI1], p2[notIndI2])
+ self.train_hist['f'][-1] = max(self.train_hist['f1'][-1], self.train_hist['f2'][-1])
+ self.train_hist['p'][-1] = max(self.train_hist['p1'][-1], self.train_hist['p2'][-1])
+
+ # add the current index and store selected p value and store the selected center
+ if one_or_two == 1:
+ self.train_hist['p selected'][-1] = p1[notIndI1[idx]]
+ indI1_.append(notIndI1[idx])
+
+ self.ind_ctrs1.append(N_ctrs_so_far + n)
+ list_ctrs.append(X1[[notIndI1[idx]], :])
+
+ elif one_or_two == 2:
+ self.train_hist['p selected'][-1] = p2[notIndI2[idx]]
+ indI2_.append(notIndI2[idx])
+
+ self.ind_ctrs2.append(N_ctrs_so_far + n)
+ list_ctrs.append(X2[[notIndI2[idx]], :])
+
+ if self.flag_val:
+ self.train_hist['p val'][-1] = max([np.max(p1_val), np.max(p2_val)])
+ self.train_hist['f val'][-1] = max(
+ [np.max(np.sum(y1_val ** 2, axis=1)), np.max(np.sum(y2_val ** 2, axis=1))])
+
+ # check if the tolerances are reacheded
+ if self.train_hist['f'][n] <= self.tol_f or self.train_hist['p'][n] <= self.tol_p:
+ n = n - 1
+ if one_or_two == 1:
+ n1 = n1 - 1
+ elif one_or_two == 2:
+ n2 = n2 - 1
+ self.print_message('end')
+ break
+
+ # compute the nth basis (including normalization), see also 14.12.21\A1
+ if one_or_two == 1:
+ # Computations (this is simply expansion in Newton basis) ..
+ Vx1[notIndI1, n] = self.kernel.d2_eval(X1[notIndI1, :], X1[indI1_[n1], :])[:, 0] \
+ - Vx_new_X1_old_ctrs[notIndI1, :] @ LVx_new_X1_old_ctrs[indI1_[n1], :].transpose() \
+ - Vx1[notIndI1, :] @ LVx1[indI1_[n1], :].transpose()
+ Vx2[notIndI2, n] = self.kernel.d2_eval(X2[notIndI2, :], X1[indI1_[n1], :])[:, 0] \
+ - Vx_new_X2_old_ctrs[notIndI2, :] @ LVx_new_X1_old_ctrs[indI1_[n1], :].transpose() \
+ - Vx2[notIndI2, :] @ LVx1[indI1_[n1], :].transpose()
+
+ LVx1[notIndI1, n] = self.kernel.dd_eval(X1[notIndI1, :], X1[indI1_[n1], :])[:, 0] \
+ - LVx_new_X1_old_ctrs[notIndI1, :] @ LVx_new_X1_old_ctrs[indI1_[n1], :].transpose() \
+ - LVx1[notIndI1, :] @ LVx1[indI1_[n1], :].transpose()
+ LVx2[notIndI2, n] = self.kernel.dd_eval(X2[notIndI2, :], X1[indI1_[n1], :])[:, 0] \
+ - LVx_new_X2_old_ctrs[notIndI2, :] @ LVx_new_X1_old_ctrs[indI1_[n1], :].transpose() \
+ - LVx2[notIndI2, :] @ LVx1[indI1_[n1], :].transpose()
+
+ Vx1[indI1_[n1], n] += self.reg_par # ToDo: This is likely wrong/incomplete/..
+
+ # Normalizations ..
+ Vx1[notIndI1, n] = Vx1[notIndI1, n] / np.sqrt(p1[indI1_[n1]])
+ Vx2[notIndI2, n] = Vx2[notIndI2, n] / np.sqrt(p1[indI1_[n1]])
+ LVx1[notIndI1, n] = LVx1[notIndI1, n] / np.sqrt(p1[indI1_[n1]])
+ LVx2[notIndI2, n] = LVx2[notIndI2, n] / np.sqrt(p1[indI1_[n1]])
+
+ # ToDo: Take care of the validation sets ..
+ # if self.flag_val:
+ # Vx_val[:, n] = self.kernel.eval(X_val, X[indI_[n], :])[:, 0]\
+ # - Vx_val_new_X_old_ctrs[:, :] @ Vx_new_X_old_ctrs[indI_[n], :].transpose()\
+ # - Vx_val[:, :n+1] @ Vx[indI_[n], 0:n+1].transpose()
+ # Vx_val[:, n] = Vx_val[:, n] / np.sqrt(p[indI_[n]])
+
+ elif one_or_two == 2:
+ # Computations (this is simply expansion in Newton basis) ..
+ Vx1[notIndI1, n] = self.kernel.eval(X1[notIndI1, :], X2[indI2_[n2], :])[:, 0] \
+ - Vx_new_X1_old_ctrs[notIndI1, :] @ Vx_new_X2_old_ctrs[indI2_[n2], :].transpose() \
+ - Vx1[notIndI1, :] @ Vx2[indI2_[n2], :].transpose()
+ Vx2[notIndI2, n] = self.kernel.eval(X2[notIndI2, :], X2[indI2_[n2], :])[:, 0] \
+ - Vx_new_X2_old_ctrs[notIndI2, :] @ Vx_new_X2_old_ctrs[indI2_[n2], :].transpose() \
+ - Vx2[notIndI2, :] @ Vx2[indI2_[n2], :].transpose()
+
+ LVx1[notIndI1, n] = self.kernel.d1_eval(X1[notIndI1, :], X2[indI2_[n2], :])[:, 0] \
+ - LVx_new_X1_old_ctrs[notIndI1, :] @ Vx_new_X2_old_ctrs[indI2_[n2], :].transpose() \
+ - LVx1[notIndI1, :] @ Vx2[indI2_[n2], :].transpose()
+ LVx2[notIndI2, n] = self.kernel.d1_eval(X2[notIndI2, :], X2[indI2_[n2], :])[:, 0] \
+ - LVx_new_X2_old_ctrs[notIndI2, :] @ Vx_new_X2_old_ctrs[indI2_[n2], :].transpose() \
+ - LVx2[notIndI2, :] @ Vx2[indI2_[n2], :].transpose()
+
+ Vx2[indI2_[n2], n] += self.reg_par # ToDo: This is likely wrong/incomplete/..
+
+ # Normalizations ..
+ Vx1[notIndI1, n] = Vx1[notIndI1, n] / np.sqrt(p2[indI2_[n2]])
+ Vx2[notIndI2, n] = Vx2[notIndI2, n] / np.sqrt(p2[indI2_[n2]])
+ LVx1[notIndI1, n] = LVx1[notIndI1, n] / np.sqrt(p2[indI2_[n2]])
+ LVx2[notIndI2, n] = LVx2[notIndI2, n] / np.sqrt(p2[indI2_[n2]])
+
+ # ToDo: Take care of the validation sets
+
+ # update the change of basis (see 14.12.21\B)
+ Cut_new_row = np.ones(N_ctrs_so_far + n + 1) # required only for the final one
+
+ if one_or_two == 1:
+ Cut_new_row[:N_ctrs_so_far + n] = \
+ -np.concatenate((LVx_new_X1_old_ctrs[indI1_[n1], :], LVx1[indI1_[n1], :n])) \
+ @ Cut_[:N_ctrs_so_far + n, :N_ctrs_so_far + n]
+
+ Cut_[N_ctrs_so_far + n, :N_ctrs_so_far + n + 1] = Cut_new_row / LVx1[indI1_[n1], n]
+ elif one_or_two == 2:
+ Cut_new_row[:N_ctrs_so_far + n] = \
+ -np.concatenate((Vx_new_X2_old_ctrs[indI2_[n2], :], Vx2[indI2_[n2], :n])) \
+ @ Cut_[:N_ctrs_so_far + n, :N_ctrs_so_far + n]
+
+ Cut_[N_ctrs_so_far + n, :N_ctrs_so_far + n + 1] = Cut_new_row / Vx2[indI2_[n2], n]
+
+ # compute the nth coefficient
+ if one_or_two == 1:
+ c[n] = y1[indI1_[n1]] / np.sqrt(p1[indI1_[n1]])
+ elif one_or_two == 2:
+ c[n] = y2[indI2_[n2]] / np.sqrt(p2[indI2_[n2]])
+
+ # update the power functions: ToDo: For Gaussian in 1D we obtain negative power vals :(
+ p1[notIndI1] = p1[notIndI1] - LVx1[notIndI1, n] ** 2
+ p2[notIndI2] = p2[notIndI2] - Vx2[notIndI2, n] ** 2
+
+ if self.flag_val:
+ p1_val = p1_val - LVx1_val[:, n] ** 2
+ p2_val = p2_val - Vx2_val[:, n] ** 2
+
+ # update the residual
+ y1[notIndI1] = y1[notIndI1] - LVx1[notIndI1, n][:, None] * c[n]
+ y2[notIndI2] = y2[notIndI2] - Vx2[notIndI2, n][:, None] * c[n]
+
+ if self.flag_val:
+ y1_val = y1_val - LVx1_val[:, n][:, None] * c[n]
+ y2_val = y2_val - Vx2_val[:, n][:, None] * c[n]
+
+ # remove the nth index from the dictionary
+ if one_or_two == 1:
+ notIndI1.pop(idx)
+ n1 += 1
+ elif one_or_two == 2:
+ notIndI2.pop(idx)
+ n2 += 1
+
+ # Report some data every now and then
+ if n % self.n_report == 0:
+ self.print_message('track')
+
+ else:
+ self.print_message('end')
+
+ # Define coefficients and centers
+ # ToDo: Check whether everything is clipped here correctly!!
+ self.c = np.concatenate((self.c, c[:n + 1]))
+ self.Cut_ = Cut_[:N_ctrs_so_far + n + 1, :N_ctrs_so_far + n + 1]
+ self.ind_ctrs1 = self.ind_ctrs1[:N_ctrs_so_far1 + n1 + 1] # ToDo: Check whether this is correct!!!!
+ self.ind_ctrs2 = self.ind_ctrs2[:N_ctrs_so_far2 + n2 + 1]
+ self.coef_ = self.Cut_.transpose() @ self.c
+
+ self.ctrs_ = np.concatenate((self.ctrs_,
+ np.concatenate(list_ctrs[:n + 1], axis=0)),
+ axis=0)
+
+ return self
+
+ def predict_s(self, X):
+ # Prediction of s, which might/should approximate the true solution u
+
+ X = check_array(X)
+
+ # Compute prediction
+ if self.ctrs_.shape[0] > 0:
+ prediction = self.kernel.d2_eval(X, self.ctrs_[self.ind_ctrs1]) @ self.coef_[self.ind_ctrs1] \
+ + self.kernel.eval(X, self.ctrs_[self.ind_ctrs2]) @ self.coef_[self.ind_ctrs2]
+ else:
+ prediction = np.zeros((X.shape[0], 1))
+
+ return prediction
+
+ def predict_Ls(self, X):
+ # Prediction of Ls, i.e. application of L to the prediction
+
+ # Validate the input
+ X = check_array(X)
+
+ # Compute prediction
+ if self.ctrs_.shape[0] > 0:
+ prediction = self.kernel.dd_eval(X, self.ctrs_[self.ind_ctrs1]) @ self.coef_[self.ind_ctrs1] \
+ + self.kernel.d1_eval(X, self.ctrs_[self.ind_ctrs2]) @ self.coef_[self.ind_ctrs2]
+ else:
+ prediction = np.zeros((X.shape[0], 1))
+
+ return prediction
+
+ def print_message(self, when):
+ if self.verbose and when == 'begin':
+ print('')
+ print('*' * 30 + ' [VKOGA] ' + '*' * 30)
+ print('Training model with')
+ print(' |_ kernel : %s' % self.kernel)
+ print(' |_ regularization par. : %2.2e' % self.reg_par)
+ print(' |_ restriction par. : %2.2e' % self.restr_par)
+ print('')
+
+ if self.verbose and when == 'end':
+ print('Training completed with')
+ print(' |_ selected points : %8d / %8d' % (
+ self.train_hist['n'][-1], self.ctrs_.shape[0] + self.maxIter))
+ if self.flag_val:
+ print(' |_ train, val residual : %2.2e / %2.2e, %2.2e' %
+ (self.train_hist['f'][-1], self.tol_f, self.train_hist['f val'][-1]))
+ print(' |_ train, val power fun: %2.2e / %2.2e, %2.2e' %
+ (self.train_hist['p'][-1], self.tol_p, self.train_hist['p val'][-1]))
+ else:
+ print(' |_ train residual : %2.2e / %2.2e' % (self.train_hist['f'][-1], self.tol_f))
+ print(' |_ train power fun : %2.2e / %2.2e' % (self.train_hist['p'][-1], self.tol_p))
+
+ if self.verbose and when == 'track':
+ print('Training ongoing with')
+ print(' |_ selected points : %8d / %8d' % (
+ self.train_hist['n'][-1], self.ctrs_.shape[0] + self.maxIter))
+ if self.flag_val:
+ print(' |_ train, val residual : %2.2e / %2.2e, %2.2e' %
+ (self.train_hist['f'][-1], self.tol_f, self.train_hist['f val'][-1]))
+ print(' |_ train, val power fun: %2.2e / %2.2e, %2.2e' %
+ (self.train_hist['p'][-1], self.tol_p, self.train_hist['p val'][-1]))
+ else:
+ print(' |_ train residual : %2.2e / %2.2e' % (self.train_hist['f'][-1], self.tol_f))
+ print(' |_ train power fun : %2.2e / %2.2e' % (self.train_hist['p'][-1], self.tol_p))
+
+
+# %% Utilities to
+import pickle
+
+
+def save_object(obj, filename):
+ with open(filename, 'wb') as output:
+ pickle.dump(obj, output, pickle.HIGHEST_PROTOCOL)
+
+
+def load_object(filename):
+ with open(filename, 'rb') as input:
+ obj = pickle.load(input)
+ return obj