Skip to content

Commit

Permalink
Add verbosity flag for MPB (#1302)
Browse files Browse the repository at this point in the history
* Return the prior verbosity level setting. Change the parameter name to `level`

* Add verbosity level global variable, and a Python function for setting it

* Switch the verbosity value and related code to a class

* switch parameter name to `level`, prep for reuse with other cvars

* Move verbosity module to meep

* Rename module to avoid name conflicts with names of instances

* Restore returning the former level from set()

* Use the Verbosity class in meep

* Add ability to set the initial level

* Make sure we use mpb's verbosity flag, and not colliding with meep's

* Move the creation of the verbosity obj so it's available when submodules are imported

* Add explicit rules for the meep.mpb .py files so they will be copied into place with a simple `make`

* Add verbosity checks in solver.py

* Add verbosity checks in libpympb/pympb.cpp

* Revert a Makefile change

* Handle the case where a cvar object does not yet have a verbosity attribute.
This will help with transitioning while the feature rolls out.

* Python2 doesn't have f-strings

* empty commit to trigger build

* Support building with an older MPB which doesn't have the mpb_verbosity global

* Set the verbosity default to 1

* match MPB

Co-authored-by: Steven G. Johnson <[email protected]>
  • Loading branch information
RobinD42 and stevengj authored Jul 31, 2020
1 parent ebaee4f commit e416477
Show file tree
Hide file tree
Showing 9 changed files with 332 additions and 150 deletions.
165 changes: 105 additions & 60 deletions libpympb/pympb.cpp

Large diffs are not rendered by default.

10 changes: 10 additions & 0 deletions libpympb/pympb.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,16 @@

namespace py_mpb {

// A global "verbosity" level.
// 0: minimal output
// 1: a little (default)
// 2: a lot
// 3: debugging
//
// This is redeclared here so SWIG will see it.
extern "C" int mpb_verbosity;


#define TWOPI 6.2831853071795864769252867665590057683943388

void map_data(mpb_real *d_in_re, int size_in_re, mpb_real *d_in_im, int size_in_im, int n_in[3],
Expand Down
4 changes: 3 additions & 1 deletion python/Makefile.am
Original file line number Diff line number Diff line change
Expand Up @@ -197,7 +197,9 @@ HL_IFACE = \
$(srcdir)/simulation.py \
$(srcdir)/source.py \
$(srcdir)/visualization.py \
$(srcdir)/materials.py
$(srcdir)/materials.py \
$(srcdir)/verbosity_mgr.py


pkgpython_PYTHON = __init__.py $(HL_IFACE)

Expand Down
56 changes: 28 additions & 28 deletions python/adjoint/optimization_problem.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ def get_gradient(self,fields_a,fields_f,frequencies,geom_list,f):
grad = np.zeros((self.num_design_params*num_freqs,)) # preallocate

# compute the gradient
mp._get_gradient(grad,fields_a,fields_f,self.volume,np.array(frequencies),geom_list,f)
mp._get_gradient(grad,fields_a,fields_f,self.volume,np.array(frequencies),geom_list,f)

return np.squeeze(grad.reshape(self.num_design_params,num_freqs,order='F'))

Expand All @@ -42,7 +42,7 @@ class OptimizationProblem(object):
"""

def __init__(self,
def __init__(self,
simulation,
objective_functions,
objective_arguments,
Expand Down Expand Up @@ -70,7 +70,7 @@ def __init__(self,
self.design_regions = design_regions
else:
self.design_regions = [design_regions]

self.num_design_params = [ni.num_design_params for ni in self.design_regions]
self.num_design_regions = len(self.design_regions)

Expand All @@ -83,7 +83,7 @@ def __init__(self,
self.nf = nf
self.frequencies = [fcen]
else:
fmax = fcen+0.5*df
fmax = fcen+0.5*df
fmin = fcen-0.5*df
dfreq = (fmax-fmin)/(nf-1)
self.frequencies = np.linspace(fmin, fmin+dfreq*nf, num=nf, endpoint=False)
Expand All @@ -98,7 +98,7 @@ def __init__(self,
self.decay_fields=decay_fields
self.decay_dt=decay_dt
self.minimum_run_time=minimum_run_time

# store sources for finite difference estimations
self.forward_sources = self.sim.sources

Expand All @@ -115,7 +115,7 @@ def __call__(self, rho_vector=None, need_value=True, need_gradient=True):
self.update_design(rho_vector=rho_vector)

# Run forward run if requested
if need_value and self.current_state == "INIT":
if need_value and self.current_state == "INIT":
print("Starting forward run...")
self.forward_run()

Expand Down Expand Up @@ -162,7 +162,7 @@ def _df(x=None):
return df

return _f, _df

def prepare_forward_run(self):
# prepare forward run
self.sim.reset_meep()
Expand Down Expand Up @@ -226,7 +226,7 @@ def adjoint_run(self,objective_idx=0):
self.sim.change_sources(self.adjoint_sources)

# register design flux
# TODO use yee grid directly
# TODO use yee grid directly
self.design_region_monitors = [self.sim.add_dft_fields([mp.Ex,mp.Ey,mp.Ez],self.frequencies,where=dr.volume,yee_grid=False) for dr in self.design_regions]

# Adjoint run
Expand All @@ -245,7 +245,7 @@ def adjoint_run(self,objective_idx=0):
def calculate_gradient(self):
# Iterate through all design regions and calculate gradient
self.gradient = [[dr.get_gradient(self.a_E[ar][dri],self.d_E[dri],self.frequencies,self.sim.geometry,self.sim.fields) for dri, dr in enumerate(self.design_regions)] for ar in range(len(self.objective_functions))]

# Cleanup list of lists
if len(self.gradient) == 1:
self.gradient = self.gradient[0] # only one objective function
Expand All @@ -256,7 +256,7 @@ def calculate_gradient(self):
self.gradient = [g[0] for g in self.gradient] # multiple objective functions bu one design region
# Return optimizer's state to initialization
self.current_state = "INIT"

def calculate_fd_gradient(self,num_gradients=1,db=1e-4,design_variables_idx=0,filter=None):
'''
Estimate central difference gradients.
Expand All @@ -272,7 +272,7 @@ def calculate_fd_gradient(self,num_gradients=1,db=1e-4,design_variables_idx=0,fi
Returns
-----------
fd_gradient ... : lists
fd_gradient ... : lists
[number of objective functions][number of gradients]
'''
Expand All @@ -292,25 +292,25 @@ def calculate_fd_gradient(self,num_gradients=1,db=1e-4,design_variables_idx=0,fi
fd_gradient_idx = np.random.choice(self.num_design_params[design_variables_idx],num_gradients,replace=False)

for k in fd_gradient_idx:

b0 = np.ones((self.num_design_params[design_variables_idx],))
b0[:] = (self.design_regions[design_variables_idx].design_parameters.design_parameters)
# -------------------------------------------- #
# left function evaluation
# -------------------------------------------- #
self.sim.reset_meep()

# assign new design vector
b0[k] -= db
self.design_regions[design_variables_idx].update_design_parameters(b0)

# initialize design monitors
self.forward_monitors = []
for m in self.objective_arguments:
self.forward_monitors.append(m.register_monitors(self.frequencies))

self.sim.run(until_after_sources=stop_when_dft_decayed(self.sim, self.forward_monitors, self.decay_dt, self.decay_fields, self.fcen_idx, self.decay_by, self.minimum_run_time))

# record final objective function value
results_list = []
for m in self.objective_arguments:
Expand All @@ -330,10 +330,10 @@ def calculate_fd_gradient(self,num_gradients=1,db=1e-4,design_variables_idx=0,fi
self.forward_monitors = []
for m in self.objective_arguments:
self.forward_monitors.append(m.register_monitors(self.frequencies))

# add monitor used to track dft convergence
self.sim.run(until_after_sources=stop_when_dft_decayed(self.sim, self.forward_monitors, self.decay_dt, self.decay_fields, self.fcen_idx, self.decay_by, self.minimum_run_time))

# record final objective function value
results_list = []
for m in self.objective_arguments:
Expand All @@ -344,13 +344,13 @@ def calculate_fd_gradient(self,num_gradients=1,db=1e-4,design_variables_idx=0,fi
# estimate derivative
# -------------------------------------------- #
fd_gradient.append( [np.squeeze((fp[fi] - fm[fi]) / (2*db)) for fi in range(len(self.objective_functions))] )

# Cleanup singleton dimensions
if len(fd_gradient) == 1:
fd_gradient = fd_gradient[0]

return fd_gradient, fd_gradient_idx

def update_design(self, rho_vector):
"""Update the design permittivity function.
Expand All @@ -360,15 +360,15 @@ def update_design(self, rho_vector):
if np.array(rho_vector[bi]).ndim > 1:
raise ValueError("Each vector of design variables must contain only one dimension.")
b.update_design_parameters(rho_vector[bi])

self.sim.reset_meep()
self.current_state = "INIT"
def get_objective_arguments(self):
'''Return list of evaluated objective arguments.
'''
objective_args_evaluation = [m.get_evaluation() for m in self.objective_arguments]
return objective_args_evaluation

def plot2D(self,init_opt=False, **kwargs):
"""Produce a graphical visualization of the geometry and/or fields,
as appropriately autodetermined based on the current state of
Expand Down Expand Up @@ -397,15 +397,15 @@ def stop_when_dft_decayed(simob, mon, dt, c, fcen_idx, decay_by, minimum_run_tim
x.append(len(xi))
y.append(len(yi))
z.append(len(zi))

# Record data in closure so that we can persitently edit
closure = {
'previous_fields': [np.ones((x[mi],y[mi],z[mi],len(c)),dtype=np.complex128) for mi, m in enumerate(mon)],
't0': 0
}

def _stop(sim):

if sim.round_time() <= dt + closure['t0']:
return False
else:
Expand All @@ -427,11 +427,11 @@ def _stop(sim):
relative_change = np.mean(relative_change) # average across monitors
closure['previous_fields'] = current_fields
closure['t0'] = sim.round_time()
if mp.cvar.verbosity > 0:

if mp.verbosity > 0:
fmt = "DFT decay(t = {0:1.1f}): {1:0.4e}"
print(fmt.format(sim.meep_time(), np.real(relative_change)))
return relative_change <= decay_by and sim.round_time() >= minimum_run_time
return relative_change <= decay_by and sim.round_time() >= minimum_run_time
return _stop


Expand All @@ -440,7 +440,7 @@ def atleast_3d(*arys):
'''
Modified version of numpy's `atleast_3d`
Keeps one dimensional array data in first dimension, as
Keeps one dimensional array data in first dimension, as
opposed to moving it to the second dimension as numpy's
version does. Keeps the meep dimensionality convention.
Expand Down
6 changes: 6 additions & 0 deletions python/mpb.i
Original file line number Diff line number Diff line change
Expand Up @@ -346,9 +346,14 @@ static mpb_real field_integral_energy_callback(mpb_real energy, mpb_real epsilon

%apply double { number };

%rename(verbosity) mpb_verbosity;
%include "pympb.hpp"


%pythoncode %{
from meep.verbosity_mgr import Verbosity
verbosity = Verbosity(_mpb.cvar, 1)

from .solver import (
MPBArray,
ModeSolver,
Expand Down Expand Up @@ -390,4 +395,5 @@ static mpb_real field_integral_energy_callback(mpb_real energy, mpb_real epsilon
from .mpb_data import (
MPBData,
)

%}
27 changes: 11 additions & 16 deletions python/simulation.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@
from meep.geom import Vector3, init_do_averaging
from meep.source import EigenModeSource, check_positive
import meep.visualization as vis

from meep.verbosity_mgr import Verbosity

try:
basestring
Expand All @@ -34,6 +34,8 @@
except ImportError:
do_progress = False

verbosity = Verbosity(mp.cvar, 1)


# Send output from Meep, ctlgeom, and MPB to Python's stdout
mp.set_meep_printf_callback(mp.py_master_printf_wrap)
Expand Down Expand Up @@ -1556,7 +1558,7 @@ def _compute_fragment_stats(self, gv):
return stats

def _init_structure(self, k=False):
if mp.cvar.verbosity > 0:
if verbosity > 0:
print('-' * 11)
print('Initializing structure...')

Expand All @@ -1579,7 +1581,7 @@ def _init_structure(self, k=False):
if self.collect_stats and isinstance(self.default_material, mp.Medium):
self.fragment_stats = self._compute_fragment_stats(gv)

if self._output_stats and isinstance(self.default_material, mp.Medium) and mp.cvar.verbosity > 0:
if self._output_stats and isinstance(self.default_material, mp.Medium) and verbosity > 0:
stats = self._compute_fragment_stats(gv)
print("STATS: aniso_eps: {}".format(stats.num_anisotropic_eps_pixels))
print("STATS: anis_mu: {}".format(stats.num_anisotropic_mu_pixels))
Expand Down Expand Up @@ -1851,7 +1853,7 @@ def use_real(self):

if use_real(self):
self.fields.use_real_fields()
elif mp.cvar.verbosity > 0:
elif verbosity > 0:
print("Meep: using complex fields.")

if self.k_point:
Expand Down Expand Up @@ -2039,7 +2041,7 @@ def use_output_directory(self, dname=''):
closure = {'trashed': False}

def hook():
if mp.cvar.verbosity > 0:
if verbosity > 0:
print("Meep: using output directory '{}'".format(dname))
self.fields.set_output_directory(dname)
if not closure['trashed']:
Expand Down Expand Up @@ -2100,7 +2102,7 @@ def stop_cond(sim):
self.progress.value = t0 + stop_time
self.progress.description = "100% done "

if mp.cvar.verbosity > 0:
if verbosity > 0:
print("run {} finished at t = {} ({} timesteps)".format(self.run_index, self.meep_time(), self.fields.t))
self.run_index += 1

Expand Down Expand Up @@ -4057,7 +4059,7 @@ def _stop(sim):
closure['cur_max'] = 0
closure['t0'] = sim.round_time()
closure['max_abs'] = max(closure['max_abs'], old_cur)
if closure['max_abs'] != 0 and mp.cvar.verbosity > 0:
if closure['max_abs'] != 0 and verbosity > 0:
fmt = "field decay(t = {}): {} / {} = {}"
print(fmt.format(sim.meep_time(), old_cur, closure['max_abs'], old_cur / closure['max_abs']))
return old_cur <= closure['max_abs'] * decay_by
Expand Down Expand Up @@ -4175,7 +4177,7 @@ def _disp(sim):
sim.progress.value = val1
sim.progress.description = "{}% done ".format(int(val2))

if mp.cvar.verbosity > 0:
if verbosity > 0:
print(msg_fmt.format(val1, t, val2, val3, val4))
closure['tlast'] = t1

Expand Down Expand Up @@ -4891,15 +4893,8 @@ def quiet(quietval=True):
output can be enabled again by passing `False`. This sets a global variable, so the
value will persist across runs within the same script.
"""
mp.cvar.verbosity = int(not quietval)
verbosity(int(not quietval))

def verbosity(v=1):
"""
Given a number `v`, specify the degree of Meep's output: `0` is quiet mode, `1` (the
default) is ordinary output, `2` is extra debugging output, and `3` is all debugging
output.
"""
mp.cvar.verbosity = v

def get_num_groups():
# Lazy import
Expand Down
Loading

0 comments on commit e416477

Please sign in to comment.