The TransportMap support team will be happy to help you out with questions ranging from theory on transport maps, to the installation and usage of the software TransportMaps.

I got an error "log_pdf() got an unexpected keyword argument 'cache'" when I ran the tutorial

0 votes
Hi,

I ran the "Laplace approximation" in tutorial but I I got an error "log_pdf() got an unexpected keyword argument 'cache'".
asked Aug 1, 2018 in usage by Aaron (12 points)

1 Answer

0 votes
 
Best answer

It appears to be a leftover bug in the documentation. Thanks for point it out.

The functions log_pdf, grad_x_log_pdf and hess_x_log_pdf now can take additional arguments and the Gumbel example provided in the tutorial is still coding an old interface. The correct code for the definition of the GumbelDistribution is as follows.

class GumbelDistribution(DIST.Distribution):
    def __init__(self, mu, beta):
        super(GumbelDistribution,self).__init__(1)
        self.mu = mu
        self.beta = beta
        self.dist = stats.gumbel_r(loc=mu, scale=beta)
    def pdf(self, x, params=None, *args, **kwargs):
        return self.dist.pdf(x).flatten()
    def log_pdf(self, x, params=None, *args, **kwargs):
        return self.dist.logpdf(x).flatten()
    def grad_x_log_pdf(self, x, params=None, *args, **kwargs):
        m = self.mu
        b = self.beta
        z = (x-m)/b
        return (np.exp(-z)-1.)/b
    def hess_x_log_pdf(self, x, params=None, *args, **kwargs):
        m = self.mu
        b = self.beta
        z = (x-m)/b
        return (-np.exp(-z)/b**2.)[:,:,np.newaxis]

I'm in the process of updating the tutorial in any place where this bug appear.

Thanks again,

 Daniele

answered Aug 1, 2018 by dabi (307 points)
selected Mar 9, 2019 by dabi
Thanks. It works now.

Best,

Aaron
Hi  Daniele,

I want to use the tutorial ' Direct transports from densities (inference)' to solve the simple groundwater flow problem. I have defined the posterior distribution and prior distribution of this problem without using the packages in the the software Transportmaps now.  Can I then construct transport maps for sampling this inverse problem directly by using the tutorial ' Direct transports from densities (inference)' like 'Linear span parametrization' and 'Sampling from the approximation'? I would very appreciate if you can give me some advice.

Best,

Aaron

Hi Aaron, yes in principle once you have defined the unnormalized density of the posterior distribution (i.e. defining at least the functions log_pdf and grad_x_log_pdf) you can construct a transport map between a standard normal and this posterior. If your actual code is not in the "TransportMaps format" you just need to define a distribution that links your old code to the TransportMaps code.

I would recommend to use the IntegratedSquared parametrization as it enforces monotonicity explicitly and thus performs often better than the LinearSpan parametrization (enforcing monotonicity only locally).

Also, if the problem you are tackling is a (high-dimensional) spatial distribution I would recommend to massage the problem first. Creating a very high dimensional map is a non-trivial task. Often we first pullback the distribution through the linear map corresponding to the Laplace approximation. Then we employ techniques such as Likelihood Informed Subspaces or Active Subspaces to reduce the dimensionality of the problem prior the construction of the transport map.

Hi  Daniele, thanks for you help. But I still have some questions to solve this problem.

1. How to define a distribution that links my old code to the TransportMaps code? My problem considers a linear elliptic PDE on the unit interval \(S=[0,1]\) which is a bit different with the examples in the tutorial.

$$ \begin{cases}
	\nabla\cdot(k(s)\nabla u) = -f(s)\\
	u|_{s=0}=0, u|_{s=1}=0
\end{cases} $$
\[ f(s) = a\exp(-\frac{1}{2b^2}(s-s_m)^2) \]
with \(a=100\), \(s_m=0.5\), and \(b=0.01\).
This is my old code and the density of the posterior distribution is 'kl_log_post' in the code.
import numpy as np 


def solvepde(k):
    n = np.size(k)
    
    h = 1/(n-1)
    s = np.linspace(0,1,n)
    b = h**2 * 100 * np.exp(-(s[1:n-1]-0.5)**2/(2*0.01**2))
    A = np.zeros((n-2,n-2))
    for i in range(n-2):
        A[i,i] = -(2*k[i+1]+k[i]+k[i+2])/2
    for i in range(n-3):
        A[i,i+1] = (k[i+1]+k[i+2])/2
        A[i+1,i] = A[i,i+1]
    tmp = np.linalg.lstsq(A,b)[0]
    u = np.zeros((n,))
    u[1:n-1] = tmp
    return u

def matern(n):
    sigma = 1.5
    lc = 0.5
    M = np.zeros((n,n))
    for i in list(range(n)):
        for j in list(range(n)):
            M[i,j] = np.abs(i-j)/(n-1)    
    return sigma**2 * np.exp(-M/lc)

def phi(u):
    n = np.size(u)
    num = 31
    result = 0
    for i in range(num):
        result = result + (true_u[1+i*(n//num)]-u[1+i*(n//num)])**2
    return result/(2*0.05**2)

def log_prior(k,M):
    k = np.mat(k).T
    n = np.size(k)
    return -0.5*k.T*np.linalg.inv(M)*k-n/2*np.log(2*np.pi)-0.5*np.log(np.linalg.det(M))

def prior(k,M):
    k = np.mat(k).T
    n = np.size(k)
    return np.exp(-0.5*k.T*np.linalg.inv(M)*k) / ( (2*np.pi)**(n/2)*np.sqrt(np.linalg.det(M)) )

def post(k):
    u = solvepde(np.exp(k)+0.5)
    return np.exp(-phi(u))*prior(k)

def log_post(k,M):
    u = solvepde(np.exp(k)+0.5)
    return log_prior(k,M)-phi(u)

def kl_log_prior(k,T):
    k = np.mat(k).T
    n = np.size(k)
    return -0.5*k.T*np.linalg.inv(T)*k-n/2*np.log(2*np.pi)-0.5*np.log(np.linalg.det(T))

def kl_log_post(x,M,A):
    k = np.dot(A,x)
    u = solvepde(np.exp(k)+0.5)
    return kl_log_prior(k,M)-phi(u)


dim = 101
kl = 10
k = np.random.randn(dim)
true_u = solvepde(k)

n = np.size(k)
M = matern(n)
a,e = np.linalg.eig(M)
A = e[:,0:kl]
T = np.diag(a[0:kl])
 
2.  Is grad_x_log_pdf necessary for Integrated Squared parametrization and optimization?
 Is it possible to set parameter 'ders = 0' in the KL-divergence minimization problem to avoid this? 
Since grad_x_log_pdf is difficult to gain in this problem, right? 
3. This is a PDE inverse problem and I need to retain about 10 Karhuen-Loeve modes. 
Besides using the likelihood informed subspaces, will the method "low-dimensional couplings" 
be useful for this problem? I also see this method in the tutorial.
Thanks again,
Aaron
Hi  Daniele,

I still haven't figured out  how to define a distribution that links my old code to the TransportMaps code. Can you help me out of the difficulty?

I really appreciate your assistance.

Best,

Aaron

Hi Aaron, 

I would write the distribution somewhat like in the following (I have not tested the code because I'm on vacation right now, but it should be along these lines):

import TransportMaps.Distributions as DIST

class EllipticPDEDistribution(DIST.Distribution):

  def __init__(self, n):
    self.M = self.matern(n)
    a, e = np.linalg.eig(M)
    self.A = e[:,0:kl]
    self.T = np.diag(a[0:kl])

  def solvepde(self, k):
    n = np.size(k)
    
    h = 1/(n-1)
    s = np.linspace(0,1,n)
    b = h**2 * 100 * np.exp(-(s[1:n-1]-0.5)**2/(2*0.01**2))
    A = np.zeros((n-2,n-2))
    for i in range(n-2):
        A[i,i] = -(2*k[i+1]+k[i]+k[i+2])/2
    for i in range(n-3):
        A[i,i+1] = (k[i+1]+k[i+2])/2
        A[i+1,i] = A[i,i+1]
    tmp = np.linalg.lstsq(A,b)[0]
    u = np.zeros((n,))
    u[1:n-1] = tmp
    return u

  def matern(self, n):
    sigma = 1.5
    lc = 0.5
    M = np.zeros((n,n))
    for i in list(range(n)):
        for j in list(range(n)):
            M[i,j] = np.abs(i-j)/(n-1)    
    return sigma**2 * np.exp(-M/lc)

  def phi(self, u):
    n = np.size(u)
    num = 31
    result = 0
    for i in range(num):
        result = result + (true_u[1+i*(n//num)]-u[1+i*(n//num)])**2
    return result/(2*0.05**2)

  def log_prior(self, k,M):
    k = np.mat(k).T
    n = np.size(k)
    return -0.5*k.T*np.linalg.inv(M)*k-n/2*np.log(2*np.pi)-0.5*np.log(np.linalg.det(M))

  def prior(self, k,M):
    k = np.mat(k).T
    n = np.size(k)
    return np.exp(-0.5*k.T*np.linalg.inv(M)*k) / ( (2*np.pi)**(n/2)*np.sqrt(np.linalg.det(M)) )

  def post(self, k):
    u = solvepde(np.exp(k)+0.5)
    return np.exp(-phi(u))*prior(k)

  def log_post(self, k,M):
    u = solvepde(np.exp(k)+0.5)
    return log_prior(k,M)-phi(u)

  def kl_log_prior(self, k,T):
    k = np.mat(k).T
    n = np.size(k)
    return -0.5*k.T*np.linalg.inv(T)*k-n/2*np.log(2*np.pi)-0.5*np.log(np.linalg.det(T))

  def kl_log_post(self, x,M,A):
    k = np.dot(A,x)
    u = solvepde(np.exp(k)+0.5)
    return kl_log_prior(k,M)-phi(u)

  def log_pdf(self, x, *args, **kwargs):
    return self.kl_log_post(x, self.M, self.A)

and then initiate it by

dim = 101
kl = 10
k = np.random.randn(dim)
true_u = solvepde(k)
n = np.size(k)
pi = EllipticPDEDistribution(n)

2. grad_x_log_pdf is not mandatory. If you set ders=0, then it will use finite difference to approximate it at each iteration. If you are in moderate dimension (say 10), it should work smoothly. If you go higher in dimension I recommend to use analitic gradients (e.g. using the adjoint method).

3. In 10 dimension you should be able to build a linear map (maybe a second order one, I don't know the number of coefficients over the top of my head). Consider that the more complexity you require, the more number of points you need to use to discretize the KL-divergence. The "low-dimensional couplings" technique is not helpfull here. It is useful when Markov structure is present, but for PDE this is a little tricky, in particular when using KL-expansion which then epresses the random field in global modes.

Since you are using the KL-expansion, it may be the case that most of the complexity is in the interactions between the leading modes. Therefore you could impose some sparsity in the map. Something like:

\[ T(x) = \left[ \begin{array}{l} T_1(x_1) \\ T_2(x_1,x_2) \\ \quad\vdots \\ T_k(x_1,\ldots ,x_k) \\ T_{k+1}(x_{k+1}) \\ \quad\vdots \\ T_n(x_n) \end{array}\right] \]

You can look at this example on how to define a map with such a sparsity structure.

Hope this helps.

Hi Daniele,

Thanks for your kind help!

I link my old code to the TransportMaps code by your distribution. But I got an error 'ValueError: shapes (101,1) and (6,1) not aligned: 1 (dim 1) != 6 (dim 0)' when I ran the code. I have tried to figure out the bug but it still can't be solved now. Can you help me out of this difficulty?

This is my python code

import numpy as np 
import TransportMaps.Distributions as DIST
import matplotlib as mpl
import matplotlib.pyplot as plt
import TransportMaps as TM
import TransportMaps.Maps as MAPS

def solvepde(k):
    n = np.size(k)
    
    h = 1/(n-1)
    s = np.linspace(0,1,n)
    b = h**2 * 100 * np.exp(-(s[1:n-1]-0.5)**2/(2*0.01**2))
    A = np.zeros((n-2,n-2))
    for i in range(n-2):
        A[i,i] = -(2*k[i+1]+k[i]+k[i+2])/2
    for i in range(n-3):
        A[i,i+1] = (k[i+1]+k[i+2])/2
        A[i+1,i] = A[i,i+1]
    tmp = np.linalg.lstsq(A,b,rcond=None)[0]
    u = np.zeros((n,))
    u[1:n-1] = tmp
    return u

def matern(n):
    sigma = 1.5
    lc = 0.5
    M = np.zeros((n,n))
    for i in list(range(n)):
        for j in list(range(n)):
            M[i,j] = np.abs(i-j)/(n-1)    
    return sigma**2 * np.exp(-M/lc)

def phi(u):
    n = np.size(u)
    num = 31
    result = 0
    for i in range(num):
        result = result + (true_u[1+i*(n//num)]-u[1+i*(n//num)])**2
    return result/(2*0.05**2)

def log_prior(k,M):
    k = np.mat(k).T
    n = np.size(k)
    return -0.5*k.T*np.linalg.inv(M)*k-n/2*np.log(2*np.pi)-0.5*np.log(np.linalg.det(M))

def prior(k,M):
    k = np.mat(k).T
    n = np.size(k)
    return np.exp(-0.5*k.T*np.linalg.inv(M)*k) / ( (2*np.pi)**(n/2)*np.sqrt(np.linalg.det(M)) )

def post(k):
    u = solvepde(np.exp(k)+0.5)
    return np.exp(-phi(u))*prior(k,M)

def log_post(k,M):
    u = solvepde(np.exp(k)+0.5)
    return log_prior(k,M)-phi(u)

def kl_log_prior(x,T):
    x = np.mat(x).T
    n = np.size(x)
    return -0.5*x.T*np.linalg.inv(T)*x-n/2*np.log(2*np.pi)-0.5*np.log(np.linalg.det(T))

def kl_log_post(x,M,A):
    k = np.dot(A,x)
    u = solvepde(np.exp(k)+0.5)
    return kl_log_prior(k,M)-phi(u)



class EllipticPDEDistribution(DIST.Distribution):

  def __init__(self, n):
    self.M = self.matern(n)
    a, e = np.linalg.eig(M)
    self.A = e[:,0:kl]
    self.T = np.diag(a[0:kl])
    self.dim = kl


  def solvepde(self, k):
    n = np.size(k)
    
    h = 1/(n-1)
    s = np.linspace(0,1,n)
    b = h**2 * 100 * np.exp(-(s[1:n-1]-0.5)**2/(2*0.01**2))
    A = np.zeros((n-2,n-2))
    for i in range(n-2):
        A[i,i] = -(2*k[i+1]+k[i]+k[i+2])/2
    for i in range(n-3):
        A[i,i+1] = (k[i+1]+k[i+2])/2
        A[i+1,i] = A[i,i+1]
    tmp = np.linalg.lstsq(A,b)[0]
    u = np.zeros((n,))
    u[1:n-1] = tmp
    return u

  def matern(self, n):
    sigma = 1.5
    lc = 0.5
    M = np.zeros((n,n))
    for i in list(range(n)):
        for j in list(range(n)):
            M[i,j] = np.abs(i-j)/(n-1)    
    return sigma**2 * np.exp(-M/lc)

  def phi(self, u):
    n = np.size(u)
    num = 31
    result = 0
    for i in range(num):
        result = result + (true_u[1+i*(n//num)]-u[1+i*(n//num)])**2
    return result/(2*0.05**2)

  def log_prior(self, k,M):
    k = np.mat(k).T
    n = np.size(k)
    return -0.5*k.T*np.linalg.inv(M)*k-n/2*np.log(2*np.pi)-0.5*np.log(np.linalg.det(M))

  def prior(self, k,M):
    k = np.mat(k).T
    n = np.size(k)
    return np.exp(-0.5*k.T*np.linalg.inv(M)*k) / ( (2*np.pi)**(n/2)*np.sqrt(np.linalg.det(M)) )

  def post(self, k):
    u = solvepde(np.exp(k)+0.5)
    return np.exp(-phi(u))*prior(k)

  def log_post(self, k,M):
    u = solvepde(np.exp(k)+0.5)
    return log_prior(k,M)-phi(u)

  def kl_log_prior(self, k,T):
    k = np.mat(k).T
    n = np.size(k)
    return -0.5*k.T*np.linalg.inv(T)*k-n/2*np.log(2*np.pi)-0.5*np.log(np.linalg.det(T))

  def kl_log_post(self, x,M,A):
    k = np.dot(A,x)
    u = solvepde(np.exp(k)+0.5)
    return kl_log_prior(k,M)-phi(u)

  def log_pdf(self, x, *args, **kwargs):
    return self.kl_log_post(x, self.M, self.A)


n = 101
kl = 1

M = matern(n)
a,e = np.linalg.eig(M)
A = e[:,0:kl]
T = np.diag(a[0:kl])

np.random.seed(4)

k = np.random.multivariate_normal(np.zeros(n),M)
true_u = solvepde(np.exp(k)+0.5)

pi = EllipticPDEDistribution(n)

dim = pi.dim
qtype = 3          # Gauss quadrature
qparams = [5]*dim  # Quadrature order
reg = None         # No regularization
tol = 1e-3         # Optimization tolerance
ders = 0           # Use gradient and Hessian

order = 2
T2 = TM.Default_IsotropicIntegratedSquaredTriangularTransportMap(dim, order, 'total')
rho = DIST.StandardNormalDistribution( dim )
push_T2_rho = DIST.PushForwardTransportMapDistribution(T2, rho)
log = push_T2_rho.minimize_kl_divergence(pi, qtype=qtype, qparams=qparams, regularization=reg, tol=tol, ders=ders)

 

I really appreciate your assistance.

Best,

Aaron

 

This is the whole Error information:

  File "<ipython-input-1-b29f62a2d2f0>", line 1, in <module>
    runfile('/home/Aaron/EllipticPDE.py', wdir='/home/Aaron')

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/spyder/utils/site/sitecustomize.py", line 705, in runfile
    execfile(filename, namespace)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/spyder/utils/site/sitecustomize.py", line 102, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "/home/Aaron/EllipticPDE.py", line 174, in <module>
    log = push_T2_rho.minimize_kl_divergence(pi, qtype=qtype, qparams=qparams, regularization=reg, tol=tol, ders=ders)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Distributions/TransportMapDistributions.py", line 527, in minimize_kl_divergence
    mpi_pool=mpi_pool, grad_check=grad_check, hess_check=hess_check)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Maps/TriangularTransportMapBase.py", line 1044, in minimize_kl_divergence
    grad_check=grad_check, hess_check=hess_check)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Maps/TriangularTransportMapBase.py", line 1168, in minimize_kl_divergence_complete
    options=options, callback=self.minimize_kl_divergence_callback)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/scipy/optimize/_minimize.py", line 597, in minimize
    return _minimize_bfgs(fun, x0, args, jac, callback, **options)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 963, in _minimize_bfgs
    gfk = myfprime(x0)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 293, in function_wrapper
    return function(*(wrapper_args + args))

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 723, in approx_fprime
    return _approx_fprime_helper(xk, f, epsilon, args=args)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 657, in _approx_fprime_helper
    f0 = f(*((xk,) + args))

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 293, in function_wrapper
    return function(*(wrapper_args + args))

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Maps/TransportMapBase.py", line 1377, in minimize_kl_divergence_objective
    mpi_pool=self.mpi_pool)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Misc.py", line 286, in mpi_map
    fval = f(**args)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Routines.py", line 148, in kl_divergence
    mpi_pool=mpi_pool_tuple[1])

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Misc.py", line 286, in mpi_map
    fval = f(**args)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Misc.py", line 654, in wrapped
    out = f(slf, *args, **kwargs)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Misc.py", line 105, in wrapped
    out = f(slf, *args, **kwargs)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Distributions/TransportMapDistributions.py", line 587, in log_pdf
    cache=cache)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Misc.py", line 654, in wrapped
    out = f(slf, *args, **kwargs)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Misc.py", line 105, in wrapped
    out = f(slf, *args, **kwargs)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Maps/MapBase.py", line 256, in log_pullback
    lpdf = pi.log_pdf(ev, params=params_pi, idxs_slice=idxs_slice, cache=pi_cache)

  File "/home/Aaron/EllipticPDE.py", line 145, in log_pdf
    return self.kl_log_post(x, self.M, self.A)

  File "/home/Aaron/EllipticPDE.py", line 140, in kl_log_post
    k = np.dot(A,x)

ValueError: shapes (101,1) and (6,1) not aligned: 1 (dim 1) != 6 (dim 0)

The method log_pdf (grad_x_log_pdf and hess_x_log_pdf) take the argument x which is \(m \times d \) where \( m \) is the number of points used to discretize the KL-divergence and \( d \) is the number of dimensions (in your case one).

I think your code should be

k = np.dot(A, x.T).T

so that \( k \) is now \( m \times n \) where n is the number of discretization points.

The method log_pdf should return a 1-dimensional array with \( m \) entries.

Thanks, it woks now...but has a new strange error 'ValueError: setting an array element with a sequence.' Do you know how it produces this error?

This is the whole error information:

 File "<ipython-input-2-b29f62a2d2f0>", line 1, in <module>
    runfile('/home/Aaron/EllipticPDE.py', wdir='/home/Aaron')

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/spyder/utils/site/sitecustomize.py", line 705, in runfile
    execfile(filename, namespace)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/spyder/utils/site/sitecustomize.py", line 102, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "/home/Aaron/EllipticPDE.py", line 174, in <module>
    log = push_T2_rho.minimize_kl_divergence(pi, qtype=qtype, qparams=qparams, regularization=reg, tol=tol, ders=ders)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Distributions/TransportMapDistributions.py", line 527, in minimize_kl_divergence
    mpi_pool=mpi_pool, grad_check=grad_check, hess_check=hess_check)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Maps/TriangularTransportMapBase.py", line 1044, in minimize_kl_divergence
    grad_check=grad_check, hess_check=hess_check)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Maps/TriangularTransportMapBase.py", line 1168, in minimize_kl_divergence_complete
    options=options, callback=self.minimize_kl_divergence_callback)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/scipy/optimize/_minimize.py", line 597, in minimize
    return _minimize_bfgs(fun, x0, args, jac, callback, **options)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 963, in _minimize_bfgs
    gfk = myfprime(x0)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 293, in function_wrapper
    return function(*(wrapper_args + args))

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 723, in approx_fprime
    return _approx_fprime_helper(xk, f, epsilon, args=args)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 657, in _approx_fprime_helper
    f0 = f(*((xk,) + args))

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 293, in function_wrapper
    return function(*(wrapper_args + args))

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Maps/TransportMapBase.py", line 1377, in minimize_kl_divergence_objective
    mpi_pool=self.mpi_pool)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Misc.py", line 286, in mpi_map
    fval = f(**args)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Routines.py", line 148, in kl_divergence
    mpi_pool=mpi_pool_tuple[1])

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Misc.py", line 286, in mpi_map
    fval = f(**args)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Misc.py", line 654, in wrapped
    out = f(slf, *args, **kwargs)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Misc.py", line 105, in wrapped
    out = f(slf, *args, **kwargs)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Distributions/TransportMapDistributions.py", line 587, in log_pdf
    cache=cache)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Misc.py", line 654, in wrapped
    out = f(slf, *args, **kwargs)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Misc.py", line 105, in wrapped
    out = f(slf, *args, **kwargs)

  File "/home/Aaron/anaconda3/lib/python3.6/site-packages/TransportMaps/Maps/MapBase.py", line 256, in log_pullback
    lpdf = pi.log_pdf(ev, params=params_pi, idxs_slice=idxs_slice, cache=pi_cache)

  File "/home/Aaron/EllipticPDE.py", line 145, in log_pdf
    return self.kl_log_post(x, self.M, self.A)

  File "/home/Aaron/EllipticPDE.py", line 141, in kl_log_post
    u = solvepde(np.exp(k)+0.5)

  File "/home/Aaron/EllipticPDE.py", line 16, in solvepde
    A[i,i] = -(2*k[i+1]+k[i]+k[i+2])/2

ValueError: setting an array element with a sequence.

This I don't know. It seems be an error happening in solvepde.

Maybe it has something to do with the dimensions of the input.

I would recommend debugging using, e.g., pudb.

Thanks, I have tried to figure out the bug but it still can't be solved now.  When I remove the class 'EllipticPDEDistribution(DIST.Distribution)', there is no error. I think the error comes from the class (or input to class) but I still don't know how to solve it. Do you have other ideas?

This code without class 'EllipticPDEDistribution(DIST.Distribution)' runs well.

import numpy as np 
import TransportMaps.Distributions as DIST
import matplotlib as mpl
import matplotlib.pyplot as plt
import TransportMaps as TM
import TransportMaps.Maps as MAPS

def solvepde(k):
    n = np.size(k)
    
    h = 1/(n-1)
    s = np.linspace(0,1,n)
    b = h**2 * 100 * np.exp(-(s[1:n-1]-0.5)**2/(2*0.01**2))
    A = np.zeros((n-2,n-2))
    for i in range(n-2):
        A[i,i] = -(2*k[i+1]+k[i]+k[i+2])/2
    for i in range(n-3):
        A[i,i+1] = (k[i+1]+k[i+2])/2
        A[i+1,i] = A[i,i+1]
    tmp = np.linalg.lstsq(A,b,rcond=None)[0]
    u = np.zeros((n,))
    u[1:n-1] = tmp
    return u

def matern(n):
    sigma = 1.5
    lc = 0.5
    M = np.zeros((n,n))
    for i in list(range(n)):
        for j in list(range(n)):
            M[i,j] = np.abs(i-j)/(n-1)    
    return sigma**2 * np.exp(-M/lc)

def phi(u):
    n = np.size(u)
    num = 31
    result = 0
    for i in range(num):
        result = result + (true_u[1+i*(n//num)]-u[1+i*(n//num)])**2
    return result/(2*0.05**2)

def log_prior(k,M):
    k = np.mat(k).T
    n = np.size(k)
    return -0.5*k.T*np.linalg.inv(M)*k-n/2*np.log(2*np.pi)-0.5*np.log(np.linalg.det(M))

def prior(k,M):
    k = np.mat(k).T
    n = np.size(k)
    return np.exp(-0.5*k.T*np.linalg.inv(M)*k) / ( (2*np.pi)**(n/2)*np.sqrt(np.linalg.det(M)) )

def post(k):
    u = solvepde(np.exp(k)+0.5)
    return np.exp(-phi(u))*prior(k,M)

def log_post(k,M):
    u = solvepde(np.exp(k)+0.5)
    return log_prior(k,M)-phi(u)

def kl_log_prior(x,T):
    x = np.mat(x).T
    n = np.size(x)
    return -0.5*x.T*np.linalg.inv(T)*x-n/2*np.log(2*np.pi)-0.5*np.log(np.linalg.det(T))

def kl_log_post(x,M,A):
    k = np.dot(A,x.T).T
    u = solvepde(np.exp(k)+0.5)
    return kl_log_prior(k,M)-phi(u)






n = 101
kl = 1

M = matern(n)
a,e = np.linalg.eig(M)
A = e[:,0:kl]
T = np.diag(a[0:kl])

np.random.seed(4)

k_1 = np.random.multivariate_normal(np.zeros(n),M)
true_u = solvepde(np.exp(k_1)+0.5)

 

This code with class 'EllipticPDEDistribution(DIST.Distribution)' has the error.

import numpy as np 
import TransportMaps.Distributions as DIST
import matplotlib as mpl
import matplotlib.pyplot as plt
import TransportMaps as TM
import TransportMaps.Maps as MAPS

def solvepde(k):
    n = np.size(k)
    
    h = 1/(n-1)
    s = np.linspace(0,1,n)
    b = h**2 * 100 * np.exp(-(s[1:n-1]-0.5)**2/(2*0.01**2))
    A = np.zeros((n-2,n-2))
    for i in range(n-2):
        A[i,i] = -(2*k[i+1]+k[i]+k[i+2])/2
    for i in range(n-3):
        A[i,i+1] = (k[i+1]+k[i+2])/2
        A[i+1,i] = A[i,i+1]
    tmp = np.linalg.lstsq(A,b,rcond=None)[0]
    u = np.zeros((n,))
    u[1:n-1] = tmp
    return u

def matern(n):
    sigma = 1.5
    lc = 0.5
    M = np.zeros((n,n))
    for i in list(range(n)):
        for j in list(range(n)):
            M[i,j] = np.abs(i-j)/(n-1)    
    return sigma**2 * np.exp(-M/lc)

def phi(u):
    n = np.size(u)
    num = 31
    result = 0
    for i in range(num):
        result = result + (true_u[1+i*(n//num)]-u[1+i*(n//num)])**2
    return result/(2*0.05**2)

def log_prior(k,M):
    k = np.mat(k).T
    n = np.size(k)
    return -0.5*k.T*np.linalg.inv(M)*k-n/2*np.log(2*np.pi)-0.5*np.log(np.linalg.det(M))

def prior(k,M):
    k = np.mat(k).T
    n = np.size(k)
    return np.exp(-0.5*k.T*np.linalg.inv(M)*k) / ( (2*np.pi)**(n/2)*np.sqrt(np.linalg.det(M)) )

def post(k):
    u = solvepde(np.exp(k)+0.5)
    return np.exp(-phi(u))*prior(k,M)

def log_post(k,M):
    u = solvepde(np.exp(k)+0.5)
    return log_prior(k,M)-phi(u)

def kl_log_prior(x,T):
    x = np.mat(x).T
    n = np.size(x)
    return -0.5*x.T*np.linalg.inv(T)*x-n/2*np.log(2*np.pi)-0.5*np.log(np.linalg.det(T))

def kl_log_post(x,M,A):
    k = np.dot(A,x.T).T
    u = solvepde(np.exp(k)+0.5)
    return kl_log_prior(k,M)-phi(u)



class EllipticPDEDistribution(DIST.Distribution):

  def __init__(self, n):
    self.M = self.matern(n)
    a, e = np.linalg.eig(M)
    self.A = e[:,0:kl]
    self.T = np.diag(a[0:kl])
    self.dim = kl


  def solvepde(self, k):
    n = np.size(k)
    
    h = 1/(n-1)
    s = np.linspace(0,1,n)
    b = h**2 * 100 * np.exp(-(s[1:n-1]-0.5)**2/(2*0.01**2))
    A = np.zeros((n-2,n-2))
    for i in range(n-2):
        A[i,i] = -(2*k[i+1]+k[i]+k[i+2])/2
    for i in range(n-3):
        A[i,i+1] = (k[i+1]+k[i+2])/2
        A[i+1,i] = A[i,i+1]
    tmp = np.linalg.lstsq(A,b)[0]
    u = np.zeros((n,))
    u[1:n-1] = tmp
    return u

  def matern(self, n):
    sigma = 1.5
    lc = 0.5
    M = np.zeros((n,n))
    for i in list(range(n)):
        for j in list(range(n)):
            M[i,j] = np.abs(i-j)/(n-1)    
    return sigma**2 * np.exp(-M/lc)

  def phi(self, u):
    n = np.size(u)
    num = 31
    result = 0
    for i in range(num):
        result = result + (true_u[1+i*(n//num)]-u[1+i*(n//num)])**2
    return result/(2*0.05**2)

  def log_prior(self, k,M):
    k = np.mat(k).T
    n = np.size(k)
    return -0.5*k.T*np.linalg.inv(M)*k-n/2*np.log(2*np.pi)-0.5*np.log(np.linalg.det(M))

  def prior(self, k,M):
    k = np.mat(k).T
    n = np.size(k)
    return np.exp(-0.5*k.T*np.linalg.inv(M)*k) / ( (2*np.pi)**(n/2)*np.sqrt(np.linalg.det(M)) )

  def post(self, k):
    u = solvepde(np.exp(k)+0.5)
    return np.exp(-phi(u))*prior(k)

  def log_post(self, k,M):
    u = solvepde(np.exp(k)+0.5)
    return log_prior(k,M)-phi(u)

  def kl_log_prior(self, k,T):
    k = np.mat(k).T
    n = np.size(k)
    return -0.5*k.T*np.linalg.inv(T)*k-n/2*np.log(2*np.pi)-0.5*np.log(np.linalg.det(T))

  def kl_log_post(self, x,M,A):
    k = np.dot(A,x.T).T
    u = solvepde(np.exp(k)+0.5)
    return kl_log_prior(k,M)-phi(u)

  def log_pdf(self, x, *args, **kwargs):
    return self.kl_log_post(x, self.M, self.A)


n = 101
kl = 2

M = matern(n)
a,e = np.linalg.eig(M)
A = e[:,0:kl]
T = np.diag(a[0:kl])

np.random.seed(4)

k_1 = np.random.multivariate_normal(np.zeros(n),M)
true_u = solvepde(np.exp(k_1)+0.5)

pi = EllipticPDEDistribution(n)

dim = pi.dim
qtype = 3          # Gauss quadrature
qparams = [5]*dim  # Quadrature order
reg = None         # No regularization
tol = 1e-3         # Optimization tolerance
ders = 0           # Use gradient and Hessian

order = 2
T2 = TM.Default_IsotropicIntegratedSquaredTriangularTransportMap(dim, order, 'total')
rho = DIST.StandardNormalDistribution( dim )
push_T2_rho = DIST.PushForwardTransportMapDistribution(T2, rho)
log = push_T2_rho.minimize_kl_divergence(pi, qtype=qtype, qparams=qparams, regularization=reg, tol=tol, ders=ders)

 Really sorry for giving you too much trouble.

Best,

Aaron

Besides,I add &lsquo;self.dim = kl&rsquo; in the class you give me. Since when I run the code without this sentence, there is an error &lsquo;EllipticPDEDistribution&rsquo; object has no attribute &lsquo;_dim&rsquo;. I am not sure if this is right.

On the init error you are right. One way to do the same thing you did is to do the following.

def __init__(self, n, kl):
  self.M = self.matern(n) 
  a, e = np.linalg.eig(M) 
  self.A = e[:,0:kl]
  self.T = np.diag(a[0:kl]) 
  super(EllipticPDEDistribution, self).__init__(kl)

 

Thanks, I understand. Do you have some ideas for the other error?
The problem is that log_pdf allows you to do everything vectorized, i.e. it takes x to be a 2d array of m points in d dimensions (in your case d=kl=1). So, when solvepde is called the argument k is 2d array of dimension m x n. For each row of this array you should call solvepde which I guess it is not vectorized. This is why you see the error. A for loop over the m rows of k should do the trick.
Thank you for your help. It works now.
...