How to define multiple initial conditions in for loop?

[jiweiqi]

Hi,

I am playing with PINN for a toy ODE problem of a reaction network model. It has five state variables, i.e., five species. Therefore, we have to set five initial conditions as

initial_values = [0.5488135, 0.71518937,0,0,0]
ic = [0]*5
ic[0] =  dde.IC(timedomain, lambda x: initial_values[0], boundary, component = 0)
ic[1] =  dde.IC(timedomain, lambda x: initial_values[1], boundary, component = 1)
ic[2] =  dde.IC(timedomain, lambda x: initial_values[2], boundary, component = 2)
ic[3] =  dde.IC(timedomain, lambda x: initial_values[3], boundary, component = 3)
ic[4] =  dde.IC(timedomain, lambda x: initial_values[4], boundary, component = 4)

We try to set ic in a for loop to make life easier, especially when you have hundreds of species. Something like

for i in range(5):
        ic[i] = dde.IC(timedomain, lambda x: initial_values[i], boundary, component = i)

However, this will not work since define a lambda function inside a for loop will make all of the ics being the value of i=5. We have tried the suggestions from https://stackoverflow.com/questions/19837486/python-lambda-in-a-loop.

for i in range(5):
        ic[i] = dde.IC(timedomain, lambda index=i: initial_values[index], boundary, component = i)

But it gives another error TypeError: only integer scalar arrays can be converted to a scalar index.

Any suggestion?

Below is the entire code

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 tensorflow.keras.layers as layers
import tensorflow.keras.activations as activations

import matplotlib.pyplot as plt

initial_values = [0.5488135, 0.71518937,0,0,0]
K = [0.1, 0.2, 0.13, 0.3]
def main():

    def ode_system(x, y):
        """
        dx_1/dt = -2 * k_1 *x_1^2 - k_2 * x_1
        dx_2/dt = -k_1 * x_1^2 - k_4 * x_2 * x_4
        dx_3/dt = -k_2 * x_1 - k_3 * x_3
        dx_4/dt = k_3 * x_3 - k_4 * x_2 * x_4
        dx_5/dt = k_4 * x_2 * x_4
        """

        x_1, x_2, x_3, x_4, x_5 = y[:, 0:1], y[:, 1:2], y[:, 2:3], y[:, 3:4], y[:, 4:5]
        dx1_t = tf.gradients(x_1, x)[0]
        dx2_t = tf.gradients(x_2, x)[0]
        dx3_t = tf.gradients(x_3, x)[0]
        dx4_t = tf.gradients(x_4, x)[0]
        dx5_t = tf.gradients(x_5, x)[0]

        return[-2.0 * K[0] * x_1**2 - K[1] * x_1 - dx1_t,
               K[0] * x_1**2 - K[3] * x_2 * x_4 - dx2_t,
               K[1] * x_1 - K[2] * x_3 - dx3_t,
               K[2] * x_3 - K[3] * x_2 * x_4 - dx4_t,
               K[3] * x_2 * x_4 - dx5_t]

    def boundary(x, on_initial):
        return x==0

    def func(x):
        pass

    def one(x):
        return 1
    def zero(x):
        return 0

    num_pathways = 5
    num_chemicals = 5

    timedomain = dde.geometry.TimeDomain(0,20)
    ic = [0]*5
    # ic[0] =  dde.IC(timedomain, lambda x: initial_values[0], boundary, component = 0)
    # ic[1] =  dde.IC(timedomain, lambda x: initial_values[1], boundary, component = 1)
    # ic[2] =  dde.IC(timedomain, lambda x: initial_values[2], boundary, component = 2)
    # ic[3] =  dde.IC(timedomain, lambda x: initial_values[3], boundary, component = 3)
    # ic[4] =  dde.IC(timedomain, lambda x: initial_values[4], boundary, component = 4)

    for i in range(5):
        ic[i] = dde.IC(timedomain, lambda x: initial_values[i], boundary, component = i)

    print(ic[0].func(0))

    data = dde.data.PDE(timedomain, ode_system, ic, 2500, 2)
    net = dde.maps.FNN([1] + [30]*4 + [5], "tanh", "Glorot normal")
    model = dde.Model(data, net)

    model.compile("adam", lr=1e-3)
    model.train(epochs=1000, display_every = 100)
    model.compile("L-BFGS-B")
    losshistory, train_state = model.train(model_save_path = "CRN model")

    dde.saveplot(losshistory, train_state, issave=True, isplot=True)

if __name__ == "__main__":
    main()

[lululxvi]

The error is due to lambda x: initial_values[0]. Here x is a Nx5 array, and the function should return a Nx1 array. So, use lambda x: np.ones((len(x), 1)) * initial_values[0])

[jiweiqi]

Xiexie Lushen. The solution works as ic = dde.IC(timedomain, lambda x: initial_values, boundary), instead of setting it up component by component. Another quick question is how to choose the number of boundary points, i.e., the last argument in data = dde.data.PDE(timedomain, ode_system, ic, 2500, 2). Is the number of boundary points related to the weighting factor for the (soft) loss at the boundary?

[lululxvi]

Are you sure that ic = dde.IC(timedomain, lambda x: initial_values, boundary) works without modifying the source code?

There is only a time domain with left and right boundary points, i.e., 2 boundary points. This is not related to the weights of the loss.

[jiweiqi]

@lululxvi No, ic = dde.IC(timedomain, lambda x: initial_values, boundary) not works. Do you know why?

For your initial comment of

The error is due to lambda x: initial_values[0]. Here x is a Nx5 array, and the function should return a Nx1 array. So, use lambda x: np.ones((len(x), 1)) * initial_values[0])

lambda x: initial_values[0] works well if you don't use it in for loop. The problem is that you can not put lambda function inside for loop, otherwise, it will be evaluated for every iteration. Such that ic[0], ic[1], ic[2], ic[3] = ic[4] =0. But ic[0] is supposed to be not 0. There are some tricks to get around it. But we didn't succeed.

[lululxvi]

• oh, yes, lambda x: initial_values[0] works without the for loop.
• DeepXDE requires that the IC of each component is defined separately. (Because in many cases, we only have partial observations of the initial conditions.)
• Does the following work?

  ic[i] = dde.IC(timedomain, lambda x, index=i: initial_values[index], boundary, component = i)

[jiweiqi]

Got it. It makes sense to define ic separately.

ic[i] = dde.IC(timedomain, lambda x, index=i: initial_values[index], boundary, component = i)

This one is not working, with the error of TypeError: only integer scalar arrays can be converted to a scalar index

[lululxvi]

What about other solutions in the stackoverflow, e.g.,

def make_func(i):
    initial_values = ...
    return lambda x: initial_values[i]

[jiweiqi]

The make_func works! Cool. I post the code here, in case someone else in interested into it.

def make_func(i):
        return lambda x: initial_values[i]

for i in range(5):
        ic[i] = dde.IC(timedomain, make_func(i), boundary, component = i)

# Print the ic for the first species to see whether it is set correctly.
print(ic[0].func(x=0))