diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..9c351a0 --- /dev/null +++ b/.clang-format @@ -0,0 +1,5 @@ +--- +Language: Cpp +AllowShortLoopsOnASingleLine: true +AllowShortFunctionsOnASingleLine: All +AllowShortIfStatementsOnASingleLine: true diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml deleted file mode 100644 index a360982..0000000 --- a/.github/workflows/release.yml +++ /dev/null @@ -1,46 +0,0 @@ -name: Build and upload priwo to PyPI. - -on: - release: - types: - - published - -concurrency: - cancel-in-progress: true - group: ${{ github.workflow }}-${{ github.ref }} - -jobs: - # Build sdist and wheels. - dist: - name: Build sdist and wheels on ${{ matrix.os }} - runs-on: ${{ matrix.os }} - - strategy: - matrix: - os: [ubuntu-latest, windows-latest, macos-latest] - fail-fast: false - - steps: - - uses: actions/checkout@v3 - - name: Build sdist and wheels. - run: pipx run build - - name: Check metadata - run: pipx run twine check dist/* - - name: Upload wheels - uses: actions/upload-artifact@v3 - with: - path: dist/* - - # Upload the wheels and source distribution to PyPI. - publish: - needs: [dist] - runs-on: ubuntu-latest - if: github.event_name == 'release' && github.event.action == 'published' - steps: - - uses: actions/download-artifact@v3 - with: - name: artifact - path: dist - - uses: pypa/gh-action-pypi-publish@v1.6.4 - with: - password: ${{ secrets.pypi_password }} diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml deleted file mode 100644 index e213441..0000000 --- a/.github/workflows/test.yml +++ /dev/null @@ -1,43 +0,0 @@ -name: Testing priwo. - -on: - push: - -concurrency: - cancel-in-progress: true - group: ${{ github.workflow }}-${{ github.ref }} - -jobs: - build: - runs-on: ${{ matrix.os }} - strategy: - matrix: - os: [ubuntu-latest, macos-latest] - python-version: ["3.7", "3.8", "3.9", "3.10", "3.11"] - include: - - os: ubuntu-20.04 - python-version: "3.6" - - os: macos-latest - python-version: "3.6" - fail-fast: false - - steps: - # Setup the testing environment. - - uses: actions/checkout@v3 - - name: Set up Python ${{ matrix.python-version }}. - uses: actions/setup-python@v4 - with: - python-version: ${{ matrix.python-version }} - - # Sanity check: check Python version. - - name: Display Python version. - run: python -c "import sys; print(sys.version)" - - # Run tests for priwo. - - name: Install dependencies. - run: | - python -m pip install --upgrade pip - pip install ward numpy - pip install -e . - - name: Run tests for priwo. - run: ward diff --git a/.gitignore b/.gitignore index f62a92a..843b830 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,7 @@ +/build +.clangd +.clang_format + dist build .egg diff --git a/.readthedocs.yaml b/.readthedocs.yaml deleted file mode 100644 index bff8be8..0000000 --- a/.readthedocs.yaml +++ /dev/null @@ -1,13 +0,0 @@ -version: 2 - -build: - os: ubuntu-20.04 - tools: - python: "3.9" - -sphinx: - configuration: docs/source/conf.py - -python: - install: - - requirements: docs/requirements.txt diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..778a3df --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,74 @@ +# Set the minimum CMake version and policies for highest tested version +cmake_minimum_required(VERSION 3.15...3.27) + +# Set up the project and ensure there is a working C++ compiler +project(priwo LANGUAGES CXX) + +# Warn if the user invokes CMake directly +if (NOT SKBUILD) + message(WARNING "\ + This CMake file is meant to be executed using 'scikit-build-core'. + Running it directly will almost certainly not produce the desired + result. If you are a user trying to install this package, use the + command below, which will install all necessary build dependencies, + compile the package in an isolated environment, and then install it. + ===================================================================== + $ pip install . + ===================================================================== + If you are a software developer, and this is your own package, then + it is usually much more efficient to install the build dependencies + in your environment once and use the following command that avoids + a costly creation of a new virtual environment at every compilation: + ===================================================================== + $ pip install nanobind scikit-build-core[pyproject] + $ pip install --no-build-isolation -ve . + ===================================================================== + You may optionally add -Ceditable.rebuild=true to auto-rebuild when + the package is imported. Otherwise, you need to rerun the above + after editing C++ files.") +endif() + +# Try to import all Python components potentially needed by nanobind +find_package(Python 3.8 + REQUIRED COMPONENTS Interpreter Development.Module + OPTIONAL_COMPONENTS Development.SABIModule) + +# Import nanobind through CMake's find_package mechanism +find_package(nanobind CONFIG REQUIRED) + +# Add the module +nanobind_add_module( + _internals + STABLE_ABI + src/priwo/bits.h + src/priwo/bits.C + src/priwo/priwo.h + src/priwo/priwo.C + src/priwo/common.h + src/priwo/common.C + + src/priwo/hdr.h + src/priwo/hdr.C + src/priwo/fil.h + src/priwo/fil.C + src/priwo/tim.h + src/priwo/tim.C + + src/priwo/inf.h + src/priwo/inf.C + src/priwo/dat.h + src/priwo/dat.C + src/priwo/fft.h + src/priwo/fft.C + src/priwo/pfd.h + src/priwo/pfd.C + src/priwo/bestprof.h + src/priwo/bestprof.C + + + src/priwo/polycos.h + src/priwo/polycos.C +) + +# Install directive for scikit-build-core +install(TARGETS _internals LIBRARY DESTINATION priwo) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md deleted file mode 100644 index e69de29..0000000 diff --git a/Dockerfile b/Dockerfile deleted file mode 100644 index b08090e..0000000 --- a/Dockerfile +++ /dev/null @@ -1,25 +0,0 @@ -FROM python:3.8 - -RUN apt update && \ - apt install -y build-essential git vim && \ - rm -rf /var/lib/apt/lists/* -RUN pip install ipython - -ENV PRIWO_PATH=/software/priwo -RUN mkdir -p ${PRIWO_PATH} -WORKDIR ${PRIWO_PATH} -COPY . ${PRIWO_PATH} -RUN python -m pip install -e . - -RUN echo "\"\e[B\":history-search-forward" >> ~/.inputrc && \ - echo "\"\e[A\":history-search-backward" >> ~/.inputrc - -RUN echo "filetype plugin indent on" >> ~/.vimrc && \ - echo "set tabstop=4" >> ~/.vimrc && \ - echo "set shiftwidth=4" >> ~/.vimrc && \ - echo "set expandtab" >> ~/.vimrc && \ - echo "set pastetoggle=" >> ~/.vimrc && \ - echo "set hlsearch" >> ~/.vimrc && \ - echo "syntax on" >> ~/.vimrc - -ENTRYPOINT [ "/bin/bash" ] diff --git a/LICENSE b/LICENSE index 4da3094..e69de29 100644 --- a/LICENSE +++ b/LICENSE @@ -1,21 +0,0 @@ -MIT License - -Copyright (c) 2021-2023 Ujjwal Panda - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. diff --git a/MANIFEST.in b/MANIFEST.in deleted file mode 100644 index 80ecce8..0000000 --- a/MANIFEST.in +++ /dev/null @@ -1,4 +0,0 @@ -graft src -include *.md -include LICENSE -include setup.cfg diff --git a/README.md b/README.md index 933e6cc..e69de29 100644 --- a/README.md +++ b/README.md @@ -1,109 +0,0 @@ -
-
-priwo: I/O for common pulsar data formats. -
-
- -![License][license-badge] -![Version][version-badge] -![Python versions][pyversions-badge] - -![Tests][tests-badge] -[![Documentation Status][docs-badge]][docs] -[![Interrogate][interrogate-badge]][interrogate] - -![Stars][stars-badge] -![Downloads][dm-badge] -[![Issues][issues-badge]][issues] - -[![Gitmoji][gitmoji-badge]][gitmoji] -[![Code style: black][black-badge]][black] - -
- -

What is this?

- -[**`priwo`**][priwo] is a library that allows you to read in and write out -pulsar data from the following data formats: - -- [**`SIGPROC`**][sigproc] headers, -- [**`PRESTO`**][presto] FFT (`*.fft`) files, -- [**`PRESTO`**][presto] infodata (`*.inf`) files, -- [**`PRESTO`**][presto] folded data (`*.pfd`) files, -- [**`PRESTO`**][presto] time series (`*.dat`) files, -- [**`SIGPROC`**][sigproc] filterbank (`*.fil`) files, -- [**`SIGPROC`**][sigproc] time series (`*.tim`) files, -- [**`PRESTO`**][presto] best pulse profile (`*.bestprof`) files, -- [**`TEMPO`**][tempo] polynomial ephemerides (`*.polycos`) files. - -`priwo`'s API is deliberately _low-level_: each function in `priwo` deals with a -single file format and takes/returns a Python dictionary. This allows users to -design arbitrary high-level APIs on top of `priwo`'s functionality. This is -unlike most other contemporary libraries, such as [**`your`**][your]. `your` (to -which this library has been frequently compared to) provides a high-level API -for reading in pulsar data, while also providing modules to help process and -analyze it. This is makes the number of dependencies it uses is a bit high (as -of 05/10/22, that is a total of 9 dependencies). On the other hand, `priwo` has -just a single dependency: [**`pabo`**][pabo][^1]. This makes -it an ideal choice to drop into your projects, without worrying about -[**dependency hell**][dependency_hell]. - -`priwo` is well-tested (via [**`ward`**][ward]) and actively maintained. No -major changes to the API are expected before `v0.1.0`. Support for many more -formats, such as [`PSRFITS`][psrfits], is on the way. If you would like to -contribute, have a look at [`CONTRIBUTING.md`](CONTRIBUTING.md), and get in -touch! If you find a bug, feel free to open an [issue][issues]. If you would -like to suggest support for any data format(s) I have missed, suggest a feature, -or just chat, feel free to jump into the [discussions][discussions]. - -

Installing

- -Installing `priwo` is as easy as: - -```bash -pip install priwo -``` - -
- -[^1]: - [**`pabo`**][pabo] is a package I made to make parsing binary data easier, - and it has just two dependencies : [**`attrs`**][attrs] and - [**`numpy`**][numpy]. - -
- -[numpy]: https://numpy.org -[attrs]: https://www.attrs.org -[gitmoji]: https://gitmoji.dev -[black]: https://github.com/psf/black -[just]: https://github.com/casey/just -[tempo]: https://tempo.sourceforge.net -[sigproc]: http://sigproc.sourceforge.net -[pabo]: https://github.com/astrogewgaw/pabo -[ward]: https://github.com/darrenburns/ward -[priwo]: https://github.com/astrogewgaw/priwo -[docs]: https://priwo.readthedocs.io/en/latest -[presto]: https://github.com/scottransom/presto -[your]: https://github.com/thepetabyteproject/your -[issues]: https://github.com/astrogewgaw/priwo/issues -[interrogate]: https://github.com/econchick/interrogate -[discussions]: https://github.com/astrogewgaw/priwo/discussions -[dependency_hell]: https://en.wikipedia.org/wiki/Dependency_hell -[psrfits]: https://www.atnf.csiro.au/research/pulsar/psrfits_definition/Psrfits.html -[dm-badge]: https://img.shields.io/pypi/dm/priwo?style=for-the-badge -[version-badge]: https://img.shields.io/pypi/v/priwo?style=for-the-badge -[wheel-badge]: https://img.shields.io/pypi/wheel/priwo?style=for-the-badge -[forks-badge]: https://img.shields.io/github/forks/astrogewgaw/priwo?style=for-the-badge -[stars-badge]: https://img.shields.io/github/stars/astrogewgaw/priwo?style=for-the-badge -[pyversions-badge]: https://img.shields.io/pypi/pyversions/priwo.svg?style=for-the-badge -[issues-badge]: https://img.shields.io/github/issues/astrogewgaw/priwo?style=for-the-badge -[license-badge]: https://img.shields.io/github/license/astrogewgaw/priwo?style=for-the-badge -[black-badge]: https://img.shields.io/badge/code%20style-black-000000.svg?style=for-the-badge -[docs-badge]: https://readthedocs.org/projects/priwo/badge/?version=latest&style=for-the-badge -[gitmoji-badge]: https://img.shields.io/badge/gitmoji-%20😜%20😍-FFDD67.svg?style=for-the-badge -[interrogate-badge]: https://raw.githubusercontent.com/astrogewgaw/priwo/main/assets/docs_cov.svg -[tests-badge]: https://img.shields.io/github/actions/workflow/status/astrogewgaw/priwo/test.yml?branch=dev&style=for-the-badge diff --git a/assets/docs_cov.svg b/assets/docs_cov.svg deleted file mode 100644 index 4143f82..0000000 --- a/assets/docs_cov.svg +++ /dev/null @@ -1,28 +0,0 @@ - - INTERROGATE: 100% - - - - - - INTERROGATE - 100% - - - - - - - - - - - - - - - - - - - diff --git a/docs/Makefile b/docs/Makefile deleted file mode 100644 index 91f74ba..0000000 --- a/docs/Makefile +++ /dev/null @@ -1,12 +0,0 @@ -SPHINXOPTS ?= -SPHINXBUILD ?= sphinx-build -SOURCEDIR = source -BUILDDIR = build - -help: - @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) - -.PHONY: help Makefile - -%: Makefile - @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/requirements.txt b/docs/requirements.txt deleted file mode 100644 index 92694c9..0000000 --- a/docs/requirements.txt +++ /dev/null @@ -1,5 +0,0 @@ -furo -myst_parser -sphinx_sitemap -sphinx_copybutton -sphinx_autodoc_typehints diff --git a/docs/source/conf.py b/docs/source/conf.py deleted file mode 100644 index 5f500e2..0000000 --- a/docs/source/conf.py +++ /dev/null @@ -1,20 +0,0 @@ -project = "priwo" -html_theme = "furo" -author = "Ujjwal Panda" -exclude_patterns = [""] -html_static_path = ["static"] -templates_path = ["templates"] -copyright = "2021-2022, Ujjwal Panda" -html_baseurl = "https://priwo.readthedocs.io" - - -extensions = extensions = [ - "myst_parser", - "sphinx_sitemap", - "sphinx_copybutton", - "sphinx.ext.autodoc", - "sphinx.ext.viewcode", - "sphinx.ext.napoleon", - "sphinx_autodoc_typehints", - "sphinx.ext.autosectionlabel", -] diff --git a/docs/source/index.md b/docs/source/index.md deleted file mode 100644 index 16805a4..0000000 --- a/docs/source/index.md +++ /dev/null @@ -1,10 +0,0 @@ -```{toctree} -:glob: -:maxdepth: 2 -:caption: "Table of Contents:" -``` - -```{include} ../../README.md -:relative-docs: docs/ -:relative-images: -``` diff --git a/justfile b/justfile index 4cebc1b..03593ad 100644 --- a/justfile +++ b/justfile @@ -1,4 +1,5 @@ -pkg := "priwo" +pkg := "priwo" +desc := "I/O for common pulsar/FRB ddata formats." alias l := loc alias t := test @@ -16,7 +17,7 @@ default: from rich.console import Console console = Console() - + grid = Table.grid(expand=True, padding=(0, 2, 0, 2)) grid.add_column(justify="left", style="bold") grid.add_column(justify="right", style="italic") @@ -33,7 +34,7 @@ default: grid, padding=2, expand=False, - title="[b]{{pkg}}[/b]: [i]I/O for common pulsar data formats.[/i]", + title="[b]{{pkg}}[/b]: [i]{{desc}}[/i]", ) ) @@ -47,7 +48,6 @@ default: echo "Cleaning..." rm -rf tmp rm -rf dist - rm -rf build rm -rf .eggs rm -rf .coverage rm -rf .mypy_cache diff --git a/pyproject.toml b/pyproject.toml index ecf145b..bdccbfe 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,38 +1,37 @@ [build-system] -requires = [ - "wheel", - "setuptools>=45", - "setuptools_scm[toml]>=6.0", -] +requires = ["scikit-build-core >=0.4.3", "nanobind >=1.3.2"] +build-backend = "scikit_build_core.build" -[tool.setuptools_scm] -write_to = "src/priwo/_version.py" +[project] +name = "priwo" +version = "0.1.0" +readme = "README.md" +dependencies = ["numpy"] +requires-python = ">=3.8" +license = { file = "LICENSE" } +description = "I/O for common pulsar/FRB data formats." +authors = [{ name = "Ujjwal Panda", email = "ujjwapanda97@gmail.com" }] -[tool.interrogate] -verbose = 0 -color = true -quiet = false -fail-under = 85 -badge-format = "svg" -whitelist-regex = [] -ignore-magic = false -ignore-module = false -ignore-private = false -ignore-setters = false -ignore-init-method = true -omit-covered-files = false -ignore-init-module = false -ignore-semiprivate = false -ignore-nested-classes = true -badge-style = "for-the-badge" -ignore-nested-functions = false -ignore-property-decorators = false -generate-badge = "assets/docs_cov.svg" -exclude = [ - "docs", - "tests", - "build", - "assets", - "setup.py", - "src/priwo/_version.py", +classifiers = [ + "Operating System :: Unix", + "Operating System :: MacOS :: MacOS X", + "Operating System :: Microsoft :: Windows", + "Programming Language :: Python :: 3.8", + "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + "License :: OSI Approved :: MIT License", + "Topic :: Scientific/Engineering :: Astronomy", ] + +[project.urls] +Source = "https://github.com/astrogewgaw/priwo" +Homepage = "https://github.com/astrogewgaw/priwo" +Issues = "https://github.com/astrogewgaw/priwo/issues" +Documentation = "https://priwo.readthedocs.io/en/latest" + +[tool.scikit-build] +wheel.py-api = "cp312" # Build stable ABI wheels for CPython 3.12+ +minimum-version = "0.4" # Protect the configuration against future changes in scikit-build-core +build-dir = "build/{wheel_tag}" # Setuptools-style build caching in a local directory diff --git a/setup.cfg b/setup.cfg deleted file mode 100644 index 2344cd4..0000000 --- a/setup.cfg +++ /dev/null @@ -1,41 +0,0 @@ -[metadata] -name = priwo -license = MIT License -license_file = LICENSE -author = "Ujjwal Panda" -author_email = "ujjwalpanda97@gmail.com" -long_description_content_type = text/markdown -description = "I/O for common pulsar data formats." -long_description = file: README.md, CONTRIBUTING.md - -project_urls = - Home Page = https://github.com/astrogewgaw/priwo - Source Code = https://github.com/astrogewgaw/priwo - Documentation = https://priwo.readthedocs.io/en/latest - Bug Reports = https://github.com/astrogewgaw/priwo/issues - -classifiers = - Operating System :: Unix - Operating System :: MacOS :: MacOS X - Operating System :: Microsoft :: Windows - Programming Language :: Python :: 3.6 - Programming Language :: Python :: 3.7 - Programming Language :: Python :: 3.8 - Programming Language :: Python :: 3.9 - Programming Language :: Python :: 3.10 - Programming Language :: Python :: 3.11 - License :: OSI Approved :: MIT License - Topic :: Scientific/Engineering :: Astronomy - -[options] -zip_safe = False -packages = find: -package_dir = =src -use_scm_version = True -python_requires = >=3.6 -include_package_data = True -install_requires = pabo>=0.1.3 -setup_requires = setuptools_scm[toml]>=6.0 - -[options.packages.find] -where=src diff --git a/setup.py b/setup.py deleted file mode 100644 index 7f1a176..0000000 --- a/setup.py +++ /dev/null @@ -1,4 +0,0 @@ -from setuptools import setup - -if __name__ == "__main__": - setup() diff --git a/src/priwo/__init__.py b/src/priwo/__init__.py index 6a65c40..68da843 100644 --- a/src/priwo/__init__.py +++ b/src/priwo/__init__.py @@ -1,28 +1,41 @@ -""" -I/O for common pulsar data formats. -""" - -from priwo.presto import * -from priwo.timing import * -from priwo.sigproc import * +from priwo._internals import ( + readhdr, + readtim, + readfil, + readinf, + readdat, + readfft, + readpfd, + readpolycos, + readbestprof, + writehdr, + writetim, + writefil, + writeinf, + writedat, + writefft, + writepfd, + writepolycos, + writebestprof, +) __all__ = [ + "readhdr", + "readtim", + "readfil", "readinf", "readdat", "readfft", - "readbpf", "readpfd", - "readhdr", - "readtim", - "readfil", "readpolycos", + "readbestprof", + "writehdr", + "writetim", + "writefil", "writeinf", "writedat", "writefft", - "writebpf", "writepfd", - "writehdr", - "writetim", - "writefil", "writepolycos", + "writebestprof", ] diff --git a/src/priwo/_bpf.py b/src/priwo/_bpf.py new file mode 100644 index 0000000..2d72afd --- /dev/null +++ b/src/priwo/_bpf.py @@ -0,0 +1,100 @@ +import re +import numpy as np + + +def _bpf(fn): + with open(fn, "r") as f: + lines = f.readlines() + + fields = {} + for line in lines: + if line.startswith("#"): + cleaned = line.strip() + cleaned = cleaned.replace("#", "") + if len(cleaned) > 0: + key, string = list(map(str.strip, re.split(r"[=<]", cleaned))) + fields[key] = string + + meta = {} + for key, string in fields.items(): + varname, vartype = { + "Input file": ("filenm", str), + "Candidate": ("candnm", str), + "Telescope": ("telescope", str), + "Epoch_topo": ("tepoch", float), + "Epoch_bary": ("bepoch", float), + "T_sample": ("dt", float), + "Data Folded": ("N", float), + "Data Avg": ("data_avg", float), + "Data StdDev": ("data_std", float), + "Profile Bins": ("proflen", int), + "Profile Avg": ("prof_avg", float), + "Profile StdDev": ("prof_std", float), + "Reduced chi-sqr": ("redchi", float), + "Prob(Noise)": (None, None), + "Best DM": ("bestdm", float), + "P_topo (ms)": ("p_topo", float), + "P'_topo (s/s)": ("pd_topo", float), + "P''_topo (s/s^2)": ("pdd_topo", float), + "P_bary (ms)": ("p_bary", float), + "P'_bary (s/s)": ("pd_bary", float), + "P''_bary (s/s^2)": ("pdd_bary", float), + "P_orb (s)": ("orb_p", float), + "asin(i)/c (s)": ("orb_x", float), + "eccentricity": ("orb_e", float), + "w (rad)": ("orb_w", float), + "T_peri": ("orb_t", float), + }[key] + + if key != "Prob(Noise)": + if key in [ + "P_topo (ms)", + "P'_topo (s/s)", + "P''_topo (s/s^2)", + "P_bary (ms)", + "P'_bary (s/s)", + "P''_bary (s/s^2)", + ]: + meta[varname], meta["".join([varname, "_err"])] = None, None + if string != "N/A": + values = string.split("+/-") + value, error = list(map(vartype, values)) + meta[varname], meta["".join([varname, "_err"])] = value, error + else: + meta[varname] = None + if string != "N/A": + meta[varname] = vartype(string) + else: + regex = re.compile(r"[+-]?(?:0|[1-9]\d*)(?:\.\d*)?(?:[eE][+\-]?\d+)?") + if (matches := re.findall(regex, string)) is not None: + meta["chi_prob"], meta["chi_sig"] = list(map(float, matches)) + + meta["topo"] = { + "p": meta.pop("p_topo"), + "pd": meta.pop("pd_topo"), + "pdd": meta.pop("pdd_topo"), + "p_err": meta.pop("p_topo_err"), + "pd_err": meta.pop("pd_topo_err"), + "pdd_err": meta.pop("pdd_topo_err"), + } + + meta["bary"] = { + "p": meta.pop("p_bary"), + "pd": meta.pop("pd_bary"), + "pdd": meta.pop("pdd_bary"), + "p_err": meta.pop("p_bary_err"), + "pd_err": meta.pop("pd_bary_err"), + "pdd_err": meta.pop("pdd_bary_err"), + } + + meta["orb"] = { + "p": meta.pop("orb_p"), + "x": meta.pop("orb_x"), + "e": meta.pop("orb_e"), + "w": meta.pop("orb_w"), + "t": meta.pop("orb_t"), + } + + data = [float(_) for _ in re.findall(r"^\s+\d+\s+(.+)$", "".join(lines), re.M)] + data = np.asarray(data) + return meta, data diff --git a/src/priwo/bestprof.C b/src/priwo/bestprof.C new file mode 100644 index 0000000..bb8c400 --- /dev/null +++ b/src/priwo/bestprof.C @@ -0,0 +1,141 @@ +#include "bestprof.h" +#include "common.h" +#include + +nb::tuple readbestprof(std::string fn) { + nb::object _bpf = nb::module_::import_("priwo._bpf").attr("_bpf"); + nb::tuple output = nb::cast(_bpf(fn)); + return output; +} + +void writebestprof(nb::dict meta, + nb::ndarray> data, + std::string fn) { + FILE *f = chkfopen(fn.c_str(), "w"); + + const char *filenm = nb::cast(meta["filenm"]); + const char *candnm = nb::cast(meta["candnm"]); + const char *telescope = nb::cast(meta["telescope"]); + + fprintf(f, "# Input file = %-s\n", filenm); + fprintf(f, "# Candidate = %-s\n", candnm); + fprintf(f, "# Telescope = %-s\n", telescope); + + double tepoch; + if (meta["tepoch"].is_none()) { + tepoch = 0.0; + fprintf(f, "# Epoch_topo = N/A\n"); + } else { + tepoch = nb::cast(meta["tepoch"]); + fprintf(f, "# Epoch_topo = %-.12f\n", tepoch); + } + + double bepoch; + if (meta["bepoch"].is_none()) { + bepoch = 0.0; + fprintf(f, "# Epoch_bary = N/A\n"); + } else { + bepoch = nb::cast(meta["bepoch"]); + fprintf(f, "# Epoch_bary (MJD) = %-.12f\n", bepoch); + } + + double N = nb::cast(meta["N"]); + double dt = nb::cast(meta["dt"]); + int proflen = nb::cast(meta["proflen"]); + double data_avg = nb::cast(meta["data_avg"]); + double data_std = nb::cast(meta["data_std"]); + double prof_avg = nb::cast(meta["prof_avg"]); + double prof_std = nb::cast(meta["prof_std"]); + + fprintf(f, "# T_sample = %.6g\n", dt); + fprintf(f, "# Data Folded = %-.0f\n", N); + fprintf(f, "# Data Avg = %-17.15g\n", data_avg); + fprintf(f, "# Data StdDev = %-17.15g\n", data_std); + fprintf(f, "# Profile Bins = %d\n", proflen); + fprintf(f, "# Profile Avg = %-17.15g\n", prof_avg); + fprintf(f, "# Profile StdDev = %-17.15g\n", prof_std); + + char tmp[80]; + + double redchi = nb::cast(meta["redchi"]); + double chi_sig = nb::cast(meta["chi_sig"]); + + sprintf(tmp, "(~%.1f sigma)", chi_sig); + fprintf(f, "# Reduced chi-sqr = %.3f\n", redchi); + fprintf(f, "# Prob(Noise) < %.3g %s\n", 0.0, tmp); + + if (!meta["bestdm"].is_none()) { + double bestdm = nb::cast(meta["bestdm"]); + fprintf(f, "# Best DM = %.3f\n", bestdm); + } + + if (tepoch != 0.0) { + double p_topo = nb::cast(meta["topo"]["p"]); + double pd_topo = nb::cast(meta["topo"]["pd"]); + double pdd_topo = nb::cast(meta["topo"]["pdd"]); + double p_topo_err = nb::cast(meta["topo"]["p_err"]); + double pd_topo_err = nb::cast(meta["topo"]["pd_err"]); + double pdd_topo_err = nb::cast(meta["topo"]["pdd_err"]); + + fprintf(f, "# P_topo (ms) = %-17.15g +/- %-.3g\n", p_topo, + p_topo_err); + fprintf(f, "# P'_topo (s/s) = %-17.15g +/- %-.3g\n", pd_topo, + pd_topo_err); + fprintf(f, "# P''_topo (s/s^2) = %-17.15g +/- %-.3g\n", pdd_topo, + pdd_topo_err); + } else { + fprintf(f, "# P_topo (ms) = N/A\n"); + fprintf(f, "# P'_topo (s/s) = N/A\n"); + fprintf(f, "# P''_topo (s/s^2) = N/A\n"); + } + + if (bepoch != 0.0) { + double p_bary = nb::cast(meta["bary"]["p"]); + double pd_bary = nb::cast(meta["bary"]["pd"]); + double pdd_bary = nb::cast(meta["bary"]["pdd"]); + double p_bary_err = nb::cast(meta["bary"]["p_err"]); + double pd_bary_err = nb::cast(meta["bary"]["pd_err"]); + double pdd_bary_err = nb::cast(meta["bary"]["pdd_err"]); + + fprintf(f, "# P_bary (ms) = %-17.15g +/- %-.3g\n", p_bary, + p_bary_err); + fprintf(f, "# P'_bary (s/s) = %-17.15g +/- %-.3g\n", pd_bary, + pd_bary_err); + fprintf(f, "# P''_bary (s/s^2) = %-17.15g +/- %-.3g\n", pdd_bary, + pdd_bary_err); + } else { + fprintf(f, "# P_bary (ms) = N/A\n"); + fprintf(f, "# P'_bary (s/s) = N/A\n"); + fprintf(f, "# P''_bary (s/s^2) = N/A\n"); + } + + if (meta["orb"]["p"].is_none()) { + fprintf(f, "# P_orb (s) = N/A\n"); + fprintf(f, "# asin(i)/c (s) = N/A\n"); + fprintf(f, "# eccentricity = N/A\n"); + fprintf(f, "# w (rad) = N/A\n"); + fprintf(f, "# T_peri = N/A\n"); + } else { + double orb_p = nb::cast(meta["orb"]["p"]); + double orb_x = nb::cast(meta["orb"]["x"]); + double orb_e = nb::cast(meta["orb"]["e"]); + double orb_w = nb::cast(meta["orb"]["w"]); + double orb_t = nb::cast(meta["orb"]["t"]); + + fprintf(f, "# P_orb (s) = %-17.15g\n", orb_p); + fprintf(f, "# asin(i)/c (s) = %-17.15g\n", orb_x); + fprintf(f, "# eccentricity = %-17.15g\n", orb_e); + fprintf(f, "# w (deg) = %-17.15g\n", orb_w); + fprintf(f, "# T_peri = %-.12f\n", orb_t); + } + + fprintf(f, "######################################################\n"); + for (int ii = 0; ii < data.size(); ii++) + fprintf(f, "%4d %.7g\n", ii, data(ii)); + fclose(f); +} + +void init_bestprof(nb::module_ m) { + m.def("readbestprof", &readbestprof); + m.def("writebestprof", &writebestprof); +} diff --git a/src/priwo/bestprof.h b/src/priwo/bestprof.h new file mode 100644 index 0000000..9b267b8 --- /dev/null +++ b/src/priwo/bestprof.h @@ -0,0 +1,12 @@ +#ifndef BESTPROF +#define BESTPROF + +#include "common.h" + +void init_bestprof(nb::module_ m); +nb::tuple readbestprof(std::string fn); +void writebestprof(nb::dict meta, + nb::ndarray> data, + std::string fn); + +#endif diff --git a/src/priwo/bits.C b/src/priwo/bits.C new file mode 100644 index 0000000..863d026 --- /dev/null +++ b/src/priwo/bits.C @@ -0,0 +1,211 @@ +#include +#include +#include + +#ifdef USE_OPENMP +#include +#endif + +#include "bits.h" +#include "common.h" + +#define LO2BITS 3 +#define LO4BITS 15 +#define HI2BITS 192 +#define HI4BITS 240 +#define LOMED2BITS 12 +#define UPMED2BITS 48 + +template +void unpack_1bit(const uint8_t *inbuffer, uint8_t *outbuffer, size_t nbytes) { +#ifdef USE_OPENMP +#pragma omp parallel for if (parallel) +#endif + for (size_t ii = 0; ii < nbytes; ii++) { + for (size_t jj = 0; jj < 8; jj++) { + if constexpr (bigEndian) { + outbuffer[(ii << 3) + (7 - jj)] = (inbuffer[ii] >> jj) & 1; + } else { + outbuffer[(ii << 3) + jj] = (inbuffer[ii] >> jj) & 1; + } + } + } +} + +template +void unpack_2bit(const uint8_t *inbuffer, uint8_t *outbuffer, size_t nbytes) { +#ifdef USE_OPENMP +#pragma omp parallel for if (parallel) +#endif + for (size_t ii = 0; ii < nbytes; ii++) { + if constexpr (bigEndian) { + outbuffer[(ii << 2) + 3] = inbuffer[ii] & LO2BITS; + outbuffer[(ii << 2) + 2] = (inbuffer[ii] & LOMED2BITS) >> 2; + outbuffer[(ii << 2) + 1] = (inbuffer[ii] & UPMED2BITS) >> 4; + outbuffer[(ii << 2) + 0] = (inbuffer[ii] & HI2BITS) >> 6; + } else { + outbuffer[(ii << 2) + 0] = inbuffer[ii] & LO2BITS; + outbuffer[(ii << 2) + 1] = (inbuffer[ii] & LOMED2BITS) >> 2; + outbuffer[(ii << 2) + 2] = (inbuffer[ii] & UPMED2BITS) >> 4; + outbuffer[(ii << 2) + 3] = (inbuffer[ii] & HI2BITS) >> 6; + } + } +} + +template +void unpack_4bit(const uint8_t *inbuffer, uint8_t *outbuffer, size_t nbytes) { +#ifdef USE_OPENMP +#pragma omp parallel for if (parallel) +#endif + for (size_t ii = 0; ii < nbytes; ii++) { + if constexpr (bigEndian) { + outbuffer[(ii << 1) + 1] = inbuffer[ii] & LO4BITS; + outbuffer[(ii << 1) + 0] = (inbuffer[ii] & HI4BITS) >> 4; + } else { + outbuffer[(ii << 1) + 0] = inbuffer[ii] & LO4BITS; + outbuffer[(ii << 1) + 1] = (inbuffer[ii] & HI4BITS) >> 4; + } + } +} + +template +void pack_1bit(const uint8_t *inbuffer, uint8_t *outbuffer, size_t nbytes) { + size_t pos; +#ifdef USE_OPENMP +#pragma omp parallel for if (parallel) +#endif + for (size_t ii = 0; ii < nbytes / 8; ii++) { + pos = ii * 8; + if constexpr (bigEndian) { + outbuffer[ii] = (inbuffer[pos + 0] << 7) | (inbuffer[pos + 1] << 6) | + (inbuffer[pos + 2] << 5) | (inbuffer[pos + 3] << 4) | + (inbuffer[pos + 4] << 3) | (inbuffer[pos + 5] << 2) | + (inbuffer[pos + 6] << 1) | inbuffer[pos + 7]; + } else { + outbuffer[ii] = inbuffer[pos + 0] | (inbuffer[pos + 1] << 1) | + (inbuffer[pos + 2] << 2) | (inbuffer[pos + 3] << 3) | + (inbuffer[pos + 4] << 4) | (inbuffer[pos + 5] << 5) | + (inbuffer[pos + 6] << 6) | (inbuffer[pos + 7] << 7); + } + } +} + +template +void pack_2bit(const uint8_t *inbuffer, uint8_t *outbuffer, size_t nbytes) { + size_t pos; +#ifdef USE_OPENMP +#pragma omp parallel for if (parallel) +#endif + for (size_t ii = 0; ii < nbytes / 4; ii++) { + pos = ii * 4; + if constexpr (bigEndian) { + outbuffer[ii] = (inbuffer[pos + 0] << 6) | (inbuffer[pos + 1] << 4) | + (inbuffer[pos + 2] << 2) | inbuffer[pos + 3]; + } else { + outbuffer[ii] = inbuffer[pos + 0] | (inbuffer[pos + 1] << 2) | + (inbuffer[pos + 2] << 4) | (inbuffer[pos + 3] << 6); + } + } +} + +template +void pack_4bit(const uint8_t *inbuffer, uint8_t *outbuffer, size_t nbytes) { + size_t pos; +#ifdef USE_OPENMP +#pragma omp parallel for if (parallel) +#endif + for (size_t ii = 0; ii < nbytes / 2; ii++) { + pos = ii * 2; + if constexpr (bigEndian) { + outbuffer[ii] = (inbuffer[pos] << 4) | inbuffer[pos + 1]; + } else { + outbuffer[ii] = inbuffer[pos] | (inbuffer[pos + 1] << 4); + } + } +} + +using PackUnpackFunc = void (*)(const uint8_t *, uint8_t *, size_t); + +constexpr std::array, 2>, 3> + unpackDispatcher = { + {{{ + {unpack_1bit, unpack_1bit}, // little + {unpack_1bit, unpack_1bit} // big + }}, + {{ + {unpack_2bit, unpack_2bit}, // little + {unpack_2bit, unpack_2bit} // big + }}, + {{ + {unpack_4bit, unpack_4bit}, // little + {unpack_4bit, unpack_4bit} // big + }}}}; + +constexpr std::array, 2>, 3> + packDispatcher = { + {{{ + {pack_1bit, pack_1bit}, // little + {pack_1bit, pack_1bit} // big + }}, + {{ + {pack_2bit, pack_2bit}, // little + {pack_2bit, pack_2bit} // big + }}, + {{ + {pack_4bit, pack_4bit}, // little + {pack_4bit, pack_4bit} // big + }}}}; + +size_t get_bitorder_index(const std::string &bitorder) { + if (bitorder.empty() || (bitorder[0] != 'l' && bitorder[0] != 'b')) { + throw std::invalid_argument( + "Invalid bitorder. Must begin with 'l' or 'b'."); + } + return (bitorder[0] == 'b') ? 1 : 0; +} + +nb::ndarray +unpack(const nb::ndarray &x, size_t nbits, + const std::string &bitorder, bool parallel = false) { + if (nbits != 1 && nbits != 2 && nbits != 4) { + throw std::invalid_argument( + "Invalid number of bits. Supported values are 1, 2, and 4."); + } + size_t bitorder_idx = get_bitorder_index(bitorder); + size_t nbits_idx = nbits >> 1; + size_t nbytes = x.size(); + auto y = EmptyArr1D(nbytes * 8 / nbits); + + PackUnpackFunc unpackFunc = + unpackDispatcher[nbits_idx][bitorder_idx][parallel ? 1 : 0]; + unpackFunc(x.data(), y.data(), nbytes); + return y; +} + +nb::ndarray +pack(const nb::ndarray &x, size_t nbits, + const std::string &bitorder, bool parallel = false) { + if (nbits != 1 && nbits != 2 && nbits != 4) { + throw std::invalid_argument( + "Invalid number of bits. Supported values are 1, 2, and 4."); + } + size_t bitorder_idx = get_bitorder_index(bitorder); + size_t nbits_idx = nbits >> 1; + size_t nbytes = x.size(); + auto y = EmptyArr1D(nbytes * nbits / 8); + + PackUnpackFunc packFunc = + packDispatcher[nbits_idx][bitorder_idx][parallel ? 1 : 0]; + packFunc(x.data(), y.data(), nbytes); + return y; +} + +void init_bits(nb::module_ m) { + m.def("unpack", &unpack, + "Unpack 1, 2 and 4-bit data from an 8-bit numpy array", nb::arg("x"), + nb::arg("nbits"), nb::arg("bitorder") = "big", + nb::arg("parallel") = false); + m.def("pack", &pack, "Pack 1, 2 and 4-bit data into an 8-bit numpy array", + nb::arg("x"), nb::arg("nbits"), nb::arg("bitorder") = "big", + nb::arg("parallel") = false); +} diff --git a/src/priwo/bits.h b/src/priwo/bits.h new file mode 100644 index 0000000..f78f38d --- /dev/null +++ b/src/priwo/bits.h @@ -0,0 +1,16 @@ +#ifndef BITS +#define BITS + +#include "common.h" + +void init_bits(nb::module_ m); + +nb::ndarray +unpack(const nb::ndarray &x, size_t nbits, + const std::string &bitorder, bool parallel); + +nb::ndarray +pack(const nb::ndarray &x, size_t nbits, + const std::string &bitorder, bool parallel); + +#endif diff --git a/src/priwo/common.C b/src/priwo/common.C new file mode 100644 index 0000000..8215876 --- /dev/null +++ b/src/priwo/common.C @@ -0,0 +1,169 @@ +#include "common.h" + +#ifndef SWAP +#define SWAP(a, b) \ + tmpswap = (a); \ + (a) = (b); \ + (b) = tmpswap; +#endif + +static unsigned char tmpswap; + +char *rmtrail(char *str) { + int i; + if (str && 0 != (i = strlen(str))) { + while (--i >= 0) { + if (!isspace(str[i])) break; + } + str[++i] = '\0'; + } + return str; +} + +char *rmlead(char *str) { + char *obuf; + if (str) { + for (obuf = str; *obuf && isspace(*obuf); ++obuf); + if (str != obuf) memmove(str, obuf, strlen(obuf) + 1); + } + return str; +} + +char *rmspace(char *str) { return rmlead(rmtrail(str)); } + +int swapint(int var) { + unsigned char *buffer; + int *iptr; + + buffer = (unsigned char *)(&var); + SWAP(buffer[0], buffer[3]); + SWAP(buffer[1], buffer[2]); + iptr = (int *)buffer; + return *iptr; +} + +unsigned int swapuint(unsigned int var) { + unsigned char *buffer; + unsigned int *iptr; + + buffer = (unsigned char *)(&var); + SWAP(buffer[0], buffer[3]); + SWAP(buffer[1], buffer[2]); + iptr = (unsigned int *)buffer; + return *iptr; +} + +short swapshort(short var) { + unsigned char *buffer; + short *sptr; + + buffer = (unsigned char *)(&var); + SWAP(buffer[0], buffer[1]); + sptr = (short *)buffer; + return *sptr; +} + +unsigned short swapushort(unsigned short var) { + unsigned char *buffer; + unsigned short *sptr; + + buffer = (unsigned char *)(&var); + SWAP(buffer[0], buffer[1]); + sptr = (unsigned short *)buffer; + return *sptr; +} + +float swapfloat(float var) { + unsigned char *buffer; + float *fptr; + + buffer = (unsigned char *)(&var); + SWAP(buffer[0], buffer[3]); + SWAP(buffer[1], buffer[2]); + fptr = (float *)buffer; + return *fptr; +} + +double swapdouble(double var) { + unsigned char *buffer; + double *dptr; + buffer = (unsigned char *)(&var); + SWAP(buffer[0], buffer[7]); + SWAP(buffer[1], buffer[6]); + SWAP(buffer[2], buffer[5]); + SWAP(buffer[3], buffer[4]); + dptr = (double *)buffer; + return *dptr; +} + +long long swaplonglong(long long var) { + unsigned char *buffer; + long long *llptr; + + buffer = (unsigned char *)(&var); + SWAP(buffer[0], buffer[7]); + SWAP(buffer[1], buffer[6]); + SWAP(buffer[2], buffer[5]); + SWAP(buffer[3], buffer[4]); + llptr = (long long *)buffer; + return *llptr; +} + +long double swaplongdouble(long double var) { + unsigned char *buffer; + long double *ldptr; + + buffer = (unsigned char *)(&var); + SWAP(buffer[0], buffer[11]); + SWAP(buffer[1], buffer[10]); + SWAP(buffer[2], buffer[9]); + SWAP(buffer[3], buffer[8]); + SWAP(buffer[4], buffer[7]); + SWAP(buffer[5], buffer[6]); + ldptr = (long double *)buffer; + return *ldptr; +} + +int parseint(FILE *infile, int byteswap) { + int itmp; + chkfread(&itmp, sizeof(int), 1, infile); + if (byteswap) itmp = swapint(itmp); + return itmp; +} + +float parsefloat(FILE *infile, int byteswap) { + float ftmp; + chkfread(&ftmp, sizeof(float), 1, infile); + if (byteswap) ftmp = swapfloat(ftmp); + return ftmp; +} + +double parsedouble(FILE *infile, int byteswap) { + double dtmp; + chkfread(&dtmp, sizeof(double), 1, infile); + if (byteswap) dtmp = swapdouble(dtmp); + return dtmp; +} + +FILE *chkfopen(const char *path, const char *mode) { + FILE *file; + if ((file = fopen(path, mode)) == NULL) + throw std::runtime_error("Cannot open file! Exiting..."); + return (file); +} + +size_t chkfread(void *data, size_t type, size_t number, FILE *stream) { + size_t num; + num = fread(data, type, number, stream); + if (num != number && ferror(stream)) + throw std::runtime_error("Cannot read from file! Exiting..."); + return num; +} + +size_t chkfwrite(void *data, size_t type, size_t number, FILE *stream) { + size_t num; + num = fwrite(data, type, number, stream); + if (num != number && ferror(stream)) + throw std::runtime_error("Cannot write to file! Exiting..."); + return num; +} diff --git a/src/priwo/common.h b/src/priwo/common.h new file mode 100644 index 0000000..827252c --- /dev/null +++ b/src/priwo/common.h @@ -0,0 +1,99 @@ +#ifndef COMMON +#define COMMON + +#include +#include +#include +#include +#include +#include +#include + +#ifdef __linux__ +#include +#endif + +#include +#include +#include +#include + +namespace nb = nanobind; + +char *rmlead(char *str); +char *rmtrail(char *str); +char *rmspace(char *str); + +int swapint(int var); +short swapshort(short var); +float swapfloat(float var); +double swapdouble(double vari); +long long swaplonglong(long long var); +unsigned int swapuint(unsigned int var); +long double swaplongdouble(long double var); +unsigned short swapushort(unsigned short var); + +int parseint(FILE *infile, int byteswap); +float parsefloat(FILE *infile, int byteswap); +double parsedouble(FILE *infile, int byteswap); + +FILE *chkfopen(const char *path, const char *mode); +size_t chkfread(void *data, size_t type, size_t number, FILE *stream); +size_t chkfwrite(void *data, size_t type, size_t number, FILE *stream); + +template +nb::ndarray EmptyArr1D(size_t N) { + T *data = new T[N]; + +#ifdef __linux__ + // Use madvise with MADV_HUGEPAGE to optimize memory usage on Linux + const size_t hugepage_threshold = 1u << 22u; // 4MB threshold + const size_t page_size = 4096u; + + if (N * sizeof(T) >= hugepage_threshold) { + uintptr_t data_addr = reinterpret_cast(data); + size_t offset = page_size - (data_addr % page_size); + size_t length = N * sizeof(T) - offset; + + // Intentionally not checking for errors, following NumPy's approach + madvise(reinterpret_cast(data_addr + offset), length, + MADV_HUGEPAGE); + } +#endif + + size_t shape[1] = {N}; + + // Create and return the ndarray with the given shape and ownership capsule + nb::capsule owner(data, [](void *p) noexcept { delete[] (T *)p; }); + return nb::ndarray(data, 1, shape, owner); +} + +template +nb::ndarray EmptyArr2D(size_t N, size_t M) { + int totalsize = N * M; + T *data = new T[totalsize]; + +#ifdef __linux__ + // Use madvise with MADV_HUGEPAGE to optimize memory usage on Linux + const size_t hugepage_threshold = 1u << 22u; // 4MB threshold + const size_t page_size = 4096u; + + if (totalsize * sizeof(T) >= hugepage_threshold) { + uintptr_t data_addr = reinterpret_cast(data); + size_t offset = page_size - (data_addr % page_size); + size_t length = totalsize * sizeof(T) - offset; + + // Intentionally not checking for errors, following NumPy's approach + madvise(reinterpret_cast(data_addr + offset), length, + MADV_HUGEPAGE); + } +#endif + + size_t shape[2] = {N, M}; + + // Create and return the ndarray with the given shape and ownership capsule + nb::capsule owner(data, [](void *p) noexcept { delete[] (T *)p; }); + return nb::ndarray(data, 2, shape, owner); +} + +#endif diff --git a/src/priwo/dat.C b/src/priwo/dat.C new file mode 100644 index 0000000..4fc431a --- /dev/null +++ b/src/priwo/dat.C @@ -0,0 +1,25 @@ +#include "dat.h" + +namespace fs = std::filesystem; +using namespace nanobind::literals; + +std::tuple>> +readdat(std::string fn) { + nb::dict meta = readinf(fs::path(fn).replace_extension("inf")); + nb::object np = nb::module_::import_("numpy"); + return std::make_tuple( + meta, nb::cast>>( + np.attr("fromfile")(fn, "dtype"_a = np.attr("float32")))); +} + +void writedat(nb::dict meta, nb::ndarray> data, + std::string fn) { + writeinf(meta, fs::path(fn).replace_extension("inf")); + nb::object np = nb::module_::import_("numpy"); + np.attr("asarray")(data).attr("tofile")(fn); +} + +void init_dat(nb::module_ m) { + m.def("readdat", &readdat); + m.def("writedat", &writedat); +} diff --git a/src/priwo/dat.h b/src/priwo/dat.h new file mode 100644 index 0000000..1894a9b --- /dev/null +++ b/src/priwo/dat.h @@ -0,0 +1,15 @@ +#ifndef DAT +#define DAT + +#include + +#include "common.h" +#include "inf.h" + +void init_dat(nb::module_ m); +std::tuple>> +readdat(std::string fn); +void writedat(nb::dict meta, nb::ndarray> data, + std::string fn); + +#endif diff --git a/src/priwo/fft.C b/src/priwo/fft.C new file mode 100644 index 0000000..079b74c --- /dev/null +++ b/src/priwo/fft.C @@ -0,0 +1,25 @@ +#include "fft.h" + +namespace fs = std::filesystem; +using namespace nanobind::literals; + +std::tuple>> +readfft(std::string fn) { + nb::dict meta = readinf(fs::path(fn).replace_extension("inf")); + nb::object np = nb::module_::import_("numpy"); + return std::make_tuple( + meta, nb::cast>>( + np.attr("fromfile")(fn, "dtype"_a = np.attr("complex64")))); +} + +void writefft(nb::dict meta, nb::ndarray> data, + std::string fn) { + writeinf(meta, fs::path(fn).replace_extension("inf")); + nb::object np = nb::module_::import_("numpy"); + np.attr("asarray")(data).attr("tofile")(fn); +} + +void init_fft(nb::module_ m) { + m.def("readfft", &readfft); + m.def("writefft", &writefft); +} diff --git a/src/priwo/fft.h b/src/priwo/fft.h new file mode 100644 index 0000000..820db5e --- /dev/null +++ b/src/priwo/fft.h @@ -0,0 +1,15 @@ +#ifndef FFT +#define FFT + +#include + +#include "common.h" +#include "inf.h" + +void init_fft(nb::module_ m); +std::tuple>> +readfft(std::string fn); +void writefft(nb::dict meta, nb::ndarray> data, + std::string fn); + +#endif diff --git a/src/priwo/fil.C b/src/priwo/fil.C new file mode 100644 index 0000000..9ce8601 --- /dev/null +++ b/src/priwo/fil.C @@ -0,0 +1,68 @@ +#include "fil.h" +#include "hdr.h" + +using namespace nanobind::literals; + +std::tuple>> +readfil(std::string fn) { + nb::dict meta = readhdr(fn); + int nf = nb::cast(meta["nchans"]); + int nbits = nb::cast(meta["nbits"]); + int nskip = nb::cast(meta["hdrlen"]); + + const char *dtype; + switch (nbits) { + case 1: + dtype = "uint8"; + break; + case 2: + dtype = "uint8"; + break; + case 4: + dtype = "uint8"; + break; + case 8: + dtype = "uint8"; + break; + case 16: + dtype = "uint16"; + break; + case 32: + dtype = "float32"; + break; + } + + nb::object np = nb::module_::import_("numpy"); + nb::object dobj = + np.attr("fromfile")(fn, "dtype"_a = np.attr(dtype), "offset"_a = nskip); + if (nbits == 1 || nbits == 2 || nbits == 4) { + nb::object priwo = nb::module_::import_("priwo"); + dobj = priwo.attr("_internals") + .attr("unpack")(dobj, "nbits"_a = nbits, "bitorder"_a = "little", + "parallel"_a = true); + } + nb::ndarray> data = + nb::cast>>(np.attr("transpose")( + np.attr("reshape")(dobj, "newshape"_a = nb::make_tuple(-1, nf)))); + return std::make_tuple(meta, data); +} + +void writefil(nb::dict meta, nb::ndarray> data, + std::string fn) { + writehdr(meta, fn); + nb::object np = nb::module_::import_("numpy"); + nb::object dobj = np.attr("asarray")(data); + int nbits = nb::cast(meta["nbits"]); + if (nbits == 1 || nbits == 2 || nbits == 4) { + nb::object priwo = nb::module_::import_("priwo"); + dobj = priwo.attr("_internals") + .attr("pack")(dobj, "nbits"_a = nbits, "bitorder"_a = "little", + "parallel"_a = true); + } + dobj.attr("T").attr("tofile")(fn); +} + +void init_fil(nb::module_ m) { + m.def("readfil", &readfil); + m.def("writefil", &writefil); +} diff --git a/src/priwo/fil.h b/src/priwo/fil.h new file mode 100644 index 0000000..eb6f2a1 --- /dev/null +++ b/src/priwo/fil.h @@ -0,0 +1,12 @@ +#ifndef FIL +#define FIL + +#include "common.h" + +void init_fil(nb::module_ m); +std::tuple>> +readfil(std::string fn); +void writefil(nb::dict meta, nb::ndarray> data, + std::string fn); + +#endif diff --git a/src/priwo/hdr.C b/src/priwo/hdr.C new file mode 100644 index 0000000..d6798e7 --- /dev/null +++ b/src/priwo/hdr.C @@ -0,0 +1,374 @@ +#include "hdr.h" + +static void _getstr(FILE *f, int *nbytes, char string[]) { + int nchar; + chkfread(&nchar, sizeof(int), 1, f); + *nbytes = sizeof(int); + if (feof(f)) throw std::runtime_error("Reached EOF! Exiting..."); + if (nchar > 80 || nchar < 1) return; + chkfread(string, nchar, 1, f); + string[nchar] = '\0'; + *nbytes += nchar; +} + +static void _sendstr(std::string string, FILE *f) { + int len; + len = strlen(string.c_str()); + chkfwrite(&len, sizeof(int), 1, f); + chkfwrite(string.data(), sizeof(char), len, f); +} + +static void _senddbl(std::string name, double double_precision, FILE *f) { + _sendstr(name, f); + chkfwrite(&double_precision, sizeof(double), 1, f); +} + +static void _sendint(std::string name, int integer, FILE *f) { + _sendstr(name, f); + chkfwrite(&integer, sizeof(int), 1, f); +} + +std::string _getteln(int telescope_id) { + char telescope[80]; + switch (telescope_id) { + case 0: + strcpy(telescope, "Fake"); + break; + case 1: + strcpy(telescope, "Arecibo"); + break; + case 2: + strcpy(telescope, "Ooty"); + break; + case 3: + strcpy(telescope, "Nancay"); + break; + case 4: + strcpy(telescope, "Parkes"); + break; + case 5: + strcpy(telescope, "Jodrell"); + break; + case 6: + strcpy(telescope, "GBT"); + break; + case 7: + strcpy(telescope, "GMRT"); + break; + case 8: + strcpy(telescope, "Effelsberg"); + break; + case 9: + strcpy(telescope, "ATA"); + break; + case 10: + strcpy(telescope, "SRT"); + break; + case 11: + strcpy(telescope, "LOFAR"); + break; + case 12: + strcpy(telescope, "VLA"); + break; + case 20: + strcpy(telescope, "CHIME"); + break; + case 21: + strcpy(telescope, "FAST"); + break; + case 30: + strcpy(telescope, "MWA"); + break; + case 64: + strcpy(telescope, "MeerKAT"); + break; + case 65: + strcpy(telescope, "KAT-7"); + break; + default: + strcpy(telescope, "Unknown"); + break; + } + return telescope; +} + +std::string _getmachn(int machine_id) { + char *backend, string[80]; + switch (machine_id) { + case 0: + strcpy(string, "FAKE"); + break; + case 1: + strcpy(string, "PSPM"); + break; + case 2: + strcpy(string, "WAPP"); + break; + case 3: + strcpy(string, "AOFTM"); + break; + case 4: + strcpy(string, "BPP"); + break; + case 5: + strcpy(string, "OOTY"); + break; + case 6: + strcpy(string, "SCAMP"); + break; + case 7: + strcpy(string, "SPIGOT"); + break; + case 11: + strcpy(string, "BG/P"); + break; + case 12: + strcpy(string, "PDEV"); + break; + case 20: + strcpy(string, "CHIME+PSR"); + break; + case 30: + strcpy(string, "MWA-VCS"); + break; + case 31: + strcpy(string, "MWAX-VCS"); + break; + case 32: + strcpy(string, "MWAX-RTB"); + break; + case 64: + strcpy(string, "KAT"); + break; + case 65: + strcpy(string, "KAT-DC2"); + break; + default: + strcpy(string, "Unknown"); + break; + } + backend = (char *)calloc(strlen(string) + 1, 1); + strcpy(backend, string); + return backend; +} + +nb::dict readhdr(std::string fn) { + nb::dict hdr; + + FILE *f = chkfopen(fn.c_str(), "rb"); + + char string[80], message[80]; + int ix, nbytes = 0, totalbytes; + int barycentric, pulsarcentric; + int expecting_rawdatafile = 0, expecting_source_name = 0; + + _getstr(f, &nbytes, string); + if (strcmp(string, "HEADER_START")) { + rewind(f); + throw std::runtime_error("SIGPROC header invalid! Exiting..."); + } + totalbytes = nbytes; + + int datatype; + char inpfile[80]; /* Input filename */ + char srcname[80]; /* Source name */ + long long N; /* Number of points (in time) in the file */ + double tstart; /* MJD start time */ + double tsamp; /* Sampling time in sec */ + double src_raj; /* Source RA (J2000) in hhmmss.ss */ + double src_dej; /* Source DEC (J2000) in ddmmss.ss */ + double az_start; /* Starting azimuth in deg */ + double za_start; /* Starting zenith angle in deg */ + double fch1; /* Highest channel frequency (MHz) */ + double foff; /* Channel stepsize (MHz) */ + double refdm; /* Reference dispersion measure (pc/cm^3) */ + int machine_id; /* Instrument ID (see backend_name() */ + int telescope_id; /* Telescope ID (see telescope_name() */ + int nchans; /* Number of filterbank channels */ + int nsamples; /* Number of filterbank samples */ + int nbits; /* Number of bits in the filterbank samples */ + int nifs; /* Number of IFs present */ + int nbeams; /* Number of beams in the observing system */ + int ibeam; /* Beam number used for this data */ + int sumifs; /* Whether the IFs are summed or not */ + int signedints; /* Whether the integer data is signed or not */ + + ibeam = 1; + sumifs = 1; + signedints = 0; + + while (1) { + _getstr(f, &nbytes, string); + if (!strcmp(string, "HEADER_END")) break; + totalbytes += nbytes; + if (!strcmp(string, "rawdatafile")) { + expecting_rawdatafile = 1; + } else if (!strcmp(string, "source_name")) { + expecting_source_name = 1; + } else if (!strcmp(string, "az_start")) { + chkfread(&(az_start), sizeof(double), 1, f); + totalbytes += sizeof(double); + hdr["az_start"] = az_start; + } else if (!strcmp(string, "za_start")) { + chkfread(&(za_start), sizeof(double), 1, f); + totalbytes += sizeof(double); + hdr["za_start"] = za_start; + } else if (!strcmp(string, "src_raj")) { + chkfread(&(src_raj), sizeof(double), 1, f); + totalbytes += sizeof(double); + hdr["src_raj"] = src_raj; + } else if (!strcmp(string, "src_dej")) { + chkfread(&(src_dej), sizeof(double), 1, f); + totalbytes += sizeof(double); + hdr["src_dej"] = src_dej; + } else if (!strcmp(string, "tstart")) { + chkfread(&(tstart), sizeof(double), 1, f); + totalbytes += sizeof(double); + hdr["tstart"] = tstart; + } else if (!strcmp(string, "tsamp")) { + chkfread(&(tsamp), sizeof(double), 1, f); + totalbytes += sizeof(double); + hdr["tsamp"] = tsamp; + } else if (!strcmp(string, "fch1")) { + chkfread(&(fch1), sizeof(double), 1, f); + totalbytes += sizeof(double); + hdr["fch1"] = fch1; + } else if (!strcmp(string, "foff")) { + chkfread(&(foff), sizeof(double), 1, f); + totalbytes += sizeof(double); + hdr["foff"] = foff; + } else if (!strcmp(string, "refdm")) { + chkfread(&(refdm), sizeof(double), 1, f); + totalbytes += sizeof(double); + hdr["refdm"] = refdm; + } else if (!strcmp(string, "nchans")) { + chkfread(&(nchans), sizeof(int), 1, f); + totalbytes += sizeof(int); + hdr["nchans"] = nchans; + } else if (!strcmp(string, "telescope_id")) { + chkfread(&(telescope_id), sizeof(int), 1, f); + totalbytes += sizeof(int); + hdr["telescope_id"] = telescope_id; + hdr["telescope"] = _getteln(telescope_id); + } else if (!strcmp(string, "machine_id")) { + chkfread(&(machine_id), sizeof(int), 1, f); + totalbytes += sizeof(int); + hdr["machine_id"] = machine_id; + hdr["machine"] = _getmachn(machine_id); + } else if (!strcmp(string, "data_type")) { + chkfread(&(datatype), sizeof(int), 1, f); + totalbytes += sizeof(int); + hdr["datatype"] = datatype; + } else if (!strcmp(string, "nbits")) { + chkfread(&(nbits), sizeof(int), 1, f); + totalbytes += sizeof(int); + hdr["nbits"] = nbits; + } else if (!strcmp(string, "barycentric")) { + chkfread(&barycentric, sizeof(int), 1, f); + totalbytes += sizeof(int); + } else if (!strcmp(string, "pulsarcentric")) { + chkfread(&pulsarcentric, sizeof(int), 1, f); + totalbytes += sizeof(int); + } else if (!strcmp(string, "nsamples")) { + chkfread(&(nsamples), sizeof(int), 1, f); + totalbytes += sizeof(int); + hdr["nsamples"] = nsamples; + } else if (!strcmp(string, "nifs")) { + chkfread(&(nifs), sizeof(int), 1, f); + if (nifs > 1) sumifs = 0; + totalbytes += sizeof(int); + hdr["nifs"] = nifs; + hdr["sumifs"] = sumifs; + } else if (!strcmp(string, "nbeams")) { + chkfread(&(nbeams), sizeof(int), 1, f); + totalbytes += sizeof(int); + hdr["nbeams"] = nbeams; + } else if (!strcmp(string, "ibeam")) { + chkfread(&(ibeam), sizeof(int), 1, f); + totalbytes += sizeof(int); + hdr["ibeam"] = ibeam; + } else if (!strcmp(string, "signed")) { + char tmp; + chkfread(&(signedints), sizeof(char), 1, f); + totalbytes += sizeof(char); + hdr["signedints"] = signedints; + } else if (expecting_rawdatafile) { + strcpy(inpfile, string); + hdr["rawdatafile"] = inpfile; + expecting_rawdatafile = 0; + } else if (expecting_source_name) { + strcpy(srcname, string); + hdr["source_name"] = srcname; + expecting_source_name = 0; + } + } + fclose(f); + totalbytes += nbytes; + hdr["hdrlen"] = totalbytes; + return hdr; +} + +void writehdr(nb::dict &hdr, std::string fn) { + FILE *f = chkfopen(fn.c_str(), "wb"); + + _sendstr("HEADER_START", f); + for (auto item : hdr) { + std::string key = nb::cast(item.first); + if (!strcmp(key.c_str(), "rawdatafile")) { + _sendstr(key, f); + _sendstr(nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "source_name")) { + _sendstr(key, f); + _sendstr(nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "az_start")) { + _senddbl(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "za_start")) { + _senddbl(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "src_raj")) { + _senddbl(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "src_dej")) { + _senddbl(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "tstart")) { + _senddbl(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "tsamp")) { + _senddbl(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "fch1")) { + _senddbl(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "foff")) { + _senddbl(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "refdm")) { + _senddbl(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "nchans")) { + _sendint(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "telescope_id")) { + _sendint(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "machine_id")) { + _sendint(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "data_type")) { + _sendint(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "nbits")) { + _sendint(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "barycentric")) { + _sendint(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "pulsarcentric")) { + _sendint(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "nsamples")) { + _sendint(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "nifs")) { + _sendint(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "nbeams")) { + _sendint(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "ibeam")) { + _sendint(key, nb::cast(item.second), f); + } else if (!strcmp(key.c_str(), "signed")) { + _sendint(key, nb::cast(item.second), f); + } + } + _sendstr("HEADER_END", f); + fclose(f); +} + +void init_hdr(nb::module_ m) { + m.def("readhdr", &readhdr); + m.def("writehdr", &writehdr); +} diff --git a/src/priwo/hdr.h b/src/priwo/hdr.h new file mode 100644 index 0000000..4c41990 --- /dev/null +++ b/src/priwo/hdr.h @@ -0,0 +1,10 @@ +#ifndef HDR +#define HDR + +#include "common.h" + +void init_hdr(nb::module_ m); +nb::dict readhdr(std::string fn); +void writehdr(nb::dict &dict, std::string fn); + +#endif diff --git a/src/priwo/inf.C b/src/priwo/inf.C new file mode 100644 index 0000000..224acfd --- /dev/null +++ b/src/priwo/inf.C @@ -0,0 +1,446 @@ +#include +#include +#include + +#include "inf.h" + +constexpr int NUMBANDS = 6; +constexpr int NUMSCOPES = 8; +constexpr int NUMFILTERS = 11; +constexpr int MAXNUMONOFF = 40; + +constexpr char BANDS[NUMBANDS][40] = {"Radio", "IR", "Optical", + "UV", "X-ray", "Gamma"}; + +constexpr char scopes[NUMSCOPES][40] = { + "None (Artificial Data Set)", "Arecibo", "Parkes", "VLA", "MMT", + "Las Campanas 2.5m", "Mt. Hopkins 48in", "Other"}; + +void _getval(FILE *f, char *val, const char *key) { + int ii, slen; + char line[250]; + char *sptr = NULL; + + sptr = fgets(line, 250, f); + if (sptr != NULL && sptr[0] != '\n' && 0 != (ii = strlen(sptr))) { + if (ii >= 40 && line[40] == '=') + sptr = line + 41; + else { + while (--ii >= 0) + if (sptr[ii] == '=') break; + if (ii + 1 == 0) { + sprintf(line, + "No '=' to separate key/value while looking for '%s'! " + "Exiting...", + key); + throw std::runtime_error(line); + } + sptr = line + ii + 1; + } + sptr = rmspace(sptr); + slen = strlen(sptr); + if (slen) { + if ((strcmp(key, "name") == 0 && slen > 199) || + (strcmp(key, "telescope") == 0 && slen > 39) || + (strcmp(key, "band") == 0 && slen > 39) || + (strcmp(key, "name") != 0 && slen > 99)) { + sprintf(line, "Value is too long (%d char) for '%s'! Exiting...", slen, + key); + throw std::runtime_error(line); + } + strcpy(val, sptr); + } else + strcpy(val, "Unknown"); + return; + } else { + if (feof(f)) + sprintf(line, "Reached EOF while looking for '%s'! Exiting...", key); + else + sprintf(line, "Found blank line while looking for '%s'! Exiting...", key); + throw std::runtime_error(line); + } +} + +double _str2dbl(char *val, const char *key) { + double retval; + + char *endptr; + char err[100]; + char *sptr = val; + + retval = strtod(sptr, &endptr); + if (retval == 0.0 && endptr == val) { + sprintf(err, "Cannot convert '%s' to double (%s)! Exiting...", val, key); + throw std::runtime_error(err); + } + return retval; +} + +long _str2long(char *val, const char *key) { + long retval; + + char *endptr; + char err[100]; + char *sptr = val; + + errno = 0; + retval = strtol(sptr, &endptr, 10); + + if ((errno == ERANGE && (retval == LONG_MAX || retval == LONG_MIN)) || + (errno != 0 && retval == 0)) { + sprintf(err, "Cannot convert '%s' to long (%s)! Exiting...", val, key); + throw std::runtime_error(err); + } + + if (endptr == val) { + sprintf(err, "No digits found in '%s' for %s! Exiting...", val, key); + throw std::runtime_error(err); + } + + return retval; +} + +void _str2radec(char *radec, int *h_or_d, int *m, double *s) { + int retval; + radec = rmspace(radec); + retval = sscanf(radec, "%d:%d:%lf\n", h_or_d, m, s); + if (retval != 3) { + char err[100]; + sprintf(err, "Cannot convert '%s' to RA or DEC! Exiting...", radec); + throw std::runtime_error(err); + } + if (radec[0] == '-' && *h_or_d == 0) { + *m = -*m; + *s = -*s; + } +} + +void _radec2str(char *radec, int h_or_d, int m, double s) { + int offset = 0; + if (h_or_d == 0 && (m < 0 || s < 0.0)) { + radec[0] = '-'; + offset = 1; + } + sprintf(radec + offset, "%.2d:%.2d:%07.4f", h_or_d, abs(m), fabs(s)); +} + +nb::dict readinf(std::string fn) { + nb::dict inf; + + char *sptr; + int ii, retval, noteslen = 0; + char tmp1[100], tmp2[100], tmp3[100]; + + FILE *f = chkfopen(fn.c_str(), "r"); + + double ra_s; /* Right ascension seconds (J2000) */ + double dec_s; /* Declination seconds (J2000) */ + double N; /* Number of bins in the time series */ + double dt; /* Width of each time series bin (sec) */ + double fov; /* Diameter of Beam or FOV in arcsec */ + double mjd_f; /* Epoch of observation (MJD) frac part */ + double dm; /* Radio -- Dispersion Measure (cm-3 pc) */ + double freq; /* Radio -- Low chan central freq (Mhz) */ + double freqband; /* Radio -- Total Bandwidth (Mhz) */ + double chan_wid; /* Radio -- Channel Bandwidth (Mhz) */ + double wavelen; /* IR,Opt,UV -- central wavelength (nm) */ + double waveband; /* IR,Opt,UV -- bandpass (nm) */ + double energy; /* x-ray,gamma -- central energy (kev) */ + double energyband; /* x-ray,gamma -- energy bandpass (kev) */ + double onoff[MAXNUMONOFF * 2]; /* Bin number pairs where obs is "on" */ + int num_chan; /* Radio -- Number Channels */ + int mjd_i; /* Epoch of observation (MJD) int part */ + int ra_h; /* Right ascension hours (J2000) */ + int ra_m; /* Right ascension minutes (J2000) */ + int dec_d; /* Declination degrees (J2000) */ + int dec_m; /* Declination minutes (J2000) */ + int bary; /* Barycentered? 1=yes, 0=no */ + int numonoff; /* The number of onoff pairs in the data */ + char notes[500]; /* Any additional notes */ + char name[200]; /* Data file name without suffix */ + char object[100]; /* Object being observed */ + char instrument[100]; /* Instrument used */ + char observer[100]; /* Observer[s] for the data set */ + char analyzer[100]; /* Who analyzed the data */ + char telescope[40]; /* Telescope used */ + char band[40]; /* Type of observation (EM band) */ + char filt[7]; /* IR,Opt,UV -- Photometric Filter */ + + _getval(f, name, "name"); + _getval(f, telescope, "telescope"); + + inf["name"] = name; + inf["telescope"] = telescope; + + int isfake = strcmp(telescope, scopes[0]) == 0; + + if (!isfake) { + _getval(f, instrument, "instrument"); + _getval(f, object, "object"); + _getval(f, tmp1, "RA string"); + _str2radec(tmp1, &ra_h, &ra_m, &ra_s); + _getval(f, tmp1, "DEC string"); + _str2radec(tmp1, &dec_d, &dec_m, &dec_s); + _getval(f, observer, "observer"); + _getval(f, tmp1, "MJD string"); + retval = sscanf(tmp1, "%d.%s", &mjd_i, tmp2); + if (retval != 2) { + sprintf(tmp3, "Cannot parse MJD string '%s'! Exiting...", tmp1); + throw std::runtime_error(tmp3); + } + sprintf(tmp3, "0.%s", tmp2); + mjd_f = _str2dbl(tmp3, "mjd_f"); + _getval(f, tmp1, "bary"); + bary = _str2long(tmp1, "bary"); + + inf["bary"] = bary; + inf["object"] = object; + inf["observer"] = observer; + inf["instrument"] = instrument; + inf["mjd"] = nb::make_tuple(mjd_i, mjd_f); + inf["ra"] = nb::make_tuple(ra_h, ra_m, ra_s); + inf["dec"] = nb::make_tuple(dec_d, dec_m, dec_s); + + } else { + strcpy(object, "fake pulsar"); + inf["object"] = object; + } + + _getval(f, tmp1, "N"); + N = _str2dbl(tmp1, "N"); + _getval(f, tmp1, "dt"); + dt = _str2dbl(tmp1, "dt"); + _getval(f, tmp1, "numonoff"); + numonoff = _str2long(tmp1, "numonoff"); + + inf["N"] = N; + inf["dt"] = dt; + inf["numonoff"] = numonoff; + + if (numonoff) { + ii = 0; + do { + _getval(f, tmp1, "on-off pairs"); + retval = sscanf(tmp1, "%lf %*[ ,] %lf", &onoff[ii], &onoff[ii + 1]); + if (retval != 2) { + sprintf(tmp3, "Cannot parse on-off pair (%d)! Exiting...", ii / 2); + throw std::runtime_error(tmp3); + } + ii += 2; + } while (onoff[ii - 1] < N - 1 && ii < 2 * MAXNUMONOFF); + numonoff = ii / 2; + if (numonoff >= MAXNUMONOFF) { + sprintf(tmp3, + "Number of onoff pairs (%d) >= MAXNUMONOFF (%d)! Exiting...", + numonoff, MAXNUMONOFF); + throw std::runtime_error(tmp3); + } + } else { + numonoff = 1; + onoff[0] = 0; + onoff[1] = N - 1; + } + inf["numonoff"] = numonoff; + + nb::list onofflist; + if (numonoff) { + for (ii = 0; ii < numonoff; ii++) { + onofflist.append(nb::make_tuple(onoff[2 * ii], onoff[2 * ii + 1])); + } + } + inf["onoff"] = onofflist; + + if (!isfake) { + _getval(f, band, "band"); + inf["band"] = band; + + if (strcmp(band, BANDS[0]) == 0) { + _getval(f, tmp1, "fov"); + fov = _str2dbl(tmp1, "fov"); + _getval(f, tmp1, "dm"); + dm = _str2dbl(tmp1, "dm"); + _getval(f, tmp1, "freq"); + freq = _str2dbl(tmp1, "freq"); + _getval(f, tmp1, "freqband"); + freqband = _str2dbl(tmp1, "freqband"); + _getval(f, tmp1, "num_chan"); + num_chan = _str2long(tmp1, "num_chan"); + _getval(f, tmp1, "chan_wid"); + chan_wid = _str2dbl(tmp1, "chan_wid"); + + inf["dm"] = dm; + inf["fov"] = fov; + inf["freq"] = freq; + inf["num_chan"] = num_chan; + inf["freqband"] = freqband; + inf["chan_wid"] = chan_wid; + + } else if ((strcmp(band, BANDS[4]) == 0) || (strcmp(band, BANDS[5]) == 0)) { + _getval(f, tmp1, "fov"); + fov = _str2dbl(tmp1, "fov"); + _getval(f, tmp1, "energy"); + energy = _str2dbl(tmp1, "energy"); + _getval(f, tmp1, "energyband"); + energyband = _str2dbl(tmp1, "energyband"); + + inf["fov"] = fov; + inf["energy"] = energy; + inf["energyband"] = energyband; + + } else { + _getval(f, filt, "filt"); + _getval(f, tmp1, "fov"); + fov = _str2dbl(tmp1, "fov"); + _getval(f, tmp1, "wavelen"); + wavelen = _str2dbl(tmp1, "wavelen"); + _getval(f, tmp1, "waveband"); + waveband = _str2dbl(tmp1, "waveband"); + + inf["fov"] = fov; + inf["filt"] = filt; + inf["wavelen"] = wavelen; + inf["waveband"] = waveband; + } + } + _getval(f, analyzer, "analyzer"); + inf["analyzer"] = analyzer; + + sptr = fgets(tmp1, 100, f); + while (1) { + sptr = fgets(tmp1, 100, f); + if (noteslen + strlen(tmp1) > 500) break; + if (sptr) { + if (noteslen == 0) + strcpy(notes + noteslen, rmlead(tmp1)); + else + strcpy(notes + noteslen, tmp1); + noteslen += strlen(notes); + } else { + if (feof(f)) break; + } + } + inf["notes"] = notes; + + fclose(f); + return inf; +} + +void writeinf(nb::dict inf, std::string fn) { + + int ii, itmp; + char tmp1[100], tmp2[100]; + + FILE *f = chkfopen(fn.c_str(), "w"); + + const char *name = nb::cast(inf["name"]); + const char *telescope = nb::cast(inf["telescope"]); + fprintf(f, " Data file name without suffix = %s\n", name); + fprintf(f, " Telescope used = %s\n", telescope); + + int isfake = strcmp(telescope, scopes[0]) == 0; + + if (!isfake) { + int bary = nb::cast(inf["bary"]); + const char *object = nb::cast(inf["object"]); + const char *observer = nb::cast(inf["observer"]); + const char *instrument = nb::cast(inf["instrument"]); + + nb::tuple ra = inf["ra"]; + nb::tuple dec = inf["dec"]; + nb::tuple mjd = inf["mjd"]; + + int ra_h = nb::cast(ra[0]); + int ra_m = nb::cast(ra[1]); + double ra_s = nb::cast(ra[2]); + + int dec_d = nb::cast(dec[0]); + int dec_m = nb::cast(dec[1]); + double dec_s = nb::cast(dec[2]); + + int mjd_i = nb::cast(mjd[0]); + double mjd_f = nb::cast(mjd[1]); + + fprintf(f, " Instrument used = %s\n", instrument); + fprintf(f, " Object being observed = %s\n", object); + _radec2str(tmp1, ra_h, ra_m, ra_s); + fprintf(f, " J2000 Right Ascension (hh:mm:ss.ssss) = %s\n", tmp1); + _radec2str(tmp1, dec_d, dec_m, dec_s); + fprintf(f, " J2000 Declination (dd:mm:ss.ssss) = %s\n", tmp1); + fprintf(f, " Data observed by = %s\n", observer); + sprintf(tmp1, "%.15f", mjd_f); + sscanf(tmp1, "%d.%s", &itmp, tmp2); + fprintf(f, " Epoch of observation (MJD) = %d.%s\n", mjd_i, + tmp2); + fprintf(f, " Barycentered? (1 yes, 0 no) = %d\n", bary); + } + nb::list onoffs = inf["onoff"]; + double N = nb::cast(inf["N"]); + double dt = nb::cast(inf["dt"]); + int numonoff = nb::cast(inf["numonoff"]); + + fprintf(f, " Number of bins in the time series = %-11.0f\n", N); + fprintf(f, " Width of each time series bin (sec) = %.15g\n", dt); + fprintf(f, " Any breaks in the data? (1 yes, 0 no) = %d\n", + numonoff > 1 ? 1 : 0); + + if (numonoff > 1) { + for (ii = 0; ii < numonoff; ii++) { + fprintf(f, + " On/Off bin pair #%3d = %-11.0f, %-11.0f\n", + ii + 1, nb::cast(onoffs[ii][0]), + nb::cast(onoffs[ii][1])); + } + } + + if (!isfake) { + const char *band = nb::cast(inf["band"]); + fprintf(f, " Type of observation (EM band) = %s\n", band); + if (strcmp(band, BANDS[0]) == 0) { + double fov = nb::cast(inf["fov"]); + double dm = nb::cast(inf["dm"]); + double freq = nb::cast(inf["freq"]); + int num_chan = nb::cast(inf["num_chan"]); + double freqband = nb::cast(inf["freqband"]); + double chan_wid = nb::cast(inf["chan_wid"]); + + fprintf(f, " Beam diameter (arcsec) = %.0f\n", fov); + fprintf(f, " Dispersion measure (cm-3 pc) = %.12g\n", dm); + fprintf(f, " Central freq of low channel (MHz) = %.12g\n", freq); + fprintf(f, " Total bandwidth (MHz) = %.12g\n", + freqband); + fprintf(f, " Number of channels = %d\n", num_chan); + fprintf(f, " Channel bandwidth (MHz) = %.12g\n", + chan_wid); + } else if ((strcmp(band, BANDS[4]) == 0) || (strcmp(band, BANDS[5]) == 0)) { + double fov = nb::cast(inf["fov"]); + double energy = nb::cast(inf["energy"]); + double energyband = nb::cast(inf["energyband"]); + + fprintf(f, " Field-of-view diameter (arcsec) = %.2f\n", fov); + fprintf(f, " Central energy (kev) = %.1f\n", energy); + fprintf(f, " Energy bandpass (kev) = %.1f\n", + energyband); + } else { + double fov = nb::cast(inf["fov"]); + double wavelen = nb::cast(inf["wavelen"]); + double waveband = nb::cast(inf["waveband"]); + const char *filt = nb::cast(inf["filt"]); + + fprintf(f, " Photometric filter used = %s\n", filt); + fprintf(f, " Field-of-view diameter (arcsec) = %.2f\n", fov); + fprintf(f, " Central wavelength (nm) = %.1f\n", wavelen); + fprintf(f, " Bandpass (nm) = %.1f\n", waveband); + } + } + const char *notes = nb::cast(inf["notes"]); + const char *analyzer = nb::cast(inf["analyzer"]); + + fprintf(f, " Data analyzed by = %s\n", analyzer); + fprintf(f, " Any additional notes:\n %s\n\n", notes); + fclose(f); +} + +void init_inf(nb::module_ m) { + m.def("readinf", &readinf); + m.def("writeinf", &writeinf); +} diff --git a/src/priwo/inf.h b/src/priwo/inf.h new file mode 100644 index 0000000..a275772 --- /dev/null +++ b/src/priwo/inf.h @@ -0,0 +1,10 @@ +#ifndef INF +#define INF + +#include "common.h" + +void init_inf(nb::module_ m); +nb::dict readinf(std::string fn); +void writeinf(nb::dict inf, std::string fn); + +#endif diff --git a/src/priwo/pfd.C b/src/priwo/pfd.C new file mode 100644 index 0000000..57fdeb1 --- /dev/null +++ b/src/priwo/pfd.C @@ -0,0 +1,400 @@ +#include "pfd.h" + +nb::dict readpfd(std::string fn) { + nb::dict pfd; + + char temp[16]; + int tmp, byteswap = 0; + FILE *f = chkfopen(fn.c_str(), "rb"); + + int numdms; /* Number of 'dms' */ + int numperiods; /* Number of 'periods' */ + int numpdots; /* Number of 'pdots' */ + int nsub; /* Number of frequency subbands folded */ + int npart; /* Number of folds in time over integration */ + int proflen; /* Number of bins per profile */ + int numchan; /* Number of channels for radio data */ + int pstep; /* Minimum period stepsize in profile phase bins */ + int pdstep; /* Minimum p-dot stepsize in profile phase bins */ + int dmstep; /* Minimum DM stepsize in profile phase bins */ + int ndmfact; /* 2*ndmfact*proflen+1 DMs to search */ + int npfact; /* 2*npfact*proflen+1 periods and p-dots to search */ + char *filenm; /* Filename of the folded data */ + char *candnm; /* String describing the candidate */ + char *telescope; /* Telescope where observation took place */ + char *pgdev; /* PGPLOT device to use */ + char rastr[16]; /* J2000 RA string in format hh:mm:ss.ssss */ + char decstr[16]; /* J2000 DEC string in format dd:mm:ss.ssss */ + double dt; /* Sampling interval of the data */ + double startT; /* Fraction of observation file to start folding */ + double endT; /* Fraction of observation file to stop folding */ + double tepoch; /* Topocentric eopch of data in MJD */ + double bepoch; /* Barycentric eopch of data in MJD */ + double avgvoverc; /* Average topocentric velocity */ + double lofreq; /* Center of low frequency radio channel */ + double chan_wid; /* Width of each radio channel in MHz */ + double bestdm; /* Best DM */ + nb::dict topo; /* Best topocentric p, pd, and pdd */ + nb::dict bary; /* Best barycentric p, pd, and pdd */ + nb::dict fold; /* f, fd, and fdd used to fold the initial data */ + nb::dict orb; /* Barycentric orbital parameters used in folds */ + double *dms; /* DMs used in the trials */ + double *periods; /* Periods used in the trials */ + double *pdots; /* P-dots used in the trials */ + double *rawfolds; /* Raw folds (nsub * npart * proflen points) */ + double *foldstats; /* Statistics for the raw folds */ + + numdms = parseint(f, byteswap); + numperiods = parseint(f, byteswap); + numpdots = parseint(f, byteswap); + nsub = parseint(f, byteswap); + npart = parseint(f, byteswap); + + if (npart < 1 || npart > 10000) { + byteswap = 1; + + nsub = swapint(nsub); + npart = swapint(npart); + numdms = swapint(numdms); + numpdots = swapint(numpdots); + numperiods = swapint(numperiods); + } + + proflen = parseint(f, byteswap); + numchan = parseint(f, byteswap); + pstep = parseint(f, byteswap); + pdstep = parseint(f, byteswap); + dmstep = parseint(f, byteswap); + ndmfact = parseint(f, byteswap); + npfact = parseint(f, byteswap); + + tmp = parseint(f, byteswap); + filenm = (char *)calloc(tmp + 1, sizeof(char)); + chkfread(filenm, sizeof(char), tmp, f); + + tmp = parseint(f, byteswap); + candnm = (char *)calloc(tmp + 1, sizeof(char)); + chkfread(candnm, sizeof(char), tmp, f); + + tmp = parseint(f, byteswap); + telescope = (char *)calloc(tmp + 1, sizeof(char)); + chkfread(telescope, sizeof(char), tmp, f); + + tmp = parseint(f, byteswap); + pgdev = (char *)calloc(tmp + 1, sizeof(char)); + chkfread(pgdev, sizeof(char), tmp, f); + + { + int ii, haspos = 1; + chkfread(temp, sizeof(char), 16, f); + for (ii = 0; ii < 16; ii++) { + if (!isdigit(temp[ii]) && temp[ii] != ':' && temp[ii] != '.' && + temp[ii] != '-' && temp[ii] != '\0') { + haspos = 0; + break; + } + } + if (haspos) { + strcpy(rastr, temp); + chkfread(decstr, sizeof(char), 16, f); + dt = parsedouble(f, byteswap); + startT = parsedouble(f, byteswap); + } else { + strcpy(rastr, "Unknown"); + strcpy(decstr, "Unknown"); + dt = *(double *)(temp + 0); + if (byteswap) dt = swapdouble(dt); + startT = *(double *)(temp + sizeof(double)); + if (byteswap) startT = swapdouble(startT); + } + } + + endT = parsedouble(f, byteswap); + tepoch = parsedouble(f, byteswap); + bepoch = parsedouble(f, byteswap); + avgvoverc = parsedouble(f, byteswap); + lofreq = parsedouble(f, byteswap); + chan_wid = parsedouble(f, byteswap); + bestdm = parsedouble(f, byteswap); + + float topopow = parsefloat(f, byteswap); + parsefloat(f, byteswap); + double topop1 = parsedouble(f, byteswap); + double topop2 = parsedouble(f, byteswap); + double topop3 = parsedouble(f, byteswap); + + float barypow = parsefloat(f, byteswap); + parsefloat(f, byteswap); + double baryp1 = parsedouble(f, byteswap); + double baryp2 = parsedouble(f, byteswap); + double baryp3 = parsedouble(f, byteswap); + + float foldpow = parsefloat(f, byteswap); + parsefloat(f, byteswap); + double foldp1 = parsedouble(f, byteswap); + double foldp2 = parsedouble(f, byteswap); + double foldp3 = parsedouble(f, byteswap); + + double orbp = parsedouble(f, byteswap); + double orbe = parsedouble(f, byteswap); + double orbx = parsedouble(f, byteswap); + double orbw = parsedouble(f, byteswap); + double orbt = parsedouble(f, byteswap); + double orbpd = parsedouble(f, byteswap); + double orbwd = parsedouble(f, byteswap); + + dms = (double *)calloc(numdms, sizeof(double)); + chkfread(dms, sizeof(double), numdms, f); + + periods = (double *)calloc(numperiods, sizeof(double)); + chkfread(periods, sizeof(double), numperiods, f); + + pdots = (double *)calloc(numpdots, sizeof(double)); + chkfread(pdots, sizeof(double), numpdots, f); + + rawfolds = (double *)calloc(nsub * npart * proflen, sizeof(double)); + chkfread(rawfolds, sizeof(double), nsub * npart * proflen, f); + + foldstats = (double *)calloc(nsub * npart * 7, sizeof(double)); + chkfread(foldstats, sizeof(double), nsub * npart * 7, f); + + if (byteswap) { + int ii; + + for (ii = 0; ii < numdms; ii++) dms[ii] = swapdouble(dms[ii]); + for (ii = 0; ii < numpdots; ii++) pdots[ii] = swapdouble(pdots[ii]); + for (ii = 0; ii < numperiods; ii++) periods[ii] = swapdouble(periods[ii]); + + for (ii = 0; ii < nsub * npart * proflen; ii++) + rawfolds[ii] = swapdouble(rawfolds[ii]); + + for (ii = 0; ii < nsub * npart * 7; ii++) { + foldstats[ii] = swapdouble(foldstats[ii]); + } + } + + pfd["numdms"] = numdms; + pfd["numperiods"] = numperiods; + pfd["numpdots"] = numpdots; + pfd["nsub"] = nsub; + pfd["npart"] = npart; + + pfd["proflen"] = proflen; + pfd["numchan"] = numchan; + pfd["pstep"] = pstep; + pfd["pdstep"] = pdstep; + pfd["dmstep"] = dmstep; + pfd["ndmfact"] = ndmfact; + pfd["npfact"] = npfact; + + pfd["filenm"] = filenm; + pfd["candnm"] = candnm; + pfd["telescope"] = telescope; + pfd["pgdev"] = pgdev; + + pfd["rastr"] = rastr; + pfd["decstr"] = decstr; + pfd["dt"] = dt; + pfd["startT"] = startT; + pfd["endT"] = endT; + pfd["tepoch"] = tepoch; + pfd["bepoch"] = bepoch; + pfd["avgvoverc"] = avgvoverc; + pfd["lofreq"] = lofreq; + pfd["chan_wid"] = chan_wid; + pfd["bestdm"] = bestdm; + + topo["pow"] = topopow; + topo["p1"] = topop1; + topo["p2"] = topop2; + topo["p3"] = topop3; + pfd["topo"] = topo; + + bary["pow"] = barypow; + bary["p1"] = baryp1; + bary["p2"] = baryp2; + bary["p3"] = baryp3; + pfd["bary"] = bary; + + fold["pow"] = foldpow; + fold["p1"] = foldp1; + fold["p2"] = foldp2; + fold["p3"] = foldp3; + pfd["fold"] = fold; + + orb["p"] = orbp; + orb["e"] = orbe; + orb["x"] = orbx; + orb["w"] = orbw; + orb["t"] = orbt; + orb["pd"] = orbpd; + orb["wd"] = orbwd; + pfd["orb"] = orb; + + size_t dmx[1] = {static_cast(numdms)}; + pfd["dms"] = nb::ndarray(dms, 1, dmx, nb::handle()); + + size_t px[1] = {static_cast(numperiods)}; + pfd["periods"] = nb::ndarray(periods, 1, px, nb::handle()); + + size_t pdotx[1] = {static_cast(numpdots)}; + pfd["pdots"] = nb::ndarray(pdots, 1, pdotx, nb::handle()); + + size_t foldx[3] = {static_cast(nsub), static_cast(npart), + static_cast(proflen)}; + pfd["rawfolds"] = + nb::ndarray(rawfolds, 3, foldx, nb::handle()); + + size_t statx[3] = {static_cast(nsub), static_cast(npart), 7}; + pfd["stats"] = + nb::ndarray(foldstats, 3, statx, nb::handle()); + + fclose(f); + return pfd; +} + +void writepfd(nb::dict pfd, std::string fn) { + FILE *f = chkfopen(fn.c_str(), "wb"); + int itmp; + + int numdms = nb::cast(pfd["numdms"]); + int numperiods = nb::cast(pfd["numperiods"]); + int numpdots = nb::cast(pfd["numpdots"]); + int nsub = nb::cast(pfd["nsub"]); + int npart = nb::cast(pfd["npart"]); + int proflen = nb::cast(pfd["proflen"]); + int numchan = nb::cast(pfd["numchan"]); + int pstep = nb::cast(pfd["pstep"]); + int pdstep = nb::cast(pfd["pdstep"]); + int dmstep = nb::cast(pfd["dmstep"]); + int ndmfact = nb::cast(pfd["ndmfact"]); + int npfact = nb::cast(pfd["npfact"]); + + std::string filenm = nb::cast(pfd["filenm"]); + std::string candnm = nb::cast(pfd["candnm"]); + std::string telescope = nb::cast(pfd["telescope"]); + std::string pgdev = nb::cast(pfd["pgdev"]); + + std::string rastr = nb::cast(pfd["rastr"]); + std::string decstr = nb::cast(pfd["decstr"]); + + double dt = nb::cast(pfd["dt"]); + double startT = nb::cast(pfd["startT"]); + double endT = nb::cast(pfd["endT"]); + double tepoch = nb::cast(pfd["tepoch"]); + double bepoch = nb::cast(pfd["bepoch"]); + double avgvoverc = nb::cast(pfd["avgvoverc"]); + double lofreq = nb::cast(pfd["lofreq"]); + double chan_wid = nb::cast(pfd["chan_wid"]); + double bestdm = nb::cast(pfd["bestdm"]); + + chkfwrite(&numdms, sizeof(int), 1, f); + chkfwrite(&numperiods, sizeof(int), 1, f); + chkfwrite(&numpdots, sizeof(int), 1, f); + chkfwrite(&nsub, sizeof(int), 1, f); + chkfwrite(&npart, sizeof(int), 1, f); + chkfwrite(&proflen, sizeof(int), 1, f); + chkfwrite(&numchan, sizeof(int), 1, f); + chkfwrite(&pstep, sizeof(int), 1, f); + chkfwrite(&pdstep, sizeof(int), 1, f); + chkfwrite(&dmstep, sizeof(int), 1, f); + chkfwrite(&ndmfact, sizeof(int), 1, f); + chkfwrite(&npfact, sizeof(int), 1, f); + + itmp = strlen(filenm.c_str()); + chkfwrite(&itmp, sizeof(int), 1, f); + chkfwrite(filenm.data(), sizeof(char), itmp, f); + + itmp = strlen(candnm.c_str()); + chkfwrite(&itmp, sizeof(int), 1, f); + chkfwrite(candnm.data(), sizeof(char), itmp, f); + + itmp = strlen(telescope.c_str()); + chkfwrite(&itmp, sizeof(int), 1, f); + chkfwrite(telescope.data(), sizeof(char), itmp, f); + + itmp = strlen(pgdev.c_str()); + chkfwrite(&itmp, sizeof(int), 1, f); + chkfwrite(pgdev.data(), sizeof(char), itmp, f); + + double topopow = nb::cast(pfd["topo"]["pow"]); + double topop1 = nb::cast(pfd["topo"]["p1"]); + double topop2 = nb::cast(pfd["topo"]["p2"]); + double topop3 = nb::cast(pfd["topo"]["p3"]); + + double barypow = nb::cast(pfd["bary"]["pow"]); + double baryp1 = nb::cast(pfd["bary"]["p1"]); + double baryp2 = nb::cast(pfd["bary"]["p2"]); + double baryp3 = nb::cast(pfd["bary"]["p3"]); + + double foldpow = nb::cast(pfd["fold"]["pow"]); + double foldp1 = nb::cast(pfd["fold"]["p1"]); + double foldp2 = nb::cast(pfd["fold"]["p2"]); + double foldp3 = nb::cast(pfd["fold"]["p3"]); + + double orbp = nb::cast(pfd["orb"]["p"]); + double orbe = nb::cast(pfd["orb"]["e"]); + double orbx = nb::cast(pfd["orb"]["x"]); + double orbw = nb::cast(pfd["orb"]["w"]); + double orbt = nb::cast(pfd["orb"]["t"]); + double orbpd = nb::cast(pfd["orb"]["pd"]); + double orbwd = nb::cast(pfd["orb"]["wd"]); + + chkfwrite(rastr.data(), sizeof(char), 16, f); + chkfwrite(decstr.data(), sizeof(char), 16, f); + chkfwrite(&dt, sizeof(double), 1, f); + chkfwrite(&startT, sizeof(double), 1, f); + chkfwrite(&endT, sizeof(double), 1, f); + chkfwrite(&tepoch, sizeof(double), 1, f); + chkfwrite(&bepoch, sizeof(double), 1, f); + chkfwrite(&avgvoverc, sizeof(double), 1, f); + chkfwrite(&lofreq, sizeof(double), 1, f); + chkfwrite(&chan_wid, sizeof(double), 1, f); + chkfwrite(&bestdm, sizeof(double), 1, f); + chkfwrite(&topopow, sizeof(double), 1, f); + chkfwrite(&topop1, sizeof(double), 1, f); + chkfwrite(&topop2, sizeof(double), 1, f); + chkfwrite(&topop3, sizeof(double), 1, f); + chkfwrite(&barypow, sizeof(double), 1, f); + chkfwrite(&baryp1, sizeof(double), 1, f); + chkfwrite(&baryp2, sizeof(double), 1, f); + chkfwrite(&baryp3, sizeof(double), 1, f); + chkfwrite(&foldpow, sizeof(double), 1, f); + chkfwrite(&foldp1, sizeof(double), 1, f); + chkfwrite(&foldp2, sizeof(double), 1, f); + chkfwrite(&foldp3, sizeof(double), 1, f); + chkfwrite(&orbp, sizeof(double), 1, f); + chkfwrite(&orbe, sizeof(double), 1, f); + chkfwrite(&orbx, sizeof(double), 1, f); + chkfwrite(&orbw, sizeof(double), 1, f); + chkfwrite(&orbt, sizeof(double), 1, f); + chkfwrite(&orbpd, sizeof(double), 1, f); + chkfwrite(&orbwd, sizeof(double), 1, f); + + nb::ndarray dms = + nb::cast>(pfd["dms"]); + + nb::ndarray periods = + nb::cast>(pfd["periods"]); + + nb::ndarray pdots = + nb::cast>(pfd["pdots"]); + + nb::ndarray rawfolds = + nb::cast>(pfd["rawfolds"]); + + nb::ndarray stats = + nb::cast>(pfd["stats"]); + + chkfwrite(dms.data(), sizeof(double), numdms, f); + chkfwrite(periods.data(), sizeof(double), numperiods, f); + chkfwrite(pdots.data(), sizeof(double), numpdots, f); + chkfwrite(rawfolds.data(), sizeof(double), nsub * npart * proflen, f); + chkfwrite(stats.data(), sizeof(double), nsub * npart * 7, f); + fclose(f); +} + +void init_pfd(nb::module_ m) { + m.def("readpfd", &readpfd); + m.def("writepfd", &writepfd); +} diff --git a/src/priwo/pfd.h b/src/priwo/pfd.h new file mode 100644 index 0000000..55f6969 --- /dev/null +++ b/src/priwo/pfd.h @@ -0,0 +1,10 @@ +#ifndef PFD +#define PFD + +#include "common.h" + +void init_pfd(nb::module_ m); +nb::dict readpfd(std::string fn); +void writepfd(nb::dict pfd, std::string fn); + +#endif diff --git a/src/priwo/polycos.C b/src/priwo/polycos.C new file mode 100644 index 0000000..a6842db --- /dev/null +++ b/src/priwo/polycos.C @@ -0,0 +1,171 @@ +#include "common.h" +#include "nanobind/nanobind.h" +#include "polycos.h" + +#include + +#define TEST_EQUAL(a, b) \ + (fabs(a) == 0.0 ? (fabs((a) - (b)) <= 2 * DBL_EPSILON ? 1 : 0) \ + : (fabs((a) - (b)) / fabs((a)) <= 2 * DBL_EPSILON ? 1 : 0)) + +double _str2dbl(char *val) { + double retval; + + char *endptr; + char err[100]; + char *sptr = val; + + retval = strtod(sptr, &endptr); + if (retval == 0.0 && endptr == val) { + sprintf(err, "Cannot convert %s to double! Exiting...", val); + throw std::runtime_error(err); + } + return retval; +} + +int readpolyco(nb::dict polyco, FILE *f) { + char PSR[11]; + char DATE[10]; + double UTC; + double TMID; + double DM; + double DOPPLER; + double LOG10RMS; + double RPHASE; + double F0; + int OBSNUM; + double DATASPAN; + int NUMCOEFF; + double OBSFREQ; + double BINPHASE = 0.0; + + ssize_t read; + size_t len = 0; + char *line = NULL; + + read = getline(&line, &len, f); + if (read == -1) return -1; + sscanf(line, "%10s %9s %12lf %20lf %21lf %6lf %7lf", PSR, DATE, &UTC, &TMID, + &DM, &DOPPLER, &LOG10RMS); + + polyco["PSR"] = PSR; + polyco["DATE"] = DATE; + polyco["UTC"] = UTC; + polyco["TMID"] = TMID; + polyco["DM"] = DM; + polyco["DOPPLER"] = DOPPLER; + polyco["LOG10RMS"] = LOG10RMS; + + read = getline(&line, &len, f); + if (read == -1) return -1; + sscanf(line, "%20lf %18lf %5d %6lf %5d %21lf %5lf", &RPHASE, &F0, &OBSNUM, + &DATASPAN, &NUMCOEFF, &OBSFREQ, &BINPHASE); + + polyco["RPHASE"] = RPHASE; + polyco["F0"] = F0; + polyco["OBSNUM"] = OBSNUM; + polyco["DATASPAN"] = DATASPAN; + polyco["NUMCOEFF"] = NUMCOEFF; + polyco["OBSFREQ"] = OBSFREQ; + if (TEST_EQUAL(BINPHASE, 0.0)) + polyco["BINPHASE"] = nb::none(); + else + polyco["BINPHASE"] = BINPHASE; + + char tmp1[30], tmp2[30], tmp3[30]; + double *COEFFS = (double *)calloc(NUMCOEFF, sizeof(double)); + + for (int ii = 0; ii < (int)(NUMCOEFF / 3); ii++) { + read = getline(&line, &len, f); + if (read == -1) return -1; + sscanf(line, "%s %s %s", tmp1, tmp2, tmp3); + COEFFS[ii * 3 + 0] = _str2dbl(tmp1); + COEFFS[ii * 3 + 1] = _str2dbl(tmp2); + COEFFS[ii * 3 + 2] = _str2dbl(tmp3); + } + + size_t shape[1] = {static_cast(NUMCOEFF)}; + polyco["COEFFS"] = nb::ndarray>( + COEFFS, 1, shape, nb::handle()); + + return 0; +} + +nb::list readpolycos(std::string fn) { + nb::list polycos; + + FILE *f = chkfopen(fn.c_str(), "r"); + + int status = 0; + while (status == 0) { + nb::dict polyco; + status = readpolyco(polyco, f); + if (!(polyco.size() == 0)) polycos.append(polyco); + } + + fclose(f); + return polycos; +} + +void writepolyco(nb::dict polyco, FILE *f) { + const char *PSR = nb::cast(polyco["PSR"]); + const char *DATE = nb::cast(polyco["DATE"]); + double UTC = nb::cast(polyco["UTC"]); + double TMID = nb::cast(polyco["TMID"]); + double DM = nb::cast(polyco["DM"]); + double DOPPLER = nb::cast(polyco["DOPPLER"]); + + fprintf(f, "%10s", PSR); + fprintf(f, "%10s", DATE); + fprintf(f, "%11lf", UTC); + fprintf(f, "%20.11lf", TMID); + fprintf(f, "%21.6lf", DM); + fprintf(f, "%7.3lf", DOPPLER); + + if (!polyco["LOG10RMS"].is_none()) { + double LOG10RMS = nb::cast(polyco["LOG10RMS"]); + fprintf(f, "%7.3lf", LOG10RMS); + } + fprintf(f, "\n"); + + double RPHASE = nb::cast(polyco["RPHASE"]); + double F0 = nb::cast(polyco["F0"]); + int OBSNUM = nb::cast(polyco["OBSNUM"]); + double DATASPAN = nb::cast(polyco["DATASPAN"]); + int NUMCOEFF = nb::cast(polyco["NUMCOEFF"]); + double OBSFREQ = nb::cast(polyco["OBSFREQ"]); + + fprintf(f, "%20.6lf", RPHASE); + fprintf(f, "%20.12lf", F0); + fprintf(f, "%5d", OBSNUM); + fprintf(f, "%5.0lf", DATASPAN); + fprintf(f, "%4d", NUMCOEFF); + fprintf(f, "%21.3lf", OBSFREQ); + + if (!polyco["BINPHASE"].is_none()) { + double BINPHASE = nb::cast(polyco["BINPHASE"]); + fprintf(f, "%4.0lf", BINPHASE); + } + fprintf(f, "\n"); + + nb::ndarray> COEFFS = + nb::cast>>(polyco["COEFFS"]); + for (int ii = 0; ii < (int)(NUMCOEFF / 3); ii++) { + fprintf(f, "%25.16E %25.16E %25.16E\n", COEFFS(ii * 3 + 0), + COEFFS(ii * 3 + 1), COEFFS(ii * 3 + 2)); + } +} + +void writepolycos(nb::list polycos, std::string fn) { + FILE *f = chkfopen(fn.c_str(), "w"); + for (auto object : polycos) { + nb::dict polyco = nb::cast(object); + writepolyco(polyco, f); + } + fclose(f); +} + +void init_polycos(nb::module_ m) { + m.def("readpolycos", &readpolycos); + m.def("writepolycos", &writepolycos); +} diff --git a/src/priwo/polycos.h b/src/priwo/polycos.h new file mode 100644 index 0000000..6091339 --- /dev/null +++ b/src/priwo/polycos.h @@ -0,0 +1,10 @@ +#ifndef POLYCOS +#define POLYCOS + +#include "common.h" + +void init_polycos(nb::module_ m); +nb::list readpolycos(std::string fn); +void writepolycos(nb::list polycos, std::string fn); + +#endif diff --git a/src/priwo/presto/__init__.py b/src/priwo/presto/__init__.py deleted file mode 100644 index 4dab4f0..0000000 --- a/src/priwo/presto/__init__.py +++ /dev/null @@ -1,22 +0,0 @@ -""" -R/W pulsar data from formats used by the PRESTO package. -""" - -from priwo.presto.inf import readinf, writeinf -from priwo.presto.dat import readdat, writedat -from priwo.presto.fft import readfft, writefft -from priwo.presto.bpf import readbpf, writebpf -from priwo.presto.pfd import readpfd, writepfd - -__all__ = [ - "readinf", - "readdat", - "readfft", - "readbpf", - "readpfd", - "writeinf", - "writedat", - "writefft", - "writebpf", - "writepfd", -] diff --git a/src/priwo/presto/bpf.py b/src/priwo/presto/bpf.py deleted file mode 100644 index b7ca18b..0000000 --- a/src/priwo/presto/bpf.py +++ /dev/null @@ -1,138 +0,0 @@ -""" -R/W PRESTO best pulse profile (*.bestprof) files. -""" - -import re -import pabo as pb - - -only = lambda _: _[0] if len(_) == 1 else _ -strings = lambda _: None if str(_) == "N/A" else str(_) -numex = re.compile(r"[+-]?(?:0|[1-9]\d*)(?:\.\d*)?(?:[eE][+\-]?\d+)?") -numeral = ( - lambda _: only( - list( - map( - float, - re.findall(numex, _), - ) - ) - ) - if re.search(numex, _) - else None -) - - -def readbpf(f): - """ - Read in a PRESTO best pulse profile (*.bestprof) file. - """ - - meta = {} - with open(f, "r") as fp: - lines = fp.read() - header = re.findall(r"^#.*", lines, re.M)[:-1] - for key, val in { - name: only([validator(_) for _ in string.strip().split("+/-")]) - for (name, validator), string in zip( - { - # fmt: off - "filename": strings, - "candidate": strings, - "telescope": strings, - "topo_epoch": numeral, - "bary_epoch": numeral, - "tsamp": numeral, - "nsamp": numeral, - "data_avg": numeral, - "data_std": numeral, - "prof_bins": numeral, - "prof_avg": numeral, - "prof_std": numeral, - "red_chi_sqr": numeral, - "noise_sigma": numeral, - "dm": numeral, - "p_topo": numeral, - "pd_topo": numeral, - "pdd_topo": numeral, - "p_bary": numeral, - "pd_bary": numeral, - "pdd_bary": numeral, - "p_orb": numeral, - "asini_by_c": numeral, - "eccentricity": numeral, - "w": numeral, - "t_peri": numeral, - # fmt: on - }.items(), - [re.split(r"\s+[=<>]\s+", line)[-1] for line in header], - ) - }.items(): - if isinstance(val, list): - if key == "noise_sigma": - meta[key] = val[-1] - else: - meta[key] = val[0] - meta["".join([key, "_err"])] = val[1] - else: - meta[key] = val - data = pb.Array(pb.Float(4)).parse( - b"".join( - [ - pb.Float(4).build(float(_)) - for _ in re.findall(r"^\s+\d+\s+(.+)$", lines, re.M) - ] - ) - ) - return meta, data - - -def writebpf(meta, data, f): - """ - Write out a PRESTO best pulse profile (*.bestprof) file. - """ - - copy = {} - for key, val in meta.items(): - val = val if val is not None else "N/A" - error = meta.get("".join([key, "_err"]), None) - val = f"{val!s}".rstrip("0").rstrip(".") if isinstance(val, float) else val - copy[key] = f"{val:17} +/- {error}" if error is not None else f"{val}" - - with open(f, "w+") as fp: - fp.write( - "\n".join( - [ - f"# Input file = {copy['filename']}", - f"# Candidate = {copy['candidate']}", - f"# Telescope = {copy['telescope']}", - f"# Epoch_topo = {copy['topo_epoch']}", - f"# Epoch_bary = {copy['bary_epoch']}", - f"# T_sample = {copy['tsamp']}", - f"# Data Folded = {copy['nsamp']}", - f"# Data Avg = {copy['data_avg']}", - f"# Data StdDev = {copy['data_std']}", - f"# Profile Bins = {copy['prof_bins']}", - f"# Profile Avg = {copy['prof_avg']}", - f"# Profile StdDev = {copy['prof_std']}", - f"# Reduced chi-sqr = {copy['red_chi_sqr']}", - f"# Prob(Noise) < 0 (~{copy['noise_sigma']} sigma)", - f"# Best DM = {copy['dm']}", - f"# P_topo (ms) = {copy['p_topo']}", - f"# P'_topo (s/s) = {copy['pd_topo']}", - f"# P''_topo (s/s^2) = {copy['pdd_topo']}", - f"# P_bary (ms) = {copy['p_bary']}", - f"# P'_bary (s/s) = {copy['pd_bary']}", - f"# P''_bary (s/s^2) = {copy['pdd_bary']}", - f"# P_orb (s) = {copy['p_orb']}", - f"# asin(i)/c (s) = {copy['asini_by_c']}", - f"# eccentricity = {copy['eccentricity']}", - f"# w (rad) = {copy['w']}", - f"# T_peri = {copy['t_peri']}", - "######################################################", - ] - ) - ) - fp.write("\n") - for i, value in enumerate(data): - fp.write(f"{i:>4} {value:e}\n") diff --git a/src/priwo/presto/dat.py b/src/priwo/presto/dat.py deleted file mode 100644 index 0844288..0000000 --- a/src/priwo/presto/dat.py +++ /dev/null @@ -1,25 +0,0 @@ -""" -R/W PRESTO time series (*.dat) files. -""" - -import pabo as pb -from pathlib import Path -from priwo.presto.inf import readinf -from priwo.presto.inf import writeinf - - -def readdat(f): - """ - Read in a PRESTO time series (*.dat) file. - """ - - return (readinf(Path(f).with_suffix(".inf")), pb.Array(pb.Float(4)).parse(f)) - - -def writedat(meta, data, f): - """ - Write out a PRESTO time series (*.dat) file. - """ - - writeinf(meta, Path(f).with_suffix(".inf")) - pb.Array(pb.Float(4)).build(data, f) diff --git a/src/priwo/presto/fft.py b/src/priwo/presto/fft.py deleted file mode 100644 index acbe3fe..0000000 --- a/src/priwo/presto/fft.py +++ /dev/null @@ -1,25 +0,0 @@ -""" -R/W PRESTO FFT (*.fft) files. -""" - -import pabo as pb -from pathlib import Path -from priwo.presto.inf import readinf -from priwo.presto.inf import writeinf - - -def readfft(f): - """ - Read in a PRESTO FFT (*.fft) file. - """ - - return (readinf(Path(f).with_suffix(".inf")), pb.Array(pb.Float(4)).parse(f)) - - -def writefft(meta, data, f): - """ - Write out a PRESTO FFT (*.fft) file. - """ - - writeinf(meta, Path(f).with_suffix(".inf")) - pb.Array(pb.Float(4)).build(data, f) diff --git a/src/priwo/presto/inf.py b/src/priwo/presto/inf.py deleted file mode 100644 index d4a8d4a..0000000 --- a/src/priwo/presto/inf.py +++ /dev/null @@ -1,127 +0,0 @@ -""" -R/W PRESTO infodata (*.inf) files. -""" - -import re - -boolean = lambda _: int(_) != 0 -splits = lambda _: [int(__) for __ in _.split(",")] -string = lambda _: {True: "1", False: "0"}.get(_, str(_)) - -# fmt: off -INFMAP = { - "Data file name without suffix": ("filename", str), - "Telescope used": ("telescope", str), - "Instrument used": ("instrument", str), - "Object being observed": ("object", str), - "J2000 Right Ascension (hh:mm:ss.ssss)": ("ra", str), - "J2000 Declination (dd:mm:ss.ssss)": ("dec", str), - "Data observed by": ("observer", str), - "Epoch of observation (MJD)": ("mjd", float), - "Barycentered? (1=yes, 0=no)": ("barycentered", boolean), - "Barycentered? (1 yes, 0 no)": ("barycentered", boolean), - "Number of bins in the time series": ("nsamples", int), - "Width of each time series bin (sec)": ("samptime", float), - "Any breaks in the data? (1=yes, 0=no)": ("breaks", boolean), - "Any breaks in the data? (1 yes, 0 no)": ("breaks", boolean), - "Type of observation (EM band)": ("emband", str), - "Beam diameter (arcsec)": ("beamdiam", float), - "Dispersion measure (cm-3 pc)": ("dm", float), - "Central freq of low channel (MHz)": ("cfreq", float), - "Total bandwidth (MHz)": ("bw", float), - "Number of channels": ("nchannels", int), - "Channel bandwidth (MHz)": ("chanwidth", float), - "Field-of-view diameter (arcsec)": ("fov", float), - "Central energy (kev)": ("cE", float), - "Energy bandpass (kev)": ("bpE", float), - "Photometric filter used": ("filter", str), - "Central wavelength (nm)": ("cwaveln", float), - "Bandpass (nm)": ("bandpass", float), - "Data analyzed by": ("analyst", str), -} -# fmt: on - - -def readinf(f): - - """ - Read in a PRESTO infodata (*.inf) file. - """ - - meta = {} - - regex = re.compile( - r""" - ^ # Beginning of string. - \s+ # Trailing whitespace, if any. - (?P.+?) # Capture key. - \s+=\s+ # Separator. - (?P.+?) # Capture value. - \s+ # Trailing whitespace, if any. - $ # End of line. - """, - re.MULTILINE | re.VERBOSE, - ) - - notes = [] - onoffs = [] - with open(f, "r") as lines: - for line in lines: - matched = re.search(regex, line) - if matched: - groups = matched.groupdict() - key = groups["key"] - val = groups["val"] - try: - nameof, typeof = INFMAP[key] - meta[nameof] = typeof(val) - except KeyError: - onoffs.append(splits(val)) - else: - notes.append(line) - notes = notes[1:] - notes = [note.strip() for note in notes] - notes = [note for note in notes if note] - notes = "\n".join(notes) - meta["onoffs"] = onoffs - meta["notes"] = notes - return meta - - -def writeinf(meta, f): - - """ - Write out a PRESTO infodata (*.inf) file. - """ - - notes = meta.pop("notes") - onoffs = meta.pop("onoffs") - frame = " {key:<37s} = {val:s}" - lines = [ - frame.format( - key={ - name: description - for ( - description, - (name, _), - ) in INFMAP.items() - }[key], - val=string(val), - ) - for key, val in meta.items() - ] - - if meta["breaks"]: - lines[12:12] = [ - frame.format( - val="{:<11d}, {:d}".format(*onoff), - key="On/Off bin pair #{npair:>3}".format(npair=i + 1), - ) - for i, onoff in enumerate(onoffs) - ] - - with open(f, "w+") as fp: - fp.write("\n".join(lines)) - fp.write("\n Any additional notes:\n") - for note in notes.split("\n"): - fp.write(" {note}\n".format(note=note)) diff --git a/src/priwo/presto/pfd.py b/src/priwo/presto/pfd.py deleted file mode 100644 index d47a526..0000000 --- a/src/priwo/presto/pfd.py +++ /dev/null @@ -1,112 +0,0 @@ -""" -R/W PRESTO folded data (*.pfd) files. -""" - -import pabo as pb - - -def prod(x): - res = 1.0 - for i in x: - res *= i - return res - - -# fmt: off -Position = pb.Spec( - fields=dict( - power = pb.Float(8) / "Power normalized with local power.", - p = pb.Float(8) / "r if rzw search, r_startfft if bin search.", - pd = pb.Float(8) / "z if rzw search, r_freqmod if bin search.", - pdd = pb.Float(8) / "w if rzw search, numfftbins if bin search.", - ) -) -# fmt: on - - -# fmt: off -PFD = pb.Spec( - fields=dict( - ndms = pb.Int(4) / "Number of trial DMs.", - nperiods = pb.Int(4) / "Number of trial periods.", - npdots = pb.Int(4) / "Number of trial period derivatives.", - nsub = pb.Int(4) / "Number of frequency subbands folded.", - npart = pb.Int(4) / "Number of folds in time over integration.", - nbin = pb.Int(4) / "Number of bins per profile.", - nchan = pb.Int(4) / "Number of channels for radio data.", - pstep = pb.Int(4) / "Minimum period stepsize in profile phase bins.", - pdstep = pb.Int(4) / "Minimum p-dot stepsize in profile phase bins.", - dmstep = pb.Int(4) / "Minimum DM stepsize in profile phase bins.", - ndmfact = pb.Int(4) / "2 * ndmfact * proflen + 1 DMs to search.", - npfact = pb.Int(4) / "2 * npfact * proflen + 1 periods and p-dots to search.", - filename = pb.PascalString(pb.Int(4), "utf8") / "Filename of the folded data.", - candname = pb.PascalString(pb.Int(4), "utf8") / "String describing the candidate.", - telescope = pb.PascalString(pb.Int(4), "utf8") / "Telescope where observation took place.", - pgdev = pb.PascalString(pb.Int(4), "utf8") / "PGPLOT device to use.", - ra = pb.PaddedString(16, "utf8") / "J2000 RA string in format hh:mm:ss.ssss.", - dec = pb.PaddedString(16, "utf8") / "J2000 DEC string in format dd:mm:ss.ssss.", - dt = pb.Float(8) / "Sampling interval of the data.", - T0 = pb.Float(8) / "Fraction of observation file to start folding.", - Tn = pb.Float(8) / "Fraction of observation file to stop folding.", - tepoch = pb.Float(8) / "Topocentric eopch of data in MJD.", - bepoch = pb.Float(8) / "Barycentric eopch of data in MJD.", - vavg = pb.Float(8) / "Average topocentric velocity.", - f0 = pb.Float(8) / "Center of low frequency radio channel.", - df = pb.Float(8) / "Width of each radio channel in MHz.", - bestdm = pb.Float(8) / "Best DM.", - topo = Position / "Best topocentric p, pd, and pdd.", - bary = Position / "Best barycentric p, pd, and pdd.", - fold = Position / "f, fd, and fdd used to fold the initial data.", - orb=pb.Spec( - fields=dict( - p = pb.Float(8) / "Orbital period (s).", - e = pb.Float(8) / "Orbital eccentricity.", - x = pb.Float(8) / "Projected semi-major axis (lt-sec).", - w = pb.Float(8) / "Longitude of periapsis (deg).", - t = pb.Float(8) / "Time since last periastron passage (s).", - pd = pb.Float(8) / "Orbital period derivative (s/yr).", - wd = pb.Float(8) / "Advance of longitude of periapsis (deg/yr).", - ) - ) / "Barycentric orbital parameters used in folds.", - ) -) -# fmt: on - - -def readpfd(f): - """ - Read in a PRESTO folded data (*.pfd) file. - """ - - with open(f, "rb") as fp: - meta = PFD.parse(fp) - for key, shape in { - "dms": (meta["ndms"],), - "periods": (meta["nperiods"],), - "pdots": (meta["npdots"],), - "profs": (meta["npart"], meta["nsub"], meta["nbin"]), - "stats": (7, meta["nsub"], meta["npart"]), - }.items(): - meta[key] = ( - pb.Array( - pb.Float(8), - count=int(prod(shape)), - ) - .parse(fp) - .reshape(shape) - ) - data = meta.pop("profs", None) - return meta, data - - -def writepfd(meta, data, f): - """ - Write out a PRESTO folded data (*.pfd) file. - """ - - with open(f, "wb+") as fp: - PFD.build(meta, fp) - for key in ["dms", "periods", "pdots"]: - meta[key].tofile(fp) - data.tofile(fp) - meta["stats"].tofile(fp) diff --git a/src/priwo/priwo.C b/src/priwo/priwo.C new file mode 100644 index 0000000..8a84f1f --- /dev/null +++ b/src/priwo/priwo.C @@ -0,0 +1,16 @@ +#include + +#include "priwo.h" + +NB_MODULE(_internals, m) { + init_bits(m); + init_hdr(m); + init_tim(m); + init_fil(m); + init_inf(m); + init_dat(m); + init_pfd(m); + init_fft(m); + init_polycos(m); + init_bestprof(m); +} diff --git a/src/priwo/priwo.h b/src/priwo/priwo.h new file mode 100644 index 0000000..2ce85c9 --- /dev/null +++ b/src/priwo/priwo.h @@ -0,0 +1,15 @@ +#ifndef PRIWO +#define PRIWO + +#include "bestprof.h" +#include "bits.h" +#include "dat.h" +#include "fft.h" +#include "fil.h" +#include "hdr.h" +#include "inf.h" +#include "pfd.h" +#include "polycos.h" +#include "tim.h" + +#endif diff --git a/src/priwo/sigproc/__init__.py b/src/priwo/sigproc/__init__.py deleted file mode 100644 index b136ed1..0000000 --- a/src/priwo/sigproc/__init__.py +++ /dev/null @@ -1,16 +0,0 @@ -""" -R/W pulsar data from formats used by the SIGPROC package. -""" - -from priwo.sigproc.hdr import readhdr, writehdr -from priwo.sigproc.tim import readtim, writetim -from priwo.sigproc.fil import readfil, writefil - -__all__ = [ - "readhdr", - "readtim", - "readfil", - "writehdr", - "writetim", - "writefil", -] diff --git a/src/priwo/sigproc/fil.py b/src/priwo/sigproc/fil.py deleted file mode 100644 index 4db8ad4..0000000 --- a/src/priwo/sigproc/fil.py +++ /dev/null @@ -1,61 +0,0 @@ -""" -R/W SIGPROC filterbank (*.fil) files. -""" - -import pabo as pb - -from priwo.sigproc.hdr import readhdr -from priwo.sigproc.hdr import writehdr - - -def readfil(f): - """ - Read in a SIGPROC filterbank (*.fil) file. - """ - - meta = readhdr(f) - size = meta["size"] - nbits = meta.get("nbits", None) - with open(f, "rb") as fp: - fp.seek(size) - data = pb.Array( - { - 1: pb.Int(1), - 2: pb.Int(1), - 4: pb.Int(1), - 8: pb.Int(1), - 16: pb.Int(2), - 32: pb.Float(4), - 64: pb.Float(8), - None: pb.Float(4), - }[nbits], - packing=(nbits if nbits in [1, 2, 4] else None), - ).parse(fp) - data = data.reshape(-1, meta["nchans"]) - data = data.T - return meta, data - - -def writefil(meta, data, f): - """ - Write out a SIGPROC filterbank (*.fil) file. - """ - - writehdr(meta, f) - size = meta["size"] - nbits = meta.get("nbits", None) - with open(f, "ab") as fp: - fp.seek(size) - pb.Array( - { - 1: pb.Int(1), - 2: pb.Int(1), - 4: pb.Int(1), - 8: pb.Int(1), - 16: pb.Int(2), - 32: pb.Float(4), - 64: pb.Float(8), - None: pb.Float(4), - }[nbits], - packing=(nbits if nbits in [1, 2, 4] else None), - ).build(data.T, fp) diff --git a/src/priwo/sigproc/hdr.py b/src/priwo/sigproc/hdr.py deleted file mode 100644 index d7bc899..0000000 --- a/src/priwo/sigproc/hdr.py +++ /dev/null @@ -1,87 +0,0 @@ -""" -R/W SIGPROC headers. -""" - -import pabo as pb - -# fmt: off -START = pb.Const("HEADER_START", pb.PascalString(pb.Int(4), "utf8")) -END = pb.Const("HEADER_END", pb.PascalString(pb.Int(4), "utf8")) -# fmt: on - -# fmt: off -HDRKEYS = { - "filename": pb.PascalString(pb.Int(4), "utf8"), - "telescope_id": pb.Int(4), - "telescope": pb.PascalString(pb.Int(4), "utf8"), - "machine_id": pb.Int(4), - "data_type": pb.Int(4), - "rawdatafile": pb.PascalString(pb.Int(4), "utf8"), - "source_name": pb.PascalString(pb.Int(4), "utf8"), - "barycentric": pb.Int(4), - "pulsarcentric": pb.Int(4), - "az_start": pb.Float(8), - "za_start": pb.Float(8), - "src_raj": pb.Float(8), - "src_dej": pb.Float(8), - "tstart": pb.Float(8), - "tsamp": pb.Float(8), - "nbits": pb.Int(4), - "nsamples": pb.Int(4), - "fch1": pb.Float(8), - "foff": pb.Float(8), - "fchannel": pb.Float(8), - "nchans": pb.Int(4), - "nifs": pb.Int(4), - "refdm": pb.Float(8), - "flux": pb.Float(8), - "period": pb.Float(8), - "nbeams": pb.Int(4), - "ibeam": pb.Int(4), - "hdrlen": pb.Int(4), - "pb": pb.Float(8), - "ecc": pb.Float(8), - "asini": pb.Float(8), - "orig_hdrlen": pb.Int(4), - "new_hdrlen": pb.Int(4), - "sampsize": pb.Int(4), - "bandwidth": pb.Float(8), - "fbottom": pb.Float(8), - "ftop": pb.Float(8), - "obs_date": pb.PascalString(pb.Int(4), "utf8"), - "obs_time": pb.PascalString(pb.Int(4), "utf8"), - "signed": pb.Int(1), - "accel": pb.Float(8), -} -# fmt: on - - -def readhdr(f): - """ - Read in a SIGPROC header. - """ - - meta = {} - with open(f, "rb") as fp: - START.parse(fp) - while True: - key = pb.PascalString(pb.Int(4), "utf8").parse(fp) - if key == "HEADER_END": - break - meta[key] = HDRKEYS[key].parse(fp) - meta["size"] = fp.tell() - return meta - - -def writehdr(meta, f): - """ - Write out a SIGPROC header. - """ - - with open(f, "wb+") as fp: - START.build(fp) - for key, val in meta.items(): - if key != "size": - pb.PascalString(pb.Int(4), "utf8").build(key, fp) - HDRKEYS[key].build(val, fp) - END.build(fp) diff --git a/src/priwo/sigproc/tim.py b/src/priwo/sigproc/tim.py deleted file mode 100644 index f27f016..0000000 --- a/src/priwo/sigproc/tim.py +++ /dev/null @@ -1,58 +0,0 @@ -""" -R/W SIGPROC time series (*.tim) files. -""" - -import pabo as pb -from priwo.sigproc.hdr import readhdr -from priwo.sigproc.hdr import writehdr - - -def readtim(f): - """ - Read in a SIGPROC time series (*.tim) file. - """ - - meta = readhdr(f) - size = meta["size"] - nbits = meta.get("nbits", 8) - with open(f, "rb") as fp: - fp.seek(size) - data = pb.Array( - { - 1: pb.Int(1), - 2: pb.Int(1), - 4: pb.Int(1), - 8: pb.Int(1), - 16: pb.Int(2), - 32: pb.Float(4), - 64: pb.Float(8), - None: pb.Float(4), - }[nbits], - packing=(nbits if nbits in [1, 2, 4] else None), - ).parse(fp) - return meta, data - - -def writetim(meta, data, f): - """ - Write out a SIGPROC time series (*.tim) file. - """ - - writehdr(meta, f) - size = meta["size"] - nbits = meta.get("nbits", None) - with open(f, "ab") as fp: - fp.seek(size) - pb.Array( - { - 1: pb.Int(1), - 2: pb.Int(1), - 4: pb.Int(1), - 8: pb.Int(1), - 16: pb.Int(2), - 32: pb.Float(4), - 64: pb.Float(8), - None: pb.Float(4), - }[nbits], - packing=(nbits if nbits in [1, 2, 4] else None), - ).build(data, fp) diff --git a/src/priwo/tim.C b/src/priwo/tim.C new file mode 100644 index 0000000..5da69f8 --- /dev/null +++ b/src/priwo/tim.C @@ -0,0 +1,66 @@ +#include "hdr.h" +#include "tim.h" + +using namespace nanobind::literals; + +std::tuple>> +readtim(std::string fn) { + nb::dict meta = readhdr(fn); + int nbits = nb::cast(meta["nbits"]); + int nskip = nb::cast(meta["hdrlen"]); + + const char *dtype; + switch (nbits) { + case 1: + dtype = "uint8"; + break; + case 2: + dtype = "uint8"; + break; + case 4: + dtype = "uint8"; + break; + case 8: + dtype = "uint8"; + break; + case 16: + dtype = "uint16"; + break; + case 32: + dtype = "float32"; + break; + } + + nb::object np = nb::module_::import_("numpy"); + nb::object dobj = + np.attr("fromfile")(fn, "dtype"_a = np.attr(dtype), "offset"_a = nskip); + if (nbits == 1 || nbits == 2 || nbits == 4) { + nb::object priwo = nb::module_::import_("priwo"); + dobj = priwo.attr("_internals") + .attr("unpack")(dobj, "nbits"_a = nbits, "bitorder"_a = "little", + "parallel"_a = true); + } + nb::ndarray> data = + nb::cast>>(dobj); + return std::make_tuple(meta, data); +} + +void writetim(nb::dict meta, nb::ndarray> data, + std::string fn) { + writehdr(meta, fn); + nb::object np = nb::module_::import_("numpy"); + nb::object dobj = np.attr("asarray")(data); + int nbits = nb::cast(meta["nbits"]); + if (nbits == 1 || nbits == 2 || nbits == 4) { + nb::object priwo = nb::module_::import_("priwo"); + dobj = priwo.attr("_internals") + .attr("pack")(dobj, "nbits"_a = nbits, "bitorder"_a = "little", + "parallel"_a = true); + } + dobj.attr("tofile")(fn); +} + +void init_tim(nb::module_ m) { + m.def("readtim", &readtim); + m.def("writetim", &writetim); +} diff --git a/src/priwo/tim.h b/src/priwo/tim.h new file mode 100644 index 0000000..c9db6fb --- /dev/null +++ b/src/priwo/tim.h @@ -0,0 +1,8 @@ +#ifndef TIM +#define TIM + +#include "common.h" + +void init_tim(nb::module_ m); + +#endif diff --git a/src/priwo/timing/__init__.py b/src/priwo/timing/__init__.py deleted file mode 100644 index 1e480f5..0000000 --- a/src/priwo/timing/__init__.py +++ /dev/null @@ -1,10 +0,0 @@ -""" -R/W pulsar data from formats used by timing packages, like TEMPO/TEMPO2. -""" - -from priwo.timing.polycos import readpolycos, writepolycos - -__all__ = [ - "readpolycos", - "writepolycos", -] diff --git a/src/priwo/timing/polycos.py b/src/priwo/timing/polycos.py deleted file mode 100644 index 3b67889..0000000 --- a/src/priwo/timing/polycos.py +++ /dev/null @@ -1,133 +0,0 @@ -""" -R/W TEMPO polynomial ephemerides (*.polycos) files. -""" - -import re -import pabo as pb - -one = lambda _: _[0] -pad = lambda _, N: (_ + [None] * N)[:N] -fx = lambda _, __: __(_) if _ is not None else _ - -# fmt: off -TEMPLATES = { - "name": "{:<10s}", - "date": "{:>10s}", - "utc": "{:>11s}", - "tmid": "{:>20.11f}", - "dm": "{:>21.6f}", - "doppler": "{:>7.3f}", - "rms": "{:>7.3f}\n", - "refphz": "{:>20.6f}", - "reff0": "{:^20.12f}", - "obsno": "{:^5d}", - "span": "{:^5.0f}", - "ncoeff": "{:^4d}", - "freq": "{:^21.3f}", - "binphz": "{:^4.0f}\n", -} -# fmt: on - - -def readpolycos(f): - """ - Read in a TEMPO polynomial ephemerides (*.polycos) file. - """ - - polycos = [] - with open(f, "r") as fp: - while True: - line = fp.readline() - if line == "": - break - - ( - utc, - tmid, - dm, - doppler, - rms, - ) = pad(re.findall(r"[+-]?[0-9]+[.][0-9]+", line), 5) - name = one(re.findall(r"[A-Z]?[0-9]{4}[+-][0-9]{2,4}", line)) - date = one(re.findall(r"[0-9]{2}[-][A-Za-z]{3}[-][0-9]{2}", line)) - - line = fp.readline() - ( - refphz, - reff0, - obsno, - span, - ncoeff, - freq, - binphz, - ) = pad(re.findall(r"[0-9]+[.]*[0-9]*", line), 7) - - polycos.append( - { - "meta": { - # fmt: off - "name": str(name), - "date": str(date), - "utc": str(utc), - "tmid": fx(tmid, float), - "dm": fx(dm, float), - "doppler": fx(doppler, float), - "rms": fx(rms, float), - "refphz": fx(refphz, float), - "reff0": fx(reff0, float), - "obsno": fx(obsno, int), - "span": fx(span, float), - "ncoeff": fx(ncoeff, int), - "freq": fx(freq, float), - "binphz": fx(binphz, float), - # fmt: on - }, - "data": pb.Array(pb.Float(8)).parse( - b"".join( - [ - __ - for _ in [ - [ - pb.Float(8).build( - float(re.sub(r"[dD]", "e", _)) if _ else 0.0 - ) - for _ in re.findall( - r"[+-]?(?:0|[1-9]\d*)(?:\.\d*)?(?:[eEdD][+\-]?\d+)?", - fp.readline(), - ) - ] - for _ in range(int(ncoeff) // 3) - ] - for __ in _ - ] - ) - ), - } - ) - return polycos - - -def writepolycos(polycos, f): - """ - Write out a TEMPO polynomial ephemerides (*.polycos) file. - """ - - with open(f, "w+") as fp: - for polyco in polycos: - fp.write( - "".join( - list( - [ - TEMPLATES[_].format(__) - if __ is not None - else TEMPLATES[_].replace("f", "s").format("") - for _, __ in polyco["meta"].items() - ] - ) - ) - ) - fp.write( - ( - "{:>25.16E}{:>25.16E}{:>25.16E}\n" * (polyco["meta"]["ncoeff"] // 3) - ).format(*list(polyco["data"])) - ) diff --git a/tests/data/test.bestprof b/tests/data/test.bestprof deleted file mode 100644 index 030fa8c..0000000 --- a/tests/data/test.bestprof +++ /dev/null @@ -1,155 +0,0 @@ -# Input file = J1646-2142_500_200_512_2.15jul2k19.raw0.fil -# Candidate = PSR_1646-2142 -# Telescope = GMRT -# Epoch_topo = 58679.685665396297 -# Epoch_bary = N/A -# T_sample = 1.024e-05 -# Data Folded = 263664000 -# Data Avg = 11042.9913182359 -# Data StdDev = 246.537468521724 -# Profile Bins = 128 -# Profile Avg = 22747152667.3743 -# Profile StdDev = 353837.098935206 -# Reduced chi-sqr = 106.001 -# Prob(Noise) < 0 (~112.9 sigma) -# Best DM = 29.741 -# P_topo (ms) = 5.85347537656615 +/- 5.16e-08 -# P'_topo (s/s) = -0 +/- 1.48e-13 -# P''_topo (s/s^2) = 0 +/- 3.55e-16 -# P_bary (ms) = N/A -# P'_bary (s/s) = N/A -# P''_bary (s/s^2) = N/A -# P_orb (s) = N/A -# asin(i)/c (s) = N/A -# eccentricity = N/A -# w (rad) = N/A -# T_peri = N/A -###################################################### - 0 2.274448e+10 - 1 2.274387e+10 - 2 2.274407e+10 - 3 2.274434e+10 - 4 2.274426e+10 - 5 2.274443e+10 - 6 2.274413e+10 - 7 2.274457e+10 - 8 2.274456e+10 - 9 2.274414e+10 - 10 2.274449e+10 - 11 2.274425e+10 - 12 2.274422e+10 - 13 2.274493e+10 - 14 2.274423e+10 - 15 2.274431e+10 - 16 2.274468e+10 - 17 2.274435e+10 - 18 2.274494e+10 - 19 2.274472e+10 - 20 2.274447e+10 - 21 2.274442e+10 - 22 2.274458e+10 - 23 2.274484e+10 - 24 2.274548e+10 - 25 2.274551e+10 - 26 2.274534e+10 - 27 2.274512e+10 - 28 2.274574e+10 - 29 2.274579e+10 - 30 2.274511e+10 - 31 2.274561e+10 - 32 2.274531e+10 - 33 2.274449e+10 - 34 2.274533e+10 - 35 2.274572e+10 - 36 2.27453e+10 - 37 2.274533e+10 - 38 2.274549e+10 - 39 2.274473e+10 - 40 2.274471e+10 - 41 2.274452e+10 - 42 2.274548e+10 - 43 2.274569e+10 - 44 2.274558e+10 - 45 2.274558e+10 - 46 2.274629e+10 - 47 2.274659e+10 - 48 2.274679e+10 - 49 2.274627e+10 - 50 2.274669e+10 - 51 2.27486e+10 - 52 2.274862e+10 - 53 2.274957e+10 - 54 2.275093e+10 - 55 2.27513e+10 - 56 2.275209e+10 - 57 2.275344e+10 - 58 2.275462e+10 - 59 2.275635e+10 - 60 2.275708e+10 - 61 2.275796e+10 - 62 2.275842e+10 - 63 2.275846e+10 - 64 2.275808e+10 - 65 2.275682e+10 - 66 2.275543e+10 - 67 2.275422e+10 - 68 2.275298e+10 - 69 2.27521e+10 - 70 2.275152e+10 - 71 2.274946e+10 - 72 2.274909e+10 - 73 2.274845e+10 - 74 2.274771e+10 - 75 2.274703e+10 - 76 2.274676e+10 - 77 2.274587e+10 - 78 2.27462e+10 - 79 2.274559e+10 - 80 2.274462e+10 - 81 2.27452e+10 - 82 2.274479e+10 - 83 2.2745e+10 - 84 2.274572e+10 - 85 2.274506e+10 - 86 2.274566e+10 - 87 2.274544e+10 - 88 2.274619e+10 - 89 2.274546e+10 - 90 2.274595e+10 - 91 2.274653e+10 - 92 2.274656e+10 - 93 2.274696e+10 - 94 2.274689e+10 - 95 2.274802e+10 - 96 2.274805e+10 - 97 2.274831e+10 - 98 2.274924e+10 - 99 2.274956e+10 - 100 2.27497e+10 - 101 2.275046e+10 - 102 2.275078e+10 - 103 2.274968e+10 - 104 2.274894e+10 - 105 2.274943e+10 - 106 2.274789e+10 - 107 2.274758e+10 - 108 2.274711e+10 - 109 2.274681e+10 - 110 2.274625e+10 - 111 2.274629e+10 - 112 2.274546e+10 - 113 2.274607e+10 - 114 2.27461e+10 - 115 2.27457e+10 - 116 2.274588e+10 - 117 2.274594e+10 - 118 2.274588e+10 - 119 2.274561e+10 - 120 2.274525e+10 - 121 2.274523e+10 - 122 2.274527e+10 - 123 2.274486e+10 - 124 2.27448e+10 - 125 2.274452e+10 - 126 2.274426e+10 - 127 2.274392e+10 diff --git a/tests/data/test.dat b/tests/data/test.dat deleted file mode 100644 index eb6a597..0000000 Binary files a/tests/data/test.dat and /dev/null differ diff --git a/tests/data/test.fft b/tests/data/test.fft deleted file mode 100644 index f86ee42..0000000 Binary files a/tests/data/test.fft and /dev/null differ diff --git a/tests/data/test.fil b/tests/data/test.fil deleted file mode 100644 index d81cd60..0000000 Binary files a/tests/data/test.fil and /dev/null differ diff --git a/tests/data/test.inf b/tests/data/test.inf deleted file mode 100644 index fc2f57b..0000000 --- a/tests/data/test.inf +++ /dev/null @@ -1,22 +0,0 @@ - Data file name without suffix = fake_presto_radio - Telescope used = Parkes - Instrument used = Multibeam - Object being observed = Pulsar - J2000 Right Ascension (hh:mm:ss.ssss) = 00:00:01.0000 - J2000 Declination (dd:mm:ss.ssss) = -00:00:01.0000 - Data observed by = Kenji Oba - Epoch of observation (MJD) = 59000.0 - Barycentered? (1 yes, 0 no) = 1 - Number of bins in the time series = 16 - Width of each time series bin (sec) = 6.4e-05 - Any breaks in the data? (1 yes, 0 no) = 0 - Type of observation (EM band) = Radio - Beam diameter (arcsec) = 981.0 - Dispersion measure (cm-3 pc) = 42.42 - Central freq of low channel (MHz) = 1182.1953125 - Total bandwidth (MHz) = 400.0 - Number of channels = 1024 - Channel bandwidth (MHz) = 0.390625 - Data analyzed by = Space Sheriff Gavan - Any additional notes: - Input filterbank samples have 2 bits. diff --git a/tests/data/test.pfd b/tests/data/test.pfd deleted file mode 100644 index 23615f6..0000000 Binary files a/tests/data/test.pfd and /dev/null differ diff --git a/tests/data/test.polycos b/tests/data/test.polycos deleted file mode 100644 index 2078c3e..0000000 --- a/tests/data/test.polycos +++ /dev/null @@ -1,432 +0,0 @@ -1646-2142 14-JUL-19 50000.00 58678.20833333330 29.741000 0.605 -6.064 - 30728623742.970901 170.849389109675 27 60 12 399.805 105 - 1.0969146311545774E-06 -6.2056431791115385E-01 2.5599064429269767E-05 - -5.2444212654562188E-09 -1.1805989251899655E-10 -9.4231511657727734E-12 - 1.7097184257639180E-13 1.9622534107682035E-14 -1.5905644515708209E-16 - -1.7482189238810964E-17 5.0468890900483809E-20 5.5993354653641204E-21 -1646-2142 14-JUL-19 60000.00 58678.25000000000 29.741000 0.602 -6.207 - 30729238763.628349 170.849389109675 27 60 12 399.805 - -2.8087855264737061E-07 -6.1760905131639865E-01 2.3319080073267456E-05 - -1.8093317152862359E-08 5.0668833920437720E-11 -2.5987694977307418E-12 - -2.7661293806924571E-13 8.8353139554091779E-15 2.9127483414166014E-16 - -9.7164823275297494E-18 -1.0266508882988786E-19 3.4568299813080993E-21 -1646-2142 14-JUL-19 70000.00 58678.29166666660 29.741000 0.600 -5.982 - 30729853784.451115 170.849389109675 27 60 12 399.805 - -6.3646801010987242E-07 -6.1504287438969962E-01 1.9134551640564395E-05 - -2.7652914152694469E-08 8.5187242034532451E-12 3.1600277125515597E-13 - -6.3754722837317780E-14 -3.4550944766961893E-15 3.7584935373312459E-17 - 4.6126764424497782E-18 -8.2798113962437442E-21 -1.7853959635168298E-21 -1646-2142 14-JUL-19 80000.00 58678.33333333330 29.741000 0.598 -6.158 - 30730468805.412231 170.849389109675 27 60 12 399.805 - -2.4185076181147680E-07 -6.1307910016514378E-01 1.3363234058091192E-05 - -3.5177994655782725E-08 -1.9152467556379766E-10 -3.4421190509884995E-12 - 4.0282176390940081E-13 5.0634753504107773E-15 -3.9594157759802094E-16 - -3.1534277521980084E-18 1.3406380875284955E-19 7.1388351344393414E-22 -1646-2142 14-JUL-19 90000.00 58678.37500000000 29.741000 0.597 -5.948 - 30731083826.468670 170.849389109675 27 60 12 399.805 - 1.5825181348865305E-06 -6.1188993702291949E-01 6.3167725782335724E-06 - -4.2573130872987786E-08 -4.4211342127064520E-11 4.9401621736761082E-12 - 3.8458709659199675E-14 -9.7677492869785195E-15 -2.4116266865904183E-17 - 8.3258759980492771E-18 5.2527643852263944E-21 -2.5453463184816893E-21 -1646-2142 14-JUL-19 100000.00 58678.41666666660 29.741000 0.597 -6.174 - 30731698847.568615 170.849389109675 27 60 12 399.805 - -4.4598444782378323E-07 -6.1159433316503886E-01 -1.4815634999753802E-06 - -4.9118994882870682E-08 3.6773782693362188E-11 1.7714149672723568E-11 - -4.8449096209430545E-14 -2.7759747029111250E-14 1.8013384463837715E-17 - 1.9793791743638711E-17 5.2905569646101498E-22 -5.2886875446160967E-21 -1646-2142 14-JUL-19 110000.00 58678.45833333330 29.741000 0.597 -5.952 - 30732313868.659325 170.849389109675 27 60 12 399.805 - 1.0703549399929124E-06 -6.1225117739537349E-01 -9.4465953390040577E-06 - -4.3156011582528771E-08 -8.1874752496277556E-11 -2.8952885325780127E-12 - 2.0328286073177856E-13 3.7603541884725031E-15 -1.8587687164428594E-16 - -1.5976138571268230E-18 5.9432557357071224E-20 4.7680615624511566E-23 -1646-2142 14-JUL-19 120000.00 58678.50000000000 29.741000 0.599 -6.017 - 30732928889.682125 170.849389109675 27 60 12 399.805 - 8.5202900608656595E-07 -6.1385172660733200E-01 -1.7091472852749631E-05 - -4.1069914992588911E-08 -1.2648022118856420E-10 1.8455599967817549E-12 - 3.0463292265255786E-13 -2.5105698667924790E-15 -2.6494863528718959E-16 - 1.2619673571505250E-18 8.2154775574142424E-20 -1.5092375132543554E-22 -1646-2142 14-JUL-19 130000.00 58678.54166666660 29.741000 0.601 -6.217 - 30733543910.580730 170.849389109675 27 60 12 399.805 - -6.1404338360521171E-08 -6.1632451070439853E-01 -2.3917305463225909E-05 - -3.5288104052355634E-08 5.0817723140936974E-11 2.0596963974551649E-12 - -3.7388477797430177E-14 -1.4929966591241536E-15 2.2177983505410851E-17 - 2.0708360956819916E-20 -2.1596255930819062E-21 2.0257452852144915E-22 -1646-2142 14-JUL-19 140000.00 58678.58333333330 29.741000 0.604 -6.149 - 30734158931.309437 170.849389109675 27 60 12 399.805 - -1.0689240059411076E-06 -6.1953764695932312E-01 -2.9406441458428863E-05 - -2.6968174368971803E-08 1.7943983265451856E-10 6.7317881889736628E-12 - -3.0618771241707744E-13 -1.5680655213756161E-14 2.7756937702746606E-16 - 1.4725113407629362E-17 -8.8651999787310088E-20 -4.7827602166092632E-21 -1646-2142 14-JUL-19 150000.00 58678.62500000000 29.741000 0.608 -6.000 - 30734773951.827602 170.849389109675 27 60 12 399.805 - 8.1718101865163595E-07 -6.2330834371944643E-01 -3.3151528203209410E-05 - -1.5895556194135656E-08 1.5338878841442409E-11 -2.7745434068210876E-13 - 9.5052475406750477E-14 1.2569363870798863E-15 -1.0179959674538156E-16 - -1.0859451194312884E-18 3.5753765490197932E-20 2.6963262714516515E-22 -1646-2142 14-JUL-19 160000.00 58678.66666666660 29.741000 0.612 -6.077 - 30735388972.106770 170.849389109675 27 60 12 399.805 - 9.3423442894258191E-07 -6.2741595638853320E-01 -3.4956813316572898E-05 - -3.2217146044164367E-09 -1.2032388855963460E-11 -6.4427792729053801E-12 - 1.3483765094387229E-13 1.3383842574005559E-14 -1.1666724815950412E-16 - -1.1679698652282618E-17 3.5322617724753012E-20 3.6241119093596227E-21 -1646-2142 14-JUL-19 170000.00 58678.70833333330 29.741000 0.616 -6.102 - 30736003992.136925 170.849389109675 27 60 12 399.805 - -3.2869193531889773E-07 -6.3161635887887158E-01 -3.4677656910797447E-05 - 7.8243739799510606E-09 -1.3076888106425182E-11 -3.6411707447291845E-12 - 9.0772402811910030E-14 9.4269600828738680E-15 -5.5583290286929290E-17 - -9.6844573699587460E-18 1.2229683053393578E-20 3.3989763161936253E-21 -1646-2142 14-JUL-19 180000.00 58678.75000000000 29.741000 0.620 -6.092 - 30736619011.918541 170.849389109675 27 60 12 399.805 - 8.6098025706832226E-07 -6.3565903176172311E-01 -3.2372637226936881E-05 - 1.5148972785578596E-08 1.3737355898266936E-10 1.5593009302295400E-11 - -1.6425513420765839E-13 -2.8924673409317005E-14 1.1938417726562335E-16 - 2.3160855085994979E-17 -3.0510794403797661E-20 -6.7065760245628923E-21 -1646-2142 14-JUL-19 190000.00 58678.79166666660 29.741000 0.624 -6.045 - 30737234031.466877 170.849389109675 27 60 12 399.805 - -1.7820899061646300E-07 -6.3930502566632941E-01 -2.8122309784185977E-05 - 2.8234849275293215E-08 7.4594668670167028E-11 2.8140451547642094E-13 - -5.3051448109788640E-14 4.0921750379183726E-16 3.6077729610532735E-17 - -1.1444101209181174E-18 -9.5877422471948200E-21 5.6009420997025799E-22 -1646-2142 14-JUL-19 200000.00 58678.83333333330 29.741000 0.627 -6.028 - 30737849050.815254 170.849389109675 27 60 12 399.805 - -1.1763476604267908E-07 -6.4234142697750896E-01 -2.2257681785471989E-05 - 3.6873526154940641E-08 1.8302142422098033E-11 -2.8867043717720552E-12 - 4.4301788989936045E-14 4.7985114669292261E-15 -5.4770226251011033E-17 - -3.0365056631021625E-18 2.1608173607564594E-20 5.9285545132269650E-22 -1646-2142 14-JUL-19 210000.00 58678.87500000000 29.741000 0.629 -6.010 - 30738464070.004124 170.849389109675 27 60 12 399.805 - -1.6715551089240693E-06 -6.4459778068667029E-01 -1.5178901956255549E-05 - 4.4252563501822627E-08 -3.6128101316231095E-11 -1.1338614426959927E-11 - 1.4284488492019625E-13 2.0563988741792526E-14 -1.4844561750022562E-16 - -1.6419932965838413E-17 5.2225630790687566E-20 4.7922395736179120E-21 -1646-2142 14-JUL-19 220000.00 58678.91666666660 29.741000 0.630 -6.024 - 30739079089.082619 170.849389109675 27 60 12 399.805 - 2.6861433244027011E-07 -6.4595668181621779E-01 -7.3984649796187242E-06 - 4.5112042832459984E-08 7.3096047705520992E-11 -3.0586387185212590E-12 - -1.2103868846109613E-13 4.6790913716636870E-15 8.2636219463303044E-17 - -3.1654539958395194E-18 -1.8552403803988747E-20 7.8998635501808997E-22 -1646-2142 14-JUL-19 230000.00 58678.95833333330 29.741000 0.631 -6.105 - 30739694108.109528 170.849389109675 27 60 12 399.805 - -1.1725700973340002E-06 -6.4636273571093694E-01 5.9856134329925191E-07 - 4.5950066734878591E-08 5.2811872237971239E-11 -8.7021171320962416E-12 - -1.2605669537837070E-13 1.5617737893223370E-14 9.9586707628047775E-17 - -1.2393819187027964E-17 -2.6670730316596304E-20 3.5698166647751534E-21 -1646-2142 15-JUL-19 0.00 58679.00000000000 29.741000 0.630 -6.042 - 30740309127.140594 170.849389109675 27 60 12 399.805 - 1.1561395545610365E-07 -6.4582502991186741E-01 8.2465906311509569E-06 - 3.9987542248782133E-08 -3.4267189189964006E-12 4.2283434344120460E-12 - -4.6546991324303836E-14 -1.1002425190238011E-14 4.0112544646256599E-17 - 1.1016177443650935E-17 -1.0106493077133440E-20 -3.7881353040332855E-21 -1646-2142 15-JUL-19 10000.00 58679.04166666660 29.741000 0.629 -6.037 - 30740924146.229046 170.849389109675 27 60 12 399.805 - -2.2384346281753564E-07 -6.4441805445767875E-01 1.5019362143360889E-05 - 3.3315100227316237E-08 -6.9994284339563690E-11 2.2358474710833351E-12 - 8.9287510390712164E-14 -1.6699347608166772E-15 -8.4059255661716572E-17 - -9.3027924782677104E-20 2.7633462262995188E-20 3.3989448691382337E-22 -1646-2142 15-JUL-19 20000.00 58679.08333333330 29.741000 0.627 -6.110 - 30741539165.426262 170.849389109675 27 60 12 399.805 - -4.6767617471743517E-07 -6.4227543588471425E-01 2.0416523290792692E-05 - 2.4854384898672157E-08 7.4303894548082820E-11 -2.0059147888935810E-12 - -2.2612745261569341E-13 1.0280172412908243E-14 1.9350084286515371E-16 - -1.2018619218804160E-17 -6.0547888862445772E-20 4.3333151086384915E-21 -1646-2142 15-JUL-19 30000.00 58679.12500000000 29.741000 0.624 -6.075 - 30742154184.769600 170.849389109675 27 60 12 399.805 - -7.0448226922176227E-08 -6.3958154416481805E-01 2.4154877213604649E-05 - 1.5718327200207438E-08 -9.4710431563989408E-11 -2.7289712607314541E-12 - 1.4536040399675340E-13 6.0012342278203575E-15 -1.5234209871301864E-16 - -5.4700459663220890E-18 5.2655295626672484E-20 1.7681077479066514E-21 -1646-2142 15-JUL-19 40000.00 58679.16666666660 29.741000 0.621 -6.021 - 30742769204.284176 170.849389109675 27 60 12 399.805 - -8.4852689519882007E-07 -6.3655776231071948E-01 2.5909933296974001E-05 - 2.4673519803688865E-09 -6.4418751010702391E-11 4.7041776131898853E-12 - -1.0081125893339214E-14 -5.8298472451191286E-15 3.4001317013506264E-17 - 2.9743766084030928E-18 -1.5810207343290627E-20 -4.7081339669507842E-22 -1646-2142 15-JUL-19 50000.00 58679.20833333330 29.741000 0.618 -6.043 - 30743384223.985458 170.849389109675 27 60 12 399.805 - 5.9526191758044189E-07 -6.3344896152002450E-01 2.5575160715323047E-05 - -8.1877257781858862E-09 -1.3060849634442140E-10 3.5291260024959481E-12 - 2.2235259950867026E-13 -7.1860577999829419E-15 -2.3383531752783705E-16 - 5.8837294606583341E-18 8.3416298492658412E-20 -1.7055124862508472E-21 -1646-2142 15-JUL-19 60000.00 58679.25000000000 29.741000 0.615 -6.220 - 30743999243.869568 170.849389109675 27 60 12 399.805 - -5.5595578056684414E-07 -6.3050512931848524E-01 2.3188370730168483E-05 - -1.8565378804897392E-08 -2.1488337873404544E-10 -1.0500936131423890E-12 - 4.0816047797414134E-13 3.7428703755669062E-15 -3.8285965529265905E-16 - -4.6533424165409427E-18 1.2387202753481428E-19 1.8257501615435813E-21 -1646-2142 15-JUL-19 70000.00 58679.29166666660 29.741000 0.613 -6.057 - 30744614263.917854 170.849389109675 27 60 12 399.805 - 1.6972896525382225E-06 -6.2796510000474837E-01 1.8847492153392598E-05 - -2.8342387864551231E-08 8.4578706043256088E-11 -1.8191059215516956E-12 - -2.3731527296482827E-13 3.8604200581267122E-15 1.9782881856132253E-16 - -3.3914583531405169E-18 -5.9499314282558531E-20 1.0747289908382919E-21 -1646-2142 15-JUL-19 80000.00 58679.33333333330 29.741000 0.611 -6.135 - 30745229284.102524 170.849389109675 27 60 12 399.805 - 6.1967347714762810E-07 -6.2604037289539061E-01 1.2988937717047779E-05 - -3.5720986598054748E-08 -1.4769765701127771E-10 -2.8536939668678808E-12 - 2.4230719617577424E-13 4.6319975358494801E-15 -1.9243157369260644E-16 - -3.3746862061107209E-18 5.2968845316191440E-20 9.0803225125831579E-22 -1646-2142 15-JUL-19 90000.00 58679.37500000000 29.741000 0.610 -6.033 - 30745844304.379871 170.849389109675 27 60 12 399.805 - 3.0221575424279820E-07 -6.2489974822544236E-01 5.8515827755317612E-06 - -4.1799753347325431E-08 1.3129018069085685E-10 1.7507814174204184E-13 - -3.2601253594303514E-13 -1.4709151861217668E-15 2.9686508241493232E-16 - 1.8430972925684359E-18 -9.5600979074072753E-20 -6.4396664265323035E-22 -1646-2142 15-JUL-19 100000.00 58679.41666666660 29.741000 0.609 -6.071 - 30746459324.697617 170.849389109675 27 60 12 399.805 - 3.5842390847841410E-09 -6.2465861088235031E-01 -1.9447506189610383E-06 - -4.5537396633899883E-08 8.4402844313115586E-11 2.3411026728598771E-12 - -2.2687426537162185E-13 -1.7697226233463004E-15 2.2343954617269333E-16 - 4.3851737543104691E-19 -7.5049030674236891E-20 2.3203188914203335E-23 -1646-2142 15-JUL-19 110000.00 58679.45833333330 29.741000 0.610 -5.985 - 30747074345.002808 170.849389109675 27 60 12 399.805 - 3.0088363055409684E-08 -6.2537070037116993E-01 -9.9150372856126729E-06 - -4.7236011041145460E-08 1.8558228479063833E-11 1.4873601564273337E-11 - -5.0166033979286580E-14 -2.7304121781741616E-14 4.9653197272599419E-17 - 2.2196818118810217E-17 -1.4443955392612225E-20 -6.5839038814086039E-21 -1646-2142 15-JUL-19 120000.00 58679.50000000000 29.741000 0.612 -6.184 - 30747689365.236801 170.849389109675 27 60 12 399.805 - 4.7517908839320167E-07 -6.2702449156620677E-01 -1.7540476335739784E-05 - -4.2449017899566646E-08 6.2473147045976901E-11 6.7670359054931419E-12 - -8.9762836902750731E-14 -7.5389886375116670E-15 7.5359456595533876E-17 - 2.8423926409327064E-18 -2.2221979841519708E-20 -4.7833779023161597E-23 -1646-2142 15-JUL-19 130000.00 58679.54166666660 29.741000 0.614 -6.091 - 30748304385.343605 170.849389109675 27 60 12 399.805 - -7.8323836993215989E-07 -6.2954374118676393E-01 -2.4261089442806999E-05 - -3.0927175139159330E-08 7.7329659006414886E-11 -1.3720580721610549E-11 - -1.4655068360211865E-13 2.5021559286827077E-14 1.6072647828023895E-16 - -2.0034308944548802E-17 -5.7604619793440438E-20 5.8697552145422275E-21 -1646-2142 15-JUL-19 140000.00 58679.58333333330 29.741000 0.617 -6.043 - 30748919405.278038 170.849389109675 27 60 12 399.805 - 8.1581445450160106E-07 -6.3279217715869496E-01 -2.9607143247795905E-05 - -2.7324585274587033E-08 -3.7895278576826918E-11 9.8682581086878423E-12 - 1.0380248247241115E-13 -1.9655322283509931E-14 -5.0775176997688925E-17 - 1.6900499700660012E-17 5.5294506129072578E-21 -5.2117932149830645E-21 -1646-2142 15-JUL-19 150000.00 58679.62500000000 29.741000 0.621 -6.075 - 30749534425.000195 170.849389109675 27 60 12 399.805 - -5.1605449525107816E-07 -6.3658536649123598E-01 -3.3265340980215377E-05 - -1.3356573718694338E-08 2.9496417948605803E-11 -5.4762638863660813E-12 - -1.4857090582486464E-14 8.8632830012836473E-15 3.8268433070384345E-17 - -6.6176178615918226E-18 -1.7030454810429034E-20 1.8408473978910805E-21 -1646-2142 15-JUL-19 160000.00 58679.66666666660 29.741000 0.625 -5.963 - 30750149444.482475 170.849389109675 27 60 12 399.805 - 5.4862682769690752E-07 -6.4069935821353985E-01 -3.4965201883191196E-05 - -3.1089064156866240E-09 9.7239498066715476E-11 -2.9528482117747812E-12 - -7.1688662022194718E-14 5.1786968104831277E-15 4.8191461282936196E-17 - -3.5602573801683093E-18 -1.2643354858035722E-20 8.5341281709338170E-22 -1646-2142 15-JUL-19 170000.00 58679.70833333330 29.741000 0.629 -6.050 - 30750764463.715832 170.849389109675 27 60 12 399.805 - 7.9875404446899667E-07 -6.4489032158328463E-01 -3.4557092684021226E-05 - 1.2139507887706574E-08 1.0438125668588218E-10 -1.7591273022507758E-11 - -8.8121809321321377E-14 3.1104202812735555E-14 6.1228055109168936E-17 - -2.4136965965464953E-17 -1.5826482428174043E-20 6.8356974468152340E-21 -1646-2142 15-JUL-19 180000.00 58679.75000000000 29.741000 0.633 -5.983 - 30751379482.701725 170.849389109675 27 60 12 399.805 - -1.6451466946024174E-06 -6.4890718038899831E-01 -3.2072033902048975E-05 - 2.0927942413448768E-08 6.7769446709332778E-11 -5.9182377261756162E-12 - -8.9841277170512237E-14 8.6762727501473067E-15 1.0039262346499895E-16 - -5.8660874444495766E-18 -3.5754594719099026E-20 1.5018085353736707E-21 -1646-2142 15-JUL-19 190000.00 58679.79166666660 29.741000 0.637 -6.121 - 30751994501.456345 170.849389109675 27 60 12 399.805 - 1.9072088087747352E-06 -6.5251186439446351E-01 -2.7723384781786524E-05 - 2.6566311503099836E-08 8.8352683032170826E-11 1.2125117071210648E-11 - -1.1781346249445668E-13 -2.3391629840773565E-14 1.0974810849043404E-16 - 1.9064381031894836E-17 -3.5472943923682681E-20 -5.5639223372963049E-21 -1646-2142 15-JUL-19 200000.00 58679.83333333330 29.741000 0.639 -6.014 - 30752609520.013851 170.849389109675 27 60 12 399.805 - 3.7318732362331963E-07 -6.5549471323232489E-01 -2.1742872525900051E-05 - 3.5041760915904177E-08 -5.6018373103474364E-11 5.3602824371016253E-12 - 1.1592332368132678E-13 -8.2892806815476240E-15 -6.9280150486045971E-17 - 6.1826623416543919E-18 1.5027070840526152E-20 -1.7541036565178784E-21 -1646-2142 15-JUL-19 210000.00 58679.87500000000 29.741000 0.642 -6.084 - 30753224538.415359 170.849389109675 27 60 12 399.805 - -1.4016768135564212E-06 -6.5768824743319876E-01 -1.4634682789258276E-05 - 4.0533719528555544E-08 2.2332019216523638E-11 6.7225499748555130E-12 - 1.4548334185511764E-14 -1.1928015306950518E-14 -3.9366741303709903E-17 - 9.2022402699955788E-18 2.0640314072439618E-20 -2.5998224978903495E-21 -1646-2142 15-JUL-19 220000.00 58679.91666666660 29.741000 0.643 -6.060 - 30753839556.710453 170.849389109675 27 60 12 399.805 - -1.9480869642950205E-07 -6.5897885889874253E-01 -6.8085169524852438E-06 - 4.5419332181347158E-08 2.9029758188099180E-11 -3.2409669793651221E-12 - -6.2550168061420550E-14 3.5511998335753297E-15 6.6410053668656195E-17 - -1.3069803580650532E-18 -2.4717045754609607E-20 1.3788675655980611E-23 -1646-2142 15-JUL-19 230000.00 58679.95833333330 29.741000 0.643 -6.064 - 30754454574.958130 170.849389109675 27 60 12 399.805 - -1.1911752430878320E-06 -6.5931489009468680E-01 1.1816511831003551E-06 - 4.3451001684556695E-08 1.6193720441171299E-11 2.1290383236375384E-12 - -8.6317790457236843E-14 -5.1652373467676290E-15 1.0883307760633830E-16 - 4.9561230256740472E-18 -4.4340479344018046E-20 -1.6448938115641944E-21 -1646-2142 16-JUL-19 0.00 58680.00000000000 29.741000 0.643 -6.037 - 30755069593.214092 170.849389109675 27 60 12 399.805 - -1.0833185315976129E-06 -6.5871041350103221E-01 8.7934686033785283E-06 - 3.9442406281038127E-08 -1.2574167348259260E-10 3.5629232503566911E-12 - 2.8056627312759339E-13 -6.5706350960709291E-15 -2.9332826692952177E-16 - 5.1684293158946724E-18 1.0471700268525053E-19 -1.4929547501454874E-21 -1646-2142 16-JUL-19 10000.00 58680.04166666660 29.741000 0.641 -6.133 - 30755684611.531261 170.849389109675 27 60 12 399.805 - 3.4507793600195569E-07 -6.5724360589562647E-01 1.5479734882597055E-05 - 3.2994019799483754E-08 -5.1281658639106370E-11 4.4900873639738416E-12 - 1.9490988369787993E-14 -9.5086952474999516E-15 -6.5409616933633422E-18 - 8.5482703636762528E-18 -2.3465157155420284E-22 -2.7791555739963690E-21 -1646-2142 16-JUL-19 20000.00 58680.08333333330 29.741000 0.639 -5.979 - 30756299629.960487 170.849389109675 27 60 12 399.805 - 6.0478298509816786E-07 -6.5505183740421047E-01 2.0812438055262596E-05 - 2.6157161724641797E-08 -1.2717248260121116E-10 -4.4265481126858104E-12 - 1.6764373403175783E-13 6.4571612281403945E-15 -1.3763305581233023E-16 - -4.0062595709497497E-18 4.0724090989232641E-20 8.7850635180967763E-22 -1646-2142 16-JUL-19 30000.00 58680.12500000000 29.741000 0.636 -6.079 - 30756914648.538403 170.849389109675 27 60 12 399.805 - 8.3072348484002079E-07 -6.5232183210919492E-01 2.4388392349982010E-05 - 1.1674197919692338E-08 -6.2508213277650651E-11 1.3105827270785112E-11 - 3.0518365330842827E-14 -2.3861358296319806E-14 -1.8882327841933992E-17 - 1.9046122574484975E-17 3.6061863548484149E-21 -5.5484899622698965E-21 -1646-2142 16-JUL-19 40000.00 58680.16666666660 29.741000 0.633 -6.165 - 30757529667.289242 170.849389109675 27 60 12 399.805 - -4.0503936789629269E-07 -6.4927831013105897E-01 2.6007055172338510E-05 - 5.8356933802628287E-09 -1.0268994931543492E-10 -8.9020925749453996E-12 - 1.4763253102851267E-13 1.3073762864513167E-14 -1.5812779367592048E-16 - -8.9561731134779658E-18 5.7986924053711889E-20 2.3546864364990219E-21 -1646-2142 16-JUL-19 50000.00 58680.20833333330 29.741000 0.630 -6.098 - 30758144686.227512 170.849389109675 27 60 12 399.805 - -1.2687853818468126E-06 -6.4616539890669800E-01 2.5515959980945999E-05 - -8.1379434569845700E-09 2.7083785173637931E-11 -1.7657780568571886E-12 - -1.4923807037190771E-13 4.6294307431546633E-15 1.1991673490193981E-16 - -4.6475524441121938E-18 -3.3996450389447711E-20 1.5783774939737859E-21 -1646-2142 16-JUL-19 60000.00 58680.25000000000 29.741000 0.627 -5.994 - 30758759705.348343 170.849389109675 27 60 12 399.805 - 7.1302303210167561E-07 -6.4323385587256643E-01 2.3012120641518147E-05 - -1.7964580317615194E-08 -1.2716514782721942E-10 -7.5135162455256370E-12 - 1.9814762175268489E-13 1.4502209295624714E-14 -1.9029866821569317E-16 - -1.1896939958536328E-17 6.3453412956578960E-20 3.4937858261643271E-21 -1646-2142 16-JUL-19 70000.00 58680.29166666660 29.741000 0.625 -6.094 - 30759374724.632160 170.849389109675 27 60 12 399.805 - 1.0190185092580351E-06 -6.4072091758068728E-01 1.8615862121575928E-05 - -3.1467829828280117E-08 -2.0349944549937411E-10 1.0202305811048701E-11 - 3.4866667634431997E-13 -1.7671636435021878E-14 -3.0053562568260006E-16 - 1.3098043269415902E-17 9.1137953371562547E-20 -3.4920337200230178E-21 -1646-2142 16-JUL-19 80000.00 58680.33333333330 29.741000 0.623 -6.125 - 30759989744.050339 170.849389109675 27 60 12 399.805 - -7.7023406835300540E-07 -6.3883574060311299E-01 1.2574678106379323E-05 - -4.0248451862470333E-08 9.4592028577544807E-11 1.4423443830213301E-11 - -2.5619567036202616E-13 -2.6797537878617611E-14 2.2787200844027115E-16 - 2.2057258051679547E-17 -7.1690444700979040E-20 -6.5984416076430977E-21 -1646-2142 16-JUL-19 90000.00 58680.37500000000 29.741000 0.622 -6.195 - 30760604763.558510 170.849389109675 27 60 12 399.805 - 2.3180999508473632E-07 -6.3774427630257668E-01 5.4462683704040882E-06 - -4.4534487313618588E-08 -5.4482032489873619E-11 8.0915685957177607E-12 - 8.5148086417058391E-14 -1.1114348483945886E-14 -7.8871607648778804E-17 - 6.6997947453988480E-18 2.5349901578397431E-20 -1.4453468910856642E-21 -1646-2142 16-JUL-19 100000.00 58680.41666666660 29.741000 0.622 -6.021 - 30761219783.103954 170.849389109675 27 60 12 399.805 - -2.5376035999390666E-07 -6.3755796009030141E-01 -2.3863144818959286E-06 - -4.5389596721810080E-08 -3.8389925873659049E-11 5.1029273193761658E-12 - 4.3984597030780682E-14 -1.0126331864067115E-14 -3.1248508576031709E-17 - 8.6662692813953835E-18 1.0410854457781258E-20 -2.6666259473221699E-21 -1646-2142 16-JUL-19 110000.00 58680.45833333330 29.741000 0.623 -6.129 - 30761834802.633518 170.849389109675 27 60 12 399.805 - 7.8154137841957017E-08 -6.3832586322925400E-01 -1.0380766455769799E-05 - -4.3999064833627125E-08 2.9070782336506895E-11 1.1420342680902956E-12 - -3.4054986408853028E-14 -2.1457628742134050E-15 2.0239916040045613E-17 - 1.9166882413382564E-18 -3.2686690912645793E-21 -6.1902955557892337E-22 -1646-2142 16-JUL-19 120000.00 58680.50000000000 29.741000 0.624 -6.013 - 30762449822.088608 170.849389109675 27 60 12 399.805 - 2.6047491756206465E-07 -6.4003242828161200E-01 -1.7958615505573728E-05 - -3.7587453827946840E-08 4.4216213724831579E-11 -9.9602846676872249E-12 - -1.3472952494364641E-15 1.6933942674134851E-14 -3.2401349643938679E-17 - -1.2552328744976534E-17 1.9464292813184483E-20 3.3997108708853653E-21 -1646-2142 16-JUL-19 130000.00 58680.54166666660 29.741000 0.627 -6.062 - 30763064841.413540 170.849389109675 27 60 12 399.805 - -2.1431161545057420E-07 -6.4259668551905358E-01 -2.4589387745277951E-05 - -3.3776810472224466E-08 -9.9663765747339841E-12 -5.8923504899482227E-13 - 9.6274624067119722E-14 3.4392270172043182E-15 -9.7861524653993575E-17 - -4.1729750824053053E-18 3.5389921953226819E-20 1.5539678792308668E-21 -1646-2142 16-JUL-19 140000.00 58680.58333333330 29.741000 0.630 -6.061 - 30763679860.563679 170.849389109675 27 60 12 399.805 - 1.2047510745266665E-07 -6.4588025844073993E-01 -2.9844395742360618E-05 - -2.5512129104921212E-08 -8.4689778136558570E-11 3.9568208967380159E-12 - 2.7654998852601330E-13 -6.4137431389582802E-15 -2.5344398384677209E-16 - 4.2821425517316729E-18 8.2249167535479223E-20 -1.0058991559733679E-21 -1646-2142 16-JUL-19 150000.00 58680.62500000000 29.741000 0.634 -6.056 - 30764294879.499847 170.849389109675 27 60 12 399.805 - -7.1695190848042323E-07 -6.4969434356281985E-01 -3.3415020054816211E-05 - -1.4028690407807038E-08 1.7816995660780527E-10 -1.8365163147539640E-12 - -2.5220229181841623E-13 4.3857266984654827E-15 2.0823421949581510E-16 - -3.9061751502828010E-18 -6.2477094963343525E-20 1.1790967088294452E-21 -1646-2142 16-JUL-19 160000.00 58680.66666666660 29.741000 0.638 -6.057 - 30764909898.195328 170.849389109675 27 60 12 399.805 - -9.1679951452761088E-07 -6.5381416866528175E-01 -3.4920732972053540E-05 - -3.1345671532416875E-09 -1.3018524884050091E-11 1.6017942405902710E-12 - 1.4196531223523833E-13 -3.0977867804354248E-15 -1.3835319121620314E-16 - 2.7899817872410146E-18 4.7685450852375216E-20 -9.5165585153929196E-22 -1646-2142 16-JUL-19 170000.00 58680.70833333330 29.741000 0.642 -6.032 - 30765524916.642040 170.849389109675 27 60 12 399.805 - -1.4787113040565236E-06 -6.5799403301006643E-01 -3.4393407248340235E-05 - 8.7724556677698359E-09 1.2243798947357433E-10 -1.2869600756153459E-12 - -1.7921040455565989E-13 4.9612298707182527E-15 1.6722163135548602E-16 - -5.7155741066117643E-18 -5.3702040261906497E-20 2.0865377452363192E-21 -1646-2142 16-JUL-19 180000.00 58680.75000000000 29.741000 0.646 -5.996 - 30766139934.842422 170.849389109675 27 60 12 399.805 - -1.4731669113580981E-06 -6.6198406562468726E-01 -3.1781075916570501E-05 - 2.0441074649857774E-08 2.3290981803247381E-11 -2.6413593383951564E-12 - 3.8019626152889040E-14 4.7212768956454758E-15 -2.7994829213005767E-17 - -3.7033201988560320E-18 7.2616147058982333E-21 1.0622893053712181E-21 -1646-2142 16-JUL-19 190000.00 58680.79166666660 29.741000 0.649 -6.027 - 30766754952.813606 170.849389109675 27 60 12 399.805 - 6.5321622508274224E-07 -6.6554722885126771E-01 -2.7306505216436970E-05 - 2.9704300281943558E-08 -3.5004843984045489E-11 -5.0426925917702495E-13 - 1.9104071032858260E-13 -3.4931083609746082E-16 -1.8646296074534826E-16 - 1.2139270009230189E-18 6.1264537900827979E-20 -5.7384996783271231E-22 -1646-2142 16-JUL-19 200000.00 58680.83333333330 29.741000 0.652 -6.085 - 30767369970.590569 170.849389109675 27 60 12 399.805 - -2.0421623957850299E-09 -6.6847620948001063E-01 -2.1284744692755857E-05 - 3.9366621775214088E-08 8.5265223429458146E-11 -7.8771059443487015E-12 - -1.1674825178136700E-13 1.3422910086993835E-14 1.1014375870212447E-16 - -1.0697333239377615E-17 -3.8443838967088908E-20 3.1714660352015652E-21 -1646-2142 16-JUL-19 210000.00 58680.87500000000 29.741000 0.654 -6.085 - 30767984988.215073 170.849389109675 27 60 12 399.805 - 1.0207051816240459E-06 -6.7060604380968358E-01 -1.4094978773637925E-05 - 4.2908898285634543E-08 9.5992688828422001E-11 -4.9258326476787876E-12 - -1.3572318907302694E-13 1.1839487928477319E-14 1.0157883130669997E-16 - -1.1303372071709341E-17 -2.8781806412350808E-20 3.7335324369433741E-21 -1646-2142 16-JUL-19 220000.00 58680.91666666660 29.741000 0.655 -6.007 - 30768600005.737148 170.849389109675 27 60 12 399.805 - 1.4019630041259279E-06 -6.7182774090495823E-01 -6.2360501618348241E-06 - 4.3595798868433080E-08 7.0454449248975137E-11 3.4316843220585464E-12 - -1.2653255276861120E-13 -5.8034399941902432E-15 1.0875524758854586E-16 - 4.3435599329146639E-18 -3.5238134010759204E-20 -1.2120634208077183E-21 -1646-2142 16-JUL-19 230000.00 58680.95833333330 29.741000 0.656 -6.040 - 30769215023.215984 170.849389109675 27 60 12 399.805 - -8.1677051358875104E-07 -6.7209412083670395E-01 1.7560965482431596E-06 - 4.3105209664863310E-08 1.7026210770236738E-11 4.0740575626097051E-12 - -7.3310417543462938E-14 -9.2390384606999324E-15 7.2994027258152459E-17 - 8.2178453732771450E-18 -2.4676210783225960E-20 -2.5614508354215359E-21 -1646-2142 17-JUL-19 0.00 58681.00000000000 29.741000 0.655 -6.013 - 30769830040.707218 170.849389109675 27 60 12 399.805 - -1.5092712970891931E-08 -6.7142323876627985E-01 9.3220482975622859E-06 - 4.1757190976232566E-08 -2.6588971124909429E-11 -4.5753957752068857E-12 - -4.9934920461462052E-14 4.2937214421701118E-15 8.1008330695105380E-17 - -1.6647672727848908E-18 -3.4556234904508793E-20 1.9270474523488295E-22 -1646-2142 17-JUL-19 10000.00 58681.04166666660 29.741000 0.653 -6.060 - 30770445058.263462 170.849389109675 27 60 12 399.805 - -1.1356303390129552E-06 -6.6989681214307595E-01 1.5941432518764502E-05 - 3.2552582172082968E-08 -1.4295903533049594E-10 2.1560075120381590E-12 - 2.7721771148208882E-13 -2.9356593070566292E-15 -2.7501981504646289E-16 - 1.8966151131337988E-18 9.3384829185689269E-20 -4.7786988887211809E-22 -1646-2142 17-JUL-19 20000.00 58681.08333333330 29.741000 0.651 -6.045 - 30771060075.935009 170.849389109675 27 60 12 399.805 - 1.2940663757180643E-06 -6.6765645376991367E-01 2.1152431667121763E-05 - 2.3906786417549138E-08 -9.2878667139130984E-11 1.4732150465670682E-13 - 1.2584354201279343E-13 1.4677483708886908E-15 -1.1480260815023348E-16 - -1.7455329446687277E-18 3.5497001601729317E-20 5.2522520882384874E-22 -1646-2142 17-JUL-19 30000.00 58681.12500000000 29.741000 0.649 -6.051 - 30771675093.757759 170.849389109675 27 60 12 399.805 - -9.9053088988346194E-07 -6.6489187478012790E-01 2.4616642389932379E-05 - 1.3004918697859133E-08 -7.3889038159512853E-11 2.1330929094158206E-12 - 9.0465234022172007E-14 -2.3550212639581387E-15 -1.0736264846013305E-16 - 1.4392788936169179E-18 4.1786581827054860E-20 -4.0102914515760781E-22 -1646-2142 17-JUL-19 40000.00 58681.16666666660 29.741000 0.646 -6.092 - 30772290111.755054 170.849389109675 27 60 12 399.805 - 9.0892705766412729E-07 -6.6182878275679768E-01 2.6093657228143540E-05 - 2.6135295208110966E-09 -2.2036567262996684E-11 -3.1765371124727722E-12 - -6.4446202833324900E-14 8.4178775889926837E-15 4.9897269770088399E-17 - -8.1328411170305371E-18 -1.2062246372087634E-20 2.6667452302623327E-21 diff --git a/tests/data/test.raw b/tests/data/test.raw deleted file mode 100644 index 03e28db..0000000 Binary files a/tests/data/test.raw and /dev/null differ diff --git a/tests/data/test_16bit.fil b/tests/data/test_16bit.fil deleted file mode 100644 index 6f71c05..0000000 Binary files a/tests/data/test_16bit.fil and /dev/null differ diff --git a/tests/data/test_1bit.fil b/tests/data/test_1bit.fil deleted file mode 100644 index f540c8e..0000000 Binary files a/tests/data/test_1bit.fil and /dev/null differ diff --git a/tests/data/test_2bit.fil b/tests/data/test_2bit.fil deleted file mode 100644 index 24151d2..0000000 Binary files a/tests/data/test_2bit.fil and /dev/null differ diff --git a/tests/data/test_32bit.fil b/tests/data/test_32bit.fil deleted file mode 100644 index b9246d4..0000000 Binary files a/tests/data/test_32bit.fil and /dev/null differ diff --git a/tests/data/test_4bit.fil b/tests/data/test_4bit.fil deleted file mode 100644 index b76332d..0000000 Binary files a/tests/data/test_4bit.fil and /dev/null differ diff --git a/tests/data/test_8bit.fil b/tests/data/test_8bit.fil deleted file mode 100644 index 81ef9a4..0000000 Binary files a/tests/data/test_8bit.fil and /dev/null differ diff --git a/tests/data/test_f32.tim b/tests/data/test_f32.tim deleted file mode 100644 index 8f3cba2..0000000 Binary files a/tests/data/test_f32.tim and /dev/null differ diff --git a/tests/data/test_i8.tim b/tests/data/test_i8.tim deleted file mode 100644 index a2a4dad..0000000 Binary files a/tests/data/test_i8.tim and /dev/null differ diff --git a/tests/data/test_radio.inf b/tests/data/test_radio.inf deleted file mode 100644 index fc2f57b..0000000 --- a/tests/data/test_radio.inf +++ /dev/null @@ -1,22 +0,0 @@ - Data file name without suffix = fake_presto_radio - Telescope used = Parkes - Instrument used = Multibeam - Object being observed = Pulsar - J2000 Right Ascension (hh:mm:ss.ssss) = 00:00:01.0000 - J2000 Declination (dd:mm:ss.ssss) = -00:00:01.0000 - Data observed by = Kenji Oba - Epoch of observation (MJD) = 59000.0 - Barycentered? (1 yes, 0 no) = 1 - Number of bins in the time series = 16 - Width of each time series bin (sec) = 6.4e-05 - Any breaks in the data? (1 yes, 0 no) = 0 - Type of observation (EM band) = Radio - Beam diameter (arcsec) = 981.0 - Dispersion measure (cm-3 pc) = 42.42 - Central freq of low channel (MHz) = 1182.1953125 - Total bandwidth (MHz) = 400.0 - Number of channels = 1024 - Channel bandwidth (MHz) = 0.390625 - Data analyzed by = Space Sheriff Gavan - Any additional notes: - Input filterbank samples have 2 bits. diff --git a/tests/data/test_radio_breaks.inf b/tests/data/test_radio_breaks.inf deleted file mode 100644 index 1c3c687..0000000 --- a/tests/data/test_radio_breaks.inf +++ /dev/null @@ -1,24 +0,0 @@ - Data file name without suffix = fake_presto_radio_breaks - Telescope used = Parkes - Instrument used = Multibeam - Object being observed = Pulsar - J2000 Right Ascension (hh:mm:ss.ssss) = 00:00:01.0000 - J2000 Declination (dd:mm:ss.ssss) = -00:00:01.0000 - Data observed by = Kenji Oba - Epoch of observation (MJD) = 59000.000000 - Barycentered? (1 yes, 0 no) = 1 - Number of bins in the time series = 16 - Width of each time series bin (sec) = 6.4e-05 - Any breaks in the data? (1 yes, 0 no) = 1 - On/Off bin pair # 1 = 0 , 14 - On/Off bin pair # 2 = 15 , 15 - Type of observation (EM band) = Radio - Beam diameter (arcsec) = 981 - Dispersion measure (cm-3 pc) = 42.42 - Central freq of low channel (MHz) = 1182.1953125 - Total bandwidth (MHz) = 400 - Number of channels = 1024 - Channel bandwidth (MHz) = 0.390625 - Data analyzed by = Space Sheriff Gavan - Any additional notes: - Input filterbank samples have 2 bits. diff --git a/tests/data/test_ui8.tim b/tests/data/test_ui8.tim deleted file mode 100644 index 56601ab..0000000 Binary files a/tests/data/test_ui8.tim and /dev/null differ diff --git a/tests/data/test_xray.inf b/tests/data/test_xray.inf deleted file mode 100644 index c2198df..0000000 --- a/tests/data/test_xray.inf +++ /dev/null @@ -1,19 +0,0 @@ - Data file name without suffix = fake_presto_xray - Telescope used = Chandra - Instrument used = HRC-S - Object being observed = Pulsar - J2000 Right Ascension (hh:mm:ss.ssss) = 00:00:01.0000 - J2000 Declination (dd:mm:ss.ssss) = -00:00:01.0000 - Data observed by = Kenji Oba - Epoch of observation (MJD) = 59000.000000 - Barycentered? (1=yes, 0=no) = 1 - Number of bins in the time series = 16 - Width of each time series bin (sec) = 6.4e-05 - Any breaks in the data? (1=yes, 0=no) = 0 - Type of observation (EM band) = X-ray - Field-of-view diameter (arcsec) = 3.000000 - Central energy (kev) = 1.000000 - Energy bandpass (kev) = 5.000000 - Data analyzed by = Space Sheriff Gavan - Any additional notes: - Full ms-resolution analysis diff --git a/tests/test_bpf.py b/tests/test_bpf.py deleted file mode 100644 index 5e2577a..0000000 --- a/tests/test_bpf.py +++ /dev/null @@ -1,57 +0,0 @@ -from ward import test -from ward import fixture -from pathlib import Path -from tempfile import NamedTemporaryFile -from priwo.presto.bpf import readbpf, writebpf - - -@fixture -def data(): - return Path(__file__).parent.joinpath("data") - - -def check(f): - meta, _ = readbpf(f) - assert meta == dict( - asini_by_c=None, - bary_epoch=None, - candidate="PSR_1646-2142", - data_avg=11042.9913182359, - data_std=246.537468521724, - dm=29.741, - eccentricity=None, - filename="J1646-2142_500_200_512_2.15jul2k19.raw0.fil", - noise_sigma=112.9, - nsamp=263664000.0, - p_bary=None, - p_orb=None, - p_topo=5.85347537656615, - p_topo_err=5.16e-08, - pd_bary=None, - pd_topo=-0.0, - pd_topo_err=1.48e-13, - pdd_bary=None, - pdd_topo=0.0, - pdd_topo_err=3.55e-16, - prof_avg=22747152667.3743, - prof_bins=128.0, - prof_std=353837.098935206, - red_chi_sqr=106.001, - telescope="GMRT", - topo_epoch=58679.6856653963, - t_peri=None, - tsamp=1.024e-05, - w=None, - ) - - -@test(f"{str(readbpf.__doc__).strip()}") -def _(f=data().joinpath("test.bestprof")): - check(f) - - -@test(f"{str(writebpf.__doc__).strip()}") -def _(f=data().joinpath("test.bestprof")): - with NamedTemporaryFile(suffix=".bestprof") as fp: - writebpf(*readbpf(f), fp.name) - check(fp.name) diff --git a/tests/test_dat.py b/tests/test_dat.py deleted file mode 100644 index f15d8c6..0000000 --- a/tests/test_dat.py +++ /dev/null @@ -1,54 +0,0 @@ -import numpy as np - -from ward import test -from ward import fixture -from pathlib import Path -from tempfile import NamedTemporaryFile -from priwo.presto.dat import readdat, writedat - - -@fixture -def data(): - return Path(__file__).parent.joinpath("data") - - -@fixture -def array(): - return np.asarray( - [ - 0.0, - 1.0, - 2.0, - 3.0, - 4.0, - 5.0, - 6.0, - 7.0, - 8.0, - 9.0, - 10.0, - 11.0, - 12.0, - 13.0, - 14.0, - 15.0, - ], - dtype=np.float32, - ) - - -def check(f): - _, data = readdat(f) - assert np.allclose(data, array()) - - -@test(f"{str(readdat.__doc__).strip()}") -def _(f=data().joinpath("test.dat")): - check(f) - - -@test(f"{str(writedat.__doc__).strip()}") -def _(f=data().joinpath("test.dat")): - with NamedTemporaryFile(suffix=".bestprof") as fp: - writedat(*readdat(f), fp.name) - check(fp.name) diff --git a/tests/test_fft.py b/tests/test_fft.py deleted file mode 100644 index 89f3a9d..0000000 --- a/tests/test_fft.py +++ /dev/null @@ -1,54 +0,0 @@ -import numpy as np - -from ward import test -from ward import fixture -from pathlib import Path -from tempfile import NamedTemporaryFile -from priwo.presto.fft import readfft, writefft - - -@fixture -def data(): - return Path(__file__).parent.joinpath("data") - - -@fixture -def array(): - return np.asarray( - [ - 120.0, - -8.0, - -7.9999995, - 40.218716, - -8.0, - 19.31371, - -8.0, - 11.972846, - -8.0, - 8.0, - -8.0, - 5.3454294, - -8.0, - 3.3137085, - -8.0, - 1.5912989, - ], - dtype=np.float32, - ) - - -def check(f): - _, data = readfft(f) - assert np.allclose(data, array()) - - -@test(f"{str(readfft.__doc__).strip()}") -def _(f=data().joinpath("test.fft")): - check(f) - - -@test(f"{str(writefft.__doc__).strip()}") -def _(f=data().joinpath("test.fft")): - with NamedTemporaryFile(suffix=".fft") as fp: - writefft(*readfft(f), fp.name) - check(fp.name) diff --git a/tests/test_fil.py b/tests/test_fil.py deleted file mode 100644 index a751e34..0000000 --- a/tests/test_fil.py +++ /dev/null @@ -1,50 +0,0 @@ -import numpy as np - -from ward import test -from ward import fixture -from pathlib import Path -from tempfile import NamedTemporaryFile -from priwo.sigproc.fil import readfil, writefil - - -@fixture -def data(): - return Path(__file__).parent.joinpath("data") - - -def check(f, array): - _, data = readfil(f) - assert np.allclose(data[64, 100:110], array) - - -for n, array in { - 1: np.asarray([1, 1, 0, 0, 0, 0, 0, 0, 1, 1], dtype=np.uint8), - 2: np.asarray([0, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=np.uint8), - 4: np.asarray([7, 4, 3, 11, 6, 10, 10, 9, 6, 7], dtype=np.uint8), - 8: np.asarray([121, 94, 94, 124, 151, 118, 132, 74, 112, 65], dtype=np.uint8), - 32: np.asarray( - [ - 1.166237, - -0.84468514, - 0.874816, - 1.4028563, - -0.98618776, - -0.80890864, - -1.6307002, - 1.1306021, - 0.50498164, - -1.6316832, - ], - dtype=np.float32, - ), -}.items(): - - @test(f"{str(readfil.__doc__).strip()} ({n}-bit data).") - def _(f=data().joinpath(f"test_{n}bit.fil"), array=array): - check(f, array) - - @test(f"{str(writefil.__doc__).strip()} ({n}-bit data).") - def _(f=data().joinpath(f"test_{n}bit.fil"), array=array): - with NamedTemporaryFile(suffix=".fil") as fp: - writefil(*readfil(f), fp.name) - check(fp.name, array) diff --git a/tests/test_hdr.py b/tests/test_hdr.py deleted file mode 100644 index 93ec73e..0000000 --- a/tests/test_hdr.py +++ /dev/null @@ -1,48 +0,0 @@ -from ward import test -from ward import fixture -from pathlib import Path -from tempfile import NamedTemporaryFile -from priwo.sigproc.hdr import readhdr, writehdr - - -@fixture -def data(): - return Path(__file__).parent.joinpath("data") - - -def check(f): - assert readhdr(f) == dict( - rawdatafile="./small.fil", - source_name="src1", - machine_id=0, - barycentric=0, - pulsarcentric=0, - telescope_id=6, - src_raj=122637.6361, - src_dej=135752.112, - az_start=-1.0, - za_start=-1.0, - data_type=0, - fch1=1465.0, - foff=-1.0, - nchans=336, - nbeams=1, - ibeam=0, - nbits=8, - tstart=58682.620316710374, - tsamp=0.00126646875, - nifs=1, - size=389, - ) - - -@test(f"{str(readhdr.__doc__).strip()}") -def _(f=data().joinpath("test.fil")): - check(f) - - -@test(f"{str(writehdr.__doc__).strip()}") -def _(f=data().joinpath("test.fil")): - with NamedTemporaryFile(suffix=".fil") as fp: - writehdr(readhdr(f), fp.name) - check(fp.name) diff --git a/tests/test_inf.py b/tests/test_inf.py deleted file mode 100644 index 815d7f1..0000000 --- a/tests/test_inf.py +++ /dev/null @@ -1,124 +0,0 @@ -from ward import test -from ward import fixture -from pathlib import Path -from tempfile import NamedTemporaryFile -from priwo.presto.inf import readinf, writeinf - - -@fixture -def data(): - return Path(__file__).parent.joinpath("data") - - -def check(f): - assert readinf(f) == dict( - analyst="Space Sheriff Gavan", - barycentered=True, - beamdiam=981.0, - breaks=False, - bw=400.0, - cfreq=1182.1953125, - chanwidth=0.390625, - dec="-00:00:01.0000", - dm=42.42, - emband="Radio", - filename="fake_presto_radio", - instrument="Multibeam", - mjd=59000.0, - nchannels=1024, - notes="Input filterbank samples have 2 bits.", - nsamples=16, - object="Pulsar", - observer="Kenji Oba", - onoffs=[], - ra="00:00:01.0000", - samptime=6.4e-05, - telescope="Parkes", - ) - - -def check_breaks(f): - assert readinf(f) == dict( - analyst="Space Sheriff Gavan", - barycentered=True, - beamdiam=981.0, - breaks=True, - bw=400.0, - cfreq=1182.1953125, - chanwidth=0.390625, - dec="-00:00:01.0000", - dm=42.42, - emband="Radio", - filename="fake_presto_radio_breaks", - instrument="Multibeam", - mjd=59000.0, - nchannels=1024, - notes="Input filterbank samples have 2 bits.", - nsamples=16, - object="Pulsar", - observer="Kenji Oba", - onoffs=[[0, 14], [15, 15]], - ra="00:00:01.0000", - samptime=6.4e-05, - telescope="Parkes", - ) - - -def check_xray(f): - assert readinf(f) == dict( - analyst="Space Sheriff Gavan", - barycentered=True, - bpE=5.0, - breaks=False, - cE=1.0, - dec="-00:00:01.0000", - emband="X-ray", - filename="fake_presto_xray", - fov=3.0, - instrument="HRC-S", - mjd=59000.0, - notes="Full ms-resolution analysis", - nsamples=16, - object="Pulsar", - observer="Kenji Oba", - onoffs=[], - ra="00:00:01.0000", - samptime=6.4e-05, - telescope="Chandra", - ) - - -@test(f"{str(readinf.__doc__).strip()}, for radio data.") -def _(f=data().joinpath("test_radio.inf")): - check(f) - - -@test(f"{str(readinf.__doc__).strip()}, for radio data, but with breaks.") -def _(f=data().joinpath("test_radio_breaks.inf")): - check_breaks(f) - - -@test(f"{str(readinf.__doc__).strip()} for xray data.") -def _(f=data().joinpath("test_xray.inf")): - check_xray(f) - - -@test(f"{str(writeinf.__doc__).strip()}, for radio data.") -def _(f=data().joinpath("test_radio.inf")): - with NamedTemporaryFile(suffix=".inf") as fp: - writeinf(readinf(f), fp.name) - check(fp.name) - - -@test(f"{str(writeinf.__doc__).strip()}, for radio data, but with breaks.") -def _(f=data().joinpath("test_radio_breaks.inf")): - with NamedTemporaryFile(suffix=".inf") as fp: - writeinf(readinf(f), fp.name) - check_breaks(fp.name) - - -@test(f"{str(writeinf.__doc__).strip()}, for xray data.") -def _(f=data().joinpath("test_xray.inf")): - with NamedTemporaryFile(suffix=".inf") as fp: - writeinf(readinf(f), fp.name) - check_xray(fp.name) diff --git a/tests/test_pfd.py b/tests/test_pfd.py deleted file mode 100644 index 52602ef..0000000 --- a/tests/test_pfd.py +++ /dev/null @@ -1,92 +0,0 @@ -import numpy as np - -from ward import test -from ward import fixture -from pathlib import Path -from tempfile import NamedTemporaryFile -from priwo.presto.pfd import readpfd, writepfd - - -@fixture -def data(): - return Path(__file__).parent.joinpath("data") - - -def check(f): - meta, _ = readpfd(f) - for _ in [ - "dms", - "pdots", - "stats", - "periods", - ]: - meta.pop(_) - - assert meta == dict( - ndms=257, - nperiods=201, - npdots=201, - nsub=128, - npart=60, - nbin=128, - nchan=512, - pstep=1, - pdstep=2, - dmstep=1, - ndmfact=1, - npfact=1, - filename="J1646-2142_500_200_512_2.15jul2k19.raw0.fil", - candname="PSR_1646-2142", - telescope="GMRT", - pgdev="J1646-2142_500_200_512_2.15jul2k19.raw0.new_freq_PSR_1646-2142.pfd.ps/CPS", - ra="16:46:18.1250", - dec="-21:42:08.9063", - dt=1.024e-05, - T0=0.0, - Tn=1.0, - tepoch=58679.6856653963, - bepoch=0.0, - vavg=0.0, - f0=300.0, - df=0.390625, - bestdm=29.741, - topo=dict( - power=0.0, - p=0.005853475376566149, - pd=-0.0, - pdd=0.0, - ), - bary=dict( - power=0.0, - p=0.0, - pd=0.0, - pdd=0.0, - ), - fold=dict( - power=6.9527192090977e-310, - p=170.8386788476822, - pd=0.0, - pdd=0.0, - ), - orb=dict( - p=0.0, - e=0.0, - x=0.0, - w=0.0, - t=0.0, - pd=0.0, - wd=0.0, - ), - ) - - -@test(f"{str(readpfd.__doc__).strip()}") -def _(f=data().joinpath("test.pfd")): - check(f) - - -@test(f"{str(writepfd.__doc__).strip()}") -def _(f=data().joinpath("test.pfd")): - with NamedTemporaryFile(suffix=".pfd") as fp: - writepfd(*readpfd(f), fp.name) - check(fp.name) diff --git a/tests/test_polycos.py b/tests/test_polycos.py deleted file mode 100644 index d33cb89..0000000 --- a/tests/test_polycos.py +++ /dev/null @@ -1,41 +0,0 @@ -from ward import test -from ward import fixture -from pathlib import Path -from tempfile import NamedTemporaryFile -from priwo.timing.polycos import readpolycos, writepolycos - - -@fixture -def data(): - return Path(__file__).parent.joinpath("data") - - -def check(f): - assert readpolycos(f)[0]["meta"] == dict( - name="1646-2142", - date="14-JUL-19", - utc="50000.00", - tmid=58678.2083333333, - dm=29.741, - doppler=0.605, - rms=-6.064, - refphz=30728623742.9709, - reff0=170.849389109675, - obsno=27, - span=60.0, - ncoeff=12, - freq=399.805, - binphz=105.0, - ) - - -@test(f"{str(readpolycos.__doc__).strip()}") -def _(f=data().joinpath("test.polycos")): - check(f) - - -@test(f"{str(writepolycos.__doc__).strip()}") -def _(f=data().joinpath("test.polycos")): - with NamedTemporaryFile(suffix=".polycos") as fp: - writepolycos(readpolycos(f), fp.name) - check(fp.name) diff --git a/tests/test_tim.py b/tests/test_tim.py deleted file mode 100644 index 6de137d..0000000 --- a/tests/test_tim.py +++ /dev/null @@ -1,59 +0,0 @@ -import numpy as np - -from ward import test -from ward import fixture -from pathlib import Path -from tempfile import NamedTemporaryFile -from priwo.sigproc.tim import readtim, writetim - - -@fixture -def data(): - return Path(__file__).parent.joinpath("data") - - -@fixture -def array(): - return np.asarray( - [ - 0.0, - 1.0, - 2.0, - 3.0, - 4.0, - 5.0, - 6.0, - 7.0, - 8.0, - 9.0, - 10.0, - 11.0, - 12.0, - 13.0, - 14.0, - 15.0, - ], - dtype=np.float32, - ) - - -def check(f): - _, data = readtim(f) - assert np.allclose(data, array()) - - -for label, fname in [ - ("Signed 8-bit integer data", "test_ui8.tim"), - ("Unsigned 8-bit integer data", "test_i8.tim"), - ("32-bit floating point data", "test_f32.tim"), -]: - - @test(f"{str(readtim.__doc__).strip()} ({label}).") - def _(f=data().joinpath(f"{fname}")): - check(f) - - @test(f"{str(writetim.__doc__).strip()} ({label}).") - def _(f=data().joinpath(f"{fname}")): - with NamedTemporaryFile(suffix=".tim") as fp: - writetim(*readtim(f), fp.name) - check(fp.name)