Commit 27dec2cf authored by Ding, Lijie's avatar Ding, Lijie
Browse files

add Rmu and rename

parent 55eca3c2

analyze/ML_analyze.py

+114 −103
Original line number Diff line number Diff line
@@ -7,22 +7,26 @@ from scipy.optimize import curve_fit 
import pickle





def get_feature_Iq2D_IqIq_af_data(folder, parameters, random=False):

def get_feature_Iq_gq_data(folder, parameters, random=False):



    all_N = []  # system size related

    all_R0 = []  # suspension size related

    all_n = []  # system size related

    all_sigma = []

    all_gxx, all_gxy, all_gyx, all_gyy = [], [], [], []  # affine transformation parameters

    all_sqrtD = []

    all_gxy = []

    # all_theta, all_Sx, all_phi = [], [], []  # affine transformation parameters

    all_Iq2D = []

    all_Iq_flatten = []

    all_Iq2D_af = []

    all_IqIq_af_flatten = []

    all_gq = []

    all_gq_flatten = []

    all_finfo = []

    qr = []

    qphi = []

    for i in range(len(parameters)):

        if random:

            N, run_num = parameters[i]

            finfo = f"N{N:.0f}_random_run{run_num:.0f}"

            run_num = parameters[i][0]

            finfo = f"random_run{run_num:.0f}"

        else:

            N, sigma, theta, Sx, phi = parameters[i]

            finfo = f"{N:.0f}_sigma{sigma:.1f}_theta{theta:.1f}_Sx{Sx:.1f}_phi{phi:.1f}"
@@ -33,46 +37,53 @@ def get_feature_Iq2D_IqIq_af_data(folder, parameters, random=False): 
            print(f"File not found: {filename}")

            continue

        data = np.genfromtxt(filename, delimiter=",", skip_header=1)

        bnum_phi = len(data[0]) - 11

        bnum_phi = len(data[0]) - 7

        bnum_r = int((len(data) - 2) / 3)



        N, sigma, gxx, gxy, gyx, gyy = data[0, 1], data[0, 6], data[0, 7], data[0, 8], data[0, 9], data[0, 10]

        qphi = data[1, 12:]

        qr = data[2 : 2 + bnum_r, 11]

        Iq2D = data[2 : 2 + bnum_r, 12:]

        Iq2D_af = data[2 + bnum_r : 2 + 2 * bnum_r, 12:]

        IqIq_af = data[2 + 2 * bnum_r : 2 + 3 * bnum_r, 12:]

        R0, n, sigma, sqrtD, gxy = data[0, 1], data[0, 2], data[0, 3], data[0, 4], data[0, 5]



        qphi = data[1, 7:]

        qr = data[2 : 2 + bnum_r, 6]

        Iq2D = data[2 : 2 + bnum_r, 7:]

        Iq2D_af = data[2 + bnum_r : 2 + 2 * bnum_r, 7:]

        IqIq_af = data[2 + 2 * bnum_r : 2 + 3 * bnum_r, 7:]



        all_Iq2D.append(Iq2D)

        all_Iq_flatten.append(Iq2D.flatten())

        all_Iq2D_af.append(Iq2D_af)

        all_IqIq_af_flatten.append(IqIq_af.flatten())

        gq = IqIq_af / (Iq2D * Iq2D)

        all_gq.append(gq)

        all_gq_flatten.append(gq.flatten())



        all_N.append(N)

        all_R0.append(R0)

        all_n.append(n)

        all_sigma.append(sigma)

        all_gxx.append(gxx)

        all_gxy.append(gxy)

        all_gyx.append(gyx)

        all_gyy.append(gyy)

        all_sqrtD.append(sqrtD)

        all_gxy.append(np.abs(gxy))



        all_finfo.append(finfo)



    all_feature = np.array([all_N, all_sigma, all_gxx, all_gxy, all_gyx, all_gyy]).T

    all_feature_name = ["N", "sigma", "gxx", "gxy", "gyx", "gyy"]

    all_feature_tex = [r"$N$", r"$\sigma$", r"$\gamma_{xx}$", r"$\gamma_{xy}$", r"$\gamma_{yx}$", r"$\gamma_{yy}$"]

    all_feature = np.array([all_R0, all_n, all_sigma, all_sqrtD, all_gxy]).T

    all_feature_name = ["R0", "n", "sigma", "sqrtD", "gxy"]

    all_feature_tex = [r"$R_0$", r"$n$", r"$\sigma$", r"$\sqrt{D}$", r"$\gamma_{xy}$"]



    qphi = np.array(qphi)

    qr = np.array(qr)

    all_Iq2D = np.array(all_Iq2D)

    all_Iq2D_af = np.array(all_Iq2D_af/(Iq2D*Iq2D))

    all_IqIq_af_flatten = np.array(all_IqIq_af_flatten)

    all_Iq_flatten = np.array(all_Iq_flatten)

    all_gq_flatten = np.array(all_gq_flatten)



    return all_feature, all_feature_name, all_feature_tex, all_Iq2D, all_Iq2D_af, all_IqIq_af_flatten, qr, qphi

    return all_feature, all_feature_name, all_feature_tex, all_Iq_flatten, all_gq_flatten, qr, qphi





def calc_svd(folder, parameters):

    all_feature, all_feature_name, all_feature_tex, all_Iq2D, all_Iq2D_af, all_IqIq_af_flatten, qr, qphi = get_feature_Iq2D_IqIq_af_data(folder, parameters, random=True)

    all_feature, all_feature_name, all_feature_tex, all_Iq_flatten, all_gq_flatten, qr, qphi = get_feature_Iq_gq_data(folder, parameters, random=True)



    F = [np.log(all_Iq_flatten), all_gq_flatten]

    Fname = ["logIq", "gq"]

    # for k in range(len(F)):

    for k in [1]:

        print("all_feature shape:", np.array(all_feature).shape)

    svd = np.linalg.svd(all_IqIq_af_flatten)

        svd = np.linalg.svd(F[k])

        print(svd.S)

        print("np.array(svd.U).shape", np.array(svd.U).shape)

        print("np.array(svd.S).shape", np.array(svd.S).shape)
@@ -89,30 +100,30 @@ def calc_svd(folder, parameters): 
        QPHI, QR = np.meshgrid(qphi, qr)

        ax01 = plt.subplot(2, 2, 2, projection="polar")

        print("np.minimum(svd.Vh[0]), np.maximum(svd.Vh[0])", svd.Vh[0].min(), svd.Vh[0].max())

    ax01.contourf(QPHI, QR, svd.Vh[0].reshape(np.shape(QPHI)), levels=np.linspace(-0.3, 0.2, 10), cmap="rainbow")

        ax01.contourf(QPHI, QR, svd.Vh[0].reshape(np.shape(QPHI)), cmap="rainbow")  # , levels=np.linspace(-0.3, 0.2, 10), cmap="rainbow")

        ax01.set_title("Left Singular Vectors (svd.Vh[0])")



        ax11 = plt.subplot(2, 2, 3, projection="polar")

        print("np.minimum(svd.Vh[1]), np.maximum(svd.Vh[1])", svd.Vh[1].min(), svd.Vh[1].max())

    ax11.contourf(QPHI, QR, svd.Vh[1].reshape(np.shape(QPHI)), levels=np.linspace(-0.3, 0.2, 10), cmap="rainbow")

        ax11.contourf(QPHI, QR, svd.Vh[1].reshape(np.shape(QPHI)),cmap="rainbow")  # , levels=np.linspace(-0.3, 0.2, 10), cmap="rainbow")

        ax11.set_title("Left Singular Vectors (svd.Vh[1])")



        ax12 = plt.subplot(2, 2, 4, projection="polar")

        print("np.minimum(svd.Vh[2]), np.maximum(svd.Vh[2])", svd.Vh[2].min(), svd.Vh[2].max())

    ax12.contourf(QPHI, QR, svd.Vh[2].reshape(np.shape(QPHI)), levels=np.linspace(-0.3, 0.2, 10), cmap="rainbow")

        ax12.contourf(QPHI, QR, svd.Vh[2].reshape(np.shape(QPHI)),cmap="rainbow" ) #, levels=np.linspace(-0.3, 0.2, 10), cmap="rainbow")

        ax12.set_title("Left Singular Vectors (svd.Vh[2])")



        plt.tight_layout()

    plt.savefig(f"{folder}/svd.png", dpi=300)

        plt.savefig(f"{folder}/{Fname[k]}_svd.png", dpi=300)

        plt.show()

        plt.close()



    SqV = np.inner(all_IqIq_af_flatten, np.transpose(svd.Vh))

        SqV = np.inner(F[k], np.transpose(svd.Vh))

        plt.figure()

        fig = plt.figure(figsize=(2 * len(all_feature_name), 8))

        axs = [fig.add_subplot(2, len(all_feature_name) // 2 + 1, i + 1, projection="3d") for i in range(len(all_feature_name))]

        for i in range(len(all_feature_name)):

        scatter = axs[i].scatter(SqV[:, 0], SqV[:, 1], SqV[:, 2], c=all_feature[:, i], cmap="jet_r", s=2)

            scatter = axs[i].scatter(SqV[:, 0], SqV[:, 1], SqV[:, 2], c=all_feature[:, i], cmap="rainbow", s=1)

            axs[i].set_xlabel("V[0]")

            axs[i].set_ylabel("V[1]")

            axs[i].set_zlabel("V[2]")
@@ -131,7 +142,7 @@ def calc_svd(folder, parameters): 
            axs[i].view_init(elev=10.0, azim=-30)



        plt.tight_layout()

    plt.savefig(f"{folder}/svd_projection_scatter_plot.png", dpi=300)

        plt.savefig(f"{folder}/{Fname[k]}_svd_projection_scatter_plot.png", dpi=300)

        plt.show()

        plt.close()
@@ -145,7 +156,7 @@ def calc_svd(folder, parameters): 
        # save svd projection data

        data = np.column_stack((all_feature, SqV[:, 0], SqV[:, 1], SqV[:, 2]))

        column_names = all_feature_name + ["sqv[0]", "sqv[1]", "sqv[2]"]

    np.savetxt(f"{folder}/data_svd_projection.txt", data, delimiter=",", header=",".join(column_names), comments="")

        np.savetxt(f"{folder}/data_{Fname[k]}_svd_projection.txt", data, delimiter=",", header=",".join(column_names), comments="")





def calc_Sq_pair_distance_distribution(all_Delta_Sq, max_z, bin_num):

analyze/main_ML_analyze.py

+5 −6
Original line number Diff line number Diff line
@@ -10,15 +10,14 @@ import time 
def main():



    print("analyzing data using ML model")

    folder = "../data/20241119_rand"

    rand_num = 1000

    rand_max = 1000

    folder = "../data/20241120_rand"

    rand_num = 4000

    rand_max = 4000

    parameters = []

    n = 100

    for i in range(rand_num):

        filename = f"{folder}/obs_N{n}_random_run{i}.csv"

        filename = f"{folder}/obs_random_run{i}.csv"

        if os.path.exists(filename):

            parameters.append([n, i])

            parameters.append([i])

        if len(parameters) >= rand_max:

            break

    print("parameters", parameters)

code/main.cpp

+11 −9
Original line number Diff line number Diff line
@@ -5,7 +5,7 @@ 
#include <string>

#include <chrono>

#include <filesystem>

#include "ideal_gas.h"

#include "suspension.h"



int main(int argc, char const *argv[])

{
@@ -27,24 +27,26 @@ int main(int argc, char const *argv[]) 
    }



    double R0 = 0.00125;

    double Rmu = 1*R0;

    double sqrtD0 = 1*R0; // for simplicity: dx, dy = sqrt(2D0)*rand_norm(gen)

    // size dependence, Einstein: D ~ 1/a, a is particle diameter



    // precision run with specified parameters

    if (argc == 6)

    if (argc == 7)

    {

        int n = std::atoi(argv[1]); // number density

        double sigma = std::atof(argv[2]);

        double sqrtD = std::atof(argv[3])*R0;

        double Rmu = std::atof(argv[2])*R0;

        double sigma = std::atof(argv[3]);

        double sqrtD = std::atof(argv[4])*R0;

        // double theta = std::atof(argv[3]);

        // double Sx = std::atof(argv[4]);

        // double phi = std::atof(argv[5]);

        double gxy = std::atof(argv[4])*R0;

        std::string folder = std::string(argv[5]);

        double gxy = std::atof(argv[5])*R0;

        std::string folder = std::string(argv[6]);



        ideal_gas gas_2d(R0, n, sigma, sqrtD, gxy, false);

        suspension gas_2d(R0, Rmu, n, sigma, sqrtD, gxy, false);



        std::string finfo = "n" + std::string(argv[1]) + "_sigma" + std::string(argv[2]) + "_sqrtD" + std::string(argv[3]) + "_gxy" + std::string(argv[4]);

        std::string finfo = "n" + std::string(argv[1]) + "_Rmu" + std::string(argv[2]) + "_sigma" + std::string(argv[3]) + "_sqrtD" + std::string(argv[4]) + "_gxy" + std::string(argv[5]);



        int number_of_config = 2000;

        int bnum_r = 100;
@@ -61,7 +63,7 @@ int main(int argc, char const *argv[]) 


        // number_of_config 2000, bin_num 101, n 200 takes 468s to run



        ideal_gas gas_2d(R0, 100, 0.0, 1.0, 0.0, true);

        suspension gas_2d(R0, 1.0*R0, 100, 0.0, 1.0, 0.0, true);

        std::string finfo = "random_run" + std::to_string(run_num);



        int number_of_config = 4000;

code/makefile

+6 −6
Original line number Diff line number Diff line 
Compile=g++ -g -std=c++20 -O3 # improve run time e.g. 10s to 4s

ideal_gas: main.o ideal_gas.o

	$(Compile) -o ideal_gas main.o ideal_gas.o

main.o: main.cpp ideal_gas.h

suspension: main.o suspension.o

	$(Compile) -o suspension main.o suspension.o

main.o: main.cpp suspension.h

	$(Compile) -o main.o -c main.cpp

ideal_gas.o: ideal_gas.cpp ideal_gas.h

	$(Compile) -o ideal_gas.o -c ideal_gas.cpp

suspension.o: suspension.cpp suspension.h

	$(Compile) -o suspension.o -c suspension.cpp

clean:

	rm -f main.o ideal_gas.o
 
 No newline at end of file
 
	rm -f main.o suspension.o
 
 No newline at end of file

code/ideal_gas→code/suspension

+32 B (94 KiB)

File changed and moved.

No diff preview for this file type.

View original file
View changed file