Skip to content
Snippets Groups Projects
Commit 2ed0814e authored by Wenzel, Tizian's avatar Wenzel, Tizian
Browse files

Initial commit.

parent 4e906c04
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
.gitignore 0 → 100644
+ 156
− 0
View file @ 2ed0814e
# 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
example.py 0 → 100644
+ 73
− 0
View file @ 2ed0814e
# 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.')
kernels_PDE.py 0 → 100644
+ 237
− 0
View file @ 2ed0814e
#!/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)
vkoga_PDE.py 0 → 100644
+ 594
− 0
View file @ 2ed0814e
# 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment