Inverse Problem: Different output sizes and recovering boundary values

[EShalom]

Hi @lululxvi,

Thank you for your wonderful framework and documentation!

I had a question on the application of the DeepXDE framework for the inverse problem which after consulting many of the FAQs and examples I still am a bit stuck on how to implement correctly.

I’m looking into systems that have spatial variable advection component (f(x)) and volume (v(x)) which acts on the evolution of c(x,t). See equation below:

Where the system has some boundary conditions:

There is measurement data for c(x,t) at discrete spatial and time points (uniformly sampled in both time/space).

The problem I’m looking at is recovering the advection component f(x) and volume component v(x) and the influx into the boundary of the system from the data sampled for c within the system.

Specifically, my questions are:

1. From issue #178 it seems I should be applying apply_output_transform, is this correct? I am unsure how to apply this within the code as the link from the original issue is no longer working.
2. Can DeepXDE function to recover boundary values or boundary influxes (e.g. c per time) or must these be defined for it to work properly?

Currently, I have this which does run but only returns all parameters at all spatial-temporal points and has boundary conditions provided.

import deepxde as dde
import matplotlib.pyplot as plt
import numpy as np
import tensorflow as t
from numpy import heaviside as H

def flow(a,num):
    b = a[0:num+1,0:1]
    return 5 + (0*b)

def volume(a,num):
    b = a[0:num+1,0:1]
    return (0.6*(0.4*np.sin(0.3*b/10)**2+0.6*np.cos(0.15*b/10)**2)+0.3)

def boundary_l(x, on_boundary):
    return on_boundary and np.isclose(x[0], 0)

def boundary_r(x, on_boundary):
    return on_boundary and np.isclose(x[0], Lx)

def C_bcl(x):
    Cl = 4*np.sin(2*np.pi*x[:,1:]/4)*H(x[:,1:],1)*H(-x[:,1:]+2,1)
    return Cl

def C_bcr(x):
    Cr = np.zeros(len(x))
    return Cr

Lx = 256
dx = 8
T = 80
Dt = 2
Nt = int(T/Dt)
Nx = int(Lx/dx) # number of sampling points from the data
xvals = np.linspace(dx/2, Lx-dx/2, Nx) # spatial points
y = np.zeros((3, xvals.size)) # network output arguement

v =(0.6*(0.4*np.sin(0.3*xvals/10)**2+0.6*np.cos(0.15*xvals/10)**2)+0.3)
f = np.ones((Nx+1))*5 # mL/s/mm^2
t=np.linspace(0,T,Nt+1)
Cp = 4*np.sin(2*np.pi*t/4)*H(t,1)*H(-t+2,1)
Cn = np.zeros(Nt+1)
Jp = Cp*((f[0]*H(f[0],0))/dx)
Jn = Cn*((-f[-1]*H(-f[-1],0))/dx)
Jp *= H(f[0],0)
Jn *= H(-f[-1],0)

C_tiss = fmod(Lx, dx, t, 0.2, f, D, v, Jp, Jn
C, tsamp = downsampler(C_tiss, T, 0.2, Dt, 0, 0)# this a function above  just produce the measrement
training_data = np.load("machine_learning_1C.npz")
t, x = training_data["t"], training_data["x"] #1D arrays
X, T = np.meshgrid(x, t)
X = np.reshape(X, (-1, 1))
T = np.reshape(T, (-1, 1))
C = np.reshape(C, (-1, 1))
observe_x, ob_C = np.hstack((X, T)), C

geom = dde.geometry.Interval(0, Lx)# creates domain mesh in space
timedomain = dde.geometry.TimeDomain(0, 80)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

def pde(x, y):
    C, f, v = y[:, 0:1], y[:, 1:2], y[:, 2:3] # separates coefficients in network output
    
    dC_x = dde.grad.jacobian(y, x, i=0, j=0)
    df_x = dde.grad.jacobian(y, x, i=1, j=0)
    dC_t = dde.grad.jacobian(y, x, i=0, j=1)
    
    return v*dC_t + f * dC_x + C*df_x


def fun_init(x):
    return 0 

observe_C = dde.icbc.PointSetBC(observe_x, ob_C, component=0)

bc_l = dde.icbc.DirichletBC(geomtime, C_bcl ,boundary_l, component=1)
ic = dde.icbc.IC(geomtime, fun_init, lambda _, on_initial: on_initial, component=0)

data = dde.data.TimePDE(
    geomtime,
    pde,
    [bc_l, ic, observe_C],
    num_domain=50,
    num_boundary=2,
    num_initial=50,
    anchors=observe_x,
    num_test=1000
)

net = dde.nn.PFNN([2]+[20, 20,20]*2 + [3], "tanh", "Glorot uniform")
model = dde.Model(data, net)
model.compile("adam", lr=1e-3)

losshistory, train_state = model.train(epochs=5000)
#%%
space_num = 50
x = geomtime.uniform_points(space_num*41)
yhat = model.predict(x)
Chat, fhat, vhat = yhat[:, 0:1], yhat[:, 1:2], yhat[:, 2:3]
#%%
ftrue = flow(x,space_num)
vtrue = volume(x,space_num)

Do you think this type of problem is compatible with DeepXDE?

If you think it is possible would you be able to provide some suggestions or direction for implementation? Any tips for this would be greatly appreciated.

Thank you very much!

Cheers, Eve

[lululxvi]

Your problem seems a standard inverse problem. Why do you think you need apply_output_transform?

You don't have to define any BC for inverse problems.

[EShalom]

Hi @lululxvi ,

Thank you for your response!

I just thought the output transform might be needed as the output fields are all different shapes e.g. C(x,t) has space and time points but f(x) and v(x) only have spatial points and any recovered BCs have time points. So I am unsure of how to write them into the output (y) so that yhat = model.predict(x) has the right shape for each type. As it seems to give errors if any y column has a different length?

Cheers, Eve

[lululxvi]

• For the boundary condition f: You can ignore the BC when you solve the problem. After you get c, then you can compute f.
• To implement v as only a function of x, but not t, see FAQ Q: A standard network doesn’t work for me, and I want to implement a network with a special structure/property.

[EShalom]

Thanks a lot for the help on this! I have put a script together - based on the FAQs you mention - which I've included here but it's having some issues I can't seem to solve. The pde has form:

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import numpy as np

import deepxde as dde
from deepxde.backend import tf

import matplotlib.pyplot as plt
 
def ground_truth(x,t,u0,D0,x0):
    C=np.zeros((len(x),len(t)))
    for n in range(0,len(t)):
        for i in range(0,len(x)):
            tv1=1/(2*x0*np.sqrt(D0*np.power(u0,2)*np.pi*t[n]))
            tv2=np.power(((1/u0)*np.log(x[i]/x0)-(1-D0*u0)*t[n]),2)
            tv3=np.exp((-tv2)/(4*D0*t[n]))
            C[i,n]=tv1*tv3
    return C

def pde(x, y):
    C, u, D = y[:, 0:1], y[:, 1:2], y[:, 2:3]
    dC_t = dde.grad.jacobian(y, x, i=0, j=1)
    dC_x = dde.grad.jacobian(y, x, i=0, j=0)
    du_x = dde.grad.jacobian(y, x, i=1, j=0)
    dD_x = dde.grad.jacobian(y, x, i=2, j=0)
    dC_xx = dde.grad.hessian(y, x, i=0, j=0,component=0)
    eq = dC_t + (C*du_x) + (u*dC_x)- (D*dC_xx)-(dC_x*dD_x)
    return eq

Lx = 10
T = 1
u0 = 5
x0 = 3
D0 = 0.02
geom = dde.geometry.Interval(0.01, Lx)
timedomain = dde.geometry.TimeDomain(0.01, T)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

x = np.linspace(0.01, Lx, num=50)
t = np.linspace(0.01, T, num=50)
training_C = ground_truth(x, t, u0, D0, x0)
ob_C = training_C.flatten()[:, None]
X, T = np.meshgrid(x, t)
X = np.reshape(X, (-1, 1))
T = np.reshape(T, (-1, 1))
observe_xt = np.hstack((X, T))
observe_C = dde.icbc.PointSetBC(observe_xt, ob_C,component=0)
#%%
data = dde.data.PDE(
    geomtime,
    pde,
    [observe_C],
    train_distribution='uniform',
    num_domain=700,
    num_boundary=200,
    anchors=observe_xt
)

net = dde.maps.FNN([2] + [150] * 3 + [3], "tanh", "Glorot uniform")

def modify_output(X, y):
    x, t = X[:, 0:1], X[:, 1:2]
    C = y[:,0:1]
    
    # u #
    u_layer_size = 60
    u = tf.layers.dense(x, u_layer_size, tf.nn.tanh)
    u = tf.layers.dense(u, u_layer_size, tf.nn.tanh)
    u = tf.layers.dense(u, u_layer_size, tf.nn.tanh)
    u = tf.layers.dense(u, u_layer_size, tf.nn.tanh)
    u = tf.layers.dense(u, u_layer_size, tf.nn.tanh)
    u_output = tf.layers.dense(u, 1, None)

    # D #
    D_layer_size = 80
    D = tf.layers.dense(x, D_layer_size, tf.nn.tanh)
    D = tf.layers.dense(D, D_layer_size, tf.nn.tanh)
    D = tf.layers.dense(D, D_layer_size, tf.nn.tanh)
    D = tf.layers.dense(D, D_layer_size, tf.nn.tanh)
    D = tf.layers.dense(D, D_layer_size, tf.nn.tanh)
    D = tf.layers.dense(D, D_layer_size, tf.nn.tanh)
    D_output = tf.layers.dense(D, 1, None)
    final_output = tf.concat([C,u_output, D_output], axis=1)
    return final_output

net.apply_output_transform(modify_output)

model = dde.Model(data, net)

model.compile("adam", lr=1e-3, metrics=["l2 relative error"])
losshistory, train_state = model.train(epochs=1000)
model.compile("L-BFGS-B", metrics=["l2 relative error"])
losshistory, train_state = model.train()
dde.saveplot(losshistory, train_state, issave=True, isplot=True)

This returns an error I'm having difficulty debugging...

  File "/anaconda3/envs/tensorflow/lib/python3.7/site-packages/deepxde/utils/internal.py", line 22, in wrapper
    result = f(*args, **kwargs)

  File "/anaconda3/envs/tensorflow/lib/python3.7/site-packages/deepxde/model.py", line 414, in train
    self._test()

  File "/anaconda3/envs/tensorflow/lib/python3.7/site-packages/deepxde/model.py", line 568, in _test
    for m in self.metrics

  File "/anaconda3/envs/tensorflow/lib/python3.7/site-packages/deepxde/model.py", line 568, in <listcomp>
    for m in self.metrics

  File "/anaconda3/envs/tensorflow/lib/python3.7/site-packages/deepxde/metrics.py", line 12, in l2_relative_error
    return np.linalg.norm(y_true - y_pred) / np.linalg.norm(y_true)

TypeError: unsupported operand type(s) for -: 'NoneType' and 'float'

The ytrue value is the None type variable metnioned in the error - ytrue should be from the training data which is defined I think? Unless it comes from the pde definition itself and that is written incorrectly?

If you would let me know what you think could be the problem that would be great.

Many thanks again, Eve

[forxltk]

The metrics in model.compile requires a reference solution in dde.data.PDE.

[EShalom]

Ah, I see - thank you for the comment @forxltk ! I am slightly concerned about this, is this not providing the answer we are hoping to determine via the inverse problem?

Maybe I am misunderstanding and the reference solutions you mention follow the same PDE but have different coefficient values and therefore a separate output from the target data set for the network?

Many thanks, Eve

[forxltk]

Ah, I see - thank you for the comment @forxltk ! I am slightly concerned about this, is this not providing the answer we are hoping to determine via the inverse problem?

Maybe I am misunderstanding and the reference solutions you mention follow the same PDE but have different coefficient values and therefore a separate output from the target data set for the network?

Many thanks, Eve

The metrics and reference solution are not necessary. They just help you evaluate the model during training. Just forget metrics.

model.compile("adam", lr=1e-3)#, metrics=["l2 relative error"])
losshistory, train_state = model.train(epochs=1000)
model.compile("L-BFGS-B")#, metrics=["l2 relative error"])

[EShalom]

Ah, I see the reference is used within the metrics! Thank you very much @forxltk for your quick response, I shall see if this will return something for me :)

[rodrigogdourado]

@EShalom hi! did you manage to solve your problem?