Skip to content

Commit

Permalink
black with 120 character limit; pyproject.toml
Browse files Browse the repository at this point in the history
  • Loading branch information
Max Imakaev committed Jan 7, 2023
1 parent df5d6a0 commit 7f9e6fb
Show file tree
Hide file tree
Showing 26 changed files with 196 additions and 607 deletions.
20 changes: 5 additions & 15 deletions examples/customIntegrators/activeBD.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,9 +27,7 @@
ids[int(N / 2) :] = 0


def run_monomer_diffusion(
gpuid, N, ids, activity_ratio, timestep=170, nblocks=10, blocksize=100
):
def run_monomer_diffusion(gpuid, N, ids, activity_ratio, timestep=170, nblocks=10, blocksize=100):
"""Run a single simulation on a GPU of a hetero-polymer with A monomers and B monomers. A monomers
have a larger diffusion coefficient than B monomers, with an activity ratio of D_A / D_B.
Expand All @@ -53,12 +51,8 @@ def run_monomer_diffusion(
"""
if len(ids) != N:
raise ValueError(
"The array of monomer identities must have length equal to the total number of monomers."
)
D = np.ones(
(N, 3)
) # Dx, Dy, Dz --> we assume the diffusion coefficient in each spatial dimension is the same
raise ValueError("The array of monomer identities must have length equal to the total number of monomers.")
D = np.ones((N, 3)) # Dx, Dy, Dz --> we assume the diffusion coefficient in each spatial dimension is the same
# by default set the average diffusion coefficient to be 1 kT/friction
# let D_A = 1 + Ddiff and D_B = 1 - Ddiff such that D_A / D_B is the given activity_ratio
Ddiff = (activity_ratio - 1) / (activity_ratio + 1)
Expand All @@ -79,9 +73,7 @@ def run_monomer_diffusion(
particleD = unit.Quantity(D, kT / friction)
integrator = ActiveBrownianIntegrator(timestep, collision_rate, particleD)
gpuid = f"{gpuid}"
reporter = HDF5Reporter(
folder="active_inactive", max_data_length=100, overwrite=True
)
reporter = HDF5Reporter(folder="active_inactive", max_data_length=100, overwrite=True)
sim = simulation.Simulation(
platform="CUDA",
# for custom integrators, feed a tuple with the integrator class reference and a string specifying type,
Expand All @@ -99,9 +91,7 @@ def run_monomer_diffusion(

polymer = starting_conformations.grow_cubic(N, int(np.ceil(r)))
sim.set_data(polymer, center=True) # loads a polymer, puts a center of mass at zero
sim.set_velocities(
v=np.zeros((N, 3))
) # initializes velocities of all monomers to zero (no inertia)
sim.set_velocities(v=np.zeros((N, 3))) # initializes velocities of all monomers to zero (no inertia)
sim.add_force(forces.spherical_confinement(sim, density=density, k=5.0))
sim.add_force(
forcekits.polymer_chains(
Expand Down
8 changes: 2 additions & 6 deletions examples/customIntegrators/corr_noise.py
Original file line number Diff line number Diff line change
Expand Up @@ -88,16 +88,12 @@ def run_correlated_diffusion(gpuid, N, rhos, timestep=170, nblocks=10, blocksize
"""
if rhos.shape[1] != N:
raise ValueError(
"The array of monomer identities must have length equal to the total number of monomers."
)
raise ValueError("The array of monomer identities must have length equal to the total number of monomers.")
# monomer density in confinement in units of monomers/volume
density = 0.224
r = (3 * N / (4 * 3.141592 * density)) ** (1 / 3)
print(f"Radius of confinement: {r}")
D = np.ones(
(N, 3)
) # Diffusion coefficients of N monomers in x,y,z spatial dimensions
D = np.ones((N, 3)) # Diffusion coefficients of N monomers in x,y,z spatial dimensions
timestep = timestep
# the monomer diffusion coefficient should be in units of kT / friction, where friction = mass*collision_rate
collision_rate = 2.0
Expand Down
12 changes: 3 additions & 9 deletions examples/loopExtrusion/directionalModel/flagshipNormLifetime.py
Original file line number Diff line number Diff line change
Expand Up @@ -70,9 +70,7 @@ def averageContacts(contactFunction, inValues, N, **kwargs):

filenameChunks = [filenames[i::nproc] for i in range(nproc)]

with closing(
mp.Pool(processes=nproc, initializer=init, initargs=sharedARrays_)
) as p:
with closing(mp.Pool(processes=nproc, initializer=init, initargs=sharedARrays_)) as p:
p.map(worker, filenameChunks)


Expand Down Expand Up @@ -201,9 +199,7 @@ def doSim(i):
hicdata = np.clip(hicdata, 0, np.percentile(hicdata, 99.99))
hicdata /= np.mean(np.sum(hicdata, axis=1))

fmap(
doSim, range(30), n=20
) # number of threads to use. On a 20-core machine I use 20.
fmap(doSim, range(30), n=20) # number of threads to use. On a 20-core machine I use 20.

arr = coarsegrain(arr, 20)
arr = np.clip(arr, 0, np.percentile(arr, 99.9))
Expand Down Expand Up @@ -322,7 +318,5 @@ def doPolymerSimulation(steps, dens, stiff, folder):
steps=5000,
stiff=2,
dens=0.2,
folder="flagshipLifetime{1}Mu{2}_try={0}".format(
sys.argv[1], sys.argv[2], sys.argv[3]
),
folder="flagshipLifetime{1}Mu{2}_try={0}".format(sys.argv[1], sys.argv[2], sys.argv[3]),
)
12 changes: 2 additions & 10 deletions examples/loopExtrusion/directionalModel/showContactmap.py
Original file line number Diff line number Diff line change
Expand Up @@ -365,16 +365,8 @@ def showCmapNew():
if end > 75000:
continue
plt.xlim([st, end])
plt.savefig(
os.path.join(
"heatmaps", "{0}_st={1}_end={2}_r=2.png".format(fname, st, end)
)
)
plt.savefig(
os.path.join(
"heatmaps", "{0}_st={1}_end={2}_r=2.pdf".format(fname, st, end)
)
)
plt.savefig(os.path.join("heatmaps", "{0}_st={1}_end={2}_r=2.png".format(fname, st, end)))
plt.savefig(os.path.join("heatmaps", "{0}_st={1}_end={2}_r=2.pdf".format(fname, st, end)))
plt.clf()

plt.show()
Expand Down
58 changes: 13 additions & 45 deletions examples/loopExtrusion/newSimulatorCode.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,9 +31,7 @@

# new parameters because some things changed
saveEveryBlocks = 10 # save every 10 blocks (saving every block is now too much almost)
skipSavedBlocksBeginning = (
20 # how many blocks (saved) to skip after you restart LEF positions
)
skipSavedBlocksBeginning = 20 # how many blocks (saved) to skip after you restart LEF positions
totalSavedBlocks = 40 # how many blocks to save (number of blocks done is totalSavedBlocks * saveEveryBlocks)
restartMilkerEveryBlocks = 100

Expand All @@ -49,16 +47,8 @@

savesPerMilker = restartMilkerEveryBlocks // saveEveryBlocks
milkerInitsSkip = saveEveryBlocks * skipSavedBlocksBeginning // restartMilkerEveryBlocks
milkerInitsTotal = (
(totalSavedBlocks + skipSavedBlocksBeginning)
* saveEveryBlocks
// restartMilkerEveryBlocks
)
print(
"Milker will be initialized {0} times, first {1} will be skipped".format(
milkerInitsTotal, milkerInitsSkip
)
)
milkerInitsTotal = (totalSavedBlocks + skipSavedBlocksBeginning) * saveEveryBlocks // restartMilkerEveryBlocks
print("Milker will be initialized {0} times, first {1} will be skipped".format(milkerInitsTotal, milkerInitsSkip))


class smcTranslocatorMilker(object):
Expand Down Expand Up @@ -93,9 +83,7 @@ def setup(self, bondForce, blocks=100, smcStepsPerBlock=1):
"""

if len(self.allBonds) != 0:
raise ValueError(
"Not all bonds were used; {0} sets left".format(len(self.allBonds))
)
raise ValueError("Not all bonds were used; {0} sets left".format(len(self.allBonds)))

self.bondForce = bondForce

Expand All @@ -115,11 +103,7 @@ def setup(self, bondForce, blocks=100, smcStepsPerBlock=1):
self.curBonds = allBonds.pop(0)

for bond in self.uniqueBonds:
paramset = (
self.activeParamDict
if (bond in self.curBonds)
else self.inactiveParamDict
)
paramset = self.activeParamDict if (bond in self.curBonds) else self.inactiveParamDict
ind = bondForce.addBond(bond[0], bond[1], **paramset)
self.bondInds.append(ind)
self.bondToInd = {i: j for i, j in zip(self.uniqueBonds, self.bondInds)}
Expand All @@ -133,9 +117,7 @@ def step(self, context, verbose=False):
:return: (current bonds, previous step bonds); just for reference
"""
if len(self.allBonds) == 0:
raise ValueError(
"No bonds left to run; you should restart simulation and run setup again"
)
raise ValueError("No bonds left to run; you should restart simulation and run setup again")

pastBonds = self.curBonds
self.curBonds = self.allBonds.pop(0) # getting current bonds
Expand All @@ -153,12 +135,8 @@ def step(self, context, verbose=False):
for bond, isAdd in zip(bondsToChange, bondsIsAdd):
ind = self.bondToInd[bond]
paramset = self.activeParamDict if isAdd else self.inactiveParamDict
self.bondForce.setBondParameters(
ind, bond[0], bond[1], **paramset
) # actually updating bonds
self.bondForce.updateParametersInContext(
context
) # now run this to update things in the context
self.bondForce.setBondParameters(ind, bond[0], bond[1], **paramset) # actually updating bonds
self.bondForce.updateParametersInContext(context) # now run this to update things in the context
return self.curBonds, pastBonds


Expand Down Expand Up @@ -195,9 +173,7 @@ def initModel():
# now polymer simulation code starts

# ------------feed smcTran to the milker---
SMCTran.steps(
1000000
) # first steps to "equilibrate" SMC dynamics. If desired of course.
SMCTran.steps(1000000) # first steps to "equilibrate" SMC dynamics. If desired of course.
milker = smcTranslocatorMilker(SMCTran) # now feed this thing to milker (do it once!)
# --------- end new code ------------

Expand All @@ -206,9 +182,7 @@ def initModel():

# simulation parameters are defined below
a = Simulation(timestep=80, thermostat=0.01)
a.setup(
platform="CUDA", PBC=True, PBCbox=[box, box, box], GPU=0, precision="mixed"
) # set up GPU here
a.setup(platform="CUDA", PBC=True, PBCbox=[box, box, box], GPU=0, precision="mixed") # set up GPU here
a.saveFolder(folder)
a.load(data)
a.addHarmonicPolymerBonds(wiggleDist=0.1)
Expand All @@ -234,13 +208,9 @@ def initModel():
) # now only one step of SMC per step
print("Restarting milker")

a.doBlock(
steps=steps, increment=False
) # do block for the first time with first set of bonds in
a.doBlock(steps=steps, increment=False) # do block for the first time with first set of bonds in
for i in range(restartMilkerEveryBlocks - 1):
curBonds, pastBonds = milker.step(
a.context
) # this updates bonds. You can do something with bonds here
curBonds, pastBonds = milker.step(a.context) # this updates bonds. You can do something with bonds here
if i % saveEveryBlocks == (saveEveryBlocks - 2):
a.doBlock(steps=steps, increment=doSave)
if doSave:
Expand All @@ -250,9 +220,7 @@ def initModel():
open(os.path.join(a.folder, "SMC{0}.dat".format(a.step)), "wb"),
)
else:
a.integrator.step(
steps
) # do steps without getting the positions from the GPU (faster)
a.integrator.step(steps) # do steps without getting the positions from the GPU (faster)

data = a.getData() # save data and step, and delete the simulation
block = a.step
Expand Down
33 changes: 8 additions & 25 deletions polychrom/cli/traj_convert.py
Original file line number Diff line number Diff line change
Expand Up @@ -103,9 +103,7 @@ def _find_matches(pat, filenames):
a = re.search(pat, filename)
if a is not None:
if len(a.groups()) != 1:
raise ValueError(
"You should have one group in regex denoting the number of the file"
)
raise ValueError("You should have one group in regex denoting the number of the file")
assert len(a.groups()) == 1
gr = int(a.groups()[0])
result[filename] = gr
Expand Down Expand Up @@ -176,12 +174,8 @@ def _find_matches(pat, filenames):
help="allow blocks to be non-consecutive (1,2,3...)",
)
@click.option("--verbose", is_flag=True)
@click.option(
"--round-to", default=2, show_default=True, help="round to this number of digits"
)
@click.option(
"--skip-files", default=1, show_default=True, help="save only every Nth file"
)
@click.option("--round-to", default=2, show_default=True, help="round to this number of digits")
@click.option("--skip-files", default=1, show_default=True, help="save only every Nth file")
@click.option(
"--HDF5-blocks-per-file",
default=100,
Expand Down Expand Up @@ -251,10 +245,7 @@ def trajcopy(
if kwargs["input_style"] == "old":
blocks = _find_matches(block_pattern, all_files)
elif kwargs["input_style"] == "new":
blocks = {
i: j
for j, i in list_URIs(in_dir, empty_error=True, return_dict=True).items()
}
blocks = {i: j for j, i in list_URIs(in_dir, empty_error=True, return_dict=True).items()}
else:
raise ValueError("input-style should be 'old' or 'new'")

Expand Down Expand Up @@ -300,15 +291,11 @@ def trajcopy(
assert len(extra_loader) == len(extra_require)

# matching patterns for extra files, populating the dataframe
for val_pat, val_name, require in zip(
extra_pattern, extra_pattern_name, extra_require
):
for val_pat, val_name, require in zip(extra_pattern, extra_pattern_name, extra_require):
datas = _find_matches(val_pat, all_files)
if require:
if len(datas) != len(blocks):
raise ValueError(
f"files missing for {val_name}: need {len(blocks)} found {len(datas)}"
)
raise ValueError(f"files missing for {val_name}: need {len(blocks)} found {len(datas)}")
if len(datas) > 0:
datas = pd.Series(data=list(datas.keys()), index=list(datas.values()))
datas.name = val_name
Expand All @@ -327,9 +314,7 @@ def trajcopy(
print(other[:: len(other) // 20 + 1])
print("Verify that none of these should be converted using extra_pattern")
print("If not, increase max_unmatched_files")
raise ValueError(
"Limit exceeded: {0} files did not match anything".format(len(other))
)
raise ValueError("Limit exceeded: {0} files did not match anything".format(len(other)))

# creating the reporter
if (len(blocks) > 0) and (not kwargs["dry_run"]):
Expand Down Expand Up @@ -361,9 +346,7 @@ def trajcopy(
cur["block"] = i
elif kwargs["input_style"] == "new":
cur = load_URI(data)
cur["pos"] = np.round(
np.asarray(cur["pos"], dtype=np.float32), kwargs["round_to"]
)
cur["pos"] = np.round(np.asarray(cur["pos"], dtype=np.float32), kwargs["round_to"])

# adding "extra" data in the dict to save
for name, ldr in zip(extra_pattern_name, extra_loader):
Expand Down
Loading

0 comments on commit 7f9e6fb

Please sign in to comment.