Skip to content

Commit

Permalink
Merge branch 'update_evaluator_modules' into build_taproom
Browse files Browse the repository at this point in the history
  • Loading branch information
jeff231li committed Sep 18, 2023
2 parents c2099de + 65e4735 commit d248c25
Show file tree
Hide file tree
Showing 12 changed files with 1,239 additions and 427 deletions.
3 changes: 3 additions & 0 deletions devtools/conda-envs/test_env.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,10 @@ dependencies:
- plumed
- intermol
- openff-interchange >=0.3.7
<<<<<<< HEAD
- openff-toolkit >=0.11
=======
>>>>>>> update_evaluator_modules
- openff-units >=0.2.0
- openff-utilities

Expand Down
237 changes: 121 additions & 116 deletions paprika/analysis/analysis.py

Large diffs are not rendered by default.

227 changes: 123 additions & 104 deletions paprika/analysis/bootstrap.py

Large diffs are not rendered by default.

30 changes: 18 additions & 12 deletions paprika/analysis/utils.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
import numpy as np
import numpy


def get_factors(n):
Expand All @@ -17,7 +17,7 @@ def get_factors(n):
"""
factors = []
sqrt_n = int(round(np.sqrt(n) + 0.5))
sqrt_n = int(round(numpy.sqrt(n) + 0.5))
i = 1
while i <= sqrt_n:
if n % i == 0:
Expand All @@ -26,6 +26,7 @@ def get_factors(n):
if j != i:
factors.append(int(j))
i += 1

return sorted(factors, key=int)


Expand All @@ -45,19 +46,24 @@ def get_nearest_max(n):
"""
max_factors = 0

if n % 2 == 0:
beg = n - 100
end = n
else:
beg = n - 101
end = n - 1

if beg < 0:
beg = 0

most_factors = 0
for i in range(beg, end + 2, 2):
num_factors = len(get_factors(i))
if num_factors >= max_factors:
max_factors = num_factors
most_factors = i

return most_factors


Expand All @@ -72,12 +78,12 @@ def get_block_sem(data_array):
Parameters
----------
data_array: :class:`np.array`
data_array: :class:`numpy.array`
Array containing data values.
Returns
-------
np.max(sems): float
numpy.max(sems): float
The maximum SEM obtained from te blocking curve.
"""
Expand All @@ -86,11 +92,11 @@ def get_block_sem(data_array):
block_sizes = get_factors(len(data_array))

# An array to store means for each block ... make it bigger than we need.
block_means = np.zeros([block_sizes[-1]], np.float64)
block_means = numpy.zeros([block_sizes[-1]], numpy.float64)

# Store the SEM for each block size, except the last two size for which
# there will only be two or one blocks total and thus very noisy.
sems = np.zeros([len(block_sizes) - 2], np.float64)
sems = numpy.zeros([len(block_sizes) - 2], numpy.float64)

# Check each block size except the last two.
for size_idx in range(len(block_sizes) - 2):
Expand All @@ -102,15 +108,15 @@ def get_block_sem(data_array):
data_beg_idx = blk_idx * block_sizes[size_idx]
data_end_idx = (blk_idx + 1) * block_sizes[size_idx]
# Compute the mean of this block and store in array
block_means[blk_idx] = np.mean(data_array[data_beg_idx:data_end_idx])
block_means[blk_idx] = numpy.mean(data_array[data_beg_idx:data_end_idx])
# Compute the standard deviation across all blocks, devide by
# num_blocks-1 for SEM
sems[size_idx] = np.std(block_means[0:num_blocks], ddof=0) / np.sqrt(
sems[size_idx] = numpy.std(block_means[0:num_blocks], ddof=0) / numpy.sqrt(
num_blocks - 1
)
# Hmm or should ddof=1? I think 0, see Flyvbjerg -----^

return np.max(sems)
return numpy.max(sems)


def get_subsampled_indices(N, g, conservative=False):
Expand Down Expand Up @@ -138,16 +144,16 @@ def get_subsampled_indices(N, g, conservative=False):

# if conservative, assume integer g and round up
if conservative:
g = np.ceil(g)
g = numpy.ceil(g)

# initialize
indices = [0]
g_idx = 1.0
int_step = int(np.round(g_idx * g))
int_step = int(numpy.round(g_idx * g))

while int_step < N:
indices.append(int_step)
g_idx += 1.0
int_step = int(np.round(g_idx * g))
int_step = int(numpy.round(g_idx * g))

return indices
54 changes: 28 additions & 26 deletions paprika/io.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,13 +6,13 @@
import traceback
from enum import Enum

import numpy as np
import pytraj as pt
import numpy
import pytraj
from openff.units import unit as openff_unit
from parmed import Structure
from parmed.amber import AmberParm

from paprika.restraints import BiasPotentialType, DAT_restraint, RestraintType
from paprika.restraints import DAT_restraint

# https://stackoverflow.com/questions/27909658/json-encoder-and-decoder-for-complex-numpy-arrays
# https://stackoverflow.com/a/24375113/901925
Expand Down Expand Up @@ -80,11 +80,11 @@ def default(self, obj):
logging.warning("Encountered Structure, which does not store filename.")
return ""

if isinstance(obj, np.ndarray):
if isinstance(obj, numpy.ndarray):
if obj.flags["C_CONTIGUOUS"]:
obj_data = obj.data
else:
cont_obj = np.ascontiguousarray(obj)
cont_obj = numpy.ascontiguousarray(obj)
assert cont_obj.flags["C_CONTIGUOUS"]
obj_data = cont_obj.data
data_b64 = base64.b64encode(obj_data)
Expand All @@ -97,23 +97,25 @@ def default(self, obj):
elif isinstance(
obj,
(
np.int_,
np.intc,
np.intp,
np.int8,
np.int16,
np.int32,
np.int64,
np.uint8,
np.uint16,
np.uint32,
np.uint64,
numpy.int_,
numpy.intc,
numpy.intp,
numpy.int8,
numpy.int16,
numpy.int32,
numpy.int64,
numpy.uint8,
numpy.uint16,
numpy.uint32,
numpy.uint64,
),
):
return int(obj)
elif isinstance(obj, (np.float_, np.float16, np.float32, np.float64)):
elif isinstance(
obj, (numpy.float_, numpy.float16, numpy.float32, numpy.float64)
):
return float(obj)
elif isinstance(obj, (np.ndarray,)):
elif isinstance(obj, (numpy.ndarray,)):
return obj.tolist()
elif isinstance(obj, openff_unit.Quantity):
return serialize_quantity(obj)
Expand All @@ -139,7 +141,7 @@ def __init__(self, *args, **kwargs):
def custom_object_hook(self, obj):
if "__ndarray__" in obj:
data = base64.b64decode(obj["__ndarray__"])
return np.frombuffer(data, obj["dtype"]).reshape(obj["shape"])
return numpy.frombuffer(data, obj["dtype"]).reshape(obj["shape"])

if "@type" in obj:
if obj["@type"] == "openff.units.unit.Quantity":
Expand Down Expand Up @@ -343,18 +345,18 @@ def load_trajectory(window, trajectory, topology, single_topology=False):
)
logger.debug(f"Loading {os.path.join(window, topology)} and {trajectory_path}")
try:
traj = pt.iterload(trajectory_path, os.path.join(window, topology))
traj = pytraj.iterload(trajectory_path, os.path.join(window, topology))
except ValueError as e:
formatted_exception = traceback.format_exception(None, e, e.__traceback__)
logger.info(
f"Failed trying to load {os.path.join(window, topology)} and {trajectory_path}: "
f"{formatted_exception}"
)
elif isinstance(topology, str) and single_topology:
traj = pt.iterload(trajectory_path, os.path.join(topology))
traj = pytraj.iterload(trajectory_path, os.path.join(topology))
else:
try:
traj = pt.iterload(trajectory_path, topology)
traj = pytraj.iterload(trajectory_path, topology)
except BaseException:
raise Exception("Tried to load `topology` object directly and failed.")

Expand Down Expand Up @@ -384,7 +386,7 @@ def read_restraint_data(
Returns
-------
data: :class:`np.array`
data: :class:`numpy.array`
The values for this restraint in this window
"""

Expand All @@ -397,7 +399,7 @@ def read_restraint_data(
and not restraint.mask4
):
data = openff_unit.Quantity(
pt.distance(
pytraj.distance(
trajectory, " ".join([restraint.mask1, restraint.mask2]), image=True
),
units=openff_unit.angstrom,
Expand All @@ -407,7 +409,7 @@ def read_restraint_data(
restraint.mask1 and restraint.mask2 and restraint.mask3 and not restraint.mask4
):
data = openff_unit.Quantity(
pt.angle(
pytraj.angle(
trajectory,
" ".join([restraint.mask1, restraint.mask2, restraint.mask3]),
),
Expand All @@ -416,7 +418,7 @@ def read_restraint_data(

elif restraint.mask1 and restraint.mask2 and restraint.mask3 and restraint.mask4:
data = openff_unit.Quantity(
pt.dihedral(
pytraj.dihedral(
trajectory,
" ".join(
[restraint.mask1, restraint.mask2, restraint.mask3, restraint.mask4]
Expand Down
8 changes: 4 additions & 4 deletions paprika/simulate/amber.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
import abc
import logging
import os
import subprocess as sp
import subprocess
from collections import OrderedDict
from enum import Enum

Expand Down Expand Up @@ -633,11 +633,11 @@ def run(self, soft_minimize=False, overwrite=False, fail_ok=False):
logger.debug("Exec line: " + " ".join(exec_list))

# Execute
amber_output = sp.Popen(
amber_output = subprocess.Popen(
exec_list,
cwd=self.path,
stdout=sp.PIPE,
stderr=sp.PIPE,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
env=os.environ,
)

Expand Down
14 changes: 7 additions & 7 deletions paprika/simulate/gromacs.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
import glob
import logging
import os
import subprocess as sp
import subprocess
from collections import OrderedDict
from enum import Enum

Expand Down Expand Up @@ -623,11 +623,11 @@ def run(self, run_grompp=True, overwrite=False, fail_ok=False):
grompp_list += ["-n", self.index_file]

# Run GROMPP
grompp_output = sp.Popen(
grompp_output = subprocess.Popen(
grompp_list,
cwd=self.path,
stdout=sp.PIPE,
stderr=sp.PIPE,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
env=os.environ,
)
grompp_stdout = grompp_output.stdout.read().splitlines()
Expand Down Expand Up @@ -707,11 +707,11 @@ def run(self, run_grompp=True, overwrite=False, fail_ok=False):
mdrun_list += ["-plumed", self.plumed_file]

# Run MDRUN
mdrun_output = sp.Popen(
mdrun_output = subprocess.Popen(
mdrun_list,
cwd=self.path,
stdout=sp.PIPE,
stderr=sp.PIPE,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
env=os.environ,
)
mdrun_out = mdrun_output.stdout.read().splitlines()
Expand Down
16 changes: 8 additions & 8 deletions paprika/simulate/namd.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
import abc
import logging
import os
import subprocess as sp
import subprocess
from collections import OrderedDict
from enum import Enum

import numpy as np
import parmed as pmd
import numpy
import parmed

from paprika.utils import get_dict_without_keys

Expand Down Expand Up @@ -584,14 +584,14 @@ def _get_cell_basis_vectors(self):
"""
Function to calculate the PBC cell basis vectors (needed when running a simulation for the first time).
"""
structure = pmd.load_file(
structure = parmed.load_file(
os.path.join(self.path, self.topology),
os.path.join(self.path, self.coordinates),
structure=True,
)
coordinates = structure.coordinates
masses = np.ones(len(coordinates))
center = pmd.geometry.center_of_mass(coordinates, masses)
masses = numpy.ones(len(coordinates))
center = parmed.geometry.center_of_mass(coordinates, masses)

self.cell_basis_vectors["cellOrigin"] = list(center)

Expand Down Expand Up @@ -878,11 +878,11 @@ def run(self, overwrite=True, fail_ok=False):
logger.debug("Exec line: " + " ".join(exec_list))

# Execute
namd_output = sp.Popen(
namd_output = subprocess.Popen(
exec_list,
cwd=self.path,
stdout=open(os.path.join(self.path, self.logfile), "w"),
stderr=sp.PIPE,
stderr=subprocess.PIPE,
env=os.environ,
)
namd_stderr = namd_output.stderr.read().splitlines()
Expand Down
Loading

0 comments on commit d248c25

Please sign in to comment.