diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..ccc3ef3 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,161 @@ +# Contributing to SymDALI + +This document is intended for collaborators who want to extend or maintain the codebase. It assumes familiarity with Wolfram Language and a basic understanding of GW parameter estimation. + +## Repository Layout + +``` +SymDALI/ +├── Kernel/ — Wolfram Language source packages +├── Rules/ — Pre-computed symbolic and numerical derivative rules +├── LibraryResources/ — Platform-specific compiled C libraries +├── Data/ASD/ — Noise curve data files +├── Documentation/ — Reference pages (Mathematica notebook format) +└── Demo/ — Usage examples +``` + +## Core Architecture + +Understanding the data flow is essential before contributing. + +### Gradient computation pipeline + +The central computation is the DALI tensor network sum over frequency bins: + +$$\mathcal{F}^{(i,j)}_{\mu_1 \cdots \mu_i \nu_1 \cdots \nu_j} = 4 \Delta f \sum_k \frac{\text{Re}\left[\partial_{\mu_1 \cdots \mu_i} \tilde{s}^*(f_k)\, \partial_{\nu_1 \cdots \nu_j} \tilde{s}(f_k)\right]}{S_n(f_k)}$$ + +where $\tilde{s}(f) = F_+ h_+(f) + F_\times h_\times(f)$ is the detector signal. The computation is factored: + +1. **Derivative rules loading** (`DerivativeTools.wl`): Pre-computed rules for amplitude $\mathcal{A}$, phase $\Phi$, and pattern functions $F_{+,\times}$ are loaded from disk as compiled C `LibraryFunction` objects. + +2. **Gradient computation** (`DALICoefficients.wl`): `iGenGrads` iteratively applies `TakeGrad` to build gradients up to the requested order, exploiting Mathematica's `SymmetrizedIndependentComponents` to work with only the linearly independent tensor components. + +3. **Tensor contraction** (`DALICoefficients.wl`): `GenDaliTerm` and `GenDaliList` perform the frequency-domain integral via the compiled `CiGenDaliTerm` function. + +4. **Post-processing** (`DALIPolynomial.wl`): `ProcessDALITensors` and `TaylorForm` convert the raw output into forms convenient for sampling. + +### Derivative rules format + +The `Rules/` directory stores pre-computed derivatives in two forms: + +- **`SymRules`**: Associations mapping `$D[{n,...}, f] -> number` for derivatives that evaluate to constants (e.g., phase derivatives that do not depend on frequency). Stored as `.wdx` or `.mx` files. +- **`NRules`**: Associations mapping `$D[{n,...}, f][args_] -> CompiledFunction[...]` for derivatives that are functions of the arguments. The `CompiledFunction` wraps a `LibraryFunction` loaded from a `.dylib` file. + +The helper `$D[{n1,n2,...}, f][x,y,...]` denotes the mixed partial derivative $\partial^{n_1+n_2+\cdots} f / \partial x_1^{n_1} \partial x_2^{n_2} \cdots$ evaluated at $(x,y,\ldots)$. + +## How to Add a New Detector + +### 1. Add detector tensor and vertex + +In `Kernel/Detectors.wl`, add entries to the `DetectorTensor` and `DetectorVertex` associations. Use `ArmDirection` and `FromGeocentricCoordinates` to compute the tensors from geographic coordinates: + +```mathematica +Module[ + {lat = ... Degree, lon = ... Degree, n1, n2}, + {n1, n2} = ArmDirection[{lat, lon}, {{omega1, psi1}, {omega2, psi2}}]; + DetectorTensor["MyDet"] = (n1⊗n1 - n2⊗n2)/2; +] + +DetectorVertex["MyDet"] = FromGeocentricCoordinates[{h, lat, lon}]; +``` + +The arm direction parameters `{omega_i, psi_i}` follow the WGS-84 convention from `LALDetectors.h`: `omega` is the altitude and `psi` is the azimuth of the arm. + +### 2. Add a noise curve + +Place a two-column text file (frequency [Hz], ASD [1/√Hz]) in `Data/ASD/`, then add an interpolation in `Kernel/Detectors.wl`: + +```mathematica +data = Import[FileNameJoin[{ASDdir, "my_detector_asd.txt"}], "Data"] // N; +ASD["MyDet"] = Interpolation[data]; +``` + +No changes to any other file are needed. + +## How to Add a New Waveform Model + +Adding a new waveform is the most involved type of contribution. The following steps are required. + +### 1. Implement the waveform function + +The waveform must be expressed as a product of an amplitude factor and a phase factor, separated so that the amplitude can be computed independently. Study `ihphcIMRPhenomD` in `Kernel/DALICoefficients.wl` for the expected structure. + +The function should take all physical parameters as separate scalar arguments (not a list or association), followed by a reference frequency and the evaluation frequency as the last two arguments. + +### 2. Generate derivative rules + +The most expensive step. Use `DerivativeRules` from `Kernel/DerivativeTools.wl` to compute all symbolic derivatives and compile them to C: + +```mathematica +(* Example: generate derivatives of the amplitude function *) +ampRules = DerivativeRules[ + {myAmp, {param1, param2, ...}, 3}, (* name, vars, max order *) + $Block[ + {{internalHead1, internalHead2}}, + internalHead1[x_] = ...; (* internal function definitions *) + myAmpExpression[param1, param2, ...] + ], + "KeepDefs" -> {localVar -> expression, ...} +] +``` + +For waveforms with many parameters, this can take hours. The demo notebooks in `Demo/` and `Documentation/English/Tutorials/` contain worked examples for IMRPhenomD. + +### 3. Save the derivative rules + +Use `SaveSymRules` and `SaveNRules` (or the serialization utilities in `DerivativeTools.wl`) to persist the rules to `Rules/YourModel/` and to `LibraryResources//DerivativeRules/YourModel/`. + +### 4. Register the model in `DALICoefficients.wl` + +Add entries in the `Module` block that loads derivative rules at package initialization, and implement the corresponding `ihphcYourModel` function. Register the model name in `headhphc`, `iihphcvars`, `isAligned`, and `iconvert` lookup tables. Add validation in `RetrieveFiducial` and `iDALITensors`. + +### 5. Add a public waveform evaluation function in `Utils.wl` + +Following the pattern of `hphcIMRPhenomD`, implement a user-facing function `hphcYourModel` with input validation and `Options` for any deviation parameters. + +## How to Generate New Derivative Rules for an Existing Model + +If you want to extend an existing model to a higher DALI order or to additional parameters, you need to generate the additional derivative rules and add them to the existing rule sets. + +The notebook `Demo/DerivativeRules_test.nb` and `Demo/manual_C_functions.nb` show examples of how this was done for IMRPhenomD. The key function is: + +```mathematica +DerivativeRules[{name, vars, n}, expr, opts] +``` + +Important options: +- `"IncludeZeroDerivative" -> True` — include the zeroth derivative (function value itself) +- `"KeepDefs" -> {var -> expr, ...}` — local variable substitutions to apply after differentiation +- `"Parallel" -> {k}` — use `k` parallel kernels + +## Platform Support + +All compiled C libraries are currently provided only for `MacOSX-ARM64`. The `DerivativeRulesLoad` function in `DerivativeTools.wl` uses `SystemInformation["Kernel", "SystemID"]` to locate the appropriate library directory. To add support for a new platform: + +1. Recompile the `.dylib` files on the target platform using the C source in `LibraryResources/MacOSX-ARM64/Polynomial_Functions/DALI.c` and by regenerating the waveform derivative rules. +2. Place the resulting libraries under `LibraryResources//`. + +The `SystemID` string for common platforms is: +- `"MacOSX-ARM64"` — macOS Apple Silicon +- `"MacOSX-x86-64"` — macOS Intel +- `"Linux-x86-64"` — Linux 64-bit + +## Code Conventions + +- **Wolfram Language style**: Follow the existing patterns. Internal (non-exported) functions are prefixed with a lowercase `i` (e.g., `iGenGrads`, `iDALITensors`). +- **Error handling**: Use `Throw[$Failed, failTag[functionName]]` for internal errors and `Throw["message"]` for user-facing errors. The top-level `DALITensors` wraps everything in `Catch`. +- **Protection**: All public symbols are `Protect`ed after definition. Use `Unprotect` before modifying, then `Protect` again. +- **Package contexts**: Each submodule lives in its own context under `FelipeBarbosa`SymDALI``. Keep package dependencies explicit via the `BeginPackage` second argument. +- **No comments on the obvious**: Comment only non-obvious derivations, physical assumptions, or known bugs/workarounds. + +## Testing + +There is no automated test suite. Validation is done by cross-checking DALI tensor components against independent numerical differentiation and against results from [gwfast](https://arxiv.org/abs/2203.02670) for the Fisher matrix limit. Before submitting changes: + +1. Run the `Demo/UserGuide.wl` examples end to end. +2. For Fisher matrix changes, verify that `DALITensors[..., 1][[1,1]]` agrees with an independent Fisher matrix calculation (e.g., numerical finite differences of the log-likelihood). +3. For waveform changes, compare `hphcYourModel` output against a reference implementation over a broad parameter range. + +## Opening Issues and Pull Requests + +This repository is currently in a pre-publication state. Please open an issue to discuss any proposed changes before starting significant development work, to avoid duplicating ongoing efforts. diff --git a/Demo/Overplot.png b/Demo/Overplot.png new file mode 100644 index 0000000..5fa6afe Binary files /dev/null and b/Demo/Overplot.png differ diff --git a/README.md b/README.md index ad3de4e..4a3029f 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,226 @@ -**SymDALI** is a Mathematica implementation of [DALI](https://arxiv.org/abs/1401.6892) with symbolic differentiation. It is actively being developed for applications in [gravitational waves](https://arxiv.org/abs/2203.02670). +# SymDALI -An example application is show in this [Figure](https://github.com/Felipe-4/SymDALI/tree/main/Demo/Overplot.pdf). It is the exact likelihood (black), Fisher approximation (blue), DALI order 2 (red) and DALI order 3 (green). The event has SNR of 109. The assumed network was 2 ET L-shaped and 1 CE detectors. The waveform is ```IMRPhenomD``` and it includes the extra parameters from the [```TIGER```](https://arxiv.org/pdf/1311.0420) framework, from inspiral to ringdown. This examples uses the parameter \delta_1. +**SymDALI** is a Wolfram Language (Mathematica) paclet that implements the [DALI](https://arxiv.org/abs/1401.6892) algorithm using symbolic differentiation. It is developed for gravitational-wave (GW) parameter estimation and is the companion code for a paper to be submitted to Physical Review D. -The DALI order 2 result takes around ~ 15 minutes (using [```emcee```](https://arxiv.org/pdf/1202.3665)) to run, while the exact likelihood takes around 13 hours (using [```nessai```](https://arxiv.org/pdf/2102.11056)). A publication with more examples is comming... . +## Overview -At the moment this repository supports the calculation of DALI tensors to order 3 for IMRPhenomD. The file ```UserGuide.wl``` in [Demo](https://github.com/Felipe-4/SymDALI/tree/main/Demo) has some examples on how to use part of the functionality of the code. Right now It is fully functional only on MacOs Apple Silicon M1. I will also be implementing Linux and Windows as soon as possible. +The Fisher information matrix gives the Gaussian (quadratic) approximation to the posterior of a GW signal, but it often fails for events with moderate signal-to-noise ratio (SNR) or when posterior distributions are non-Gaussian. DALI (Derivative Approximation for LIkelihoods) systematically extends the Fisher approximation to higher orders in the Taylor expansion of the log-likelihood: + +$$\ln \mathcal{L}(\theta) \approx -\frac{1}{2} \Gamma_{ij} \Delta\theta^i \Delta\theta^j - \frac{1}{6} \mathcal{F}_{ijk} \Delta\theta^i \Delta\theta^j \Delta\theta^k - \ldots$$ + +where $\Gamma_{ij}$ is the Fisher matrix, $\mathcal{F}_{ijk}$ and higher tensors are the DALI coefficients, and $\Delta\theta^i = \theta^i - \theta^i_0$ is the deviation from the fiducial point. The result is a polynomial approximation to the likelihood that can be sampled orders of magnitude faster than the exact likelihood. + +The figure below shows an example: exact likelihood (black), Fisher approximation (blue), DALI order 2 (red), and DALI order 3 (green) for a GW event with SNR 109, observed by a network of 2 Einstein Telescope L-shaped detectors and 1 Cosmic Explorer. The DALI order 2 result takes ~15 minutes with [emcee](https://arxiv.org/pdf/1202.3665), while the exact likelihood takes ~13 hours with [nessai](https://arxiv.org/pdf/2102.11056). + +![Overplot](Demo/Overplot.png) + +## Features + +- DALI tensor computation up to order 3 for GW signals +- Waveform models: + - **IMRPhenomD** (aligned-spin, fully supported) + - **IMRPhenomHM** (aligned-spin with higher modes, fully supported) + - **IMRPhenomPv2** (precessing, partially supported) +- GR deviation parameters from the [TIGER](https://arxiv.org/pdf/1311.0420) framework (inspiral δφ_{-2},...,δφ_7; intermediate δβ_2, δβ_3; ringdown δα_2, δα_3, δα_4) +- Detector network support: + - Current detectors: LIGO H1, L1; Virgo V1; KAGRA K; LIGO India I1 + - Next-generation: Einstein Telescope (triangle and L-shaped at Sardinia and Meuse-Rhine, 10/15/20 km arm lengths), Cosmic Explorer (Idaho and New Mexico) + - Noise curves: ET-D, CE-20, CE-40, CE-40-lf, CE-20-pm, A+/O5, KAGRA 80 Mpc, Virgo O5 +- Signal-to-noise ratio (SNR) computation +- Waveform evaluation (h+ and hx) for IMRPhenomD and IMRPhenomHM +- Pattern function computation including time delay from detector to Earth center +- Population model utilities: Madau-Dickinson star formation, Power-Law+Peak mass distribution (GWTC-3), spin distribution +- Cosmological utilities: comoving distance, differential comoving volume (Planck 2018) + +## Requirements + +- **Wolfram Mathematica 14.0+** +- **Platform**: macOS ARM64 (Apple Silicon, M1 or later) — Linux and Windows support is planned + +## Installation + +1. Clone or download this repository. +2. In a Mathematica notebook, load the paclet by pointing to the directory containing `PacletInfo.wl`: + +```mathematica +PacletDirectoryLoad["/path/to/SymDALI"]; +<< FelipeBarbosa`SymDALI` +``` + +Alternatively, install it as a paclet: + +```mathematica +PacletInstall["/path/to/SymDALI"]; +<< FelipeBarbosa`SymDALI` +``` + +## Quick Start + +### DALI Tensors + +The main function is `DALITensors`. It takes an approximant name, a fiducial point association, a list of detector associations, and the expansion order. + +```mathematica +(* Define detectors *) +det1 = <| + "Position" -> DetectorVertex["ET-S"], + "DetectorTensor" -> DetectorTensor["ET-S-L"], + "ASD" -> ASD["ET-10-lfhf"] +|>; + +det2 = <| + "Position" -> DetectorVertex["CE-I"], + "DetectorTensor" -> DetectorTensor["CE-I"], + "ASD" -> ASD["CE-40"] +|>; + +(* Define the fiducial point — all keys are strings *) +fp = <| + "\[ScriptCapitalM]c" -> 30.0, (* chirp mass [solar masses] *) + "\[Delta]" -> 0.0, (* mass asymmetry (m1-m2)/(m1+m2), range [0,1) *) + "\[Chi]s" -> 0.0, (* symmetric spin *) + "\[Chi]a" -> 0.0, (* antisymmetric spin *) + "\[Iota]" -> 0.4, (* inclination [radians] *) + "\[Theta]" -> 1.2, (* sky polar angle [radians] *) + "\[Phi]" -> 0.8, (* sky azimuthal angle [radians] *) + "\[Psi]" -> 0.3, (* polarization angle [radians] *) + "1/dL" -> 0.2, (* inverse luminosity distance [1/Gpc] *) + "tc" -> 0.0, (* coalescence time [s] *) + "\[Phi]ref" -> 0.0 (* reference phase [radians] *) +|>; + +(* Compute DALI tensors up to order 2 *) +dali = DALITensors["IMRPhenomD", fp, {det1, det2}, 2]; +``` + +The result is a nested list `dali[[i, j]]` where `i` runs from 1 to the expansion order and `j` from 1 to `i`. The diagonal entries `dali[[i, i]]` contain the purely symmetric DALI tensors of order `2i`. The full tensor at order `i+j` mixes gradients of order `i` and `j`. + +**Options:** + +| Option | Default | Description | +|--------|---------|-------------| +| `"fmin"` | `10` | Minimum frequency [Hz] | +| `"fmax"` | `1024` | Maximum frequency [Hz] (capped at the ISCO) | +| `"res"` | `1000` | Number of frequency points | +| `"AllFisherMatrices"` | `False` | If `True`, return per-detector contributions instead of the network sum | + +### Post-processing + +```mathematica +(* ProcessDALITensors: LI components with multiplicities, ready for contraction *) +processed = ProcessDALITensors[dali]; + +(* TaylorForm: totally symmetric tensors organized by rank *) +taylor = TaylorForm[dali]; +``` + +### SNR + +```mathematica +snr = SNR["IMRPhenomD", fp, {det1, det2}]; +``` + +### Waveform Evaluation + +```mathematica +fvec = Range[10, 512, 0.1]; (* frequency vector [Hz] *) + +(* hphcIMRPhenomD returns {{h+}, {hx}} as a matrix with dimensions {2, Length[fvec]} *) +hphc = hphcIMRPhenomD[fvec, Mc, delta, chis, chia, iota, tc, phiref, invdL]; + +(* With TIGER parameters *) +hphcTIGER = hphcIMRPhenomD[fvec, Mc, delta, chis, chia, iota, tc, phiref, invdL, + "\[Delta]\[CurlyPhi]-2" -> 0.5, "\[Delta]\[Beta]3" -> 0.1]; +``` + +### Detector Utilities + +```mathematica +(* Detector tensor and vertex position *) +Dij = DetectorTensor["ET-S-L"]; +pi = DetectorVertex["ET-S"]; + +(* Noise curve as InterpolatingFunction *) +Sn = ASD["ET-10-lfhf"]; +Sn[100] (* ASD at 100 Hz *) + +(* Pattern functions *) +{Fp, Fx} = PatternFunctions[fvec, delta, alpha, psi, GMST, pi, Dij]; +``` + +### Population Utilities + +```mathematica +(* Comoving distance (InterpolatingFunction) *) +Dc = ComovingDistance[]; + +(* Differential comoving volume (InterpolatingFunction) *) +dVcdz = DifferentialComovingVolume[]; + +(* Power-Law+Peak mass distribution (GWTC-3 hyperparameters) *) +p = PowerLawPlusPeak[m1, q]; + +(* Madau-Dickinson star formation rate *) +psi = MadauDickinsonProfile[z]; + +(* Spin distribution *) +p = SpinDistribution[chi1, chi2, costheta1, costheta2]; +``` + +## Package Structure + +``` +SymDALI/ +├── PacletInfo.wl # Paclet metadata +├── Kernel/ +│ ├── SymDALI.wl # Entry point, loads subpackages +│ ├── DALICoefficients.wl # Core DALI tensor computation (DALITensors) +│ ├── DerivativeTools.wl # Symbolic differentiation engine +│ ├── Detectors.wl # Detector tensors, vertices, and ASDs +│ ├── DALIPolynomial.wl # Post-processing (TaylorForm, ProcessDALITensors) +│ ├── Population.wl # Population distributions and cosmology +│ └── Utils.wl # SNR, waveform evaluation, pattern functions +├── Rules/ +│ ├── IMRPhenomD/ # Pre-computed derivative rules for IMRPhenomD +│ ├── IMRPhenomHM/ # Pre-computed derivative rules for IMRPhenomHM +│ └── IMRPhenomPV2/ # Pre-computed derivative rules for IMRPhenomPv2 +├── LibraryResources/ +│ └── MacOSX-ARM64/ +│ ├── DerivativeRules/ # Compiled C libraries (.dylib) for fast evaluation +│ └── Polynomial_Functions/ # Compiled DALI polynomial evaluator +├── Data/ +│ └── ASD/ # Noise curves for current and next-gen detectors +├── Documentation/ +│ └── English/ # Reference pages for public symbols +└── Demo/ + ├── UserGuide.wl # Usage examples + └── Overplot.pdf # Example figure +``` + +## Key Design Decisions + +**Symbolic differentiation with compiled evaluation.** The derivatives of the waveform amplitude and phase with respect to all parameters are computed once symbolically (using `DerivativeRules`) and saved to disk. At runtime, these are loaded as compiled C library functions (`LibraryLink`), allowing fast numerical evaluation without repeating the expensive symbolic computation. + +**Factored gradient structure.** The signal in detector `d` is `s_d(f) = F_+^d h_+(f) + F_x^d h_x(f)`. The gradient of `s_d` is computed by separately differentiating the pattern functions `(F_+, F_x)` with respect to sky and orientation angles, and the polarizations `(h_+, h_x)` with respect to intrinsic parameters. These are then combined via the chain rule. + +**SymmetrizedArray tensors.** All DALI tensors are inherently totally symmetric. Mathematica's `SymmetrizedArray` objects are used to store and manipulate them efficiently, keeping only the linearly independent (LI) components. + +## Citation + +If you use SymDALI in your research, please cite the companion paper (reference to be added upon publication) and the original DALI paper: + +- Sellentin, Heavens & Jasche (2014), [arXiv:1401.6892](https://arxiv.org/abs/1401.6892) + +## License + +MIT License. See [LICENSE](LICENSE) for details. + + +## Acknowledgements + +- [gwfast](https://arxiv.org/abs/2203.02670) — inspiration and cross-validation +- [emcee](https://arxiv.org/pdf/1202.3665) — MCMC sampling used in examples +- [nessai](https://arxiv.org/pdf/2102.11056) — nested sampling used in examples +- Detector sensitivity curves from the [Einstein Telescope](https://www.et-gw.eu) and [Cosmic Explorer](https://cosmicexplorer.org) projects + +This work is supported by FAPES-Brazil \ No newline at end of file