Width diagrams now computes using ordered calculations, and parallelizes it. Also, small fixes for other scripts.

This commit is contained in:
Kenneth Jao 2021-09-27 04:59:45 -04:00
parent d31029d95d
commit edf6ad9659
6 changed files with 144 additions and 167 deletions

View File

@ -1,16 +1,22 @@
from pathlib import Path
import sys, numpy as np
import numpy as np, pickle, sys
from squish import Simulation
def main():
n = int(sys.argv[1])
all_widths = set(np.round(np.arange(3, 10.05, 0.05), 2))
for file in Path(f"simulations/Radial[T]T - N{n}R4.0").iterdir():
i = file.name.index("x")
all_widths.remove(float(file.name[i-4:i]))
for file in Path(f"squish_output/Radial[T]Search - N{n} - 500").iterdir():
sim, frames = Simulation.load(file / 'data.squish')
if sim.domain.n == n:
try:
all_widths.remove(next(frames)["domain"][1])
except StopIteration:
pass
remain_widths = sorted(list(all_widths))[::-1]
print(remain_widths)
print([int(round((10-w)/.05)) for w in remain_widths])
print("Remaining:", remain_widths)
if __name__ == "__main__":

View File

@ -4,6 +4,7 @@ import argparse, pickle, numpy as np, os
from pathlib import Path
import matplotlib.pyplot as plt
from squish import Simulation
def main():
parser = argparse.ArgumentParser("Graphs convergence graphs for a collection of simulations.")
@ -16,12 +17,14 @@ def main():
data = {}
for file in Path(args.sims_path).iterdir():
with open(file, "rb") as f:
all_info, _ = pickle.load(f)
sim = Simulation.load(file / "data.squish")
sim_info = next(sim)
step = sim_info["step_size"]
for frame in sim:
step = float(file.name[:-4].split("-")[1])
data[step] = {"times": [], "values": [], "diffs": []}
for i, frame_info in enumerate(all_info):
for i, frame_info in enumerate(sim):
data[step]["times"].append(step*i)
data[step]["values"].append(np.linalg.norm(frame_info["arr"]))
data[step]["diffs"].append(np.linalg.norm(all_info[-1]["arr"] - frame_info["arr"]))

View File

@ -1,101 +1,93 @@
from __future__ import annotations
from typing import List
import os, math, argparse, numpy as np, pickle
from typing import List, Tuple, Dict
import argparse, math, numpy as np, os
import matplotlib.pyplot as plt
import matplotlib.ticker as mtick
from multiprocessing import Pool, cpu_count
from pathlib import Path
from ..simulation import Diagram, Simulation
from .._squish import AreaEnergy, RadialALEnergy, RadialTEnergy
import squish.ordered as order
from squish import Simulation, DomainParams
from squish.common import OUTPUT_DIR
ENERGY_R_STR = {AreaEnergy: "Area", RadialALEnergy: "Radial[AL]", RadialTEnergy: "Radial[T]"}
ENERGY_I_STR = {AreaEnergy: "area", RadialALEnergy: "radial-al", RadialTEnergy: "radial-t"}
I_TO_R = {"area": "Area","radial-t": "Radial[AL]", "radial-t": "Radial[T]"}
def order_process(domain: DomainParams) -> Tuple[float, float, float]:
energies = []
configs = order.configurations(domain)
for config in configs:
energies.append(2*domain.w*domain.h + \
2*math.pi*domain.n*(domain.r**2 - 2*domain.r*order.avg_radius(domain, config)))
def get_torus_config_energies(n: int, widths: np.ndarray, h: float, r: float,
energy: str) -> Tuple:
sim_file = SIM_FOLDER / f"{I_TO_R[energy]} - TorusConfigEnergy - N{n}.data"
if sim_file.is_file():
with open(sim_file, "rb") as data:
return pickle.load(data)
torus_min_energies, torus_max_energies = np.empty(widths.shape), np.empty(widths.shape)
torus_min_configs, torus_max_configs = [None]*len(widths), [None]*len(widths)
for i, w in enumerate(widths):
sim = Simulation(n, w, h, r, energy)
configs = []
for j in range(2):
for c in range(1,n): # Ignore 0, tends to error.
config = (1,c) if j == 0 else (c,1)
sim.add_frame(torus=config)
configs.append(config)
# eigs = np.sort(np.linalg.eig(sim.frames[-1].hessian(10e-5))[0])
# if eigs[0] > 1e-4:
# del sim.frames[-1]
# del config[-1]
hashes = int(21*i/len(widths))
print(f'Generating at width {w:.02f}... ' + \
f'|{"#"*hashes}{" "*(20-hashes)}| {i+1}/{len(widths)}, ' + \
f'{c + (n-1)*j}/{2*(n-1)} completed.', flush=True, end='\r')
pair = list(zip(configs,[frame.energy for frame in sim.frames]))
torus_min_configs[i], torus_min_energies[i] = min(pair, key=lambda x: x[1])
torus_max_configs[i], torus_max_energies[i] = max(pair, key=lambda x: x[1])
print(flush=True)
out_tup = (torus_min_energies, torus_max_energies, torus_min_configs, torus_max_configs)
with open(sim_file, "wb") as output:
pickle.dump(out_tup, output)
return out_tup
return domain.w, min(energies), max(energies)
def get_equilibria_data(filepath: Path):
if filepath.is_file():
with open(filepath, "rb") as data:
return pickle.load(data)
def get_ordered_energies(orig_domain: DomainParams, widths: np.ndarray) -> Dict:
data = {}
domains = [DomainParams(orig_domain.n, w, orig_domain.h, orig_domain.r) for w in widths]
with Pool(cpu_count()) as pool:
mins, maxes = {}, {}
for i, res in enumerate(pool.imap_unordered(order_process, domains)):
mins[res[0]] = res[1]
maxes[res[0]] = res[2]
hashes = int(21*i/len(widths))
print(f'Generating at width {res[0]:.02f}... |{"#"*hashes}{" "*(20-hashes)}|' + \
f' {i+1}/{len(widths)} completed.', flush=True, end='\r')
print(flush=True)
data["min"] = list([x[1] for x in sorted(mins.items())])
data["max"] = list([x[1] for x in sorted(maxes.items())])
return data
def eq_file_process(file: Path) -> Tuple[float, List[float], List[float]]:
sim, frames = Simulation.load(file / 'data.squish')
alls = []
for frame_info in frames:
alls.append([
frame_info["energy"],
np.var(frame_info["stats"]["avg_radius"]) <= 1e-8,
np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6)
])
sim, frames = Simulation.load(file / 'data.squish')
sim.frames = list(frames)
counts = sim.get_distinct()
distincts = []
for j, frame_info in enumerate(sim.frames):
distincts.append([
frame_info["energy"],
np.var(frame_info["stats"]["avg_radius"]) <= 1e-8,
np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6),
counts[j]
])
return sim.domain.w, alls, distincts
def get_equilibria_data(filepath: Path) -> Tuple[Dict, numpy.ndarray, DomainParams]:
data = {"all": {}, "distinct": {}}
files = list(Path(filepath).iterdir())
for i, file in enumerate(files):
sim = Simulation.load(file)
data["all"][sim.w] = []
for frame in sim.frames:
data["all"][sim.w].append([
frame.energy,
np.var(frame.stats["avg_radius"]) <= 1e-8,
np.count_nonzero(frame.stats["site_edge_count"] != 6)
])
with Pool(cpu_count()) as pool:
for i, res in enumerate(pool.imap_unordered(eq_file_process, files)):
data["all"][res[0]] = res[1]
data["distinct"][res[0]] = res[2]
sim.get_distinct()
data["distinct"][sim.w] = []
for frame in sim.frames:
data["distinct"][sim.w].append([
frame.energy,
np.var(frame.stats["avg_radius"]) <= 1e-8,
np.count_nonzero(frame.stats["site_edge_count"] != 6)
])
hashes = int(21*i/len(files))
print(f'Loading simulations... |{"#"*hashes}{" "*(20-hashes)}|' + \
hashes = int(21*i/len(files))
print(f'Loading simulations... |{"#"*hashes}{" "*(20-hashes)}|' + \
f' {i+1}/{len(files)} simulations loaded.', flush=True, end='\r')
print(flush=True)
print(flush=True)
sim, frames = Simulation.load(files[0] / 'data.squish')
widths = np.asarray(sorted(data["all"]))
n, h, r, energy = sim.n, sim.h, sim.r, sim.energy
domain = DomainParams(sim.domain.n, widths[-1], sim.domain.h, sim.domain.r)
return data, widths, domain
sim_file = SIM_FOLDER / f"{ENERGY_R_STR[energy]} - EquilibriaData - N{n}.data"
out_tup = (widths, data, n, h, r, ENERGY_I_STR[energy])
with open(sim_file, "wb") as output:
pickle.dump(out_tup, output)
return out_tup
def axis_settings(ax, widths):
ax.invert_xaxis()
@ -107,7 +99,7 @@ def axis_settings(ax, widths):
def main():
# Loading arguments.
parser = argparse.ArgumentParser("Compiles the equilibriums for each width into a diagram.")
parser = argparse.ArgumentParser("Outputs width search data into diagrams")
parser.add_argument('sims_path', metavar='path/to/data',
help="folder that contains simulation files, or cached data file.")
parser.add_argument('-q', '--quiet', dest='quiet', action='store_true', default=False,
@ -115,13 +107,10 @@ def main():
args = parser.parse_args()
widths, data, n, h, r, energy = get_equilibria_data(Path(args.sims_path))
data, widths, domain = get_equilibria_data(Path(args.sims_path))
order_data = get_ordered_energies(domain, widths)
torus_min_energies, torus_max_energies, _, _ = get_torus_config_energies(
n, widths, h, r, energy
)
fig_folder = Path(f"figures/ShrinkEnergyComparison - N{n}")
fig_folder = OUTPUT_DIR / Path(f"ShrinkEnergyComparison - N{domain.n}")
fig_folder.mkdir(exist_ok=True)
# Torus minimum energies used as reference.
@ -136,7 +125,7 @@ def main():
ax.plot(widths, all_disorder_count)
axis_settings(ax, widths)
ax.yaxis.set_major_formatter(mtick.PercentFormatter())
ax.title.set_text('Probability of Disorder')
ax.title.set_text(f"Probability of Disorder - N{domain.n}")
ax.set_xlabel("Width")
ax.set_ylabel("Disordered Equilibria")
boa_y_min = round(min(all_disorder_count)/20)*20 - 5
@ -152,15 +141,22 @@ def main():
distinct_ordered.append(equal_shape.count(True))
distinct_unordered.append(equal_shape.count(False))
ax.plot(widths, distinct_unordered, label="Unordered Equilibria")
ax.plot(widths, distinct_ordered, label="Ordered Equilibria")
ax.legend()
ax2 = ax.twinx()
ax.plot(widths, distinct_unordered, label="Unordered Equilibria", color='C0')
ax2.plot(widths, distinct_ordered, label="Ordered Equilibria", color='C1')
axis_settings(ax, widths)
ax.title.set_text('Density of States')
ax.title.set_text(f"Density of States - N{domain.n}")
ax.set_xlabel("Width")
ax.set_ylabel("Number of States")
dos_y_max = 1.05*max(distinct_ordered + distinct_unordered)
ax.set_yticks(np.arange(0, dos_y_max, round(dos_y_max/200, 1)*10))
ax.set_ylabel("Number of States (Disordered)", color='C0')
ax2.set_ylabel("Number of States (Ordered)", color='C1')
dos_y_max_unorder = 1.05*max(distinct_unordered)
dos_y_max_order = 1.05*max(distinct_ordered)
ax.set_yticks(np.linspace(0, dos_y_max_unorder, 20).astype(int))
#ax.set_yticks(np.arange(0, dos_y_max_unorder, round(dos_y_max_unorder/200, 1)*10))
ax2.set_yticks(np.arange(0, dos_y_max_order))
fig.savefig(fig_folder / "Density Of States.png")
# Defect density diagram
@ -172,7 +168,7 @@ def main():
ax.plot(widths, defects)
axis_settings(ax, widths)
ax.title.set_text('Average Defects')
ax.title.set_text(f"Average Defects - N{domain.n}")
ax.set_xlabel("Width")
ax.set_ylabel("Defects")
ax.set_yticks(np.arange(0, 1+max(defects), 0.5))
@ -187,53 +183,37 @@ def main():
ordered_energies.append([c[0] for c in data["distinct"][width] if c[1]])
unordered_energies.append([c[0] for c in data["distinct"][width] if not c[1]])
for i in range(len(torus_min_energies)):
ordered_energies[i].append(torus_min_energies[i])
ordered_energies[i].append(torus_max_energies[i])
for i in range(len(order_data["min"])):
ordered_energies[i].append(order_data["min"][i])
ordered_energies[i].append(order_data["max"][i])
null_unorder = []
for i, energies in enumerate(unordered_energies):
if len(energies) == 0:
null_unorder.append(i)
energies.append(torus_min_energies[i])
energies.append(order_data["min"][i])
min_order = np.asarray([min(width) for width in ordered_energies])
max_order = np.asarray([max(width) for width in ordered_energies])
min_unorder = np.asarray([min(width) for width in unordered_energies])
max_unorder = np.asarray([max(width) for width in unordered_energies])
min_unorder_off = min_unorder - torus_min_energies
max_unorder_off = max_unorder - torus_min_energies
ax.plot(widths, min_order - torus_min_energies, color='C1')
#ax.plot(widths, max_order - torus_min_energies, color='C1', linestyle='dotted')
min_unorder_off = min_unorder - min_order
max_unorder_off = max_unorder - min_order
ax.plot(widths, min_order - min_order, color='C1')
#ax.plot(widths, max_order - offset, color='C1', linestyle='dotted')
ax.plot(widths, min_unorder_off, color='C0')
ax.plot(widths, max_unorder_off, color='C0', linestyle='dotted')
axis_settings(ax, widths)
for i in null_unorder:
ax.scatter(widths[i], min_unorder[i] - torus_min_energies[i],
ax.scatter(widths[i], 0,
marker='X', color="blue", s=50, zorder=4)
# ax.scatter(widths[i], max_unorder[i] - torus_min_energies[i],
# ax.scatter(widths[i], max_unorder[i] - offset[i],
# marker='X', edgecolors="blue", facecolors='none', s=100, zorder=4)
# for i, marker in enumerate(min_markers):
# if marker:
# ax.scatter(widths[i], min_energies[i]-torus_min_energies[i],
# marker='H', color="orange", s=20, zorder=4)
# else:
# ax.scatter(widths[i], min_energies[i]-torus_min_energies[i],
# marker='d', color="blue", s=20, zorder=4)
# for i, marker in enumerate(max_markers):
# if marker:
# ax.scatter(widths[i], max_energies[i]-torus_min_energies[i],
# marker='H', edgecolors="orange", s=20, facecolors='none', zorder=4)
# else:
# ax.scatter(widths[i], max_energies[i]-torus_min_energies[i],
# marker='d', edgecolors="blue", s=20, facecolors='none', zorder=4)
ax.title.set_text('Reduced Energy vs. Width')
ax.title.set_text(f"Reduced Energy vs. Width - N{domain.n}")
ax.set_xlabel("Width")
ax.set_ylabel("Reduced Energy")
bif_y_max = np.max(np.abs(np.concatenate((min_unorder_off, max_unorder_off))))
@ -245,8 +225,6 @@ def main():
if __name__ == '__main__':
os.environ["QT_LOGGING_RULES"] = "*=false"
SIM_FOLDER = Path(f"simulations/ShrinkEnergyComparison")
SIM_FOLDER.mkdir(exist_ok=True)
try:
main()
except KeyboardInterrupt:

View File

@ -3,7 +3,7 @@ from typing import List, Tuple, Union, Optional, Iterator, Generator
import pickle, numpy as np
from math import gcd
from pathlib import Path
from ._squish import AreaEnergy, RadialALEnergy, RadialTEnergy
from ._squish import VoronoiContainer, AreaEnergy, RadialALEnergy, RadialTEnergy
OUTPUT_DIR = Path("squish_output")
OUTPUT_DIR.mkdir(exist_ok=True)
@ -84,7 +84,7 @@ class Energy:
except KeyError:
raise ValueError(f"\'{mode}\' is not a valid energy!")
else:
if mode is not VoronoiContainer and issubclaass(mode, VoronoiContainer):
if mode is not VoronoiContainer and issubclass(mode, VoronoiContainer):
raise ValueError("Provided class is not a valid energy!")
self.mode = mode

View File

@ -36,8 +36,9 @@ def get_config_generators(domain: DomainParams, config: Config) -> Tuple[Config,
all_sites = np.concatenate((q1, q1-[w,0], q1-[w,h], q1-[0,h]))[2:]
# Checking 0 < ax + by < v*v to make the sites are within the region.
tol = 1e-3
vdot = np.matmul(all_sites, v)
in_box = all_sites[np.where((0 <= vdot) & (vdot <= v.dot(v)))[0]]
in_box = all_sites[np.where((-tol <= vdot) & (vdot <= (v.dot(v)+tol)))[0]]
in_box = np.expand_dims(in_box, 0).swapaxes(0,1) # Used for the next step, getting site*site
w = in_box[np.argmin(np.squeeze(np.matmul(in_box, in_box.transpose(0,2,1))))].flatten()

View File

@ -62,8 +62,14 @@ class Simulation:
distinct_avg_radii, distinct_count, new_frames = [], [], []
for frame in self.frames:
avg_radii = np.sort(frame.stats["avg_radius"])
try:
stats = frame.stats
except AttributeError:
stats = frame["stats"] # When we have a loaded simulations.
avg_radii = np.sort(stats["avg_radius"])
is_in = False
for i, dist_radii in enumerate(distinct_avg_radii):
if np.allclose(avg_radii, dist_radii, atol=1e-5):
@ -73,6 +79,7 @@ class Simulation:
if not is_in:
distinct_avg_radii.append(avg_radii)
distinct_count.append(1)
new_frames.append(frame)
self.frames = new_frames
@ -97,7 +104,7 @@ class Simulation:
f = self[index]
info = {
"arr": f.site_arr,
"domain": (f.n, f.h, f.w, f.r),
"domain": (f.n, f.w, f.h, f.r),
"energy": f.energy,
"stats": f.stats
}
@ -123,8 +130,13 @@ class Simulation:
def load(path: str) -> Tuple[Simulation, Generator]:
def frames() -> Dict:
with open(path, 'rb') as infile:
first = True
while True:
try:
if first:
pickle.load(infile)
first = False
continue
yield pickle.load(infile)
except EOFError:
break
@ -139,31 +151,6 @@ class Simulation:
return sim, frames()
@staticmethod
def load_old(filename: str) -> Simulation:
"""
Loads the points at every point into a file.
:param filename: [str] name of the file
"""
frames = []
with open(filename, 'rb') as data:
all_info, sim_class = pickle.load(data)
if type(sim_class) == str:
sim_class = {"flow": Flow, "search": Search, "shrink": Shrink}[sim_class]
if all_info[0]["params"][0] in [53, 59, 61, 83, 131]:
thres = 1e-5
else:
thres = 1e-6
sim = sim_class(DomainParams(*all_info[0]["params"]), Energy("radial-t"), 0.05, thres, True, 0.1, 500)
for frame_info in all_info:
frames.append(sim.energy.mode(*frame_info["params"], frame_info["arr"]))
#frames[-1].stats = frame_info["stats"]
sim.frames = frames
return sim
class Flow(Simulation):
"""Finds an equilibrium from initial sites.
@ -409,3 +396,5 @@ STR_TO_SIM = {
"search": Search,
"shrink": Shrink
}
simulation = Simulation