Unit test for constraint functions of a periodic design region

[oskooi]

#2465 added support for periodic design regions to the filters used in the adjoint solver. To verify that this feature is working correctly, a possible test that we discussed involves evaluating the constraint function for the solid and void features of the design weights and checking that the results are unchanged whether or not a feature with known lengthscale (i.e., a circle with a given diameter) crosses the edge of the design region. (The test in #2465 used a design region with random weights which may not be appropriate since it can contain multiple features with the same lengthscale.)

The test described above is set up in this gist. However, the results of this test are different for the two cases which seems to suggest that #2465 is not working correctly.

Case 1: circle does not cross the edge of the design region

For this test, we conduct two checks:

1. Given a minimum length scale larger than the diameter of the circle, the constraint function for the solid and void regions should be < ~0 (i.e., constraint is not violated).
2. Given a minimum length scale smaller than the diameter of the circle, the constraint function for the solid and void regions should be > 0 (i.e., constraint is violated).

The results from this test match the expected values:

test1-solid:, 1.2413e-09 (should be < ~0)
test1-void:,  1.2413e-09 (should be < ~0)
test2-solid:, 0.000902384 (should be > 0)
test2-void:,  0.000902384 (should be > 0)

Also, note that the value of the constraint function for the solid and void regions is identical (as expected based on the single test feature of a circle).

Case 2: circle crosses the edge of the design region

The results from this test do not match those in Case 1.

test1-solid:, 0.000205845 (should be < ~0)
test1-void:,  0.000205845 (should be < ~0)
test2-solid:, 0.000229071 (should be > 0)
test2-void:,  0.000229071 (should be > 0)

The value for all checks is identical and positive which means the constraint is violated.

[mawc2019]

The argument periodic_axes needs to specified in both constraint functions and filter functions, as demonstrated in the test. In the above example, this argument is not specified in the filter function, so the test fails as expected.
That being said, specifying periodic_axes twice is inconvenient. Before #2465, only resolution needs to be specified twice when a constraint function is used.

[mawc2019]

In addition, in Case 2, under periodic boundary conditions, the pattern is repeated quadrants instead of repeated discs. To make repeated discs, we need to add three quadrants at the other three corners.

[oskooi]

The argument periodic_axes needs to specified in both constraint functions and filter functions, as demonstrated in the test. In the above example, this argument is not specified in the filter function, so the test fails as expected.

In addition, in Case 2, under periodic boundary conditions, the pattern is repeated quadrants instead of repeated discs. To make repeated discs, we need to add three quadrants at the other three corners.

Good catch. I added the periodic_axes parameter to the function filter_f and also updated the design weights in Case 2 to ensure the circle is preserved under translations by the lattice vectors.

However, even with these changes, the Case 2 results do not match Case 1 (which are unaffected by these changes).

Case 2: circle crosses the edge of the design region

test1-solid:, 1.72534e-13 (should be < ~0)
test1-void:,  1.72534e-13 (should be < ~0)
test2-solid:, 2.60596e-10 (should be > 0)
test2-void:,  2.60596e-10 (should be > 0)

[mawc2019]

We need to make more revisions on the test code:

1. This line needs to be deleted. If not, the maximum value in the test array is 0.25 instead of 1.
2. The argument periodic_axes should be assigned as (0,1) instead of (1,1). This argument specifies the axes along which the input array is periodic, as described in the docstring. It follows the same convention as some numpy and scipy functions, e.g., numpy.fft.fft2.
3. The four quadrants in Case 2 composes a disc that is roughly the same as the complete disc in Case 1, but not exactly identical. To make them exactly identical, we can shift the complete disc as done in the test. To be specific, we can replace these lines with the following lines.

    idx_row_shift, idx_col_shift = int(Nx/2), int(Ny/2)
    d = np.vstack((d[idx_row_shift:, :], d[0:idx_row_shift, :]))
    d = np.hstack((d[:, idx_col_shift:], d[:, 0:idx_col_shift]))

[oskooi]

Applying fixes (1) and (2) seems to finally produce the expected agreement for Case 1 and 2. Fix 3 is not necessary because the results for Case 1 and Case 2 (second check in which the constraint is violated) converge to the same value with resolution: a log-log plot of the relative error in the constraint-function value (for which Case 1 is the reference result) vs. resolution shows linear convergence.

I think these results would be useful as an additional unit test which will be added shortly.

resolution = 100

Case 1

test1-solid:, 1.24131e-09 (should be < ~0)
test1-void:,  1.24131e-09 (should be < ~0)
test2-solid:, 0.000902384 (should be > 0)
test2-void:,  0.000902384 (should be > 0)

Case 2

test1-solid:, 1.72768e-09 (should be < ~0)
test1-void:,  1.72768e-09 (should be < ~0)
test2-solid:, 0.000916146 (should be > 0)
test2-void:,  0.000916146 (should be > 0)

resolution = 200

Case 1

test1-solid:, 1.29696e-09 (should be < ~0)
test1-void:,  1.29696e-09 (should be < ~0)
test2-solid:, 0.000904221 (should be > 0)
test2-void:,  0.000904221 (should be > 0)

Case 2

test1-solid:, 1.51548e-09 (should be < ~0)
test1-void:,  1.51548e-09 (should be < ~0)
test2-solid:, 0.000910708 (should be > 0)
test2-void:,  0.000910708 (should be > 0)

resolution = 400

Case 1

test1-solid:, 1.3198e-09 (should be < ~0)
test1-void:,  1.3198e-09 (should be < ~0)
test2-solid:, 0.00090495 (should be > 0)
test2-void:,  0.00090495 (should be > 0)

Case 2

test1-solid:, 1.42613e-09 (should be < ~0)
test1-void:,  1.42613e-09 (should be < ~0)
test2-solid:, 0.000908179 (should be > 0)
test2-void:,  0.000908179 (should be > 0)

[stevengj]

Maybe plot the constraint function as a function of the lengthscale constraint L — it should become > 0 for L > 1. (For L < 1, the current algorithm should give something ≈ 0.)

Should also be pretty easy to go back to old versions and see whether it changed, since it doesn't require any optimization etc.

[oskooi]

The constraint functions constraint_solid and constraint_void contain a resolution parameter that is not documented:

meep/python/adjoint/filters.py

Lines 870 to 893 in 79e8d17

	def constraint_solid(
	x, c, eta_e, filter_f, threshold_f, resolution, periodic_axes=None
	):
	"""Calculates the constraint function of the solid phase needed for minimum length optimization [1].
	Parameters
	----------
	x : array_like
	Design parameters
	c : float
	Decay rate parameter (1e0 - 1e8)
	eta_e : float
	Erosion threshold limit (0-1)
	filter_f : function_handle
	Filter function. Must be differntiable by autograd.
	threshold_f : function_handle
	Threshold function. Must be differntiable by autograd.
	periodic_axes: array_like (1D)
	List of axes (x, y = 0, 1) that are to be treated as periodic (default is none: all axes are non-periodic)
	Returns
	-------
	float
	Constraint value

This resolution argument is then passed from the constraint functions to indicator_solid and indicator_void which also do not document this parameter:

meep/python/adjoint/filters.py

Lines 918 to 934 in 79e8d17

	def indicator_void(x, c, filter_f, threshold_f, resolution, periodic_axes=None):
	"""Calculates the indicator function for the void phase needed for minimum length optimization [1].
	Parameters
	----------
	x : array_like
	Design parameters
	c : float
	Decay rate parameter (1e0 - 1e8)
	eta_d : float
	Dilation threshold limit (0-1)
	filter_f : function_handle
	Filter function. Must be differntiable by autograd.
	threshold_f : function_handle
	Threshold function. Must be differntiable by autograd.
	periodic_axes: array_like (1D)
	List of axes (x, y = 0, 1) that are to be treated as periodic (default is none: all axes are non-periodic)

@smartalecH: is resolution the resolution of the design region?

[smartalecH]

Yes.

[oskooi]

Thanks for clarifying. In my test (results to be posted shortly), mpa.constraint_solid (and mpa.constraint_void) return nearly the same result (i.e., up to the sixth decimal digit) whether resolution is 1 or the actual resolution of the design region (i.e., 100). The resolution parameter therefore does not seem to have much of an impact on the final results. Is this expected?

[smartalecH]

IIRC it should absolutely have an impact. It scales the numerical spatial gradient within the indicator function. Not including it is analogous to scaling c (The indicator function parameter).

[oskooi]

Maybe plot the constraint function as a function of the lengthscale constraint L — it should become > 0 for L > 1. (For L < 1, the current algorithm should give something ≈ 0.)

Here is a plot of the constraint function (mpa.constraint_solid) vs. the lengthscale constraint $L$ for a circle with diameter $d$. This test does not involve the period_axes feature. The results demonstrate that the constraint-function values are opposite to convention:

1. When the constraint is violated ($L &gt; d$), the constraint function returns a value of 0.
2. When the constraint is not violated ($L \le d$), the constraint functions returns a value $\geq 0$.

Note, however, that it seems that as $L \rightarrow 0$, the constraint function $\rightarrow 0$ for all four test structures.

The test used to generated these results is available here.

[smartalecH]

Ya this doesn't look good.

Can you take a few cases where the constraint should obviously be failing (and it's not) and plot

• the indicator function
• the projected field

[oskooi]

Here are plots of (1) the original design weights, (2) projected field (design_field of mpa.indicator_solid), and (3) indicator function, for a circle with diameter d = 1.0 and a lengthscale constraint of L = 0.8367 for which mpa.constraint_solid returns a value of 0.000503. This is the largest value from the data in the plot from my previous comment. In this case, L < d (no violation) and thus the constraint function should return a value of ~0.

The plot of (2) and (3) shows that the circle has a diameter which is clearly smaller than (1).

[smartalecH]

Ok there are multiple issues with the test you proposed. But here's what's killing things: you are trying to place a lengthscale on the latent field, and ignoring what the projected field looks like. In reality, you need to place a lengthscale on the projected field. Which means you need to figure out what the latent field for a particular $\beta$ and filter radius looks like such that it yields a circle with a radius of 0.5. In other words, every point on your graph will have a slightly different latent field.

Let's look at example, in particular, the extreme case when l_c=2.0. This corresponds to a filter radius of 2.0 um as well, which actually smears your circle quite a bit:

The maximum value is 0.15, well below $\eta_i=0.5$! If we project, we may think things are fine... until we add a colorbar:

Wow, $10^{-10}$ is the largest value! In other words, your constraint function thinks it's looking at empty space:

The indicator function is practically 0 everywhere (which means the constraint function should be 0).

So. If we instead determine what the latent field is supposed to look like (doing a mini optimization problem) we can get much better results:

Where we can start to see the constraint increasing around $L=1$, which is the diameter of our projected (not the latent) circle.

How does the c0 parameter affect all of this? Let's now increase our c0 by a factor of 100 (from 800 to 80,000). We now get the following:

We see that the exponential is "steeper". It looks like the curve doesn't start increasing until 1.25 now... but in reality it's hard to tell in linear scale. So in log scale:

We see that the threshold point (at $L=1$) is actually the point that first dips below the numerical noise floor on my system (about 1e-15).

Ideally we would use plots like this to determine what our threshold should be.* The plot is good for all optimizations for the given resolution, minimum length scale, and c0 parameter.

Here's some other things I noticed:

• Don't bother flattening your design parameters... we have a spatial gradient function in the indicator function that depends on array dimensionality. While you account for this with an explicit reshape, it's just asking for an error.
• Don't pass in the complement of the design parameters into your void constraint. The complementary math is already factored into the function (and its indicator function) where needed.
• Your c0 prefac seems waaay too large to me.

[oskooi]

Thanks for the detailed analysis. However, even with these modifications, the constraint function value does not produce the expected behavior. In your plot of the constraint function value vs. $L$, there is no discernible change at $L = 1.0$ where the constraint function should transition from no violation ($L \leq 1$) to violation ($L &gt; 1$). The user would somehow need to define, in advance, a threshold value for the constraint function to be able to identify the minimum $L$ for which a violation occurs as you point out:

Ideally we would use plots like this to determine what our threshold should be.

However, this does not seem to be practical for typical applications where essentially nothing is known about the design weights a priori.

What we would like is for the constraint function to change discontinuously at the minimum lengthscale. This does not seem to be possible with our existing approach, unfortunately.

Your c0 prefac seems waaay too large to me.

These values for the hyperparameters are from the tutorial Broadband Waveguide Mode Converter with Minimum Feature Size:

meep/python/examples/adjoint_optimization/mode_converter.py

Lines 273 to 274 in b3f42f2

	# hyper parameter (constant factor and exponent)
	c0 = 1e7 * (filter_radius * 1 / resolution) ** 4

The values were obtained through some experimentation and introduced in NanoComp/photonics-opt-testbed#32 (comment) because they seemed to produce decent results.

[smartalecH]

there is no discernible change at $L=1$

Not true. Please look again at the log plot, posted again below for convenience.

The user would somehow need to define, in advance, a threshold value for the constraint function to be able to identify the minimum

Yes, this is always the case for constraint functions susceptible to numerical noise (and numerical solvers in general!)

However, this does not seem to be practical for typical applications where essentially nothing is known about the design weights a priori.

... but it doesn't depend on the weights themselves. As I said earlier, it depends on the minimum length scale and resolution, and c0 parameter.

What we would like is for the constraint function to change discontinuously at the minimum lengthscale.

Wrong for two reasons. For one, we do have that. And, that's not ideal. Ideally, we'd have something with a continuous derivative (and a continuous second derivative etc) so that the optimizer can more easily traverse (or in the case of ccsa, approximate it with an expansion).

These values for the hyperparameters are from the tutorial

Again, they seem too high (even if they heuristically worked for one problem). I would go back to the original paper. Things break down quickly when c0 gets too high. The best way to check, is to look at the indicator function.

[smartalecH]

@oskooi just to add a bit more signal to the discussion, I decided to run a different test. The current test keeps the (projected) geometry fixed, and sweeps the lengthscale constraint. This isn't what the optimizer does, of course. In reality, we want to keep the lengthscale constraint fixed, and sweep the geometry size (the diameter in this case).

So I modified the above, fixing the linewidth constraint to $L=1.0$. I ran this test for various values of $c$:

Again we see that discontinuity that you were looking for (although, remember, that's not necessarily a good thing... it depends on the optimizer and where it's at). We see different thresholds for each $c$. If we changed the resolution, we would similarly see different thresholds (although it'd be good to check).

But again, these thresholds shouldn't depend on the densities.

[stevengj]

As commented above and in #2517, the lengthscale constraints need to (a) be applied to filtered/projected densities and (b) must be calibrated for a particular filter radius.

To calibrate it, you could calibrate it on a disc. Then, as a test case we could use the same constraint on a different shape with a known lengthscale (e.g. stripes or annuli).