Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Pade approximant based spectral extrapolation #2440

Merged
merged 25 commits into from
Mar 31, 2023

Conversation

Arjun-Khurana
Copy link
Contributor

After too much delay, this is an initial PR for Pade-approximant based spectral extrapolation, hopefully closing #2192. It essentially works in the same way as Harminv (thanks @Alex-Kaylor for help with the initial implementation) by aggregating field and computing an analytic frequency response using the rational Pade Approximant.

Like Harminv, Pade is constructed as a class that holds as a member variable (data) the aggregated field data, parameters for sampling the data, and the computed analytic function.

To initialize, one need only provide a field component c and a point pt at which to perform the computation.

Optional arguments include sample_rate, start_time and stop_time, which sample the collected field data according to samples = self.data[self.start_time : self.stop_time : self.sample_rate].

Additionally, there are m and n, which respectively set the orders of the numerator and denominator. Alternatively, one can specify m_frac or n_frac, which specify m and n as fractions of the length of the collected time series. If none of these are specified, the default behavior is to use half the length of the time series as the order of the numerator m, and N-m-2 as the order of the denominator.

Speaking of the algorithm, the code is currently implemented by calling the builtin scipy pade routine. Separately, I have investigated the behavior of the Trefethen SVD-based algorithm and found it to be much harder to get it to converge (though in theory it is more robust to noise in the input data). Thus, the existing implementation only uses the builtin scipy algorithm, but it would be nice to at least have the option of using either algorithm. The way I have implemented the Trefethen algorithm here provides straightforward potential interchange of Trefethen and scipy.

I have also written a test case that compares the center frequency of a peak in the Pade frequency response to the first Harminv mode, included in test_ring.py. This case passes when I run it locally with python3 -m pytest python/tests/test_ring.py but fails on my fork's CI workflow, I'm sure I'm doing something incorrectly there.

@Arjun-Khurana Arjun-Khurana marked this pull request as draft March 14, 2023 14:02
@Arjun-Khurana
Copy link
Contributor Author

Things to do before I would consider this PR "complete":

  • Understand and correctly implement the Trefethen algorithm as an option
  • Create a stop_when_pade_converged style step function that computes the Pade frequency response every some number of time steps and stops when it converges within some tolerance

Suggestions and improvements are welcome.

@smartalecH
Copy link
Collaborator

I would punt on the SVD approach. You can always add more PRs with additional functionality later.

For now, focus on the simplest working routine that's easy to augment later.

@smartalecH
Copy link
Collaborator

smartalecH commented Mar 14, 2023

This is why the CI dies (you can download the log as an artifact):

ERROR: test_pade (__main__.TestRing)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/runner/work/meep/meep/build/meep-1.27.0-beta/_build/sub/python/../../../python/tests/test_ring.py", line 77, in test_pade
    self.init()
  File "/home/runner/work/meep/meep/build/meep-1.27.0-beta/_build/sub/python/../../../python/tests/test_ring.py", line 49, in init
    self.p = mp.Pade(mp.Ez, mp.Vector3(r + 0.1), sample_rate=4)
AttributeError: module 'meep' has no attribute 'Pade'

Be sure to include your new module/function in meep.i (The swig file):

from .simulation import (

python/simulation.py Outdated Show resolved Hide resolved
Copy link
Collaborator

@oskooi oskooi left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good work. We should add a tutorial for this feature after this is merged.

python/tests/test_ring.py Outdated Show resolved Hide resolved
python/tests/test_ring.py Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
python/simulation.py Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
@codecov-commenter
Copy link

codecov-commenter commented Mar 15, 2023

⚠️ Please install the 'codecov app svg image' to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

Attention: Patch coverage is 98.03922% with 1 line in your changes missing coverage. Please review.

Project coverage is 72.90%. Comparing base (d5d172c) to head (32284d7).
Report is 140 commits behind head on master.

Files with missing lines Patch % Lines
python/simulation.py 98.03% 1 Missing ⚠️

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files
@@            Coverage Diff             @@
##           master    #2440      +/-   ##
==========================================
+ Coverage   72.66%   72.90%   +0.24%     
==========================================
  Files          17       17              
  Lines        5162     5212      +50     
==========================================
+ Hits         3751     3800      +49     
- Misses       1411     1412       +1     
Files with missing lines Coverage Δ
python/simulation.py 77.06% <98.03%> (+0.44%) ⬆️

Copy link
Collaborator

@oskooi oskooi left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Getting there. Some more suggested changes.

python/simulation.py Outdated Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
python/tests/test_ring.py Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
@oskooi
Copy link
Collaborator

oskooi commented Mar 16, 2023

You will also need to create a new entry for the Pade class in doc/docs/Python_User_Interface.md.in (similar to Harminv):

#### Harminv Step Function
The following step function collects field data from a given point and runs [Harminv](https://github.com/NanoComp/harminv) on that data to extract the frequencies, decay rates, and other information.
* [Harminv class](#Harminv)

This is necessary for converting its docstrings into Markdown for the readthedocs manual.

@stevengj
Copy link
Collaborator

In a separate PR, it would be good to show how to use this to calculate transmission and reflection spectra in problems with single-mode ports and sharp resonances. For example, to compute the transmission/reflection spectrum through N layers of a 2d photonic crystal at normal incidence.

python/simulation.py Outdated Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
python/simulation.py Outdated Show resolved Hide resolved
Copy link
Collaborator

@oskooi oskooi left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just two final comments/suggestions.

python/simulation.py Show resolved Hide resolved
python/simulation.py Show resolved Hide resolved
Copy link
Collaborator

@oskooi oskooi left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM.

@Arjun-Khurana Arjun-Khurana marked this pull request as ready for review March 31, 2023 14:56
@stevengj stevengj merged commit 1fe3899 into NanoComp:master Mar 31, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants