Core
PASCal.core
¶
This module defines the core functionality of PASCal: the fit function and the results container.
Classes¶
PASCalResults
dataclass
¶
A container for the results of a PASCal run, providing convenience wrappers to all data objects and plotting functions.
Source code in PASCal/core.py
@dataclass
class PASCalResults:
"""A container for the results of a PASCal run, providing
convenience wrappers to all data objects and plotting functions.
"""
options: Options
"""The provided fit options."""
x: Union[Pressure, Temperature, Charge]
"""The input control variable $X$, either pressure, temperature or charge."""
x_errors: Union[Pressure, Temperature, Charge]
"""The input errors on the control variable $X$, either pressure, temperature or charge."""
unit_cells: np.ndarray
"""The input unit cells."""
diagonal_strain: Strain
"""The diagonalised strain."""
cell_volumes: Volume
"""The computed cell volumes."""
principal_components: np.ndarray
"""The eigenvalues of the strain."""
indicatrix: Tuple[Any, Any, Any, Any, Any]
"""The indicatrix data used to construct the surface plot."""
norm_crax: np.ndarray
"""The normalised crystal axes."""
median_x: int
"""The index of the middle value of the control variable $X$."""
median_principal_axis_crys: np.ndarray
"""The median principal axis of the strain."""
principal_axis_crys: np.ndarray
"""The principal axes of the strain."""
strain_fits: Dict[str, List[Any]]
"""A dictionary containing the strain fits. Under each key is a list of fit models for each principal axis,
containing either a statsmodel fit result (temperature data), a tuple of `popt`, `pcov`, `callable` results from SciPy's `curve_fit` (pressure data),
or an array of Chebyshev coefficients from NumPy's `chebfit` (electrochemical data).
"""
volume_fits: Dict[str, Any]
"""A dictionary containing the volume fits. Under each key is a model for the cell volume,
containing either a statsmodel fit result, a tuple of `popt`, `pcov`, `callable` results from SciPy's `curve_fit`,
or an array of Chebyshev coefficients from NumPy's `chebfit`.
"""
compressibility: Optional[np.ndarray]
"""The computed compressibility."""
compressibility_errors: Optional[np.ndarray]
"""The computed compressibility errors."""
warning: List[str]
"""Any warnings generated during the fit."""
named_coefficients: Dict[str, Any] = field(default_factory=dict)
"""Any additional named coefficients to render in the table."""
def plot_strain(
self, return_json: bool = False, show_errors: bool = False
) -> Union[str, plotly.graph_objs.Figure]:
"""Plots the strain fits."""
return plot_strain(
self.x,
self.x_errors,
self.diagonal_strain,
self.strain_fits,
self.options.data_type,
return_json=return_json,
show_errors=show_errors,
)
def plot_volume(
self, return_json: bool = False, show_errors: bool = False
) -> Union[str, plotly.graph_objs.Figure]:
return plot_volume(
self.x,
self.x_errors,
self.cell_volumes,
self.volume_fits,
self.options.data_type,
return_json=return_json,
show_errors=show_errors,
)
def plot_indicatrix(
self,
return_json: bool = False,
plot_size: int = 800,
) -> Union[str, plotly.graph_objs.Figure]:
return plot_indicatrix(
self.norm_crax,
*self.indicatrix,
self.options.data_type,
return_json=return_json,
plot_size=plot_size,
)
def plot_compressibility(
self, return_json: bool = False
) -> Union[str, plotly.graph_objs.Figure]:
assert (
self.compressibility is not None
), "Cannot plot before compresibility calculation."
assert (
self.compressibility_errors is not None
), "Cannot plot before compresibility calculation."
return plot_compressibility(
self.x,
self.compressibility,
self.compressibility_errors,
self.strain_fits,
self.options.data_type,
return_json=return_json,
)
def plot_charge_derivative(
self, return_json: bool = False
) -> Union[str, plotly.graph_objs.Figure]:
return plot_charge_derivative(
self.x,
self.diagonal_strain,
self.strain_fits,
self.options.data_type,
return_json=return_json,
)
def plot_residual(
self, return_json: bool = False
) -> Union[str, plotly.graph_objs.Figure]:
return plot_residual(
self.strain_fits,
self.volume_fits,
self.options.data_type,
return_json=return_json,
)
def _set_named_coeffs(self):
"""Compute a series of named coefficients for displaying in the app."""
self.named_coefficients = {}
if self.options.data_type == PASCalDataType.TEMPERATURE:
self.named_coefficients["CalAlphaErr"] = np.array(
[self.strain_fits["linear"][i].HC0_se[1] * K_to_MK for i in range(3)]
)
self.named_coefficients["VolLin"] = np.array(
self.volume_fits["linear"].params[1] * self.x
+ self.volume_fits["linear"].params[0],
)
self.named_coefficients["VolCoef"] = (
self.volume_fits["linear"].params[1] / self.cell_volumes[0] * K_to_MK
)
self.named_coefficients["VolCoefErr"] = (
self.volume_fits["linear"].HC0_se[1] / self.cell_volumes[0] * K_to_MK
)
self.named_coefficients["XCal"] = np.array(
[
PERCENT
* (
self.strain_fits["linear"][i].params[1] * self.x
+ self.strain_fits["linear"][i].params[0]
)
for i in range(3)
]
)
if self.options.data_type == PASCalDataType.PRESSURE:
self.named_coefficients["VolCoef"] = -(
self.volume_fits["linear"].params[1] / self.cell_volumes[0] * GPa_to_TPa
)
self.named_coefficients["VolCoefErr"] = (
self.volume_fits["linear"].HC0_se[1] / self.cell_volumes[0] * GPa_to_TPa
)
self.named_coefficients["CalEpsilon0"] = np.array(
[self.strain_fits["empirical"][0][i][0] for i in range(3)]
)
self.named_coefficients["CalLambda"] = np.array(
[self.strain_fits["empirical"][0][i][1] for i in range(3)]
)
self.named_coefficients["CalPc"] = np.array(
[self.strain_fits["empirical"][0][i][2] for i in range(3)]
)
self.named_coefficients["CalNu"] = np.array(
[self.strain_fits["empirical"][0][i][3] for i in range(3)]
)
self.named_coefficients["XCal"] = np.array(
[
PERCENT
* (
empirical_pressure_strain_relation(
self.x, *self.strain_fits["empirical"][0][i]
)
)
for i in range(3)
]
)
self.named_coefficients["B0"] = np.array(
[self.volume_fits[k][0][1] for k in self.volume_fits if k != "linear"]
)
sigmas = [
np.sqrt(np.diag(self.volume_fits[k][1]))
for k in self.volume_fits
if k != "linear"
]
self.named_coefficients["SigB0"] = np.array([d[1] for d in sigmas])
self.named_coefficients["V0"] = np.array(
[self.volume_fits[k][0][0] for k in self.volume_fits if k != "linear"]
)
self.named_coefficients["SigV0"] = np.array([d[0] for d in sigmas])
self.named_coefficients["BPrime"] = np.array(
[
self.volume_fits[k][0][2] if len(self.volume_fits[k][0]) > 2 else 4
for k in self.volume_fits
if k != "linear"
]
)
self.named_coefficients["SigBPrime"] = np.array(
[d[2] if len(d) > 2 else "n/a" for d in sigmas], dtype=object
)
self.named_coefficients["PcCoef"] = np.array(
[0.0, 0.0, self.options.pc_val if self.options.pc_val else 0.0]
)
self.named_coefficients["CalPress"] = np.zeros(
(len(self.volume_fits), len(self.cell_volumes))
)
self.named_coefficients["CalPress"][0][:] = (
self.cell_volumes - self.volume_fits["linear"].params[0]
) / self.volume_fits["linear"].params[1]
axis = 1
for fn in self.volume_fits:
if fn != "linear":
self.named_coefficients["CalPress"][axis][:] = fn(
self.cell_volumes, *self.volume_fits[fn][0]
)
axis += 1
if self.options.data_type == PASCalDataType.ELECTROCHEMICAL:
best_degrees, _ = get_best_chebyshev_strain_fit(
self.strain_fits["chebyshev"]
)
self.named_coefficients["XCal"] = np.array(
[
np.polynomial.chebyshev.chebval(
self.x, self.strain_fits["chebyshev"][best_degrees[i]][0][i]
)
for i in range(3)
]
)
cheby_deriv = [
np.polynomial.chebyshev.chebder(
self.strain_fits["chebyshev"][best_degrees[i]][0][i],
m=1,
scl=1,
axis=0,
)
for i in range(3)
]
self.named_coefficients["Deriv"] = np.array(
[
mAhg_to_kAhg
* np.polynomial.chebyshev.chebval(self.x, cheby_deriv[i])
for i in range(3)
]
)
best_degree, vol_coeff = get_best_chebyshev_volume_fit(
self.volume_fits["chebyshev"]
)
self.named_coefficients["VolCheb"] = np.polynomial.chebyshev.chebval(
self.x, self.volume_fits["chebyshev"][best_degree][0]
)
vol_der = np.polynomial.chebyshev.chebder(vol_coeff, m=1, scl=1, axis=0)
self.named_coefficients["VolCoef"] = (
np.polynomial.chebyshev.chebval(self.x, vol_der)[self.median_x]
* mAhg_to_kAhg
)
Attributes¶
options: Options
dataclass-field
¶
The provided fit options.
x: ndarray
dataclass-field
¶
The input control variable \(X\), either pressure, temperature or charge.
x_errors: ndarray
dataclass-field
¶
The input errors on the control variable \(X\), either pressure, temperature or charge.
unit_cells: ndarray
dataclass-field
¶
The input unit cells.
diagonal_strain: ndarray
dataclass-field
¶
The diagonalised strain.
cell_volumes: ndarray
dataclass-field
¶
The computed cell volumes.
principal_components: ndarray
dataclass-field
¶
The eigenvalues of the strain.
indicatrix: Tuple[Any, Any, Any, Any, Any]
dataclass-field
¶
The indicatrix data used to construct the surface plot.
norm_crax: ndarray
dataclass-field
¶
The normalised crystal axes.
median_x: int
dataclass-field
¶
The index of the middle value of the control variable \(X\).
median_principal_axis_crys: ndarray
dataclass-field
¶
The median principal axis of the strain.
principal_axis_crys: ndarray
dataclass-field
¶
The principal axes of the strain.
strain_fits: Dict[str, List[Any]]
dataclass-field
¶
A dictionary containing the strain fits. Under each key is a list of fit models for each principal axis,
containing either a statsmodel fit result (temperature data), a tuple of popt
, pcov
, callable
results from SciPy's curve_fit
(pressure data),
or an array of Chebyshev coefficients from NumPy's chebfit
(electrochemical data).
volume_fits: Dict[str, Any]
dataclass-field
¶
A dictionary containing the volume fits. Under each key is a model for the cell volume,
containing either a statsmodel fit result, a tuple of popt
, pcov
, callable
results from SciPy's curve_fit
,
or an array of Chebyshev coefficients from NumPy's chebfit
.
compressibility: Optional[numpy.ndarray]
dataclass-field
¶
The computed compressibility.
compressibility_errors: Optional[numpy.ndarray]
dataclass-field
¶
The computed compressibility errors.
warning: List[str]
dataclass-field
¶
Any warnings generated during the fit.
named_coefficients: Dict[str, Any]
dataclass-field
¶
Any additional named coefficients to render in the table.
Methods¶
plot_strain(self, return_json: bool = False, show_errors: bool = False) -> Union[str, plotly.graph_objs._figure.Figure]
¶
Plots the strain fits.
Source code in PASCal/core.py
def plot_strain(
self, return_json: bool = False, show_errors: bool = False
) -> Union[str, plotly.graph_objs.Figure]:
"""Plots the strain fits."""
return plot_strain(
self.x,
self.x_errors,
self.diagonal_strain,
self.strain_fits,
self.options.data_type,
return_json=return_json,
show_errors=show_errors,
)
plot_volume(self, return_json: bool = False, show_errors: bool = False) -> Union[str, plotly.graph_objs._figure.Figure]
¶
Source code in PASCal/core.py
def plot_volume(
self, return_json: bool = False, show_errors: bool = False
) -> Union[str, plotly.graph_objs.Figure]:
return plot_volume(
self.x,
self.x_errors,
self.cell_volumes,
self.volume_fits,
self.options.data_type,
return_json=return_json,
show_errors=show_errors,
)
plot_indicatrix(self, return_json: bool = False, plot_size: int = 800) -> Union[str, plotly.graph_objs._figure.Figure]
¶
Source code in PASCal/core.py
def plot_indicatrix(
self,
return_json: bool = False,
plot_size: int = 800,
) -> Union[str, plotly.graph_objs.Figure]:
return plot_indicatrix(
self.norm_crax,
*self.indicatrix,
self.options.data_type,
return_json=return_json,
plot_size=plot_size,
)
plot_compressibility(self, return_json: bool = False) -> Union[str, plotly.graph_objs._figure.Figure]
¶
Source code in PASCal/core.py
def plot_compressibility(
self, return_json: bool = False
) -> Union[str, plotly.graph_objs.Figure]:
assert (
self.compressibility is not None
), "Cannot plot before compresibility calculation."
assert (
self.compressibility_errors is not None
), "Cannot plot before compresibility calculation."
return plot_compressibility(
self.x,
self.compressibility,
self.compressibility_errors,
self.strain_fits,
self.options.data_type,
return_json=return_json,
)
plot_charge_derivative(self, return_json: bool = False) -> Union[str, plotly.graph_objs._figure.Figure]
¶
Source code in PASCal/core.py
def plot_charge_derivative(
self, return_json: bool = False
) -> Union[str, plotly.graph_objs.Figure]:
return plot_charge_derivative(
self.x,
self.diagonal_strain,
self.strain_fits,
self.options.data_type,
return_json=return_json,
)
plot_residual(self, return_json: bool = False) -> Union[str, plotly.graph_objs._figure.Figure]
¶
Source code in PASCal/core.py
def plot_residual(
self, return_json: bool = False
) -> Union[str, plotly.graph_objs.Figure]:
return plot_residual(
self.strain_fits,
self.volume_fits,
self.options.data_type,
return_json=return_json,
)
Functions¶
fit(x, x_errors, unit_cells, options: Union[PASCal.options.Options, dict]) -> PASCalResults
¶
Perform the PASCal fits for the given data.
For temperature data, linear models are fitted for strain vs temperature and volume vs temperature. For pressure data, an empirical fit is made for strain vs pressure, and Birch-Murnaghan fits are made for volume vs pressure. For electrochemical data, Chebyshev polynomials are fitted for strain vs state of charge and volume vs state of charge.
Parameters:
Name | Type | Description | Default |
---|---|---|---|
x |
The independent variable (pressure, temperature or charge). |
required | |
x_errors |
The errors on the independent variable. |
required | |
unit_cells |
The unit cell volumes. |
required | |
options |
Union[PASCal.options.Options, dict] |
The options for the fit. |
required |
Returns:
Type | Description |
---|---|
PASCalResults |
A PASCalResults object containing the results of the fit. |
Source code in PASCal/core.py
def fit(x, x_errors, unit_cells, options: Union[Options, dict]) -> PASCalResults:
"""Perform the PASCal fits for the given data.
For temperature data, linear models are fitted for strain vs temperature and volume
vs temperature.
For pressure data, an empirical fit is made for strain vs pressure, and
Birch-Murnaghan fits are made for volume vs pressure.
For electrochemical data, Chebyshev polynomials are fitted for strain vs state
of charge and volume vs state of charge.
Parameters:
x: The independent variable (pressure, temperature or charge).
x_errors: The errors on the independent variable.
unit_cells: The unit cell volumes.
options: The options for the fit.
Returns:
A PASCalResults object containing the results of the fit.
"""
if not isinstance(options, Options):
options = Options.from_dict(options)
warning = options.precheck_inputs(x)
print(f"Performing fit with {options=}")
cell_volumes = PASCal.utils.cell_vols(unit_cells)
orthonormed_cells = PASCal.utils.orthomat(unit_cells) # cell in orthogonal axes
strain = PASCal.utils.calculate_strain(
orthonormed_cells,
mode="eulerian" if options.eulerian_strain else "lagrangian",
finite=options.finite_strain,
)
# Diagonalising to get eigenvectors and values, principal_axes are in orthogonal coordinates
diagonal_strain, principal_axes = np.linalg.eigh(strain)
median_x = int(np.ceil(len(x) / 2)) - 1 # median
### Axes matching
principal_axes, diagonal_strain = PASCal.utils.match_axes(
principal_axes, orthonormed_cells, diagonal_strain
)
### Calculating Eigenvectors and Cells in different coordinate systems
(
median_principal_axis_crys,
principal_axis_crys,
crys_prin_ax,
) = PASCal.utils.rotate_to_principal_axes(
orthonormed_cells, principal_axes, median_x
)
strain_fits: dict = {}
volume_fits: dict = {}
strain_fits["linear"] = fit_linear_wls(diagonal_strain, x, x_errors)
volume_fits["linear"] = fit_linear_wls(cell_volumes, x, x_errors)[0]
if options.data_type == PASCalDataType.TEMPERATURE:
principal_components = np.array(
[strain_fits["linear"][i].params[1] * K_to_MK for i in range(3)]
)
elif options.data_type == PASCalDataType.PRESSURE:
# do empirical fits
strain_fits["empirical"] = fit_empirical_strain_pressure(
diagonal_strain,
x,
x_errors,
strain_fits["linear"],
)
compressibility = np.zeros((3, len(x)))
compressibility_errors = np.zeros((3, len(x)))
popt, pcov, _ = strain_fits["empirical"]
for i in range(3):
empirical_popts = popt[i][1:]
compressibility[i][:] = (
PASCal.utils.get_compressibility(x, *empirical_popts) * GPa_to_TPa
)
empirical_pcovs = pcov[i]
compressibility_errors[i][:] = (
PASCal.utils.get_compressibility_errors(
empirical_pcovs, x, *empirical_popts
)
* GPa_to_TPa
)
principal_components = np.array(
[compressibility[i][median_x] for i in range(3)]
)
bm_popts, bm_pcovs = fit_birch_murnaghan_volume_pressure(
cell_volumes,
x,
x_errors,
critical_pressure=options.pc_val if options.use_pc else None,
)
for k in bm_popts:
volume_fits[k] = bm_popts[k], bm_pcovs[k]
elif options.data_type == PASCalDataType.ELECTROCHEMICAL:
strain_fits["chebyshev"] = {}
for deg in range(1, options.deg_poly_strain + 1):
strain_fits["chebyshev"][deg] = fit_chebyshev(
diagonal_strain,
x,
deg,
)
volume_fits["chebyshev"] = {}
for deg in range(1, options.deg_poly_vol + 1):
volume_fits["chebyshev"][deg] = fit_chebyshev(cell_volumes, x, deg)
best_degrees, best_coeffs = get_best_chebyshev_strain_fit(
strain_fits["chebyshev"]
)
cheby_deriv_coeffs = [
np.polynomial.chebyshev.chebder(
strain_fits["chebyshev"][best_degrees[i]][0][i],
m=1,
scl=1,
axis=0,
)
for i in range(3)
]
deriv = [
mAhg_to_kAhg * np.polynomial.chebyshev.chebval(x, cheby_deriv_coeffs[i])
for i in range(3)
]
principal_components = np.array([deriv[i][median_x] for i in range(3)])
norm_crax = PASCal.utils.normalise_crys_axes(
crys_prin_ax[median_x, :, :], principal_components
)
indicatrix = PASCal.utils.indicatrix(principal_components)
results = PASCalResults(
options=options,
x=x,
x_errors=x_errors,
unit_cells=unit_cells,
principal_components=principal_components,
median_x=median_x,
diagonal_strain=diagonal_strain,
cell_volumes=cell_volumes,
strain_fits=strain_fits,
volume_fits=volume_fits,
indicatrix=indicatrix,
norm_crax=norm_crax,
median_principal_axis_crys=median_principal_axis_crys,
principal_axis_crys=principal_axis_crys,
warning=warning,
compressibility=compressibility
if options.data_type == PASCalDataType.PRESSURE
else None,
compressibility_errors=compressibility_errors
if options.data_type == PASCalDataType.PRESSURE
else None,
)
results._set_named_coeffs()
return results