<jupyter_start><jupyter_text>1.3) Take vectors and matrices from numpy to torch  Deep Prior Distribution of Relaxation Times  In this tutorial we will reproduce Figure 2 in <NAME>., & <NAME>. (2020). The Deep-Prior Distribution of Relaxation Times. Journal of The Electrochemical Society, 167(2), 026506 https://iopscience.iop.org/article/10.1149/1945-7111/ab631a/metaThe DP-DRT method is our next newly developed deep learning based approach to obtain the DRT from the EIS data. The DP-DRT is trained on a single electrochemical impedance spectrum. A single random input is given to the nerural network underlying the DP-DRT.<jupyter_code>import math
import os
import random as rnd
from math import cos, pi, sin

import compute_DRT
import matplotlib.pyplot as plt
import numpy as np
import torch
import torch.nn.functional as F

%matplotlib inline

# check the device
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
print("Using device:", device)
if device.type == "cuda":
    print(torch.cuda.get_device_name(0))
    print("Memory Usage:")
    print("Allocated:", round(torch.cuda.memory_allocated(0) / 1024 ** 2, 1), "MB")
    print("Cached:   ", round(torch.cuda.memory_cached(0) / 1024 ** 2, 1), "MB")

# we will assume you have a cpu
# if you want to use a GPU, you will need to use cuda<jupyter_output><empty_output><jupyter_text>1) Problem setup 1.1) Generate a single stochastic experiment  note: the exact circuit is a ZARCThe impedance of a ZARC can be written as$$Z^{\rm exact}(f) = R_\infty + \displaystyle \frac{1}{\displaystyle \frac{1}{R_{\rm ct}}+C \left(i 2\pi f\right)^\phi}$$where $\displaystyle C = \frac{\tau_0^\phi}{R_{\rm ct}}$.The analytical DRT can be computed analytically as$$\gamma(\log \tau) =  \displaystyle \frac{\displaystyle R_{\rm ct}}{\displaystyle 2\pi} \displaystyle \frac{\displaystyle \sin\left((1-\phi)\pi\right)}{\displaystyle \cosh(\phi \log(\tau/\tau_0))-\cos(\pi(1-\phi))}$$<jupyter_code># set the seed for the random number generators
rng = rnd.seed(214975)
rng_np = np.random.seed(213912)
torch.manual_seed(213912)

# define frequency range, from 1E-4 to 1E4 with 10 ppd
N_freqs = 81
freq_vec = np.logspace(-4.0, 4.0, num=N_freqs, endpoint=True)
tau_vec = 1.0 / freq_vec

# define parameters for ZARC model and calculate the impedance and gamma following the above equations
R_inf = 10
R_ct = 50
phi = 0.8
tau_0 = 1
C = tau_0 ** phi / R_ct

# exact Z and gamma
Z = R_inf + 1.0 / (1.0 / R_ct + C * (1j * 2.0 * pi * freq_vec) ** phi)
gamma_exact = (
    (R_ct)
    / (2.0 * pi)
    * sin((1.0 - phi) * pi)
    / (np.cosh(phi * np.log(tau_vec / tau_0)) - cos((1.0 - phi) * pi))
)

# adding noise to the impedance data
sigma_n_exp = 0.1
Z_exp = Z + sigma_n_exp * (
    np.random.normal(0, 1, N_freqs) + 1j * np.random.normal(0, 1, N_freqs)
)<jupyter_output><empty_output><jupyter_text>1.2) Build $\mathbf A_{\rm re}$ and $\mathbf A_{\rm im}$ matrices<jupyter_code># define the matrices that calculate the impedace from DRT, i.e., Z_re = A_re * gamma, Z_im = A_im * gamma
A_re = compute_DRT.A_re(freq_vec)
A_im = compute_DRT.A_im(freq_vec)
# size of the arbitrary zeta input
N_zeta = 1

# define the neural network
# N is batch size, D_in is input dimension, H is hidden dimension, D_out is output dimension.
N = 1
D_in = N_zeta
H = max(N_freqs, 10 * N_zeta)
# the output also includes the R_inf, so it has dimension N_freq+1
# note that
# 1) there is no inductance (in this specific example - the DP-DRT can include inductive features, see article)
# 2) R_inf is stored as the last item in the NN output

D_out = N_freqs + 1

# Construct the neural network structure
class vanilla_model(torch.nn.Module):
    def __init__(self):
        super(vanilla_model, self).__init__()
        self.fct_1 = torch.nn.Linear(D_in, H)
        self.fct_2 = torch.nn.Linear(H, H)
        self.fct_3 = torch.nn.Linear(H, H)
        self.fct_4 = torch.nn.Linear(H, D_out)

        # initialize the weight parameters
        torch.nn.init.zeros_(self.fct_1.weight)
        torch.nn.init.zeros_(self.fct_2.weight)
        torch.nn.init.zeros_(self.fct_3.weight)
        torch.nn.init.zeros_(self.fct_4.weight)

    # forward
    def forward(self, zeta):
        h = F.elu(self.fct_1(zeta))
        h = F.elu(self.fct_2(h))
        h = F.elu(self.fct_3(h))
        gamma_pred = F.softplus(self.fct_4(h), beta=5)

        return gamma_pred
# transform impedance variables to tensors
Z_exp_re_torch = (
    torch.from_numpy(np.real(Z_exp)).type(torch.FloatTensor).reshape(1, N_freqs)
)
Z_exp_im_torch = (
    torch.from_numpy(np.imag(Z_exp)).type(torch.FloatTensor).reshape(1, N_freqs)
)
# tranform gamma
gamma_exact_torch = torch.from_numpy(gamma_exact).type(torch.FloatTensor)

# transform these matrices into tensors
A_re_torch = torch.from_numpy(A_re.T).type(torch.FloatTensor)
A_im_torch = torch.from_numpy(A_im.T).type(torch.FloatTensor)<jupyter_output><empty_output><jupyter_text>2) Setup DP-DRT model 2.1) Deep network   2.2) Loss function<jupyter_code>def loss_fn(output, Z_exp_re_torch, Z_exp_im_torch, A_re_torch, A_im_torch):

    # we assume no inductance and the R_inf is stored as the last item in the NN output

    MSE_re = torch.sum(
        (output[:, -1] + torch.mm(output[:, 0:-1], A_re_torch) - Z_exp_re_torch) ** 2
    )
    MSE_im = torch.sum((torch.mm(output[:, 0:-1], A_im_torch) - Z_exp_im_torch) ** 2)
    MSE = MSE_re + MSE_im

    return MSE<jupyter_output><empty_output><jupyter_text>3) Train the model<jupyter_code>model = vanilla_model()

# initialize following variables
zeta = torch.randn(N, N_zeta)
loss_vec = np.array([])
distance_vec = np.array([])
lambda_vec = np.array([])

# optimize the neural network
learning_rate = 1e-5
optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)

# max iterations
max_iters = 100001
gamma_NN_store = torch.zeros((max_iters, N_freqs))
R_inf_NN_store = torch.zeros((max_iters, 1))

for t in range(max_iters):
    # Forward pass: compute predicted y by passing x to the model.
    gamma = model(zeta)

    # Compute the loss
    loss = loss_fn(gamma, Z_exp_re_torch, Z_exp_im_torch, A_re_torch, A_im_torch)
    # save it
    loss_vec = np.append(loss_vec, loss.item())

    # store gamma
    gamma_NN = gamma[:, 0:-1].detach().reshape(-1)
    gamma_NN_store[t, :] = gamma_NN

    # store R_inf
    R_inf_NN_store[t, :] = gamma[:, -1].detach().reshape(-1)

    # Compute the distance
    distance = math.sqrt(torch.sum((gamma_NN - gamma_exact_torch) ** 2).item())
    # save it
    distance_vec = np.append(distance_vec, distance)

    # and print it
    if not t % 100:
        print("iter=", t, "; loss=", loss.item(), "; distance=", distance)

    # zero all gradients (purge any cache)
    optimizer.zero_grad()

    # compute the gradient of the loss with respect to model parameters
    loss.backward()

    # Update the optimizer
    optimizer.step()<jupyter_output><empty_output><jupyter_text>4) Analyze results 4.1) Find early stopping value<jupyter_code>index_opt = np.argmin(distance_vec)
index_early_stop = np.flatnonzero(np.abs(np.diff(loss_vec)) < 1e-8)

gamma_DIP_torch_opt = gamma_NN_store[index_opt, :]
R_inf_DIP_torch_opt = R_inf_NN_store[index_opt, :]

gamma_DIP_opt = gamma_DIP_torch_opt.detach().numpy()
R_DIP_opt = R_inf_DIP_torch_opt.detach().numpy()

if len(index_early_stop):
    gamma_DIP_torch_early_stop = gamma_NN_store[index_early_stop[0], :]
    gamma_DIP = gamma_DIP_torch_early_stop.detach().numpy()
    R_DIP = R_inf_NN_store[index_early_stop[0], :]
    R_DIP = R_DIP.detach().numpy()
else:
    gamma_DIP = gamma_DIP_opt
    R_DIP = R_DIP_opt<jupyter_output><empty_output><jupyter_text>4.2) Plot the loss<jupyter_code>plt.semilogy(loss_vec, linewidth=4, color="black")
plt.semilogy(
    np.array([index_early_stop[0], index_early_stop[0]]),
    np.array([1e-3, 1e7]),
    ":",
    linewidth=3,
    color="red",
)
plt.semilogy(
    np.array([index_opt, index_opt]),
    np.array([1e-3, 1e7]),
    ":",
    linewidth=3,
    color="blue",
)
plt.text(
    30000,
    1e2,
    r"early stop",
    {
        "color": "red",
        "fontsize": 20,
        "ha": "center",
        "va": "center",
        "rotation": 90,
        "bbox": dict(boxstyle="round", fc="white", ec="red", pad=0.2),
    },
)
plt.text(
    0.93e5,
    1e2,
    r"optimal",
    {
        "color": "blue",
        "fontsize": 20,
        "ha": "center",
        "va": "center",
        "rotation": 90,
        "bbox": dict(boxstyle="round", fc="white", ec="blue", pad=0.2),
    },
)
plt.rc("text", usetex=False)
plt.rc("font", family="serif", size=15)
plt.rc("xtick", labelsize=15)
plt.rc("ytick", labelsize=15)
plt.xlabel(r"iter", fontsize=20)
plt.ylabel(r"loss", fontsize=20)
plt.axis([0, 1.01e5, 0.9e-2, 1.1e6])
fig = plt.gcf()
fig.set_size_inches(5, 4)
plt.show()<jupyter_output><empty_output><jupyter_text>4.3) Plot the error curve vs. iterationThe error is defined as the distance between predicted DRT and exact DRT, i.e.,$ \rm error = ||\mathbf \gamma_{\rm exact} - \mathbf \gamma_{\rm DP-DRT}||$<jupyter_code>plt.semilogy(distance_vec, linewidth=4, color="black")
plt.semilogy(
    np.array([index_early_stop[0], index_early_stop[0]]),
    np.array([1e-3, 1e7]),
    ":",
    linewidth=4,
    color="red",
)
plt.semilogy(
    np.array([index_opt, index_opt]),
    np.array([1e-3, 1e7]),
    ":",
    linewidth=4,
    color="blue",
)
plt.text(
    30000,
    2e1,
    r"early stop",
    {
        "color": "red",
        "fontsize": 20,
        "ha": "center",
        "va": "center",
        "rotation": 90,
        "bbox": dict(boxstyle="round", fc="white", ec="red", pad=0.2),
    },
)
plt.text(
    0.93e5,
    2e1,
    r"optimal",
    {
        "color": "blue",
        "fontsize": 20,
        "ha": "center",
        "va": "center",
        "rotation": 90,
        "bbox": dict(boxstyle="round", fc="white", ec="blue", pad=0.2),
    },
)
plt.rc("text", usetex=False)
plt.rc("font", family="serif", size=15)
plt.rc("xtick", labelsize=15)
plt.rc("ytick", labelsize=15)
plt.xlabel(r"iter", fontsize=20)
plt.ylabel(r"error", fontsize=20)
plt.axis([0, 1.01e5, 0.9e0, 1.1e2])
fig = plt.gcf()
fig.set_size_inches(5, 4)
plt.show()<jupyter_output><empty_output><jupyter_text>4.4) Plot the impedanceWe compare the DP-DRT EIS spectrum against the one from the stochastic experiment<jupyter_code>Z_DIP = R_DIP + np.matmul(A_re, gamma_DIP) + 1j * np.matmul(A_im, gamma_DIP)

plt.plot(
    np.real(Z_exp),
    -np.imag(Z_exp),
    "o",
    markersize=10,
    color="black",
    label="synth exp",
)
plt.plot(np.real(Z_DIP), -np.imag(Z_DIP), linewidth=4, color="red", label="DP-DRT")
plt.rc("text", usetex=False)
plt.rc("font", family="serif", size=20)
plt.annotate(
    r"$10^{-2}$",
    xy=(np.real(Z_exp[20]), -np.imag(Z_exp[20])),
    xytext=(np.real(Z_exp[20]) - 2, 10 - np.imag(Z_exp[20])),
    arrowprops=dict(arrowstyle="-", connectionstyle="arc"),
)
plt.annotate(
    r"$10^{-1}$",
    xy=(np.real(Z_exp[30]), -np.imag(Z_exp[30])),
    xytext=(np.real(Z_exp[30]) - 2, 6 - np.imag(Z_exp[30])),
    arrowprops=dict(arrowstyle="-", connectionstyle="arc"),
)
plt.annotate(
    r"$1$",
    xy=(np.real(Z_exp[40]), -np.imag(Z_exp[40])),
    xytext=(np.real(Z_exp[40]), 10 - np.imag(Z_exp[40])),
    arrowprops=dict(arrowstyle="-", connectionstyle="arc"),
)
plt.annotate(
    r"$10$",
    xy=(np.real(Z_exp[50]), -np.imag(Z_exp[50])),
    xytext=(np.real(Z_exp[50]) - 1, 10 - np.imag(Z_exp[50])),
    arrowprops=dict(arrowstyle="-", connectionstyle="arc"),
)
plt.rc("xtick", labelsize=15)
plt.rc("ytick", labelsize=15)
plt.legend(frameon=False, fontsize=15)
plt.xlim(10, 65)
plt.ylim(0, 55)
plt.xticks(range(0, 70, 10))
plt.yticks(range(0, 60, 10))
plt.gca().set_aspect("equal", adjustable="box")
plt.xlabel(r"$Z_{\rm re}/\Omega$", fontsize=20)
plt.ylabel(r"$-Z_{\rm im}/\Omega$", fontsize=20)
fig = plt.gcf()
size = fig.get_size_inches()
plt.show()<jupyter_output><empty_output><jupyter_text>4.5) Plot the DRTWe compare the $\gamma$ from the DP-DRT model against the exact one<jupyter_code>plt.semilogx(tau_vec, gamma_exact, linewidth=4, color="black", label="exact")
plt.semilogx(tau_vec, gamma_DIP, linewidth=4, color="red", label="early stop")
plt.semilogx(
    tau_vec, gamma_DIP_opt, linestyle="None", marker="o", color="blue", label="optimal"
)
plt.rc("text", usetex=False)
plt.rc("font", family="serif", size=15)
plt.rc("xtick", labelsize=15)
plt.rc("ytick", labelsize=15)
plt.axis([1e-4, 1e4, -0.4, 25])
plt.legend(frameon=False, fontsize=15)
plt.xlabel(r"$\tau/{\rm s}$", fontsize=20)
plt.ylabel(r"$\gamma/\Omega$", fontsize=20)
fig = plt.gcf()
fig.set_size_inches(5, 4)
plt.show()<jupyter_output><empty_output><jupyter_text>4.6) Ancillary data<jupyter_code>print("total number parameters = ", compute_DRT.count_parameters(model))
print("distance_early_stop = ", distance_vec[index_early_stop[0]])
print("distance_opt= ", distance_vec[index_opt])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1619674670797831, 'fraction_numerical': 0.032997676219984506, 'mean_word_length': 4.093096646942801, 'pattern_counts': {'":': 29, '<': 48, '<?xml version=': 0, '>': 47, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '96982', 'license': 'MIT', 'n_tokens_mistral': 5387, 'n_tokens_neox': 5044, 'n_words': 1072, 'path': 'src/elchempy/experiments/_EIS/DP-DRT-master/example.ipynb', 'repo_name': 'MyPyDavid/ECpy', 'size': '12913'}