# Fileset

[3D-CGAN.zip](https://mdr.nims.go.jp/filesets/a2174a43-8720-436d-b7c4-a665b640b6f1/download)

## Creator

[Xiaoyang Zheng](https://orcid.org/0000-0003-1452-5855), [Ikumu Watanabe](https://orcid.org/0000-0002-7693-1675)

## Rights

[Creative Commons BY-NC Attribution-NonCommercial 4.0 International](https://creativecommons.org/licenses/by-nc/4.0/)

## Other metadata

[Generating 3D voxelized architected materials using 3D conditional generative adversarial network](https://mdr.nims.go.jp/datasets/0172ab31-42a0-4ff6-980e-4340b4a6c9fc)

## Fulltext

3d_plot_voxel.pyimport matplotlib.pyplot as pltimport numpy as npfrom matplotlib import pyplot as pltfrom mpl_toolkits.mplot3d import Axes3Ddef make_ax(grid=False):    fig = plt.figure()    ax = fig.add_subplot(projection='3d')    # ax.set_xlabel("x")    # ax.set_ylabel("y")    # ax.set_zlabel("z")    ax.grid(grid)    ax.set_axis_off()    return axfilled = np.array([    [[1, 0, 1], [0, 0, 1], [0, 1, 0]],    [[0, 1, 1], [1, 0, 0], [1, 0, 1]],    [[1, 1, 0], [1, 1, 1], [0, 0, 0]]])ax = make_ax()filedirection = "/home/xiaoyang/PycharmProjects/pythonProject30_3d_foam/cgan20220421/fake_image/154.npy"data = np.squeeze(np.load(filedirection)).astype("int32")data_pro = np.reshape(data,[16,-1])data_pro = np.sum(data_pro,-1)/64**3print(data_pro)ax.voxels(data[1], facecolors='#1f77b430', edgecolors='gray')plt.show()CGAN_main.pyimport osos.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'import numpy as npimport tensorflow as tffrom PIL import Imageimport glob, datetimefrom pure_gan_kernel4 import Generator, Discriminatorfrom regression_with_ckpt import Solverfrom load_data import load_dataimport matplotlib.pyplot as pltimport mat73import mathcurrent_time = datetime.datetime.now().strftime("%Y%m%d-%H%M%S")log_dir = 'logs/' + current_timesummary_writer = tf.summary.create_file_writer(log_dir)def new_tri_plot():    lims_rho = [0.2, 0.8]    lims_young = [0, 1]    lims_diffusion = [0.1, 0.4]    fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=[12, 4], constrained_layout=True)    ax1.plot(lims_rho, lims_rho, 'gray')    ax1.set_xlim(lims_rho)    ax1.set_ylim(lims_rho)    ax1.set_aspect(1)    ax1.set_xlabel("Input")    ax1.set_ylabel("Output")    ax1.set_title("rho")    ax2.plot(lims_young, lims_young, 'gray')    ax2.set_xlim(lims_young)    ax2.set_ylim(lims_young)    ax2.set_aspect(1)    ax2.set_xlabel("Input")    ax2.set_ylabel("Output")    ax2.set_title("Young's modulus")    ax3.plot(lims_diffusion, lims_diffusion, 'gray')    ax3.set_xlim(lims_diffusion)    ax3.set_ylim(lims_diffusion)    ax3.set_aspect(1)    ax3.set_xlabel("Input")    ax3.set_ylabel("Output")    ax3.set_title("diffusion")    return ax1, ax2, ax3def tri_plot(y, pred, ax1, ax2, ax3):    x1 = y[:, 0]    x2 = y[:, 1]    x3 = y[:, 2]    y1 = pred[:, 0]    y2 = pred[:, 1]    y3 = pred[:, 2]    ax1.scatter(x1, y1, c='C0', alpha=0.2)    ax2.scatter(x2, y2, c='C1', alpha=0.2)    ax3.scatter(x3, y3, c='C2', alpha=0.2)def celoss_ones(logits):    # [b, 1]    # [b] = [1, 1, 1, 1,]    # Label Smoothing, replace the label with a random number between 0.7 and 1.2    loss = tf.reduce_mean(tf.nn.sigmoid_cross_entropy_with_logits                          (logits=logits, labels=tf.ones_like(logits)-0.3 + np.random.uniform(size=logits.shape) * 0.5))    # loss = tf.keras.losses.categorical_crossentropy(y_pred=logits,    #                                                 y_true=tf.ones_like(logits))    return tf.reduce_mean(loss)def celoss_zeros(logits):    # [b, 1]    # [b] = [1, 1, 1, 1,]    loss = tf.reduce_mean(tf.nn.sigmoid_cross_entropy_with_logits                          (logits=logits, labels=tf.zeros_like(logits)+np.random.uniform(size=logits.shape) * 0.3))    # loss = tf.keras.losses.categorical_crossentropy(y_pred=logits,    #                                                 y_true=tf.zeros_like(logits))    return tf.reduce_mean(loss)def gradient_penalty(discriminator, real_seeds, fake_seeds):    alpha = tf.random.uniform(shape=real_seeds.get_shape(), minval=0., maxval=1.)    differences = fake_seeds - real_seeds  # This is different from MAGAN    interpolates = real_seeds + (alpha * differences)    with tf.GradientTape() as tape:        tape.watch([interpolates])        d_interplote_logits = discriminator(interpolates, training=True)    grads = tape.gradient(d_interplote_logits, interpolates)    # grads:[b, 64, 2] => [b, -1]    grads = tf.reshape(grads, [grads.shape[0], -1])    gp = tf.norm(grads, axis=1)  # [b]    gp = tf.reduce_mean((gp - 1) ** 2)    return gpdef d_loss_fn(generator, discriminator, batch_z, batch_x, condition, is_training):    # 1. treat real image as real    # 2. treat generated image as fake    fake_image = generator(batch_z, condition, is_training)    d_fake_logits = discriminator(fake_image, is_training)    d_real_logits = discriminator(batch_x, is_training)    d_loss_real = celoss_zeros(d_real_logits)    d_loss_fake = celoss_ones(d_fake_logits)    gp = gradient_penalty(discriminator, batch_x, fake_image)    loss = d_loss_real + d_loss_fake + 10 * gp    return loss, gpdef g_loss_fn(generator, discriminator, solver, batch_z, condition, is_training):    fake_image = generator(batch_z, condition, is_training)    d_fake_logits = discriminator(fake_image, is_training)    loss = celoss_zeros(d_fake_logits)    reference = solver(fake_image)    mse = tf.reduce_mean(tf.losses.mean_squared_error(condition, reference))    return loss+0.1*mse, msedef preprocess(x):    x = tf.cast(x, dtype=tf.float32)    return xdef get_data_space(batch_size):    rho = np.random.uniform(low=0.2, high=0.8, size=[batch_size, 1])    young = (250.9/140) * rho**2.845    diffusion = rho * 0.5    return np.concatenate((rho, young, diffusion), axis=-1)def main():    tf.random.set_seed(222)    np.random.seed(222)    assert tf.__version__.startswith('2.')    # hyper parameters    z_dim = 128    epochs = 200    batch_size = 32    learning_rate = 0.0001    is_training = True    generator = Generator()    # generator.load_weights(r'/home/xiaoyang/Desktop/pycharm_project/project10_matrix64/ckpt/20210527-181242_smoothing lable_3rd/generator_133.ckpt')    discriminator = Discriminator()    # discriminator.load_weights(r'/home/xiaoyang/Desktop/pycharm_project/project10_matrix64/ckpt/20210527-181242_smoothing lable_3rd/discriminator_133.ckpt')    solver = Solver()    solver.load_weights(r"ckpt/solver_199.ckpt")    g_optimizer = tf.optimizers.Adam(learning_rate=learning_rate, beta_1=0.5)    d_optimizer = tf.optimizers.Adam(learning_rate=learning_rate, beta_1=0.5)    for epoch in range(epochs):        mat = mat73.loadmat(r'/home/xiaoyang/PycharmProjects/pythonProject30_3d_foam/gan20220415/mat/geo.mat')        data_npy = mat["geo"]        data_npy = np.expand_dims(data_npy,-1)        print(data_npy.shape)        dataset = tf.data.Dataset.from_tensor_slices(data_npy)        dataset = dataset.map(preprocess).shuffle(1000)        dataset = dataset.batch(batch_size)        # train G for 2 times        for step, real_images in enumerate(dataset):            batch = len(real_images)            for _ in range(1):                noise = tf.random.normal([batch, z_dim])  # (b, 128)                condition = get_data_space(batch)                # E and nu => E and G => [-1~1]                with tf.GradientTape() as tape:                    d_loss, gp = d_loss_fn(generator, discriminator, batch_z=noise, batch_x=real_images,                                           condition=condition, is_training=is_training)                if d_loss>0.7:                    grads = tape.gradient(d_loss, discriminator.trainable_variables)                    d_optimizer.apply_gradients(zip(grads, discriminator.trainable_variables))            for _ in range(1):                # E and nu => E and G => [-1~1]                noise = tf.random.normal([batch, z_dim])                condition = get_data_space(batch)                with tf.GradientTape() as tape:                    g_loss, mse = g_loss_fn(generator, discriminator, solver, batch_z=noise,                                            condition=condition, is_training=is_training)                grads = tape.gradient(g_loss, generator.trainable_variables)                g_optimizer.apply_gradients(zip(grads, generator.trainable_variables))            for _ in range(epoch//7+1):                condition = get_data_space(batch)                noise = tf.random.normal([batch, z_dim])                with tf.GradientTape() as tape:                    fake_image = generator(noise, condition, is_training)                    reference = solver(fake_image)                    mse = tf.reduce_mean(tf.losses.mean_squared_error(condition, reference))                grads = tape.gradient(mse, generator.trainable_variables)                g_optimizer.apply_gradients(zip(grads, generator.trainable_variables))        # plot test        ax1, ax2, ax3 = new_tri_plot()        total_sum = 0        total_error = 0        for i in range(10):            noise_test = tf.random.normal([16, z_dim])            condition = get_data_space(16)            fake_image = generator(noise_test, condition, training=False)            pred = solver(fake_image)            loss_test = tf.reduce_mean(tf.losses.mean_squared_error(condition, pred))            total_sum += 1            total_error += loss_test            tri_plot(condition, pred, ax1, ax2, ax3)        # save_generated_voxels        np.save("fake_image/%d.npy" % epoch, np.round(fake_image).astype("bool"))        mse = total_error / total_sum        plt.savefig('results/test/%d_test.png' % epoch)        plt.close()        print(epoch, 'd-loss: ', float(d_loss), 'gp:', float(gp), 'g-loss', float(g_loss), "mse: ", float(mse.numpy()))        # vusulize it on tensorboard: http://localhost:6006/        with summary_writer.as_default():            tf.summary.scalar('d-loss: ', float(d_loss), step=epoch)            tf.summary.scalar('gp: ', float(gp), step=epoch)            tf.summary.scalar('g-loss: ', float(g_loss), step=epoch)            tf.summary.scalar('mse: ', float(mse.numpy()), step=epoch)        # save weights        generator.save_weights('ckpt/generator_%d.ckpt' % (epoch))        discriminator.save_weights('ckpt/discriminator_%d.ckpt' % epoch)if __name__ == '__main__':    main()load_data.pyimport numpy as npimport osdef load_data(loadfile, voxel=66):    files = os.listdir(loadfile)    # print(len(files))    output = []    epoch = 0    for file in files:        # load file        file_dir = loadfile + '/' + file        # num = [int(index) for index in file.split('.') if index.isdigit()][0]        with open(file_dir, "r") as f:            x = np.loadtxt(f, dtype="int")            if len(x) == voxel**3:                x = np.reshape(x, [voxel, voxel, voxel, 1])[1:-1,1:-1,1:-1].tolist()                output.append(x)        if len(output):            np.save(r"/home/xiaoyang/PycharmProjects/pythonProject30_3d_foam/models/voxel/npy_0713_27b/%d.npy" % epoch, output)            epoch += 1            output = []    # return np.array(output)def main():    loadfile = r"/home/xiaoyang/PycharmProjects/pythonProject30_3d_foam/models/voxel/0713_27b"    x = load_data(loadfile)    print(x.shape)if __name__ == "__main__":    main()pure_gan_kernel4.pyimport osos.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'import tensorflow as tffrom tensorflow import kerasfrom tensorflow.keras import layersimport matplotlib.pyplot as pltimport numpy as npdef pad_dim(x, n=1):    x = tf.concat((x[:, :, :, -n:, :], x, x[:, :, :, :n, :]), axis=-2)    x = tf.concat((x[:,:,-n:,:, :], x, x[:,:,:n,:,:]), axis=-3)    x = tf.concat((x[:,-n:,:,:,:], x, x[:,:n,:,:,:]), axis=-4)    return xdef delete_dim(x):    x = x[:, 3:-3, 3:-3, 3:-3, :]    return xclass Generator(keras.Model):    def __init__(self):        super(Generator, self).__init__()        self.fc = layers.Dense(2*2*2*512)        self.bn0 = layers.BatchNormalization(momentum=0.9, epsilon=1e-5)        self.conv1 = layers.Conv3DTranspose(256, kernel_size=[4,4,4], strides=[2,2,2], padding='valid')        self.bn1 = layers.BatchNormalization(momentum=0.9, epsilon=1e-5)        self.conv2 = layers.Conv3DTranspose(128, kernel_size=4, strides=2, padding='valid')        self.bn2 = layers.BatchNormalization(momentum=0.9, epsilon=1e-5)        self.conv3 = layers.Conv3DTranspose(64, kernel_size=4, strides=2, padding='valid')        self.bn3 = layers.BatchNormalization(momentum=0.9, epsilon=1e-5)        self.conv4 = layers.Conv3DTranspose(32, kernel_size=4, strides=2, padding='valid')        self.bn4 = layers.BatchNormalization(momentum=0.9, epsilon=1e-5)        self.conv5 = layers.Conv3DTranspose(16, kernel_size=4, strides=2, padding='valid')        self.bn5 = layers.BatchNormalization(momentum=0.9, epsilon=1e-5)        self.conv6 = layers.Conv3DTranspose(1, kernel_size=4, strides=2, padding='valid')    def call(self, inputs_noise, input_condition, training=None):        # inputs_noise: (b, 128), inputs_condition: (b, 2)        inputs_noise = tf.cast(inputs_noise, dtype=tf.float32)        input_condition = tf.cast(input_condition, dtype=tf.float32)        net = tf.concat((inputs_noise, input_condition), axis=-1)        net = self.fc(net)  # (b, 4*4*512)        net = self.bn0(net, training=training)        net = tf.reshape(net, [-1, 2, 2, 2, 512])  # (b, 4, 4, 512)        net = pad_dim(net, n=1)        net = tf.nn.leaky_relu(self.bn1(self.conv1(net), training=training))        net = delete_dim(net)        net = pad_dim(net)        net = tf.nn.leaky_relu(self.bn2(self.conv2(net), training=training))        net = delete_dim(net)        net = pad_dim(net)        net = tf.nn.leaky_relu(self.bn3(self.conv3(net), training=training))        net = delete_dim(net)        net = pad_dim(net)        net = tf.nn.leaky_relu(self.bn4(self.conv4(net), training=training))        net = delete_dim(net)        # net = pad_dim(net)        # net = tf.nn.leaky_relu(self.bn5(self.conv5(net), training=training))        # net = delete_dim(net)        net = pad_dim(net)        net = self.conv6(net)        net = delete_dim(net)        net = tf.sigmoid(net)  # (b, 256, 256, 1)        return netclass Discriminator(keras.Model):    def __init__(self):        super(Discriminator, self).__init__()        self.conv1 = layers.Conv3D(16, kernel_size=4, strides=2, padding='valid')  # => (b, 128, 128, 16)        self.conv2 = layers.Conv3D(32, kernel_size=4, strides=2, padding='valid')  # => (b, 64, 64, 32)        self.conv3 = layers.Conv3D(64, kernel_size=4, strides=2, padding='valid')  # => (b, 32, 32, 64)        self.conv4 = layers.Conv3D(128, kernel_size=4, strides=2, padding='valid')  # => (b, 16, 16, 128)        self.conv5 = layers.Conv3D(256, kernel_size=4, strides=2, padding='valid')  # => (b, 8, 8, 256)        self.conv6 = layers.Conv3D(512, kernel_size=4, strides=2, padding='valid')  # => (b, 2, 2, 512)        self.flatten = layers.Flatten()        self.fc = layers.Dense(1)    def call(self, inputs_img, training=None):        inputs_img = tf.cast(inputs_img, dtype=tf.float32)        x = pad_dim(inputs_img)        # inputs_img: (b, 256, 256, 1) => (b, 4, 4, 384)        x = layers.Dropout(0.3)(tf.nn.leaky_relu(self.conv1(x)))        x = pad_dim(x)        x = layers.Dropout(0.3)(tf.nn.leaky_relu(self.conv2(x)))        x = pad_dim(x)        x = layers.Dropout(0.3)(tf.nn.leaky_relu(self.conv3(x)))        x = pad_dim(x)        x = layers.Dropout(0.3)(tf.nn.leaky_relu(self.conv4(x)))        x = pad_dim(x)        x = layers.Dropout(0.3)(tf.nn.leaky_relu(self.conv5(x)))        x = pad_dim(x)        x = layers.Dropout(0.3)(tf.nn.leaky_relu(self.conv6(x)))        net = self.flatten(x)  # (b, 4*4*512)        net = self.fc(net)        return netdef main():    g = Generator()    d = Discriminator()    # x = tf.random.uniform([128, 256, 256, 1])    z = tf.random.normal([2, 128])    image_fake = g(z, training=False)    # image_fake = tf.cast((image_fake*127.5)+127.5, dtype=tf.int32)    # image_fake = image_fake.numpy()    # print(image_fake.shape)    sss = d(image_fake)    g.summary()    d.summary()    print(image_fake)    print(sss)    print(image_fake.shape)    # print(image_fake.shape)    # print(image_fake[0])    # print(np.max(image_fake[0]), np.min(image_fake[0]))    # print(np.max(image_fake[1]), np.min(image_fake[1]))    # plt.imshow(image_fake[0], cmap='gray')    # plt.figure()    # plt.imshow(image_fake[1], cmap='gray')    # plt.show()if __name__ == '__main__':    main()readme.docxGenerating 3D voxelized architected materials using 3D conditional generative adversarial networkXiaoyang Zheng 1,2, Ikumu Watanabe 1,21 Center for Basic Research on Materials, National Institute for Materials Science, 1-2-1 Sengen, Tsukuba 305-0047, Japan2 Graduate School of Pure and Applied Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba 305-8573, JapanEmails: ZHENG.Xiaoyang@nims.go.jp (X.Z.); WATANABE.Ikumu@nims.go.jp (I.W.)AbstractThis tutorial aims to give an introduction of how to use a deep generative model, 3D conditional generative adversarial network (3D-CGAN). The 3D-CGAN can be used for the inverse design of 3D voxelized microstructures with target properties. The 3D-CGAN is trained with supervised learning using a labeled dataset. The dataset consists of a large number of geometries (3D arrays) and their corresponding properties (e.g., elastic moduli). After training, the 3D-CGAN can generate a batch of geometries using target properties at inputs. In our previous tutorial, we have demonstrated how to use CGAN for the inverse design of 2D microstructures.[1] This work is based on our previous publication for the inverse design of 3D architected materials.[2] We hope this tutorial can be useful for those who are interested in the inverse design problems of microstructures. 1. IntroductionFigure 1. Deep learning in mechanical metamaterials. In traditional forward-design approach, a structure is created by experienced designers with CAD modeling tools or computational algorithms, and its mechanical properties are then theoretically predicted by FEM simulations and experimentally verified by mechanical testing. On the other hand, deep learning can speed up the forward-design approach by replacing the geometry generation and property prediction processes using ANNs. Moreover, deep learning can invert the forward-design approach with an inverse-design approach, in which a structure can be directly generated using a trained ANN taking target properties as input.[3]Forward design is a conventional approach to design microstructures, such as architected materials, mechanical metamaterials, lattices, etc. This forward design approach follows a general process: a structure is created firstly and then its mechanical properties are investigated by finite element simulation or mechanical testing (Figure 1). The mechanical properties of designed materials will be only known after time-consumingly simulations or experiments. In contrast, deep generative models, such as GAN, enable inverse design of microstructures. In inverse design, microstructures can be automatically generated by inputting target properties to a deep generative model, which outputs corresponding geometries of microstructures. In addition, deep learning can also be used for the property prediction and geometry generation of mechanical metamaterials.[3] In our previous tutorial, we have introduced the basics of CGAN, and how to use the CGAN for the inverse design of 2D auxetic metamaterial. As CGAN and 3D-CGAN share the similar architecture, we here only give a brief introduction of 3D-CGAN. The detailed fundamentals, usage, and guidance can be referred in the previous tutorial.[1] In the 3D-CGAN, it consists of three neural network structures: a generator, a discriminator, and a calculator (Figure 2). The generator is trained to generate the realistic geometries from latent variables (multivariate normal distribution) and user-defined labels (target properties, e.g., elastic moduli, porosity, etc.), and simultaneously aims to deceive the discriminator and the calculator. The discriminator is trained to distinguish geometries produced by the generator from the real dataset. The calculator is trained to predict the properties of a given geometry. Figure 2. Architecture of 3D-CGANThe training process is based on supervised learning, where the dataset consists of paired datapoints. In each datapoint, it is composed of a 3D geometry (with a shape of 64×64×64) and its corresponding properties (with a shape of c, where c is the total numbers of properties). In our case, we used 10,000 datapoints for the training process. The geometries were generated using Voronoi tessellation and the properties were calculated using homogenization method.2. Procedures of CGAN training 2.1 Solver trainingAs the solver is independent of the generator and discriminator, we firstly train the solver with supervised learning. Follow the steps below:a. Open regression_with_ckpt.py with, e.g., PyCharmb. Change the source file in line 22 loaddata(filepath. It returns “geo” with shape of b×64×64×64 and “labels” with shape of b×c. b means the total number of geometries, and 64 means the dimension a geometry voxel. The “geo” consists of 0s and 1s, where 0 represents void part and 1 represents solid part. In “labels”, c means the number of properties for a geometry.Note that the dataset will be divided with 80% for training and 20% for testing. (lines 197-200) The batch size is set to 32, and training epoch is 200.The comparison between the real and predicted values is plotted using line 257.The performance is investigated using the mean square error (MSE) of testing dataset using line 274.The check points (trained weights and biases) are saved using line 285.The check points will be used to train the generator and discriminator.2.2 Generator and discriminator trainingAfter training the solver, we can use the saved weights of the solver to train the generator and discriminator. Follow the steps below:a. Open CGAN_main.py with, e.g., PyCharmb. Change the source file in line 172 mat = mat73.loadmat().c. Change the source file of saved weights in line 166 solver.load_weights(r"/ckpt/solver_199.ckpt")d. Change your sampling condition in line 187, condition = get_data_space(batch)Note that the loss functions and MSE will be saved using lines 248-251.The weighs of the generator and the discriminator will be saved using lines 254 and 255. After training, the saved weights can be used for the inverse design.Some generated geometries will be saved using line 237, which can be read using 3d_plot_voxel.npy.3. ConclusionsWe give a tutorial of how to use 3D-CGAN for the inverse design of 3D geometries. After suitable training, the CGAN can take target properties as inputs and output corresponding 3D geometries. Also, you can find the 2D case in our previous tutorial and paper.[1,4] If you meet any problems or have any comments, just feel free to contact us via emails. Hope you can enjoy our codes.Reference1. Zheng X, Watanabe I. Tutorial for conditional generative adversarial network. https://doi.org/10.48505/nims.3869.2. Zheng X, Chen TT, Jiang X, Naito M, Watanabe I. Deep-learning-based inverse design of three-dimensional architected cellular materials with the target porosity and stiffness using voxelized Voronoi lattices. Science and Technology of Advanced Materials. 2023 Dec 31;24(1):2157682.3. Zheng X, Zhang X, Chen TT, Watanabe I. Deep Learning in Mechanical Metamaterials: From Prediction and Generation to Inverse Design. Advanced Materials. 2023 Jun 18:2302530.4. Zheng X, Chen TT, Guo X, Samitsu S, Watanabe I. Controllable inverse design of auxetic metamaterials using deep learning. Materials & Design. 2021 Dec 1;211:110178.image1.pngimage2.pngregression_with_ckpt.pyimport osimport numpy as npos.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'import tensorflow as tffrom tensorflow import kerasfrom tensorflow.keras import datasets, layers, optimizers, Sequential, metricsimport glob, datetimeimport matplotlib.pyplot as pltfrom scipy.io import loadmatimport mat73physical_devices = tf.config.list_physical_devices('GPU')tf.config.experimental.set_memory_growth(physical_devices[0], True)current_time = datetime.datetime.now().strftime("%Y%m%d-%H%M%S")log_dir = 'logs/' + current_timesummary_writer = tf.summary.create_file_writer(log_dir)def loaddata(filepath):    geo=mat73.loadmat(filepath+"/geo.mat")    geo=geo["geo"]    rho = loadmat(filepath+"/rho.mat")    rho = rho["rho"]    young = loadmat(filepath+"/Young.mat")    young = young["Young"]/150    diffusion = loadmat(filepath+"/diffusion_average.mat")    diffusion = diffusion["diffusion_average"]    # Random shuffling images and labels    np.random.seed(9487)    indice = np.array(range(len(geo)))    np.random.shuffle(indice)    geo = geo[indice]    geo = tf.cast(tf.expand_dims(geo, -1), tf.float64)    labels = np.concatenate((rho, young, diffusion), -1)    labels = tf.cast(labels[indice], tf.float64)    return geo, labelsdef new_tri_plot():    # lims_young = [0, 12]    # lims_poi = [-0.5, 0.5]    fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=[12, 4], constrained_layout=True)    # ax1.plot(lims_young, lims_young, 'gray')    # ax1.set_xlim(lims_young)    # ax1.set_ylim(lims_young)    ax1.set_aspect(1)    ax1.set_xlabel("True values")    ax1.set_ylabel("Predictions")    ax1.set_title("rho")    # ax2.plot(lims_poi, lims_poi, 'gray')    # ax2.set_xlim(lims_poi)    # ax2.set_ylim(lims_poi)    ax2.set_aspect(1)    ax2.set_xlabel("True values")    ax2.set_ylabel("Predictions")    ax2.set_title("Young's modulus")    ax3.set_aspect(1)    ax3.set_xlabel("True values")    ax3.set_ylabel("Predictions")    ax3.set_title("diffusion")    return ax1, ax2, ax3def tri_plot(y, pred, ax1, ax2, ax3):    x1 = y[:, 0]    x2 = y[:, 1]    x3 = y[:, 2]    y1 = pred[:, 0]    y2 = pred[:, 1]    y3 = pred[:, 2]    ax1.scatter(x1, y1, c='C0', alpha=0.2)    ax2.scatter(x2, y2, c='C1', alpha=0.2)    ax3.scatter(x3, y3, c='C2', alpha=0.2)def pad_dim(x, n=1):    x = tf.concat((x[:, :, :, -n:, :], x, x[:, :, :, :n, :]), axis=-2)    x = tf.concat((x[:, :, -n:, :, :], x, x[:, :, :n, :, :]), axis=-3)    x = tf.concat((x[:, -n:, :, :, :], x, x[:, :n, :, :, :]), axis=-4)    return xclass Solver(keras.Model):    def __init__(self):        super(Solver, self).__init__()        # unit 1, [64,64,64,1] => [32,32,32,16]        self.conv1a = layers.Conv3D(16, kernel_size=3, padding='valid', activation=tf.nn.relu)        self.conv1b = layers.Conv3D(16, kernel_size=3, padding='valid', activation=tf.nn.relu)        self.max1 = layers.MaxPool3D(pool_size=2, strides=2, padding='valid')        # unit 2, => [16,16,16,32]        self.conv2a = layers.Conv3D(32, kernel_size=3, padding='valid', activation=tf.nn.relu)        self.conv2b = layers.Conv3D(32, kernel_size=3, padding='valid', activation=tf.nn.relu)        self.max2 = layers.MaxPool3D(pool_size=2, strides=2, padding='valid')        # unit 3, => [8,8,8,64]        self.conv3a = layers.Conv3D(64, kernel_size=3, padding='valid', activation=tf.nn.relu)        self.conv3b = layers.Conv3D(64, kernel_size=3, padding='valid', activation=tf.nn.relu)        self.max3 = layers.MaxPool3D(pool_size=2, strides=2, padding='valid')        # unit 4, => [4,4,4,128]        self.conv4a = layers.Conv3D(128, kernel_size=3, padding='valid', activation=tf.nn.relu)        self.conv4b = layers.Conv3D(128, kernel_size=3, padding='valid', activation=tf.nn.relu)        self.max4 = layers.MaxPool3D(pool_size=2, strides=2, padding='valid')        # unit 5, => [2,2,2,256]        self.conv5a = layers.Conv3D(256, kernel_size=3, padding='valid', activation=tf.nn.relu)        self.conv5b = layers.Conv3D(256, kernel_size=3, padding='valid', activation=tf.nn.relu)        self.max5 = layers.MaxPool3D(pool_size=2, strides=2, padding='valid')        # unit 6, => [1,1,1,512]        self.conv6a = layers.Conv3D(512, kernel_size=3, padding='valid', activation=tf.nn.relu)        self.conv6b = layers.Conv3D(512, kernel_size=3, padding='valid', activation=tf.nn.relu)        self.max6 = layers.MaxPool3D(pool_size=2, strides=2, padding='valid')        # unit 7, => [2]        self.fc1 = layers.Dense(256, activation=tf.nn.relu)        self.fc2 = layers.Dense(128, activation=tf.nn.relu)        self.fc3 = layers.Dense(3, activation=None)    def call(self, x):        # inputs_noise: (b, 64), inputs_condition: (b, 3)        x = pad_dim(x)        x = self.conv1a(x)        x = pad_dim(x)        x = self.conv1b(x)        x = self.max1(x)        x = pad_dim(x)        x = self.conv2a(x)        x = pad_dim(x)        x = self.conv2b(x)        x = self.max2(x)        x = pad_dim(x)        x = self.conv3a(x)        x = pad_dim(x)        x = self.conv3b(x)        x = self.max3(x)        x = pad_dim(x)        x = self.conv4a(x)        x = pad_dim(x)        x = self.conv4b(x)        x = self.max4(x)        x = pad_dim(x)        x = self.conv5a(x)        x = pad_dim(x)        x = self.conv5b(x)        x = self.max5(x)        x = pad_dim(x)        x = self.conv6a(x)        x = pad_dim(x)        x = self.conv6b(x)        x = self.max6(x)        x = tf.keras.layers.Flatten()(x)        x = self.fc1(x)        x = self.fc2(x)        x = self.fc3(x)        return xdef preprocess(x, y):    x = tf.cast(x, dtype=tf.float32)    # x = tf.expand_dims(x, -1)    y = tf.cast(y, dtype=tf.float32)    return x, ydef main():    # tf.random.set_seed(2345)    filepath = r'/home/xiaoyang/PycharmProjects/pythonProject30_3d_foam/cgan20220421/dataset'    dataset_matrixes, dataset_labels = loaddata(filepath)    data_length = dataset_labels.shape[0]    x_test = dataset_matrixes[int(0.8 * len(dataset_matrixes)):-1]    y_test = dataset_labels[int(0.8 * len(dataset_labels)):-1]    x = dataset_matrixes[0:int(0.8 * len(dataset_matrixes))]    y = dataset_labels[0:int(0.8 * len(dataset_labels))]    train_db = tf.data.Dataset.from_tensor_slices((x, y))    train_db = train_db.shuffle(10000).map(preprocess).batch(32)    test_db = tf.data.Dataset.from_tensor_slices((x_test, y_test))    test_db = test_db.map(preprocess).batch(32)    solver = Solver()    solver.build(input_shape=[None, 64, 64, 64, 1])    optimizer = optimizers.Adam(learning_rate=1e-4, beta_1=0.5)  ## lr = 2e-4, beta_1 = 0.9    # plot test    ax1, ax2, ax3 = new_tri_plot()    total_sum = 0    total_error = 0    for x, y in test_db:        pred = solver(x)        loss_test = tf.losses.mean_squared_error(y, pred)        loss_test = tf.reduce_mean(loss_test)        total_sum += 1        total_error += loss_test        tri_plot(y, pred, ax1, ax2, ax3)    mse = total_error / total_sum    print(0, "mse: ", mse.numpy())    plt.savefig('results/test/%d_test.png' % 0)    plt.close()    with summary_writer.as_default():        tf.summary.scalar('loss:', float(mse.numpy()), step=0)        tf.summary.scalar('mse', float(mse.numpy()), step=0)    for epoch in range(200):        # plot train        ax1, ax2, ax3 = new_tri_plot()        epoch += 1        for step, (x, y) in enumerate(train_db):            with tf.GradientTape() as tape:                logits = solver(x)                loss = tf.losses.mean_squared_error(y, logits)                loss = tf.reduce_mean(loss)            grads = tape.gradient(loss, solver.trainable_variables)            optimizer.apply_gradients(zip(grads, solver.trainable_variables))            if step < 300:                tri_plot(y, logits, ax1, ax2, ax3)        plt.savefig('results/train/%d_train.png' % epoch)        plt.close()        print(epoch, "loss: ", loss.numpy())        # plot test        ax1, ax2, ax3 = new_tri_plot()        total_sum = 0        total_error = 0        for x, y in test_db:            pred = solver(x)            loss_test = tf.losses.mean_squared_error(y, pred)            loss_test = tf.reduce_mean(loss_test)            total_sum += 1            total_error += loss_test            tri_plot(y, pred, ax1, ax2, ax3)        mse = total_error / total_sum        print(epoch, "mse: ", mse.numpy())        plt.savefig('results/test/%d_test.png' % epoch)        plt.close()        with summary_writer.as_default():            tf.summary.scalar('loss:', float(loss.numpy()), step=epoch)            tf.summary.scalar('mse', float(mse.numpy()), step=epoch)        solver.save_weights('ckpt/solver_%d.ckpt' % epoch)if __name__ == "__main__":    main()