Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • wenxuan.jia/pygwinc
  • sean-leavey/pygwinc
  • sebastian.steinlechner/pygwinc
  • nicholas.demos/pygwinc
  • chris.whittle/pygwinc
  • raymond.robie/pygwinc
  • mateusz.bawaj/pygwinc
  • anchal.gupta/pygwinc
  • 40m/pygwinc
  • evan.hall/pygwinc
  • kevin.kuns/pygwinc
  • geoffrey-lovelace/pygwinc
  • brittany.kamai/pygwinc
  • daniel-brown/pygwinc
  • lee-mcculler/pygwinc
  • jameson.rollins/pygwinc
  • gwinc/pygwinc
17 results
Show changes
Commits on Source (35)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1049 additions and 623 deletions
.gitignore
View file @ a5f2d45e
#for docs and setup.py outputs
build/
# test cache
gwinc/test/cache
# Byte-compiled / optimized / DLL files
__pycache__/
......
.gitlab-ci.yml
View file @ a5f2d45e
image: igwn/base:buster
stages:
- test
- deploy
- dist
- test
- review
- deploy
# have to specify this so that all jobs execute for all commits
# including merge requests
workflow:
rules:
- if: $CI_MERGE_REQUEST_ID
- if: $CI_COMMIT_BRANCH
variables:
GIT_STRATEGY: clone
test:
# build the docker image we will use in all the jobs, with all
# dependencies pre-installed/configured.
gwinc/base:
stage: dist
variables:
IMAGE_TAG: $CI_REGISTRY_IMAGE/$CI_JOB_NAME:$CI_COMMIT_REF_NAME
GIT_STRATEGY: none
script:
- docker login -u gitlab-ci-token -p $CI_JOB_TOKEN $CI_REGISTRY
- |
cat <<EOF > Dockerfile
FROM igwn/base:buster
RUN apt-get update -qq
RUN apt-get -y install --no-install-recommends git python3-gitlab python3 python3-yaml python3-scipy python3-matplotlib python3-lalsimulation python3-pypdf2 python3-h5py
RUN git clone https://gitlab-ci-token:ci_token@git.ligo.org/gwinc/inspiral_range.git
EOF
- docker build -t $IMAGE_TAG .
- docker push $IMAGE_TAG
# create plots for the canonical IFOs
generate_budgets:
stage: test
before_script:
- echo $CI_COMMIT_SHA | cut -b1-8 > gitID.txt
image: $CI_REGISTRY_IMAGE/gwinc/base:$CI_COMMIT_REF_NAME
script:
- rm -rf ifo gwinc_test_report.pdf
- mkdir ifo
- apt-get update -qq
- apt-get install -y -qq git python3-yaml python3-scipy python3-matplotlib python3-ipython lalsimulation-python3 python3-pypdf2 python3-h5py
- git clone https://gitlab-ci-token:ci_token@git.ligo.org/gwinc/inspiral_range.git
- export PYTHONPATH=inspiral_range
- export MPLBACKEND=agg
- python3 -m gwinc.test -r gwinc_test_report.pdf
- for ifo in aLIGO Aplus Voyager CE1 CE2; do
- python3 -m gwinc $ifo -s ifo/$ifo.png
- python3 -m gwinc $ifo -s ifo/$ifo.h5
- done
- python3 -m gwinc.ifo -s ifo/all_compare.png
after_script:
- rm gitID.txt
cache:
key: "$CI_PROJECT_NAMESPACE:$CI_PROJECT_NAME:$CI_JOB_NAME"
untracked: true
- mkdir -p ifo
- export PYTHONPATH=/inspiral_range
- for ifo in $(python3 -c "import gwinc; print(' '.join(gwinc.IFOS))"); do
- python3 -m gwinc $ifo -s ifo/$ifo.png -s ifo/$ifo.h5
- done
- python3 -m gwinc.ifo -s ifo/all_compare.png
artifacts:
when: always
expire_in: 4w
paths:
- ifo
- gwinc_test_report.pdf
- ifo
# this is a special job intended to run only for merge requests.
# budgets are compared against those from the target branch. if the
# merge request has not yet been approved and noise changes are found,
# the job will fail. once the merge request is approved the job can
# be re-run, at which point the pipeline should succeed allowing the
# merge to be merged.
noise_change_approval:
stage: review
rules:
# - if: '$CI_MERGE_REQUEST_TARGET_BRANCH_NAME == "master"'
- if: $CI_MERGE_REQUEST_ID
image: $CI_REGISTRY_IMAGE/gwinc/base:$CI_COMMIT_REF_NAME
script:
- export PYTHONPATH=/inspiral_range
- |
cat <<EOF > check_approved.py
import sys
import gitlab
project_id = sys.argv[1]
mr_iid = sys.argv[2]
# this only works for public repos, otherwise need to specify
# private_token=
gl = gitlab.Gitlab('https://git.ligo.org')
project = gl.projects.get(project_id)
mr = project.mergerequests.get(mr_iid)
approvals = mr.approvals.get()
print(approvals.approved)
EOF
- echo CI_MERGE_REQUEST_PROJECT_ID=$CI_MERGE_REQUEST_PROJECT_ID
- echo CI_MERGE_REQUEST_IID=$CI_MERGE_REQUEST_IID
- echo CI_MERGE_REQUEST_TARGET_BRANCH_NAME=$CI_MERGE_REQUEST_TARGET_BRANCH_NAME
- approved=$(python3 check_approved.py $CI_MERGE_REQUEST_PROJECT_ID $CI_MERGE_REQUEST_IID )
- if [[ $approved != True ]] ; then
- echo "Approval not yet given, checking for noise changes..."
- target=origin/$CI_MERGE_REQUEST_TARGET_BRANCH_NAME
- if ! python3 -m gwinc.test --git-rev $target -r gwinc_test_report.pdf ; then
- echo "NOISE CHANGES RELATIVE TO $CI_MERGE_REQUEST_TARGET_BRANCH_NAME."
- echo "Approval required to merge this branch."
- /bin/false
- else
- echo "No noise changes detected, merge may proceed."
- fi
- else
- echo "Merge request approved, noise change accepted."
- fi
artifacts:
when: on_failure
paths:
- gwinc_test_report.pdf
expose_as: 'PDF report of noise changes relative to target branch'
# generate the html doc web pages. the "pages" job has special
# meaning, as it's "public" artifact becomes the directory served
# through gitlab static pages
pages:
stage: deploy
dependencies:
- test
only:
- master
needs:
- job: generate_budgets
artifacts: true
image: $CI_REGISTRY_IMAGE/gwinc/base:$CI_COMMIT_REF_NAME
script:
- rm -rf public
- mv ifo public
- apt-get install -y -qq python3-sphinx-rtd-theme
- cd docs
- make html
- cd ..
- mv ./build/sphinx/html/* public/
- rm -rf public
- apt-get install -y -qq python3-sphinx-rtd-theme
- cd docs
- make html
- cd ..
- mv ./build/sphinx/html public
- mv ifo public/
artifacts:
when: always
paths:
- public
expire_in: 4w
only:
- master
- public
CONTRIBUTING.md
View file @ a5f2d45e
......@@ -2,25 +2,48 @@
The pygwinc project welcomes your contributions. Our policy is that
all contributions should be peer-reviewed. To facilitate the review,
please do not commit your work to this repository yourself. Instead,
fork the repository and send a merge request.
please do not commit your work to this repository directly. Instead,
please fork the repository and create a [merge
request](https://git.ligo.org/gwinc/pygwinc/-/merge_requests/new)
against the main pygwinc master branch.
`pygwinc` comes with a test suite that will calculate all canonical
IFO budgets with the current state of the code, and compare them
against cached hdf5 traces (stored in the repo with
[git-lfs](https://git-lfs.github.com/) ( [see
also](https://wiki.ligo.org/Computing/GitLFS)). Contributors should
run the tests before submitting any merge requests:
When submitting code for merge, please follow good coding practice.
Respect the existing coding style, which for `pygwinc` is standard
[PEP8](https://www.python.org/dev/peps/pep-0008/) (with some
exceptions). Make individual commits as logically distinct and atomic
as possible, and provide a complete, descriptive log of the changes
(including a top summary line). Review efficiency follows code
legibility.
`pygwinc` comes with a validation command that can compare budgets
from the current code against those produced from different versions
in git (by default it compares against the current HEAD). The command
can be run with:
```shell
$ python3 -m gwinc.test
```
The test suite will be run automatically upon push to git.ligo.org as
part of the continuous integration, and code that fails the CI will
not be merged, unless failures are well justified.
Use the '--plot' or '--report' options to produce visual comparisons
of the noise differences. The comparison can be done against an
arbitrary commit using the '--git-ref' option. Traces for referenced
commits are cached, which speeds up subsequent comparison runs
significantly.
Once you submit your merge request a special CI job will determine if
there are budgets differences between your code and the master branch.
If there are, explicit approval from reviewers will be required before
your changes can be merged (see "approving noise" below).
## For reviewers: approving noise curve changes
If your change affects the shape of a noise curve, your commit message
should make note of that, and provide a justification. It will also
be necessary to update the reference curves stored in cache, but that
should be done in a separate commit not a part of the original merge
request, so that all reviewers can see the changes introduced in the
CI test failure report.
As discussed above, merge requests that generate noise changes will
cause a pipeline failure in the `review:noise_change_approval` CI job.
The job will generate a report comparing the new noise traces against
those from master, which can be found under the 'View exposed
artifacts' menu item in the pipeline report. Once you have reviewed
the report and the code, and understand and accept the noise changes,
click the 'Approve' button in the MR. Once sufficient approval has
been given, `review:noise_change_approval` job can be re-run, which
should now pick up that approval has been given and allow the pipeline
to succeed. Once the pipeline succeeds the merge request can be
merged. Click the 'Merge' button to finally merge the code.
IFO.md
View file @ a5f2d45e
......@@ -2,43 +2,43 @@
CI-generated plots and data for all IFOs included in pygwinc.
![IFO compare](https://gwinc.docs.ligo.org/pygwinc/all_compare.png)
![IFO compare](https://gwinc.docs.ligo.org/pygwinc/ifo/all_compare.png)
## aLIGO
* [ifo.yaml](gwinc/ifo/aLIGO/ifo.yaml)
* [aLIGO.h5](https://gwinc.docs.ligo.org/pygwinc/aLIGO.h5)
* [aLIGO.h5](https://gwinc.docs.ligo.org/pygwinc/ifo/aLIGO.h5)
![aLIGO](https://gwinc.docs.ligo.org/pygwinc/aLIGO.png)
![aLIGO](https://gwinc.docs.ligo.org/pygwinc/ifo/aLIGO.png)
## A+
* [ifo.yaml](gwinc/ifo/Aplus/ifo.yaml)
* [Aplus.h5](https://gwinc.docs.ligo.org/pygwinc/Aplus.h5)
* [Aplus.h5](https://gwinc.docs.ligo.org/pygwinc/ifo/Aplus.h5)
![Aplus](https://gwinc.docs.ligo.org/pygwinc/Aplus.png)
![Aplus](https://gwinc.docs.ligo.org/pygwinc/ifo/Aplus.png)
## Voyager
* [ifo.yaml](gwinc/ifo/Voyager/ifo.yaml)
* [Voyager.h5](https://gwinc.docs.ligo.org/pygwinc/Voyager.h5)
* [Voyager.h5](https://gwinc.docs.ligo.org/pygwinc/ifo/Voyager.h5)
![Voyager](https://gwinc.docs.ligo.org/pygwinc/Voyager.png)
![Voyager](https://gwinc.docs.ligo.org/pygwinc/ifo/Voyager.png)
## Cosmic Explorer 1
* [ifo.yaml](gwinc/ifo/CE1/ifo.yaml)
* [CE1.h5](https://gwinc.docs.ligo.org/pygwinc/CE1.h5)
* [CE1.h5](https://gwinc.docs.ligo.org/pygwinc/ifo/CE1.h5)
![CE1](https://gwinc.docs.ligo.org/pygwinc/CE1.png)
![CE1](https://gwinc.docs.ligo.org/pygwinc/ifo/CE1.png)
## Cosmic Explorer 2
* [ifo.yaml](gwinc/ifo/CE2/ifo.yaml)
* [CE2.h5](https://gwinc.docs.ligo.org/pygwinc/CE2.h5)
* [CE2.h5](https://gwinc.docs.ligo.org/pygwinc/ifo/CE2.h5)
![CE2](https://gwinc.docs.ligo.org/pygwinc/CE2.png)
![CE2](https://gwinc.docs.ligo.org/pygwinc/ifo/CE2.png)
README.md
View file @ a5f2d45e
......@@ -2,12 +2,14 @@
# Python Gravitational Wave Interferometer Noise Calculator
![gwinc](https://gwinc.docs.ligo.org/pygwinc/aLIGO.png)
[![aLIGO](https://gwinc.docs.ligo.org/pygwinc/ifo/aLIGO.png "Canonical
IFOs")](IFO.md)
`pygwinc` is a multi-faceted tool for processing and plotting noise
budgets for ground-based gravitational wave detectors. It's primary
feature is a collection of mostly analytic noise calculation functions
for various sources of noise affecting detectors (`gwinc.noise`):
feature is a collection of mostly analytic [noise calculation
functions](#noise-functions) for various sources of noise affecting
detectors (`gwinc.noise`):
* quantum noise
* mirror coating thermal noise
......@@ -17,22 +19,20 @@ for various sources of noise affecting detectors (`gwinc.noise`):
* Newtonian/gravity-gradient noise
* residual gas noise
See [noise functions](#noise-functions) below.
`pygwinc` also includes a generalized noise budgeting tool
(`gwinc.nb`) that allows users to create arbitrary noise budgets (for
any experiment, not just ground-based GW detectors) using measured or
analytically calculated data. See the [budget
interface](#Budget-interface) section below.
`pygwinc` is also a generalized noise budgeting tool (`gwinc.nb`) that
allows users to create arbitrary noise budgets (for any experiment,
not just ground-based GW detectors) using measured or analytically
calculated data. See the [budget interface](#Budget-interface)
section below.
`pygwinc` includes canonical budgets for various well-known current
and future detectors (`gwinc.ifo`):
and future GW detectors (`gwinc.ifo`):
* [aLIGO](https://gwinc.docs.ligo.org/pygwinc/aLIGO.png)
* [A+](https://gwinc.docs.ligo.org/pygwinc/Aplus.png)
* [Voyager](https://gwinc.docs.ligo.org/pygwinc/Voyager.png)
* [Cosmic Explorer 1](https://gwinc.docs.ligo.org/pygwinc/CE1.png)
* [Cosmic Explorer 2](https://gwinc.docs.ligo.org/pygwinc/CE2.png)
* [aLIGO](https://gwinc.docs.ligo.org/pygwinc/ifo/aLIGO.png)
* [A+](https://gwinc.docs.ligo.org/pygwinc/ifo/Aplus.png)
* [Voyager](https://gwinc.docs.ligo.org/pygwinc/ifo/Voyager.png)
* [Cosmic Explorer 1](https://gwinc.docs.ligo.org/pygwinc/ifo/CE1.png)
* [Cosmic Explorer 2](https://gwinc.docs.ligo.org/pygwinc/ifo/CE2.png)
See [IFO.md](IFO.md) for the latest CI-generated plots and hdf5 cached
data.
......@@ -68,6 +68,12 @@ Materials.Coating.Philown 5e-05 3e-05
$ python3 -m gwinc my_aLIGO.yaml
```
You can also use the `--ifo` option to change parameters from the
command line:
```shell
$ python3 -m gwinc aLIGO --ifo Optics.SRM.Tunephase=3.14
```
Stand-alone YAML files will always assume the nominal ['aLIGO' budget
description](gwinc/ifo/aLIGO).
......@@ -325,13 +331,7 @@ traces = budget.run()
```
The IFOs included in `gwinc.ifo` provide examples of the use of the
budget interface:
* [aLIGO](gwinc/ifo/aLIGO)
* [Aplus](gwinc/ifo/Aplus)
* [Voyager](gwinc/ifo/Voyager)
* [CE1](master/gwinc/ifo/CE1)
* [CE2](master/gwinc/ifo/CE2)
budget interface (e.g. [gwinc.ifo.aLIGO](gwinc/ifo/aLIGO)).
### extracting single noise terms
......
gwinc/__main__.py
View file @ a5f2d45e
......@@ -3,7 +3,6 @@ import os
import signal
import argparse
import numpy as np
from IPython.terminal.embed import InteractiveShellEmbed
import logging
logging.basicConfig(
......@@ -34,8 +33,12 @@ input file (IFO) can be an HDF5 file saved from a previous call, in
which case all noise traces and IFO parameters will be loaded from
that file.
Individual IFO parameters can be overriden with the --ifo option:
gwinc --ifo Optics.SRM.Tunephase=3.14 ...
If the inspiral_range package is installed, various figures of merit
can be calculated for the resultant spectrum with the --fom argument,
can be calculated for the resultant spectrum with the --fom option,
e.g.:
gwinc --fom horizon ...
......@@ -46,39 +49,60 @@ See documentation for inspiral_range package for details.
"""
IFO = 'aLIGO'
FLO = 5
FHI = 6000
NPOINTS = 3000
FREQ = '5:3000:6000'
DATA_SAVE_FORMATS = ['.hdf5', '.h5']
parser = argparse.ArgumentParser(
prog='gwinc',
description=description,
formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument('--flo', '-fl', default=FLO, type=float,
help="lower frequency bound in Hz [{}]".format(FLO))
parser.add_argument('--fhi', '--fh', default=FHI, type=float,
help="upper frequency bound in Hz [{}]".format(FHI))
parser.add_argument('--npoints', '-n', default=NPOINTS,
help="number of frequency points [{}]".format(NPOINTS))
parser.add_argument('--title', '-t',
help="plot title")
parser.add_argument('--fom',
help="calculate inspiral range for resultant spectrum ('func[:param=val,param=val]')")
parser.add_argument(
'--freq', '-f', metavar='FSPEC',
help="frequency array specification in Hz, either as 'flo:fhi' or 'flo:npoints:fhi' [{}]".format(FREQ))
parser.add_argument(
'--ifo', '-o',
#nargs='+', action='extend',
action='append',
help="override budget IFO parameter (may be specified multiple times)")
parser.add_argument(
'--title', '-t',
help="plot title")
parser.add_argument(
'--fom',
help="calculate inspiral range for resultant spectrum ('func[:param=val,param=val]')")
group = parser.add_mutually_exclusive_group()
group.add_argument('--interactive', '-i', action='store_true',
help="interactive plot with interactive shell")
group.add_argument('--save', '-s',
help="save budget traces (.hdf5/.h5) or plot (.pdf/.png/.svg) to file")
group.add_argument('--yaml', '-y', action='store_true',
help="print IFO as yaml to stdout and exit")
group.add_argument('--text', '-x', action='store_true',
help="print IFO as text table to stdout and exit")
group.add_argument('--diff', '-d', metavar='IFO',
help="show differences table between another IFO description")
group.add_argument('--no-plot', '-np', action='store_false', dest='plot',
help="supress plotting")
parser.add_argument('IFO',
help="IFO name, description file path (.yaml, .mat, .m), budget module (.py), or HDF5 data file (.hdf5, .h5)")
group.add_argument(
'--interactive', '-i', action='store_true',
help="interactive plot with interactive shell")
group.add_argument(
'--save', '-s', action='append',
help="save plot (.png/.pdf/.svg) or budget traces (.hdf5/.h5) to file (may be specified multiple times)")
group.add_argument(
'--yaml', '-y', action='store_true',
help="print IFO as yaml to stdout and exit (budget not calculated)")
group.add_argument(
'--text', '-x', action='store_true',
help="print IFO as text table to stdout and exit (budget not calculated)")
group.add_argument(
'--diff', '-d', metavar='IFO',
help="show differences table between another IFO description and exit (budget not calculated)")
group.add_argument(
'--no-plot', '-np', action='store_false', dest='plot',
help="supress plotting")
parser.add_argument(
'IFO',
help="IFO name, or path to budget module (.py), description file (.yaml/.mat/.m), or HDF5 data file (.hdf5/.h5)")
def freq_from_spec(spec):
fspec = spec.split(':')
if len(fspec) == 2:
fspec = fspec[0], FREQ.split(':')[1], fspec[1]
return np.logspace(
np.log10(float(fspec[0])),
np.log10(float(fspec[2])),
int(fspec[1]),
)
def main():
......@@ -89,22 +113,41 @@ def main():
##########
# initial arg processing
if os.path.splitext(os.path.basename(args.IFO))[1] in ['.hdf5', '.h5']:
if os.path.splitext(os.path.basename(args.IFO))[1] in DATA_SAVE_FORMATS:
from .io import load_hdf5
Budget = None
freq, traces, attrs = load_hdf5(args.IFO)
ifo = getattr(attrs, 'IFO', None)
plot_style = attrs
if args.freq:
logging.warning("ignoring frequency specification for frequencies defined in HDF5...")
else:
Budget = load_budget(args.IFO)
ifo = Budget.ifo
# FIXME: this should be done only if specified, to allow for
# using any FREQ specified in the Budget
freq = np.logspace(np.log10(args.flo), np.log10(args.fhi), args.npoints)
if args.freq:
try:
freq = freq_from_spec(args.freq)
except IndexError:
parser.error("improper frequency specification '{}'".format(args.freq))
else:
freq = getattr(Budget, 'freq', freq_from_spec(FREQ))
plot_style = getattr(Budget, 'plot_style', {})
traces = None
out_data_files = set()
out_plot_files = set()
if args.save:
for path in args.save:
if os.path.splitext(path)[1] in DATA_SAVE_FORMATS:
out_data_files.add(path)
out_plot_files = set(args.save) - out_data_files
if args.ifo:
for paramval in args.ifo:
param, val = paramval.split('=', 1)
ifo[param] = float(val)
if args.yaml:
if not ifo:
parser.exit(2, "no IFO structure available.")
......@@ -199,31 +242,12 @@ def main():
plot_style['title'] += '\n{}'.format(fom_title)
##########
# output
# save noise traces to HDF5 file
if args.save and os.path.splitext(args.save)[1] in ['.hdf5', '.h5']:
from .io import save_hdf5
logging.info("saving budget traces {}...".format(args.save))
if ifo:
plot_style['IFO'] = ifo.to_yaml()
save_hdf5(
path=args.save,
freq=freq,
traces=traces,
**plot_style
)
# interactive
# interactive shell plotting
elif args.interactive:
if args.interactive:
from IPython.terminal.embed import InteractiveShellEmbed
ipshell = InteractiveShellEmbed(
user_ns={
'freq': freq,
'traces': traces,
'ifo': ifo,
'plot_style': plot_style,
'plot_noise': plot_noise,
},
banner1='''
GWINC interactive plotter
......@@ -231,14 +255,39 @@ You may interact with plot using "plt." methods, e.g.:
>>> plt.title("foo")
>>> plt.savefig("foo.pdf")
''')
''',
user_ns={
'freq': freq,
'traces': traces,
'ifo': ifo,
'plot_style': plot_style,
'plot_noise': plot_noise,
},
)
ipshell.enable_pylab()
ipshell.run_code("plot_noise(freq, traces, **plot_style)")
ipshell.run_code("plt.title('{}')".format(plot_style['title']))
ipshell.ex("fig = plot_noise(freq, traces, **plot_style)")
ipshell.ex("plt.title('{}')".format(plot_style['title']))
ipshell()
##########
# output
# save noise traces to HDF5 file
if out_data_files:
from .io import save_hdf5
attrs = dict(plot_style)
attrs['IFO'] = ifo.to_yaml()
for path in out_data_files:
logging.info("saving budget traces {}...".format(path))
save_hdf5(
path=path,
freq=freq,
traces=traces,
**attrs,
)
# standard plotting
elif args.plot:
if args.plot or out_plot_files:
logging.info("plotting noises...")
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
......@@ -249,10 +298,9 @@ You may interact with plot using "plt." methods, e.g.:
**plot_style
)
fig.tight_layout()
if args.save:
fig.savefig(
args.save,
)
if out_plot_files:
for path in out_plot_files:
fig.savefig(path)
else:
plt.show()
......
gwinc/ifo/noises.py
View file @ a5f2d45e
......@@ -141,21 +141,9 @@ def precomp_suspension(f, ifo):
ifo.Suspension.VHCoupling = Struct()
ifo.Suspension.VHCoupling.theta = ifo.Infrastructure.Length / const.R_earth
if ifo.Suspension.Type == 'BQuad':
material = 'Silicon'
else:
material = 'Silica'
hForce, vForce, hTable, vTable = suspension.suspQuad(f, ifo.Suspension, material)
try:
# full TF (conventional)
ifo.Suspension.hForce = hForce.fullylossy
ifo.Suspension.vForce = vForce.fullylossy
# TFs with each stage lossy
ifo.Suspension.hForce_singlylossy = hForce.singlylossy
ifo.Suspension.vForce_singlylossy = vForce.singlylossy
except:
ifo.Suspension.hForce = hForce
ifo.Suspension.vForce = vForce
hForce, vForce, hTable, vTable = suspension.suspQuad(f, ifo.Suspension)
ifo.Suspension.hForce = hForce
ifo.Suspension.vForce = vForce
ifo.Suspension.hTable = hTable
ifo.Suspension.vTable = vTable
......
gwinc/noise/substratethermal.py
View file @ a5f2d45e
......@@ -61,8 +61,6 @@ def substrate_carrierdensity(f, materials, wBeam, exact=False):
def substrate_thermorefractive(f, materials, wBeam, exact=False):
"""Substrate thermal displacement noise spectrum from thermorefractive fluctuations
For semiconductor substrates.
:f: frequency array in Hz
:materials: gwinc optic materials structure
:wBeam: beam radius (at 1 / e^2 power)
......
gwinc/noise/suspensionthermal.py
View file @ a5f2d45e
......@@ -5,11 +5,12 @@ from __future__ import division
import numpy as np
from numpy import pi, imag
from .. import const
from ..const import kB
from ..suspension import getJointParams
def suspension_thermal(f, sus):
"""Suspention thermal noise for a single suspended test mass
"""Suspension thermal noise.
:f: frequency array in Hz
:sus: gwinc suspension structure
......@@ -17,67 +18,45 @@ def suspension_thermal(f, sus):
:returns: displacement noise power spectrum at :f:
:Temp: must either be set for each stage individually, or globally
in :sus.Temp:. The latter will be preferred if specified, so if
you wish to use per-stage tempurature you must remove :sus.Temp:.
in :sus.Temp:. If both are specified, :sus.Temp: is interpreted as
the temperature of the upper joint of the top stage.
Assumes suspension transfer functions and V-H coupling have been
pre-calculated and populated into the relevant `sus` fields.
"""
# Assign Physical Constants
kB = const.kB
# and vertical to beamline coupling angle
theta = sus.VHCoupling.theta
noise = np.zeros((1, f.size))
noise = np.zeros_like(f)
# if the temperature is uniform along the suspension
if 'Temp' in sus:
##########################################################
# Suspension TFs
##########################################################
##########################################################
# Suspension TFs
##########################################################
hForce = sus.hForce
vForce = sus.vForce
hForce = sus.hForce
vForce = sus.vForce
##########################################################
# Thermal Noise Calculation
##########################################################
# convert to beam line motion
# theta is squared because we rotate by theta into the suspension
# basis, and by theta to rotate back to the beam line basis
dxdF = hForce + theta**2 * vForce
# thermal noise (m^2/Hz) for one suspension
w = 2*pi*f
noise = 4 * kB * sus.Temp * abs(imag(dxdF)) / w
# if the temperature is set for each suspension stage
else:
##########################################################
# Suspension TFs
##########################################################
hForce = sus.hForce_singlylossy[:, :]
vForce = sus.vForce_singlylossy[:, :]
##########################################################
# Thermal Noise Calculation
##########################################################
dxdF = np.zeros(hForce.shape, dtype=complex)
for n, stage in enumerate(reversed(sus.Stage)):
# add up the contribution from each stage
##########################################################
# Thermal Noise Calculation
##########################################################
# Here the suspension stage list is reversed so that the last stage
# in the list is the test mass
stages = sus.Stage[::-1]
w = 2*pi*f
dxdF = np.zeros_like(hForce) # complex valued
for n, stage in enumerate(stages):
# add up the contribution from each joint
jointParams = getJointParams(sus, n)
for (isLower, Temp, wireMat, bladeMat) in jointParams:
# convert to beam line motion. theta is squared because
# we rotate by theta into the suspension basis, and by
# theta to rotate back to the beam line basis
dxdF[n, :] = hForce[n, :] + theta**2 * vForce[n, :]
# thermal noise (m^2/Hz) for one suspension
w = 2*pi*f
noise += 4 * kB * stage.Temp * abs(imag(dxdF[n, :])) / w
m = 2*n + isLower
dxdF[m, :] = hForce[m, :] + theta**2 * vForce[m, :]
noise += 4 * kB * Temp * abs(imag(dxdF[m, :])) / w
return np.squeeze(noise)
# thermal noise (m^2/Hz) for one suspension
return noise
gwinc/suspension.py
View file @ a5f2d45e
This diff is collapsed.
gwinc/test/__main__.py
View file @ a5f2d45e
import os
import sys
import glob
import shutil
import signal
import logging
import tempfile
import argparse
import subprocess
import numpy as np
import matplotlib.pyplot as plt
from collections import OrderedDict
from collections.abc import Mapping
from PyPDF2 import PdfFileReader, PdfFileWriter
logging.basicConfig(format='%(message)s',
level=os.getenv('LOG_LEVEL', logging.INFO))
from .. import IFOS, load_budget
from ..io import load_hdf5, save_hdf5
from ..io import load_hdf5
logging.basicConfig(
format='%(message)s',
level=os.getenv('LOG_LEVEL', logging.INFO))
try:
import inspiral_range
except ImportError:
logging.warning("inspiral_range package not found, range will not be calculated")
inspiral_range = None
TOLERANCE = 1e-6
CACHE_PATH = os.path.join(os.path.dirname(__file__), 'cache')
CACHE_LIMIT = 5
def test_path(*args):
"""Return path to package file."""
return os.path.join(os.path.dirname(__file__), *args)
def git_find_upstream_name():
try:
remotes = subprocess.run(
['git', 'remote', '-v'],
capture_output=True, universal_newlines=True,
check=True,
).stdout
except subprocess.CalledProcessError as e:
logging.error(e.stderr.split('\n')[0])
for remote in remotes.split('\n'):
name, url, fp = remote.split()
if 'gwinc/pygwinc.git' in url:
return name
def git_rev_resolve_hash(git_rev):
"""Resolve a git revision into its hash string."""
try:
return subprocess.run(
['git', 'show', '-s', '--format=format:%H', git_rev],
capture_output=True, universal_newlines=True,
check=True,
).stdout
except subprocess.CalledProcessError as e:
logging.error(e.stderr.split('\n')[0])
def prune_cache_dir():
"""Prune all but the N most recently accessed caches.
"""
cache_dir = test_path('cache')
if not os.path.exists(cache_dir):
return
expired_paths = sorted(
[os.path.join(cache_dir, path) for path in os.listdir(cache_dir)],
key=lambda path: os.stat(path).st_atime, reverse=True,
)[CACHE_LIMIT:]
if not expired_paths:
return
for path in expired_paths:
logging.info("pruning old cache: {}".format(path))
shutil.rmtree(path)
def gen_cache(git_hash, path):
"""generate cache for specified git hash at the specified path
The included shell script is used to extract the gwinc code from
the appropriate git commit, and invoke a new python instance to
generate the noise curves.
"""
logging.info("creating new cache for hash {}...".format(git_hash))
subprocess.run(
[test_path('gen_cache.sh'), git_hash, path],
check=True,
)
def load_cache(path):
"""load a cache from the specified path
returns a "cache" dictionary with 'git_hash' and 'ifos' keys.
"""
logging.info("loading cache {}...".format(path))
cache = {}
git_hash_path = os.path.join(path, 'git_hash')
if os.path.exists(git_hash_path):
with open(git_hash_path) as f:
git_hash = f.read().strip()
else:
git_hash = None
logging.debug("cache hash: {}".format(git_hash))
cache['git_hash'] = git_hash
cache['ifos'] = {}
for f in sorted(os.listdir(path)):
name, ext = os.path.splitext(f)
if ext != '.h5':
continue
cache['ifos'][name] = os.path.join(path, f)
return cache
def walk_traces(traces, root=()):
......@@ -105,17 +199,15 @@ def compare_traces(tracesA, tracesB, tolerance=TOLERANCE, skip=None):
return diffs
def plot_diffs(freq, diffs, tolerance,
name, labelA, labelB, fom_title='',
save=None):
def plot_diffs(freq, diffs, styleA, styleB):
spec = (len(diffs)+1, 2)
sharex = None
for i, nname in enumerate(diffs):
noiseA, noiseB, frac = diffs[nname]
axl = plt.subplot2grid(spec, (i, 0), sharex=None)
axl.loglog(freq, np.sqrt(noiseA), label=labelA)
axl.loglog(freq, np.sqrt(noiseB), label=labelB)
axl.loglog(freq, np.sqrt(noiseA), **styleA)
axl.loglog(freq, np.sqrt(noiseB), **styleB)
axl.grid()
axl.legend(loc='upper right')
axl.set_ylabel(nname)
......@@ -134,133 +226,167 @@ def plot_diffs(freq, diffs, tolerance,
if i == 0:
axr.set_title("fractional difference")
plt.suptitle('''{} {}/{} noise comparison
(noises that differ by more than {} ppm)
{}'''.format(name, labelA, labelB, tolerance*1e6, fom_title))
axl.set_xlabel("frequency [Hz]")
axr.set_xlabel("frequency [Hz]")
plt.subplots_adjust(top=0.8, right=0.85, wspace=0.3)
if save:
pwidth = 10
pheight = (len(diffs) * 5) + 2
plt.gcf().set_size_inches(pwidth, pheight)
plt.savefig(save)
else:
plt.show()
##################################################
def main():
parser = argparse.ArgumentParser()
parser.add_argument('--tolerance', '-t', type=float, default=TOLERANCE,
help='fractional tolerance [{}]'.format(TOLERANCE))
parser.add_argument('--skip', '-k', metavar='NOISE', action='append',
help='traces to skip in comparison (multiple may be specified)')
parser.add_argument('--cache', '-c', metavar='PATH', default=CACHE_PATH,
help='specify alternate IFO traces cache path')
parser = argparse.ArgumentParser(
description="""GWINC noise validation
This command calculates the canonical noise budgets with the current
code and compares them against those calculated with code from a
specified git revision. You must be running from a git checkout of
the source for this to work. The command will fail if it detects any
noise differences. Plots or a PDF report of differences can be
generated with the '--plot' or '--report' commands respectively.
By default it will attempt to determine the git reference for upstream
master for your current configuration (usually 'origin/master' or
'upstream/master'). You may specify an arbitrary git revision with
the --git-rev command. For example, to compare against another
remote/branch use:
$ python3 -m gwinc.test --git-rev remote/dev-branch
or if you have uncommitted changes compare against the current head
with:
$ python3 -m gwinc.test -g HEAD
See gitrevisions(7) for various ways to refer to git revisions.
A cache of traces from reference git revisions will be stored in
gwinc/test/cache/<SHA1>. Old caches are automatically pruned.""",
formatter_class=argparse.RawDescriptionHelpFormatter,
)
parser.add_argument(
'--tolerance', '-t', type=float, default=TOLERANCE,
help="fractional tolerance of comparison [{}]".format(TOLERANCE))
parser.add_argument(
'--skip', '-k', metavar='NOISE', action='append',
help="traces to skip in comparison (multiple may be specified)")
rgroup = parser.add_mutually_exclusive_group()
rgroup.add_argument('--plot', '-p', action='store_true',
help='plot differences')
rgroup.add_argument('--report', '-r', metavar='REPORT.pdf',
help='create PDF report of test results (only created if differences found)')
rgroup.add_argument('--gen-cache', action='store_true',
help='update/create IFO traces cache directory')
parser.add_argument('ifo', metavar='IFO', nargs='*',
help='specific ifos to test (default all)')
rgroup.add_argument(
'--git-rev', '-g', metavar='REV',
help="specify specific git revision to compare against")
rgroup.add_argument(
'--head', '-gh', action='store_const', dest='git_rev', const='HEAD',
help="shortcut for '--git-rev HEAD'")
ogroup = parser.add_mutually_exclusive_group()
ogroup.add_argument(
'--plot', '-p', action='store_true',
help="show interactive plot differences")
ogroup.add_argument(
'--report', '-r', metavar='REPORT.pdf',
help="create PDF report of test results (only created if differences found)")
parser.add_argument(
'ifo', metavar='IFO', nargs='*',
help="specific ifos to test (default all)")
args = parser.parse_args()
if args.gen_cache:
try:
os.makedirs(args.cache)
except FileExistsError:
pass
freq = np.logspace(np.log10(5), np.log10(6000), 3000)
for name in IFOS:
Budget = load_budget(name)
traces = Budget(freq).run()
path = os.path.join(args.cache, name+'.h5')
save_hdf5(path, freq, traces)
return
# get the reference hash
if args.git_rev:
git_rev = args.git_rev
else:
remote = git_find_upstream_name()
if not remote:
sys.exit("Could not resolve upstream remote name")
git_rev = '{}/master'.format(remote)
logging.info("git rev: {}".format(git_rev))
git_hash = git_rev_resolve_hash(git_rev)
if not git_hash:
sys.exit("Could not resolve reference, could not run test.")
logging.info("git hash: {}".format(git_hash))
# load the cache
cache_path = test_path('cache', git_hash)
if not os.path.exists(cache_path):
prune_cache_dir()
gen_cache(git_hash, cache_path)
cache = load_cache(cache_path)
if args.report:
base, ext = os.path.splitext(args.report)
if ext != '.pdf':
parser.error("Test reports only support PDF format.")
outdir = tempfile.TemporaryDirectory()
# find all cached IFOs
logging.info("loading cache {}...".format(args.cache))
cached_ifos = {}
for f in sorted(os.listdir(args.cache)):
name, ext = os.path.splitext(f)
if ext != '.h5':
continue
cached_ifos[name] = os.path.join(args.cache, f)
# select
if args.ifo:
ifos = {name:cached_ifos[name] for name in args.ifo}
ifos = args.ifo
else:
ifos = cached_ifos
ifos = IFOS
labelA = 'cache'
labelB = 'head'
style_cache = dict(label='reference', linestyle='-')
style_head = dict(label='head', linestyle='--')
fail = False
# compare
for name, path in ifos.items():
for name in ifos:
logging.info("{} tests...".format(name))
freq, tracesA, attrs = load_hdf5(path)
try:
path = cache['ifos'][name]
except KeyError:
logging.warning("IFO {} not found in cache")
fail |= True
continue
freq, traces_cache, attrs = load_hdf5(path)
Budget = load_budget(name)
tracesB = Budget(freq).run()
traces_head = Budget(freq).run()
if inspiral_range:
totalA = tracesA['Total'][0]
totalB = tracesB['Total'][0]
total_cache = traces_cache['Total'][0]
total_head = traces_head['Total'][0]
range_func = inspiral_range.range
H = inspiral_range.waveform.CBCWaveform(freq)
fomA = range_func(freq, totalA, H=H)
tracesA['int73'] = inspiral_range.int73(freq, totalA)[1], None
fomB = range_func(freq, totalB, H=H)
tracesB['int73'] = inspiral_range.int73(freq, totalB)[1], None
fom_cache = range_func(freq, total_cache, H=H)
traces_cache['int73'] = inspiral_range.int73(freq, total_cache)[1], None
fom_head = range_func(freq, total_head, H=H)
traces_head['int73'] = inspiral_range.int73(freq, total_head)[1], None
fom_summary = """
inspiral {func} {m1}/{m2} Msol:
{labelA}: {fomA:.2f} Mpc
{labelB}: {fomB:.2f} Mpc
{label_cache}: {fom_cache:.2f} Mpc
{label_head}: {fom_head:.2f} Mpc
""".format(
func=range_func.__name__,
m1=H.params['m1'],
m2=H.params['m2'],
labelA=labelA,
fomA=fomA,
labelB=labelB,
fomB=fomB,
label_cache=style_cache['label'],
fom_cache=fom_cache,
label_head=style_head['label'],
fom_head=fom_head,
)
else:
fom_summary = ''
diffs = compare_traces(tracesA, tracesB, args.tolerance, args.skip)
diffs = compare_traces(traces_cache, traces_head, args.tolerance, args.skip)
if diffs:
logging.warning("{} tests FAIL".format(name))
fail |= True
if args.plot or args.report:
if args.report:
save = os.path.join(outdir.name, name+'.pdf')
else:
save = None
plot_diffs(
freq, diffs, args.tolerance,
name, labelA, labelB, fom_summary,
save=save,
)
else:
if not diffs:
logging.info("{} tests pass.".format(name))
continue
logging.warning("{} tests FAIL".format(name))
fail |= True
if args.plot or args.report:
plot_diffs(freq, diffs, style_cache, style_head)
plt.suptitle('''{} {}/{} noise comparison
(noises that differ by more than {} ppm)
reference git hash: {}
{}'''.format(name, style_cache['label'], style_head['label'],
args.tolerance*1e6, cache['git_hash'], fom_summary))
if args.report:
pwidth = 10
pheight = (len(diffs) * 5) + 2
plt.gcf().set_size_inches(pwidth, pheight)
plt.savefig(os.path.join(outdir.name, name+'.pdf'))
else:
plt.show()
if not fail:
logging.info("all tests pass.")
......
gwinc/test/cache/Aplus.h5 deleted 100644 → 0 LFS
View file @ 94de6c42
File deleted
gwinc/test/cache/CE1.h5 deleted 100644 → 0 LFS
View file @ 94de6c42
File deleted
gwinc/test/cache/CE2.h5 deleted 100644 → 0 LFS
View file @ 94de6c42
File deleted
gwinc/test/cache/Voyager.h5 deleted 100644 → 0 LFS
View file @ 94de6c42
File deleted
gwinc/test/cache/aLIGO.h5 deleted 100644 → 0 LFS
View file @ 94de6c42
File deleted
gwinc/test/gen_cache.sh 0 → 100755
View file @ a5f2d45e
#!/bin/bash -e
if [ -z "$1" ] || [ -z "$2" ] ; then
echo "usage: $(basename $0) git_hash cache_dir_path"
echo "generate a cache of IFO budget traces from a particular git commit"
exit 1
fi
git_hash="$1"
cache_dir="$2"
mkdir -p $cache_dir
cache_dir=$(cd $cache_dir && pwd)
gwinc_dir=$cache_dir/gwinc
mkdir -p $gwinc_dir
git archive $git_hash | tar -x -C $gwinc_dir
cd $gwinc_dir
export LOG_LEVEL=INFO
for ifo in $(python3 -c "import gwinc; print(' '.join(gwinc.IFOS))") ; do
python3 -m gwinc --save $cache_dir/${ifo}.h5 $ifo
done
echo $git_hash > $cache_dir/git_hash