Commit 55eca3c2 authored by Ding, Lijie's avatar Ding, Lijie
Browse files

only have shear along y

parent 20854332

.gitignore

+2 −0
Original line number Diff line number Diff line 
/*

!code

code/*.o

code/local_run.sh

code/nohup.out

!analyze/

analyze/__pycache__

analyze/.DS_Store

analyze/ML_analyze.py

+52 −60
Original line number Diff line number Diff line
@@ -7,20 +7,22 @@ from scipy.optimize import curve_fit 
import pickle





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

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



    all_N = []  # system size related

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

    all_gamma_xx, all_gamma_xy, all_gamma_yx, all_gamma_yy = [], [], [], []  # affine transformation parameters

    #all_Sq2D = []

    #all_Sq2D_af = []

    all_SqSq2D_flatten = []

    all_sigma = []

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

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

    all_Iq2D = []

    all_Iq2D_af = []

    all_IqIq_af_flatten = []

    all_finfo = []

    qD = []

    qr = []

    qphi = []

    for i in range(len(parameters)):

        if random:

            n, run_num = parameters[i]

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

            N, run_num = parameters[i]

            finfo = f"N{N:.0f}_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}"
@@ -31,56 +33,46 @@ def get_feature_SqSq2D_data(folder, parameters, random=False): 
            print(f"File not found: {filename}")

            continue

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

        bin_num = len(data[0]) - 10

        print(f"bin_num: {bin_num}")



        qD = data[1, 10:]

        Sq2D = data[2 : bin_num + 2, 10:]

        Sq2D_af = data[bin_num + 2 : 2 * bin_num + 2, 10:]

        SqSq2D = data[2 * bin_num + 2 : 3 * bin_num + 2, 10:]



        center = bin_num // 2

        mask_range = 5

        Sq2D[center - mask_range : center + mask_range + 1, center - mask_range : center + mask_range + 1] = 0

        Sq2D_af[center - mask_range : center + mask_range + 1, center - mask_range : center + mask_range + 1] = 0

        SqSq2D[center - mask_range : center + mask_range + 1, center - mask_range : center + mask_range +1] = 0



        all_SqSq2D_flatten.append(SqSq2D.flatten())



        n = data[0, 1]

        theta = data[0, 2]

        Sx = data[0, 3]

        phi = data[0, 5]

        gamma_xx = data[0, 6]

        gamma_xy = data[0, 7]

        gamma_yx = data[0, 8]

        gamma_yy = data[0, 9]



        all_N.append(n)

        all_theta.append(theta)

        all_Sx.append(Sx)

        all_phi.append(phi)



        all_gamma_xx.append(gamma_xx)

        all_gamma_xy.append(gamma_xy)

        all_gamma_yx.append(gamma_yx)

        all_gamma_yy.append(gamma_yy)

        bnum_phi = len(data[0]) - 11

        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:]



        all_Iq2D.append(Iq2D)

        all_Iq2D_af.append(Iq2D_af)

        all_IqIq_af_flatten.append(IqIq_af.flatten())



        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_finfo.append(finfo)



    all_feature = np.array([all_theta, all_Sx, all_phi, all_gamma_xx, all_gamma_xy, all_gamma_yx, all_gamma_yy]).T

    all_feature_name = ["theta", "Sx", "phi", "gamma_xx", "gamma_xy", "gamma_yx", "gamma_yy"]

    all_feature_tex = [r"$\theta$", r"$S_x$", r"$\phi$", r"$\gamma_{xx}$", r"$\gamma_{xy}$", r"$\gamma_{yx}$", r"$\gamma_{yy}$"]

    qD = np.array(qD)

    all_SqSq2D_flatten = np.array(all_SqSq2D_flatten)

    return all_feature, all_feature_name, all_feature_tex, all_SqSq2D_flatten, qD

    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}$"]

    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)



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





def calc_svd(folder, parameters):

    all_feature, all_feature_name, all_feature_tex, all_SqSq2D_flatten, qD = get_feature_SqSq2D_data(folder, parameters, random=True)

    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)



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

    svd = np.linalg.svd(all_SqSq2D_flatten)

    svd = np.linalg.svd(all_IqIq_af_flatten)

    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)
@@ -94,20 +86,20 @@ def calc_svd(folder, parameters): 
    ax00.set_title("Singular Values (svd.S)")



    # Subplot for svd.U

    qDx, qDy = np.meshgrid(qD, qD)

    ax01 = plt.subplot(2, 2, 2)

    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(qDx, qDy, svd.Vh[0].reshape(np.shape(qDx)), levels=np.linspace(-0.3, 0.2, 10), cmap="rainbow")

    ax01.contourf(QPHI, QR, svd.Vh[0].reshape(np.shape(QPHI)), 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)

    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(qDx, qDy, svd.Vh[1].reshape(np.shape(qDx)), levels=np.linspace(-0.3, 0.2, 10), cmap="rainbow")

    ax11.contourf(QPHI, QR, svd.Vh[1].reshape(np.shape(QPHI)), 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)

    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(qDx, qDy, svd.Vh[2].reshape(np.shape(qDx)), levels=np.linspace(-0.3, 0.2, 10), cmap="rainbow")

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

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



    plt.tight_layout()
@@ -115,7 +107,7 @@ def calc_svd(folder, parameters): 
    plt.show()

    plt.close()



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

    SqV = np.inner(all_IqIq_af_flatten, 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))]
@@ -124,7 +116,7 @@ def calc_svd(folder, parameters): 
        axs[i].set_xlabel("V[0]")

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

        axs[i].set_zlabel("V[2]")

        axs[i].set_title(all_feature_name[i])

        axs[i].set_title(all_feature_tex[i])

        axs[i].set_box_aspect([1, 1, 1])  # Set the aspect ratio of the plot

        # Set the same range for each axis

        max_range = np.array([SqV[:, 0].max() - SqV[:, 0].min(), SqV[:, 1].max() - SqV[:, 1].min(), SqV[:, 2].max() - SqV[:, 2].min()]).max() / 2.0

analyze/main_ML_analyze.py

+3 −3
Original line number Diff line number Diff line
@@ -10,13 +10,13 @@ import time 
def main():



    print("analyzing data using ML model")

    folder = "../data/20241101"

    folder = "../data/20241119_rand"

    rand_num = 1000

    rand_max = 1000

    parameters = []

    n = 200

    n = 100

    for i in range(rand_num):

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

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

        if os.path.exists(filename):

            parameters.append([n, i])

        if len(parameters) >= rand_max:

analyze/main_analyze.py

+7 −14
Original line number Diff line number Diff line
@@ -9,30 +9,23 @@ def main(): 


    if 1:

        folder = "../data/data_local/data_pool"

        N = 20

        n, sigma, sqrtD, gxy = 100, 0.0, 1.0, 5.0

        parameters = [

            [N, 0.0, 1.00, 0.00, 0.00, 1.00],

            [N, 0.0, 1.00, 0.10, 0.00, 1.00],

            [N, 0.0, 1.00, 0.20, 0.00, 1.00],

            #[N, 0.0, 1.00, 0.20, 0.00, 1.00],

            #[N, 0.0, 1.00, 0.30, 0.00, 1.00],

            #[N, 0.0, 1.00, 1.00, 0.00, 1.00],

            #[N, 0.1, 1.00, 0.10, 0.00, 1.00],

            #[N, 1.0, 1.00, 0.10, 0.00, 1.00],

            [100, 0.0, 1.0, 10.0],

            [150, 0.0, 1.0, 10.0],

        ]

        for parameter in parameters:

            N, sigma, gxx, gxy, gyx, gyy = parameter

            finfo = f"N{N:.0f}_sigma{sigma:.1f}_gxx{gxx:.2f}_gxy{gxy:.2f}_gyx{gyx:.2f}_gyy{gyy:.2f}"

            n, sigma, sqrtD, gxy = parameter

            finfo = f"n{n:.0f}_sigma{sigma:.1f}_sqrtD{sqrtD:.1f}_gxy{gxy:.1f}"

            filename = folder + f"/config_{finfo}.csv"

            plot_gas_config(filename, finfo, show=True)

            # plot_gas_Sq_SqSq(folder, parameter, show=True)

            plot_gas_Iq_IqIq(folder, finfo, show=True)

    if 0:

        folder = "../data/data_local/data_pool"

        N = 200

        rand_num = 2

        for rnum in range(rand_num):

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

            finfo = f"random_run{rnum}"

            filename = folder + f"/config_{finfo}.csv"

            plot_gas_config(filename, finfo, show=True)

            # plot_gas_Sq_SqSq(folder, parameter, show=True)

analyze/plot_analyze.py

+37 −30
Original line number Diff line number Diff line
@@ -7,13 +7,16 @@ import os 
def plot_gas_config(filename, finfo, show=False):

    # Load the data

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

    R, x, y, x_af, y_af = data[:, 0], data[:, 1], data[:, 2], data[:, 3], data[:, 4]

    R0, R, x, y, x_af, y_af = data[:, 0], data[:, 1], data[:, 2], data[:, 3], data[:, 4], data[:, 5]

    # Plot the data

    fig = plt.figure()

    ax = fig.add_subplot(111)

    ax.plot([-0.5, 0.5, 0.5, -0.5, -0.5], [-0.5, -0.5, 0.5, 0.5, -0.5], color="black")  # frame

    ax.scatter(x, y, s=R * 3e3, color="royalblue", label="initial")

    ax.scatter(x_af, y_af, s=R * 3e3, color="tomato", label="affine-transformed")

    L = 1.0 / R0[0]

    ax.plot(np.array([-0.5, 0.5, 0.5, -0.5, -0.5]) * L, np.array([-0.5, -0.5, 0.5, 0.5, -0.5]) * L, color="black", lw=0.1)  # frame



    s = (1 * ax.get_window_extent().width / (1.1 * L + 1.0) * 72.0 / fig.dpi) ** 2

    ax.scatter(x / R0, y / R0, s=s * R / R0, facecolors="royalblue", edgecolors="none", label="initial")

    ax.scatter(x_af / R0, y_af / R0, s=s * R / R0, facecolors="tomato", edgecolors="none", label="affine-transformed")

    ax.legend()

    ax.set_xlabel("x")

    ax.set_ylabel("y")
@@ -23,7 +26,7 @@ def plot_gas_config(filename, finfo, show=False): 
    ax.set_axis_off()



    ax.set_title(finfo)

    plt.savefig(filename.replace(".csv", ".png"))

    plt.savefig(filename.replace(".csv", ".pdf"), format="pdf")

    if show:

        plt.show()

    plt.close()
@@ -148,10 +151,12 @@ def plot_gas_Sq_SqSq(folder, parameter, show=False): 


def get_feature_Iq2D_IqIq_af_data(folder, finfos):



    all_N = []  # system size related

    all_R0 = []  # system size related

    all_n = []  # system size related

    all_sigma = []

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

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

    all_sqrtD = []

    all_gxy = []



    all_Iq2D = []

    all_Iq2D_af = []

    all_IqIq_af = []
@@ -159,52 +164,52 @@ def get_feature_Iq2D_IqIq_af_data(folder, finfos): 
    qr = []

    qphi = []

    for i in range(len(finfos)):

        # N, sigma, gxx, gxy, gyx, gyy = parameters[i]

        # finfo = f"N{N:.0f}_sigma{sigma:.1f}_gxx{gxx:.2f}_gxy{gxy:.2f}_gyx{gyx:.2f}_gyy{gyy:.2f}"

        filename = f"{folder}/obs_{finfos[i]}.csv"



        if not os.path.exists(filename):

            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

        print(f"bnum_phi: {bnum_phi}")

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

        print(f"bnum_r: {bnum_r}")



        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_Iq2D_af.append(Iq2D_af)

        all_IqIq_af.append(IqIq_af)



        all_N.append(N)

        all_R0.append(R0)

        all_n.append(n)

        all_sigma.append(sigma)

        all_gxx.append(gxx)

        all_sqrtD.append(sqrtD)

        all_gxy.append(gxy)

        all_gyx.append(gyx)

        all_gyy.append(gyy)



    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 = 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)

    all_IqIq_af = np.array(all_IqIq_af)



    return all_feature, all_feature_name, all_Iq2D, all_Iq2D_af, all_IqIq_af, qr, qphi

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





def plot_gas_Iq_IqIq(folder, finfo, show=True):



    # get single data

    all_feature, all_feature_name, all_Iq2D, all_Iq2D_af, all_IqIq_af, qr, qphi = get_feature_Iq2D_IqIq_af_data(folder, [finfo])

    all_feature, all_feature_name, all_feature_tex, all_Iq2D, all_Iq2D_af, all_IqIq_af, qr, qphi = get_feature_Iq2D_IqIq_af_data(folder, [finfo])



    plt.figure(figsize=(9, 3))

    ax1 = plt.subplot(231, projection="polar")
@@ -274,14 +279,15 @@ def plot_gas_Iq_IqIq(folder, finfo, show=True): 
    ax5.set_xlabel(r"$q_\phi$")

    ax5.set_ylabel(r"$\sum_{r}g(q)$")





    for ax in [ax1, ax2, ax3]:

        ax.set_thetamin(0)

        ax.set_thetamax(180)

        ax.grid(False)

        ax.set_axis_off()



    plt.title(f"N={all_feature[0, 0]:.2f}, sigma={all_feature[0, 1]:.2f}, gxx={all_feature[0, 2]:.2f}, gxy={all_feature[0, 3]:.2f}, gyx={all_feature[0, 4]:.2f}, gyy={all_feature[0, 5]:.2f}")

    ax = plt.subplot(111)

    ax.set_axis_off()

    ax.text(0.75, 0.3, f"{finfo}", fontsize=12, ha="center", va="center", transform=ax.transAxes)

    plt.tight_layout()

    plt.savefig(f"{folder}/Iq_IqIq_{finfo}.png")

    if show:
@@ -292,6 +298,7 @@ def plot_gas_Iq_IqIq(folder, finfo, show=True): 
def calc_single_sphere_Iq(qr, R):

    return 3.0 / (qr * R) ** 3 * (np.sin(qr * R) - qr * R * np.cos(qr * R))





def test_plot():

    r = np.linspace(400, 4000, 40)

    phi = np.linspace(0, 2 * np.pi / 20 * 19, 20)