alirezayazdani1 / sbinns Goto Github PK

View Code? Open in Web Editor NEW

36.0 36.0 25.0 103 KB

License: Apache License 2.0

Python 65.49% Jupyter Notebook 34.51%

sbinns's People

Contributors

Stargazers

Watchers

sbinns's Issues

Setting search range for a trainable parameter in inverse ODE problem

Hello @alirezayazdani1 @lululxvi

I have an inverse ODE problem and I am trying to infer the value of one variable (based on Lorenz inverse problem) but I don´t get the true value of A (I created the synthetic data using the true value of A) , I wanted to set a search range like you did in the glucose-insulin code and using this part of the code :
def get_variable(v): low, up = v * 0.2, v * 1.8 l = (up - low) / 2 return l * tf.tanh(tf.Variable(0, trainable=True, dtype=tf.float32)) + l + low

and using deepxde I get this error message :
File "/usr/local/lib/python3.7/dist-packages/keras/optimizer_v2/utils.py", line 80, in filter_empty_gradients
([v.name for v in vars_with_empty_grads]))
File "/usr/local/lib/python3.7/dist-packages/keras/optimizer_v2/utils.py", line 80, in
([v.name for v in vars_with_empty_grads]))

AttributeError: Tensor.name is meaningless when eager execution is enabled.

Do you know why this error comes up ? Thanks in advance !

Errors because of inside function in PointSet class

Dear Prof Lu,

Appreciate your excellent work. I am wondering work like yours for a long time.

DeepXDE is in 0.11 version now. The api for PointSet class has updated, compared with 0.9 version, especially inside function. All the codes will not work for 0.11 version because of changing in inside function.
Maybe, you should mention that.

Thank you.

Running the Code on DeepXDE v 0.14.1

Hi,
I was trying to run your code using DeepXDE v 0.14.1 and replaced dde.bc.PointSet with dde.PointSetBC. However, I am getting the following error while running the code:

TypeError                                 Traceback (most recent call last)
<ipython-input-17-39c3ba0a95b5> in <module>()
      1 # Train
----> 2 var_list = pinn(t, y, noise)

<ipython-input-16-9d5866f87c8f> in pinn(data_t, data_y, noise)
     56         np.random.choice(np.arange(1, n - 1), size=n // 4, replace=False), [0, n - 1]
     57     )
---> 58     ptset = dde.PointSetBC(data_t[idx])
     59     inside = lambda x, _: ptset.inside(x)
     60     observe_y4 = dde.DirichletBC(

TypeError: __init__() missing 1 required positional argument: 'values'

Can you please propose how to solve this issue?

Best,
Akash.

How to setup "ODE_weights" and “data_weights” ？

Hi ,Thank you for your contribyution to this work.
When I run a inverse ODE problem，I am confused about the setting of ODE weights. What rules or experience do you use to determine ODE_weights and data_weights?

I failed to try other models with glycolysis code

I am trying to solve the ODE equation based on a model that modifies the glycolysis reaction, but the loss function keeps reporting nan and I can't get the result, even if I set the loss_weights of the ODE to 0, it doesn't help, the first 15 of the 17 state variables have a loss of nan. I hope someone can help me find the cause of the problem if possible.

The ODE equation of the system is as follows.

The code for data generation is as follows：
Biodiesel_ODE.py

import numpy as np
from scipy.integrate import odeint
from data_config import *
from k_config import *

K = [Vreactf, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k_1, k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9, k_10, Fa, CoAlco, af, Ae]


def enzymatic_biodiesel_model(t, k=None):
    if k is None:
        k = K
    CoAlco = k[22]
    af = k[23]
    Ae = k[24]
    Vreactf = k[0]

    def ode_func(x, t):
        T, D, M, BD, FFA, G, W, CH, E, EX, ET, ED, EM, ECH, Ef, Vp, V = x
        if t <= 120:
            Fa = k[21][0]
        else:
            Fa = k[21][1]
        if x[16] >= Vreactf:
            Fa = 0

        aT = af * (Vp / V)
        Af = (aT - Ae * (E + EX + ET + ED + EM + ECH)) / Ae

        r1 = k[1] * Ef * Af - k[11] * E
        r2 = k[2] * T * E - k[12] * ET
        r3 = k[3] * ET - k[13] * EX * D
        r4 = k[4] * D * E - k[14] * ED
        r5 = k[5] * ED - k[15] * EX * M
        r6 = k[6] * M * E - k[16] * EM
        r7 = k[7] * EM - k[17] * EX * G
        r8 = k[8] * EX * W - k[18] * FFA * E
        r9 = k[9] * EX * CH - k[19] * BD * E
        r10 = k[10] * CH * E - k[20] * ECH

        rg = (r7 * V * 92 / 1261)
        rw = (-r8 * V * 18 / 1000)

        return [
            -T * Fa / V - r2,
            -D * Fa / V + (r3 - r4),
            -M * Fa / V + (r5 - r6),
            -BD * Fa / V + r9,
            -FFA * Fa / V + r8,
            -G * Fa / V + r7,
            -W * Fa / V - r8,
            (CoAlco - CH) * Fa / V - (r9 + r10),
            -E * Fa / V + (r1 - r2 - r4 - r6 + r8 + r9 - r10),
            -EX * Fa / V + (r3 + r5 + r7 - r8 - r9),
            -ET * Fa / V + (r2 - r3),
            -ED * Fa / V + (r4 - r5),
            -EM * Fa / V + (r6 - r7),
            -ECH * Fa / V + r10,
            -Ef * Fa / V - r1,
            rg + rw,
            Fa
        ]

    # Experimental data in [mole/L]
    st = 0
    tmp = data[32][st:]
    tmp[0] = 1e-4

    FAME = tmp
    FFA = data[33][st:]
    TAG = data[34][st:]
    DAG = data[35][st:]
    MAG = data[36][st:]

    # Initial Condition
    T = TAG[0]
    D = DAG[0]
    M = MAG[0]
    B = FAME[0]
    FA = FFA[0]
    G = 1e-6
    W = mH2O
    CH = Alco
    E = 0.0
    EX = 0.0
    ET = 0.0
    ED = 0.0
    EM = 0.0
    ECH = 0.0
    Ef = Enzyme
    Vp = Vpo
    V = Vo
    x0 = [T, D, M, B, FA, G, W, CH, E, EX, ET, ED, EM, ECH, Ef, Vp, V]

    return odeint(ode_func, x0, t)


def main():
    t = np.arange(0, 1500, 1)[:, None]
    np.ravel(t)

    y, info = enzymatic_biodiesel_model(np.ravel(t))
    # y = enzymatic_biodiesel_model(np.ravel(t))
    np.savetxt("enzymatic_biodiesel.dat", np.hstack((t, y)))


if __name__ == "__main__":
    main()

The neural network code is as follows：
Biodisel_NN.py

import torch
import numpy as np
import deepxde as dde
from Biodiesel_ODE import enzymatic_biodiesel_model
from data_config import *


def enzymatic_biodiesel_sbinn(data_t, data_y, noise):
   
    k1_ = dde.Variable(0.0)
    k2_ = dde.Variable(0.0)
    k3_ = dde.Variable(0.0)
    k4_ = dde.Variable(0.0)
    k5_ = dde.Variable(0.0)
    k6_ = dde.Variable(0.0)
    k7_ = dde.Variable(0.0)
    k8_ = dde.Variable(0.0)
    k9_ = dde.Variable(0.0)
    k10_ = dde.Variable(0.0)
    k_1_ = dde.Variable(0.0)
    k_2_ = dde.Variable(0.0)
    k_3_ = dde.Variable(0.0)
    k_4_ = dde.Variable(0.0)
    k_5_ = dde.Variable(0.0)
    k_6_ = dde.Variable(0.0)
    k_7_ = dde.Variable(0.0)
    k_8_ = dde.Variable(0.0)
    k_9_ = dde.Variable(0.0)
    k_10_ = dde.Variable(0.0)
    var_list = [k1_, k2_, k3_, k4_, k5_, k6_, k7_, k8_, k9_, k10_, k_1_, k_2_, k_3_, k_4_, k_5_, k_6_, k_7_, k_8_, k_9_, k_10_]

    def ODE(t, y):
        row, col = y.shape
        T, D, M, BD, FFA, G, W, CH, E, EX, ET, ED, EM, ECH, Ef, Vp, V = [y[:, i:i + 1] for i in range(col)]

        Fa_ = torch.tensor(Fa)[(t > 120).type(torch.long)]
        Fa_ = Fa_ * ((V < Vreactf).type(torch.long).reshape(len(V), -1))

        k1 = torch.nn.functional.softplus(k1_)
        k2 = torch.nn.functional.softplus(k2_)
        k3 = torch.nn.functional.softplus(k3_)
        k4 = torch.nn.functional.softplus(k4_)
        k5 = torch.nn.functional.softplus(k5_)
        k6 = torch.nn.functional.softplus(k6_)
        k7 = torch.nn.functional.softplus(k7_)
        k8 = torch.nn.functional.softplus(k8_)
        k9 = torch.nn.functional.softplus(k9_)
        k10 = torch.nn.functional.softplus(k10_)
        k_1 = torch.nn.functional.softplus(k_1_)
        k_2 = torch.nn.functional.softplus(k_2_)
        k_3 = torch.nn.functional.softplus(k_3_)
        k_4 = torch.nn.functional.softplus(k_4_)
        k_5 = torch.nn.functional.softplus(k_5_)
        k_6 = torch.nn.functional.softplus(k_6_)
        k_7 = torch.nn.functional.softplus(k_7_)
        k_8 = torch.nn.functional.softplus(k_8_)
        k_9 = torch.nn.functional.softplus(k_9_)
        k_10 = torch.nn.functional.softplus(k_10_)

        aT = af * (Vp / V)
        Af = (aT - Ae * (E + EX + ET + ED + EM + ECH)) / Ae

        r1 = k1 * Ef * Af - k_1 * E
        r2 = k2 * T * E - k_2 * ET
        r3 = k3 * ET - k_3 * EX * D
        r4 = k4 * D * E - k_4 * ED
        r5 = k5 * ED - k_5 * EX * M
        r6 = k6 * M * E - k_6 * EM
        r7 = k7 * EM - k_7 * EX * G
        r8 = k8 * EX * W - k_8 * FFA * E
        r9 = k9 * EX * CH - k_9 * BD * E
        r10 = k10 * CH * E - k_10 * ECH

        rg = (r7 * V * 92 / 1261)
        rw = (-r8 * V * 18 / 1000)

        return [
            dde.grad.jacobian(y, t, i=0, j=0) - (-T * Fa_ / V - r2),
            dde.grad.jacobian(y, t, i=1, j=0) - (-D * Fa_ / V + (r3 - r4)),
            dde.grad.jacobian(y, t, i=2, j=0) - (-M * Fa_ / V + (r5 - r6)),
            dde.grad.jacobian(y, t, i=3, j=0) - (-BD * Fa_ / V + r9),
            dde.grad.jacobian(y, t, i=4, j=0) - (-FFA * Fa_ / V + r8),
            dde.grad.jacobian(y, t, i=5, j=0) - (-G * Fa_ / V + r7),
            dde.grad.jacobian(y, t, i=6, j=0) - (-W * Fa_ / V - r8),
            dde.grad.jacobian(y, t, i=7, j=0) - ((CoAlco - CH) * Fa_ / V - (r9 + r10)),
            dde.grad.jacobian(y, t, i=8, j=0) - (-E * Fa_ / V + (r1 - r2 - r4 - r6 + r8 + r9 - r10)),
            dde.grad.jacobian(y, t, i=9, j=0) - (-EX * Fa_ / V + (r3 + r5 + r7 - r8 - r9)),
            dde.grad.jacobian(y, t, i=10, j=0) - (-ET * Fa_ / V + (r2 - r3)),
            dde.grad.jacobian(y, t, i=11, j=0) - (-ED * Fa_ / V + (r4 - r5)),
            dde.grad.jacobian(y, t, i=12, j=0) - (-EM * Fa_ / V + (r6 - r7)),
            dde.grad.jacobian(y, t, i=13, j=0) - (-ECH * Fa_ / V + r10),
            dde.grad.jacobian(y, t, i=14, j=0) - (-Ef * Fa_ / V - r1),
            dde.grad.jacobian(y, t, i=15, j=0) - (rg + rw),
            dde.grad.jacobian(y, t, i=16, j=0) - Fa_
        ]

    
    geom = dde.geometry.TimeDomain(data_t[0, 0], data_t[-1, 0])

    def boundary(x, _):  
        return np.isclose(x[0], data_t[-1, 0])  

    y1 = data_y[-1]
    bc = [dde.DirichletBC(geom, lambda X: y1[i], boundary, component=i) for i in range(0, 17)]

    n = len(data_t)
    idx = np.append(
        np.random.choice(np.arange(1, n - 1), size=n // 4, replace=False), [0, n - 1]
    )
    ic = [dde.PointSetBC(data_t[idx], data_y[idx, i:i + 1], component=0) for i in range(5)]
    np.savetxt("biodiesel_input.data", np.hstack((data_t[idx], data_y[idx, 0:1], data_y[idx, 1:2], data_y[idx, 2:3], data_y[idx, 3:4], data_y[idx, 4:5])))

    data = dde.data.PDE(geom, ODE, bc + ic, anchors=data_t)

    net = dde.maps.FNN([1] + [17] + [128] * 3 + [17], "relu", "Glorot normal")

    # def feature_transform(t):
    #     return torch.concat(
    #         (
    #             t,
    #             torch.exp(-t),
    #             torch.sin(2 * t),
    #             torch.sin(3 * t),
    #             torch.sin(4 * t),
    #             torch.sin(5 * t),
    #             torch.sin(6 * t),
    #         ),
    #         axis=1,
    #     )
    #
    # net.apply_feature_transform(feature_transform)
    # ？？？？？？TODO
    # def output_transform(t, y):
    #     return (
    #             torch.as_tensor(data_y[0]) + torch.tanh(t) * torch.tensor([1, 1, 0.1, 0.1, 0.1, 1, 0.1]) * y
    #     )
    #
    # net.apply_output_transform(output_transform)

    model = dde.Model(data, net)
    checkpointer = dde.callbacks.ModelCheckpoint(
        "./model/model.ckpt", verbose=1, save_better_only=True, period=1000
    )
    variable = dde.callbacks.VariableValue(
        var_list, period=1000, filename="variables.dat", precision=3,
    )
    callbacks = [checkpointer, variable]
    # TODO
    bc_weights = [1] * 17
    if noise >= 0.1:
        bc_weights = [w * 10 for w in bc_weights]
    # TODO
    data_weights = [1] * 5
    # Large noise requires small data_weights
    if noise >= 0.1:
        data_weights = [w / 10 for w in data_weights]
    model.compile("adam", lr=1e-3, loss_weights=[0] * 17 + bc_weights + data_weights, external_trainable_variables=var_list)
    model.train(iterations=1000, display_every=1000)
    # TODO
    ode_weights = [1] * 17
    # Large noise requires large ode_weights
    if noise > 0:
        ode_weights = [10 * w for w in ode_weights]
    model.compile("adam", lr=1e-3, loss_weights=ode_weights + bc_weights + data_weights)

    losshistory, train_state = model.train(
        epochs=900000 if noise == 0 else 2000000,
        display_every=1000,
        callbacks=callbacks,
        disregard_previous_best=True,
        # model_restore_path="./model/model.ckpt-"
    )
    dde.saveplot(losshistory, train_state, issave=True, isplot=True)
    var_list = [model.sess.run(v) for v in var_list]
    return var_list


def main():
    data_set = np.loadtxt("enzymatic_biodiesel.dat")
    t = data_set[:, 0].reshape(1500, -1)
    y = data_set[:, 1:]
    noise = 0
    # Add noise
    # if noise > 0:
    #     std = noise * y.std(0)
    #     y[1:-1, :] += np.random.normal(0, std, (y.shape[0] - 2, y.shape[1]))
    #     np.savetxt("enzymatic_biodiesel_noise.dat", np.hstack((t, y)))

    # Train
    var_list = enzymatic_biodiesel_sbinn(t, y, noise)

    # Prediction
    y = enzymatic_biodiesel_model(np.ravel(t), *var_list)
    np.savetxt("glycolysis_pred.dat", np.hstack((t, y)))


if __name__ == "__main__":
    main()

I can sort out the parameters in k_config and data_config later if needed. One thing to note is that Fa will change after 120 minutes and will equal 0 when Fa>Vreactf, and this mutation may have some impact. I'm not sure.

Recommend Projects

React

A declarative, efficient, and flexible JavaScript library for building user interfaces.
Vue.js

🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
Typescript

TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
TensorFlow

An Open Source Machine Learning Framework for Everyone
Django

The Web framework for perfectionists with deadlines.
Laravel

A PHP framework for web artisans
D3

Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

javascript

JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
web

Some thing interesting about web. New door for the world.
server

A server is a program made to process requests and deliver data to clients.
Machine learning

Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Visualization

Some thing interesting about visualization, use data art
Game

Some thing interesting about game, make everyone happy.

Recommend Org

Facebook

We are working to build community through open source technology. NB: members must have two-factor auth.
Microsoft

Open source projects and samples from Microsoft.
Google

Google ❤️ Open Source for everyone.
Alibaba

Alibaba Open Source for everyone
D3

Data-Driven Documents codes.
Tencent

China tencent open source team.