Skip to content

Commit

Permalink
[pre-commit.ci] auto fixes from pre-commit.com hooks
Browse files Browse the repository at this point in the history
for more information, see https://pre-commit.ci
  • Loading branch information
pre-commit-ci[bot] committed Jun 20, 2024
1 parent 4ce7ac6 commit 24fcd39
Show file tree
Hide file tree
Showing 3 changed files with 72 additions and 48 deletions.
5 changes: 1 addition & 4 deletions ipsuite/analysis/model/dynamics_checks.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,7 @@
import numpy as np
import zntrack
from ase.geometry import conditional_find_mic
from ase.neighborlist import build_neighbor_list
from ase.neighborlist import natural_cutoffs
from ase.neighborlist import build_neighbor_list, natural_cutoffs

from ipsuite import base
from ipsuite.utils.ase_sim import get_energy
Expand Down Expand Up @@ -61,7 +60,6 @@ def initialize(self, atoms):
self.is_initialized = True

def check(self, atoms: ase.Atoms) -> bool:

p1 = atoms.positions[self.idx_i]
p2 = atoms.positions[self.idx_j]
_, dists = conditional_find_mic(p1 - p2, atoms.cell, atoms.pbc)
Expand All @@ -80,7 +78,6 @@ def check(self, atoms: ase.Atoms) -> bool:
atoms.numbers[first_atom] = 3
atoms.numbers[second_atom] = 3


if self.bonded_max_dist:
max_dist = np.max(dists)
too_far = max_dist > self.bonded_max_dist
Expand Down
92 changes: 55 additions & 37 deletions ipsuite/analysis/model/plots.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,8 @@
import numpy as np
import seaborn as sns
from scipy.interpolate import interpn
from scipy.stats import gaussian_kde, foldnorm
from scipy.optimize import curve_fit
from scipy.stats import foldnorm, gaussian_kde


def density_scatter(ax, x, y, bins, **kwargs) -> None:
Expand Down Expand Up @@ -48,7 +48,13 @@ def density_scatter(ax, x, y, bins, **kwargs) -> None:


def get_figure(
true, prediction, datalabel: str, xlabel: str, ylabel: str, figsize: tuple = (10, 7), density=True,
true,
prediction,
datalabel: str,
xlabel: str,
ylabel: str,
figsize: tuple = (10, 7),
density=True,
) -> plt.Figure:
"""Create a correlation plot for true, prediction values.
Expand Down Expand Up @@ -104,31 +110,42 @@ def get_cdf_figure(x, y, figsize: tuple = (10, 7)):
return fig


def get_calibration_figure(error, std, markersize:float = 3.0, datalabel="", figsize: tuple = (10, 7)):
fig, ax = plt.subplots(1,1,figsize=figsize, dpi=300)
def get_calibration_figure(
error, std, markersize: float = 3.0, datalabel="", figsize: tuple = (10, 7)
):
fig, ax = plt.subplots(1, 1, figsize=figsize, dpi=300)

x = np.linspace(1e-6, 5e3, 5)
noise_level_2 = x

quantiles_lower_01 = [foldnorm.ppf(0.15, 0.,0.,i) for i in noise_level_2]
quantiles_upper_01 = [foldnorm.ppf(0.85, 0.,0.,i) for i in noise_level_2]
quantiles_lower_05 = [foldnorm.ppf(0.05, 0.,0.,i) for i in noise_level_2]
quantiles_upper_05 = [foldnorm.ppf(0.95, 0.,0.,i) for i in noise_level_2]
quantiles_lower_005 = [foldnorm.ppf(0.005, 0.,0.,i) for i in noise_level_2]
quantiles_upper_005 = [foldnorm.ppf(0.995, 0.,0.,i) for i in noise_level_2]

ax.scatter(std, error, s=markersize, alpha=0.3, color="tab:blue", rasterized=True, linewidth=0., label=datalabel)
quantiles_lower_01 = [foldnorm.ppf(0.15, 0.0, 0.0, i) for i in noise_level_2]
quantiles_upper_01 = [foldnorm.ppf(0.85, 0.0, 0.0, i) for i in noise_level_2]
quantiles_lower_05 = [foldnorm.ppf(0.05, 0.0, 0.0, i) for i in noise_level_2]
quantiles_upper_05 = [foldnorm.ppf(0.95, 0.0, 0.0, i) for i in noise_level_2]
quantiles_lower_005 = [foldnorm.ppf(0.005, 0.0, 0.0, i) for i in noise_level_2]
quantiles_upper_005 = [foldnorm.ppf(0.995, 0.0, 0.0, i) for i in noise_level_2]

ax.scatter(
std,
error,
s=markersize,
alpha=0.3,
color="tab:blue",
rasterized=True,
linewidth=0.0,
label=datalabel,
)
ax.loglog()
ax.plot(x, quantiles_upper_05, color='gray', alpha=0.5)
ax.plot(x, quantiles_lower_05, color='gray', alpha=0.5)
ax.plot(x, quantiles_upper_01, color='gray', alpha=0.5)
ax.plot(x, quantiles_lower_01, color='gray', alpha=0.5)
ax.plot(x, quantiles_upper_005, color='gray', alpha=0.5)
ax.plot(x, quantiles_lower_005, color='gray', alpha=0.5)
ax.plot(x, quantiles_upper_05, color="gray", alpha=0.5)
ax.plot(x, quantiles_lower_05, color="gray", alpha=0.5)
ax.plot(x, quantiles_upper_01, color="gray", alpha=0.5)
ax.plot(x, quantiles_lower_01, color="gray", alpha=0.5)
ax.plot(x, quantiles_upper_005, color="gray", alpha=0.5)
ax.plot(x, quantiles_lower_005, color="gray", alpha=0.5)

ax.plot(np.logspace(-3,100.0),np.logspace(-3,100.0),linestyle="--", color="grey")
ax.set_xlim(np.min(std)/1.5, np.max(std)*1.5)
ax.set_ylim(np.min(error)/1.5, np.max(error)*1.5)
ax.plot(np.logspace(-3, 100.0), np.logspace(-3, 100.0), linestyle="--", color="grey")
ax.set_xlim(np.min(std) / 1.5, np.max(std) * 1.5)
ax.set_ylim(np.min(error) / 1.5, np.max(error) * 1.5)

ax.set_xlabel(r"$\sigma_{f_{i\alpha}}(A)$ [eV/$\AA$] ")
ax.set_ylabel(r"$|\Delta f_{i\alpha}(A)|$ [eV/$\AA$] ")
Expand All @@ -138,34 +155,35 @@ def get_calibration_figure(error, std, markersize:float = 3.0, datalabel="", fig

def gauss(x, *p):
m, s = p
return np.exp(-((x-m)/s)**2*0.5)/np.sqrt(2*np.pi*s**2)
return np.exp(-(((x - m) / s) ** 2) * 0.5) / np.sqrt(2 * np.pi * s**2)


def slice_ensemble_uncertainty(true, pred_ens, slice_start, slice_end):
pred_mean = np.mean(pred_ens, axis=1)
pred_std = np.std(pred_ens, axis=1)
isel = np.where((slice_start<pred_std) & (pred_std<slice_end))[0]
error_true = np.reshape(true[isel,0]-pred_mean[isel,0], -1)
error_pred = np.reshape(pred_ens[isel,:,0]-pred_mean[isel,np.newaxis,0], -1)

isel = np.where((slice_start < pred_std) & (pred_std < slice_end))[0]

error_true = np.reshape(true[isel, 0] - pred_mean[isel, 0], -1)
error_pred = np.reshape(pred_ens[isel, :, 0] - pred_mean[isel, np.newaxis, 0], -1)
return error_true, error_pred

def slice_uncertainty(true, pred_mean, pred_std, slice_start, slice_end):
isel = np.where((slice_start<pred_std) & (pred_std<slice_end))[0]

error_true = np.reshape(true[isel]-pred_mean[isel], -1)

def slice_uncertainty(true, pred_mean, pred_std, slice_start, slice_end):
isel = np.where((slice_start < pred_std) & (pred_std < slice_end))[0]

error_true = np.reshape(true[isel] - pred_mean[isel], -1)
error_pred = pred_std[isel]
return error_true, error_pred


def get_gaussianicity_figure(error_true, error_pred):
true_kde_sel = gaussian_kde(error_true)
ens_kde_sel = gaussian_kde(error_pred)
ens_kde_sel = gaussian_kde(error_pred)

bounds = 1.5 * max(np.max(np.abs(error_true)), np.max(np.abs(error_pred)))

xgrid = np.linspace(-bounds,bounds,400)
xgrid = np.linspace(-bounds, bounds, 400)

ens_sel = ens_kde_sel(xgrid)
true_sel = true_kde_sel(xgrid)
Expand All @@ -175,11 +193,11 @@ def get_gaussianicity_figure(error_true, error_pred):

fig, ax = plt.subplots()

ax.semilogy(xgrid,gauss(xgrid,0,std), 'k--', label="Gaussian")
ax.semilogy(xgrid, true_sel, 'r-', label="empirical")
ax.semilogy(xgrid, ens_sel, 'b-', label="predicted")
ax.semilogy(xgrid, gauss(xgrid, 0, std), "k--", label="Gaussian")
ax.semilogy(xgrid, true_sel, "r-", label="empirical")
ax.semilogy(xgrid, ens_sel, "b-", label="predicted")
ymax = 5 * max(np.max(true_sel), np.max(ens_sel))
ax.set_ylim(1e-6,ymax)
ax.set_ylim(1e-6, ymax)
ax.set_yscale("log")

ax.set_xlabel(r"$\Delta (S)$ / meV/Ang")
Expand Down
23 changes: 16 additions & 7 deletions ipsuite/analysis/model/predict.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,14 @@

from ipsuite import base, models, utils
from ipsuite.analysis.model.math import decompose_stress_tensor, force_decomposition
from ipsuite.analysis.model.plots import get_cdf_figure, get_figure, get_hist, get_calibration_figure, get_gaussianicity_figure, slice_ensemble_uncertainty
from ipsuite.analysis.model.plots import (
get_calibration_figure,
get_cdf_figure,
get_figure,
get_gaussianicity_figure,
get_hist,
slice_ensemble_uncertainty,
)
from ipsuite.geometry import BarycenterMapping
from ipsuite.utils.ase_sim import freeze_copy_atoms

Expand Down Expand Up @@ -316,20 +323,22 @@ def get_plots(self, save=False):
)

gaussianicy_figures = []
for (start, end) in self.force_dist_slices:
error_true, error_pred = slice_ensemble_uncertainty(self.content["forces_true"], self.content["forces_ensemble"], start, end)
fig = get_gaussianicity_figure(
error_true, error_pred
for start, end in self.force_dist_slices:
error_true, error_pred = slice_ensemble_uncertainty(
self.content["forces_true"],
self.content["forces_ensemble"],
start,
end,
)
fig = get_gaussianicity_figure(error_true, error_pred)
gaussianicy_figures.append(fig)
if save:
forces_plot.savefig(self.plots_dir / "forces.png")
forces_cdf_plot.savefig(self.plots_dir / "forces_cdf.png")

for ii, fig in enumerate(gaussianicy_figures):
fig.savefig(self.plots_dir / f"forces_gaussianicity_{ii}.png")


def run(self):
self.nwd.mkdir(exist_ok=True, parents=True)
self.get_data()
Expand Down

0 comments on commit 24fcd39

Please sign in to comment.