automated DFT decimation for adjoint sources #1753

oskooi · 2021-09-05T22:40:20Z

Closes #1751.

This is an incomplete attempt to add automated DFT decimation for the adjoint sources. It follows the outline in #1751 by removing the tol parameter of get_fwidth and adding a new function set_fwidth(double fw) to the src_time class and its two variants gaussian_src_time and custom_src_time.

The piece that has yet to be added involves modifying the Python adjoint source FilteredSource (which is a wrapper around CustomSource) to call its member function set_fwidth for each frequency in its list frequencies after the CustomSource instance has been created:

meep/python/adjoint/filter_source.py

Lines 57 to 61 in e2739e4

    
           # initialize super 
        
           super(FilteredSource, self).__init__(src_func=f, 
        
                                                center_frequency=center_frequency, 
        
                                                is_integrated=False, 
        
                                                end_time=self.T)

The problem though is that the bandwidth of the Nuttall window-function implementation of python/adjoint/filter_source.py does not seem to be well defined. To demonstrate this, here is a plot of the DTFT of an actual Nuttall window (real and imaginary parts) using arbitrary but realistic values:

import numpy as np
import matplotlib.pyplot as plt


res = 40
dt = 0.5/res
fcen = 1.0
df = 0.05
nfrqs = 50
frqs = np.linspace(fcen-0.5*df,fcen+0.5*df,nfrqs)
T = np.max(np.abs(1 / np.diff(frqs)))
N = np.rint(T / dt)

def sinc(f, f0):
    num = np.where(f == f0, N + 1,
                   (1 - np.exp(1j * (N + 1) * (2 * np.pi) *
                               (f - f0) * dt)))
    den = np.where(f == f0, 1,
                   (1 - np.exp(1j * (2 * np.pi) * (f - f0) * dt)))
    return num / den

def cos_window_fd(a, f, f0):
    df = 1 / (N * dt)
    cos_sum = a[0] * sinc(f, f0)
    for k in range(1, len(a)):
        cos_sum += (-1)**k * a[k] / 2 * sinc(f, f0 - k * df) + (-1)**k * a[k] / 2 * sinc(f, f0 + k * df)
    return cos_sum

def nuttall_dtft(f, f0):
    a = [0.355768, 0.4873960, 0.144232, 0.012604]
    return cos_window_fd(a, f, f0)

if __name__ == '__main__':
    plt.figure()
    dtft = nuttall_dtft(frqs,fcen)
    plt.subplot(1,2,1)
    plt.plot(frqs,np.real(dtft),'bo-')
    plt.xlabel('frequencies')
    plt.ylabel('DTFT of nuttall window (real part)')
    plt.subplot(1,2,2)
    plt.plot(frqs,np.imag(dtft),'ro-')
    plt.xlabel('frequencies')
    plt.ylabel('DTFT of nuttall window (imaginary part)')
    plt.subplots_adjust(hspace=0.5)
    plt.tight_layout(pad=1.0)
    plt.savefig('nuttall_dtft.png',dpi=150)

(Note: the imaginary part is several orders of magnitude larger than the real part.)

One possible definition for the bandwidth of this Nuttall window is the smallest interval in which its magnitude is non-zero but this may be too conservative.

For comparison, the Nuttall window from the scipy signal package has well-defined sidelobes which obviously makes the definition of the bandwidth unambiguous.

codecov-commenter · 2021-09-05T22:55:45Z

Codecov Report

Merging #1753 (46e99b5) into master (cc626e2) will increase coverage by 0.06%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master    #1753      +/-   ##
==========================================
+ Coverage   74.43%   74.50%   +0.06%     
==========================================
  Files          13       13              
  Lines        4589     4600      +11     
==========================================
+ Hits         3416     3427      +11     
  Misses       1173     1173

Impacted Files	Coverage Δ
python/adjoint/filter_source.py	`88.46% <100.00%> (+1.69%)`	⬆️
python/source.py	`95.45% <100.00%> (+0.03%)`	⬆️

smartalecH · 2021-09-06T00:51:40Z

When defining bandwidth, we typically look at the power spectral density (PSD), not the real and imaginary parts.

oskooi · 2021-09-06T00:52:46Z

When defining bandwidth, we typically look at the power spectral density (PSD), not the real and imaginary parts.

Plotting the magnitude of the DTFT of the Nuttall window shows the side lobes more clearly as well as some of its other features including an (unexpected) discontinuity. In this specific example, it looks like the frequency response decays by roughly five orders of magnitude from its peak before reaching some kind of (slowly-decaying) floor. For comparison, the Nuttall spectrum in the scipy example decays by roughly 10 orders of magnitude from its peak value before reaching a constant floor. The spectrum also does not contain a discontinuity.

    plt.figure()
    dtft = nuttall_dtft(frqs,fcen)
    plt.semilogy(frqs,np.abs(dtft),'bo-')
    plt.xlabel('frequency')
    plt.ylabel('DTFT of Nuttall window (magnitude)')
    plt.savefig('nuttall_dtft_abs.png',dpi=150,bbox_inches='tight')

smartalecH · 2021-09-06T01:04:18Z

The scipy example is in dB (multiplied by 20, not 10) which, again, is PSD. That's why the rolloff differs by a factor of 2.

stevengj · 2021-09-15T01:23:13Z

You'll just want to add fwidth as a constructor parameter in CustomSource (in sources.py) and in custom_py_src_time (in meep-python.hpp).

Basically, handle fwidth in exactly the same way that we handle the center frequency now.

stevengj · 2021-09-15T01:32:36Z

The scipy plot looks odd because the amplitude is not decaying outside of the center lobe — probably that is because they are using fft and not the exact DTFT. The Wikipedia plot is more like what we'd expect: https://en.wikipedia.org/wiki/Window_function#/media/File:Window_function_and_frequency_response_-_Nuttall_(continuous_first_derivative).svg

Try plotting it just from 0.99 to 1.01, but sample the points 10x more finely (and plot with b- not bo-) to see if you can reproduce the oscillations of the side lobes.

The fact that scipy doesn't provide the exact DTFT means we have to implement that part on our own anyway, so we should implement our own time-domain version also to ensure that they are consistent.

stevengj · 2021-09-22T01:24:37Z

The Nutall window W(ω) (centered at ω=0) has a complicated formula, but we only care about its asymptotic behavior (for large |ω|), and that should approximately be a simple power law (proportional to #/ω^n for some exponent n and some coefficient #, both of which can be calculated analytically). This asymptotic power law should be accurate enough for computing the bandwidth.

Mathematica (available online via WolframAlpha) should be able to do this asymptotic expansion for you. Since we expect it to be a polynomial in x = 1/ω, probably you can give Mathematica W(1/x) and tell it to Taylor expand around x=0.

smartalecH · 2021-09-23T17:30:06Z

The Nutall window W(ω) (centered at ω=0) has a complicated formula, but we only care about its asymptotic behavior (for large |ω|),

Don't we already know its asymptotic behavior? We specify the basis width of each window to be the smallest Δω. Why not just use that value to determine the bandwidth?

Specifically, as long as the Nyquist rate is greater than ω_max+Δω/2, we should be good.

stevengj · 2021-09-29T01:28:00Z

You can also just get the asymptotics numerically to decent accuracy pretty easily:

For a Nutall window centered at ω = 0, plot the DFT amplitude |W(ω)| vs. ω on a log–log scale.
From the log–log scale, read off the asymptotic power law 1/ω^p (e.g 1/ω^4)
Now you just need the coefficient #/ω^p: compute # approximately by setting # = ω^p |W(ω)| for a sufficiently large ω. (If you get enough digits for #, you might even guess the formula for it from the inverse symbolic calculator.)

Given this, now you just need to work out the rescaling of the asymptotic when you change the bandwidth to Δω — I'm guessing it's something like #(Δω/ω)^p, but you should double-check.

oskooi · 2021-09-30T17:34:54Z

Following the outline above, I was able to fit the asymptotic limit of ω→∞ for the Nuttall window function to a power law with exponent p=3. Including a bandwidth Δf in the fitting function does not seem to change the results.

1. semilog plot of the DTFT window function centered at frequency f=0

2. loglog plots of the power-law fit without (left) and with (right) a bandwidth Δf

import numpy as np
import matplotlib.pyplot as plt

res = 40
dt = 0.5/res
fcen = 0
df = 4.0
nfrqs = 400
frqs = np.linspace(fcen-0.5*df,fcen+0.5*df,nfrqs)
T = np.max(np.abs(1 / np.diff(frqs)))
N = np.rint(T / dt)


def sinc(f, f0):
    num = np.where(f == f0, N + 1,
                   (1 - np.exp(1j * (N + 1) * (2 * np.pi) *
                               (f - f0) * dt)))
    den = np.where(f == f0, 1,
                   (1 - np.exp(1j * (2 * np.pi) * (f - f0) * dt)))
    return num / den

def cos_window_fd(a, f, f0):
    df = 1 / (N * dt)
    cos_sum = a[0] * sinc(f, f0)
    for k in range(1, len(a)):
        cos_sum += (-1)**k * a[k] / 2 * sinc(f, f0 - k * df) + (-1)**k * a[k] / 2 * sinc(f, f0 + k * df)
    return cos_sum

def nuttall_dtft(f, f0):
    a = [0.355768, 0.4873960, 0.144232, 0.012604]
    return cos_window_fd(a, f, f0)


if __name__ == '__main__':
    plt.figure()
    plt.subplot(1,2,1)
    dtft = nuttall_dtft(frqs,fcen)
    coeff = ((df/2)**3)*np.abs(dtft[-1])
    print("coeff:, {}".format(coeff))
    plt.loglog(frqs,np.abs(dtft),'b-',label='|W(f)|')
    plt.loglog(frqs,coeff/np.power(frqs,3),'k--',label='{:.6e}/f$^3$'.format(coeff))
    plt.xlabel('frequency, f')
    plt.ylabel('DTFT of Nuttall window, |W(f)|')
    plt.title('fit to C/f$^p$')
    plt.legend(loc='lower left')

    fwidth = 0.1
    coeff = (((df/2)/fwidth)**3)*np.abs(dtft[-1])
    print("coeff:, {}".format(coeff))
    plt.subplot(1,2,2)
    plt.loglog(frqs,np.abs(dtft),'b-',label='|W(f)|')
    plt.loglog(frqs,coeff*np.power(fwidth/frqs,3),'k--',label='{:.6e}(Δf/f)$^3$'.format(coeff))
    plt.xlabel('frequency, f')
    plt.ylabel('DTFT of Nuttall window, |W(f)|')
    plt.title('fit to C*(Δf/f)$^p$, Δf={:.1f}'.format(fwidth))
    plt.legend(loc='lower left')
    plt.savefig('nuttall_dtft_abs_fwidth.png',dpi=150,bbox_inches='tight')

stevengj · 2021-09-30T20:11:15Z

You need to make sure that this formula correctly predicts the behavior when you change the bandwidth Δf = 1/(NΔt).

In your tests above, you have N=7980 and dt=0.0125.

So, if we want 6.112498e-5 / f^3 = C * (Δf / f)^3, then we should use C = 6.112498e-5 * (7980 * 0.0125)**3 = 60.667688.

So hopefully 60.667688 * (Δf / f)^3 works, but you should check this in cases there is some additional scale factor missing if you change N or dt.

oskooi · 2021-10-02T16:22:37Z

This seems to be finally working: the frequency bandwidth of the Nuttall window used in the adjoint sources is computed analytically by fitting to an asymptotic power law. Similar to the GaussianSource, the frequency bandwidth of the Nuttall window is defined as a property of the CustomSource object (the base class of the FilteredSource of the adjoint solver) which is used to determine the optimal decimation factor of the DFT fields monitors.

I verified that the fitting procedure for C * (Δf / f)^3 works for different values of N and Δt. As an example, the figures below show the fit to N=6225 and Δt=0.01 without (left) and with (right) the bandwidth Δf=0.0161.

import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt


res = 50
dt = 0.5/res
fcen = 0
df = 4.0
nfrqs = 250
frqs = np.linspace(fcen-0.5*df,fcen+0.5*df,nfrqs)
T = np.max(np.abs(1 / np.diff(frqs)))
N = np.rint(T / dt)
fwidth = 1/(N*dt)

def sinc(f, f0):
    num = np.where(f == f0, N+1,
                   (1 - np.exp(1j * (N + 1) * (2 * np.pi) *
                               (f - f0) * dt)))
    den = np.where(f == f0, 1,
                   (1 - np.exp(1j * (2 * np.pi) * (f - f0) * dt)))
    return num / den

def cos_window_fd(a, f, f0):
    df = 1 / (N * dt)
    cos_sum = a[0] * sinc(f, f0)
    for k in range(1, len(a)):
        cos_sum += (-1)**k * a[k] / 2 * sinc(f, f0 - k * df) + (-1)**k * a[k] / 2 * sinc(f, f0 + k * df)
    return cos_sum

def nuttall_dtft(f, f0):
    a = [0.355768, 0.4873960, 0.144232, 0.012604]
    return cos_window_fd(a, f, f0)

if __name__ == '__main__':
    plt.figure()
    plt.subplot(1,2,1)
    dtft = nuttall_dtft(frqs,fcen)
    coeff1 = ((df/2)**3)*np.abs(dtft[-1])
    print("coeff1:, {}".format(coeff1))
    plt.loglog(frqs,np.abs(dtft),'b-',label='|W(f)|')
    plt.loglog(frqs,coeff1/np.power(frqs,3),'k--',label='{:.6e}/f$^3$'.format(coeff1))
    plt.xlabel('frequency, f')
    plt.ylabel('DTFT of Nuttall window, |W(f)|')
    plt.title('fit to C/f$^p$')
    plt.legend(loc='lower left')

    coeff2 = coeff1 * 1/fwidth**3
    print("coeff2:, {}".format(coeff2))
    plt.subplot(1,2,2)
    plt.loglog(frqs,np.abs(dtft),'b-',label='|W(f)|')
    plt.loglog(frqs,coeff2*np.power(fwidth/frqs,3),'k--',label='{:.6e}(Δf/f)$^3$'.format(coeff2))
    plt.xlabel('frequency, f')
    plt.ylabel('DTFT of Nuttall window, |W(f)|')
    plt.title('fit to C*(Δf/f)$^p$, Δf={:.4f}'.format(fwidth))
    plt.legend(loc='lower left')
    plt.savefig('nuttall_dtft.png',dpi=150,bbox_inches='tight')

stevengj · 2021-10-02T21:39:51Z

Note that coefficient C in C * (Δf/f)^3 should not change when you change Δf (i.e. N or dt); otherwise, you have the wrong scaling.

python/adjoint/filter_source.py

stevengj · 2021-10-07T19:17:15Z

You should also check that this gives the same answer (to many digits) as setting the decimation factor to 1.

…wer law

oskooi · 2021-10-08T05:50:12Z

After switching the fit to the function C / f^3 from C * (Δf / f)^3 and using a larger f to obtain a better approximation of the asymptotic limit, there was no longer any need to adjust the tolerances for the unit test test_adjoint_jax.py. Only a single subtest in test_adjoint_solver.py involving the DFT fields and just for the single-precision case required a less than 1% change in its tolerance but this specific test is actually flaky as previously reported in #1593 (comment). In fact, running the tests on my local i7-7700k machine using this branch with the commit 66d55e2 (which does not include any changes to the tolerances), everything passes with single and double precision so clearly there is something about this one test that is strange and can be dealt with as a separate issue.

I think therefore that this PR is ready to be merged.

oskooi · 2021-10-08T05:52:53Z

python/tests/test_adjoint_solver.py

@@ -213,7 +208,7 @@ def test_adjoint_solver_DFT_fields(self):
            adj_scale = (dp[None,:]@adjsol_grad).flatten()
            fd_grad = S12_perturbed-S12_unperturbed
            print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad))
-            tol = 0.04 if mp.is_single_precision() else 0.005
+            tol = 0.0461 if mp.is_single_precision() else 0.005


This is the one place where the tolerance needed to be slightly adjusted.

oskooi · 2021-10-13T01:42:38Z

Results from running test_adjoint_cyl.py with and without automatic DFT decimation using this branch:

Objective function

automatic DFT decimation

adjsol_obj:, 0.5175003706296474

no DFT decimation

adjsol_obj:, 0.5175003706296474

Directional derivative

automatic DFT decimation

adj_scale:, [-8.58436192e-05]

no DFT decimation

adj_scale:, [-8.58433114e-05]

edit: these results were generated with Meep compiled using single-precision floating point.

* remove tol parameter from get_fwidth and add new set_fwidth function * add fwidth parameter to custom_src_time class * compute bandwidth of Nuttall DTFT window function using asymptotic power law * fixes and docs * add set_fwidth and get_fwidth functions to custom_py_src_time class * adjust tolerances of failing unit tests * remove fwidth from formula for fitting coefficient * revert changes to tolerances in unit tests * slightly increase tolerance for single-precision case of DFT fields test

oskooi added the enhancement label Sep 5, 2021

oskooi force-pushed the adjoint_decimation branch from 720f541 to 6fc1b00 Compare October 1, 2021 17:26

oskooi changed the title ~~WIP: automated DFT decimation for adjoint sources~~ automated DFT decimation for adjoint sources Oct 2, 2021

oskooi commented Oct 3, 2021

View reviewed changes

python/adjoint/filter_source.py Outdated Show resolved Hide resolved

oskooi added 8 commits October 7, 2021 17:45

remove tol parameter from get_fwidth and add new set_fwidth function

52ed6fc

add fwidth parameter to custom_src_time class

d58e4ae

compute bandwidth of Nuttall DTFT window function using asymptotic po…

24a538c

…wer law

fixes and docs

63600df

add set_fwidth and get_fwidth functions to custom_py_src_time class

80d7c8b

adjust tolerances of failing unit tests

794bb33

remove fwidth from formula for fitting coefficient

804ef13

revert changes to tolerances in unit tests

66d55e2

oskooi force-pushed the adjoint_decimation branch from c1498ae to 66d55e2 Compare October 8, 2021 00:48

slightly increase tolerance for single-precision case of DFT fields test

46e99b5

oskooi commented Oct 8, 2021

View reviewed changes

stevengj merged commit 8f0deee into NanoComp:master Oct 13, 2021

oskooi deleted the adjoint_decimation branch October 13, 2021 01:44

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

automated DFT decimation for adjoint sources #1753

automated DFT decimation for adjoint sources #1753

oskooi commented Sep 5, 2021

codecov-commenter commented Sep 5, 2021 •

edited

Loading

smartalecH commented Sep 6, 2021

oskooi commented Sep 6, 2021

smartalecH commented Sep 6, 2021

stevengj commented Sep 15, 2021 •

edited

Loading

stevengj commented Sep 15, 2021 •

edited

Loading

stevengj commented Sep 22, 2021 •

edited

Loading

smartalecH commented Sep 23, 2021

stevengj commented Sep 29, 2021 •

edited

Loading

oskooi commented Sep 30, 2021

stevengj commented Sep 30, 2021

oskooi commented Oct 2, 2021

stevengj commented Oct 2, 2021 •

edited

Loading

stevengj commented Oct 7, 2021 •

edited

Loading

oskooi commented Oct 8, 2021

oskooi Oct 8, 2021

oskooi commented Oct 13, 2021 •

edited

Loading

	# initialize super
	super(FilteredSource, self).__init__(src_func=f,
	center_frequency=center_frequency,
	is_integrated=False,
	end_time=self.T)

automated DFT decimation for adjoint sources #1753

automated DFT decimation for adjoint sources #1753

Conversation

oskooi commented Sep 5, 2021

codecov-commenter commented Sep 5, 2021 • edited Loading

Codecov Report

smartalecH commented Sep 6, 2021

oskooi commented Sep 6, 2021

smartalecH commented Sep 6, 2021

stevengj commented Sep 15, 2021 • edited Loading

stevengj commented Sep 15, 2021 • edited Loading

stevengj commented Sep 22, 2021 • edited Loading

smartalecH commented Sep 23, 2021

stevengj commented Sep 29, 2021 • edited Loading

oskooi commented Sep 30, 2021

stevengj commented Sep 30, 2021

oskooi commented Oct 2, 2021

stevengj commented Oct 2, 2021 • edited Loading

stevengj commented Oct 7, 2021 • edited Loading

oskooi commented Oct 8, 2021

oskooi Oct 8, 2021

Choose a reason for hiding this comment

oskooi commented Oct 13, 2021 • edited Loading

codecov-commenter commented Sep 5, 2021 •

edited

Loading

stevengj commented Sep 15, 2021 •

edited

Loading

stevengj commented Sep 15, 2021 •

edited

Loading

stevengj commented Sep 22, 2021 •

edited

Loading

stevengj commented Sep 29, 2021 •

edited

Loading

stevengj commented Oct 2, 2021 •

edited

Loading

stevengj commented Oct 7, 2021 •

edited

Loading

oskooi commented Oct 13, 2021 •

edited

Loading