{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nCreated on Mon Nov 8 12:15:09 2021\n\n@author: maout\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "###Systematic multi-N inhomogeneous Kuramoto\n\nimport numpy as np\n\n\n#from matplotlib import pyplot as plt\n\nfrom scipy.spatial.distance import cdist\nimport time\nimport ot\nimport numba\nimport math\nimport random\nimport sys\nimport pickle\nfrom functools import reduce\n\n@numba.njit(parallel=True,fastmath=True)\ndef my_cdist(r,y, output,dist='euclidean'):\n \"\"\" \n Fast computation of pairwise distances between data points in r and y matrices.\n Stores the distances in the output array.\n Available distances: 'euclidean' and 'seucledian'\n Parameters\n ----------\n r : NxM array\n First set of N points of dimension M.\n y : N2xM array\n Second set of N2 points of dimension M.\n output : NxN2 array\n Placeholder for storing the output of the computed distances.\n dist : type of distance, optional\n Select 'euclidian' or 'sqeuclidian' for Euclidian or squared Euclidian\n distances. The default is 'euclidean'.\n\n Returns\n -------\n None. (The result is stored in place in the input array output).\n\n \"\"\"\n N, M = r.shape\n N2, M2 = y.shape\n #assert( M == M2, 'The two inpus have different second dimention! Input should be N1xM and N2xM')\n if dist == 'euclidean':\n for i in numba.prange(N):\n for j in numba.prange(N2):\n tmp = 0.0\n for k in range(M):\n tmp += (r[i, k] - y[j, k])**2 \n output[i,j] = math.sqrt(tmp)\n elif dist == 'sqeuclidean':\n for i in numba.prange(N):\n for j in numba.prange(N2):\n tmp = 0.0\n for k in range(M):\n tmp += (r[i, k] - y[j, k])**2 \n output[i,j] = tmp \n elif dist == 'l1':\n for i in numba.prange(N):\n for j in numba.prange(N2):\n tmp = 0.0\n for k in range(M):\n tmp += (r[i, k] - y[j, k])**2 \n output[i,j] = math.sqrt(tmp) \n return 0\n\ndef score_function_multid_seperate(X,Z,func_out=False, C=0.001,kern ='RBF',l=1,which=1,which_dim=1):\n \"\"\"\n returns function psi(z)\n Input: X: N observations\n Z: sparse points\n func_out : Boolean, True returns function if False return grad-log-p on data points \n l: lengthscale of rbf kernel\n C: weighting constant \n which: return 1: grad log p(x) \n which_dim: which gradient of log density we want to compute (starts from 1 for the 0-th dimension)\n Output: psi: array with density along the given dimension N or N_s x 1\n \n \"\"\"\n \n if kern=='RBF':\n #l = 1 # lengthscale of RBF kernel\n #@numba.njit(parallel=True,fastmath=True)\n def Knumba(x,y,l,res,multil=False): #version of kernel in the numba form when the call already includes the output matrix\n if multil: \n #print('here')\n #res = np.ones((x.shape[0],y.shape[0])) \n for ii in range(len(l)): \n tempi = np.zeros((x[:,ii].size, y[:,ii].size ), dtype=np.float64)\n ##puts into tempi the cdist result\n #print(x[:,ii:ii+1].shape)\n my_cdist(x[:,ii:ii+1], y[:,ii:ii+1],tempi,'sqeuclidean')\n \n res = np.multiply(res,np.exp(-tempi/(2*l[ii]*l[ii]))) \n ##res = np.multiply(res,np.exp(-cdist(x[:,ii].reshape(-1,1), y[:,ii].reshape(-1,1),'sqeuclidean')/(2*l[ii]*l[ii])))\n #return res\n else:\n tempi = np.zeros((x.shape[0], y.shape[0] ), dtype=np.float64)\n #return np.exp(-cdist(x, y,'sqeuclidean')/(2*l*l))\n my_cdist(x, y,tempi,'sqeuclidean') #this sets into the array tempi the cdist result\n res = np.exp(-tempi/(2*l*l))\n return 0\n \n def K(x,y,l,multil=False):\n if multil: \n #print('here')\n res = np.ones((x.shape[0],y.shape[0])) \n for ii in range(len(l)): \n tempi = np.zeros((x[:,ii].size, y[:,ii].size ))\n ##puts into tempi the cdist result\n my_cdist(x[:,ii].reshape(-1,1), y[:,ii].reshape(-1,1),tempi,'sqeuclidean')\n res = np.multiply(res,np.exp(-tempi/(2*l[ii]*l[ii]))) \n ##res = np.multiply(res,np.exp(-cdist(x[:,ii].reshape(-1,1), y[:,ii].reshape(-1,1),'sqeuclidean')/(2*l[ii]*l[ii])))\n return res\n else:\n tempi = np.zeros((x.shape[0], y.shape[0] ))\n #return np.exp(-cdist(x, y,'sqeuclidean')/(2*l*l))\n my_cdist(x, y,tempi,'sqeuclidean') #this sets into the array tempi the cdist result\n return np.exp(-tempi/(2*l*l))\n #return np.exp(-(x-y.T)**2/(2*l*l))\n #return np.exp(np.linalg.norm(x-y.T, 2)**2)/(2*l*l) \n #@njit\n def grdx_K(x,y,l,which_dim=1,multil=False): #gradient with respect to the 1st argument - only which_dim\n N,dim = x.shape \n diffs = x[:,None]-y \n redifs = np.zeros((1*N,N))\n ii = which_dim -1\n #print('diffs:',diffs)\n if multil:\n redifs = np.multiply(diffs[:,:,ii],K(x,y,l,True))/(l[ii]*l[ii]) \n else:\n redifs = np.multiply(diffs[:,:,ii],K(x,y,l))/(l*l) \n return redifs\n #return -(1./(l*l))*(x-y.T)*K(x,y)\n \n def grdy_K(x,y): # gradient with respect to the second argument\n N,dim = x.shape\n diffs = x[:,None]-y \n redifs = np.zeros((N,N))\n ii = which_dim -1 \n redifs = np.multiply(diffs[:,:,ii],K(x,y,l))/(l*l) \n return -redifs\n #return (1./(l*l))*(x-y.T)*K(x,y)\n #@njit \n def ggrdxy_K(x,y):\n N,dim = Z.shape\n diffs = x[:,None]-y \n \n redifs = np.zeros((N,N))\n for ii in range(which_dim-1,which_dim): \n for jj in range(which_dim-1,which_dim):\n redifs[ii, jj ] = np.multiply(np.multiply(diffs[:,:,ii],diffs[:,:,jj])+(l*l)*(ii==jj),K(x,y))/(l**4) \n return -redifs\n #return np.multiply((K(x,y)),(np.power(x[:,None]-y,2)-l**2))/l**4\n #############################################################################\n elif kern=='periodic': ###############################################################################################\n ###periodic kernel\n ## K(x,y) = exp( -2 * sin^2( pi*| x-y |/ (2*pi) ) /l^2)\n \n ## Kx(x,y) = (K(x,y)* (x - y) cos(abs(x - y)/2) sin(abs(x - y)/2))/(l^2 abs(x - y))\n ## -(2 K(x,y) \u03c0 (x - y) sin((2 \u03c0 abs(x - y))/per))/(l^2 s abs(x - y))\n per = 2*np.pi ##period of the kernel\n #l = 0.5\n def K(x,y,l,multil=False):\n \n if multil: \n #print('here')\n res = np.ones((x.shape[0],y.shape[0])) \n for ii in range(len(l)): \n tempi = np.zeros((x[:,ii].size, y[:,ii].size ))\n ##puts into tempi the cdist result\n #my_cdist(x[:,ii].reshape(-1,1), y[:,ii].reshape(-1,1),tempi, 'l1') \n #res = np.multiply(res, np.exp(- 2* (np.sin(tempi/ 2 )**2) /(l[ii]*l[ii])) )\n res = np.multiply(res, np.exp(- 2* (np.sin(cdist(x[:,ii].reshape(-1,1), y[:,ii].reshape(-1,1),'minkowski', p=1)/ 2 )**2) /(l[ii]*l[ii])) )\n return -res\n else:\n tempi = np.zeros((x.shape[0], y.shape[0] ))\n ##puts into tempi the cdist result\n #my_cdist(x, y, tempi,'l1')\n #res = np.exp(-2* ( np.sin( tempi / 2 )**2 ) /(l*l) )\n res = np.exp(-2* ( np.sin( cdist(x, y,'minkowski', p=1) / 2 )**2 ) /(l*l) )\n return res\n \n def grdx_K(x,y,l,which_dim=1,multil=False): #gradient with respect to the 1st argument - only which_dim\n N,dim = x.shape \n diffs = x[:,None]-y \n #print('diffs:',diffs)\n redifs = np.zeros((1*N,N))\n ii = which_dim -1\n #print(ii)\n if multil:\n redifs = np.divide( np.multiply( np.multiply( np.multiply( -2*K(x,y,l,True),diffs[:,:,ii] ),np.sin( np.abs(diffs[:,:,ii]) / 2) ) ,np.cos( np.abs(diffs[:,:,ii]) / 2) ) , (l[ii]*l[ii]* np.abs(diffs[:,:,ii])) ) \n else:\n redifs = np.divide( np.multiply( np.multiply( np.multiply( -2*diffs[:,:,ii],np.sin( np.abs(diffs[:,:,ii]) / 2) ) ,K(x,y,l) ),np.cos( np.abs(diffs[:,:,ii]) / 2) ) ,(l*l* np.abs(diffs[:,:,ii])) ) \n return -redifs\n\n\n\n if isinstance(l, (list, tuple, np.ndarray)):\n ### for different lengthscales for each dimension \n #K_xz = np.ones((X.shape[0],Z.shape[0]), dtype=np.float64) \n #Knumba(X,Z,l,K_xz,multil=True) \n K_xz = K(X,Z,l,multil=True) \n #Ks = np.ones((Z.shape[0],Z.shape[0]), dtype=np.float64) \n #Knumba(Z,Z,l,Ks,multil=True) \n Ks = K(Z,Z,l,multil=True) \n multil = True\n #print(Z.shape)\n Ksinv = np.linalg.inv(Ks+ 1e-3 * np.eye(Z.shape[0]))\n A = K_xz.T @ K_xz \n gradx_K = -grdx_K(X,Z,l,which_dim=which_dim,multil=True) #-\n # if not(Test_p == 'None'):\n # K_sz = K(Test_p,Z,l,multil=True)\n \n else:\n multil = False\n \n K_xz = K(X,Z,l,multil=False) \n \n Ks = K(Z,Z,l,multil=False) \n \n Ksinv = np.linalg.inv(Ks+ 1e-3 * np.eye(Z.shape[0]))\n A = K_xz.T @ K_xz \n \n gradx_K = -grdx_K(X,Z,l,which_dim=which_dim,multil=False)\n sumgradx_K = np.sum(gradx_K ,axis=0)\n #print( sumgradx_K.shape )\n if func_out==False: #if output wanted is evaluation at data points\n ### evaluatiion at data points\n res1 = -K_xz @ np.linalg.inv( C*np.eye(Z.shape[0], Z.shape[0]) + Ksinv @ A + 1e-3 * np.eye(Z.shape[0]))@ Ksinv@sumgradx_K\n else: \n #### for function output \n if multil: \n #res = np.ones((x.shape[0],y.shape[0])) \n #for ii in range(len(l)): \n if kern=='RBF': \n K_sz = lambda x: reduce(np.multiply, [ np.exp(-cdist(x[:,iii].reshape(-1,1), Z[:,iii].reshape(-1,1),'sqeuclidean')/(2*l[iii]*l[iii])) for iii in range(x.shape[1]) ])\n \n \n elif kern=='periodic':\n K_sz = lambda x: np.multiply(np.exp(-2*(np.sin( cdist(x[:,0].reshape(-1,1), Z[:,0].reshape(-1,1), 'minkowski', p=2)/(l[0]*l[0])))),np.exp(-2*(np.sin( cdist(x[:,1].reshape(-1,1), Z[:,1].reshape(-1,1),'sqeuclidean')/(l[1]*l[1])))))\n \n \n #return K_sz\n else:\n if kern=='RBF':\n K_sz = lambda x: np.exp(-cdist(x, Z,'sqeuclidean')/(2*l*l))\n elif kern=='periodic':\n K_sz = lambda x: np.exp(-2* ( np.sin( cdist(x, Z,'minkowski', p=1) / 2 )**2 ) /(l*l) )\n #return K_sz\n\n res1 = lambda x: K_sz(x) @ ( -np.linalg.inv( C*np.eye(Z.shape[0], Z.shape[0]) + Ksinv @ A + 1e-3 * np.eye(Z.shape[0])) ) @ Ksinv@sumgradx_K\n\n\n \n return res1\n\n\n\n\n\ndef score_function_multid_seperate_all_dims(X,Z,func_out=False, C=0.001,kern ='RBF',l=1,which=1):\n \"\"\"\n returns function psi(z)\n Input: X: N observations\n Z: sparse points\n func_out : Boolean, True returns function if False return grad-log-p on data points \n l: lengthscale of rbf kernel\n C: weighting constant \n which: return 1: grad log p(x) \n \n Output: psi: array with density along the given dimension N or N_s x 1\n \n \"\"\"\n \n if kern=='RBF':\n #l = 1 # lengthscale of RBF kernel\n #@numba.njit(parallel=True,fastmath=True)\n def Knumba(x,y,l,res,multil=False): #version of kernel in the numba form when the call already includes the output matrix\n if multil: \n #print('here')\n #res = np.ones((x.shape[0],y.shape[0])) \n for ii in range(len(l)): \n tempi = np.zeros((x[:,ii].size, y[:,ii].size ), dtype=np.float64)\n ##puts into tempi the cdist result\n #print(x[:,ii:ii+1].shape)\n my_cdist(x[:,ii:ii+1], y[:,ii:ii+1],tempi,'sqeuclidean')\n \n res = np.multiply(res,np.exp(-tempi/(2*l[ii]*l[ii]))) \n ##res = np.multiply(res,np.exp(-cdist(x[:,ii].reshape(-1,1), y[:,ii].reshape(-1,1),'sqeuclidean')/(2*l[ii]*l[ii])))\n #return res\n else:\n tempi = np.zeros((x.shape[0], y.shape[0] ), dtype=np.float64)\n #return np.exp(-cdist(x, y,'sqeuclidean')/(2*l*l))\n my_cdist(x, y,tempi,'sqeuclidean') #this sets into the array tempi the cdist result\n res = np.exp(-tempi/(2*l*l))\n return 0\n \n def K(x,y,l,multil=False):\n if multil: \n #print('here')\n res = np.ones((x.shape[0],y.shape[0])) \n for ii in range(len(l)): \n tempi = np.zeros((x[:,ii].size, y[:,ii].size ))\n ##puts into tempi the cdist result\n my_cdist(x[:,ii].reshape(-1,1), y[:,ii].reshape(-1,1),tempi,'sqeuclidean')\n res = np.multiply(res,np.exp(-tempi/(2*l[ii]*l[ii]))) \n ##res = np.multiply(res,np.exp(-cdist(x[:,ii].reshape(-1,1), y[:,ii].reshape(-1,1),'sqeuclidean')/(2*l[ii]*l[ii])))\n return res\n else:\n tempi = np.zeros((x.shape[0], y.shape[0] ))\n #return np.exp(-cdist(x, y,'sqeuclidean')/(2*l*l))\n my_cdist(x, y,tempi,'sqeuclidean') #this sets into the array tempi the cdist result\n return np.exp(-tempi/(2*l*l))\n #return np.exp(-(x-y.T)**2/(2*l*l))\n #return np.exp(np.linalg.norm(x-y.T, 2)**2)/(2*l*l) \n #@njit\n def grdx_K_all(x,y,l,multil=False): #gradient with respect to the 1st argument - only which_dim\n N,dim = x.shape \n M,_ = y.shape\n diffs = x[:,None]-y \n redifs = np.zeros((1*N,M,dim))\n for ii in range(dim): \n \n if multil:\n redifs[:,:,ii] = np.multiply(diffs[:,:,ii],K(x,y,l,True))/(l[ii]*l[ii]) \n else:\n redifs[:,:,ii] = np.multiply(diffs[:,:,ii],K(x,y,l))/(l*l) \n return redifs\n #return -(1./(l*l))*(x-y.T)*K(x,y)\n \n def grdx_K(x,y,l,which_dim=1,multil=False): #gradient with respect to the 1st argument - only which_dim\n N,dim = x.shape \n M,_ = y.shape\n diffs = x[:,None]-y \n redifs = np.zeros((1*N,M))\n ii = which_dim -1\n #print('diffs:',diffs)\n if multil:\n redifs = np.multiply(diffs[:,:,ii],K(x,y,l,True))/(l[ii]*l[ii]) \n #print(redifs.shape)\n else:\n redifs = np.multiply(diffs[:,:,ii],K(x,y,l))/(l*l) \n return redifs\n \n \n #############################################################################\n elif kern=='periodic': ###############################################################################################\n ###periodic kernel\n ## K(x,y) = exp( -2 * sin^2( pi*| x-y |/ (2*pi) ) /l^2)\n \n ## Kx(x,y) = (K(x,y)* (x - y) cos(abs(x - y)/2) sin(abs(x - y)/2))/(l^2 abs(x - y))\n ## -(2 K(x,y) \u03c0 (x - y) sin((2 \u03c0 abs(x - y))/per))/(l^2 s abs(x - y))\n per = 2*np.pi ##period of the kernel\n #l = 0.5\n def K(x,y,l,multil=False):\n \n if multil: \n #print('here')\n res = np.ones((x.shape[0],y.shape[0])) \n for ii in range(len(l)): \n tempi = np.zeros((x[:,ii].size, y[:,ii].size ))\n ##puts into tempi the cdist result\n #my_cdist(x[:,ii].reshape(-1,1), y[:,ii].reshape(-1,1),tempi, 'l1') \n #res = np.multiply(res, np.exp(- 2* (np.sin(tempi/ 2 )**2) /(l[ii]*l[ii])) )\n res = np.multiply(res, np.exp(- 2* (np.sin(cdist(x[:,ii].reshape(-1,1), y[:,ii].reshape(-1,1),'minkowski', p=1)/ 2 )**2) /(l[ii]*l[ii])) )\n return -res\n else:\n tempi = np.zeros((x.shape[0], y.shape[0] ))\n ##puts into tempi the cdist result\n #my_cdist(x, y, tempi,'l1')\n #res = np.exp(-2* ( np.sin( tempi / 2 )**2 ) /(l*l) )\n res = np.exp(-2* ( np.sin( cdist(x, y,'minkowski', p=1) / 2 )**2 ) /(l*l) )\n return res\n \n def grdx_K(x,y,l,which_dim=1,multil=False): #gradient with respect to the 1st argument - only which_dim\n N,dim = x.shape \n diffs = x[:,None]-y \n #print('diffs:',diffs)\n redifs = np.zeros((1*N,N))\n ii = which_dim -1\n #print(ii)\n if multil:\n redifs = np.divide( np.multiply( np.multiply( np.multiply( -2*K(x,y,l,True),diffs[:,:,ii] ),np.sin( np.abs(diffs[:,:,ii]) / 2) ) ,np.cos( np.abs(diffs[:,:,ii]) / 2) ) , (l[ii]*l[ii]* np.abs(diffs[:,:,ii])) ) \n else:\n redifs = np.divide( np.multiply( np.multiply( np.multiply( -2*diffs[:,:,ii],np.sin( np.abs(diffs[:,:,ii]) / 2) ) ,K(x,y,l) ),np.cos( np.abs(diffs[:,:,ii]) / 2) ) ,(l*l* np.abs(diffs[:,:,ii])) ) \n return -redifs\n\n dim = X.shape[1]\n\n if isinstance(l, (list, tuple, np.ndarray)):\n multil = True\n ### for different lengthscales for each dimension \n #K_xz = np.ones((X.shape[0],Z.shape[0]), dtype=np.float64) \n #Knumba(X,Z,l,K_xz,multil=True) \n K_xz = K(X,Z,l,multil=True) \n #Ks = np.ones((Z.shape[0],Z.shape[0]), dtype=np.float64) \n #Knumba(Z,Z,l,Ks,multil=True) \n Ks = K(Z,Z,l,multil=True) \n \n #print(Z.shape)\n Ksinv = np.linalg.inv(Ks+ 1e-3 * np.eye(Z.shape[0]))\n A = K_xz.T @ K_xz \n \n gradx_K = -grdx_K_all(X,Z,l,multil=True) #-\n gradxK = np.zeros((X.shape[0],Z.shape[0],dim))\n for ii in range(dim):\n gradxK[:,:,ii] = -grdx_K(X,Z,l,multil=True,which_dim=ii+1)\n # if not(Test_p == 'None'):\n # K_sz = K(Test_p,Z,l,multil=True)\n np.testing.assert_allclose(gradxK, gradx_K) \n else:\n multil = False\n \n K_xz = K(X,Z,l,multil=False) \n \n Ks = K(Z,Z,l,multil=False) \n \n Ksinv = np.linalg.inv(Ks+ 1e-3 * np.eye(Z.shape[0]))\n A = K_xz.T @ K_xz \n \n gradx_K = -grdx_K_all(X,Z,l,multil=False) #shape: (N,M,dim)\n sumgradx_K = np.sum(gradx_K ,axis=0) ##last axis will have the gradient for each dimension ### shape (M, dim)\n #print( sumgradx_K.shape )\n if func_out==False: #if output wanted is evaluation at data points\n \n # res1 = np.zeros((N, dim)) \n # ### evaluatiion at data points\n # for di in range(dim):\n # res1[:,di] = -K_xz @ np.linalg.inv( C*np.eye(Z.shape[0], Z.shape[0]) + Ksinv @ A + 1e-3 * np.eye(Z.shape[0]))@ Ksinv@sumgradx_K[:,di]\n \n \n res1 = -K_xz @ np.linalg.inv( C*np.eye(Z.shape[0], Z.shape[0]) + Ksinv @ A + 1e-3 * np.eye(Z.shape[0]))@ Ksinv@sumgradx_K\n \n \n #res1 = np.einsum('ik,kj->ij', -K_xz @ np.linalg.inv( C*np.eye(Z.shape[0], Z.shape[0]) + Ksinv @ A + 1e-3 * np.eye(Z.shape[0]))@ Ksinv, sumgradx_K)\n \n \n else: \n #### for function output \n if multil: \n #res = np.ones((x.shape[0],y.shape[0])) \n #for ii in range(len(l)): \n if kern=='RBF': \n K_sz = lambda x: reduce(np.multiply, [ np.exp(-cdist(x[:,iii].reshape(-1,1), Z[:,iii].reshape(-1,1),'sqeuclidean')/(2*l[iii]*l[iii])) for iii in range(x.shape[1]) ])\n \n \n elif kern=='periodic':\n K_sz = lambda x: np.multiply(np.exp(-2*(np.sin( cdist(x[:,0].reshape(-1,1), Z[:,0].reshape(-1,1), 'minkowski', p=2)/(l[0]*l[0])))),np.exp(-2*(np.sin( cdist(x[:,1].reshape(-1,1), Z[:,1].reshape(-1,1),'sqeuclidean')/(l[1]*l[1])))))\n \n \n #return K_sz\n else:\n if kern=='RBF':\n K_sz = lambda x: np.exp(-cdist(x, Z,'sqeuclidean')/(2*l*l))\n elif kern=='periodic':\n K_sz = lambda x: np.exp(-2* ( np.sin( cdist(x, Z,'minkowski', p=1) / 2 )**2 ) /(l*l) )\n #return K_sz\n\n res1 = lambda x: K_sz(x) @ ( -np.linalg.inv( C*np.eye(Z.shape[0], Z.shape[0]) + Ksinv @ A + 1e-3 * np.eye(Z.shape[0])) ) @ Ksinv@sumgradx_K\n # res1 = np.zeros((N, dim))\n # for di in range(dim):\n # res1 = lambda x: K_sz(x) @ ( -np.linalg.inv( C*np.eye(Z.shape[0], Z.shape[0]) + Ksinv @ A + 1e-3 * np.eye(Z.shape[0])) ) @ Ksinv@sumgradx_K[:,di]\n \n \n #np.testing.assert_allclose(res2, res1)\n \n return res1 ### shape out N x dim" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class BRIDGE_ND_reweight:\n def __init__(self,t1,t2,y1,y2,f,g,N,k,M,reweight=False, U=1,dens_est='nonparametric',reject=True,plotting=True,kern='RBF',dt=0.001):\n \"\"\"\n Bridge initialising function\n t1: starting time point\n t2: end time point\n y1: initial observation/position\n y2: end observation/position\n f: drift function handler\n g: diffusion coefficient or function handler \n N: number of particles/trajectories\n k: discretisation steps within bridge \n M: number of sparse points for grad log density estimation\n reweight: boolean - determines if reweighting will follow\n U: function, reweighting function to be employed during reweighting: dim_y1 \\to 1\n dens_est: density estimation function\n > 'nonparametric' : non parametric density estimation\n > 'hermit1' : parametic density estimation empoying hermite polynomials (physiscist's)\n > 'hermit2' : parametic density estimation empoying hermite polynomials (probabilists's)\n > 'poly' : parametic density estimation empoying simple polynomials\n > 'rbf' : parametric density estimation employing radial basis functions\n kern: type of kernel: 'RBF' or 'periodic'\n reject: boolean parameter indicating whether non valid bridge trajectories will be rejected\n plotting: boolean parameter indicating whether bridge statistics will be plotted\n dt: integration time step [default: 0.001]\n \"\"\"\n self.dim = y1.shape[0] # dimensionality of the problem\n self.t1 = t1\n self.t2 = t2\n self.y1 = y1\n self.y2 = y2\n\n \n ##density estimation stuff\n self.kern = kern\n # DRIFT /DIFFUSION\n self.f = f\n self.g = g #scalar or array\n \n ### PARTICLES DISCRETISATION\n self.N = N \n self.k = k\n self.N_sparse = M\n \n self.dt = dt#0.001 #((t2-t1)/k)\n \n ### reweighting\n self.reweight = reweight\n if self.reweight:\n self.U = U\n\n ### reject\n self.reject = reject\n\n \n self.finer = 1#200 #discetasation ratio between numerical BW solution and particle bridge solution\n self.timegrid = np.arange(self.t1,self.t2+self.dt/2,self.dt)\n #self.timegrid_fine = np.arange(self.t1, self.t2+self.dt*(1./self.finer)/2, self.dt*(1./self.finer) )\n \n # print(self.k == self.timegrid.size)\n # print(self.timegrid)\n \n self.Z = np.zeros((self.dim,self.N,self.k)) #storage for forward trajectories\n self.B = np.zeros((self.dim,self.N,self.k)) #storage for backward trajectories\n self.ln_roD = [] \n self.BPWE = np.zeros((self.dim,self.N,self.timegrid.size))\n self.BPWEmean = np.zeros((self.dim,self.k*self.finer))\n self.BPWEstd = np.zeros((self.dim,self.k*self.finer))\n self.BPWEskew = np.zeros((self.dim,self.k*self.finer))\n self.BPWEkurt = np.zeros((self.dim,self.k*self.finer))\n \n print('start forward')\n #self.forward_sampling()\n self.forward_sampling_Otto()\n # plt.figure(figsize=(6,4)),plt.plot(self.Z[0].T,self.Z[1].T,alpha=0.3);\n # plt.plot(self.y1[0],self.y1[1],'go')\n # plt.plot(self.y2[0],self.y2[1],'ro')\n # plt.show()\n\n \n #self.density_estimation()\n self.backward_simulation()\n #self.reject_trajectories() \n #self.calculate_true_statistics()\n #if plotting:\n # self.plot_statistics()\n \n def forward_sampling(self): \n print('Sampling forward...')\n W = np.ones((self.N,1))/self.N\n for ti,tt in enumerate(self.timegrid):\n\n if ti == 0:\n self.Z[:,:,0] = self.y1\n #self.Z[1,:,0] = self.y1[1]\n else:\n for i in range(self.N):\n #self.Z[:,i,:] = sdeint.itoint(self.f, self.g, self.Z[i,0], self.timegrid)[:,0] \n self.Z[:,i,ti] = ( self.Z[:,i,ti-1] + self.dt* self.f(self.Z[:,i,ti-1]) + \\\n (self.g)*np.random.normal(loc = 0.0, scale = np.sqrt(self.dt),size=(self.dim,)) )% (2*np.pi)\n \n ###WEIGHT\n if self.reweight == True:\n if ti>0:\n W[:,0] = np.exp(1*self.dt*self.U(self.Z[:,:,ti])) \n W = W/np.sum(W)\n \n ###REWEIGHT \n #Tstar = reweight_optimal_transport_multidim(self.Z[:,:,ti].T,W)\n #P = Tstar *N\n # print(Tstar.shape)\n # print(X.shape)\n #self.Z[:,:,ti] = ( (self.Z[:,:,ti])@Tstar )% (2*np.pi)\n M = ot.dist(self.Z[:,:,ti].T, self.Z[:,:,ti].T)\n M /= M.max()\n a = W[:,0]\n b = np.ones_like(W[:,0])/self.N\n T2 = ot.emd(a, b, M)\n #T2 = ot.sinkhorn(a, b, M, 0.1)\n #T2 = ot.bregman.sinkhorn_epsilon_scaling(a, b, M, 0.01)\n #T2 = ot.bregman.sinkhorn_stabilized(a, b, M, 0.001)\n #T2 = ot.optim.cg(a, b, M, reg, fi, dfi, verbose=False)\n self.Z[:,:,ti] = (self.N*self.Z[:,:,ti]@T2) % (2*np.pi)\n \n #for di in range(self.dim):\n # self.Z[di,:,-1] = self.y2[di]\n print('Forward sampling done!')\n return 0\n \n ### effective forward drift - estimated seperatelly for each dimension\n def f_seperate(self,x,t):#plain GP prior\n \n dimi, N = x.shape \n bnds = np.zeros((dimi,2))\n for ii in range(dimi):\n bnds[ii] = [np.min(x[ii,:]),np.max(x[ii,:])]\n sum_bnds = np.sum(bnds)\n \n\n Sxx = np.array([ np.random.uniform(low=bnd[0],high=bnd[1],size=(self.N_sparse)) for bnd in bnds ] ) \n \n lnthsc = 2*np.std(x,axis=1) \n # gpsi2 = np.zeros((dimi, N)) \n # for ii in range(dimi): \n # gpsi2[ii,:]= score_function_multid_seperate(x.T,Sxx.T,False,C=0.001,which=1,l=lnthsc,which_dim=ii+1, kern=self.kern) \n gpsi = score_function_multid_seperate_all_dims(x.T,Sxx.T,False,C=0.001,which=1,l=lnthsc, kern=self.kern).T \n #\n #np.testing.assert_allclose(gpsi, gpsi2)\n return (self.f(x,t)-0.5* self.g**2* gpsi)\n \n \n def forward_sampling_Otto(self):\n print('Sampling forward with deterministic particles...')\n W = np.ones((self.N,1))/self.N\n for ti,tt in enumerate(self.timegrid): \n print(ti) \n if ti == 0:\n # for di in range(self.dim):\n # self.Z[di,:,0] = self.y1[di] \n self.Z[:,:,0] = np.random.multivariate_normal(self.y1, np.eye(self.dim)*self.g/2, self.N).T\n elif ti==1: #propagate one step with stochastic to avoid the delta function\n #for i in range(self.N): #substract dt because I want the time at t-1\n self.Z[:,:,ti] = (self.Z[:,:,ti-1] + self.dt*self.f(self.Z[:,:,ti-1],tt-self.dt)+\\\n (self.g)*np.random.normal(loc = 0.0, scale = np.sqrt(self.dt),size=(self.dim,self.N)) )% (2*np.pi)\n else:\n \n self.Z[:,:,ti] = ( self.Z[:,:,ti-1] + self.dt* self.f_seperate(self.Z[:,:,ti-1],tt-self.dt) )% (2*np.pi)\n ###WEIGHT\n if self.reweight == True:\n if ti>0:\n W[:,0] = np.exp(self.U(self.Z[:,:,ti])) #-1 \n #Neff = 1/np.sum(W[:,0]**2)\n #print('Neff:', Neff)\n #print(W[:,0])\n W = W/np.sum(W) \n #Neff = 1/np.sum(W[:,0]**2)\n #print('Neff:', Neff)\n #print(W[:,0])\n #print('-----')\n ###REWEIGHT \n #start = time.time()\n #Tstar = reweight_optimal_transport_multidim(self.Z[:,:,ti].T,W)\n #print(Tstar)\n # if ti ==3:\n # stop = time.time()\n # print('Timepoint: %d needed '%ti, stop-start)\n #P = Tstar *N\n # print(Tstar.shape)\n # print(X.shape)\n #self.Z[:,:,ti] = ((self.Z[:,:,ti])@Tstar )% (2*np.pi) ##### \n M = ot.dist(self.Z[:,:,ti].T, self.Z[:,:,ti].T)\n M /= M.max()\n a = W[:,0]\n b = np.ones_like(W[:,0])/self.N\n T2 = ot.emd(a, b, M,numItermax=1000000)\n #T2 = ot.sinkhorn(a, b, M, 0.1)\n #T2 = ot.bregman.sinkhorn_epsilon_scaling(a, b, M, 0.01)\n #T2 = ot.bregman.sinkhorn_stabilized(a, b, M, 0.001)\n #T2 = ot.optim.cg(a, b, M, reg, fi, dfi, verbose=False)\n self.Z[:,:,ti] = (self.N*self.Z[:,:,ti]@T2) % (2*np.pi)\n \n print('Forward sampling with Otto is ready!')\n \n return 0\n \n def density_estimation(self, ti,rev_ti):\n rev_t = rev_ti-1#########################################################-1\n grad_ln_ro = np.zeros((self.dim,self.N))\n lnthsc = 2*np.std(self.Z[:,:,rev_t],axis=1)\n \n bnds = np.zeros((self.dim,2))\n for ii in range(self.dim):\n bnds[ii] = [max(np.min(self.Z[ii,:,rev_t]),np.min(self.B[ii,:,rev_ti])),min(np.max(self.Z[ii,:,rev_t]),np.max(self.B[ii,:,rev_ti]))]\n sum_bnds = np.sum(bnds)\n #print(bnds)\n #print(np.min(self.B[:,:,rev_ti]))\n #print(np.max(self.B[:,:,rev_ti]))\n #print('-------')\n \n #sparse points\n Sxx = np.array([ np.random.uniform(low=bnd[0],high=bnd[1],size=(self.N_sparse)) for bnd in bnds ] )\n \n # for di in range(self.dim): \n # #estimate density from forward (Z) and evaluate at current postitions of backward particles (B) \n # grad_ln_ro[di,:] = score_function_multid_seperate(self.Z[:,:,rev_t].T,Sxx.T,func_out= True,C=0.001,which=1,l=lnthsc,which_dim=di+1, kern=self.kern)(self.B[:,:,rev_ti].T)\n #print(grad_ln_ro[:,:3]) \n grad_ln_ro = (score_function_multid_seperate_all_dims(self.Z[:,:,rev_t].T,Sxx.T,func_out= True,C=0.001,which=1,l=lnthsc, kern=self.kern)(self.B[:,:,rev_ti].T) ).T \n #np.testing.assert_allclose(grad_ln_ro2.T, grad_ln_ro)\n return grad_ln_ro \n\n\n def bw_density_estimation(self, ti, rev_ti):\n #grad_ln_b = np.zeros((self.dim,self.N))\n lnthsc = 2*np.std(self.B[:,:,rev_ti],axis=1)\n #print(ti, rev_ti, rev_ti-1)\n bnds = np.zeros((self.dim,2))\n for ii in range(self.dim):\n bnds[ii] = [max(np.min(self.Z[ii,:,rev_ti-1]),np.min(self.B[ii,:,rev_ti])),min(np.max(self.Z[ii,:,rev_ti-1]),np.max(self.B[ii,:,rev_ti]))]\n #sparse points\n #print(bnds)\n sum_bnds = np.sum(bnds)\n #if np.isnan(sum_bnds) or np.isinf(sum_bnds):\n # plt.figure(figsize=(6,4)),plt.plot(self.B[0].T,self.B[1].T,alpha=0.3);\n # plt.plot(self.y1[0],self.y1[1],'go')\n # plt.plot(self.y2[0],self.y2[1],'ro')\n # plt.show()\n Sxx = np.array([ np.random.uniform(low=bnd[0],high=bnd[1],size=(self.N_sparse)) for bnd in bnds ] )\n \n # for di in range(self.dim): \n # grad_ln_b[di,:] = score_function_multid_seperate(self.B[:,:,rev_ti].T,Sxx.T,func_out= False,C=0.001,which=1,l=lnthsc,which_dim=di+1, kern=self.kern)#(self.B[:,:,-ti].T)\n grad_ln_b = score_function_multid_seperate_all_dims(self.B[:,:,rev_ti].T,Sxx.T,func_out= False,C=0.001,which=1,l=lnthsc, kern=self.kern).T#(self.B[:,:,-ti].T)\n \n return grad_ln_b # this should be function\n \n \n def backward_simulation(self): \n \n for ti,tt in enumerate(self.timegrid[:-1]): \n W = np.ones((N,1))/N \n if ti==0:\n for di in range(self.dim):\n self.B[di,:,-1] = self.Z[di,:,-1]#self.y2[di] \n else:\n #tti = np.where(self.timegrid==tt)[0][0] \n Ti = self.timegrid.size\n rev_ti = Ti- ti \n #print(rev_ti) \n grad_ln_ro = self.density_estimation(ti,rev_ti) #density estimation of forward particles \n \n if ti==1: \n #print(rev_ti,rev_ti-1)\n self.B[:,:,rev_ti-1] = (self.B[:,:,rev_ti] - self.f(self.B[:,:,rev_ti], self.timegrid[rev_ti])*self.dt + self.dt*self.g**2*grad_ln_ro \\\n + (self.g)*np.random.normal(loc = 0.0, scale = np.sqrt(self.dt),size=(self.dim,self.N)) )% (2*np.pi)\n else:\n grad_ln_b = self.bw_density_estimation(ti,rev_ti)\n self.B[:,:,rev_ti-1] = (self.B[:,:,rev_ti] -\\\n ( self.f(self.B[:,:,rev_ti], self.timegrid[rev_ti])- self.g**2*grad_ln_ro +0.5*self.g**2 * grad_ln_b )*self.dt)% (2*np.pi)\n \n \n \n for di in range(self.dim):\n self.B[di,:,0] = self.y1[di]\n \n return 0 \n\n\n\n def reject_trajectories(self):\n fplus = self.y1+self.f(self.y1,self.t1)*self.dt+4*self.g**2 *np.sqrt(self.dt)\n fminus = self.y1+self.f(self.y1,self.t1) *self.dt-4*self.g**2 *np.sqrt(self.dt)\n for iii in range(2):\n if fplus[iii] < fminus[iii]:\n temp = fminus[iii]\n fminus[iii] = fplus[iii]\n fplus[iii] = temp\n\n sinx = np.where( np.logical_or(np.logical_not(np.logical_and( self.B[0,:,1]fminus[0])) , np.logical_not( np.logical_and(self.B[0,:,1]fminus[0])) ) )[0]\n #((self.B[1,:,-2]fminus[1]) ) ))[0]\n print(sinx)\n temp = len(sinx)\n print(\"Identified %d invalid bridge trajectories \"%len(sinx))\n if self.reject:\n print(\"Deleting invalid trajectories...\")\n sinx = sinx[::-1]\n for element in sinx:\n self.B = np.delete(self.B, element, axis=1)\n return 0\n\n def calculate_u(self,grid_x,ti):\n \"\"\"\n \n\n Parameters\n ----------\n grid_x : array of size d x number of points on the grid\n ti : time index in timegrid for the computation of u\n Computes the control u on the grid or on a the point .\n \n\n Returns\n -------\n The control u(grid_x, t), where t=timegrid[ti].\n\n \"\"\"\n #a = 0.001\n #grad_dirac = lambda x,di: - 2*(x[di] -self.y2[di])*np.exp(- (1/a**2)* (x[0]- self.y2[0])**2)/(a**3 *np.sqrt(np.pi)) \n \n lnthsc = 2*np.std(self.B[:,:,ti],axis=1)\n \n bnds = np.zeros((self.dim,2))\n for ii in range(self.dim):\n bnds[ii] = [max(np.min(self.Z[ii,:,ti]),np.min(self.B[ii,:,ti])),min(np.max(self.Z[ii,:,ti]),np.max(self.B[ii,:,ti]))]\n \n #sum_bnds = np.sum(bnds)\n \n Sxx = np.array([ np.random.uniform(low=bnd[0],high=bnd[1],size=(self.N_sparse)) for bnd in bnds ] )\n #u_t = np.zeros(grid_x.shape)\n # for di in range(self.dim): \n # u_t[di] = score_function_multid_seperate(self.B[:,:,ti].T,Sxx.T,func_out= True,C=0.001,which=1,l=lnthsc,which_dim=di+1, kern=self.kern)(grid_x.T) \\\n # - score_function_multid_seperate(self.Z[:,:,ti].T,Sxx.T,func_out= True,C=0.001,which=1,l=lnthsc,which_dim=di+1, kern=self.kern)(grid_x.T)\n \n u_t = (score_function_multid_seperate_all_dims(self.B[:,:,ti].T,Sxx.T,func_out= True,C=0.001,which=1,l=lnthsc, kern=self.kern)(grid_x.T) \\\n - score_function_multid_seperate_all_dims(self.Z[:,:,ti].T,Sxx.T,func_out= True,C=0.001,which=1,l=lnthsc, kern=self.kern)(grid_x.T) ).T\n #np.testing.assert_allclose(u_t2, u_t)\n return u_t\n \n \n def check_if_covered(self, X, ti):\n \"\"\"\n Checks if test point X falls within forward and backward densities at timepoint timegrid[ti]\n\n Parameters\n ----------\n X : TYPE\n DESCRIPTION.\n ti : TYPE\n DESCRIPTION.\n\n Returns\n -------\n Boolean variable - True if the text point X falls within the densities.\n\n \"\"\"\n covered = True\n bnds = np.zeros((self.dim,2))\n for ii in range(self.dim):\n bnds[ii] = [max(np.min(self.Z[ii,:,ti]),np.min(self.B[ii,:,ti])),min(np.max(self.Z[ii,:,ti]),np.max(self.B[ii,:,ti]))]\n #bnds[ii] = [np.min(self.B[ii,:,ti]),np.max(self.B[ii,:,ti])]\n \n covered = covered * ( (X[ii] >= bnds[ii][0]) and (X[ii] <= bnds[ii][1]) )\n \n return covered\n\n\n\n\nscript_name = sys.argv[0] \ndim = int(sys.argv[1]) #number of oscillators\ninput_num = int(sys.argv[2]) ##task id\nnoise = int(sys.argv[3])\nrepetition = int(sys.argv[4]) ##here I repeat from outside different instances with same parameters\nfre = int(sys.argv[5])\n\nfreqs = [0.25, 0.5, 1]\nfreq = freqs[fre-1]\n\nKs = np.linspace(0,10,11) \n\nK = Ks[input_num-1]\n\n\n\n\n\n\n\n\ndef fsquare(x): ##apply dot product for every 2dim vector in the array f^2\n x = np.atleast_2d( x )\n #return np.apply_along_axis(lambda xi: np.dot(f(xi),f(xi.T) ) , 1, x)\n return np.linalg.norm(f(x.T),axis=0)** 2\n\n\n\n###order parameter\ndef R_t(x):\n return (1/dim)* np.abs( np.sum(np.exp( 1j* x)) )\n\n@numba.njit(parallel=False,fastmath=True)\ndef U(x): #constraint\n ### 1- order parameter\n #difs = (x-x[:,None] ) %(2*np.pi)\n #print(difs.shape)\n #return -(1/dim)* np.array([np.sum(np.triu(difs[:,:,ii])) for ii in range(difs.shape[-1]) ])\n return -1*(1-(1/dim)* np.abs( np.sum(np.exp( 1j* x),axis=0) ) )\n\n\n\n\n\nh = 0.001 #sim_prec\nt_start = 0.\n\n \nt1 = 0\nt2 = 0.5\nT = t2-t1\ntimegrid = np.arange(0,T+h/2,h)\nif dim==6:\n N = 3000#0#0#0#0#2000\nelif dim==10:\n N=4000\nif noise==0:\n g = 0.5\nelif noise==1:\n g = 1\nk = timegrid.size\nM = 300#0\n\n\ny2 = np.ones(dim) \n\nsigmas = np.array([0.1,0.5,1])\nrep_bridge = 1 #10 different bridge instances for every setting\nreps = 20 ##instanses for stochastic path evaluation of each bridge\n \n### first run in kura had scale 0.5\n### in kura_wide0 I have scale 1\nrandom.seed(input_num*100+repetition)\ny0 = np.random.normal( loc=3, scale=0.5,size=dim ) \nws = np.ones(dim)\nws[:round(dim/2)] = freq\nws[round(dim/2):] = -freq\n\n\n\ndef f(x,t=0):\n #print(x.shape)\n difs = -np.sin(x-x[:,None] )\n #print(difs)\n #print( (np.ones((1,dim)) @ difs).shape )\n dn = x.shape ## shape of input state\n if len(dn) ==1: ##if function is called for sinle point, set ni to 1\n \n xout = ws + (K/dim)* np.sum( difs, axis=0)\n else:\n ni = dn[1] ## else if functionn is called for an ensemble of points, set ni to that number \n \n xout = np.tile(ws,(ni,1)).T + (K/dim)* np.sum( difs, axis=0)\n \n return xout\n\nsave_dir = '/work/maoutsa/kura_wide0/'\nnaming = 'Delta%d_N_syst_Kuramoto_inh_k_%d_gi_%d_N_%d_M_%d_repetition_%d_fr_%.3f'%(dim, K,noise,N, M,repetition, freq)\nbridg2d = BRIDGE_ND_reweight(t1,t2,y0,y2,f,g,N,k,M,reweight=True, U=U,dens_est='nonparametric',reject=False,plotting=True,kern='RBF',dt=h)\n\nFnon = np.zeros((dim,timegrid.size,reps)) * np.nan\nFcont = np.zeros((dim,timegrid.size,reps)) * np.nan\nused_us = np.zeros((dim,timegrid.size,reps))\nRttcont = np.zeros((timegrid.size,reps))\nRttnon = np.zeros((timegrid.size,reps))\nfor repi in range(reps):\n with open(save_dir+naming+\"output.txt\", \"a\") as fi:\n print('>>>>>>>>>>>>>>>>>>>',repi, file=fi)\n for ti,tt in enumerate(timegrid[:-1]):\n ### ti is local time, tti is global time - both are time indices\n tti = ti ## index of timepoint in the initial timegrid- the real time axis\n #print(tti)\n if ti==0:\n Fcont[:,tti,repi] = y0 \n Fnon[:,tti,repi] = y0 \n else: \n ###use previous grad log for current step\n if ti== bridg2d.timegrid.size-1: ##this is the timegrid_sub size\n uu = bridg2d.calculate_u(np.atleast_2d(Fcont[:,tti-1,repi]).T,ti-1).reshape(-1)\n else:\n uu = bridg2d.calculate_u(np.atleast_2d(Fcont[:,tti-1,repi]).T,ti).reshape(-1)\n \n used_us[:,tti,repi] = uu#.T\n Fcont[:,tti,repi] = ( Fcont[:,tti-1,repi]+ h* f(Fcont[:,tti-1,repi])+h*g**2 *uu+(g)*np.random.normal(loc = 0.0, scale = np.sqrt(h),size=(dim)) )% (2*np.pi)\n Fnon[:,tti,repi] = ( Fnon[:,tti-1,repi]+ h* f(Fnon[:,tti-1,repi])+(g)*np.random.normal(loc = 0.0, scale = np.sqrt(h),size=(dim)) )% (2*np.pi)\n Rttcont[ti,repi] = R_t(Fcont[:,ti,repi]) \n Rttnon[ti,repi] = R_t(Fnon[:,ti,repi])\n \n with open(save_dir+naming+\"output.txt\", \"a\") as fi:\n print(Rttcont[:,repi], file=fi)\n \n \n \nto_save = dict()\nto_save['Fcont'] = Fcont\nto_save['Rttcont'] = Rttcont\nto_save['Fnon'] = Fnon\nto_save['Rttnon'] = Rttnon\nto_save['timegrid'] = timegrid\nto_save['B'] = bridg2d.B\nto_save['Z'] = bridg2d.Z\nto_save['K'] = K\nto_save['ws'] = ws\nto_save['g'] = g\nto_save['N'] = N\nto_save['M'] = M\nto_save['used_us'] = used_us\nto_save['y0'] = y0\nto_save['repetition'] = repetition\nto_save['freq'] = freq\npickle.dump(to_save, open(save_dir+naming+'.dat', \"wb\"))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.15" } }, "nbformat": 4, "nbformat_minor": 0 }