Patching Dec 9, 2021 6-7a CST- All GitLab services may be unavailable for 5-10 minutes

Commit 0ac59b9e authored by Wei's avatar Wei
Browse files

Initial commit (with rompsline python2)

parents
# Byte-compiled / optimized / DLL files
__pycache__/
*.py[cod]
*$py.class
# C extensions
*.so
# Distribution / packaging
.Python
build/
develop-eggs/
dist/
downloads/
eggs/
.eggs/
lib/
lib64/
parts/
sdist/
var/
wheels/
*.egg-info/
.installed.cfg
*.egg
# PyInstaller
# Usually these files are written by a python script from a template
# before PyInstaller builds the exe, so as to inject date/other infos into it.
*.manifest
*.spec
# Installer logs
pip-log.txt
pip-delete-this-directory.txt
# Unit test / coverage reports
htmlcov/
.tox/
.coverage
.coverage.*
.cache
nosetests.xml
coverage.xml
*.cover
.hypothesis/
# Translations
*.mo
*.pot
# Django stuff:
*.log
local_settings.py
# Flask stuff:
instance/
.webassets-cache
# Scrapy stuff:
.scrapy
# Sphinx documentation
docs/_build/
# PyBuilder
target/
# Jupyter Notebook
.ipynb_checkpoints
# pyenv
.python-version
# celery beat schedule file
celerybeat-schedule
# SageMath parsed files
*.sage.py
# Environments
.env
.venv
env/
venv/
ENV/
# Spyder project settings
.spyderproject
.spyproject
# Rope project settings
.ropeproject
# mkdocs documentation
/site
# mypy
.mypy_cache/
#!/usr/bin/env python3
import errno
import subprocess
import mmap
# import romspline
import h5py
import numpy as np
import re
import sys
import os
import datetime
import romspline
try:
import ConfigParser # renamed to configparser in Python3
except ImportError:
import configparser as ConfigParser
class CactusConverter:
''' A class that converts output format from cactus to LIGO format. Internally it uses Mathematica script (with SimulationTools package) to generate metadata and waveform (hdf5), and then it uses python script to parse and combine metadata with waveform in LIGO compatible format.
@author Wei Ren
@param simulationPath - path to simulations being converted
outputPath - path to output directory
simulations - list of simulation names
nrGroup - Name of the group (default to "NCSA").
'''
def __init__(self, simulationPath, outputPath, simulationList, nrGroup = "NCSA"):
self.simulationPath = simulationPath
self.outputPath = outputPath
self.simulationList = simulationList
self.nrGroup = nrGroup
def toSXS(self):
'''
Export Cactus BBH simulations to SXS format.
@return - None.
@output - output files in the outputPath.
'''
procList = []
for sim in self.simulationList:
simPath = os.path.join(self.simulationPath, sim)
simfactoryPath = os.path.join(simPath, "SIMFACTORY")
if not os.path.isdir(simfactoryPath):
os.mkdir(simfactoryPath, 0o755)
proc = subprocess.Popen(["wolfram", "-nohup", "-script", "CactusToSXS.m", sim, self.outputPath, self.simulationPath])
procList.append(proc)
# Simple barrier
for proc in procList:
proc.wait()
self.__hackEccentricity()
print("Export to SXS format done!")
def toLIGO(self, exportSXS=True):
'''
Export Cactus BBH simulations to LIGO format. A boolean flag is also provided. It's useful when the user already exported SXS format, but the user must ensure SXS output files are in the outputPath otherwise the function will fail.
@param exportSXS - a boolean flag to export SXS (or not).
@return - None.
@output - output files in outputPath.
'''
if exportSXS:
self.toSXS()
for sim in self.simulationList:
self.__toLIGO(sim)
def exportStrainData(self, r_0 = 500, omega = 0.01):
'''
Export strain information (phase, frequency, amplitude) in *.dat files.
@param r_0 - default to 500.
@param omega - default to 0.01.
@return - None.
@output - output files in outputPath.
'''
procList = []
for sim in self.simulationList:
proc = subprocess.Popen(["wolfram", "-nohup", "-script", "exportStrainData.m", str(omega), str(r_0), sim, self.outputPath, self.simulationPath])
procList.append(proc)
for proc in procList:
proc.wait()
def __toLIGO(self, simName):
'''
Convert to LIGO format.
@param - Name of the one simulation.
@return - None.
@output - output hdf5 file with simulation information in outputPath.
'''
def getfloats(self, section, option):
# helper function to parse floats
return np.array([float(elem) for elem in re.split(", *", self.get(section, option))])
# see https://stackoverflow.com/a/2982
simPath = os.path.join(self.outputPath, simName)
metafile = os.path.join(simPath, "metadata.txt")
SXSfile = os.path.join(simPath, "rhOverM_Asymptotic_GeometricUnits.h5")
if not os.path.isfile(metafile):
raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), metafile)
ConfigParser.RawConfigParser.getfloats = getfloats
metadata = ConfigParser.RawConfigParser()
metadata.read(metafile)
outputLIGO = os.path.join(simPath, simName+".h5")
with h5py.File(outputLIGO, "w") as output:
attr = output.attrs
attr['Format'] = 1
attr['type'] = 'NRinjection'
# TODO: determine type (and name) based on information in metadata file
attr['name'] = '{group}:BBH:{simulation}'.format(
group=self.nrGroup, simulation=simName)
attr['alternative-names'] = '{group}:Binary Black Hole:{simulation}'.format(
group=self.nrGroup, simulation=simName)
attr['NR-group'] = self.nrGroup
attr['NR-code'] = 'Einstein Toolkit'
attr['modification-date'] = datetime.date.today().strftime('%Y/%m/%d')
attr['point-of-contact-email'] = metadata.get("metadata", "point-of-contact-email")
# TODO: check simulation-type based on spin values
attr['simulation-type'] = 'non-spinning'
attr['INSPIRE-bibtex-keys'] = ''
attr['license'] = 'LVC-internal' # TODO: choose one.
attr['Lmax'] = 2
attr['NR-techniques'] = '' # TODO: refer to the updated document
attr['files-in-error-series'] = ''
attr['comparable-simulation'] = ''
attr['production-run'] = 1
attr['object1'] = 'BH'
attr['object2'] = 'BH'
attr['PN_approximant'] = 'None'
mass1 = metadata.getfloat('metadata', 'relaxed-mass1')
attr['mass1'] = mass1
mass2 = metadata.getfloat('metadata', 'relaxed-mass2')
attr['mass2'] = mass2
spin1 = metadata.getfloats('metadata', 'relaxed-spin1')
attr['spin1x'] = spin1[0]
attr['spin1y'] = spin1[1]
attr['spin1z'] = spin1[2]
spin2 = metadata.getfloats('metadata', 'relaxed-spin2')
attr['spin2x'] = spin2[0]
attr['spin2y'] = spin2[1]
attr['spin2z'] = spin2[2]
position1 = metadata.getfloats('metadata', 'relaxed-position1')
position2 = metadata.getfloats('metadata', 'relaxed-position2')
n = position1 - position2 # from 2 to 1
nhat = n / np.linalg.norm(n)
attr['nhatx'] = nhat[0]
attr['nhaty'] = nhat[1]
attr['nhatz'] = nhat[2]
attr['eccentricity'] = metadata.getfloat(
'metadata', 'relaxed-eccentricity')
attr['mean_anomaly'] = metadata.getfloat(
'metadata', 'relaxed-mean-anomaly')
# hard-codes z axis as orbital axis
orbit_f = metadata.getfloats('metadata', 'relaxed-orbital-frequency')
attr['Omega'] = orbit_f[2]
f_lower_at_1MSUN = orbit_f[2] / (np.pi * 4.92549095e-6)
attr['f_lower_at_1MSUN'] = f_lower_at_1MSUN
attr['eta'] = mass1 * mass2 / (mass1 + mass2)**2
attr['LNhatx'] = 0.
attr['LNhaty'] = 0.
# NOTE: check based on spin vector instead of hardcode.
attr['LNhatz'] = 1.
relaxed_measurement_time = metadata.getfloat("metadata", "relaxed-measurement-time")
output.create_group('auxiliary-info')
with h5py.File(SXSfile, 'r') as fSXS:
extrapolatedGroups = list(filter(lambda n: re.match(
"Extrapolated_N[0-9]+\\.dir", n), fSXS.keys()))
print(type(extrapolatedGroups))
if(len(extrapolatedGroups) != 1):
raise IndexError(
"Could not identify unique data group in %s" % ",".join(fSXS.keys()))
extrapolatedGroup = extrapolatedGroups[0]
ds_regexp = "Y_l([0-9]+)_m(-?[0-9]+)\\.dat"
modes = [re.match(ds_regexp,ds).groups()[0:2]
for ds in fSXS[extrapolatedGroup] if re.match(ds_regexp, ds)]
for (l,m) in modes:
print("Exporting l=%s, m=%s" % (l, m))
ampGroup = 'amp_l%s_m%s' % (l, m)
phaseGroup = 'phase_l%s_m%s' % (l, m)
output.create_group(ampGroup)
output.create_group(phaseGroup)
strain = fSXS[extrapolatedGroup + "/Y_l%s_m%s.dat" % (l, m)][:]
# we need to enforce that the first data point is actually at t=relaxed_measurement_time
# TODO: if it's not there, use a spline interpolation to create it?
if not (relaxed_measurement_time >= strain[0, 0] and relaxed_measurement_time <= strain[-1, 0]) :
raise IndexError(
"relaxed_measurement_time %g is not a node of the time series for strain" %
relaxed_measurement_time)
mask = strain[:, 0] > relaxed_measurement_time # Truncate data before relaxed_measurement_time
strain = strain[mask, :]
t = strain[:, 0]
amplitude = np.sqrt(strain[:, 1]**2 + strain[:, 2]**2)
phase = np.unwrap(np.arctan2(strain[:, 2], strain[:, 1]))
spline = romspline.ReducedOrderSpline(t, amplitude, verbose=False)
spline.write(output[ampGroup])
spline = romspline.ReducedOrderSpline(t, phase, verbose=False)
spline.write(output[phaseGroup])
output.create_dataset('NRtimes', data=t)
def simulationPath():
doc = "The simulationPath property."
def fget(self):
return self._simulationPath
def fset(self, value):
if not os.path.isdir(value):
raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), value)
self._simulationPath = value
def fdel(self):
del self._simulationPath
return locals()
simulationPath = property(**simulationPath())
def outputPath():
doc = "The outputPath property."
def fget(self):
return self._outputPath
def fset(self, value):
if not os.path.isdir(value):
os.mkdir(value, 0o755)
self._outputPath = value
def fdel(self):
del self._outputPath
return locals()
outputPath = property(**outputPath())
def simulationList():
doc = "The simulationList property."
def fget(self):
return self._simulationList
def fset(self, value):
try:
if isinstance(value, str):
raise TypeError("'simulationList' should be a **list** of simulation names (strings).")
for elem in value:
sim = os.path.join(self.simulationPath, elem)
if not os.path.isdir(sim):
raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), sim)
except Exception as e:
print("'simulationList' is not a list of simulation names or it contains simulations that does not exists in 'simulationPath'.\n")
raise
self._simulationList = value
def fdel(self):
del self._simulationList
return locals()
simulationList = property(**simulationList())
def __hackEccentricity(self):
# Eccentricity is complex in the output metadata.txt for some reason.
# A hacky function to keep the real part.
for sim in self.simulationList:
simOutputPath = os.path.join(self.outputPath, sim)
metafile = os.path.join(simOutputPath, "metadata.txt")
with open(metafile, "r") as f:
tempString = f.read()
newString = tempString
match = re.search("(^.*?)(Complex[(])([-+]?\d*\.\d+|\d+)([,])([-+]?\d*\.\d+|\d+)([)])(.*?)$", tempString, re.DOTALL)
if match:
newString = match.group(1) + match.group(3) + match.group(7)
with open(metafile, "w") as f:
f.write(newString)
# Testing
if __name__ == '__main__':
a = CactusConverter("/home/wei/numericalRelativity/simulations", "/home/wei/numericalRelativity/temp", ["E0017_N32", "E0017_N36", "E0017_N40"], "NCSA")
print(a.simulationPath)
print(a.outputPath)
print(a.simulationList)
#!/usr/bin/env wolframscript
# Mathematica script to export Cactus BBH simulation to SXS format.
Get["SimulationTools`"];
# Parse command line argument
SimulationTools`$SimulationPath={$CommandLine[[-1]]};
sim=SimulationNames[$CommandLine[[-3]]];
outDir=$CommandLine[[-2]];
# Export hdf5 and metadata file
ExportSXSSimulation[sim[[1]],FileNameJoin[{outDir,sim[[1]]}]];
# Quit
LinkClose[#] & /@ Links[];
CloseKernels[];
Quit[];
# Cactus To LIGO
Convert the output format from cactus to LIGO so that it can be directly used by LIGO Algorithm Library (LAL).
# Prerequisite
- **python 3.5+**
- numpy
- h5py
- [romSpline](https://bitbucket.org/chadgalley/romspline) (no need to install since it's already packed in this repository for python3 compatibility)
- **Mathematica**
- [SimulationTools](http://www.simulationtools.org/)
- h5mma
## Installing Prerequisite
### Python 3.5+
Make sure `numpy` and `h5py` packages are installed in the python environment that will be running.
### Mathematica
#### `SimulationTools` and `h5mma`
The official website provides information how to [download](http://simulationtools.org/download.shtml) and install `SimulationTools` and `h5mma`. Below is a summary of instructions:
1. Locate the path to Mathematica Applications. To find out the path, you can run the following command in Mathematica:
```
FileNameJoin[{$UserBaseDirectory, "Applications"}]
```
*By default Mathematica will install the Applications Path to `$HOME/Library/Mathematica/Applications` (macOS) or `$HOME/.Mathematica/Applications` (Linux).*
2. Change directory to the Mathematica Applications path. For example, you can do the following if on Linux with default Mathematica installation:
```
cd $HOME/.Mathematica/Applications
```
*Note: this directory might not be created after fresh installation. In that case, you can manually create this directory or just open Mathematica at least once to make sure the path exists.*
3. Download and install `SimulationTools` and `h5mma`
```bash
# Check out SimulationTools
git clone --recursive https://bitbucket.org/simulationtools/SimulationTools.git
# Download h5mma
curl -L --remote-name "https://bitbucket.org/simulationtools/h5mma/downloads/h5mma-1.2.1.tar.gz"
tar -xzvf h5mma-1.2.1.tar.gz
```
4. Check out the NCSA-converter
```bash
git clone https://git.ncsa.illinois.edu/weiren2/NCSA-Converter.git
cd NCSA-Converter
```
# Running the Script
The NCSA-converter is already encapsulated into a **CactusConverter** class in **CactusConverter.py**. To use it, simply construct an **CactusConverter** object with correct parameters (most importantly the input and output path). Check comments in **CactusConverter.py** for more details on the functions.
The python script **convertToLIGO.py** shows an example that you can use to convert to LIGO format and extract relevant information (such as amplitude, phase...) from Cactus output. It also added support for command line inputs. Also see comments in **convertToLIGO.py** for some example uses.
Usage:
```bash
./convertToLIGO.py [group name] [simDir] [outdir] [simulation names]
```
## Some concrete examples
- Run from command line:
```bash
./convertToLIGO.py /path/to/simulations /path/to/output E0001_N32 E0001_N36 E0001_N40
```
- Or use it in a shell script or PBS job script
```bash
# Define input and output path below
INPUT_SIMULATION_PATH=/path/to/simulations
OUTPUT_PATH=/path/to/output
# Make sure qsub script is in the same directory with CactusToLIGO
cd $PBS_O_WORKDIR
# Remember to change the list of simulations at the end.
./convertToLIGO.py $INPUT_SIMULATION_PATH $OUTPUT_PATH E0001_N32 E0001_N36 E0001_N40
```
#### Author information ####
Authors: Wei Ren (*weiren2 "at" illinois "dot" edu*), Eliu Huerta (*elihu "at" illinois "dot" edu*), Roland Haas (*rhaas "at" illinois "dot" edu*), Ian Hinder (*ian.hinder "at" aei "dot" mpg "dot"de*)
Copyright (C) 2016 Trustees of the University of Illinois
Included Romspline
Copyright (C) 2016 Trustees of the University of Illinois
Copyright (C) 2015 Chad Galley (*crgalley "at" tapir "dot" caltech "dot" edu*).
# Issues
Please contact [ncsarelativitygroup@ncsa.illinois.edu](mailto:ncsarelativitygroup@ncsa.illinois.edu) for any bug or improvement.
## Improvements
- [ ] Remove the hardcoded part when filling in attributes.
- [x] ~~Make python script use command line argument.~~
#!/usr/bin/env python3
import argparse
import sys
import os
from CactusConverter import CactusConverter
def main(simulationPath, outputPath, simulationList):
converter = CactusConverter(simulationPath, outputPath, simulationList)
# SXS files are generated in outputPath.
converter.toSXS()
# Set to False if you can make sure SXS is already generated.
# (By default toLIGO() fuction will run toSXS anyway)
# (To avoid duplicating the work, we can set exportSXS flag to False.)
converter.toLIGO(exportSXS = False)
# exportStrainData() function exports amplitudes (real and imaginary part),
# phase and frequency raw data in separate files (.dat).
converter.exportStrainData()
if __name__ == '__main__':
if(len(sys.argv) < 4):
sys.stderr.write("usage: %s [group name] [simDir] [outdir] [simulation names]\n" % sys.argv[0])
sys.exit(1)
simDir = sys.argv[1]
outDir = sys.argv[2]
simNames = sys.argv[3:]
if not os.path.exists(simDir):
sys.stderr.write("%s: File %s not found\n" % (sys.argv[0], simDir))
sys.exit(1)
simDir = os.path.abspath(simDir)
if not os.path.exists(outDir):
sys.stderr.write("%s: File %s not found\n" % (sys.argv[0], outDir))
sys.exit(1)
outDir = os.path.abspath(outDir)
print(simDir)
print(outDir)
print(simNames)
main(simDir, outDir, simNames)
#!/usr/bin/env wolframscript
# A Mathematica script to export data in [2,2] mode.
Get["SimulationTools`"];
# Parse command line argument
SimulationTools`$SimulationPath={$CommandLine[[-1]]};
sim=$CommandLine[[-3]];
outDir=$CommandLine[[-2]];
r0 = ToExpression[$CommandLine[[-4]]];
\[Omega]0=ToExpression[$CommandLine[[-5]]];
# Export data
MADM = ReadADMMass[sim];
\[Psi]4[2,2]=ReadPsi4[sim,2,2,r0];
r\[Psi]4[2,2]=Shifted[r0 \[Psi]4[2,2],-RadialToTortoise[r0,MADM]];
dest = outDir<>"/"<>sim<>"/";
If[!DirectoryQ[dest],CreateDirectory[dest]];
Export[dest<>sim<>"_omega.dat", -Frequency[r\[Psi]4[2, 2]]];
Export[dest<>sim<>"_phase.dat",-Phase[r\[Psi]4[2,2]]];
rh[2,2]=Psi4ToStrain[r\[Psi]4[2,2],\[Omega]0];
Export[dest<>sim<>"_real.dat",Re[rh[2,2]]];
Export[dest<>sim<>"_imag.dat", Im[rh[2,2]]];
Export[dest<>sim<>"_amp.dat", Abs[r\[Psi]4[2,2]]]
# Quit
Print["Export strain data of " <> sim <> " done."];
LinkClose[#] & /@ Links[];
CloseKernels[];
Quit[];
# Byte-compiled / optimized / DLL files
__pycache__/
*.py[cod]
*$py.class
# C extensions
*.so
# Distribution / packaging
.Python
build/
develop-eggs/
dist/
downloads/
eggs/
.eggs/
lib/
lib64/
parts/
sdist/
var/
wheels/
*.egg-info/
.installed.cfg
*.egg
# PyInstaller
# Usually these files are written by a python script from a template
# before PyInstaller builds the exe, so as to inject date/other infos into it.
*.manifest
*.spec
# Installer logs
pip-log.txt
pip-delete-this-directory.txt
# Unit test / coverage reports
htmlcov/
.tox/
.coverage
.coverage.*
.cache
nosetests.xml
coverage.xml
*.cover
.hypothesis/
# Translations
*.mo
*.pot
# Django stuff:
*.log
local_settings.py
# Flask stuff:
instance/
.webassets-cache
# Scrapy stuff:
.scrapy
# Sphinx documentation
docs/_build/
# PyBuilder
target/