scaling_plotter.py

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm 
from matplotlib.widgets import Slider

# these need to first be defined outide of any function since we are treating them as global later 
nu = 1.; beta = 1/8.; gamma = 1.75; z = 1. # ising values 
# nu = 1.; beta = .2; gamma = 1.6 # s3 values 

def scaling_plotter(data,plot,xc,nu0=1.,gamma0=1.75,beta0=.125,z0=2.1667,raw=False,d=2,title="",cmap=cm.coolwarm): 
    """
    makes interactive scaling plots for various critical exponents. 
    data: a dictionary with the following entries 
        * "Ls": system sizes 
        * "xs": raw error rates (not reduced to pc) for each system size (#Ls x #(x data points) matrix)
        * "chis": susceptibilities 
        * "mags": magnetizations 
        * "binds": binder cumulants 
        * "trels": memory lifetimes 
    plot: quantitiy to plot, one of ["binds" "mags" "chis" "trels"]
    xc: value of critical point 
    nu,gamma,beta: estimates for critical exponents 
    raw: if true, just plots the raw data without scaling
    d: dimension; used only for checking hyperscaling relations 
    title: plot title 
    critical exponents tuned using key presses. left/right tunes nu, right/left tunes gamma, and ,/. tunes beta 
    """
    # cmap = cm.coolwarm 
    ### asthetics for plots ## 
    plt.rcParams['font.family'] = 'serif'
    plt.rcParams['font.serif'] = ['Computer Modern Roman'] + plt.rcParams['font.serif']
    mew = 1.
    ms = 6.5 / 1.
    lw = 2.2 #* 0
    a = .8

    fig,ax = plt.subplots(figsize=(2.5,2.5))
    ax.minorticks_off()

    Ls = data["Ls"]; chis = data["chis"]; mags = data["mags"]; binds = data["binds"]; trels = data["trels"]

    ### plot raw unscaled data ### 
    if raw: 
        for l in range(len(Ls)): 
            col = cmap((l+1) / len(Ls))
            L = Ls[l]

            xs = data["xs"][l]; ts = (xs - xc)/xc 
            ys = data[plot][l] # stuff to plot at this system size 
            if plot == "trels":
                ys = L**0/ys # want to plot the inverse memory lifetime 
                ax.set_yscale('log'); ax.set_xscale('log')

            ax.plot(xs,ys,c=col,label=r'$%d$'%Ls[l],marker='o',mec='k',mew=mew,ms=ms,lw=lw)
        ax.set_xlabel(r'$p$'); ax.set_ylabel(r'$1/t_{\sf mem}$')
        if plot == "mags": 
            ax.set_ylabel(r'$M$')
        elif plot == "chis": 
            ax.set_ylabel(r'$\chi$')
        elif plot == "binds": 
            ax.set_ylabel(r'$B$')

        ax.legend(title=r'$L$',title_fontsize=12)
        ax.set_title(title)

    ### plot scaled data (in an interactive way) ### 
    else:  
        lines0 = [] # will hold the plot lines; dynamically updated 

        # draw lines with guessed values of critical exponents 
        for l in range(len(Ls)): 
            col = cmap((l+1) / len(Ls))
            L = Ls[l]

            xs = data["xs"][l]
            ts = (xs-xc)/xc # reduced temperatures 
            ys = data[plot][l] # stuff to plot at this system size 

            if plot == "binds": 
                line0, = ax.plot(ts * Ls[l]**(1/nu0),ys,c=col,label=r'$%d$'%Ls[l],marker='o',mfc=col,mew=mew,ms=ms,lw=lw*0,alpha=a)
            elif plot == "chis": 
                line0, = ax.plot(ts * Ls[l]**(1/nu0),ys * L**(-gamma0 / nu0),c=col,label=r'$%d$'%Ls[l],marker='o',mfc=col,mew=mew,ms=ms,lw=lw*0,alpha=a)
            elif plot == "mags": 
                line0, = ax.plot(ts * Ls[l]**(1/nu0),ys * L**(beta0 / nu0),c=col,label=r'$%d$'%Ls[l],marker='o',mfc=col,mew=mew,ms=ms,lw=lw*0,alpha=a)
            elif plot == "trels": 
                ys = 1/ys # want to plot the inverse memory lifetime 
                line0, = ax.plot(ts * Ls[l]**(1/nu0),ys * L**(z0),c=col,label=r'$%d$'%Ls[l],marker='o',mfc=col,mew=mew,ms=ms,lw=lw*0,alpha=a)
            lines0.append(line0)

        def update(nu,beta,gamma,z):
            print("nu       = ",nu)
            print("gamma    = ",gamma)
            print("beta     = ",beta)
            print("z        = ",z)
            print("beta(hs) = ",(d*nu-gamma)/2)
            print("eta (hs) = ",2-gamma/nu)
            print("Delta_ep = ",d-1/nu) # dimension of lightest Z2 neutral field
            print(" ")

            # will need to dynamically update the plot ranges 
            xmin = 1e10 
            xmax = -1e10
            ymin = 1e10 
            ymax = -1e10

            for l in range(len(Ls)): 
                L = Ls[l]
                xdat = (data["xs"][l]-xc)/xc * L**(1/nu)

                if plot == "chis":
                    ydat = data["chis"][l] * L**(-gamma/nu)
                elif plot == "mags": 
                    ydat = data["mags"][l] * L**(beta/nu)
                elif plot == "binds": 
                    ydat = data["binds"][l] 
                elif plot == "trels": 
                    ydat = (1/data["trels"][l]) * L**(z)
                xmin = min(xmin,np.min(xdat))
                xmax = max(xmax,np.max(xdat))
                ymin = min(ymin,np.min(ydat))
                ymax = max(ymax,np.max(ydat))

                # update the line data 
                lines0[l].set_xdata(xdat) 
                lines0[l].set_ydata(ydat)
            
            ax.set_xlim(xmin*1.05,xmax*1.05)
            ax.set_ylim(ymin*.95,ymax*1.05)
            fig.canvas.draw_idle()
            
            # Create dynamic title with current exponent values
            if plot == "binds": 
                dynamic_title = title + r"$\, \, | \, \, \nu = %.3f$"%nu
            elif plot == "chis": 
                dynamic_title = title + r"$\, \, | \, \, \nu = %.3f$"%nu
            elif plot == "mags": 
                dynamic_title = title + r"$\, \, | \, \, \beta = %.3f$"%beta
            elif plot == "trels": 
                dynamic_title = title + r"$\, \, | \, \, z = %.3f$"%z
            else:
                dynamic_title = title  # fallback for any other plot types
            ax.set_title(dynamic_title)


        # define a function that updates values of critical exponents based on arrow key presses 
        def on_key(event):
            global nu,beta,gamma,z
            d = .005 # incriment by which exponents are tuned 
            numin = 0; numax = 4
            gammamin = 0; gammamax = 4
            betamin = 0; betamax = 4  
            zmin = .75; zmax = 4
            if event.key == 'l':
                z = max(zmin,z-d)
            elif event.key == ':':
                z = min(zmax,z+d)
            if event.key == 'left':
                nu = max(numin,nu - d) 
            elif event.key == 'right':
                nu = min(numax,nu + d) 
            elif event.key == 'up':
                gamma = min(gammamax,gamma+d)
            elif event.key == 'down':
                gamma = max(gammamin,gamma-d)
            elif event.key == '.': 
                beta = max(betamin,beta-d)
            elif event.key == ',': 
                beta = min(betamax,beta+d)

            update(nu,beta,gamma,z)
            
            
        fig.canvas.mpl_connect('key_press_event', lambda event : on_key(event))

        ax.legend(title=r'$L$',title_fontsize=14)
        if plot == "binds": 
            ylab = r'$B$'
        elif plot == "chis": 
            ylab = r'$L^{-\gamma/\nu}\chi$'
        elif plot == "mags": 
            ylab = r'$L^{\beta/\nu} m$'
        elif plot == "trels": 
            ylab = r'$L^z / t_{\sf mem}$'
        ax.set(xlabel=r'$\overline{p}L^{1/\nu}$',ylabel = ylab)
        if plot == "trels": 
            ax.set_yscale('log')
    plt.show()