Submitted version (SympGPR)

.gitignore

@@ -0,0 +1,13 @@
+__pycache__
+*.dll
+*.so
+*.pyc
+*.pyd
+*~
+*.mod
+outcmaes
+*.o
+*.pyf
+*.dylib
+*.npz
+*.png

python/01_pendulum/explicit/Makefile

@@ -0,0 +1,13 @@
+FC = gfortran
+FFLAGS = -Wall -march=native -O2 -g -fbacktrace
+PYTHON = python3
+NAME = kernels_sum
+
+all: $(NAME).f90
+	$(PYTHON) -m numpy.f2py -m $(NAME) -c $(NAME).f90 --f90flags='$(FFLAGS)' -lgomp
+
+$(NAME).f90: init_func.py
+	$(PYTHON) init_func.py
+
+clean:
+	rm -f *.x *.so *.o *.mod *.pyf *.pyd *.dylib *.dll $(NAME).f90

python/01_pendulum/explicit/func_expl.py

@@ -0,0 +1,144 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Apr 21 16:47:00 2020
+
+@author: Katharina Rath
+"""
+
+import numpy as np
+from scipy.linalg import solve_triangular
+import scipy
+from scipy.integrate import solve_ivp
+
+
+from kernels_sum import *
+
+
+def f_kern(x, y, x0, y0, l):
+    return kern_num(x,y,x0,y0,l[0], l[1])
+
+def dkdx(x, y, x0, y0, l):
+    return dkdx_num(x,y,x0,y0,l[0], l[1])
+
+def dkdy(x, y, x0, y0, l):
+    return dkdy_num(x,y,x0,y0,l[0], l[1])
+
+def dkdx0(x, y, x0, y0, l):
+    return dkdx0_num(x,y,x0,y0,l[0], l[1])
+
+def dkdy0(x, y, x0, y0, l):
+    return dkdy0_num(x,y,x0,y0,l[0], l[1])
+
+def d2kdxdx0(x, y, x0, y0, l):
+    return d2kdxdx0_num(x,y,x0,y0,l[0], l[1])
+
+def d2kdydy0(x, y, x0, y0, l):
+    return d2kdydy0_num(x,y,x0,y0,l[0], l[1])
+
+def d2kdxdy0(x, y, x0, y0, l):
+    return d2kdxdy0_num(x,y,x0,y0,l[0], l[1])
+
+def d2kdydx0(x, y, x0, y0, l):
+    return d2kdxdy0(x, y, x0, y0, l)
+
+
+def intode(y, t, dtsymp):
+    ys = np.zeros([2, len(t)])
+    ys[0,0] = y[0]
+    ys[1,0] = y[1]            
+
+    for kt in range(1,len(t)):
+            ys[1, kt]=ys[1, kt-1] - dtsymp*(np.sin(ys[0, kt-1] + np.pi))
+            ys[0, kt]=ys[0, kt-1] + dtsymp*ys[1, kt]
+    return ys.T
+
+def build_K(xin, x0in, hyp, K):
+    # set up covariance matrix with derivative observations, Eq. (38)
+    l = hyp[:-1]
+    sig = hyp[-1]
+    N = K.shape[0]//2
+    N0 = K.shape[1]//2
+    x0 = x0in[0:N0]
+    x = xin[0:N]
+    y0 = x0in[N0:2*N0]
+    y = xin[N:2*N]
+    for k in range(N):
+        for lk in range(N0):
+            K[k,lk] = d2kdxdx0(
+                x0[lk], y0[lk], x[k], y[k], l) 
+            K[N+k,lk] = d2kdxdy0(
+                 x0[lk], y0[lk], x[k], y[k], l) 
+            K[k,N0+lk] = d2kdydx0(
+                 x0[lk], y0[lk], x[k], y[k], l) 
+            K[N+k,N0+lk] = d2kdydy0(
+                x0[lk], y0[lk], x[k], y[k], l) 
+    K[:,:] = sig*K[:,:]
+
+def energy(x, U0): 
+    return x[1]**2/2 + U0*(1 + np.cos(x[0]))
+
+
+# compute log-likelihood according to RW, p.19
+def solve_cholesky(L, b):
+    return solve_triangular(
+        L.T, solve_triangular(L, b, lower=True, check_finite=False), 
+        lower=False, check_finite=False)
+
+# negative log-posterior
+def nll_chol(hyp, x, y):
+    K = np.empty((len(x), len(x)))
+    build_K(x, x, hyp[:-1], K)
+    Ky = K + np.abs(hyp[-1])*np.diag(np.ones(len(x)))
+    L = scipy.linalg.cholesky(Ky, lower = True)
+    alpha = solve_cholesky(L, y)
+    ret = 0.5*y.T.dot(alpha) + np.sum(np.log(L.diagonal()))
+    return ret
+
+
+def calcQ(x,y, xtrain, l, Kyinv, ztrain, Ntest):
+    Kstar = np.empty((len(xtrain), 2))
+    build_K(xtrain, np.hstack(([x], [y])), l, Kstar)
+    qGP = Kstar.T.dot(Kyinv.dot(ztrain))
+    f = qGP[1]
+
+    return f, qGP[0]
+
+def calcP(x,y, l, xtrain, ztrain, Kyinv, Ntest):
+    Kstar = np.empty((len(xtrain), 2))
+    build_K(xtrain, np.hstack(([x], [y])), l, Kstar)
+    pGP = Kstar.T.dot(Kyinv.dot(ztrain))
+    res = -pGP[0] 
+    return res
+
+def applymap(l, Q0map, P0map, xtrain, ztrain, Kyinv, Ntest, nm):
+    pmap = np.zeros([nm, Ntest])
+    qmap = np.zeros([nm, Ntest])
+    #set initial conditions
+    pmap[0,:] = P0map
+    qmap[0,:] = Q0map
+    for i in range(0,nm-1):
+        for k in range(0, Ntest): 
+            pmap[i+1, k] = pmap[i, k] + calcP(qmap[i,k], pmap[i, k], l, xtrain, ztrain, Kyinv, Ntest)
+        for k in range(0, Ntest):
+            if np.isnan(pmap[i+1, k]):
+                qmap[i+1,k] = np.nan
+            else: 
+                qmap[i+1, k], dpmap = calcQ(qmap[i,k], pmap[i+1,k],xtrain, l, Kyinv, ztrain, Ntest) 
+                qmap[i+1, k] = np.mod(qmap[i+1,k] + qmap[i, k], 2.0*np.pi)
+    return qmap, pmap
+
+def dydt_ivp(t, y):
+    ydot = np.zeros([2])
+    ydot[0] = y[1] # dx/dt
+    ydot[1] = -np.sin(y[0]+np.pi)# - eps*np.sin(y_[0] - om*(t))#- 0.5*np.pi) # dpx/dt
+    return ydot
+
+def integrate_pendulum(q0, p0, t):
+    ysint = np.zeros([len(t), 2, len(q0)]) # initial values for y
+    ysint = []
+    for ik in range(len(q0)):
+        res_int = solve_ivp(dydt_ivp, [t[0], t[-1]], np.array((q0[ik], p0[ik])), t_eval = t, method='LSODA', rtol = 1e-13, atol = 1e-16)
+        temp = res_int.y
+        ysint.append(temp)
+    return ysint
+

python/01_pendulum/explicit/init_func.py

@@ -0,0 +1,81 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed May  6 14:50:32 2020
+
+@author: ert
+"""
+#%%
+from sympy import *
+from sympy.utilities.autowrap import autowrap, ufuncify
+from sympy.utilities.lambdify import lambdastr, lambdify, implemented_function
+import shutil
+
+tmp = '.'
+#try:
+#    shutil.rmtree(tmp)
+#except:
+#    pass
+
+#%% # Kernel functions (kernel, 1st and second derivatives)
+x, y, xa, ya, xb, yb, l, lx, ly = symbols(r'x y x_a y_a x_b, y_b, l, lx, ly')
+
+
+kern_sqexp = exp(-x**2/(2.0*l**2))
+kern_sqexp_sin = exp(-sin(0.5*x)**2/(2.0*l**2))
+
+kern_y = kern_sqexp.subs(x, ya-yb).subs(l, ly)
+kern_x = kern_sqexp_sin.subs(x, xa-xb).subs(l, lx)
+kern = (kern_x+kern_y).simplify()
+
+dkdxa = diff(kern, xa).simplify()
+dkdxb = diff(kern, xb).simplify()
+dkdya = diff(kern, ya).simplify()
+dkdyb = diff(kern, yb).simplify()
+dkdxadxb = diff(kern, xa, xb).simplify()
+dkdyadyb = diff(kern, ya, yb).simplify()
+dkdxadyb = diff(kern, xa, yb).simplify()
+
+d3kdxdx0dy0 = diff(kern, xa, xb, yb).simplify()
+d3kdydy0dy0 = diff(kern, ya, yb, yb).simplify()
+d3kdxdy0dy0 = diff(kern, xa, yb, yb).simplify()
+
+dkdlx = diff(kern, lx).simplify()
+dkdly = diff(kern, ly).simplify()
+
+d3kdxdx0dlx = diff(kern, xa, xb, lx).simplify()
+d3kdydy0dlx = diff(kern, ya, yb, lx).simplify()
+d3kdxdy0dlx = diff(kern, xa, yb, lx).simplify()
+
+d3kdxdx0dly = diff(kern, xa, xb, ly).simplify()
+d3kdydy0dly = diff(kern, ya, yb, ly).simplify()
+d3kdxdy0dly = diff(kern, xa, yb, ly).simplify()
+
+#%%
+from sympy.utilities.codegen import codegen
+seq = (xa, ya, xb, yb, lx, ly)
+[(name, code), (h_name, header)] = \
+codegen([('kern_num', kern),
+         ('dkdx_num', dkdxa),
+         ('dkdy_num', dkdya),
+         ('dkdx0_num', dkdxb),
+         ('dkdy0_num', dkdyb),
+         ('d2kdxdx0_num', dkdxadxb),
+         ('d2kdydy0_num', dkdyadyb),
+         ('d2kdxdy0_num', dkdxadyb),
+         ('d3kdxdx0dy0_num', d3kdxdx0dy0),
+         ('d3kdydy0dy0_num', d3kdydy0dy0),
+         ('d3kdxdy0dy0_num', d3kdxdy0dy0),
+         ('dkdlx_num', dkdlx),
+         ('dkdly_num', dkdly),
+         ('d3kdxdx0dlx_num', d3kdxdx0dlx),
+         ('d3kdydy0dlx_num', d3kdydy0dlx),
+         ('d3kdxdy0dlx_num', d3kdxdy0dlx),
+         ('d3kdxdx0dly_num', d3kdxdx0dly),
+         ('d3kdydy0dly_num', d3kdydy0dly),
+         ('d3kdxdy0dly_num', d3kdxdy0dly),
+         ], "F95", "kernels_sum",
+        argument_sequence=seq,
+        header=False, empty=False)
+with open(name, 'w') as f:
+    f.write(code)

python/01_pendulum/explicit/kernels.py

@@ -0,0 +1,15 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Aug  9 22:30:35 2020
+
+@author: manal
+"""
+
+def __bootstrap__():
+   global __bootstrap__, __loader__, __file__
+   import sys, pkg_resources, imp
+   __file__ = pkg_resources.resource_filename(__name__,'kernels_sum.cpython-36m-x86_64-linux-gnu.so')
+   __loader__ = None; del __bootstrap__, __loader__
+   imp.load_dynamic(__name__,__file__)
+__bootstrap__()