Skip to content

Commit

Permalink
Merge pull request #33 from elbeejay/kwarg-particleiterationlimit
Browse files Browse the repository at this point in the history
make max_iter kwarg
  • Loading branch information
elbeejay authored Apr 12, 2022
2 parents 9aa3ed7 + 6eeca40 commit 179105d
Show file tree
Hide file tree
Showing 4 changed files with 96 additions and 38 deletions.
2 changes: 1 addition & 1 deletion dorado/__init__.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
__version__ = "2.5.0"
__version__ = "2.5.1"


from . import lagrangian_walker
Expand Down
22 changes: 11 additions & 11 deletions dorado/lagrangian_walker.py
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ def big_sliding_window(raster):
index, in the order [NW, N, NE, W, 0, E, SW, S, SE]. Outputs are
ordered to match np.ravel(), so it functions similarly to a loop
applying ravel to the elements around each index.
For example, the neighboring values in raster indexed at (i,j) are
For example, the neighboring values in raster indexed at (i,j) are
raster(i-1:i+2, j-1:j+2).ravel(). These 9 values have been mapped to
big_ravel(i,j,:) for ease of computations. Helper function for make_weight.
Expand All @@ -59,7 +59,7 @@ def big_sliding_window(raster):
**Outputs** :
big_ravel : `ndarray`
3D array which sorts the D8 neighbors at index (i,j) in
3D array which sorts the D8 neighbors at index (i,j) in
raster into the 3rd dimension at (i,j,:)
"""
Expand Down Expand Up @@ -123,7 +123,7 @@ def tile_domain_array(raster):


def clear_borders(tiled_array):
"""Set to zero all the edge elements of a vertical stack
"""Set to zero all the edge elements of a vertical stack
of 2D arrays. Helper function for make_weight.
**Inputs** :
Expand Down Expand Up @@ -165,45 +165,45 @@ def make_weight(Particles):
"""
L, W = Particles.stage.shape

# calculate surface slope weights
weight_sfc = (tile_domain_array(Particles.stage) \
- big_sliding_window(Particles.stage))
weight_sfc /= tile_local_array(Particles.distances, L, W)
weight_sfc[weight_sfc <= 0] = 0
clear_borders(weight_sfc)

# calculate inertial component weights
weight_int = (tile_domain_array(Particles.qx)*tile_local_array(Particles.jvec, L, W)) \
+ (tile_domain_array(Particles.qy)*tile_local_array(Particles.ivec, L, W))
weight_int /= tile_local_array(Particles.distances, L, W)
weight_int[weight_int <= 0] = 0
clear_borders(weight_int)

# get depth and cell types for neighboring cells
depth_ind = big_sliding_window(Particles.depth)
ct_ind = big_sliding_window(Particles.cell_type)

# set weights for cells that are too shallow, or invalid 0
weight_sfc[(depth_ind <= Particles.dry_depth) | (ct_ind == 2)] = 0
weight_int[(depth_ind <= Particles.dry_depth) | (ct_ind == 2)] = 0

# if sum of weights is above 0 normalize by sum of weights
norm_sfc = np.nansum(weight_sfc, axis=2)
idx_sfc = tile_domain_array((norm_sfc > 0))
weight_sfc[idx_sfc] /= tile_domain_array(norm_sfc)[idx_sfc]

norm_int = np.nansum(weight_int, axis=2)
idx_int = tile_domain_array((norm_int > 0))
weight_int[idx_int] /= tile_domain_array(norm_int)[idx_int]

# define actual weight by using gamma, and weight components
weight = Particles.gamma * weight_sfc + \
(1 - Particles.gamma) * weight_int

# modify the weight by the depth and theta weighting parameter
weight = depth_ind ** Particles.theta * weight

# if the depth is below the minimum depth then set weight to 0
weight[(depth_ind <= Particles.dry_depth) | (ct_ind == 2)] = 0

Expand Down Expand Up @@ -506,7 +506,7 @@ def particle_stepper(Particles, current_inds, travel_times):
Particles.dx),
new_cells))
# put new indices into array
new_inds = np.array(new_inds, dtype=np.int)
new_inds = np.array(new_inds, dtype=np.int32)
new_inds[np.array(dist) == 0] = 0
# see if the indices are at boundaries
new_inds = check_for_boundary(new_inds, inds, Particles.cell_type)
Expand Down
28 changes: 20 additions & 8 deletions dorado/particle_track.py
Original file line number Diff line number Diff line change
Expand Up @@ -559,8 +559,7 @@ def generate_particles(self,
# run an iteration where particles are moved
# have option of specifying the particle start locations
# otherwise they are randomly placed within x and y seed locations
def run_iteration(self,
target_time=None):
def run_iteration(self, target_time=None, max_iter=1e4):
"""Run an iteration of the particle routing.
Runs an iteration of the particle routing.
Expand All @@ -573,7 +572,13 @@ def run_iteration(self,
end of this iteration. If left undefined, then just one
iteration is run and the particles will be out of sync in time.
Note that this loop will terminate before the target_time if
the particle exceeds the hard-coded limit of 1e4 steps
the particle exceeds the maximum number of permitted
iterations (max_iter)
max_iter : `int`, optional
The maximum number of iterations a particle is allowed to
take in order to try to reach a specified target_time.
This is an optional parameter with a default value of 10,000.
**Outputs** :
Expand Down Expand Up @@ -677,7 +682,7 @@ def run_iteration(self,
est_next_dt = max(0.1, all_times[ii][-1] -
all_times[ii][-2])
count += 1
if count > 1e4:
if count > max_iter:
_iter_particles.append(ii)
break

Expand Down Expand Up @@ -883,7 +888,8 @@ def ind2coord(walk_data, raster_origin, raster_size, cellsize):


def exposure_time(walk_data,
region_of_interest):
region_of_interest,
verbose=True):
"""Measure exposure time distribution of particles in a specified region.
Function to measure the exposure time distribution (ETD) of particles to
Expand All @@ -901,6 +907,11 @@ def exposure_time(walk_data,
with 1's everywhere inside the region in which we want to
measure exposure time, and 0's everywhere else.
verbose : `bool`, optional
Toggles whether or not messages print to the
console. If False, nothing is output, if True, messages are output.
Default value is True. Errors are always raised.
**Outputs** :
exposure_times : `list`
Expand Down Expand Up @@ -950,9 +961,10 @@ def exposure_time(walk_data,

# single print statement
if len(_short_list) > 0:
print(str(len(_short_list)) + ' Particles within ROI at final'
' timestep.\n' + 'Particles are: ' + str(_short_list) +
'\nRun more iterations to get full tail of ETD.')
if verbose:
print(str(len(_short_list)) + ' Particles within ROI at final'
' timestep.\n' + 'Particles are: ' + str(_short_list) +
'\nRun more iterations to get full tail of ETD.')

return exposure_times.tolist()

Expand Down
82 changes: 64 additions & 18 deletions dorado/routines.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,8 @@
# Functions for running random walk model
# --------------------------------------------------------
def steady_plots(particle, num_iter,
folder_name=None, save_output=True):
folder_name=None, save_output=True,
verbose=True):
"""Automated particle movement in steady flow field.
Function to automate plotting of particle movement over a steady flow
Expand All @@ -45,6 +46,11 @@ def steady_plots(particle, num_iter,
Controls whether or not the output images/data are saved to disk.
Default value is True.
verbose : `bool`, optional
Toggles whether or not messages print to the
console. If False, nothing is output, if True, messages are output.
Default value is True. Errors are always raised.
**Outputs** :
walk_data : `dict`
Expand All @@ -61,7 +67,8 @@ def steady_plots(particle, num_iter,
if folder_name is None:
folder_name = os.getcwd()
if os.path.exists(folder_name):
print('Saving files in existing directory')
if verbose:
print('Saving files in existing directory')
else:
os.makedirs(folder_name)
if not os.path.exists(folder_name+os.sep+'figs'):
Expand Down Expand Up @@ -111,7 +118,7 @@ def steady_plots(particle, num_iter,

def unsteady_plots(dx, Np_tracer, seed_xloc, seed_yloc, num_steps, timestep,
output_base, output_type,
folder_name=None):
folder_name=None, verbose=True):
"""Automated particle movement in unsteady flow.
Function to automate plotting of particle movement in an unsteady flow
Expand Down Expand Up @@ -156,6 +163,11 @@ def unsteady_plots(dx, Np_tracer, seed_xloc, seed_yloc, num_steps, timestep,
folder_name : `str`, optional
Path to folder in which to save output plots.
verbose : `bool`, optional
Toggles whether or not messages print to the
console. If False, nothing is output, if True, messages are output.
Default value is True. Errors are always raised.
**Outputs** :
walk_data : `dict`
Expand All @@ -174,7 +186,8 @@ def unsteady_plots(dx, Np_tracer, seed_xloc, seed_yloc, num_steps, timestep,
if folder_name is None:
folder_name = os.getcwd()
if os.path.exists(folder_name):
print('Saving files in existing directory')
if verbose:
print('Saving files in existing directory')
else:
os.makedirs(folder_name)
if not os.path.exists(folder_name+os.sep+'figs'):
Expand All @@ -195,9 +208,10 @@ def unsteady_plots(dx, Np_tracer, seed_xloc, seed_yloc, num_steps, timestep,
qylist = sorted(qylist)
datalist = sorted(datalist)
if num_steps > max(len(depthlist), len(datalist)):
print('Warning: num_steps exceeds number of model outputs in'
' output_base')
print('Setting num_steps to equal number of model outputs')
if verbose:
print('Warning: num_steps exceeds number of model outputs in'
' output_base')
print('Setting num_steps to equal number of model outputs')
num_steps = max(len(depthlist), len(datalist))

# Create vector of target times
Expand Down Expand Up @@ -277,7 +291,7 @@ def unsteady_plots(dx, Np_tracer, seed_xloc, seed_yloc, num_steps, timestep,
return walk_data


def time_plots(particle, num_iter, folder_name=None):
def time_plots(particle, num_iter, folder_name=None, verbose=True):
"""Steady flow plots with particle travel times visualized.
Make plots with each particle's travel time visualized.
Expand All @@ -296,6 +310,11 @@ def time_plots(particle, num_iter, folder_name=None):
folder_name : `str`, optional
Path to folder in which to save output plots.
verbose : `bool`, optional
Toggles whether or not messages print to the
console. If False, nothing is output, if True, messages are output.
Default value is True. Errors are always raised.
**Outputs** :
walk_data : `dict`
Expand All @@ -308,7 +327,8 @@ def time_plots(particle, num_iter, folder_name=None):
if folder_name is None:
folder_name = os.getcwd()
if os.path.exists(folder_name):
print('Saving files in existing directory')
if verbose:
print('Saving files in existing directory')
else:
os.makedirs(folder_name)
if not os.path.exists(folder_name+os.sep+'figs'):
Expand Down Expand Up @@ -366,7 +386,7 @@ def time_plots(particle, num_iter, folder_name=None):
# --------------------------------------------------------
# Functions for plotting/interpreting outputs
# --------------------------------------------------------
def get_state(walk_data, iteration=-1):
def get_state(walk_data, iteration=-1, verbose=True):
"""Pull walk_data values from a specific iteration.
Routine to return slices of the walk_data dict at a given iteration #.
Expand All @@ -382,6 +402,11 @@ def get_state(walk_data, iteration=-1):
Iteration number at which to slice the dictionary. Defaults
to -1, i.e. the most recent step
verbose : `bool`, optional
Toggles whether or not messages print to the
console. If False, nothing is output, if True, messages are output.
Default value is True. Errors are always raised.
**Outputs** :
xinds : `list`
Expand Down Expand Up @@ -416,13 +441,14 @@ def get_state(walk_data, iteration=-1):
iter_exceeds_warning += 1

if iter_exceeds_warning > 0:
print('Note: %s particles have not reached %s iterations' % \
(iter_exceeds_warning, iteration))
if verbose:
print('Note: %s particles have not reached %s iterations' % \
(iter_exceeds_warning, iteration))

return xinds, yinds, times


def get_time_state(walk_data, target_time):
def get_time_state(walk_data, target_time, verbose=True):
"""Pull walk_data values nearest to a specific time.
Routine to return slices of the walk_data dict at a given travel time.
Expand All @@ -436,6 +462,11 @@ def get_time_state(walk_data, target_time):
target_time : `float`
Travel time at which to slice the dictionary.
verbose : `bool`, optional
Toggles whether or not messages print to the
console. If False, nothing is output, if True, messages are output.
Default value is True. Errors are always raised.
**Outputs** :
xinds : `list`
Expand Down Expand Up @@ -469,7 +500,8 @@ def get_time_state(walk_data, target_time):
times.append(walk_data['travel_times'][ii][tt])

if times_ii[-1] < target_time:
print('Note: Particle '+str(ii)+' never reached target_time')
if verbose:
print('Note: Particle '+str(ii)+' never reached target_time')

return xinds, yinds, times

Expand All @@ -479,7 +511,8 @@ def plot_exposure_time(walk_data,
folder_name=None,
timedelta=1,
nbins=100,
save_output=True):
save_output=True,
verbose=True):
"""Plot exposure time distribution of particles in a specified region.
Function to plot the exposure time distribution (ETD) of particles to
Expand Down Expand Up @@ -509,6 +542,11 @@ def plot_exposure_time(walk_data,
Controls whether or not the output images/data are saved to disk.
Default value is True.
verbose : `bool`, optional
Toggles whether or not messages print to the
console. If False, nothing is output, if True, messages are output.
Default value is True. Errors are always raised.
**Outputs** :
If `save_output` is set to True, script saves plots of the cumulative
Expand Down Expand Up @@ -537,7 +575,8 @@ def plot_exposure_time(walk_data,
if folder_name is None:
folder_name = os.getcwd()
if os.path.exists(folder_name):
print('Saving files in existing directory')
if verbose:
print('Saving files in existing directory')
else:
os.makedirs(folder_name)
if not os.path.exists(folder_name+os.sep+'figs'):
Expand Down Expand Up @@ -835,7 +874,8 @@ def snake_plots(grid,
interval=4,
tail_length=12,
rgba_start=[1, 0.4, 0.2, 1],
rgba_end=[1, 0.3, 0.1, 0]):
rgba_end=[1, 0.3, 0.1, 0],
verbose=True):
"""Plot particle positions with a trailing tail.
Loops through existing walk_data history and creates a series of
Expand Down Expand Up @@ -883,6 +923,11 @@ def snake_plots(grid,
path, specified as [red, green, blue, alpha] in the range 0-1.
Default slowly fades to red-ish with an alpha of zero
verbose : `bool`, optional
Toggles whether or not messages print to the
console. If False, nothing is output, if True, messages are output.
Default value is True. Errors are always raised.
**Outputs** :
Saves a plot of the particle travel history at each step as a png
Expand All @@ -895,7 +940,8 @@ def snake_plots(grid,
if folder_name is None:
folder_name = os.getcwd()
if os.path.exists(folder_name):
print('Saving files in existing directory')
if verbose:
print('Saving files in existing directory')
else:
os.makedirs(folder_name)
if not os.path.exists(folder_name+os.sep+'figs'):
Expand Down

0 comments on commit 179105d

Please sign in to comment.