quadratic_cv_bias.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from gwbench import snr\n",
    "from gwbench import waveform as wfc\n",
    "from gwbench import injections\n",
    "from gwbench import network\n",
    "from gwbench import basic_relations as br\n",
    "from gwbench import detector_response_derivatives as drds\n",
    "\n",
    "from pycbc import conversions as conv\n",
    "from  pycbc.filter.matchedfilter import optimized_match\n",
    "\n",
    "from astropy.cosmology import Planck18, z_at_value\n",
    "\n",
    "import gwbench_network_funcs as gwnet\n",
    "\n",
    "from pycbc.types import FrequencySeries\n",
    "from pycbc.filter import match\n",
    "from pycbc.filter import matchedfilter\n",
    "\n",
    "from pycbc import types, fft, waveform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_file = 'smooth_powerlaw_pop.npz'\n",
    "with np.load(input_file, mmap_mode='r') as binaries:\n",
    "    Mcs = np.array(binaries['Mcs'])\n",
    "    etas = np.array(binaries['etas'])\n",
    "    chi1z = np.array(binaries['chi1z'])\n",
    "    chi2z = np.array(binaries['chi2z'])\n",
    "    DLs = np.array(binaries['DLs'])\n",
    "    iotas = np.array(binaries['iotas'])\n",
    "    ras = np.array(binaries['ras'])\n",
    "    decs = np.array(binaries['decs'])\n",
    "    psis = np.array(binaries['psis'])\n",
    "\n",
    "m1s = conv.mass1_from_mchirp_eta(Mcs, etas)\n",
    "mtots = conv.mtotal_from_mchirp_eta(Mcs, etas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate waveforms and derivatives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-11-20 12:49:21,134 - Network - INFO : Polarizations calculated.\n",
      "2023-11-20 12:49:21,135 - Network - INFO : Calculate numeric derivatives of polarizations.\n",
      "2023-11-20 12:49:21,143 - Network - INFO : Polarizations calculated.\n",
      "2023-11-20 12:49:22,413 - Network - INFO : Numeric derivatives of polarizations calculated.\n",
      "2023-11-20 12:49:22,488 - Network - INFO : PSDs loaded.\n",
      "2023-11-20 12:49:22,491 - Network - INFO : Antenna patterns and LPFs loaded.\n",
      "2023-11-20 12:49:22,520 - Network - INFO : Detector responses calculated.\n",
      "2023-11-20 12:49:22,521 - Network - INFO : Calculate numeric derivatives of detector responses.\n",
      "2023-11-20 12:49:22,522 - Network - INFO :    CE-40_C\n",
      "2023-11-20 12:49:23,416 - Network - INFO :    CE-20_S\n",
      "2023-11-20 12:49:24,331 - Network - INFO :    ET_ET1\n",
      "2023-11-20 12:49:25,326 - Network - INFO :    ET_ET2\n",
      "2023-11-20 12:49:26,179 - Network - INFO :    ET_ET3\n",
      "2023-11-20 12:49:27,391 - Network - INFO : Numeric derivatives of detector responses calculated.\n",
      "2023-11-20 12:49:27,394 - Network - INFO : SNRs calculated.\n",
      "2023-11-20 12:49:27,395 - Network - INFO : Calculate errors (Fisher & cov matrices).\n",
      "2023-11-20 12:49:27,395 - Network - INFO :    CE-40_C\n",
      "2023-11-20 12:49:27,456 - Network - INFO :    CE-20_S\n",
      "2023-11-20 12:49:27,515 - Network - INFO :    ET_ET1\n",
      "2023-11-20 12:49:27,584 - Network - INFO :    ET_ET2\n",
      "2023-11-20 12:49:27,663 - Network - INFO :    ET_ET3\n",
      "2023-11-20 12:49:27,761 - Network - INFO : Errors calculated.\n"
     ]
    }
   ],
   "source": [
    "i = 54837\n",
    "\n",
    "inj_params = {\n",
    "    'Mc':    Mcs[i],\n",
    "    'eta':   etas[i],\n",
    "    'chi1x': 0.,\n",
    "    'chi2x': 0.,\n",
    "    'chi1y': 0.,\n",
    "    'chi2y': 0.,\n",
    "    'chi1z': chi1z[i],\n",
    "    'chi2z': chi2z[i],\n",
    "    'DL':    DLs[i],\n",
    "    'tc':    0,\n",
    "    'phic':  0,\n",
    "    'iota':  iotas[i],\n",
    "    'ra':    ras[i],\n",
    "    'dec':   decs[i],\n",
    "    'psi':   psis[i],\n",
    "    'gmst0': 0,\n",
    "    'hybr': 0.0\n",
    "} \n",
    "\n",
    "net_key = '3G'\n",
    "\n",
    "approx1 = 'IMRPhenomXAS'\n",
    "approx2 = 'IMRPhenomD'\n",
    "\n",
    "d_f = 2**-4\n",
    "f_low = 5.0\n",
    "mtotal = conv.mtotal_from_mchirp_eta(inj_params['Mc'], inj_params['eta'])\n",
    "f_high = np.round(4*br.f_isco_Msolar(mtotal))\n",
    "\n",
    "net_ap = gwnet.get_hybrid_network_response(inj_params=inj_params, network_key=net_key, f_max=f_high,\n",
    "                            approximant1=approx1, approximant2=approx2, cond_num=1e25)\n",
    "\n",
    "# net_true = gwnet.get_network_response(inj_params=inj_params, f_max=f_high, approximant=approx1, network_key=net_key, calc_detector_responses=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  5.    ,   5.0625,   5.125 , ..., 194.8125, 194.875 , 194.9375])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "net_ap.hfc\n",
    "net_ap.f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numeric derivatives of detector responses calculated.\n"
     ]
    }
   ],
   "source": [
    "step=1e-9\n",
    "method='central' \n",
    "order=2\n",
    "n=1\n",
    "\n",
    "for det in net_ap.detectors:\n",
    "    det.calc_det_responses_derivs_num(net_ap.inj_params, net_ap.deriv_variables, net_ap.wf, net_ap.deriv_symbs_string, net_ap.conv_cos, net_ap.conv_log, net_ap.use_rot, step, method, order, n, user_locs=None)\n",
    "print('Numeric derivatives of detector responses calculated.')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "'numpy.ndarray' object is not callable",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m/Users/veomekapil/Library/CloudStorage/OneDrive-JohnsHopkins/Research/Waveform Systematics/hybrid_waveforms/quadratic_cv_bias.ipynb Cell 7\u001b[0m line \u001b[0;36m2\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/veomekapil/Library/CloudStorage/OneDrive-JohnsHopkins/Research/Waveform%20Systematics/hybrid_waveforms/quadratic_cv_bias.ipynb#X10sZmlsZQ%3D%3D?line=0'>1</a>\u001b[0m \u001b[39mfrom\u001b[39;00m \u001b[39mgwbench\u001b[39;00m \u001b[39mimport\u001b[39;00m wf_derivatives_num \u001b[39mas\u001b[39;00m wdn\n\u001b[0;32m----> <a href='vscode-notebook-cell:/Users/veomekapil/Library/CloudStorage/OneDrive-JohnsHopkins/Research/Waveform%20Systematics/hybrid_waveforms/quadratic_cv_bias.ipynb#X10sZmlsZQ%3D%3D?line=1'>2</a>\u001b[0m wdn\u001b[39m.\u001b[39;49mpart_deriv(net_ap\u001b[39m.\u001b[39;49mdetectors[\u001b[39m0\u001b[39;49m]\u001b[39m.\u001b[39;49mdel_hf[\u001b[39m'\u001b[39;49m\u001b[39mdel_Mc_hf\u001b[39;49m\u001b[39m'\u001b[39;49m], net_ap\u001b[39m.\u001b[39;49mf, net_ap\u001b[39m.\u001b[39;49mderiv_variables, pl_cr\u001b[39m=\u001b[39;49m\u001b[39m0\u001b[39;49m, compl\u001b[39m=\u001b[39;49m\u001b[39mNone\u001b[39;49;00m, step\u001b[39m=\u001b[39;49m\u001b[39m1e-9\u001b[39;49m, method\u001b[39m=\u001b[39;49m\u001b[39m'\u001b[39;49m\u001b[39mcentral\u001b[39;49m\u001b[39m'\u001b[39;49m, order\u001b[39m=\u001b[39;49m\u001b[39m2\u001b[39;49m, n\u001b[39m=\u001b[39;49m\u001b[39m1\u001b[39;49m)\n",
      "File \u001b[0;32m/opt/anaconda3/envs/gwbench_dev/lib/python3.9/site-packages/gwbench/wf_derivatives_num.py:122\u001b[0m, in \u001b[0;36mpart_deriv\u001b[0;34m(func, f, params_list, pl_cr, compl, step, method, order, n)\u001b[0m\n\u001b[1;32m    119\u001b[0m     \u001b[39mreturn\u001b[39;00m utils\u001b[39m.\u001b[39mz_deriv_from_amp_pha(amp, pha, del_amp, del_pha)\n\u001b[1;32m    121\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[0;32m--> 122\u001b[0m     \u001b[39mreturn\u001b[39;00m part_deriv_ndGradient(func, f, params_list, \u001b[39mNone\u001b[39;49;00m, step, method, order, n)\n",
      "File \u001b[0;32m/opt/anaconda3/envs/gwbench_dev/lib/python3.9/site-packages/gwbench/wf_derivatives_num.py:131\u001b[0m, in \u001b[0;36mpart_deriv_ndGradient\u001b[0;34m(func, f, params_list, funcid, step, method, order, n)\u001b[0m\n\u001b[1;32m    128\u001b[0m     \u001b[39melse\u001b[39;00m:              \u001b[39mreturn\u001b[39;00m func(f,\u001b[39m*\u001b[39mx)[funcid]\n\u001b[1;32m    130\u001b[0m \u001b[39mif\u001b[39;00m params_list \u001b[39m==\u001b[39m \u001b[39mNone\u001b[39;00m: \u001b[39mreturn\u001b[39;00m nd\u001b[39m.\u001b[39mGradient(wraps, step\u001b[39m=\u001b[39mstep, method\u001b[39m=\u001b[39mmethod, order\u001b[39m=\u001b[39morder, n\u001b[39m=\u001b[39mn)\n\u001b[0;32m--> 131\u001b[0m \u001b[39melse\u001b[39;00m:                   \u001b[39mreturn\u001b[39;00m nd\u001b[39m.\u001b[39;49mGradient(wraps, step\u001b[39m=\u001b[39;49mstep, method\u001b[39m=\u001b[39;49mmethod, order\u001b[39m=\u001b[39;49morder, n\u001b[39m=\u001b[39;49mn)(params_list)\n",
      "File \u001b[0;32m/opt/anaconda3/envs/gwbench_dev/lib/python3.9/site-packages/numdifftools/core.py:490\u001b[0m, in \u001b[0;36mGradient.__call__\u001b[0;34m(self, x, *args, **kwds)\u001b[0m\n\u001b[1;32m    489\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39m__call__\u001b[39m(\u001b[39mself\u001b[39m, x, \u001b[39m*\u001b[39margs, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwds):\n\u001b[0;32m--> 490\u001b[0m     result \u001b[39m=\u001b[39m \u001b[39msuper\u001b[39;49m(Gradient, \u001b[39mself\u001b[39;49m)\u001b[39m.\u001b[39;49m\u001b[39m__call__\u001b[39;49m(np\u001b[39m.\u001b[39;49matleast_1d(x)\u001b[39m.\u001b[39;49mravel(), \u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwds)\n\u001b[1;32m    491\u001b[0m     \u001b[39mif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mfull_output:\n\u001b[1;32m    492\u001b[0m         \u001b[39mreturn\u001b[39;00m result[\u001b[39m0\u001b[39m]\u001b[39m.\u001b[39msqueeze(), result[\u001b[39m1\u001b[39m]\n",
      "File \u001b[0;32m/opt/anaconda3/envs/gwbench_dev/lib/python3.9/site-packages/numdifftools/core.py:431\u001b[0m, in \u001b[0;36mJacobian.__call__\u001b[0;34m(self, x, *args, **kwds)\u001b[0m\n\u001b[1;32m    430\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39m__call__\u001b[39m(\u001b[39mself\u001b[39m, x, \u001b[39m*\u001b[39margs, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwds):\n\u001b[0;32m--> 431\u001b[0m     \u001b[39mreturn\u001b[39;00m \u001b[39msuper\u001b[39;49m(Jacobian, \u001b[39mself\u001b[39;49m)\u001b[39m.\u001b[39;49m\u001b[39m__call__\u001b[39;49m(np\u001b[39m.\u001b[39;49matleast_1d(x), \u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwds)\n",
      "File \u001b[0;32m/opt/anaconda3/envs/gwbench_dev/lib/python3.9/site-packages/numdifftools/core.py:288\u001b[0m, in \u001b[0;36mDerivative.__call__\u001b[0;34m(self, x, *args, **kwds)\u001b[0m\n\u001b[1;32m    286\u001b[0m x_i \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39masarray(x)\n\u001b[1;32m    287\u001b[0m \u001b[39mwith\u001b[39;00m np\u001b[39m.\u001b[39merrstate(divide\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mignore\u001b[39m\u001b[39m'\u001b[39m, invalid\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mignore\u001b[39m\u001b[39m'\u001b[39m):\n\u001b[0;32m--> 288\u001b[0m     results, f_xi \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_derivative(x_i, args, kwds)\n\u001b[1;32m    289\u001b[0m     derivative, info \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_extrapolate(\u001b[39m*\u001b[39mresults)\n\u001b[1;32m    290\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mfull_output:\n",
      "File \u001b[0;32m/opt/anaconda3/envs/gwbench_dev/lib/python3.9/site-packages/numdifftools/core.py:422\u001b[0m, in \u001b[0;36mJacobian._derivative_nonzero_order\u001b[0;34m(self, x_i, args, kwds)\u001b[0m\n\u001b[1;32m    420\u001b[0m diff, f \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_get_functions(args, kwds)\n\u001b[1;32m    421\u001b[0m steps, step_ratio \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_get_steps(x_i)\n\u001b[0;32m--> 422\u001b[0m fxi \u001b[39m=\u001b[39m f(x_i)\n\u001b[1;32m    423\u001b[0m results \u001b[39m=\u001b[39m [diff(f, fxi, x_i, h) \u001b[39mfor\u001b[39;00m h \u001b[39min\u001b[39;00m steps]\n\u001b[1;32m    425\u001b[0m steps2 \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_expand_steps(steps, x_i, fxi)\n",
      "File \u001b[0;32m/opt/anaconda3/envs/gwbench_dev/lib/python3.9/site-packages/numdifftools/core.py:257\u001b[0m, in \u001b[0;36mDerivative._get_functions.<locals>.export_fun\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m    256\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mexport_fun\u001b[39m(x):\n\u001b[0;32m--> 257\u001b[0m     \u001b[39mreturn\u001b[39;00m fun(x, \u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwds)\n",
      "File \u001b[0;32m/opt/anaconda3/envs/gwbench_dev/lib/python3.9/site-packages/gwbench/wf_derivatives_num.py:127\u001b[0m, in \u001b[0;36mpart_deriv_ndGradient.<locals>.wraps\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m    126\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mwraps\u001b[39m(x):\n\u001b[0;32m--> 127\u001b[0m     \u001b[39mif\u001b[39;00m funcid \u001b[39m==\u001b[39m \u001b[39mNone\u001b[39;00m: \u001b[39mreturn\u001b[39;00m func(f,\u001b[39m*\u001b[39;49mx)\n\u001b[1;32m    128\u001b[0m     \u001b[39melse\u001b[39;00m:              \u001b[39mreturn\u001b[39;00m func(f,\u001b[39m*\u001b[39mx)[funcid]\n",
      "\u001b[0;31mTypeError\u001b[0m: 'numpy.ndarray' object is not callable"
     ]
    }
   ],
   "source": [
    "from gwbench import wf_derivatives_num as wdn\n",
    "wdn.part_deriv(net_ap.detectors[0].del_hf['del_Mc_hf'], net_ap.f, net_ap.deriv_variables, pl_cr=0, compl=None, step=1e-9, method='central', order=2, n=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gwbench_dev",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}